THIRD-ORDER PERTURBATION THEORY WITH NONLINEAR PRESSURE

Pressure plays an important role for the structure formation in the universe. Pressure determines the Jeans scale, λ_J, below which the growth of structure slows down, and eventually stops and oscillates: while fluctuations in the cold dark matter (CDM) and the pressured component evolve in the same way above the Jeans scale, their evolutions are significantly different below the Jeans scale.

The dominant source of gravity is CDM, which is cold and its velocity dispersion is negligible before the collapse of halos. However, the subdominant matter components—baryons and neutrinos—have significant velocity dispersions, which should be included in the calculation when precision is required. While the accurate calculations have been done for the linear perturbations, the effects of the pressure on the nonlinear evolution of matter fluctuations on cosmological scales (∼10–100 Mpc) have not been studied very much in the literature.

We address this issue by calculating the nonlinear matter power spectrum using the third-order perturbation theory (3PT; see Bernardeau et al. 2002, for a review), with the pressure gradient term in the Euler equation explicitly included. This enables us to study the effects of the pressure on the nonlinear evolution of matter fluctuations in a self-consistent manner.

The rest of this paper is organized as follows. In Section 2, we find the linear, second-order, and third-order solutions of the coupled continuity, Euler, and Poisson equations for two matter components with and without the pressure gradient. In Section 3, we calculate the nonlinear matter power spectrum from the solutions obtained in Section 2. In Section 4, we compare our full 3PT calculation with the approximation used by Saito et al. (2008) for the effects of massive neutrinos on the matter power spectrum. Finally, in Section 5, we discuss the implications of our results for a few practical astrophysical and cosmological applications. In the appendices we give the detailed derivations of the 3PT results used in the main body of the paper.

2.1. Basic Equations

The main goal of this paper is to find the perturbative solutions for the CDM density contrast, δ_c, for which the pressure gradient is ignored, and the density contrast of another matter component, δ_b, for which the pressure gradient is retained. This component may be identified with baryons (hence the subscript "b") or neutrinos, depending on the sound speed one uses in the Euler equation.¹

The equations that we are going to solve include two continuity equations

$$\begin{eqnarray} &&\dot{\delta }_c(\mathbf {x},\tau)+\mathbf {\nabla }\cdot [(1+\delta _c(\mathbf {x},\tau))\mathbf {v}_c(\mathbf {x},\tau)]=0, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} &&\dot{\delta }_b(\mathbf {x},\tau)+\mathbf {\nabla }\cdot [(1+\delta _b(\mathbf {x},\tau))\mathbf {v}_b(\mathbf {x},\tau)]=0, \end{eqnarray} \tag{ 2 }$$

two Euler equations

$$\begin{eqnarray} &&\dot{\mathbf {v}}_c(\mathbf {x},\tau)+[\mathbf {v}_c(\mathbf {x},\tau)\cdot \mathbf {\nabla }]\mathbf {v}_c(\mathbf {x},\tau) =-\frac{\dot{a}}{a}\mathbf {v}_c(\mathbf {x},\tau)-\mathbf {\nabla }\phi (\mathbf {x},\tau),\nonumber\\ \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} &&\dot{\mathbf {v}}_b(\mathbf {x},\tau)+[\mathbf {v}_b(\mathbf {x},\tau)\cdot \mathbf {\nabla }]\mathbf {v}_b(\mathbf {x},\tau) \nonumber \\ &&=-\frac{\dot{a}}{a}\mathbf {v}_b(\mathbf {x},\tau)-\mathbf {\nabla }\phi (\mathbf {x},\tau) - \frac{c_s^2(\mathbf {x},\tau)\mathbf {\nabla }\delta _b(\mathbf {x},\tau)}{1+\delta _b(\mathbf {x},\tau)}, \end{eqnarray} \tag{ 4 }$$

and one Poisson equation

$$\begin{eqnarray} &&\nabla ^2\phi (\mathbf {x},\tau)=4\pi Ga^2[\bar{\rho }_c(\tau)\delta _c(\mathbf {x},\tau)+\bar{\rho }_b(\tau)\delta _b(\mathbf {x},\tau)], \end{eqnarray} \tag{ 5 }$$

where $\delta _i\equiv (\rho _i-\bar{\rho }_i)/\bar{\rho }_i$ is the density contrast of a matter component i = (c, b), $\bar{\rho }$ is the background matter density, a is the scale factor, v_i is the peculiar velocity field of a matter component i, ϕ is the gravitational potential, and c_s is the sound speed of the matter component with pressure. Here, the dots denote the partial derivatives with respect to the conformal time, τ, i.e., $\dot{\delta }=\partial \delta /\partial \tau$, and ∇ denotes the partial derivatives with respect to the comoving coordinates.

We rewrite the Poisson equation as

$$\begin{eqnarray} &&\nabla ^2\phi (\mathbf {x},\tau)=\frac{6}{\tau ^2}\delta (\mathbf {x},\tau), \end{eqnarray} \tag{ 6 }$$

where we have assumed an Einstein–de Sitter (EdS) universe (we shall generalize the results to other cosmological models later), for which the energy density of the universe is dominated entirely by the matter density, and a ∝ τ². The background Friedmann equation is given by

$$\begin{equation} \frac{8\pi G}{3}[\bar{\rho }_c(\tau)+\bar{\rho }_b(\tau)]a^2 = \frac{4}{\tau ^2}. \end{equation} \tag{ 7 }$$

We have also defined the total matter fluctuation, δ, which is given by

$$\begin{eqnarray} \delta (\mathbf {x},\tau)&&\equiv \frac{\bar{\rho }_c(\tau)\delta _c(\mathbf {x},\tau)+\bar{\rho }_b(\tau)\delta _b(\mathbf {x},\tau)}{\bar{\rho }_c(\tau)+\bar{\rho }_b(\tau)}\nonumber\\ &&=f_c\delta _c(\mathbf {x},\tau)+f_b\delta _b(\mathbf {x},\tau), \end{eqnarray} \tag{ 8 }$$

where $f_c\equiv \bar{\rho }_c/(\bar{\rho }_c+\bar{\rho }_b)=\Omega _c/\Omega _m$, and $f_b\equiv \bar{\rho }_b/(\bar{\rho }_c+\bar{\rho }_b)=\Omega _b/\Omega _m$. For an EdS universe, Ω_m = 1.

Taking the divergence of the Euler equations, we obtain the equations for the velocity divergence fields, θ_i ≡ ∇ · v_i. Moving nonlinear terms to the right-hand side of the equations and using the Poisson equation, we obtain

$$\begin{eqnarray} &&\dot{\delta }_c(\mathbf {x},\tau)+\theta _c(\mathbf {x},\tau) =-\mathbf \nabla \cdot [\delta _c(\mathbf {x},\tau)\mathbf {v}_c(\mathbf {x},\tau)], \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} &&\dot{\delta }_b(\mathbf {x},\tau)+\theta _b(\mathbf {x},\tau) =-\mathbf \nabla \cdot [\delta _b(\mathbf {x},\tau)\mathbf {v}_b(\mathbf {x},\tau)], \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} && \dot{\theta }_c(\mathbf {x},\tau) \!+\!\frac{2}{\tau }\theta _c(\mathbf {x},\tau)\!+\!\frac{6}{\tau ^2}\delta (\mathbf {x},\tau)=-\mathbf {\nabla }\!\cdot \!\lbrace [\mathbf {v}_c(\mathbf {x},\tau)\cdot \mathbf {\nabla }]\mathbf {v}_c(\mathbf {x},\tau)\rbrace,\nonumber\\ \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} &&\fl \dot{\theta }_b(\mathbf {x},\tau) \!+\!\frac{2}{\tau }\theta _b(\mathbf {x},\tau)\!+\!\frac{6}{\tau ^2}\delta (\mathbf {x},\tau)\!\hphantom{aaaaaaaaaaaaaaaaaaaaaaaa} \nonumber \\ &&\quad =\!-\mathbf {\nabla }\cdot \lbrace [\mathbf {v}_b(\mathbf {x},\tau)\cdot \mathbf {\nabla }]\mathbf {v}_b(\mathbf {x},\tau)\rbrace \!-\!\mathbf {\nabla }\!\cdot \bigg[\!\frac{c_s^2(\mathbf {x},\tau)\mathbf {\nabla }\delta _b(\mathbf {x},\tau)}{1+\delta _b(\mathbf {x},\tau)}\!\bigg].\nonumber\\ \end{eqnarray} \tag{ 12 }$$

Note that the second term on the right-hand side of Equation (12) still contains the linear order term. All the other terms on the right-hand side of the above equations are nonlinear.

We shall simplify the pressure term, the second term on the right-hand side of Equation (12), as follows. First, we shall assume that the sound speed is homogeneous, i.e., ∇c²_s = 0. See Naoz & Barkana (2005) for the analysis of linear perturbations with ∇c²_s ≠ 0. Second, we expand the pressure term to the third order in perturbations:

$$\begin{eqnarray} &&\frac{\mathbf {\nabla }\delta \rho _b}{\rho _b} =\frac{\mathbf {\nabla }\delta _b}{1+\delta _b} \simeq \mathbf {\nabla }\delta _b -\delta _b\mathbf {\nabla }\delta _b+\delta _b^2\mathbf {\nabla }\delta _b +\mathcal {O}\big(\delta _b^4\big). \end{eqnarray} \tag{ 13 }$$

Going to Fourier space, we obtain

$$\begin{eqnarray} \dot{\tilde{\delta }}_c(\mathbf {k},\tau)+\tilde{\theta }_c(\mathbf {k},\tau) &=& -\!\frac{1}{(2\pi)^3}\int \!\int d\mathbf {q}_1d\mathbf {q}_2 \delta _D(\mathbf {q}_1+\mathbf {q}_2-\mathbf {k})\nonumber\\ && \times \frac{\mathbf {k}\,{\bm \cdot}\, \mathbf {q}_1}{q_1^2} \tilde{\theta }_c(\mathbf {q}_1,\tau)\tilde{\delta }_c(\mathbf {q}_2,\tau), \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} \dot{\tilde{\delta }}_b(\mathbf {k},\tau)+\tilde{\theta }_b(\mathbf {k},\tau) &=& -\!\frac{1}{(2\pi)^3}\int \!\int d\mathbf {q}_1d\mathbf {q}_2 \delta _D(\mathbf {q}_1+\mathbf {q}_2-\mathbf {k})\nonumber\\ && \times\frac{\mathbf {k}\,{\bm \cdot}\, \mathbf {q}_1}{q_1^2} \tilde{\theta }_b(\mathbf {q}_1,\tau)\tilde{\delta }_b(\mathbf {q}_2,\tau), \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} && \dot{\tilde{\theta }}_c(\mathbf {k},\tau)+ \frac{2}{\tau }\tilde{\theta }_c(\mathbf {k},\tau)+\frac{6}{\tau ^2} \tilde{\delta }(\mathbf {k},\tau) \nonumber\\ &&\quad= -\frac{1}{(2\pi)^3}\int \!\int d\mathbf {q}_1d\mathbf {q}_2 \delta _D(\mathbf {q}_1+\mathbf {q}_2-\mathbf {k}) \nonumber\\ && \,\,\,\,\,\,\quad\times\frac{k^2(\mathbf {q}_1 \,{\bm \cdot}\, \mathbf {q}_2)}{2q_1^2q_2^2} \tilde{\theta }_c(\mathbf {q}_1,\tau)\tilde{\theta }_c(\mathbf {q}_2,\tau), \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} && \dot{\tilde{\theta }}_b(\mathbf {k},\tau) +\frac{2}{\tau }\tilde{\theta }_b(\mathbf {k},\tau) +\frac{6}{\tau ^2}\tilde{\delta }(\mathbf {k},\tau) \nonumber\\ && \quad = -\frac{1}{(2\pi)^3}\int \!\int d\mathbf {q}_1d\mathbf {q}_2 \delta _D(\mathbf {q}_1+\mathbf {q}_2-\mathbf {k})\nonumber\\ &&\,\,\,\,\,\,\quad\times \frac{k^2(\mathbf {q}_1\cdot \mathbf {q}_2)}{2q_1^2q_2^2} \tilde{\theta }_b(\mathbf {q}_1,\tau)\tilde{\theta }_b(\mathbf {q}_2,\tau) \nonumber \\ && \,\,\,\,\,\,\quad-\mathcal {F} \bigg[\mathbf {\nabla }\cdot \bigg(\frac{c_s^2(\tau)\mathbf {\nabla }\delta _b(\mathbf {x},\tau)}{1+\delta _b(\mathbf {x},\tau)}\bigg)\bigg](\mathbf {k}), \end{eqnarray} \tag{ 17 }$$

where

$$\begin{eqnarray} &&\!\!\!\!\!\!\!\!\!\!\!\!\mathcal {F}\bigg[\mathbf {\nabla }\cdot \bigg(\frac{c_s^2(\tau)\mathbf {\nabla }\delta _b(\mathbf {x},\tau)}{1+\delta _b(\mathbf {x},\tau)}\bigg) \bigg](\mathbf {k}) \nonumber\\ && \!\!\!\!\equiv -k^2c_s^2(\tau)\bigg[ \tilde{\delta }_b(\mathbf {k}) -\frac{1}{2(2\pi)^3}\int \!\int d\mathbf {q}_1d\mathbf {q}_2 \tilde{\delta }_b(\mathbf {q}_1,\tau)\tilde{\delta }_b \nonumber\\ && \,\,\times\,(\mathbf {q}_2,\tau) \delta _D(\mathbf {q}_1+\mathbf {q}_2-\mathbf {k}) +\frac{1}{3(2\pi)^6}\int \!\int \!\int d\mathbf {q}_1 d\mathbf {q}_2d\mathbf {q}_3 \tilde{\delta }_b\nonumber\\ && \,\,\times\,(\mathbf {q}_1,\tau)\tilde{\delta }_b (\mathbf {q}_2,\tau) \tilde{\delta }_b(\mathbf {q}_3,\tau) \delta _D(\mathbf {q}_1+\mathbf {q}_2+\mathbf {q}_3-\mathbf {k}) \bigg]. \end{eqnarray} \tag{ 18 }$$

In the subsequent subsections, we shall solve these coupled equations perturbatively. Hereafter, we shall omit the tildes on the perturbation variables in Fourier space.

2.2. Linear Order Solution: Jeans Filtering Scale

In the linear order, one finds

$$\begin{eqnarray} &&\dot{\delta }_{1,c}(\mathbf {k},\tau)+\theta _{1,c}(\mathbf {k},\tau)=0, \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} &&\dot{\delta }_{1,b}(\mathbf {k},\tau)+\theta _{1,b}(\mathbf {k},\tau)=0, \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} &&\dot{\theta }_{1,c}(\mathbf {k},\tau)+\frac{2}{\tau }\theta _{1,c}(\mathbf {k},\tau)+\frac{6}{\tau ^2} \delta _1(\mathbf {k},\tau)=0, \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} &&\dot{\theta }_{1,b}(\mathbf {k},\tau)+\frac{2}{\tau }\theta _{1,b}(\mathbf {k},\tau)+\frac{6}{\tau ^2} \delta _1(\mathbf {k},\tau) \nonumber \\ &&\quad-k^2c_s^2(\tau)\delta _{1,b}(\mathbf {k},\tau)=0, \end{eqnarray} \tag{ 22 }$$

where the subscripts "1" mean that these quantities denote the first-order perturbations, and δ₁ = f_cδ_1,c + f_bδ_1,b. We rewrite Equation (22) as

