= Solution
In <Feynman gauge>, use the mode expansion
$$
\widehat A^\mu(x)=
\int\frac{d^3\mathbf k}{(2\pi)^3\,2|\mathbf k|}
\sum_{\lambda=0}^3
\left[
\epsilon^\mu_\lambda(\mathbf k)a_\lambda(\mathbf k)e^{-ikx}
+\epsilon^{\mu*}_\lambda(\mathbf k)a_\lambda^\dagger(\mathbf k)e^{ikx}
\right],
$$
with $k^0=|\mathbf k|$ and the covariant polarization completeness relation. For $x^0>y^0$, time ordering retains the annihilation-creation contraction and gives the positive-frequency term; for $x^0<y^0$, it gives the negative-frequency term. Combining them by a $k^0$ contour integral gives
$$
\boxed{
\langle0|\mathcal T\widehat A^\mu(x)\widehat A^\nu(y)|0\rangle
=\int_{C_F}\frac{d^4k}{(2\pi)^4}
\frac{-i\eta^{\mu\nu}}{k^2}
e^{-ik\cdot(x-y)}}.
$$
The <Feynman i-epsilon prescription> places the positive-energy pole $k^0=|\mathbf k|-i\epsilon$ below the real axis and the negative-energy pole $k^0=-|\mathbf k|+i\epsilon$ above it. Equivalently the denominator is $k^2+i\epsilon$ with the contour along the real axis.
Back to article page