text:
# High-Order Variational Perturbation Theory for the Free Energy

## I Introduction

Most physical problems can only be solved by approximation methods. One of them is perturbation theory which yields weak-coupling expansions. Unfortunately, they often do not converge.
The ground state energy of the anharmonic oscillator is the simplest example where this phenomenon can be studied. Algebraic recursion relations à la Bender and Wu Bender/Wu  yield perturbation series for the eigenvalues (energy) and eigenfunctions (wave functions) of the time-independent Schrödinger equation up to arbitrarily high orders. In Ref. exp  the calculation was extended to 250th order. The Bender-Wu recursion relation yields a power series for the anharmonic part of the wave function both in the coupling strength $`g`$ and in the coordinate $`x`$. The power series in $`x`$ can be cut off naturally by comparing the recursive results with those obtained from evaluating Feynman diagrams. The resulting weak-coupling series for the ground state energy does not converge for any value of the coupling strength. This paper deals with both problems: Obtaining high-order perturbation expressions and making them converge for all values of the coupling strength. It is organized as follows:
In Sec. II we perturbatively evaluate the path integral representation for the imaginary-time evolution amplitude of the anharmonic oscillator by means of a generalized Wick’s theorem weissbach ; bachmann . In Sec. III we represent the first-order results diagrammatically. Doing so, we demonstrate that the algebraic computational cost is very high for the diagrammatic approach. We also obtain a cross check for the results which are derived from a differential recursion relation in the subsequent Sec. IV. In order to cut down on the algebraic computational cost we introduce a strategy to exploit the symmetry property of the imaginary-time evolution amplitude in Sec. V, thus laying the foundation for our high-order results. In Sec. VI we combine the resulting algebraic recursion relation with the original differential recursion relation, thus generalizing the Bender-Wu approach Bender/Wu . From our perturbative results for the imaginary-time evolution amplitude we then gain a perturbation expression for the free energy of the anharmonic oscillator in Sec. VII which we check again diagrammatically in Sec. VIII. The perturbative results are then re-summed in Sec. IX by means of variational perturbation theory Kleinert  for intermediate coupling $`g=1`$ for which the usual weak-coupling series would diverge. This theory is a systematic extension of a simple variational approach, first developed by Feynman and Kleinert in the path integral formalism. Feynman introduced the path integral formalism as a quantization regulation, that represents the operator properties of quantum physics by fluctuations of the dynamical variables Feynman1 ; Feynmanhibbs . By extending analytically real time to imaginary time, also quantum statistical quantities can be obtained by summing over quantum mechanical and thermal fluctuations with the help of path integrals Feynmanhibbs ; Feynman2 . In order to evaluate the path integral for the free energy approximatively, Feynman and Kleinert developed a variational method in 1986 Feynman . It replaces the relevant system by the exactly solvable harmonic oscillator whose frequency becomes a variational parameter which has to be optimized. Starting with Ref. Kleinert2 , this method has been systematically extended by Kleinert to higher orders Kleinert ; Frohlinde . It is now known as variational perturbation theory and yields results for all temperatures and all coupling strengths.
In Sec. X we extend this procedure to higher orders of the free energy and cross check the results in Sec. XI.

## II Path Integral Representation

The path integral representation for the imaginary-time evolution amplitude of a particle of mass $`M`$ moving in a one dimensional potential $`V(x)`$ reads Kleinert

$`(x_b\mathrm{}\beta |x_a\mathrm{\hspace{0.17em}0})`$ $`=`$ $`{\displaystyle _{x(0)=x_a}^{x(\mathrm{}\beta )=x_b}}𝒟x\mathrm{exp}\{{\displaystyle \frac{1}{\mathrm{}}}{\displaystyle _0^\mathrm{}\beta }d\tau `$ (1)
$`[{\displaystyle \frac{M}{2}}\dot{x}^2(\tau )+V(x(\tau ))]\}.`$
For the anharmonic oscillator potential

$`V(x)={\displaystyle \frac{M}{2}}\omega ^2x^2+gx^4`$ (2)
the imaginary-time evolution amplitude (1) can be expanded in powers of the coupling constant $`g`$. Thus we obtain the perturbation series

$`(x_b\mathrm{}\beta |x_a\mathrm{\hspace{0.17em}0})=(x_b\mathrm{}\beta |x_a\mathrm{\hspace{0.17em}0})_\omega `$
$`\times \left[1{\displaystyle \frac{g}{\mathrm{}}}{\displaystyle _0^\mathrm{}\beta }𝑑\tau _1x^4(\tau _1)_\omega +\mathrm{}\right],`$
where we have introduced the harmonic imaginary-time evolution amplitude

$`(x_b\mathrm{}\beta |x_a\mathrm{\hspace{0.17em}0})_\omega {\displaystyle _{x(0)=x_a}^{x(\mathrm{}\beta )=x_b}}𝒟x`$
$`\times \mathrm{exp}\left\{{\displaystyle \frac{1}{\mathrm{}}}{\displaystyle _0^\mathrm{}\beta }𝑑\tau \left[{\displaystyle \frac{M}{2}}\dot{x}^2(\tau )+{\displaystyle \frac{M}{2}}\omega ^2x^2(\tau )\right]\right\},`$
and the harmonic path expectation value for an arbitrary functional $`F[x]`$:

$`F[x]_\omega {\displaystyle \frac{1}{(x_b\mathrm{}\beta |x_a\mathrm{\hspace{0.17em}0})_\omega }}{\displaystyle _{x(0)=x_a}^{x(\mathrm{}\beta )=x_b}}𝒟xF[x]`$
$`\times \mathrm{exp}\left\{{\displaystyle \frac{1}{\mathrm{}}}{\displaystyle _0^\mathrm{}\beta }𝑑\tau \left[{\displaystyle \frac{M}{2}}\dot{x}^2(\tau )+{\displaystyle \frac{M}{2}}\omega ^2x^2(\tau )\right]\right\}.`$
The latter is evaluated with the help of the generating functional for the harmonic oscillator, whose path integral representation reads

$`(x_b\mathrm{}\beta |x_a\mathrm{\hspace{0.17em}0})_\omega [j]={\displaystyle _{x(0)=x_a}^{x(\mathrm{}\beta )=x_b}}𝒟x\mathrm{exp}\{{\displaystyle \frac{1}{\mathrm{}}}{\displaystyle _0^\mathrm{}\beta }d\tau `$
$`\times [{\displaystyle \frac{M}{2}}\dot{x}^2(\tau )+{\displaystyle \frac{M}{2}}\omega ^2x^2(\tau )j(\tau )x(\tau )]\},`$ (6)
leading to Kleinert

$`(x_b\mathrm{}\beta |x_a\mathrm{\hspace{0.17em}0})_\omega [j]=(x_b\mathrm{}\beta |x_a\mathrm{\hspace{0.17em}0})_\omega `$ (7)
$`\times \mathrm{exp}[{\displaystyle \frac{1}{\mathrm{}}}{\displaystyle _0^\mathrm{}\beta }d\tau _1x_{\mathrm{cl}}(\tau _1)j(\tau _1)`$
$`+{\displaystyle \frac{1}{2\mathrm{}^2}}{\displaystyle _0^\mathrm{}\beta }d\tau _1{\displaystyle _0^\mathrm{}\beta }d\tau _2G^{(\mathrm{D})}(\tau _1,\tau _2)j(\tau _1)j(\tau _2)]`$
with the harmonic imaginary-time evolution amplitude

$`(x_b\mathrm{}\beta |x_a\mathrm{\hspace{0.17em}0})_\omega =\sqrt{{\displaystyle \frac{M\omega }{2\pi \mathrm{}\mathrm{sinh}\mathrm{}\beta \omega }}}`$
$`\mathrm{exp}\left\{{\displaystyle \frac{M\omega }{2\mathrm{}\mathrm{sinh}\mathrm{}\beta \omega }}[(x_a^2+x_b^2)\mathrm{cosh}\mathrm{}\beta \omega 2x_ax_b]\right\}.`$
In equation (7) we have introduced the classical path

$`x_{\mathrm{cl}}(\tau ){\displaystyle \frac{x_a\mathrm{sinh}(\mathrm{}\beta \tau )\omega +x_b\mathrm{sinh}\omega \tau }{\mathrm{sinh}\mathrm{}\beta \omega }},`$ (9)
and the Dirichlet Green’s function

$`G^{(\mathrm{D})}(\tau _1,\tau _2)`$ $``$ $`{\displaystyle \frac{\mathrm{}}{M\omega }}{\displaystyle \frac{1}{\mathrm{sinh}\mathrm{}\beta \omega }}`$
$`[\theta (\tau _1\tau _2)\mathrm{sinh}(\mathrm{}\beta \tau _1)\omega \mathrm{sinh}\omega \tau _2`$
$`+\theta (\tau _2\tau _1)\mathrm{sinh}(\mathrm{}\beta \tau _2)\omega \mathrm{sinh}\omega \tau _1].`$
We follow Ref. weissbach ; bachmann  and evaluate harmonic path expectation values of polynomials in $`x`$ arising from the generating functional (7) according to Wick’s theorem. Let us illustrate the procedure to reduce the power of polynomials by the example of the harmonic path expectation value $`x^n(\tau _1)x^m(\tau _2)_\omega `$:

