This was spun off my "Einstein-Maxwell Conservation Law" worksheet, which started as this but became what it is.

Combining Maxwell and fluids

primitives:
$U^I = \left[\matrix{ \rho \\ \rho v^i \\ E^i \\ B^i \\ E_{hydro} }\right]$

conservatives:
$U^I = \downarrow I(i) \left[\matrix{ \rho \\ \rho v^i \\ E_{hydro} \\ \epsilon E^i \\ B^i }\right]$
$E_{hydro} = E_{kin} + E_{int}$
$= \frac{1}{2} \rho v^2 + \frac{P}{\gamma - 1} =$ total hydrodynamic energy
$E_{EM} = \frac{1}{2} ( \epsilon_0 E^2 + \frac{1}{\mu_0} B^2 ) = $ electomagnetic energy, but I'm not incorporating it into the conservative energy just yet...

Conservative PDEs:
$\rho_{,t} + (\rho v^i)_{,i} = 0$
$(\rho v^i)_{,t} + (\rho v^i v^j + \delta^{ij} P)_{,j} = 0$
$(E_{hydro})_{,t} + (v^j H_{hydro})_{,j} = 0$
$(\epsilon_0 E^i)_{,t} - (\frac{1}{\mu_0} {\epsilon^{ij}}_k B^k)_{,j} = -J^i$
${B^i}_{,t} + ({\epsilon^{ij}}_k E^k)_{,j} = 0$

Election vs Ion

using Abgrall, Kumar "Robust Finite Volume Schemes for Two-Fluid Plasma Equations

variables:
$q_\alpha = $ charge of species $\alpha$
$m_\alpha = $ mass of species $\alpha$
$r_\alpha = q_\alpha / m_\alpha$ = ratio of charge to mass
$\rho_\alpha = n_\alpha m_\alpha$
$P_\alpha = n_\alpha k_B T_\alpha$

Flux derivative wrt conservative variables:
$\left[\matrix{ \rho_{ion} \\ \rho_{ion} v_{ion}^i \\ E_{hydro,ion} \\ \rho_{elec} \\ \rho_{elec} v_{elec}^i \\ E_{hydro,elec} \\ \epsilon_0 E^i \\ B^i }\right]_{,t} + \left[\matrix{ 0 & \delta^{jk} & 0 & 0 & 0 & 0 & 0 & 0 \\ -v_{ion}^i v_{ion}^j + \frac{1}{2} \delta^{ij} \tilde\gamma v_{ion}^2 & \delta^i_k v_{ion}^j + \delta^j_k v_{ion}^i - \delta^{ij} \tilde\gamma v_{ion,k} & \delta^{ij} \tilde\gamma & 0 & 0 & 0 & 0 & 0 \\ v_{ion}^j (\frac{1}{2} \tilde\gamma v_{ion}^2 - h_{hydro,ion}) & -\tilde\gamma v_{ion}^j v_{ion,k} + \delta^j_k h_{hydro,ion} & \gamma v_{ion}^j & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 &0 & \delta^{jk} & 0 & 0 & 0 \\ 0 & 0 & 0 & -v_{elec}^i v_{elec}^j + \frac{1}{2} \delta^{ij} \tilde\gamma v_{elec}^2 & \delta^i_k v_{elec}^j + \delta^j_k v_{elec}^i - \delta^{ij} \tilde\gamma v_{elec,k} & \delta^{ij} \tilde\gamma & 0 & 0 \\ 0 & 0 & 0 & v_{elec}^j (\frac{1}{2} \tilde\gamma v_{elec}^2 - h_{hydro,elec}) & -\tilde\gamma v_{elec}^j v_{elec,k} + \delta^j_k h_{hydro,elec} & \gamma v_{elec}^j & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -\frac{1}{\mu_0} {\epsilon^{ij}}_k \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{\epsilon_0} {\epsilon^{ij}}_k & 0 \\ }\right] \left[\matrix{ \rho_{ion} \\ \rho_{ion} v_{ion}^k \\ E_{hydro,ion} \\ \rho_{elec} \\ \rho_{elec} v_{elec}^k \\ E_{hydro,elec} \\ \epsilon_0 E^k \\ B^k }\right]_{,j} = \left[\matrix{ 0 \\ r_{ion} \rho_{ion} (E^i + {\epsilon^i}_{jk} B^k v_{ion}^j) \\ r_{ion} \rho_{ion} (E_j v_{ion}^j) \\ 0 \\ r_{elec} \rho_{elec} (E^i + {\epsilon^i}_{jk} B^k v_{elec}^j) \\ r_{elec} \rho_{elec} (E_j v_{elec}^j) \\ - \Sigma_\alpha r_\alpha \rho_\alpha v^i_\alpha 0 \\ }\right]$

with constraints:
${\epsilon_0 E^i}_{,i} = \Sigma_\alpha r_\alpha \rho_\alpha$
${B^i}_{,i} = 0$