I am still struggling to get my head around QFT and whilst I think I understand the method of generating functionals to compute correlation functions (as in my question here), my course notes often use Wick's theorem for calculations, and I'd really like to understand how that works.
So if I return to the standard example of scalar Yukawa theory with interaction Lagrangian:
$$\mathcal{L}_{\text{int}}=-g\psi^{\dagger}\psi\phi$$
I would like to know how to compute quantities like:
$$\langle \Omega|T[\psi^{\dagger}(x)\psi(y)\phi(z)]|\Omega\rangle$$
Where $T[...]$ represents the time-ordering operator. My understanding is that this can be expressed as:
$$\langle \Omega|T[\psi^{\dagger}(x)\psi(y)\phi(z)]|\Omega\rangle = \lim_{t\to\infty}\frac{\langle0|T[\psi_{I}^{\dagger}(x)\psi_{I}(y)\phi_{I}(z)\exp\left(-i\int_{-t}^{t}\mathcal{H}_{\text{int}}\:\mathrm{d}^{4}s\right)]|0\rangle}{\langle0|T[\exp\left(-i\int_{-t}^{t}\mathcal{H}_{\text{int}}\:\mathrm{d}^{4}s\right)]|0\rangle}$$
We can expand the exponential to (arbitrarily) first-order in $g$:
$$\langle0|T[\psi_{I}^{\dagger}(x)\psi_{I}(y)\phi_{I}(z)\exp\left(-i\int_{-t}^{t}\mathcal{H}_{\text{int}}\:\mathrm{d}^{4}s\right)]|0\rangle = \langle 0|T\left[\psi_{I}^{\dagger}(x)\psi_{I}(y)\phi_{I}(z)(1-ig\int_{-t}^{t}\psi_{I}^{\dagger}(s)\psi_{I}(s)\phi_{I}(s)\:\mathrm{d}^{4}s)\right]|0\rangle$$
By Wick's theorem, the zeroth order in $g$ has only 3 operators and so will inevitably yield a contribution of 0. The second gives:
$$\langle0|\int_{-t}^{t}T[\psi_{I}^{\dagger}(x)\psi_{I}(y)\phi_{I}(z)\psi_{I}^{\dagger}(s)\psi_{I}(s)\phi_{}(s)]\:\mathrm{d}^{4}s|0\rangle$$
However, I'm not sure what the best way to proceed is. I can use Wick's theorem to expand the time-ordered string of operators in terms of products of operators, but then where do I go?