Coordinate basis:
$dx^\mu = (dx^0, dx^1, dx^2, dx^3) = (c dt, dx, dy, dz)$ has units $m$
$\partial_\mu = (\partial_0, \partial_1, \partial_2, \partial_3) = (\frac{1}{c} \partial_t, \partial_x, \partial_y, \partial_z)$ has units $\frac{1}{m}$

Non-coordinate basis:
$e^\mu = {e^\mu}_{\bar{\mu}} dx^{\bar\mu}$
$e_\mu = {e_\mu}^{\bar{\mu}} \partial_{\bar\mu}$

speed of light: $c = 299792458 \cdot \frac{m}{s}$

$\epsilon_0 =$ vacuum permittivity, in units $[\frac{C^2 \cdot s^2}{kg \cdot m^3}]$
$\mu_0 =$ vacuum permeability, in units $[\frac{kg \cdot m}{C^2}]$
$\frac{1}{\sqrt{\epsilon_0 \mu_0}} = c$


$E^i$ = electric field, in units $[\frac{kg \cdot m}{C \cdot s^2}]$
$B^i$ = magnetic field, in units $[\frac{kg}{C \cdot s}]$
$D^i = \epsilon E^i =$ electric displacement field, in units $[\frac{C}{m^2}]$
$B^i = \mu H^i$, for $H^i =$ magnetizing field, in units $[\frac{C}{m \cdot s}]$

$\epsilon =$ permittivity, in units $[\frac{C^2 \cdot s^2}{kg \cdot m^3}]$
$\mu =$ permeability, in units $[\frac{kg \cdot m}{C^2}]$
$v_p = \frac{1}{\sqrt{\epsilon \mu}} = $ phase velocity, in units $[\frac{m}{s}]$

Wedge product:
$a \wedge b = a \otimes b - b \otimes a$

Exterior derivative in a non-coordinate basis:
$dT = e^\mu \wedge e_\mu (T)$
$= (d {T^{\alpha_1 ... \alpha_p}}_{\beta_1 ... \beta_q} e_{\alpha_1} \otimes ... \otimes e_{\alpha_p} \otimes e^{\beta_1} \otimes ... \otimes e^{\beta_q}) = e_\mu ({T^{\alpha_1 ... \alpha_p}}_{\beta_1 ... \beta_q}) \cdot e^\mu \wedge ( e_{\alpha_1} \otimes ... \otimes e_{\alpha_p} \otimes e^{\beta_1} \otimes ... \otimes e^{\beta_q} )$

So for $T = T_{a_1 ... a_p}$,
$(dT)_{a_0 a_1 ... a_p} = (p+1) \partial_{[a_0} T_{a_1 ... a_p]}$

Electromagnetic four-potential:
$A = A_\mu e^\mu = A_0 e^0 + A_i e^i$
$A_\mu$ is in units of $\frac{kg \cdot m}{C \cdot s}$
$A_0 = \frac{1}{c} \phi$
$A = \phi dt + A_i e^i$
$\phi$ is in units of $\frac{kg \cdot m^2}{C \cdot s^2}$

Faraday tensor:
$F_{uv}$ has units of $\frac{kg}{C \cdot s}$
$F = F_{uv} e^u \otimes e^v = \frac{1}{2} F_{uv} e^u \wedge e^v$
$F = dA$
$= e^\mu \wedge e_\mu (A_\nu e^\nu)$
$= e_\mu (A_\nu) \cdot e^\mu \wedge e^\nu$
$= {e_\mu}^{\bar\mu} \partial_{\bar\mu} (A_\nu) \cdot e^\mu \wedge e^\nu$
$= 2 {e_{[\mu|}}^{\bar\mu} \partial_{\bar\mu} (A_{|\nu]}) \cdot e^\mu \otimes e^\nu$
$= ( {e_\mu}^{\bar\mu} \partial_{\bar\mu} (A_\nu) - {e_\nu}^{\bar\nu} \partial_{\bar\nu} (A_\mu) )\cdot e^\mu \otimes e^\nu$
$F_{\mu\nu} = {e_\mu}^{\bar\mu} \partial_{\bar\mu} (A_\nu) - {e_\nu}^{\bar\nu} \partial_{\bar\nu} (A_\mu) $

From 2009 Alcubierre, Degollado, Salgado:
Let $F^{ab} = n^a E^b - n^b E^a + \epsilon^{abcd} n_c B_d$
So $E^a = F^{ab} n_b$ for the hypersurface normal $n^a$
So $B^a = -\frac{1}{2} \epsilon^{abcd} n_b F_{cd}$
So $(\star F)_{uv} = \frac{1}{2} F^{ab} \epsilon_{abuv} = \frac{1}{2} (n^a E^b - n^b E^a + \epsilon^{abcd} n_c B_d) \epsilon_{abuv}$
i.e. $(\star F)^{ab} = -n^a B^b + n^b B^a + \epsilon^{abcd} n_c E_d$
So $B^a = -(\star F)^{ab} n_b$
So $E^a = -\frac{1}{2} \epsilon^{abcd} n_b (\star F)_{cd}$

From the 'relativistic electromagnetism metric in medium' worksheet:
(Notice that $\alpha, \beta^i, \gamma_{ij}$ are ADM metric parameters)
Let $\epsilon_{abcd} = \sqrt{|g|} [abcd] = \sqrt{|g|} \bar\epsilon_{abcd} $ is the spacetime Levi-Civita permutation tensor
So $\epsilon_{0ijk} = \alpha \sqrt\gamma \bar\epsilon_{0ijk} = \alpha \sqrt\gamma \bar\epsilon_{ijk}$ (see the ADM section for more details)
$\epsilon_{ijk} = \sqrt\gamma [ijk] = \sqrt\gamma \bar\epsilon_{ijk}$ is the spatial Levi-Civita permutation tensor
So $\epsilon_{0ijk} = \alpha \epsilon_{ijk}$
And ${\epsilon^0}_{ijk} = g^{0a} \epsilon_{aijk} = g^{00} \epsilon_{0ijk} + g^{0\Sigma} \epsilon_{\Sigma ijk}$, but since $ijk \in \Sigma$ we know $\epsilon_{\Sigma ijk} = 0$
So ${\epsilon^0}_{ijk} = -\frac{1}{\alpha^2} \epsilon_{0ijk} = -\frac{1}{\alpha} \epsilon_{ijk} = -\frac{1}{\alpha} \sqrt\gamma \bar\epsilon_{0ijk} = -\frac{1}{\alpha} \sqrt\gamma \bar\epsilon_{ijk}$
So ${\epsilon^{0i}}_{jk} = -\frac{1}{\alpha^2} \gamma^{im} \epsilon_{0mjk}$

NOTICE if we are going to decompose things in ADM formalism, then we must ensure only a coordinate basis, as per the conditions of using ADM formalism.

Write the Maxwell laws in macro form, using D instead of F.
$D_{ab} = \frac{1}{2} {Z_{ab}}^{cd} F_{cd}$
...where Z is antisymmetric between a and b, and between c and d, and symmatric between ab and cd, just like the Riemann tensor ... hmm ...

$d^2 A = 0$
$A_{[a,bc]} = 0$ ... yes, partials, even in curved spacetime. If you want to convert them to covariant derivatives then you have to add an extra $A_b {R^b}_d$ term.
$F_{[ab,c]} = 0$
$\star d^2 A = 0$
$\epsilon^{abcd} F_{[ab,c]} = 0$
$\epsilon^{abcd} F_{ab,c} = 0$

In vacuum:
$\star d \star d A = J$
$\epsilon^{abcd} (\star F)_{[ab,c]} = J^a$
macroscopic:
$\star d \star Z(d A) = J$
$\epsilon^{abcd} (\star D)_{[ab,c]} = J^a$
$\epsilon^{abcd} (\frac{1}{2} D_{uv} \epsilon^{uvab})_{[ab,c]} = J^a$

Formulate them as a first-order system wrt D and B and constraints wrt H and E.
For $\star d^2 A = 0$:
$\epsilon^{abcd} (n_a E_b - n_b E_a + \epsilon_{abuv} n^u B^v)_{,c} = 0$
$\epsilon^{abcd} ( n_{a,c} E_b + n_a E_{b,c} - n_{b,c} E_a - n_b E_{a,c} + (\sqrt{|g|})_{,c} \bar\epsilon_{abuv} n^u B^v + \epsilon_{abuv} {n^u}_{,c} B^v + \epsilon_{abuv} n^u {B^v}_{,c} ) = 0$
$\epsilon^{abcd} ( n_{a,c} E_b + n_a E_{b,c} - n_{b,c} E_a - n_b E_{a,c} ) + \epsilon^{abcd} \epsilon_{abuv} ( \frac{1}{\sqrt{|g|}} (\sqrt{|g|})_{,c} n^u B^v + {n^u}_{,c} B^v + n^u {B^v}_{,c} ) = 0$
$\epsilon^{abcd} ( n_{a,c} E_b + n_a E_{b,c} - n_{b,c} E_a - n_b E_{a,c} ) - 2 (\delta^c_u \delta^d_v - \delta^c_v \delta^d_u) ( {n^u}_{,c} B^v + n^u {B^v}_{,c} + \frac{1}{\sqrt{|g|}} (\sqrt{|g|})_{,c} n^u B^v ) = 0$
$2 \epsilon^{abcd} ( n_{a,c} E_b + n_a E_{b,c} ) - 2 ( {n^c}_{,c} B^d + n^c {B^d}_{,c} - {n^d}_{,c} B^c - n^d {B^c}_{,c} + \frac{1}{\sqrt{|g|}} (\sqrt{|g|})_{,c} ( n^c B^d - n^d B^c ) ) = 0$

