Boundary conditions / uniqueness of the propagators / Green's functions My question(s) concern the interpretation and uniqueness of the propagators / Green's functions for both classical and quantum fields.
It is well known that the Green's function for the Laplace equation
$$ \Delta_x G(x,x') = \delta^{(3)}(x-x') $$
with boundary conditions
$$ G(x,x') = 0 \text{ for } x\in\partial D, x'\in D$$
is determined uniquely. This is because any harmonic function that vanishes on the boundary of domain must vanish inside as well.
In contrast, I am confused about the uniqueness of the Green's function for the wave equation
$$ (∂_t^2 - \Delta)G(x,t;x',t') = \delta^{(3)}(x-x')\delta(t-t') $$
Without specifying boundary conditions, there exist many homogenuous solutions to the wave equation and thus many different Green's functions.
Physicists usually offer me several canonical choices, like the retarded, the advanced and also the Feynman propagator, but I don't understand what makes each choice unique compared to the others, or which boundary conditions correspond to which choice. Hence, my question is

Which boundary conditions correspond to retarded, the advanced and the Feynman propagator? What other possibilities for boundary conditions and propagators are there?

I am also confused by the situation in quantum field theory, where we have various conventions for propagators and time-ordering like
$$ \langle0|T_t\lbrace a(x,t)a^\dagger(x',t')\rbrace |0\rangle $$
Apparently, the ground state is very important for selecting the right Green's function (see also my previous question), but I still don't understand why that is.

How does the vacuum state act as a boundary condition in the spirit of the previous question?

There is also the issue of imaginary time. The point is that imaginary time turns the wave equation $(\partial_t^2 -\Delta)\phi=0$ into the Laplacian $(\partial_\tau^2 + \Delta)\phi = 0$, but I don't understand how the usual analytic continuation $\tau \to it \pm i\eta$ for various propagators depends on the boundary conditions. For instance, Tsvelik writes in equation (22.4) that the imaginary time Green's function for the Laplacian in 2D is "obviously" the one that vanishes when $\tau$ and $x$ go to infinity, but I don't understand the reasoning behind this choice.

What are the right boundary conditions for Green's functions in imaginary time? Is it clear that the path integral formalism selects the right ones?

 A: Retarded propagators are those with $G(\dots, t,t')=0$ for all $t<t'$. They're vanishing before $t=t'$, the delta-function "stimulates" the field at $t=t'$, and the Green's function for positive $t-t'$ measures the response of the field. One may view this description as a construction of the Green's function which also proves that it's unique.
Advanced propagators are the time reversal of the retarded ones, with the sign of $t-t'$ interchanged in the description above. (Recheck whether I haven't permuted the two propagators by a mistake: beware of sign errors.) Feynman propagators are the arithmetic averages of the corresponding retarded and advanced ones.
The (non-thermal) Green's functions are defined as correlation functions in the vacuum state. This is useful because the ground state is the simplest state and all other multiparticle states may be built from the ground state by adding linear combinations of fields (check the LSZ formula etc.). For this reason, the structure of the vacuum including all the correlations in it "knows" about all dynamical questions. 
In thermal field theory, the expectation values in the vacuum are replaced by ${\rm Tr}(\rho \dots)$ where $\rho=\exp(-\beta H)$ is the thermal density matrix.
The latter point also clarifies why the vacuum statet acts like a "boundary condition": the $T=0$ (non-thermal) correlators may be calculated for $t=\infty (1+i\epsilon)$ which has an infinite real part but a relatively smaller imaginary part, too. By analyticity, the small imaginary "angle" doesn't change the Green's function much. However, it simplifies the calculation because the evolution operator will contain an extra factor of $\exp(-\epsilon \infty H)$ which is still enough to exponentially suppress all excited states relatively to the ground state. When $\infty$ is really sent to infinity, only the matrix elements evaluated relatively to the ground state survive. The same comments and choices apply to the infinite future. For this reason, the vacuum matrix elements are automatically picked from the simplest boundary conditions one may assume at $t=\pm\infty$.
The small but nonzero value of $\epsilon$ gets translated to the $i\epsilon$ in Feynman's propagator: a procedure like this one also shows why Feynman's propagators are the correct ones to use in the Feynman diagrams.
