© ESO 2022

1. Introduction

The accepted orthodoxy in cosmology is that the Universe prior to a redshift of 1500 consisted partly of fully ionised baryonic plasma at a temperature above 4000° K in thermal equilibrium, a consequent black body radiation field with a radiation pressure far exceeding the thermal pressure of the plasma, and a dark matter component with a mass an order of magnitude greater than that of the baryons. Due to the radiation pressure acting on the baryons, the fluctuations in baryon density, inferred from the anisotropies in the cosmic microwave background (CMB) radiation, would be unable to collapse to form galaxies and galaxy clusters, since the pressure was such that the Jeans length (i.e. the minimum scale for a density enhancement to collapse under gravity) would be greater than the scale of the Universe.

Subsequent to redshift 1500, the temperature of the expanding plasma would fall to below the temperature needed for hydrogen recombination to occur in thermal equilibrium, and baryonic matter would then decouple from the radiation field after this era of recombination. Consequently, radiation pressure would no longer prevent the baryon density perturbations from collapsing under gravity, and the era of galaxy and cluster formation would begin.

Large-scale cosmological computer simulations start their initial conditions after the era of full recombination at redshift ∼100 and assume that there is coincidence between dark matter and the baryonic perturbations that are revealed by the CMB anisotropies; in other words, they assume that the positions and amplitudes of peaks and troughs of density perturbations in terms of a percentage of the ambient density are identical for baryonic and dark matter (Vogelsberger et al. 2014; Schaye et al. 2015; Kaviraj et al. 2017). However, there is a serious issue with this assumption: though we do not know the specific nature of dark matter, it seems certain that it does not absorb or emit electromagnetic radiation and hence cannot be affected at all by the radiation pressure. This means that prior to the era of recombination, perturbations of dark matter cannot be supported by radiation pressure, while baryonic perturbations will be supported. Furthermore, standard cosmology postulates cold dark matter (ΛCDM), so there is also no internal kinetic energy that would be able to prevent the collapse. Therefore, dark matter should be subject to gravitational collapse on all scales. This would leave any assumption of coincidence at recombination between baryonic matter perturbations, as evidenced by the electron density perturbations deduced from the CMB, and dark matter perturbations open to question.

To quantify this idea, it is useful to carry out a simple plane wave analysis of the perturbation of a system consisting of baryonic plasma, black body radiation, and dark matter. This analysis considers the period between horizon crossing of the perturbation wavelength and recombination. A similar analysis has been the basis of the baryon acoustic model for the interpretation of the CMB, as exemplified by Dodelson (2003). In this paper we follow a similar physical model for the baryons and the dark matter but use a different physical paradigm for the radiation. Dodelson uses a kinetic theory model for the photons of the radiation field, in the form of the Einstein Boltzmann equations, leading to a hierarchy of moment equations for the radiation temperature perturbation. We assume strict thermodynamic equilibrium, with the radiation having a black body spectrum and the energy density given by the Stefan–Boltzmann law based on the baryon temperature.

It should be noted that the fully ionised baryon plasma does not support simple acoustic waves, but, in common with contemporary plasmas, it supports a variety of waves, such as Langmuir waves (also known as plasma waves, where the electrons oscillate rapidly, independent of the ions), ion-acoustic waves (Boyd & Sanderson 1969), and, if a magnetic field is present in the primordial plasma (Ruis-Granados et al. 2016), a variety of magnetohydrodynamic waves. However, ion acoustic waves are the nearest mode to a simple acoustic wave, and the plasma frequency at recombination is many magnitudes higher than the frequency of ion acoustic waves at the relevant wavelengths, so the electrostatic binding of electrons and ions is very strong. Therefore, in common with Dodelson’s analysis, we assume that a simple acoustic-like wave is the origin of the CMB anisotropies and follow a similar set of equations for baryons and dark matter as that on page 111 of Dodelson (2003).

In this paper we perform an eigenmode analysis of the plane waves supported by this model. This analysis highlights the very weak coupling between the baryon and dark matter perturbations and calls into question the assumptions behind the initial conditions for simulations starting at redshift ∼100, and indeed the analysis of the CMB anisotropies in terms of coincident baryon and dark matter perturbations.

2. Basic equations with black body radiation pressure

