two-fluid plasma model coupled with GEM
Going off of 2014 Abgrall, Kumar
and 2003 Loverich
with some extras from Wikipedia

constants:
$c =$ $\frac{m}{s} =$ speed of light.
$k_e = $ $\cdot \frac{kg \cdot m^3}{C^2 \cdot s^2} =$ Coulomb's constant (typically $\frac{1}{4 \pi \epsilon_0}$).
$k_B = $ $\cdot \frac{kg \cdot m^2}{K \cdot s^2} =$ Boltzmann's constant.
$\epsilon_0 = \frac{1}{4 \pi k_e} = $ {{vacuum_permittivity_in_C2_s2_per_kg_m3 = 1 / (4 * Math.PI * Coulomb_constant_in_kg_m3_per_C2_s2)}} $\cdot \frac{C^2 \cdot s^2}{kg \cdot m^3}$
$\mu_0 = \frac{1}{c^2 \cdot \epsilon_0} = $ {{vacuum_permeability_in_kg_m_per_C2 = 1 / (speed_of_light_in_m_per_s * speed_of_light_in_m_per_s * vacuum_permittivity_in_C2_s2_per_kg_m3)}} $\cdot \frac{kg \cdot m}{C^2} =$ vacuum permeability
such that $\frac{1}{\mu_0 \epsilon_0} = c^2 =$ {{speed_of_light_in_m_per_s * speed_of_light_in_m_per_s}} $\frac{m^2}{s^2}$

material properties:
$\epsilon =$ material permittivity, in $\frac{C^2 \cdot s^2}{kg \cdot m^3}$
$\mu =$ material permeability, in $\frac{kg \cdot m}{C^2}$
$v_p = \frac{1}{\sqrt{\mu \epsilon} } =$ phase velocity through medium

electromagnetic properties:
$E_i =$ electric field, in $\frac{kg \cdot m}{C \cdot s^2}$
$D_i = \epsilon E =$ displacement field, in $\frac{C}{m^2}$
$B_i =$ magnetic field, in $\frac{kg}{C \cdot s}$
$H_i = \frac{1}{\mu} B = $ magnetizing field, in $\frac{C}{m \cdot s}$

$\rho_{charge} = (\hat{\lambda}_D^2 \hat{r}_L)^{-1} (r_{ion} \rho_{ion} + r_{elec} \rho_{elec}) =$ charge density, in $\frac{C}{m^3}$
$J_i = ( \hat{\lambda}_D^2 \hat{r}_L )^{-1} (r_{ion} \rho_{ion} v_{ion}^i + r_{elec} \rho_{elec} v_{elec}^i) =$ current density, in $\frac{C}{m^2 \cdot s}$
(TODO why the extra Larmor / Debye variables? Do the units of $( \hat{\lambda}_D^2 \hat{r}_L )^{-1}$ cancel?)

GLM variables:
$\phi =$ electric divergence potential, in $\frac{C}{m^2}$
$\xi =$ speed of propagation of electric divergence, in $\frac{m}{s}$
$\psi =$ magnetic divergence potential, in $\frac{kg}{C \cdot s}$
$\kappa =$ speed of propagation of magnetic divergence, in $\frac{m}{s}$

GLM Maxwell equations:
$[\frac{C}{m^2 \cdot s}]: \partial_t D_i - \mu^{-1} {\epsilon_i}^{jk} \nabla_j B_k + \xi \nabla_i \phi = -J_i$
$[\frac{kg}{C \cdot s^2}]: \partial_t B_i + \epsilon^{-1} {\epsilon_i}^{jk} \nabla_j D_k + \kappa \nabla_i \psi = 0$
$[\frac{C}{m^2 \cdot s}]: \partial_t \phi + \xi \nabla_j D^j = \xi \rho_{charge}$
$[\frac{kg}{C \cdot s^2}]: \partial_t \psi + \kappa \nabla_j B^j = 0$

As a linear system:
$\partial_t \left[\begin{matrix} D_i \\ B_i \\ \phi \\ \psi \end{matrix}\right] = \left[\begin{matrix} 0 & -\frac{1}{\mu} {\epsilon_i}^{jk} & \xi & 0 \\ \frac{1}{\epsilon} {\epsilon_i}^{jk} & 0 & 0 & \kappa \\ \delta^{jk} \xi & 0 & 0 & 0 \\ 0 & \delta^{jk} \kappa & 0 & 0 \end{matrix}\right] \nabla_j \left[\begin{matrix} D_k \\ B_k \\ \phi \\ \psi \end{matrix}\right] = \left[\begin{matrix} -J_i \\ 0 \\ \xi \rho_{charge} \\ 0 \end{matrix}\right]$

