<figure><img src="image_1.png" /></figure>

HAL
open science

# Path integral for the quartic oscillator: An accurate analytic formula for the partition function

Michel Caffarel

<figure><img src="image_2.png" /></figure>

## To cite this version:

Michel Caffarel. Path integral for the quartic oscillator: An accurate analytic formula for the partition function. Journal of Mathematical Chemistry, 2024, 10.1007/s10910-024-01671-z. http://dx.doi.org/10.1007/s10910-024-01671-z

HAL Id: hal-04663791

https://hal.science/hal-04663791v1

Submitted on 29 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Path integral for the quartic oscillator: An accurate analytic formula for the partition function

Michel Caffarel¹, *

¹Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, France

In this work an approximate analytic expression for the quantum partition function of the quartic oscillator described by the potential $V(x) = \frac{1}{2}\omega^2x^2 + gx^4$ is presented. Using a path integral formalism, the exact partition function is approximated by the partition function of a harmonic oscillator with an effective frequency depending both on the temperature and coupling constant $g$. By invoking a Principle of Minimal Sensitivity (PMS) of the path integral to the effective frequency, we derive a mathematically well-defined analytic formula for the partition function. Quite remarkably, the formula reproduces qualitatively and quantitatively the key features of the exact partition function. The free energy is accurate to a few percent over the entire range of temperatures and coupling strengths $g$. Both the harmonic ($g \to 0$) and classical (high-temperature) limits are exactly recovered. The divergence of the power series of the ground-state energy at weak coupling, characterized by a factorial growth of the perturbational energies, is reproduced as well as the functional form of the strong-coupling expansion along with accurate coefficients. Explicit accurate expressions for the ground- and first-excited state energies, $E_0(g)$ and $E_1(g)$ are also presented.

## I. INTRODUCTION

In this work we consider the familiar one-dimensional anharmonic quartic oscillator described by the Hamiltonian

$$H = -\frac{1}{2}\frac{d^2}{dx^2} + \frac{\omega^2}{2}x^2 + gx^4, \quad (1)$$

where $\omega^2$ and $g \ge 0$ denote the harmonic force and coupling constants, respectively. The quartic oscillator is one of the simplest yet non-trivial system of quantum mechanics. It is widely used as a theoretical model in many fields including quantum chemistry (anharmonic vibrational effects in molecular spectroscopy), solid-state physics, laser theory, nuclear physics, and quantum field theory. The quartic oscillator has led to a vast literature from a mathematical, numerical, and physical point of view. A particularly interesting feature of this elementary model is the divergence of the Rayleigh-Schrödinger perturbation series of the ground-state energy for all $g > 0$. Understanding the origin of this divergence and establishing efficient resummation techniques for this paradigmatic divergent series is at the heart of most studies. Among these, probably the most influential one is that of Bender and Wu[1, 2] who made a detailed mathematical/numerical analysis of the model. The authors were able to answer a number of important questions and, in particular, to shed some light on the origin of the divergence by investigating the location of the singularities of the energy in the complex $g$-plane. Many others works have followed, as referenced below. With respect to the divergence issue, virtually all known resummation methods have been applied to this series. The two predominant approaches

are the use of Padé approximants[3–5] and Borel’s integral summation method.[6] However, various alternative approaches have been developed, including the use of hypergeometric functions ([7] and references therein), nonlinear sequence of transformations,[8, 9] Borel transformation with a conformal mapping,[10] or an order-dependent mapping,[11], as well as sequences of analytic approximations[12], among others.

In this study, we are interested in evaluating the partition function (PF) of the quartic oscillator defined as

$$Z = \mathrm{Tr}e^{-\beta H} = \sum_{n=0}^{+\infty} e^{-\beta E_n} \quad (2)$$

where $E_n$ are the discretized energies of the system and $\beta$ the inverse temperature. A straightforward and direct method for evaluating the partition function consists in adding up the exponential components (Boltzmann weights) of Eq.(2), see e.g. [13, 14]. However, the exact energies being not known, approximate numerical energies are required, for example by diagonalizing the Hamiltonian matrix within a sufficiently large basis set of Gaussian functions.[15] In view of the simplicity of the model, highly accurate energies can be obtained, enabling the computation of an "exact" numerical PF, at least for not too small values of $\beta$. In the following, we shall employ these precise numerical values as our "exact" reference. Among the analytic methods for evaluating the PF, the oldest one is certainly the Wigner-Kirkwood perturbation expansion of the PF in powers of $\hbar$ (or inverse temperature $\beta$), a method that systematically evaluates the quantum corrections to the classical partition function.[16–18] This approach has been improved and extended in different ways, e.g. [19],[20]. Another approach is the thermodynamic variation perturbation method[21, 22] based on a Schwinger-type expansion of the partition function.[23] Furthermore, exact upper and lower bounds for the PF have been obtained.[24–26] A natural framework to work out approximations for the

<footnote>* caffarel@irsamc.ups-tlse.fr</footnote>
PF is the path integral formalism as used in the present work. The most prominent approach following this route was pioneered by Feynman[27] and then refined by Feynman and Kleinert[28] (also, independently, by Giachetti and Tognetti[29]). Over the years, the method of Feynman and Kleinert has been systematically improved and extended by Kleinert and collaborators[30–33].

The main result of this work is to derive a closed-form expression for the partition function using a path integral formalism. Although approximate, it captures some of the key features of the exact partition function for all temperatures and coupling constants (from the weak- to the strong- coupling regimes). From the partition function we also derive explicit expressions for the ground- and first-excited state energies $E_0(g)$ and $E_1(g)$, respectively. From a general point of view, our simple partition function reproducing the essential features of the exact partition function provides an interesting model to investigate the properties of the quartic oscillator.

This paper is organized as follows. In Section II we first derive an exact path integral expression for the PF. Then, it is shown that at zero coupling the partition function of the harmonic oscillator is recovered. In Section III a gaussian approximation for the probability density appearing in the path integral is introduced. In practice, this amounts to approximating the exact partition function by that of a harmonic oscillator with an effective frequency depending both on the temperature and coupling constant $g$. Then, we propose in Section IV to invoke a Principle of Minimal Sensitivity (PMS) of the PF to the effective frequency. Remarkably, this principle leads to a mathematical constraint allowing to define a mathematically well-founded partition function, which is our final expression. Section V is devoted to the derivation of explicit expressions for $E_0(g)$ and $E_1(g)$, and to the calculation of the coefficients of the weak- and strong-coupling expansions of the ground-state energy. In Section VI we present a comparative study between our formula and the partition functions of Feynman and Kleinert[28] and of Büttner and Flytzanis[34]. Finally, in Section VII, a summary of our main results and a few perspectives are presented.

## II. PARTITION FUNCTION AS A PATH INTEGRAL

In the position representation the partition function writes

$$Z = \int dx \langle x | e^{-\beta H} | x \rangle. \quad (3)$$

To obtain its path integral representation, we follow the standard route (see, e.g., [35]): The exponential operator is broken into a product of $n$ exponential operators as $e^{-\beta H} = \prod_{i=1}^{n} e^{-\tau H}$ with $\tau = \frac{\beta}{n}$, and the spectral

resolution of the identity operator, $\mathbb{1} = \int dx_i |x_i\rangle\langle x_i|$, is introduced between each operator giving

$$Z = \int dx_1 \dots \int dx_n \prod_{i=1}^{n} \langle x_i | e^{-\tau H} | x_{i+1} \rangle \quad (4)$$

where the initial and final points are identified to $x$, $x_1 = x_{n+1} = x$. To proceed, we introduce a high-temperature (or small $\tau$) approximation of the quantity $\langle x_i | e^{-\tau H} | x_{i+1} \rangle$ (known as the propagator or Green’s function). At the lowest order in $\tau$, we have

$$\langle x_i | e^{-\tau H} | x_{i+1} \rangle = \frac{1}{\sqrt{2\pi\tau}} e^{-\frac{(x_{i+1}-x_i)^2}{2\tau}} e^{-\tau V(x_i)} + O(\tau^2) \quad (5)$$

where $V(x)$ is the potential, here $V(x) = \frac{\omega^2}{2}x^2 + gx^4$. By taking the large-$n$ limit the contribution of the quadratic error vanishes, leading to the following exact path integral expression of the partition function

$$Z = \lim_{n \to \infty} Z_n \quad (6)$$

where the discretized partition function $Z_n$ writes

$$Z_n = \left( \frac{1}{\sqrt{2\pi\tau}} \right)^n \int dx_1 \dots \int dx_n \quad (7)$$

$$e^{-\frac{1}{\tau} \sum_{i,j=1}^{n} x_i A_{ij} x_j - \tau \sum_{i=1}^{n} V(x_i)}. \quad (8)$$

Here, the matrix $A$ is given by

$$A_{ij} = \delta_{ij} - \frac{1}{2}(\delta_{i i+1} + \delta_{i-1 i}) \quad (9)$$

with the boundary conditions $A_{1n} = A_{n1} = -\frac{1}{2}$. In the following, for the sake of convenience, the $n$-dimensional integrals will be denoted as $\int dx$, where $x = (x_1, \dots, x_n)$.

The matrix $A$ can be diagonalized by using a Fourier transform, we get

$$Z_n = \left( \frac{1}{\sqrt{2\pi\tau}} \right)^n \int dx e^{-\frac{1}{\tau} \sum_{i=1}^{n} \lambda_i \tilde{x}_i^2} e^{-\tau \sum_{i=1}^{n} V(x_i)} \quad (10)$$