Our assumption that the radiation is in strict thermodynamic equilibrium with the baryonic plasma is based on the fact that we considered only scales much greater than the mean free path of a photon, meaning that the radiation and plasma move together; in other words, we used the tight coupling approximation. This neglects the effect of Silk damping at short wavelengths (Silk 1967) but yields a simple, but plausible, physical model that identifies the different baryon and dark matter wave modes. The energy density of the radiation is then given by the Stefan–Boltzmann law, leading to the associated isotropic, but spatially varying, radiation pressure.

We note that under the tight coupling approximation, discussed by Dodelson (2003) on page 226, the radiation dipole term in the equation of motion of the baryons in his formalism, proportional to Θ₁, can be shown, using the equation of motion for Θ₁, to reduce to a gradient of the monopole term (i.e. ikΘ₀) when the frequency of a wave is much lower than the photon collision frequency. Consequently, under these assumptions, the two radiation formalisms are equivalent.

The total thermal energy density is then

$$\begin{array}{c}\hfill U=U_R+U_B,\end{array}$$

where

$$\begin{array}{c}\hfill U_B=\frac{3}{2}{n}_{B}KT=3\frac{{\rho }_B}{m_p}KT.\end{array}$$

Here, ρ_B is the density of baryons, m_p is the proton mass, T is the baryonic temperature, and K is Boltzmann’s constant, while

$$\begin{array}{c}\hfill U_R=4\frac{\sigma T^4}{c},\end{array}$$

where σ is Stefan’s constant, and c is the speed of light.

With T = 4000° K, and ρ_B ∼ 10⁻¹⁷ kg m⁻³ at z = 1500 we have U_R ≫ U_B.

Hence, we deemed the work done by radiation pressure on baryonic compressions to go wholly into U_R since when T increases radiation will be emitted, increasing U_R. We therefore used a conservation equation for the U_R effectively carried by the baryons, with the addition of a source term that represents adiabatic compressional heating together with a relativistic gravitational correction term.

The perturbed, expanding flat Universe metric in the Newtonian conformal gauge is given by

$$\begin{array}{cc}\hfill g_{00}(\mathbf{x},t)& =-1-2\psi (\mathbf{x},t)\hfill \\ \hfill g_{0i}(\mathbf{x},t)& =0\hfill \\ \hfill g_{\mathit{ij}}(\mathbf{x},t)& =a{\left(t\right)}^2{\delta }_{\mathit{ij}}(1+2\varphi (\mathbf{x},t)),\hfill \end{array}$$

where a(t) is the cosmological scale factor, ψ is the Newtonian potential, and ϕ is the perturbation of the spatial curvature. We ignored any neutrino contribution since such a weakly interacting component will have little effect on the relationship between baryons and dark matter, and, since the radiation field is isotropic, we have (see Eq. (5.33) of Dodelson 2003)

$$\begin{array}{c}\hfill \varphi =-\psi .\end{array}$$

Hence, the set of equations governing a combined baryonic and dark matter fluid accompanied by a black body radiation field with dominant radiation pressure, replacing Eqs. (4.100) to (4.106) of Dodelson, is

$$\begin{array}{cc}\hfill \frac{\partial {\rho }_B}{\partial t}+▿·\left({\rho }_B{\mathbf{v}}_B\right)& =-\frac{3}{c^2}{\rho }_B\frac{\partial \varphi }{\partial t}\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill {\rho }_B(\frac{\partial }{\partial t}+{\mathbf{v}}_B·▿){\mathbf{v}}_B+▿p+\frac{\mathrm{d}a/\mathrm{d}t}{a}{\rho }_B{\mathbf{v}}_B& =-{\rho }_B▿\varphi \hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \frac{\partial {\rho }_D}{\partial t}+▿·\left({\rho }_D{\mathbf{v}}_D\right)& =-\frac{3}{c^2}{\rho }_D\frac{\partial \varphi }{\partial t}\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill {\rho }_D(\frac{\partial }{\partial t}+{\mathbf{v}}_D·▿){\mathbf{v}}_D+\frac{\mathrm{d}a/\mathrm{d}t}{a}{\rho }_D{\mathbf{v}}_D& =-{\rho }_D▿\varphi \hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill \frac{\partial U_R}{\partial t}+▿·\left(U_R{\mathbf{v}}_B\right)+p▿·{\mathbf{v}}_B& =-\frac{3}{c^2}{U}_R\frac{\partial \varphi }{\partial t}\hfill \end{array}$$(5)

and finally