$$\begin{eqnarray} & &\dot{\theta }_{1,b}(\mathbf {k},\tau)+\frac{2}{\tau }\theta _{1,b}(\mathbf {k},\tau) +\frac{6}{\tau ^2}\bigg[ \delta _1(\mathbf {k},\tau)-\frac{k^2c_s(\tau)^2\tau ^2}{6}\delta _{1,b}(\mathbf {k},\tau) \bigg]\nonumber \\ &&\quad{=}\,0, \nonumber \\ &&\dot{\theta }_{1,b}(\mathbf {k},\tau)\,{+}\,\frac{2}{\tau }\theta _{1,b}(\mathbf {k},\tau)\,{+}\,\frac{6}{\tau ^2} \bigg[\delta _1(\mathbf {k},\tau)\,{-}\,\frac{k^2}{k^2_J}\delta _{1,b}(\mathbf {k},\tau)\bigg]\,{=}\,0,\nonumber\\ \end{eqnarray} \tag{ 23 }$$

where we have used the usual definition of the Jeans wavenumber, k_J:

$$\begin{equation} k_J(\tau)\equiv \frac{\sqrt{6}}{c_s(\tau) \tau }. \end{equation} \tag{ 24 }$$

The Jeans wavenumber divides the solutions for δ_1,b into two classes: the growing solution for k ≪ k_J, and the oscillatory solution for k ≫ k_J, when there is no CDM, i.e., f_b = 1 and δ₁ = δ_1,b. When δ₁ ≠ δ_1,b, the Jeans wavenumber does not provide a dividing scale for the solutions of δ_1,b.

The Jeans wavenumber depends on the temperature of the matter component "b" as k_J ∝ T^−1/2_bτ⁻¹; thus, k_J depends on time in general, k_J = k_J(τ). However, in order to simplify the problem and obtain physical insights into the effects of pressure on the nonlinear growth of structure, we shall assume that k_J is independent of time, which requires that the matter temperature evolves as if the matter were coupled to radiation, T_b ∝ 1/a ∝ 1/τ². This is not a realistic assumption especially in a low-redshift universe where baryons are decoupled from the radiation background and neutrinos are nonrelativistic—in both cases the temperature evolves as T_b ∝ 1/a² ∝ 1/τ⁴ and thus k_J evolves as k_J ∝ τ ∝ a^1/2, for the adiabatic evolution.

We shall solve the above coupled linear equations iteratively: as CDM is always the most dominant source of gravity, the zeroth-order iterative solution may be found by setting δ₁ → δ_1,c (i.e., f_c → 1). We find the solution for the ratio of the density contrasts, which is often called the "Jeans filtering function" (Gnedin & Hui 1998)

$$\begin{equation} g_1(\mathbf {k},\tau)\equiv \frac{\delta _{1,b}(\mathbf {k},\tau)}{\delta _{1,c}(\mathbf {k},\tau)}, \end{equation} \tag{ 25 }$$

which should be a decreasing function of k due to the effect of pressure. At the zeroth order of iteration, the CDM density contrast grows as

$$\begin{equation} \delta _{1,c}^{(0)}(\mathbf {k},\tau)\propto a\propto \tau ^2, \end{equation} \tag{ 26 }$$

and thus the equation for g₁ simplifies to

$$\begin{eqnarray} &&\ddot{g}_1^{(0)}(\mathbf {k},\tau)+\frac{6}{\tau }\dot{g}_1^{(0)}(\mathbf {k},\tau) +\frac{6}{\tau ^2}\bigg(1+\frac{k^2}{k^2_J}\bigg)g_1^{(0)}(\mathbf {k},\tau) =\frac{6}{\tau ^2}.\qquad \end{eqnarray} \tag{ 27 }$$

The solution for g₁(k, τ) must be normalized such that g₁(k, τ) → 1 as k → 0. We find

$$\begin{eqnarray} &&g_1^{(0)}(k,\tau)=\frac{1}{1+\frac{k^2}{k^2_J}}+\mathcal {O}(\tau ^{m(k)}), \end{eqnarray} \tag{ 28 }$$

where

$$\begin{equation} m(k)\equiv -\frac{5}{2}\bigg[1\pm \sqrt{1-\frac{24}{25}\bigg(1+\frac{k^2}{k^2_J}\bigg)}\bigg]. \end{equation} \tag{ 29 }$$

The second term is a decaying mode, whose amplitude is set by the initial condition, e.g., at the epoch when the baryon temperature was raised (by, say, cosmic re-ionization) to the point where the pressure became important, or at the epoch when the neutrinos became nonrelativistic.

Ignoring the decaying mode (although we shall come back to this later), we have the zeroth-order solution:

$$\begin{eqnarray} &&g_1^{(0)}(k)=\frac{1}{1+\frac{k^2}{k^2_J}}. \end{eqnarray} \tag{ 30 }$$

At the first-order iteration, we have the pressure feedback on the growth of CDM. The evolution of δ⁽¹⁾_1,c depends on k, and is given by

$$\begin{eqnarray} &&\delta _{1,c}^{(1)}(k,\tau)\propto \tau ^{n(k)}, \end{eqnarray} \tag{ 31 }$$

where

$$\begin{eqnarray} n(k)&\equiv & \frac{1}{2}\bigg[-\!1\pm 5\sqrt{1-\frac{24}{25}f_b\big(1-g_1^{(0)}(k)\big)}\bigg] \nonumber \\ &\simeq &\bigg\lbrace\!\!\! \begin{array}{@{}rr@{}}2- \frac{6}{5} f_b\big[1-g^{(0)}_1(k)\big]\\ -3+\frac{6}{5} f_b\big[1-g^{(0)}_1(k)\big]\end{array}. \end{eqnarray} \tag{ 32 }$$

The second equality is valid for f_b[1 − g⁽⁰⁾₁(k)] ≪ 1. The growing mode solution is given by

$$\begin{eqnarray} &&n_+(k)=2-\case{6}{5}f_b\big[1-g^{(0)}_1(k)\big]. \end{eqnarray} \tag{ 33 }$$

As g⁽⁰⁾(k) → 1 and 0 for k → 0 and ∞, respectively, the large-scale and small-scale limits of the growing mode solution is (see, e.g., Section 8.3 of Weinberg 2008 for a recent review)

$$\begin{eqnarray} \delta _{1,c+}^{(1)}(k,\tau)&\propto &\tau ^2\propto a, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ k\ll k_J, \end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray} \delta _{1,c+}^{(1)}(k,\tau)&\propto &\tau ^{2-\frac{6}{5}f_b}\propto a^{1-\frac{3}{5}f_b},\ \ \ \ \ \ \ \ k\gg k_J. \end{eqnarray} \tag{ 35 }$$

The growth of δ_1,c on the spatial scales below the Jeans scale is suppressed relative to that of the large-scale modes.

Taking the first-order iteration solution for δ⁽¹⁾_1,c+ into account, the first-order iteration equation for g⁽¹⁾₁ is

$$\begin{eqnarray} &&\ddot{g}_1^{(1)}(k,\tau)+\frac{1}{\tau }\bigg[1+5\sqrt{1-\frac{24}{25}f_b\big(1-g^{(0)}_1(k)\big)}\bigg] \dot{g}_1^{(1)}(k,\tau) \nonumber \\ &&+\frac{6}{\tau ^2}\bigg[1+\frac{k^2}{k^2_J}-f_b\big(2-g_1^{(0)}(k)\big)\bigg]g_1^{(1)}(k,\tau) =\frac{6(1-f_b)}{\tau ^2},\qquad \end{eqnarray} \tag{ 36 }$$

whose growing mode solution (with the normalization that g⁽¹⁾₁ → 1 for k → 0) is

$$\begin{eqnarray} \nonumber g^{(1)}_1(k)&=&\frac{1-f_b}{1+\frac{k^2}{k^2_J}-f_b\big[2-g^{(0)}_1(k)\big]}\\ &=&\frac{1-f_b}{1-f_b+\frac{k^2}{k^2_J}\big(1-\frac{f_b}{1\,+\,k^2/k_J^2}\big)}. \end{eqnarray} \tag{ 37 }$$

This iteration converges quickly for f_b < 0.5, and further iterations are not necessary. The largest difference between g⁽⁰⁾₁(k) and g⁽¹⁾₁(k) occurs as k/k_J → ∞, and is 100% for f_b = 0.5. If the component "b" is identified with baryons, f_b ≃ 1/6, and the difference is reduced to ∼20%. The difference between g⁽¹⁾₁(k) and g⁽²⁾₁(k) occurs at k ∼ k_J, and is ∼4% for f_b = 0.5, and 0.2% for f_b ≃ 1/6. The difference is much smaller for neutrinos, whose f_b are smaller for the modest choices of the neutrino masses (m_nu,i < 1 eV).

To simplify the subsequent analysis, we shall adopt the zeroth-order iterative solution for the filtering function, g⁽⁰⁾₁ = 1/(1 + k²/k²_J), and the first-order iterative solution for the CDM growth factor, Equation (33), as the solution at the first order in perturbations. This solution is sufficiently accurate for our obtaining the physical insights.

Let us comment on the decaying mode that we have ignored in obtaining Equation (30). This decaying mode is an oscillatory function at $k/k_J> 1/(2\sqrt{6})\simeq 0.2$, representing the acoustic oscillation of the pressured component (Nusser 2000). While this term is a decaying mode, it decays slowly, and is not quite negligible even at low redshift. We show the decaying mode at the zeroth-order iterative solution in Figure 1,

$$\begin{equation} \Delta g_1^{(0)}(k,\tau)\equiv g_1^{(0)}(k,\tau)-\frac{1}{1+\frac{k^2}{k_J^2}}, \end{equation} \tag{ 38 }$$

assuming that the pressure became important at z_* = 10. This figure shows that the decaying mode remains important even until z ∼ 0; thus, technically speaking, ignoring the decaying mode results in an inaccurate form of the filtering function. Nevertheless, we shall ignore it and adopt g₁(k) = 1/(1 + k²/k²_J).

Zoom In Zoom Out Reset image size

The exact form of g₁(k, τ) is not so important for our purposes. The main goal of this paper is to study how nonlinearities affect this function. In other words, we are interested in how the higher order filtering functions, g_n(k, τ), are related to the linear one, g₁(k, τ). One may use any forms of g₁(k, τ) for a better accuracy, depending on the problem (baryons or neutrinos).

2.3. Second- and Third-order Solutions

For the higher order (nth order) density perturbations and velocity-divergence fields, we define the Jeans filtering functions such that

$$\begin{eqnarray} &&g_n(\mathbf {k},\tau)\equiv \frac{\delta _{n,b}(\mathbf {k},\tau)}{\delta _{n,c}(\mathbf {k},\tau)}, \end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray} &&h_n(\mathbf {k},\tau)\equiv \frac{\theta _{n,b}(\mathbf {k},\tau)}{\theta _{n,c}(\mathbf {k},\tau)}. \end{eqnarray} \tag{ 40 }$$

Assuming that CDM dominates the gravitational potential, we find the zeroth-order iteration ansatz in an EdS universe:

$$\begin{eqnarray} \delta _b(\mathbf {k},\tau)&=&\sum ^{\infty }_{n=1}a^n(\tau)\delta _{n,c}(\mathbf {k})g_n(\mathbf {k},\tau), \end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray} \theta _b(\mathbf {k},\tau)&=&\sum ^{\infty }_{n=1}\dot{a}(\tau)a^{n-1}(\tau)\theta _{n,c}(\mathbf {k})h_n(\mathbf {k},\tau). \end{eqnarray} \tag{ 42 }$$

Detailed derivations of the nonlinear filtering functions at the second order, g₂(k, τ), and the third order, g₃(k, τ), are given in Appendix B. The second-order solution is

$$\begin{equation} g_2(\mathbf {k},\tau)=\frac{\frac{10}{3} -\frac{7}{3}\big[1-\frac{\delta _{2,c}^{\prime }(\mathbf {k})}{\delta _{2,c}(\mathbf {k})} \big] }{\frac{10}{3}+\frac{k^2}{k_J^2}} +\mathcal {O}(\tau ^{-9/2}), \end{equation} \tag{ 43 }$$

where

$$\begin{eqnarray} \delta _{2,c}(\mathbf {k})&=&\frac{1}{(2\pi)^3} \int d\mathbf {q}F_2^{(s)}(\mathbf {q},\mathbf {k}-\mathbf {q}) \delta _{1,c}(\mathbf {q})\delta _{1,c}(\mathbf {k}-\mathbf {q}),\quad \end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray} \delta ^{\prime }_{2,c}(\mathbf {k})&=&\frac{1}{(2\pi)^3} \int d\mathbf {q}\bigg[F_2^{(s)}(\mathbf {q},\mathbf {k}-\mathbf {q})+\frac{3}{14}\frac{k^2}{k_J^2}\bigg]\nonumber \\ & &\times g_1(\mathbf {q})g_1(\mathbf {k}-\mathbf {q}) \delta _{1,c}(\mathbf {q})\delta _{1,c}(\mathbf {k}-\mathbf {q}), \end{eqnarray} \tag{ 45 }$$

and F^(s)₂ is a mathematical function given by Equation (A28). The third-order solution is

$$\begin{equation} g_3(\mathbf {k},\tau)=\frac{7-6\big[1-\frac{\delta _{3,c}^{\prime }(\mathbf {k})}{\delta _{3,c}(\mathbf {k})}\big]}{7+\frac{k^2}{k^2_J}}+\mathcal {O}(\tau ^{-13/2}), \end{equation} \tag{ 46 }$$

where

$$\begin{eqnarray} \delta _{3,c}(\mathbf {k})&=& \frac{1}{(2\pi)^6}\int \!\int \!\int d\mathbf {q}_1d\mathbf {q}_2d\mathbf {q}_3 \delta _D(\mathbf {k}-\mathbf {q}_1-\mathbf {q}_2-\mathbf {q}_3)\nonumber \\ & &\times F_3^{(s)}(\mathbf {q}_1,\mathbf {q}_2,\mathbf {q}_3) \delta _{1,c}(\mathbf {q}_1)\delta _{1,c}(\mathbf {q}_2) \delta _{1,c}(\mathbf {q}_3), \end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray} \delta ^{\prime }_{3,c}(\mathbf {k})&=& \frac{1}{(2\pi)^6}\int \!\int \!\int d\mathbf {q}_1d\mathbf {q}_2d\mathbf {q}_3 \delta _D(\mathbf {k}-\mathbf {q}_1-\mathbf {q}_2-\mathbf {q}_3)\nonumber \\ & &\times {\mathcal F}_3^{(s)}(\mathbf {q}_1,\mathbf {q}_2,\mathbf {q}_3) \delta _{1,c}(\mathbf {q}_1)\delta _{1,c}(\mathbf {q}_2) \delta _{1,c}(\mathbf {q}_3), \end{eqnarray} \tag{ 48 }$$

and F^(s)₃ and ${\mathcal F}_3^{(s)}$ are mathematical functions given by Equations (A30) and (B27), respectively. One may check that these functions are properly normalized, i.e., g_n → 1 as k → 0, using δ'_2,c → δ_2,c and δ'_3,c → δ_3,c as k → 0.

Ignoring the decaying modes, let us rewrite g₂ and g₃ as

$$\begin{eqnarray} g_2(\mathbf {k})&=&\frac{1 -\frac{7}{10}\big[1-\frac{\delta _{2,c}^{\prime }(\mathbf {k})}{\delta _{2,c}(\mathbf {k})} \big] }{1+\frac{3}{10}\frac{k^2}{k_J^2}}, \end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray} g_3(\mathbf {k})&=& \frac{1-\frac{6}{7}\big[1-\frac{\delta _{3,c}^{\prime }(\mathbf {k})}{\delta _{3,c}(\mathbf {k})}\big] }{1+\frac{1}{7}\frac{k^2}{k^2_J}}. \end{eqnarray} \tag{ 50 }$$

