Wick's theorem for non-equilibrium steady state I am working on a grand canonical Hamiltonian which has the form:
$$
\hat{K}=\hat{H}_{SC}+\hat{H}_{tip}+\hat{H}_{T}-\mu\hat{N}_{SC}-(\mu+eV)\hat{N}_{tip}
$$
where $\hat{H}_{T}=-t_0\sum_{\sigma}(c^{\dagger}_{r\sigma}\hat{d}_{\sigma}+h.c.)$ is the tunnelling term between the superconducting system and the tip of a scanning microscope that is coupled to a site $r$. For the following, I let
\begin{align*}
\hat{H_0}&=\hat{H}_{SC}+\hat{H}_{tip}-\mu\hat{N}_{SC}-\mu\hat{N}_{tip}\\
\hat{H}&=\hat{H}_{SC}+\hat{H}_{tip}+\hat{H}_{T}-\mu\hat{N}_{SC}-\mu\hat{N}_{tip}\\
\hat{K}_0&=\hat{H}_{SC}+\hat{H}_{tip}-\mu\hat{N}_{SC}-(\mu+eV)\hat{N}_{tip}\\
\end{align*}
My question is when I want to find the two particles contour-ordered Green's function, e.g.
$$
G_{cdcd}(\sigma \tau,\sigma'\tau')=i^2\left\langle\mathcal{T}_{C}\left[ c^{\dagger}_{r\sigma}(\tau)d_{\sigma}(\tau)c^{\dagger}_{r\sigma'}(\tau')d_{\sigma}(\tau')\right]\right\rangle
$$
where the operators evolve with the Hamiltonian $\hat{H}$. Can I just use wick's theorem to split it into one particle Green's function given that $\hat{K}$ is quadratic?
Since the expectation value is taken with the state
$$
\mid\rangle=\lim_{\eta\to 0^{+}}\mathcal{T}\exp\left[-i\int_{-\infty}^{0}\mathrm{d}t' e^{i\hat{H}_0 t'}e^{\eta t'}\hat{H}_{T}e^{-i\hat{H}_0 t'}\right]\mid\rangle_{0,V\neq 0}
$$
where $\mid\rangle_{0,V\neq 0}$ is the state that was in individual equilibrium at remote past, i.e. the ground state of $\hat{K}_0$.
If there is no bias, it is then obvious that $\hat{K}_0=\hat{H}_0$ and assuming no level crossing, adiabatic theorem states that $\mid\rangle$ is the ground state of $\hat{H}$. In this case, Wick's theorem can be used.
But if there is a voltage bias, $\hat{K}_0\neq \hat{H}_0$, though $\left[\hat{H}_0,\hat{K}_{0}\right]=0$ suggests that $\mid\rangle_{0,V\neq 0}$ is an eigenstate of $\hat{H}_0$, most likely not the ground state. By adiabatic theorem, $\mid\rangle$ will not be a ground state of $\hat{H}$ anymore. But I am not sure if it will still be a ground state of some weird looking quadratic Hamiltonian, thus making the Wick's theorem applicable.
 A: As you have pointed out, since the Hamiltonian is quadratic, Wick's theorem can be applied but the formalism must be suitably altered to account for the fact that the system is inherently a non-equilibrium one. The most important change being that the final and the initial states, which are used to contract the fermionic operators while calculating the expectation values, are different. Probably a way to see it is that the number operator doesn't commute with the SC and the tip parts separately, which is a consequence of the tunneling term. We typically avoid the calculation of an exact ground state which in itself is a highly non-trivial issue, instead connecting all expectation values to the non-interacting equilibrium states. Anyway, the Wick's theorem still applies, but in the sense of the Keldysh diagrammatic technique where one takes care of the Keldysh contour dependent times entering each correlator. Consequently one gets four kinds of two-point Green's functions on the standard contour basis depending on the choice of the contour (two for each field).
However, at the level of linear response, I guess one can still use the standard equilibrium Green's functions and the associated field theory. Probably because at the first order in perturbation, the energy levels are altered while the states remain unchanged.
See http://www.physics.arizona.edu/~stafford/Courses/560A/nonequilibrium.pdf, in particular Fig. 3.1 and the associated text on page 48. It mentions that the zero bias result can be obtained by linear response theory while the finite bias results require the full non-equilibrium theory.
Further, one of the initial papers and a seminal one on this topic is https://iopscience.iop.org/article/10.1088/0022-3719/4/8/018. See section 2.2 and specifically the text below Eq. (20), which too confirms the applicability of linear response theory near zero-bias.
