The axial anomaly.
In this answer we calculate the axial anomaly, following Fujikawa's method, in an arbitrary number of spacetime dimensions $d$. We assume $d\in 2\mathbb N$ inasmuch as there is no $\gamma_5$ for odd $d$ (and thus no anomlay). We will mainly follow Ref.1.
OP is interested in the particular case $d=2$. This case is analysed, by means of a different method, in Ref.2, which OP might find useful.
The setup.
We work in euclidean spacetime to avoid cumbersome factors of $i$ here and there.
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} 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} 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.
The anomaly function.
Recall that the axial transformation is a global symmetry of the classical action. Therefore, if we let the gauge parameter $\epsilon$ 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} and this relation defines the axial current $j_A$. In the classical theory, this implies that the current is conserved, $\partial\cdot j_A\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} where $\mathcal A=\mathcal 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} from which we obtain an explicit formula for the anomaly: \begin{equation}\tag{8} \partial\cdot j_A\equiv\mathcal A \end{equation}
With this, all we have to do is to calculate $\mathcal A(x)$. To this end, we need to calculate the Jacobian of the transformation.
Fermionic Jacobian.
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} 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$. 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} 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}
Regularisation.
The expression above is not very useful as written, because $\text{tr}(\gamma_5)$ vanishes and $\delta(0)$ diverges, so we must introduce a regulator to be able to make sense out of this indeterminacy. For definiteness we introduce the gauge-invariant regulator \begin{equation}\tag{12} \delta(x-y)\gamma_5\to \gamma_5f(-\not D{}^2/\Lambda^2)\delta(x-y) \end{equation} where $\Lambda\to\infty$ is a cutoff, and $f$ is a function such that $f(0)\equiv 1$. We also assume that it is perfectly flat at the origin (so that the regulator does not affect the low-energy physics), and that it vanishes at infinity fast enough (so that it is able to kill the high-energy modes). For example, we could use a hard cutoff $f(s)=\Theta(1-s)$ or some smooth modification thereof (smoothness seems more sensible from the physical point of view). In practice, it suffices to assume that $f$ satisfies \begin{equation} \begin{aligned} s^jf^{(j)}(s)&\stackrel{s\to0}\longrightarrow\delta_{0j}\quad\, j=0,\dots,d/2-1\\ s^jf^{(j)}(s)&\stackrel{s\to\infty}\longrightarrow0 \qquad j=0,\dots,d/2-1 \end{aligned}\tag{13} \end{equation}
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}$. 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}
We note that when expanding $f$ in series in $\Lambda^{-2}$, we only need to keep the term $\not D^d$, because the rest either vanish on account of the gamma matrices properties, or vanishes in the limit $\Lambda\to\infty$. We therefore need to compute \begin{equation}\tag{15} \frac{1}{(d/2)!}\text{tr}\left(\gamma_5\not D^d\right)f^{(d/2)}(k^2) \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}
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}
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} as required.
Alternative method.
The calculation above is based on Feynman's functional integral formulation of quantum mechanics. For completeness, we also mention that it is possible to formulate the analysis directly in terms of operators. We merely sketch the calculation, and leave it to the reader to fill in the details. Our presentation is inspired by the point-splitting method of Ref.2, although we generalise their analysis to arbitrary $d$ (and to non-abelian algebras).
In the operator formalism, the axial current is formally defined as \begin{equation}\tag{21} j^\mu_A\sim \bar\psi\gamma^\mu \gamma_5t\psi \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$, 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$ 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} 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$ defined above is meaningless as written.
With this, we see that the naïve approach leads to an $0\times\infty$ indeterminacy, and we must proceed, once again, by introducing a regulator.
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$. 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} 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} 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 with $\delta_{\mu_1\mu_2}\cdots\delta_{\mu_{n-1}\mu_n}$.
With this, and taking the divergence of $j_A$ and 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, 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)$. This is, in principle, a straightforward computation.
Further reading.
Ref.1 does also analyse the calculation of the anomaly function for an arbitrary gauge theory.
References.
- Weinberg S. - Quantum theory of fields, Vol.2, chapter 22.
- Peskin, Schroeder - An introduction To Quantum Field Theory, chapter 19.