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=a_R\left(\frac{e}{\rho c_v}\right)^4

and

e=\rho c_v T=\rho c_v \left(\frac{E}{a_R}\right)^{1/4}

and we just consider the implicit terms, we can combine the gas and radiation diffusion equations to arrive at:

\frac{\partial \left (e+E\right) }{\partial t}=\color{red}{ \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E}

which simplifies to

\left(1+\frac{\rho c_v}{4 E}\left(\frac{E}{a_R}\right )^{1/4}\right) \frac{\partial E}{\partial t}= \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E

The second term in parenthesis represents the extra 'inertia' the radiation field has due to its coupling with the gas. It is non-linear and this limits the time step that can be taken.

\Delta t \approx \frac{E}{\frac{\partial E}{\partial t}} = \frac{\left(E+\frac{\rho c_v}{4}\left(\frac{E}{a_R}\right )^{1/4}\right)}{\nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E}

Changes to the discretization