Solution (source code)

= Solution

Use <suffix notation>, with repeated indices summed. Contracting the two <Levi-Civita symbols> gives
$$
[\mathbf u\times(\nabla\times\mathbf A)]_i
=\epsilon_{ijk}u_j\epsilon_{klm}\partial_lA_m
=u_j(\partial_iA_j-\partial_jA_i).
$$
Applying the <product rule> to $h=A_iB_i$ and substituting the two evolution equations yields
$$
\begin{aligned}
\partial_t h
&=B_iu_j\partial_iA_j-B_iu_j\partial_jA_i+B_i\partial_i\psi
+A_iB_j\partial_ju_i-A_iu_j\partial_jB_i\\
&=B_i\partial_i(u_jA_j+\psi)-u_j\partial_j(A_iB_i).
\end{aligned}
$$
In the second line we renamed dummy indices in $A_iB_j\partial_ju_i$ and then combined product derivatives. Hence
$$
\boxed{f=\mathbf u\cdot\mathbf A+\psi,\qquad
\partial_t h=\mathbf B\cdot\nabla f-\mathbf u\cdot\nabla h.}
$$
No zero-divergence assumption was used in this local identity.

When both fields are <solenoidal vector fields>, the <product rule for divergence> rewrites this as a <conservation law>:
$$
\partial_t h=\nabla\cdot(f\mathbf B-h\mathbf u).
$$
Because the volume is fixed, differentiation under its integral and the <divergence theorem> give
$$
\frac{dH}{dt}=\int_{\partial V}(f\mathbf B-h\mathbf u)\cdot\mathbf n\,dS=0
$$
under the two tangency conditions. Thus \b[the volume integral is conserved]. This is <helicity conservation by tangent boundary conditions>; when $\mathbf B=\nabla\times\mathbf A$, the density is the <magnetic helicity> density.

For the given axisymmetric field, the <curl in spherical coordinates> gives
$$
\begin{aligned}
B_r&=\frac1{r\sin\theta}\partial_\theta(\sin\theta A_\phi)
=2(a^2-r^2)\cos\theta,\\
B_\theta&=-\frac1r\partial_r(rA_\phi)
=(4r^2-2a^2)\sin\theta,\\
B_\phi&=\frac1r\partial_r(rA_\theta)=3ar\sin\theta.
\end{aligned}
$$
There are no $\phi$ derivatives and $A_r=0$. Consequently
$$
\begin{aligned}
h&=A_\theta B_\theta+A_\phi B_\phi\\
&=ar^2(4r^2-2a^2)\sin^2\theta+3ar^2(a^2-r^2)\sin^2\theta\\
&=\boxed{ar^2(a^2+r^2)\sin^2\theta.}
\end{aligned}
$$
The <volume form> in <spherical polar coordinates> is $r^2\sin\theta\,dr\,d\theta\,d\phi$, so
$$
\begin{aligned}
H&=2\pi a\int_0^a r^4(a^2+r^2)\,dr\int_0^\pi\sin^3\theta\,d\theta\\
&=\frac{8\pi a}{3}\left(\frac{a^7}{5}+\frac{a^7}{7}\right)
=\boxed{\frac{32\pi a^8}{35}.}
\end{aligned}
$$
Here $B_r=0$ at $r=a$, as required for the magnetic-field tangency condition in the conservation argument; tangency of $\mathbf u$ is a separate hypothesis.