Solution (source code)

= Solution

Dropping gravity and expansion from the <photon Boltzmann hierarchy> gives
$$
\boxed{\dot\Theta_0=-\frac{k}{3}\Theta_1},
\qquad
\boxed{\dot\Theta_1=k\Theta_0-\frac{2k}{5}\Theta_2+
\dot\tau(\Theta_1+v_b)}.
$$
The first equation is photon energy continuity: the <photon dipole> transports the <photon monopole>. The second is the <photon Euler equation>: a pressure gradient accelerates the radiation, its quadrupole contributes shear stress, and <Thomson scattering> drags the velocity toward the baryon velocity.

With $v_\gamma=-\Theta_1$, the scattering contribution to $\dot v_\gamma$ is $\dot\tau(v_\gamma-v_b)$. Radiation momentum density is $(4\bar\rho_\gamma/3)\mathbf v_\gamma$, whereas baryon momentum density is $\bar\rho_b\mathbf v_b$. Conservation of total momentum therefore requires
$$
\frac43\bar\rho_\gamma\,\dot v_\gamma\big|_{\rm sc}
+\bar\rho_b\,\dot v_b\big|_{\rm sc}=0,
$$
or
$$
\boxed{\dot v_b=-\frac{\dot\tau}{R}(v_\gamma-v_b),
\qquad R=\frac{3\bar\rho_b}{4\bar\rho_\gamma}}.
$$
Since the <cosmological optical depth> is measured to the observation time, $\dot\tau<0$; the two drag terms reduce their relative velocity.