1. Contracting $`x(\tau _1)`$ with $`x^{n1}(\tau _1)`$ and $`x^m(\tau _2)`$ leads to Green’s functions $`G^{(\mathrm{D})}(\tau _1,\tau _1)`$ and $`G^{(\mathrm{D})}(\tau _1,\tau _2)`$ with multiplicity $`n1`$ and $`m`$, respectively. The rest of the polynomial remains within the harmonic path expectation value, leading to $`x^{n2}(\tau _1)x^m(\tau _2)_\omega `$ and $`x^{n1}(\tau _1)x^{m1}(\tau _2)_\omega `$.
2. If $`n>1`$, extract one $`x(\tau _1)`$ from the path expectation value giving $`x_{\mathrm{cl}}(\tau _1)`$ multiplied by $`x^{n1}(\tau _1)x^m(\tau _2)_\omega `$.
3. Add the terms from (i) and (ii).
4. Repeat the previous steps until only products of path expectation values $`x(\tau _1)_\omega =x_{\mathrm{cl}}(\tau _1)`$ remain.

With the help of this procedure, we obtain to first order

$`x^4(\tau _1)_\omega =x_{\mathrm{cl}}^4(\tau _1)`$ (11)
$`+\mathrm{\hspace{0.17em}6}x_{\mathrm{cl}}^2(\tau _1)G^{(\mathrm{D})}(\tau _1,\tau _1)+3G^{(\mathrm{D})^2}(\tau _1,\tau _1).`$
## III Feynman Diagrams

These contractions can be illustrated by Feynman diagrams with the following rules: A vertex represents the integration over $`\tau `$

$`\text{ }\text{}={\displaystyle _0^\mathrm{}\beta }𝑑\tau ,`$ (12)
a line denotes the Dirichlet Green’s function

$`\text{ }\text{}=G^{(\mathrm{D})}(\tau _1,\tau _2),`$ (13)
and a cross or a “current” pictures a classical path

$`\text{ }\text{}=x_{\mathrm{cl}}(\tau _1).`$ (14)
Inserting the harmonic path expectation value (11) into the perturbation expansion (II) leads in first order to the diagrams

$`{\displaystyle _0^\mathrm{}\beta }𝑑\tau _1x^4(\tau _1)_\omega `$
$`=\text{ }\text{}+6\text{ }\text{}+3\text{ }\text{}.`$
We now evaluate the first-order Feynman diagrams in (III) for finite temperatures and arbitrary $`x_a,x_b`$. Thus we will get a first-order result for the imaginary-time evolution amplitude in (II). The first diagram leads to

$`\text{ }\text{}={\displaystyle \frac{1}{32\omega \mathrm{sinh}^4\mathrm{}\beta \omega }}[(x_a^4+x_b^4)(\mathrm{sinh}4\mathrm{}\beta \omega `$
$`8\mathrm{sinh}2\mathrm{}\beta \omega +12\mathrm{}\beta \omega )`$
$`+(x_a^3x_b+x_ax_b^3)(4\mathrm{sinh}3\mathrm{}\beta \omega +36\mathrm{sinh}\mathrm{}\beta \omega `$
$`48\mathrm{}\beta \omega \mathrm{cosh}\mathrm{}\beta \omega )`$
$`+x_a^2x_b^2(36\mathrm{sinh}2\mathrm{}\beta \omega +48\mathrm{}\beta \omega +24\mathrm{}\beta \omega \mathrm{cosh}2\mathrm{}\beta \omega )],`$
and the second diagram reduces to

$`\text{ }\text{}={\displaystyle \frac{\mathrm{}}{32M\omega ^2\mathrm{sinh}^3\mathrm{}\beta \omega }}[(x_a^2+x_b^2)(\mathrm{sinh}3\mathrm{}\beta \omega `$
$`+9\mathrm{sinh}\mathrm{}\beta \omega 12\mathrm{}\beta \omega \mathrm{cosh}\mathrm{}\beta \omega )`$
$`+x_ax_b(12\mathrm{sinh}2\mathrm{}\beta \omega +16\mathrm{}\beta \omega `$
$`+8\mathrm{}\beta \omega \mathrm{cosh}2\mathrm{}\beta \omega )],`$ (17)
whereas the last diagram turns out to be

$`\text{ }\text{}={\displaystyle \frac{\mathrm{}^2}{16M^2\omega ^3\mathrm{sinh}^2\mathrm{}\beta \omega }}(3\mathrm{sinh}2\mathrm{}\beta \omega `$
$`+4\mathrm{}\beta \omega +2\mathrm{}\beta \omega \mathrm{cosh}2\mathrm{}\beta \omega ).`$ (18)
So all in all we get the following first-order result for the imaginary-time evolution amplitude

$`(x_b\mathrm{}\beta |x_a\mathrm{\hspace{0.17em}0})=(x_b\mathrm{}\beta |x_a\mathrm{\hspace{0.17em}0})_\omega `$ (19)
$`\times (1{\displaystyle \frac{g}{\mathrm{}}}\{{\displaystyle \frac{\mathrm{}^2}{M^2\omega ^3\mathrm{sinh}^2\mathrm{}\beta \omega }}[{\displaystyle \frac{9}{16}}\mathrm{sinh}2\mathrm{}\beta \omega +{\displaystyle \frac{3}{4}}\mathrm{}\beta \omega +{\displaystyle \frac{3}{8}}\mathrm{}\beta \omega \mathrm{cosh}2\mathrm{}\beta \omega ]`$
$`+{\displaystyle \frac{\mathrm{}}{M\omega ^2\mathrm{sinh}^3\mathrm{}\beta \omega }}[(x_a^2+x_b^2)({\displaystyle \frac{3}{16}}\mathrm{sinh}3\mathrm{}\beta \omega +{\displaystyle \frac{27}{16}}\mathrm{sinh}\mathrm{}\beta \omega {\displaystyle \frac{9}{4}}\mathrm{}\beta \omega \mathrm{cosh}\mathrm{}\beta \omega )`$
$`+x_ax_b({\displaystyle \frac{9}{4}}\mathrm{sinh}2\mathrm{}\beta \omega +3\mathrm{}\beta \omega +{\displaystyle \frac{3}{2}}\mathrm{}\beta \omega \mathrm{cosh}2\mathrm{}\beta \omega )]`$
$`+{\displaystyle \frac{1}{\omega \mathrm{sinh}^4\mathrm{}\beta \omega }}[(x_a^4+x_b^4)({\displaystyle \frac{1}{32}}\mathrm{sinh}4\mathrm{}\beta \omega {\displaystyle \frac{1}{4}}\mathrm{sinh}2\mathrm{}\beta \omega +{\displaystyle \frac{3}{8}}\mathrm{}\beta \omega )`$
$`+(x_a^3x_b+x_ax_b^3)\left({\displaystyle \frac{1}{8}}\mathrm{sinh}3\mathrm{}\beta \omega +{\displaystyle \frac{9}{8}}\mathrm{sinh}\mathrm{}\beta \omega {\displaystyle \frac{3}{2}}\mathrm{}\beta \omega \mathrm{cosh}\mathrm{}\beta \omega \right)`$
$`+x_a^2x_b^2({\displaystyle \frac{9}{8}}\mathrm{sinh}2\mathrm{}\beta \omega +{\displaystyle \frac{3}{2}}\mathrm{}\beta \omega +{\displaystyle \frac{3}{4}}\mathrm{}\beta \omega \mathrm{cosh}2\mathrm{}\beta \omega )]\}+\mathrm{}).`$
The imaginary-time evolution amplitude thus has the time reversal behaviour

$`(x_b\mathrm{}\beta |x_a\mathrm{\hspace{0.17em}0})=(x_a\mathrm{}\beta |x_b\mathrm{\hspace{0.17em}0})^{},`$ (20)
while it is known that the imaginary-time evolution amplitude is real for one-dimensional problems.

## IV Partial Differential Equation

Consider the Schrödinger equation for the real-time evolution amplitude

$`i\mathrm{}{\displaystyle \frac{}{t}}(x_bt|x_a\mathrm{\hspace{0.17em}0})`$
$`={\displaystyle \frac{\mathrm{}^2}{2M}}{\displaystyle \frac{^2}{x_b^2}}(x_bt|x_a\mathrm{\hspace{0.17em}0})+V(x_b)(x_bt|x_a\mathrm{\hspace{0.17em}0}).`$
In order to get a corresponding quantum statistical Schrödinger equation we now have to change from real time to imaginary time, i.e. we have to perform the Wick rotation $`ti\tau `$. Thus the Schrödinger equation (IV) becomes