First look at normal component:
$n_d ( \epsilon^{abcd} ( n_{a,c} E_b + n_a E_{b,c} ) - {n^c}_{,c} B^d - n^c {B^d}_{,c} + n^d {B^c}_{,c} + {n^d}_{,c} B^c - \frac{1}{\sqrt{|g|}} (\sqrt{|g|})_{,c} ( n^c B^d - n^d B^c ) ) = 0$
$\epsilon^{abcd} n_{a,c} E_b n_d - {B^c}_{,c} - n_d {B^d}_{,c} n^c + n_d {n^d}_{,c} B^c - \frac{1}{\sqrt{|g|}} (\sqrt{|g|})_{,c} B^c = 0$
$n_{a,c} \epsilon^{abcd} E_b n_d - {B^c}_{,c} + {}^4 {\Gamma^d}_{cd} B^c - n_d {B^d}_{,c} n^c + n_d {n^d}_{,c} B^c = 0$
Using $B^d n_d = 0$ therefore ${B^d}_{,c} n_d = -B^d n_{d,c}$
Using $n^d n_d = -1$ therefore ${n^d}_{,c} n_d = -n^d n_{d,c}$
$n_{a,b} (-n^a B^b + n^b B^a + \epsilon^{abcd} n_c E_d) - {B^c}_{,c} - {}^4 {\Gamma^d}_{cd} B^c = 0$
$n_{a,b} (\star F)^{ab} - {B^c}_{,c} - {}^4 {\Gamma^d}_{cd} B^c = 0$
Using $n_a (\star F)^{ab} = -B^b$, so $-{B^b}_{,b} = (n_a (\star F)^{ab})_{,b} = n_{a,b} (\star F)^{ab} + n_a {(\star F)^{ab}}_{,b}$
so $n_{a,b} (\star F)^{ab} = -{B^b}_{,b} - n_a {(\star F)^{ab}}_{,b}$
$- {B^c}_{,c} - n_a {(\star F)^{ab}}_{,b} - {B^c}_{,c} - {}^4 {\Gamma^d}_{cd} B^c = 0$
$-n_a {(\star F)^{ab}}_{,b} - 2 {B^c}_{,c} - {}^4 {\Gamma^d}_{cd} B^c = 0$
BLEH.
I'm just using the 2009 paper.
$\nabla^\perp_a D^a = \frac{1}{\sqrt\gamma} (\sqrt\gamma D^i)_{,i} = \rho = J^a n_a$
$\nabla^\perp_a B^a = \frac{1}{\sqrt\gamma} (\sqrt\gamma B^i)_{,i} = 0$
${D^i}_{,t} - \mathcal{L}_\beta D^i = (\nabla^\perp \times \alpha B)^i + \alpha K D^i - \alpha (\perp J)^i$
${B^i}_{,t} - \mathcal{L}_\beta B^i = -(\nabla^\perp \times \alpha D)^i + \alpha K B^i$

$(\sqrt\gamma D^i)_{,i} = \sqrt\gamma J^a n_a$

$(\sqrt\gamma B^i)_{,i} = 0$

Then find eigenvectors/eigenvalues.

$F = dA = (\partial_\mu A_\nu - \partial_\nu A_\mu) \cdot dx^\mu \otimes dx^\nu$
$F_{\mu\nu} = 2 \partial_{[\mu} A_{\nu]}$

for $D_{ab} = {T_{ab}}^{cd} F_{cd}$

$F_{uv} = \downarrow u(i) \overset{\rightarrow v(j)}{ \left[\begin{matrix} 0 & -\frac{1}{c} \alpha E_j \\ \frac{1}{c} \alpha E_i & {\epsilon_{ij}}^k B_k \end{matrix}\right] }$
$F_{uv} = \downarrow u(i) \overset{\rightarrow v(j)}{ \left[\begin{matrix} 0 & -\frac{1}{c} \alpha E_j \\ \frac{1}{c} \alpha E_i & \gamma \bar\epsilon_{ijk} B^k \end{matrix}\right] }$
In 3+1:
$F_{uv} = \left[\begin{matrix} 0 & -\frac{1}{c} \alpha E_1 & -\frac{1}{c} \alpha E_2 & -\frac{1}{c} \alpha E_3 \\ \frac{1}{c} \alpha E_1 & 0 & \gamma B^3 & -\gamma B^2 \\ \frac{1}{c} \alpha E_2 & -\gamma B^3 & 0 & \gamma B^1 \\ \frac{1}{c} \alpha E_3 & \gamma B^2 & -\gamma B^1 & 0 \end{matrix}\right]$

$\star F = \star (F_{ab} dx^a \wedge dx^b) = F_{ab} \star (dx^a \wedge dx^b) = (\star F)_{ab} dx^a \wedge dx^b$
(Notice that $(\star F)_{ab}$ is a variable denotation, and not an operator, while $\star dx^a$ is an operator.)
TODO don't forget to add in the influence from $g_{ab}$
$(\star F)_{ab} = \frac{1}{2} F_{uv} {\epsilon^{uv}}_{ab}$
$(\star F)_{0i} = \frac{1}{2} F_{jk} {\epsilon^{jk}}_{0i} = \frac{1}{2} F^{jk} \epsilon_{jk0i} = \frac{1}{2} F^{jk} \epsilon_{0ijk}$
$(\star F)_{0i} = \frac{1}{2} \alpha \epsilon_{ijk} g^{ja} g^{kb} F_{ab}$
$(\star F)_{0i} = \frac{1}{2} \alpha \epsilon_{ijk} (g^{j0} g^{k0} F_{00} + g^{j0} g^{km} F_{0m} + g^{jn} g^{k0} F_{n0} + g^{jm} g^{kn} F_{mn})$ $(\star F)_{0i} = \frac{1}{2} \alpha \epsilon_{ijk} (-\frac{1}{\alpha^2} \beta^j (\gamma^{km} - \frac{1}{\alpha^2} \beta^k \beta^m) \frac{1}{c} \alpha E_m + \frac{1}{\alpha^2} (\gamma^{jn} - \frac{1}{\alpha^2} \beta^j \beta^n) \beta^k \frac{1}{c} \alpha E_n + (\gamma^{jm} - \frac{1}{\alpha^2} \beta^j \beta^m) (\gamma^{kn} - \frac{1}{\alpha^2} \beta^k \beta^n) \epsilon_{mnp} B^p)$
$(\star F)_{0i} = \frac{1}{2} \alpha \epsilon_{ijk} (-2 \frac{1}{\alpha^2} \beta^j (\gamma^{km} - \frac{1}{\alpha^2} \beta^k \beta^m) \frac{1}{c} \alpha E_m + (\gamma^{jm} - \frac{1}{\alpha^2} \beta^j \beta^m) (\gamma^{kn} - \frac{1}{\alpha^2} \beta^k \beta^n) \epsilon_{mnp} B^p)$
$(\star F)_{0i} = \frac{1}{2} ( - 2 \frac{1}{c} \epsilon_{ijk} \beta^j E^k + \alpha \epsilon_{ijk} \epsilon^{mjk} B^m + 2 \frac{1}{\alpha} \epsilon_{ijk} \epsilon^{mnk} \beta^j \beta_m B_n )$
$(\star F)_{0i} = \frac{1}{2} ( - 2 \frac{1}{c} \epsilon_{ijk} \beta^j E^k + \alpha \delta^{mj}_{ij} B^m + 2 \frac{1}{\alpha} \delta^{mn}_{ij} \beta^j \beta_m B_n )$
$(\star F)_{0i} = \frac{1}{2} ( - 2 \frac{1}{c} \epsilon_{ijk} \beta^j E^k + \alpha (\delta^m_i \delta^j_j - \delta^m_j \delta^j_i) B^m + 2 \frac{1}{\alpha} (\delta^m_i \delta^n_j - \delta^m_j \delta^n_i) \beta^j \beta_m B_n )$
$(\star F)_{0i} = \alpha B^i + \frac{1}{\alpha} \beta_i \beta^j B_j - \frac{1}{\alpha} \beta^j \beta_j B_i - \frac{1}{c} \epsilon_{ijk} \beta^j E^k $
TODO here I think my sign is backwards.

$(\star F)_{jk} = \frac{1}{2} F_{0i} {\epsilon^{0i}}_{jk}$
$(\star F)_{jk} = \frac{1}{2} F_{0i} g^{0a} g^{ib} \epsilon_{abjk}$

$(\star F)_{ab} = \left[\begin{matrix} 0 & -\alpha B_1 & -\alpha B_2 & -\alpha B_3 \\ \alpha B_1 & 0 & -\frac{1}{c} E_3 & \frac{1}{c} E_2 \\ \alpha B_2 & \frac{1}{c} E_3 & 0 & -\frac{1}{c} E_1 \\ \alpha B_3 & -\frac{1}{c} E_2 & \frac{1}{c} E_1 & 0 \end{matrix}\right]$

Deriving Maxwell equations:
$dF = d^2 A = 0$
$0 = \partial_a F_{bc} dx^a \wedge dx^b \wedge dx^c$
becomes:
$0 = (F_{23,1} + F_{31,2} + F_{12,3} - F_{32,1} - F_{13,2} - F_{21,3}) dx^1 \wedge dx^2 \wedge dx^3$
$0 = (F_{23,0} + F_{30,2} + F_{02,3} - F_{32,0} - F_{03,2} - F_{20,3}) dx^0 \wedge dx^2 \wedge dx^3$
$0 = (F_{13,0} + F_{30,1} + F_{01,3} - F_{31,0} - F_{03,1} - F_{10,3}) dx^0 \wedge dx^1 \wedge dx^3$
$0 = (F_{12,0} + F_{20,1} + F_{01,2} - F_{21,0} - F_{02,1} - F_{10,2}) dx^0 \wedge dx^1 \wedge dx^2$
becomes:
$B_{1,1} + B_{2,2} + B_{3,3} = 0$
$B_{1,t} + E_{3,2} - E_{2,3} = 0$
$B_{2,t} + E_{1,3} - E_{3,1} = 0$
$B_{3,t} + E_{2,1} - E_{1,2} = 0$

$D_{uv} = {T_{uv}}^{ab} F_{ab}$
$D_{0i} = c^2 \epsilon F_{0i}$
$D_{i0} = c^2 \epsilon F_{i0}$
$D_{ij} = \frac{1}{\mu} F_{ij}$
Therefore
$c^2 \epsilon = \frac{\epsilon}{\epsilon_0} \frac{1}{\mu_0} = {T_{0i}}^{0i} = {T_{i0}}^{i0} = -{T_{0i}}^{i0} = -{T_{i0}}^{0i}$
$\frac{1}{\mu} = \frac{\mu}{\mu_0} \frac{1}{\mu_0} = {T_{jk}}^{jk} = {T_{kj}}^{kj} = -{T_{jk}}^{kj} = -{T_{kj}}^{jk}$
${T_{aa}}^{bc} = {T_{bc}}^{aa} = 0$

$D_{uv} = \downarrow u(i) \overset{\rightarrow v(i)}{ \left[\begin{matrix} 0 & -c D_j \\ c D_i & {\bar\epsilon_{ij}}^k H_k \end{matrix}\right]}$