$$\begin{array}{c}\hfill {▿}^2\varphi -\frac{3}{c^2}\frac{\frac{\mathrm{d}a}{\mathrm{d}t}}{a}(\frac{\partial \varphi }{\partial t}+\frac{\frac{\mathrm{d}a}{\mathrm{d}t}}{a}\varphi )=4\pi G({\rho }_D+{\rho }_B+\frac{U_R}{c^2})·\end{array}$$(6)

Here, v is velocity and

$$\begin{array}{c}\hfill p=\frac{4}{3}\frac{\sigma T^4}{c}\gg n_{B}KT.\end{array}$$

3. Plane wave perturbation solutions with frequency $\gg \frac{\mathrm{d}a/\mathrm{d}t}{a}$

We firstly ignore the rate of change of the cosmological scale factor terms, under the assumption that we are considering only perturbations with timescales much smaller than the expansion time, so the expansion rate is negligible compared to the rate of change of the perturbations. The terms proportional to $\frac{\partial \varphi }{\partial t}$ on the right-hand side of the conservation equations, Eqs. (1), (3), and (5), are the only relativistic correction terms required under this assumption and will only be non-negligible in the case of a wave mode where the phase velocity is comparable to c.

Linearising with respect to the small perturbations of the principal physical quantities (i.e. ${\overline{\rho }}_B$, ${\overline{\mathbf{v}}}_B$, ${\overline{\rho }}_D$, ${\overline{\mathbf{v}}}_D$, $\overline{\varphi }$, and $\overline{T}$) and then assuming longitudinal plane waves propagating in the x direction (i.e. all perturbations ∝exp(iωt−ikx) with ${\overline{\mathbf{v}}}_{B,D}=({\overline{v}}_{B,D},o,o)$), we have from Eqs. (1) to (6)

$$\begin{array}{cc}\hfill i\omega {\overline{\rho }}_B-ik{\rho }_B{\overline{v}}_B& =-\frac{3}{c^2}{\rho }_{B}i\omega \overline{\varphi }\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill i\omega {\rho }_B{\overline{v}}_B-ik\frac{16\sigma T^3}{3c}\overline{T}& ={\rho }_{B}ik\overline{\varphi }\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill i\omega {\overline{\rho }}_D-ik{\rho }_D{\overline{v}}_D& =-\frac{3}{c^2}{\rho }_{D}i\omega \overline{\varphi }\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill i\omega {\rho }_D{\overline{v}}_D& ={\rho }_{D}ik\overline{\varphi }\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill i\omega \frac{16\sigma T^3}{c}\overline{T}-ik\frac{16\sigma T^4}{3c}{\overline{v}}_B& =-\frac{12\sigma T^4}{c^3}i\omega \overline{\varphi }\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill {\text{k}}^2\overline{\varphi }& =-4\pi G({\overline{\rho }}_D+{\overline{\rho }}_B+\frac{16\sigma T^3}{c^3}\overline{T}).\hfill \end{array}$$(12)

We need to now reduce these equations to a single equation for one of the six physical quantities in which the multiplying expression will yield a dispersion relation for frequency, ω, as a function of wavenumber, k .

From Eq. (10),

$$\begin{array}{c}\hfill \overline{\varphi }=\frac{\omega }{k}{\overline{v}}_D,\end{array}$$(13)

and from Eq. (7),

$$\begin{array}{c}\hfill {\overline{v}}_B=\frac{\omega }{k}\frac{{\overline{\rho }}_B}{{\rho }_B}+\frac{3}{c^2}\frac{\omega }{k}\overline{\varphi },\end{array}$$(14)

while from Eq. (11)

$$\begin{array}{c}\hfill \overline{T}=\frac{k}{3\omega }T{\overline{v}}_B-\frac{3}{4c^2}T\overline{\varphi }.\end{array}$$

Using Eq. (14), this becomes

$$\begin{array}{c}\hfill \overline{T}=\frac{T}{3}\frac{{\overline{\rho }}_B}{{\rho }_B}+\frac{1}{4c^2}T\overline{\varphi }.\end{array}$$(15)

Therefore, substituting for $\overline{T}$ in Eq. (8) we have

$$\begin{array}{c}\hfill \omega {\rho }_B{\overline{v}}_B-k\frac{16\sigma T^4}{9c}\frac{{\overline{\rho }}_B}{{\rho }_B}-k\frac{4\sigma T^4}{3c^3}\overline{\varphi }={\rho }_{B}k\overline{\varphi },\end{array}$$

