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

or

\left(1+\frac{e}{4E}\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{E+\frac{e}{4}}{\nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E}

Changes to the discretization

For the coupled system of equations we had the following:

\color{purple}{\left [ 1 + \psi \left( \alpha_{i+1/2} + \alpha_{i-1/2} + \frac{\epsilon_i}{ 1 +\psi \phi_i}\right ) \right ] E^{n+1}_i - \left ( \psi \alpha_{i+1/2} \right ) E^{n+1}_{i+1} - \left ( \psi \alpha_{i-1/2} \right ) E^{n+1}_{i-1} =\left [ 1 - \bar{\psi} \left( \alpha_{i+1/2} + \alpha_{i-1/2} +\frac{\epsilon_i }{ 1 +\psi \phi_i} \right ) \right ] E^n_i + \left ( \bar{\psi} \alpha_{i+1/2} \right ) E^{n}_{i+1} + \left ( \bar{\psi} \alpha_{i-1/2} \right ) E^{n}_{i-1} + \frac{\theta_i}{ 1 +\psi \phi_i}}

If the gas and radiation are in thermal equilibrium, then we have \theta_i = \epsilon_i E_i and we also have that in the limit that \kappa_P \rightarrow \infty, we have \epsilon \rightarrow \infty and \phi \rightarrow \infty

This simplifies the above equation to

\color{purple}{\left [ \left(1 + \frac{\epsilon_i}{\phi_i} \right) + \psi \left( \alpha_{i+1/2} + \alpha_{i-1/2} \right ) \right ] E^{n+1}_i - \left ( \psi \alpha_{i+1/2} \right ) E^{n+1}_{i+1} - \left ( \psi \alpha_{i-1/2} \right ) E^{n+1}_{i-1} =\left [ \left(1 + \frac{\epsilon_i}{\phi_i}\right) - \bar{\psi} \left( \alpha_{i+1/2} + \alpha_{i-1/2} \right ) \right ] E^n_i + \left ( \bar{\psi} \alpha_{i+1/2} \right ) E^{n}_{i+1} + \left ( \bar{\psi} \alpha_{i-1/2} \right ) E^{n}_{i-1} }

And \frac{\epsilon_i}{\phi_i} = \frac{\frac{T_i}{4\Gamma}}{E_i} = \frac{\frac{T_i}{4\frac{\partial T}{\partial e}}}{E_i} which if we use our equation of state where e \propto T gives \frac{\epsilon_i}{\phi_i} = \frac{e_i}{4E_i}

Now if we go back and calculate \frac{\partial e}{\partial E} = \frac{\partial e}{\partial T}\frac{\partial T}{\partial E} = \frac{\frac{1}{\Gamma}}{\frac{4 E}{T}} we arrive at \frac{\epsilon_i}{\phi_i} instead of \frac{e_i}{4E_i} which is consistent with our derivation above.

Also our time equation should be

\Delta t \approx \frac{E}{\frac{\partial E}{\partial t}} = \frac{E+\frac{T_i}{4\Gamma}}{\nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E}

So in principle combining the gas and energy equations in the limit that of high planck opacity, does not change the matrix or rhs vector, however it does limit the ability for there to be strong source terms on the right in regions where the gas and radiation have gotten out of equilibrium. It is not clear how this effects the ability of the elliptic solver to converge to given tolerances.

Modifications to time steps

More importantly is the recognition of the time scales over which the internal energy can change.

Previously we looked at the decoupled equation for the gas energy density

\frac{\partial e}{\partial t}=-\kappa_{0P}(4 \pi B-cE)

\Delta e = \Delta t \kappa_{0P} \left | 4 \pi B_0 -cE \right | < \xi \frac{T_0}{4 \Gamma}

which gives \Delta t < \xi \frac{T_0}{4 \Gamma \kappa_{0P} \left | 4 \pi B_0 - cE \right | }

however, if the gas is in equilibrium with the radiation, this does not limit the time step at all - even though diffusion may quickly move the gas out of equilibrium with the radiation.

We can account for this by expanding our equation

\Delta t < \xi \frac{T_0}{4 \Gamma \kappa_{0P} \left | 4 \pi B_0 - c\left(E_0+\partial_tE\Delta t\right ) \right | } which gives us a quadratic for \Delta t

-4 \Gamma \kappa_{0P} c \partial_tE \Delta t^2 + \left(4 \Gamma \kappa_{0P} \left ( 4 \pi B_0 - cE_0 \right ) \right ) \Delta t - \xi T_0 < 0

4 \Gamma \kappa_{0P} c \partial_tE \Delta t^2 - \left(4 \Gamma \kappa_{0P} \left ( 4 \pi B_0 - cE_0 \right ) \right ) \Delta t - \xi T_0 < 0

yucky…

How about we assume that the gas is close to equilibrium? Then we get

\Delta t < \xi \frac{T_0}{4 \Gamma \kappa_{0P} \left | \partial_tE\Delta t \right | }

\Delta t < \sqrt{\xi \frac{T_0}{4 \Gamma \kappa_{0P} c \left | \nabla \cdot \frac{c\lambda}{\kappa_{0R}} \nabla E \right |} } which is just the geometric average of the coupling time and the diffusion time