# Analytic Model for the Energy Spectrum of the Anharmonic Oscillator

Michel Caffarel

► To cite this version:

Michel Caffarel. Analytic Model for the Energy Spectrum of the Anharmonic Oscillator. Physics Letters A, 2024, 525 (129925), 10.1016/j.physleta.2024.129925. hal-04705550

HAL Id: hal-04705550

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

Submitted on 23 Sep 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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Analytic Model for the Energy Spectrum of the Anharmonic Oscillator

Michel Caffarel¹, *

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

In a recent work we have proposed an original analytic expression for the partition function of the quartic oscillator. This partition function, which has a simple and compact form with *no adjustable parameters*, reproduces some key mathematical properties of the exact partition function and provides free energies accurate to a few percent over a wide range of temperatures and coupling constants. In this work, we present the derivation of the energy spectrum of this model. We also generalize our previous study limited to the quartic oscillator to the case of a general anharmonic oscillator. Numerical application for a potential of the form $V(x) = \frac{\omega^2}{2}x^2 + gx^{2m}$ show that the energy levels are obtained with a relative error of about a few percent, a precision which we consider to be quite satisfactory given the simplicity of the model, the absence of adjustable parameters, and the negligible computational cost.

The one-dimensional quantum anharmonic oscillator[1] plays an important role in quantum mechanics as a simple yet nontrivial model for describing nonlinear (anharmonic) effects. As such, it is commonly used in various scientific fields including molecular physics (rovibrational spectra), condensed matter physics (phonons), nuclear physics (collective vibrational motions of nuclei), and quantum field theory (e.g., $\phi^4$-theory, Higgs mechanism), to cite the main ones. No exact solution for a general anharmonic potential has been found so far, although some partial solutions have been developed for specific potentials, for example in the case of the so-called quasi-exactly-solvable problems; see Turbiner[2] and references therein. From a numerical point of view, the simplicity of the ordinary differential equation to solve allows very accurate solutions. To cite one numerical approach among many, let us mention the recently developed Lagrange Mesh method of Baye [3] as recently implemented by del Valle[4]. However, as for any physical model, having analytical solutions -even approximate ones- is important since it may lead a deeper insight into the nature of the problem, reveal underlying hidden structures that may not be evident from numerical solutions and, also, guide further analysis or generalizations. A great variety of approximate analytical approaches have been proposed, making it very difficult to provide an exhaustive account of the literature. In this work, we are more specifically interested in evaluating the energy levels of the anharmonic oscillator. Among the main approaches developed for this purpose, let us mention the semi-classical approaches, such as WKB[5, 6] or the phase-integral method based on the generalized Bohr-Sommerfeld quantization condition ([7] and references therein) and the methods based on the use of a variational approach either using parameterized excited wavefunctions[8] or path integrals[9, 10]. Other interesting approaches include the design of simple analytic expressions for the energy levels built from the weak- and/or strong-coupling expansions[11, 12] and

methods based on the tuning of boundary conditions[13] or eigenvector continuation [14].

The aim of this work is to present novel analytical expressions for the energy levels of the anharmonic oscillator. For that, we take advantage of a recently proposed partition function for the quartic oscillator,[15] which we extend here to an arbitrary potential. As shown in our previous study[15], this partition function provides an appealing *simple* and *physically meaningful* model for the exact solution. Indeed, it has a simple and compact form with *no adjustable parameters*, a desirable property for an analytical model. Furthermore, it reproduces some key mathematical properties of the exact partition function, thus supporting the idea that the important features of the non-trivial mathematical structure of the solution are, at least partially, accounted for.

These properties are the following: i) The harmonic and classical limits are exactly recovered, ii) the well-known divergence of the weak-coupling (Rayleigh-Schrödinger) expansion of the energy is reproduced. As for the exact solution, the energy corrections are found to be rational numbers and display a factorial-like growth in terms of the perturbational order, and iii) the functional form of the strong-coupling expansion is recovered.

From a quantitative point of view, the free energy is found to be accurate to a few percent over a wide range of temperatures and coupling constants. A similar precision is also obtained for the ground- and first-excited state energies.

To the best of our knowledge, no partition function in closed-form proposed so far is capable of simultaneously reproducing all these features (see, for example, references [9, 16–18] and a comparative study with other models in [15]).

Having at our disposal a simple and physically meaningful model for the partition function, it is of interest to derive the full energy spectrum of the model. This is the purpose of the present work. As we shall see, the energy levels obtained turn out to be accurate with a relative error of about a few percent, a precision which we consider to be quite satisfactory for such a simple model. Of course, a higher precision can be obtained without

<footnote>* caffarel@irsamc.ups-tlse.fr</footnote>
difficulty by using some of the numerical/analytical approaches cited above. However, we emphasize that the aim of this work is not to achieve an ultimate precision in the spectrum but, instead, to propose a simple yet faithful analytical model for the spectrum of a general anharmonic oscillator that can be readily applied in various scientific contexts with essentially no computational cost.

The Hamiltonian considered here is as follows

$$H = -\frac{1}{2} \frac{d^2}{dx^2} + V(x) \quad (1)$$

where $V(x)$ is a rather general potential function bounded from below and verifying $\lim_{|x| \to \infty} V(x) = +\infty$ (in other words, there are only discrete energies). In the numerical applications presented below, we will particu- larize the potential in the form

$$V(x) = \frac{\omega^2}{2} x^2 + g x^{2m} \quad (2)$$

where $\omega^2$ and $g$ denote the harmonic force constant and coupling constant, respectively.

As a first step let us briefly recap the main steps leading to the model partition function.

