Solution (source code)

= Solution

For a massless Dirac fermion,
$$
\mathcal L=\bar\psi i\gamma^\mu D_\mu\psi,
\qquad D_\mu=\partial_\mu+ieA_\mu.
$$
The local gauge symmetry gives the vector current $J^\mu=e\bar\psi\gamma^\mu\psi$. The global <chiral transformation> $\psi\mapsto e^{i\beta\gamma^5}\psi$ gives the classical <axial current>
$$
J_{\rm ax}^\mu=\bar\psi\gamma^\mu\gamma^5\psi.
$$
Using the massless <Dirac equation> and $\{\gamma^5,\gamma^\mu\}=0$ gives $\partial_\mu J_{\rm ax}^\mu=0$ classically.

Quantum mechanically, a gauge-invariant regulator for the fermion measure is not invariant under the axial rotation. In the Fujikawa form, its infinitesimal Jacobian contains
$$
-2i\beta\lim_{\Lambda\to\infty}
\operatorname{tr}\left[\gamma^5e^{-\not D^2/\Lambda^2}\right].
$$
The first nonzero term in the <heat kernel expansion> is quadratic in $F_{\mu\nu}$. Using
$$
\operatorname{tr}\{\gamma^5[\gamma^\mu,\gamma^\nu]
[\gamma^\rho,\gamma^\sigma]\}=16i\epsilon^{\mu\nu\rho\sigma}
$$
and the Gaussian momentum integral gives the <chiral anomaly>
$$
\boxed{\partial_\mu J_{\rm ax}^\mu
=-\frac{e^2}{16\pi^2}\epsilon^{\alpha\beta\gamma\delta}
F_{\alpha\beta}F_{\gamma\delta}}.
$$
The overall sign follows the gamma-matrix, charge, and Levi-Civita conventions stated in the question. The same coefficient is obtained from the one-loop axial-vector-vector triangle diagram.