These results may be interpreted as, roughly speaking, the nonlinear filtering functions having smaller effective filtering scales (larger filtering wavenumbers): $k_J\rightarrow \tilde{k}_J=\sqrt{\frac{10}{3}}k_J$ for the second order, $k_J\rightarrow \tilde{k}_J=\sqrt{7}k_J$ for the third order, and $k_J\rightarrow \tilde{k}_J= \sqrt{\frac{2}{3}n(n+\frac{1}{2})}k_J$ for the nth order perturbations. In other words, the higher order solutions for δ_n,b are less suppressed relative to the CDM solutions. In the following section, we shall quantify this effect in more detail.

In this section, we calculate the nonlinear matter power spectrum using the results obtained in the previous section. The total matter fluctuation, δ, is given by δ = f_cδ_c + f_bδ_b, and thus the total matter power spectrum, P_tot(k), is given by the sum of three contributions:

$$\begin{equation} P_{\rm tot}(k,\tau) = f_c^2 P_c(k,\tau) + f_cf_bP_{bc}(k,\tau) + f_b^2P_b(k,\tau), \end{equation} \tag{ 51 }$$

where P_c(k) and P_b(k) are the power spectra of the CDM and another matter component with pressure, respectively, and P_bc(k) is the cross-correlation power spectrum. Each term is the sum of the linear part, P₁₁(k, τ), and the nonlinear parts, P₂₂(k, τ) and P₁₃(k, τ):

$$\begin{equation} P_i(k,\tau) = P_{11,i}(k,\tau) + P_{22,i}(k,\tau) + 2P_{13,i}(k,\tau), \end{equation} \tag{ 52 }$$

where i = (c, b, bc).

The 3PT power spectrum of CDM has been found in the literature (see Bernardeau et al. 2002, for a review)

$$\begin{eqnarray} P_{22,c}(k,\tau) &=& 2\int \frac{d\mathbf {q}}{(2\pi)^3} P_{11,c}(q,\tau)P_{11,c}(|\mathbf {k}-\mathbf {q}|,\tau) \nonumber \\ & &\times \big[F^{(s)}_2(\mathbf {q},\mathbf {k}-\mathbf {q})\big]^2, \end{eqnarray} \tag{ 53 }$$

where F^(s)₂ is a mathematical function given by Equation (A28), and

$$\begin{eqnarray} P_{13,c}(k,\tau)\!\!\!&=& \!\!\!\frac{2\pi }{252}k^2P_{11,c}(k,\tau)\int _0^{\infty }\frac{dq}{(2\pi)^3} P_{11,c}(q,\tau)\nonumber \\ && \!\!\!\times \bigg[ 50\frac{q^2}{k^2} -21\frac{q^4}{k^4}-79+6\frac{k^2}{q^2} \nonumber \\ &&\!\!\!\,+\frac{3}{2} \frac{(q^2-k^2)^3(2k^2+7q^2)}{k^5q^3}\ln {\frac{k+q}{|k-q|}} \bigg]. \end{eqnarray} \tag{ 54 }$$

See Appendix A for the detailed derivations.

Here, we have implicitly generalized the results from an EdS universe to general cosmological models, by writing

$$\begin{eqnarray} &&\frac{a^2(\tau)}{a^2(\tau _i)} P_{11}(k,\tau _i) \rightarrow P_{11}(k,\tau)\! \nonumber \\ &&\quad=\! \frac{D^2(\tau)}{D^2(\tau _i)}\! \Bigg(\frac{\delta _{1,c+}^{(1)}(k,\tau)/\delta _{1,c+}^{(0)}(k,\tau)}{\delta _{1,c+}^{(1)}(k,\tau _*)/\delta _{1,c+}^{(0)}(k,\tau _*)}\Bigg)^2 \!\!P_{11}(k,\tau _i), \end{eqnarray} \tag{ 55 }$$

where τ_i is some arbitrary epoch, τ_* is the epoch where the pressure effect becomes non-negligible (i.e., re-ionization epoch for baryons and the relativistic to nonrelativistic transition epoch for massive neutrinos), and D(τ) is the linear growth factor appropriate to a given cosmological model. This simple generalization has been shown to provide an excellent approximation to the full calculation: see Bernardeau et al. (2002) for models with nonzero curvature and/or a cosmological constant, and Takahashi (2008) for dynamical dark energy models with a constant equation of state of dark energy.

The linear spectra of the other contributions, P_11,bc and P_11,b, are given by

$$\begin{eqnarray} P_{11,bc}(k,\tau) &=& g_1(k)P_{11,c}(k,\tau), \end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray} P_{11,b}(k,\tau) &=& g_1^2(k)P_{11,c}(k,\tau). \end{eqnarray} \tag{ 57 }$$

The nonlinear terms, the main results of this paper, are given by

$$\begin{eqnarray} P_{22,bc}(k,\tau) &=&\frac{1}{\frac{10}{3}+\frac{k^2}{k^2_J}} \bigg[P_{22,c}(k,\tau)\nonumber \\ &&+\frac{14}{3} \int \frac{d\mathbf {q}}{(2\pi)^3} P_{11,c}(q,\tau)P_{11,c}(|\mathbf {k}-\mathbf {q}|,\tau) \nonumber \\ & &\times F^{(s)}_2(\mathbf {q},\mathbf {k}-\mathbf {q}) {\mathcal F}^{(s)}_2(\mathbf {q},\mathbf {k}-\mathbf {q}) \bigg], \end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray} P_{22,b}(k,\tau) &=&\frac{1}{\big(\frac{10}{3}+\frac{k^2}{k^2_J}\big)^2} \bigg[P_{22,c}(k,\tau)\nonumber \\ &&+\frac{28}{3} \int \frac{d\mathbf {q}}{(2\pi)^3} P_{11,c}(q,\tau)P_{11,c}(|\mathbf {k}-\mathbf {q}|,\tau) \nonumber \\ & &\times F^{(s)}_2(\mathbf {q},\mathbf {k}-\mathbf {q}) {\mathcal F}^{(s)}_2(\mathbf {q},\mathbf {k}-\mathbf {q})\nonumber \\ &&+\frac{98}{9} \int \frac{d\mathbf {q}}{(2\pi)^3} P_{11,c}(q,\tau)P_{11,c}(|\mathbf {k}-\mathbf {q}|,\tau) \nonumber \\ & &\times \big({\mathcal F}^{(s)}_2(\mathbf {q},\mathbf {k}-\mathbf {q})\big)^2\bigg], \end{eqnarray} \tag{ 59 }$$

$$\begin{eqnarray} P_{13,bc}(k,\tau) &=&\frac{1}{2}\bigg[\bigg(g_1(k) +\frac{1}{7+\frac{k^2}{k^2_J}}\bigg)P_{13,c}(k,\tau)+\frac{18}{7+\frac{k^2}{k^2_J}}P_{11,c}\nonumber \\ && \times(k,\tau)\int \frac{d\mathbf {q}}{(2\pi)^3} {\mathcal F}^{(s)}_3(\mathbf {q},-\mathbf {q},\mathbf {k}) P_{11,c}(q,\tau)\bigg],\nonumber \\ \end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray} P_{13,b}(k,\tau)\!\!\! &=& \!\!\!\frac{g_1(k)}{7+\frac{k^2}{k^2_J}}\bigg[P_{13,c}(k,\tau)+18P_{11,c}(k,\tau)\nonumber \\ && \times \int \frac{d\mathbf {q}}{(2\pi)^3}{\mathcal F}^{(s)}_3(\mathbf {q},-\mathbf {q},\mathbf {k})P_{11,c}(q,\tau) \bigg]. \end{eqnarray} \tag{ 61 }$$

See Appendix C for the detailed derivations.

How would P_tot(k) compare with the CDM part, P_c(k)?

• 1.  
  In the linear limit, we should recover P_tot(k)/P_c(k) → [f_c + f_bg₁(k)]², which approaches unity as k → 0.
• 2.  
  In the very small scale limit (k → ∞), the pressured component is completely smooth (δ_b(k) → 0) because g₁(k) → 0; thus, P_tot(k)/P_c(k) approaches a constant value, f²_c.
• 3.  
  In the intermediate regime, especially at the transition scale between the super-Jeans scale (k < k_J) and the sub-Jeans scale (k > k_J), the shape of P_tot(k)/P_c(k) is significantly distorted away from the linear prediction. Nonlinear clustering of the pressured component adds power at k ∼ k_J, which shifts the effective filtering scale to smaller spatial scales as we go to lower redshifts.

In Figure 2 we show the ratio, P_tot(k, z)/P_c(k, z) (solid lines in the online version, thick lines in the print version), for different redshifts (z = 0.1, 1, 3, 5, 10, and 30), and different k_J (k_J = 1 and 3 h Mpc⁻¹ for the left and right panels, respectively). In the linear regime (see the bottom lines, z = 30), the ratio agrees with the linear prediction shown by the dashed lines in the online version (thin lines in the print version). As we go to lower redshifts, we find that the filtering wavenumbers continue to shift to larger values, i.e., the filtering scales continue to shift to smaller spatial scales as we go to lower redshifts. This effect cannot be predicted from the linear theory, where all the modes evolve in the same way.

Zoom In Zoom Out Reset image size

The nonlinear power spectrum with a significant contribution from a pressured component has not been studied very much in the literature, with one exception. Saito et al. (2008; hereafter STT) have studied effects of massive neutrinos on the nonlinear matter power spectrum using 3PT (also see Wong 2008; Lesgourgues et al. 2009). However, their treatment is not satisfactory: they have entirely ignored nonlinearities in neutrinos, but approximated the neutrino perturbations as linear perturbations. More precisely, they calculated the nonlinear matter power spectrum as

$$\begin{equation} P^{\rm STT}_{\rm tot}(k,z)=f_c^2P_c(k,z)+2f_cf_{\nu }P_{11,\nu c}(k,z)+f_{\nu }^2P_{11,\nu }(k,z). \end{equation} \tag{ 62 }$$

In our language this leads to

$$\begin{equation} P^{\rm STT}_{\rm tot}(k,z)= f_c^2P_c(k)+\big[2f_cf_{\nu }g_1(k)+f_{\nu }^2g_1^2(k)\big]P_{11,c}(k,z). \end{equation} \tag{ 63 }$$

Here, we have replaced the subscript "b" with "ν" to avoid confusion in notation.

How accurate is the STT approximation? To study this, we compare Equation (63) to the full calculation given in the previous section. Figure 3 shows the fractional difference between our full calculation and STT's approximation, [P_tot(k) − P^STT_tot(k)]/P_tot(k), for Ω_ν/Ω_m = 1/10, 1/20, and 1/100, which correspond to the sum of neutrino masses of ∑_im_ν,i ≃ 1.3, 0.64, and 0.13 eV, respectively, where i = (e, μ, τ). We find that STT's approximation clearly underestimates the power at k ≈ k_FS, where k_FS is the neutrino free-streaming scale, or it is the Jeans wavenumber computed with the velocity dispersion of the neutrinos. More precisely,

$$\begin{equation} k_{{\rm FS},i}(\tau) \equiv \frac{\sqrt{6}}{\sigma _{\nu,i}(\tau)\tau }, \end{equation} \tag{ 64 }$$

in an EdS universe, where σ²_ν,i(τ) is the velocity dispersion of neutrino species i (see, e.g., Appendix A.3 of Takada et al. 2006).

Zoom In Zoom Out Reset image size

One may argue that STT's approximation should be better for a smaller neutrino mass: the errors in the total matter power spectrum are 3.5%, 0.6%, and 0.003% ∑_im_ν,i = 1.3, 0.64, and 0.13 eV, respectively, at z = 0.1; however, our results indicate that their approximation is conceptually not correct: neutrinos should not be treated as linear perturbations, as the neutrino velocity dispersion has no effect in suppressing the neutrino perturbations at and above the free-streaming scale. In other words, the errors may happen to be small in the total matter power spectrum for small neutrino masses because neutrinos contribute only a tiny fraction of the total matter density anyway, but the errors in the neutrino power spectrum are large. Figure 4 shows the fractional difference between the nonlinear neutrino power spectrum, P_ν(k), and the linear power spectrum, P^lin_ν(k), i.e., ΔP/P = [P_ν(k) − P^lin_ν(k)]/P_ν(k). It is clear that neutrinos are significantly nonlinear, even well below the free-streaming scale, k ≫ k_FS. Nevertheless, the STT approximation may still provide a convenient phenomenological tool for calculating the nonlinear total matter power spectrum in the presence of massive neutrinos.

Zoom In Zoom Out Reset image size

In this paper, we have obtained the second- and third-order solutions for the density perturbations in a system consisting of two matter components with and without the pressure gradient. This is the first self-consistent analytical calculation, with nonlinearities in the pressured component fully retained up to the third order in perturbations.

As our study is focused on understanding the physics of the nonlinear pressure effect on the matter power spectrum, we have studied a toy model in which the Jeans wavenumber, k_J, is independent of time. This is equivalent to the temperature of the pressured component following that of radiation, i.e., T ∝ 1/a.

Nevertheless, we have found several results that have qualitative implications for the practical applications. We have found that nonlinearities in the pressured component shift the filtering scale from the well known linear filtering scale (Gnedin & Hui 1998) to a smaller spatial scale (larger wavenumber) by a factor depending on the redshift and the Jeans scale. In other words, the actual filtering scale for a given sound speed (or temperature) is smaller than the linear scale. Therefore, if one used the linear filtering scale to interpret the fall-off of, e.g., the flux power spectrum of the Lyα forests (Zaldarriaga et al. 2001), one would underestimate the temperature of the pressured component.

How important is this effect? For example, when the Jeans wavenumber is k_J = 10 h Mpc^-1, our calculation predicts that the effective filtering wavenumber is ≃10, 12, 13, 13, and 14 h Mpc^-1 at z = 30, 10, 5, 3, and 1, respectively. While we do not expect 3PT to be valid at such high wavenumbers, our results clearly indicate that the expected changes in the filtering scale cannot be ignored. Table 1 summarizes the ratios of the effective (actual) and the linear filtering wavenumbers. Note that the linear filtering wavenumber is the same as the Jeans wavenumber in our model; thus, we show k_F,eff/k_J in Table 1. We extracted the effective filtering wavenumber, k_F,eff, by fitting [f_c + f_b/(1 + k²/k²_F,eff)]² to P_tot(k, z)/P_c(k, z). We find that a factor of 1.4 change in the filtering scale is quite common over a wide range of redshifts and k_J.

A factor of 1.4 change in the filtering scale changes the inferred temperature by a factor of 2; thus, one implication of our result is that the temperature of the intergalactic medium (IGM) obtained from the Lyα forests at z = 3 by Zaldarriaga et al. (2001) might have been underestimated by a factor of 2.

A factor of 1.4 change in the filtering scale gives a factor of ∼3 change in the filtering mass. Our calculation shows that the actual filtering mass is similar to the linear one only in high redshifts, while the former is significantly smaller than the latter in low redshift. This result is qualitatively similar to those found in Okamoto et al. (2008) and Hoeft et al. (2006); however, a quantitative comparison is not possible, as our results apply only to the system with a constant Jeans wavenumber.

What is next? As for baryons, we need to extend our formalism for incorporating a realistic thermal history of the universe with a proper time dependence of k_J. As for neutrinos, we need to incorporate not only the pressure gradient but also the anisotropic stress in the Euler equation. To do this we need to solve the Boltzmann equation. Nevertheless, our results presented in this paper show that neutrinos are significantly nonlinear, even well below the free-streaming scale.

This material is based in part upon work supported by the Texas Advanced Research Program under Grant 003658-0005-2006, by NASA grants NNX08AM29G and NNX08AL43G, and by NSF grant AST-0807649. E.K. acknowledges support from an Alfred P. Sloan Research Fellowship.

The continuity, Euler, and Poisson equations of CDM are given by