i) The partition function is first expressed as a path integral in a standard way (see, e.g., [19])

$$Z = \text{Tr} e^{-\beta H} = \lim_{n \to \infty} Z_n \quad (3)$$

with

$$Z_n = \left( \frac{1}{\sqrt{2\pi\tau}} \right)^n \int_{-\infty}^{\infty} dx_1 \dots \int_{-\infty}^{\infty} dx_n \\ e^{-\frac{1}{2\tau} \sum_{i=1}^{n} (x_{i+1} - x_i)^2} e^{-\tau \sum_{i=1}^{n} V(x_i)}. \quad (4)$$

where $\tau = \frac{\beta}{n}$ is the time-step and periodic conditions are used, $x_{n+1} = x_1$.

ii) Second, the short-time anharmonic contribution, $e^{-\tau V(x)}$, is approximated by a gaussian distribution cen- tered on some position $x^*$ (typically, the position of the lowest minimum of the potential; here, $x^* = 0$ in our applications) with an effective frequency $\omega_g(\tau)$, that is

$$\frac{e^{-\tau V(x)}}{\int_{-\infty}^{\infty} dx e^{-\tau V(x)}} \sim \frac{e^{-\tau \frac{1}{2} \omega_g^2(\tau) (x-x^*)^2}}{\int_{-\infty}^{\infty} dx e^{-\tau \frac{1}{2} \omega_g^2(\tau) (x-x^*)^2}}. \quad (5)$$

To set the frequency $\omega_g(\tau)$ we have proposed to impose to the two distributions to have the same variance. In the case of the quartic oscillator treated previously,[15] it leads to $\omega_g = \omega \sqrt{B[\frac{4g}{\tau\omega^4}]}$ where B is some function. In the more general case $V(x) = \frac{\omega^2}{2} x^2 + g x^{2m}$, we get

$$\omega_g(\tau) = \omega \sqrt{B^{(m)} \left[ \frac{2mg}{\tau^{m-1}\omega^{2m}} \right]} \quad (6)$$

where the parameter-free function $B^{(m)}(x)$ is given by

$$B^{(m)}(x) = \frac{1}{2} \frac{\int_{-\infty}^{\infty} dy e^{-y^2 - xy^{2m}}}{\int_{-\infty}^{\infty} dy y^2 e^{-y^2 - xy^{2m}}} \quad (7)$$

For a general potential V the formula writes

$$\omega_g(\tau) = \sqrt{\frac{1}{\tau} \frac{\int_{-\infty}^{\infty} dx e^{-\tau V(x)}}{\int_{-\infty}^{\infty} dx (x-x^*)^2 e^{-\tau V(x)}}} \quad (8)$$

iii) The gaussian approximation being made, the infinite-$n$ limit of $Z_n$ is no longer defined, $\lim_{n \to \infty} Z_n = +\infty$. To circumvent this problem, we have proposed to introduce a Principle of Minimal Sensitivity (PMS) for the path integral. More precisely, we impose to the path integral to minimally depend on the effective frequency used in the gaussian approximation, that is, $\frac{\partial Z_n}{\partial \omega_g(\tau)} = 0$. For a given temperature, the PMS condition holds only for a unique value of $n$ denoted as $n_c(\beta)$ [and, thus, a unique time-step, $\tau_c(\beta) = \frac{\beta}{n_c(\beta)}$]. We have then proposed to define the model partition function as the value of $Z_n$ at this "optimal" value of $n$, $Z = Z_{n_c(\beta)}$. In particular, the ill-defined $n \to \infty$-limit is avoided. The nonlinear implicit equation for $n_c(\beta)$ resulting from the PMS condition writes

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

iv) Finally, the analytical expression of the model partition function is given by

$$Z = \frac{C(\beta)^{n_c(\beta)}}{e^{\frac{\beta \omega_g[\tau_c(\beta)]}{2}} - e^{-\frac{\beta \omega_g[\tau_c(\beta)]}{2}}} \quad (10)$$

with

$$C(\beta) = \sqrt{\frac{\omega_g[\tau_c(\beta)]}{\pi \coth \frac{\beta \omega_g[\tau_c(\beta)]}{2}}} I_g(\beta) \quad (11)$$

and

$$I_g(\beta) = \int_{-\infty}^{\infty} dx e^{-\tau_c(\beta)V(x)}, \quad (12)$$

where the time-step $\tau_c(\beta) = \frac{\beta}{n_c(\beta)}$ is obtained from Eq.(9).

At this point, an important remark is in order. Introducing a gaussian approximation of a non-gaussian quantity is a standard practice in physics (harmonic phonons, Gaussian approximation for the Ginzburg-Landau action, etc.) In short, it is done by restricting to the second-order a Taylor expansion of some Hamiltonian or action. Here, our gaussian approximation is very different in nature. Instead of approximating the potential $V(x)$ by a quadratic potential, we approximate the quantity $e^{-\tau V(x)}$ by a
gaussian distribution with an effective frequency which depends explicitly on the time-step. It is this dependence on the time-step that makes $Z_n$ to diverge in the large-$n$ limit, but which, after application of the PMS condition, allows to fix $n$ at a finite value $n_c(\beta)$ and, then, leads to an accurate model for the exact partition function. Approximating $V(x)$ using a quadratic potential independent on the time-step would merely lead to the partition function of a simple harmonic oscillator and, thus, to a very poor model for the exact partition function.

A. Ground-state energy. The ground-state energy is obtained from the large-$\beta$ behavior of $Z$

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