and using Eq. (14) we have

$$\begin{array}{c}\hfill (\frac{{\omega }^2}{k^2}-\alpha c^2){\overline{\rho }}_B=-(\frac{3}{c^2}\frac{{\omega }^2}{k^2}-1-\frac{3}{4}\alpha ){\rho }_B\overline{\varphi },\end{array}$$(16)

where

$$\begin{array}{c}\hfill \alpha =\frac{16\sigma T^4}{9{\rho }_B{c}^3}·\end{array}$$

From Eq. (9),

$$\begin{array}{c}\hfill {\overline{\rho }}_D=\frac{k}{\omega }{\rho }_D{\overline{v}}_D-\frac{3}{c^2}{\rho }_D\overline{\varphi }\end{array}$$(17)

and, using Eq. (13),

$$\begin{array}{c}\hfill {\overline{\rho }}_D=(\frac{k^2}{{\omega }^2}-\frac{3}{c^2}){\rho }_D\overline{\varphi .}\end{array}$$(18)

From Eq. (15), Eq. (12) becomes

$$\begin{array}{c}\hfill k^2\overline{\varphi }=-4\pi G({\overline{\rho }}_D+{\overline{\rho }}_B(1+3\alpha )+\frac{4\sigma T^4}{c^5}\overline{\varphi }).\end{array}$$

Therefore, substituting for ${\overline{\rho }}_B$ and ${\overline{\rho }}_D$ from Eqs. (16) and (18), respectively, we obtain

$$\begin{array}{cc}& \{k^2+4\pi G(\frac{k^2}{{\omega }^2}-\frac{3}{c^2}){\rho }_D-4\pi G(1+3\alpha )\frac{(\frac{3}{c^2}\frac{{\omega }^2}{k^2}-1-\frac{3\alpha }{4})}{(\frac{{\omega }^2}{k^2}-\alpha c^2)}{\rho }_B\hfill \\ \hfill & +\frac{9\pi G\alpha {\rho }_B}{c^2}\}\overline{\varphi }=0.\hfill \end{array}$$

The expression inside the curly brackets must therefore be zero and represents the dispersion relation for the perturbations. Multiplying it by ${\omega }^2(\frac{{\omega }^2}{k^2}-\alpha c^2)$ and defining the characteristic frequencies (which are actually inverse timescales),

$$\begin{array}{cc}& {\omega }_D^2=4\pi G{\rho }_D\mathrm{and}{\omega }_B^2=4\pi G{\rho }_B.\hfill \end{array}$$

We then obtain the following quadratic equation for ω²:

$$\begin{array}{cc}& [1-\frac{3{\omega }_D^2}{k^2{c}^2}-3(1+\frac{9\alpha }{4})\frac{{\omega }_B^2}{k^2{c}^2}]{\omega }^4\hfill \\ \hfill & -[\alpha k^2{c}^2-(1+3\alpha ){\omega }_D^2-(1+\frac{15\alpha }{4}){\omega }_B^2]{\omega }^2-\alpha k^2{c}^2{\omega }_D^2=0.\hfill \end{array}$$

This has the solutions

$$\begin{array}{c}\hfill {\omega }^2=\left\{\frac{[1-(1+3\alpha )\frac{{\omega }_D^2}{\alpha k^2{c}^2}-(1+\frac{15\alpha }{4})\frac{{\omega }_B^2}{\alpha k^2{c}^2}]\pm \sqrt[]{A+B}}{2[1-\frac{3{\omega }_D^2}{k^2{c}^2}-3(1+\frac{9\alpha }{4})\frac{{\omega }_B^2}{k^2{c}^2}]}\right\}\alpha k^2{c}^2,\end{array}$$(19)

where

$$\begin{array}{c}\hfill A={[1-(1+3\alpha )\frac{{\omega }_D^2}{\alpha k^2{c}^2}-(1+\frac{15\alpha }{4})\frac{{\omega }_B^2}{\alpha k^2{c}^2}]}^2\end{array}$$

and

$$\begin{array}{c}\hfill B=\frac{4{\omega }_D^2}{\alpha k^2{c}^2}[1-\frac{3{\omega }_D^2}{k^2{c}^2}-3(1+\frac{9\alpha }{4})\frac{{\omega }_B^2}{k^2{c}^2}]·\end{array}$$