$`\mathrm{}{\displaystyle \frac{}{\tau }}(x_b\tau |x_a\mathrm{\hspace{0.17em}0})`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2M}}{\displaystyle \frac{^2}{x_b^2}}(x_b\tau |x_a\mathrm{\hspace{0.17em}0})`$ (22)
$`+V(x_b)(x_b\tau |x_a\mathrm{\hspace{0.17em}0}).`$
For both the real and the imaginary-time evolution amplitude the initial condition reads

$`(x_b\mathrm{\hspace{0.17em}0}|x_a\mathrm{\hspace{0.17em}0})=\delta (x_bx_a).`$ (23)
Plugging the anharmonic oscillator potential (2) into the Schrödinger equation (22) we finally get

$`\left\{\mathrm{}{\displaystyle \frac{}{\tau }}+{\displaystyle \frac{\mathrm{}^2}{2M}}{\displaystyle \frac{^2}{x_b^2}}{\displaystyle \frac{M}{2}}\omega ^2x_b^2gx_b^4\right\}`$
$`(x_b\tau |x_a\mathrm{\hspace{0.17em}0})=0.`$ (24)
Making the ansatz

$`(x_b\tau |x_a\mathrm{\hspace{0.17em}0})=(x_b\tau |x_a\mathrm{\hspace{0.17em}0})_\omega A(x_b,x_a,\tau ),`$ (25)
where $`(x_b\tau |x_a\mathrm{\hspace{0.17em}0})_\omega `$ is the harmonic imaginary-time evolution amplitude (II), we conclude from (IV) a partial differential equation for $`A(x_b,x_a,\tau )`$:

$`\left\{{\displaystyle \frac{}{\tau }}{\displaystyle \frac{\mathrm{}}{2M}}{\displaystyle \frac{^2}{x_b^2}}+\omega {\displaystyle \frac{x_b\mathrm{cosh}\omega \tau x_a}{\mathrm{sinh}\omega \tau }}{\displaystyle \frac{}{x_b}}+{\displaystyle \frac{g}{\mathrm{}}}x_b^4\right\}`$
$`A(x_b,x_a,\tau )=0.`$ (26)
We now choose our ansatz for $`A(x_b,x_a,\tau )`$ by introducing three expansions in $`g`$, in $`x_a`$ and in $`x_b`$, respectively. Also we take out the factor $`\mathrm{sinh}^l\omega \tau `$, such that the ordinary differential equations for the expansion coefficients become as simple as possible:

$`A(x_b,x_a,\tau )={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{k=0}{\overset{2n}{}}}{\displaystyle \underset{l=0}{\overset{2k}{}}}g^n{\displaystyle \frac{c_{2k|l}^{(n)}(\tau )}{\mathrm{sinh}^l\omega \tau }}x_a^{2kl}x_b^l.`$ (27)
In order to obtain the unperturbed result $`A(x_b,x_a,\tau )=1`$ for $`g=0`$ we need $`c_{0|0}^{(0)}(\tau )=1`$. The superscript $`n`$ in equation (27) denotes the perturbative order, whereas $`2k`$ counts the (even) powers of the various products $`x_a^ix_b^j`$. The summations over the coordinates $`x_a`$, $`x_b`$ can be truncated at $`k=4n`$, because we learn from Feynman diagrammatic considerations that the diagram with the most currents $`x`$ in the $`n`$th order looks like

$`\text{ }\text{}\text{ }\text{}\mathrm{}\text{ }\text{}.`$ (28)
Inserting the new ansatz (27) into the Schrödinger equation (IV) and arranging the indices in such a way that each term is proportional to $`x_a^{2kl}x_b^l`$, we get for the different powers of $`g`$ and for $`n>0`$:

$`{\displaystyle \underset{k=0}{\overset{2n}{}}}{\displaystyle \underset{l=0}{\overset{2k}{}}}{\displaystyle \frac{x_a^{2kl}x_b^l}{\mathrm{sinh}^l\omega \tau }}{\displaystyle \frac{c_{2k|l}^{(n)}(\tau )}{\tau }}{\displaystyle \frac{\mathrm{}}{2M}}{\displaystyle \underset{k=1}{\overset{2n1}{}}}{\displaystyle \underset{l=2}{\overset{2k2}{}}}(l+2)(l+1){\displaystyle \frac{c_{2k+2|l+2}^{(n)}(\tau )}{\mathrm{sinh}^{l+2}\omega \tau }}x_a^{2kl}x_b^l`$ (29)
$`\omega {\displaystyle \underset{k=0}{\overset{2n}{}}}{\displaystyle \underset{l=1}{\overset{2k1}{}}}(l+1){\displaystyle \frac{c_{2k|l+1}^{(n)}(\tau )}{\mathrm{sinh}^{l+2}\omega \tau }}x_a^{2kl}x_b^l+{\displaystyle \frac{1}{\mathrm{}}}{\displaystyle \underset{k=2}{\overset{2n}{}}}{\displaystyle \underset{l=4}{\overset{2k+4}{}}}{\displaystyle \frac{c_{2k4|l4}^{(n1)}(\tau )}{\mathrm{sinh}^{l4}\omega \tau }}x_a^{2kl}x_b^l=0.`$
Thus the sums over $`k`$ and over $`l`$ collapse and we determine the master equation for our coefficients $`c_{2k|l}^{(n)}(\tau )`$

$`{\displaystyle \frac{c_{2k|l}^{(n)}(\tau )}{\tau }}=(l+2)(l+1){\displaystyle \frac{\mathrm{}}{2M}}{\displaystyle \frac{c_{2k+2|l+2}^{(n)}(\tau )}{\mathrm{sinh}^2\omega \tau }}`$
$`+(l+1)\omega {\displaystyle \frac{c_{2k|l+1}^{(n)}(\tau )}{\mathrm{sinh}^2\omega \tau }}{\displaystyle \frac{1}{\mathrm{}}}c_{2k4|l4}^{(n1)}(\tau )\mathrm{sinh}^4\omega \tau ,`$ (30)
which is solved by

$`c_{2k|l}^{(n)}(\tau )=(l+2)(l+1){\displaystyle \frac{\mathrm{}}{2M}}{\displaystyle 𝑑\tau \frac{c_{2k+2|l+2}^{(n)}(\tau )}{\mathrm{sinh}^2\omega \tau }}`$ (31)
$`+(l+1)\omega {\displaystyle 𝑑\tau \frac{c_{2k|l+1}^{(n)}(\tau )}{\mathrm{sinh}^2\omega \tau }}`$
$`{\displaystyle \frac{1}{\mathrm{}}}{\displaystyle 𝑑\tau c_{2k4|l4}^{(n1)}(\tau )\mathrm{sinh}^4\omega \tau }+d_{2k|l}^{(n)}.`$
Here the $`d_{2k|l}^{(n)}`$ denote the integration constants which are fixed by applying the initial condition

$`\underset{\tau 0}{lim}\left|{\displaystyle \frac{c_{2k|l}^{(n)}(\tau )}{\mathrm{sinh}^l\omega \tau }}\right|<\mathrm{}.`$ (32)
However the above master equation (IV) is not valid for all $`k`$ and $`l`$. Therefore we now introduce a set of empirical rules telling us which of the coefficients $`c_{2k|l}^{(n)}(\tau )`$ have to be dropped once we write down (31) for any order $`n`$:

1. Drop all terms containing a $`c_{2k|l}^{(n)}(\tau )`$ where $`2k>4n`$.
2. Drop all terms containing a $`c_{2k|l}^{(n)}(\tau )`$ with $`l>2k`$.
3. Neglect all terms containing a $`c_{2k|l}^{(n)}(\tau )`$ with any negative indices $`k`$ and $`l`$.

To convince the reader that equation (31) together with this procedure leads to the correct results we now reobtain our first-order result from (19). To that end we set $`n=1`$, such that $`k`$ runs from 0 to 2 and $`l`$ from 0 to 4. Fixing $`k=2`$ and counting down from $`l=4`$ to $`l=0`$ we get

$`c_{4|4}^{(1)}(\tau )`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}}}{\displaystyle 𝑑\tau c_{0|0}^{(0)}(\tau )\mathrm{sinh}^4\omega \tau }+d_{4|4}^{(1)}={\displaystyle \frac{1}{\mathrm{}\omega }}\left({\displaystyle \frac{1}{32}}\mathrm{sinh}4\omega \tau {\displaystyle \frac{1}{4}}\mathrm{sinh}2\omega \tau +{\displaystyle \frac{3}{8}}\mathrm{sinh}\omega \tau \right),`$ (33)
$`c_{4|3}^{(1)}(\tau )`$ $`=`$ $`4\omega {\displaystyle 𝑑\tau \frac{c_{4|4}^{(1)}(\tau )}{\mathrm{sinh}^2\omega \tau }}+d_{4|3}^{(1)}={\displaystyle \frac{1}{\mathrm{}\omega \mathrm{sinh}\omega \tau }}\left({\displaystyle \frac{1}{8}}\mathrm{sinh}3\omega \tau +{\displaystyle \frac{9}{8}}\mathrm{sinh}\omega \tau {\displaystyle \frac{3}{2}}\mathrm{cosh}\omega \tau \right),`$ (34)
$`c_{4|2}^{(1)}(\tau )`$ $`=`$ $`3\omega {\displaystyle 𝑑\tau \frac{c_{4|3}^{(1)}(\tau )}{\mathrm{sinh}^2\omega \tau }}+d_{4|2}^{(1)}={\displaystyle \frac{1}{\mathrm{}\omega \mathrm{sinh}^2\omega \tau }}\left({\displaystyle \frac{9}{8}}\mathrm{sinh}2\omega \tau +{\displaystyle \frac{3}{2}}\omega \tau +{\displaystyle \frac{3}{4}}\omega \tau \mathrm{cosh}2\omega \tau \right),`$ (35)
$`c_{4|1}^{(1)}(\tau )`$ $`=`$ $`2\omega {\displaystyle 𝑑\tau \frac{c_{4|2}^{(1)}(\tau )}{\mathrm{sinh}^2\omega \tau }}+d_{4|1}^{(1)}={\displaystyle \frac{1}{\mathrm{}\omega \mathrm{sinh}^3\omega \tau }}\left({\displaystyle \frac{1}{8}}\mathrm{sinh}3\omega \tau +{\displaystyle \frac{9}{8}}\mathrm{sinh}\omega \tau {\displaystyle \frac{3}{2}}\omega \tau \mathrm{cosh}\omega \tau \right),`$ (36)
$`c_{4|0}^{(1)}(\tau )`$ $`=`$ $`\omega {\displaystyle 𝑑\tau \frac{c_{4|1}^{(1)}(\tau )}{\mathrm{sinh}^2\omega \tau }}+d_{4|0}^{(1)}={\displaystyle \frac{1}{\mathrm{}\omega \mathrm{sinh}^4\omega \tau }}\left({\displaystyle \frac{1}{32}}\mathrm{sinh}4\omega \tau {\displaystyle \frac{1}{4}}\mathrm{sinh}2\omega \tau +{\displaystyle \frac{3}{8}}\mathrm{sinh}\omega \tau \right).`$ (37)
Correspondingly, for $`k=1`$ we obtain

