Solution (source code)

= Solution

Since
$$
\alpha n=\partial_t-\beta
$$
and $j_\mu n^\mu=0$, the scaling rule for a covector Lie derivative has no extra term:
$$
\alpha\mathcal L_nj_\mu
=\mathcal L_{\partial_t-\beta}j_\mu
=\partial_tj_\mu-\mathcal L_\beta j_\mu.
$$
The acceleration is the spatial lapse gradient,
$$
a_\mu=D_\mu\log\alpha=\frac{D_\mu\alpha}{\alpha}.
$$
Multiplying the momentum equation by $\alpha$ therefore gives
$$
\boxed{
\partial_tj_\mu
=\mathcal L_\beta j_\mu
-\alpha D_\nu S^\nu{}_\mu
-S^\nu{}_\mu D_\nu\alpha
+\alpha Kj_\mu-\rho D_\mu\alpha}.
$$
For spatial components, the shift term can be written
$$
\mathcal L_\beta j_i
=\beta^kD_kj_i+j_kD_i\beta^k,
$$
so explicitly
$$
\boxed{
\partial_tj_i
=\beta^kD_kj_i+j_kD_i\beta^k
-\alpha D_kS^k{}_i-S^k{}_iD_k\alpha
+\alpha Kj_i-\rho D_i\alpha}.
$$