1 Introduction↩︎

The QCD axion remains one of the dark matter (DM) candidates with the most significant motivation, and the most motivated such candidate with a non-thermal primary production mechanism. Since the intuitive picture of DM as a thermal relic is quite prevalent in the literature and pedagogical sources, the precise behavior through cosmic time of a non-thermal DM candidate such as the axion can sometimes be elusive. In particular, the process by which the axion excites adiabatic field perturbations in the early universe, arriving at dynamics nearly identical to cold dark matter (CDM), is not as straightforward as the equivalent process for a thermal candidate.¹ Similarly, the long-term evolution of initial isocurvature fluctuations is not as straightforward as in the CDM case. For these reasons, this work expounds upon the generation mechanism and early superhorizon evolution of field fluctuations for the axion and the corresponding density perturbations. The formulae and procedures laid out in this work are generalizable to other scalar DM scenarios beyond the QCD axion, in particular scenarios where the DM does not thermalize with the standard model (SM) bath.

The QCD axion is a compact pseudo-scalar field \(\theta\), postulated to explain the observed non-CP violating nature of QCD. It can arise as a pseudo-Nambu-Goldstone boson for a spontaneously broken, approximate global \(U(1)\) Peccei-Quinn (PQ) symmetry [1]–[4], or as a mode of an extra-dimensional gauge field [5]. The axion DM population is primarily produced via an initial vacuum misalignment in which the scalar is frozen in its potential until the time when the Hubble rate of expansion drops below the effective mass of the field, at which point the field begins to undergo oscillations with an energy density that redshifts like matter [6]–[8]. This misalignment mechanism proceeds differently depending on whether the axion field is already present as a massless degree of freedom during inflation (pre-inflationary), or emerges from a PQ-breaking phase transition afterwards (post-inflationary). In the pre-inflationary case, the initial phase \(\theta_\text{ini}\) is approximately uniform in the observable universe, and thus the axion field is (initially) homogeneous. The post-inflationary case has a complicated cosmology involving production of topological defects, but because the PQ scalar is initially in thermal equilibrium with the SM, the origin of adiabatic fluctuations is less mysterious. For reviews of both cosmological scenarios, see, e.g., [9]–[11].

Let us focus on the pre-inflationary axion case. In simple inflationary scenarios, the inflaton exhibits only adiabatic fluctuations, which are then inherited by SM fields following reheating as the universe attains local thermal equilibrium [12]–[16]. Then, with all fields in equilibrium, the perturbations remain adiabatic and their evolution proceeds according to the usual cosmological perturbation theory and decoupling histories, setting the stage for the cosmic background radiation and evolution of structure that follows [17]–[20]. However, since the axion is present as a massless field during inflation, entropy-type (isocurvature) inhomogeneities in the axion field are induced, uncorrelated with the adiabatic (curvature) perturbations of the inflaton. For the axion to be a viable DM candidate, these isocurvature perturbations should be small so as not to exceed constraints on DM isocurvature in the cosmic microwave background (CMB) [21]–[23]. Assuming these inhomogeneities are suppressed, and in the absence of large couplings to the inflaton or the SM particles, the axion field would remain homogeneous in the early universe as it does not come into thermal equilibrium with the rest of the plasma. On the other hand, after the axion begins oscillations, it will behave like DM and thus should have the same inhomogeneities as the adiabatic DM fluctuations (which should be related to the adiabatic perturbations of all other fluids). In fact, as made clear by a theorem of Steven Weinberg regarding the existence of an adiabatic solution to the perturbed Einstein equations at superhorizon scales [24]–[26], the initial homogeneous state of the axion is already adiabatic, and it will proceed to remain adiabatic into the regime of its DM-like behavior. Several works in the literature have elaborated on the behavior of the axion and axion-like-particle (ALP) perturbations and the adiabatic nature of the perturbations, with various assumptions about the cosmic background evolution and axion potential [10], [27]–[31].

The QCD axion, in addition to being non-thermal DM, has other complications to its early-time dynamics. The axion acquires a potential from QCD dynamics, which gives it an effective mass and higher order couplings. This potential, however, depends on temperature until the completion of the QCD phase transition at a temperature \(T_\text{QCD} \sim \mathcal{O}(100 MeV)\). While the evolution of the homogeneous background field remains similar to DM in the presence of this temperature-dependent potential, the perturbations of the field deviate significantly from the DM-like behavior with new contributions arising from temperature fluctuations in the plasma (see, e.g., [32]–[35] for other sources pointing out this contribution to axion perturbations). This deviation only contributes to the axion field fluctuations in the era prior to the QCD phase transition, but it is not immediately obvious that its effects do not appear at late times by, for instance, spoiling the adiabatic nature of the axion perturbations.

Both the non-thermal nature of the production of axions and the temperature dependence of the QCD axion potential prior to the QCD phase transition seem at first troublesome for guaranteeing that the axion eventually exhibits an adiabatic fluctuation equivalent to DM. We will see that the adiabaticity of the inhomogeneities survives both of these hurdles, spectacularly guaranteed by Weinberg’s theorem. Beyond the adiabaticity of the scalar perturbations, the isocurvature perturbations of the axion also relate to the background evolution. In the context of the QCD axion, the temperature dependent axion potential will give rise to deviations from an otherwise constant entropy perturbation. We will also underscore the influence of the background evolution on the isocurvature spectrum and non-Gaussianity (via the bispectrum). Though various elements of this discussion are present throughout the literature, a comprehensive explanation of the excitation and evolution of the axion’s DM-like perturbations in a general cosmological background is lacking. In the rest of this work, we will show explicitly that the axion indeed exhibits perturbations which are automatically adiabatic, regardless of its non-thermal production and complicated potential, by first reviewing the theorem by Weinberg followed by an explicit example calculation. We will also comment on the behavior of isocurvature perturbations, which may be excited during inflation, throughout the complex early evolution of the axion field.

2 Adiabatic and Entropy Perturbations of Axion↩︎

Fields in thermal equilibrium with each other acquire correlated perturbations, so thermally produced dark matter automatically has adiabatic perturbations. The axion is instead produced via the misalignment mechanism and is therefore not immediately required to be adiabatic. Indeed, the standard pre-inflationary axion naturally exhibits fluctuations during inflation which are uncorrelated with the inflaton, resulting in late-time DM isocurvature. In this context, we aim to demonstrate clearly how the axion DM fluid naturally excites an adiabatic mode, despite the complications of the post-inflationary universe such as the QCD phase transition, and that the late-time isocurvature could differ from what is generated by the end of inflation.