$`c_{2|2}^{(1)}(\tau )`$ $`=`$ $`{\displaystyle \frac{6\mathrm{}}{M}}{\displaystyle 𝑑\tau \frac{c_{4|4}^{(1)}(\tau )}{\mathrm{sinh}^2\omega \tau }}+d_{2|2}^{(1)}={\displaystyle \frac{1}{M\omega ^2\mathrm{sinh}\omega \tau }}\left({\displaystyle \frac{3}{16}}\mathrm{sinh}3\omega \tau +{\displaystyle \frac{27}{16}}\mathrm{sinh}\omega \tau {\displaystyle \frac{9}{4}}\omega \tau \mathrm{cosh}\omega \tau \right),`$ (38)
$`c_{2|1}^{(1)}(\tau )`$ $`=`$ $`{\displaystyle \frac{3\mathrm{}}{M}}{\displaystyle 𝑑\tau \frac{c_{4|3}^{(1)}(\tau )}{\mathrm{sinh}^2\omega \tau }}+2\omega {\displaystyle 𝑑\tau \frac{c_{2|2}^{(1)}(\tau )}{\mathrm{sinh}^2(\tau )}}+d_{2|1}^{(1)}`$ (39)
$`=`$ $`{\displaystyle \frac{1}{M\omega ^2\mathrm{sinh}^2\omega \tau }}\left({\displaystyle \frac{9}{4}}\mathrm{sinh}2\omega \tau +3\omega \tau +{\displaystyle \frac{3}{2}}\omega \tau \mathrm{cosh}2\omega \tau \right),`$
$`c_{2|0}^{(1)}(\tau )`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{M}}{\displaystyle 𝑑\tau \frac{c_{4|2}^{(1)}(\tau )}{\mathrm{sinh}^2\omega \tau }}+\omega {\displaystyle 𝑑\tau \frac{c_{2|1}^{(1)}(\tau )}{\mathrm{sinh}^2(\tau )}}+d_{2|0}^{(1)}`$ (40)
$`=`$ $`{\displaystyle \frac{1}{\mathrm{}\omega ^2\mathrm{sinh}^3\omega \tau }}\left({\displaystyle \frac{3}{16}}\mathrm{sinh}3\omega \tau +{\displaystyle \frac{27}{16}}\mathrm{sinh}\omega \tau {\displaystyle \frac{9}{4}}\omega \tau \mathrm{cosh}\omega \tau \right).`$
Finally for $`k=0`$ we get the equation

$`c_{0|0}^{(1)}(\tau )`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{M}}{\displaystyle 𝑑\tau \frac{c_{2|2}^{(1)}(\tau )}{\mathrm{sinh}^2\omega \tau }}+d_{0|0}^{(1)}={\displaystyle \frac{\mathrm{}}{M^2\omega ^3\mathrm{sinh}^2\omega \tau }}\left({\displaystyle \frac{9}{16}}\mathrm{sinh}2\omega \tau +{\displaystyle \frac{3}{4}}\omega \tau +{\displaystyle \frac{3}{8}}\omega \tau \mathrm{cosh}2\omega \tau \right).`$ (41)
The path of recursion which follows from this procedure is shown in Fig. 1.

## V Exploiting the symmetries

As seen above we already have to solve nine ordinary differential equations for the first-order imaginary-time evolution amplitude. For any order $`n`$ the number $`p`$ of integrals to solve is

$`p={\displaystyle \underset{j=1}{\overset{2n+1}{}}}(2j1)=4n^2+4n+1.`$ (42)
Due to the time reversal behaviour (20), the coefficients $`c_{2k|l}^{(n)}(\tau )`$ show a symmetry, namely:

$`{\displaystyle \frac{c_{2k|l}^{(n)}(\tau )}{\mathrm{sinh}^l\omega \tau }}={\displaystyle \frac{c_{2k|2kl}^{(n)}(\tau )}{\mathrm{sinh}^{2kl}\omega \tau }}.`$ (43)
Exploiting the symmetry (43), we can cut down the number (42) considerably. At first sight it is reduced to

$`p^{}={\displaystyle \underset{j=1}{\overset{2n+1}{}}}j=2n^2+3n+1,`$ (44)
so there are only six integrals left for the first order. But we can go even further. Employing these symmetries we can eventually change almost all recursive differential equations into purely algebraic ones leaving only $`p^{\prime \prime }=(2n+1)`$ integrations. So for the first order we are left with three integrations only, namely with equations (33), (38), and (41). These coefficients $`c_{4|4}^{(1)}(\tau )`$, $`c_{2|2}^{(1)}(\tau )`$, and $`c_{0|0}^{(1)}(\tau )`$ are integrated recursively. The other coefficients can then be obtained algebraically: Once we have $`c_{4|4}^{(1)}(\tau )`$ we also know $`c_{4|0}^{(1)}(\tau )`$ because of the symmetry (43). Comparing equation (IV) for $`k=2,l=4`$ and $`k=2,l=0`$ we then obtain an algebraic equation for $`c_{4|1}^{(1)}(\tau )`$. The knowledge of $`c_{4|1}^{(1)}(\tau )`$ gives us $`c_{4|3}^{(1)}(\tau )`$ because of the symmetry (43) and by comparing (IV) this time for $`k=2,l=3`$ on the one hand and $`k=2,l=1`$ on the other hand we are left with an algebraic equation for $`c_{4|2}^{(1)}(\tau )`$. Thus we get all the coefficients for $`k=2`$ only by solving one differential equation, namely the one for $`c_{4|4}^{(1)}(\tau )`$. For $`k=1`$ the procedure is similar, $`k=0`$ only generates one coefficient anyway, namely $`c_{0|0}^{(1)}(\tau )`$, which still has to be solved by evaluating one integral. The new path of recursion is shown in Fig. 2.
So finally three out of the nine first-order coefficients are obtained by integration, three more are clear for symmetry reasons and three come from an algebraic recursion.
We now generalize the algebraic part of our recursion. Consider again the symmetry property (43). Differentiation on both sides yields

$`{\displaystyle \frac{c_{2k|l}^{(n)}(\tau )}{\tau }}={\displaystyle \frac{1}{\mathrm{sinh}^{2k2l}\omega \tau }}{\displaystyle \frac{c_{2k|2kl}^{(n)}(\tau )}{\tau }}`$
$`2(kl)\omega \mathrm{cosh}\omega \tau {\displaystyle \frac{c_{2k|2kl}^{(n)}(\tau )}{\mathrm{sinh}^{2k2l+1}\omega \tau }}.`$ (45)
Now we substitute for the two partial derivatives according to equation (IV). Solving for the $`(l+1)`$-st coefficient and shifting the index $`l`$ down by one we obtain

$`c_{2k|l}^{(n)}(\tau )`$ $`=`$ $`{\displaystyle \frac{(l+1)\mathrm{}}{2M\omega }}c_{2k+2|l+1}^{(n)}(\tau )+{\displaystyle \frac{c_{2k4|l5}^{(n1)}(\tau )}{\mathrm{}\omega l}}\mathrm{sinh}^6\omega \tau `$ (46)
$`+{\displaystyle \frac{(2kl+3)(2kl+2)\mathrm{}}{2M\omega l}}{\displaystyle \frac{c_{2k+2|2kl+3}^{(n)}(\tau )}{\mathrm{sinh}^{2k2l+2}\omega \tau }}+{\displaystyle \frac{2kl+2}{l}}{\displaystyle \frac{c_{2k|2kl+2}^{(n)}(\tau )}{\mathrm{sinh}^{2k2l+2}\omega \tau }}`$
$`{\displaystyle \frac{1}{\mathrm{}\omega l}}{\displaystyle \frac{c_{2k4|2kl3}^{(n1)}(\tau )}{\mathrm{sinh}^{2k2l4}\omega \tau }}{\displaystyle \frac{(2k2l+2)\mathrm{cosh}\omega \tau }{l}}{\displaystyle \frac{c_{2k|2kl+1}^{(n)}(\tau )}{\mathrm{sinh}^{2k2l+1}\omega \tau }},`$
which is the algebraic recursion relation for any non-diagonal coefficient $`c_{2k|l}^{(n)}(\tau )`$ with $`0<lk`$. (The coefficients with $`k<l<2k`$ are then clear for symmetry reasons.) The diagonal coefficients $`c_{2k|2k}^{(n)}(\tau )`$ still have to be integrated.

