 Solution

Preserve the primary constraints using the corrected Hamiltonian from solution. After a spatial integration by parts, the $B_0$ terms in its density are $B_0{\partial }_i{\pi }_i+B_0^{\ast }{\partial }_i{\pi }_i^{\ast }-M^2{B}_0^{\ast }{B}_0$. Hamilton's equations giveWith the definitions of $E^i,E^{\ast i}$ in the paper, the actual massive Gauss law is thereforeThis is also the $\nu =0$ component of the Proca equation. The zero right-hand sides printed in the source apply only to the massless theory or to special configurations with $B_0=0$.

For a concrete counterexample, take $k\ne 0$, ${\omega }^2=k^2+M^2$, and the longitudinal on-shell plane waveIt satisfies ${\partial }_{\mu }{B}^{\mu }=0$, but $E^3=-iM{e}^{-i\omega t+ikz}$ and ${\partial }_i{E}^i=kM{e}^{-i\omega t+ikz}\ne 0$. Imposing the printed zero-divergence condition would incorrectly remove the physical longitudinal particle polarization.