• 1.  
  Continuity equation:

  $$\begin{equation} \dot{\delta }(\mathbf {x},\tau)+\mathbf {\nabla }\cdot [(1+\delta (\mathbf {x},\tau))\mathbf {v}(\mathbf {x},\tau)]=0. \end{equation} \tag{ A1 }$$
• 2.  
  Euler equations:

  $$\begin{equation} \dot{\mathbf {v}}(\mathbf {x},\tau)+[\mathbf {v}(\mathbf {x},\tau)\cdot \mathbf {\nabla }]\mathbf {v}(\mathbf {x},\tau) =-\frac{\dot{a}}{a}\mathbf {v}(\mathbf {x},\tau)-\mathbf {\nabla }\phi (\mathbf {x},\tau). \end{equation} \tag{ A2 }$$
• 3.  
  Poisson equation (for an EdS universe):

  $$\begin{equation} \nabla ^2\phi (\mathbf {x},\tau)=\frac{6}{\tau ^2}\delta (\mathbf {x},\tau). \end{equation} \tag{ A3 }$$

First, we take the divergence of Equation (A2) and substitute Equation (A3). Moving all the nonlinear terms to the right-hand side of the equations, we find

$$\begin{eqnarray} &&\dot{\delta }(\mathbf {x},\tau)+\mathbf {\nabla }\cdot \mathbf {v}(\mathbf {x},\tau) =-\mathbf \nabla \cdot [\delta (\mathbf {x},\tau)\mathbf {v}(\mathbf {x},\tau)], \end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray} &\displaystyle\frac{\partial}{\partial \tau }[\mathbf {\nabla }\cdot \mathbf {v}(\mathbf {x},\tau)] +\frac{\dot{a}}{a}[\mathbf {\nabla }\cdot \mathbf {v}(\mathbf {x},\tau)]+\frac{6}{\tau ^2}\delta (\mathbf {x},\tau)& \nonumber \\ &\!\!\!\!\!\!\!=-\mathbf {\nabla }\cdot \lbrace [\mathbf {v}(\mathbf {x},\tau)\cdot \mathbf {\nabla }]\mathbf {v}(\mathbf {x},\tau)\rbrace.& \end{eqnarray} \tag{ A5 }$$

Let us take the Fourier transform of Equations (A4) and (A5)

$$\begin{eqnarray} &&\dot{\tilde{\delta }}(\mathbf {k},\tau)+\tilde{\theta }(\mathbf {k},\tau) \,{=}\,\,{-}\,\frac{1}{(2\pi)^3}\nonumber \\ &&\quad\,{\times}\,\int \!\!\!\int d\mathbf {q}_1d\mathbf {q}_2 \delta _D(\mathbf {q}_1\,{+}\,\mathbf {q}_2\,{-}\,\mathbf {k}) \frac{\mathbf {k}\mathbf \,{\bm \cdot}\, \mathbf {q}_1}{q_1^2} \tilde{\theta }(\mathbf {q}_1,\tau)\tilde{\delta }(\mathbf {q}_2,\tau),\nonumber \\ \end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray} &&\dot{\tilde{\theta }}(\mathbf {k},\tau)+\frac{\dot{a}}{a}\tilde{\theta }(\mathbf {k},\tau)+\frac{6}{\tau ^2} \tilde{\delta }(\mathbf {k},\tau) =-\frac{1}{(2\pi)^3}\nonumber \\ &&\quad\,{\times}\,\int \!\!\!\int d\mathbf {q}_1d\mathbf {q}_2 \delta _D(\mathbf {q}_1\,{+}\,\mathbf {q}_2\,{-}\,\mathbf {k}) \frac{k^2(\mathbf {q}_1\cdot \mathbf {q}_2)}{2q_1^2q_2^2} \tilde{\theta }(\mathbf {q}_1,\tau)\tilde{\theta }(\mathbf {q}_2,\tau),\nonumber \\ \end{eqnarray} \tag{ A7 }$$

where we have defined θ ≡ ∇ · v, and its Fourier transform is given by

$$\begin{equation} \tilde{\mathbf {v}}(\mathbf {k},\tau)=-i\frac{\mathbf {k}}{k^2}\tilde{\theta }(\mathbf {k},\tau). \end{equation} \tag{ A8 }$$

One can decompose the solutions of the nonlinear continuity and Euler equations, $\tilde{\delta }$ and $\tilde{\theta }$, into the sum of infinite series of nth order perturbations of density and velocity divergence fields:

$$\begin{eqnarray} \tilde{\delta }(\mathbf {k},\tau)&=&\sum ^{\infty }_{n=1}a^n(\tau)\delta _n(\mathbf {k}), \end{eqnarray} \tag{ A9 }$$

$$\begin{eqnarray} \tilde{\theta }(\mathbf {k},\tau)&=&\sum ^{\infty }_{n=1}\dot{a}(\tau)a^{n-1}(\tau)\theta _n(\mathbf {k}), \end{eqnarray} \tag{ A10 }$$

respectively. Note that, strictly speaking, this particular decomposition, a decomposition into a series with powers of a(τ), is valid only for an EdS universe. However, generalization to arbitrary cosmological models can be done in the end by replacing a(τ) with the appropriate growth factor, D(τ) (Bernardeau et al. 2002; Takahashi 2008).

Now, let us solve Equations (A6) and (A7) at each order of perturbations. The nth (n > 1) term of Equation (A6) is given by

$$\begin{eqnarray} &&\dot{a}(\tau)a^{n-1}(\tau) [n\delta _n(\mathbf {k})+\theta _n(\mathbf {k})] =-\frac{1}{(2\pi)^3}\int \!\int d\mathbf {q}_1d\mathbf {q}_2 \delta _D\nonumber \\ &&\quad\times(\mathbf {q}_1\,{+}\,\mathbf {q}_2-\mathbf {k}) \frac{\mathbf {k}\cdot \mathbf {q}_1}{q_1^2} \sum ^{n-1}_{m\,{=}\,1}\dot{a}(\tau)a^{n-1}(\tau) \theta _m(\mathbf {q}_1)\delta _{n-m}(\mathbf {q}_2).\nonumber \\ \end{eqnarray} \tag{ A11 }$$

Dividing both sides by $\dot{a}(\tau)a^{n-1}(\tau)$, one obtains

$$\begin{eqnarray} &&n\delta _n(\mathbf {k})+\theta _n(\mathbf {k})=A_n(\mathbf {k}), \end{eqnarray} \tag{ A12 }$$

where

$$\begin{eqnarray} A_n(\mathbf {k})&=&-\frac{1}{(2\pi)^3}\int \!\int d\mathbf {q}_1d\mathbf {q}_2 \delta _D(\mathbf {q}_1+\mathbf {q}_2-\mathbf {k}) \nonumber \\ &&\times\frac{\mathbf {k}\cdot \mathbf {q}_1}{q_1^2} \sum ^{n-1}_{m=1}\theta _m(\mathbf {q}_1)\delta _{n-m}(\mathbf {q}_2). \end{eqnarray} \tag{ A13 }$$

Similarly, from the Euler equation, Equation (A7), one obtains

$$\begin{eqnarray} &&3\delta _n(\mathbf {k})+(1+2n)\theta _n(\mathbf {k})=B_n(\mathbf {k}), \end{eqnarray} \tag{ A14 }$$

where

$$\begin{eqnarray} B_n(\mathbf {k})&=&-\frac{1}{(2\pi)^3}\int \!\int d\mathbf {q}_1d\mathbf {q}_2 \delta _D(\mathbf {q}_1+\mathbf {q}_2-\mathbf {k}) \nonumber \\ &&\times\frac{k^2(\mathbf {q}_1\cdot \mathbf {q}_2)}{q_1^2q_2^2} \sum ^{n-1}_{m=1}\theta _m(\mathbf {q}_1)\theta _{n-m}(\mathbf {q}_2). \end{eqnarray} \tag{ A15 }$$

The forms of Equations (A12) and (A14) indicate that the nth order solutions are written in terms of the sum of first to (n − 1)th order solutions, with δ₁(k) = −θ₁(k). By solving Equations (A12) and (A14) for δ_n and θ_n, one obtains

$$\begin{eqnarray} &&\delta _n(\mathbf {k})=\frac{(1+2n)A_n(\mathbf {k})-B_n(\mathbf {k})}{(2n+3)(n-1)}, \end{eqnarray} \tag{ A16 }$$

$$\begin{eqnarray} &&\theta _n(\mathbf {k})=\frac{-3A_n(\mathbf {k})+nB_n(\mathbf {k})}{(2n+3)(n-1)}, \end{eqnarray} \tag{ A17 }$$

which can be rewritten as

$$\begin{eqnarray} \delta _n(\mathbf {k})&=&\frac{1}{(2\pi)^{3n-3}} \int d\mathbf {q}_1 \ldots d\mathbf {q}_n \delta _D(\mathbf {q}_1+ \cdots +\mathbf {q}_n-\mathbf {k}) \nonumber \\ && \times F_n(\mathbf {q}_1, \ldots,\mathbf {q}_n) \delta _1(\mathbf {q}_1) \ldots \delta _1(\mathbf {q}_n), \end{eqnarray} \tag{ A18 }$$

$$\begin{eqnarray} \theta _n(\mathbf {k})&=&-\frac{1}{(2\pi)^{3n-3}} \int d\mathbf {q}_1 \ldots d\mathbf {q}_n \delta _D(\mathbf {q}_1+ \cdots +\mathbf {q}_n-\mathbf {k}) \nonumber \\ &&\times G_n(\mathbf {q}_1, \ldots,\mathbf {q}_n) \delta _1(\mathbf {q}_1) \ldots \delta _1(\mathbf {q}_n). \end{eqnarray} \tag{ A19 }$$

Here, the newly defined kernels, F_n and G_n, can be found from the following recursion relations:

$$\begin{eqnarray} F_n(\mathbf {q}_1, \ldots,\mathbf {q}_n)&=&\sum ^{n-1}_{m=1} \frac{G_m(\mathbf {q}_1, \ldots,\mathbf {q}_m)}{(2n+3)(n-1)} \nonumber \\ &&\times \bigg[(1+2n)\frac{\mathbf {k}\cdot \mathbf {q}_1}{q_1^2} F_{n-m}(\mathbf {q}_{m+1}, \ldots,\mathbf {q}_n)\nonumber \\ &&\times \frac{k^2(\mathbf {q}_1\cdot \mathbf {q}_2)}{q_1^2q_2^2} G_{n-m}(\mathbf {q}_{m+1}, \ldots,\mathbf {q}_n) \bigg],\nonumber \\ \end{eqnarray} \tag{ A20 }$$

and

$$\begin{eqnarray} G_n(\mathbf {q}_1, \ldots,\mathbf {q}_n)&=&\sum ^{n-1}_{m=1} \frac{G_m(\mathbf {q}_1, \ldots,\mathbf {q}_m)}{(2n+3)(n-1)}\nonumber \\ &&\times \bigg[3\frac{\mathbf {k}\cdot \mathbf {q}_1}{q_1^2} F_{n-m}(\mathbf {q}_{m+1}, \ldots,\mathbf {q}_n) \nonumber \\ &&+n\frac{k^2(\mathbf {q}_1\cdot \mathbf {q}_2)}{q_1^2q_2^2} G_{n-m}(\mathbf {q}_{m+1}, \ldots,\mathbf {q}_n) \bigg],\nonumber \\ \end{eqnarray} \tag{ A21 }$$

with the boundary conditions of F₁ = 1 = G₁. The second-order solutions are

$$\begin{eqnarray} F_2(\mathbf {q}_1,\mathbf {q}_2)&=&\frac{5}{7}\frac{\mathbf {k}\cdot \mathbf {q}_1}{q^2_1} +\frac{k^2(\mathbf {q}_1\cdot \mathbf {q}_2)}{7q^2_1q^2_2}, \end{eqnarray} \tag{ A22 }$$

$$\begin{eqnarray} G_2(\mathbf {q}_1,\mathbf {q}_2)&=&\frac{3}{7}\frac{\mathbf {k}\cdot \mathbf {q}_1}{q^2_1} +\frac{2k^2(\mathbf {q}_1\cdot \mathbf {q}_2)}{7q^2_1q^2_2}, \end{eqnarray} \tag{ A23 }$$

where k = q₁ + q₂. The third-order solutions are

$$\begin{eqnarray} F_3(\mathbf {q}_1,\mathbf {q}_2,\mathbf {q}_3) &=& \frac{1}{18} \bigg[\frac{7\mathbf {k}\cdot \mathbf {q}_1}{q^2_1}F_2(\mathbf {q}_2,\mathbf {q}_3) \nonumber \\ &&+\frac{k^2(\mathbf {q}_1\cdot \mathbf {q}_{23})}{q^2_1q^2_{23}}G_2(\mathbf {q}_2,\mathbf {q}_3) \bigg] +\frac{G_2(\mathbf {q}_1,\mathbf {q}_2)}{18}\nonumber \\ &&\times \bigg[\frac{7\mathbf {k}\cdot \mathbf {q}_{12}}{q^2_{12}} +\frac{k^2(\mathbf {q}_{12}\cdot \mathbf {q}_3)}{q^2_{12}q^2_3} \bigg], \end{eqnarray} \tag{ A24 }$$

where q_ij ≡ q_i + q_j and k = ∑q_i.

It is often convenient to have the symmetrized forms of the above kernels. They are

$$\begin{eqnarray} F^{(s)}_2(\mathbf {q}_1,\mathbf {q}_2)&=&\case{1}{2}[ F_2(\mathbf {q}_1,\mathbf {q}_2)+F_2(\mathbf {q}_2,\mathbf {q}_1)], \end{eqnarray} \tag{ A25 }$$

$$\begin{eqnarray} G^{(s)}_2(\mathbf {q}_1,\mathbf {q}_2)&=&\case{1}{2}[ G_2(\mathbf {q}_1,\mathbf {q}_2)+G_2(\mathbf {q}_2,\mathbf {q}_1)], \end{eqnarray} \tag{ A26 }$$

$$\begin{eqnarray} F^{(s)}_3(\mathbf {q}_1,\mathbf {q}_2,\mathbf {q}_3)&=&\case{1}{6}[ F_3(\mathbf {q}_1,\mathbf {q}_2,\mathbf {q}_3)+F_3(\mathbf {q}_1,\mathbf {q}_3,\mathbf {q}_2) \nonumber \\ & &+F_3(\mathbf {q}_2,\mathbf {q}_1,\mathbf {q}_3) +F_3(\mathbf {q}_2,\mathbf {q}_3,\mathbf {q}_1) \nonumber \\ & &+F_3(\mathbf {q}_3,\mathbf {q}_1,\mathbf {q}_2) +F_3(\mathbf {q}_3,\mathbf {q}_2,\mathbf {q}_1) ].\qquad \end{eqnarray} \tag{ A27 }$$

The explicit forms are

$$\begin{eqnarray} F^{(s)}_2(\mathbf {q}_1,\mathbf {q}_2)&=& \frac{5}{7}+\frac{2}{7}\frac{(\mathbf {q}_1\cdot \mathbf {q}_2)^2}{q_1^2q_2^2} +\frac{1}{2}\frac{(\mathbf {q}_1\cdot \mathbf {q}_2)\big(q_1^2+q_2^2\big)}{q_1^2q_2^2},\qquad \end{eqnarray} \tag{ A28 }$$

$$\begin{eqnarray} G^{(s)}_2(\mathbf {q}_1,\mathbf {q}_2)&=& \frac{3}{7}+\frac{4}{7}\frac{(\mathbf {q}_1\cdot \mathbf {q}_2)^2}{q_1^2q_2^2} +\frac{1}{2}\frac{(\mathbf {q}_1\cdot \mathbf {q}_2)\big(q_1^2+q_2^2\big)}{q_1^2q_2^2},\qquad \end{eqnarray} \tag{ A29 }$$

