\begin{equation} \partial\cdot j_A=\frac{2}{(4 \pi )^{d/2} (d/2)!}\varepsilon^{\mu_1\mu_2\cdots\mu_d}\text{tr}\left(tF_{\mu_1\mu_2}\cdots F_{\mu_{d-1}\mu_d}\right) \end{equation}\begin{equation} \partial\cdot j^5_a=\frac{2}{(4 \pi )^{d/2} (d/2)!}\varepsilon^{\mu_1\mu_2\cdots\mu_d}\text{tr}\left(tt_aF_{\mu_1\mu_2}\cdots F_{\mu_{d-1}\mu_d}\right) \end{equation}
which is nothing but the Euler class associated to the gauge field. Here $t$ is a certain matrix in flavour space, and $a$ is an index in colour space (with $t_a$ the generator of the representation under which the fermions transform).
To simplify the notation, we use DeWitt generalised indices, where an index is actually a pair of indices, the first one ranging over the spacetime manifold, and the second one over a finite-dimensional vector space (a fibre). Contraction of indices corresponds to both integration and summation. Moreover, when there is no ambiguity, all contracted indices will be left implicit. For example, the classical action will be denoted by \begin{equation}\tag{1} S=\bar\psi\not D\psi:=\int_{\mathbb R^{2d}} \bar\psi_\alpha(x) (\not D)^\alpha{}_\beta(x,y) \psi^\beta(y)\ \mathrm dx\,\mathrm dy \end{equation} where $\alpha,\beta$ are indices over flavour, colour and spin, and $\not D$ is the (operator-valued) matrix \begin{equation}\tag{2} \not D(x,y):=-\gamma^\mu\left(1_\mathrm{flavour}\oplus1_\mathrm{colour} \partial_\mu+ A_\mu(x)\right)\delta(x-y) \end{equation}\begin{equation}\tag{2} \not D(x,y):=-\gamma^\mu1_\mathrm{flavour}\left(1_\mathrm{colour} \partial_\mu+ A_\mu(x)\right)\delta(x-y) \end{equation} with $A(x)$ the (external) Yang-Mills field.
Let \begin{equation}\tag{3} Z=\int \mathrm e^{-\bar\psi\not D\psi}\ \mathrm d\psi\,\mathrm d\bar\psi \end{equation}\begin{equation}\tag{3} Z[A]=\int \mathrm e^{-\bar\psi\not D\psi}\ \mathrm d\psi\,\mathrm d\bar\psi \end{equation} be the partition function. The strategy to derive an expression for the axial anomaly is to study how the integrand responds to the axial transformation.
Recall that the axial transformation is a global symmetry of the classical action. Therefore, if we let the gauge parameterparameters $\epsilon$$\epsilon^a$ become a function of spacetime, the change in the action must be a total derivative (cf. this PSE post): \begin{equation}\tag{5} S\to S+\epsilon\partial\cdot j_A \end{equation}\begin{equation}\tag{5} S\to S+\epsilon^a\partial\cdot j^5_a \end{equation} and this relation defines the axial current $j_A$$j^5_a$. In the classical theory, this implies that the current is conserved, $\partial\cdot j_A\equiv 0$$\partial\cdot j_a^5\equiv 0$.
If the integration measure $\mathrm d\psi\,\mathrm d\bar\psi$ were invariant as well, we would conclude that the current is also conserved in the quantum theory. We shall argue below that the measure is actually not invariant. We shall find that \begin{equation}\tag{6} \mathrm d\psi\,\mathrm d\bar\psi\to \mathrm e^{\epsilon\cdot \mathcal A}\mathrm d\psi\,\mathrm d\bar\psi \end{equation}\begin{equation}\tag{6} \mathrm d\psi\,\mathrm d\bar\psi\to \mathrm e^{\epsilon^a\mathcal A_a}\mathrm d\psi\,\mathrm d\bar\psi \end{equation} where $\mathcal A=\mathcal A(x)$$\mathcal A=\mathcal A_a(x)$ is the so-called anomaly function.
As the axial transformation is nothing but a change of variables, the partition function remains unchanged, so \begin{equation} \begin{aligned} 0&=\int \mathrm e^{-\bar\psi\not D\psi}\left(\mathrm e^{-\epsilon\partial\cdot j_A+\epsilon\cdot\mathcal A}-1\right)\ \mathrm d\psi\,\mathrm d\bar\psi\\ &=\epsilon \left\langle -\partial\cdot j_A+\mathcal A\right\rangle+\mathcal O(\epsilon^2) \end{aligned}\tag{7} \end{equation}\begin{equation} \begin{aligned} 0&=\int \mathrm e^{-\bar\psi\not D\psi}\left(\mathrm e^{-\epsilon\partial\cdot j^5+\epsilon\mathcal A}-1\right)\ \mathrm d\psi\,\mathrm d\bar\psi\\ &=\epsilon \left\langle -\partial\cdot j^5+\mathcal A\right\rangle+\mathcal O(\epsilon^2) \end{aligned}\tag{7} \end{equation} from which we obtain an explicit formula for the anomaly: \begin{equation}\tag{8} \partial\cdot j_A\equiv\mathcal A \end{equation}\begin{equation}\tag{8} \langle\partial\cdot j^5_a\rangle\equiv\mathcal A_a \end{equation}
We now calculate the change in the integration measure under $\psi\to U\psi$, where \begin{equation}\tag{9} U(x,y)=\delta(x-y)\mathrm e^{t\gamma_5\epsilon(x)} \end{equation}\begin{equation}\tag{9} U(x,y)=\delta(x-y)\mathrm e^{tt_a\gamma_5\epsilon^a(x)} \end{equation} is the chiral transformation matrix, with $t$ an hermitian matrix in flavour space. In practice one is usually interested in the case where different flavours transform independently, i.e., $t=0$$t=1$. We leave this matrix arbitrary to include a more general situation.
As $\psi,\bar\psi$ are fermionic coordinates, the Jacobian appears in the denominator: \begin{equation} \begin{aligned} \mathrm d\psi\,\mathrm d\bar\psi&\to (\det U\det\bar U)^{-1}\mathrm d\psi\,\mathrm d\bar\psi\\ &=\mathrm e^{2\epsilon\,\text{tr}(\delta t\gamma_5)}\mathrm d\psi\,\mathrm d\bar\psi \end{aligned}\tag{10} \end{equation}\begin{equation} \begin{aligned} \mathrm d\psi\,\mathrm d\bar\psi&\to (\det U\det\bar U)^{-1}\mathrm d\psi\,\mathrm d\bar\psi\\ &=\mathrm e^{2\epsilon^a\,\text{tr}(\delta tt_a\gamma_5)}\mathrm d\psi\,\mathrm d\bar\psi \end{aligned}\tag{10} \end{equation} where we have used that $\bar U=U$ and Jacobi's formula $\det \mathrm e^X=\mathrm e^{\text{tr}X}$.
From this we learn that \begin{equation} \begin{aligned} \mathcal A(x)&=2\text{tr}(\delta t\gamma_5)\\ &=2\delta(0)\text{tr}(t\gamma_5) \end{aligned}\tag{11} \end{equation}\begin{equation} \begin{aligned} \mathcal A_a(x)&=2\text{tr}(\delta tt_a\gamma_5)\\ &=2\delta(0)\text{tr}(tt_a\gamma_5) \end{aligned}\tag{11} \end{equation} which is, apparently, both independent of the gauge field $A$ and the spacetime point $x$. This is, in fact, not really the case.
These conditions can be obtained by reverse-engineering the final result and making sure the anomaly function remains finite in the $\Lambda\to\infty$ limit (which is a sensible requirement). For the record, Fujikawa chose $f(s)=\mathrm e^{-s^2}$$f(s)=\mathrm e^{s}$. We leave $f$ arbitrary to stress the fact that the result is, to some extent, regulator-independent, as long as it satisfies the conditions above. We could consider even weaker conditions on $f$, but they would lead to some finite renormalisations in the final result, which are not really inadmissible but they carry no relevant information for our purposes.
With this, and changing into momentum space, the anomaly becomes \begin{equation}\tag{14} \mathcal A=2\lim_{\Lambda\to\infty}\Lambda^d\int_{\mathbb R^d} \text{tr}\left(t\gamma_5f\left[-(i\not k+\not D/\Lambda)^2\right]\right)\frac{\mathrm dk}{(2\pi)^d} \end{equation}\begin{equation}\tag{14} \mathcal A_a=2\lim_{\Lambda\to\infty}\Lambda^d\int_{\mathbb R^d} \text{tr}\left(tt_a\gamma_5f\left[-(i\not k+\not D/\Lambda)^2\right]\right)\frac{\mathrm dk}{(2\pi)^d} \end{equation}
Noting that \begin{equation}\tag{16} \text{tr}(\gamma_5\gamma^{\mu_1}\gamma^{\mu_2}\cdots\gamma^{\mu_d})=2^{d/2}\varepsilon^{\mu_1\mu_2\cdots\mu_d} \end{equation} we get \begin{equation} \begin{aligned} \text{tr}\left(t\gamma_5\not D^d\right)&=2^{d/2}\varepsilon^{\mu_1\mu_2\cdots\mu_d}\text{tr}\left(tD_{\mu_1}D_ {\mu_2}\cdots D_{\mu_d}\right)\\ &=\varepsilon^{\mu_1\mu_2\cdots\mu_d}\text{tr}\left(t[D_{\mu_1},D_ {\mu_2}]\cdots [D_{\mu_{d-1}},D_{\mu_d}]\right)\\ &=\varepsilon^{\mu_1\mu_2\cdots\mu_d}\text{tr}\left(tF_{\mu_1\mu_2}\cdots F_{\mu_{d-1}\mu_d}\right) \end{aligned}\tag{17} \end{equation}\begin{equation} \begin{aligned} \text{tr}\left(tt_a\gamma_5\not D^d\right)&=2^{d/2}\varepsilon^{\mu_1\mu_2\cdots\mu_d}\text{tr}\left(tt_aD_{\mu_1}D_ {\mu_2}\cdots D_{\mu_d}\right)\\ &=\varepsilon^{\mu_1\mu_2\cdots\mu_d}\text{tr}\left(tt_a[D_{\mu_1},D_ {\mu_2}]\cdots [D_{\mu_{d-1}},D_{\mu_d}]\right)\\ &=\varepsilon^{\mu_1\mu_2\cdots\mu_d}\text{tr}\left(tt_aF_{\mu_1\mu_2}\cdots F_{\mu_{d-1}\mu_d}\right) \end{aligned}\tag{17} \end{equation}
Finally, using \begin{equation} \begin{aligned} \int_{\mathbb R^d}f^{(d/2)}(k^2)\frac{\mathrm dk}{(2\pi)^d}&=\frac{1}{(2\pi)^d}\frac{2\pi^{d/2}}{\left(\frac d2-1\right)!}\int_{\mathbb R^+} k^{d-1}f^{(d/2)}(k^2)\ \mathrm dk\\ &\overset{\mathrm{IBP}}=\frac{1}{(4\pi)^{d/2}}(f(0)-f(\infty)) \end{aligned}\tag{18} \end{equation} we get \begin{equation}\tag{19} \partial\cdot j_A=\frac{2}{(4 \pi )^{d/2} (d/2)!}\varepsilon^{\mu_1\mu_2\cdots\mu_d}\text{tr}\left(tF_{\mu_1\mu_2}\cdots F_{\mu_{d-1}\mu_d}\right) \end{equation}\begin{equation}\tag{19} \partial\cdot j_a^5=\frac{2}{(4 \pi )^{d/2} (d/2)!}\varepsilon^{\mu_1\mu_2\cdots\mu_d}\text{tr}\left(tt_aF_{\mu_1\mu_2}\cdots F_{\mu_{d-1}\mu_d}\right) \end{equation} as promised.
For $d=4$ we recover the standard expression, and for $d=2$ we get \begin{equation}\tag{20} \partial\cdot j_A\bigg|_{d=2}=\frac{1}{2\pi}\varepsilon^{\mu_1\mu_2}\text{tr}\left(tF_{\mu_1\mu_2}\right) \end{equation}\begin{equation}\tag{20} \partial\cdot j_a^5\bigg|_{d=2}=\frac{1}{2\pi}\varepsilon^{\mu_1\mu_2}\text{tr}\left(tt_aF_{\mu_1\mu_2}\right) \end{equation} as requiredOP suspected.
If the chiral transformation is trivial in flavour space, we have $t=1_\mathrm{flavour}$, and its trace is just $N_\mathrm{flavour}$, the number of flavours.
In the operator formalism, the operators $\psi,\bar\psi$ satisfy the Euler-Lagrange equations of $S$, that is, \begin{equation}\tag{21} \overset{\rightarrow}{\not D}\psi=0=\bar\psi\overset{\leftarrow}{\not D} \end{equation}
Moreover, the axial current is formally defined as \begin{equation}\tag{21} j^\mu_A\sim \bar\psi\gamma^\mu \gamma_5t\psi \end{equation}\begin{equation}\tag{22} j^\mu_a(x)\sim \bar\psi(x)\gamma^\mu \gamma_5tt_a\psi(x) \end{equation}
If we formally take the divergence of this current, and use the equations of motion, we get $\partial\cdot j_A\sim 0$$\partial\cdot j_a^5\sim 0$, so it seems that the axial current is conserved. This argument fails because $\bar\psi(x) M\psi(y)$ diverges as $y\to x$, and therefore the operator $j_A$$j_a^5$ defined above is meaningless as written. Indeed, if we let $D(x)=\frac{\Gamma\left(\frac{d-2}{2}\right)}{4\pi^{d/2}}(1/x^2)^{\frac{d-2}{2}}$ be the massless bosonic propagator, the OPE of two fermionic fields is easily seen to be \begin{equation} \begin{aligned} \bar\psi(x)M\psi(y)&=\text{tr}(M\!\not\!\partial D(x-y))+\cdots\\ &=\frac{\Gamma(d/2)}{2\pi^{d/2}}\left(\frac{1}{(x-y)^2}\right)^{\frac{d-2}{2}}\text{tr}\left(\frac{\not\!\epsilon M}{\epsilon^2}\right)+\cdots \end{aligned}\tag{22} \end{equation}\begin{equation} \begin{aligned} \bar\psi(x)M\psi(y)&=\text{tr}(M\!\not\!\partial D(x-y))+\cdots\\ &=\frac{\Gamma(d/2)}{2\pi^{d/2}}\left(\frac{1}{(x-y)^2}\right)^{\frac{d-2}{2}}\text{tr}\left(\frac{(\not\!x-\not\!y) M}{(x-y)^2}\right)+\cdots \end{aligned}\tag{23} \end{equation} where $M$ is an arbitrary matrix, and $\cdots$ denote terms of lower order in $x-y$. This bilinear form is clearly divergent as $y\to x$, and therefore the current $j_A$$j_a^5$ defined above is meaningless as written.
The short distance divergence is easily avoided by smearing out the coincidence points, $\bar\psi(x)\psi(x+\epsilon)$, for some $\epsilon\in\mathbb R^d$ (unrelated to the gauge parameter $\epsilon^a$). As we are dealing with gauge-dependent objects, we must introduce a Wilson line so as to be able to compare them at different space-time points. In other words, in order to obtain a gauge-invariant current we are lead to define \begin{equation}\tag{23} j_A^\mu\equiv \lim_{\epsilon\to0} \bar\psi(x+\epsilon/2)\gamma^\mu\gamma_5t\, \mathrm{Pexp}\left[\int_{x-\epsilon/2}^{x+\epsilon/2} A\cdot\mathrm ds\right]\psi(x-\epsilon/2) \end{equation}\begin{equation}\tag{24} j_a^\mu\equiv \lim_{\epsilon\to0} \bar\psi(x+\epsilon/2)\gamma^\mu\gamma_5tt_a\, \mathrm{Pexp}\left[\int_{x-\epsilon/2}^{x+\epsilon/2} A(s)\cdot\mathrm ds\right]\psi(x-\epsilon/2) \end{equation} where $\mathrm{Pexp}$ is the path-ordered exponential. The limit $\epsilon\to0$ is highly formal, and is defined by the so-called symmetric prescription (cf. this PSE post) \begin{equation}\tag{24} \lim_{\epsilon\to0}\frac{\epsilon^{\mu_1}\cdots\epsilon^{\mu_n}}{(\epsilon^2)^m}\equiv \frac{\Gamma(d/2)}{2^m\Gamma(m+d/2)}\begin{cases} \delta^{(\mu_1\mu_2}\cdots\delta^{\mu_{n-1}\mu_n)} & n=2m\\0&\text{otherwise}\end{cases} \end{equation}\begin{equation}\tag{25} \lim_{\epsilon\to0}\frac{\epsilon^{\mu_1}\cdots\epsilon^{\mu_n}}{(\epsilon^2)^m}\equiv \frac{\Gamma(d/2)}{2^m\Gamma(m+d/2)}\begin{cases} \delta^{(\mu_1\mu_2}\cdots\delta^{\mu_{n-1}\mu_n)} & n=2m\\0&\text{otherwise}\end{cases} \end{equation} where the parentheses denote symmetrisation (without the sometimes conventional $1/n!$ factor). The $n\neq 2m$ case is fixed by dimensional analysis, and the $n=2m$ case is enforced by requiring that both sides agree when we contract both sides withcontracted with $\delta_{\mu_1\mu_2}\cdots\delta_{\mu_{n-1}\mu_n}$.
With this, and taking the divergence of $j_A$ and$j_a$, collecting the terms that survive the $\epsilon\to0$ limit, and using the known value of the OPE $\bar\psi(y)M\psi(x)$ and the symmetric prescription for the limit $\epsilon\to0$, we should recover the expression for the axial anomaly. This is a messy calculation that we shall not reproduce here, although we stress that all we need to do is to use the equations of motion together with equations $(22),(24)$$(21)-(25)$. This is, in principle, a straightforward computation.