In 3+1:
$D_{uv} = \left[\begin{matrix} 0 & -c D_1 & -c D_2 & -c D_3 \\ c D_1 & 0 & -H_3 & H_2 \\ c D_2 & H_3 & 0 & -H_1 \\ c D_3 & -H_2 & H_1 & 0 \end{matrix}\right]$

$(\star D)_{ab} = \left[\begin{matrix} 0 & -H_1 & -H_2 & -H_3 \\ H_1 & 0 & -c D_3 & c D_2 \\ H_2 & c D_3 & 0 & -c D_1 \\ H_3 & -c D_2 & c D_1 & 0 \end{matrix}\right]$

$J = J_u dx^u = c \rho_{free} dx^0 + J_i dx^i = c^2 \rho_{free} dt + J_i dx^i$
$J_i$ is of units $[\frac{C}{m^2 \cdot s}]$.
$\rho_{charge}$ is of units $[\frac{C}{m^3}]$.
$\rho_{free}$ is of units $[\frac{C}{m^3}]$.
TODO I'm using $\rho_{charge}$ in the expression for the Lorentz force law. Are these expressions supposed to use the free or bound or total charge?
TODO The displacemet field below uses the free charge only. That's why I'm assigning $\delta D = J$, but is this true in general, or should the source term of this equation be modified?
TODO Notice, the Wikipedia 'Four-Current' page defines the four-current in contravariant terms, which means it will differ from my definition by the metric lapse and shift.

Lorentz force law:
$\ddot{x}^a = \frac{q}{m} {F^a}_b \dot{x}^b$
...where m = mass of particle at x, q = charge of particle at x, and $F_{ab}$ is the Faraday tensor equal to $F = dA$, made up of E and B,
where B is constrained by $\nabla \cdot B = 0$ because of the property $dF = d^2 A = 0$,
and $E = \epsilon^{-1}(D)$ is constrained by $\nabla \cdot D = \rho_{free}$ because of the property $\delta D = J$,
and ${F^a}_b = g^{ac} F_{cb}$ for $g^{ab}$ the inverse of the metric tensor of the manifold that this takes place on.
The metric itself $g_{ab}$ is derived from the mass (of the other particles?), $D^2 \omega $'s contracted trace reversal $= R_{uv} = 8 \pi T_{uv}$


${(\star F)^{\alpha\beta}}_{;\beta} - A_\beta {R^\beta}_\alpha = 0$

${F^{\alpha\beta}}_{;\beta} - A_\beta {R^\beta}_\alpha = \mu_0 J^\alpha$

Maxwell Equations of charge

$B_{x,x} + B_{y,y} + B_{z,z} = 0$
$B_{x,t} + E_{z,y} - E_{y,z} = 0$
$B_{y,t} + E_{x,z} - E_{z,x} = 0$
$B_{z,t} + E_{y,x} - E_{x,y} = 0$
$D_{x,x} + D_{y,y} + D_{z,z} = \rho_{free}$
$D_{x,t} - H_{z,y} + H_{y,z} = -J_x$
$D_{y,t} - H_{x,z} + H_{z,x} = -J_y$
$D_{z,t} - H_{y,x} + H_{x,y} = -J_z$

in vector notation:

$\partial_t \left[ \begin{matrix} D_x \\ D_y \\ D_z \\ B_x \\ B_y \\ B_z \end{matrix} \right] + \vec\partial \cdot \left[\begin{matrix} 0 & -H_z & H_y \\ H_z & 0 & -H_x \\ -H_y & H_x & 0 \\ 0 & E_z & -E_y \\ -E_z & 0 & E_x \\ E_y & -E_x & 0 \end{matrix}\right] = \left[\begin{matrix} -J_x \\ -J_y \\ -J_z \\ 0 \\ 0 \\ 0 \end{matrix}\right]$

in index notation:

$D_{i,i} = \rho_{free}$
$B_{i,i} = 0$

$D_{i,t} - {\bar\epsilon_i}^{jk} H_{k,j} = -J_i$
$B_{i,t} + {\bar\epsilon_i}^{jk} E_{k,j} = 0$

$J_i = \sigma E_i = $ current density
$\sigma = $ electric conductivity

conservation form:
$\int (\partial_t U_I + \nabla_n F_I) dV = S_I$

conservation variables:
$U_I = \left[\begin{matrix} D_i \\ B_i \end{matrix}\right]$

flux:
$F_I(n) = \left[\begin{matrix} -{\bar\epsilon_i}^{kj} H_j n_k \\ {\bar\epsilon_i}^{kj} E_j n_k \end{matrix}\right]$

source:
$S_I = \left[\begin{matrix} -J_i \\ 0 \end{matrix}\right]$

conservation form:
$U_I + \frac{\partial F_I (n^k)}{\partial U_J} \frac{\partial U_J}{\partial x^k} = S_I$

index/matrix form:
$\left[\matrix{ D_i \\ B_i }\right] + \left[\matrix{ 0 & - \frac{1}{\mu} {\bar\epsilon_i}^{jk} \\ \frac{1}{\epsilon} {\bar\epsilon_i}^{jk} & 0 }\right] \left[\matrix{ D_k \\ B_k }\right]_{,j} = \left[\matrix{ -J_i \\ 0 }\right]$

$\frac{\partial F(n)}{\partial U} = \left[\begin{matrix} 0 & 0 & 0 & 0 & \frac{1}{\mu} n^z & -\frac{1}{\mu} n^y \\ 0 & 0 & 0 & -\frac{1}{\mu} n^z & 0 & \frac{1}{\mu} n^x \\ 0 & 0 & 0 & \frac{1}{\mu} n^y & -\frac{1}{\mu} n^x & 0 \\ 0 & -\frac{1}{\epsilon} n^z & \frac{1}{\epsilon} n^y & 0 & 0 & 0 \\ \frac{1}{\epsilon} n^z & 0 & -\frac{1}{\epsilon} n^x & 0 & 0 & 0 \\ -\frac{1}{\epsilon} n^y & \frac{1}{\epsilon} n^x & 0 & 0 & 0 & 0 \end{matrix}\right]$

For $a = a_i dx^i = a_x dx + a_y dy + a_z dz$,
Let the one-form $a$ be represented in matrix form as a row vector: $a^\flat = \left[\begin{matrix} a_x & a_y & a_z \end{matrix}\right]$.

Let $\star a = a_i \star dx^i = a_x dy \wedge dz + a_y dz \wedge dx + a_z dx \wedge dy = a_x (dy \otimes dz - dz \otimes dy) + a_y (dz \otimes dx - dx \otimes dz) + a_z (dx \otimes dy - dy \otimes dz)$.
$\star a = \left[\begin{matrix} 0 & a_z & -a_y \\ -a_z & 0 & a_x \\ a_y & -a_x & 0 \end{matrix}\right]$,
such that $a \cdot (\star b) = \star (a \wedge b) = a \times b$
Then $\frac{\partial F(n)}{\partial U} = \left[\begin{matrix} 0 & \star n \frac{1}{\mu} \\ -\star n \frac{1}{\epsilon} & 0 \end{matrix}\right]$

Let $\{ n, n_2, n_3 \}$ form an orthonormal basis, such that $g_{ij} (n_a)^i (n_b)^j = \delta_{ab}$ and $\bar\epsilon_{ijk} n^i (n_2)^j (n_3)^k = 1$

eigen decomposition:
$\frac{\partial F}{\partial U} = R \Lambda L$



solving for eigenvalues:
$A v = \lambda v$ $(A - \lambda I) v = 0$
$\left|\begin{matrix} -\lambda & 0 & 0 & 0 & \frac{1}{\mu} n^z & -\frac{1}{\mu} n^y \\ 0 & -\lambda & 0 & -\frac{1}{\mu} n^z & 0 & \frac{1}{\mu} n^x \\ 0 & 0 & -\lambda & \frac{1}{\mu} n^y & -\frac{1}{\mu} n^x & 0 \\ 0 & -\frac{1}{\epsilon} n^z & \frac{1}{\epsilon} n^y & -\lambda & 0 & 0 \\ \frac{1}{\epsilon} n^z & 0 & -\frac{1}{\epsilon} n^x & 0 & -\lambda & 0 \\ -\frac{1}{\epsilon} n^y & \frac{1}{\epsilon} n^x & 0 & 0 & 0 & -\lambda \end{matrix}\right| = 0$

$-\lambda \left|\begin{matrix} -\lambda & 0 & -\frac{1}{\mu} n^z & 0 & \frac{1}{\mu} n^x \\ 0 & -\lambda & \frac{1}{\mu} n^y & -\frac{1}{\mu} n^x & 0 \\ -\frac{1}{\epsilon} n^z & \frac{1}{\epsilon} n^y & -\lambda & 0 & 0 \\ 0 & -\frac{1}{\epsilon} n^x & 0 & -\lambda & 0 \\ \frac{1}{\epsilon} n^x & 0 & 0 & 0 & -\lambda \end{matrix}\right| + \frac{1}{\mu} n^z \left|\begin{matrix} 0 & -\lambda & 0 & -\frac{1}{\mu} n^z & \frac{1}{\mu} n^x \\ 0 & 0 & -\lambda & \frac{1}{\mu} n^y & 0 \\ 0 & -\frac{1}{\epsilon} n^z & \frac{1}{\epsilon} n^y & -\lambda & 0 \\ \frac{1}{\epsilon} n^z & 0 & -\frac{1}{\epsilon} n^x & 0 & 0 \\ -\frac{1}{\epsilon} n^y & \frac{1}{\epsilon} n^x & 0 & 0 & -\lambda \end{matrix}\right| + \frac{1}{\mu} n^y \left|\begin{matrix} 0 & -\lambda & 0 & -\frac{1}{\mu} n^z & 0 \\ 0 & 0 & -\lambda & \frac{1}{\mu} n^y & -\frac{1}{\mu} n^x \\ 0 & -\frac{1}{\epsilon} n^z & \frac{1}{\epsilon} n^y & -\lambda & 0 \\ \frac{1}{\epsilon} n^z & 0 & -\frac{1}{\epsilon} n^x & 0 & -\lambda \\ -\frac{1}{\epsilon} n^y & \frac{1}{\epsilon} n^x & 0 & 0 & 0 \end{matrix}\right| = 0$