where $\lambda_i$ are the eigenvalues of $A$. For $n > 2$ the $\lambda_i$'s are given by

$$\lambda_i = 1 - \cos \frac{2\pi}{n}(i-1) \quad i=1 \text{ to } n. \quad (11)$$

The eigenvectors will be denoted by $\tilde{x}_i$ and decomposed as

$$\tilde{x}_i = \sum_{j=1}^{n} O_{ij} x_j \quad (12)$$

where $O_{ij}$ is the orthogonal matrix diagonalizing $A$. The orthogonality condition writes

$$\sum_{k=1}^{n} O_{ik} O_{jk} = \delta_{ij}. \quad (13)$$
Let us now introduce the probability density function
defined as

$$
\pi(x) = \frac{e^{-\tau V(x)}}{I_v} \qquad (14)
$$

where $I_v$ is the normalization factor

$$
I_v = \int dx e^{-\tau V(x)}. \qquad (15)
$$

The discretized partition function can be written as

$$
Z_n = \left( \frac{I_v}{\sqrt{2\pi\tau}} \right)^n \int dx \prod_{i=1}^{n} \pi(x_i) e^{-\frac{\lambda_i}{\tau} \hat{x}_i^2(x)}. \quad (16)
$$

At this point, no approximation has been made. As $n \to \infty$, $Z_n$ converges to the exact partition function. Unfortunately, for an arbitrary density $\pi(x)$ [equivalently, an arbitrary potential $V(x)$], the multi-dimensional integrals cannot be performed analytically and approximations are to be introduced. This will be the subject of the two following sections.

Before doing this, let us first verify that the harmonic
partition function, denoted here as $Z^{(0)}$, is recovered at
zero coupling. $Z^{(0)}$ is given by

$$
Z^{(0)} = \sum_{k=0}^{+\infty} e^{-\beta(k+\frac{1}{2})\omega} = \frac{1}{e^{\frac{\beta\omega}{2}} - e^{-\frac{\beta\omega}{2}}}. \quad (17)
$$

Besides checking the validity of the expression for $Z_n$ in
a particular case, the derivation of $Z^{(0)}$ to follow will
also be of interest for the next section.

At $g=0$ the probability density $\pi(x)$, Eq.(14), is gaus-
sian and writes

$$
\pi(x) = \frac{e^{-\tau \frac{1}{2} \omega^2 x^2}}{\sqrt{\frac{2\pi}{\tau \omega^2}}} \quad (18)
$$

Using the orthogonality of the matrix O, Eq.(13), leading
to $\sum_{i=1}^n \hat{x}_i^2 = \sum_{i=1}^n x_i^2$. $Z_n$ can be written as a product
of one-dimensional gaussian integrals. Performing the
gaussian x-integration we get

$$
Z_n^{(0)} = \left( \frac{\omega \sqrt{\tau} I_v^0}{\sqrt{2\pi}} \right)^n \left( \frac{1}{\omega \tau} \right) \left[ \prod_{k=2}^n \frac{1}{\sqrt{2\lambda_k}} \right] \prod_{i=2}^n \frac{1}{\sqrt{1 + \frac{\omega^2 \beta^2}{2\lambda_i n^2}}} \quad (19)
$$

Here, $I_v^{(0)}$ is equal to $\sqrt{\frac{2\pi}{\tau\omega^2}}$, leading the first factor of the RHS of the equation to be one. The quantity given in brackets is equal to $\frac{1}{n}$, a result which follows from the relation (see, its derivation in Appendix A)

$$
\prod_{k=1}^{n-1} \sin \frac{k\pi}{n} = \frac{n}{2^{n-1}} \quad (20)
$$

which leads to

$$
\prod_{k=2}^{n} \frac{1}{\sqrt{2\lambda_k}} = \prod_{k=1}^{n-1} \frac{1}{2 \sin \frac{k\pi}{n}} = \frac{1}{n}. \quad (21)
$$

Now, let us introduce the function $P_n(x)$ defined as

$$
P_n(x) = \prod_{i=2}^{n} \frac{1}{\sqrt{1 + \frac{x^2}{2\lambda_i n^2}}}. \qquad (22)
$$

Using $P_n$ and the preceding results, the partition function
can be written as

$$
Z_n^{(0)} = \frac{P_n(\beta\omega)}{\beta\omega}. \qquad (23)
$$

Evaluating $P_n(x)$ as $n \to \infty$ is not straightforward and
requires some algebra. In Appendix B we show that

$$
\lim_{n \to \infty} P_n(x) = \frac{x}{e^{\frac{x}{2}} - e^{-\frac{x}{2}}}. \qquad (24)
$$

Finally,

$$
Z^{(0)} = \lim_{n \to \infty} Z_n^{(0)} = \frac{1}{e^{\frac{\beta\omega}{2}} - e^{-\frac{\beta\omega}{2}}}, \quad (25)
$$