Now at z = 1500, ρ_B ∼ 10⁻¹⁷ kg m⁻³ and ρ_D ∼ 10⁻¹⁶ kg m⁻³; therefore, ${\omega }_B^2\sim {10}^{-26}$ and ${\omega }_D^2\sim {10}^{-24}$. For wavelengths ∼10 Kpc at recombination (and therefore within the horizon), k²c² ∼ 10⁻²².

Therefore, both $\frac{{\omega }_D^2}{k^2{c}^2}$ and $\frac{{\omega }_B^2}{k^2{c}^2}\ll 1$, and, to first order in these small quantities, the argument of the square root here becomes

$$\begin{array}{c}\hfill A+B=1+2(1-3\alpha )\frac{{\omega }_D^2}{\alpha k^2{c}^2}-2(1+\frac{15\alpha }{4})\frac{{\omega }_B^2}{\alpha k^2{c}^2},\end{array}$$

and the top line of the fraction becomes

$$\begin{array}{cc}& 1-(1+3\alpha )\frac{{\omega }_D^2}{\alpha k^2{c}^2}-(1+\frac{15\alpha }{4})\frac{{\omega }_B^2}{\alpha k^2{c}^2}\hfill \\ \hfill & \pm [1+(1-3\alpha )\frac{{\omega }_D^2}{\alpha k^2{c}^2}-(1+\frac{15\alpha }{4})\frac{{\omega }_B^2}{\alpha k^2{c}^2}]·\hfill \end{array}$$

For the minus sign, this reduces to $-2\frac{{\omega }_D^2}{\alpha k^2{c}^2}$, meaning the solution for ω² is

$$\begin{array}{c}\hfill {\omega }^2=-{\omega }_D^2.\end{array}$$

Therefore, the time dependence of all perturbations is exp(±ω_Dt) (i.e. there is one decaying mode and one unstable mode). For both these modes, Eqs. (16) and (18) imply

$$\begin{array}{c}\hfill (\frac{{\omega }_D^2}{k^2}+\alpha c^2){\overline{\rho }}_B=-(\frac{3{\omega }_D^2}{c^2{k}^2}+1+\frac{3}{4}\alpha ){\rho }_B\overline{\varphi }\end{array}$$

and

$$\begin{array}{c}\hfill {\overline{\rho }}_D=-(\frac{k^2}{{\omega }_D^2}+\frac{3}{c^2}){\rho }_D\overline{\varphi }.\end{array}$$

Therefore,

$$\begin{array}{c}\hfill \frac{{\overline{\rho }}_D}{{\overline{\rho }}_B}=\frac{1+\frac{\alpha k^2{c}^2}{{\omega }_D^2}+\frac{3{\omega }_D^2}{c^2{k}^2}+3\alpha }{1+\frac{3{\omega }_D^2}{c^2{k}^2}+\frac{3\alpha }{4}}\frac{{\rho }_D}{{\rho }_B}·\end{array}$$

The dominant term here is the second one on the top line, and hence

$$\begin{array}{c}\hfill \frac{{\overline{\rho }}_D}{{\overline{\rho }}_B}\approx \frac{\alpha k^2{c}^2}{{\omega }_D^2(1+\frac{3\alpha }{4})}\frac{{\rho }_D}{{\rho }_B}\gg 1,\end{array}$$

that is, the dark matter and baryonic matter density perturbations are in phase but with the dark matter perturbation very much greater than the baryonic perturbation. Although the baryonic matter reacts to the potential well created by the dark matter, radiation pressure prevents the baryonic matter from falling significantly into that well. Now from Eq. (17), and substituting for $\overline{\varphi }$ from Eq. (18), we obtain

$$\begin{array}{c}\hfill {\overline{v}}_D=\frac{k}{\omega (\frac{k^2}{{\omega }^2}-\frac{3}{c^2})}\frac{{\overline{\rho }}_D}{{\rho }_D},\end{array}$$

which for these two modes yields

$$\begin{array}{c}\hfill {\overline{v}}_D=\pm \frac{\mathit{ik}}{{\omega }_D(\frac{k^2}{{\omega }_D^2}+\frac{3}{c^2})}\frac{{\overline{\rho }}_D}{{\rho }_D},\end{array}$$

