*"My problem here is the following: to find G we need boundary conditions of the problem. I can't understand, though, what boundary conditions we should impose here.*

*So to solve scattering problems using Green's functions like that, what are the boundary conditions we need to impose to compute the Green's function?"*

Why do you want to set a boundary condition? Clearly here your problem has no *boundary* since you can detect some particles very far away from the place where the interaction with the potential takes place. However, surely you will need an *intitial condition* provided by $\psi^0_\textbf{k}(\textbf{r})$.

I feel (tell me if I'm wrong) that what you mean by *"boundary condition"* refers actually to what kind of solution you are interested in. As you pointed out, the equation
$$
(\nabla^2+k^2)G(\textbf{r},\textbf{r}')=-4\pi\delta(\textbf{r}-\textbf{r}')
$$
has two solution called *retarded* $G^+$ and *advanced* $G^-$ Green functions. This gives rise to two different kind of solutions to your differential equation :

The causal solution :
$$
\psi_\textbf{k}(\textbf{r})=\psi^{\text{(in)}}_\textbf{k}(\textbf{r})-\frac{1}{4\pi}\int\mathrm{d}\textbf{r}'\,G^+(\textbf{r}-\textbf{r}')\,U(\textbf{r}')\,\psi_\textbf{k}(\textbf{r}')
$$
where $\psi^{\text{(in)}}_\textbf{k}\equiv \psi^{\text{0}}_\textbf{k}$ is interpreted as the *incoming field*, *i.e.* the asymptotic field you obtain when taking the limit $t\rightarrow -\infty$ and which evolves freely with time. This is usually the solution people are interested in because it describes the scattering event as a result of the interaction with the potential.

The anti-causal solution :
$$
\psi_\textbf{k}(\textbf{r})=\psi^{\text{(out)}}_\textbf{k}(\textbf{r})-\frac{1}{4\pi}\int\mathrm{d}\textbf{r}'\,G^-(\textbf{r}-\textbf{r}')\,U(\textbf{r}')\,\psi_\textbf{k}(\textbf{r}')
$$
where $\psi^{\text{(out)}}_\textbf{k}$ is interpreted as the *outcoming field*, *i.e.* the asymptotic field you obtain when taking the limit $t\rightarrow +\infty$ and which had evolved freely with time from the past.

In both of these cases, the limits $t\rightarrow\pm\infty$ allow you to get rid of the integral terms and give rise to different interpretations of the solutions.

- You may also be only interested in the radiated field which is defined as :
$$
\psi^{\text{(rad)}}_\textbf{k}(\textbf{r})=\psi^{\text{(out)}}_\textbf{k}(\textbf{r})-\psi^{\text{(in)}}_\textbf{k}(\textbf{r})=\frac{1}{4\pi}\int\mathrm{d}\textbf{r}'\,\mathcal{G}(\textbf{r}-\textbf{r}')\,U(\textbf{r}')\,\psi_\textbf{k}(\textbf{r}')
$$
with
$$
\mathcal{G}(\textbf{r})=G^+(\textbf{r})-G^-(\textbf{r})
$$

*"I've seem this method in some notes and as things are written it seems the only imposed condition is that $G(\textbf{r},\textbf{r}')=G(\textbf{r}-\textbf{r}')$."*

It seems to me that this is not directly related to your initial problem. One can show independently that the property
$$
G(\textbf{r},\textbf{r}')=G(\textbf{r}-\textbf{r}')
$$
is just a consequence of the fact that the dispersion relation $\langle\textbf{k}| H|\textbf{k}\rangle =E(\textbf{k})$ only depends on the *magnitude* of the impulsion, *i.e.* $E(\textbf{k})\equiv E(k)$, which is a consequence of the fact that your system is invariant under translation.

Any solution to the differential equation for the Green's function can
be used. If you solve for the Green's function using contour integration
over the frequencies, the solutions for the portion of the contour
near the pole are solutions of the homogeneous wave equation. Choosing
different contours corresponds to adding different solutions to the
homogeneous wave equation and therefore give different boundary conditions.
You could also find these solutions using other methods to
solve the differential equation without contour integration.

The reason to choose various boundary conditions is that when you use
Green's theorem to write the desired solution, you have terms where
you integrate over the source, and terms which come from the
space-time surface. If you can choose the boundary conditions for your
Green's function so that the surface terms are zero, your calculation
becomes simpler.
Every Green's function gives one particular solution.
Choosing the
retarded Green's function will make the time surface terms zero if the
fields are all zero before the source turns on. Often this is the solution
you want. Other boundary conditions solve other cases directly.

## Best Answer

Consider the case of a free scalar field, governed by the usual Lagrangian,

$$\mathcal{L} = \frac{1}{2}\partial_\mu \phi \partial^\mu \phi - \frac{1}{2}m^2 \phi^2$$

The propagator, or equivalently Green's function for the theory is a function which can be though of as a response when we use a delta function as an input in the equations of motion, i.e.

$$\square \Delta_F(x-y) \sim \delta^{(4)}(x-y)$$

Explicitly, it is given by a Fourier integral over four-momentum,

$$\Delta_F (x-y)=\int \frac{d^4 p}{(2\pi)^4} \frac{i}{p^2-m^2} e^{-ip \cdot (x-y)}$$

We encounter a singularity at $p^0 = \pm \sqrt{p^2+m^2} = \pm E_{p}$. We can choose a contour which avoids these by dipping below the first, and then above the other. However, to apply the residue theorem, it must be closed. For $x^0 > y^0$, we close it in an anti-clockwise direction in the upper-half plane, and the opposite if $y^0 > x^0$. Alternatively, we can define the Feynman propagator,

$$\Delta_F = \int \frac{d^4 p}{(2\pi)^4} \, \frac{i}{p^2-m^2 + i\epsilon} e^{-ip\cdot (x-y)}$$

The $i\epsilon$ prescription due to Feynman shifts the poles by an infinitesimal amount away from the real axis; as a result the integral going straight through the real line is equivalent to the aforementioned contour.

Depending on your purpose, it may be useful to pick a particular contour, in which case we can define the retarded propagator $\Delta_R$ as the one which chooses to go over each pole on the real line, and the advanced contour going under. See the depiction below:

To understand what they mean physically, consider in the context of response theory, the response function $\chi$ which determines how a system changes under the addition of a source $\phi$, i.e.

$$\delta \langle \mathcal{O}_i(t) \rangle = \int dt' \chi(t-t')\phi(t')$$

It's clear the aforementioned is in fact a convolution, $\chi \ast \phi$, and $\chi$ also has the interpretation of a Green's function. But we can't affect the past, so clearly,

$$\chi(t) = 0, \quad t < 0$$

For the Fourier transform, this is equivalent to stating,

$$\chi(\omega) \quad \mathrm{analytic}, \quad \mathrm{Im} \, \omega > 0$$

In other words, $\chi(\omega)$ has no poles in the upper-half plane. This object $\chi(t)$ is in fact our retarded Green's function, and is also called the causal Green's function, precisely because the above requirement is imposed.