in agreement with Eq.(17). Note that the derivation of
the partition function of the harmonic oscillator and,
more generally, of the Green’s function $K(x, x', \beta) =$
$\langle x|e^{-\beta\hat{H}}|x' \rangle$, has been done in the literature in many dif-
ferent ways, see for example [35–41]. The present deriva-
tion is related to that of Cohen[38].

III. GAUSSIAN APPROXIMATION

The first approximation introduced here for evaluating
$Z_n$, Eq.(16), consists in approximating the general den-
sity $\pi(x) = \frac{e^{-\tau V(x)}}{\int dx e^{-\tau V(x)}}$ by a gaussian density, $\pi_G$, corre-
sponding to an effective harmonic oscillator of frequency
$\omega_g(\tau)$

$$
\pi(x) \sim \pi_G(x) = \frac{e^{-\tau \frac{1}{2} \omega_g^2(\tau) x^2}}{\int dx e^{-\tau \frac{1}{2} \omega_g^2(\tau) x^2}} \quad (26)
$$

Note that when the coupling constant *g* goes to zero,
the effective frequency reduces to *ω*. Different criteria
can be chosen to define the optimal effective frequency
minimizing the error in the approximation. Here, we
propose to impose the variance of *π*<sub>G</sub> to be equal to the
exact one, that is

$$
\sigma^2(\tau) = \langle x^2 \rangle_{\pi_G} = \langle x^2 \rangle_{\pi}. \qquad (27)
$$

After some manipulations, this equality leads to

$$
\omega_g(\tau) = \omega \sqrt{B\left(\frac{4g}{\tau\omega^4}\right)} \quad (28)
$$
where the parameter-free function $B(x)$ is defined as

$$B(x) = \frac{1}{2} \int dy e^{-y^2 - xy^4} \qquad (29)$$

a function which can be expressed as

$$B(x) = \frac{4x K_{\frac{1}{4}}\left(\frac{1}{8x}\right)}{K_{-\frac{3}{4}}\left(\frac{1}{8x}\right) + K_{\frac{5}{4}}\left(\frac{1}{8x}\right) - 2(1+2x)K_{\frac{1}{4}}\left(\frac{1}{8x}\right)} \quad (30)$$

where $K_ν(x)$ is the modified Bessel function of the second kind. The function $B(x)$ is positive and monotonically increasing. It starts at $x = 0$ with $B(0) = 1$, then increases linearly at small $x$'s and, finally, behaves as $\sim \sqrt{x}$ at large $x$.

The gaussian approximation being made, $Z_n$ becomes the partition function of a harmonic oscillator

$$Z_n = \left( \frac{I_v}{\sqrt{2\pi\tau}} \right)^n \int dx \prod_{j=1}^{n} \frac{e^{-\frac{\tau\omega_g^2(\tau)x_j^2}{2}}}{\int dx_j e^{-\frac{\tau\omega_g^2(\tau)x_j^2}{2}}} \quad (31)$$

which can be evaluated as done in the previous section for the usual harmonic oscillator. After the $x$ gaussian integration, $Z_n$ takes the form of Eq.(19), with the replacements

$$I_v^0 \to I_v \quad \text{and} \quad \omega \to \omega_g(\tau). \qquad (32)$$

The partition function is then given by

$$Z_n = \left( \frac{\omega_g(\tau)\sqrt{\tau}I_v}{\sqrt{2\pi}} \right)^n \frac{1}{\beta\omega_g(\tau)} P_n[\beta\omega_g(\tau)]. \quad (33)$$

Finally, by performing the limit $n \to \infty$ for $P_n[\beta\omega_g(\tau)]$ only, we are led to the following expression for $Z_n$

$$Z_n = \left( \frac{\omega_g(\tau)\sqrt{\tau}I_v}{\sqrt{2\pi}} \right)^n \frac{1}{e^{\frac{\beta\omega_g(\tau)}{2}} - e^{-\frac{\beta\omega_g(\tau)}{2}}} \quad (34)$$

Note that, since $P_n$ has been replaced by $P_\infty$, this formula is now only valid in the large-$n$ limit. In the harmonic case, the quantity $C(\tau) \equiv \frac{\omega_g(\tau)\sqrt{\tau}I_v}{\sqrt{2\pi}}$ is equal to one and $\omega_g(\tau) = \omega$. $Z_n$ thus becomes independent of $n$ and the infinite-$n$ limit becomes trivial. At non-zero coupling, this is no longer the case. In the large-$n$ limit, we have $\omega_g(\tau) \sim \tau^{-1/4}$ and $I_v \sim \tau^{-1/4}$. As a consequence, in the small-$\tau$ limit $C(\tau)$ converges to a finite constant independent of $\omega$ and $g$, given by $C(0) = 2\sqrt{\frac{2}{\pi}}\frac{\Gamma(\frac{5}{4})}{\Gamma(\frac{1}{4})} = 1.24397...$. This constant being greater than 1, the partition function diverges geometrically as $Z_n \sim C(0)^n e^{-\alpha n^{1/4}}$ with $\alpha > 0$

$$Z_n \to +\infty \quad \text{as} \quad n \to +\infty. \qquad (35)$$

In the next section we propose to overcome this problem by introducing an additional constraint to $Z_n$.

## IV. PRINCIPLE OF MINIMAL SENSITIVITY

In the preceding section, the exact partition function has been approximated by the partition function of a harmonic oscillator with an effective frequency $\omega_g(\tau)$. Now, the exact partition function being independent of $\omega_g(\tau)$ it is natural to impose a minimal sensitivity of $Z_n$ to the choice of this frequency. Mathematically, this condition writes

$$\frac{\partial Z_n}{\partial \omega_g(\tau)} = 0. \qquad (36)$$

Using the expression of $Z_n$ given by Eq.(34), the stationarity condition leads to the following implicit equation for the variable $n$

$$n = \frac{\beta\omega_g\left(\frac{\beta}{n}\right)}{2}\coth\left(\frac{\beta\omega_g\left(\frac{\beta}{n}\right)}{2}\right). \qquad (37)$$

As shown below, this equation actually admits a unique solution denoted, here, as $n = n_c(\beta)$ (in general, not an integer). We thus see that imposing to the gaussian approximation of the exact $Z_n$ both a minimal error in the approximation and a minimal sensitivity to the choice of the effective frequency turns out to be possible only for a single value of $n = n_c(\beta)$. We thus propose to define $Z$ as

$$Z = Z_{n_c(\beta)}. \qquad (38)$$

The unicity of the solution is shown as follows. Let us rewrite Eq.(37) as $f(n) = g(n)$ where $f(n) = n$ and $g(n) = \frac{\beta\omega_g(\frac{\beta}{n})}{2}\coth\frac{\beta\omega_g(\frac{\beta}{n})}{2}$. Both functions are monotonically increasing functions of $n$. In the case of the function $g(n)$ this is true because $\omega_g(\tau)$ is an increasing function of $n$. At large $n$, $f(n)$ increases linearly, while $g(n)$ increases slower as $\sim n^{1/4}$. Now, the function $f$ starting with a smaller value than $g$, $f(0) = 0 < g(0)$, and increasing more rapidly, the two functions necessarily cross at a unique value of $n$.

To summarize, our final formula for the partition function of the quartic oscillator is given as follows. Introducing the convenient $\beta$-dependent effective $\tau$ denoted as $\tau_c(\beta)$ and defined as

$$\tau_c(\beta) = \frac{\beta}{n_c(\beta)}, \qquad (39)$$

the partition function writes

$$Z = \left( \frac{\omega_g[\tau_c(\beta)]\sqrt{\tau_c(\beta)}I_v[\tau_c(\beta)]}{\sqrt{2\pi}} \right)^{n_c(\beta)} \frac{1}{e^{\frac{\beta\omega_g[\tau_c(\beta)]}{2}} - e^{-\frac{\beta\omega_g[\tau_c(\beta)]}{2}}} \quad (40)$$

with

$$I_v[\tau_c(\beta)] = \int dx e^{-\tau_c(\beta)V(x)}. \qquad (41)$$
The quantities $n_c(\beta)$ and $\tau_c(\beta) = \frac{\beta}{n_c(\beta)}$ are obtained by solving (for example, iteratively) the pair of equations

$$n_c(\beta) = \frac{\beta \omega_g[\tau_c(\beta)]}{2} \coth \frac{\beta \omega_g[\tau_c(\beta)]}{2} \quad (42)$$

and

$$\omega_g[\tau_c(\beta)] = \omega \sqrt{B \left[ \frac{4g}{\tau_c(\beta)\omega^4} \right]} \quad (43)$$

Note that, besides reducing to the partition function of the harmonic oscillator at zero coupling, our expression also reduces to the classical partition function in the high-temperature limit, $\beta \to 0$. Indeed, in this limit we have $n_c(\beta) \to 1$, $\tau_c(\beta) \to \beta$, and $\omega_g(\beta) \to \beta^{-1/4}$ leading to the exact classical partition function

$$Z \to \frac{1}{\sqrt{2\pi\beta}} \int dx e^{-\beta V(x)} \quad (44)$$

Introducing the principle of minimal sensitivity being a critical step of this work, it is worth presenting some numerical calculations illustrating this principle. In figure 1 the partition function $Z_n$ given by Eq.(34), that is, before the principle of minimal sensitivity has been introduced, is shown as a function of $n$ for three different values of the inverse temperature, $\beta = 5, 7.5$, and $10$ and for $\omega = g = 1$. For each $\beta$ and $n$, $Z_n$ is evaluated for thirteen different values of the effective frequency $\omega_g$. For each $n$, the $\omega_g$'s have been chosen to be uniformly distributed around $\omega_g(\tau)$, the optimal effective frequency leading to a minimal error in the gaussian approximation, as given by Eq.(28). Figure 1 illustrates the fact that $Z_n$ diverges at large $n$, Eq.(35). We also see that, for each $n$, $Z_n$ is strongly sensitive to the value of $\omega_g$ used, except at a unique value of $n [=n_c(\beta)]$ where the stationarity condition, $\frac{\partial Z_n}{\partial \omega_g(\tau)} = 0$ holds. The values of $n_c$ observed on the figure coincide with those obtained by solving the pair of equations, Eqs.(42) and (43), as it should. These various results, discussed here only for $g = 1$ and three different $\beta$'s, have been found to be valid for any value of $\beta$ and $g$ (data not shown here). In figure 2 we compare the partition function $Z$ as a function of $\beta$ and the exact one obtained by explicit summation of the exponential components of $Z$ using highly accurate numerical energies (solid line in the figure). As seen, at the scale of the figure, both sets of data are nearly indistinguishable. To increase the resolution, the inset presents the difference $Z(\beta) - Z_{ex}(\beta)$. The maximum error is about 0.01. As expected, the error vanishes in the high-temperature limit, $\beta \to 0$.

## V. ENERGY

In this section the formulas for the ground- and first-excited energies, $E_0(g)$ and $E_1(g)$ are derived. Next, the coefficients of the weak- and strong-coupling expansions of the ground-state energy are evaluated and compared to the exact ones.

<figure><img src="image_3.png" /><figcaption>Figure 1. $Z_n$, Eq.(34), as a function of $n$ for $\omega = 1$, $g = 1$, and $\beta = 5, 7.5$, and $10$. For each $n$, $Z_n$ is evaluated for thirteen different values of $\omega_g$'s uniformly distributed around $\omega_g(\tau)$ as given by Eq.(28).</figcaption></figure>

### A. Ground- and first-excited energies

The ground-state energy is obtained from the zero-temperature limit of the free energy as follows

$$E_0(g) = \lim_{\beta \to \infty} -\frac{1}{\beta} \ln Z. \quad (45)$$

At large $\beta$, the (unique) solution of Eqs.(42) and (43) is found to be proportional to $\beta$

$$n_c(\beta) = \frac{\bar{\omega}_g}{2} \beta, \quad (46)$$

where $\bar{\omega}_g$ denotes the solution of the implicit equation given by

$$\bar{\omega}_g = \omega \sqrt{B \left( \frac{2g\bar{\omega}_g}{\omega^4} \right)} \quad (47)$$

Note that the proportionality of $n_c(\beta)$ with $\beta$ is consistent with the data presented in Fig.1. Using the zero-temperature limit of the free energy and the expression
<figure><img src="image_4.png" /><figcaption>Figure 2. Comparison between $Z(\beta)$ and the exact numerical partition function $Z_{\text{ex}}(\beta)$. The inset shows the error. Here, $\omega = 1$ and $g = 1$.</figcaption></figure>

for the partition function, we get

$$E_0(g) = \frac{\bar{\omega}_g}{2} \left[ 1 - \ln \left( \sqrt{\frac{\bar{\omega}_g}{\pi}} \bar{I}(g) \right) \right] \quad (48)$$

where

$$\bar{I}(g) = \int dx e^{-\frac{2}{\bar{\omega}_g} V(x)}. \quad (49)$$

To get the first-excited energy $E_1(g)$, the first subleading contribution to $Z(\beta)$ at large $\beta$ has to be evaluated. For that, we first rewrite the equation obeyed by $\omega_g(\beta)$ under the form

$$\omega_g(\beta) = \omega \sqrt{B} \left[ \frac{2g}{\omega^4} \omega_g(\beta) \coth \left( \frac{\beta \omega_g(\beta)}{2} \right) \right]. \quad (50)$$

Introducing the variable $y$ defined as

$$y = e^{-\beta \omega_g} \quad (51)$$

the hyperbolic cotangent can be expanded for all values of $y < 1$ as

$$\coth \frac{\beta \omega_g(\beta)}{2} = 1 + 2 \sum_{n=1}^{\infty} y^n. \quad (52)$$

As a consequence of Eq. (50) and of the fact that the function $B(x)$ is infinitely differentiable, the effective frequency can also be expanded in powers of $y$. In addition, at the lowest order $y \sim e^{-\beta \bar{\omega}_g}$ up to an exponentially small correction. By performing the Taylor expansion of $B(x)$ at the value $\frac{2g}{\bar{\omega}_g^4} \bar{\omega}_g$ we get at the lowest order

$$\omega_g(\beta) = \bar{\omega}_g + \bar{c}_g e^{-\beta \bar{\omega}_g} + O[e^{-2\beta \bar{\omega}_g}] \quad (53)$$

with

$$\bar{c}_g = \frac{2u_g}{1 - \frac{u_g}{\bar{\omega}_g}} \quad (54)$$

and

$$u_g = \frac{g}{\omega^2} B' \left( \frac{2g}{\omega^4} \bar{\omega}_g \right). \quad (55)$$

Incorporating the first-order expression for $\omega_g(\beta)$ into the expression of the partition function, it is also possible to expand the PF with respect to $y$. After some algebra the two first leading contributions to the PF write

$$Z(\beta) = e^{-\beta E_0} + e^{-\beta(E_0 + \bar{\omega}_g)} [1 - \beta \Delta(g)] + O\left[e^{-\beta(E_0 + 2\bar{\omega}_g)}\right] \quad (56)$$

with

$$\Delta(g) = \left(1 + \frac{\bar{c}_g}{2\bar{\omega}_g}\right) \left[1 - \frac{4\bar{V}(g)}{\bar{\omega}_g \bar{I}(g)} - 2 \ln\left(\sqrt{\frac{\bar{\omega}_g}{\pi}} \bar{I}(g)\right)\right] \quad (57)$$

and

$$\bar{V}(g) = \int dx V(x) e^{-\frac{2}{\bar{\omega}_g} V(x)}. \quad (58)$$

Due to the presence of the linear term $[1 - \beta \Delta(g)]$ in $Z(\beta)$, our approximate partition function cannot be expressed as an infinite sum of Boltzmann factors as it should be for the exact one. As a consequence, in order to define what is meant by excited-state energies in our model, we need to introduce some prescription. Here, we propose to exponentiate the linear contribution $[1 - \beta \Delta(g)]$ to give $e^{-\beta \Delta(g)}$. By doing this, the first subleading contribution to $Z(\beta)$ becomes a pure exponential and the first excited-state energy is defined as

$$E_1(g) = E_0(g) + \bar{\omega}_g + \Delta(g). \quad (59)$$

Note that, in practice, the exponentiation has a marginal quantitative impact since the product $\beta \Delta(g)$ remains small for typical values of $\beta$. Preliminary investigations show that a generalization of such a strategy for evaluating higher excited-state energies is possible. However, the general structure is not so simple. This study is left for future research.

In figure 3 the variation of $E_0$ and $E_1$ as a function of $g$ for $\omega = 1$ is shown. The exact curves are represented by solid lines. At the scale of the figure, $E_0(g)$ is in excellent agreement with the exact energy over the full range of $g$. The first excited-energy is slightly less accurate and has also the correct overall behavior. Note that both $E_0$ and $E_1$ are systematically smaller than the exact ones for all values of the coupling constant.
<figure><img src="image_5.png" /><figcaption>Figure 3. $E_0$ and $E_1$ as a function of $g$. Exact results given by the solid lines. $\omega = 1$.</figcaption></figure>

## B. Weak coupling expansion

For small values of $g$ (weak-coupling) the Rayleigh-Schrödinger perturbation theory expresses the ground-state energy as a power series

$$E_0(g) = \sum_{n=0}^{+\infty} E^{(n)} g^n. \qquad (60)$$

The $n$-th order perturbational energies $E^{(n)}$ are known to be rational numbers. For $\omega = 1$ they are given by $\frac{1}{2}, \frac{3}{4}, -\frac{21}{8}, \frac{333}{16}, -\frac{30885}{128}, ...$ In the large-$n$ regime they have the following asymptotic behavior[1]

$$E^{(n)} \sim_{n \to \infty} -\frac{\sqrt{6}}{\pi} \left(\frac{-3}{4}\right)^n \frac{(2n)!}{n!}. \qquad (61)$$

Because of the factorial increase of the coefficients the series has a zero radius of convergence. In practice, by summing a finite number of components of this divergent series the ground-state energy can be accurately obtained only for (very) small values of $g$. For larger values, it is necessary to use one of the methods to sum up divergent series, see references in the introduction.

Let us now evaluate the weak-coupling expansion of our approximate ground-state energy $E_0(g)$, Eqs.(48) and (49). At zero coupling, $\bar{\omega}_g = \omega$, and the zero-th order of the expansion reduces to the ground-state energy of the harmonic oscillator,

$$E^{(0)} = \frac{\omega}{2}. \qquad (62)$$

The ground-state energy depending only on $\bar{\omega}_g$, let us expand it in powers of $g$

$$\bar{\omega}_g = \sum_{n=0}^{+\infty} u_n g^n \quad \text{with} \quad u_0 = \omega. \qquad (63)$$

Using the implicit equation determining $\bar{\omega}_g$, Eq.(47), the coefficients $u_n$ can be expressed in terms of the coefficients $B_n$ of the polynomial expansion of $B(x)$, Eq. (29). Now, let us define the function $J(g)$ as

$$J(g) = \sqrt{\frac{\bar{\omega}_g}{\pi}} \bar{I}(g) = \frac{1}{\sqrt{\pi}} \frac{\bar{\omega}_g}{\omega} \int dx e^{-x^2} e^{-\frac{2g\bar{\omega}_g}{\omega^4}x^4} \qquad (64)$$

and expand it in power series of $g$

$$J(g) = \sum_{n=0}^{+\infty} J_n g^n \quad \text{with} \quad J_0 = 1. \qquad (65)$$

The coefficients are given by

$$J_n = \frac{1}{\sqrt{\pi}} \frac{\bar{\omega}_g^{n+1}}{\omega^{4n+1}} \frac{(-2)^n [4n]}{n!} \qquad (66)$$

where the symbol $[4n]$ is defined as

$$[4n] = \int dx e^{-x^2} x^{4n} = \frac{(4n-1)!!}{2^{2n}} \sqrt{\pi}. \qquad (67)$$

Using the expansion of $\bar{\omega}_g$ the coefficients $J_n$ can be expressed in terms of the coefficients $\mathbf{u} = (u_1, u_2, ...)$.

$$E_0(g) = \frac{1}{2} \left( \sum_{n=0}^{+\infty} u_n g^n \right) \left[ 1 - \ln \left( 1 + \sum_{n=1}^{+\infty} J_n[\mathbf{u}] g^n \right) \right]. \qquad (68)$$

As seen, the power expansion of the energy can be expressed only in terms of the coefficients $\mathbf{u}$. The evaluation of these coefficients is rather tedious but can be easily performed using a symbolic computation program such as Mathematica[42]. In Table I we report the exact and approximate coefficients of the power series up to $n = 5$. Quite interestingly, the approximate coefficients are found to be rational numbers like in the exact case, the two first coefficients being identical. Finally, we note that the two series display a similar mathematical pattern. Indeed, the magnitude of the coefficients increases dramatically with the perturbational order in both cases (factorially, for the exact series) with an even more rapid rate for the approximate series and, furthermore, the sign pattern of the coefficients is identical.
<table><caption>Table I. Coefficients of the power series of $E_0(g)$ up to $n=5$. The first column reports the exact coefficients. The second column the coefficients obtained in the present work.</caption><thead><tr><th>n</th><th>E<sub>ex</sub><sup>(n)</sup></th><th>E<sup>(n)</sup></th></tr></thead><tbody><tr><td>0</td><td>1/2= 0.5</td><td>1/2= 0.5</td></tr><tr><td>1</td><td>3/4= 0.75</td><td>3/4= 0.75</td></tr><tr><td>2</td><td>-21/8= -2.625</td><td>-15/4= -3.75</td></tr><tr><td>3</td><td>333/16= 20.8125</td><td>54</td></tr><tr><td>4</td><td>-30885/128= -241.289...</td><td>-20817/16= -1301.0625</td></tr><tr><td>5</td><td>916731/256= 3580.98...</td><td>216243/5= 43248.6</td></tr></tbody></table>

### C. Strong coupling expansion

In the strong-coupling regime (large $g$) it has been shown that the expansion in terms of $1/g$ takes the form[43]

$$E_0(g) = g^{\frac{1}{3}} \sum_{n=0}^{+\infty} \alpha_n \left(\frac{1}{g}\right)^{\frac{2n}{3}}. \qquad (69)$$

Accurate values for the coefficients $\alpha_n$ have been reported.[43]

Here, to build the strong-coupling expansion of $E_0(g)$ we first need to evaluate the large-$x$ expansion of $B(x)$. From the definition of $B(x)$, Eq.(29), we easily show that (see, Appendix C)

$$B(x) = \sqrt{x} \sum_{n=0}^{+\infty} B_n \left( \frac{1}{\sqrt{x}} \right)^n . \qquad (70)$$

The three first values of $B_n$ are given in Appendix C. Now, by using the implicit equation for $\bar{\omega}_g$, Eq.(47), we find by simple inspection that

$$\bar{\omega}_g(g) = g^{\frac{1}{3}} \sum_{n=0}^{+\infty} \omega_n \left(\frac{1}{g}\right)^{\frac{2n}{3}}. \qquad (71)$$

The leading coefficient $\omega_0$ is given by

$$\omega_0 = \left(\sqrt{2}B_0\right)^{\frac{2}{3}} g^{\frac{1}{3}} \qquad (72)$$

with

$$B_0 = 2 \frac{\Gamma\left(\frac{5}{4}\right)}{\Gamma\left(\frac{3}{4}\right)}. \qquad (73)$$

The coefficients $\omega_1$ and $\omega_2$ are reported in Appendix C. Now, the strong coupling expansion of the energy may be obtained from the expansion of $\bar{\omega}_g(g)$ and the expression of $E_0(g)$, Eq.(48). We get

$$E_0(g) = g^{\frac{1}{3}} \sum_{n=0}^{+\infty} \alpha_n \left(\frac{1}{g}\right)^{\frac{2n}{3}}. \qquad (74)$$

Remarkably, the functional form of the exact energy, Eq.(69), is recovered by our model. The explicit expression of the first three coefficients $\alpha_n$ are given in Appendix C. Their numerical values are reported in Table

II and compared to the exact ones calculated in [43]. The relative error on the leading coefficient $\alpha_0$ is only about 4%, a remarkable result in view of the simplicity of the model. The accuracy of the next coefficients decreases as a function of $n$ but the values are still reasonable.

<table><caption>Table II. Strong-coupling coefficients $\alpha_n$ compared to the exact values of [44]. Here, $\omega = 1$.</caption><thead><tr><th>n</th><th>α<sub>n</sub></th><th>α<sub>n</sub><sup>ex</sup></th></tr></thead><tbody><tr><td>0</td><td>0.6393...</td><td>0.6679...</td></tr><tr><td>1</td><td>0.1576...</td><td>0.1436...</td></tr><tr><td>2</td><td>-0.0152...</td><td>-0.0086...</td></tr></tbody></table>

## VI. COMPARISON WITH THE FEYNMAN-KLEINERT AND BÜTTNER-FLYTZANIS APPROXIMATE PARTITION FUNCTIONS

Two simple yet accurate partition functions for the anharmonic oscillator proposed in the literature are those of Feynman and Kleinert[28] and of Büttner and Flytzanis[34]. In this section, the respective quality of the three approximate formulas is evaluated.

### 1. Feynman-Kleinert (FK) approach

In the same way as in the present work, the authors make use of the path integral formalism. However, the route followed is very different. In short (for details, see [28]) the starting idea is to rewrite $Z$ as a classical partition function involving an effective classical potential $V_{\text{eff}}(x)$,

$$Z = \frac{1}{\sqrt{2\pi\beta}} \int dx_0 e^{-\beta V_{\text{eff}}(x_0)}. \qquad (75)$$

This representation is exact and particularly well-suited to the high-temperature limit where the formula reduces to the classical partition function, Eq.(44). In this limit the only path contributing to the partition function is the constant path connecting the initial and final points, that is, $x(0) = x(\beta) = x_0$. At small $\beta$, it is useful to introduce the collective (centroïd) variable, $\bar{x} = \int_0^\beta dsx(s)$, since the paths contributing the most to the path integral are those for which $x(s)$ does not fluctuate too much around $\bar{x}$. The deviations from $\bar{x}$ are then locally approximated in a harmonic way, thus generating an approximation for the effective classical potential that can be explicitly obtained by solving a pair of coupled equations. In Appendix D the working formulae giving the effective potential are reported. From the partition function, the approximate ground-state energy can be obtained (see the derivation in Appendix D). We have

$$E_0(g) = \frac{\Omega_0(g)}{4} + \frac{\omega^2}{4\Omega_0(g)} + \frac{3g}{4\Omega_0^2(g)} \qquad (76)$$
with

$$ \Omega_0(g) = \frac{\omega^2}{3^{\frac{1}{3}} \Delta(g)^{\frac{1}{3}}} + \frac{\Delta(g)^{\frac{1}{3}}}{3^{\frac{2}{3}}} \quad (77) $$

and

$$ \Delta(g) = 27g + \sqrt{3\sqrt{243g^2 - \omega^6}}. \quad (78) $$

To obtain the first excited-state energy, $E_1(g)$, reported in the figure 4 below, the subleading component of the FK partition function, $Z_{FK}$, is extracted as follows

$$ E_1(g) = \lim_{\beta \to \infty} -\frac{1}{\beta} \log \left( Z_{FK} - e^{-\beta E_0(g)} \right). \quad (79) $$

## 2. Büttner-Flytzanis (BF) approach

The Büttner-Flytzanis approach is based on a variational approach for the free energy, $F = -\frac{1}{\beta} \ln Z$. As noted by Feynman,[27] an approximation of the free energy is the following upper limit for $F$

$$ \tilde{F} = F_0 + \langle H - H_0 \rangle_0 \quad (80) $$

where $F_0$ is the free energy of some reference Hamiltonian $H_0$ and the average is taken with respect to $H_0$. Using as reference Hamiltonian a harmonic oscillator of frequency $\tilde{\omega}$ with its center displaced by a quantity $b$ we have

$$ \tilde{F} = F_0 + \langle V(x) - \frac{1}{2}\tilde{\omega}^2(x-b)^2 \rangle_0. \quad (81) $$

$\tilde{F}$ can be evaluated analytically and its minimization with respect to the variational parameters $\tilde{\omega}$ and $b$ leads to $b=0$ and to a third-order polynomial equation giving the optimal value for $\tilde{\omega}$. For completeness, the equations are presented in the Appendix E.

## 3. Comparison between models

In Table III a comparison between the free energies obtained with the three models (BF, FK, and our model) is presented. Results are given at some selected values of $\beta$ and coupling constants $g$. For each approach, the relative errors $\epsilon$ (in %) are reported. A first remark is that the BF and FK free energies have positive errors, a result expected due to the variational character of both theories. In constrast, as already noted above, our free energies are found to be systematically smaller than the exact ones, a property which is a priori not expected. Unfortunately, we have not been able to show whether or not this property is true, in particular when other types of potentials $V(x)$ are considered. Now, taking this property for granted, it is tempting to maximize the energy with respect to the effective frequency, $\tilde{\omega}_g$, considered as a variational parameter [note that $E_0(g)$ and $E_1(g)$ depends only on $\tilde{\omega}_g$]. Quite surprisingly,

<table><caption>Table III. Comparison between the free energies obtained with the three models (BF, FK, and the present approach). Relative errors (in %) are reported for different values of $\beta$ in the small- ($g=0.01$), intermediate- ($g=1$), and strong-coupling ($g=10$) regimes. $\omega=1$.</caption><thead><tr><th></th><th>β</th><th>ε(BF)(%)</th><th>ε(FK)(%)</th><th>ε(This work)(%)</th></tr></thead><tbody><tr><td colspan="5">g = 0.01</td></tr><tr><td>10</td><td>0.006</td><td>0.004</td><td>-0.02</td></tr><tr><td>5</td><td>0.007</td><td>0.002</td><td>-0.02</td></tr><tr><td>2</td><td>0.03</td><td>0.0002</td><td>-0.03</td></tr><tr><td>1</td><td>0.9</td><td>0.0002</td><td>-0.3</td></tr><tr><td>0.1</td><td>0.5</td><td>0.01</td><td>-0.002</td></tr><tr><td colspan="5">g = 1</td></tr><tr><td>10</td><td>1.</td><td>1.</td><td>-2.</td></tr><tr><td>5</td><td>1.</td><td>0.8</td><td>-2.</td></tr><tr><td>2</td><td>1.</td><td>0.2</td><td>-2.</td></tr><tr><td>1</td><td>3.</td><td>0.04</td><td>-3.</td></tr><tr><td>0.1</td><td>3.</td><td>0.000001</td><td>-0.08</td></tr><tr><td colspan="5">g = 10</td></tr><tr><td>10</td><td>2.</td><td>2.</td><td>-4.</td></tr><tr><td>5</td><td>2.</td><td>2.</td><td>-4.</td></tr><tr><td>2</td><td>2.</td><td>1.</td><td>-4.</td></tr><tr><td>1</td><td>2.2</td><td>0.4</td><td>-4.</td></tr><tr><td>0.1</td><td>5.</td><td>0.00006</td><td>-0.4</td></tr></tbody></table>

no improvement has been observed at all couplings, a result which would indicate that our parameter-free partition function is already optimal if the property is true. A second remark is that, as already noted by Srivastava and Vishwamittar[45], the quality of the FK free energies is superior to that of the BF approach. This is particularly true in the high-temperature limit where the BF approach fails to converge to the classical limit. The present method, which has both the exact harmonic and classical limits, like in the FK approach, leads also to small errors over the full range of temperatures and couplings (here, maximum error of about 4%). However, the FK approach remains superior (maximum error of about 2%), particularly at high temperatures where the convergence of the FK partition function to the classical limit is particularly rapid.

An important issue regarding approximate partition functions is their ability to reproduce the gap in energy, $\Delta E = E_1(g) - E_0(g)$, a critical quantity for the low-energy properties of the model. In Fig.4 we present the temperature-dependent energy gap defined as follows

$$ \Delta E(\beta) = E_1(\beta) - E_0(\beta) \quad (82) $$

with

$$ E_0(\beta) = -\frac{1}{\beta} \log Z(\beta) \quad (83) $$

and

$$ E_1(\beta) = -\frac{1}{\beta} \log \left[ Z(\beta) - e^{-\beta E_0(\beta \to \infty)} \right]. \quad (84) $$
At $\beta \to \infty$, $\Delta E(\beta)$ is expected to converge to the energy gap. In Fig.4 the gap in energy as a function of $\beta$ is presented for $g=1$ (intermediate coupling). As seen, because of its very construction, the FK approach is not at all suited for reproducing the energy gap. In contrast, this is not the case for the two other methods. In the case of our model, the partition function has been modified by adding to it the small correction
$$\delta Z = e^{-\beta(E_0+\bar{\omega}_g)} [e^{-\beta\Delta} - (1-\beta\Delta)] \text{ in order to stabilize its zero-temperature limit (see, the discussion above about the excited-state energies). Both approaches converge to accurate values with an error of about 3.4\% for the BF approach and a significantly smaller error of about 0.9\% for the present approach.}$$
Finally, it is important to emphasize that, in this section, we have only compared *simple* formulas for the partition function. For the FK and BF methods the simple formulae used here correspond only to the lowest order of approximation in their formalism. In the case of FK, systematic corrections to the calculation of the effective classical potential can be performed,[30] thus leading to systematically better energies and gaps. It is also true in the variational framework for free energies, see for example[46]. However, in both cases, the formalism and formulas needed to evaluate the systematic corrections to the simple formulas become much more involved and the advantage of having a simple model is lost.

Finally, it is important to emphasize that, in this section, we have only compared *simple* formulas for the partition function. For the FK and BF methods the simple formulae used here correspond only to the lowest order of approximation in their formalism. In the case of FK, systematic corrections to the calculation of the effective classical potential can be performed,[30] thus leading to systematically better energies and gaps. It is also true in the variational framework for free energies, see for example[46]. However, in both cases, the formalism and formulas needed to evaluate the systematic corrections to the simple formulas become much more involved and the advantage of having a simple model is lost.

## VII. SUMMARY AND PERSPECTIVES

Let us first briefly summarize the derivation of our model partition function. The first step is to express the partition function as a path integral, Eq.(16). Next, a gaussian approximation for the quartic contribution is introduced, Eq.(26). This approximation amounts to reformulating $Z$ as the partition function of a harmonic oscillator with an effective frequency depending both on the temperature and coupling constant $g$, Eq.(28). However, because of the gaussian approximation the limit $n \to \infty$ for $Z_n$ is no longer defined, Eq.(35). To address this problem, we introduce a Principle of Minimal Sensitivity of the PF to the effective frequency, Eq.(36). Thanks to this principle, we are able to propose our model partition function, Eq.(40) with Eqs.(41,42,43).

Despite its apparent simplicity, this partition function encapsulates a number of important properties which can be summarized as follows.

i) Both the harmonic ($g \to 0$) and classical ($\beta \to 0$) limits are recovered.

