Following the paper:

Alcubierre, Degollado, Salgado, "The Einstein-Maxwell System in 3+1 form and initial data for multiple charged black holes" 2009.

variables:
$\mu_0 =$ vacuum permeability
$\epsilon_0 =$ vacuum permissivity
$c = \frac{1}{\sqrt{\mu_0 \epsilon_0}} =$ vacuum speed-of-light
$\rho =$ density
$v^i = $ 3-velocity
$u^a = \left[u^t, u^i\right] = \left[u^t, u^t v^i\right] = $ 4-velocity
$u^t = (1 - v^2)^{-1/2} =$ Lorentz factor
$p^a = \rho u^a =$ 4-momentum

Special relativistic basis:
$\partial_t = \frac{1}{c} \frac{\partial}{\partial t}$
$\partial^t = c dt$
$\partial_i = \frac{\partial}{\partial x^i}$
$\partial^i = dx^i$

Special relativistic metric:
$\eta_{ab} = \left[\matrix{ -1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 }\right]$
$\eta^{ab} = \left[\matrix{ -1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 }\right]$

Faraday tensor:
$\phi = $ electric potential
$A^i = $ magnetic vector potential
$A^u = \left[ -\frac{1}{c} \phi, A^i \right] $ electromagnetic 4-potential
$A_u = \left[ \frac{1}{c} \phi, A_i \right] $
$F_{uv} = 2 \partial_{[u} A_{v]}$ in flat space
$F_{uv} $$ = \downarrow u \overset{\rightarrow v}{ \left[\matrix{ 0 & \partial_t A_j - \partial_j A_t \\ \partial_i A_t - \partial_t A_i & \partial_i A_j - \partial_j A_i }\right] }$$ = \downarrow u \overset{\rightarrow v}{ \left[\matrix{ 0 & \frac{1}{c} (\frac{\partial}{\partial t} A_j - \partial_j \phi) \\ \frac{1}{c} (\partial_i \phi - \frac{\partial}{\partial t} A_i) & \partial_i A_j - \partial_j A_i }\right] }$$ = \downarrow u \overset{\rightarrow v}{ \left[\matrix{ 0 & -\frac{1}{c} E_j \\ \frac{1}{c} E_i & \epsilon_{ijk} B^k }\right] }$
${F^u}_v = \downarrow u \overset{\rightarrow v}{ \left[\matrix{ 0 & \frac{1}{c} E_j \\ \frac{1}{c} E^i & {\epsilon^i}_{jk} B^k }\right] }$

$F^{uv} = \downarrow u \overset{\rightarrow v}{ \left[\matrix{ 0 & \frac{1}{c} E^j \\ -\frac{1}{c} E^i & \epsilon^{ijk} B_k }\right] }$
$(\star F)_{uv} = \downarrow u \overset{\rightarrow v}{ \left[\matrix{ 0 & -B_j \\ B_i & \frac{1}{c} \epsilon_{ijk} E^k }\right] }$
$(\star F)^{uv} = \downarrow u \overset{\rightarrow v}{ \left[\matrix{ 0 & -B^j \\ B^i & \frac{1}{c} \epsilon^{ijk} E_k }\right] }$

Lorentz force law:
$\dot{p}^a = q {F^a}_b u^b$
$\frac{1}{c} \frac{\partial}{\partial t} (\rho u^t v^i) = q u^t (\frac{1}{c} E^i + {\epsilon^i}_{jk} B^k v^j)$
in flat space:
$\frac{1}{c} \frac{\partial}{\partial t} (\rho v^i) = q (\frac{1}{c} E^i + {\epsilon^i}_{jk} B^k v^j)$

Maxwell in flat space:
$\partial_b (\star F)^{ab} = \partial_b (\frac{1}{2} F^{uv} {\epsilon_{uv}}^{ab}) = 0$
$\partial_i B^i = 0$
$\frac{1}{c} \frac{\partial}{\partial t} B^i + \frac{1}{c} {\epsilon_i}^{jk} \partial_j E_k = 0$
$\frac{\partial}{\partial t} B^i + {\epsilon_i}^{jk} \partial_j E_k = 0$

$\sigma = $ electric conductivity
$j^i = \sigma E^i = $ current density
$J^a = \left[ c \rho_{charge}, j^i \right] =$ 4-current
$\partial_b F^{ab} = \mu_0 J^a$
$\frac{1}{c} \partial_i E^i = \mu_0 c \rho_{charge}$
$\mu_0 \epsilon_0 \partial_i E^i = \mu_0 \rho_{charge}$
$\partial_i E^i = \frac{1}{\epsilon_0} \rho_{charge}$

