Some questions about density perturbations in the early universe This is really a set of 3 questions in total, if that's fine. They are concerned with the theory of the evolution of small density perturbations in the early universe.
Question 1: Consider the following two equations (see source below):


Source: https://www.astro.rug.nl/~weygaert/tim1publication/lss2009/lss2009.linperturb.pdf
The first equation is the Poisson equation, and the second is the Euler equation. My first question is how these two equations are derived. In particular, what's the reason for the pressure term $3P/c^2$ in the first equation (the density term of course coming from Newton's gravitational law), and the term $P_{rad}/c^2$ in the denominator under the pressure gradient in the second equation?
Question 2: Combining the above equations with the continuity equation (not shown here), you may get a differential equation like the following, describing the evolution of the Fourier components of the (matter) density perturbations:

where k is the wavenumber of the Fourier component. Source: https://www.physicsforums.com/insights/poor-mans-cmb-primer-part-4-cosmic-acoustics/
Ignoring the friction term (the second term on the LHS), it's concluded in the article above that this equation either describes (acoustic) oscillations or gravitational collapse, depending on the value of k. What the article does not address, however, is that the factor $c_s^2k^2/a^2-4\pi G\rho$ in front of the density perturbation in the third term is dependent on time; the average matter density ($\rho$) is proportional to $a^{-3}$, and the other term contains $a^{-2}$. Since the scale factor depends on time, then how can the differential equation above describe harmonic oscillations (since we won't get a simple harmonic differential equation when this factor is time dependent)?
Question 3: Later on in the article above, it says that when you have a Fourier component with a wavelength equal to twice the sound horizon distance, then the perturbation will have got just the right amount of time at decoupling to fully compress. This is shown with the following diagram:

But doesn't this conclusion depend on the exact phase of the Fourier component? Why couldn't the perturbation have started fully compressed, only to reach maximum rarefaction at decoupling instead?
 A: Question 1
The pressure appears because General Relativity generalizes Newtonian gravity to include relativity. Instead of mass sourcing a gravitational field, in GR the source of gravity is the energy-momentum tensor $T_{\mu\nu}$. This is a matrix of numbers. One of the components of this tensor is the density, and in the non-relativistic limit, that component dominates Einstein's equations. However, other components include things like pressure and shear stresses, which will also contribute to the gravitational field if they are sufficiently strong.
Just to give an example of how expressions like the one you wrote can arise, the trace of the energy-momentum tensor is given by
\begin{equation}
T^\mu_{\ \ \mu} = -\rho + \frac{3P}{c^2}
\end{equation}
Note that the second term disappears in the non-relativistic limit, $c\rightarrow \infty$.
Question 2
The connection to the harmonic oscillator only applies when you look at solutions over a short enough time interval that $a$ is approximately constant. Note that you need to make this approximation to be able to ignore the Hubble friction term, which amounts to ignoring the time derivative of $a$.
In that approximation, the solution will oscillate for short wavelength modes with $k/a \gg 4 \pi G \bar{\rho}/c_s^2$, and will decay for long wavelength modes with $k/a \ll 4 \pi G \bar{\rho}/c_s^2$. Modes with $k \sim 4 \pi G \bar{\rho}/c_s^2$ experience a transition from the "decay" regime to the "oscillate" regime; to understand this in detail you need to do a more careful analysis, for example by using the WKB approximation to match the solutions in the two different regions, or numerically solving the equation.
Question 3
The equation has two solutions, a growing mode and a decaying mode. The phase you refer to amounts to an initial condition that tells you the relative strength of these two modes. In the current understanding of structure formation, the initial conditions are set by inflation (or some other process), which sets the growing mode to zero. In the past, people considered alternative scenarios, for example where initial perturbations where created by cosmic strings, which would excite both the growing and decaying mode, and lead to a different conclusion. However, these scenarios are now disfavored, because the phases of the peaks observed in the CMB agree with what you would expect with the initial conditions/phases generated naturally by inflation.
A: This is not a complete answer; it is a partial answer which justifies the location of the first acoustic peak but not the others.
Suppose that at early times (first minutes or so) there is a localized over-density in the dark matter (and of course most of the matter is dark). A spherical wave of over-density in the baryonic matter will propagate outwards from there, going at the speed of sound. It stops propagating at recombination when the pressure drops. So now you have a sphere of over-density, with radius equal to the sound horizon $s$, with an over-density at its centre. Many such spheres at random locations result in a distribution where over-densities are slightly more likely to be separated by $s$ than by some other distance.
That argument only applies to the first acoustic peak, but it gives a flavour of what is going on, and by bringing in the idea of oscillation one may be able to get some intuition for the other dips and peaks.