ii) The free energy is accurate to a few percent over the entire range of temperatures and coupling strengths.

iii) The ground-state and first-excited energies are accurately reproduced for all coupling strengths, the energy gap having a typical error of about one percent.

iv) The divergence of the weak-coupling expansion with a factorial-like growth of the coefficients is recovered. In addition, the coefficients are found to be rational numbers as the exact ones.

<figure><img src="image_6.png" /><figcaption>Figure 4. Comparison between the temperature-dependent energy gaps as a function of $\beta$. The energy gap, $\Delta E = E_1(g) - E_0(g)$, is obtained as $\beta \to \infty$. The solid line gives the exact values. $\omega = 1$ and $g = 1$.</figcaption></figure>

v.) The functional form of the strong-coupling expansion is reproduced. The leading coefficient is obtained with an error of a few percent.

The quartic oscillator being one of the simplest paradigmatic examples of non-trivial quantum system, it has been and still will be the subject of numerous studies in quantum, statistical, and quantum field theory domains to test new ideas and algorithms. Then, we believe that to have at our disposal a simple analytic model for the partition function reproducing the key features of the exact partition function can be of interest to the community. Note that the approach presented here can be generalized without difficulty to one-dimensional oscillators with an arbitrary anharmonicity (for example, $x^p$ with $p$ even). Finally, it will be interesting to explore the applicability of the approach to multi-dimensional systems.
# ACKNOWLEDGMENTS

