If we plug the expressions for the radiation 4-momentum back into the gas equations and keep terms necessary to maintain accuracy we get:

\frac{\partial }{\partial t} \left(\rho\mathbf{v}\right)+\nabla\cdot\left(\rho\mathbf{vv}\right)=\nabla P\color{green}{-\lambda \nabla E}

\frac{\partial e}{\partial t}+\nabla\cdot\left[\left(e+P\right)\mathbf{v}\right]=\color{red}{-\kappa_{0P}(4 \pi B-cE)} \color{green}{+\lambda \left ( 2 \frac{\kappa_{0P}}{\kappa_{0R}}-1\right)\mathbf{v}\cdot\nabla E}\color{blue}{-\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E}

\frac{\partial E}{\partial t} \color{red}{ - \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E}=\color{red}{\kappa_{0P} (4 \pi B-cE)} \color{green}{-\lambda \left(2\frac{\kappa_{0P}}{\kappa_{0R}}-1\right)\mathbf{v}\cdot \nabla E} \color{green}{-\nabla \cdot \left ( \frac{3-R_2}{2}\mathbf{v}E\right )}\color{blue}{+\frac{3-R_2}{2}\kappa_{0P}\frac{v^2}{c}E}

Now if

E=a_RT^4

and

e=\rho c_v T

we can combine the gas and radiation diffusion equations to arrive at:

\frac{\partial \left (e+E\right) }{\partial t}+\nabla\cdot\left[\left(e+P\right)\mathbf{v}\right]=\color{red}{ \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E} \color{green}{-\nabla \cdot \left ( \frac{3-R_2}{2}\mathbf{v}E\right )}