Suppose you have a momentum distribution of some decoupled $X$ particles in the early universe $f(\mathbf{p})$ that is injected in (well above the electroweak scale so that degrees of freedom for all species are relativistic). Assume the universe to be the flat FRW with the usual particle content.
This additional particle content contributes to the total energy density as $$\rho_X=\dfrac{g_{X}}{(2\pi)^3}\int E(\mathbf{p})~ f(\mathbf{p})~d^3\mathbf{p}.$$
In general, since $f(\mathbf{p})$ contains both relativistic and non-relativistic momenta (i.e. both $|\mathbf{p}|\ll m_{X}$ and $|\mathbf{p}|\gg m_{X}$) is seems oversimplified to just write down that either $\rho_{X}\propto a^{-3,-4}$ i.e., matter-like or radiation-like, because some of the distribution is relativistic, whilst some is not.
Is there a generalisation to describe the evolution of this energy density component at any given future time? or is there a reason why the energy density should always scale as exactly matter or radiation?
Edit: So I am thinking we want to figure out how it evolves with time, so of course this means considering the time derivative of this energy density. Using that $$\dfrac{d}{dt} \equiv \dfrac{\partial}{\partial t} + \dot{\mathbf{p}}\cdot \nabla_{\mathbf{p}}$$, and assuming there is no time dependence, then
$$\dfrac{d}{dt} \sim \dot{\mathbf{p}}\cdot \nabla_{\mathbf{p}}$$
If we make the assumption of homogeneity and isotropy, then the energy density only depends on the magnitude of the momentum too, so we can simplify further to
$$\dfrac{d}{dt}\sim \dfrac{d|\mathbf{p}|}{dt}\cdot \dfrac{\partial}{\partial |\mathbf{p}|}$$, and also
$$\rho_{X} \sim \dfrac{g_X}{2\pi^2}\int_{|\mathbf{p}|=0}^{\infty} |\mathbf{p}|^{2}~E(|\mathbf{p}|)~f(|\mathbf{p}|)~d|\mathbf{p}|$$
Therefore, the derivative of the energy density with respect to cosmic time should satisfy
$$\dfrac{d\rho_{X}}{dt}\sim \dfrac{g_{X}}{2\pi^2}~\dfrac{d|\mathbf{p}|}{dt}~\int 2|\mathbf{p}|E(|\mathbf{p}|)f(|\mathbf{p}|) + |\mathbf{p}|^2 \dfrac{|\mathbf{p}|}{E(|\mathbf{p}|)}f(|\mathbf{p}|) + |\mathbf{p}|^2 E(|\mathbf{p}|)~\dfrac{\partial f (|\mathbf{p}|)}{\partial |\mathbf{p}|}$$
I guess that energy conservation needs to somehow present itself here...
