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. We shall argue that
\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).
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 those cumbersome factors of $i$. Beware: it is still quite possible that I got some signs wrong here and there. If you find any sign that I missed, feel free to fix them.
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^\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[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.
The anomaly function.
Recall that the axial transformation is a global symmetry of the classical action. Therefore, if we let the gauge parameters $\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^a\partial\cdot j^5_a
\end{equation}
and this relation defines the axial current $j^5_a$. In the classical theory, this implies that the current is conserved, $\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^a\mathcal A_a}\mathrm d\psi\,\mathrm d\bar\psi
\end{equation}
where $\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^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}
\langle\partial\cdot j^5_a\rangle\equiv\mathcal A_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^{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=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^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_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.
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}$. 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_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}
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(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^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^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 OP 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.
Alternative method 1.
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 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{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^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^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\!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^5$ 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$ (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{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{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 contracted with $\delta_{\mu_1\mu_2}\cdots\delta_{\mu_{n-1}\mu_n}$.
With this, taking the divergence of $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 $(21)-(25)$. This is, in principle, a straightforward computation.
Alternative method 2.
So far we have discussed the anomaly in both the functional integral and the operator formalisms. One can also define a QFT directly by its Feynman rules (i.e., we regard the perturbative expansion as the definition of the QFT). It is reassuring to note that we may derive the axial anomaly directly from the computation of a certain Feynman diagrams (the $(d/2)!$ inequivalent polygons with $\frac{d}{2}+1$ sides), and the result turns out to be exact (no higher order corrections). This is known as the Adler-Bell-Jackiw theorem (cf. this PSE post). We shall not discuss this here any further.
Further reading.
In Ref.1 the author argues that the anomaly is associated to the eta invariant of $\not D$ (which is formally defined as the number of its positive eigenvalues minus the number of its negative eigenvalues; note that $i\not D$ is hermitian in the Minkowski setting, cf. this PSE post). To make this sum $\infty-\infty$ meaningful, the author proceeds by heat-kernel regularisation, which makes the analysis completely rigorous (at least, up to physicists standards).
Ref.1 does also analyse the calculation of the anomaly function for an arbitrary gauge theory, and discusses the relation between the eta invariant and the celebrated Atiyah-Singer index theorem.
References.
- Weinberg S. - Quantum theory of fields, Vol.2, chapter 22.
- Peskin, Schroeder - An introduction To Quantum Field Theory, chapter 19.