Cosmological Perturbations of Extreme Axion in the Radiation Era

TL;DR

This paper investigates sub-horizon perturbations in the extreme axion model during the radiation era, revealing a spectral excess that could significantly impact early universe structures and galaxy formation.

Contribution

It introduces a novel analysis of extreme axion perturbations, showing they can cause a spectral excess and are described by Mathieu's equation, highlighting features absent in free-particle models.

Findings

Spectral excess in certain wave number ranges compared to CDM.

Perturbations described by Mathieu's equation with parametric instability.

Potential impact on high-redshift universe and early galaxy formation.

Abstract

Sub-horizon perturbations under the extreme initial condition of the axion model are investigated, where initial axion angles start near the potential maximum. This work focuses on a few new features found in the extreme axion model but absent in the free-particle model. A particularly novel new feature is the spectral excess relative to the CDM model in some wave number range, where the excess may be so large that landscapes of high-redshift universe beyond $z=10$ can be significantly altered. For axions of particle mass $10^{-22}$ eV, this range of wave number corresponds to first galaxies of few times $10^9-10^{10} M_\odot$. We demonstrate that sub-horizon perturbations are accurately described by Mathieu's equation and subject to parametric instability, which explains this novel feature. Actually the axion model is not a special one; perturbations in a wide range of scalar field models can share the similar characteristic.

Equations

$${\theta }^{{}^{\prime \prime }}+2H{\theta }^{{}^{\prime }}+m^2{a}^2\sin (\theta )=0,$$

$$\delta {\theta }^{{}^{\prime \prime }}+2H\delta {\theta }^{{}^{\prime }}+[k^2+m^2{a}^2\cos (\theta )]\delta \theta =4{\varphi }^{{}^{\prime }}{\theta }^{{}^{\prime }}-2m^2{a}^2\sin (\theta )\varphi ,$$

$${\Delta }_{\theta }\equiv \frac{{\theta }^{{}^{\prime }}\delta {\theta }^{{}^{\prime }}+m^2{a}^2\sin (\theta )\delta \theta -({\theta }^{{}^{\prime }}{)}^2\varphi +3H{\theta }^{{}^{\prime }}\delta \theta }{(1/2)({\theta }^{{}^{\prime }}{)}^2+m^2{a}^2[1-\cos (\theta )]},$$

$$\ddot{Q}+{\omega }_{\theta }^2[1+\eta +ϵ\cos (2{\omega }_{\theta }t)]Q=0,$$

$$\left.\begin{aligned} &\ddot{q}+\Big{\{}{3\over 4}\mathscr{H}^{2}+{{k^{2}}\over{a^{2}}}+m^{2}\Big{[}1-{{\langle\theta^{2}\rangle}\over 4}\Big{(}1+\cos\Big{(}2\int_{t_{m}}^{t}\omega_{\theta}dt\Big{)}\Big{)}\Big{]}\Big{\}}q\\ &={Sr\over{a^{{1\over 2}}}},\end{aligned}\right.$$

$$\omega_{\theta}^{2}=m^{2}\Big{[}1-{{\theta_{0}^{2}}\over 8}\Big{(}{t\over{t_{m}}}\Big{)}^{-{3\over 2}}\Big{]}.$$

$$\omega_{\delta}^{2}=\omega_{\theta}^{2}\Big{[}1-{{\theta_{0}^{2}}\over 8}\Big{(}{t\over{t_{m}}}\Big{)}^{-{3\over 2}}\Big{]},$$

$$\left.\begin{aligned} &\ddot{q}+\Big{\{}{{k^{2}}\over{a^{2}}}+\omega_{\theta}^{2}\Big{[}1-{{\theta_{0}^{2}}\over 8}\Big{(}{t\over{t_{m}}}\Big{)}^{-{3\over 2}}\Big{(}1+2\beta\cos\Big{(}2\int_{t_{m}}^{t}\omega_{\theta}dt\Big{)}\Big{)}\Big{]}\Big{\}}q\\ &=0,\end{aligned}\right.$$

$$\left.\begin{aligned} &\ddot{q}+\Big{\{}{{k^{2}}\over{a^{2}}}+\omega_{\theta}^{2}\Big{[}1-{\theta_{0}^{2}\over 8}\Big{(}{t\over{t_{m}}}\Big{)}^{-{3\over 2}}\Big{(}1+2\beta\cos(2\omega_{\theta}(t-t_{m}))\Big{)}\Big{]}\Big{\}}q\\ &=0,\end{aligned}\right.$$

$$i\dot{\hat{q}}={1\over{2\omega_{\theta}}}\Big{[}\Big{(}\Big{(}{k\over a}\Big{)}^{2}-\alpha\Big{)}\hat{q}-\beta\alpha\hat{q}^{*}\Big{]},$$

$$\omega={k\over 2\omega_{\theta}a}\Big{[}\Big{(}{k\over a}\Big{)}^{2}-2\alpha\Big{]}^{1/2},$$

$$\left.\begin{aligned} A(a_{\kappa})={{2k^{2}}\over{k_{c}^{2}}}\Big{\{}&\Big{[}{{a_{m0}}\over{a_{\kappa}}}\Big{(}{{a_{m}}\over{a_{m0}}}\Big{)}^{3}{{k_{c}^{2}}\over{k^{2}}}{{\pi^{2}}\over 4}-1\Big{]}^{1\over 2}-\\ &\tan^{-1}\Big{[}\Big{(}{{a_{m0}}\over{a_{\kappa}}}\Big{(}{{a_{m}}\over{a_{m0}}}\Big{)}^{3}{{k_{c}^{2}}\over{k^{2}}}{{\pi^{2}}\over 4}-1\Big{)}^{1\over 2}\Big{]}\Big{\}},\end{aligned}\right.$$

$$\left.\begin{aligned} B(\eta)=2{{k^{2}}\over{k_{c}^{2}}}\Big{\{}&-[1-r\exp(-\eta)]^{1/2}\\ &+{1\over 2}\ln\Big{[}{{1+[1-r\exp(-\eta)]^{1\over 2}}\over{1-[1-r\exp(-\eta)]^{1\over 2}}}\Big{]}\Big{\}}.\end{aligned}\right.$$

$$\left.\begin{aligned} \eta_{peak}=&{\pi\over 4}\Big{(}{{k_{c}}\over k}\Big{)}^{2}+2(1-\ln(2))+\\ &\Big{[}2\ln\Big{(}{{k_{c}}\over{k}}\Big{)}+3\ln\Big{(}{{a_{m}}\over{a_{m0}}}\Big{)}+2\ln\Big{(}{\pi\over 2}\Big{)}\Big{]},\end{aligned}\right.$$

$$\left.\begin{aligned} \Delta\eta_{peak}\approx&-{\pi\over 4}\Big{(}{{k_{c}}\over k}\Big{)}^{2}+(1+\kappa)\ln\Big{(}{{k_{c}}\over{k}}\Big{)}+3\ln\Big{(}{{a_{m}}\over{a_{m0}}}\Big{)}\\ &+2+2\ln\Big{(}{\pi\over 4}\Big{)}+\ln\Big{(}{2\over\sqrt{3}}\Big{)},\end{aligned}\right.$$

$$\Big{(}{{a_{m}}\over{a_{m0}}}\Big{)}^{2}\approx 3.5-{2\over 3}\ln\Big{(}{{\delta\theta_{0}}\over{\delta\theta_{01}}}\Big{)},$$

$$i{{\partial}\over{\partial t}}\Phi=\Big{[}-{{\nabla^{2}}\over{2m}}+g|\Phi|^{2}\Big{]}\Phi.$$

$$i\dot{\delta}p+g|\Phi_{0}|^{2}\delta p=\Big{[}{{k^{2}}\over{2m}}+2g|\Phi_{0}|^{2}\Big{]}\delta p+g|\Phi_{0}|^{2}\delta p^{*}.$$

$$\left.\begin{aligned} &{{d^{2}}\over{d\eta^{2}}}\Big{(}{{\delta n}\over n}\Big{)}+\Big{(}{k\over{k_{c}}}\Big{)}^{4}\Big{[}1-\Big{(}{{k_{c}}\over k}\Big{)}^{2}{{\theta_{0}^{2}}\over 4}\Big{(}{{a_{m}}\over{a_{m0}}}\Big{)}^{2}\exp(-\eta)\Big{]}\Big{(}{{\delta n}\over n}\Big{)}\\ &=0,\end{aligned}\right.$$

