Real-space correlations of massless Majorana fermions Consider an action of free, massless Majorana fermions in real time and 1+1 dimensions of the form
$$
S[\psi] = \frac{1}{2} \int d^2x \ \psi^{T}\gamma^0 (i \gamma^{\mu} \partial_{\mu}) \psi
$$
Here, $\psi$ is a two-component spinor. I would like to compute real space correlations of the field $\psi$. I've already done this problem in imaginary time, where the result is relatively simple and can be checked in many CFT books. I'm trying to check my understanding by doing it in real time as well, but I'm running into some issues.
For simplicity, I'll work in the chiral basis: set $\gamma^0 = \sigma^x$ and $\gamma^1 = -i\sigma^y$, so that the action reads (with $\psi = (\psi_+, \psi_-)^T$)
$$
S = \frac{1}{2} \int d^2 x \ \left\{ \psi_+ (i \partial_t + i \partial_x) \psi_+ + \psi_-(i \partial_t - i \partial_x) \psi_- \right\}
$$
From the action, I would guess that the correlator for $\psi_+$ is given by
$$
\langle \psi_+(t,x) \psi_+(0) \rangle = \int \frac{dp d\omega}{(2\pi)^2} \frac{i}{\omega-p+i\varepsilon} e^{ipx-i\omega t}
$$
However, this correlation function seems quite problematic to me, for a few reasons. For example, suppose I evaluate the $\omega$ integral first. Then I find that the propagator is zero for $t<0$, since the contour must be closed in the upper half plane. Similarly, the correlator is zero for $x<0$. How can these be? I would have thought that since $\psi_+$ is a Grassmann field inside the path integral,
$$
\langle \psi_+(t,x) \psi_+(t',x') \rangle = -\langle \psi_+(t',x') \psi_+(t,x) \rangle
$$
But this is clearly inconsistent with the above.
My suspicion is that this has something to do with an incorrect use of the $i\varepsilon$ prescription, since the $+i\varepsilon$ in the denominator seems to be the issue for all three problems. But I'm not sure how to handle these issues. In the bosonic case typically one might argue that the $i\varepsilon$ must be added to the action to aid in convergence, which is at least a nice mnemonic. But in the Grassmann case there should be no issue with unbounded Gaussian integrals.
So, here's my question: is my expression for the $\psi_+$ correlator correct? If so, what is wrong with my sanity checks following the expression? And if not, how should it be corrected?
 A: I believe I've found the answer to my own question. As discussed in the comments, the issue can be traced back to the $i\varepsilon$ prescription -- in my question, the $i\varepsilon$ prescription assumed essentially imposes a retarded response function. But as I stated, the path integral should give time-ordered correlation functions. So, how do I see what the correct $i\varepsilon$ prescription is?
To make the answer more transparent, it's worth noting that the same issue as in the Majorana action above occurs in a much more familiar setting as well. Consider the path integral for nonrelativistic free fermions:
$$
\hat{H} = \sum_{k} \xi_k c^{\dagger}_k c_k, \quad S[\bar{c},c] = \int dt \sum_k \bar{c}_k[i \hbar \partial_t - \xi_k]c_k
$$
Clearly, the same problem as before arises here: evaluating the correlation function $\langle Tc_k(t) c^{\dagger}_{k}(t') \rangle$ using the path integral also requires picking an $i\varepsilon$ prescription, and picking the naive $+i\varepsilon$ in the denominator gives the wrong answer.
The solution is to remember that correlation functions in the ground state are actually defined in real time by
$$
\langle c_k(t) \bar{c}_{k}(t') \rangle = \lim_{T \to \infty(1-i\varepsilon)} \frac{1}{Z} \int D[\bar{c},c] \ c_k(t) \bar{c}_k(t') \ \exp \left\{ \frac{i}{\hbar} \int_{-T}^T dt \sum_k \bar{c}_k[i \hbar \partial_t - \xi_k]c_k \right\}
$$
where $Z$ is the same expression without $c_k(t) \bar{c}_{k}(t')$. In other words, you have to take $t$ to lie along a slightly imaginary direction to project onto the ground state. With this in mind, the action should really be written
$$
S[\bar{c},c] = \int_{-\infty}^{\infty} dt \ (1-i\varepsilon) \sum_k \bar{c}_k [i \hbar (1+i\varepsilon) \partial_t - \xi_k] c_k
$$
With this in mind, the correlation function is really
$$
G(k,\omega) = \frac{1}{\hbar \omega -\xi_k(1-i\varepsilon)} = \frac{1}{\hbar \omega - \xi_k + i \varepsilon \text{sgn}(\xi_k)}
$$
which is the correct answer for the ground state momentum-frequency space Green's function.
So, back to my initial question: with the $T \rightarrow \infty (1-i\varepsilon)$ prescription in mind, the correct Green's function is
$$
\langle \psi_+(t,x) \psi_+(0,0) \rangle = \int \frac{dpd\omega}{(2\pi)^2} \frac{i}{\omega - p + i \varepsilon \text{sgn}(p)} e^{ipx-i\omega t}
$$