I would like to thank the Centre National de la Recherche Scientifique (CNRS) for its continued support. I also acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 863481).

Appendix A: Derivation of $\prod_{k=1}^{n-1} \sin \frac{\pi k}{n} = \frac{n}{2^{n-1}}$

The derivation of this equality can be found at several places in the literature. However, for completeness, we provide the following derivation.

For $1 \le k \le n - 1$ we have

$$|1 - e^{i2\pi \frac{k}{n}}| = 2 \sin \frac{\pi k}{n}. \quad (A1)$$

The product of sine can thus be written as

$$\prod_{k=1}^{n-1} \sin \frac{\pi k}{n} = \frac{1}{2^{n-1}} \prod_{k=1}^{n-1} |1 - e^{2\pi i \frac{k}{n}}| = \frac{1}{2^{n-1}} \prod_{k=1}^{n-1} (1 - e^{2\pi i \frac{k}{n}}). \quad (A2)$$

Let us introduce the roots of the unity, $w_k = e^{i \frac{2\pi k}{n}}$, solutions of $x^n - 1 = 0$. We have

$$x^n - 1 = \prod_{k=0}^{n-1} (x - w_k). \quad (A3)$$

Defining $Q(x)$ as

$$Q(x) = \prod_{k=1}^{n-1} (x - w_k) \quad (A4)$$