$$\left.\begin{aligned} &A(\eta)=\\ &2{{k^{2}}\over{k_{c}^{2}}}\Big{\{}[r\exp(-\eta)-1]^{1\over 2}-\tan^{-1}[((r\exp(-\eta)-1)^{1\over 2}]\Big{\}},\\ &B(\eta)=\\ &2{{k^{2}}\over{k_{c}^{2}}}\Big{\{}-[1-r\exp(-\eta)]^{1\over 2}+{1\over 2}\ln\Big{[}{{1+[1-r\exp(-\eta)]^{1\over 2}}\over{1-[1-r\exp(-\eta)]^{1\over 2}}}\Big{]}\Big{\}},\end{aligned}\right.$$

$$\left.\begin{aligned} &{\pi\over 4}\Big{(}{{k_{c}}\over{k_{peak}}}\Big{)}^{2}+\ln\Big{[}{\pi\over 4}\Big{(}{{k_{c}}\over{k_{peak}}}\Big{)}^{2}\Big{]}\\ &=\ln\Big{(}{{a_{eq}}\over{a_{m0}}}\Big{)}-3\ln\Big{(}{{a_{m}}\over{a_{m0}}}\Big{)}-2-\ln\Big{(}{\pi\over 4}\Big{)},\end{aligned}\right.$$

Full text

I Introduction

Scalar fields as dark matter candidates have a long history of development Dehnen and Rose (1993); Sin (1994); Ji and Sin (1994); Dehnen et al. (1995); Lee and Koh (1996); Hu et al. (2000); Guzmán and Matos (2000); Matos et al. (2000); Matos and Ureña-López (2000); Sahni and Wang (2000), but most works were addressing the matter-dominated era where comparisons with observations can be made. While these models have some degrees of freedom to accommodate a suite of observational oddities, it is inevitable that they must introduce one or more energy or mass scales, in sharp contrast to the CDM model, which is extremely insensitive to the particle mass. For example, wave dark matter ($\psi$DM)Schive et al. (2014), ultra-light bosonic dark matter Woo, and Chiueh (2009) or fuzzy dark matter Hu et al. (2000) introduces one boson mass $m$, which is found to be around $10^{-22}$ eV to explain kilo-parsec-scale core structures in dwarf spheroidal galaxies Schive et al. (2014); Chen et al. (2017). If $m\gg 10^{-22}$ eV, $\psi$ dark matter becomes indistinguishable from CDM observationally. The axion model, being a nonlinear field model, introduces a second energy scale $f$, in addition to $m$, where $f$ is the axion decay constant that is above the GUT scale to explain the cosmic background dark matter mass density to be so close to the cosmic critical density in a non-QCD axion model involving the dark sectorChiueh (2014); Davoudiasl, and Murphy (2017). Indeed, recent developments of string theories also favor extremely light axions with a large axion decay constant $f$ much greater than the electroweak scale Svrček, and Witten (2006); Arvanitaki et al. (2010); Hui et al. (2017); Diez-Tejedor, and Marsh (2017).

However, there is one more free degree of freedom, i.e., the initial field amplitude, which is a dimensionless parameter not present in the field Lagrangian but is able to control the solution. Whether the initial field is located in a linear regime or in a nonlinear regime may make a difference in the solution space and affects the observable. In the context of cosmology, as the universe expands the field amplitude quickly decreases due to Hubble friction, and soon the field samples only the quadratic part of the potential to become free particle. Hence the free-particle model ($\psi$DM) is the ultimate time asymptotic attractor for the axion model and for many other nonlinear scalar field models. One therefore hopes that the Hubble friction may erase the memory of the initial condition, and the solution converges to the free-particle solution.

We therefore studied linear perturbations of the free-particle model in the radiation dominant era in the previous work Zhang, and Chiueh (2017) (Paper (@slowromancapi@)). Four phases of evolution are identified. Central to the four phases is the critical wavenumber $k_{c}$, for which the mode enters horizon when the horizon size equals the Compton wavelength. This critical wavenumber lies at the boundary of the four phases and gives rise to a sharp spectral transition. We have also numerically investigated perturbations of the axion model to investigate the attractor aspect of the problem, and indeed found that the time-asymptotic solution depends very weakly on the initial angles, except when the axion field starts from very close to the top of the field potential, a highly nonlinear initial field. In such an exceptional case, the perturbation begins to behave quite unexpectedly from when the field starts elsewhere. We call this singular case the extreme axion model. This narrow window of new degree of freedom is interesting, and may allows for accommodating the tension concerning the particle mass of $\psi$DM determined by the high-redshift Lyman-$\alpha$ forests Armengaud et al. (2017); Iršič et al. (2017) and by the flat cores of nearby dwarf spheroid galaxies Chen et al. (2017); Lora and Magaña (2014); Calabrese and Spergel (2016). In this work, we follow up this finding of Paper (@slowromancapi@) for the extreme axion model and analyze the perturbation evolution in details. Particular emphasis is placed on sub-horizon modes after the onset of mass oscillation, as it holds the key to the unexpected. The analysis developed in this work can be extended to other nonlinear scalar field models with a finite potential barrier.

For the fiducial boson mass $m$ as small as $10^{-22}$ eV, the particle number density is extremely high, yielding a critical temperature so high that any conceivable background temperature is way below the critical. These bosons therefore form a Bose-Einstein condensate and many-particle wave functions collapse to a single wave function. To acquire phase coherence for many-particle wave functions, nonlinearity is essential to couple these wave functions and locks their phases. The nonlinearity of the scalar field for a Bose-Einstein condensate is just a manifestation of the microscopic two-body scattering that correlate wave functions. To the leading order of interactions, the simplest S-wave scattering of a dilute boson gas results in the Gross-Pitaevskii equation as an effective macroscopic theory Gross (1961); Pitaevskii (1961) via the well-known Bogoliubov’s reduction formalism in the non-relativistic limit Bogoliubov (1947). We will show that a general class of relativistic scalar field models can be reduced to the Gross-Pitaevskii equation to the leading-order nonlinearity in the non-relativistic limit, and the axion model is no exception. Therefore the axion model gains a microscopic support for being a Bose-Einstein condensate, and this finding will be elaborated later.

Recently there have been numerical works attempting to address axion perturbations Linares Cedeño et al. (2017), partly motivated by string theories Svrček, and Witten (2006); Arvanitaki et al. (2010); Diez-Tejedor, and Marsh (2017) and partly by the emerging interest in wave dark matterSchive et al. (2014); Chen et al. (2017); Armengaud et al. (2017); Iršič et al. (2017); Hui et al. (2017); Marsh and Silk (2014); Schive et al. (2014); Marsh and Pop (2015); Marsh (2016). The difficulty of computing axion perturbations arises from that the equation demands high numerical accuracy to solve, as it must stably track two near-by-frequency oscillations for thousands of, or even much more, periods to determine the precise relative phase shift between the two oscillations. It is therefore essential that numerical results have a support from detailed analyses of the solution. One of the aims of this work is to serve for such a purpose. We find for the free-particle case, as in Paper (@slowromancapi@), all numerical works largely agree. For the extreme axion case, we find the numerical results of other works Diez-Tejedor, and Marsh (2017); Linares Cedeño et al. (2017) deviate from the result of this work. It remains to be seen at what numerical bottleneck these differences arise.

This paper is structured as follows. In Sec. (II), we briefly review the previous work and pose problems pertinent to the three unexpected features to be understood. Sec. (III) addresses the first feature. The remaining two features require a new mechanism involving the parametric drive and amplification, and we elaborate this mechanism in Sec. (IV). The treatment can generalized to other scalar field models as shown in Sec. (V). The matter power spectrum of the extreme axion model is also shown in this section. In Sec. (VI), we make contact of the axion model to the Gross-Pitaevskii equation and also discuss the concern about quantum tunneling when the initial field assumes a classically unstable value. We conclude this work in Sec. (VII). The particle mass dependence of our results is discussed in Appendix. This work is confined to the radiation-dominant era unless otherwise mentioned, the fiducial boson mass is chosen $10^{-22}$ eV, and standard cosmological parameters of the concordance model are adopted, i.e., $H_{0}=70km/sec/Mpc$, $\Omega_{DM}=0.24$, $\Omega_{b}=0.06$. We also set the speed of light $c$ and the Planck constant $\hbar$ equal to $1$. Throughout the analysis we adopt the Newtonian gauge for perturbations.

II Review on Unique Features of the Extreme Axion Model

In the following, we shall first adopt the axion model with the $m^{2}a^{2}(1-\cos(\theta))$ field potential as a working template, and in a later section we will extend the analysis to more general nonlinear potentials. The equations of motion for the axion background field $\theta$ and the perturbed field $\delta\theta$ are respectively

[TABLE]

and

[TABLE]

where the prime denotes the derivative with respect to the conformal time $\tau(\equiv\int_{0}dt/a\propto a)$, $a$ is the scaling factor, $H(\equiv da/dt\propto a^{-1})$ the conformal Hubble parameter, and $k$ the wavenumber. The right-hand side of the perturbed field equation is the source $Sr$, which contains the metric perturbation $\phi$ and is contributed from perturbations of all species through the Poisson equation. The metric perturbation is dominated by the photon perturbation in the radiation era, and to a good approximation we can regarded $\phi$ as independent of $\delta\theta$ till near radiation-matter equality. This approximation is called the passive evolution in Paper (@slowromancapi@). The full treatment of $\phi$ must include additionally the dark matter, baryons and neutrinos, as elucidated in Paper (@slowromancapi@), which shows that passive evolution provides a good approximation before the epoch of radiation-matter equality.

We need also define one more quantity, the dimensionless gauge-covariant energy density of the perturbed field

[TABLE]

which is the normalized physical energy density perturbation and the denominator is the background field denisty $\epsilon_{\theta}$.

A short summary of Paper (@slowromancapi@) is in order. In the small $\theta$ limit, which most parts of Paper (@slowromancapi@) addresses, is the free-particle limit as the field nonlinearity vanishes. In this limit, we have well-defined four phases, (a) before mass oscillation and superhorizon, (b) before mass oscillation and subhorizon, (c) after mass oscillation and super-horizon, and (d) after mass oscillation and subhorizon. Long wave go through phases (a),(c) and (d), but short waves through (a),(b) and (d). The division of long and short waves is the critical wavenumber $k=k_{c}$, where $k_{c}=ma=2H$, for which all terms on the left-hand side of Eq. (2) are equally important and mass oscillation just begins to set in. In phase (a), $\Delta_{\theta}$ grows as $a^{6}$, phase (b) as $a^{2}\cos((ka/\sqrt{3})-\chi_{1})$, phase (c) as $a^{2}$ and phase (d) as $\sin((k^{2}/2maH)\ln(a)-\chi_{2})$, where $\chi_{1}$ and $\chi_{2}$ are oscillation phases associated with photon and matter wave, respectively.

When the initial angle $\theta_{0}$ is not small, the $\theta_{0}=\pi/2$ case is of very little difference from the free-particle model found in Paper (@slowromancapi@). In fact, Fig. (1) sums up nicely the features to be discussed, which illustrates the evolution of $\Delta_{\theta}(a)$ for long, medium and short waves with initial angles $\theta_{0}=179^{o}$ and $90^{o}$. When comparing $\theta_{0}=179^{o}$ and $90^{o}$ cases, the first feature common to all wave numbers is (1) a steep rise in amplitude at the onset of mass oscillation for $\theta_{0}=179^{o}$ not present for $\theta_{0}=90^{o}$. A second feature of the $\theta_{0}=179^{o}$ case is (2) a substantially longer duration of the first half cycle of matter-wave oscillation for some $k\sim k_{c}$ than that of the $\theta_{0}=\pi/2$ case for the same $k$. Associated with the second feature is a third feature that (3) the perturbation amplitude of the $\theta_{0}=179^{o}$ case is higher than that of the CDM model during a certain period also for some $k\sim k_{c}$, which has never been observed in the free-particle model. Clearly, these three features are not caused by ordinary nonlinear mass oscillation of $\theta$, but associated with the extreme condition where $\theta_{0}\to\pi$. Note that the second and third features do not show up prominently for $k\ll k_{c}$ and $k\gg k_{c}$ in Fig. (1). This requires an explanation.

Finally, Fig. (1) demonstrates that the passive evolution approximates the evolution of full treatment quite well till near the radiation-matter equality. The focus of our analysis in this work is placed upon after the onset of mass oscillation but still far away from the epoch of radiation-matter equality. Hence, passive evolution provides a fair simplification for understanding the above three features; however, our numerical solutions will include the full treatment.

III Abrupt Growth of $\delta\theta$

Prior to the onset of mass oscillation, the perturbation grows as $a^{6}$, corresponding to the earliest phase (a) in the evolution. When $\theta_{0}$ is near the top of the field potential, it delays the mass oscillation, and substantial delay makes the friction $2H\delta\theta^{{}^{\prime}}$ negligible in Eq. (2) at the onset of mass oscillation. This creates an almost frictionless background for perturbed field dynamics. This rapid growth occurs only in a short time when $\theta$ first rolls down from the potential top. The duration of exponential growth is independent of the exact location of the initial angle $\theta_{0}$ from the potential top as long as $\theta_{0}$ is close to the top. More importantly the abrupt growth is insensitive to wave number $k$, as evidenced from the same abrupt growth for $k=0.1k_{c}$ mode and $k=10k_{c}$ mode in the case of $\theta_{0}=179^{o}$ in Fig. (1). This provides a crucial clue for the growth mechanism.

The restoring forces of long-wave modes and short-wave modes are very different, with the former being negative and the latter being positive. Hence the same growth for all $k$ modes indicates that the cause of the growth should be from the source, the right-hand side of Eq. (2). Unlike the free-particle model, well before the onset of nonlinear mass oscillation the source is almost zero, where $\theta^{\prime}\to 0$ and $m^{2}a^{2}\sin(\theta)\to 0$ as $\theta\to\pi$. The weak constant source yields a small coefficient in the $a^{6}$ growth (phase (a)) before the abrupt growth, more so for $\theta_{0}$ closer to the potential top, as opposed to a much larger coefficient due to a much larger source in the $\theta_{0}=\pi/2$ case. Just at the onset the source suddenly rises to its full strength when the field $\theta$ rolls down the hill on its first pass. Such a drive is so abrupt that the perturbed field gets amplified regardless of the nature of its restoring force, since the restoring force has no time to respond. One may analogize this mechanism as the ”direct current (DC)” drive, as opposed to the ”alternative current (AC)” drive of the parametric instability to be discussed in the next section. After this short period of time, the source strength either stays in full strength or declines depending on whether the mode has entered horizon. For super-horizon modes, the source stays in full strength and the modes enter phase (c) of a slow $a^{2}$ growth, and for sub-horizon modes, they enter a new phase of parametric instability or matter-wave oscillation, phase (d).

The exact location of $\theta_{0}$ from the top would, however, affect the onset time of mass oscillation. For a given $k$, the more delay of the mass oscillation, the longer the duration of the $a^{6}$ growth, and the perturbation can grow to a greater amplitude. On the other hand, the closer $\theta_{0}$ is to the field potential top, the smaller is the source, and the smaller the coefficient of the $a^{6}$ growth as mentioned above. These two opposite trends almost cancel, and by the end of the abrupt growth, $\delta\theta$ is brought to nearly the same amplitude as the free-particle model. An alternative way to understand this is that once the source becomes at its full strength, it drives the perturbed field to a level comparable to the photon perturbation before the perturbed field becomes decoupling from the source shortly after the onset of mass oscillation. Such a driving mechanism applies to all adiabatic perturbations.

IV Parametric Instability

Parametric instability refers to the presence of an oscillating restoring force of almost twice the natural frequency for an oscillator, described by Mathieu’s equation Abramowitz and Stegun (1964):

[TABLE]

where $Q$ is the oscillator solution and the overdot denotes $d/dt$. The parameters $\epsilon$ and $\eta$ are the driver strength and the detuning squared frequency. The $(\eta,\epsilon)$ phase diagram at small $\epsilon$ and $\eta$ marks the marginally stability curve as $|\eta|=(1/2)|\epsilon|$. In the limit $\eta=0$, the oscillator is unconditionally unstable even for an tiny but finite $\epsilon$.

To make a comparison with Mathieu’s equation, we change the variable of Eq. (2) from the scaling factor $a$ to the ordinary time $t$. A straightforward algebra shows that Eq. (2), Taylor-expanded up to the first-order nonlinearity in $\theta$, can be cast into the equation:

[TABLE]

where $q\equiv\delta\theta/\sqrt{\epsilon_{\theta}}$ with $\epsilon_{\theta}$ being the background energy density, $\mathscr{H}=H/a$, $t_{m}$ is the time for the onset of mass oscillation, $\omega_{\theta}$ the frequency of $\theta$ containing a nonlinear frequency shift, and $Sr$ is the right-hand side of Eq. (2). The short-time average $\langle\theta^{2}\rangle$ decays as $a^{-3}$, which we model as $\langle\theta^{2}\rangle=(\theta_{0}^{2}/2)(t/t_{m})^{-3/2}$.

The frequency $\omega_{\theta}$ of the nonlinear oscillation of $\theta$ can be derived from Eq. (1), where the restoring force $\sin(\theta)=\theta(1-\theta^{2}/6)+O(\theta^{5})\sim\theta_{0}(t/t_{m})^{-3/4}[1-(t/t_{m})^{-3/2}(\theta_{0}^{2}/8)]\cos(\int\omega_{\theta}dt))-[\theta_{0}^{3}(t/t_{m})^{-9/4}/24]\cos(3\int\omega_{\theta}dt)$, assuming $\theta=\theta_{0}(t/t_{m})^{-3/4}\cos(\int\omega_{\theta}dt)$. Ignoring the triple frequency term and retaining the coefficient of $\cos(\int\omega_{\theta}dt)$, we have the driving frequency