$$\begin{eqnarray} &&F^{(s)}_3(\mathbf {q}_1,\mathbf {q}_2,\mathbf {q}_3)= \frac{7}{54}\mathbf {k}\cdot \bigg[F^{(s)}_2(\mathbf {q}_2,\mathbf {q}_3) \frac{\mathbf {q}_1}{q_1^2} +F^{(s)}_2(\mathbf {q}_1,\mathbf {q}_3)\frac{\mathbf {q}_2}{q_2^2} \nonumber \\ &&\quad+F^{(s)}_2(\mathbf {q}_1,\mathbf {q}_2)\frac{\mathbf {q}_3}{q_3^2}\bigg] +\frac{1}{27}k^2\bigg[ G^{(s)}_2(\mathbf {q}_2,\mathbf {q}_3)\frac{\mathbf {q}_1\cdot \mathbf {q}_{23}}{q_1^2q_{23}^2}\nonumber \\ &&\quad+G^{(s)}_2(\mathbf {q}_1,\mathbf {q}_3)\frac{\mathbf {q}_2\cdot \mathbf {q}_{13}}{q_2^2q_{13}^2} +G^{(s)}_2(\mathbf {q}_1,\mathbf {q}_2)\frac{\mathbf {q}_3\cdot \mathbf {q}_{12}}{q_3^2q_{12}^2} \bigg] +\frac{7}{54}\mathbf {k}\cdot \nonumber \\ &&\quad\times\bigg[G^{(s)}_2(\mathbf {q}_2,\mathbf {q}_3) \frac{\mathbf {q}_{23}}{q_{23}^2} +G^{(s)}_2(\mathbf {q}_1,\mathbf {q}_3)\frac{\mathbf {q}_{13}}{q_{13}^2} +G^{(s)}_2(\mathbf {q}_1,\mathbf {q}_2)\frac{\mathbf {q}_{12}}{q_{12}^2}\bigg].\nonumber \\ \end{eqnarray} \tag{ A30 }$$

In order to calculate the next-to-linear-order density power spectrum, one needs to use the solutions of the density fluctuations up to the third order:

$$\begin{eqnarray} &&(2\pi)^3P(k,\tau)\delta _D(\mathbf {k}+\mathbf {k}^{\prime })= \langle \tilde{\delta }(\mathbf {k},\tau)\tilde{\delta }(\mathbf {k}^{\prime },\tau)\rangle \nonumber \\ &&\quad\,{=}\,\bigg<\bigg(\sum ^{\infty }_{m\,{=}\,1}a^m(\tau)\tilde{\delta }_m(\mathbf {k})\bigg) \bigg(\sum ^{\infty }_{l\,{=}\,1}a^l(\tau)\tilde{\delta }_l(\mathbf {k}^{\prime })\bigg)\bigg>\,{\simeq}\, a^2(\tau)\nonumber \\ &&\quad\quad \,{\times}\,\langle \delta _1(\mathbf {k})\delta _1(\mathbf {k}^{\prime })\rangle +a^4(\tau)\langle \delta _1(\mathbf {k})\delta _3(\mathbf {k}^{\prime }) \nonumber \\ &&\quad\quad \,{+}\,\delta _2(\mathbf {k})\delta _2(\mathbf {k}^{\prime }) +\delta _3(\mathbf {k})\delta _1(\mathbf {k}^{\prime })\rangle, \end{eqnarray} \tag{ A31 }$$

which yields

$$\begin{equation} P(k,\tau)= a^2(\tau)P_{11}(k)+a^4(\tau)[P_{22}(k)+2P_{13}(k)]+\mathcal {O}(\delta ^6). \end{equation} \tag{ A32 }$$

Here, we have defined the quantity, P_ij(k), given by

$$\begin{eqnarray} &&(2\pi)^3P_{ij}(k)\delta _D(\mathbf {k}+\mathbf {k}^{\prime }) =\langle \delta _i(\mathbf {k})\delta _j(\mathbf {k}^{\prime })\rangle. \end{eqnarray} \tag{ A33 }$$

The nonlinear corrections, P₂₂(k) and P₁₃(k), are

$$\begin{eqnarray} &&P_{22}(k)= 2\int \frac{d\mathbf {q}}{(2\pi)^3} P_{11}(q)P_{11}(|\mathbf {k}-\mathbf {q}|) \big[F^{(s)}_2(\mathbf {q},\mathbf {k}-\mathbf {q})\big]^2,\nonumber \\ \end{eqnarray} \tag{ A34 }$$

where

$$\begin{eqnarray} &&F^{(s)}_2(\mathbf {q},\mathbf {k}-\mathbf {q}) =\frac{5}{7}\nonumber \\ &&\quad+\frac{1}{14}\bigg[ \frac{-10q^4+20kq^3\mu -10k^2q^2\mu ^2-7k^2q^2+7k^3q\mu }{q^2(k^2+q^2-2kq\mu)}\bigg],\nonumber\\ \end{eqnarray} \tag{ A35 }$$

and $\mu \equiv \hat{\mathbf {k}}\cdot \hat{\mathbf {q}}$, and

$$\begin{eqnarray} &&P_{13}(k) = 3P_{11}(k) \int \frac{d\mathbf {q}}{(2\pi)^3}F^{(s)}_3(\mathbf {q},-\mathbf {q},\mathbf {k}) P_{11}(q). \end{eqnarray} \tag{ A36 }$$

Using

$$\begin{eqnarray*} \int _{-1}^1d\mu F^{(s)}_3(\mathbf {q},-\mathbf {q},\mathbf {k})&=& \frac{1}{756}\bigg[ 50-21\frac{q^2}{k^2}-79\frac{k^2}{q^2}+6\frac{k^4}{q^4}+\frac{3}{2} \nonumber\\ &&\,{\times}\,\frac{(q^2\,{-}\,k^2)^3(2k^2\,{+}\,7q^2)}{k^3q^5}\ln {\frac{k+q}{|k\,{-}\,q|}} \bigg], \nonumber\\ \end{eqnarray*}$$

one obtains (Makino et al. 1992)

$$\begin{eqnarray} P_{13}(k)&=& \frac{2\pi }{252}k^2P_{11}(k)\int _0^{\infty }\frac{dq}{(2\pi)^3}P_{11}(q) \bigg[ 50\frac{q^2}{k^2}-21\frac{q^4}{k^4}-79\nonumber\\ &&+6\frac{k^2}{q^2}+\frac{3}{2} \frac{(q^2-k^2)^3(2k^2+7q^2)}{k^5q^3}\ln {\frac{k+q}{|k-q|}} \bigg]. \end{eqnarray} \tag{ A37 }$$

In this appendix, we shall derive the higher order filtering functions. We shall solve Equations (14)–(17) perturbatively, up to the third order in perturbations. The density contrasts and velocity divergence fields of CDM and the matter with pressure are all expanded into the infinite sum of nth order perturbations as

$$\begin{eqnarray} &&\tilde{\delta }_c(\mathbf {k},\tau)=\sum ^{\infty }_{n=1}a^n(\tau)\delta _{n,c}(\mathbf {k}), \end{eqnarray} \tag{ B1 }$$

$$\begin{eqnarray} &&\tilde{\theta }_c(\mathbf {k},\tau)=\sum ^{\infty }_{n=1}\dot{a}(\tau)a^{n-1}(\tau)\theta _{n,c}(\mathbf {k}), \end{eqnarray} \tag{ B2 }$$

$$\begin{eqnarray} \tilde{\delta }_b(\mathbf {k},\tau)&=&\sum ^{\infty }_{n=1}a^n(\tau)\delta _{n,c}(\mathbf {k})g_n(\mathbf {k},\tau), \end{eqnarray} \tag{ B3 }$$

$$\begin{eqnarray} \tilde{\theta }_b(\mathbf {k},\tau)&=&\sum ^{\infty }_{n=1}\dot{a}(\tau)a^{n-1}(\tau)\theta _{n,c}(\mathbf {k})h_n(\mathbf {k},\tau), \end{eqnarray} \tag{ B4 }$$

where g_n(k, τ) and h_n(k, τ) are the filtering functions for the density and velocity divergence fields, respectively, at the nth order.

With the above series expansion, Equations (15) and (17) yield

$$\begin{eqnarray} &&\sum _{n=1}^{\infty }[ (n\dot{a}(\tau)a^{n-1}(\tau)g_n(\mathbf {k},\tau)+a^n(\tau)\dot{g}_n(\mathbf {k},\tau))\nonumber\\ &&\quad\quad\times\delta _{n,c}(\mathbf {k}) +\dot{a}(\tau)a^{n-1}(\tau)h_n(\mathbf {k},\tau)\theta _{n,c}(\mathbf {k})] \nonumber \\ &&\quad=-\frac{1}{(2\pi)^3}\int \!\int d\mathbf {q}_1d\mathbf {q}_2 \delta _D(\mathbf {q}_1+\mathbf {q}_2-\mathbf {k}) \frac{\mathbf {k}\cdot \mathbf {q}_1}{q_1^2}\sum _{m=1}^{\infty }\nonumber\\ &&\quad\quad\times \sum _{l=1}^{\infty }\dot{a}a^{m+l-1} h_m(\mathbf {q}_1,\tau)g_l(\mathbf {q}_2,\tau) \theta _{m,c}(\mathbf {q}_1)\delta _{l,c}(\mathbf {q}_2),\qquad \end{eqnarray} \tag{ B5 }$$

$$\begin{eqnarray} &&\sum _{n=1}^{\infty }[ (\ddot{a}(\tau)a^{n-1}(\tau)+\dot{a}^2(\tau)a^{n-2}(\tau)(n-1))h_n(\mathbf {k},\tau)\theta _{n,c}(\mathbf {k}) \nonumber \\ &&\quad+\dot{a}(\tau)a^{n-1}(\tau)\dot{h}_n(\mathbf {k},\tau)\theta _{n,c}(\mathbf {k}) +\frac{2}{\tau }\dot{a}(\tau)a^{n-1}(\tau)h_n(\mathbf {k})(\mathbf {k},\tau)\theta _{n,c} \nonumber \\ &&\quad\times(\mathbf {k})+\frac{6}{\tau ^2}a^n(\tau)(f_c+f_bg_n(\mathbf {k},\tau))\delta _{n,c}(\mathbf {k}) ] =-\frac{1}{(2\pi)^3}\nonumber \\ &&\quad\times\int \!\int d\mathbf {q}_1d\mathbf {q}_2 \delta _D(\mathbf {q}_1+\mathbf {q}_2-\mathbf {k}) \frac{k^2(\mathbf {q}_1\cdot \mathbf {q}_2)}{2q_1^2q_2^2} \sum _{m=1}^{\infty }\sum _{l=1}^{\infty }\dot{a}^2(\tau)\nonumber \\ &&\quad\times a^{m+l-2}(\tau) h_m(\mathbf {q}_1,\tau)h_l(\mathbf {q}_2,\tau) \theta _{m,c}(\mathbf {q}_1)\theta _{l,c}(\mathbf {q}_2) +k^2c_s^2(\tau)\nonumber \\ &&\quad\,{\times}\,\sum _{n\,{=}\,1}^{\infty }a^n(\tau)g_n(\mathbf {k},\tau)\delta _{n,c}(\mathbf {k})\,{-}\,\frac{1}{2(2\pi)^3}k^2c_s^2(\tau)\int \!\!\int\! d\mathbf {q}_1d\mathbf {q}_2 \delta _D\nonumber \\ &&\quad\times(\mathbf {q}_1+\mathbf {q}_2-\mathbf {k}) \sum _{m=1}^{\infty }\sum _{l=1}^{\infty }a^{m+l}(\tau) g_m(\mathbf {q}_1,\tau)g_l(\mathbf {q}_2,\tau) \delta _{m,c}\nonumber \\ &&\quad\times(\mathbf {q}_1)\delta _{l,c}(\mathbf {q}_2) +\frac{1}{3(2\pi)^6}k^2c_s^2(\tau)\int \!\int \!\int d\mathbf {q}_1 d\mathbf {q}_2d\mathbf {q}_3 \delta _D\nonumber \\ &&\quad\times(\mathbf {q}_1+\mathbf {q}_2+\mathbf {q}_3-\mathbf {k}) \sum _{m=1}^{\infty }\sum _{l=1}^{\infty }\sum _{p=1}^{\infty }a^{m+l+p}(\tau) g_m(\mathbf {q}_1,\tau)g_l\nonumber \\ &&\quad\times(\mathbf {q}_2,\tau) g_p(\mathbf {q}_3,\tau) \delta _{m,c}(\mathbf {q}_1)\delta _{l,c}(\mathbf {q}_2) \delta _{p,c}(\mathbf {q}_3). \end{eqnarray} \tag{ B6 }$$

From now on, we shall write the sound speed, c_s, in terms of the usual Jeans wavenumber, k_J, as $c_s=\sqrt{6}/(k_J\tau)$. We shall ignore the inhomogeneity in c_s (i.e., spatial dependence of c_s) throughout this paper. For the linear analysis for ∇c_s ≠ 0, see Naoz & Barkana (2005).

B.1. Second-order Solutions

We have derived the linear filtering function, g₁(k), in Equation (30). For n = 2, the continuity and Euler equations are given by

$$\begin{eqnarray} &&\delta _{2,c}(\mathbf {k})\dot{g}_2(\mathbf {k},\tau)+\frac{4}{\tau }\delta _{2,c}(\mathbf {k})g_2(\mathbf {k},\tau) +\frac{2}{\tau }\theta _{2,c}(\mathbf {k})h_2(\mathbf {k},\tau) \nonumber \\ &&=\frac{2}{\tau }\frac{1}{(2\pi)^3} \int \!\int d\mathbf {q}_1d\mathbf {q}_2 \delta _D(\mathbf {q}_1+\mathbf {q}_2-\mathbf {k}) \nonumber \\ &&\quad\times\frac{\mathbf {k}\cdot \mathbf {q}_1}{q_1^2} \delta _{1,c}(\mathbf {q}_1)\delta _{1,c}(\mathbf {q}_2) g_1(\mathbf {q}_1)g_1(\mathbf {q}_2) \equiv \frac{2}{\tau }A_2(\mathbf {k}), \end{eqnarray} \tag{ B7 }$$

$$\begin{eqnarray} &&\frac{10}{\tau ^2}\theta _{2,c}(\mathbf {k})h_2(\mathbf {k},\tau) +\frac{2}{\tau }\theta _{2,c}(\mathbf {k})\dot{h}_2(\mathbf {k},\tau)+\frac{6}{\tau ^2}\delta _{2,c}(\mathbf {k}) \nonumber \\ &&\quad-\frac{6}{\tau ^2}\frac{k^2}{k_J^2}\delta _{2,c}(\mathbf {k})g_2(\mathbf {k},\tau) =\frac{4}{\tau ^2}\frac{1}{(2\pi)^3} \int \!\int d\mathbf {q}_1d\mathbf {q}_2 \delta _D\nonumber \\ &&\quad\times(\mathbf {q}_1+\mathbf {q}_2-\mathbf {k}) \bigg[ -\frac{3}{4}\frac{k^2}{k_J^2} -\frac{k^2(\mathbf {q}_1\cdot \mathbf {q}_2)}{2q_1^2q_2^2} \bigg]\nonumber \\ &&\quad\times \delta _{1,c}(\mathbf {q}_1)\delta _{1,c}(\mathbf {q}_2)g_1(\mathbf {q}_1)g_1(\mathbf {q}_2) \equiv \frac{4}{\tau ^2}B_2(\mathbf {k}). \end{eqnarray} \tag{ B8 }$$

