Solution (source code)

= Solution

Since $D\Delta=i\delta$, the source-induced one-point function is
$$
\langle A_\mu(x)\rangle_j
=-ie\int d^4z\,\Delta_{\mu\rho}(x-z)j^\rho(z).
$$
The action is quadratic and the source is linear, so completing the square gives the exact full two-point function
$$
G_{\mu\nu}(x,y)=\Delta_{\mu\nu}(x-y)
-e^2\int d^4z\,d^4w\,
\Delta_{\mu\rho}(x-z)j^\rho(z)
\Delta_{\nu\sigma}(y-w)j^\sigma(w).
$$
Diagrammatically these are a free line joining $x$ to $y$ and a disconnected pair of lines, each joining one external insertion to one current cross. The connected two-point function remains exactly $\Delta_{\mu\nu}$. The expansion truncates at $e^2$ because a Gaussian integral has no interaction vertices and its mean is linear in $e$.