At large $\beta$, the (unique) solution $n_c(\beta)$ of Eqs.(6) and (9) is proportional to $\beta$ and given by

$$n_c(\beta) = \frac{\beta \bar{\omega}_g}{2}, \qquad (14)$$

where $\bar{\omega}_g = \omega_g[\tau_c(+\infty)]$ is the solution of the implicit equation given by

$$\bar{\omega}_g = \omega \sqrt{B^{(m)} \left( \frac{2g\bar{\omega}_g^{m-1}}{\omega^{2m}} \right)}. \qquad (15)$$

Using the zero-temperature limit of the free energy and the expression 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] \qquad (16)$$

where

$$\bar{I}_g = \int_{-\infty}^{\infty} dx e^{-\frac{2}{\bar{\omega}_g} V(x)}. \qquad (17)$$

B. Excited-state energies. The exact partition function of the quartic oscillator decomposes as a discrete sum of exponentials

$$Z = \sum_{n=0}^{\infty} e^{-\beta E_n} \qquad (18)$$

where $E_n$ are the excited-state energies. As just seen, the ground-state energy $E_0$ is obtained by extracting the leading exponential component of the PF at large $\beta$. To get excited-state energies subleading components are to be evaluated.

In the following we will show that the model partition function actually does not write as a sum of simple exponentials, as it should be for the exact PF, but, instead, as a sum of exponentials multiplied by a polynomial term as follows

$$Z = \sum_{n=0}^{\infty} P_n(\beta) e^{-\beta(E_0 + n\bar{\omega}_g)} \qquad (19)$$

where $P_n(\beta)$ is a polynomial of degree $n$ in $\beta$ with $P_n(0) = 1$. The partition function being no longer expressed as a sum of simple exponentials, the definition of what is meant by excited-state energies becomes problematic. Comparisons with the "exact" numerical spectrum show that defining the excited energies as the exponents of the exponential contributions, that is,

$$E_n = E_0 + n\bar{\omega}_g \qquad (20)$$

gives a very poor approximation of the exact spectrum. Note that this is not surprising since, in such a case, the energy differences between successive states would remain constant, a property which is clearly wrong for the exact spectrum. In sharp contrast, we have found that modifying the partition function by exponentiating the linear contribution of the polynomials and incorporating it into the exponential part leads to remarkably good energy levels. Precisely, we propose to replace the polynomial

$$P_n(\beta) = 1 + P_{n1}\beta + \dots + P_{nn}\beta^n \qquad (21)$$

by $e^{\beta P_{n1}}$ and to define the excited-state energies of our model as

$$E_n = E_0 + n\bar{\omega}_g - P_{n1}. \qquad (22)$$

Let us emphasize that this replacement should not be considered as a quantitative approximation, but rather, as a "minimal" modification of the model to impose to the partition function to have the exact expansion, Eq.(18). Unfortunately, we have not been able to understand why this simple additional prescription to our model is so effective in giving accurate energy levels (see, the figures to follow). However, it should be considered as a salient result of this work.

To derive the functional form of the partition function, Eq.(19), we first need to expand the effective frequency $\omega_g[\tau_c(\beta)]$ as a power series of exponential-like contributions.

Let us first consider the case $\omega \neq 0$. Using Eqs.(6) and (9), and introducing the variable $y$ defined as

$$y = e^{-\beta\omega_g} \qquad (23)$$

the implicit equation obeyed by $\omega_g$ is rewritten as

$$\omega_g = \omega \sqrt{B \left[ \frac{2g}{\omega^{2m}} \omega_g^{m-1} \left( 1 + 2 \sum_{n=1}^{\infty} y^n \right)^{m-1} \right]} \qquad (24)$$

where the coth function has been expanded in power series of $y$. Note that the superscript (m) has been removed from $B^{(m)}(x)$ to simplify the notation. This will also be the case in the following for most quantities depending on $m$ when no confusion is possible. The function $B$ being infinitely differentiable, $\omega_g$ can be expanded in powers of $y$

$$\omega_g = \bar{\omega}_g + \sum_{n=1}^{\infty} \omega_n y^n. \qquad (25)$$
Remark that the equation obeyed by $\omega_g$ has no explicit dependence on $\beta$ (the dependence on the temperature is only through the variable $y$). Accordingly, the coefficients $\omega_n$ are independent of $\beta$. Their evaluation is done i) by introducing the expansion of $\omega_g$ in Eq.(24), ii) by Taylor-expanding the function $B$ at $x_0 = \frac{2g}{\omega^{2m}}\bar{\omega}_g^{m-1}$, and, finally, by identifying the contributions corresponding to a given power of $y$. After some algebra, an explicit expression for the $\omega_n$'s can be derived. For example, the first coefficient, $\omega_1$'s is given by

$$ \omega_1 = \bar{\omega}_g \frac{\frac{B_1}{B_0} x_0 (m-1)}{1 - \frac{B_1}{B_0} x_0 (m-1)} \quad (26) $$

where $B_0 = B^{(m)}(x_0)$ and $B_1 = \frac{dB^{(m)}}{dx}(x_0)$. For $n \ge 2$, $\omega_n$ can be expressed as a function of the preceding coefficients, $\omega_p$ with $1 \le p \le n-1$. The explicit formula is given in Appendix A.

When $\omega = 0$ the coefficients $\omega_k$ are more easily derived. From the general expression of the effective frequency, Eq.(8), we get $\omega_g(\tau) = c_0(\frac{g}{\tau})^{\frac{1}{4}}$ with $c_0 = \sqrt{\frac{\int dxe^{-x^4}}{\int dxx^2e^{-x^4}}}$. Using the PMS condition, Eq.(9), we obtain