Here, θ_1,c(k) = −δ_1,c(k). Combining Equations (B7) and (B8), we get the second-order inhomogeneous partial differential equation:

$$\begin{eqnarray} \ddot{g}_2(\mathbf {k},\tau)&&+\frac{10}{\tau ^2}\dot{g}_2(\mathbf {k},\tau) +\frac{1}{\tau ^2}\bigg[20+6\frac{k^2}{k_J^2}\bigg]g_2(\mathbf {k},\tau) \nonumber \\ &&+\frac{1}{\tau ^2}\bigg[-6-\frac{10A_2(\mathbf {k})}{\delta _{2,c}(\mathbf {k})} +\frac{4B_2(\mathbf {k})}{\delta _{2,c}(\mathbf {k})}\bigg]=0, \end{eqnarray} \tag{ B9 }$$

where δ_2,c(k) is given by

$$\begin{equation} \delta _{2,c}(\mathbf {k})=\frac{1}{(2\pi)^3} \int d\mathbf {q}F_2^{(s)}(\mathbf {q},\mathbf {k}-\mathbf {q}) \delta _{1,c}(\mathbf {q})\delta _{1,c}(\mathbf {k}-\mathbf {q}). \end{equation} \tag{ B10 }$$

Solving the above differential equation, we have

$$\begin{equation} g_2(\mathbf {k},\tau)=\frac{6 +\frac{10A_2(\mathbf {k})}{\delta _{2,c}(\mathbf {k})} -\frac{4B_2(\mathbf {k})}{\delta _{2,c}(\mathbf {k})}}{20+6\frac{k^2}{k_J^2}} +\mathcal {O}(\tau ^{-9/2}), \end{equation} \tag{ B11 }$$

where the oscillation component,

$$\begin{equation} \mathcal {O}(\tau ^{-9/2})\propto \tau ^{-\frac{9}{2}\bigg(1\pm \sqrt{1-\frac{4}{81}\left(20+6\frac{k^2}{k_J^2}\right)} \bigg)}, \end{equation} \tag{ B12 }$$

decays for any choice of 0 ⩽ k/k_J. The second-order filtering function for the velocity divergence field, h₂(h, τ), is given by

$$\begin{equation} h_2(\mathbf {k})=\case{1}{\theta _{2,c}(\mathbf {k})} [A_2(\mathbf {k})-2\delta _{2,c}(\mathbf {k})g_2(\mathbf {k})], \end{equation} \tag{ B13 }$$

where we have ignored the decaying term.

Using the explicit forms of A₂(k) and B₂(k) given by Equations (B6) and (B8), respectively, we obtain

$$\begin{equation} g_2(\mathbf {k},\tau)=\frac{\frac{10}{3} -\frac{7}{3}\big[1-\frac{\delta _{2,c}^{\prime }(\mathbf {k})}{\delta _{2,c}(\mathbf {k})} \big] }{\frac{10}{3}+\frac{k^2}{k_J^2}} +\mathcal {O}(\tau ^{-9/2}), \end{equation} \tag{ B14 }$$

where δ'_2,c is

$$\begin{eqnarray} &&\delta ^{\prime }_{2,c}(\mathbf {k})= \frac{1}{(2\pi)^3} \int d\mathbf {q} \mathcal {F}_2^{(s)}(\mathbf {q},\mathbf {k}-\mathbf {q}) \delta _{1,c}(\mathbf {q})\delta _{1,c}(\mathbf {k}-\mathbf {q}),\qquad \end{eqnarray} \tag{ B15 }$$

where

$$\begin{equation} \mathcal {F}_2^{(s)}(\mathbf {q}_1,\mathbf {q}_2)\equiv \bigg[ F_2^{(s)}(\mathbf {q}_1,\mathbf {q}_2) +\frac{3}{14}\frac{k^2}{k^2_J}\bigg] g_1(\mathbf {q}_1)g_1(\mathbf {q}_2). \end{equation} \tag{ B16 }$$

In the limit where k_J → ∞, $\mathcal {F}_2^{(s)}(\mathbf {q}_1,\mathbf {q}_2)=F_2^{(s)}(\mathbf {q}_1,\mathbf {q}_2)$, and thus g₂ → 1. For the velocity divergence filtering function, we find

$$\begin{eqnarray} h_2(\mathbf {k})\!&=&\!\frac{1}{\theta _{2,c}(\mathbf {k})}\! \bigg[ \frac{1}{(2\pi)^3}\!\int \!\!\!\int \! d\mathbf {q}_1d\mathbf {q}_2 \delta _D(\mathbf {q}_1+\mathbf {q}_2-\mathbf {k}) \delta _{1,c}(\mathbf {q}_1)\delta _{1,c}\nonumber \\ &&\times(\mathbf {q}_2) (2F_2^{(s)}\!(\mathbf {q}_1,\!\mathbf {q}_2)\!-G_2^{(s)}\!(\mathbf {q}_1,\!\mathbf {q}_2)) g_1(\mathbf {q}_1)g_1(\mathbf {q}_2)\nonumber \\ &&\!-2\delta _{2,c}(\mathbf {k})g_2(\mathbf {k})\bigg]= \!\frac{1}{\theta _{2,c}(\mathbf {k})}\! \bigg[ \frac{1}{(2\pi)^3}\!\int \!\!\!\int \! d\mathbf {q}_1d\mathbf {q}_2 \delta _D\nonumber \\ &&\times(\mathbf {q}_1+\mathbf {q}_2-\mathbf {k}) \delta _{1,c}(\mathbf {q}_1)\delta _{1,c}(\mathbf {q}_2)\! \bigg(\!1\!+\!\frac{(\mathbf {q}_1\cdot \mathbf {q}_2)(q_1^2+q_2^2)}{2q_1^2q_2^2} \bigg)\nonumber \\ &&\times g_1(\mathbf {q}_1)g_1(\mathbf {q}_2)\bigg] \!-2\frac{\delta _{2,c}(\mathbf {k})}{\theta _{2,c}(\mathbf {k})}g_2(\mathbf {k}), \end{eqnarray} \tag{ B17 }$$

where we have used $2F_2(\mathbf {q}_1,\mathbf {q}_2)-G_2(\mathbf {q}_1,\mathbf {q}_2)=\frac{\mathbf {k}\,\cdot\, \mathbf {q}_1}{q_1^2}$. This expression also converges to h₂ = 1 as we take the limit of k_J → ∞.

B.2. Third-order Solutions

For n = 3, the continuity and Euler equations are given by

$$\begin{eqnarray} &&3\dot{a}(\tau)a^2(\tau)g_3(\mathbf {k},\tau)\delta _{3,c}(\mathbf {k}) +a^3(\tau)\dot{g}_3\nonumber \\ &&\quad\times(\mathbf {k},\tau)\delta _{3,c}(\mathbf {k}) +\dot{a}(\tau)a^2(\tau)h_3(\mathbf {k},\tau)\theta _{3,c}(\mathbf {k}) \nonumber \\ &&=\dot{a}(\tau)a^2(\tau)\frac{1}{(2\pi)^6}\int \!\int \!\int d\mathbf {q}_1d\mathbf {q}_2d\mathbf {q}_3 \delta _D\nonumber \\ &&\quad\times(\mathbf {q}_1+\mathbf {q}_2+\mathbf {q}_3-\mathbf {k}) \delta _{1,c}(\mathbf {q}_1)\delta _{1,c}(\mathbf {q}_2) \delta _{1,c}(\mathbf {q}_3) \nonumber \\ &&\quad\times \bigg[ \frac{\mathbf {k}\cdot \mathbf {q}_1}{q_1^2} g_1(\mathbf {q}_1)g_2(\mathbf {q}_{23}) F_2^{(s)}(\mathbf {q}_2,\mathbf {q}_3) \nonumber \\ && \,\quad+\frac{\mathbf {k}\cdot \mathbf {q}_{12}}{q_{12}^2} h_2(\mathbf {q}_{12})g_1(\mathbf {q}_3) G_2^{(s)}(\mathbf {q}_1,\mathbf {q}_2) \bigg] \equiv \dot{a}(\tau)a^2(\tau)A_3(\mathbf {k}),\nonumber \\ \end{eqnarray} \tag{ B18 }$$

$$\begin{eqnarray} &&[\ddot{a}(\tau)a^2(\tau)+2\dot{a}^2(\tau)a(\tau)] h_3(\mathbf {k},\tau)\theta _{3,c}(\mathbf {k}) +\dot{a}(\tau)a^2(\tau)\nonumber \\ &&\quad\times \dot{h}_3(\mathbf {k},\tau)\theta _{3,c}(\mathbf {k}) +\frac{2}{\tau }\dot{a}(\tau)a^2(\tau)h_3(\mathbf {k},\tau)\theta _{3,c}(\mathbf {k}) \nonumber \\ &&\quad+\frac{6}{\tau ^2}a^3(\tau)\delta _{3,c}(\mathbf {k}) -\frac{6}{\tau ^2}\frac{k^2}{k^2_J}a^3(\tau)\delta _{3,c}(\mathbf {k}) \nonumber \\ &&=\dot{a}^2(\tau)a(\tau)\frac{1}{(2\pi)^6}\int \!\!\int \!\!\int\! d\mathbf {q}_1 d\mathbf {q}_2d\mathbf {q}_3 \delta _D\nonumber\\ && \quad\times(\mathbf {q}_1+\mathbf {q}_2+\mathbf {q}_3-\mathbf {k}) \delta _{1,c} (\mathbf {q}_1)\delta _{1,c}(\mathbf {q}_2)\delta _{1,c}(\mathbf {q}_3) \nonumber \\ &&\quad\times \bigg[ -\frac{k^2(\mathbf {q}_1\cdot \mathbf {q}_{23})}{2q_1^2q_{23}^2} g_1(\mathbf {q}_1)h_2(\mathbf {q}_{23}) G_2^{(s)}(\mathbf {q}_2,\mathbf {q}_3) \nonumber \\ &&\quad-\frac{k^2(\mathbf {q}_{12}\cdot \mathbf {q}_3)}{2q_{12}^2q_3^2} h_2 (\mathbf {q}_{12})g_1(\mathbf {q}_3) G_2^{(s)}(\mathbf {q}_1,\mathbf {q}_2) \nonumber \\ &&\quad-\frac{3}{4}\frac{k^2}{k_J^2} g_1(\mathbf {q}_1)g_2(\mathbf {q}_{23}) F_2^{(s)}(\mathbf {q}_2,\mathbf {q}_3) -\frac{3}{4}\frac{k^2}{k_J^2} g_2(\mathbf {q}_{12})g_1(\mathbf {q}_3) F_2^{(s)}\nonumber \\ &&\quad\times(\mathbf {q}_1,\mathbf {q}_2) +\frac{1}{2}\frac{k^2}{k_J^2} g_1(\mathbf {q}_1)g_1(\mathbf {q}_2)g_1(\mathbf {q}_3) \bigg] \nonumber \\ &&\equiv \dot{a}^2(\tau)a(\tau)B_3(\mathbf {k}). \end{eqnarray} \tag{ B19 }$$

In an EdS universe, $a(\tau)=\frac{\tau ^2}{9}$, we have

$$\begin{eqnarray} && \delta _{3,c}(\mathbf {k})\dot{g}_3(\mathbf {k},\tau)+\frac{6}{\tau }\delta _{3,c}(\mathbf {k})g_3(\mathbf {k},\tau) +\frac{2}{\tau }\theta _{3,c}(\mathbf {k})h_3(\mathbf {k},\tau)\nonumber \\ &&\quad=\frac{2}{\tau }A_3(\mathbf {k}), \end{eqnarray} \tag{ B20 }$$

$$\begin{eqnarray} &&\frac{14}{\tau ^2}h_3(\mathbf {k},\tau)\theta _{3,c}(\mathbf {k})+\frac{2}{\tau }\dot{h}_3(\mathbf {k},\tau)\theta _{3,c}(\mathbf {k}) +\frac{6}{\tau ^2}\delta _{3,c}(\mathbf {k})\nonumber \\ &&\quad-\frac{6}{\tau ^2}\frac{k^2}{k_J^2}\delta _{3,c}(\mathbf {k})g_3(\mathbf {k},\tau) =\frac{4}{\tau ^2}B_3(\mathbf {k}). \end{eqnarray} \tag{ B21 }$$

Combining Equations (B20) and (B21), we have the second-order differential equation:

$$\begin{eqnarray} &&\ddot{g}_3(\mathbf {k},\tau)+\frac{14}{\tau }\dot{g}_3(\mathbf {k},\tau)+\frac{1}{\tau ^2} \bigg(42+6\frac{k^2}{k^2_J}\bigg)g_3(\mathbf {k},\tau)\nonumber \\ &&\quad+\frac{1}{\tau ^2}\bigg(-6-\frac{14A_3(\mathbf {k})}{\delta _{3,c}(\mathbf {k})}+\frac{4B_3(\mathbf {k})}{\delta _{3,c}(\mathbf {k})} \bigg)=0. \end{eqnarray} \tag{ B22 }$$

Solving this, we obtain

$$\begin{equation} g_3(\mathbf {k},\tau)= \frac{1+\frac{7A_3(\mathbf {k})}{3\delta _{3,c}(\mathbf {k})}-\frac{2B_3(\mathbf {k})}{3\delta _{3,c}(\mathbf {k})}}{7+\frac{k^2}{k^2_J}} +\mathcal {O}(\tau ^{-13/2}), \end{equation} \tag{ B23 }$$

where the oscillation component,

$$\begin{equation} \mathcal {O}(\tau ^{-13/2})\propto \tau ^{-13/2\bigg(1\pm \sqrt{1-\frac{24}{169}\big(7+\frac{k^2}{k^2_J}\big)}\bigg)}, \end{equation} \tag{ B24 }$$

decays for any $0\le \frac{k}{k_J}$. The velocity divergence filtering function at the third order is

$$\begin{equation} h_3(\mathbf {k})=\frac{1}{\theta _{3,c}(\mathbf {k})} [A_3(\mathbf {k})-3\delta _{3,c}(\mathbf {k})g_3(\mathbf {k})], \end{equation} \tag{ B25 }$$

where we have ignored the decaying term.

Let us rewrite 7A₃(k) − 2B₃(k) in Equation (B23) as