meaning the velocity perturbation is 90° out of phase with the density perturbation and, for the decaying mode (shown with the plus sign), involves expansion of the dark matter away from the density peak. On the other hand, the unstable mode (minus sign) involves gravitational collapse into the density peak. Any initial dark matter perturbation will potentially involve both of these modes; for instance, an initial density perturbation at rest involves equal amplitude contributions for both modes – one mode decays away, while the other mode grows exponentially in amplitude. The cold dark matter is therefore unconditionally unstable to gravitational collapse on all scales.

For the plus sign the solution (19) for ω² to first order in $\frac{{\omega }_D^2}{k^2{c}^2}$ and $\frac{{\omega }_B^2}{k^2{c}^2}$ is

$$\begin{array}{c}\hfill {\omega }^2=\alpha k^2{c}^2+(\frac{27}{4}{\alpha }^2-\frac{3}{4}\alpha -1){\omega }_B^2\approx \alpha k^2{c}^2.\end{array}$$

Again, for these modes, Eqs. (16) and (18) imply

$$\begin{array}{cc}& (\frac{27}{4}{\alpha }^2-\frac{3}{4}\alpha -1)\frac{{\omega }_B^2}{k^2}{\overline{\rho }}_B\hfill \\ \hfill & =(1-\frac{9}{4}\alpha -3(\frac{27}{4}{\alpha }^2-\frac{3}{4}\alpha -1)\frac{{\omega }_B^2}{k^2{c}^2}){\rho }_B\overline{\varphi }\hfill \end{array}$$

and

$$\begin{array}{c}\hfill {\overline{\rho }}_D=(\frac{k^2}{\alpha k^2{c}^2+(\frac{27}{4}{\alpha }^2-\frac{3}{4}\alpha -1){\omega }_B^2}-\frac{3}{c^2}){\rho }_D\overline{\varphi }.\end{array}$$

That is, to first order,

$$\begin{array}{c}\hfill \frac{{\overline{\rho }}_D}{{\overline{\rho }}_B}=\frac{(1-3\alpha )(\frac{27}{4}{\alpha }^2-\frac{3}{4}\alpha -1)}{\alpha (1-\frac{9}{4}\alpha )}\frac{{\omega }_B^2}{k^2{c}^2}\frac{{\rho }_D}{{\rho }_B}\ll 1\end{array}$$

at the era of recombination α ≈ 0.098, and hence

$$\begin{array}{c}\hfill \frac{{\overline{\rho }}_D}{{\overline{\rho }}_B}\approx 9.29\frac{{\omega }_B^2}{k^2{c}^2}\frac{{\rho }_D}{{\rho }_B},\end{array}$$

meaning ${\overline{\rho }}_D$ is in phase with ${\overline{\rho }}_B$ but is very small in comparison. This mode is the postulated baryon acoustic oscillations, with a phase velocity close to 0.3c due to being driven by radiation pressure rather than baryonic thermal pressure.

For the dark matter gravitational collapse mode, the assumption that $\omega \gg \frac{\mathrm{d}a/\mathrm{d}t}{a}$ turns out to be untenable, and hence neglecting the scale factor terms is not justifiable – prior to recombination there would only be a few e-folding times for the gravitational collapse to proceed. However, for the baryon acoustic oscillations, for wavelengths ∼10 Kpc at recombination (and therefore within the horizon),

$$\begin{array}{c}\hfill \omega \sim kc\gg \frac{\mathrm{d}a/\mathrm{d}t}{a},\end{array}$$

so the analysis holds for baryon acoustic oscillations.

4. Analysis including the scale factor terms for the dark matter collapse perturbations

We sought a solution of Eqs. (1)–(6) where the perturbations of ρ_B, v_B, and U_R are negligible, but we retained the scale factor terms in Eq. (4). We can still ignore the scale factor terms in Eq. (6) since $k^2{c}^2\gg \frac{\mathrm{d}a/\mathrm{d}t}{a}{\omega }_D$ and ${\left(\frac{\frac{\mathrm{d}a}{\mathrm{d}t}}{a}\right)}^2$. We can also drop the $\frac{\partial \varphi }{\partial t}$ term on the right-hand side of Eq. (3) since we expect the characteristic timescale to be ≫1/kc.

Equations (3), (4), and (6) for the Fourier component exp(−ikx) then become

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\overline{\rho }}_D}{\mathrm{d}t}& =ik{\rho }_D{\overline{v}}_D,\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{\overline{v}}_D}{\mathrm{d}t}+\frac{\mathrm{d}a/\mathrm{d}t}{a}{\overline{v}}_D& =ik\overline{\varphi ,}\hfill \end{array}$$(21)