$$ \omega_k = \bar{\omega}_g \sum_{l=1}^{k} 2^l \binom{\frac{1}{3}}{l} \binom{k-1}{l-1} \quad k \ge 1 \quad (27) $$

with

$$ \bar{\omega}_g = c_0^{\frac{4}{3}} \left(\frac{g}{2}\right)^{\frac{1}{3}}. \quad (28) $$

In this $\omega = 0$-case, the ground-state energy can be explicitly written as

$$ E_0 = \left[ \frac{\Gamma(\frac{5}{4})}{\Gamma(\frac{3}{4})} \right]^{\frac{2}{3}} \left( 1 - \log \sqrt{\frac{[2\Gamma(\frac{5}{4})]^3}{\pi\Gamma(\frac{3}{4})}} \right) g^{\frac{1}{3}}. \quad (29) $$

Now, to proceed, let us introduce the following new variable $y_0$

$$ y_0 = e^{-\beta \bar{\omega}_g}. \quad (30) $$

The next step consists in expressing the partition function as a power series in $y_0$, thus leading to Eq.(19). For that, we first need to derive the expansion of the effective frequency in terms of $y_0$. We have

$$ \omega_g = \bar{\omega}_g + \sum_{n=1}^{\infty} \omega_n y_0^n e^{-\beta n [\omega_g - \bar{\omega}_g]}. \quad (31) $$

By expanding the exponential in power series and by simple inspection, the form of the solution is

$$ \omega_g = \bar{\omega}_g + \sum_{n=1}^{\infty} Q_n(\beta) y_0^n \quad (32) $$

where $Q_n(\beta)$ are polynomials of degree $n-1$ in $\beta$. The polynomials can be evaluated by differentiation

$$ Q_n(\beta) = \frac{1}{n!} \frac{\partial^n \omega_g}{\partial y_0^n} (y_0 = 0). \quad (33) $$

Using Eq.(31) and the Leibniz formula for derivatives, the n-th derivative of $\omega_g$ writes

$$ \frac{1}{n!} \frac{\partial^n \omega_g}{\partial y_0^n} (y_0 = 0) = \sum_{k=1}^{n} \omega_k \frac{1}{(n-k)!} \frac{\partial^{n-k}}{\partial y_0^{n-k}} \left[ e^{-\beta k (\omega_g - \bar{\omega}_g)} \right] (y_0 = 0). \quad (34) $$

To proceed we make use of the Faà di Bruno formula

$$ \frac{\partial^n}{\partial x^n} e^{f(x)} = e^{f(x)} \sum_{m_1, m_2, \dots, m_n}^{\prime} \frac{n!}{m_1! m_2! \dots m_n!} \prod_{j=1}^{n} \left[ \frac{f^{(j)}(x)}{j!} \right]^{m_j} \quad (35) $$

where the prime on the sum indicates summing with the constraint

$$ 1m_1 + 2m_2 + \dots nm_n = n. \quad (36) $$

We then have

$$ \frac{1}{n!} \frac{\partial^n \omega_g(\beta)}{\partial y_0^n} (y_0 = 0) = \sum_{k=1}^{n} \omega_k \frac{1}{(n-k)!} $$

$$ \times \sum_{m_1, m_2, \dots, m_{n-k}}^{\prime} \frac{(n-k)!}{m_1! m_2! \dots m_{n-k}!} \prod_{j=1}^{n-k} \left[ \frac{(-k\beta)(\omega_g)^{(j)}(x)}{j!} \right]^{m_j} \quad (y_0 = 0) $$

and, finally

$$ Q_n(\beta) = \omega_n + \sum_{k=1}^{n-1} \omega_k $$

$$ \times \sum_{m_1, m_2, \dots, m_{n-k}}^{\prime} \frac{[-k\beta Q_1(\beta)]^{m_1}}{m_1!} \dots \frac{[-k\beta Q_{n-k}(\beta)]^{m_{n-k}}}{m_{n-k}!} $$

a relation which allows to evaluate iteratively the polynomials $Q_n(\beta)$. From this expression, we see that $Q_n(\beta)$ are indeed polynomials of degree $n-1$ as stated above. Let us give the first five polynomials

$$ Q_1(\beta) = \omega_1 $$

$$ Q_2(\beta) = \omega_2 - \beta \omega_1^2 $$

$$ Q_3(\beta) = \omega_3 - 3\beta\omega_1\omega_2 + \frac{3}{2}\beta^2\omega_1^3 $$

$$ Q_4(\beta) = \omega_4 - 2\beta(2\omega_1\omega_3 + \omega_2^2) + 8\beta^2\omega_2\omega_1^2 - \frac{8}{3}\omega_1^4\beta^3 $$
$$Q_5(\beta) = \omega_5 - 5\beta(\omega_2\omega_3 + \omega_1\omega_4) + \frac{25}{2}\beta^2(\omega_3\omega_1^2 + \omega_1\omega_2^2)$$

$$-\frac{125\beta^3}{6}\omega_1^3\omega_2 + \frac{125}{24}\omega_1^5\beta^4$$

We are now ready to derive the expansion of the partition function in the variable $y_0$. It is convenient to introduce the following three quantities

$$\Delta\omega_g \equiv \omega_g - \bar{\omega}_g, \qquad (37)$$

$$\Delta T \equiv \ln C(\beta) - \ln \sqrt{\frac{\bar{\omega}_g}{\pi}} \bar{I}_g, \qquad (38)$$

and