$$\begin{eqnarray} &&\hspace{-1.5pc}7A_3(\mathbf {k})-2B_3(\mathbf {k}) \nonumber \\ &&\hspace{-0.5pc}= \frac{1}{(2\pi)^6}\int \!\int \!\int d\mathbf {q}_1d\mathbf {q}_2d\mathbf {q}_3 \delta _D(\mathbf {q}_1+\mathbf {q}_2+\mathbf {q}_3-\mathbf {k}) \delta _{1,c}\nonumber \\ &&\quad(\mathbf {q}_1)\delta _{1,c}(\mathbf {q}_2) \delta _{1,c}(\mathbf {q}_3) \nonumber \\ &&\;\,{\times}\, \bigg[ \frac{7\mathbf {k}\cdot \mathbf {q}_1}{q_1^2} g_1(\mathbf {q}_1)g_2(\mathbf {q}_{23}) F_2^{(s)}(\mathbf {q}_2,\mathbf {q}_3)\nonumber \\ && \;\,{+}\,\frac{7\mathbf {k}\cdot \mathbf {q}_{12}}{q_{12}^2} h_2(\mathbf {q}_{12})g_1(\mathbf {q}_3) G_2^{(s)}(\mathbf {q}_1,\mathbf {q}_2) \nonumber \\ &&\;\,{+}\,\frac{k^2(\mathbf {q}_1\cdot \mathbf {q}_{23})}{q_1^2q_{23}^2} g_1(\mathbf {q}_1)h_2(\mathbf {q}_{23}) G_2^{(s)}(\mathbf {q}_2,\mathbf {q}_3)\nonumber \\ &&\;\,{+}\,\frac{k^2(\mathbf {q}_{12}\cdot \mathbf {q}_3)}{q_{12}^2q_3^2} h_2(\mathbf {q}_{12})g_1(\mathbf {q}_3) G_2^{(s)}(\mathbf {q}_1,\mathbf {q}_2) \nonumber \\ &&\;\,{+}\,\frac{3}{2}\frac{k^2}{k_J^2} g_1(\mathbf {q}_1)g_2(\mathbf {q}_{23}) F_2^{(s)}(\mathbf {q}_2,\mathbf {q}_3)\nonumber \\ &&\;\,{+}\,\frac{3}{2}\frac{k^2}{k_J^2} g_2(\mathbf {q}_{12})g_1(\mathbf {q}_3) F_2^{(s)}(\mathbf {q}_1,\mathbf {q}_2)\nonumber \\ && \;\,{-}\,\frac{k^2}{k_J^2} g_1(\mathbf {q}_1)g_1(\mathbf {q}_2)g_1(\mathbf {q}_3) \bigg] \nonumber \\ &&\hspace{-0.5pc}\equiv \frac{18}{(2\pi)^6}\int \!\int \!\int d\mathbf {q}_1d\mathbf {q}_2d\mathbf {q}_3 \delta _D(\mathbf {q}_1+\mathbf {q}_2+\mathbf {q}_3-\mathbf {k}) \mathcal {F}_3\nonumber \\ &&\;\,{\times}\,(\mathbf {q}_1,\mathbf {q}_2,\mathbf {q}_3) \delta _{1,c}(\mathbf {q}_1)\delta _{1,c}(\mathbf {q}_2) \delta _{1,c}(\mathbf {q}_3) \nonumber \\ &&\hspace{-0.5pc}\equiv 18\delta ^{\prime }_{3,c}(\mathbf {k}). \end{eqnarray} \tag{ B26 }$$

The new kernel, ${\cal F}_3(\mathbf {q}_1,\mathbf {q}_2,\mathbf {q}_3)$, can be symmetrized as

$$\begin{eqnarray} &&\mathcal {F}^{(s)}_3(\mathbf {q}_1,\mathbf {q}_2,\mathbf {q}_3) \nonumber \\ &&=\case{1}{6} [ \mathcal {F}_3(\mathbf {q}_1,\mathbf {q}_2,\mathbf {q}_3) +\mathcal {F}_3(\mathbf {q}_1,\mathbf {q}_3,\mathbf {q}_2)\nonumber \\ &&\quad +\mathcal {F}_3(\mathbf {q}_2,\mathbf {q}_1,\mathbf {q}_3) +\mathcal {F}_3(\mathbf {q}_2,\mathbf {q}_3,\mathbf {q}_1)\nonumber \\ &&\quad +\mathcal {F}_3(\mathbf {q}_3,\mathbf {q}_1,\mathbf {q}_2) +\mathcal {F}_3(\mathbf {q}_3,\mathbf {q}_2,\mathbf {q}_1) ] \nonumber \\ &&=\frac{7}{54}\mathbf {k}\cdot \bigg[ F_2^{(s)}(\mathbf {q}_2,\mathbf {q}_3)\frac{\mathbf {q}_1}{q_1^2} g_1(\mathbf {q}_1)g_2(\mathbf {q}_{23})\nonumber \\ &&\quad +F_2^{(s)}(\mathbf {q}_1,\mathbf {q}_3)\frac{\mathbf {q}_2}{q_2^2} g_1(\mathbf {q}_2)g_2(\mathbf {q}_{13})\nonumber \\ &&\quad +F_2^{(s)}(\mathbf {q}_1,\mathbf {q}_2)\frac{\mathbf {q}_3}{q_3^2} g_1(\mathbf {q}_3)g_2(\mathbf {q}_{12}) \bigg] \nonumber \\ &&\quad+\frac{1}{27}k^2\bigg[ G_2^{(s)}(\mathbf {q}_2,\mathbf {q}_3)\frac{\mathbf {q}_1\cdot \mathbf {q}_{23}}{q_1^2q_{23}^2} g_1(\mathbf {q}_1)h_2(\mathbf {q}_{23})\nonumber \\ &&\quad +G_2^{(s)}(\mathbf {q}_1,\mathbf {q}_3)\frac{\mathbf {q}_2\cdot \mathbf {q}_{13}}{q_2^2q_{13}^2} g_1(\mathbf {q}_2)h_2(\mathbf {q}_{13})\nonumber \\ &&\quad +G_2^{(s)}(\mathbf {q}_1,\mathbf {q}_2)\frac{\mathbf {q}_3\cdot \mathbf {q}_{12}}{q_3^2q_{12}^2} g_1(\mathbf {q}_3)h_2(\mathbf {q}_{12}) \bigg] \nonumber \\ &&\quad+\frac{7}{54}\mathbf {k}\cdot \bigg[ G_2^{(s)}(\mathbf {q}_2,\mathbf {q}_3)\frac{\mathbf {q}_{23}}{q_{23}^2} g_1(\mathbf {q}_1)h_2(\mathbf {q}_{23})\nonumber \\ &&\quad +G_2^{(s)}(\mathbf {q}_1,\mathbf {q}_3)\frac{\mathbf {q}_{13}}{q_{13}^2} g_1(\mathbf {q}_2)h_2(\mathbf {q}_{13})\nonumber \\ &&\quad +G_2^{(s)}(\mathbf {q}_1,\mathbf {q}_2)\frac{\mathbf {q}_{12}}{q_{12}^2} g_1(\mathbf {q}_3)h_2(\mathbf {q}_{12}) \bigg] \nonumber \\ &&\quad+\frac{1}{18}\frac{k^2}{k^2_J}\bigg[ g_1(\mathbf {q}_1)g_2(\mathbf {q}_{23}) F_2^{(s)}(\mathbf {q}_2,\mathbf {q}_3)\nonumber \\ &&\quad +g_1(\mathbf {q}_2)g_2(\mathbf {q}_{13}) F_2^{(s)}(\mathbf {q}_1,\mathbf {q}_3)\nonumber \\ &&\quad +g_1(\mathbf {q}_3)g_2(\mathbf {q}_{12}) F_2^{(s)}(\mathbf {q}_1,\mathbf {q}_2)\nonumber \\ &&\quad -g_1(\mathbf {q}_1)g_1(\mathbf {q}_2)g_1(\mathbf {q}_3) \bigg]. \end{eqnarray} \tag{ B27 }$$

In the limit of k_J → ∞, $\mathcal {F}_3\rightarrow F_3$, and g₃(k) = 1. Using δ'_3,c(k) introduced above, we write g₃ as

$$\begin{equation} g_3(\mathbf {k})=\frac{7-6\bigg[1-\frac{\delta _{3,c}^{\prime }(\mathbf {k})}{\delta _{3,c}(\mathbf {k})}\bigg]}{7+\frac{k^2}{k^2_J}}. \end{equation} \tag{ B28 }$$

We calculate the power spectrum of the total matter fluctuations, δ = f_cδ_c + f_bδ_b = f_cδ_c + (1 − f_c)δ_b, which is given, up to the third order in perturbations, by

$$\begin{eqnarray} \delta (\mathbf {k},\tau)&=&f_c\delta _c(\mathbf {k},\tau)+f_b\delta _b(\mathbf {k},\tau) \nonumber \\ &=&f_c[\delta _{1,c}(\mathbf {k},\tau)+\delta _{2,c}(\mathbf {k},\tau)+\delta _{3,c}(\mathbf {k},\tau)]\nonumber \\ && +(1-f_c)[\delta _{1,b}(\mathbf {k},\tau)+\delta _{2,b}(\mathbf {k},\tau)+\delta _{3,b}(\mathbf {k},\tau)] \nonumber \\ &=&f_c[\delta _{1,c}(\mathbf {k},\tau)+\delta _{2,c}(\mathbf {k},\tau)+\delta _{3,c}(\mathbf {k},\tau)]\nonumber \\ && +(1-f_c)[g_1(k)\delta _{1,c}(\mathbf {k},\tau)+g_2(\mathbf {k})\delta _{2,c}(\mathbf {k},\tau)\nonumber \\ &&+g_3(\mathbf {k})\delta _{3,c}(\mathbf {k},\tau)]. \end{eqnarray} \tag{ C1 }$$

The power spectrum is

$$\begin{eqnarray} (2\pi)^3P_{\rm tot}(k)\delta _D(\mathbf {k}+\mathbf {k}^{\prime }) &=&\langle \delta (\mathbf {k})\delta (\mathbf {k}^{\prime })\rangle \nonumber \\ &=&\langle \lbrace f_c\delta _c(\mathbf {k})+(1-f_c)\delta _b(\mathbf {k})\rbrace \lbrace f_c\delta _c(\mathbf {k}^{\prime })\nonumber \\ &&+(1-f_c)\delta _b(\mathbf {k}^{\prime })\rbrace \rangle \nonumber \\ &=&f_c^2\langle \delta _c(\mathbf {k})\delta _c(\mathbf {k}^{\prime })\rangle +2f_c(1\!-\!f_c)\langle \delta _b(\mathbf {k})\delta _c\nonumber \\ &&(\mathbf {k}^{\prime })\rangle +(1-f_c)^2\langle \delta _b(\mathbf {k})\delta _b(\mathbf {k}^{\prime })\rangle \nonumber \\ &\equiv & (2\pi)^3[ f_c^2P_c(k)+2f_c(1-f_c)P_{b,c}(k)\nonumber \\ &&+(1-f_c)^2P_b(k)]\delta _D(\mathbf {k}+\mathbf {k}^{\prime }), \end{eqnarray} \tag{ C2 }$$

where P_c, P_b,c, and P_b are

$$\begin{eqnarray} (2\pi)^3P_{c}(k)\delta _D(\mathbf {k}+\mathbf {k}^{\prime }) &=&\langle \delta _c(\mathbf {k})\delta _c(\mathbf {k}^{\prime })\rangle \nonumber \\ &=&\langle \lbrace \delta _{1,c}(\mathbf {k}) +\delta _{2,c}(\mathbf {k}) +\delta _{3,c}(\mathbf {k})\rbrace \lbrace \delta _{1,c}(\mathbf {k}^{\prime })\nonumber \\ && +\delta _{2,c}(\mathbf {k}^{\prime }) +\delta _{3,c}(\mathbf {k}^{\prime })\rbrace \rangle \nonumber \\ &=&\langle \delta _{1,c}(\mathbf {k})\delta _{1,c}(\mathbf {k}^{\prime })\rangle +2\langle \delta _{1,c}(\mathbf {k})\delta _{3,c}(\mathbf {k}^{\prime })\rangle\nonumber \\ && +\langle \delta _{2,c}(\mathbf {k})\delta _{2,c}(\mathbf {k}^{\prime })\rangle \nonumber \\ &\equiv &(2\pi)^3[P_{11,c}(k)+2P_{13,c}(k)+P_{22,c}(k)]\nonumber \\ &&\times \delta _D(\mathbf {k}+\mathbf {k}^{\prime }), \end{eqnarray} \tag{ C3 }$$

$$\begin{eqnarray} (2\pi)^3P_{bc}(k)\delta _D(\mathbf {k}+\mathbf {k}^{\prime }) &=&\langle \delta _b(\mathbf {k})\delta _c(\mathbf {k}^{\prime })\rangle \nonumber \\ &=&\langle \lbrace g_1(k)\delta _{1,c}(\mathbf {k}) +g_2(\mathbf {k})\delta _{2,c}(\mathbf {k})\nonumber \\ && +g_3(\mathbf {k})\delta _{3,c}(\mathbf {k})\rbrace \lbrace \delta _{1,c}(\mathbf {k}^{\prime })\nonumber \\ && +\delta _{2,c}(\mathbf {k}^{\prime }) +\delta _{3,c}(\mathbf {k}^{\prime })\rbrace \rangle \nonumber \\ &=&g_1(k)\langle \delta _{1,c}(\mathbf {k})\delta _{1,c}(\mathbf {k}^{\prime })\rangle\nonumber \\ && +g_1(k)\langle \delta _{1,c}(\mathbf {k})\delta _{3,c}(\mathbf {k}^{\prime })\rangle\nonumber \\ && +\langle g_3(\mathbf {k})\delta _{3,c}(\mathbf {k})\delta _{1,c}(\mathbf {k}^{\prime })\rangle\nonumber \\ && +\langle g_2(\mathbf {k})\delta _{2,c}(\mathbf {k})\delta _{2,c}(\mathbf {k}^{\prime })\rangle \nonumber \\ &\equiv &(2\pi)^3[P_{11,bc}(k)+2P_{13,bc}(k)\nonumber \\ &&+P_{22,bc}(k)] \delta _D(\mathbf {k}+\mathbf {k}^{\prime }), \end{eqnarray} \tag{ C4 }$$

$$\begin{eqnarray} (2\pi)^3P_{b}(k)\delta _D(\mathbf {k}+\mathbf {k}^{\prime }) &=&\langle \delta _b(\mathbf {k})\delta _b(\mathbf {k}^{\prime })\rangle \nonumber \\ &=&\langle \lbrace g_1(k)\delta _{1,c}(\mathbf {k}) +g_2(\mathbf {k})\delta _{2,c}(\mathbf {k})\nonumber \\ && +g_3(\mathbf {k})\delta _{3,c}(\mathbf {k})\rbrace \times \lbrace g_1(k^{\prime })\delta _{1,c}(\mathbf {k}^{\prime })\nonumber \\ && +g_2(\mathbf {k}^{\prime })\delta _{2,c}(\mathbf {k}^{\prime }) +g_3(\mathbf {k}^{\prime })\delta _{3,c}(\mathbf {k}^{\prime })\rbrace \rangle \nonumber \\ &=&g_1^2(k)\langle \delta _{1,c}(\mathbf {k})\delta _{1,c}(\mathbf {k}^{\prime })\rangle\nonumber \\ && +2g_1(k)\langle \delta _{1,c}(\mathbf {k})g_3(\mathbf {k}^{\prime })\delta _{3,c}(\mathbf {k}^{\prime })\rangle\nonumber \\ && +\langle g_2(\mathbf {k})\delta _{2,c}(\mathbf {k})g_2(\mathbf {k}^{\prime })\delta _{2,c}(\mathbf {k}^{\prime })\rangle \nonumber \\ &\equiv &(2\pi)^3[P_{11,b}(k)+2P_{13,b}(k)+P_{22,b}(k) ]\nonumber \\ &&\times \delta _D(\mathbf {k}+\mathbf {k}^{\prime }), \end{eqnarray} \tag{ C5 }$$

respectively.

Now, P_11,c(k), P_13,c(k), and P_22,c(k) can be numerically calculated with the corresponding kernels, F^(s)₂ and F^(s)₃:

$$\begin{eqnarray} (2\pi)^3P_{11,bc}(k)\delta _D(\mathbf {k}+\mathbf {k}^{\prime }) &=&\langle \delta _{1,b}(\mathbf {k})\delta _{1,c}(\mathbf {k}^{\prime })\rangle \nonumber \\ &=&g_1(k)\langle \delta _{1,c}(\mathbf {k})\delta _{1,c}(\mathbf {k}^{\prime })\rangle, \end{eqnarray} \tag{ C6 }$$

