Maxwell Gauss constraints propagate by form transport

Put $U:=dF$ and $V:=d{\star }_{g}F$. They are continuous three-forms. By [F1], ${\iota }_{n}U={\iota }_{n}V=0$ throughout the evolution. In a local frame with timelike vector $n$, a three-form annihilated by $n$ has only its spatial volume component, so there are unique continuous scalars $q,r$ with $U=q{\iota }_n{\mathrm{v}\mathrm{o}\mathrm{l}}_g$ and $V=r{\iota }_n{\mathrm{v}\mathrm{o}\mathrm{l}}_g$. On ${\mathrm{\Sigma }}_0$, the spatial components of $dF$ and $d{\star }_{g}F$ are respectively $(D_i{b}^i){\mathrm{v}\mathrm{o}\mathrm{l}}_h$ and, up to the fixed orientation sign, $(D_i{e}^i){\mathrm{v}\mathrm{o}\mathrm{l}}_h$. Thus $q=r=0$ on the initial slice.

Since $dU=d^{2}F=0$ and $dV=d^2{\star }_{g}F=0$ distributionally, the coordinate formula for exterior derivative gives ${\mathrm{\partial }}_{\mu }(\sqrt{\mathrm{\mid }\det g\mathrm{\mid }}q{n}^{\mu })=0$ and the same equation for $r$. The identity ${\mathrm{\partial }}_{\mu }(\sqrt{\mathrm{\mid }\det g\mathrm{\mid }}{n}^{\mu })=\sqrt{\mathrm{\mid }\det g\mathrm{\mid }}{\mathrm{\nabla }}_{\mu }{n}^{\mu }$ follows from ${\mathrm{\Gamma }}^{\mu }{}_{\mu \nu }={\mathrm{\partial }}_{\nu }\log \sqrt{\mathrm{\mid }\det g\mathrm{\mid }}$, itself obtained by differentiating the determinant and using the Levi-Civita formula. Therefore $n(q)+q{\mathrm{div}}_{g}n=0$ and $n(r)+r{\mathrm{div}}_{g}n=0$ in the distributional sense.

The field $n$ is $C^1$ and transverse to ${\mathrm{\Sigma }}_0$ by hypothesis. Its local flow map $\mathrm{\Phi }(\tau ,x)$ is $C^1$ and the differential at $\tau =0$ maps the time direction to $n$ and the slice directions identically, hence is invertible. By [F3], after shrinking, $(\tau ,x)$ are flow coordinates. Pullback preserves distributional closedness of these continuous three-forms under the $C^1$ flow chart: test against compactly supported smooth functions, use the change-of-variables formula, and approximate the resulting compactly supported $C^1$ test functions in their $C^1$ norm. In these coordinates $n={\mathrm{\partial }}_{\tau }$ and ${\mathrm{v}\mathrm{o}\mathrm{l}}_g=J(\tau ,x)d\tau \wedge d{x}^1\wedge d{x}^2\wedge d{x}^3$ with continuous $J>0$. Thus $dU=0$ reads ${\mathrm{\partial }}_{\tau }(Jq)=0$ distributionally, and $dV=0$ reads ${\mathrm{\partial }}_{\tau }(Jr)=0$. A continuous function whose distributional derivative in $\tau$ vanishes is constant in $\tau$: convolving in $\tau$ gives smooth functions with zero derivative, and uniform convergence on compact subcylinders passes that constancy to the original function; testing against spatial test functions first gives the same conclusion when both variables are present. Their initial values vanish by step 1.1, so $q=r=0$ throughout the flow neighborhood.

The identities $U=q{\iota }_n{\mathrm{v}\mathrm{o}\mathrm{l}}_g$ and $V=r{\iota }_n{\mathrm{v}\mathrm{o}\mathrm{l}}_g$ now give $dF=0$ and $d{\star }_{g}F=0$. The argument uses only the Maxwell evolution, its initial Gauss constraints, and $d^2=0$.