## VI Combined differential and algebraic equation

We now combine the differential recursion with the algebraic one. As only the diagonal coefficients have to be evaluated by integrating the differential recursive equation, we can even further simplify the solution (31) to our master equation (IV). We only need it for the diagonal coefficients, for which $`l+1=2k+1`$ is always greater than $`2k`$. And according to our index rule (ii), coefficients of the shape $`c_{2k|2k+1}^{(n)}`$ have to be neglected. We get

$`c_{2k|2k}^{(n)}(\tau )=(2k+2)(2k+1){\displaystyle \frac{\mathrm{}}{2M}}{\displaystyle 𝑑\tau \frac{c_{2k+2|2k+2}^{(n)}(\tau )}{\mathrm{sinh}^2\omega \tau }}`$
$`{\displaystyle \frac{1}{\mathrm{}}}{\displaystyle 𝑑\tau c_{2k4|2k4}^{(n1)}(\tau )\mathrm{sinh}^4\omega \tau }+d_{2k|2k}^{(n)}.`$ (47)
Index rules (i) and (iii) still have to be applied, $`k`$ runs from 0 to $`2n`$.
Let us quickly summarize the combined differential and algebraic recursion relation considering the first order as an example. Fig. 2 shows all first-order coefficients for the imaginary time evolution amplitude. Each coefficient is represented by a little circle. Now the coefficients on the diagonal line $`2k=l`$ have to be obtained by referring to equation (VI) together with rules (i) and (iii). These two rules tell us which of the coefficients either from the the same order $`n`$ or from the previous order $`n1`$ have to be integrated and which ones can be put to zero.
Once we have the diagonal coefficients $`c_{2k|2k}^{(n)}(\tau )`$ we can calculate the off-diagonal ones with $`lk`$ with the help of equation (46). The coefficients with $`k<l<2k`$ are then clear for symmetry reasons.
Using the computer algebra programme Maple V R7 we managed to calculate seven perturbative orders of the imaginary-time evolution amplitude which can be found at internet .

## VII Free Energy

In this section we obtain perturbative results for the partition function by integrating the diagonal elements of our perturbative expression for the imaginary-time evolution amplitude from the previous sections:

$`Z={\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑x(x\mathrm{}\beta |x\mathrm{\hspace{0.17em}0}).`$ (48)
From the partition function we then compute the free energy perturbatively:

$`F={\displaystyle \frac{1}{\beta }}\mathrm{log}Z.`$ (49)
We have to expand the logarithm in order to obtain a perturbation expansion for the free energy $`F`$. For the first order we insert (19) together with (II) into (48) and evaluate the integral and for the second order we use, correspondingly, the data from Ref. internet . By taking the logarithm we get with (49) and with the expansion for the logarithm for the free energy to second order

$`F^{(2)}(\beta )={\displaystyle \frac{1}{\beta }}\mathrm{log}2\mathrm{sinh}{\displaystyle \frac{\mathrm{}\beta \omega }{2}}`$ (50)
$`+{\displaystyle \frac{3g\mathrm{}^2}{4M^2\omega ^2}}\mathrm{coth}^2{\displaystyle \frac{\mathrm{}\beta \omega }{2}}{\displaystyle \frac{g^2\mathrm{}^3}{64M^4\omega ^5}}({\displaystyle \frac{54\mathrm{}\beta \omega }{\mathrm{sinh}^4\frac{\mathrm{}\beta \omega }{2}}}`$
$`+{\displaystyle \frac{36\mathrm{}\beta \omega \mathrm{cosh}\mathrm{}\beta \omega +60\mathrm{sinh}\mathrm{}\beta \omega +21\mathrm{sinh}2\mathrm{}\beta \omega }{\mathrm{sinh}^4\frac{\mathrm{}\beta \omega }{2}}}).`$
The higher orders are omitted for the sake of keeping the type face clear.
With Maple we came as high as the fifth variational order which is two orders more than what has been obtained in previous work Meyer .

## VIII Diagrammatical Check

It is possible to check the perturbative results for the free energy for all temperatures. Namely, we can expand $`Z`$ in terms of harmonic expectations in a similar way as for the imaginary-time evolution amplitude in (II). To that end we need the generating functional