$$\Delta R \equiv \frac{1}{2}(\coth \frac{\beta\omega_g}{2} - 1). \qquad (39)$$

The partition function can be written as

$$Z = e^{-\beta E_0} P. \qquad (40)$$

Simple algebra shows that $P$ can be expressed as

$$P = e^{\beta S}[1 + \Delta R] \qquad (41)$$

where

$$S = c_1\Delta\omega_g + c_2\Delta T + c_3\Delta R + c_4\Delta\omega_g\Delta T + c_5\Delta\omega_g\Delta R + c_6\Delta T\Delta R + c_7\Delta\omega_g\Delta T\Delta R \qquad (42)$$

with the following coefficients

$$c_1 = \frac{1}{2} \left( \ln \sqrt{\frac{\bar{\omega}_g}{\pi}} \bar{I}_g - 1 \right)$$

$$c_2 = \frac{\bar{\omega}_g}{2}$$

$$c_3 = \bar{\omega}_g \ln \sqrt{\frac{\bar{\omega}_g}{\pi}} \bar{I}_g$$

$$c_4 = \frac{1}{2}$$

Let us introduce the following form for the expansion of a quantity $X$

$$X = \sum_{n=1}^{\infty} X_n(\beta) y_0^n \qquad (43)$$

where $X_n$ is a polynomial of order $n-1$ in $\beta$. Note that this form is stable by multiplication by a scalar, addition, multiplication, and exponentiation. We have already seen that $\Delta\omega_g$ admits such a representation. From the relation

$$y = y_0 e^{-\beta \Delta \omega_g} \qquad (44)$$

we see that it is also the case for $y$. The quantities $\Delta T$ and $\Delta R$ can be written as a function of the variable $y$ only. By expanding these quantities in terms of $y$ we thus find that they also both admit this representation. From (42) it is also true for $S$. Finally, the partition function, Eqs.(40,41), can be written as

$$Z = e^{-\beta E_0} \left[ 1 + \sum_{n=1}^{\infty} P_n(\beta) y_0^n \right] \qquad (45)$$

where $P_n(\beta)$ is a polynomial of order $n$ in $\beta$. This is the form given in Eq.(19).

To get the excited-state energies the quantities $P_{n1}$ need to be evaluated, see Eq.(22). Using the formula

$$P_{n1} = \frac{1}{n!} \frac{\partial^{n+1} P(\beta)}{\partial \beta \partial y_0^n} (\beta = 0, y_0 = 0)$$

we get

$$P_{n1} = S_{n0} + \Delta R_{n1} + \sum_{k+l=n} S_{k0} \Delta R_{l0}. \qquad (46)$$

where $S_{nk}$ and $\Delta_{nk}$ are the coefficients of the $n$-th order polynomial in the representation, Eq.(43), for $S$ and $\Delta R$, respectively. Using Eq.(42) and $Q_{n0} = \omega_n$ we have

$$S_{n0} = c_1\omega_n + c_2\Delta T_{n0} + c_3\Delta R_{n0} + c_4\sum_{k_1+k_2=n} \omega_{k_1}\Delta T_{k_20} + c_5\sum_{k_1+k_2=n} \omega_{k_1}\Delta R_{k_20} + c_6\sum_{k_1+k_2=n} \Delta T_{k_10}\Delta R_{k_20} + c_7\sum_{k_1+k_2=n} \left( \sum_{k_3+k_4=k_1} \omega_{k_3}\Delta T_{k_40} \right) \Delta R_{k_20} \qquad (47)$$

where $\Delta T_{nk}$ are the coefficients of the polynomial for $\Delta T$. In each sum above, $k_i \ge 1$. Let us evaluate the
quantities $\Delta R_{n0}$, $\Delta R_{n1}$, and $\Delta T_{n0}$.

Starting from

$$ \Delta R = \sum_{n=1}^{\infty} y_0^n e^{-\beta n \Delta \omega_g} $$

we have

$$ \Delta R = \sum_{n=1}^{\infty} y_0^n \sum_{k=1}^{\infty} \frac{(-\beta n)^k}{k!} \left( \sum_{l=1}^{\infty} Q_l(\beta) y_0^l \right)^k = \sum_{n=1}^{\infty} y_0^n \sum_{k=1}^{\infty} \alpha_k y_0^k $$

with

$$ \alpha_k = \sum_{l=1}^{k} \frac{(-\beta n)^l}{l!} \sum_{k_1 + \dots + k_l = k} Q_{k_1}(\beta) \dots Q_{k_l}(\beta) $$

which gives

$$ \Delta R_{n0} = 1, \Delta R_{11} = 0, \text{ and } \Delta R_{n1} = - \sum_{k=1}^{n-1} k \omega_{n-k} \text{ for } n \ge 2 \quad (48) $$

Let us now calculate $\Delta T_{n0}$. $\Delta T$, as defined by Eq.(38), is decomposed as

$$ \Delta T = \Delta A + \Delta B $$

with

$$ \Delta A = \ln \sqrt{\frac{\omega_g(\beta)}{\pi \coth \frac{\beta \omega_g[\tau_c(\beta)]}{2}}} - \ln \sqrt{\frac{\bar{\omega}_g}{\pi}} $$

and

$$ \Delta B = \ln I_g(\beta) - \ln \bar{I}_g $$

Let us begin with $\Delta A$ which can be written as

$$ \Delta A = \frac{1}{2} \ln \left[ \frac{1 + \sum_{n=1}^{\infty} \frac{Q_n(\beta)}{\bar{\omega}_g} y_0^n}{1 + \sum_{n=1}^{\infty} 2 \Delta R_n(\beta) y_0^n} \right]. $$

