Green's function for free scalar field theory as an inverse of $-(\partial^2+m^2)$ In his book [QFT in a Nutshell], Zee argues that the functional integral
$$Z= \int D\varphi \; e^{i\int d^4x[ -\frac{1}{2}\varphi(\partial^2+m^2)\varphi + J\varphi]}$$
can be evaluated by discretizing the spacetime. To be specific, we imagine a grid of spacing $a$ and replace the function $\varphi(na)$ by the vector $\varphi_n = \varphi(na)$, the derivative $\partial\varphi(na)$ by $(1/a)(\varphi_{n+1}-\varphi_{n})$, and so on.
Then the discrete analog of the path integral is just a Gaussian integral
$$\int dq_1\cdots dq_N e^{(i/2)q\cdot Aq + iJ\cdot q} = \left( \frac{(2\pi i)^N}{\det[A]}\right)^{1/2} e^{-(i/2)J\cdot A^{-1}J}.$$
Motivated by this, Zee writes the final answer,
$$Z = C e^{-(i/2) \int d^4x d^4y J(x)D(x-y)J(y)},$$
where $D(x-y)$ is the Green function $-(\partial^2+m^2)D(x-y)=\delta^{(4)}(x-y).$
However, since the differential operator $(\partial^2+m^2)$ has a nonzero kernel, I'm curious about how the operator can be replaced by an invertible matrix $A$. This might be answered by the line of thoughts in this note. That is, by encoding appropriate boundary conditions into $A$. What should be the boundary condition here? Would $\varphi(x)\rightarrow 0$ as $\|x\|\rightarrow\infty$ be the answer?
Furthermore, how do I figure out the explicit form of $D(x-y)$, which is the continuum limit of $A^{-1}_{ij}?$ What is the reasoning behind Zee's assertion that
$$D(x-y)= \int \frac{d^4k}{(2\pi)^4} \frac{e^{ik(x-y)}}{k^2-m^2+i\epsilon}
= -i\int \frac{d^3k}{(2\pi)^3 2\omega_k}[e^{-i(\omega_k t-\vec{k}\cdot\vec{x})}\theta(x^0)+e^{i(\omega_k t-\vec{k}\cdot\vec{x})}\theta(-x^0)],$$
and not one of the other solutions, e.g. the retarded and advanced propagators $D_{ret}(x-y), D_{adv}(x-y)$? This question is closely related to this one, which was answered comparing the canonical and path integral formulations. Could it be done by imposing a boundary condition without referring to the canonical approach?
 A: In general you may choose whatever boundary conditions you like for your fields. The condition you mention is standard, but certainly not the only kind of condition that appears in application. It really depends on what you're looking to do.
To compute $D$, the standard procedure is to being Fourier transforming everything into momentum space. Once you do so, you'll note that the Fourier transform of $D$ may be solved for algebraically. Then simply Fourier transform back. This is precisely what Zee what done to obtain the expression you've noted for $D$. This will always work for Poincare invariant quadratic terms in the Lagrangian because explicit spacetime coordinate dependence is prohibited by translation symmetry, and hence the Fourier transform of the kinetic term will always be algebraic in momentum space.
You are correct that the operator $k^2-m^2$ is generically not invertible and indeed will not be invertible on-shell since then $k^2=m^2$. This is why the $i\epsilon$ is there. The operator $k^2-m^2+i\epsilon$ is invertible (for real-valued momenta) because it has no zeros on the real line as they have been shifted off it by the addition of the small complex number $i\epsilon$.
The proper way of introducing this $\epsilon$ into the action is to note that the partition function produces vacuum expectation values, and hence there is technically an insertion of $\Psi_0\Psi_0^*$ into the path integral, where $\Psi_0$ is the vacuum wave functional. Using canonical methods, you can show that the free theory vacuum wave functional is a functional Gaussian and hence this insertion would contribute a term like $\epsilon \int\phi^2$ to the action (there are some fine points about how this should be done. Not sure of a good source that shows this.). This is usually left out, but should be included in a careful analysis. So essentially it is the vacuum wave functional which is ensuring that our kinetic term is invertible.
The assertion about which propagator this represents is dependent only on the direction in which the two zeros in the kinetic term have been moved off the real line. The way these zeros are shifted is fixed by the vacuum wave functionals so there is no ambiguity if the calculation is carried out carefully.
As a final point, I will note that while there are many very interesting things that can be said about boundary conditions in QFT, essentially all of them lead to non-perturbative considerations and can be ignored almost entirely for the sake of a perturbative treatment of the theory (which is what I believe Zee, and essentially all introductory texts on the topic, do). By the way, the demand that the field vanish at infinity can be seen as equivalent to the demand that the action be finite. If you were to Wick rotate the action, you would very quickly see that configurations for which the rotated action (slightly different than the Lorentzian action) diverges have zero contribution to the path integral anyway (since the integrand is $e^{-S_E}$).
