Solution (source code)

= Solution

For a smooth map $\mathbf u$, the <Jacobian determinant> is the determinant of its <Jacobian matrix> $A_{ai}=\partial_i u_a$. In <Einstein summation convention>,
$$
J[\mathbf u]=\det A=\frac16\epsilon_{ijk}\epsilon_{abc}(\partial_i u_a)(\partial_j u_b)(\partial_k u_c).
$$
Differentiate the supplied expression for $V_i$ by the <product rule>. The term differentiating $\partial_j u_a$ vanishes because $\partial_i\partial_j u_a$ is symmetric in $i,j$, while $\epsilon_{ijk}$ is antisymmetric; the term differentiating $\partial_k u_b$ vanishes similarly. Only differentiation of $u_c$ remains:
$$
\partial_iV_i=\frac16\epsilon_{ijk}\epsilon_{abc}(\partial_j u_a)(\partial_k u_b)(\partial_i u_c)=J[\mathbf u].
$$
The last equality follows by the cyclic relabelling $(i,j,k)\mapsto(j,k,i)$, which has positive sign. This proves the <Jacobian determinant as a divergence> identity.

For a composition $\mathbf w=\mathbf u\circ\mathbf v$, the <chain rule> gives
$$
\partial_iw_a(\mathbf x)=(\partial_bu_a)(\mathbf v(\mathbf x))\,\partial_iv_b(\mathbf x),\qquad D\mathbf w(\mathbf x)=D\mathbf u(\mathbf v(\mathbf x))D\mathbf v(\mathbf x).
$$
Taking <determinants> yields
$$
\boxed{J[\mathbf w](\mathbf x)=J[\mathbf u](\mathbf v(\mathbf x))J[\mathbf v](\mathbf x).}
$$
For the first-order perturbation, use <spatial derivative control in perturbations of the identity map>: the usual smooth-family interpretation gives $D\mathbf u_t=I+tD\mathbf F+o(t)$ locally. Multilinearity of the <determinant> shows that its linear term comes from choosing one column of $tD\mathbf F$ and all other columns from $I$. This gives the <trace>:
$$
\boxed{J[\mathbf u_t]=1+t\operatorname{tr}(D\mathbf F)+o(t)=1+t\nabla\cdot\mathbf F+o(t),\qquad Q=\nabla\cdot\mathbf F.}
$$
The regularity qualification matters: a pointwise $o(t)$ remainder for maps does not by itself give an $o(t)$ remainder for spatial derivatives. For instance $(x+t^2\sin(x/t^2),y,z)$ has a uniform $o(t)$ displacement and $\mathbf F=0$, but its <Jacobian determinant> at $x=0$ is $2$ for $t\ne0$. Thus the first expansion needs the remainder in $C^1$ locally, or an equivalent smoothness assumption on the family.

Under the subsequent group law, $\mathbf u_0$ is the identity and $\mathbf u_{-t}$ is the inverse of $\mathbf u_t$. Write $j(t,\mathbf x)=J[\mathbf u_t](\mathbf x)$. Applying the composition formula to $\mathbf u_{t+h}=\mathbf u_h\circ\mathbf u_t$ and the short-time expansion gives the <Jacobian evolution of a smooth flow>:
$$
\partial_tj(t,\mathbf x)=(\nabla\cdot\mathbf F)(\mathbf u_t(\mathbf x))j(t,\mathbf x),\qquad j(0,\mathbf x)=1.
$$
Hence
$$
j(t,\mathbf x)=\exp\left(\int_0^t(\nabla\cdot\mathbf F)(\mathbf u_s(\mathbf x))\,ds\right).
$$
If only the pointwise generator expansion is assumed initially, the group law still gives $\partial_t\mathbf u_t(\mathbf x)=\mathbf F(\mathbf u_t(\mathbf x))$. Uniqueness and smooth dependence for this smooth <ordinary differential equation> make it a smooth <flow map>, thereby justifying the differentiated expansion used above. Thus \b[a divergence-free generator has]
$$
\boxed{J[\mathbf u_t]=1\quad\text{for every }t.}
$$
By the <change of variables formula>, these maps preserve ordinary Euclidean volume and orientation. This is the Euclidean instance of a <volume-preserving vector field>.