Decomposed into individual matrices:
Let $U = \left[\begin{matrix} D_i \\ B_i \\ \phi \\ \psi \end{matrix}\right], S = \left[\begin{matrix} -J_i \\ 0 \\ \xi \rho_{charge} \\ 0 \end{matrix}\right]$

the system becomes:
$\partial_t U + A^j \nabla_j U = S$

for flux jacobians:
$A^x = \left[\begin{matrix} 0 & 0 & 0 & 0 & 0 & 0 & {\chi_\phi} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{\mu} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{-1}{\mu} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & {\chi_\psi} \\ 0 & 0 & \frac{-1}{\epsilon} & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{\epsilon} & 0 & 0 & 0 & 0 & 0 & 0 \\ {\chi_\phi} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & {\chi_\psi} & 0 & 0 & 0 & 0 \end{matrix}\right]$, $A^y = \left[\begin{matrix} 0 & 0 & 0 & 0 & 0 & \frac{-1}{\mu} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & {\chi_\phi} & 0 \\ 0 & 0 & 0 & \frac{1}{\mu} & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{\epsilon} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & {\chi_\psi} \\ \frac{-1}{\epsilon} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & {\chi_\phi} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & {\chi_\psi} & 0 & 0 & 0 \end{matrix}\right]$, $A^z = \left[\begin{matrix} 0 & 0 & 0 & 0 & \frac{1}{\mu} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{-1}{\mu} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & {\chi_\phi} & 0 \\ 0 & \frac{-1}{\epsilon} & 0 & 0 & 0 & 0 & 0 & 0 \\ \frac{1}{\epsilon} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & {\chi_\psi} \\ 0 & 0 & {\chi_\phi} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & {\chi_\psi} & 0 & 0 \end{matrix}\right]$