[TABLE]

On the other hand, the perturbed field $q$ has a natural frequency $\omega_{\delta}$ different from the driving frequecy $\omega_{\theta}$ and related by

[TABLE]

to the leading order.

We shall address the sub-horizon regime where $k/a\gg\mathscr{H}$. Hence we can ignore both the weak source term $Sr$ as the driver $\phi$ declines as $a^{-2}$ shown in Paper (@slowromancapi@) and the $\mathscr{H}^{2}$ term in Eq. (5), thus arriving at a simplified equation that describes the homogeneous solution of $q$,

[TABLE]

with $\beta=1$. An addtional parameter $\beta$ is introduced so as to make a close contact with the Mathieu’s equation which has two parameters $\epsilon$ and $\eta$.

Now, Eq. (8) is the Mathieu’s equation with time-dependent coefficients, where the detuning squared frequency $\eta=(k^{2}/a^{2})-\delta\omega^{2}$ and the driver strength $\epsilon=2\delta\omega^{2}$ with $\delta\omega^{2}\equiv\omega_{\theta}^{2}(\theta_{0}^{2}/8)(t/t_{m})^{-3/2}$. Note that when $(k/a)^{2}\to 0$, we have $|\eta|=(1/2)|\epsilon|$ and it satisfies the marginally stable condition of Mathieu’s equation. Worth noting is that the squared detuning frequency $\eta$ has zero-crossing for a range of $k$ and $a$, and these $k$-modes can temporarily be parametrically unstable. This provides a crude explanation as to why in some range of $k$ the matter-wave oscillation appears to be amplified and has a relatively high amplitude, i.e., feature (3), but more details will follow.