$`Z[j(\tau )]`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑x(x\mathrm{}\beta |x\mathrm{\hspace{0.17em}0})_\omega [j]`$ (51)
which we get from (7)-(II). It is of the form

$`Z[j(\tau )]=Z[0]`$ (52)
$`\times \mathrm{exp}\left[{\displaystyle \frac{1}{2\mathrm{}^2}}{\displaystyle _0^\mathrm{}\beta }𝑑\tau _1{\displaystyle _0^\mathrm{}\beta }𝑑\tau _2G^{(\mathrm{p})}(\tau _1,\tau _2)j(\tau _1)j(\tau _2)\right],`$
where the harmonic partition function reads

$`Z[0]={\displaystyle \frac{1}{2\mathrm{sinh}\frac{\mathrm{}\beta \omega }{2}}}`$ (53)
and

$`G^{(\mathrm{p})}(\tau _1,\tau _2)={\displaystyle \frac{\mathrm{}}{2M\omega }}{\displaystyle \frac{\mathrm{cosh}\left(\frac{\mathrm{}\beta \omega }{2}|\tau _1\tau _2|\omega \right)}{\mathrm{sinh}\frac{\mathrm{}\beta \omega }{2}}}`$ (54)
denotes the periodic Green’s function of the harmonic oscillator. We now obtain the partition function $`Z`$ of the anharmonic oscillator from the generating functional $`Z[j(\tau )]`$ by differentiating with respect to the current $`j(\tau )`$ while setting $`j(\tau )=0`$ afterwards:

$`Z=\mathrm{exp}\left\{{\displaystyle \frac{1}{\mathrm{}}}{\displaystyle _0^\mathrm{}\beta }𝑑\tau g\left[{\displaystyle \frac{\mathrm{}\delta }{\delta j(\tau )}}\right]^4\right\}Z[j(\tau )]|_{j=0}.`$ (55)
Thus we get

$`Z=Z[0]\{1{\displaystyle \frac{3g}{\mathrm{}}}{\displaystyle _0^\mathrm{}\beta }d\tau _1G^{(\mathrm{p})^2}(\tau _1,\tau _1)`$ (56)
$`+{\displaystyle \frac{g^2}{2\mathrm{}^2}}{\displaystyle _0^\mathrm{}\beta }d\tau _1{\displaystyle _0^\mathrm{}\beta }d\tau _2[9G^{(\mathrm{p})^2}(\tau _1,\tau _1)G^{(\mathrm{p})^2}(\tau _2,\tau _2)`$
$`+72G^{(\mathrm{p})}(\tau _1,\tau _1)G^{(\mathrm{p})^2}(\tau _1,\tau _2)G^{(\mathrm{p})}(\tau _2,\tau _2)`$
$`+24G^{(\mathrm{p})^4}(\tau _1,\tau _2)]+\mathrm{}\}.`$
In terms of Feynman diagrams this reads

$`Z=Z[0][1{\displaystyle \frac{3g}{\mathrm{}}}\text{ }\text{}`$ (57)
$`+{\displaystyle \frac{g^2}{2\mathrm{}^2}}(9\text{ }\text{}\text{ }\text{}`$
$`+\mathrm{\hspace{0.17em}72}\text{ }\text{}+24\text{ }\text{})+\mathrm{}]`$
$`=\mathrm{exp}[{\displaystyle \frac{1}{2}}\text{ }\text{}{\displaystyle \frac{3}{\mathrm{}}}\text{ }\text{}`$
$`+{\displaystyle \frac{g^2}{2\mathrm{}^2}}(72\text{ }\text{}+24\text{ }\text{})+\mathrm{}],`$ (58)
where we have introduced the symbol

$`{\displaystyle \frac{1}{2}}\text{ }\text{}\mathrm{log}Z[0].`$ (59)
Once we rewrite the partition function $`Z`$ in the form of the cumulant expansion as on the right hand side of equation (57), the disconnected Feynman diagrams disappear Kleinert . Now we can easily take the logarithm. Following (49) we obtain for the free energy

$`F={\displaystyle \frac{1}{\beta }}[{\displaystyle \frac{1}{2}}\text{ }\text{}{\displaystyle \frac{3g}{\mathrm{}}}\text{ }\text{}`$ (60)
$`+{\displaystyle \frac{g^2}{2\mathrm{}^2}}(72\text{ }\text{}+24\text{ }\text{})+\mathrm{}].`$
The above Feynman diagrams are of course constructed with the help of the same rules as for the imaginary-time evolution amplitude (12), (13), and (14), but instead of the Dirichlet’s Green’s function (II) we have to use the periodic Green’s function (54). We now want to evaluate the four diagrams in (60) so that we get a second-order expression for the free energy for finite temperatures. According to (53) and (59) we get for the zeroth-order contribution

$`{\displaystyle \frac{1}{2}}\text{ }\text{}=\mathrm{log}\left[{\displaystyle \frac{1}{2\mathrm{sinh}\frac{\mathrm{}\beta \omega }{2}}}\right],`$ (61)
whereas the first-order diagram becomes

$`\text{ }\text{}={\displaystyle \frac{\mathrm{}^3\beta }{4M^2\omega ^2}}\mathrm{coth}^2{\displaystyle \frac{\mathrm{}\beta \omega }{2}}.`$ (62)
The integration in (62) is trivial, because $`G^{(\mathrm{p})}(\tau ,\tau )`$ does not depend on $`\tau `$ any more according to (54). For the second order the integrations become more sophisticated:

$`\text{ }\text{}={\displaystyle \frac{\mathrm{}^5\beta \mathrm{coth}^2\frac{\mathrm{}\beta \omega }{2}}{32M^4\omega ^5\mathrm{sinh}^2\frac{\mathrm{}\beta \omega }{2}}}`$
$`\times (\mathrm{}\beta \omega +\mathrm{sinh}\mathrm{}\beta \omega ).`$
The other contribution to the second order yields

$`\text{ }\text{}={\displaystyle \frac{\mathrm{}^5\beta }{256M^4\omega ^5\mathrm{sinh}^4\frac{\mathrm{}\beta \omega }{2}}}`$ (65)
$`\times (\mathrm{sinh}2\mathrm{}\beta \omega +8\mathrm{sinh}\mathrm{}\beta \omega +6\mathrm{}\beta \omega ).`$
So all in all we get for the free energy (60) up to second order in the coupling constant $`g`$ the result (50). It shows the correct low-temperature behaviour

$`\underset{\beta \mathrm{}}{lim}F^{(2)}(\beta )={\displaystyle \frac{\mathrm{}\omega }{2}}+{\displaystyle \frac{3g\mathrm{}^2}{4M^2\omega ^2}}{\displaystyle \frac{21g^2\mathrm{}^3}{8M^4\omega ^5}},`$ (66)
which is the ground state energy and can be found for instance in Bender/Wu ; Kleinert .

## IX Variational Perturbation Theory

Variational perturbation theory is a method that enables us to resum divergent Borel-type perturbation series in such a way that they converge even for infinitely large values of the perturbative coupling Kleinert ; Frohlinde . To this end we add and subtract a trial harmonic oscillator with trial frequency $`\mathrm{\Omega }`$ to our anharmonic oscillator (2):

$`V(x)={\displaystyle \frac{M}{2}}\mathrm{\Omega }^2x^2+g{\displaystyle \frac{M}{2}}{\displaystyle \frac{\omega ^2\mathrm{\Omega }^2}{g}}x^2+gx^4.`$ (67)
Now we treat the second term as if it was of the order of the coupling constant $`g`$. The result is obtained most simply by substituting for the frequency $`\omega `$ in the original anharmonic oscillator potential (2) according to Kleinert’s square-root trick Kleinert

$`\omega \mathrm{\Omega }\sqrt{1+gr},`$ (68)
where

$`r{\displaystyle \frac{\omega ^2\mathrm{\Omega }^2}{g\mathrm{\Omega }^2}}.`$ (69)
These substitutions are not the most general ones. The square root is just a special case for the anharmonic oscillator.
We now apply the trick (68) to our first-order series representation for the free energy $`F`$ found in (50). Substituting for the frequency $`\omega `$ according to (68), expanding for fixed $`r`$ up to the first order in $`g`$ and re-substituting for $`r`$ according to (69) we get

$`F^{(1)}(\beta ,\mathrm{\Omega })={\displaystyle \frac{1}{\beta }}\mathrm{log}{\displaystyle \frac{1}{2\mathrm{sinh}\frac{\mathrm{}\beta \mathrm{\Omega }}{2}}}+{\displaystyle \frac{3g\mathrm{}^2}{4M^2\mathrm{\Omega }^2}}\mathrm{coth}^2{\displaystyle \frac{\mathrm{}\beta \mathrm{\Omega }}{2}}`$
$`+{\displaystyle \frac{\mathrm{}\mathrm{\Omega }}{4}}\left({\displaystyle \frac{\omega ^2}{\mathrm{\Omega }^2}}1\right)\mathrm{coth}{\displaystyle \frac{\mathrm{}\beta \mathrm{\Omega }}{2}}.`$ (70)
So the free energy (IX) now depends on the trial frequency $`\mathrm{\Omega }`$ which is of no physical relevance. In order to get rid of it, we have to minimize its effect by employing the principle of least sensitivity Stevenson . This principle suggests to search for local extrema of $`F(\beta ,\mathrm{\Omega })`$ with respect to $`\mathrm{\Omega }`$:

$`{\displaystyle \frac{F^{(1)}(\beta ,\mathrm{\Omega })}{\mathrm{\Omega }}}=0.`$ (71)
For the first order $`F^{(1)}(\beta ,\mathrm{\Omega })`$ it turns out that there are several extrema for each $`\beta `$. As we seek a curve $`\mathrm{\Omega }^{(1)}(\beta )`$ that is as smooth as possible the choice is easy — we take the lowest branch for the others are not bounded (see Fig. 3). Moreover the other branches lead to diverging results.
To second order, we proceed in a similar way and we find that there are no extrema at all for $`F^{(2)}(\beta ,\mathrm{\Omega })`$. In accordance with the principle of least sensitivity we look for inflection points instead, i.e. we look for solutions to the equation

$`{\displaystyle \frac{^2F^{(2)}(\beta ,\mathrm{\Omega })}{\mathrm{\Omega }^2}}=0.`$ (72)
In general we try to solve the equation

$`{\displaystyle \frac{^nF^{(N)}(\beta ,\mathrm{\Omega })}{\mathrm{\Omega }^n}}=0`$ (73)
for the smallest possible $`n`$. Plugging $`\mathrm{\Omega }^{(N)}(\beta )`$ into $`F^{(N)}(\beta ,\mathrm{\Omega })`$, we finally get back a re-summed expression for the physical quantity $`F(\beta )`$. The results for the first three orders are given in Fig. 4.

In order to check our results we have to compare them to the numerically evaluated free energy $`F_{\mathrm{num}}^{(N)}(\beta )`$ which is discussed in Sec. XI.

## X Higher Orders

We now evaluate the convergence behaviour for the variational perturbative results for the free energy $`F^{(N)}(\beta )`$ up to the fifth order. However, in order to reduce the computational cost we restrict ourselves to a certain value of the inverse temperature $`\beta `$. Results are shown in Fig. 5. For odd variational perturbation orders we optimized the free energy according to (71), i.e. we determined $`\mathrm{\Omega }`$ by setting the first derivative of $`F^{(N)}(\beta )`$ with respect to $`\mathrm{\Omega }`$ to zero. For even orders we had to go for inflection points, instead, so we had to solve equation (72).
It turns out that odd and even orders oscillate about an exponential best fit curve. For $`\beta =1`$ the numerical result, $`F_{\mathrm{num}}^{(9)}(1.0)=0.6571`$, turns out to lie within the interval obtained by fitting the five perturbative orders of the free energy with origin, as shown in Fig. 5. This interval is $`[0.657,0.659]`$, and

clearly, the variational perturbative results converge exponentially.

## XI Checking our Results

The spectral representation of the partition function reads

$`Z={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}e^{\beta E_n},`$ (74)
where the $`E_n`$ are the energy eigenvalues. Let us define the numerical approximants

$`Z_{\mathrm{num}}^{(N)}={\displaystyle \underset{n=0}{\overset{N}{}}}e^{\beta E_n}`$ (75)
and

$`F_{\mathrm{num}}^{(N)}={\displaystyle \frac{1}{\beta }}\mathrm{log}Z^{(N)},`$ (76)
respectively. One possibility to obtain numerical results for the eigenvalues $`E_n`$ is provided for by the so called “shooting method”. We integrate the Schrödinger equation numerically for the potential (2) and for a particular value of the coupling strength $`g`$. If the energy $`E`$ which we plug into the program does not coincide with one of the energy eigenvalues $`E_n`$, the solution to the Schrödinger equation explodes already for relatively small values of the coordinate $`x`$. If the energy eigenvalue is close to the exact answer, we have $`|\mathrm{\Psi }(x)|<\mathrm{}`$ also for larger values of $`x`$. This method yields the unnormalized eigenfunctions (the wave functions which still have to be normalized) and the energy eigenvalues to very high accuracy (see Tab. 1). Renormalization is necessary, for the computer algebra programme needs an initial value $`\mathrm{\Psi }(0)`$ which we set to one.

Plugging the first ten numeric energy eigenvalues

into equation (75) and evaluating (76) up to $`N=9`$, we obtain a function $`F_{\mathrm{num}}^{(N)}(\beta )`$. It converges rapidly for low temperatures, corresponding to high values of $`\beta `$. For high temperatures more terms should be taken into account.
Alternatively one can also use classical results as a high-temperature cross check: High temperatures correspond to classical statistical distributions such that we can evaluate the classical partition function according to

$`Z_{\mathrm{cl}}={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{dx}{\lambda _{\mathrm{th}}}}\mathrm{exp}\left[\beta V(x)\right],`$ (77)
with the potential (2) and $`\lambda _{\mathrm{th}}=\sqrt{2\pi \mathrm{}^2/Mk_\mathrm{B}T}`$. This integral reduces to

$`Z_{\mathrm{cl}}={\displaystyle \frac{1}{2\lambda _{\mathrm{th}}}}\sqrt{{\displaystyle \frac{M\omega ^2}{2g}}}`$
$`\times \mathrm{exp}\left({\displaystyle \frac{\beta M^2\omega ^4}{32g}}\right)K_{1/4}\left({\displaystyle \frac{\beta M^2\omega ^4}{32g}}\right),`$ (78)
where $`K_{1/4}(z)`$ is a modified Bessel function. The classical partition function (XI) can be evaluated for high temperatures which corresponds to small values of $`\beta `$. Consequently, when we test our variational perturbative results, we compare them to the classical free energy for low values of $`\beta `$, namely $`\beta <1/4`$. And for high values of $`\beta `$ we use the numerical approximation, $`F_{\mathrm{num}}^{(9)}(\beta )`$, for comparison (see Fig. 6).
In natural units $`\mathrm{}=k_\mathrm{B}=1`$ a value of $`\beta =1/4`$ corresponds to a physical temperature of $`T=1.26\times 10^6K`$.

## XII Conclusion and Outlook

The recursive technique that has been developed throughout Sections IVVI definitely out classes all diagrammatical perturbative calculations. Using the conventional evaluation of Feynman diagrams, the partition function and the free energy have been evaluated up to third order Meyer , here we obtained the fifth-order result.
For the free energy the convergence of variational perturbation theory was found to be exponential. The fact that the principle of least sensitivity Stevenson  produces extrema for the odd variational orders and inflection points for even orders is reflected in the respective convergence behaviours: Odd and even orders can best be fitted separately by exponentials as emphasized in Ref. Frohlinde . Thus we obtained intervals of convergence for certain values of the free energy which always turned out to contain the exact numerical result when taking into account the statistical errors associated with the boundaries of the intervals. For the free energy, the numerical results were obtained using its spectral representation reverting on the first ten energy eigenvalues obtained with the “shooting method”, sketched in Sec. XI.
Finally, we note that our high-order perturbative results for the anharmonic imaginary-time evolution amplitude are useful for calculating other thermodynamic quantities as the correlation function or the ground state wave function weissbach ; ref2 ; Kunihiro . Furthermore, it remains to compare these perturbative results with the semiclassical approximation ref1 .

###### Acknowledgements.

The authors wish to thank Prof. Kleinert for fruitful discussions on variational perturbation theory. FW is especially grateful to Hagen Kleinert who enthusiastically supervised his diploma thesis.


format: text

text:
# Superconductivity in hole-doped C60 from electronic correlations.

## Perturbation Theory:

Let us start by considering the following Hubbard model on a single C60 molecule,

$$H_{C_{60}}=\underset{ij,\sigma }{}t_{ij}(c_{i\sigma }^{}c_{j\sigma }+\text{h.c.})+\frac{U}{2}\underset{i,\sigma }{}n_{i\sigma }n_{i\sigma },$$
(1)
where the only non-vanishing hopping integrals are $`t_{ij}=t`$ for nearest neighbor (nn) hopping on pentagons and $`t_{ij}=t^{}`$ for nn hopping between pentagons. We use $`t=2eV`$ and $`t^{}/t=1.2`$ in accordance with the values used in  and allow $`U`$ to vary. Values of $`U512eV`$ have been suggested in the literature .

We do standard second order perturbation theory in Hubbard U. Since the Hamiltonian is spin rotationally invariant the states fall into degenerate sets corresponding to irreducible representations of the icosahedral group $`I_h`$ and spin. The states are well characterized by angular momentum and only weakly split by the icosahedral symmetry.

The validity of second order perturbation theory for the large-U Hubbard model and the neglect of longer range Coulomb interaction for this problem has been under debate. It has been shown by exact diagonalization that for small Hubbard clusters (e.g. the 12-site truncated tetrahedron) the second order perturbation theory is qualitatively correct giving positive pair-binding energies for moderately large U . In addition, longer ranged repulsions are more effectively screened by the metallic environment than the high energy second order processes that give rise to the attraction . Here, we explore the consequences of this model assuming that it gives a reasonable estimate of the molecular spectrum and the resulting pair attraction.

## Effective Interactions:

Using the perturbation theory results we derive a set of interactions for the crystal. Although the perturbative result contains terms to higher order in fermion operators, we take the effective Hamiltonian

$$H_{eff}=e_0e_1\underset{l,\sigma }{}n_{l\sigma }+t_iT_{klmn\alpha \beta \gamma \delta }^ic_{k\alpha }^{}c_{l\beta }^{}c_{m\gamma }c_{n\delta },$$
(2)
which acts on a space of one particle states in the five-dimensional $`H_u`$ representation of $`I_h`$. Here $`c_{k\alpha }^{}`$ creates a hole with quantum number $`k=1,\mathrm{},5`$ and spin $`\alpha `$, $`n_{l\sigma }`$ is the number operator and alike indices are henceforth summed over. The $`e_0`$, $`e_1`$ and $`t_i`$ are parameters and the $`T_{klmn\alpha \beta \gamma \delta }^i`$ are tensors chosen such to make the four-fermion terms $`T^i`$ invariant independently under spin and icosahedral symmetry.

Group theory reveals that there are nine such independent four-fermion terms. These can be derived by constructing all two-fermion terms $`c_{k\alpha }c_{l\beta }`$ transforming in a particular representation of spin and angular momentum and taking the tracing with their hermitian conjugates. We can write for the product of two fermions in the representation $`H`$ of $`I_h`$ and spin$`1/2`$, $`HH=(AG2H)_s+(T_1T_2G)_a`$, $`1/21/2=0_a1_s`$ where $`s`$ and $`a`$ mean the symmetric and antisymmetric parts of the tensor products and where $`A,T_1,T_2,G`$ and $`H`$ are the 1,3,3,4 and 5 dimensional representations of $`I_h`$ respectively.

The product of two anticommuting fermion operators thus reduces into seven irreducible parts, given by finding the antisymmetric part of the product of angular momentum and spin. We then construct the invariant four-fermion operators $`T^i`$, with corresponding coupling constants $`t_i`$, labeled according to from which two-fermion operators they are constructed using the composite index

$`i`$ $`=`$ $`\{(A,0,0),(H,2,0),(G,4,0),(H,4,0),(T_1,1,1),`$ (4)
$`(T_2,3,1),(G,3,1)\}.`$
Here $`(R_i,L_i,S_i)`$ indicates icosahedral representation, corresponding angular momentum in the case of full rotational symmetry, $`O(3)`$, and spin, respectively. The tensors $`T_{klmn\alpha \beta \gamma \delta }^i`$ are normalized such that they are projection operators into the $`i`$th irreducible subspace of the two-fermion products. We neglect the additional two invariants that can be constructed from tracing the two different realizations of $`H`$ together since these are not allowed under $`O(3)`$ and since the deviance from full rotational symmetry is small.

The effective Hamiltonian (2) is then used to match the spectrum found from the perturbation theory. $`e_0`$ is given by the energy of the neutral molecule (filled HOMO), $`e_1`$ by the 10-fold degenerate 1-hole states and the 2-hole spectrum is in one-one correspondence to the seven four-fermion terms with energies $`e_02e_1+2t_i`$, which fixes $`t_i=t_i(U,U^2)`$ as shown in figure 1.

In spite of its relative simplicity this effective model reproduces the trend found in the perturbative calculation, namely that for moderately large $`U`$ states with low spin and low “angular momentum” have lower energy, i.e. Hund’s rule is not valid. To make this statement more transparent we can consider the conventional four-fermion operators $`n^2`$, $`\stackrel{}{S}^2`$ and $`\stackrel{}{L}^2`$, invariant under $`O(3)\times SU(2)`$. To fit these to the normal ordered four-fermion operators $`T^i`$ we also need the tensor invariants of angular momentum $`𝒬_L^a𝒬_L^a`$ (no sum over $`L`$) where $`𝒬_L^a=Q_{L,ij}^ac_{i\alpha }^{}c_{j\alpha }`$ transform as the $`L=`$2,3 or 4 representation of $`SO(3)`$ ($`a=1,\mathrm{},2L+1`$). We find, $`n^2=n+_iT^i`$, $`\stackrel{}{S}^2=\frac{3}{4}n+_i(\frac{S_i(S_i+1)}{2}\frac{3}{4})T^i`$ and $`\stackrel{}{L}^2=6n+_i(\frac{L_i(L_i+1)}{2}6)T^i`$, with $`i`$, $`S_i`$ and $`L_i`$ as defined in (4) and with similar expressions for the other operators. Apart from the small split due to icosahedral symmetry, which we average over, we get new coupling constants as shown in fig. (2). (The $`L=2,3`$ invariants are normalized as $`\stackrel{}{L}^2`$ and due to overcompleteness we choose not to include the $`L=4`$ tensor invariant.)

We find that for large $`U`$, $`n^2`$, $`\stackrel{}{L}^2`$ and $`\stackrel{}{S}^2`$ dominate the energetics and we will subsequently only use these as a minimal model expected to capture the important physics of the interactions.

## Lattice Hamiltonian:

We can now write down the lattice Hamiltonian

$`H`$ $`=`$ $`{\displaystyle \underset{<\stackrel{}{r}\stackrel{}{r}^{}>}{}}t_{\stackrel{}{r}\stackrel{}{r}^{},kl}(c_{\stackrel{}{r},k\sigma }^{}c_{\stackrel{}{r}^{},l\sigma }+\text{h.c.})`$ (5)
$`+`$ $`{\displaystyle \underset{\stackrel{}{r}}{}}g\stackrel{}{S}_\stackrel{}{r}\stackrel{}{S}_\stackrel{}{r}+\lambda \stackrel{}{L}_\stackrel{}{r}\stackrel{}{L}_\stackrel{}{r}+un_\stackrel{}{r}^2,`$ (6)
where $`\stackrel{}{r}`$ are the sites of the lattice and $`<\stackrel{}{r}\stackrel{}{r}^{}>`$ runs over the range of intermolecular hopping. This model is quite general and could be used also for electron doped systems in two or three dimensions . Here we consider a model where the charge is confined to the  surface of the fcc crystal, i.e. a two dimensional system where $`<\stackrel{}{r}\stackrel{}{r}^{}>`$ correspond to nn on a square lattice. From the perturbative calculation (fig. 2) we take $`g,\lambda ,u>0`$. Note that the confinement of electron propagation to the surface completely breaks the five-fold degeneracy of the HOMO states. This will be manifest in the tight-binding part of the Hamiltonian which reflects the symmetry group of the surface.

Typical hopping integrals are of the order of 0.1eV  which is comparable to the interactions (fig. 2), the problem is in an intermediate coupling regime. Nevertheless, we do a standard BCS/Hartree-Fock construction, replacing the Hamiltonian (5) by a reduced non-interacting Hamiltonian. We keep only spatially uniform superconducting mean-fields $`b_{kl\alpha \beta }=\frac{1}{V}_\stackrel{}{r}<c_{\stackrel{}{r},k\alpha }c_{\stackrel{}{r},l\beta }>`$ and mean-fields of the number operators $`n_{l\alpha }=\frac{1}{V}_\stackrel{}{r}<c_{\stackrel{}{r},l\alpha }^{}c_{\stackrel{}{r},l\alpha }>`$. We can then derive the following effective Hamiltonian in momentum ($`p`$) space

$`H_{MF}`$ $`=H_0+H_{pair}+H_{HF}={\displaystyle \underset{p}{}}(t(p)_{kl}\mu \delta _{kl})c_{p,k\sigma }^{}c_{p,l\sigma }`$ (7)
$``$ $`{\displaystyle \underset{L,S}{}}(c_{p,l\beta }^{}c_{p,k\alpha }^{}𝒪_{kl\alpha \beta }^{L,S}+\text{h.c.})+h_{HF,kl}c_{p,k\alpha }^{}c_{p,l\alpha },`$ (8)
where we included a chemical potential $`\mu `$. We define the components of the order parameter

$$𝒪_{kl\alpha \beta }^{L,S}=V_{L,S}\underset{i:L_i=L,S_i=S}{}T_{klmn\alpha \beta \gamma \delta }^ib_{mn\gamma \delta },$$
(9)
where $`V_{L,S}=(\frac{3}{4}\frac{S(S+1)}{2})g+(6\frac{L(L+1)}{2})\lambda u`$ and $`V_{L,S}>0`$ corresponds to attraction.

Assuming no net magnetization the Hartee-Fock terms can be written (no sum over $`l`$)

$$h_{HF,ll}=\frac{3}{4}gn_l\lambda \tau _l+u(n_l+2\underset{kl}{}n_k),$$
(10)
where $`n_l=n_l+n_l`$ is the total particle number with angular momentum component $`l`$ and $`\stackrel{}{\tau }=(3n_4+3n_5,4n_3+n_4+n_5,4n_2+n_4+n_5,3n_1+n_2+n_3+n_5,3n_1+n_2+n_3+n_4)`$. (In addition there is one off-diagonal component $`\sqrt{3}\lambda (n_5n_4)(c_{p,1\alpha }^{}c_{p,2\alpha }+\text{h.c.})`$, which is included in our calculations but which will in general be small.)

The tight-binding part of the Hamiltonian for the simplest (unidirectional) crystal structure takes the form

$$t(p)_{kl}\left[\begin{array}{ccccc}t_1f_p& t_{12}f_p& t_{13}g_p& 0& 0\\ t_{12}f_p& t_2f_p& t_{23}g_p& 0& 0\\ t_{13}g_p& t_{23}g_p& t_3f_p& 0& 0\\ 0& 0& 0& t_4f_p& t_{45}g_p\\ 0& 0& 0& t_{45}g_p& t_5f_p\end{array}\right]$$
(11)
where $`f_p=\mathrm{cos}(p_x)\mathrm{cos}(p_y)`$ and $`g_p=\mathrm{sin}(p_x)\mathrm{sin}(p_y)`$.

We have taken hopping parameters from . For hole hopping we have, $`t_1=.107,t_2=.198,t_3=.134,t_4=.032,t_5=.170,t_{12}=.087,t_{13}=.073,t_{23}=.162,t_{45}=.115`$ eV. The Hamiltonian has the symmetry $`t(p)=t(p+(\pi ,0))`$ implying a symmetric band structure around zero energy where all bands will be half-filled and there are Van-Hove singularities at zero energy at $`(p_x,p_y)=(\pm \pi /2,\pm \pi /2)`$.

The Hartree-Fock terms (10) have interesting properties related to their multiband nature. For positive $`g`$ and $`u`$ it is energetically favorable to fill up as few bands as possible for a given particle number. The term $`\stackrel{}{S}_\stackrel{}{r}\stackrel{}{S}_\stackrel{}{r}`$ gives on-site spin-triplet states with higher energy than singlets, so that by putting particles in a single angular momentum state the energy can be lowered by exclusion, and similarly for the on-site charging energy $`n_\stackrel{}{r}^2`$, where two particles with the same spin and angular momentum cannot occupy the same site. The $`\stackrel{}{L}_\stackrel{}{r}\stackrel{}{L}_\stackrel{}{r}`$ term on the other hand gives rise to an anisotropic attraction between the components.

For positive parameters $`g,u`$ and $`\lambda `$ the $`L=0`$, $`S=0`$ pairing channel of (7) is strongest and we can expect this to dominate. But, since the rotational invariance of the C60 molecules is broken by the lattice, subdominant order parameters with non-zero angular momentum appear and in general all three $`S=0`$ order parameters may be non-zero.

What kind of physics can we expect from this model? For large $`g`$ or $`u`$ there may be Mott insulating states at even integer $`2n`$ filling where $`n`$ bands of angular momentum states will fill up completely. If $`u`$ is not very large compared to $`g`$ the insulating state will be non-magnetic due to the low energy of molecular singlets. (For the regime $`ug`$, not realized here, there may also be magnetic insulating states at odd integer filling.) The pairing terms compete with a putative insulating state due to the Hartree-Fock terms so that even for large $`g`$ or $`\lambda `$ there may be s.c. ground states also at even integer fillings, although a suppression of Tc is likely due to the low DOS when the bands are nearly filled or empty. In general we can expect the Hartree-Fock terms to completely recast the DOS compared to that given by the band structure and consequently also Tc’s.

## Results:

By numerical iteration of the mean particle number in the five angular momentum components and the s.c. mean-fields at fixed chemical potential we arrive at self-consistent solutions. For all plots the system size is $`100\times 100`$ with at least ten sampling points per unit shift in particle number. Figure (3) shows the energy gap $`2\mathrm{\Delta }`$, the norm of the s.c order parameter (defined as $`\sqrt{Tr𝒪𝒪^{}}`$ with $`𝒪=_{L,S}𝒪^{L,S}`$) at T=0 and Tc as a function of doped holes (up to 7 holes) for parameters $`u=.09,g=.06,\lambda =.02`$ eV, both with and without Hartree-Fock terms. We find that Tc scales roughly linearly with $`\mathrm{\Delta }`$ and the reduced gap $`\frac{2\mathrm{\Delta }}{T_c}3.2`$ is close to the weak coupling BCS value (3.53). For the calculation with H-F terms the magnitude of the order parameter (not in the figure) fits very well with $`2\mathrm{\Delta }`$. Without H-F terms there is a deviation from this fit around 5 holes, due to a momentum dependent gap.

Fig. 4 shows the density of states, $`\frac{n}{\mu }`$, for the same parameters and with Hartree-Fock terms and the DOS from the noninteracting band structure (11), i.e. the DOS without H-F terms. Since we find self-consistent solutions of both the band fillings and the gaps we calculate the DOS at finite temperature, above Tc , for the realization with H-F terms.

The values of $`\lambda `$ and $`u`$ chosen here correspond roughly to the perturbative spectrum (fig. 2) at $`U=5`$, but we have reduced $`g`$ significantly in order for the pairing attraction not to dominate completely. One could argue that the second order perturbation theory can be expected to overestimate the core polarization effect that gives the $`\stackrel{}{S}^2`$ term. Of course, the actual magnitude of Tc and the gap that we find should not be taken too literally since parameters are only rough estimates and we are doing mean-field theory at relatively strong coupling. The important result is the qualitative behavior of Tc as a function of doping.

Without Hartree-Fock terms, Tc follows the DOS of the band structure with a corresponding maximum at the half-filled band. This should be contrasted with the results for the full Hamiltonian which has a Tc that is maximized close to three holes and which is strongly suppressed at five holes. The DOS is suppressed at four and six holes due to the commensurate lock-in discussed above. In fact, the influence of these special fillings is such that the DOS is low in the whole region between four and six holes where it is the highest without interactions.

We thank S.A. Kivelson for valuable comments. This research was supported by STINT and the Swedish Research Council.