I'd like to know how to solve the dirac equation with some small gauge potential $\epsilon \gamma^\mu{A}_\mu(x)$ by applying perturbation theory. The equations reads as $$(\gamma^\mu\partial_\mu-m+\epsilon\gamma^\mu A_\mu(x))\psi(x) = 0.$$
Now I use the ansatz $\psi(x) = \psi_0(x)+\epsilon\psi_1(x)$ where $\psi_0$ is the free solution of the dirac equation. Putting the ansatz into the dirac equation we get for different orders:
$$\epsilon^0: (\gamma^\mu\partial_\mu-m)\psi_0(x) = 0,\\ \epsilon^1: (\gamma^\mu\partial_\mu-m)\psi_1(x) = -\gamma^\mu A_\mu(x) \psi_0(x),\\ \epsilon^2: \gamma^\mu A_\mu (x) \psi_1(x) = 0.$$
At this point I'm not sure how to solve the second equation such that $\psi_1 \in \ker(\gamma^\mu A_\mu)$.
The usual, conventional way we apply perturbation theory to the equation
$(\gamma^\mu \partial_\mu - m + \epsilon \gamma^\mu A_\mu(x)) \psi(x) = 0 \tag 1$
is to assume $\psi(x)$ can be expanded in powers of $\epsilon$:
$\psi(x) = \displaystyle \sum_0^\infty \psi_k(x) \epsilon^k; \tag 2$
we observe we allow an infinite number of term in (2), one for each non-negative power $\epsilon^k$ of $\epsilon$; when this series is inserted into (1) we find
$\displaystyle \sum_0^\infty \epsilon^k (\gamma^\mu \partial_\mu - m) \psi_k(x) + \sum_1^\infty e^k \gamma_\mu A_\mu(x) \psi_{k - 1}(x) = 0; \tag 3$
we equate terms with like powers of $\epsilon$:
$\epsilon^0: \; \; (\gamma^\mu \partial_\mu - m) \psi_0(x) = 0; \tag 4$
$\epsilon^1: \; \; (\gamma^\mu \partial_\mu - m) \psi_1(x) = -\gamma_\mu A_\mu(x) \psi_0(x); \tag 5$
$\epsilon^2: \; \; (\gamma^\mu \partial_\mu - m) \psi_2(x) = -\gamma_\mu A_\mu(x) \psi_1(x); \tag 6$
$\vdots \tag 7$
$\epsilon^k: \; \; (\gamma^\mu \partial_\mu - m) \psi_k(x) = -\gamma_\mu A_\mu(x) \psi_{k - 1}(x); \tag 8$
$\vdots \tag 9$
Now in principle we can find $\psi(x)$ to any order in $\epsilon$ by successively solving these equations (4)-(9); but of course in practice we usually in fact halt the process after some relatively small, finite number of iterations. From this point of view, that is, granting we are going to truncate the series (2) at some point, the key question becomes how best to do it?
The most commonly utilized method of which I am aware is to simply sequentially solve (4)-(9), taking as many iterations, as many coefficients $\psi_k(x)$ of the series (2), as is deemed necessary to obtain an approximate solution of sufficient accuracy. That is, we consecutively find the $\psi_k(x)$ out to some value of $k$, and take our solution to be
$\psi(x) = \displaystyle \sum_0^k \epsilon^j \psi_j(x); \tag{10}$
we also often take $\epsilon = 1$ after the iterative scheme (4)-(9) has been set up; then to order $k$ we have
$\psi(x) = \displaystyle \sum_0^k \psi_j(x); \tag{11}$
the choice $\epsilon = 1$ is often favored since it allows the perturbing term $\gamma^\mu A_\mu(x)$ to be $\epsilon$-independent, desirable since $\epsilon$ has no essential physical interpretation; taking $\epsilon = 1$ allows the operator $\gamma^\mu A_\mu(x)$ to retain its original physical significance.
It should be observed that the above scheme never requires setting
$\psi_k(x) = 0 \tag{12}$
a priori, as our OP Alpha001 has done in his question, where the assumption (12) with $k =2$ is made in concert with the foreshortened $\epsilon$-expansion
$\psi(x) = \psi_0(x) + \epsilon \psi_1(x); \tag{13}$
this as has been seen leads to the additional condition
$\gamma^\mu A_\mu(x) \psi_1(x) = 0, \tag{14}$
which must also be satisfied by $\psi_1(x)$; this presents difficulties of its own: for example, it seems quite likely that the operator $\gamma^\mu A_\mu(x)$ may be non-singular for many choices of $A_\mu(x)$, leading to
$\psi_1(x) = 0, \tag{15}$
in which case the entire perturbation scheme breaks down.
The upshot of these remarks is that it appears the the expansion (2) cannot be a priori truncated to a shortened series for some finite value of $k$ as in (10), (11); no, we must allow terms of all orders in $\epsilon$ at the beginning, and then truncate the solution after solving for as many $\psi_k(x)$ as required. Of course, such truncation is a simple operation in this case, since we merely need ignore terms beyond a certain order; but we do this "after the fact" (the fact of calculating solutions, as it were); nevertheless, it is an easy operation to perform since it is tantamount to simply terminating the iterative process at some value of $k$; conceptually, $\psi_{k + 1}(x)$ is available to us if we want it; we just don't want it.