$$\begin{eqnarray} (2\pi)^3P_{13,bc}(k)\delta _D(\mathbf {k}+\mathbf {k}^{\prime }) &=&\case{1}{2}[\langle \delta _{1,b}(\mathbf {k})\delta _{3,c}(\mathbf {k}^{\prime })\rangle\nonumber \\ && +\langle \delta _{1,c}(\mathbf {k})\delta _{3,b}(\mathbf {k}^{\prime })\rangle] \nonumber \\ &=&\case{1}{2} [ g_1(k)\langle \delta _{1,c}(\mathbf {k})\delta _{3,c}(\mathbf {k}^{\prime })\rangle\nonumber \\ && +\langle \delta _{1,c}(\mathbf {k})g_3(\mathbf {k}^{\prime })\delta _{3,c}(\mathbf {k}^{\prime })\rangle] \nonumber \\ &=&\frac{1}{2}\bigg[\bigg(g_1(k)+\frac{1}{7+\frac{k^2}{k^2_J}}\bigg)\nonumber \\ &&\times \langle \delta _{1,c}(\mathbf {k})\delta _{3,c}(\mathbf {k}^{\prime })\rangle +\frac{6}{7+\frac{k^2}{k^2_J}}\nonumber \\ &&\times \langle \delta _{1,c}(\mathbf {k})\delta ^{\prime }_{3,c}(\mathbf {k}^{\prime })\rangle \bigg], \end{eqnarray} \tag{ C7 }$$

$$\begin{eqnarray} (2\pi)^3P_{22,bc}(k)\delta _D(\mathbf {k}+\mathbf {k}^{\prime }) &=&\langle \delta _{2,b}(\mathbf {k})\delta _{2,c}(\mathbf {k}^{\prime })\rangle \nonumber \\ &=&\langle g_2(\mathbf {k})\delta _{2,c}(\mathbf {k})\delta _{2,c}(\mathbf {k}^{\prime })\rangle \nonumber \\ &=&\frac{1}{\frac{10}{3}+\frac{k^2}{k^2_J}} \bigg[\langle \delta _{2,c}(\mathbf {k})\delta _{2,c}(\mathbf {k}^{\prime })\rangle\nonumber \\ && +\frac{7}{3}\langle \delta ^{\prime }_{2,c}(\mathbf {k})\delta _{2,c}(\mathbf {k}^{\prime })\rangle \bigg], \end{eqnarray} \tag{ C8 }$$

$$\begin{eqnarray} (2\pi)^3P_{11,b}(k)\delta _D(\mathbf {k}+\mathbf {k}^{\prime }) &=&\langle \delta _{1,b}(\mathbf {k})\delta _{1,b}(\mathbf {k}^{\prime })\rangle \nonumber \\ &=&g_1^2(k)\langle \delta _{1,c}(\mathbf {k})\delta _{1,c}(\mathbf {k}^{\prime })\rangle, \end{eqnarray} \tag{ C9 }$$

$$\begin{eqnarray} (2\pi)^3P_{13,b}(k)\delta _D(\mathbf {k}+\mathbf {k}^{\prime }) &=&\langle \delta _{1,b}(\mathbf {k})\delta _{3,b}(\mathbf {k}^{\prime })\rangle \nonumber \\ &=&g_1(k)\langle \delta _{1,c}(\mathbf {k})g_3(\mathbf {k}^{\prime })\delta _{3,c}(\mathbf {k}^{\prime })\rangle \nonumber \\ &=&\frac{g_1(k)}{7+\frac{k^2}{k^2_J}} [\langle \delta _{1,c}(\mathbf {k})\delta _{3,c}(\mathbf {k}^{\prime })\rangle\nonumber \\ && +6\langle \delta _{1,c}(\mathbf {k})\delta ^{\prime }_{3,c}(\mathbf {k}^{\prime })\rangle ], \end{eqnarray} \tag{ C10 }$$

$$\begin{eqnarray} (2\pi)^3P_{22,b}(k)\delta _D(\mathbf {k}+\mathbf {k}^{\prime }) &=&\langle \delta _{2,b}(\mathbf {k})\delta _{2,b}(\mathbf {k}^{\prime })\rangle \nonumber \\ &=&\langle g_2(\mathbf {k})\delta _{2,c}(\mathbf {k})g_2(\mathbf {k}^{\prime })\delta _{2,c}(\mathbf {k}^{\prime })\rangle \nonumber \\ &=&\frac{1}{\big(\frac{10}{3}+\frac{k^2}{k^2_J}\big)^2} \bigg[ \langle \delta _{2,c}(\mathbf {k})\delta _{2,c}(\mathbf {k}^{\prime })\rangle\nonumber \\ && +\frac{14}{3}\langle \delta _{2,c}(\mathbf {k})\delta ^{\prime }_{2,c}(\mathbf {k}^{\prime })\rangle\nonumber \\ && +\frac{49}{9}\langle \delta ^{\prime }_{2,c}(\mathbf {k})\delta ^{\prime }_{2,c}(\mathbf {k}^{\prime })\rangle\bigg]. \end{eqnarray} \tag{ C11 }$$

The ensemble averages of the products involving δ'_n,c(k) are given by

$$\begin{eqnarray} \langle \delta _{1,c}(\mathbf {k})\delta ^{\prime }_{3,c}(\mathbf {k}^{\prime })\rangle &=&3\delta _D(\mathbf {k}+\mathbf {k}^{\prime })P_{11,c}(k)\nonumber \\ && \times\int d\mathbf {q}\mathcal {F}_3^{(s)}(\mathbf {q},-\mathbf {q},\mathbf {k}) P_{11,c}(q)\nonumber \\ &=&6\pi \delta _D(\mathbf {k}+\mathbf {k}^{\prime })P_{11,c}(k) \int ^{\infty }_0 dq\ q^2P_{11,c}(q)\nonumber \\ &&\times\int ^1_{-1}d\mu \mathcal {F}_3^{(s)}(\mathbf {q},-\mathbf {q},\mathbf {k}), \end{eqnarray} \tag{ C12 }$$

$$\begin{eqnarray} \langle \delta _{2,c}(\mathbf {k})\delta ^{\prime }_{2,c}(\mathbf {k}^{\prime })\rangle &=& 2\delta _D(\mathbf {k}+\mathbf {k}^{\prime })\nonumber \\ &&\times\int d\mathbf {q} P_{11,c}(q)P_{11,c}(|\mathbf {k}-\mathbf {q}|) F_2^{(s)}\nonumber \\ &&\times(\mathbf {q},\mathbf {k}-\mathbf {q}) \mathcal {F}_2^{(s)}(\mathbf {q},\mathbf {k}-\mathbf {q}), \end{eqnarray} \tag{ C13 }$$

$$\begin{eqnarray} \langle \delta ^{\prime }_{2,c}(\mathbf {k})\delta ^{\prime }_{2,c}(\mathbf {k}^{\prime })\rangle &=&2\delta _D(\mathbf {k}+\mathbf {k}^{\prime })\int d\mathbf {q} P_{11,c}(q)P_{11,c}(|\mathbf {k}-\mathbf {q}|)\nonumber \\ &&\times \big[ \mathcal {F}_2^{(s)}(\mathbf {q},\mathbf {k}-\mathbf {q}) \big]^2. \end{eqnarray} \tag{ C14 }$$

Here, the term, $\int \frac{d\mathbf {q}}{(2\pi)^3}\mathcal {F}_3(\mathbf {q},-\mathbf {q},\mathbf {k})P_{11,c}(q)$, in Equation (C12) is given by

$$\begin{eqnarray} &&\!\!\!\!\!\!\int \frac{d\mathbf {q}}{(2\pi)^3}\mathcal {F}^{(s)}_3(\mathbf {q},-\mathbf {q},\mathbf {k})P_{11,c}(q) \nonumber \\ &&=\int \frac{d\mathbf {q}}{(2\pi)^3}\bigg\lbrace \frac{7}{54}\mathbf {k}\cdot \bigg[ F_2^{(s)}(-\mathbf {q},\mathbf {k})\frac{\mathbf {q}}{q^2} g_1(\mathbf {q})g_2(\mathbf {k}-\mathbf {q}) \nonumber \\ &&\quad-F_2^{(s)}(\mathbf {q},\mathbf {k})\frac{\mathbf {q}}{q^2} g_1(\mathbf {q})g_2(\mathbf {k}+\mathbf {q}) \bigg] \nonumber \\ &&\quad\,{+}\,\frac{2}{27}k^2\bigg[ F_2^{(s)}(-\mathbf {q},\mathbf {k})\frac{\mathbf {q}\cdot (\mathbf {k}-\mathbf {q})}{q^2(\mathbf {k}-\mathbf {q})^2} g_1(\mathbf {q})g_2(\mathbf {k}-\mathbf {q})\nonumber \\ && \quad-F_2^{(s)}(\mathbf {q},\mathbf {k})\frac{\mathbf {q}\cdot (\mathbf {k}+\mathbf {q})}{q^2(\mathbf {k}+\mathbf {q})^2} g_1(\mathbf {q})g_2(\mathbf {k}+\mathbf {q}) \bigg] \nonumber \\ &&\quad\,{+}\,\frac{14}{54}\mathbf {k}\cdot \bigg[ F_2^{(s)}(-\mathbf {q},\mathbf {k})\frac{\mathbf {k}-\mathbf {q}}{(\mathbf {k}-\mathbf {q})^2} g_1(\mathbf {q})g_2(\mathbf {k}-\mathbf {q})\nonumber \\ && \quad\,{+}\,F_2^{(s)}(\mathbf {q},\mathbf {k})\frac{\mathbf {k}+\mathbf {q}}{(\mathbf {k}+\mathbf {q})^2} g_1(\mathbf {q})g_2(\mathbf {k}+\mathbf {q}) \bigg] \nonumber \\ &&\quad-\frac{1}{27}k^2\bigg[ \bigg(1+\frac{(-\mathbf {q}\cdot \mathbf {k})(q^2+k^2)}{2q^2k^2} \bigg) \frac{\mathbf {q}\cdot (\mathbf {k}-\mathbf {q})}{q^2(\mathbf {k}-\mathbf {q})^2} g_1^2(\mathbf {q})g_1(\mathbf {k})\nonumber \\ && \quad-\bigg(1+\frac{(\mathbf {q}\cdot \mathbf {k})(q^2+k^2)}{2q^2k^2} \bigg) \frac{\mathbf {q}\cdot (\mathbf {k}+\mathbf {q})}{q^2(\mathbf {k}+\mathbf {q})^2} g_1^2(\mathbf {q})g_1(\mathbf {k}) \bigg] \nonumber \\ &&\quad\,{-}\,\frac{7}{54}\mathbf {k}\cdot \bigg[ \bigg(1+\frac{(-\mathbf {q}\cdot \mathbf {k})(q^2+k^2)}{2q^2k^2} \bigg) \frac{\mathbf {k}-\mathbf {q}}{(\mathbf {k}-\mathbf {q})^2} g_1^2(\mathbf {q})g_1(\mathbf {k})\nonumber \\ && \quad\,{+}\,\bigg(1+\frac{(\mathbf {q}\cdot \mathbf {k})(q^2+k^2)}{2q^2k^2} \bigg) \frac{\mathbf {k}+\mathbf {q}}{(\mathbf {k}+\mathbf {q})^2} g_1^2(\mathbf {q})g_1(\mathbf {k}) \bigg] \nonumber \\ &&\quad\,{+}\,\frac{1}{18}\frac{k^2}{k^2_J}\bigg[ g_1(\mathbf {q})g_2(\mathbf {k}-\mathbf {q}) F_2^{(s)}(-\mathbf {q},\mathbf {k})\nonumber \\ && \quad\,{+}\,g_1(\mathbf {q})g_2(\mathbf {k}+\mathbf {q}) F_2^{(s)}(\mathbf {q},\mathbf {k})\nonumber \\ && \quad\,{-}\,g_1^2(\mathbf {q})g_1(\mathbf {k}) \bigg]\bigg\rbrace P_{11,c}(q), \end{eqnarray} \tag{ C15 }$$

where we have used Equation (B17) and F^(s)₂(q, − q) = G^(s)₂(q, − q) = 0. We then calculate the angular average of $\mathcal {F}_3^{(s)}$, i.e., $\int d\mu \mathcal {F}_3^{(s)}$, for the linear filtering function of g₁(k) = 1/(1 + k²/k²_J):

$$\begin{eqnarray} \int ^1_{-1}d\mu \mathcal {F}^{(s)}_3&=& \frac{1}{612360r^8s(1+r^2)(r^2+s^2)^2}\bigg[\bigg[30r^2s^3[-14000s^6\nonumber \\ &&+810r^{10}(1+s^2) +900r^2s^4(-7+5s^2) \nonumber \\ &&+60r^4s^2(105-125s^2+78s^4)\nonumber \\ &&+9r^8(321-248s^2+159s^4)\nonumber \\ &&+27r^6(126-87s^2+70s^4+9s^6)] \nonumber \\ &&-243r^8(-7+5s^2+2s^4)[5(r^4+s^2)(-1+s^2)^2\nonumber \\ &&+r^2(5-5s^2-19s^4+5s^6)]\ln \frac{1+s}{|1-s|} \nonumber \\ &&+[10s^2+3r^2(1+s^2)][-35s^2\nonumber \\ &&+3r^2(-7+s^2)] [-2000s^6+135r^8(-1+s^2)^2 \nonumber \\ &&+240r^4s^2(3-4s^2+3s^4)+300r^2(s^4+s^6)\nonumber \\ &&+27r^6(5+5s^2-9s^4+5s^6)]\nonumber \\ &&\times \frac{1}{2}\ln \bigg[\frac{10s^2+3r^2(1+s)^2}{10s^2+3r^2(1-s)^2}\bigg]\bigg]\bigg], \end{eqnarray} \tag{ C16 }$$

where r ≡ k/k_J and s ≡ k/q. We find that the calculation of $\mathcal {F}_3$ is numerically unstable as k/k_J → 0 (r → 0). The exact limit of $\mathcal {F}_3$ is $\lim _{k/k_J\rightarrow 0}\mathcal {F}_3\rightarrow F_3$, and thus one may replace $\mathcal {F}_3$ with F₃ for a sufficiently small value of k/k_J.

Finally, we generalize the above results from an EdS universe to general cosmological models, by writing

$$\begin{eqnarray} &&\hspace{-1.4pc}\frac{a^2(\tau)}{a^2(\tau _i)} P_{11}(k,\tau _i) \rightarrow P_{11}(k,\tau) = \frac{D^2(\tau)}{D^2(\tau _i)}\nonumber \\ &&\hspace{-0.3pc} \times\bigg(\frac{\delta _{1,c+}^{(1)}(k,\tau)/\delta _{1,c+}^{(0)}(k,\tau)}{\delta _{1,c+}^{(1)}(k,\tau _*)/\delta _{1,c+}^{(0)}(k,\tau _*)}\bigg)^2 P_{11}(k,\tau _i), \end{eqnarray} \tag{ C17 }$$

where τ_i is some arbitrary epoch, τ_* is the epoch where the pressure effect becomes non-negligible (i.e., re-ionization epoch for baryons and nonrelativistic transition for massive neutrinos), and D(τ) is the linear growth factor appropriate to a given cosmological model. We obtain Equation (58) from combining Equations (C8), (C13), and P_22,c given by Equation (A34). Similarly, we obtain Equation (59) from combining Equations (C11), (C13), (C14), and P_22,c, Equation (60) from combining Equations (C7), (C12), and P_13,c given by Equation (A36), and Equation (61) from combining Equations (C10), (C12), and P_13,c.

Figure 5 shows the dimensionless 3PT and linear power spectra, Δ²(k) = k³P(k)/(2π²), for a matter component with pressure at different redshifts (z = 0.1, 1.0, 3.0, 5.0, 10, and 30) with k_J = 1.0 and 3.0 h Mpc^-1. The 3PT and linear power spectra are similar at the highest redshift, whereas the 3PT has significantly more power than the linear spectrum at larger wavenumbers as we go to lower redshifts. As a result, the filtering scale for a given linear filtering scale migrates toward larger wavenumbers in lower redshifts.

Zoom In Zoom Out Reset image size