and

$$\begin{array}{c}\hfill k^2\overline{\varphi }=-4\pi G{\overline{\rho }}_D.\end{array}$$(22)

Combining Eqs. (21) and (22) and multiplying both sides by a, we have

$$\begin{array}{c}\hfill \frac{\mathrm{d}\left(a{\overline{v}}_D\right)}{\mathrm{d}t}=-i\frac{4\pi Ga}{k}{\overline{\rho }}_D.\end{array}$$

This yields the set of equations

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\overline{\rho }}_D}{\mathrm{d}t}=\frac{k{\rho }_D}{a}{\overline{u}}_D\end{array}$$

and

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\overline{u}}_D}{\mathrm{d}t}=\frac{4\pi Ga}{k}{\overline{\rho }}_D,\end{array}$$

where $ia{\overline{v}}_D={\overline{u}}_D$.

Solving these numerically, starting with $\frac{{\overline{\rho }}_D}{{\rho }_D}={10}^{-6}$ and ${\overline{u}}_D=0$ at a  = 10⁻⁵, ${\overline{\rho }}_D$ increases almost linearly by a factor of just under 3 at a  = 10⁻³, corresponding to recombination. This solution has a potential field given by

$$\begin{array}{c}\hfill \overline{\varphi }=-\frac{4\pi G}{k^2}{\overline{\rho }}_D.\end{array}$$

This can be used to estimate the response of the baryons driven by this dark matter potential by equating the radiation pressure, which will oppose the collapse of the baryons to the gravitational pull of the potential well, yielding

$$\begin{array}{c}\hfill \frac{16\sigma T^3}{3c}\overline{T}=-{\rho }_B\overline{\varphi ,}\end{array}$$

meaning

$$\begin{array}{c}\hfill \frac{16\sigma T^3}{3c}\overline{T}=\frac{{\omega }_B^2}{k^2}{\overline{\rho }}_D.\end{array}$$(23)

Now, dropping the relativistic $\frac{\partial \varphi }{\partial t}$ term, Eq. (1) is

$$\begin{array}{c}\hfill \frac{\mathrm{d}{\rho }_B}{\mathrm{d}t}+{\rho }_B▿·v_B=0\end{array}$$(24)

and Eq. (5) is

$$\begin{array}{c}\hfill \frac{\mathrm{d}{U}_R}{\mathrm{d}t}+(U_R+p)▿·v_B=0.\end{array}$$(25)

Substituting for ▿ ⋅ v_B from Eqs. (24) to (25) yields

$$\begin{array}{c}\hfill \frac{1}{T}\frac{\mathrm{d}T}{\mathrm{d}t}=\frac{4}{3{\rho }_B}\frac{\mathrm{d}{\rho }_B}{\mathrm{d}t},\end{array}$$

and therefore

$$\begin{array}{c}\hfill \frac{\overline{T}}{T}=\frac{4}{3}\frac{{\overline{\rho }}_B}{{\rho }_B}·\end{array}$$

And from Eq. (23),

$$\begin{array}{c}\hfill \frac{48\sigma T^4}{9c}\frac{{\overline{\rho }}_B}{{\rho }_B}=\frac{{\omega }_B^2}{k^2}{\overline{\rho }}_D.\end{array}$$

Therefore,

$$\begin{array}{c}\hfill {\overline{\rho }}_B=\frac{1}{4\alpha }\frac{{\omega }_B^2}{k^2{c}^2}{\overline{\rho }}_D,\end{array}$$

that is,

$$\begin{array}{c}\hfill {\overline{\rho }}_B\ll {\overline{\rho }}_D.\end{array}$$

The baryons would be virtually unaffected by the collapse of the dark matter.

5. Conclusions

The eigenmode analysis here demonstrates the distinction between the baryon acoustic oscillation mode and the dark matter perturbation mode in the epoch between horizon crossing and recombination. Baryon acoustic waves are fast moving, with a high frequency oscillating potential to which the dark matter barely reacts; on the other hand, the dark matter perturbation mode is a static, relatively slow gravitational instability. The two modes are uncorrelated, and the distinction between them is ignored in modern analyses, where a constraint that couples the two modes is universally applied.