$-\lambda \left( -\lambda \left|\begin{matrix} -\lambda & \frac{1}{\mu} n^y & -\frac{1}{\mu} n^x & 0 \\ \frac{1}{\epsilon} n^y & -\lambda & 0 & 0 \\ -\frac{1}{\epsilon} n^x & 0 & -\lambda & 0 \\ 0 & 0 & 0 & -\lambda \end{matrix}\right| - \frac{1}{\mu} n^z \left|\begin{matrix} 0 & -\lambda & -\frac{1}{\mu} n^x & 0 \\ -\frac{1}{\epsilon} n^z & \frac{1}{\epsilon} n^y & 0 & 0 \\ 0 & -\frac{1}{\epsilon} n^x & -\lambda & 0 \\ \frac{1}{\epsilon} n^x & 0 & 0 & -\lambda \end{matrix}\right| + \frac{1}{\mu} n^x \left|\begin{matrix} 0 & -\lambda & \frac{1}{\mu} n^y & -\frac{1}{\mu} n^x \\ -\frac{1}{\epsilon} n^z & \frac{1}{\epsilon} n^y & -\lambda & 0 \\ 0 & -\frac{1}{\epsilon} n^x & 0 & -\lambda \\ \frac{1}{\epsilon} n^x & 0 & 0 & 0 \end{matrix}\right| \right) + \frac{1}{\mu} n^z \left( \lambda \left|\begin{matrix} 0 & -\lambda & \frac{1}{\mu} n^y & 0 \\ 0 & \frac{1}{\epsilon} n^y & -\lambda & 0 \\ \frac{1}{\epsilon} n^z & -\frac{1}{\epsilon} n^x & 0 & 0 \\ -\frac{1}{\epsilon} n^y & 0 & 0 & -\lambda \end{matrix}\right| + \frac{1}{\mu} n^z \left|\begin{matrix} 0 & 0 & -\lambda & 0 \\ 0 & -\frac{1}{\epsilon} n^z & \frac{1}{\epsilon} n^y & 0 \\ \frac{1}{\epsilon} n^z & 0 & -\frac{1}{\epsilon} n^x & 0 \\ -\frac{1}{\epsilon} n^y & \frac{1}{\epsilon} n^x & 0 & -\lambda \end{matrix}\right| + \frac{1}{\mu} n^x \left|\begin{matrix} 0 & 0 & -\lambda & \frac{1}{\mu} n^y \\ 0 & -\frac{1}{\epsilon} n^z & \frac{1}{\epsilon} n^y & -\lambda \\ \frac{1}{\epsilon} n^z & 0 & -\frac{1}{\epsilon} n^x & 0 \\ -\frac{1}{\epsilon} n^y & \frac{1}{\epsilon} n^x & 0 & 0 \end{matrix}\right| \right) + \frac{1}{\mu} n^y \left( \lambda \left|\begin{matrix} 0 & -\lambda & \frac{1}{\mu} n^y & -\frac{1}{\mu} n^x \\ 0 & \frac{1}{\epsilon} n^y & -\lambda & 0 \\ \frac{1}{\epsilon} n^z & -\frac{1}{\epsilon} n^x & 0 & -\lambda \\ -\frac{1}{\epsilon} n^y & 0 & 0 & 0 \end{matrix}\right| + \frac{1}{\mu} n^z \left|\begin{matrix} 0 & 0 & -\lambda & -\frac{1}{\mu} n^x \\ 0 & -\frac{1}{\epsilon} n^z & \frac{1}{\epsilon} n^y & 0 \\ \frac{1}{\epsilon} n^z & 0 & -\frac{1}{\epsilon} n^x & -\lambda \\ -\frac{1}{\epsilon} n^y & \frac{1}{\epsilon} n^x & 0 & 0 \end{matrix}\right| \right) = 0$

$-\lambda \left( \lambda^2 \left|\begin{matrix} -\lambda & \frac{1}{\mu} n^y & -\frac{1}{\mu} n^x \\ \frac{1}{\epsilon} n^y & -\lambda & 0 \\ -\frac{1}{\epsilon} n^x & 0 & -\lambda \end{matrix}\right| + \frac{1}{\mu} n^z \lambda \left|\begin{matrix} 0 & -\lambda & -\frac{1}{\mu} n^x \\ -\frac{1}{\epsilon} n^z & \frac{1}{\epsilon} n^y & 0 \\ 0 & -\frac{1}{\epsilon} n^x & -\lambda \\ \end{matrix}\right| - \frac{1}{\mu \epsilon} (n^x)^2 \left|\begin{matrix} -\lambda & \frac{1}{\mu} n^y & -\frac{1}{\mu} n^x \\ \frac{1}{\epsilon} n^y & -\lambda & 0 \\ -\frac{1}{\epsilon} n^x & 0 & -\lambda \end{matrix}\right| \right) + \frac{1}{\mu} n^z \left( -\lambda^2 \left|\begin{matrix} 0 & \frac{1}{\epsilon} n^y & -\lambda \\ \frac{1}{\epsilon} n^z & -\frac{1}{\epsilon} n^x & 0 \\ -\frac{1}{\epsilon} n^y & 0 & 0 \end{matrix}\right| - \frac{1}{\mu} n^z \lambda \left|\begin{matrix} 0 & 0 & -\lambda \\ 0 & -\frac{1}{\epsilon} n^z & \frac{1}{\epsilon} n^y \\ \frac{1}{\epsilon} n^z & 0 & -\frac{1}{\epsilon} n^x \end{matrix}\right| + -\lambda \frac{1}{\mu} n^x \left|\begin{matrix} 0 & -\frac{1}{\epsilon} n^z & -\lambda \\ \frac{1}{\epsilon} n^z & 0 & 0 \\ -\frac{1}{\epsilon} n^y & \frac{1}{\epsilon} n^x & 0 \end{matrix}\right| - \frac{1}{\mu^2} n^x n^y \left|\begin{matrix} 0 & -\frac{1}{\epsilon} n^z & \frac{1}{\epsilon} n^y \\ \frac{1}{\epsilon} n^z & 0 & -\frac{1}{\epsilon} n^x \\ -\frac{1}{\epsilon} n^y & \frac{1}{\epsilon} n^x & 0 \end{matrix}\right| \right) + \frac{1}{\mu} n^y \left( \lambda \frac{1}{\epsilon} n^y \left|\begin{matrix} -\lambda & \frac{1}{\mu} n^y & -\frac{1}{\mu} n^x \\ \frac{1}{\epsilon} n^y & -\lambda & 0 \\ -\frac{1}{\epsilon} n^x & 0 & -\lambda \end{matrix}\right| -\lambda \frac{1}{\mu} n^z \left|\begin{matrix} 0 & -\frac{1}{\epsilon} n^z & 0 \\ \frac{1}{\epsilon} n^z & 0 & -\lambda \\ -\frac{1}{\epsilon} n^y & \frac{1}{\epsilon} n^x & 0 \end{matrix}\right| + \frac{1}{\mu^2} n^z n^x \left|\begin{matrix} 0 & -\frac{1}{\epsilon} n^z & \frac{1}{\epsilon} n^y \\ \frac{1}{\epsilon} n^z & 0 & -\frac{1}{\epsilon} n^x \\ -\frac{1}{\epsilon} n^y & \frac{1}{\epsilon} n^x & 0 \end{matrix}\right| \right) = 0$

$-\lambda \left( -\frac{1}{\mu} n^x \lambda^2 \left|\begin{matrix} \frac{1}{\epsilon} n^y & -\lambda \\ -\frac{1}{\epsilon} n^x & 0 \end{matrix}\right| - \lambda^3 \left|\begin{matrix} -\lambda & \frac{1}{\mu} n^y \\ \frac{1}{\epsilon} n^y & -\lambda \end{matrix}\right| + \frac{1}{\mu} n^z \lambda \frac{1}{\epsilon} n^z \left|\begin{matrix} -\lambda & -\frac{1}{\mu} n^x \\ -\frac{1}{\epsilon} n^x & -\lambda \end{matrix}\right| + \frac{1}{\mu \epsilon} (n^x)^2 \frac{1}{\mu} n^x \left|\begin{matrix} \frac{1}{\epsilon} n^y & -\lambda \\ -\frac{1}{\epsilon} n^x & 0 \end{matrix}\right| + \frac{1}{\mu \epsilon} (n^x)^2 \lambda \left|\begin{matrix} -\lambda & \frac{1}{\mu} n^y \\ \frac{1}{\epsilon} n^y & -\lambda \end{matrix}\right| \right) + \frac{1}{\mu} n^z \left( + \lambda^3 \left|\begin{matrix} \frac{1}{\epsilon} n^z & -\frac{1}{\epsilon} n^x \\ -\frac{1}{\epsilon} n^y & 0 \end{matrix}\right| - \frac{1}{\mu} n^z \lambda \frac{1}{\epsilon} n^z \left|\begin{matrix} 0 & -\lambda \\ -\frac{1}{\epsilon} n^z & \frac{1}{\epsilon} n^y \\ \end{matrix}\right| + \lambda^2 \frac{1}{\mu} n^x \left|\begin{matrix} \frac{1}{\epsilon} n^z & 0 \\ -\frac{1}{\epsilon} n^y & \frac{1}{\epsilon} n^x \end{matrix}\right| - \frac{1}{\mu^2} n^x n^y \frac{1}{\epsilon} n^z \left|\begin{matrix} \frac{1}{\epsilon} n^z & -\frac{1}{\epsilon} n^x \\ -\frac{1}{\epsilon} n^y & 0 \end{matrix}\right| - \frac{1}{\mu^2} n^x n^y \frac{1}{\epsilon} n^y \left|\begin{matrix} \frac{1}{\epsilon} n^z & 0 \\ -\frac{1}{\epsilon} n^y & \frac{1}{\epsilon} n^x \end{matrix}\right| \right) + \frac{1}{\mu} n^y \left( -\lambda \frac{1}{\epsilon} n^y \frac{1}{\mu} n^x \left|\begin{matrix} \frac{1}{\epsilon} n^y & -\lambda \\ -\frac{1}{\epsilon} n^x & 0 \end{matrix}\right| - \lambda^2 \frac{1}{\epsilon} n^y \left|\begin{matrix} -\lambda & \frac{1}{\mu} n^y \\ \frac{1}{\epsilon} n^y & -\lambda \end{matrix}\right| -\lambda^2 \frac{1}{\mu} n^z \left|\begin{matrix} 0 & -\frac{1}{\epsilon} n^z \\ -\frac{1}{\epsilon} n^y & \frac{1}{\epsilon} n^x \end{matrix}\right| + \frac{1}{\mu^2} n^z n^x \frac{1}{\epsilon} n^z \left|\begin{matrix} \frac{1}{\epsilon} n^z & -\frac{1}{\epsilon} n^x \\ -\frac{1}{\epsilon} n^y & 0 \end{matrix}\right| + \frac{1}{\mu^2} n^z n^x \frac{1}{\epsilon} n^y \left|\begin{matrix} \frac{1}{\epsilon} n^z & 0 \\ -\frac{1}{\epsilon} n^y & \frac{1}{\epsilon} n^x \end{matrix}\right| \right) = 0$