with flux jacobian eigensystems:
$A^x = \left[\begin{matrix} -{\frac{1}{\sqrt{2} \sqrt{\epsilon}}} & 0 & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\epsilon}} & 0 \\ 0 & 0 & 0 & \frac{-{\sqrt{\epsilon}}}{\sqrt{2}} & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} & 0 & 0 \\ 0 & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} & 0 & \frac{-{\sqrt{\epsilon}}}{\sqrt{2}} & 0 & 0 & 0 \\ 0 & -{\frac{1}{\sqrt{2} \sqrt{\epsilon}}} & 0 & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\epsilon}} \\ 0 & 0 & \frac{\sqrt\mu}{\sqrt{2}} & 0 & \frac{\sqrt\mu}{\sqrt{2}} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{\sqrt\mu}{\sqrt{2}} & 0 & \frac{\sqrt\mu}{\sqrt{2}} & 0 & 0 \\ \frac{1}{\sqrt{2} \sqrt{\epsilon}} & 0 & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\epsilon}} & 0 \\ 0 & \frac{1}{\sqrt{2} \sqrt{\epsilon}} & 0 & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\epsilon}} \end{matrix}\right] \cdot \left[\begin{matrix} -\chi_\phi & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -\chi_\psi & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{-1}{\sqrt{\mu \cdot \epsilon}} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{-1}{\sqrt{\mu \cdot \epsilon}} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{\sqrt{\mu \cdot \epsilon}} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{\mu \cdot \epsilon}} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & {\chi_\phi} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & {\chi_\psi} \end{matrix}\right] \cdot \left[\begin{matrix} \frac{-{\sqrt{\epsilon}}}{\sqrt{2}} & 0 & 0 & 0 & 0 & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} & 0 \\ 0 & 0 & 0 & \frac{-{\sqrt{\epsilon}}}{\sqrt{2}} & 0 & 0 & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} \\ 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\epsilon}} & 0 & \frac{1}{\sqrt{2} \sqrt{\mu}} & 0 & 0 & 0 \\ 0 & \frac{1}{-{\sqrt{2} \sqrt{\epsilon}}} & 0 & 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\mu}} & 0 & 0 \\ 0 & 0 & -{\frac{1}{\sqrt{\epsilon} \sqrt{2}}} & 0 & \frac{1}{\sqrt{2} \sqrt{\mu}} & 0 & 0 & 0 \\ 0 & \frac{1}{\sqrt{\epsilon} \sqrt{2}} & 0 & 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\mu}} & 0 & 0 \\ \frac{\sqrt{\epsilon}}{\sqrt{2}} & 0 & 0 & 0 & 0 & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} & 0 \\ 0 & 0 & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} & 0 & 0 & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} \end{matrix}\right]$
$A^y = \left[\begin{matrix} 0 & 0 & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} & 0 & \frac{-{\sqrt{\epsilon}}}{\sqrt{2}} & 0 & 0 \\ -{\frac{1}{\sqrt{2} \sqrt{\epsilon}}} & 0 & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\epsilon}} & 0 \\ 0 & 0 & \frac{-{\sqrt{\epsilon}}}{\sqrt{2}} & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} & 0 & 0 & 0 \\ 0 & 0 & \frac{\sqrt\mu}{\sqrt{2}} & 0 & \frac{\sqrt\mu}{\sqrt{2}} & 0 & 0 & 0 \\ 0 & -{\frac{1}{\sqrt{2} \sqrt{\epsilon}}} & 0 & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\epsilon}} \\ 0 & 0 & 0 & \frac{\sqrt\mu}{\sqrt{2}} & 0 & \frac{\sqrt\mu}{\sqrt{2}} & 0 & 0 \\ \frac{1}{\sqrt{2} \sqrt{\epsilon}} & 0 & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\epsilon}} & 0 \\ 0 & \frac{1}{\sqrt{2} \sqrt{\epsilon}} & 0 & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\epsilon}} \end{matrix}\right] \cdot \left[\begin{matrix} -\chi_\phi & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -\chi_\psi & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{-1}{\sqrt{\mu \cdot \epsilon}} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{-1}{\sqrt{\mu \cdot \epsilon}} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{\sqrt{\mu \cdot \epsilon}} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{\mu \cdot \epsilon}} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & {\chi_\phi} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & {\chi_\psi} \end{matrix}\right] \cdot \left[\begin{matrix} 0 & \frac{-{\sqrt{\epsilon}}}{\sqrt{2}} & 0 & 0 & 0 & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} & 0 \\ 0 & 0 & 0 & 0 & \frac{-{\sqrt{\epsilon}}}{\sqrt{2}} & 0 & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} \\ 0 & 0 & \frac{1}{-{\sqrt{2} \sqrt{\epsilon}}} & \frac{1}{\sqrt{2} \sqrt{\mu}} & 0 & 0 & 0 & 0 \\ \frac{1}{\sqrt{2} \sqrt{\epsilon}} & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\mu}} & 0 & 0 \\ 0 & 0 & \frac{1}{\sqrt{\epsilon} \sqrt{2}} & \frac{1}{\sqrt{2} \sqrt{\mu}} & 0 & 0 & 0 & 0 \\ -{\frac{1}{\sqrt{\epsilon} \sqrt{2}}} & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\mu}} & 0 & 0 \\ 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} & 0 & 0 & 0 & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} & 0 \\ 0 & 0 & 0 & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} & 0 & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} \end{matrix}\right]$
$A^z = \left[\begin{matrix} 0 & 0 & 0 & \frac{-{\sqrt{\epsilon}}}{\sqrt{2}} & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} & 0 & 0 \\ 0 & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} & 0 & \frac{-{\sqrt{\epsilon}}}{\sqrt{2}} & 0 & 0 & 0 \\ -{\frac{1}{\sqrt{2} \sqrt{\epsilon}}} & 0 & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\epsilon}} & 0 \\ 0 & 0 & \frac{\sqrt\mu}{\sqrt{2}} & 0 & \frac{\sqrt\mu}{\sqrt{2}} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{\sqrt\mu}{\sqrt{2}} & 0 & \frac{\sqrt\mu}{\sqrt{2}} & 0 & 0 \\ 0 & -{\frac{1}{\sqrt{2} \sqrt{\epsilon}}} & 0 & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\epsilon}} \\ \frac{1}{\sqrt{2} \sqrt{\epsilon}} & 0 & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\epsilon}} & 0 \\ 0 & \frac{1}{\sqrt{2} \sqrt{\epsilon}} & 0 & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\epsilon}} \end{matrix}\right] \cdot \left[\begin{matrix} -\chi_\phi & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -\chi_\psi & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{-1}{\sqrt{\mu \cdot \epsilon}} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{-1}{\sqrt{\mu \cdot \epsilon}} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{\sqrt{\mu \cdot \epsilon}} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{\sqrt{\mu \cdot \epsilon}} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & {\chi_\phi} & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & {\chi_\psi} \end{matrix}\right] \cdot \left[\begin{matrix} 0 & 0 & \frac{-{\sqrt{\epsilon}}}{\sqrt{2}} & 0 & 0 & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{-{\sqrt{\epsilon}}}{\sqrt{2}} & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} \\ 0 & \frac{1}{\sqrt{2} \sqrt{\epsilon}} & 0 & \frac{1}{\sqrt{2} \sqrt{\mu}} & 0 & 0 & 0 & 0 \\ \frac{1}{-{\sqrt{2} \sqrt{\epsilon}}} & 0 & 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\mu}} & 0 & 0 & 0 \\ 0 & \frac{1}{-{\sqrt{2} \sqrt{\epsilon}}} & 0 & \frac{1}{\sqrt{\mu} \sqrt{2}} & 0 & 0 & 0 & 0 \\ \frac{1}{\sqrt{\epsilon} \sqrt{2}} & 0 & 0 & 0 & \frac{1}{\sqrt{2} \sqrt{\mu}} & 0 & 0 & 0 \\ 0 & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} & 0 & 0 & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} & 0 & \frac{\sqrt{\epsilon}}{\sqrt{2}} \end{matrix}\right]$

