QM anharmonic oscillator - Feynman diagrams, calculating free energy I have the quantum anharmonic oscillator:
$$H = \frac{\hat{p}^2}{2m} + \frac{1}{2}m\omega^2 \hat{x}^2 + \frac{\lambda}{4!}\hat{x}^4$$
and I want to find the free energy (so essentially $\log Z(\beta)$) up to 2nd order in $\lambda$ using Feynman diagrams.
I don't know how to do this. I know that the logarithm means one only needs to look at connected diagrams. The first order would be
$$A_1 = 3 \left(-\frac{\lambda}{\hbar}\right) \frac{1}{4!} \int_0^{\hbar \beta} d\tau \left< x^2(\tau)\right>_0 \left< x^2(\tau)\right>_0 $$
and since $\left< x^2(\tau)\right>_0  = \frac{\hbar}{m}G(0)$, where $G$ is the propagator, then
$$A_1 = 3 \left(-\frac{\lambda}{\hbar}\right) \frac{1}{4!} \hbar \beta \left(\frac{\hbar}{m}G(0)\right)^2 $$
For the second order there's two kinds of connected diagrams, with combinatorial factors $2{4\choose 2}{4\choose 2}=72$ and $4\cdot 3 \cdot 2 =24$. So there's two terms in this order:
$$A_2^1 = 72 \frac{1}{2!}\left(-\frac{\lambda}{\hbar}\right)^2 \frac{1}{4!^2} \int_0^{\hbar \beta} d\tau \int_0^{\hbar \beta} d\tau' \left[ \left< x(\tau)x(\tau') \right>_0 \right]^2\left< x^2(\tau)\right>_0 \left< x^2(\tau')\right>_0 $$
$$A_2^2 = 24 \frac{1}{2!}\left(-\frac{\lambda}{\hbar}\right)^2 \frac{1}{4!^2} \int_0^{\hbar \beta} d\tau \int_0^{\hbar \beta} d\tau' \left[ \left< x(\tau)x(\tau') \right>_0 \right]^4$$
which involve double integrals of powers of the propagator. For the harmonic oscillator, I found that it's given by a very complicated expression:
$$G(\tau - \tau') = \frac{\cosh(\hbar\beta\omega/2 - \omega |\tau - \tau'|)}{2\omega \sinh (\hbar \omega \beta /2)}$$
All in all, the free energy to second order would be $$F = -kT(1+A_1 + A_2^1 + A_2^2)$$
Now it seems that to get the $A_2^i$s I need to integrate $G^2(\tau - \tau')$ and $G^4(\tau - \tau')$ in order to get the result, and I'd really like to know that this is the right way to go before I even attempt it. The integrals look absolutely hellish so even if this is right, how do I calculate them? I don't have access to mathematica so I don't even know if they give a reasonable result.
 A: If the expression for the propagator you provided is correct, the integral is not daunting as it may at first seem. We essentially are faced with integrating,
$$\int_0^c d\tau \int_0^c d\tau' \, \cosh^n(a-b|\tau-\tau'|)$$
for $a,b, c>0$ and $n \in \mathbb Z$. I will demonstrate the $n=2$ case. We can expand the integrand as,
$$\cosh^2(a-b|\tau-\tau'|) = \frac12 + \frac14 e^{2a} e^{-2b|\tau-\tau'|} + \frac14e^{-2a}e^{2b|\tau-\tau'|}.$$
Thus, the only tricky part is figuring out how to integrate,
$$\int_0^c d\tau \int_0^c d\tau' \, e^{\gamma |\tau-\tau'|}, \quad \gamma \in \mathbb R.$$
We see we are integrating over a square $[0,c]\times[0,c]$. We can split this square diagonally into two regions, the top triangle with $\tau-\tau' <0$ and the bottom half with $\tau-\tau' >0$, and then we can sum the contributions. For the bottom half, we would have,
$$\int_0^c d\tau \int_0^\tau d\tau' \, e^{\gamma(\tau-\tau')} = \frac{1}{\gamma}\int_0^c d\tau \, (e^{\gamma \tau} -1) = \frac{1}{\gamma^2}(e^{\gamma c} - \gamma c - 1).$$
Now for the top region, $|\tau-\tau'| = \tau'-\tau$ and we have,
$$\int_0^c d\tau' \int_0^{\tau'} d\tau \, e^{\gamma(\tau'-\tau)} = \frac{1}{\gamma^2}(e^{\gamma c} - \gamma c - 1)$$
and so we have that,
$$\int_0^c d\tau \int_0^c d\tau' \, e^{\gamma |\tau-\tau'|} = \frac{2}{\gamma^2}(e^{\gamma c} - \gamma c - 1).$$
Going back to the original $n=2$ integral, using this result and simplifying yields,
$$\int_0^c d\tau \int_0^c d\tau' \, \cosh^2(a-b|\tau-\tau'|) = \frac{1}{4b^2} \left[ \cosh(2(a-bc)) - \cosh 2a + 2bc(bc+\sinh 2a)\right].$$

It is tempting to now generalise the result for all even powers. Using the generalised binomial theorem and the expansion of $\cosh$ in terms of exponentials, we have,
$$\cosh^{2n}(x) = \frac{1}{2^{2n}}\sum_{k=0}^{2n} \binom{2n}{k} e^{(2n-2k)x} = \frac{1}{2^{2n}} \binom{2n}{n} + \sum_{k\neq n} \binom{2n}{k} e^{(2n-2k)x}.$$
Hence the integral can be written as,
$$\int_0^c d\tau \int_0^c d\tau' \, \cosh^{2n}(a-b|\tau-\tau'|) = \frac{c^2}{2^{2n}}\binom{2n}{n}+\sum_{k\neq n} \binom{2n}{k} \int_0^c d\tau \int_0^c d\tau' \, e^{2(n-k)(a-b|\tau-\tau'|)}.$$
This integral we already know how to compute:
$$\int_0^c d\tau \int_0^c d\tau' \, e^{2(n-k)(a-b|\tau-\tau'|)} = e^{2a(n-k)}\int_0^c d\tau \int_0^c d\tau' \, e^{2b(k-n)|\tau-\tau'|}$$
$$=\frac{1}{2b^2(k-n)^2}e^{2a(n-k)} \left[ e^{2bc(k-n)} + 2bc(n-k)-1\right].$$
Thus, the integral reduces to the finite sum from $k=0,\dots, 2n$ (without $k=n$), as
$$\frac{c^2(2n)!}{2^{2n} (n!)^2} + \frac{1}{2b^2}\sum_{k \neq n} \frac{1}{(k-n)^2}\binom{2n}{k}e^{2a(n-k)}\left[ e^{2bc(k-n)} + 2bc(n-k)-1\right].$$
Returning the propagator, plugging in the constants, the integral of $G^{2n}(\tau-\tau')$ is the finite sum,
$$\boxed{\frac{(2n)!}{(2\omega)^{2n}\sinh^{2n} \hbar \omega \beta /2} \left[ \frac{\hbar^2 \beta^2}{2^{2n}(n!)^2 } + \frac{1}{2\omega^2} \sum_{k\neq n}\frac{1}{(k-n)^2(2n-k)! k!} e^{\hbar\omega \beta (n-k)} \left( e^{2\hbar\omega\beta(k-n)}+2\hbar\omega\beta(n-k)-1\right)\right]}$$