Expanding the two logarithmic terms, we get after some algebra

$$ \Delta A_{n0} = \frac{1}{2} \sum_{k=1}^{n} \frac{(-1)^{k-1}}{k} \sum_{l_1 + \dots + l_k = n} \left[ \left( \frac{\omega_{l_1}}{\bar{\omega}_g} \right) \dots \left( \frac{\omega_{l_k}}{\bar{\omega}_g} \right) - 2^k \right]. $$

We are now left with the calculation of $\Delta B_{n0}$. We have

$$ I_g(\beta) = \sum_{n=0}^{\infty} \frac{(-1)^n \bar{I}_n}{n!} \Delta \tau_c^n $$

with

$$ \Delta \tau_c \equiv \tau_c - \bar{\tau}_c $$

where $\bar{\tau}_c = \tau_c(\infty) = \frac{2}{\bar{\omega}_g}$ and

$$ \bar{I}_n \equiv \int_{-\infty}^{\infty} dx x^n e^{-\bar{\tau}_c V(x)}. $$

After some algebra

$$ \Delta \tau_c = \frac{2}{\bar{\omega}_g} \sum_{n=1}^{\infty} (-1)^n \left[ \sum_{k=1}^{\infty} X_k(\beta) y_0^k \right]^n $$

with

$$ X_k(\beta) = 2\Delta R_k(\beta) + \frac{Q_k(\beta)}{\bar{\omega}_g} + \frac{2}{\bar{\omega}_g} \sum_{l+m=k} \sum_{l \ge 1} \sum_{m \ge 1} Q_l(\beta) Q_m(\beta). $$

We then have

$$ \Delta \tau_c = \sum_{n=1}^{\infty} Y_n(\beta) y_0^n $$

with

$$ Y_n = \frac{2}{\bar{\omega}_g} \sum_{k=1}^{n} (-1)^k \sum_{l_1 + \dots + l_k = n} X_{l_1} \dots X_{l_k} $$

Then

$$ I_g = \bar{I}_g + \sum_{n=1}^{\infty} \frac{(-1)^n \bar{I}_n}{n!} \left( \sum_{k=1}^{\infty} Y_k(\beta) y_0^k \right)^n $$

$$ I_g = \bar{I}_g + \sum_{n=1}^{\infty} Z_n(\beta) y_0^n $$

$$ Z_n = \sum_{k=1}^{n} (-1)^k \frac{\bar{I}_k}{k!} \sum_{l_1 + \dots + l_k = n} Y_{l_1} \dots Y_{l_k} $$

Thus

$$ \ln I_g = \ln \bar{I}_g + \ln \left( 1 + \sum_{n=1}^{\infty} \frac{Z_n}{\bar{I}_g} y_0^n \right). $$

Finally,

$$ \ln I_g = \ln \bar{I}_g + \sum_{n=1}^{\infty} \Delta B_n(\beta) y_0^n $$

with

$$ \Delta B_n = \sum_{k=1}^{n} \frac{(-1)^k}{k \bar{I}_g^k} \sum_{l_1 + \dots + l_k = n} Z_{l_1} \dots Z_{l_k}. $$

The quantities $\Delta B_{n0}$ are then evaluated by taking $\beta = 0$ in the preceding expressions.
Figure 1 shows the ground-state energy and the first eight excited state energies of the quartic oscillator with $\omega = 1$ versus the coupling constant $g$. The "exact" numerical energies, obtained by diagonalization of the Hamiltonian in a sufficiently large Gaussian basis set, are shown as solid lines. At the scale of the figure, the energy levels and their overall behavior as functions of $g$ and $n$ are well reproduced and appear satisfactory. Quantitative results are presented in Table I, which reports the relative errors $\epsilon$ in the computed energies. We present results for three different potentials: the first two are the quartic ($m = 2, \omega = 1$) and sextic ($m = 3, \omega = 1$) anharmonic potentials, respectively. The third is the "pure" quartic potential, defined by the absence of a quadratic term ($\omega = 0$). Results are given for five values of $g$ ($g = 0.1, 1, 10, 40, \text{ and } 400$), spanning the weak- to strong-coupling regimes. A first remark is that relative errors are all negative. In other words, the computed energies are always smaller than the exact ones. Unfortunately, we were not able to understand the origin of this interesting observed property. A second remark is that relative errors are all of the order of a few percent (up to 10% in the worst case of the sextic potential at large $g$). We consider this level of accuracy as quite satisfactory in view of the simplicity of the model. Let us insist on the fact that the model partition function has a particularly simple and compact form, Eq.(10) with Eqs.(11) and (12); and, most importantly, no adjustable parameters. In addition, the computational cost for calculating the energies is negligible. Another remark is that for each potential the relative error is maximal for the two lowest energies. Quite remarkably, the errors on the higher energies ($n \ge 2$) are almost constant and, also, nearly independent on the value of $g$. This is a quite interesting feature of the model. The comparison between the three different potentials is instructive. As expected, the accuracy reached for the sextic oscillator is inferior to that obtained for the quartic oscillator, a consequence of the greater anharmonicity of the former potential. This is also the case at small values of $g$ for the pure quartic oscillator in which the quadratic term has been removed. However, for large $g$'s the errors become nearly identical. A striking feature of the energies of the pure quartic potential is that the relative errors are independent on $g$ (actually, all digits of the errors are identical for the different $g$'s, a result not shown here). This result is explained as follows. By a simple rescaling of the Schrödinger equation, the exact energies $E_n(g)$ of the pure quartic oscillator can be shown to scale as $g^{1/3}$. Quite satisfactorily, it is also the case for our model. This can be proved by noting that $\bar{\omega}_g$, Eq.(28), and the $\omega_k$'s, Eq.(27) scale also as $g^{1/3}$ and, then, by invoking the series of equations, (22),(46),(47) and (48). The exact and model energies having the same scaling in $g$, the relative errors are thus independent on $g$.