Other than the aforementioned growth due to the parametric drive, the frequency of the solution of Eq. (8) actually deviates from its natural frequency and is locked to near half of the driving frequency for some period; therefore the solution becomes nearly phase locked to the driver during this period. As shown in Paper (@slowromancapi@), we may let $q=\Re[\hat{q}\exp(-i\int\omega_{\theta}dt)]=\Re[\hat{q}]\cos(\int\omega_{\theta}dt)+\Im[\hat{q}]\sin(\int\omega_{\theta}dt)]$ while the background field $\theta=\sqrt{\epsilon_{\theta}}\cos(\int\omega_{\theta}dt)$, where $\hat{q}$ is a slowly varying complex amplitude. When $q$ and $\theta$ are phase locked, the amplitudes of $\Re[\hat{q}]$ and $\Im[\hat{q}]$ will remain fixed and do not oscillate until the nonlinearity dies out, after which the perturbation assumes free-particle matter-wave oscillation. This picture provides a rough baseline as to why the first half cycle of matter-wave oscillation in $\Delta_{\theta}$ has a long duration. Again, more details are to come.

We shall analyze an even more simplified version of Eq. (8) below, which bears more resemblance to Mathieu’s equation, in order to bring out the aforementioned frequency locking and the amplitude excess in a quantitative manner. We assume the background field oscillates at a fixed frequency, $\omega_{\theta}=m$, ignoring the nonlinear contribution to the driving frequency which is a high-order effect for our purpose. Equation (8) thus becomes

[TABLE]

Using this $\hat{q}$ representation for sub-horizon modes after mass oscillation, one can show that the normalized energy density $\Delta_{\theta}\approx 2\Re[\hat{q}]$. Aside from the coefficient, the interaction ($\beta$) term in Eq. (9) yields $(\Re[\hat{q}]/2)[\cos(\omega_{\theta}(t-t_{m}))+\cos(3\omega_{\theta}(t-t_{m}))]-(\Im[\hat{q}]/2)[\sin(\omega_{\theta}(t-t_{m}))-\sin(3\omega_{\theta}(t-t_{m}))]$. Again ignoring the triple frequency contribution, the interaction term is then proportional to $(1/2)[\Re[\hat{q}]\cos(\omega_{\theta}(t-t_{m}))-\Im[\hat{q}]\sin(\omega_{\theta}(t-t_{m}))]=(1/2)\Re[\hat{q}^{*}\exp(-i\omega_{\theta}(t-t_{m}))]$. Substituting this result into Eq. (9), we have a reduced perturbation equation satisfying

[TABLE]

where $\alpha=\omega_{\theta}^{2}(\theta_{0}^{2}/8)(a_{m}/a)^{3}$, and $a_{m}$ is the scaling factor at the onset of nonlinear mass oscillation. Separating the real and imaginary parts of $\hat{q}$, one can straightforwardly show that the dispersion relation for this equation is $\omega=(1/2\omega_{\theta})[[(k/a)^{2}-\alpha]^{2}-(\beta\alpha)^{2}]^{1/2}$ with $\omega$ being the matter-wave frequency. This dispersion relation yields the characteristics of parametric instability. For the $\beta=1$ axion case, the dispersion relation becomes

[TABLE]

and the mode is unstable when $k^{2}/2a^{2}\alpha<1$, and stable with $\omega\to k^{2}/2\omega_{\theta}a^{2}$ when $k^{2}/2a^{2}\alpha\gg 1$. This dispersion relation is valid even when $\omega_{\theta}\neq m$ where Eq. (6) holds. For simplicity we shall continue to ignore the nonlinear correction to the driving frequency and assume $\omega_{\theta}=m$.

We first note the factor $k^{2}/(2a^{2}\alpha)\propto k^{2}(a/a_{m})(a_{m})^{-2}$, the greater $a_{m}$ of nonlinear mass oscillation is, or the closer $\theta_{0}$ is to $\pi$, the smaller the magnitude of this factor at a given $a/a_{m}$, thus mimicking a smaller $k$ for the free-particle model that has a longer matter-wave oscillation period and accounts for feature (2). On the other hand, for a given $\theta_{0}$ near $\pi$ and in the limit $k\to 0$, the parametric instability is weak, as by the time when the mode enters horizon where Eq. (11) becomes valid, the nonlinearity is already small. So the only range of $k$ exhibiting a strong parametric growth is when $k$ is on the same order of $k_{c}$. This explains why the amplitude excess occurs for $k$ on the same order of $k_{c}$, feature (3).