$-\lambda \left( \frac{1}{\epsilon\mu} \lambda^3 (n^x)^2 - \lambda^3 (\lambda^2 - \frac{1}{\epsilon\mu} (n^y)^2) + \frac{1}{\epsilon\mu} (n^z)^2 \lambda (\lambda^2 - \frac{1}{\epsilon\mu} (n^x)^2 ) - \frac{1}{\mu^2 \epsilon^2} (n^x)^4 \lambda + \frac{1}{\mu \epsilon} (n^x)^2 \lambda (\lambda^2 - \frac{1}{\epsilon\mu} (n^y)^2 ) \right) + \frac{1}{\mu} n^z \left( + \lambda^3 \frac{1}{\epsilon^2} n^x n^y + \frac{1}{\mu} n^z \lambda \frac{1}{\epsilon} n^z \lambda \frac{1}{\epsilon} n^z + \lambda^2 \frac{1}{\mu} n^x \frac{1}{\epsilon} n^z \frac{1}{\epsilon} n^x + \frac{1}{\mu^2} n^x n^y \frac{1}{\epsilon} n^z \frac{1}{\epsilon} n^x \frac{1}{\epsilon} n^y - \frac{1}{\mu^2} n^x n^y \frac{1}{\epsilon} n^y \frac{1}{\epsilon} n^z \frac{1}{\epsilon} n^x \right) + \frac{1}{\mu} n^y \left( +\lambda^2 \frac{1}{\epsilon} n^y \frac{1}{\epsilon\mu} (n^x)^2 - \lambda^2 \frac{1}{\epsilon} n^y (\lambda^2 - \frac{1}{\epsilon\mu} (n^y)^2 ) + \lambda^2 \frac{1}{\mu} n^z \frac{1}{\epsilon^2} n^z n^y - \frac{1}{\mu^2} n^z n^x \frac{1}{\epsilon} n^z \frac{1}{\epsilon^2} n^x n^y + \frac{1}{\epsilon^3 \mu^2} (n^x)^2 n^y (n^z)^2 \right) = 0$

$ - \lambda^4 \frac{1}{\epsilon\mu} (n^x)^2 + \lambda^6 - \lambda^4 \frac{1}{\epsilon\mu} (n^y)^2 - \lambda^4 \frac{1}{\epsilon\mu} (n^z)^2 + \lambda^2 \frac{1}{\epsilon^2 \mu^2} (n^z)^2 (n^x)^2 + \lambda^2 \frac{1}{\mu^2 \epsilon^2} (n^x)^4 - \lambda^4 \frac{1}{\epsilon\mu} (n^x)^2 + \lambda^2 \frac{1}{\epsilon^2 \mu^2} (n^x)^2 (n^y)^2 + \lambda^3 \frac{1}{\epsilon^2} \frac{1}{\mu} n^x n^y n^z + \lambda^2 \frac{1}{\epsilon^2 \mu^2} (n^z)^4 + \lambda^2 \frac{1}{\epsilon^2 \mu^2} (n^x)^2 (n^z)^2 + \lambda^2 \frac{1}{\epsilon^2 \mu^2} (n^x)^2 (n^y)^2 - \lambda^4 \frac{1}{\epsilon \mu} (n^y)^2 + \lambda^2 \frac{1}{\epsilon^2 \mu^2} (n^y)^4 + \lambda^2 \frac{1}{\epsilon^2 \mu^2} n^z n^z n^y n^y + \frac{1}{\epsilon^3 \mu^3} (n^x)^2 (n^y)^2 (n^z)^2 - \frac{1}{\epsilon^3 \mu^3} (n^x)^2 (n^y)^2 (n^z)^2 - \frac{1}{\epsilon^3 \mu^3} (n^x)^2 (n^y)^2 (n^z)^2 + \frac{1}{\epsilon^3 \mu^3} (n^x)^2 (n^y)^2 (n^z)^2 = 0$

$ \lambda^2 ( + \lambda^4 - \lambda^2 \frac{1}{\epsilon\mu} (2 (n^x)^2 + 2 (n^y)^2 + (n^z)^2) + \lambda \frac{1}{\epsilon^2} \frac{1}{\mu} n^x n^y n^z + \frac{1}{\epsilon^2 \mu^2} ( + (n^x)^4 + (n^y)^4 + (n^z)^4 + 2 (n^x)^2 (n^y)^2 + 2 (n^x)^2 (n^z)^2 + (n^z)^2 (n^y)^2 ) ) = 0$

TODO should be...

$ \lambda^2 ( + \lambda^4 - \lambda^2 \frac{1}{\epsilon\mu} (2 (n^x)^2 + 2 (n^y)^2 + 2 (n^z)^2) + \frac{1}{\epsilon^2 \mu^2} ( + (n^x)^4 + (n^y)^4 + (n^z)^4 + 2 (n^x)^2 (n^y)^2 + 2 (n^x)^2 (n^z)^2 + 2 (n^z)^2 (n^y)^2 ) ) = 0$

Let $|n|^2 = (n^x)^2 + (n^y)^2 + (n^z)^2$

$ \lambda^2 ( \lambda^4 - 2 \lambda^2 \frac{1}{\epsilon\mu} |n|^2 + \frac{1}{\epsilon^2 \mu^2} |n|^4 ) = 0$

$ \lambda^2 ( \lambda^2 - \frac{1}{\epsilon\mu} |n|^2 )^2 = 0$

$ \lambda^2 (\lambda + \frac{1}{\sqrt{\epsilon\mu}} |n|)^2 (\lambda - \frac{1}{\sqrt{\epsilon\mu}} |n|)^2 = 0$

Eigenvalues:
$\lambda = 0$ multiplicity 2
$\lambda = -\frac{|n|}{\sqrt{\epsilon\mu}}$ multiplicity 2
$\lambda = \frac{|n|}{\sqrt{\epsilon\mu}}$ multiplicity 2

Eigenvectors of $\lambda = 0$:

$\left[\begin{matrix} 0 & 0 & 0 & 0 & \frac{1}{\mu} n^z & -\frac{1}{\mu} n^y \\ 0 & 0 & 0 & -\frac{1}{\mu} n^z & 0 & \frac{1}{\mu} n^x \\ 0 & 0 & 0 & \frac{1}{\mu} n^y & -\frac{1}{\mu} n^x & 0 \\ 0 & -\frac{1}{\epsilon} n^z & \frac{1}{\epsilon} n^y & 0 & 0 & 0 \\ \frac{1}{\epsilon} n^z & 0 & -\frac{1}{\epsilon} n^x & 0 & 0 & 0 \\ -\frac{1}{\epsilon} n^y & \frac{1}{\epsilon} n^x & 0 & 0 & 0 & 0 \end{matrix}\right] \left[\begin{matrix} X_{D_x} \\ X_{D_y} \\ X_{D_z} \\ X_{B_x} \\ X_{B_y} \\ X_{B_z} \end{matrix}\right] = 0$

$\left[\begin{matrix} n^y & -n^x & 0 & 0 & 0 & 0 \\ n^z & 0 & -n^x & 0 & 0 & 0 \\ 0 & n^z & -n^y & 0 & 0 & 0 \\ 0 & 0 & 0 & n^y & -n^x & 0 \\ 0 & 0 & 0 & n^z & 0 & -n^x \\ 0 & 0 & 0 & 0 & n^z & -n^y \end{matrix}\right] \left[\begin{matrix} X_{D_x} \\ X_{D_y} \\ X_{D_z} \\ X_{B_x} \\ X_{B_y} \\ X_{B_z} \end{matrix}\right] = 0$

$\left[\begin{matrix} n^y & -n^x & 0 & 0 & 0 & 0 \\ 0 & n^z & -n^y & 0 & 0 & 0 \\ 0 & n^z & -n^y & 0 & 0 & 0 \\ 0 & 0 & 0 & n^y & -n^x & 0 \\ 0 & 0 & 0 & n^z & 0 & -n^x \\ 0 & 0 & 0 & 0 & n^z & -n^y \end{matrix}\right] \left[\begin{matrix} X_{D_x} \\ X_{D_y} \\ X_{D_z} \\ X_{B_x} \\ X_{B_y} \\ X_{B_z} \end{matrix}\right] = 0$

$\left[\begin{matrix} n^y & -n^x & 0 & 0 & 0 & 0 \\ 0 & n^z & -n^y & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & n^y & -n^x & 0 \\ 0 & 0 & 0 & 0 & n^z & -n^y \\ 0 & 0 & 0 & 0 & 0 & 0 \end{matrix}\right] \left[\begin{matrix} X_{D_x} \\ X_{D_y} \\ X_{D_z} \\ X_{B_x} \\ X_{B_y} \\ X_{B_z} \end{matrix}\right] = 0$

Let $u_1 = X_{D_z}, u_2 = X_{B_z}$

$X_{D_y} n^z - X_{D_z} n^y = 0$
$X_{D_y} n^z - u_1 n^y = 0$
$X_{D_y} = u_1 \frac{n^y}{n^z}$

$X_{D_x} n^y - X_{D_y} n^x = 0$
$X_{D_x} n^y - u_1 \frac{n^y}{n^z} n^x = 0$
$X_{D_x} = u_1 \frac{n^x}{n^z}$
$\left[\begin{matrix} X_{D_x} \\ X_{D_y} \\ X_{D_z} \\ X_{B_x} \\ X_{B_y} \\ X_{B_z} \end{matrix}\right] = \left[\begin{matrix} \frac{n^x}{n^z} & 0 \\ \frac{n^y}{n^z} & 0 \\ 1 & 0 \\ 0 & \frac{n^x}{n^z} \\ 0 & \frac{n^y}{n^z} \\ 0 & 1 \end{matrix}\right] \left[\begin{matrix} u_1 \\ u_2 \end{matrix}\right]$