<figure><img src="image_2.png" /><figcaption>Figure 1. Energy levels of the quartic oscillator as a function of $g$ for $n = 0$ to $n = 8$. Exact results given by the solid lines. Harmonic frequency, $\omega = 1$</figcaption></figure>

## 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).
<table><caption>Table I. Relative errors $\epsilon$ (in %) on the energy levels $E_n$ for different values of $g$ and potentials $V_{[a,n]}(x) = ax^2 + gx^n$; $\epsilon = (E_{\text{model}} - E_{\text{ex}})/E_{\text{ex}}$.</caption><thead><tr><th></th><th>E<sub>0</sub></th><th>E<sub>1</sub></th><th>E<sub>2</sub></th><th>E<sub>3</sub></th><th>E<sub>4</sub></th><th>E<sub>5</sub></th><th>E<sub>6</sub></th><th>E<sub>7</sub></th><th>E<sub>8</sub></th></tr></thead><tbody><tr><td>g=0.1</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>V<sub>(1/2),4</sub></td><td>-0.5</td><td>-0.4</td><td>-0.5</td><td>-0.6</td><td>-0.7</td><td>-0.8</td><td>-0.8</td><td>-0.9</td><td>-0.9</td></tr><tr><td>V<sub>(1/2),6</sub></td><td>-2.5</td><td>-2.6</td><td>-2.7</td><td>-2.9</td><td>-3.1</td><td>-3.2</td><td>-3.3</td><td>-3.4</td><td>-3.4</td></tr><tr><td>V<sub>(0,4)</sub></td><td>-4.2</td><td>-2.0</td><td>-1.6</td><td>-1.7</td><td>-1.7</td><td>-1.7</td><td>-1.7</td><td>-1.7</td><td>-1.7</td></tr><tr><td>g=1</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>V<sub>(1/2),4</sub></td><td>-2.3</td><td>-1.4</td><td>-1.3</td><td>-1.4</td><td>-1.4</td><td>-1.4</td><td>-1.4</td><td>-1.5</td><td>-1.5</td></tr><tr><td>V<sub>(1/2),6</sub></td><td>-6.3</td><td>-4.6</td><td>-3.7</td><td>-3.7</td><td>-3.7</td><td>-3.7</td><td>-3.8</td><td>-3.8</td><td>-3.8</td></tr><tr><td>V<sub>(0,4)</sub></td><td>-4.2</td><td>-2.0</td><td>-1.6</td><td>-1.7</td><td>-1.7</td><td>-1.7</td><td>-1.7</td><td>-1.7</td><td>-1.7</td></tr><tr><td>g=10</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>V<sub>(1/2),4</sub></td><td>-3.7</td><td>-1.9</td><td>-1.6</td><td>-1.6</td><td>-1.6</td><td>-1.6</td><td>-1.6</td><td>-1.6</td><td>-1.6</td></tr><tr><td>V<sub>(1/2),6</sub></td><td>-9.0</td><td>-5.4</td><td>-4.0</td><td>-4.0</td><td>-4.0</td><td>-3.9</td><td>-3.9</td><td>-3.9</td><td>-3.9</td></tr><tr><td>V<sub>(0,4)</sub></td><td>-4.2</td><td>-2.0</td><td>-1.6</td><td>-1.7</td><td>-1.7</td><td>-1.7</td><td>-1.7</td><td>-1.7</td><td>-1.7</td></tr><tr><td>g=40</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>V<sub>(1/2),4</sub></td><td>-4.0</td><td>-2.0</td><td>-1.6</td><td>-1.7</td><td>-1.6</td><td>-1.6</td><td>-1.6</td><td>-1.6</td><td>-1.6</td></tr><tr><td>V<sub>(1/2),6</sub></td><td>-9.0</td><td>-5.5</td><td>-4.1</td><td>-4.1</td><td>-4.1</td><td>-4.0</td><td>-4.0</td><td>-4.0</td><td>-4.0</td></tr><tr><td>V<sub>(0,4)</sub></td><td>-4.2</td><td>-2.0</td><td>-1.6</td><td>-1.7</td><td>-1.7</td><td>-1.7</td><td>-1.7</td><td>-1.7</td><td>-1.7</td></tr><tr><td>g=400</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>V<sub>(1/2),4</sub></td><td>-4.2</td><td>-2.0</td><td>-1.6</td><td>-1.7</td><td>-1.7</td><td>-1.7</td><td>-1.7</td><td>-1.7</td><td>-1.7</td></tr><tr><td>V<sub>(1/2),6</sub></td><td>-10.</td><td>-5.7</td><td>-4.2</td><td>-4.1</td><td>-4.0</td><td>-4.0</td><td>-4.0</td><td>-4.0</td><td>-4.0</td></tr><tr><td>V<sub>(0,4)</sub></td><td>-4.2</td><td>-2.0</td><td>-1.6</td><td>-1.7</td><td>-1.7</td><td>-1.7</td><td>-1.7</td><td>-1.7</td><td>-1.7</td></tr></tbody></table>

## Appendix A: Formula for evaluating $\omega_n$'s

In this appendix, we give the explicit formulas for the coefficients of the expansion of $\omega_g$ in powers of $y$ written as

