Briefly, I would like to have a reference to a clear detailed exposition of the computation of the multi-instanton contribution to the path integral while computing the energy levels splitting of the Hamiltonian with symmetric double well potential on the line.

In more detal, consider the symmetric double well potential on the line $$V(x)=\frac{\omega^2m}{8x_0^2}(x^2-x_0^2)^2,$$ where $\omega,m,x_0>0$.

Consider the following matrix element of the time evolution operator with imaginary time $$\langle x_0|e^{-HT/\hbar}|x_0\rangle,$$ where $H=-\frac{d^2}{dx^2}+V(x)$, $T\to +\infty,\hbar\to +0$. Using the path integral method one can compute this quantity as a sum of contributions of multi-instantons $$\langle x_0|e^{-HT/\hbar}|x_0\rangle=\sqrt{\frac{m\hbar \omega}{\pi}}e^{-\omega T/2}\cdot\sum_{N\mbox{ even}}\frac{(KT)^N}{N!},$$ where the constant $K$ is independent of $T$ and $\hbar$ and can be computed explicitly. The $N$th summand is interpreted as the constribution of $N/2$ pairs instanton - anti-instanton.

I would like to understand in all the details how to compute the contribution of configurations with $N>1$ (which is equal to $\sqrt{\frac{m\hbar \omega}{\pi}}e^{-\omega T/2}\cdot \frac{(KT)^N}{N!}$). Say the case $N=2$ would be interesting enough. A reference would be most helpful.

I have tried to read several expositions of the topic. I do understand computation of contribution with $N=0$ and with $N=1$ (the latter is for the similar problem of computation of $\langle x_0|e^{-HT/\hbar}|-x_0\rangle$). However the generalization to multi-instantons was not clear to me in all the literature I looked at.