Let $u_1 \rightarrow u_1 n^z, u_2 \rightarrow u_2 n^z$

$\left[\begin{matrix} X_{D_x} \\ X_{D_y} \\ X_{D_z} \\ X_{B_x} \\ X_{B_y} \\ X_{B_z} \end{matrix}\right] = \left[\begin{matrix} n^x & 0 \\ n^y & 0 \\ n^z & 0 \\ 0 & n^x \\ 0 & n^y \\ 0 & n^z \end{matrix}\right] \left[\begin{matrix} u_1 \\ u_2 \end{matrix}\right]$

$\left[\begin{matrix} \vec{X}_D \\ \vec{X}_B \end{matrix}\right] = \left[\begin{matrix} \vec{n} & 0 \\ 0 & \vec{n} \\ \end{matrix}\right] \left[\begin{matrix} u_1 \\ u_2 \end{matrix}\right]$

Eigenvectors of $\lambda = -\frac{|n|}{\sqrt{\epsilon\mu}}$:

$\left[\begin{matrix} \frac{|n|}{\sqrt{\epsilon\mu}} & 0 & 0 & 0 & \frac{1}{\mu} n^z & -\frac{1}{\mu} n^y \\ 0 & \frac{|n|}{\sqrt{\epsilon\mu}} & 0 & -\frac{1}{\mu} n^z & 0 & \frac{1}{\mu} n^x \\ 0 & 0 & \frac{|n|}{\sqrt{\epsilon\mu}} & \frac{1}{\mu} n^y & -\frac{1}{\mu} n^x & 0 \\ 0 & -\frac{1}{\epsilon} n^z & \frac{1}{\epsilon} n^y & \frac{|n|}{\sqrt{\epsilon\mu}} & 0 & 0 \\ \frac{1}{\epsilon} n^z & 0 & -\frac{1}{\epsilon} n^x & 0 & \frac{|n|}{\sqrt{\epsilon\mu}} & 0 \\ -\frac{1}{\epsilon} n^y & \frac{1}{\epsilon} n^x & 0 & 0 & 0 & \frac{|n|}{\sqrt{\epsilon\mu}} \end{matrix}\right] \left[\begin{matrix} X_{D_x} \\ X_{D_y} \\ X_{D_z} \\ X_{B_x} \\ X_{B_y} \\ X_{B_z} \end{matrix}\right] = 0$

$\left[\begin{matrix} 1 & 0 & 0 & 0 & \frac{1}{|n|} \sqrt{\frac{\epsilon}{\mu}} n^z & -\frac{1}{|n|} \sqrt{\frac{\epsilon}{\mu}} n^y \\ 0 & \frac{|n|}{\sqrt{\epsilon\mu}} & 0 & -\frac{1}{\mu} n^z & 0 & \frac{1}{\mu} n^x \\ 0 & 0 & \frac{|n|}{\sqrt{\epsilon\mu}} & \frac{1}{\mu} n^y & -\frac{1}{\mu} n^x & 0 \\ 0 & -\frac{1}{\epsilon} n^z & \frac{1}{\epsilon} n^y & \frac{|n|}{\sqrt{\epsilon\mu}} & 0 & 0 \\ 0 & 0 & -\frac{1}{\epsilon} n^x & 0 & \frac{1}{\sqrt{\epsilon\mu}} (|n| - \frac{1}{|n|} (n^z)^2) & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} n^y n^z \\ 0 & \frac{1}{\epsilon} n^x & 0 & 0 & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} n^y n^z & \frac{1}{\sqrt{\epsilon\mu}} (|n| - \frac{1}{|n|} (n^y)^2) \end{matrix}\right] \left[\begin{matrix} X_{D_x} \\ X_{D_y} \\ X_{D_z} \\ X_{B_x} \\ X_{B_y} \\ X_{B_z} \end{matrix}\right] = 0$

$\left[\begin{matrix} |n| \sqrt{\frac{\mu}{\epsilon}} & 0 & 0 & 0 & n^z & -n^y \\ 0 & 1 & 0 & -\frac{1}{|n|} \sqrt{\frac{\epsilon}{\mu}} n^z & 0 & \frac{1}{|n|} \sqrt{\frac{\epsilon}{\mu}} n^x \\ 0 & 0 & \frac{|n|}{\sqrt{\epsilon\mu}} & \frac{1}{\mu} n^y & -\frac{1}{\mu} n^x & 0 \\ 0 & 0 & \frac{1}{\epsilon} n^y & \frac{1}{\sqrt{\epsilon\mu}} (|n| - \frac{1}{|n|} (n^z)^2) & 0 & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} n^x n^z \\ 0 & 0 & -\frac{1}{\epsilon} n^x & 0 & \frac{1}{\sqrt{\epsilon\mu}} (|n| - \frac{1}{|n|} (n^z)^2) & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} n^y n^z \\ 0 & 0 & 0 & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} n^x n^z & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} n^y n^z & \frac{1}{\sqrt{\epsilon\mu}} (|n| - \frac{1}{|n|} (n^x)^2 - \frac{1}{|n|} (n^y)^2) \end{matrix}\right] \left[\begin{matrix} X_{D_x} \\ X_{D_y} \\ X_{D_z} \\ X_{B_x} \\ X_{B_y} \\ X_{B_z} \end{matrix}\right] = 0$

$\left[\begin{matrix} |n| \sqrt{\frac{\mu}{\epsilon}} & 0 & 0 & 0 & n^z & -n^y \\ 0 & |n| \sqrt{\frac{\mu}{\epsilon}} & 0 & -n^z & 0 & n^x \\ 0 & 0 & 1 & \frac{1}{|n|} \sqrt{\frac{\epsilon}{\mu}} n^y & -\frac{1}{|n|} \sqrt{\frac{\epsilon}{\mu}} n^x & 0 \\ 0 & 0 & 0 & \frac{1}{\sqrt{\epsilon\mu}} (|n| - \frac{1}{|n|} (n^y)^2 - \frac{1}{|n|} (n^z)^2) & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} n^x n^y & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} n^x n^z \\ 0 & 0 & 0 & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} n^x n^y & \frac{1}{\sqrt{\epsilon\mu}} (|n| - \frac{1}{|n|} (n^x)^2 - \frac{1}{|n|} (n^z)^2) & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} n^y n^z \\ 0 & 0 & 0 & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} n^x n^z & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} n^y n^z & \frac{1}{\sqrt{\epsilon\mu}} (|n| - \frac{1}{|n|} (n^x)^2 - \frac{1}{|n|} (n^y)^2) \end{matrix}\right] \left[\begin{matrix} X_{D_x} \\ X_{D_y} \\ X_{D_z} \\ X_{B_x} \\ X_{B_y} \\ X_{B_z} \end{matrix}\right] = 0$

$\left[\begin{matrix} |n| \sqrt{\frac{\mu}{\epsilon}} & 0 & 0 & 0 & n^z & -n^y \\ 0 & |n| \sqrt{\frac{\mu}{\epsilon}} & 0 & -n^z & 0 & n^x \\ 0 & 0 & |n| \sqrt{\frac{\mu}{\epsilon}} & n^y & -n^x & 0 \\ 0 & 0 & 0 & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} (n^x)^2 & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} n^x n^y & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} n^x n^z \\ 0 & 0 & 0 & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} n^x n^y & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} (n^y)^2 & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} n^y n^z \\ 0 & 0 & 0 & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} n^x n^z & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} n^y n^z & \frac{1}{|n|} \frac{1}{\sqrt{\epsilon\mu}} (n^z)^2 \end{matrix}\right] \left[\begin{matrix} X_{D_x} \\ X_{D_y} \\ X_{D_z} \\ X_{B_x} \\ X_{B_y} \\ X_{B_z} \end{matrix}\right] = 0$

$\left[\begin{matrix} |n| \sqrt{\frac{\mu}{\epsilon}} & 0 & 0 & 0 & n^z & -n^y \\ 0 & |n| \sqrt{\frac{\mu}{\epsilon}} & 0 & -n^z & 0 & n^x \\ 0 & 0 & |n| \sqrt{\frac{\mu}{\epsilon}} & n^y & -n^x & 0 \\ 0 & 0 & 0 & n^x & n^y & n^z \\ 0 & 0 & 0 & n^x & n^y & n^z \\ 0 & 0 & 0 & n^x & n^y & n^z \end{matrix}\right] \left[\begin{matrix} X_{D_x} \\ X_{D_y} \\ X_{D_z} \\ X_{B_x} \\ X_{B_y} \\ X_{B_z} \end{matrix}\right] = 0$

$\left[\begin{matrix} |n| \sqrt{\frac{\mu}{\epsilon}} & 0 & 0 & 0 & n^z & -n^y \\ 0 & |n| \sqrt{\frac{\mu}{\epsilon}} & 0 & 0 & \frac{n^y n^z}{n^x} & n^x + \frac{(n^z)^2}{n^x} \\ 0 & 0 & |n| \sqrt{\frac{\mu}{\epsilon}} & 0 & -n^x - \frac{(n^y)^2}{n^x} & -\frac{n^y n^z}{n^x} \\ 0 & 0 & 0 & n^x & n^y & n^z \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{matrix}\right] \left[\begin{matrix} X_{D_x} \\ X_{D_y} \\ X_{D_z} \\ X_{B_x} \\ X_{B_y} \\ X_{B_z} \end{matrix}\right] = 0$

Let $u_2 = X_{B_y}, u_1 = X_{B_z}$

$|n| \sqrt{\frac{\mu}{\epsilon}} X_{D_x} = -n^z u_2 + n^y u_1$
$|n| \sqrt{\frac{\mu}{\epsilon}} X_{D_y} = -\frac{n^y n^z}{n^x} u_2 - (n^x + \frac{(n^z)^2}{n^x}) u_1$
$|n| \sqrt{\frac{\mu}{\epsilon}} X_{D_z} = (n^x + \frac{(n^y)^2}{n^x}) u_2 + \frac{n^y n^z}{n^x} u_1$
$n^x X_{B_x} = -n^y u_2 - n^z u_1$