$$ \omega_g = \bar{\omega}_g + \sum_{n=1}^{\infty} \omega_n y^n. \qquad (A1) $$

Note that to simplify the notation the dependence on $m$ of the various quantities considered here will be omitted.

Let us define now the following quantities

$$ B_n \equiv \frac{d^n B^{(m)}}{dx^n}(x_0) \qquad (A2) $$

$$ x_0 = \frac{2g}{\omega^{2m}} \bar{\omega}_g^{m-1} \qquad (A3) $$

$$ \alpha_n = \frac{1}{n!} \frac{B_n}{B_0} x_0^n \qquad (A4) $$

and

$$ \beta_n = \frac{\frac{1}{2}(\frac{1}{2}-1)\dots(\frac{1}{2}-(n-1))}{n!} \qquad (A5) $$

$\bar{\omega}_g$ is given by

$$ \bar{\omega}_g = \omega \sqrt{B_0} \qquad (A6) $$

The first coefficient $\omega_1$ is given by

$$ \omega_1 = A\alpha_1(m-1) \qquad (A7) $$

with

$$ A = \frac{\bar{\omega}_g}{1 - \frac{\alpha_1(m-1)}{2}} \qquad (A8) $$

For $n \ge 2$ the $\omega_n$'s are function of the previous coefficients $\omega_p$ with $1 \le p \le n-1$ as follows

$$ \omega_n = A \left[ \alpha_1(m-1)S_n + \frac{1}{2}\alpha_1 T_n + U_n + V_n \right] \qquad (A9) $$

where

$$ S_n = 1 + \sum_{k+l=n} \frac{\omega_k}{\bar{\omega}_g} \qquad (A10) $$

$$ T_n = \sum_{k=2}^{m-1} \binom{m-1}{k} \sum_{l_1+\dots+l_k=n} a_{l_1}\dots a_{l_k} \qquad (A11) $$

with

$$ a_n = \frac{\omega_n}{\bar{\omega}_g} + 2S_n \quad n \ge 1 \qquad (A12) $$

$$ U_n = \frac{1}{2} \sum_{k=2}^{n} \alpha_k \sum_{l_1+\dots+l_k=n} b_{l_1} \dots b_{l_k} \qquad (A13) $$

with

$$ b_n = (m-1)a_n + T_n \quad n \ge 1 \qquad (A14) $$

and

$$ V_n = \sum_{k=2}^{n} \beta_k \sum_{l_1+\dots+l_k=n} c_{l_1}\dots c_{l_k} \qquad (A15) $$

with

$$ c_n = \alpha_1 b_n + 2U_n \quad n \ge 1 \qquad (A16) $$

When no index fulfills the constraint in the sum, the corresponding quantity is equal to zero. Here, it means that $S_1 = 1$, $T_1 = 0$, $U_1 = 0$, $V_1 = 0$, and $T_n = 0$ for $m=2$.
[1] A. Turbiner and J.C. Del Valle Rosale. *Quantum Anharmonic Oscillator*. World Scientific, 2023.

[2] Alexander V. Turbiner. One-dimensional quasi-exactly solvable schrödinger equations. *Physics Reports*, 642:1-71, 2016.

[3] Daniel Baye. The lagrange-mesh method. *Physics Reports*, 565:1-107, 2015. The Lagrange-mesh method.

[4] J. C. del Valle. Solving the one-dimensional time-independent schrödinger equation with high accuracy: The lagrangemesh mathematica® package. *International Journal of Modern Physics C*, 35(01):2450011, 2024.

[5] Carl M. Bender, Kaare Olaussen, and Paul S. Wang. Numerological analysis of the wkb approximation in large order. *Phys. Rev. D*, 16:1740-1748, Sep 1977.

[6] R. N. Kesarwani and Y. P. Varshni. Eigenvalues of an anharmonic oscillator. *Journal of Mathematical Physics*, 22(9):1983-1989, 09 1981.

[7] M. Lakshmanan, F. Karlsson, and P. O. Fröman. Phase-integral calculation of the energy levels of a quantal anharmonic oscillator. *Phys. Rev. D*, 24:2586-2598, Nov 1981.

[8] C.C. Yao and K.S. Cheng. A simple variational method for excited states of coupled anharmonic oscillators. *Nuov Cim B*, 111:645-652, 1996.

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

[10] H. Kleinert. Improving the variational approach to path integrals. *Physics Letters B*, 280(3):251-260, 1992.

[11] Ranjan Bhattacharya, Dhiranjan Roy, and Siddhartha Bhowmick. Simple systematics in the ground state energies of quantum anharmonic oscillators. *Physics Letters A*, 244(1):9-12, 1998.

[12] Ananda Dasgupta, Dhiranjan Roy, and Ranjan Bhattacharya. Simple systematics in the energy eigenvalues of quantum anharmonic oscillators. *Journal of Physics A: Mathematical and Theoretical*, 40(4):773, jan 2007.

[13] David Leonard and Paul Mansfield. Solving the anharmonic oscillator: tuning the boundary condition. *Journal of Physics A: Mathematical and Theoretical*, 40(33):10291, aug 2007.

[14] M. Companys Franzke, A. Tichai, K. Hebeler, and A. Schwenk. Excited states from eigenvector continuation: The anharmonic oscillator. *Physics Letters B*, 830:137101, 2022.

[15] M. Caffarel. Path integral for the quartic oscillator: An accurate analytic expression for the partition function, arxiv quant-ph arxiv:2312.09859v3., 2023.

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

[17] 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.

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

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