Proper interpretation of $\rho_\textrm{rad} \, a^4$ in cosmology Radiation fluid is usually represented by the equation of state $p = \frac{1}{3} \, \rho$, where $p$ is the pressure and $\rho$ is the energy density.  The local conservation of energy states that
\begin{equation}\tag{1}
\rho_\textrm{rad} = \rho_0 \, \frac{a_0^4}{a^4},
\end{equation}
where $\rho_0$ and $a_0$ are constants and $a$ is the cosmological scale factor with units of length.  This is a well known subject in standard cosmology.  What puzzles me is the constant $\rho_0 \, a_0^4$.  What is its rigorous interpretation ?  For cold dust-like matter, it's easy :  $\rho_\textrm{mat} \propto a^{- 3}$ and $M = \rho_0 \, a_0^3$ is interpreted as the conserved total proper mass of matter inside the volume $a^3$.
In case of radiation, energy $E_\textrm{rad} = \rho \, a^3$ is not conserved, since there is pressure and the photons wavelength is redshifted while the universe expands.  But then what is the constant $\rho \, a^4$ ?  I strongly suspect it's related to the total radiation entropy in the volume $a^3$, or maybe the number of ultra-relativistic particles, but I can't find any reliable source on this.
I know that in standard RWFL (Robertson-Walker-Friedmann-Lemaître) models, entropy is conserved because of the local conservation of energy-momentum, which is equivalent of saying
\begin{equation}\tag{2}
\mathrm dE = T \, \mathrm dS - p \, \mathrm dV = -\, p \, \mathrm dV.
\end{equation}
I also know that Stefan-Boltzmann law states that $\rho_\textrm{rad} \propto T^4$, and classical radiation entropy is $S \propto T^3 V$ within a volume $V$.  Thus $\rho \, a^4 \propto T^4 \, V^{4/3} \propto S^{4/3}$, which is indeed conserved since entropy doesn't change.  But I'm not totally satisfied by this reasoning and never saw this in my books on General Relativity and cosmology.
Anyone has a convincing argument that this indeed should be true ?

EDIT :  The following is one aspect of the problem that puzzles me.  The thermodynamics argument above gives
\begin{equation}\tag{3}
\rho_\textrm{rad} \, a^4 \propto S^{4/3}.
\end{equation}
The exponent $\frac{4}{3}$ coincides with the adiabatic index $\gamma = \frac{4}{3}$ of the ultra-relativistic gas.  What is it doing here ?  If I introduce the entropy density $\sigma = S/V$, then we could apparently write this :
\begin{equation}\tag{4}
\rho_\textrm{rad} \propto \sigma^{\gamma},
\end{equation}
which recalls the polytropic gaz equation of state $p \propto \rho_\textrm{mass}^{\gamma}$.  Why is that ?  Can this be made more general (for $\gamma \ne \frac{4}{3}$) ?  I never saw the relations (3) and (4) in my thermodynamics books, as far as I can tell.  I need a more rigorous proof of these, and references if possible.
If $\rho \, a^4$ is related to the number of particles instead of entropy, then the polytropic gas state would make more sense, since $p \propto \rho_\textrm{rad} \propto n^\gamma$, where $n = N/V$ is the particles density.
 A: I think I got it!
I would like to have some confirmations that all the reasoning below is sound, and don't include some mistakes that I can't see.
In general, the macroscopic state of a perfect fluid can be described by a polytropic relation like this :
    \begin{equation}\tag{1}
  p = n \, k_B \, T = \kappa \, n^{\gamma},
 \end{equation}
where $n = N / V$ is the particles density, $k_B$ is the Boltzmann constant, $T$ is the fluid's temperature and $\gamma$ is the adiabatic index of the perfect fluid (note that the temperature is usually a function of the particles density ; $T(n) \propto n^{\gamma \,-\, 1}$, which physically makes sense when $\gamma \ge 1$).  The proportionality constant $\kappa$ is arbitrary and can be ignored for the rest of this answer.  The usual state relation
\begin{equation}
p = w \, \rho
\end{equation}
then let us express the energy density $\rho$ in terms of the particles density $n$.  If the state parameter $w$ is a constant ($w \ne 0$), then
    \begin{equation}\tag{2}
  \rho \propto n^{\gamma}
 \end{equation}
The fluid's energy inside a volume $V$ is then
    \begin{equation}\tag{3}
  E = \rho \, V \propto N^{\gamma} \, V^{1 \,-\, \gamma}.
 \end{equation}
Its differential gives $dE = -\: p \; dV$ if the mean number $N$ is fixed.  It also gives
    \begin{equation}\tag{4}
  w = \gamma - 1,
 \end{equation}
which would be usefull below.  Because of the first law of thermodynamics ; $dE = T \; dS - p \; dV$, this implies that the entropy is conserved.  Because of equation (2), we can write the following relation :
    \begin{equation}\tag{5}
  \rho \, V^{\gamma} \propto N^{\gamma}.
 \end{equation}
In the special case of a perfect fluid made of non-interacting dust-like particles, the temperature is negligible ; $p = n \, k_B \, T \approx 0$, which is a good approximation if the particles density is very low.  The total proper mass of the dust inside a volume $V \propto a^3$ is $M = \rho \, V \propto N$.  This suggest to define the adiabatic index $\gamma = 1$ for the dust-like matter (of course, we have $w \equiv 0$ since pressure is negligible).  Thus, relation (5) stay good even for the fluid made of "dust" ($w = 0$).
In the case of ultra-relativistic particles (radiation), we have $w = \frac{1}{3}$ and equation (4) gives $\gamma = \frac{4}{3}$, so the relation (5) gives our final result, using $V \propto a^3$ :
    \begin{equation}\tag{6}
  \rho \, a^4 \propto N^{\frac{4}{3}}.
 \end{equation}
A: It's actually much simpler than that. The number of photons is a conserved quantity in vacuum when gravitation is weak enough to avoid the Unruh effect, so the photon number density scales the same way the matter number density does:
$$n = n_0 \frac{a_0^3}{a^3}.$$
To get the total energy density of radiation you take the spectral number density (density of photons per unit space volume per unit frequency, $n_\nu$) and calculate:
$$\rho_{\mathrm{rad}} = \int_0^\infty h\nu n_\nu \operatorname{d} \nu.$$
Note that because $n_\nu \operatorname{d}\nu$, a differential form, is an ordinary number density it will scale the same way as any other number density, taken as a unit, as long as we keep track of the limits of the integral correctly. Since they're $0$ and $\infty$, though, that's not a problem. The result is:
$$\begin{align}\rho_{\mathrm{rad}} & = \int_0^\infty h \left(\nu_0 \frac{a_0}{a}\right) \left( n_{\nu_0} \operatorname{d} \nu_0 \frac{a_0^3}{a^3} \right) \\
& = \rho_0 \frac{a_0^4}{a^4} .\end{align}$$
Put another way, when you spread out photons by scaling up the volume by a factor of $a^3$ and increase their wavelenth by $a$, their energy density drops by $a^4$.
Note that the photon spectrum in this derivation can take any form, so it is more general than an entropy based derivation that assumes a gas in thermodynamic equilibrium at some temperature.
Edit to add:

What puzzles me is the constant $ρ_0a_0^4$. What is its rigorous interpretation ?

The specific interpretation of it is that it's the radiation density when the scale factor, $a$, is 1. Usually $a$ is taken as the present day, by convention, but any other convention is allowable. The likely reason for this convention is that the present radiation density is the quantity we can most easily measure directly when we measure the temperature of the cosmic microwave background.