$\left[\begin{matrix} X_{D_x} \\ X_{D_y} \\ X_{D_z} \\ X_{B_x} \\ X_{B_y} \\ X_{B_z} \end{matrix}\right] = \left[\begin{matrix} \frac{1}{|n|} \sqrt{\frac{\epsilon}{\mu}} n^y & -\frac{1}{|n|} \sqrt{\frac{\epsilon}{\mu}} n^z \\ -\frac{1}{|n|} \sqrt{\frac{\epsilon}{\mu}} (n^x + \frac{(n^z)^2}{n^x}) & -\frac{1}{|n|} \sqrt{\frac{\epsilon}{\mu}} \frac{n^y n^z}{n^x} \\ \frac{1}{|n|} \sqrt{\frac{\epsilon}{\mu}} \frac{n^y n^z}{n^x} & \frac{1}{|n|} \sqrt{\frac{\epsilon}{\mu}} (n^x + \frac{(n^y)^2}{n^x}) \\ -\frac{n^z}{n^x} & -\frac{n^y}{n^x} \\ 0 & 1 \\ 1 & 0 \end{matrix}\right] \left[\begin{matrix} u_1 \\ u_2 \end{matrix}\right]$

Let $u_1 \rightarrow \frac{|n|}{n^x} u_1, -u_2 \rightarrow \frac{|n|}{n^x} u_2$

$\left[\begin{matrix} X_{D_x} \\ X_{D_y} \\ X_{D_z} \\ X_{B_x} \\ X_{B_y} \\ X_{B_z} \end{matrix}\right] = \left[\begin{matrix} \sqrt{\frac{\epsilon}{\mu}} n^x n^y & \sqrt{\frac{\epsilon}{\mu}} n^x n^z \\ \sqrt{\frac{\epsilon}{\mu}} ((n^y)^2 - |n|^2) & \sqrt{\frac{\epsilon}{\mu}} n^y n^z \\ \sqrt{\frac{\epsilon}{\mu}} n^y n^z & \sqrt{\frac{\epsilon}{\mu}} ((n^z)^2 - |n|^2) \\ -|n| n^z & |n| n^y \\ 0 & -|n| n^x \\ |n| n^x & 0 \end{matrix}\right] \left[\begin{matrix} u_1 \\ u_2 \end{matrix}\right]$

Rewritten in vector form:

$\left[\begin{matrix} \vec{X}_D \\ \vec{X}_B \end{matrix}\right] = \left[\begin{matrix} \sqrt{\frac{\epsilon}{\mu}} ((\vec{n} \cdot \vec{e}_y) \vec{n} - (\vec{n} \cdot \vec{n}) \vec{e}_y) & \sqrt{\frac{\epsilon}{\mu}} ((\vec{n} \cdot \vec{e}_z) \vec{n} - (\vec{n} \cdot \vec{n}) \vec{e}_z) \\ |n| (\vec{n} \times \vec{e}_y) & |n| (\vec{n} \times \vec{e}_z) \end{matrix}\right] \left[\begin{matrix} u_1 \\ u_2 \end{matrix}\right]$

$\left[\begin{matrix} \vec{X}_D \\ \vec{X}_B \end{matrix}\right] = \left[\begin{matrix} \sqrt{\frac{\epsilon}{\mu}} \vec{n} \times (\vec{n} \times \vec{e}_y) & \sqrt{\frac{\epsilon}{\mu}} \vec{n} \times (\vec{n} \times \vec{e}_z) \\ |n| (\vec{n} \times \vec{e}_y) & |n| (\vec{n} \times \vec{e}_z) \end{matrix}\right] \left[\begin{matrix} u_1 \\ u_2 \end{matrix}\right]$

TODO notice that I'm missing the metric from all of this.
Also notice below that I'm pulling the metric out on the left hand of the flux Jacobian.
If I don't pull the metric out then I bet I will find it in the inner product of n above.



generalization of the explicitly stated terms below:

$R = \left[\begin{matrix} \sqrt{\epsilon} n_2 & \sqrt{\epsilon} n_3 & \sqrt{\epsilon} n & 0 & \sqrt{\epsilon} n_2 & -\sqrt{\epsilon} n_3 \\ -\sqrt{\mu} n_3 & \sqrt{\mu} n_2 & 0 & \sqrt{\mu} n & \sqrt{\mu} n_3 & -\sqrt{\mu} n_2 \end{matrix}\right]$

$L = \left[\begin{matrix} n_2 & -n_3 \\ n_3 & n_2 \\ n & 0 \\ 0 & n \\ n_2 & n_3 \\ n_3 & -n_2 \end{matrix}\right]$



for the x-direction:
$R = \left[\matrix{ 0 & 0 & \sqrt{\epsilon} & 0 & 0 & 0 \\ \sqrt{\epsilon} & 0 & 0 & 0 & \sqrt{\epsilon} & 0 \\ 0 & \sqrt{\epsilon} & 0 & 0 & 0 & -\sqrt{\epsilon} \\ 0 & 0 & 0 & \sqrt{\mu} & 0 & 0 \\ 0 & \sqrt{\mu} & 0 & 0 & 0 & -\sqrt{\mu} \\ -\sqrt{\mu} & 0 & 0 & 0 & \sqrt{\mu} & 0 }\right]$
...Trangenstein says...
$R = \left[\matrix{ 0 & 0 & -\sqrt{\epsilon} & \sqrt{\epsilon} & 0 & 0 \\ 0 & -\sqrt{\epsilon} & 0 & 0 & \sqrt{\epsilon} & 0 \\ -\sqrt{\epsilon} & 0 & 0 & 0 & 0 & -\sqrt{\epsilon} \\ 0 & 0 & \sqrt{\mu} & \sqrt{\mu} & 0 & 0 \\ \sqrt{\mu} & 0 & 0 & 0 & 0 & \sqrt{\mu} \\ 0 & \sqrt{\mu} & 0 & 0 & \sqrt{\mu} & 0 }\right]$

$L = \left[\matrix{ 0 & 1 & 0 & 0 & 0 & -1 \\ 0 & 0 & 1 & 0 & 1 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 & -1 &0 }\right]$

$\Lambda = \left[\matrix{ -v_p & 0 & 0 & 0 & 0 & 0 \\ 0 & -v_p & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & v_p & 0 \\ 0 & 0 & 0 & 0 & 0 & v_p }\right]$

curvilinear coordinates using vectors:

conserved quantities:
$U^I = \left[\matrix{ D^i \\ B^i }\right]$

Flux:

$F^{Ij} = \left[\matrix{ -{\bar\epsilon^{ij}}_k H^k \\ {\bar\epsilon^{ij}}_k E^k }\right]$
$= \left[\matrix{ -\frac{1}{\sqrt{g}} \bar\epsilon^{ijk} g_{kl} H^l \\ \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} g_{kl} E^l }\right]$

Flux partial derivative divergence:

${F^{Ij}}_{,j} = \left[\matrix{ -\frac{1}{\sqrt{g}} \bar\epsilon^{ijk} g_{kl} H^l \\ \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} g_{kl} E^l }\right]_{,j}$
$= \left[\matrix{ -\bar\epsilon^{ijk} ((\frac{1}{\sqrt{g}})_{,j} g_{kl} H^l + \frac{1}{\sqrt{g}} g_{kl,j} H^l + \frac{1}{\sqrt{g}} g_{kl} {H^l}_{,j}) \\ \bar\epsilon^{ijk} ((\frac{1}{\sqrt{g}})_{,j} g_{kl} E^l + \frac{1}{\sqrt{g}} g_{kl,j} E^l + \frac{1}{\sqrt{g}} g_{kl} {E^l}_{,j}) }\right]$
$= \left[\matrix{ -\frac{1}{\sqrt{g}} \bar\epsilon^{ijk} ( -\frac{1}{\sqrt{g}} (\sqrt{g})_{,j} H_k + g_{kl,j} H^l + g_{kl} {H^l}_{,j}) \\ \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} ( -\frac{1}{\sqrt{g}} (\sqrt{g})_{,j} E_k + g_{kl,j} E^l + g_{kl} {E^l}_{,j}) }\right]$
$= \left[\matrix{ -\frac{1}{\sqrt{g}} \bar\epsilon^{ijk} g_{kl} {H^l}_{,j} \\ \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} g_{kl} {E^l}_{,j} }\right] + \left[\matrix{ -\frac{1}{\sqrt{g}} \bar\epsilon^{ijk} (-\frac{1}{\sqrt{g}} (\sqrt{g})_{,j} H_k + g_{kl,j} H^l) \\ \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} (-\frac{1}{\sqrt{g}} (\sqrt{g})_{,j} E_k + g_{kl,j} E^l) }\right]$

flux covariant derivative:
${F^{Ij}}_{;j} = \left[\matrix{ \bar\epsilon^{ijk} \frac{1}{\sqrt{g}} g_{kl} {H^l}_{;j} \\ -\bar\epsilon^{ijk} \frac{1}{\sqrt{g}} g_{kl} {E^l}_{;j} }\right]$
$= \left[\matrix{ \bar\epsilon^{ijk} \frac{1}{\sqrt{g}} g_{kl} ({H^l}_{,j} + {\Gamma^l}_{mj} H^m) \\ -\bar\epsilon^{ijk} \frac{1}{\sqrt{g}} g_{kl} ({E^l}_{,j} + {\Gamma^l}_{mj} E^m) }\right]$
$= \left[\matrix{ \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} g_{kl} {H^l}_{,j} \\ -\frac{1}{\sqrt{g}} \bar\epsilon^{ijk} g_{kl} {E^l}_{,j} }\right] + \left[\matrix{ -\frac{1}{\sqrt{g}} \bar\epsilon^{ijk} \Gamma_{jkl} H^l \\ \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} \Gamma_{jkl} E^l }\right]$

Add and subtract the flux weighted by the volume gradient and the metric partial terms.
This is what the code in my curvilinear solver does by default. But maybe this should be enabled on a per-equation-system basis. It seems to be unnecessary for the Maxwell equations.

