Which inverse of $-(\partial^2 + m^2)$ should be used in the path integral? The partition functional for free scalar field is 
$$Z=\int D\varphi e^{i\int d^4x[-\frac{1}{2}\varphi (\partial^2+m^2)\varphi+J\varphi]}.\tag{1}$$
To evaluate this functional integral, we usually compare it with the discrete case
$$\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}...\int_{-\infty}^{\infty} dq_1 dq_2...dq_N e^{(i/2)q\cdot A\cdot q+iJ\cdot q}=\left(\frac{(2\pi i)^N}{\det[A]}\right)^{\frac{1}{2}} e^{-(i/2)J\cdot A^{-1}\cdot J},$$
where $A^{-1}$ is the inverse of matrix $A$.
Now by analogy, we usually take the result of Eq.(1) as 
$$Z=\mathcal{C}e^{-(1/2)\int\int d^4x d^4y J(x)D(x-y)J(y)},\tag{2}$$
where $D(x-y)$ is said to be the inverse of the operator $-(\partial^2+m^2)$ (I think it should be $-(\partial_x^2+m^2)\delta(x-y)$ in a more rigorous sense). Namely,
$$\int d^4y -(\partial_x^2+m^2)\delta(x-y)D(y-z)=-(\partial_x^2+m^2)D(x-z)=\delta(x-z),\tag{3}$$
an analogue of $A^i_j{A^{-1}}^j_k=\delta^i_k$ in the discrete case.
But now I have a question. Since there is no unique inverse of the operator $-(\partial_x^2+m^2)\delta(x-y)$ (recall that we need extra boundary conditions to obtain a specific solution to Eq.(3), see my another question here), then which Green's function $D(x-y)$ should be used in Eq.(2)?
And what is the physical meaning that we choose the Green's function 
$$D(x-y)=\int \frac{d^4k}{(2\pi)^4}\frac{e^{ik(x-y)}}{k^2-m^2},\tag{4}$$
where I neglect the $i\epsilon$? Apparently, we can in principle add any function $g(x-y)$ to Eq.(4) provided 
$$-(\partial_x^2+m^2)g(x-z)=0.$$
 A: At some point you would want to compute your scalar field propagator from your generating functional 
$$Z=\mathcal{C}e^{-(1/2)\int\int d^4x d^4y J(x)D(x-y)J(y)},\tag{1}$$ and would obviously obtain $$\langle0|T \phi(x_1) \phi(x_2)|0\rangle = D(x_1-x_2). \tag{2}$$
The result of the physical propagator in the path integral approach should be identical to any other method of computing it. One could, for example use the canonical approach to compute the scalar field propagator:
\begin{align} \langle0|T \phi(x_1) \phi(x_2)|0\rangle= & \langle 0| \int \frac{{\rm d}^3 p \,{\rm d}^3 q}{(2\pi)^6} \frac{1}{\sqrt{4 E_p E_q}}  a_p e^{-ipx_1} a_q^\dagger e^{iqx_2}|0\rangle \\
=&  \int \frac{{\rm d}^3 p\, {\rm d}^3 q}{(2\pi)^6} \frac{1}{\sqrt{4 E_p E_q}} \langle 0| a_p  a_q^\dagger |0\rangle e^{-ipx_1 +iqx_2}\\
=& \int \frac{{\rm d}^3  p}{(2\pi)^3} \frac{1}{2 E_p } e^{-ip(x_1 -x_2)}.\tag{3}
\end{align}
Your Green's function $$D(x_1-x_2)=\int \frac{d^4k}{(2\pi)^4}\frac{i e^{ik(x_1-x_2)}}{k^2-m^2+i\epsilon},\tag{4}$$ can be shown to be equivalent to the above expression (3) by considering contour integration around the poles with the infinitesimal $i\epsilon$ in the denominator (cf. 44-46). By taking another form of Green's function $D(x_1 - x_2)$ the results for the different calculation methods would not match.