and writing

$$x^n - 1 = (x - 1)(1 + x + \dots + x^{n-1}) \quad (A5)$$

we then get

$$Q(x) = 1 + x + \dots + x^{n-1} \quad (A6)$$

and

$$|\prod_{k=1}^{n-1} (1 - e^{2\pi i \frac{k}{n}})| = |Q(1)| = n. \quad (A7)$$

Finally,

$$\prod_{k=1}^{n-1} \sin \frac{\pi k}{n} = \frac{n}{2^{n-1}}. \quad (A8)$$

# Appendix B: Derivation of $\lim_{n \to \infty} P_n(x)$

$P_n(x)$ is defined as

$$P_n(x) = \prod_{k=2}^{n} \frac{1}{\sqrt{1 + \frac{x^2}{2\lambda_k n^2}}} \quad (B1)$$

where

$$\lambda_k = 1 - \cos \frac{2\pi}{n}(k-1) \quad k=1, n. \quad (B2)$$

Introducing the logarithm of $P_n$

$$\ln P_n = -\frac{1}{2} \sum_{k=2}^{n} \ln \left( 1 + \frac{x^2}{2\lambda_k n^2} \right) \quad (B3)$$

and the Taylor expansion of $\ln(1+x)$

$$\ln (1 + x) = \sum_{l=1}^{+\infty} \frac{(-1)^{l-1}}{l} x^l \quad (B4)$$

we get

$$\ln P_n = \frac{1}{2} \sum_{l=1}^{+\infty} \frac{(-1)^l}{l} \left( \frac{x^2}{2n^2} \right)^l S_n^l \quad (B5)$$

where

$$S_n^l = \sum_{k=2}^{n} \left( \frac{1}{\lambda_k} \right)^l = \sum_{k=2}^{n} \left( \frac{1}{1 - \cos \left( 2\pi \frac{k-1}{n} \right)} \right)^l. \quad (B6)$$

$S_n^l$ can be rewritten as

$$S_n^l = \frac{1}{2l} \sum_{k=1}^{n-1} \frac{1}{\sin^{2l} \left( \pi \frac{k}{n} \right)}. \quad (B7)$$

As shown by Fisher[47], the large-$n$ behavior of $S_n^l$ can be related to the Riemann zeta function as follows

$$\lim_{n \to \infty} n^{-2l} \sum_{k=1}^{n-1} \frac{1}{\sin^{2l} \left( \pi \frac{k}{n} \right)} = 2 \frac{\zeta(2l)}{\pi^{2l}}. \quad (B8)$$

Now, using the well-known equality

$$\zeta(2l) = (-1)^{l+1} \frac{B_{2l}(2\pi)^{2l}}{2(2l)!} \quad (B9)$$

we get the asymptotic behavior at large $n$ of $S_n^l$

$$S_n^l \sim_{n \to \infty} \frac{|B_{2l}| 2^l}{(2l)!} n^{2l} \quad (B10)$$

where $B_{2l}$ are the Bernoulli numbers. From now on, all following expressions must be understood as valid only in the large-$n$ regime. Using Eqs. (B5) and (B10) we get

$$\ln P_n = \frac{1}{2} \sum_{l=1}^{+\infty} (-1)^l x^{2l} \frac{|B_{2l}|}{l(2l)!}. \quad (B11)$$
Let us define

$$u(x) = \frac{1}{2} \sum_{l=1}^{+\infty} (-1)^l x^{2l} \frac{|B_{2l}|}{l(2l)!} \quad (B12)$$

so that

$$P_n = e^{u(x)}. \quad (B13)$$

For $l \ge 1$ we have

$$|B_{2l}| = (-1)^{l-1} B_{2l}. \quad (B14)$$

Then, $u(x)$ can be written as

$$u(x) = -\frac{1}{2} \sum_{l=1}^{+\infty} x^{2l} \frac{B_{2l}}{l(2l)!}. \quad (B15)$$

Now, taking the derivative of $u(x)$

$$u'(x) = -\frac{1}{x} \sum_{l=1}^{+\infty} x^{2l} \frac{B_{2l}}{(2l)!} \quad (B16)$$

and using the fact that