To put the above into quantitative perspectives, one readily sees from the dispersion relation, Eq. (11), that the frequency’s being either imaginary and small compared with the free-particle frequency is to contribute to higher perturbation amplitudes and a longer duration in the first half cycle of matter-wave oscillation. The unstable phase takes place during a time interval $\ln(a_{0}/a_{m0})=2\ln(\theta_{0}/2)+3\ln(a_{m}/a_{m0})+2\ln(k_{c}/k)]$ or $\ln(a_{0}/a_{k})=2\ln(\theta_{0}/2)+3\ln(a_{m}/a_{m0})+\ln(k_{c}/k)$ depending on whether the mode has entered horizon or not, respectively, at the onset of nonlinear mass oscillation, where $a_{0}$ is the scaling factor at the end of growth $\omega=0$, and $a_{k}$ that at the horizon entry $k=2H$. Here $a_{m0}$ corresponds to $a$ at the onset of free-particle mass oscillation, i.e., $2H(a_{m0})=ma_{m0}$, and the critical wave number for the free particle model, $k_{c}=(2maH)^{1/2}=ma_{m0}$; by the same token, we have $a_{k}/a_{m0}=k_{c}/k$ using the definition of horizon entry that $a_{k}k=2aH$. (The quantity $aH=a_{m0}H(a_{m0})$ since it is redshift-independent in the radiation era, and we thus have $k_{c}\propto m^{1/2}\propto a_{m0}^{-1}$.)

As a supplementary remark, the above estimate for the duration of unstable phase has taken into account that prior to this parametric growth, low-$k$ super-horizon mode must go through the $a^{2}$ growth of phase (c) even after the onset of nonlinear mass oscillation, where the driving source $Sr$ is still strong and Eq. (8) is not valid; for such modes, only after horizon entry, $a=a_{k}$, does Eq. (8) become valid and hence the solution of this equation starts at $a_{k}$.

For these sub-horizon $k$ modes subject to parametric instabilities, the amplitudes increase by a growth factor proportional to $\exp[A(a_{\kappa})-A(a)]$, and the exponent of the growth factor at the end of parametric growth can be shown to be

[TABLE]

using WKB approximation, where $a_{\kappa}=a_{k}$ and $a_{\kappa}=a_{m}$ correspond to $a_{k}>a_{m}$ and $a_{k}<a_{m}$, respectively, $A(a_{0})$ is defined to be [math], and $\theta_{0}$ takes the value $\pi$. (See Appendix for derivation.) This growth factor is responsible for the power excess. When $(a_{m}/a_{\kappa})[(a_{m}/a_{m0})(k_{c}/k)(\pi/2)]^{2}\gg 1$, the growth factor $\exp[A(a_{\kappa})]$ becomes $\exp[\pi(k/k_{c})(a_{m}/a_{m0})^{3/2}(a_{m0}/a_{\kappa})^{1/2}]$. Therefore for long waves we have small amplification, $A(a_{k})\propto(k/k_{c})^{3/2}$ as $a_{m0}/a_{k}=k/k_{c}$. On the opposite limit for short waves, $k>k_{c}(a_{m}/a_{m0})(\pi/2)$, the parametric growth would never occur. This explains why we do not see the power excess for long waves and short waves in Fig. (1).

This unstable phase is followed by a matter-wave oscillation phase but with a lower frequency than normal. The solution in this phase has a form $\sin(B(\eta)+\pi/4)$, where the detail is also given in Appendix using WKB approximation, and we have

[TABLE]

Here $\eta\equiv\ln(a/a_{m})$ and $r=(k_{c}/k)^{2}(\pi^{2}/4)(a_{m}/a_{m0})^{2}$. This oscillation has an initial $B(\eta)=0$ at $\eta=\eta_{0}$ where $\omega=0$, but otherwise $B(\eta)\geq 0$ for $\eta\geq\eta_{0}$.

The peak of the power excess (feature 3) should be located in this oscillation phase since the solution $\sin(B(\eta)+\pi/4)$ is still on the rise at $\eta_{0}$. If one is to assume that solutions have resumed free-particle matter-wave oscillations when they reach the peaks, i.e., $r\exp[-\eta]\ll 1$, then one may find the timing $\eta_{peak}(\equiv\ln(a_{peak}/a_{m0}))$ of the solution peaks as a function of $k$, $m$ and nonlinearity given by $\sin[B(\eta_{peak})+\pi/4]=1$, or $B(\eta_{peak})+\pi/4=\pi/2$. We thus have

[TABLE]

where terms in the squared bracket are contributed from the growing phase.

