How to justifiy that $\rho_{\text{rad}} \approx \rho_{\text{mat}}$ at recombination time? In standard cosmology, the recombination time is estimated to be $t_{\text{rec}} \approx 380~000~\mathrm{years}$ after the Big Bang, when matter and electromagnetic radiation becomes decoupled.  It's the time the CMB was created (or "liberated" from matter).
Since $\rho = \rho_{\text{mat}} + \rho_{\text{rad}}$ and $p = \frac{1}{3} \, \rho_{\text{rad}}$, the sound velocity in the cosmic fluid is (I'm using $c = 1$ here):
\begin{equation}\tag{1}
v_{\text{sound}} = \sqrt{\frac{dp}{d\rho}} = \frac{1}{\sqrt{3}} \, \frac{1}{\sqrt{1 + \frac{3 \rho_{\text{mat}}}{4 \rho_{\text{rad}}}}}.
\end{equation}
At time $t_{\text{rec}}$, it is often said (written) that $\frac{3 \rho_{\text{mat}}}{4 \rho_{\text{rad}}} \approx 1$ (so the sound velocity is $v_{\text{sound}} \approx 1/\sqrt{6}$ at recombination time).  How can we justify this?
In other words, how do we calculate (or predict) the recombination time?  I suspect that it's from the Saha equation (which I don't know much), and not from the energy density relation which feels a bit arbitrary:
\begin{equation}\tag{2}
3 \rho_{\text{mat}} \approx 4 \rho_{\text{rad}}.
\end{equation}
 A: The matter density scales as $a^{-3}$ whilst the radiation density scales as $a^{-4}$.
The parameters $\rho_{\rm rad}$ and $\rho_{\rm mat}$ in your equation correspond to the densities in the photon-baryon fluid. i.e. They are the densities of only the interacting photons and baryons and do not include neutrinos or dark matter. At the present epoch $\rho_{{\rm rad},0} \sim 5\times 10^{-5}$ of the critical density, whereas $\rho_{{\rm mat},0} \sim 0.047$ of the critical density.
At the epoch of recombination $z_{\rm rec} \sim 1100$. This does indeed come from the Saha equation and is a standard piece of bookwork which calculates the redshift at which some fraction of the hydrogen becomes atomic, though getting the exact answer is a bit more complex than that (see for example here).
The scale factor $a$ changes as $(1+z)^{-1}$, so, if we go back to the epoch of recombination
$$ \frac{\rho_{\rm mat}}{\rho_{\rm rad}} = \frac{\rho_{{\rm mat},0}}{ \rho_{{\rm rad},0}} \left(1+z_{\rm rec}\right)^{-1} \simeq 0.85 \ .$$
In terms of how you work out the epoch of recombination; in brief, the process is that the Saha equation, combined with the baryon to photon ratio, tells us at what temperature there are insufficient photons in the tail of the Planck function with $E>13.6$ eV to keep the (relatively small) number of electrons from recombining with protons. Once the temperature falls below this, recombination begins and we can get an expression for the ionisation fraction as a function of $T$ or $z$.
From the Saha equation:
$$\frac{n_e^2}{n_H} = \left(\frac{m_e k_b T}{2\pi \hbar^2} \right)^{3/2} \exp\left(\frac{-13.6{\rm eV}}{k_bT}\right)\ , $$
where $n_H$ is the number density of hydrogen atoms. If we then let $x$ be the ratio of free electrons to baryons, we can say $n_H + n_e \simeq n_b$ and $n_e = xn_b$(there is a correction for helium of course) and thus
$$\frac{x^2}{1-x}= \left(\frac{m_e k_b T}{2\pi \hbar^2 n_b} \right)^{3/2} \exp\left(\frac{-13.6{\rm eV}}{k_bT}\right)\ .$$
$n_b$ is estimated from the current value of $\Omega_b$ and critical density, scaled by $a^{-3}= (1+z)^3$ and then a plot of $x$ vs $z$ can be produced. The redshift of recombination depends on the adopted value of $x$ which varies from 1 to almost zero, between $1500 > z > 900$.
What is commonly termed the epoch of recombination is actually the redshift of decoupling. One calculates the rate of photon interactions with matter (mainly through Thomson scattering) as the speed of light divided by the mean free path. The rate falls rapidly when recombination begins, since the mean free path is inversely proportional to electron density. When the interaction rate falls below the Hubble parameter, then the radiation and matter are decoupled and this occurs at $z\sim 1100$.