$$\coth x = \sum_{n=0}^{+\infty} \frac{2^{2n} B_{2n}}{(2n)!} x^{2n-1} = \frac{1}{x} \sum_{n=0}^{+\infty} \frac{B_{2n}}{(2n)!} (2x)^{2n} \quad (B17)$$

the derivative of $u(x)$ is found to be

$$u'(x) = -\frac{1}{2} \coth \frac{x}{2} + \frac{1}{x}. \quad (B18)$$

Using Eq.(B13) and $P_n(0) = 1$, we can write

$$P_n(x) = e^{\int_0^x dy yu'(y)}. \quad (B19)$$

Finally, after integrating we get

$$P_n(x) = \frac{x}{e^{-\frac{x}{2}} - e^{-\frac{x}{2}}}. \quad (B20)$$

## Appendix C: Strong coupling

The function $B(x)$ is defined as

$$B(x) = \frac{1}{2} \frac{\int dy e^{-y^2 - xy^4}}{\int dy y^2 e^{-y^2 - xy^4}} \quad \text{for } x \ge 0. \quad (C1)$$

After a change of variable, it can be rewritten as

$$B(x) = \frac{\sqrt{x}}{2} \frac{\int dz e^{-\frac{1}{\sqrt{x}}z^2 - z^4}}{\int dz z^2 e^{-\frac{1}{\sqrt{x}}z^2 - z^4}}. \quad (C2)$$

The strong coupling expansion of $B(x)$ at large $x$ has thus the following form

$$B(x) = \sqrt{x} \sum_{n=0}^{+\infty} B_n \left( \frac{1}{\sqrt{x}} \right)^n. \quad (C3)$$

The three first coefficients are given by

$$B_0 = 2 \frac{\Gamma(\frac{5}{4})}{\Gamma(\frac{3}{4})} \quad (C4)$$

$$B_1 = \frac{4\Gamma^2(\frac{5}{4}) - \Gamma^2(\frac{3}{4})}{2\Gamma^2(\frac{3}{4})} \quad (C5)$$

and

$$B_2 = - \frac{\Gamma(\frac{5}{4})[\Gamma^2(\frac{3}{4}) - 8\Gamma^2(\frac{5}{4}) + 4\Gamma(\frac{3}{4})\Gamma(\frac{7}{4})]}{4\Gamma^3(\frac{3}{4})}. \quad (C6)$$

As shown in the text, the expansion of $\bar{\omega}_g(g)$ writes

$$\bar{\omega}_g(g) = (2g)^{\frac{1}{3}} \sum_{n=0}^{+\infty} \omega_n y^n \quad (C7)$$

where the variable $y$ is defined as

$$y = g^{-\frac{2}{3}}. \quad (C8)$$

The coefficients $\omega_n$ are expressed in terms of the coefficients $B_p$'s ($p \le n$) as follows

$$\omega_0 = B_0^{\frac{2}{3}} \quad (C9)$$

$$\omega_1 = \frac{2^{\frac{2}{3}} \omega B_1}{3 B_0^{\frac{2}{3}}} \quad (C10)$$

$$\omega_2 = \frac{\omega^2}{2^{\frac{2}{3}}} \left( \frac{2B_0 B_2 - (B_1)^2}{3(B_0)^2} \right) \quad \text{etc.} \quad (C11)$$

The expression of the ground-state energy as a function of $\bar{\omega}_g$ writes

$$E_0(g) = \frac{\bar{\omega}_g}{2} \left[ 1 - \ln \left( \sqrt{\frac{\bar{\omega}_g}{\pi}} \bar{I}(g) \right) \right] \quad (C12)$$

where

$$\bar{I}(g) = \int dx e^{-\frac{2}{\bar{\omega}_g} \left( \frac{1}{2} \omega^2 x^2 + g x^4 \right)}. \quad (C13)$$

Defining $J(g)$ as

$$J(g) = \sqrt{\frac{\bar{\omega}_g}{\pi}} \bar{I}(g) = \sum_{n=0}^{+\infty} J_n y^n \quad (C14)$$

the three first coefficients $J_n$ are found to be

$$J_0 = \frac{\omega_0^{\frac{3}{4}} \Gamma(\frac{1}{4})}{2\sqrt{\pi}} \quad (C15)$$

$$J_1 = \frac{2^{\frac{1}{3}} \omega^2 \sqrt{\omega_0 \Gamma(-\frac{1}{4}) + 6 \omega_1 \Gamma(\frac{1}{4})}}{16 \sqrt{\pi} \omega_0^{\frac{1}{4}}} \quad (C16)$$
and

$$J_2 = \frac{2^{\frac{1}{3}} \sqrt{\omega_0 \omega_1 \omega^2 \Gamma(-\frac{1}{4}) + (2^{\frac{2}{3}} \omega^4 \omega_0 - 3 \omega_1^2 + 24 \omega_0 \omega_2) \Gamma(\frac{1}{4})}}{64 \sqrt{\pi} \omega_0^{\frac{5}{4}}}. \quad (C17)$$

Writing $E_0(g)$ under the form

$$E_0(g) = \frac{(2g)^{\frac{1}{3}}}{2} \left( \sum_{n=0}^{+\infty} \omega_n y^n \right) \left[ 1 - \ln \left( \sum_{n=0}^{+\infty} J_n y^n \right) \right] \quad (C18)$$

we have

$$E_0(y) = g^{\frac{1}{3}} \sum_{n=0}^{+\infty} \alpha_n y^n. \quad (C19)$$

The first coefficients are given by

$$\alpha_0 = \omega_0(1 - \ln J_0) \quad (C20)$$

$$\alpha_1 = -\frac{J_1 \omega_0}{J_0} + \omega_1(1 - \ln J_0) \quad (C21)$$

and

$$\alpha_2 = \frac{(J_1^2 - 2J_0J_2)\omega_0}{2J_0^2} - \frac{J_1\omega_1}{J_0} + \omega_2(1 - \ln J_0). \quad (C22)$$

## Appendix D: Feynman-Kleinert variational approach

Using the notations of [28] (in particular, $t \equiv \beta$), the basic equations of the Feynman-Kleinert variational approach are

$$Z = \frac{1}{\sqrt{2\pi t}} \int dx_0 e^{-tW_1(x_0)} \quad (D1)$$

where $W_1$ is an effective classical potential obtained as

$$W_1(x_0) = \min_{a^2, \Omega} \tilde{W}_1(x_0, a^2, \Omega) \quad (D2)$$

and $\tilde{W}_1$ is given by

$$\tilde{W}_1(x_0, a^2, \Omega) = \frac{1}{t} \ln \frac{\sinh \frac{\Omega t}{2}}{\frac{\Omega t}{2}} - \frac{\Omega^2}{2} a^2 + V_{a^2}(x_0). \quad (D3)$$

Here, the potential $V_{a^2}$ is defined by

$$V_{a^2}(x_0) = \int \frac{dx'}{\sqrt{2\pi a^2}} e^{-\frac{1}{2a^2}(x_0 - x')^2} V(x'). \quad (D4)$$

After minimization of $\tilde{W}_1$ a pair of optimal parameters depending on $x_0$, $(a^2(x_0), \Omega(x_0))$ is obtained. The minimization of $\tilde{W}_1$ with respect to $\Omega$ gives the equation

$$a^2 = \frac{1}{\Omega^2 t} \left( \frac{\Omega t}{2} \coth \frac{\Omega t}{2} - 1 \right) \quad (D5)$$

and the minimization with respect to $a^2$ gives

$$\Omega^2 = \frac{\partial^2 V_{a^2}(x_0)}{\partial x_0^2}. \quad (D6)$$

Let us explicit these equations in the case of the quartic oscillator, $V(x) = \frac{1}{2}\omega^2 x^2 + gx^4$. We get

$$V_{a^2}(x) = \frac{1}{2} a^2 \omega^2 + 3a^4 g + \frac{1}{2} (12a^2 g + \omega^2) x^2 + gx^4 \quad (D7)$$

and

$$\Omega^2(x) = 12a^2g + \omega^2 + 12gx^2. \quad (D8)$$

Now, to get the ground-state energy, large-time limit has to be considered. As $t \to \infty$ we have

$$\frac{1}{t} \ln \frac{\sinh \frac{\Omega(x_0)t}{2}}{\frac{\Omega(x_0)t}{2}} \to \frac{\Omega(x_0)}{2} \quad (D9)$$

and then

$$a^2(x_0) = \frac{1}{2\Omega(x_0)}. \quad (D10)$$

The equation determining $\Omega(x_0)$ is

$$\Omega(x_0)^3 - \Omega(x_0)(\omega^2 + 12gx_0^2) - 6g = 0 \quad (D11)$$

and the partition function writes

$$Z = \frac{1}{\sqrt{2\pi t}} \int dx_0 e^{-t f(x_0)} \quad (D12)$$

with

$$f(x_0) = \frac{\Omega(x_0)}{4} + V_{a^2}(x_0). \quad (D13)$$

Now, in the large-time limit only the constant part of $f(x)$ will give a contribution to $-\frac{1}{t} \ln Z$, we then have

$$E_0(g) = f(0). \quad (D14)$$

The real solution of the third-order polynomial equation, Eq. (D11), with $x_0 = 0$ is

$$\Omega_0(g) = \frac{\omega^2}{3^{\frac{1}{3}} \Delta(g)^{\frac{1}{3}}} + \frac{\Delta(g)^{\frac{1}{3}}}{3^{\frac{2}{3}}} \quad (D15)$$

where

$$\Delta(g) = 27g + \sqrt{3\sqrt{243g^2 - \omega^6}}. \quad (D16)$$

Finally, the ground-state energy of the Feynman-Kleinert model writes

$$E_0(g) = \frac{\Omega_0(g)}{4} + \frac{\omega^2}{4\Omega_0(g)} + \frac{3g}{4\Omega_0^2(g)}. \quad (D17)$$
Appendix E: Büttner-Flytzanis variational approach

In the Büttner-Flytzanis variational approach,[34] the free energy is written as

$$F(\beta) = F_0(\beta) + \frac{\omega^2 - \tilde{\omega}^2}{4} \left[ \frac{\coth \frac{\tilde{\omega}\beta}{2}}{\tilde{\omega}} \right] + \frac{3g}{4} \left[ \frac{\coth \frac{\tilde{\omega}\beta}{2}}{\tilde{\omega}} \right]^2 \quad (E1)$$

where the frequency $\tilde{\omega}$ is the variational parameter and $F_0(\beta)$ the free energy of the harmonic oscillator with fre-

quency $\tilde{\omega}$

$$F_0(\beta) = \beta^{-1} \ln \left[ 2 \sinh \left( \frac{\tilde{\omega}\beta}{2} \right) \right]. \quad (E2)$$

Minimizing Eq. (E1) with respect to $\tilde{\omega}$ gives the equation determining $\tilde{\omega}$

$$\tilde{\omega}^3 - \omega^2 \tilde{\omega} - 6g \coth \left( \frac{\tilde{\omega}\beta}{2} \right) = 0. \quad (E3)$$

[1] Carl M. Bender and Tai Tsun Wu. Anharmonic oscillator. Phys. Rev., 184:1231–1260, Aug 1969.

[2] Carl M. Bender and Tai Tsun Wu. Large-order behavior of perturbation theory. Phys. Rev. Lett., 27:461–465, Aug 1971.

[3] Jr. G. A. Baker. Essentials of Padé Approximants. Academic, New-York, 1975.

[4] Barry Simon and A Dicke. Coupling constant acidity for the anharmonic oscillator. Annals of Physics, 58(1):76–136, 1970.

[5] Barry Simon. Large orders and summability of eigenvalue perturbation theory: A mathematical overview. International Journal of Quantum Chemistry, 21(1):3–25, 1982.

[6] E. Borel. Leçons sur les Séries Divergentes. Gautier-Villars, Paris, 1928.

[7] Héctor Mera, Thomas G. Pedersen, and Branislav K. Nikolić. Fast summation of divergent series and resurgent transseries from meijer-g approximants. Phys. Rev. D, 97:105027, May 2018.

[8] Ernst Joachim Weniger. Nonlinear sequence transformations for the acceleration of convergence and the summation of divergent series. Computer Physics Reports, 10(5):189–371, 1989.

[9] Ernst Joachim Weniger, Jiří Čížek, and Francine Vinette. Very accurate summation for the infinite coupling limit of the perturbation series expansions of anharmonic oscillators. Physics Letters A, 156(3):169–174, 1991.

[10] J. C. Le Guillou and J. Zinn-Justin. Critical exponents from field theory. Phys. Rev. B, 21:3976–3998, May 1980.

[11] J. Zinn-Justin R. Seznec. Summation of divergent series by order dependent mappings: Application to the anharmonic oscillator and critical exponents in field theory. J. Math. Phys., pages 1398–1408, 1979.

[12] I. A. Ivanov. Reconstruction of the exact ground-state energy of the quartic anharmonic oscillator from the coefficients of its divergent perturbation expansion. Phys. Rev. A, 54:81–86, Jul 1996.

[13] Martin Schwarz. Statistical thermodynamics of an anharmonic oscillator. Journal of Statistical Physics, 15(9):255–261, 1976.

[14] M.M. Pant and S.K. Mitra. Thermodynamics of anharmonic oscillator. Chemical Physics Letters, 62(3):495–498, 1979.

[15] F. T. Hioe and E. W. Montroll. Quantum theory of anharmonic oscillators. I. Energy levels of oscillators with positive quartic anharmonicity. Journal of Mathematical Physics, 16(9):1945–1955, 09 1975.

[16] E. Wigner. On the quantum correction for thermodynamic equilibrium. Phys. Rev., 40:749–759, Jun 1932.

[17] John G. Kirkwood. Quantum statistics of almost classical assemblies. Phys. Rev., 44:31–37, Jul 1933.

[18] L.D. Landau and EM. Lifshitz. Statistical physics. Pergamon Press, London, 1959.

[19] W. Witschel. The semiclassical and the quasiclassical partition function of a quartic anharmonic oscillator. Chemical Physics, 50(2):265–271, 1980.

[20] Petr Jizba and Václav Zatloukal. Path-integral approach to the wigner-kirkwood expansion. Phys. Rev. E, 89:012135, Jan 2014.

[21] W. Witschel. An improved variation-perturbation technique for the partition function of anharmonic oscillators. Physics Letters A, 97(8):315–319, 1983.

[22] W. Witschel, K. Murmann, and G. Taubmann. A critical test of the diagonal variation-perturbation technique for the partition function of anharmonic oscillators. The Journal of Chemical Physics, 97(11):8637–8643, 12 1992.

[23] A. W. Sáenz and R. C. O'Rourke. Number of states and the magnetic properties of an electron gas. Rev. Mod. Phys., 27:381–398, Oct 1955.

[24] F. Vinette and J. Čížek. Upper and lower bounds of the ground state energy of anharmonic oscillators using renormalized inner projection. Journal of Mathematical Physics, 32(12):3392–3404, 12 1991.

[25] M. D. Girardeau and R. M. Mazo. Variational methods in statistical mechanics. in: Advances in Chemical Physics, volume 24. I. Prigogine and S. A. Rice Academic Press, New York, 1974.

[26] L. de Falco, R. Mignani, and R. Scipioni. Approximation method for the partition function of an anharmonic-oscillator gas. Europhysics Letters, 29(9):659, mar 1995.

[27] R.P. Feynman. Statistical mechanics. Benjamin, Reading, 1972.

[28] R. P. Feynman and H. Kleinert. Effective classical partition functions. Phys. Rev. A, 34:5080–5084, Dec 1986.

[29] Riccardo Giachetti and Valerio Tognetti. Variational approach to quantum statistical mechanics of nonlinear systems with application to sine-gordon chains. Phys. Rev. Lett., 55:912–915, Aug 1985.

[30] H. Kleinert. Systematic corrections to the variational calculation of the effective classical potential. Physics Letters A, 173(4):332–342, 1993.

[31] H. Kleinert and H. Meyer. Variational calculation of effective classical potential at non-zero temperature to higher orders. Physics Letters A, 184(4):319–327, 1994.
[32] R. Karrlein and H. Kleinert. Precise variational tunneling rates of anharmonic oscillator for g negative. *Physics Letters A*, 187(2):133-139, 1994.

[33] Hagen Kleinert. *Path integrals in quantum mechanics, statistics, polymer physics, and financial markets; 4th ed.* World Scientific, River Edge, NJ, 2006. Based on a Course on Path Integrals, Freie Univ. Berlin, 1989/1990.

[34] H. Büttner and N. Flytzanis. Effective free energies. *Phys. Rev. A*, 36:3443-3445, Oct 1987.

[35] L.S. Schulman. *Techniques And Applications of Path Integration*. Dover Publications Inc., 2005.

[36] R.P. Feynman and A.R. Hibbs. *Quantum Mechanics and Paths Integrals*. McGraw-Hill, New York, 1965.

[37] John T. Marshall and John L. Pell. Path-integral evaluation of the space-time propagator for quadratic Hamiltonian systems. *Journal of Mathematical Physics*, 20(7):1297-1302, 07 2008.

[38] S. M. Cohen. Path integral for the quantum harmonic oscillator using elementary methods. *American Journal of Physics*, 66(6):537-540, 06 1998.

[39] Barry R. Holstein. The harmonic oscillator propagator. *American Journal of Physics*, 66(7):583-589, 07 1998.

[40] F. A. Barone, H. Boschi-Filho, and C. Farina. Three methods for calculating the Feynman propagator. *American Journal of Physics*, 71(5):483-491, 05 2003.

[41] Kiyoto Hira. Derivation of the harmonic oscillator propagator using the feynman path integral and recursive relations. *European Journal of Physics*, 34(3):777, apr 2013.

[42] Wolfram Research, Inc. Mathematica, Version 13.3. Champaign, IL, 2023.

[43] W. Janke and H. Kleinert. Convergent strong-coupling expansions from divergent weak-coupling perturbation theory. *Phys. Rev. Lett.*, 75:2787-2791, Oct 1995.

[44] W. Janke and H. Kleinert. Convergent strong-coupling expansions from divergent weak-coupling perturbation theory. *Phys. Rev. Lett.*, 75:2787-2791, Oct 1995.

[45] Sunita Srivastava and Vishwamittar. Determination of free energies of an oscillator with mixed quartic and sexual anharmonicities. *Phys. Rev. A*, 44:8006-8011, Dec 1991.

[46] K. Vlachos. Variational method for the free-energy approximation of generalized anharmonic oscillators. *Phys. Rev. A*, 47:838-846, Feb 1993.

[47] L.A. Gardner and M.L. Fisher. Sum of inverse powers of cosines. *SIAM Review*, 13:116, 1971.