---

GEM constants:
$G =$ $\frac{m^3}{kg \cdot s^2} =$ gravitational constant.
$\epsilon_g = \frac{1}{4 \pi G} = $ {{GEM_permittivity_in_kg_s2_per_m3 = 1 / (4 * Math.PI * gravitational_constant_in_m3_per_kg_s2)}} $\cdot \frac{kg \cdot s^2}{m^3} =$ GEM vacuum permittivity.
$\mu_g = \frac{1}{c^2 \cdot \epsilon_g} = $ {{GEM_permeability_in_m_per_kg = 1 / (speed_of_light_in_m_per_s * speed_of_light_in_m_per_s * GEM_permittivity_in_kg_s2_per_m3)}} $\cdot \frac{m}{kg} =$ GEM vacuum permeability.
such that $\frac{1}{\mu_g \epsilon_g} = c^2$

GEM properties:
$(E_g)^i =$ gravito-electric field, in $\frac{m}{s^2}$
$(D_g)^i = \epsilon_g E_g =$ gravito-displacement field, in $\frac{kg}{m^2}$
$(B_g)^i =$ gravito-magnetic field, in $\frac{1}{s}$
$(H_g)^i = \frac{1}{\mu_g} B_g =$ gravito-magnetic field, in $\frac{kg}{m \cdot s}$

stress-energy terms:
$T_{00} = c^2 (\rho_{ion} + \rho_{elec}) + \frac{1}{2} (\frac{1}{\epsilon} D^2 + \frac{1}{\mu} B^2) = $ energy from stress-energy
$T_{0i} = -c (\rho_{ion} v_{ion,i} + \rho_{elec} v_{elec,i}) + \frac{1}{2 c} (n^2 + 1) S_i = $ momentum from stress-energy
(does gravitation self-interact? should we also consider the energy of $\frac{1}{2} (\frac{1}{\epsilon_g} (D_g)^2 + \frac{1}{\mu_g} (B_g)^2)$? How significant would that be?)

$(\rho_g) = \frac{1}{c^2} T_{00} = \rho_{ion} + \rho_{elec} + \frac{1}{2 c^2} (\frac{1}{\epsilon} D^2 + \frac{1}{\mu} B^2) =$ matter density, in $\frac{kg}{m^3}$
$(J_g)_i = \frac{1}{c} T_{0i} = -\rho_{ion} v_{ion,i} - \rho_{elec} v_{elec,i} + \frac{1}{2 c^2} (n^2 + 1) S_i$ matter momentum, in $\frac{kg}{m^2 \cdot s}$

GLM-GEM variables:
$\phi_g =$ gravito-electric divergence potential, in $\frac{kg}{m^2}$
$\xi_g =$ speed of propagation of gravito-electric divergence, in $\frac{m}{s}$
$\psi_g =$ gravito-magnetic divergence potential, in $\frac{1}{s}$
$\kappa_g =$ speed of propagation of gravito-magnetic divergence, in $\frac{m}{s}$

