Wick theorem on correlator in quantum mechanics I have an exercise to calculate the following one dimensional integral explicitly using the Wick theorem:
$$\langle q(t_1)q(t_2)q(t_3)q(t_4)q(t_5)q(t_6)\rangle=\frac{\int Dq \,q(t_1)q(t_2)q(t_3)q(t_4)q(t_5)q(t_6)e^{-S_E} }{\int Dq\,e^{-S_E}}$$
where
$$S_E=\int_{-\infty}^{+\infty}\Big(\frac{1}{2}\dot{q}^2+\frac{m^2q^2}{2}  \Big)dt. $$
I know what the answer is going to be as it is just using Wick theorem. But as the problem is to explicitly calculate, I assume that I actually need to calculate the integral myself.
This is what I have done:
$$\int Dq \,q(t_1)q(t_2)q(t_3)q(t_4)q(t_5)q(t_6)e^{-S_E}
=\left. \frac{\delta}{\delta f(t_1)}\cdots\frac{\delta}{\delta f(t_6)}\int e^{-S_E+\int\! dt~ fq}\right|_{f=0} . $$
I need some help in how to proceed from here.
 A: Sorry for the delay. 
I will work with integrals defined on $\mathbf{R}^{n}$.
A path integral can be thought as the limit $n \to  \infty$.
Generalizations to complex numbers are also straightforward.
Notation.
I will use the summation convention. Every repeated index is summed over.
For example,
$$
\sum_{i,j=1}^{n} x_{i} A_{ij} x_{j} \text{ is written as }
x_{i} A_{ij} x_{j}.
$$
Gaussian Integrals.
Consider the integral
$$
Z(A) = \int d^{n} x \, \exp\left( - \frac{1}{2} x_{i} A_{ij} x_{j} \right).
$$
It converges if the eigenvalues of A are non-negative and non-vanishing. One can also prove that
$$ 
Z(A) = ( 2 \pi )^{n/2} \left( \det A \right)^{-1/2}.
$$ 
by several methods. A possible way is to consider an orthogonal transformation which diagonalizes $A$, such that
$$
A = O D O^{T} , \quad O^{T} O = 1 \quad \text{and} \quad D_{ij} = a_{i} \delta_{ij}
$$
where $a_{i}$ are the eigenvalues of $A$. The integral factorizes because
$$
x_{i} A_{ij} x_{j} = x_{i} O_{ik} a_{k} O_{jk} x_{j} = a_{i} y_{i}^{2}
$$
for $y_{i} = O_{ij} x_{j}$. And it converges to
$$
Z(A) = \Pi_{i=1}^{n} \int dy_{i} e^{-a_{i} y_{i}^{2}/2} 
= (2 \pi)^{n/2} (a_{1} a_{2} \cdots a_{n})^{-1/2} = (2 \pi)^{n/2} (\det A)^{-1/2}.
$$
Of special interest is the case
$$
Z(A,b) = \int d^{n}x \exp \left( - \frac{1}{2} x_{i} A_{ij} x_{j} + b_{i} x_{i} \right)
$$
which gives
$$
Z(A, b) = (2 \pi)^{n/2} (\det A)^{-1/2} \exp \left[ \frac{1}{2} b_{i} \Delta_{ij} b_{j} \right]
$$
where $\Delta$ is the inverse of $A$ (prove this result yourself. Hint: Consider the change of variables 
$$
x_{i} = \Delta_{ij} b_{j} + y_i
$$
from $x_{i}$ to $y_{i}$).
Wick Theorem.
A Gaussian integrand can be considered a probability distribution. We can use it to calculate expectation values:
$$
\langle F(x) \rangle = N \int d^{n} x \, F(x) \, \exp \left( - \frac{1}{2}
x_{i} A_{ij} x_{j} \right)
$$
and the constant $N$ is determined from the condition $\langle 1 \rangle = 1$.
It is
$$
N = Z^{-1}(A, 0) = (2 \pi)^{-n/2} (\det A)^{1/2}. 
$$
The function
$$
Z(A,b)/Z(A,0) = \langle e^{b_{i} x_{i}} \rangle    
$$
is then the generating function of the moments of distribution. Expectation values are obtained by differentiating with respect to $b_{i}$:
$$
\langle x_{k_{1}} \cdots x_{k_{l}} \rangle =
\left( \frac{\partial}{\partial b_{k_{1}}} \cdots \frac{\partial}{\partial b_{k_{l}}}
\right) \exp\left[ \frac{1}{2} b_{i} \Delta_{ij} b_{j} \right] |_{b=0}.
$$
This can be inspected by expanding both sides of the generating function in powers of $b_{i}$. 
If $F(x)$ is polynomial in $x$, then
$$
\langle F(x) \rangle = F\left( \frac{\partial}{\partial b} \right) 
\exp \left[ \frac{1}{2} b_{i} \Delta_{ij} b_{j} \right]|_{b=0}.
$$
This is known as Wick's theorem. Each time one differentiates the exponential, one gets a factor of $b$. Since $b$ is set to zero at the end, one must differentiate this factor of $b$ later, otherwise the corresponding contribution vanishes. 
Thus, the expectation value, in a Gaussian theory, is given by all the possible ways of pairing the indices $k_{1}$ to $k_{l}$ in $\Delta_{k_{i} k_{j}}$. One finds, successively,
$$
\langle x_{i_1} x_{i_2} \rangle = \Delta_{i_{1} i_{2}} ,
$$
$$
\langle x_{i_{1}} x_{i_{2}} x_{i_{3}} x_{i_{4}} \rangle =
\Delta_{ i_{1} i_{2} } \Delta_{i_{3} i_{4}} + 
\Delta_{i_{1} i_{3} } \Delta_{i_{2} i_{4}} + 
\Delta_{i_{1} i_{4}} \Delta_{i_{2} i_{3}}.
$$
EDIT.
Quantum Mechanics.
You can now make the substitution $x_{i} \mapsto x(t)$. So instead of using vectors and partial derivatives, you will need functions and functional differentiation.
Of ultimate importance is $\Delta$, which is now called propagator. It was defined via
$$
A_{ij} \Delta_{jk}  = \delta_{ik}
$$ 
and now is interpreted in the continuum limit:
$$ 
A \Delta(t-t') = \delta (t - t'),
$$
where $A$ is the quadratic piece in the Lagrangian. In your case
$$
A = \frac{d^{2}}{dt^{2}} + \omega^{2}.
$$
In the literature, this $\Delta$ is known as the Green function.
I have given you all ingredients so that you can do this exercise.
Hope this gives you some insight.
A: I find useful the following alternative proof to that of OkThen. I learned it from https://arxiv.org/abs/1202.1554. In the same finite dimensional setting we have
$$0=\int_{\mathbb{R}^n}\text{d}^nx\, \frac{\partial}{\partial x^i}\left(e^{-\frac{1}{2}x^kA_{kl}x^l}x^{i_1}\cdots x^{i_n}\right)\\=-A_{ij}\int_{\mathbb{R}^n}\text{d}^nx\, e^{-\frac{1}{2}x^kA_{kl}x^l}x^jx^{i_1}\cdots x^{i_n}+\sum_{r=1}^n\int_{\mathbb{R}^n}\text{d}^nx\, e^{-\frac{1}{2}x^kA_{kl}x^l}x^{i_1}\cdots\hat{x^{i_r}}\cdots x^{i_n}\\=-A_{ij}\langle x^jx^{i_1}\cdots x^{i_n}\rangle+\sum_{r=1}^n\langle x^{i_1}\cdots\hat{x^{i_r}}\cdots x^{i_n}\rangle.$$
Denoting the inverse matrix by upper indices, we obtain
$$\langle x^ix^{i_1}\cdots x^{i_n}\rangle=\sum_{r=1}^nA^{ii_r}\langle x^{i_1}\cdots\hat{x^{i_r}}\cdots x^{i_n}\rangle.$$
Wick's theorem can then be obtained by induction from this formula: "Contract the first term of the correlation function with all the other terms, rinse and repeat."
