In the first place, I am struggling when trying to derive the path integral formulation of the Green function for non-interacting particles
$$G_{ij}(\tau)=-\frac{1}{Z}\int D(\bar{\psi},\psi) \psi_i(\tau)\bar{\psi}_j(0)e^{-iS(\bar{\psi},\psi)}$$
with the action $$S=\sum_i\int^{\beta}_{0}d\tau \bar{\psi}_i(\tau)(i\partial_{\tau}+\epsilon_i-\mu)\psi_i(\tau),$$ from the definition
$$G_{ij}(\tau)=-\langle \text{T}_\tau \psi_i(\tau)\psi_j^\dagger(0)\rangle$$
with $\text{T}_\tau$ being the time-ordering operator and $\psi_i(\tau)=e^{\hat{H}\tau}\psi_i(0)e^{-\hat{H}\tau}$. I am told to use the fact that $\langle \hat{A} \rangle =-\frac{1}{Z}\text{tr}(e^{-\beta\hat{H}}\hat{A})$ and I know the general construction allowing to express traces in terms of path integrals, but the latter requires the operators in the argument of the trace to be normal ordered. If I try to rearrenge them in order to achieve this, the calculations get really messy. Do you have any ideas on what I am missing?
Secondly, how can I rewrite the path integral formulation of $G_{ij}(\tau)$ as a Matsubara summation via Gaussian integrals? In fact, the integral in the exponent prevents me to have a prototypical Gaussian integral which I can rewrite right away. It seems I am missing something here as well.