Solution (source code)

= Solution

Let $\boldsymbol\sigma$ be the stated <Maxwell stress tensor>. Taking its divergence and using
$$
(\mathbf A\mathbin\cdot\nabla)\mathbf A
-\frac12\nabla|\mathbf A|^2
=-\mathbf A\times(\nabla\times\mathbf A)
$$
gives
$$
\begin{aligned}
\partial_j\sigma_{ij}
={}&\epsilon_0[\mathbf E\times(\nabla\times\mathbf E)]_i
-\epsilon_0E_i\nabla\mathbin\cdot\mathbf E\\
&+\frac1{\mu_0}[\mathbf B\times(\nabla\times\mathbf B)]_i
-\frac1{\mu_0}B_i\nabla\mathbin\cdot\mathbf B.
\end{aligned}
$$
The <Maxwell equations> in vacuum with sources are
$$
\nabla\mathbin\cdot\mathbf E=\frac\rho{\epsilon_0},
\qquad
\nabla\mathbin\cdot\mathbf B=0,
\qquad
\nabla\times\mathbf E=-\partial_t\mathbf B,
\qquad
\nabla\times\mathbf B=\mu_0\mathbf J+\mu_0\epsilon_0\partial_t\mathbf E.
$$
Substitution, together with antisymmetry of the <cross product>, yields
$$
\begin{aligned}
\partial_j\sigma_{ij}
&=-\rho E_i-(\mathbf J\times\mathbf B)_i
-\epsilon_0[\mathbf E\times\partial_t\mathbf B
+\partial_t\mathbf E\times\mathbf B]_i\\
&=-[\rho\mathbf E+\mathbf J\times\mathbf B]_i
-\partial_t[\epsilon_0\mathbf E\times\mathbf B]_i.
\end{aligned}
$$
Therefore the required vector is the <electromagnetic momentum density>
$$
\boxed{\mathbf g=\epsilon_0\mathbf E\times\mathbf B
=\frac{\mathbf S}{c^2}},
$$
and
$$
\boxed{\sum_{j=1}^3\frac{\partial\sigma_{ij}}{\partial x_j}
+\frac{\partial g_i}{\partial t}
=-(\rho\mathbf E+\mathbf J\times\mathbf B)_i}.
$$
This is <local conservation of electromagnetic momentum>. The tensor component $\sigma_{ij}$ is the flux of $i$-momentum in the $j$ direction in the convention of the question, $\mathbf g$ is field momentum per unit volume, and $\rho\mathbf E+\mathbf J\times\mathbf B$ is the <Lorentz force density> transferred to matter.