Subtracting $\eta_{peak0}$ of the free-particle model from this $\eta_{peak}$ , we can determine the total delay in the first quarter cycle of nonlinear mass oscllation. As shown in Appendix (A) of Paper (@slowromancapi@), the free-particle model has a oscillating solution $\propto\sin[(k/k_{c})^{2}[\eta-\ln((\sqrt{3}/2)(a_{\kappa}/a_{m0})+c_{0}(k)]]$ in phase (d). Here, the phase $c_{0}(k)=\cos^{-1}[1/((1+(k/k_{c})(\gamma-0.5))^{2})^{1/2}]$ with $\gamma$ being the Euler number $\sim 0.577$, and thus $c_{0}(k)$ is nearly [math] for a wide range of $k/k_{c}$. Now, using $(k/k_{c})^{2}(\eta_{peak0}-\ln((\sqrt{3}/2)(a_{\kappa}/a_{m0}))=\pi/2$ to fix $\eta_{peak0}$, we obtain the total delay $\Delta\eta(=\eta_{peak}-\eta_{peak0})$ as

[TABLE]

where $a_{k}/a_{m0}=k_{c}/k$ has been used, and $\kappa=0$ for long waves where $a_{\kappa}=a_{k}$ and $\kappa=1$ for short waves where $a_{\kappa}=a_{m0}$.

This is an interesting prediction, in that the strong negative $(k_{c}/k)^{2}$ dependence of $\Delta\eta_{peak}$ can make the delay be negative. The cause of the reverse effect is that the growing phase of parametric instability brings the amplitude to $1/\sqrt{2}$ of the peak in a time weakly dependent on $k_{c}/k$. This period can be short compared to the free matter-wave oscillation to bring the amplitude to a similar level for long waves, which take a time $\delta\eta\sim(4/\pi)(k_{c}/k)^{2}$. The maximum $\Delta\eta_{peak}$ can be found by taking a derivative of it with respect to $k_{c}/k$ and the maximum delay is found to be near $k/k_{c}\sim 1$. This explains why the delay in the first half cycle of nonlinear mass oscillation is prominent only around $k\sim k_{c}$ in Fig. (1).

Finally, since our results above depend on $a_{m}$, it is useful to pin down the relation between $\delta\theta_{0}(\equiv\pi-\theta_{0})$ and $a_{m}$. One can Taylor expand the field potential gradient near $\theta_{0}=\pi$ in Eq. (1) where $\sin(\theta)\approx\theta-\pi=\delta\theta$. Since the results, Eqs. (12), (14) and (15), are derived using Mathieu’s equation, Eq. (9), to be consistent with these results, one should define $a_{m}$ in accordance with the power-law-amplitude assumption. We extrapolate the asymptotic power-law solution, $\langle\theta^{2}\rangle\propto t^{-3/2}$, backward in time till it intercepts the actual background solution at $\pi^{2}$ to define the onset time $t_{m}$. In so doing, we find the following analytical formula provides the best fit:

[TABLE]

where $\delta\theta_{01}=1^{o}$. This expression works fairly fine; using it to compute $\eta_{peak}$ of Eq. (14) gives $<10\%$ errors against the measured $\eta_{peak}$.

V Numerical Solution and General Nonlinear Model

In Fig. (2), we plot $\Delta_{\theta}$’s constructed from the numerical solutions of Eq. (2) for passive evolution. The fiducial particle mass $m=10^{-22}$ eV is chosen. A comparison to $\Delta_{\theta}$’s constructed from Eq. (8) is also shown here, where we take $\beta=1$ and the metric fluctuation $\phi=0$ (c.f., Eq. (3)); the initial solution slope is set to $\dot{q}/q=[\omega_{\theta}^{2}(3\pi^{2}/8-1)-k^{2}/a_{m}^{2}]^{1/2}$ (c.f., Eq. (8)). Clearly seen in Fig. (2) is good agreement between the two solutions, except in the early time where our leading-order Taylor expansion of the nonlinearity fails. This plot demonstrates that peculiar features (2) and (3) appearing in the solution of Eq. (2) indeed arise from the parametric drive. These two features show strongly for $k\sim k_{c}$ modes than for $k\ll k_{c}$ modes as explained in the last section. In Fig. (2), we also plot a third solution of a fluid equation derived from the Gross-Pitaevskii equation, which is the subject of the next section, and in Appendix, we present the fluid equation, Eq. (A1). This fluid equation filters out the high-freqeuncy mass oscillation and is therefore an equation for the slowly varying amplitude. One can see that this third solution agrees with the solution of Mathieu’s equation, Eq. (8), extremely well.

The matter spectra not long after the radiation-matter equality is particularly interesting since by then the primary spectral feature can hardly evolve and is frozen throughout the later evolution. We plot in Fig. (3) the transfer function, $|\Delta_{\theta}/\Delta_{cdm}|^{2}$, of several initial field angles $\Delta\theta_{0}$’s at $a=a_{eq},2.5a_{eq}$ and $5a_{eq}$ using the full treatment, where $a_{eq}$ is the scaling factor at the radiation-matter equality. One can clearly see the broad spectral bumps in all extreme initial angles. For $k$ at and smaller than the spectral peak the transfer function barely evolves after $a_{eq}$, but are opposite for larger $k$. This is due to the Jeans wave number around $a_{eq}$ is close to $k_{c}$ [Paper (@slowromancapi@)]; above the Jeans wavelength, perturbations grow self-similarly as those of CDM but below, perturbations become neutrally stable oscillating matter waves Hu et al. (2000); Woo, and Chiueh (2009).

In evolving toward $a=5a_{eq}$, the photon-electron decoupling occurs around $a=3a_{eq}$ and the photon perturbation contributes very little to the metric perturbation since then. In Appendix B of Paper (@slowromancapi@), we showed that the drag between baryons and photons already damps out the photon perturbation prior to photon-electron decoupling for $k>k_{c}$, and hence the metric perturbation indeed has no contribution from photons regardless of whether or not we have considered the electron-baryon recombination physics. But for $k\ll k_{c}$, the drag is ineffective and hence metric perturbations are affected by photons at $a>3a_{eq}$ if the silk damping is not properly accounted for, and can produce some errors in the matter power spectrum. Our full treatment does not take in to account the silk damping. However, this error for matter perturbations is practically small since matters are cold and photons are hot, and gravity responds to cold matter. Moreover, these small errors in both CDM and axion perturbations are the same in our full treatment since at long waves the two perturbations are almost identical. Hence the transfer function is insensitive to such errors present in their respective spectra. It is based on this rationale the transfer functions at $a=5a_{eq}$ are presented in Fig. (3).

In Fig. (3), we have made sure that the dark matter energy density $m^{2}f^{2}\theta^{2}$, together with baryon energy density, equals the radiation energy density at $a_{eq}$. As $\theta_{0}$ approaches $\pi$ that delays the onset of nonlinear mass oscillation, the value of $f$ must decrease to satisfy the above condition. Normally the field strength is characterized by a dimensionless parameter $f^{2}/m_{p}^{2}(\equiv 4\pi Gf^{2})$, where $m_{p}$ is the Planck mass. The free-particle case corresponds to $f^{2}/m_{p}^{2}\to\infty$. The appropriate parameter regime for the extreme axion model has values $f^{2}/m^{2}\sim O(10^{-5})$, and so $f$ is on the order of GUT scale.

To further demonstrate the general validity of Eq. (8) approximating the original perturbed field equation, we consider the potential $(m^{2}a^{2})[1-\operatorname{sech}(x)]$, where $x$ is the field. We replace mass terms in the field equations, Eqs. (1) and (2), by $(m^{2}a^{2})\tanh(x)\operatorname{sech}(x)$ and $(m^{2}a^{2})\operatorname{sech}(x)[-1+2\tanh^{2}(x)]\delta x$, respectively, which yield $\beta=1$, $\alpha=(5x_{0}^{2}/8)(t/t_{m})^{-3/2}$ and $\omega_{x}^{2}=m^{2}-\alpha$. Here we also choose the initial field value $x_{0}=\pi$. Plotted also in Fig. (2) is the comparison of $\Delta_{x}$’s constructed from the passive evolution and Eq. (8). Again, excellent agreement is found when $t\gg t_{m}$, reinforcing our claim for the parametric drive of the original perturbed field equation.

To end this section, we notice that $\beta=1$ is quite generic to all symmetric field potentials, and this can be shown as follows. Let the Taylor expansion of the field potential be $V(x)=(m^{2}a^{2}/2)(x^{2}-(b/2)x^{4}+O(x^{6}))$, and the potential gradient $V^{{}^{\prime}}=(ma)^{2}(1-bx^{2})x\sim m^{2}(1-(3b/4)x_{0}^{2})x\equiv\omega_{x}^{2}x$, where $x_{0}$ is the oscillation amplitude, and $\omega_{x}$ is the nonlinear driving frequency adopting the technique used for the axion case. The coefficient of the restoring force in the perturbation equation is $V^{{}^{\prime\prime}}=(ma)^{2}(1-3bx^{2})$, which can be reduced to $\omega_{x}^{2}[1-(3b/4)x_{0}^{2}(1+2\cos(2\omega_{x}t))]$. In Sec. (IV), we have parametrized the last factor as $(1+2\beta\cos(2\omega_{x}t))$, and hence $\beta=1$ for all nonlinear models with symmetric potentials with a finite mass.

VI Axion connection to Gross-Pitaevskii Equation

The Gross-Pitaevskii equation reads

[TABLE]

This non-relativistic equation can be deduced from the Bogoliubov’s formalism in the zero-temperature limit of a dilute interacting Bose-Einstein condensate Bogoliubov (1947); Rogel-Salazar (2013). The interaction energy density is proportional to $g|\Phi|^{4}$, where $g$ is a coupling constant related to the microscopic scattering length $a_{0}$ as $g\equiv(4\pi/m)a_{0}$. The dilute gas approximation demands that $|a_{0}|^{3}|\Phi|^{2}\ll 1$, and the interactions are repulsive (attractive) when the coupling constant $g$, or the scattering length $a_{0}$, is positive (negative).

The dynamics of the linear perturbation $\delta\Phi=\Phi-\Phi_{0}$, with the background field $\Phi_{0}$ chosen to be real, can be best seen using a fluid approach. The field $\Phi$ is expressed in the polar coordinate as $\Phi=\xi\exp[i\chi]$, the fluid density $n\equiv\xi^{2}$ and the fluid velocity ${\bf v}\equiv\nabla\chi/m$, the quantum specific pressure $-(\nabla^{2}\xi)/(2m\xi)$ and the fluid specific pressure $gn$. The total specific pressure is $-\nabla^{2}\xi/(2m\xi)+gn$, and when $g$ is negative the total specific pressure can become negative in the long-wave limit and it is straightforward to show that the dispersion relation is $\omega^{2}=(k^{2}/2m)[k^{2}/(2m)+2mgn]$. (See Appendix.) If we identify $\alpha=-2mgn$ with a negative $g$ and $k/a$ replaces $k$, this dispersion relation is identical to Eq. (11) for the axion model.

To compare with Mathieu’s equation Eq. (4) in detail, we examine the excitations of the Gross-Pitaevskii equation, Eq. (17). Note that the background filed ($k=0$) has a frequency $\omega_{0}=g|\Phi_{0}|^{2}$, and $\delta(|\Phi|^{2}\Phi)=2|\Phi_{0}|(\Re[\Phi_{0}]\Re[\delta\Phi]+\Im[\Phi_{0}]\Im[\delta\Phi])+|\Phi_{0}|^{2}\delta\Phi$. We can remove the above complication by defining a frequency-shift wave function $p=\Phi\exp[ig|\Phi_{0}|^{2}t]$ for all modes including the background field, then $\delta(|\Phi|^{2}\Phi)=\exp[-ig|\Phi_{0}|^{2}t]|\Phi_{0}|^{2}(2\delta p+\delta p^{*})$. Therefore the linearized equation becomes

[TABLE]

This equation is identical to Eq. (10) for the axion model ($\beta=1$) after replacement of coefficients: $\alpha=-2mg|\Phi_{0}|^{2}$ and $a=1$. Thus, aside from the source due to metric perturbations, the linear perturbations of Eq. (5) and the excitations of Gross-Pitaevskii’s equation are the same. Not only that, the background field of Eq. (5) and the ground state of the Gross-Pitaevskii equation are also the same. The ground state is uniform with a frequency $g|\Phi_{0}|^{2}$, and the background field is also uniform having a frequency $m-(\alpha/2m)$. Subtracting off the leading order mass oscillation frequency to get to the non-relativistic regime, we find the ground state and the background field oscillating at the same frequency. Therefore the axion model can therefore be fit into Bogoliubov’s framework of dilute interacting boson gas. (We have noted that there have been previous works on the connection between the perturbed scalar-field equation, Eq. (2), and the Gross-Pitaevskii equation, but unfortunately in a rather ad hoc manner where the nonlinear shift of the driving frequency was not taken into account thus yielding a higher nonlinear coupling strength $g$ by a factor $12$ than it actually is Magaña and Matos (2012).)

The interaction potential energy for axion, to the leading order, is about $(mf)^{2}\theta^{4}=(m/f)^{2}\Psi^{4}$, and so $|g|=f^{-2}$, since the mass (energy) density $2m^{2}\langle\Psi^{2}\rangle=m|\Phi|^{2}$. It implies that the scattering length $|a_{0}|=m/(4\pi f^{2})$, and the dilute gas condition, $|a_{0}|^{3}|\Phi|^{2}\sim(4\pi)^{-3}(m/f)^{4}\ll 1$, valid to an excellent degree. However, for all practical purposes, this naive estimation gives too small a $|a_{0}|^{3}|\Phi|^{2}$ by many orders of magnitude for $m$ as small as $10^{-22}$eV, due to the fact that the scattering length $|a_{0}|$ is many orders of magnitude smaller than the Planck length. This problem arises from the fact that the gravity has been ignored in Bogoliubov’s formulation. If the scattering length $|a|$ is limited to the smallest possible length, the Planck length $l_{p}$, the dilute gas condition becomes $l_{p}^{3}f^{2}m$ and is still much less than unity even when $f$ takes the largest possible value, the Planck scale $l_{p}^{-1}$, therefore justifying the dilute gas approximation.

Having the microscopic physics in place, we can now ask how reasonable it is for the initial angle so close to the top of the field potential, as quantum tunneling may have made its way to render the system unstable. Note that the macroscopic axion field $\Psi=\sqrt{N}\psi_{q}$, where $\psi_{q}$ is the quantum field of individual particle and $N$ is the number of particles overlapped in some macroscopic coherent length. The background field has an infinite coherent length in the zero-temperature limit. Let us take a conserved position where the perturbed field has a coherent length comparable to one wavelength $2\pi/k$; the number $N$ is still a huge number due to the extremely small particle mass. As an example, the fiducial mode $k=k_{c}\sim 10$ Mpc*-1* encloses a particle number $N=n_{m}\lambda_{c}^{3}\sim 10^{98}$, where $n_{m}$ is the number density at the onset of free-particle mass oscillation and $\lambda_{c}$ is the Compton wave length of the fiducial particle mass $m_{22}=10^{-22}$ eV. Therefore the particle number amounts to $N\sim 10^{98}(m/m_{22})^{-5/2}(k/k_{c})^{-3}$. The quantum tunneling is suppressed by the factor $\prod_{i}\exp[-S_{i}/\hbar]$, where $S_{i}$ is the Eulerian action of an individual particle $i$. The suppression factor $\exp[-S_{i}/\hbar]$ is usually taken to be a Gaussian around the classical field $\psi_{c}$. Though individual $S_{i}/\hbar$ can be small enough to permit quantum tunneling, especially near a classical bifurcation point, the coherent N particles share the same phase space coordinates, thus $S_{i}=S$, and the suppression factor becomes $\exp[-NS/\hbar]$. It greatly narrows the variance of the Gaussian around the classical field. Hence, tunneling through quantum fluctuations is impossible even when the initial field angle is very close to the potential top.

In Bogoliubov theory, there is another quantity, the healing length $l_{h}$ defined to be $(2m|g||\Phi|^{2})^{1/2}l_{h}\equiv 1$, which characterizes the relative strength of destabilizing nonlinear to stabilizing linear terms of Eq. (8). In the absence of self-gravity, structures beyond $l_{h}$ grow due to the weak instability of negative $g$ until nonlinear structures, such as vortex filaments Carr and Clark (2014), form on the scale of $l_{h}$. But with self-gravity, smaller nonlinear structures can form through much stronger gravitational instabilities. The healing length $l_{h}$ is about the Compton wavelength at the onset of mass oscillation, and since $l_{h}\propto a^{3/2}$ in physical coordinate, it becomes about $150$ Mpc in the present universe for particle mass $m=10^{-22}$ eV. It is interesting to note that this healing length coincides with the observed correlation length of the baryon acoustic oscillation at present Eisenstein et al. (2005). Will the healing length have cosmological footprints on very large-scale structures? As the baryon acoustic oscillation scale depends linearly on $a$ and the healing length on $a^{3/2}$, future observations at high redshift will tell.

VII Conclusion

In this work, we analyze the three unexpected features of the extreme axion and explain their underlying mechanisms. Among them, the parametric drive and amplification mechanism accounts for two non-trivial features. To illustrate of the mechanisms, we show the original perturbed field equation can be made equivalent to Mathieu’s equation, which is able to faithfully recover the two features. We also disclose that in the non-relativistic limit and to the leading-order nonlinearity, the equation of motion for the axion model is identical to the Gross-Pitaevskii equation, a macroscopic manifestation of a zero-temperature interacting Bose-Einstein condensate. Based on this connection, we explain why the quantum tunneling of this system is impossible.

The two nontrivial features, i.e., extension of spectral cutoff to higher wave number and spectral excess of the extreme axion model can have important impacts in structure formation of the high-$z$, matter-dominant universe, due to the fact that most first-generation galaxies formed out of perturbations near the spectral cutoff. The spectral cutoff is determined solely by particle mass in the free-particle model. Its extension to higher $k$ for the extreme axion model mimics the effect of higher particles mass for free particle. Therefore, the high-$k$ power spectrum may not be a good indicator accurately reflecting the true particle mass in the extreme axion model. Recent simulations addressing the high-redshift Lyman-$\alpha$ absorption features indicate that substantially higher particle mass than $10^{-22}$eV is required or implied Armengaud et al. (2017); Iršič et al. (2017). On the other hand, approximately $10^{-22}$eV particle mass is needed to account for the flat cores of dwarf spheroidal galaxies Schive et al. (2014); Chen et al. (2017). The tension in particle mass may be lessened with the extreme axion model.

However we have found a limit to the high-$k$ spectral extensions, no matter how extreme a condition the initial angle assumes. The spectral extensions are all confined to wave numbers less than a factor $2$ higher from that of the free-particle model, i.e., approximately $<k_{c}$. That is, the spectral excess peaks around $0.6-0.8k_{c}$ and immediately following the spectral peak is a sharp cutoff. This spectral shape renders the first collapsed halo of mass $[(4\pi/3)(\pi/k)^{3}]\rho_{0}\sim 10^{10}M_{\odot}$, where $k$ is near the peak of the spectral excess and $\rho_{0}$ is the background mass density. As a reference, the first galaxies in the free-particle axion model of $m=10^{-22}$ eV have masses several $\times 10^{10}M_{\odot}$ Schive et al. (2016).

The spectral excess is perhaps our most surprising finding, since conventional dark matter candidates proposed so far are unable to produce power excess over the CDM model across the perturbation spectrum. When $\delta\theta_{0}<1$ degree, the spectral excess can be so distinct that may completely revise the standard scenario of first galaxy formation. First of all, the spectral excess leads to earlier formation of first generation galaxies and push the reionization epoch Planck Collaboration et al. (2016) earlier than the free-particle model Schive et al. (2016). Second, taking the more extreme case $\theta_{0}=0.2$ degree as an example, c.f., Fig. (3), the broad spectral peak yields first collapsed halo of mass $10^{9}-10^{10}M_{\odot}$, and frequent mergers of these over-abundant first halos than the conventional are to quickly build up more massive halos. Furthermore, busy mergers are prone to sustain intense star bursts and even rapid super-massive black hole growths. Finding quasar at $z=7$ Mortlock et al. (2011) and recent discovery of galaxies more massive than the Milky Way inferred to already form at $z=7$ Glazebrook et al. (2017) have posed challenges for the CDM model. Given the aforementioned possible outcomes, the extreme axion model may stand a better chance to meet such a challenges. Whether successful or not, only future simulations can tell.

To place our results in a concrete ground, we provide a formula for the wave number of spectral peak $k_{peak}$ and a procedure for the peak height to be calculated, as functions of the initial angle and the particle mass in Appendix. Aside from that, our identification of the perturbed field equation to the Gross-Pitaevskii equation in the non-relativistic limit is of practical relevance. It permits calculations of pertubation dynamics using a fluid approach outlined in the beginning of Sec. (VI) and carried out in Appendix.

Finally, we must stress that the extreme axion model is not a special model capable of producing the three peculiar features studied in this work. A wide range of scalar field models have the same characteristics, as demonstrated in Sec. (V) by an example. In any of these scalar field models, the particle mass $m$ has to be extremely light, not far from our fiducial mass $10^{-22}$ eV, to produce astronomical observable effects. As to the nonlinear effect of general scalar field models, if the coupling constant $g$ in the Gross-Pitaevskii equation, proportional to the microscopic scattering length, increases with time as the scaling factor $a$, then the parametric instability explored in this work will have a long-lasting, but weak, effect for perturbations beyond a particular length scale.

Appendix A Particle Mass Dependence

In Sec. (IV) we consider how solutions change with a changing $a_{m}$ for a fixed particle mass $m$, and concluded that when $a_{m}$ becomes large, the effect is equivalent to make $k$ appear smaller. When the particle mass $m$ is allowed to change, Eq. (10) is invariant to a transformation, which we call the mass-tempo-size trasformation. We normalize time $t$ to $t_{m}$, and note that $\alpha\propto m^{2}(a_{m}/a)^{3}$. It is straightforward to show that the solution is invariant, up to a shift in $\ln(a)$ space, to the changing $a_{m}$, $k$ and $m$ so long as $ka_{m}$ and $ma_{m}^{2}$ are kept fixed. This transformation for a changing $a_{m}$ is in fact a transformation only of the particle mass $m$. We note that $\ln(a_{m})=\ln(a_{m0})+\ln K(\delta\theta_{0})$ for a $K(\delta\theta_{0})$ given in Eq. (16); the quantity $\ln(a_{m0})$ depends only on the particle mass $m$ and another quantity $K(\delta\theta_{0})$ on how nonlinear the background field is, and hence the particle mass and the nonlinearity appear to be able to vary independently. However, the transformation requires a fixed $ma_{m}^{2}$, which becomes $ma_{m0}^{2}K^{2}(\delta\theta_{0})=(aH)K^{2}(\delta\theta_{0})$ where $aH$ is independent of the particle mass $m$, the scaling factor $a$ and the nonlinearity. As a result, a fixed $ma_{m}^{2}$ implies fixed nonlinearity $K(\delta\theta_{0})$ in the tranformation, and therefore this transformation involves only a varying particle mass $m$.

Realizing it, the $m$ dependence of the spectral peak $k_{peak}$ can be straightforwardly obtained. This can be carried out since an analytical solution of the perturbed field can be found to a good approximation. Adopting the fluid approach to the Gross-Pitaevskii equation, Eq. (17), in the comoving frame with proper replacement of $\nabla\to ik/a$ and others explained in last section, Sec. (VI), we find that the sub-horzon normalized fluid density $\delta n/n$ satisfies

[TABLE]

where $\eta\equiv\ln(a/a_{m})$ and $\delta n/n=2\Re[\hat{q}]\approx\Delta_{\theta}$.

As a check of the accuracy of this fluid equation to Eq.(8), we solve Eq.(A1) numerically. The initial $\eta_{i}$ starts from $a_{\kappa}/a_{m0}$ with the initial slope $S(\eta_{i})\equiv(k/k_{c})^{2}[1-(k_{c}/k)^{2}(\theta_{0}^{2}/4)(a_{m}/a_{m0})^{2}\exp[-\eta_{i}]]^{1/2}$. Solutions are plotted in Fig. (2), and one can see that the solutions excellently agree with those of Mathieu’s equation.

One may adopt the WKB approximation to analyze the solution, for which the phase, $\int d\eta S(\eta)$, has an analytical expression. When the integrand is imaginary, it represents a growing solution $\exp(A(\eta_{i})-A(\eta))$, as $A(\eta)$ is a decreasing function of $\eta$. When the integrand is real, it has an oscillating solution, $\sin(B(\eta)+d_{0})$, where $B(\eta)\geq 0$ and is an increasing function of $eta$, and $d_{0}$ is the phase. In between, the integrand crosses zero, the WKB approximation fails and we have an Airy function that connects solutions on the two sides, i.e., $\exp[A_{(}\eta_{i})-A(\eta)]\to 2\sin(B(\eta)+\pi/4)$. Here, the analytical expressions of $A$ and $B$ are

[TABLE]

where $r=(k_{c}/k)^{2}(\pi^{2}/4)(a_{m}/a_{m0})^{2}$, and both $A(\eta)$ and $B(\eta)$ equal zero when the integrand crosses zero. The expression of $A(\eta)$ has been used in Sec.(IV) to find the growth factor.

When $\eta\gg 1$, we recover that free-particle oscillation, i.e., $B\to(k^{2}/k_{c}^{2})(\ln(a/a_{m})+d)$ with a phase $d$. Note that the peak of the solution is located at the oscillating side of the solution. The relation between $k_{peak}$ and $m$ is just the solution of a transcendental equation $B(\eta_{eq})+\pi/4=\pi/2$, where the sine oscillation phase is equal to $\pi/2$, where $\eta_{eq}$ is the value of $\eta$ at radiation-matter equality. We make a further approximation to simplfy the matter. The nonlinear contribution to $B$ is negligible at the peak of the solution, i.e., $r\exp^{(}-\eta)\ll 1$, so that $B(\eta)\approx(k/k_{c})^{2}(\ln(a/a_{m0})+[\ln(a_{m0}/a_{m})+2(\ln(2)-1)-\ln(r)])$. This approximation is better for a large $m$ than for a small $m$. We thus obtain a simpler transcentdental equation for $k_{peak}$:

[TABLE]

which can be solved numerically rather easily. It can be easily seen that $k_{peak}/k_{c}<1$ for $m=10^{-22}$ eV, and more so for a larger particle mass because $\ln(a_{eq}/a_{m0})$ gets larger. Aside from the mass dependence of $k_{c}\propto m^{1/2}$, $k_{peak}$ has another weak mass dependence on $\ln(a_{m0}/a_{eq})$, with $a_{m0}\propto m^{-1/2}$. Another term, $\ln(a_{m}/a_{m0})$, is a measure of $\delta\theta_{0}$, which is given in Eq. (16) and has no mass dependence.

Though $k_{peak}$ is derived here from the passive evolution, $k_{peak}$ of the full treatment deviates only slightly from this formula; thus, to a good approximation Eq. (A3) provides an analytical expression for $k_{peak}$ of the full treatment, and we find this expression is accurate within $10\%$ of the peak of $k^{3}|\Delta_{\theta}(k)|^{2}$ of the full treatment. Moreover, $k_{peak}$ is largely frozen after $a=a_{eq}$ shown in Fig. (3), as $k_{peak}$ is smaller than the Jeans wave number $k_{J}$ in the matter-dominated regime, and therefore this spectral peak persists in the linear matter power spectrum throughout the later epoch.

The particle mass dependence of the quantity $k_{peak}/k_{c}$ is mild, shown in the right-hand side of Eq. (A3). When the particle mass $m\to\infty$, we have $\ln(a_{eq}/a_{m0})\to 0$, the peak $k_{peak}/k_{c}\to 0$, and the growth factor of Eq. (12) approaches zero. The particle mass scaling of the growth factor can be calculated by determining $k_{peak}$ and substituting into the growth factor Eq. (12). Comparing the growth factor for $m=10^{-22}$ eV to find the ratio, one is then able to determine the spectral peak height for any $m$ by referring to the peak height of Fig. (3) for $m=10^{-22}$ eV, which can be approximated to be $\ln|\Delta_{\theta}/\Delta_{cdm}|^{2}|_{peak}=17(k_{peak}/k_{c})-10.2$ as a fit.