${F^{Ij}}_{;j} = \left[\matrix{ \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} g_{kl} {H^l}_{,j} \\ -\frac{1}{\sqrt{g}} \bar\epsilon^{ijk} g_{kl} {E^l}_{,j} }\right] + \left[\matrix{ -\frac{1}{\sqrt{g}} \bar\epsilon^{ijk} (-\frac{1}{\sqrt{g}} (\sqrt{g})_{,j} H_k + g_{kl,j} H^l) \\ \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} (-\frac{1}{\sqrt{g}} (\sqrt{g})_{,j} E_k + g_{kl,j} E^l) }\right] - \left[\matrix{ -\frac{1}{\sqrt{g}} \bar\epsilon^{ijk} (-\frac{1}{\sqrt{g}} (\sqrt{g})_{,j} H_k + g_{kl,j} H^l) \\ \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} (-\frac{1}{\sqrt{g}} (\sqrt{g})_{,j} E_k + g_{kl,j} E^l) }\right] + \frac{\sqrt{g}_{,j}}{\sqrt{g}} \left[\matrix{ -\frac{1}{\sqrt{g}} \bar\epsilon^{ijk} g_{kl} H^l \\ \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} g_{kl} E^l }\right] - \frac{\sqrt{g}_{,j}}{\sqrt{g}} \left[\matrix{ -\frac{1}{\sqrt{g}} \bar\epsilon^{ijk} g_{kl} H^l \\ \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} g_{kl} E^l }\right] + \left[\matrix{ -\frac{1}{\sqrt{g}} \bar\epsilon^{ijk} \Gamma_{jkl} H^l \\ \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} \Gamma_{jkl} E^l }\right]$
$ = \frac{\sqrt{g}}{\sqrt{g}} {F^{Ij}}_{,j} + \frac{\sqrt{g}_{,j}}{\sqrt{g}} F^{Ij} + \left[\matrix{ -\frac{1}{\sqrt{g}} \bar\epsilon^{ijk} \Gamma_{jkl} H^l \\ \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} \Gamma_{jkl} E^l }\right] - \left[\matrix{ -\frac{1}{\sqrt{g}} \bar\epsilon^{ijk} (-{\Gamma^m}_{jm} H_k + (\Gamma_{klj} + \Gamma_{lkj}) H^l) \\ \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} (-{\Gamma^m}_{jm} E_k + (\Gamma_{klj} + \Gamma_{lkj}) E^l) }\right] - \left[\matrix{ -{\Gamma^m}_{jm} \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} g_{kl} H^l \\ {\Gamma^m}_{jm} \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} g_{kl} E^l }\right] $
$ = \frac{1}{\sqrt{g}} \left( \sqrt{g} F^{Ij} \right)_{,j} + \left[\matrix{ -2 \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} \Gamma_{jkl} H^l \\ 2 \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} \Gamma_{jkl} E^l }\right] $

hyperbolic conservation law system using vectors:

${U^i}_{,t} + {F^{ij}}_{;j} = S^i$
${U^i}_{,t} + \frac{1}{\sqrt{g}} \left( \sqrt{g} F^{ij} \right)_{,j} = S^i + \left[\matrix{ 2 \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} \Gamma_{jkl} H^l \\ -2 \frac{1}{\sqrt{g}} \bar\epsilon^{ijk} \Gamma_{jkl} E^l }\right]$

curvilinear coordinates using one-forms:

$\left[\matrix{ D_i \\ B_i }\right] + \left[\matrix{ -{\epsilon_i}^{jk} H_k \\ {\epsilon_i}^{jk} E_k }\right]_{;j} = \left[\matrix{ -J_i \\ 0 }\right]$

Raise all indexes of the Levi-Civita tensor. Replace it with the associated Levi-Civita tensor of flat space times the inverse (square root optional) metric determinant.

$\left[\matrix{ D_i \\ B_i }\right] + \left[\matrix{ -\frac{1}{\sqrt{g}} g_{il} \bar\epsilon^{ljk} H_k \\ \frac{1}{\sqrt{g}} g_{il} \bar\epsilon^{ljk} E_k }\right]_{;j} = \left[\matrix{ -J_i \\ 0 }\right]$

Factor out metric terms since their covariant derivative value is zero.

$\left[\matrix{ D_i \\ B_i }\right] + \frac{1}{\sqrt{g}} \left[\matrix{ g_{il} & 0 \\ 0 & g_{il} }\right] \left[\matrix{ -\bar\epsilon^{ljk} H_k \\ \bar\epsilon^{ljk} E_k }\right]_{;j} = \left[\matrix{ -J_i \\ 0 }\right]$

Apply covariant derivative. Replace one-forms with partial derivatives and connection products.

$\left[\matrix{ D_i \\ B_i }\right] + \frac{1}{\sqrt{g}} \left[\matrix{ g_{il} & 0 \\ 0 & g_{il} }\right] \left[\matrix{ -\bar\epsilon^{ljk} (H_{k,j} - {\Gamma^m}_{kj} H_m) \\ \bar\epsilon^{ljk} (E_{k,j} - {\Gamma^m}_{kj} E_m) }\right] = \left[\matrix{ -J_i \\ 0 }\right]$

Multiply the antisymmetric Levi-Civita tensor with the symmetric (only in holonomic torsion-free case) connection coefficients to cancel them.
Factor out the remaining partial derivative terms.

$\left[\matrix{ D_i \\ B_i }\right] + \frac{1}{\sqrt{g}} \left[\matrix{ g_{il} & 0 \\ 0 & g_{il} }\right] \underbrace{ \left[\matrix{ 0 & -\frac{1}{\mu} \bar\epsilon^{ljk} \\ \frac{1}{\epsilon} \bar\epsilon^{ljk} & 0 }\right] \left[\matrix{ D_k \\ B_k }\right]_{,j} }{} = \left[\matrix{ -J_i \\ 0 }\right]$

Now the underbraced terms have an eigen-decomposition. We can compute the flux using this, and transform it afterwards using the metric.

with GLM variables

removing divergence as a PDE:
$D_{i,t} = -D_{i,i} = -\phi_{,i}$ for $\phi_{,t} = -E_{i,i}$
$B_{i,t} = -B_{i,i} = -\psi_{,i}$ for $\psi_{,t} = -B_{i,i}$

combined with the original system:
(from http://www.ammar-hakim.org/maxwell-eigensystem.html)
Notice that I assume $\epsilon, \mu$ are constant, otherwise I would have to pick out their derivatives (as I am doing with the curvilinear coordinate system).
$B_{i,t} + {\epsilon_i}^{jk} E_{k,j} + \chi_\psi \psi_{,i} = 0$
$\mu \epsilon E_{i,t} - {\epsilon_i}^{jk} B_{k,j} + \chi_\phi \phi_{,i} = -\mu J_i$
$\frac{1}{\chi_\phi} \phi_{,t} + E_{i,i} = \epsilon^{-1} \rho_{free}$
$\epsilon \mu \frac{1}{\chi_\psi} \psi_{,t} + B_{i,i} = 0$

Rewritten in terms of $E \rightarrow D = \epsilon E$
$\mu D_{i,t} - {\epsilon_i}^{jk} B_{k,j} + \chi_\phi \phi_{,i} = -\mu J_i$
$B_{i,t} + \epsilon^{-1} {\epsilon_i}^{jk} D_{k,j} + \chi_\psi \psi_{,i} = 0$
$\frac{1}{\chi_\phi} \phi_{,t} + \epsilon^{-1} D_{i,i} = \epsilon^{-1} \rho_{free}$
$\epsilon \mu \frac{1}{\chi_\psi} \psi_{,t} + B_{i,i} = 0$

Scale line #1 by $\mu^{-1}$, line #3 by $\chi_\phi$, line #4 by $\epsilon^{-1} \mu^{-1}$
$D_{i,t} - \mu^{-1} {\epsilon_i}^{jk} B_{k,j} + \mu^{-1} \chi_\phi \phi_{,i} = -J_i$
$B_{i,t} + \epsilon^{-1} {\epsilon_i}^{jk} D_{k,j} + \chi_\psi \psi_{,i} = 0$
$\phi_{,t} + \epsilon^{-1} \chi_\phi D_{i,i} = \epsilon^{-1} \chi_\phi \rho_{free}$
$\psi_{,t} + \epsilon^{-1} \mu^{-1} \chi_\psi B_{i,i} = 0$

change of variables $\phi \rightarrow \phi' = \mu \phi$
$D_{i,t} - \mu^{-1} {\epsilon_i}^{jk} B_{k,j} + \chi_\phi \phi_{,i} = -J_i$
$B_{i,t} + \epsilon^{-1} {\epsilon_i}^{jk} D_{k,j} + \chi_\psi \psi_{,i} = 0$
$\mu \phi_{,t} + \epsilon^{-1} \chi_\phi D_{i,i} = \epsilon^{-1} \chi_\phi \rho_{free}$
$\psi_{,t} + \epsilon^{-1} \mu^{-1} \chi_\psi B_{i,i} = 0$

substitute $v_p^2 = \epsilon^{-1} \mu^{-1}$
$D_{i,t} - \mu^{-1} {\epsilon_i}^{jk} B_{k,j} + \chi_\phi \phi_{,i} = -J_i$
$B_{i,t} + \epsilon^{-1} {\epsilon_i}^{jk} D_{k,j} + \chi_\psi \psi_{,i} = 0$
$\phi_{,t} + (v_p)^2 \chi_\phi D_{i,i} = (v_p)^2 \chi_\phi \rho_{free}$
$\psi_{,t} + (v_p)^2 \chi_\psi B_{i,i} = 0$

scale $\phi \rightarrow \phi' = v_p \phi, \psi \rightarrow \psi' = v_p \psi$, $\chi \rightarrow \chi' = \frac{1}{v_p} \chi$
$D_{i,t} - \mu^{-1} {\epsilon_i}^{jk} B_{k,j} + \chi_\phi \phi_{,i} = -J_i$
$B_{i,t} + \epsilon^{-1} {\epsilon_i}^{jk} D_{k,j} + \chi_\psi \psi_{,i} = 0$
$v_p \phi_{,t} + v_p \chi_\phi D_{i,i} = v_p \chi_\phi \rho_{free}$
$v_p \psi_{,t} + v_p \chi_\psi B_{i,i} = 0$

Scale line #3 and line #4 by $(v_p)^{-1}$
$D_{i,t} - \mu^{-1} {\epsilon_i}^{jk} B_{k,j} + \chi_\phi \phi_{,i} = -J_i$
$B_{i,t} + \epsilon^{-1} {\epsilon_i}^{jk} D_{k,j} + \chi_\psi \psi_{,i} = 0$
$\phi_{,t} + \chi_\phi D_{i,i} = \chi_\phi \rho_{free}$
$\psi_{,t} + \chi_\psi B_{i,i} = 0$

conserved quantities:
$U^I = \left[\matrix{ D^i \\ B^i \\ \phi \\ \psi }\right]$

Sources:
https://en.wikipedia.org/wiki/Electric_displacement_field
Trangenstein, "Numeric Solutions of Hyperbolic Partial Differential Equations," 2009