$-\frac{1}{c^2} \frac{\partial}{\partial t} E^i + \epsilon^{ijk} \partial_j B_k = \mu_0 j^i$
$-\mu_0 \epsilon_0 \frac{\partial}{\partial t} E^i + \epsilon^{ijk} \partial_j B_k = \mu_0 j^i$
$\epsilon_0 \frac{\partial}{\partial t} E^i - \epsilon^{ijk} \partial_j (\frac{1}{\mu_0} B)_k = -j^i$

Curved space

$\alpha =$ lapse
$\beta^i =$ shift
$\gamma_{ij} =$ spatial metric
$g_{uv} = \downarrow u \overset{\rightarrow v}{\left[\matrix{ -\alpha^2 + \beta^2 & \beta_j \\ \beta_i & \gamma_{ij} }\right]}$
$g^{uv} = \downarrow u \overset{\rightarrow v}{\left[\matrix{ -\alpha^{-2} & \alpha^{-2} \beta^j \\ \alpha^{-2} \beta^i & \gamma^{ij} - \alpha^{-2} \beta^i \beta^j }\right]}$

$n_a = (-\alpha, 0)$ and $n^a = (1 / \alpha, -\beta^i / \alpha)$.

$K_{ab} = -{\gamma_a}^c {\gamma_b}^d \nabla_{(c} n_{d)}$
Eqn 2.10 in the paper says
$K = g^{ab} K_{ab} = -g^{ab} {\gamma_a}^c {\gamma_b}^d \nabla_{(c} n_{d)}$
$= -g^{ab} \nabla_a n_b$
$= -\nabla_a n^a$
Note that $K$ is the 4D trace, not the 3D trace that it is in numerical relativity texts. (Right? those extra $g^{tt}$ and $g^{ti}$ terms don't all eliminate, do they?)
Some math later, hopefully without any errors...
$K = \frac{1}{2 \alpha^2} (\partial_t \alpha + 2 \alpha \partial_i \beta^i - \beta^i \partial_i \alpha - \frac{\alpha}{\gamma} ( \partial_t \gamma - \beta^i \partial_i \gamma) )$
So the $c$'s appear inside $K$ on a per-term basis...

$\nabla_u F^{au} = 2 \nabla_u \nabla^{[a} A^{u]} + {R^a}_u A^u$ in curved space

$^{(3)} \epsilon^{\perp afg} = n_b {\gamma^a}_e {\gamma^f}_c {\gamma^g}_d \epsilon^{bacd} = n_b \epsilon^{bafg}$
$\epsilon^{0123} = \frac{-1}{\sqrt{-g}} = \frac{-1}{\alpha \sqrt\gamma}$
$^{(3)} \epsilon^{\perp afg} = -\alpha {\gamma^a}_e {\gamma^f}_c {\gamma^g}_d \epsilon^{0ecd}$
$^{(3)} \epsilon^{\perp 123} = -\alpha \epsilon^{0123} = \frac{1}{\sqrt\gamma}$
$^{(3)} \epsilon^{abc} = ^{(3)} \epsilon^{\perp abc}$
such that $^{(3)} \epsilon^{123} = \frac{1}{\sqrt\gamma}$
also note $\epsilon^{abc}_F = \sqrt\gamma ^{(3)} \epsilon^{abc}$ and $\epsilon^F_{abc} = \frac{1}{\sqrt\gamma} ^{(3)} \epsilon_{abc}$
for $\epsilon^{abc}_F$ and $\epsilon^F_{abc}$ the flat-space undensitized permutation tensor (typically denoted as $[abc]$ right?)

Faraday tensor projection:
$F^{ab} = ^{(3)} F^{ab} + n^a {}^{(3)} F^{\perp b} + {}^{(3)} F^{a \perp} n^b$

Measure electric and magnetic fields by Eulerian observers:
$E^a = -n_b F^{ba} = ^{(3)}F^{\perp a}$
$B^a = -n_b (\star F)^{ba} = ^{(3)} (\star F)^{\perp a}$

Note $n_a E^a = 0$ and $n_a B^a = 0$, therefore $E^a$ and $B^a$ are vectors in the spacelike hypersurface.

Combine the definitions to find:
$F^{ab} = ^{(3)} F^{ab} + n^a E^b - E^a n^b$

Combine further to define $B^a$:
$B^a = \frac{1}{2} n_b {\epsilon^{ba}}_{cd} F^{cd}$
$= \frac{1}{2} n_b \epsilon^{bacd} (^{(3)} F_{cd} n_c E_d - E_c n_d)$
$= \frac{1}{2} n_b \epsilon^{bacd} {}^{(3)} F_{cd}$
Jumping to the end:
$\frac{1}{c} D_a E^a = \frac{1}{c} \frac{1}{\sqrt\gamma} \partial_i (\sqrt\gamma E^i) = \mu_0 c \rho$ for $\rho = -n_a j^a$
$D_a B^a = \frac{1}{\sqrt\gamma} \partial_i (\sqrt\gamma B^i) = 0$

using $\frac{d}{dt} = \frac{1}{c} \frac{\partial}{\partial t} - \mathcal{L}_\vec\beta$
$\frac{1}{c} \frac{d}{dt} E^i = ^{(3)} \epsilon^{ijk} \partial_j (\alpha B_k) + \frac{1}{c} \alpha K E^i - \mu_0 \alpha j^i$
$\frac{d}{dt} B^i = -\frac{1}{c} ^{(3)} \epsilon^{ijk} \partial_j (\alpha E_k) + \alpha K B^i$

with flat-space permutation tensors, expanded partials, in hyperbolic conservation form:
$\frac{1}{c^2} \frac{\partial}{\partial t} E^i - \frac{1}{c} \beta^j \partial_j E^i - \frac{1}{\sqrt\gamma} \alpha \epsilon_F^{ijk} \gamma_{kl} \partial_j B^l = \frac{1}{\sqrt\gamma} \epsilon_F^{ijk} (\partial_j \alpha B_k + \alpha \partial_j \gamma_{kl} B^l) + \frac{1}{c} \alpha K E^i - \frac{1}{c} E^j \partial_j \beta^i - \mu_0 \alpha (j^i + \beta^i j^t) $
$\frac{1}{c} \frac{\partial}{\partial t} B^i + \frac{1}{c} \frac{1}{\sqrt\gamma} \alpha \epsilon_F^{ijk} \gamma_{kl} \partial_j E^l - \partial_j B^i \beta^j = -\frac{1}{c} \frac{1}{\sqrt\gamma} \epsilon_F^{ijk} (\partial_j \alpha E_k + \alpha \partial_j \gamma_{kl} E^l) + \alpha K B^i - B^j \partial_j \beta^i $

fix units:
$\frac{\partial}{\partial t} (\epsilon_0 E^i) - c \beta^j \partial_j (\epsilon_0 E^i) - \frac{1}{\mu_0} \frac{1}{\sqrt\gamma} \alpha \epsilon_F^{ijk} \gamma_{kl} \partial_j B^l = \frac{1}{\mu_0} \frac{1}{\sqrt\gamma} \epsilon_F^{ijk} (\partial_j \alpha B_k + \alpha \partial_j \gamma_{kl} B^l) + c \alpha K \epsilon_0 E^i - c \epsilon_0 E^j \partial_j \beta^i - \alpha (j^i + \beta^i j^t) $
$\frac{\partial}{\partial t} B^i + \frac{1}{\sqrt\gamma} \alpha \epsilon_F^{ijk} \gamma_{kl} \partial_j E^l - c \partial_j B^i \beta^j = -\frac{1}{\epsilon_0} \frac{1}{\sqrt\gamma} \epsilon_F^{ijk} (\partial_j \alpha \epsilon_0 E_k + \alpha \partial_j \gamma_{kl} \epsilon_0 E^l) + c \alpha K B^i - c B^j \partial_j \beta^i $

as a matrix:
$\frac{\partial}{\partial t} \left[\matrix{ \epsilon_0 E^i \\ B^i }\right] + \partial_j \left[\matrix{ - c \beta^j \delta^i_l & -\frac{1}{\mu_0} \frac{\alpha}{\sqrt\gamma} \epsilon_F^{ijk} \gamma_{kl} \\ \frac{1}{\epsilon_0} \frac{\alpha}{\sqrt\gamma} \epsilon_F^{ijk} \gamma_{kl} & - c \beta^j \delta^i_l }\right] \left[\matrix{ \epsilon_0 E^l \\ B^l }\right] = \left[\matrix{ \frac{1}{\mu_0} \frac{1}{\sqrt\gamma} \epsilon_F^{ijk} (\partial_j \alpha B_k + \alpha \partial_j \gamma_{kl} B^l) + c \epsilon_0 \alpha K E^i - c \epsilon_0 E^j \partial_j \beta^i - \alpha (j^i + \beta^i j^t) \\ -\frac{1}{\epsilon_0} \frac{1}{\sqrt\gamma} \epsilon_F^{ijk} (\partial_j \alpha \epsilon_0 E_k + \alpha \partial_j \gamma_{kl} \epsilon_0 E^l) + c \alpha K B^i - c B^j \partial_j \beta^i }\right]$

I'm leaving the $\partial \alpha$ and $\partial \gamma_{ij}$ terms on the source side because they are replaced with state variables when converting the system to 1st order.

charge conservation:
$\nabla_c j^c = 0$
in 3+1 form:
$\frac{1}{c} \frac{\partial}{\partial t} \rho - \mathcal{L}_\vec\beta \rho = -D_a (\alpha j^a) + \alpha \rho K$
$\frac{1}{c} \frac{\partial}{\partial t} \rho - \beta^i \partial_i \rho - \rho \partial_i \beta^i = -\partial_i (\alpha j^i) - {\Gamma^i}_{ki} j^k + \alpha \rho K$

But what does the system look like if we rewrite it in terms of $\alpha E^i$ and $\alpha B^i$?
That should make some of the Levi-Civita tensor source terms disappear.
$\frac{d}{dt} (\alpha E^i) = E^i \frac{d}{dt} \alpha + \alpha \frac{d}{dt} E^i$
$\frac{d}{dt} E^i = \frac{1}{\alpha} (\frac{d}{dt} (\alpha E^i) - E^i \frac{d}{dt} \alpha)$

$\frac{1}{c} \frac{1}{\alpha} (\frac{d}{dt} (\alpha E^i) - E^i \frac{d}{dt} \alpha) = ^{(3)} \epsilon^{ijk} \partial_j (\alpha B_k) + \frac{1}{c} \alpha K E^i - \mu_0 \alpha j^i$
$\frac{1}{\alpha} (\frac{d}{dt} (\alpha B^i) - B^i \frac{d}{dt} \alpha) = -\frac{1}{c} ^{(3)} \epsilon^{ijk} \partial_j (\alpha E_k) + \alpha K B^i$

$\partial_t (\alpha E^i) - \beta^j \delta^i_l \partial_j (\alpha E^l) - c \alpha ^{(3)} \epsilon^{ijk} \partial_j (\alpha B^l \gamma_{kl}) = \alpha K (\alpha E^i) - c \mu_0 \alpha^2 \cdot j^i + \frac{1}{\alpha} (\alpha E^i) ( \partial_t \alpha - \beta^i \partial_i \alpha - \alpha \partial_i \beta^i ) - (\alpha E^j) \partial_j \beta^i $
$\partial_t (\alpha B^i) - \beta^j \delta^i_l \partial_j (\alpha B^l) + \frac{1}{c} \alpha ^{(3)} \epsilon^{ijk} \partial_j (\alpha E^l \gamma_{kl}) = \alpha K (\alpha B^i) + \frac{1}{\alpha} (\alpha B^i) ( \partial_t \alpha - \beta^i \partial_i \alpha - \alpha \partial_i \beta^i ) - (\alpha B^j) \partial_j \beta^i $

Hmm, there's still the $\epsilon^{ijk} (\alpha E^l \gamma_{kl})$ ...
So there will always be an extra $\epsilon^{ijk}$ in the source term...
Unless we make the state vars equal to $\alpha E_i$ lowered ...
Then I'll have to add a $\frac{d}{dt} \gamma_{ij}$ to the source term ...
...which is at least equal to $-2 \alpha K_{ij}$, which isn't so bad ...

$\partial_t (\alpha E^i) = \partial_t (\alpha E_j \gamma^{ij})$
$ = \partial_t (\alpha E_j) \gamma^{ij} + \alpha E_j \partial_t \gamma^{ij}$
$ = \partial_t (\alpha E_j) \gamma^{ij} - \alpha E^k \gamma^{ij} \partial_t \gamma_{jk}$
$ = \partial_t (\alpha E_j) \gamma^{ij} + 2 \alpha K^{ij} (\alpha E_j)$

Hmm but we still have the $\beta^i$'s ...
How about just re-derive $\frac{d}{dt} (\alpha E_i)$...
$= \partial_t (\alpha E_i) - \beta^j \partial_j (\alpha E_i) - (\alpha E_j) \partial_i \beta^j$

But really, whether I use the upper or lower vector form, either has a flux term of its dual, which means a Levi-Civita times its opposite, which means the partial of the gamma will always be there. So there's no point in rewriting it in covariant form over contravariant form.

$\partial_t (\alpha E_j) \gamma^{ij} - \beta^j \delta^i_l \partial_j (\alpha E^l) - c \alpha ^{(3)} \epsilon^{ijk} \partial_j (\alpha B_k) = \alpha K (\alpha E^i) + \frac{1}{\alpha} (\alpha E^i) ( \partial_t \alpha - \beta^i \partial_i \alpha - \alpha \partial_i \beta^i ) - (\alpha E^j) \partial_j \beta^i - 2 \alpha K^{ij} (\alpha E_j) - c \mu_0 \alpha^2 j^i $
$\partial_t (\alpha B_j) \gamma^{ij} - \beta^j \delta^i_l \partial_j (\alpha B^l) + \frac{1}{c} \alpha ^{(3)} \epsilon^{ijk} \partial_j (\alpha E_k) = \alpha K (\alpha B^i) + \frac{1}{\alpha} (\alpha B^i) ( \partial_t \alpha - \beta^i \partial_i \alpha - \alpha \partial_i \beta^i ) - (\alpha B^j) \partial_j \beta^i - 2 \alpha K^{ij} (\alpha B_j) $