GEM equations:
$[\frac{kg}{m^2 \cdot s}]: \partial_t (D_g)_i - (\mu_g)^{-1} {\epsilon_i}^{jk} \nabla_j (B_g)_k + \xi_g (\phi_g)_{,i} = -(J_g)_i$
$[\frac{1}{s^2}]: \partial_t (B_g)_i + (\epsilon_g)^{-1} {\epsilon_i}^{jk} \nabla_j (D_g)_k + \kappa_g (\psi_g)_{,i} = 0$
$[\frac{kg}{m^2 \cdot s}]: \partial_t (\phi_g) + \xi_g \nabla_j (D_g)^j = \xi_g \rho_g$
$[\frac{1}{s^2}]: \partial_t (\psi_g) + \kappa_g \nabla_j (B_g)^j = 0$

Decomposed into individual matrices:
Let $U = \left[\begin{matrix} (D_g)_i \\ (B_g)_i \\ \phi \\ \psi \end{matrix}\right], S = \left[\begin{matrix} -(J_g)_i \\ 0 \\ \xi_g \rho_g \\ 0 \end{matrix}\right]$

the system becomes:
$\partial_t U + A^j \nabla_j U = S$

The eigen decomposition of the GLM GEM system is similar to that of the GLM Maxwell system above.

---

fluid properties (for $\alpha = \{ ion, elec\}$):
$\rho_\alpha = $ ion and electron density, in $\frac{kg}{m^3}$
$v_\alpha = $ ion and electron velocity, in $\frac{m}{s}$
$P_\alpha = $ ion and electron pressure, in $\frac{kg}{m \cdot s^2}$
$\mathcal{E}_{\alpha,int} = \frac{P_\alpha}{\Gamma - 1} = $ densitized internal energy, related to the pressure by the ideal gas law, in $\frac{kg}{m \cdot s^2}$
for $\Gamma = \frac{5}{3} = $ the heat capacity ratio, in units $[1]$
$\mathcal{E}_{\alpha,kin} = \frac{1}{2} \rho_\alpha |v_\alpha|^2 = $ densitized kinetic energy, in $\frac{kg}{m \cdot s^2}$
$\mathcal{E}_{\alpha,total} = \mathcal{E}_{\alpha,int} + \mathcal{E}_{\alpha,kin} = $ densitized total energy, in $\frac{kg}{m \cdot s^2}$

$r_\alpha = \frac{q_\alpha}{m_\alpha}$ = charge/mass ratio of the particle of the species, in $\frac{C}{kg}$
$m_\alpha = $ ion or electron mass, in $kg$. In the simulation we set $m_{ion} = 1$. (TODO more on normalizing units.)

$v_{th,\alpha} = $ reference thermal velocity, in units of $[\frac{m}{s}]$.
Defined in 2003 Loverich as $v_{th,\alpha} = \sqrt{\frac{2 P_\alpha}{\rho_\alpha}}$
No mention of it in 2014 Abgrall, Kumar ... is it set to a constant? 2014 Abgrall, Kumar do use it for normalizing speed of light, etc (whereas 2003 Loverich uses a separate $U_0$). If it isn't constant in 2014 Abgrall, Kumar, then neither would other constants (like the speed of light) be constant.

Reference lengths:
$x_0 = $ reference length.
$v_0 = $ reference speed, in units of $[\frac{m}{s}]$. 2014 Abgrall, Kumar uses $v_{th,\alpha}$ in place of this (which I think it fixes to be a constant?).
$B_0 = $ reference magnetic field, in units of $[\frac{kg}{C \cdot s}]$. This is not the 0'th component of B, the magnetic field (which has a value of zero).

$W = $ Lorentz boost, in units of $[1]$. 2014 Abgrall, Kumar sets $W = 1$