The existence and conservation of the adiabatic solution regardless of background cosmological dynamics was shown by Steven Weinberg in a theorem first published in a series of foundational works [24], [25] and later expounded upon in his textbook Cosmology [26]. Weinberg’s theorem guarantees that there always exists a solution to the perturbed scalar field equations in Newtonian gauge (generalizable to other gauges) containing constant scalar curvature \(\zeta\) on superhorizon scales with a wavenumber \(k\) below the conformal Hubble scale, \(k\ll R H\). The perturbed metric in Newtonian gauge is given by \[{\rm d}s^2 = -(1+2 \Psi){\rm d}t^2 + R^2(t)(1-2 \Phi) {\rm d}\boldsymbol{x}^2 \, ,\] where \(R(t)\) is the scale factor, \(H(t)\equiv\dot{R}/R\) is the Hubble rate, and the metric perturbations are observed to be small \(\Psi(t, \boldsymbol{x}), \Phi(t, \boldsymbol{x}) \ll 1\). The curvature \(\zeta\) is \[\zeta \equiv -\Psi- H \frac{\delta\rho}{\dot{\overline{\rho}}} \, ,\] with \(\delta\rho\) and \(\overline{\rho}\) the total perturbed and background energy densities, and a dot denoting a derivative with respect to \(t\).² Weinberg’s argument constructs new inhomogeneous solutions to the equations of motion (up to small-\(k\) corrections) by first finding exact new homogeneous solutions and then promoting them to nonzero \(k\). Given a background solution to the Einstein field equations, the new homogeneous solution is constructed by performing a gauge transformation \(x_\mu \mapsto x_\mu + \epsilon_\mu\) which preserves the Newtonian gauge, \(\epsilon_\mu=(\epsilon(t), R^2(t)\omega_{ij}x^j)\), on the background solution to yield a scalar fluctuation \(\delta s = -\epsilon(t) \dot{\bar{s}}(t)\), where \(\delta s\) can be any four-scalar (\(\delta\rho\), \(\delta p\), \(\delta\theta\), etc.). This can always be done by virtue of diffeomorphism invariance. ³ Demanding that this solution is a physical solution to the field equations fixes \(\epsilon_\mu\), and therefore all of the metric and scalar fluctuations.⁴ In fact, there are two such solutions, one of which decays in an expanding universe and thus is ignored, while the other is often referred to as the (growing) adiabatic mode (this solution is also referred to as a curvature perturbation). Upon solving for \(\epsilon(t)\) and \(\omega_{ij}\), one finds that \[\begin{align} \label{eq:epssolution} \epsilon(t) &= - \zeta \mu(t), \nonumber \\ \text{where}\quad\mu(t) &\equiv \frac{1}{R(t)} \int^t R(t')\,\mathrm{d}t'. \end{align}\tag{1}\] Here \(\zeta\) is constant, and different choices of the lower limit of integration give rise to solutions differing by a multiple of the decaying mode (and hence are negligibly different at sufficiently late times). This mode can be promoted to an inhomogeneous mode of long wavelength by extending \(\epsilon(t)\) into a function \(\epsilon(t,\boldsymbol{x})\). Its Fourier modes, \(\epsilon_{\boldsymbol{k}}(t)\), are then related to those of the metric perturbation, \(\zeta_{\boldsymbol{k}}\), by the same proportionality factor \(\mu(t)\) as in 1 . They solve the equations of motion up to corrections of order \(k^2/(R^2 H^2)\).

The (non-decaying) adiabatic solution amounts to the following Fourier modes of scalar metric perturbations \[\label{eq:weinbergAdiabMetric} \Phi_{\boldsymbol{k}}(t)=\Psi_{\boldsymbol{k}}(t)=-\zeta_{\boldsymbol{k}} {\dot{\mu}}(t) = \zeta_{\boldsymbol{k}}\left(-1+\frac{H(t)}{R(t)}\int^t R(t')dt'\right).\tag{2}\] The characteristic distinguishing the adiabaticity of this solution is that each species has the same ratio of its energy perturbation to the rate of change of its background energy density \(\delta\rho_i/\dot{{\overline{\rho}}}_i\), regardless of whether each species’ energy is independently conserved. In fact, the proof of Weinberg’s theorem guarantees this relationship for any scalar quantity \(s\) [24], such that \[\label{eq:weinbergAdiabScalar} \delta s_{\boldsymbol{k}}(t)= -\zeta_{\boldsymbol{k}} \mu(t)\dot{\overline{s}}(t) = -\frac{\zeta_{\boldsymbol{k}}\dot{\overline{s}}(t)}{R(t)}\int^t R(t') dt' \, .\tag{3}\] Weinberg’s proof guarantees the existence of this adiabatic mode. Data informs us that this mode is a good description of our universe [21]. This is an important clue about initial conditions. Single-field inflation always excites the adiabatic mode, but more general physics in the early universe could have led to other solutions. Any deviation from this behavior, then, is non-adiabatic and is labeled an entropy perturbation, or equivalently an isocurvature mode. Let us define a generalization of scalar curvature \(\zeta_i\) from the perturbed energy of the species \(i\) as \[\zeta_i \equiv -\Psi- H \frac{\delta\rho_i}{\dot{\overline{\rho}}_i} \, .\] Then, the adiabatic mode exhibits equal scalar curvature for each species \(\zeta_i = \zeta_j = \ldots\), and any deviation from this behavior contributes to an isocurvature mode. We can define the relative entropy perturbation between two species as \[S_{ij} \equiv 3 \left[\zeta_i - \zeta_j\right] \, .\] In the absence of isocurvature, \(S_{ij}=0\), the arguments above guarantee that the curvature remains conserved, and thus the isocurvature perturbation remains absent. On the other hand, in the presence of nonzero isocurvature, Weinberg’s theorem does not guarantee that the isocurvature necessarily remain constant.

The power and elegance of this theorem is that it applies for arbitrary dynamics of the background cosmology. Examining our case of interest, the QCD axion, we will explicitly show that nontrivial adiabatic perturbations are generated when the axion becomes DM, despite the nontrivial dynamics of the axion. Weinberg’s theorem ensures this. At early times, the axion field is frozen, with \({\dot{\overline{\theta}}}(t) = 0\). The general result 3 , then, tells us that the adiabatic mode is the one with zero perturbations, \(\delta \theta_{\boldsymbol{k}} = 0\). (This point has been emphasized before, e.g., in section 4.4.1 of the review [10].) The theorem then guarantees that if the initial condition is in the adiabatic mode, at late times when the axion field oscillates and \({\dot{\overline{\theta}}}(t) \neq 0\), the solution will still track the adiabatic mode, now with \(\delta \theta_{\boldsymbol{k}}\) nonzero and proportional to the photon perturbations. Prior to the QCD phase transition, the instanton-induced potential of the axion depends on temperature, modifying the otherwise matter-like behavior of the scalar field and subsequently changing the behavior of the perturbed fields and energies. The remarkable power of Weinberg’s theorem is to guarantee that these complications lead to the same late-time solution (on superhorizon scales), provided \(\delta \theta_{\boldsymbol{k}} = 0\) at sufficiently early times.

In addition, for a pre-inflationary axion some small initial perturbation \(\delta\theta_{\rm ini}\) is expected to be excited during inflation, which will persist as an isocurvature (entropy) perturbation in the post-inflationary universe. In this case, Weinberg’s theorem cannot be straightforwardly applied. Especially upon accounting for the QCD phase transition, which generates a coupling between radiation and axions on large-scales, the nonzero entropy mode may evolve nontrivially. We show that while the entropy mode is temporarily excited during the phase transition, the adiabatic mode evolves unscathed, owing primarily to the fact that the axion is an extremely subdominant component of the universe during this period.

In this section, we will present our main derivation in detail. We first derive the equations of motion for both the background and perturbation of the QCD axion DM assuming a general temperature dependent potential and a cosmology dominated by a single fluid with a constant equation of state. In this general setup, we show the exact solution for the axion field perturbation on large scales and identify the adiabatic and isocurvature perturbations. We also comment on non-Gaussianity by computing the bispectrum.

2.1 Equations of Motion↩︎

The action for the dimensionless axion field \(\theta\) with decay constant \(f_a\) is \[S = -\int {\rm d}t\,{\rm d}^3 x \, \sqrt{-g}\left[\frac{1}{2} f_a^2 g^{\mu \nu}\partial_\mu \theta \partial_\nu\theta + V(\theta,T)\right]\,,\] where \(g^{\mu\nu}\) is the metric tensor with determinant \(g\), and \(V(\theta, T)\) is the axion potential. \(V(\theta, T)\) could be quite complicated and depends on the temperature, \(T\), prior to the QCD phase transition. When \(T\) is above the characteristic temperature \(T_{\rm QCD}\) of the QCD phase transition, the mass takes an approximately power-law form \[\label{eq:32axion32m40T41} \tilde{m}(T) \approx m_a\left(\frac{T_\textrm{\tiny QCD}}{T}\right)^n,\tag{4}\] where \(m_a\) is the zero-temperature mass and \(n \approx 4\). When \(T < T_{\rm QCD}\), \(\tilde{m}(T) = m_a\). Lattice gauge theory calculations have determined the form of \(m(T)\) more precisely, confirming that the exponent \(n\) at high temperatures is close to that predicted by the dilute instanton gas approximation [38]. Beyond the quadratic level, the temperature-dependent interaction terms remain to be calculated. These higher-order terms could be important in particular when the initial misalignment angle is large and the axion oscillations explore a wide range of the potential [29], [39].

Let us write the axion field in terms of a (space-independent) background \(\overline{\theta}(t)\) and a perturbation \(\delta\theta(t,\boldsymbol{x})\), \[\theta(t,\boldsymbol{x}) = \overline{\theta}(t)+\delta\theta(t,\boldsymbol{x}) \, .\] We also take the temperature to be perturbed as \(T (t, \boldsymbol{x}) = \overline{T}(t) + \delta T(t, \boldsymbol{x})\) with \(\overline{T}\) as the background value and \(\delta T\) as the perturbation. Then, one can expand the action to quadratic order in the complete set of perturbations \(\{\delta\theta,\delta T,\Phi,\Psi\}\) and obtain the equations of motion for \(\theta\). Collecting terms by order, we arrive at the zeroth-order evolution equation for the background field \(\overline{\theta}(t)\) as \[\ddot{\overline{\theta}} + 3 H \dot{\overline{\theta}} + \frac{1}{f_a^2}\frac{\partial \overline{V}}{\partial\overline{\theta}} = 0 \, , \label{eq:background}\tag{5}\] where \(\overline{V}=V(\overline{\theta},\overline{T})\). It is well known that the background evolution of the QCD axion is sensitive to the detailed form of \(\overline{V}\). The temperature dependence and anharmonic effects could modify the solution of \(\overline{\theta}\) and consequently the predicted relic abundance.

The first-order perturbation equation for the field perturbation \(\delta\theta(t,\boldsymbol{x})\) is \[\delta\ddot{\theta}+3H\delta\dot{\theta}+\left(-\nabla^2 + \frac{1}{f_a^2}\frac{\partial^2 \overline{V}}{\partial{\overline{\theta}}^2}\right)\delta\theta -\dot{\overline{\theta}}(3\dot{\Phi}+\dot{\Psi}) + 2\Psi\frac{1}{f_a^2}\frac{\partial \overline{V}}{\partial\overline{\theta}}+\frac{1}{f_a^2}\frac{\partial^2 \overline{V}}{\partial \overline{\theta}\partial \overline{T}}\delta T=0 \, .\] Alternatively, we can write the Fourier transform of the perturbed equation for \(\delta\theta_{\boldsymbol{k}} (t) = \int d^3 \boldsymbol{x}\, \mathrm{e}^{-i{\boldsymbol{k}} \cdot\boldsymbol{x}} \delta\theta(t,\boldsymbol{x})\). In Fourier space, the perturbation equation becomes \[\delta\ddot{\theta}_{\boldsymbol{k}}+3H \delta\dot{\theta}_{\boldsymbol{k}}+\left(k^2 + \frac{1}{f_a^2}\frac{\partial^2 \overline{V}}{\partial\overline{\theta}^2}\right)\delta\theta_{\boldsymbol{k}} \,-\dot{\overline{\theta}}(3\dot{\Phi}_{\boldsymbol{k}}+\dot{\Psi}_{\boldsymbol{k}})+2\Psi_{\boldsymbol{k}}\frac{1}{f_a^2}\frac{\partial \overline{V}}{\partial\overline{\theta}}+\frac{1}{f_a^2}\frac{\partial^2 \overline{V}}{\partial \overline{\theta}\partial \overline{T}}\delta T_{\boldsymbol{k}} =0 \, , \label{eq:perteq}\tag{6}\] where \(k^2\equiv \boldsymbol{k}^2\); and \(\Phi_{\boldsymbol{k}}\), \(\Psi_{\boldsymbol{k}}\), and \(\delta T_{\boldsymbol{k}}\) are the Fourier modes of \(\Phi\), \(\Psi\) and \(\delta T\) with wavenumber \(k\). The last three terms in the equation above are all source terms for the axion perturbation. The first two source terms, depending on \(\Phi\), \(\Psi\) and their derivatives, originate from the axion field coupling to the metric. They are universally present in perturbation equations of general ALPs, no matter whether the potential is a constant or temperature-dependent. The last term proportional to \(\delta T\) takes into account the interaction between the QCD axion and radiation (e.g., gluons). This term is important and has to be included, when the axion potential is temperature dependent, to get the right adiabatic initial conditions, as will be discussed in the next section.⁵

The linear-order equations are sufficient to predict the adiabatic and isocurvature power spectrum on large scales. However, we will be interested in determining the axion isocurvature bispectrum as well, for which we will need to account for higher order terms in \(\delta\theta_{\boldsymbol{k}}\). We will discuss this in detail below.

2.2 Solution of the Axion Field Perturbation↩︎

We can understand the small-\(k\) solution to the axion perturbation equation 6 quite generally in terms of the family of background solutions \({\bar \theta}(t, \theta_\text{ini}, {\dot{\theta}}_\text{ini})\), which are parametrized in terms of the initial homogeneous axion value \(\theta_\text{ini}\) and its time derivative. Specifically, we find that, to good approximation, the general solution to linear order in perturbations at small \(k\) has the form \[\delta\theta_{\boldsymbol{k}}(t) = A_{\boldsymbol{k}} \frac{\partial\overline{\theta}(t,\theta_\text{ini},{\dot{\theta}}_\text{ini})}{\partial\theta_\text{ini}} + B_{\boldsymbol{k}} \frac{\partial\overline{\theta}(t,\theta_\text{ini},{\dot{\theta}}_\text{ini})}{\partial{\dot{\theta}}_\text{ini}} - \zeta_{\boldsymbol{k}} {\dot{\overline{\theta}}}(t) \frac{1}{R(t)} \int^tR(t')dt' \, . \label{eq:generalanalytic}\tag{7}\] Here \(A_{\boldsymbol{k}}\), \(B_{\boldsymbol{k}}\), and \(\zeta_{\boldsymbol{k}}\) are constant in time. The last term, proportional to \(\zeta_{\boldsymbol{k}}\), follows from Weinberg’s theorem and the corresponding general formula 3 . The first two terms, then, can be viewed as isocurvature contributions. They are solutions to the homogeneous perturbation equation, that is, to equation 6 in the small \(k\) limit and neglecting the sources \(\Phi_{\boldsymbol{k}}\), \(\Psi_{\boldsymbol{k}}\), and \(\delta T_{\boldsymbol{k}}\). Weinberg’s solution, on the other hand, depends crucially on these source terms. These source terms, in turn, are not simply constant but depend on \(\delta\theta_{\boldsymbol{k}}\) through Einstein’s equations. As a result, it is not obvious that the use of Weinberg’s solution 7 remains valid when the isocurvature contributions are nonzero. We argue in Appendix 4 that this is, in fact, a good approximation provided that the axion is a subdominant contribution to the total energy density of the universe.

In the remainder of this subsection, we will first elaborate on the isocurvature solution and its extension to higher order in perturbations (which is necessary to compute non-Gaussianities in the isocurvature modes). For typical cosmological assumptions, we may neglect the \(B_{\boldsymbol{k}}\) terms. Then, we will give a direct demonstration that Weinberg’s adiabatic ansatz solves the equations of motion under simplifying assumptions. For this explicit calculation below, we will assume that the universe is dominated by a single fluid and thus evolves according to a power law growth \[\label{eq:Rpowerlaw} R(t)\propto t^p\tag{8}\] (where for instance, \(p=1/2\) corresponds to radiation domination and \(p=2/3\) to matter domination). This assumption goes beyond the usual assumption of radiation domination and simplifies the calculation, but it is otherwise unnecessary; we comment later on generalizing these results. In addition, though similarly unnecessary, we assume the background is free of anisotropic stress such that \(\Phi=\Psi\), working with \(\Psi\) from now on.

Under these simplifying assumptions, we will examine the following solution for \(\delta\theta_{\boldsymbol{k}}\) as a special case of 7 : \[\delta\theta_{\boldsymbol{k}}(t) = A_{\boldsymbol{k}} \frac{\partial\overline{\theta}(t,\theta_\text{ini})}{\partial\theta_\text{ini}} + t \dot{\overline{\theta}}(t) \Psi_{\boldsymbol{k}} \, . \label{deltathetasoln}\tag{9}\] We will discuss these two terms separately below.

2.2.1 Isocurvature Solution↩︎

First, we will obtain the general homogeneous solution of the axion fluctuations. We will make minimal assumptions about the background solution, which allow our final result to be applicable beyond the standard misalignment mechanism and include, e.g. kinetic misalignment [33], [40].

In general, the axion fluctuations evolve non-linearly due to the non-linear potential, which in turn will produce non-Gaussian statistics for the axion density fluctuations, and therefore also for the isocurvature mode. Note that these non-Gaussianities would arise even if the axion obeyed purely Gaussian statistics by the end of inflation and are therefore distinct from inflationary non-Gaussianities.

In order to obtain the bispectrum, we must first solve the homogeneous equation of motion for the fluctuations (without source terms) beyond linear order. Up to quadratic order we have \[\delta\ddot{\theta} + 3H \delta\dot{\theta}+\frac{1}{f_a^2}\bar{V}'' \delta\theta+ \frac{1}{2 f_a^2} \bar{V}''' \delta\theta^2 =0 \, .\] We can attempt to solve this system perturbatively. A good way to organize the perturbative series is through the initial conditions. To that end, we can write \(\delta\theta= \delta\theta_{\rm ini} \hat{\delta\theta}\), with \(\delta\theta_{\rm ini}\) the primordial axion fluctuation generated during inflation. This parametrization makes it explicit that the non-linear terms are all suppressed due to the small initial fluctuations. Therefore we can write \[\hat{\delta\theta} = \hat{\delta\theta}^{(1)}+ \delta\theta_{\rm ini} \hat{\delta\theta}^{(2)} + \cdots \, ,\] where the initial conditions are enforced at linear order by setting \(\hat{\delta\theta}^{(1)}(t=0)=1\) and \(\dot{\hat{\delta\theta}}^{(1)}(t=0)=0\). The initial conditions for the higher order terms may then be trivial.

At linear order the equation of motion admits a very simple solution in terms of the background solution. It can be verified as \(\hat{\delta\theta}_{\boldsymbol{k}}^{(1)}(t)= \frac{\partial \overline{\theta}}{\partial \theta_{\rm ini}}\). In fact it is the unique solution which satisfies the initial conditions \(\hat{\delta\theta}_{\boldsymbol{k}}(0)=1\) and \(\dot{\hat{\delta\theta}}_{\boldsymbol{k}}(0)=0\).⁶

At second order the equations of motion (in position space) are \[\ddot{\hat{\delta\theta}}^{(2)} + 3H \dot{\hat{\delta\theta}}^{(2)}+\frac{1}{f_a^2}\bar{V}'' \hat{\delta\theta}^{(2)} + \frac{1}{2 f_a^2} \bar{V}''' (\hat{\delta\theta}^{(1)})^2 = 0 \, .\] It can be similarly verified by differentiating (5 ) twice that the unique solution to this equation is \(\frac{1}{2}\frac{\partial^2\overline{\theta}}{\partial \theta_{\rm ini}^2}\). Therefore the complete solution at second order is \[\delta\theta(\boldsymbol{x},t) = \delta\theta_{\rm ini}(\boldsymbol{x})\left[\frac{\partial \overline{\theta}}{\partial \theta_{\rm ini}} + \frac{1}{2}\delta\theta_{\rm ini}(\boldsymbol{x}) \frac{\partial^2 \overline{\theta}}{\partial \theta_{\rm ini}^2}\right] \, .\] Indeed, we can recognize this to be a series expansion in the initial data, and we can therefore obtain the general solution to \(n^{\rm th}\)-order \[\delta\theta(\boldsymbol{x},t) = \sum_{n=1}^\infty \left[\frac{1}{n!}\frac{\partial^n \overline{\theta}}{\partial \theta_{\rm ini}^n}\right] \delta\theta_{\rm ini}^n(\boldsymbol{x}) \, .\] Therefore the complete non-linear solution for the homogeneous, or isocurvature, field perturbations at all orders of the initial fluctuations \(\delta\theta_{\rm ini}(\boldsymbol{x})\) are entirely determined by the background solution. At leading order, we recover the first term in eq. 9 by taking the Fourier transform and identifying \(A_{\boldsymbol{k}} = \delta\theta_{\rm ini}(\boldsymbol{k})\).

2.2.2 Adiabatic Solution↩︎

We now specialize to the power-law scale factor \(R(t) \propto t^p\). The second term in eq. 9 is, in fact, exactly Weinberg’s adiabatic solution. From 2 , we obtain \[\Psi_{\boldsymbol{k}}(t) = \Phi_{\boldsymbol{k}}(t) = -\frac{\zeta_{\boldsymbol{k}}}{1+p} \,\] and in eq. (3 ) \[\delta\theta_{\boldsymbol{k}}(t) = -\frac{\zeta_{\boldsymbol{k}}}{1+p}t \,\dot{\overline{\theta}}(t) = t \dot{\overline{\theta}}(t)\Psi_{\boldsymbol{k}}.\] (Apart from Weinberg’s general result, such solutions have appeared in past literature, e.g., [41], in specific contexts like a quadratic and temperature-independent potential.) Let us verify that the second term, \(\delta\theta_{\boldsymbol{k}}= t\dot{\overline{\theta}}\Psi_{\boldsymbol{k}}\), is indeed a (particular) solution to the full perturbation equation. For superhorizon modes, we neglect the \(k^2\) term and take the metric potentials to be constant \(\dot{\Phi}_{\boldsymbol{k}}=\dot{\Psi}_{\boldsymbol{k}}=0\) (though once more, the latter assumption may be relaxed). Plugging this solution in to the left side of eq. 6 , we have \[(t\dddot{\overline{\theta}}+2\ddot{\overline{\theta}})\Psi_{\boldsymbol{k}} + 3 H\Psi_{\boldsymbol{k}} (t\ddot{\overline{\theta}}+\dot{\overline{\theta}}) + \frac{1}{f_a^2}\frac{\partial^2 {\overline{V}}}{\partial\overline{\theta}^2} (t \dot{\overline{\theta}}\Psi_{\boldsymbol{k}}) +2\Psi_{\boldsymbol{k}}\frac{1}{f_a^2}\frac{\partial {\overline{V}}}{\partial\overline{\theta}}+\frac{1}{f_a^2}\frac{\partial^2 \overline{V}}{\partial \overline{\theta}\partial \overline{T}}\delta T_{\boldsymbol{k}} = 0\, .\] Taking the time derivative of the background equation, eq. 5 , we have \[\dddot{\overline{\theta}}+3H\ddot{\overline{\theta}}+3\dot{H}\dot{\overline{\theta}}+\frac{1}{f_a^2}\left(\frac{\partial^2 \overline{V}}{\partial \overline{\theta}^2}\dot{\overline{\theta}}+\frac{\partial^2 \overline{V}}{\partial \overline{\theta}\partial \overline{T}}\dot{\overline{T}}\right)=0 \, .\] Combining these two expressions above, we find the left-hand side of Eq. 6 to be \[t\Psi_{\boldsymbol{k}}\left(-3\dot{H}\dot{\overline{\theta}}-\frac{1}{f_a^2}\frac{\partial^2 \overline{V}}{\partial \overline{\theta}\partial \overline{T}}\dot{\overline{T}}\right)+2 \Psi_{\boldsymbol{k}} \ddot{\overline{\theta}}+3H \Psi_{\boldsymbol{k}} \dot{\overline{\theta}} +2\Psi_{\boldsymbol{k}}\frac{1}{f_a^2}\frac{\partial \overline{V}}{\partial\overline{\theta}}+\frac{1}{f_a^2}\frac{\partial^2 \overline{V}}{\partial \overline{\theta}\partial \overline{T}}\delta T_{\boldsymbol{k}} \, = 0 \,.\] For \(R(t)\propto t^p\), \(t\dot{H}=-H\). Thus, using the equation of motion for \(\overline{\theta}\) in 5 , we are left only with \[- \frac{1}{f_a^2}\frac{\partial^2 \overline{V}}{\partial \overline{\theta}\partial \overline{T}}\left(t\Psi_{\boldsymbol{k}}\dot{\overline{T}}-\delta T_{\boldsymbol{k}} \right)\, = 0 \,. \label{eq:laststep}\tag{10}\] However, in the adiabatic mode, the perturbation \(\delta T_{\boldsymbol{k}}\) is related to \(\dot{\overline{T}}\) by 3 , which for our background 8 implies that \(\delta T_{\boldsymbol{k}}/\dot{\overline{T}} = t \Psi_{\boldsymbol{k}}\). Thus eq. 10 vanishes and indeed \(t\dot{\overline{\theta}} \Psi_{\boldsymbol{k}}\) is a solution to the full axion field perturbation equation.

Transforming back from the Fourier mode, we could rewrite the adiabatic solution as \(\delta \theta (t, \boldsymbol{x})= t \dot{\overline{\theta}}(t) \Psi(t, \boldsymbol{x})\). We have a few comments on this exact solution at superhorizon scales:

• Since the evolution of \(\overline{\theta}\) depends on \(V\), this particular solution tells us that the axion field perturbation \(\delta \theta\) also varies for different \(V\)’s. However, the proportionality relation between \(\delta \theta\) and \(\dot{\overline{\theta}}\) is independent of \(V\) and serves as an example of Weinberg’s general adiabaticity theorem [24]–[26].

• The more intuitive picture suggested by the solution is as follows. Initially, the QCD axion is at rest with \(\dot{\overline{\theta}}=0\), and no field perturbation is generated. Once the axion field starts moving (even before entering the full coherent oscillation phase), the axion field perturbation starts to see the metric and temperature perturbations (correlated as they both arise from the primordial inflaton fluctuations), via the motion of the field.

• This solution holds in scenarios beyond the standard misalignment mechanism happening during radiation domination. In particular, it is also a solution in a matter dominated universe. Indeed, one highly motivated scenario to achieve high-scale QCD axion DM with \(f_a > 10^{12}\) GeV is to have it start oscillating in an early matter domination phase sourced by a scalar (see, e.g., [42]–[46] for such constructions).

• While our main physics target is QCD axion DM, the derivation above applies to arbitrary ALP and non-thermal scalar DM, no matter what the potential is.

• We assume a single-fluid-dominated universe with a simple power law expansion and constant gravitational perturbations at the superhorizon scale. These are not necessary assumptions. Based on Weinberg’s theorem, the adiabatic solution always guarantees the proportionality between \(\delta \theta_k\) and \(\dot{\overline{\theta}} \Psi\). The proportionality factor could be more complicated than \(t\), if the background evolution is more complicated. We will discuss this further in Appendix 4.

2.3 Adiabatic Energy Density Perturbation of Axion DM↩︎

Let us explore the adiabatic mode in terms of energy and pressure. We turn now to the QCD axion as an example with a concrete potential which depends on the field and the temperature. The energy-momentum tensor for the QCD axion is given by \[T^\mu_\nu = f_a^2\partial^\mu\theta\partial_\nu\theta - \frac{1}{2}\left[f_a^2\partial^\lambda\theta\partial_\lambda\theta+2V(\theta,T)\right]\delta^\mu_\nu.\] From the \(00^{\rm th}\) component, we can expand the energy density of the fluid to the leading two orders and find the background (\(\overline{\rho}_\theta\)) and first-order perturbation (\(\delta \rho_\theta\)) of the energy density for axion DM to be \[\begin{align} \overline{\rho}_\theta &= \frac{f_a^2}{2} \dot{\overline{\theta}}^2 + \overline{V}, \label{eq:rho32theta} \\ \delta\rho_\theta &= f_a^2(\dot{\overline{\theta}}\delta\dot{\theta}-\dot{\overline{\theta}}^2\Psi)+\frac{\partial \overline{V}}{\partial\overline{\theta}}\delta\theta+\frac{\partial \overline{V}}{\partial\overline{T}}\delta T. \end{align}\tag{11}\] Then, using the adiabatic solution for the axion field perturbation \(\delta\theta=t\Psi\dot{\overline{\theta}}\), we find that the energy density perturbation \(\delta \rho_\theta\) is \[\label{eq:adiabatic32delta32rho} \delta \rho_\theta=-3Hf_a^2\dot{\overline{\theta}}^2 t \Psi+\frac{\partial \overline{V}}{\partial\overline{T}}\delta T\,.\tag{12}\] The time derivative of the axion background energy density is \[\dot{\overline{\rho}}_\theta = f_a^2 \dot{\overline{\theta}} \,\ddot{\overline{\theta}} + \frac{{\rm d}\overline{V}}{{\rm d}t} = f_a^2 \dot{\overline{\theta}} \left[-3H \dot{\overline{\theta}} - \frac{1}{f_a^2}\frac{\partial \overline{V}}{\partial \overline{\theta}}\right] + \dot{\overline{\theta}} \frac{\partial \overline{V}}{\partial \overline{\theta}} + \dot{\overline{T}} \frac{\partial \overline{V}}{\partial \overline{T}}\,,\] where using \(\dot{\overline{T}}=-p \overline{T}/t = - H \overline{T}\), we have \[\dot{\overline{\rho}}_\theta = -3Hf_a^2\dot{\overline{\theta}}^2- H \overline{T}\frac{\partial \overline{V}}{\partial \overline{T}} \, .\] Since \(H=p/t\) and \(\delta T/\overline{T}= - p \Psi\), the generalized curvature perturbation for the axion is \[\zeta_\theta=-\Psi- H\frac{\delta \rho_\theta}{ \dot{\overline{\rho}}_\theta} =-\Psi- H \Psi \frac{-3p f_a^2 \dot{\overline{\theta}}^2 - p \overline{T}\frac{\partial \overline{V}}{\partial \overline{T}}}{-3Hf_a^2\dot{\overline{\theta}}^2- H \overline{T}\frac{\partial \overline{V}}{\partial \overline{T}}} = -(p+1) \Psi \, .\] On the other hand, for radiation \(\delta \rho_\gamma/\overline{\rho}_\gamma = -4 p \Psi\) and \(\dot{\overline{\rho}}_\gamma = -4 H \overline{\rho}_\gamma\). Thus \[\zeta_\gamma = - \Psi- H\frac{\delta \rho_\gamma}{\dot{\overline{\rho}}_\gamma} = -(p+1) \Psi= \zeta_\theta \, .\] With \(\zeta_\gamma=\zeta_\theta\), we have \(S_{\theta\gamma}=0\), conclusively showing that \(\delta \theta = t \Psi\dot{\overline{\theta}}\) is indeed the adiabatic mode.

Note that for many standard cosmological fluids, \(\dot{\overline{\rho}}\propto H \overline{\rho}\), such that the ratio \(\delta\rho/\dot{\overline{\rho}}\) is proportional to \(\delta\rho/\overline{\rho}\), and thus often the adiabatic nature of initial conditions for a given fluid is stated as a condition on \(\delta\rho/\overline{\rho}\). In particular, if a component \(\rho_i\) is independently conserved and has a definite equation of state \(w_i\), then we have \(\dot{\overline{\rho}}_i = -3H(1+w_i) \overline{\rho}_i\). For example, in standard \(\Lambda\)CDM cosmology, adiabatic fluctuations for CDM are initially given by the density perturbation \(\delta_\text{cdm} \equiv \delta\rho_\text{cdm}/\overline{\rho}_\text{cdm}=-3\Psi/2\), while the fluctuation in photon energy density is \(\delta_\gamma=-2\Psi\) during radiation domination. However, when the background energy density evolves nontrivially such that \(\dot{\overline{\rho}} \not\propto H \overline{\rho}\), then one must examine the relationship between \(\delta\rho\) and \(\dot{\overline{\rho}}\) to determine the adiabaticity of a solution.

In the preceding discussion, we have shown analytically that the adiabaticity of the solution is preserved even when the QCD axion has a temperature-dependent potential. It is useful to examine the behavior of the axion density perturbation \(\delta_\theta\equiv \delta \rho_\theta/\overline{\rho}_\theta\), to see that while it remains adiabatic, its behavior differs from thermal CDM at early times. 1 shows both a numerical and analytical solution for \(\delta_\theta\). For both cases, high-frequency oscillations are filtered out such that the figure displays oscillation-averaged \(\langle\delta_\theta\rangle\). To implement a numerical computation, we assume the QCD axion potential to be quadratic with a temperature-dependent mass given in eq. 4 , until \(t=t_{\rm QCD}\) when the mass becomes constant. Further, we assume radiation domination (\(p=1/2\)) and superhorizon conditions (\(k \ll R H |_{t_{\rm QCD}}\)) for the duration of the numerical solution. The analytical solution in the figure corresponds to \(\delta\theta_{\boldsymbol{k}}(t) = t \dot{\overline{\theta}}(t) \Psi_{\boldsymbol{k}} \,\), with \(\dot{\overline{\theta}}(t)\) obtained using the numerical solution for the background \(\overline{\theta}(t)\).

As shown in 1, initially the density perturbation is a constant during the vacuum-dominated phase of the axion’s evolution. This differs from a vacuum-dominated scalar with constant potential energy density (which would have \(\delta_\theta = 0\)), as the axion’s mass is temperature dependent and thus the potential energy density is time dependent. Combining 12 and 11 , and taking \(\dot{\overline{\theta}} = 0\), the initial constant is simply \(\delta_\theta = \frac{\partial V}{\partial \overline{T}} \frac{\delta T}{V} =n\Psi_{\boldsymbol{k}}\), where the last equality is valid for the power-law in temperature corresponding to 4 , quadratic axion potential, and radiation domination. During this phase, \(\delta_\theta\) need not be a constant if the temperature-derivative of the potential is not proportional to the potential itself; indeed for the QCD axion, the potential’s temperature dependence could change at early times. Following the onset of oscillations, \(\delta_\theta\) oscillates about a constant which differs from the initial condition for \(\delta_\text{cdm}\). Again from 12 and 11 , we can see that the initial constant should be \(\langle \delta_\theta\rangle=-(3-n)\Psi/2\). Then, after the QCD phase transition is over at \(t/t_\text{QCD}=1\) in the figure, the axion potential becomes temperature-independent and it begins to behave as CDM, with the averaged density perturbation \(\langle\delta_\theta\rangle=\delta_\text{cdm}=-3\Psi/2\).

2.4 Isocurvature Energy Density Perturbation of Axion DM↩︎

Let us now return to the homogeneous solution for \(\delta\theta\), which is responsible for the isocurvature component of the DM density fluctuations. Since we will be interested in determining both the isocurvature power spectrum and bispectrum, we will need to determine the axion density fluctuations to second order in the initial conditions.

One way to obtain the density fluctuations beyond linear order is to use the solution of the axion field perturbations to higher order in the initial data. A shortcut is to recognize that the total energy density is a function of \(\theta_{\rm ini}+ \delta\theta_{\rm ini}\), which means that non-zero initial fluctuations provide a different misalignment angle at each point in space. Therefore one may simply Taylor expand in \(\delta\theta_{\rm ini}(\boldsymbol{x})\) to obtain \[\label{eq:axion95iso95fluct} \delta\rho_\theta^{\rm iso}(\boldsymbol{x},t) = \frac{\partial \bar{\rho}_\theta}{\partial \theta_{\rm ini}} \delta\theta_{\rm ini}(\boldsymbol{x}) + \frac{1}{2} \frac{\partial^2 \bar{\rho}_\theta}{\partial \theta_{\rm ini}^2} \delta\theta^2_{\rm ini}(\boldsymbol{x}) + \cdots \, .\tag{13}\] One may verify that the field perturbation \[\delta\theta(\boldsymbol{x},t) = \delta\theta_{\rm ini}(\boldsymbol{x})\left[\frac{\partial \overline{\theta}(t)}{\partial \theta_{\rm ini}} + \frac{1}{2} \delta\theta_{\rm ini}(\boldsymbol{x})\frac{\partial^2 \overline{\theta}(t)}{\partial \theta_{\rm ini}^2}\right]\,\] reproduces (13 ). To check this, we can use the fact that at linear order \[\delta \rho_\theta^\text{iso,(1)} = f_a^2\dot{\overline{\theta}}\left(A_{\boldsymbol{k}} \frac{\partial^2 \overline{\theta}}{\partial \theta_{\rm ini} \partial t}\right) + \frac{\partial \overline{V}}{\partial \overline{\theta}}A_{\boldsymbol{k}} \frac{\partial \overline{\theta}}{\partial \theta_{\rm ini}} = A_{\boldsymbol{k}} \frac{\partial \overline{\rho}_\theta}{\partial \theta_{\rm ini} } = \frac{\partial \overline{\rho}_\theta}{\partial \theta_{\rm ini} }\delta\theta_\text{ini} \, .\] Likewise, one can compute the density fluctuation to quadratic order as well.

Now let us determine the entropy perturbations. At linear order we have \[\begin{align} S_{\theta \gamma}(\boldsymbol{k}) &\equiv 3\left[\zeta_{\theta}-\zeta_{\gamma}\right] = -3 H \frac{\delta \rho_{\theta}^{\rm iso}}{\dot{\bar{\rho}}_\theta} \\ &=-\frac{3H}{\dot{\bar{\rho}}_\theta}\frac{\partial \bar{\rho}_\theta}{\partial \theta_{\rm ini}} \delta\theta_{\rm ini}(\boldsymbol{k}) \, . \end{align}\] In the absence of temperature dependent effects, after the background oscillations of the axion have begun we have \(\dot{\bar{\rho}}_\theta = -3 H \bar{\rho}_\theta\) which yields the simple expression \[\label{eq:energyisopert} S_{\theta \gamma} = \frac{\partial \ln \overline{\rho}_\theta}{\partial \theta_{\rm ini} }\delta\theta_\text{ini} = \frac{\partial \ln \Omega_\theta}{\partial \theta_{\rm ini} }\frac{H_I}{2\pi f_a} \, ,\tag{14}\] where in the last step, we use \(\frac{\partial \ln \overline{\rho}_\theta}{\partial \theta_{\rm ini} } =\frac{\partial \ln \Omega_\theta}{\partial \theta_{\rm ini} }\), with \(\Omega_\theta\) being the relic abundance of axion dark matter. Strictly speaking, this equality holds only after the axion starts to oscillate in the quadratic part of the potential near the minimum. That is indeed the case for the cosmologically observable modes at large scales. We also use that the primordial axion field fluctuation is set by the inflationary Hubble scale \(H_I\): \(\delta \theta_{\rm ini} = \frac{H_I}{2\pi f_a}\). To convert this further to CDM isocurvature, one needs to multiply it by the fraction of axion dark matter in the total amount of dark matter, \(r S_{\theta \gamma}\) with \(r=\Omega_\theta/\Omega_c\).

Note that eq. 14 had been derived before [29], [47], using the \(\Delta N\) formalism and during radiation domination, among other assumptions. We have shown that this form of the isocurvature mode is more general and follows immediately from the background evolution of the axion. As pointed out in [47], in the large misalignment limit with large \(\theta_{\rm ini} \to \pi\), \(\delta \rho_\theta^{\rm iso}\) will be enhanced since the axion energy density is more sensitive to the initial misalignment angle in this hilltop region.

However, accounting for the temperature dependent potential and keeping the axion equation of state arbitrary complicates this expression. Considering the quadratic potential of the QCD axion, we have that \(V \propto T^{-2n}\), which yields \[\begin{align} S_{\theta \gamma} &= \left(f_a^2 \dot{\overline{\theta}}^2 + \frac{1}{3}\overline{T}\frac{\partial \overline{V}}{\partial \overline{T}}\right)^{-1}\left(\frac{\partial \bar{\rho}_\theta}{\partial \theta_{\rm ini}} \right)\delta\theta_{\rm ini} \\ &= \frac{1}{(1-\tfrac{n}{3})+w_\theta(1+\tfrac{n}{3})}\left(\frac{\partial \log\bar{\rho}_\theta}{\partial \theta_{\rm ini}} \right)\delta\theta_{\rm ini} \\ &\approx \frac{2}{1-\tfrac{n}{3}} \frac{\delta\theta_{\rm ini}}{\theta_{\rm ini}} \, , \end{align}\] where in the last line we have used the fact that for a quadratic potential \(\bar{\rho}_\theta \propto \theta_{\rm ini}^2\), and taken \(w_\theta \approx 0\).⁷ Choosing, for instance, \(n=4\) gives an entropy perturbation \(S_{\theta\gamma}=-6\frac{\delta \theta_{\rm ini}}{\theta_{\rm ini}}\), three times larger (and negative) prior to the QCD phase transition when the temperature-dependence is significant, compared to later times when \(S_{\theta\gamma}=2\frac{\delta \theta_{\rm ini}}{\theta_{\rm ini}}\).

While this time-dependent shift in the entropy perturbation is of theoretical interest, this behavior is only present at such early times that render it irrelevant for most scales of observational interest (for example the scales that impact the CMB or large-scale-structure observations). Instead, one could readily conceive of scenarios where similar behavior can become observationally relevant. For instance, a dark sector with a phase transition which occurs at a time significantly later than the QCD phase transition might exhibit this shift in entropy perturbations which may impact the CMB power spectra differently than typical constant isocurvature.

2.5 Non-Gaussianity↩︎

Lastly, we will briefly comment on the isocurvature bispectrum. In general, the axion fluctuations evolve non-linearly due to the non-linear potential, which in turn will produce non-Gaussian statistics for the axion density fluctuations, and therefore also for the isocurvature mode [48]. These non-Gaussianities would arise even if the axion obeyed purely Gaussian statistics by the end of inflation and are therefore distinct from inflationary non-Gaussianities. For this calculation, we will turn off the temperature and metric fluctuations as our goal is to calculate the pure isocurvature bispectrum.⁸

To compute the bispectrum we must first determine the comoving curvature of the axion fluid to second order in the density fluctuations which is⁹ \[\zeta_\theta^{(2)} = - H \frac{\delta \rho_{\theta}^{(2)}}{\dot{\bar{\rho}}_\theta} +\frac{1}{2}\left(\frac{\dot{H}}{\dot{\bar{\rho}}_{\theta}^2}- H \frac{\ddot{\bar{\rho}}_{\theta}}{\dot{\bar{\rho}}_\theta^3}\right)\left(\delta\rho_\theta^{(1)}\right)^{2} + \frac{H}{\dot{\bar{\rho}}_\theta} \delta\rho_\theta^{(1)}\delta\dot{\rho}_\theta^{(1)} \, ,\] where the superscripts indicate the order at which \(\delta\theta_{\rm ini}\) appears. The second and third terms are more complicated, but can be simplified after some manipulations to yield \[\frac{1}{2}\left(\frac{\dot{H}}{\dot{\bar{\rho}}_{\theta}^2}- H \frac{\ddot{\bar{\rho}}_{\theta}}{\dot{\bar{\rho}}_\theta^3}\right)\left(\delta\rho_\theta^{(1)}\right)^{2} + \frac{H}{\dot{\bar{\rho}}_\theta} \delta\rho_\theta^{(1)}\delta\dot{\rho}_\theta^{(1)} = \frac{1}{2 \dot{\bar{\rho}}_\theta}\partial_t\left[\frac{H}{\dot{\bar{\rho}}_\theta}\left(\frac{\partial \bar{\rho}_\theta}{\partial \theta_{\rm ini}}\right)^2\right]\delta\theta_{\rm ini}^2 \, .\]

In general, as we have previously seen due to temperature dependent effects and the changing axion equation of state, \(\dot{\bar{\rho}}_\theta \neq -3 H \bar{\rho}_\theta\), so this expression must be evaluated carefully. In total the curvature perturbation is then \[\zeta_\theta = -\frac{H}{\dot{\bar{\rho}}_\theta}\left(\frac{\partial \bar{\rho}_\theta}{\partial \theta_{\rm ini}} \right)\delta\theta_{\rm ini} - \frac{1}{2}\left[\frac{H}{\dot{\bar{\rho}}_\theta}\frac{\partial^2 \bar{\rho}_\theta}{\partial \theta_{\rm ini}^2} - \frac{1}{\dot{\bar{\rho}}_\theta} \partial_t\left[\frac{H}{\dot{\bar{\rho}}_\theta}\left(\frac{\partial \bar{\rho}_\theta}{\partial \theta_{\rm ini}}\right)^2\right]\right]\delta\theta_{\rm ini}^2 \, .\] During QCD phase transition, the bispectrum will also acquire corrections due to the temperature dependent potential, which quickly become irrelevant as the phase transition ends. Therefore, in order to predict the bispectrum relevant for CMB observations, let us assume that the potential is temperature independent and the axion has begun oscillating about the minimum of its potential, and therefore has equation of state \(w_\theta=0\). This expression could then be simplified to \[\frac{1}{2 \dot{\bar{\rho}}_\theta}\partial_t\left[\frac{H}{\dot{\bar{\rho}}_\theta}\left(\frac{\partial \bar{\rho}_\theta}{\partial \theta_{\rm ini}}\right)^2\right] =- \frac{1}{6}\left(\frac{1}{\bar{\rho}_\theta}\frac{\partial \bar{\rho}_\theta}{\partial \theta_{\rm ini}}\right)^2 \, ,\] and the comoving curvature is \[\label{eq:zeta95axion952nd95gen} \zeta_\theta = \frac{1}{3}\frac{\partial \ln \bar{\rho}_\theta}{\partial \theta_{\rm ini}} \delta\theta_{\rm ini} + \frac{1}{6}\left[\frac{1}{\bar{\rho}_\theta}\frac{\partial^2 \bar{\rho}_\theta}{\partial \theta_{\rm ini}^2}-\left(\frac{\partial \ln \bar{\rho}_\theta}{\partial \theta_{\rm ini}}\right)^2 \right]\delta\theta_{\rm ini}^2 \, .\tag{15}\] The DM isocurvature mode \(S_{c\gamma}\) up to second order in the axion fluctuations, upon multiplying with the appropriate factor of the relative axion density \(r\) is then \[S_{c \gamma} = 3 r \zeta_\theta = r \frac{\partial \ln \bar{\rho}_\theta}{\partial \theta_{\rm ini}} \delta\theta_{\rm ini} + \frac{r}{2}\left[\frac{1}{\bar{\rho}_\theta}\frac{\partial^2 \bar{\rho}_\theta}{\partial \theta_{\rm ini}^2}-\left(\frac{\partial \ln \bar{\rho}_\theta}{\partial \theta_{\rm ini}}\right)^2 \right]\delta\theta_{\rm ini}^2 \, ,\] The bispectrum can then be readily calculated to be of the local shape \[B_S(k_1, k_2, k_3) = f_{\rm NL}^{(S)}\left[P_S(k_1)P_S(k_2) + {\rm cyc.}\right] \, ,\] where \(P_s\)’s are power spectra, cyc. stand for all possible cyclic permutations, and the amplitude is \[f_{\rm NL}^{(S)} = -1 + \frac{1}{r \bar{\rho}_\theta}\left(\frac{\partial \ln \bar{\rho}_\theta}{\partial \theta_{\rm ini}}\right)^{-2}\frac{\partial^2 \bar{\rho}_\theta}{\partial \theta_{\rm ini}^2} \, ,\] thus reproducing the expression obtained in [47].

3 Conclusions↩︎

In this article, we aim at presenting a clear general derivation for the generation of cosmic perturbations of QCD axion DM in the pre-inflationary scenario. The perturbations consist of two parts, the adiabatic one and the isocurvature one. The most striking feature is that irrespective of the complexity of the axion potential and the cosmic background evolution, the two types of perturbations in the axion field are entirely determined by the family of axion zero-mode solutions with different initial conditions.

For the adiabatic perturbation, we apply Weinberg’s theorem on the adiabatic mode to show that there is always an exact solution in which the axion field (or energy density) perturbation is proportional to the time derivative of axion background field (or energy density). For the isocurvature part, we show that the axion field perturbation could be different from the initial perturbation generated during inflation. This results in an enhancement of the late-time energy density fluctuation, as well as non-Gaussianity, by the sensitivity of the DM abundance to the initial conditions, which can be large when the axion starts oscillating near the hilltop of its potential. This generalizes earlier results [29], [47] to generic background cosmologies (e.g., early matter domination).

While we focus on the QCD axion, our formalism applies to generic ALP or scalar DM, and could aid in understanding the cosmic evolution of non-thermal DM in more complicated scenarios. In addition, some results such as the non-trivial evolution of isocurvature perturbations for a temperature-dependent potential may motivate new cosmological signatures from a dark QCD phase transition.

Acknowledgements↩︎

We thank Adrienne Erickcek, Paddy Fox, Akshay Ghalsasi, Keisuke Harigaya, Eiichiro Komatsu, Lingfeng Li, Andrew Long, Doddy Marsh, and Sarunas Verner for useful discussions. IJA and JF are supported by the NASA grant 80NSSC22K081 and the DOE grant DE-SC-0010010. PC and MR are supported in part by the DOE Grant DE-SC0013607. This work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-2210452. JF also acknowledges support of the Institut Henri Poincaré (UAR 839 CNRS-Sorbonne Université), and LabEx CARMIN (ANR-10-LABX-59-01), where part of the work was performed. PC receives support from the Max Planck–IAS–NTU Center, partly funded by NSTC Grant No. 1142923-M-002-011-MY5. Part of this work was conducted using computational resources and services at the Center for Computation and Visualization, Brown University.

4 Validity of a Constant \(\Psi\) at Superhorizon Scale↩︎

In the main text, we have identified the adiabatic solution of the axion-radiation system and verified that it holds through non-trivial effects, such as the QCD phase transition, under the mild assumption of a power-law background cosmology. If the number of effective degrees of freedom \(g_*\) varies, the power-law background assumption does not hold, and then one needs to rely on the general adiabatic solution in eq. 7 . Another assumption we have made is that the radiation fluid is unaffected by the axion fluid, and the Newton potential is primarily set by the radiation temperature fluctuations. Here we will verify that this is a good assumption during the QCD phase transition, due to the fact that the axion is an extremely subdominant component of the cosmological fluid and because axion isocurvature fluctuations are observed to be suppressed compared to the adiabatic mode.

In Sec. 2, we work with a power-law expansion history, \(R(t) \propto t^p\), and constant gravitational potentials at superhorizon scales. These are valid approximations in a single-fluid dominated universe [26]. Despite corrections from \(g_*\), and neglecting the effect of axion isocurvature, we have seen that there is still an exact adiabatic solution for both the gravitational potential and the axion field perturbation, in (2 ) and (3 ) respectively. Recall that to leading order, in a single-fluid dominated universe with \(R(t) \propto t^p\), the solution above just reduces to \(\Phi_{\boldsymbol{k}}=\Psi_{\boldsymbol{k}}=- \zeta_{\boldsymbol{k}}/(1+p)\), i.e. a constant and \(\delta \theta_{\boldsymbol{k}}/\dot{\overline{\theta}} = t \Psi_{\boldsymbol{k}}\), where \(\zeta_{\boldsymbol{k}}\) is the constant curvature perturbation, as in Sec. 2.

Now we want to comment on the correction to a constant gravitational potential from the isocurvature perturbation of QCD axion dark matter. Take \(\Phi_{\boldsymbol{k}}\) as an example. The equation of motion for \(\Phi_{\boldsymbol{k}}\) at superhorizon scales is \[6 H \dot{\Phi}_{\boldsymbol{k}} + 6 H^2 \Phi_{\boldsymbol{k}} = -8 \pi G \left[\overline{\rho}_\gamma(t)\delta_{\gamma,\boldsymbol{k}} + \overline{\rho}_\theta(t)\delta_{\theta,\boldsymbol{k}}\right] \, ,\] where \(G\) is the Newton constant, \(\delta_\gamma \equiv \delta \rho_\gamma/\overline{\rho}_\gamma\) and \(\delta_\theta = \delta \rho_\theta/\overline{\rho}_\theta\), and \(\mathcal{O}(k^2)\) terms are ignored as usual. We will write the solution as \(\Phi_{\boldsymbol{k}} = - \frac{1}{2}\delta_{\gamma,\boldsymbol{k}} + \delta\Phi_{\boldsymbol{k}}\), where \(- \frac{1}{2}\delta_{\gamma,\boldsymbol{k}}\) is the leading solution and a constant at superhorizon scale in Newtonian gauge, and \(\delta\Phi\) is the correction due to the axion as a subdominant component during radiation domination. At leading order, \(3H^2(t) \approx 8\pi G \overline{\rho}_\gamma(t)\). This yields \[\delta\dot{\Phi}_{\boldsymbol{k}} = -H(t) \delta\Phi_{\boldsymbol{k}} -\frac{1}{2}H(t) \frac{\overline{\rho}_{\theta}}{\overline{\rho}_\gamma} \delta_{\theta, \boldsymbol{k}} \Rightarrow \frac{{\rm d}}{{\rm d}t}(R\, \delta \Phi_{\boldsymbol{k}}) = -\frac{1}{2}\dot{R} \frac{\overline{\rho}_\theta}{\overline{\rho}_\gamma} \delta_{\theta,\boldsymbol{k}}\, .\] The ratio of axion energy density and radiation energy density at an early time, when the scale factor \(R \lesssim R_{\rm osc}\), the scale factor when axion oscillations become cold dark matter, is \[\frac{\overline{\rho}_\theta(R)}{\overline{\rho}_\gamma(R)} = \frac{\overline{\rho}_\theta(R_{\rm osc})}{\overline{\rho}_{\gamma}(R_{\rm osc})} \frac{\overline{\rho}_\gamma (R_{\rm osc})}{\overline{\rho}_\gamma(R)} \frac{\overline{\rho}_\theta(R)}{\overline{\rho}_\theta(R_{\rm osc})} \approx \left(\frac{\Omega_\theta}{\Omega_\gamma}\right)R_{\rm osc} \left(\frac{R}{R_{\rm osc}}\right)^{4} \ll 1 \, ,\] where \(\Omega_{i=\theta, r}\) is the abundance of each species today and we take \(\frac{\overline{\rho}_\theta(R)}{\overline{\rho}_\theta(R_{\rm osc})} \sim 1\). For \(R > R_{\rm osc}\), the ratio becomes \[\frac{\overline{\rho}_\theta(R)}{\overline{\rho}_\gamma(R)} = \left(\frac{\Omega_\theta}{\Omega_\gamma}\right)R \, ,\] which is due to the scaling of energy density for matter and radiation. Thus we could estimate that if \(R \lesssim R_{\rm osc}\), \[\delta\Phi_{\boldsymbol{k}} \lesssim \left(\frac{\Omega_{\theta}}{\Omega_\gamma}\right) \left(\frac{R}{R_{\rm osc}}\right)^3 R \, \delta_{\theta,\boldsymbol{k}} \, ,\] and if \(R>R_{\rm osc}\) we have \[\delta\Phi_{\boldsymbol{k}} \sim \left(\frac{\Omega_{\theta}}{\Omega_\gamma}\right) R \, \delta_{\theta, \boldsymbol{k}} \, .\] Let us consider the latter case since it sees a larger correction of \(\delta \Phi_{\boldsymbol{k}}\). The scale factor can be approximated as the ratio of the temperatures at recombination and QCD phase transition, \(R = T_{\rm rec}/T_{\rm QCD} \sim (1\,{\rm eV}/100 \,{\rm MeV})\sim 10^{-8}\), and the ratio of abundances today is at most \(\Omega_{\theta}/\Omega_\gamma \sim 10^{3}\). When \(\delta_\theta\) is dominated by the inflationary isocurvature contribution, we have \[\frac{\delta\Phi}{\Phi} \sim \left(\frac{\Omega_{\theta}}{\Omega_\gamma}\right) 10^{-8} \frac{\delta\theta_{\rm ini}}{\Phi \theta_{\rm ini}} \sim \left(\frac{\Omega_{c}}{\Omega_\gamma}\right) 10^{-8} \left(\frac{\Omega_{\theta}}{\Omega_c}\right)\frac{\delta\theta_{\rm ini}}{\Phi \theta_{\rm ini}}\sim 10^{-5} \sqrt{\frac{A_{\rm iso}}{A_s}} \, ,\] where we dropped the subscript \(\boldsymbol{k}\), \(\left(\frac{\Omega_{\theta}}{\Omega_c}\right)\delta\theta_i/\theta_i \sim \sqrt{A_{\rm iso}}\) with \(\Omega_c\) the cold dark matter relic abundance today and \(A_{\rm iso}\) the amplitude of primordial isocurvature, and \(\Phi \sim \sqrt{A_s}\) with \(A_s\) the amplitude of the primordial curvature perturbation. Since the isocurvature mode is constrained to be much smaller than the adiabatic mode [21], we have roughly \(\delta \Phi/\Phi \lesssim 10^{-6}\). This shows that it is a very good approximation to neglect the effects of axion isocurvature onto the adiabatic mode.

References↩︎