Naoz & Barkana (2005), for instance, couple the two wave modes together by applying a constraint that the baryon and dark matter perturbations have the same spatial structure (i.e. a Fourier component with the same wavenumber and phase). They then developed and numerically solved the resulting set of time-dependent differential equations for the perturbations. This dictates that the result will be a standing wave. This technique is used by many papers – another example, using the Einstein-Boltzmann equations for photons, is Pan et al. (2016), where their Eq. (8) explicitly shows the standing wave.

Figure 3 of Naoz & Barkana indicates that under this constraint the fractional perturbations of the baryons and dark matter tend to equality after recombination due to the baryons falling into the coincident dark matter potential wells. This justifies the initial conditions assumed by large-scale cosmological simulations.

However, the analysis here shows that the baryon element of a dark matter perturbation wave is very small, and not oscillatory. Hence, the standing wave in these calculations actually involves the coincidence of three separate waves – a dark matter perturbation wave and two baryon acoustic waves of equal magnitude travelling in opposite directions to create an oscillatory standing wave in phase with the dark matter wave. This situation would arise if the quantum fluctuations present in the inflationary epoch, which are surmised to be responsible for the perturbations considered here, consisted solely of gravitational potential perturbations. The effect of these gravitational perturbations would be to induce equal velocity perturbations in both the dark matter and baryons, leading to exactly this scenario. However, while the density perturbation in baryons and dark matter would be initially coincident, before recombination the density of the baryons would oscillate but that of the dark matter would not, so coincidence would occur only periodically. Authors often make a statement to the effect that the potential generated by a dark matter density perturbation drives the oscillations in the baryons (Li et al. 2008 and Eq. (8.18) of Dodelson 2003), but in fact an oscillatory baryon acoustic standing wave in phase with a dark matter perturbation would be a consequence of initial conditions set up by gravitational fluctuations during inflation.

If the quantum fluctuations during inflation actually consist of density perturbations, rather than pure gravitational perturbations, then there is no reason why perturbations in the dark matter and baryons would be related. Indeed, the differences highlighted here in the respective wave modes would make that highly unlikely.

On the assumption that all galaxies, and indeed galactic clusters, have their dark matter haloes, it could be argued that this justifies the assumption that dark matter and baryons should have coincident initial conditions for simulations. However, this would have difficulty accounting for the steeply falling rotation curve galaxies observed at z >  ∼1 (Lang et al. 2017; Genzel et al. 2017). These early galaxies seem to lack dark matter haloes. There are, however, very few, if any, steeply falling rotation curve galaxies in the local Universe – though recent stellar observations of a satellite galaxy of NGC 1052 indicate that it may contain no dark matter (van Dokkum et al. 2018). If steeply falling rotation curves turn out to be present in all early galaxies, and observations by the James Webb telescope may soon shed light on this, there would remain the requirement to explain how flat rotation curves develop as galaxies age (see for instance Nelson 1988). Initial conditions for simulations should at least consider the possibility of dark matter and baryons not being coincident at z ∼ 100. The possibility of propagating baryon acoustic oscillation (BAO) waves such as those that are deemed to create the ∼150 Mpc signature at z = 0 should also be considered. Tseliakhovich & Hirata (2010) propose that a relative velocity between dark matter and baryons due to the decrease in the sound velocity at recombination, going from radiation pressure to gas pressure, will have a significant effect on the subsequent development of structure. This effect is, however, second order in the perturbation amplitudes; the velocity and structural differences that would arise from the different evolution of dark matter perturbations and baryonic waves would be first order in the amplitudes.

Lastly, although the peaks in the power spectrum of the CMB anisotropies can be interpreted as standing waves in the finite length scale given by the sound horizon at the last scattering surface, for wavelengths much smaller than the sound horizon the evidence that there are standing waves is not so clear. Indeed, the idea that BAO waves propagate away from any initial inflation-induced perturbation is backed up by the BAO signature at z = 0 in the distribution of galaxies (Eisenstein et al. 2005). Purely standing waves at smaller scales are not consistent with this. This, and the analysis here, calls into question the assumption of standing baryonic waves in phase with dark matter in the analysis of the CMB anisotropies, as seen in Fig. 4 of Hu & Dodelson (2002), where fluctuations in the baryons are assumed to be driven by a potential created by perturbations in the dark matter. It could be argued that the success of this model in reproducing the peaks in the power spectrum justifies this assumption, but Lopez-Corredoira (2017) points out that matching the peaks can be achieved in a variety of ways not related to ΛCDM.

Acknowledgments

I thank Simon White for comments on the first draft of this paper.

References