Solution (source code)

= Solution

Differentiate $E=\rho c^2a^3V_0$ and use $dE=-P\,dV$, with $V=a^3V_0$. Division by $c^2a^3V_0$ gives
$$
\boxed{\dot\rho=-3\frac{\dot a}{a}\left(\rho+\frac{P}{c^2}\right).}
$$
This is homogeneous <adiabatic expansion> without <entropy production> or exchange across the comoving element.