We start by considering the Klein-Gordon operator $\square +m^2$ acting in a space of functions $\phi(x)$ decaying sufficiently fast for $|x^\mu| \to \infty$. The general solution of the inhomogeneous Klein-Gordon equation
$$
(\square+m^2)\phi(x) = f(x) \tag{1} \label{1}
$$
can be written in the form
$$
\phi(x) = \phi_{\rm s}(x) + \phi_{\rm h}(x), \tag{2} \label{2}
$$
where $\phi_{\rm s}(x)$ is an arbitrary special solution of eq. \eqref{1} and $\phi_{\rm h}(x)$ denotes the general solution of the homogeneous equation $(\square +m^2)\phi_{\rm h}(x)=0$. Clearly, the inverse of $\square +m^2$ does not exist in this space of functions.
As a consequence, the Green function $G(x)$ of the Klein-Gordon operator, defined by
$$
(\square + m^2) G(x) = \delta^{(4)} (x) \tag{3} \label{3}
$$
is not uniquely determined by eq. \eqref{3} as all functions of the form $G^\prime(x)=G(x)+\phi_h(x)$ satisfy eq. \eqref{3} as well.
Writing $G(x)$ as a Fourier integral,
$$
G(x)= \int \! \frac{d^4k}{(2\pi)^4} \, e^{-ikx} \, \tilde{G}(k), \tag{4} \label{4}$$
eq. \eqref{3} becomes
$$
(-k^2+m^2)\, \tilde{G}(k) = 1 \tag{5} \label{5}
$$
and the ambiguity of the Green function is reflected by the fact that the integral
$$
\int \! \frac{d^4k}{(2\pi)^4} \, \frac{1}{m^2-k^2}= \int \! \frac{d^4k}{(2\pi)^4} \, \frac{1}{m^2+ \mathbf{k}^2- (k^0)^2} \tag{6} \label{6}
$$
does not exist because of the presence of the poles at $k^0 = \pm (\mathbf{k}^2+m^2)^{1/2}=: \pm \omega(\mathbf{k})$ on the $k^0$-axis.
However, by deforming the path of the $k^0$-integration into the complex $k^0$-plane, the poles on the real $k^0$ can be avoided, making sense of the integral given in \eqref{6}. Depending on the choice of the deformed path, one obtains e.g. the retarded Green function (with $G_{\rm R}(x)= 0$ for $x^0 < 0$), the advanced Green function (with $G_{\rm A}(x)=0$ for $x^0 >0$) or the Feynman propagator (in position space) $G_{\rm F}(x)$, characterized by containing only positive frequency plane waves $\sim \exp(- i \omega(\mathbf{k}) x^0)$ for $x^0 \to \infty$ and only negative frequency plane waves $\sim \exp(i \omega(\mathbf{k}) x^0)$ for $x^0 \to -\infty$. This procedure is equivalent to imposing certain boundary conditions on the Green function, singling out a specific solution of \eqref{3}.
Instead of deforming the path of integration in the complex $k^0$-plane, it is often slightly more convenient to shift the poles away from the real axis, while leaving the $k^0$-integration on the real axis.
Concentrating on the Green function $G_{\rm F}(x)$ with Feynman boundary conditions, shifting the poles to $k^0 = \pm \omega(\mathbf{k}) \mp i \epsilon$ with $\epsilon >0$ is equivalent to considering the differential equation
$$
(\square +m^2 - i \epsilon) \, G_{\rm F}(x) =\delta^{(4)}(x). \tag{7} \label{7}
$$
In contrast to eq. \eqref{3}, the solution of eq. \eqref{7} is now uniquely determined. The solutions of the corresponding homogeneous equation, being superpositions of plane waves $\sim \exp(\mp i \omega(\mathbf{k}) x^0) \exp(\mp \epsilon x^0)$, explode either for $x^0 \to \infty$ or $x^0 \to - \infty$ and are thus forbidden by our requirement $\phi(x) \to 0$ for $|x^\mu| \to \infty$.
In other words, the inverse of the operator $\square +m^2-i \epsilon$ exists in the space of functions under consideration and the Fourier representation of the corresponding Green function is given by
\begin{align}
G_{\rm F}(x) &= \int \! \frac{d^4k}{(2\pi)^4} \, \frac{e^{-ikx}}{m^2-k^2-i \epsilon} \tag{8} \label{8} \\[5pt]
&=
i \theta(x^0) \! \int \!\frac{d^3 k}{(2\pi)^3 2 \omega(\mathbf{k})} \, e^{-i \omega(\mathbf{k}) x^0} e^{i \mathbf{k} \cdot \mathbf{x}} \\[5pt] &+i \theta(-x^0) \! \int \! \frac{d^3 k}{(2\pi)^3 2 \omega(\mathbf{k})} e^{i \omega(\mathbf{k}) x^0 } e^{-i \mathbf{k} \cdot \mathbf{x}}.
\end{align}
For the discussion of the field equation of an $U(1)$ gauge field $A^\mu(x)$ we start from the Lagrangian
$$
\mathcal{L}= - \frac{1}{4} F_{\mu \nu}F^{\mu \nu} -j^\mu A_\mu \tag{9} \label{9}
$$
with the field strength tensor $F_{\mu \nu}= \partial_\mu A_\nu - \partial_\nu A_\mu$ and a conserved current $j^\mu$. The corresponding equation of motion
is found as $$
(g_{\mu \nu} \square - \partial_\mu \partial_\nu )A^\nu =j_\mu. \tag{10} \label{10}
$$
Again, the solution of eq. \eqref{10} is not unique. As in the previous case, we may always add a solution of the corresponding homogeneous equation (with $j^\mu=0$) to an arbitrary solution of the inhomogeneous equation \eqref{10}.
But we encounter an additional problem. The fields $A^\mu$ and $A^\mu+ \partial^\mu \Lambda$ (with an arbitrary scalar field $\Lambda(x)$) describe the same electromagnetic field given by the gauge invariant (physical) field $F_{\mu \nu}$. Indeed,
$$
(g_{\mu \nu} \square - \partial_\mu \partial_\nu) \partial^\nu \Lambda
= \square \partial_\mu \Lambda - \partial_\mu \square \Lambda =0 \tag{11}$$
and the $i \epsilon$ prescription alone is obviously not sufficient to cure also this second "disease".
As already well-known from classical electrodynamics, the way out is to impose a gauge condition on the gauge field $A_\mu$. In principle, this can be done in infinitely many ways. A convenient choice is provided by Fermi's trick, where a gauge fixing term is introduced in the Lagrangian,
$$
\mathcal{L} \to \mathcal{L}_{\rm eff} = \mathcal{L}- \frac{\xi}{2}(\partial_\mu A_\mu)^2 \tag{12} \label{12}
$$
with the resulting field equation
$$
\left[ g_{\mu \nu} \square -(1-\xi) \partial_\mu \partial_\nu \right] A^\nu = j_\mu. \tag{13} \label{13}
$$
Using again the $i \epsilon$ prescription, the operator
$$
g_{\mu \nu} (\square-i \epsilon) -(1-\xi) \partial_\mu \partial_\nu \tag{14}
$$
can be inverted, where the corresponding Green function is given by
$$
D_{\rm F}^{\mu \nu} (x)= \int \! \frac{d^4k}{(2\pi)^4} e^{-ikx} \frac{-1}{k^2+i \epsilon} \left[ g^{\mu \nu} -\left(1-\frac{1}{\xi}\right) \frac{k^\mu k^\nu}{k^2} \right]. \tag{15} \label{15}
$$
Note that observable quantities are, of course, independent of the the (arbitrary) gauge parameter $\xi$, being a good check in actual calculations.