$\omega_{L,\alpha} = r_\alpha B = \frac{q_\alpha B}{m_\alpha} =$ gyrofrequency, or Larmor frequency, in units of $[\frac{1}{s}]$. (from https://en.wikipedia.org/wiki/Gyroradius).

$r_L =$ Larmor radius, in units of $[m]$.
$r_L = \frac{W m_{ion} v_{th,ion}}{q_{ion} B_0} =$ according to 2014 Abgrall, Kumar.
$r_L = \frac{v_{th,ion} m c^2}{e |B|}$ according to this page: https://www.encyclopediaofmath.org/index.php/Larmor_radius.
$r_L = \frac{m_\alpha v_\perp}{|q_\alpha| B} = \frac{v_\perp}{\omega_{L,\alpha}}$ according to wikipedia: https://en.wikipedia.org/wiki/Gyroradius

$\lambda_D = $ ion Debye length
$\lambda_D = \sqrt{\frac{\epsilon_0 (v_{th,ion})^2}{ n_0 q_{ion}}}$ in 2014 Abgrall, Kumar... and I don't think they define $n_0$, but other papers do as number density ... and the units are messed up, in $m \cdot \sqrt{\frac{C}{kg}}$... if we square the charge (like everyone else does) then all we seem to be missing is the mass of the particle ... inserting that brings us to units of $m$.
$\lambda_D = \sqrt{\frac{\epsilon_0 v_0^2 m_0}{n_0 q_0^2}}$ in another source (?), which is in units of $m$. Seems to be the corrected source from above.
From https://en.wikipedia.org/wiki/Debye_length:
$\lambda_D = \sqrt{\frac{\epsilon_0 k_B / q_e^2}{n_e / T_e + \Sigma_j z_j^2 n_j / T_i}}$, for $z_j = q_j / q_{electron}$, in units of $m$
$\lambda_D = \sqrt{\frac{\epsilon_0 k_B T_e}{n_e q_e^2}}$ for isothermal cold plasma ... if $q_e$ is in $C$ and $n_e$ is in $\frac{1}{m^3}$ then this is in units of $m$.

Normalized units:
$x_0 =$ reference length, in $m$.
$v_0 =$ reference velocity, in $\frac{m}{s}$. 2014 Abgrall, Kumar uses $v_0 = v_{th,ion}$.

$\hat{r}_L = \frac{r_L}{x_0} =$ normalized ion Larmor radius, in units of $[1]$
$\hat{c} = \frac{c}{v_0}$ = normalized speed of light, in units of $[1]$
$\hat{\lambda}_D = \frac{\lambda_D}{r_L} =$ ion Debye length in units of the Larmor radius.

So, derived,
$L = \hat{\lambda}_D^2 \hat{r}_L = \frac{\epsilon (v_{th,ion})^2 }{r_L x_0 n_0 q_{ion}}$
$= \frac{\epsilon B_0 v_{th,ion}}{W x_0 n_0 m_i}$

ions:
$[\frac{kg}{m^3 \cdot s}]: \partial_t (\rho_{ion}) + \nabla_k (\rho_{ion} v_{ion}^k) = 0$
$[\frac{kg}{m^2 \cdot s^2}]: \partial_t ( \rho_{ion} v_{ion}^j ) + \nabla_k ( \rho_{ion} v_{ion}^j v_{ion}^k + P_{ion} g^{jk}) = \rho_{ion} ( (\hat{r}_L)^{-1} (E^j + \epsilon^{jkl} (v_{ion})_k B_l) + ((E_g)^j + 4 \epsilon^{jkl} (v_{ion})_k (B_g)_l) )$
$[\frac{kg}{m \cdot s^3}]: \partial_t ( \mathcal{E}_{ion,total} ) + \nabla_k (( E_{ion,total} + P_{ion}) v_{ion}^k) = \rho_{ion} v_{ion}^j ( (\hat{r}_L)^{-1} E_j + (E_g)^j )$

electrons:
$[\frac{kg}{m^3 \cdot s}]: \partial_t (\rho_{elec}) + \nabla_k (\rho_{elec} v_{elec}^k ) = 0$
$[\frac{kg}{m^2 \cdot s^2}]: \partial_t ( \rho_{elec} v_{elec}^j ) + \nabla_k ( \rho_{elec} v_{elec}^j v_{elec}^k + P_{elec} g^{jk}) = \rho_{elec} ( (\hat{r}_L)^{-1} (E^j + \epsilon^{jkl} (v_{elec})_k B_l) + ((E_g)^j + 4 \epsilon^{jkl} (v_{elec})_k (B_g)_l) )$
$[\frac{kg}{m \cdot s^3}]: \partial_t ( E_{elec,total} ) + \nabla_k (( E_{elec,total} + P_{elec} ) v_{elec}^k ) = \rho_{elec} v_{elec}^j ( (\hat{r}_L)^{-1} E_j + (E_g)^j )$