# Asymptotics of eigenvalue sums when some turning points are complex

Pavel Okun Department of Chemistry, University of California, Irvine, CA
92617, USA Kieron Burke Departments of Physics and Astronomy and of Chemistry,
University of California, Irvine, CA 92617

###### Abstract

Recent work has shown a deep connection between semilocal approximations in
density functional theory and the asymptotics of the sum of the WKB
semiclassical expansion for the eigenvalues. However, all examples studied to
date have potentials with only real classical turning points. But systems with
complex turning points generate subdominant terms beyond those in the WKB
series. The simplest case is a pure quartic oscillator. We show how to
generalize the asymptotics of eigenvalue sums to include subdominant
contributions to the sums, if they are known for the eigenvalues. These
corrections to WKB greatly improve accuracy for eigenvalue sums, especially
for many levels. We obtain further improvements to the sums through
hyperasymptotics. For the lowest level, our summation method has error below
$2\times 10^{-4}$. For the sum of the lowest 10 levels, our error is less than
$10^{-22}$. We report all results to many digits and include copious details
of the asymptotic expansions and their derivation.

##  1 Introduction

This work grew out of an effort to derive systematic corrections to local
density functional approximations. The ultimate goal is to put density
functional theory (DFT) on a similar theoretical foundation to wavefunction-
based quantum chemistry, where it is well understood how to achieve greater
accuracy (often starting from the Hartree-Fock approximation [1]) with
increasing computational cost [2]. The analagous starting point for density
functional approximations are local density approximations but, because the
expansion is an asymptotic semiclassical expansion, the mathematics required
to perform this expansion is much less straightforward.

The connection between local density and semiclassical approximations has a
long history [3, 4]. Work over the past 16 years has attempted to find that
connection for the exchange-correlation (XC) energy of Kohn-Sham DFT [5, 6, 7,
8], to provide insight for XC approximations [9, 10]. Recently, in simple one-
dimensional systems, it has been possible to derive the leading corrections to
local approximations from semiclassical expansions [11, 12, 13]. In typical
cases, this expansion is asymptotic and so requires careful treatment [14,
15]. The basic idea is that local density functional approximations are the
leading-order terms in an asymptotic, semiclassical expansion and today’s
approximate functionals, such as generalized gradient approximations, are
crude attempts to approximate higher order corrections.

We work directly with the total energy as a functional of the potential, but
approximations can ultimately be converted into density functionals [16, 17]
for the 1D kinetic energy (whose three-dimensional counterpart is vital for
orbital-free DFT [18]). The key is to approximate the sum of the first $N$
eigenvalues, $E_{N}$ (in chemistry parlance this is the ground state energy of
$N$ non-interacting same-spin fermions). The connection between this sum and
density functionals is explained in Refs. [12, 19].

The key idea is to generate an asymptotic expansion for the sum of eigenvalues
directly from the known WKB series for individual eigenvalues. Such
asymptotics of sums are much more accurate than summing the asymptotic
approximations for the individual eigenvalues. Since this method for summing
eigenvalues uses the WKB series [20], we call it SWKB, where the ”S” stands
for sum. It was developed in Refs. [13, 11, 12]. In simple cases (particle in
a box and harmonic oscillator), the semiclassical expansion truncates at low
order, yielding exact results for both the eigenvalues and their sums [12].
Reference [11] analyzed the Pöschl-Teller well, where the WKB and SWKB series
are convergent. Reference [13] tested the SWKB approximation on a linear
potential with a hard wall (the linear half-well), which is truly asymptotic.

But the linear half-well has only real classical turning points, whereas many
other simple analytic forms allow for complex turning points. The quartic
oscillator is the simplest, and it has long been known to have sub-dominant
contributions to its semiclassical eigenvalue expansion, which reduce errors
significantly relative to WKB to the same order. So here we apply the
asymptotic summation method to the quartic oscillator. We show how the
subdominant corrections to the WKB eigenvalues lead to subdominant corrections
to the SWKB sums, and generate an asymptotic expansion for these corrections.
We find that the subdominant corrections to the sums likewise greatly improve
accuracy. This methodology can be easily applied to any such case.

More broadly, the quartic oscillator is of interest in particle physics for
testing model field theories [21, 22]. In mathematical physics it is a model
asymptotic system [23]. The quartic oscillator Schrödinger equation is a
special case of Huen’s differential equation (Sec. 31.2 of Ref. [24]) and its
exact solution was studied in this context in Refs. [25, 26]. While we focus
on the asymptotic approximation to all eigenvalues Ref. [27] uses a different
technique to approximate the ground state wavefunctions and energies.

In Sec. 3 we thoroughly analyze the semiclassical expansion of the quartic
oscillator eigenvalues to many orders. While the basic structure of this
expansion (WKB plus subdominant corrections) has long been known [28, 29, 30],
we examine it more thoroughly and systematically than typical studies in the
literature. We find asymptotic expansions beyond 20 orders, subdominant
corrections to such expansions, 20 eigenvalues given to 40 digits, and in some
cases errors less than $10^{-20}$. We also give explicit formulas for the
exact analytic expressions for the lower-order asymptotic coefficients, and
the formulas to generate the exact WKB series and its inversion to all orders,
making it easy for others to use our work.

But all this is a prelude to the main purpose of our work, which is to apply
newly developed methods of finding sums of the lowest $N$ eigenvalues of
quantum problems. In Sec. 4 we asymptotically sum over the explicit WKB
expansion of Sec. 3, finding its asymptotic coefficients. The crucial point is
that, from the subdominant corrections to WKB, we can easily find the
subdominant corrections to the WKB asymptotic expansion for the sums, which
are far more accurate than their WKB counterparts. We also use tricks of
hyperasymptotics for the sums, to further improve accuracy. This is the most
accurate approximation which we derive. Finally, we take differences between
approximations for sums to approximate the individual eigenvalues, which we
can again compare to the WKB series including subdominant corrections.

While our results should be of general interest to anyone studying
semiclassical expansions of eigenenergies in any context, our driving interest
is in density functional theory (DFT). Thus in Sec. 5 we summarize our
findings and explain their relevance to DFT.

##  2 Notation and conventions

We choose units so that our Hamiltonian is

$\hat{\mathscr{H}}=-\frac{1}{2}\frac{d^{2}}{dx^{2}}+\frac{x^{4}}{2}.$ (1)

Many benchmark results to many digits were reported in Ref. [19], with a
different constant. So we give the first 20 eigenvalues to 41 decimal places
in Table S1. References [31, 20, 32] also report eigenvalues for this system,
but to less digits.

We adopt the following notation for asymptotic series. Define the function

$B_{M}(y)=\sum_{m=0}^{M}b_{m}(y)=\sum_{m=0}^{M}\frac{b_{m}}{y^{m}},$ (2)

i.e., given a set of coefficients $b_{m}$, $B_{M}(y)$ is the $M$-th order
polynomial in $1/y$. We choose $b_{0}=1$ always, so that $b_{0}(y)=1$ and
$B_{M}(y)\to 1$ as $y\to\infty$. Thus the dominant asymptotic behavior is
always given by setting the polynomial to $1$. We refer to each $b_{m}(y)$ as
the $m$-th addition.

We will often use optimal truncation to create approximations that are more
accurate than any of fixed order (superasymptotics) [15]. We consider three
different choices. Empirical ($E$) optimal truncation implies truncation at
the order yielding the lowest error. This requires knowing which function is
being expanded. Least-addition ($L$) optimal truncation means truncating at
that value of $m$ for which $|b_{m}(y)|$ is least. Finally, asymptotic ($A$)
optimal truncation uses the (usually simple but not always available) formula
for least addition as $m\to\infty$, using the asymptotically dominant
expression for $b_{m}(y)$ as $m\to\infty$. All three must agree for
sufficiently large $y$, but may disagree for moderate values of $y$. We
generate at least 7 different asymptotic expansions in this paper, and for
ease of use the coefficients of the most important expansions are given in
Table 9 of Appendix D to 13 digits, but 41 digits are given in the
supplemental info.

Throughout, we define error as the difference between approximate and exact.

##  3 Semiclassical expansion of quartic oscillator eigenvalues

In this section, we thoroughly review the known WKB series for the quartic
oscillator, and its leading subdominant corrections. All the essential facts
already appear in various places in the literature [20, 28, 28, 29, 30], but
we collect the relevant ones here and are more thorough.

The WKB series for the quartic oscillator eigenvalues is well-known:

$\left(\frac{\epsilon}{\alpha}\right)^{3/4}A_{M}[(2\epsilon)^{3/2}]=z,$ (3)

where $\gamma=\Gamma(1/4)\approx 3.62561$,
$\alpha=3\pi^{2}[3/(2\gamma^{8})]^{1/3}\approx 1.09253$ and eigenvalues
$\epsilon_{j}$ occur at $z=j+1/2$, $j=0,1,2,\cdots$. We derive Eq. (3) in
Appendix A. The first 6 $a_{n}$ are reported in Bender and Orszag (where
$A_{2n}=A_{0}a_{n}$). We give the $a_{n}$ exactly up to $n=28$ in Table S3.
Like many of the series in this paper, Eq. (3) is an asymptotic expansion and
never converges. The coefficients at first become smaller, but after $m=2$,
they grow rapidly, with asymptotic behavior [28]

$a_{n}^{(0)}=-(-1)^{{\lfloor}n/2{\rfloor}}\frac{2}{\pi}\left(\frac{9\pi}{\gamma^{4}}\right)^{n}(2n-2)!,\qquad
n\rightarrow\infty.$ (4)

Their signs also oscillate, with pairs of positive values followed by pairs of
negative values. The $a_{n}$ appear in Table 9 and to more digits in Table S2.
This table also compares the asymptotic approximation in Eq. (4) with the
exact $a_{n}$.

$\begin{array}[]{|c|c|c|r|r|r|r|r|}\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Error}\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr
j&M&L&\lx@intercol\hfil\epsilon\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{WKB}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{BCWKB}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\epsilon^{\text{WKB}}_{M}+\Delta\epsilon_{L}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{CWKB}\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr 0&1&1&0.5301810&-3.5270398\times 10^{-02}&1.8285750\times
10^{-02}&1.7091440\times 10^{-02}&1.4571908\times 10^{-02}\\\
1&2&3&1.8998365&4.3044942\times 10^{-03}&5.6300300\times
10^{-05}&5.2731082\times 10^{-05}&5.2724811\times 10^{-05}\\\
2&6&6&3.7278490&-2.2822571\times 10^{-04}&-1.9645769\times
10^{-07}&-1.9824144\times 10^{-07}&-1.9979145\times 10^{-07}\\\
3&8&8&5.8223728&1.1248072\times 10^{-05}&1.2933375\times
10^{-09}&1.2905883\times 10^{-09}&1.2914117\times 10^{-09}\\\
4&10&11&8.1309130&-5.3436485\times 10^{-07}&-4.1292807\times
10^{-12}&-4.1333145\times 10^{-12}&-4.1333154\times 10^{-12}\\\
5&12&13&10.6191865&2.4862367\times 10^{-08}&-3.3798630\times
10^{-14}&-3.3805755\times 10^{-14}&-3.3805755\times 10^{-14}\\\
6&14&15&13.2642356&-1.1414056\times 10^{-09}&1.0688951\times
10^{-15}&1.0688829\times 10^{-15}&1.0688829\times 10^{-15}\\\
7&16&17&16.0492989&5.1916633\times 10^{-11}&-1.8636440\times
10^{-17}&-1.8636461\times 10^{-17}&-1.8636461\times 10^{-17}\\\
8&18&19&18.9615005&-2.3453802\times 10^{-12}&2.7529489\times
10^{-19}&2.7529485\times 10^{-19}&2.7529485\times 10^{-19}\\\
9&21&21&21.9905790&1.0540356\times 10^{-13}&7.6509652\times
10^{-20}&7.6509652\times 10^{-20}&7.6509652\times 10^{-20}\\\
10&23&23&25.1281273&-4.7174532\times 10^{-15}&-8.9559882\times
10^{-22}&-8.9559882\times 10^{-22}&-8.9559882\times 10^{-22}\\\
\hline\cr\end{array}$

Table 1: Errors in WKB series, including various approximations for
subdominant contributions: WKB [Eq. (5)], BCWKB [Eq. (11)],
$\epsilon^{\text{WKB}}_{M}+\Delta\epsilon_{L}$ [Eq. (13)], and CWKB [Eqs. (14)
& (15)]. $M$ is the order at which we truncate, and is one less than the order
of least addition. $L$ is the order of least addition for Eq. (13).

Equation (3) yields an implicit definition of $\epsilon(z)$ at each order $M$.
We can invert the expansion order-by-order to find

$\epsilon^{\text{WKB}}_{M}(z)=\alpha\,z^{4/3}\,B_{M}(z^{2}),$ (5)

where the coefficients $b_{n}$ are also listed numerically in Table 9 and
given to more digits in Table S5. Since Eq. (5) is the inversion of Eq. (3)
the $\\{b_{n}\\}$ are linear combinations of the $\\{a_{n^{\prime}}\\}$ with
$n^{\prime}\leq n$. We discuss how to generate the $b_{n}$ and their exact
forms in Appendix B. Assuming that the dominant contribution to the $b_{n}$ is
the term with the highest order $a_{n}$ in Eq. (38) yields

$b_{n}^{(0)}=(-1)^{{\lfloor}n/2{\rfloor}}\frac{8(2n-2)!}{3\pi(2\pi^{2})^{n}},\qquad
n\rightarrow\infty.$ (6)

Table S5 gives the ratio of $b_{n}$ coefficients to this asymptotic form, and
strongly suggests it is correct.

Figure 1: Errors (approximate minus exact) of the WKB eigenvalue series as a
function of order, $M$, for the ground state (black) to the sixth excited
state (gray). The sign of the error is denoted by the symbol (X’s positive,
full circles negative).

Figure 1 plots the errors for Eq. (5) as a function of order for the first 7
eigenvalues. While the WKB series is extremely accurate and improves as $n$
increases, plateaus develop of (approximately) fixed error, and these widen as
$n$ increases. Thus there is a range of values of $M$ that all produce
comparable error. On the other hand the WKB series for the linear half well in
Refs. [12, 13] smoothly decreases with order down to some minium value before
smoothly increasing with increasing order.

Figure 2: WKB expansion errors (blue) and additions (black) for $j=4$, with
conventions of Fig. 1.

To see this in more detail we consider the fourth excited state ($j=4$) in
Fig. 2 and include the magnitude of the additions
[$\epsilon^{\text{WKB}}_{M}(z)-\epsilon^{\text{WKB}}_{M-1}(z)$]. On the
plateau, the additions become much smaller than the error. To the right, the
additions grow and dominate the error. This strongly suggests that something
is missing from the WKB series itself.

For this problem, the WKB series is known not to give the complete
semiclassical expansion of the eigenvalues. The WKB series accounts for only
the real turning points, but there are also two imaginary turning points in
the complex plane. For the linear half well [13] the WKB series gave the full
semiclassical expansion of the eigenvalues because there was only one real
turning point (the hard wall only kills off the even states of the linear well
and does not count as a regular turning point [19]). For the quartic
oscillator, the true semiclassical expansion has additional subdominant pieces
which the WKB series neglects. Voros et. al [28] found a quantization rule
that, unlike Eq. (3), which considers only real turning points, takes all four
turning points (both real and imaginary) into account:

$z=(\epsilon/\alpha)^{3/4}A_{M}[(2\epsilon)^{3/2}]-g_{jM}(\epsilon),$ (7)

where

$g_{jM}(\epsilon)=\frac{(-1)^{j}}{\pi}\text{atan}\exp\left[-\pi(\epsilon/\alpha)^{3/4}A_{M}[-(2\epsilon)^{3/2}]\right].$
(8)

This more general quantization rule [28, 29, 30, 33, 34] comes from the
complex WKB method [35] which extends the WKB expansion into the complex
plane.

$\begin{array}[]{|c|c|r|}\hline\cr j&L&\lx@intercol\hfil\epsilon_{L}^{\rm
SD}(j+1/2)\hfil\lx@intercol\vrule\lx@intercol\\\ \hline\cr
0&2&4.98423055821\times 10^{-02}\\\ 1&7&-4.25176939432\times 10^{-03}\\\
2&11&2.28025923533\times 10^{-04}\\\ 3&19&-1.12467805687\times 10^{-05}\\\
4&23&5.34360714255\times 10^{-07}\\\ 5&27&-2.48624007875\times 10^{-08}\\\
6&31&1.14140669578\times 10^{-09}\\\ 7&35&-5.19166517060\times 10^{-11}\\\
8&39&2.34538047637\times 10^{-12}\\\ 9&43&-1.05403481095\times 10^{-13}\\\
10&47&4.71745228714\times 10^{-15}\\\ \hline\cr\end{array}$

Table 2: Subdominant corrections to WKB eigenvalues of Eq. (14), truncated at
least addition, L. Figure 3: WKB (black) and corrected WKB errors (red) after
optimal truncation, from Table 1. All corrected approximations are
indistinguishable, conventions as in Fig. 1.

Noting that $g$ is much smaller than any power of $z$ in the large $z$ limit,
the net effect of $g$ is simply to shift the $z$ value. Then repeating the
inversion which gave us Eq. (5) from Eq. (3), we find

$\epsilon_{MM^{\prime}}(z)=\epsilon^{\text{WKB}}_{M}[z+g_{M^{\prime}}(z)].$
(9)

Recalling that $j=z-1/2$, we rewrite Eq. (8) as

$g_{M^{\prime}}(z)=\frac{(-1)^{j}}{\pi}\text{atan}~{}\exp[-\pi\,z\,D_{M^{\prime}}(z^{2})],$
(10)

where we find the $d_{n}$ by plugging the WKB series in Eq. (5) into Eq. (8)
and rearranging the double sum into a single sum by grouping terms of like
powers in $z$. The $d_{n}$ are reported in Table 9 and more precisely in Table
S8. Analytic expressions for the $d_{n}$ are given in Appendix B.

Figure 4: Same as Fig. 2, but adding the BCWKB (red) approximation of Eq.
(11), using subdominant corrections truncated via least addition

.

There is no necessary connection between $M$ and $M^{\prime}$ in Eq. (9): each
can be chosen independently. Here $M$ is the order at which we truncate the
WKB series, and has a large effect on the accuracy of our approximations
beyond WKB. Thus we choose $M$ very carefully. The order $M^{\prime}$
determines the subdominant (SD) correction to the WKB values. This is the
first of several series we will derive for the SD values. We use least
addition to choose the order for all these SD series. Varying $M^{\prime}$
around the order of least addition has little effect on accuracy. We label
this BCWKB, for bare corrected WKB:

$\epsilon^{\rm BCWKB}(z)=\epsilon^{\text{WKB}}_{M}[z+g_{L}(z)],$ (11)

where $g_{L}(z)$ is optimally truncated via least addition. Equation 11 yields
our best estimate for the individual eigenvalues, Table S18 shows the
$D_{L}(z^{2})$ and the values of $L$, which we find to be $L=2j+4$ for all
$j>2$. We plug these $D_{L}$ into Eq. (10) to estimate $g$ for each $j$, and
then evaluate $\epsilon^{\rm BCWKB}(z)$. This yields the corrections shown in
Fig. 4.

This figure shows extraordinary improvement over the bare WKB series. To
determine the asymptotic ($j\rightarrow\infty)$ order of least addition for
the WKB series, we approximate the $b_{n}$ with Eq. (6) and minimize the
terms, the $b_{n}(z^{2})$, in Eq. (5) with respect to $M$:

$A(j)=\bigg{\lfloor}1+\frac{\pi}{\sqrt{2}}\left(j+\frac{1}{2}\right)\bigg{\rfloor}.$
(12)

The asymptotic and least addition orders of optimal truncation are identical
for $j\geq 5$. However, the empirical least error is one order less, for
reasons we have been unable to fathom. The figure, made for $j=4$, has $L=11$,
but the actual least error is at $M=10$, by almost two orders of magnitude.
Table 1 documents the enormous improvement of Eq. (11) over WKB. For the first
level, $L=2$, and the improvement is only a factor of 2. For $j=10$, the
improvement is 7 orders of magnitude. The prefactor of $j$ in Eq. (12),
$\pi/\sqrt{2}$, explains why our accuracies are lower than those of Ref. [13].
For both eigenvalues and their sums in the next section, logs of errors are
$1/\sqrt{2}$ smaller than those of Ref. [13] as $j\rightarrow\infty$.

Last in this section, we find the asymptotic behavior of the subdominant
corrections to WKB. We expand around $g=0$, i.e., the WKB series, to find the
difference with WKB to be

$\Delta\epsilon_{M}(z)\approx
g_{L}(z)\frac{d\epsilon^{\text{WKB}}_{M}}{dz}=g_{L}(z)\alpha
z^{1/3}B_{M}^{\prime}(z^{2}),$ (13)

where $b_{m}^{\prime}=(4/3-2m)b_{m}$. In Table 1, we show that adding this to
the WKB series yields results about as accurate as those of Eq. (11). Now our
expression is linear in $g(z)$. We next expand all functions in $g(z)$ in
inverse powers of $1/z$, except for the dominant exponential, $\exp(-\pi z)$,
yielding

$\epsilon^{\rm SD}_{M}(z)=(-1)^{j}\frac{4\alpha}{3\pi}z^{1/3}e^{-\pi
z}H_{M}(z),$ (14)

where the superscript SD refers to the leading subdominant contribution. We
define the corrected-WKB (CWKB) approximation:

$\epsilon^{\text{CWKB}}_{M}(j)=\epsilon^{\text{WKB}}_{M}(j+1/2)+\epsilon^{\rm
SD}_{L}(j+1/2),$ (15)

where the SD corrections are least-addition optimally truncated. For $j>2$, we
find $L=4j+7$ (Tables 2 and S19). In Table 1 we show that Eqs. (11), (13), and
(14) all produce results of comparable accuracy, with differences becoming
negligible with increasing level. We give these subdominant energies in Table
2 and more fully in Table S19. We give the $h_{n}$ numerically in Table 9,
more precisely in Table S11, and analytically in Appendix B. Equation (15)
yields about the same errors as Eq. (13) in Table 1 to the number of digits
shown in that table. In Table S21 we show that WKB is relatively insensitive
to the order of truncation, but CWKB is extremely sensitive, so we choose
$M=L-1$ in Eq. (15).

Our treatment of the subdominant corrections to the WKB series is simpler than
that employed in modern particle theory [36, 37], but suffices for our study
of the asymptotic approximation to the eigenvalue sum. For those interested in
the techniques of modern asymptotic theory and it’s application to quantum
oscillators we refer them to Ref. [38].

##  4 Asymptotics of sums

We now turn to the main topic of this work, the asymptotics of sums of
eigenvalues, and whether or not subdominant contributions can also be
incorporated in such approximate sums. Given the extraordinary accuracy of
modern DFT approximations (energy differences can be accurate to within one
part in 109 of total energies in Kohn-Sham (KS) calculations), the asymptotic
DFT for model systems developed in Refs. [11, 12, 13] might not be sufficient
in the presence of subdominant corrections. It is also of intellectual
interest to see if they can be handled.

Figure 5: Errors in WKB approximation for sums (SWKB) from $N=1$ (black) to
$N=7$ (gray), conventions as in Fig. 1.

The sum of the first $N$ eigenvalues is

$E(N)=\sum_{j=0}^{N-1}\epsilon_{j}.$ (16)

We can derive the SWKB approximation to this sum using any of the methods in
Refs. [12, 11, 13]. We choose to start with Eq. (2.4) of Ref. [13]:

$E_{R}(N)=E(\infty)-\int_{N}^{\infty}dj\,\epsilon_{j}-\frac{\epsilon_{N}}{2}+\sum_{r=1}^{R}\frac{B_{2r}}{(2r)!}\left(\frac{d}{dN}\right)^{2r-1}\epsilon_{N},$
(17)

where the $B_{2m}$ are the Bernoulli numbers and $E(\infty)$ is the carefully
regularized sum over all eigenvalues (here $E(\infty)=0$). We insert the WKB
expansion for $\epsilon_{j}$, Eq. (5), into Eq. (17). Equation (2.6) in Ref.
[13] explains how to evaluate the integral in Eq. (17):

$-\int_{N}^{\infty}dz\,\epsilon^{\text{WKB}}_{M}(z)=\alpha\sum_{m=0}^{M}\frac{b_{m}}{7/3-2m}Z^{7/3-2m},$
(18)

where $Z=N+1/2$. Noting that

$\left(\frac{d}{dN}\right)^{q}Z^{p+q}=Z^{p}(p+1)_{q},$ (19)

where $(x)_{n}=\Gamma(x+n)/\Gamma(x)$ is the Pochhammer symbol (Sec. 5.2 of
Ref. [24]), we find

$E_{M}(N)=\alpha\sum_{m=0}^{M}\left(\frac{3}{7}q_{m}Z-\frac{b_{m}}{2}\right)Z^{4/3-2m},$
(20)

where

$q_{n}=\frac{7/3}{7/3-2n}\sum_{m=0}^{n}B_{2(n-m)}\binom{4/3-2m}{2(n-m)}b_{m},$
(21)

where we use the binomial coefficient: $\binom{x}{r}=x...(x-r+1)/r!$. The
structure of Eq. (20) is analyzed in Appendix C. Expanding in powers of
$N^{-1}$ yields

$E_{M}(N)=\frac{3}{7}\alpha N^{7/3}C_{M}(N^{2}),$ (22)

where

$c_{n}=\frac{7}{7-6n}\sum_{m=0}^{n}\binom{2n-7/3}{2m}K_{m}b_{n-m},$ (23)

with $K_{0}=1$ and

$K_{m}=\frac{1-4^{1-m}m}{2}+\sum_{r=1}^{m}\left(\frac{1}{4^{r}}-\frac{1}{1+2r}\right)\binom{2m}{2r}B_{2(m-r)},$
(24)

for $m>0$. We give the $c_{n}$ to 41 decimal places in Table S14. The errors
of Eq. (22) appear in Fig. 5, showing missing subdominant corrections.

$\begin{array}[]{|c|c|r|c|r|r|r|r|}\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Error}\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr
N&M&\lx@intercol\hfil\text{Exact}\hfil\lx@intercol\vrule\lx@intercol&\Sigma\text{WKB}&\lx@intercol\hfil\text{SWKB}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{CSWKB}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Asym.
CSWKB}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{HYP}\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr 1&3&0.53018&-3.52704\times 10^{-2}&-3.72955\times
10^{-03}&3.04929\times 10^{-04}&3.04280\times 10^{-04}&-1.68851\times
10^{-04}\\\ 2&5&2.43002&-3.09659\times 10^{-2}&2.13376\times
10^{-04}&-3.91377\times 10^{-06}&-3.91161\times 10^{-06}&4.31969\times
10^{-07}\\\ 3&7&6.15787&-3.11941\times 10^{-2}&-1.06904\times
10^{-05}&4.58054\times 10^{-08}&4.57841\times 10^{-08}&-1.40070\times
10^{-09}\\\ 4&9&11.98024&-3.11829\times 10^{-2}&5.10050\times
10^{-07}&-5.40474\times 10^{-10}&-5.40187\times 10^{-10}&-1.09301\times
10^{-11}\\\ 5&11&20.11115&-3.11834\times 10^{-2}&-2.37643\times
10^{-08}&6.38318\times 10^{-12}&6.37835\times 10^{-12}&3.15656\times
10^{-13}\\\ 6&13&30.73034&-3.11834\times 10^{-2}&1.09166\times
10^{-09}&-7.54597\times 10^{-14}&-7.53647\times 10^{-14}&-5.18497\times
10^{-15}\\\ 7&15&43.99457&-3.11834\times 10^{-2}&-4.96713\times
10^{-11}&8.93476\times 10^{-16}&8.91385\times 10^{-16}&7.31437\times
10^{-17}\\\ 8&19&60.04387&-3.11834\times 10^{-2}&2.24450\times
10^{-12}&8.58161\times 10^{-18}&8.63167\times 10^{-18}&-8.89917\times
10^{-19}\\\ 9&21&79.00537&-3.11834\times 10^{-2}&-1.00888\times
10^{-13}&-9.97495\times 10^{-20}&-1.01030\times 10^{-19}&6.14869\times
10^{-21}\\\ 10&23&100.99595&-3.11834\times 10^{-2}&4.51598\times
10^{-15}&1.16502\times 10^{-21}&1.19954\times 10^{-21}&-2.03952\times
10^{-23}\\\ \hline\cr\end{array}$

Table 3: Errors in various approximations for sums: $\Sigma$WKB (sum of WKB
eigenvalues from Table 1), SWKB from Eq. (22), CSWKB from Eq. (28) with
subdominant corrections of Eq. (25), Asym. CSWKB using Eq. (27), and the
hyperasymptotic (HYP) approximation from Eq. (29).

$\begin{array}[]{|c|r|r|}\hline\cr
N&\lx@intercol\hfil\text{SD}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Asym.
SD}\hfil\lx@intercol\vrule\lx@intercol\\\ \hline\cr 1&4.03447966131\times
10^{-03}&4.03383052423\times 10^{-03}\\\ 2&-2.17289733012\times
10^{-04}&-2.17287572833\times 10^{-04}\\\ 3&1.07361905207\times
10^{-05}&1.07361692575\times 10^{-05}\\\ 4&-5.10590048004\times
10^{-07}&-5.10589761517\times 10^{-07}\\\ 5&2.37706662505\times
10^{-08}&2.37706614235\times 10^{-08}\\\ 6&-1.09173453706\times
10^{-09}&-1.09173444207\times 10^{-09}\\\ 7&4.96721587271\times
10^{-11}&4.96721566364\times 10^{-11}\\\ 8&-2.24449297891\times
10^{-12}&-2.24449292884\times 10^{-12}\\\ 9&1.00887497460\times
10^{-13}&1.00887496180\times 10^{-13}\\\ 10&-4.51598363452\times
10^{-15}&-4.51598360000\times 10^{-15}\\\ \hline\cr\end{array}$

Table 4: Subdominant corrections to sums, $E^{\rm SD}$, using Eq. (25) and its
asymptotic approximation, Eq. (27). Figure 6: Errors in approximations of
sums: SWKB (black), CSWKB (purple), and HYP (red), using Table 3, conventions
as in Fig. 1.

Table 3 shows results for various approximations to the sum of eigenvalues. We
plot the errors in Table 3 in Fig. 6. The performance of SWKB for the sum is
noticeably better than that of WKB for the eigenvalues. Since
$E(1)=\epsilon_{0}$, they are directly comparable for the lowest eigenvalue,
and here SWKB is an order of magnitude better than WKB, because the optimal
truncation occurs at $M$ = 3 for SWKB but at $M$ = 1 for WKB. More remarkably,
when subdominant corrections are included, the hyperasymptotic (HYP)
approximation has an error about 100 times smaller than CWKB. Thus, because it
increases the optimal truncation order, the summation method produces much
better results in even the most difficult case. We can also compare the SWKB
error for $N=10$ with the error of the WKB eigenvalue for $j=9$: The absolute
error is about 20 times smaller (and the relative error is even better).
Finally, we mention that the sum of WKB eigenvalues can never compete with the
asymptotic sums, as it is always dominated by its worst case, the ground-state
value, as shown by the fourth column of the table.

Figure 7: Errors and additions of approximate sums for $N=7$: SWKB (blue),
CSWKB (red), even additions(black), odd additions (grey), with conventions as
in Fig. 1.

Finally, we turn to the subdominant contributions to the asymptotics of the
sums. As $E(\infty)=0$, so $E^{\rm SD}(\infty)=0$, yielding

$E^{\rm SD}(N)=-\sum_{j=N}^{\infty}\epsilon^{\rm SD}_{L}(j+1/2),$ (25)

where we use Eq. (14) for $\epsilon^{\rm SD}_{L}(z)$. The results are shown
for $N=7$ in Fig. 7. The corrections greatly improve over the SWKB results, as
in the case of the eigenvalues (see Fig. 4). Here the even SWKB additions
closely follow the CSWKB error [Eq. (28)], while the odd additions are two
orders of magnitude smaller. We report the exact SD values and their
asymptotic approximations (Eqs. (25) & (27) respectively), used in Table 3, in
Table 4. We report the asymptotic values to more digits in Table S20.

Our last step is to turn the SD corrections to the sum into an asymptotic
series. We allow the order to vary in Eq. (25) and rewrite it as

$E^{\rm SD}_{M}(N)=-\epsilon_{M}^{\rm
SD}(Z)\sum_{j=0}^{\infty}\frac{(-1)^{j}}{e^{\pi
j}}\left(1+\frac{j}{Z}\right)^{1/3}\frac{H_{M}(j+Z)}{H_{M}(Z)}.$ (26)

Expanding the above summand in powers of $Z$ and summing each term from $j=0$
to infinity (in this case we can evaluate the sums analytically but in a more
general case we can still generate an asymptotic series using the Euler-
Maclaurin formula, see Refs. [12, 39] and Sec. 24.17 of Ref. [24]) yields

$E^{\rm SD}(N)=-\frac{\epsilon^{\rm
SD}(N)}{1+e^{-\pi}}\sum_{n=0}^{4}\frac{f_{n}}{Z^{n}},$ (27)

where we found that we only need the $f_{n}$ up to $n=4$ to recover the same
accuracy as Eq. (25). As always $f_{0}=1$ and we give the other $f_{n}$ in Eq.
(48). Even for $j=0$, these are always decreasing in magnitude, so we have not
reached the least addition.

We can now define a corrected SWKB (CSWKB) approximation:

$E^{\text{CSWKB}}_{M}(N)=E^{\rm SWKB}_{M}(N)+E^{\rm SD}(N),$ (28)

where we can use either Eq. (25) or (27) for $E^{\rm SD}(N)$. Table 3 shows
that both choices yield comparable accuracy. Since the even additions
determine the shape of the error, they determine the best truncation order. We
find the best results when we evaluate the SWKB contribution at the odd order
immediately before the even least addition ( given in Table 3). This is true
from $N=8$ onward, but we do not have enough $c_{n}$ coefficients to see how
far this trend goes. Table S22 shows that the errors for orders near the least
addition are similar.

Since the even additions are much larger in magnitude than the odd additions,
we ignore the odd additions and note that the remaining even additions
alternate in sign. We can therefore apply hyperasymptotics in exactly the same
fashion as for the linear half-well [13], adding half of the optimal SWKB term
plus half of the term two orders larger. Then we add this quantity to the
asymptotic approximation to $E^{\rm SD}(N)$ in Eq. (27). The end result is

$E^{\text{HYP}}_{M}(N)=\frac{E^{\rm SWKB}_{M}(N)+E^{\rm
SWKB}_{M+2}(N)}{2}+E^{\rm SD}(N).$ (29)

This hyperasymptotic approximation improves results by up to two orders of
magnitude over CSWKB, as shown in Table 3. Table S22 shows the hyperasymptotic
energies are more sensitive to changes in order than the CSWKB energies.

$\begin{array}[]{|c|r|r|r|r|r|r|}\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Error}\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr
j&\lx@intercol\hfil\text{Exact}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{WKB}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{SWKB}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{CWKB}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{CSWKB}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{HYP}\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr 0&0.53018&-3.52704\times 10^{-02}&-3.72955\times
10^{-03}&1.45719\times 10^{-02}&3.04280\times 10^{-04}&-1.68851\times
10^{-04}\\\ 1&1.89984&4.30449\times 10^{-03}&3.94293\times
10^{-03}&5.27248\times 10^{-05}&-3.08191\times 10^{-04}&1.69283\times
10^{-04}\\\ 2&3.72785&-2.28226\times 10^{-04}&-2.24066\times
10^{-04}&-1.99791\times 10^{-07}&3.95739\times 10^{-06}&-4.33369\times
10^{-07}\\\ 3&5.82237&1.12481\times 10^{-05}&1.12004\times
10^{-05}&1.29141\times 10^{-09}&-4.63243\times 10^{-08}&1.38977\times
10^{-09}\\\ 4&8.13091&-5.34365\times 10^{-07}&-5.33814\times
10^{-07}&-4.13332\times 10^{-12}&5.46566\times 10^{-10}&1.12457\times
10^{-11}\\\ 5&10.61919&2.48624\times 10^{-08}&2.48559\times
10^{-08}&-3.38058\times 10^{-14}&-6.45372\times 10^{-12}&-3.20841\times
10^{-13}\\\ 6&13.26424&-1.14141\times 10^{-09}&-1.14133\times
10^{-09}&1.06888\times 10^{-15}&7.62561\times 10^{-14}&5.25811\times
10^{-15}\\\ 7&16.04930&5.19166\times 10^{-11}&5.19158\times
10^{-11}&-1.86365\times 10^{-17}&-8.82753\times 10^{-16}&-7.40337\times
10^{-17}\\\ 8&18.96150&-2.34538\times 10^{-12}&-2.34539\times
10^{-12}&2.75295\times 10^{-19}&-8.73270\times 10^{-18}&8.96066\times
10^{-19}\\\ 9&21.99058&1.05404\times 10^{-13}&1.05404\times
10^{-13}&7.65097\times 10^{-20}&1.02230\times 10^{-19}&-6.16909\times
10^{-21}\\\ \hline\cr\end{array}$

Table 5: Errors in approximate eigenvalues, including differences of
approximate sums: WKB and CWKB from Table 1, SWKB, CSWKB, and hyperasymptotic
CSWKB of Table 3 plugged into Eq. (30), with the asymptotic subdominant
corrections of Eq. (27). Figure 8: Errors in individual eigenvalues, including
differences of approximate sums: WKB (black), SWKB (blue), CWKB (purple),
CSWKB (orange), and HYP (red), using Table 5 and the conventions of Fig. 1.

Next we use the various approximations for sums developed here to find new
approximations for the eigenvalues, by taking differences of sums:

$\epsilon_{j}=E(j+1)-E(j).$ (30)

We compare our sum approximations with the eigenvalue approximations of Sec. 3
in Table 5 and plot the errors in Fig. 8. As mentioned above, the SWKB result
is an order of magnitude better for $j=0$, but otherwise their errors are
almost identical. However, when subdominant corrections are included, the CWKB
results generally cannot be beaten, even with hyperasymptotic corrections to
the sums. Equation (17) also yields a simple approximation for
$E(N)+\epsilon_{N}/2$. In Appendix C, we check that such midway sums do not
yield sums of higher accuracy than those presented here.

##  5 Conclusion and Relevance to DFT

In the first half of this paper, we analyzed the semiclassical expansion of
the quartic oscillator eigenvalues. Although the formal structure of this
expansion has been known since the late 1970s, we have provided a thorough
analysis of the performance of the semiclassical expansion and reported the
coefficients of the relevant asymptotic series to many orders and many digits.
Moreover, in Eqs. (33), (36), and (38), we give the explicit formulas to
construct the implicit and explicit WKB series to any order. This should prove
useful for chemists, mathematicians, or physicists developing methods in the
many different contexts in which the quartic oscillator is studied.

But our primary results occur in the second half, where we report the first
application of the asymptotics of sums to a semiclassical expansion with
subdominant contributions (beyond the WKB series), due to the presence of
complex turning points. The method of asymptotics of sums, first demonstrated
in Ref. [13], can be easily adapted (at least in this case) to produce
subdominant corrections to the sums of eigenvalues. When these are included,
the performance for sums is better than the original series for the individual
eigenvalues, often by more than an order of magnitude. We can employ
hyperasymptotics to improve performance even further, yielding errors less
than $10^{-22}$ for the sum of the first 10 terms. Unsurprisingly, taking
differences of the asymptotic sums does not in general yield better results
than the original series, except for the ground state.

Local density functionals become relatively exact in a well-defined
semiclassical limit [19], but the leading corrections to these local
functionals in this limit are unknown. Modern GGAs are approximations to these
corrections. We do not have a general procedure for deriving these corrections
for the unknown XC energy of Kohn-Sham DFT, or even for the non-interacting
kinetic energy in three dimensions of orbital-free DFT. But the asymptotics of
sums of eigenvalues is directly related to orbital-free DFT in one dimension
because one considers the levels as filled by same-spin fermions. Using the
spin-scaling relation for the kinetic energy [40], and going from potential
functionals to density functionals [16, 17], one can relate the asymptotic
expansion for the sum of eigenvalues to the expansion of the non-interacting
kinetic energy around the semiclassical limit. Thus, the zero-order
contribution to that sum is a simple local functional of the potential. In the
full semiclassical limit, this result is identical to that of non-interacting
Thomas-Fermi theory, i.e., the result given by a local density approximation
to the non-interacting kinetic energy. Likewise, the leading asymptotic
correction to the sum yields the leading correction to the local
approximation, and is determined both by the derivatives of the potential (up
to the second) and by the Maslov index for the problem. Our results here show
that, in principle, far more accurate results can be extracted from these
semiclassical expansions than is currently possible with density functional
approximations. Inclusion of even the exact leading correction ($M=1$ here)
can greatly improve results [11].

The subdominant corrections appearing in the current work are largest for
$N=1$ (one particle in the lowest eigenvalue), where they would be non-
negligible in an electronic structure calculation. The present work shows that
such terms can be handled in terms of asymptotics of sums. Of course, this
would also require the further step of finding a general expression for their
form and coefficients, not just the value for a quartic oscillator, but that
step is beyond the scope of the current work.

##  6 Acknowledgments

We thank the NSF (CHE-1856165) for funding and Micheal Berry for useful
discussions.

## References

  * Hartree and Hartree [1935] D. R. Hartree and W. Hartree, Self-consistent field, with exchange, for beryllium, Proceedings of the Royal Society of London. Series A - Mathematical and Physical Sciences 150, 9 (1935).
  * Szabo and Ostlund [1996] A. Szabo and N. S. Ostlund,  _Modern quantum chemistry_ (Dover Publications, 1996).
  * March and Plaskett [1956] N. H. March and J. S. Plaskett, The relation between the wentzel-kramers-brillouin and the thomas-fermi approximations, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 235, 419 (1956).
  * Englert [1988] B.-G. Englert,  _Semiclassical Theory of Atoms_ (Springer, Verlag Berlin Heidelberg, 1988).
  * Perdew  _et al._ [2006] J. P. Perdew, L. A. Constantin, E. Sagvolden, and K. Burke, Relevance of the slowly varying electron gas to atoms, molecules, and solids, Phys. Rev. Lett. 97, 223002 (2006).
  * Elliott and Burke [2009] P. Elliott and K. Burke, Non-empirical derivation of the parameter in the b88 exchange functional, Can. J. Chem. Ecol. 87, 1485 (2009).
  * Burke  _et al._ [2016] K. Burke, A. Cancio, T. Gould, and S. Pittalis, Locality of correlation in density functional theory, The Journal of Chemical Physics 145, 054112 (2016).
  * Cancio  _et al._ [2018] A. Cancio, G. P. Chen, B. T. Krull, and K. Burke, Fitting a round peg into a round hole: Asymptotically correcting the generalized gradient approximation for correlation, The Journal of Chemical Physics 149, 084116 (2018).
  * Sun  _et al._ [2015] J. Sun, A. Ruzsinszky, and J. P. Perdew, Strongly constrained and appropriately normed semilocal density functional, Phys. Rev. Lett. 115, 036402 (2015).
  * Perdew  _et al._ [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Restoring the density-gradient expansion for exchange in solids and surfaces, Phys. Rev. Lett. 100, 136406 (2008).
  * Burke [2020a] K. Burke, Leading correction to the local density approximation of the kinetic energy in one dimension, The Journal of Chemical Physics 152, 081102 (2020a).
  * Burke [2020b] K. Burke, Deriving approximate density functionals with asymptotics, Faraday Discuss. 224, 98 (2020b).
  * Berry and Burke [2020] M. V. Berry and K. Burke, Exact and approximate energy sums in potential wells, Journal of Physics A: Mathematical and Theoretical 53, 095203 (2020).
  * Boyd [1999] J. P. Boyd, The devil’s invention: Asymptotic, superasymptotic and hyperasymptotic series, Acta Applicandae Mathematica 56, 1 (1999).
  * Berry [1991] M. Berry, Asymptotics, superasymptotics, hyperasymptotics…, in  _Asymptotics beyond All Orders_ , edited by H. Segur, S. Tanveer, and H. Levine (Springer, Boston, 1991) pp. 1–14.
  * Cangi  _et al._ [2011] A. Cangi, D. Lee, P. Elliott, K. Burke, and E. K. U. Gross, Electronic structure via potential functional approximations, Phys. Rev. Lett. 106, 236404 (2011).
  * Cangi  _et al._ [2013] A. Cangi, E. K. U. Gross, and K. Burke, Potential functionals versus density functionals, Phys. Rev. A 88, 062505 (2013).
  * Lignères and Carter [2005] V. L. Lignères and E. A. Carter, An introduction to orbital-free density functional theory, in  _Handbook of Materials Modeling: Methods_ , edited by S. Yip (Springer Netherlands, Dordrecht, 2005) pp. 137–148.
  * Okun and Burke [2021a] P. Okun and K. Burke, Semiclassics: The hidden theory behind the success of dft (2021a), arXiv:2105.04384.
  * Bender and Orszag [1999] C. M. Bender and S. A. Orszag,  _Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory_ (Springer, Verlag New York, 1999).
  * Dutt  _et al._ [1993] R. Dutt, A. Gangopadhyaya, A. Khare, A. Pagnamenta, and U. Sukhatme, Semiclassical approach to quantum-mechanical problems with broken supersymmetry, Phys. Rev. A 48, 1845 (1993).
  * Cooper  _et al._ [1995] F. Cooper, A. Khare, and U. Sukhatme, Supersymmetry and quantum mechanics, Physics Reports 251, 267 (1995).
  * Aniceto  _et al._ [2019] I. Aniceto, G. Basar, and R. Schiappa, A primer on resurgent transseries and their asymptotics, Physics Reports 809, 1 (2019), a primer on resurgent transseries and their asymptotics.
  * [24] DLMF, NIST Digital Library of Mathematical Functions, http://dlmf.nist.gov/, Release 1.1.3 of 2021-09-15, f. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
  * Lay [1997] W. Lay, The quartic oscillator, Journal of Mathematical Physics 38, 639 (1997), https://doi.org/10.1063/1.531857 .
  * Bay and Lay [1997] K. Bay and W. Lay, The spectrum of the quartic oscillator, Journal of Mathematical Physics 38, 2127 (1997).
  * Liverts  _et al._ [2006] E. Z. Liverts, V. B. Mandelzweig, and F. Tabakin, Analytic calculation of energies and wave functions of the quartic and pure quartic oscillators, Journal of Mathematical Physics 47, 062109 (2006).
  * Balian  _et al._ [1979] R. Balian, G. Parisi, and A. Voros, Quartic oscillator, in  _Feynman Path Integrals_ , edited by S. Albeverio, P. Combe, R. Høegh-Krohn, G. Rideau, M. Sirugue-Collin, M. Sirugue, and R. Stora (Springer Berlin Heidelberg, Berlin, Heidelberg, 1979) pp. 337–360.
  * Voros [1983] A. Voros, The return of the quartic oscillator. the complex wkb method, Annales de l’I.H.P. Physique théorique 39, 211 (1983).
  * Balian  _et al._ [1978] R. Balian, G. Parisi, and A. Voros, Discrepancies from asymptotic series and their relation to complex classical trajectories, Phys. Rev. Lett. 41, 1141 (1978).
  * Blinder [2019] S. M. Blinder, Eigenvalues for a pure quartic oscillator (2019), arXiv:1903.07471.
  * Reid [1970] C. E. Reid, Energy eigenvalues and matrix elements for the quartic oscillator, Journal of Molecular Spectroscopy 36, 183 (1970).
  * Voros [1980] A. Voros, The zeta function of the quartic oscillator, Nuclear Physics B 165, 209 (1980).
  * Delabaere and Pham [1997] E. Delabaere and F. Pham, Unfolding the quartic oscillator, Annals of Physics 261, 180 (1997).
  * White and Kutlin [2017] R. B. White and A. G. Kutlin, Bound state energies using phase integral analysis (2017), arXiv:1704.01170.
  * Dunne and Ünsal [2014] G. V. Dunne and M. Ünsal, Uniform wkb, multi-instantons, and resurgent trans-series, Phys. Rev. D 89, 105009 (2014).
  * Dunne and Ünsal [2017] G. V. Dunne and M. Ünsal, Wkb and resurgence in the mathieu equation, in  _Resurgence, Physics and Numbers_ , edited by F. Fauvet, D. Manchon, S. Marmi, and D. Sauzin (Scuola Normale Superiore, Pisa, 2017) pp. 249–298.
  * Delabaere  _et al._ [1997] E. Delabaere, H. Dillinger, and F. Pham, Exact semiclassical expansions for one-dimensional quantum oscillators, Journal of Mathematical Physics 38, 6126 (1997).
  * Hua [2012] L.-K. Hua,  _Introduction to Higher Mathematics_ (Cambridge University Press, Cambridge, UK, 2012).
  * Burke [2007] K. Burke, The abc of dft (2007), https://dft.uci.edu/doc/g1.pdf.
  * Okun and Burke [2021b] P. Okun and K. Burke, Uncommonly accurate energies for the general quartic oscillator, International Journal of Quantum Chemistry 121, e26554 (2021b).

##  Appendix A Deriving the Implicit WKB Series

Following Bender & Orszag [20] we write the WKB wavefunction as
$\psi(x)=\exp[iS(x)]$ and expand the action in a semiclassical series
$S(x)=S_{0}(x)+S_{1}(x)+\cdots$. For the pure quartic oscillator

$S_{n}^{\prime}(x)=(-1)^{{\lfloor}(n+1)/2{\rfloor}}i^{n+1}\sum_{m=0}^{{\lfloor}3n/4{\rfloor}}\frac{k_{n}^{(m)}x^{3n-4m}}{[p(x)]^{3n-2m-1}},$
(31)

where $p(x)=\sqrt{2\epsilon-x^{4}}$ is the classical momentum and the prime
indicates differentiation, and $k_{0}^{(0)}=k_{1}^{(0)}=1$ [20]. The higher
order ($n>1$) action coefficients are determined by

$2S_{0}^{\prime}S_{n}^{\prime}+S_{n-1}^{\prime\prime}+\sum_{m=1}^{n-1}S_{m}^{\prime}S_{n-m}^{\prime}=0.$
(32)

This recursion relation, derived in Ref. [20], yields

$k_{n}^{(m)}=(-1)^{n}\bigg{[}(3n-2m-4)k_{n-1}^{(m)}+\\\
\frac{3n-4m+1}{2}k_{n-1}^{(m-1)}\bigg{]}-\frac{1}{2}\sum_{l=1}^{n-1}(-1)^{l(n+1)}\sum_{r=0}^{m}k_{l}^{(r)}k_{n-l}^{(m-r)},\\\
$ (33)

for $n>1$. To use Eq. (33), we set all coefficients $k_{n}^{(m)}$ with
$m>{\lfloor}3n/4{\rfloor}$ and $m<0$ to 0. In general each $k_{n}^{(m)}$ is a
function of the set $\\{k_{n^{\prime}}^{(m^{\prime})}\\}$ with $m^{\prime}\leq
m$ and $n^{\prime}<n$. We give the $k_{n}^{(m)}$ up to $n=21$ in Table S4. The
eigenvalues are determined only by the $S_{n}(x)$ with even $n$ so we drop the
odd $n$ values and integrate Eq. (31):

$S_{2n}(\epsilon)=2i\sum_{m=0}^{{\lfloor}3n/2{\rfloor}}k_{2n}^{(m)}\int_{0}^{x_{0}}dx\,\frac{x^{6n-4m}}{[p(x)]^{6n-2m-1}},$
(34)

where $x_{0}=x_{0}(\epsilon)$ is the classical turning point defined by
$p[x_{0}(\epsilon)]=0$. Formally there are divergences at $x_{0}$ but we
remove them by integrating under the integral sign with respect to $\epsilon$
[41]. We know from Ref. [20] that

$S_{2n}(\epsilon)=\frac{i\,a_{n}\gamma^{2}(2\epsilon)^{3/4-3n/2}}{3\sqrt{2\pi}},$
(35)

so after some algebra we get

$a_{n}=\frac{3}{\gamma^{2}}\sqrt{\frac{\pi}{2}}\sum_{m=0}^{{\lfloor}3n/2{\rfloor}}k_{2n}^{(m)}B\left(m-3n+\frac{3}{2},\frac{3n}{2}-m+\frac{1}{4}\right),$
(36)

where $B(x,y)$ is the incomplete Euler beta function
$\Gamma(x)\Gamma(y)/\Gamma(x+y)$ (Sec. 8.17 of Ref. [24]). Equation (36) was
found by Voros, Balian, and Parisi in Ref. [28]. The first three nontrivial
($n=1,2,3$) $a_{n}$ are

$-\frac{3\pi^{2}}{2\gamma^{4}},\qquad\frac{11}{512},\qquad\frac{4697\pi^{2}}{5120\gamma^{4}}.$
(37)

##  Appendix B Analytic Expressions for Coefficients

We generate the $b_{n}$ in Eq. (5) using the $a_{n}$ and lower order $b_{n}$:

$b_{n}=\sum_{m=0}^{n}a_{m}\sum_{r=0}^{n-m}G_{nmr}\mathscr{P}_{r}^{(n-m)}(\mathbf{b}),\qquad
n\geq 1,$ (38)

where $\mathscr{P}_{0}^{(n)}=\delta_{n,0}$ and

$\mathscr{P}_{m}^{(n)}(\mathbf{b})=\sum_{\begin{subarray}{c}k_{1},k_{2},...,k_{m}=1\\\
k_{1}+k_{2}+...+k_{m}=n\end{subarray}}^{n+1-m}b_{k_{1}}b_{k_{2}}...b_{k_{m}},$
(39)

e.g. $\mathscr{P}_{3}^{(5)}(\mathbf{b})=3b_{1}b_{2}^{2}+3b_{1}^{2}b_{3}$, and
$G_{n01}=0$ and

$G_{nmr}=-\frac{4}{3(2\alpha)^{3m/2}}\binom{3/4-3m/2}{r}.$ (40)

The coefficients we have derived in this paper ($b_{n}$, $c_{n}$, $d_{n}$,
$h_{n}$, $k_{n}$, $f_{n}$) can each be expressed in terms of the $a_{n}$
coefficients. In general each set of coefficients is composed of more
fundamental coefficients. Here we give both the relation to simpler
coefficients and explicit exact formulas. The supplemental info contains all
tables needed to reconstruct these, as well as numerical values to 40 digits.
We remind the reader that our convention is to set the $n=0$ coefficient to 1.

Eliminating lower order $b_{n}$’s from Eq. (38) yields the $b_{n}$ in terms of
the $a_{n}$ alone:

$b_{n}=-\sum_{m=0}^{n-1}L_{nm}\mathscr{P}^{(n)}_{n-m}(\mathbf{a}),~{}~{}n\geq
1,$ (41)

where we give the $L_{nm}$ in Table S6. Inserting the $a_{n}$ from Eq. (36)
yields exact values for the $b_{n}$ explicitly given in Table S7. The first 3
nontrivial $b_{n}$ are

$\frac{1}{9\pi},~{}~{}-\frac{1}{648\pi^{2}}\left(5+\frac{11\gamma^{8}}{192\pi^{4}}\right),~{}~{}\frac{11}{2916\pi^{3}}\left(\frac{1}{3}-\frac{31\gamma^{8}}{640\pi^{4}}\right).$
(42)

Similarly, we can express the $d_{n}$ in terms of the $a_{n}$ and $b_{n}$:

$d_{n}=\sum_{m=0}^{n}a_{n-m}\sum_{r=0}^{m}F_{nmr}\mathscr{P}^{(m)}_{r}(\mathbf{b}),$
(43)

where we give the $F_{nmr}$ in Table S9. Inserting the $a_{n}$ from Eq. (36)
and the $b_{n}$ from Eq. (41) yields explicit expressions for the $d_{n}$
which we give in Table S10. The first 3 nontrivial $d_{n}$ are

$\frac{1}{6\pi},\quad-\frac{1}{72\pi^{2}},\qquad\frac{1}{432\pi^{3}}\left(1-\frac{1133\gamma^{8}}{8640\pi^{4}}\right).$
(44)

For our next set, it is best to represent the $h_{n}$ directly in terms of the
$a_{n}$:

$h_{n}=\sum_{m=0}^{{\lfloor}n/2{\rfloor}}P_{nm}(\mathbf{a}),$ (45)

where the $P_{nm}(\mathbf{a})$ are polynomials of $a_{n}$, given in Table S12,
for which we have not been able to find a general pattern. We give the $h_{n}$
themselves in Table S13. The first 3 nontrivial $h_{n}$ are

$-\frac{1}{6},\qquad\frac{1}{72}\left(1-\frac{4}{\pi}\right),\qquad\frac{1}{216}\left(\frac{5}{\pi}-\frac{1}{6}\right).$
(46)

We give the $c_{n}$ in terms of the $b_{n}$ in Eqs. (23) & (24). We give the
exact forms of the $c_{n}$ up to $n=10$ in Table S15. The first 3 nontrivial
$c_{n}$ are

$\frac{7}{9\pi}\left(1-\frac{\pi}{6}\right),\qquad\frac{7}{324}\left(\frac{1}{3\pi}+\frac{1}{2\pi^{2}}-\frac{7}{180}+\frac{11\gamma^{8}}{1920\pi^{6}}\right),\\\
\frac{7}{2916}\left[\frac{31}{756}-\frac{7}{18\pi}-\frac{5}{6\pi^{2}}-\frac{1}{3\pi^{3}}-\frac{\gamma^{8}}{128\pi^{6}}\left(\frac{11}{9}-\frac{31}{5\pi}\right)\right].\\\
$ (47)

We give the $f_{n}$ in terms of the $h_{n}$ using $f_{n}=k_{n}/(1+e^{\pi})$
where

$\displaystyle\begin{split}k_{1}&=-1/3,\\\
k_{2}&=\frac{\text{tanh}(\pi/2)}{9}+h_{1},\\\
k_{3}&=-\frac{5(1-4e^{\pi}+e^{2\pi})}{81(1+e^{\pi})^{2}}-\frac{2h_{1}}{3}\text{tanh}\left(\frac{\pi}{2}\right)-h_{1}^{2}+2h_{2},\\\
k_{4}&=\frac{20e^{\pi}(\cosh\pi-5)}{243(1+e^{\pi})^{2}}\text{tanh}\left(\frac{\pi}{2}\right)+\frac{5(1-4e^{\pi}+e^{2\pi})h_{1}}{9(1+e^{\pi})^{2}}\\\
&+\frac{2h_{1}^{2}-7h_{2}}{3}\text{tanh}\left(\frac{\pi}{2}\right)+h_{1}^{3}+3(h_{3}-h_{1}h_{2}).\\\
\end{split}$ (48)

And the first three nontrivial $h_{n}$ are given exactly in Eq. (46).

##  Appendix C Midway sums

$\begin{array}[]{|c|r|c|c|c|r|r|r|}\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Order}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Error}\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr
N&\lx@intercol\hfil\text{Exact}\hfil\lx@intercol\vrule\lx@intercol&\text{SWKB}&\text{CSWKB}&\text{HYP}&\lx@intercol\hfil\text{SWKB}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{CSWKB}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{HYP}\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr 0&0.26509&3&4&2&3.71006\times 10^{-02}&-1.24217\times
10^{-03}&-2.01413\times 10^{-03}\\\ 1&1.48010&5&4&3&-2.03890\times
10^{-03}&5.66343\times 10^{-05}&-1.97554\times 10^{-05}\\\
2&4.29394&7&6&6&1.04352\times 10^{-04}&-8.70496\times 10^{-07}&3.55482\times
10^{-08}\\\ 3&9.06905&9&8&8&-5.12434\times 10^{-06}&1.05384\times
10^{-08}&-4.83320\times 10^{-11}\\\ 4&16.04570&11&12&10&2.43535\times
10^{-07}&1.17658\times 10^{-10}&-4.52413\times 10^{-12}\\\
5&25.42075&13&14&12&-1.13409\times 10^{-08}&-1.33613\times
10^{-12}&9.01092\times 10^{-14}\\\ 6&37.36246&15&16&14&5.21047\times
10^{-10}&1.53515\times 10^{-14}&-1.35614\times 10^{-15}\\\
7&52.01922&17&18&16&-2.37140\times 10^{-11}&-1.78113\times
10^{-16}&1.83652\times 10^{-17}\\\ 8&69.52462&21&20&20&1.07180\times
10^{-12}&2.08591\times 10^{-18}&-1.62795\times 10^{-19}\\\
9&90.00066&23&22&22&-4.81857\times 10^{-14}&-2.47297\times
10^{-20}&8.23128\times 10^{-22}\\\ 10&113.56002&25&26&24&2.15726\times
10^{-15}&-2.86261\times 10^{-22}&6.70277\times 10^{-24}\\\
\hline\cr\end{array}$

Table 6: The exact $G(N)$ and the errors of various approximations. We use
least odd order optimal truncation for SWKB and empirical optimal truncation
for the other approximations and report the orders.

An observant reader will have noticed that Eq. (20) yields a simpler
expression for the following midway sum:

$G(N)=E(N)+\frac{\epsilon_{N}}{2}.$ (49)

By rearranging Eq. (20) we find an SWKB approximation for $G(N)$:

$G^{\rm SWKB}_{M}(N)=\frac{3}{7}\alpha Z^{7/3}Q_{M}(Z^{2}).$ (50)

While this is not the sum of primary interest, we explore here if this might
yield more accurate approximations for $E(N)$. The first three nontrivial
$q_{n}$ are

$\frac{7}{9}\left(\frac{1}{3}+\frac{1}{\pi}\right),\qquad\frac{7}{162}\left(\frac{1}{45}-\frac{1}{3\pi}+\frac{1}{4\pi^{2}}\right)+\frac{77\gamma^{8}}{622080\pi^{6}},\\\
-\frac{1}{2187}\left(\frac{2}{9}-\frac{7}{3\pi}-\frac{35}{4\pi^{2}}+\frac{7}{4\pi^{3}}\right)+\frac{7\gamma^{8}}{186624\pi^{6}}\left(\frac{11}{9}+\frac{31}{10\pi}\right).\\\
$ (51)

We give the $q_{n}$ exactly up to $n=10$ in Table S17 and numerically in Table
S16. Equation (21) gives the $q_{n}$ in terms of the $b_{n}$.

It is easy to extract the SD contribution to $G$ from that of $E$:

$G^{\rm SD}(N)=\frac{E^{\rm SD}(N+1)+E^{\rm SD}(N)}{2},$ (52)

where we use Eq. (27) for $E^{\rm SD}$. Thus

$G^{\text{CSWKB}}_{M}(N)=G^{\rm SWKB}_{M}(N)+G^{\rm SD}(N).$ (53)

Figure 9: The even (gray) and odd (black) SWKB additions with the SWKB (blue)
and CSWKB (red) errors for $N=5$. Figure 10: The errors in Table 6: SWKB
(blue), CSWKB (red), and HYP (purple).

First, we consider how accurately we can approximate $G(N)$. In Fig. 9 we show
the SWKB additions and compare the SWKB and CSWKB errors for $N=5$. As always,
the SWKB errors plateau but the SD corrections increase the accuracy by many
orders of magnitude. Since the even SWKB additions are much smaller in
magnitude than the odd additions, a hyperasymptotic approximation to $G$ is

$G^{\text{HYP}}_{M}(N)=\frac{G^{\rm SWKB}_{M}(N)+G^{\rm
SWKB}_{M+2}(N)}{2}+G^{\rm SD}(N).$ (54)

We compare the errors of the SWKB, CSWKB, and hyperasymptotic approximations
to $G(N)$ in Table 6, while Fig. 10 shows the errors. Comparing Tables 3 and 6
shows that the CSWKB and hyperasymptotic approximations are more accurate for
$G(N)$ than for $E(N)$.

It would be nice if we could make use of this accuracy to turn the $G(N)$
approximation into an $E(N)$ approximation. We can do this by subtracting half
of the Nth eigenvalue from our $G(N)$ approximation. Table 5 shows that CWKB
is our most accurate eigenvalue approximation. We show in Table 7 that
combining this with our hyperasymptotic $G(N)$ approximation gives a more
accurate $E(N)$ approximation than the hyperasymptotic $E(N)$ (our best $E(N)$
approximation in Table 3). Thus we need to mix the $G(N)$ approximation with
an eigenvalue approximation to generate an $E(N)$ approximation. We prefer to
work with pure $E(N)$ approximations instead. Thus we relegate our discussion
of the $G(N)$ approximations to this appendix.

$\begin{array}[]{|c|r|c|r|r|}\hline\cr\vrule\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{Error}\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr
N&\lx@intercol\hfil\text{Exact}\hfil\lx@intercol\vrule\lx@intercol&M&\lx@intercol\hfil\text{HYP
}E\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil\text{HYP
}G\hfil\lx@intercol\vrule\lx@intercol\\\ \hline\cr 1&0.53018&3&-1.68851\times
10^{-04}&-4.61178\times 10^{-05}\\\ 2&2.43002&6&4.31969\times
10^{-07}&1.35444\times 10^{-07}\\\ 3&6.15787&7&-1.40070\times
10^{-09}&-4.73745\times 10^{-10}\\\ 4&11.98024&10&-1.09301\times
10^{-11}&-2.45748\times 10^{-12}\\\ 5&20.11115&12&3.15656\times
10^{-13}&1.07012\times 10^{-13}\\\ 6&30.73034&14&-5.18497\times
10^{-15}&-1.89058\times 10^{-15}\\\ 7&43.99457&16&7.31437\times
10^{-17}&2.76835\times 10^{-17}\\\ 8&60.04387&20&-8.89917\times
10^{-19}&-3.00443\times 10^{-19}\\\ 9&79.00537&16&6.14869\times
10^{-21}&1.47877\times 10^{-20}\\\ \hline\cr\end{array}$

Table 7: The hyperasymptotic $E(N)$ approximation of Table 3 compared to the
hyperasymptotic $G(N)$ approximation of Table 7 turned into an $E(N)$
approximation by subtracting away half of our best estimate for
$\epsilon_{N}$: $\epsilon^{\text{CWKB}}(N)$ in Eq. (15). We give the
hyperasymptotic $G(N)$ approximation in Eq. (54) and use empirical optimal
truncation to get the order $M$.

For $N=0$, $2G(0)$ is an estimate of the ground state eigenvalue. In Table 8
we test various approximations to the ground state against the various $G(N)$
approximations developed here.

$\begin{array}[]{|c|c|r|}\hline\cr\text{Approx.}&\text{Order}&\lx@intercol\hfil\text{Error}\hfil\lx@intercol\vrule\lx@intercol\\\
\hline\cr\text{WKB}&1&-0.0352704\\\ \text{SWKB}&3&-0.0037296\\\ \text{G
SWKB}&1&0.0323487\\\ \text{CWKB}&1&0.0145719\\\ \text{CSWKB}&3&0.0003043\\\
\text{G CSWKB}&4&-0.0024843\\\ \text{E HYP}&3&-0.0001689\\\ \text{G
HYP}&2&-0.0040283\\\ \hline\cr\end{array}$

Table 8: Errors in various approximations to $\epsilon_{0}=0.53018$. We take
the G versions of the SWKB, CSWKB, and HYP approximations at the order of
empirical optimal truncation. The other approximations come from Tables 1 and
3.

##  Appendix D Table of Constants

$\begin{array}[]{|c|r|r|r|r|r|r|}\hline\cr n&\lx@intercol\hfil
a_{n}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil
b_{n}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil
d_{n}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil
h_{n}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil
h_{n+20}\hfil\lx@intercol\vrule\lx@intercol&\lx@intercol\hfil
c_{n}\hfil\lx@intercol\vrule\lx@intercol\\\ \hline\cr 1&-8.567748394584\times
10^{-02}&3.536776513153\times 10^{-02}&5.305164769730\times
10^{-02}&-1.666666666667\times 10^{-1}&5.194304272632\times
10^{04}&1.179447262911\times 10^{-01}\\\ 2&2.148437500000\times
10^{-02}&-3.527579746032\times 10^{-03}&-1.407238661699\times
10^{-03}&-3.794993676877\times 10^{-3}&1.622194035444\times
10^{05}&6.390774747964\times 10^{-03}\\\ 3&5.239936746010\times
10^{-02}&-1.765731366249\times 10^{-03}&-2.926108161501\times
10^{-03}&6.596679464131\times 10^{-3}&-5.982812586891\times
10^{04}&1.036757028064\times 10^{-05}\\\ 4&-3.188528333391\times
10^{-01}&4.173466791769\times 10^{-03}&4.726357002028\times
10^{-04}&6.104871777506\times 10^{-3}&-4.628688644053\times
10^{06}&-1.972265116005\times 10^{-03}\\\ 5&-2.906253581665\times
10^{+00}&1.016471233994\times 10^{-02}&1.628857392426\times
10^{-02}&8.012575730396\times 10^{-3}&-3.434326216453\times
10^{07}&-4.866940504878\times 10^{-05}\\\ 6&4.149581824006\times
10^{+01}&-5.136197691366\times 10^{-02}&-4.352082661762\times
10^{-03}&4.764793921301\times 10^{-3}&-1.291129925755\times
10^{08}&1.509652803690\times 10^{-02}\\\ 7&9.106051111453\times
10^{+02}&-3.151839648774\times 10^{-01}&-4.935113066975\times
10^{-01}&-2.526532163520\times 10^{-3}&4.682761458754\times
10^{07}&-6.970540753456\times 10^{-05}\\\ 8&-2.759178233251\times
10^{+04}&3.221164312655\times 10^{+00}&1.837863615021\times
10^{-01}&-2.047834105176\times 10^{-2}&5.015568432482\times
10^{09}&-6.683710532285\times 10^{-01}\\\ 9&-1.086497889469\times
10^{+06}&3.653193762830\times 10^{+01}&5.661119950435\times
10^{+01}&-4.760940814614\times 10^{-2}&4.359779213987\times
10^{10}&-1.864942552345\times 10^{-04}\\\ 10&5.444198816876\times
10^{+07}&-6.046509882828\times 10^{+02}&-2.705683388100\times
10^{+01}&-5.781448172999\times 10^{-2}&1.916566221097\times
10^{11}&9.795335354117\times 10^{+01}\\\ 11&3.396513559809\times
10^{+09}&-1.103195365766\times 10^{+04}&-1.698417456463\times
10^{+04}&2.484442856411\times 10^{-2}&-6.868413739548\times
10^{10}&-1.398393112909\times 10^{+00}\\\ 12&-2.575567562301\times
10^{+11}&2.727413882117\times 10^{+05}&9.916174888450\times
10^{+03}&4.073552890690\times 10^{-1}&-9.740013068265\times
10^{12}&-3.613621803518\times 10^{+04}\\\ 13&-2.330666020371\times
10^{+13}&7.291976034200\times 10^{+06}&1.117814650219\times
10^{+07}&1.481124311328\times 10^{+0}&-9.705364249021\times
10^{13}&7.386689584732\times 10^{+02}\\\ 14&2.482592926162\times
10^{+15}&-2.510854389322\times 10^{+08}&-7.711376112560\times
10^{+06}&2.741314325143\times 10^{+0}&-4.884300202923\times
10^{14}&2.816398059630\times 10^{+07}\\\ 15&3.075586665318\times
10^{+17}&-9.255480436013\times 10^{+09}&-1.414428669791\times
10^{+10}&-1.081551161491\times 10^{+0}&1.734725923736\times
10^{14}&-8.442791957114\times 10^{+05}\\\ 16&-4.384018949045\times
10^{+19}&4.238678892387\times 10^{+11}&1.125769296301\times
10^{+10}&-3.528471843003\times 10^{+1}&3.145775486985\times
10^{16}&-4.121307364131\times 10^{+10}\\\ 17&-7.123931840874\times
10^{+21}&2.060168174148\times 10^{+13}&3.141070779556\times
10^{+13}&-1.719162719101\times 10^{+2}&3.534910367892\times
10^{17}&1.580982245534\times 10^{+09}\\\ 18&1.309156462631\times
10^{+24}&-1.210778641220\times 10^{+15}&-2.833222450361\times
10^{+13}&-4.274928378220\times 10^{+2}&2.003996635297\times
10^{18}&1.038956384366\times 10^{+14}\\\ 19&2.701473930086\times
10^{+26}&-7.503088536990\times 10^{+16}&-1.141906853721\times
10^{+17}&1.616892923721\times 10^{+2}&-7.067807271510\times
10^{17}&-4.901331150371\times 10^{+12}\\\ 20&-6.219943513829\times
10^{+28}&5.505039532830\times 10^{+18}&1.151129625648\times
10^{+17}&8.443442098119\times 10^{+3}&-1.595548697681\times
10^{20}&-4.227171913689\times 10^{+17}\\\ \hline\cr\end{array}$

Table 9: We present to 13 digits the first 21 examples from each set of
coefficients used in this paper (remember that the zeroth order coefficient
for each set is 1). We give more coefficients and more digits in the
supplemental info (Tables S2, S5, S8, S11, and S14).


##  0.1  The Enigma of the Mass.

V.G. Plekhanov

Computer Science College, Erika Street 7a, Tallinn, 10416, Estonia

Abstract.

The different manifestations of the mass effects in the microphysics (isotope
effect) are presented for the first time. The bright effects observe in all
branches of physics: nuclear, atomic, and molecular as well as solid state
physics. Charge symmetry breaking in the strong interaction occurs because of
the difference between the masses of the up and down quarks. At present the
Standard Model can’t explain the observed mass pattern (Mn, Mp, mu, md etc.)
and their hierarchy. The last one doesn’t permit us to find the origin of the
isotope effect. The origin of the mass of the matter will be clarified when
the mechanism of chiral symmetry breaking in QCD is established.

Mass is a one of the fundamental properties of matter. It relates to classical
as well as modern physics (quantum mechanics or modern theory of gravitation
(see, e.g. [1]). Although the physical meaning of mass was discovered by
Einstein more than a century ago, when he introduced in physics the concept of
rest energy (E0) [2], the concept of mass still doesn’t have strict
mathematical determination. Indeed, according to the notion of the
relativistical physics (see, e.g. [3]) mass is determined by the next
expression

m2 = $\frac{\text{E}^{2}}{\text{c}^{4}}$ \-
$\frac{\overrightarrow{\text{p}^{2}}}{\text{c}^{2}}$ (1).

And in the case of resting body ($\overrightarrow{\text{p}}$ = 0) we have

m = $\frac{\text{E}_{0}}{\text{c}^{2}}$ (2).

From equation (2) it can be seen that the mass is proportioned to the rest
energy. If we put c = 1, in that case we see that the mass of body equals its
rest energy. The mass of a body is not a constant, it varies with changes in
its energy. Namely, rest energy ”slumbering” in massive bodies partly is
released in chemical and especially nuclear reactions. In spite of equivalence
of the mass of the body and rest energy, especially nuclear physics and
physics of elementary particles, the task of mass has not been solved. Until
present time the spectrum of the discrete hierarchy of elementary particles
mass hasn’t had a successful theoretical explantion [4,6]. As is well-known on
the boundary of the 19 and 20 centuries there was an opinion that the mass of
the electron has the electromagnetic origin [1,4]. However, later
investigations showed that the electromagnetic part of the mass of the
electron has a small contribution to its full mass [3]. Nevertheless, the
modern view connects the origin of the mass with nonlocal gravitational
fields, which nature is due to electromagnetic interaction [8 -11].This
conclusion reflects those fact, that the space between separated particles in
essence isn’t empty, it is filled with the material medium - the physical
fields. The space inside the atom is filled with electromagnetic field, and
inside nucleus - more densier and stronger field which is called sometimes
meson one.

The present letter is devoted to the elucidation of the origin of mass, so far
as only its nature closely connected with the origin of the isotope effect,
the experimental manifestation of which more persuasively testified in the
last fifty years in all branches of physics (nuclear, atomic, molecular as
well as solid state (see, e.g. reviews [12-14])). On the other hand it is
necessary to underline that only isotope effect is a direct manifestation of
the mass effect in microphysics. It should be added that the direct
measurments of the energy of zero-point vibrations owing to isotope effect in
solids shows the good agreement of the experimental values with the results of
the calculation of quantum electrodynamics in solids [13, 14].

Below we describe shortly the manifestations of the isotope effect in
molecular as well as solid state physics (more details see [14]). The
discovery [15] of the new fullerene allotropes of carbon, exemplified by C60
and soon followed by an efficient method for their synthesis [14], led to a
burst of theoretical and experimental activity on their physical properties.
Much of this activity concentrated on the vibrational properties of C60 and
their elucidation by Raman scattering [15]. Comparison between theory and
experiment was greatly simplified by the high symmetry (Ih), resulting in only
ten Raman active modes for the isolated molecule and the relative weakness of
solid state effect [15], causing the crystalline C60 (c - C60) Raman spectrum
at low resolution to deviate only slightly from that expected for the isolated
molecule [15]. Since the natural abundance of 13C is 1.11% (see, e.g. [12]),
almost half of all C60 molecules made from natural graphite contain one or
more 13C isotopes. If the squared frequency of a vibrational mode in a C60
molecule with n13C isotopes is written as a series $\ \omega^{2}$ =
$\omega_{(0)}^{2}$ \+ $\omega_{(1)}^{2}$ \+ $\omega_{(2)}^{2}$ \+
$\omega_{(3)}^{2}$ \+ …… in the mass perturbation (where
$\omega_{\left(0\right)}$ is an eigenmode frequency in a C60 molecule with 60
12C atoms), nondegenerate perturbation theory predicts for the two totally
symmetric Ag modes a first - order correction given

$\frac{\omega_{(1)}^{2}}{\omega_{(0)}^{2}}$ = - $\frac{\text{n}}{\text{720}}$.
(3)

This remarkable result, independent of the relative position of the isotopes
within the molecule and equally independent of the unperturbed eigenvector, is
a direct consequence of the equivalence of all carbon atoms in icosahedral
C60. To the same order of accuracy within nondegenerate perturbation theory,
the Raman polarizability derivatives corresponding to the perturbed modes are
equal to their unperturbed counterparts, since the mode eigenvectors remain
unchanged. These results lead to the following conclusion [15]: The Ag Raman
spectrum from a set of noninteracting C60 molecules will mimic their mass
spectrum if the isotope effect on these vibrations can be described in terms
of first - order nondegenerate perturbation theory. It is no means obvious
that C60 will meet the requirements for the validity of this simple theorem. A
nondegenerate perturbation expansion is only valid if the Agmode is
sufficiently isolated in frequency from its neighboring modes. Such isolation
is not, of course, required by symmetry. Even if a perturbation expansion
converges, there is no a priori reason why second - and higher - order
correction to Eq. (3) should be negligible. As was shown in cited paper the
experimental Raman spectrum (see below) of C60 does agree with the prediction
of Eq. (3). Moreover, as was shown in quoted paper, experiments with
isotopically enriched samples display the striking correlation between mass
and Raman spectra predicted by the above simple theorem. Fig. 1 shows a high -
resolution Raman spectrum at 30 K in an energy range close to the high -
energy pentagonal - pinch Ag(2) vibration according to [15]. Three peaks are
resolved, with integrated intensity of 1.00; 0.95; and 0.35 relative to the
strongest peak. The insert of this figure shows the evolution of this spectrum
as the sample is heated. The peaks cannot be resolved beyond the melting
temperature of CS2 at 150 K. The theoretical fit yields a separation of 0.98
$\pm$ 0.01 cm-1 between two main peaks and 1.02 $\pm$ 0.02 cm-1 between the
second and third peaks. The fit also yields full widths at half maximum
(FDWHM) of 0.64; 0.70 and 0.90 cm-1, respectively. The mass spectrum of this
solution shows three strong peaks (Fig. 1b) corresponding to mass numbers 720;
721 and 722, with intensities of 1.00; 0.67 and 0.22 respectively as predicted
from the known isotopic abundance of 13C. The authors [15] assign the highest
- energy peak at 1471 cm-1 to the Ag(2) mode of isotopically pure C60 (60 12C
atoms). The second peak at 1470 cm-1 is assigned to C60 molecules with one 13C
isotope, and the third peak at 1469 cm-1 to C60 molecules with two 13C
isotopes. The separation between the peaks is in excellent agreement with the
prediction from Eq. (3), which gives 1.02 cm-1. In addition, the width of the
Raman peak at 1469 cm-1, assigned to a C60 molecule with two 13C atoms, is
only 30 % larger than the width of the other peaks. This is consistent with
the prediction of Eq. (3) too, that the frequency of the mode will be
independent of the relative position of the 13C isotopes within the molecule.
The relative intensity between two isotope and one isotope Raman lines agrees
well with the mass spectrum ratios. Concluding this part we stress that the
Raman spectra of C60 molecules show remarkable correlation with their mass
spectra. Thus the study of isotope - related shift offers a sensitive means to
probe the vibrational dynamics of C60.

Next examples of the dependence of the exciton spectra in solds on the isotope
effect demonstrate below. Isotopic substitution only affects the wavefunction
of phonons; therefore, the energy values of electron levels in the Schrödinger
equation ought to have remained the same. This, however, is not so, since
isotopic substitution modifies not only the phonon spectrum, but also the
constant of electron-phonon interaction (see [12]). It is for this reason that
the energy values of purely electron transition in molecules of hydride and
deuteride are found to be different. This effect is even more prominent when
we are dealing with a solid [16]. Intercomparison of absorption spectra for
thin films of LiH and LiD at room temperature revealed that the longwave
maximum (as we know now, the exciton peak ) moves 64.5 meV towards the shorter
wavelengths when H is replaced with D [17].

The mirror reflection spectra of mixed and pure LiD crystals cleaved in liquid
helium are presented in Fig. 2. For comparison, on the same diagram we have
also plotted the reflection spectrum of LiH crystals with clean surface. All
spectra have been measured with the same apparatus under the same conditions.
As the deuterium concentration increases, the long-wave maximum broadens and
shifts towards the shorter wavelengths. As can clearly be seen in Fig. 2, all
spectra exhibit a similar long-wave structure. This circumstance allows us to
attribute this structure to the excitation of the ground (Is) and the first
excited (2s) exciton states. The energy values of exciton maxima for pure and
mixed crystals at 2 K are presented in Table 22 of ref. [12]. The binding
energies of excitons E${}_{\text{b}}$, calculated by the hydrogen-like
formula, and the energies of interband transitions E${}_{\text{g}}$ are also
given in Table 22.

Going back to Fig. 2, it is hard to miss the growth of $\Delta_{\text{12}}$,
which in the hydrogen-like model causes an increase of the exciton Rydberg
with the replacement of isotopes . When hydrogen is completely replaced with
deuterium, the exciton Rydberg (in the Wannier-Mott model) increases by 20%
from 40 to 50 meV, whereas E${}_{\text{g}}$ exhibits a 2% increase, and at 2
$\div$ 4.2 K is $\Delta$E${}_{\text{g}}$ = 103 meV. This quantity depends on
the temperature, and at room temperature is 73 meV, which agrees well enough
with $\Delta$E${}_{\text{g}}$ = 64.5 meV as found in the paper of Kapustinsky
et al. [17]. The single-mode nature of exciton reflection spectra of mixed
crystals LiH${}_{\text{x}}$D${}_{\text{1-x}}$ agrees qualitatively with the
results obtained with the virtual crystal model (see e.g. Elliott et al. [18];
Onodera and Toyozawa [19]), being at the same time its extreme realization,
since the difference between ionization potentials ($\Delta\zeta$) for this
compound is zero. According to the virtual crystal model, $\Delta\zeta$ = 0
implies that $\Delta$E${}_{\text{g}}$ = 0, which is in contradiction with the
experimental results for LiH${}_{\text{x}}$D${}_{\text{1}}$-${}_{\text{x}}$
crystals. The change in E${}_{\text{g}}$ caused by isotopic substitution has
been observed for many broad-gap and narrow-gap semiconductor compounds (see
also [12]).

All of these results are documented in Table 22 of Ref.[12], where the
variation of E${}_{\text{g}}$, E${}_{\text{b}}$, are shown at the isotope
effect. We should highlighted here that the most prominent isotope effect is
observed in LiH crystals, where the dependence of E${}_{\text{b}}$ = f
(C${}_{\text{H}}$) is also observed and investigated. To end this section, let
us note that E${}_{\text{g}}$ decreases by 97 cm${}^{\text{-1}}$ when
${}^{\text{7}}$Li is replaced with ${}^{\text{6}}$Li.

Detailed investigations of the exciton reflectance spectrum in CdS crystals
were done by Zhang et al. [20]. Zhang et al. studied only the effects of Cd
substitutions, and were able to explain the observed shifts in the band gap
energies, together with the overall temperature dependence of the band gap
energies in terms of a two-oscillator model provided that they interpreted the
energy shifts of the bound excitons and n = 1 polaritons as a function of
average S mass reported earlier by Kreingol’d et al. [21] as shifts in the
band gap energies. However, Kreingol’d et al. [21] had interpreted these
shifts as resulting from isotopic shifts of the free exciton binding energies
, and not the band gap energies, based on their observation of different
energy shifts of features which they identified as the n = 2 free exciton
states (for details see [21]). The observations and interpretations, according
Meyer at al. [22], presented by Kreingol’d et al. [21] are difficult to
understand, since on the one hand a significant band gap shift as a function
of the S mass is expected , whereas it is difficult to understand the origin
of the relatively huge change in the free exciton binding energies which they
claimed. Very recently Meyer et al. [22] reexamine the optical spectra of CdS
as function of average S mass, using samples grown with natural Cd and either
natural S ($\sim$ 95% 32S), or highly enriched (99% 34S). These author
observed shifts of the bound excitons and the n = 1 free exciton edges
consistent with those reported by Kreingol’d et al. [21], but, contrary to
their results, Meyer et al. observed essentially identical shifts of the free
exciton excited states, as seen in both reflection and luminescence
spectroscopy. The reflectivity and photoluminescence spectra in polarized
light ($\overrightarrow{E}$ $\bot$ $\overrightarrow{C}$) over the A and B
exciton energy regions for the two samples depicted on the Fig. 3. For the
$\overrightarrow{E}$ $\bot$ $\overrightarrow{C}$ polarization used in Fig. 3
both A and B excitons have allowed transitions, and therefore reflectivity
signatures. Fig. 3 reveals both reflectivity signatures of the n = 2 and 3
states of the A exciton as well that of the n = 2 state of the B exciton.

In Table 18 of Ref. [14] the results of Meyer et al. summarized the energy
differences $\Delta$E = E (Cd34S) - E (CdnatS), of a large number of bound
exciton and free exciton transitions, measured using photoluminescence,
absorption, and reflectivity spectroscopy, in CdS made from natural S (CdnatS,
95% 32S) and from highly isotopically enriched 34S (Cd34S, 99% 34S). As we can
see from Fig. 3, all of the observed shifts are consistent with a single
value, 10.8$\pm$0.2 cm-1. Several of the donor bound exciton photoluminescence
transitions, which in paper [22] can be measured with high accuracy, reveal
shifts which differ from each other by more than the relevant uncertainties,
although all agree with the 10.8$\pm$0.2 cm-1 average shift. These small
differences in the shift energies for donor bound exciton transitions may
reflect a small isotopic dependence of the donor binding energy in CdS (see,
also [12]). This value of 10.8$\pm$0.2 cm-1 shift agrees well with the value
of 11.8 cm-1 reported early by Kreingol’d et al. [21] for the Bn=1 transition,
particularly when one takes into account the fact that enriched 32S was used
in that earlier study, whereas Meyer et al. have used natural S in place of an
isotopically enriched Cd32S (for details see [22]).

Authors [21] conclude that all of the observed shifts arise predominantly from
an isotopic dependence of the band gap energies, and that the contribution
from any isotopic dependence of the free exciton binding energies is much
smaller. On the basis of the observed temperature dependencies of the
excitonic transitions energies, together with a simple two-oscillator model,
Zhang et al. [20] earlier calculated such a difference, predicting a shift
with the S isotopic mass of 950 $\mu$eV/amu for the A exciton and 724
$\mu$eV/amu for the B exciton. Reflectivity and photoluminescence study of
natCd32S and natCd34S performed by Kreingol’d et al. [21] shows that for anion
isotope substitution the ground state (n = 1) energies of both A and B
excitons have a positive energy shifts with rate of $\partial$E/$\partial$MS =
740 $\mu$eV/amu. Results of Meyer et al. [22] are consistent with a shift of
$\sim$710 $\mu$eV/amu for both A and B excitons. Finally, it is interesting to
note that the shift of the exciton energies with Cd mass is 56 $\mu$eV/amu
[20], an order of magnitude less than found for the S mass (more details see
[12, 13]).

The brought examples clearly indicate mass dependence of the electron and
phonon states (see more details [14]) but on the other side it simply
underlines the primary importance in microphysics the difference of mass
between neutron (Mn) and proton (Mp). Really small difference in their masses
Mn \- Mp = 1.2333317 MeV leads to the bright effects in microphysics.
According to the last data [9], the experimental neutron-proton mass
difference of Mn \- Mp = 1.2333317 MeV is received as estimated
electromagnetic contribution Mn \- M${}_{p}\mid^{\text{em}}=$ -0.76 $\pm$ 0.30
MeV, and the remaining mass difference is determined to a strong isospin
breaking contribution of Mn \- M${}_{p}\mid^{\text{d-u}}$= 2.05 $\pm$ 0.30
MeV. In other words the last contribution is a result of difference in mass of
d - and u - quarks (see, also [10, 11]).

As we all know, the observed world - stars, planets, galaxy as well as
surrounding objects consist from the nuclei, neutrons, protons and electrons.
The mass of electrons has a small contribution to the total mass ( less than
0.1% (see, e.g. [1]). Therefore, that we knew that the origin of the mass of
the observed worlds needs to be elucidated the origin of nuclear mass. As we
know the nucleon consists from u - and d - quarks. But the mass of u - and d -
quarks is so small, that is their sum is a small part of the nucleon mass (1 -
2 % [6]). In modern physics of elementary particles it is considered that the
mass of nucleon arises from the spontaneous breaking of a chiral symmetry in
quantum chromodynamics (QCD) [23] and may be expressed over vacuum condensate
(see [5] and references therein).This model has an approximate formula which
expresses the mass of nucleon over quarks condensate [5]

m = [-2(2$\pi$)${}^{2}\langle$0$\mid\overline{\text{q}}$q$\mid$0$\rangle$]1/3
(4),

where m is nucleon mass, $\langle$0$\mid\overline{\text{q}}$q$\mid$0$\rangle$
is quarks condensate, q is the field of u - or d - quarks. The chiral symmetry
in QCD tresult in the expression for the quarks condensate (so called Gell -
Mann - Oakes - Renner formula [24])

$\langle$0$\mid\overline{\text{q}}$q$\mid$0$\rangle$ =
-$\frac{\text{1}}{\text{2}}\frac{\text{m}_{\pi}\text{
f}_{\pi}}{\text{m}_{u}\text{ + m}_{\text{d}}}$ (5).

Here mπ and fπ are the mass and decay constant of $\pi$ \- meson. The defined
value of quarks condensate on the ground of $\tau$ \- decay [5,6] equals

$\langle$0$\mid\overline{\text{q}}$q$\mid$0$\rangle$ = - (254 MeV)3 $\pm$ 10
%. (6)

Put the last value into the expression (4) it gives the nuclon’s mass m = 1.08
GeV, when the experimental value of nucleon’s mass equals m = 0.94 MeV. From
comparison of these values we see that the difference between experimental
value of m and theoretical estimation is 0.15 GeV, that surpasses the
experimental value of the difference Mn \- Mp = 1.2333317 MeV much order. The
last one means that in such model (as well as in the model of constituent
quarks) we have neither the mass difference of the nucleons nor its number in
nuclei and, consequently, isotope effect. But the experimental manifestations
of the isotope effect was demonstrated above in the different branches of
microphysics. Considering the quarks structure of nucleon (the wavefunction of
the neutron is udd, and for proton one is uud) that is the quark strucure
indicates the different construction of the neutron and proton, but this model
doesn’t quantative describe the mass of nucleons.

Thus, the origin of the isotope effect is closely connected with the different
origin of u - and d - quarks and with solution the spectrum and hierarhy of
the elementary particles mass and more common with the solution of the nature
of mass (see, also [25]).

Acknowledgements. I deeply thank to Prof. B.L. Ioffe for the enlighting
discussion on the origin of mass, and Prof. P. Sneider for improving my
English.

Figure captions.

Fig. 1. a - unpolarized Raman spectrum in the frequency region of the
pentagonal - pinch mode, for a frozen sample of nonisotopically enriched C60
in CS2 at 30 K. The points are the experimental data, and the solid curve is a
three - Lorentzian fit. The highest - frequency peak is assigned to the
totally symmetric pentagonal - pinch Ag mode in isotopically pure 12C60. The
other two peaks are assigned to the perturbed pentagonal - pinch mode in
molecules having one and two 13C - enriched C60, respectively. The insert
shows the evolution of these peaks as the solution is heated. b - the points
give the measured unpolarized raman spectrum in the pentagonal - pinch region
for a frozen solution of 13C - enriched C60 in CS2 at 30 K. The solid line is
a theoretical spectrum computed using the sample’s mass spectrum, as described
in the text (after [15]).

Fig. 2. Mirror reflection spectra of crystals: 1 - LiH; 2 - LiHxD1-x; 3 - LiD;
at 4.2 K. 4 - source of light without crystal. Spectral resolution of the
instrument is indicated on the diagram (after [12]).

Fig. 3. a - Reflection spectra in the A and B excitonic polaritons region of
CdnatS and Cd34S at 1.3K with incident light in the
$\overrightarrow{\text{E}}$ $\perp$ $\overrightarrow{C}$. The broken vertical
lines connecting peaks indicate measured enrgy shifts reported in Table 18 of
Ref. [14]. In this polarization, the n = 2 and 3 excited states of the A
exciton, and the n = 2 excited state of the B exciton, can be observed. b -
Polarized photoluminescence spectra in the region of the A${}_{\text{n = 2}}$
and A${}_{\text{n = 3}}$ free exciton recombination lines of CdnatS and Cd34S
taken at 1.3 K with the $\overrightarrow{\text{E}}$ $\perp$
$\overrightarrow{C}$. The broken vertical lines connecting peaks indicate
measured enrgy shifts reported in Table 18 of Ref. [14] (after [22]).

References.

1\. M. Jammer, Concepts of mass in classical and modern physics, Harvard
University Press, Cambridge - Massachsets (1961).

2\. A. Einstein, Ann. Phys. (Leipzig) 20, 371 (1906).

3\. L.D. Landau, E.M. Lifshitz, The classical theory of fields, Pergamon, New
York (1958).

4\. L.B. Okun, Physics Today, June 1989; Uspekhi Fiz. Nauk 158, 512 (1989) (in
Russian).

5\. B.L. Ioffe, Uspekhi Fiz. Nauk 171, 1273 (2001) (in Russian); Progr. Part.
Nucl. Phys. 56, 232 (2006).

6\. C.D. Frogatt, Surveys High Energy Physics 18, 77 (2003); The Problem of
Mass , ArX:hep - ph/0312220.

7\. A. Dobado and A.L. Maroto, Phys. Rev. D60, 104045-9 (1999).

8\. J.J. Kelly, Phys. Rev. C70, 068202 (2004).

9\. S.R. Beane, K. Originas and M.J. Savage, Nucl. Phys. B768, 38, (2007).

10\. G.A. Miller, A.K. Opper, E.J. Stephenson, Annual Review of Nuclear
Science 56, 253 (2006).

11\. G.A. Miller, Phys. Rev. Lett. 99, 112001 (2007); The Neutron Negative
Central Charge Density: an Inclusive - Exclusive Connection, ArXiv 0806.3977.

12\. V.G. Plekhanov, Phys. Reports 410, 1 (2005).

13\. M.Cardona, M.L.W. Thewalt, Rev. Mod. Phys. 77, 1173 (2005).

14\. V.G. Plekhanov, will be published.

15\. J. Menendez and J.B. Page, Vibrational spectroscopy of C60, in, M.
Cardona and G. Guntherodt, eds., Light Scattering in Solids VIII, Springer,
Berlin - Heidelberg (2000) (Vol. 76 in Topics in Applied Physics).

16\. V.G. Plekhanov, Isotope effects in solid state physics, Academic Press,
San Diego (2001).

17\. A.F. Kapustinsky, L.M. Shamovsky, K.S. Bayushkina, Acta Physicochim.
(USSR) 7, 799 (1937).

18\. R.J. Elliott, J.AA. Krumhansl, P.L. Leath, Rev. Mod. Phys. 46, 465
(1974).

19\. Y.Onodera and Y. Toyozawa, J. Phys. Soc. Japan 24, 341 (1968).

20\. M. Zhang, M. Ghieler, T. Ruf, Phys. Rev. B57, 9716 (1998).

21\. F.I. Kreingol’d, K.F. Lider, M.B. Shabaeva, Fiz. Tverd. Tela 26, 3940
(1984) (in Russian).

22\. T.A. Meyer, M.L.W. Thewalt and R. Lauck, Phys. Rev. B69, 115214-5 (2004).

23.J. Grasser and H. Leutwyller, Phys. Reports 87, 77 (1982); H. Leutwyller,
Insights and Puzzles in Light Quark Physics, ZrXiv:hep - ph/070063138.

24\. M. Gell-Mann, R.J. Oakes, B. Renner, Phys. Rev. 175, 2195 (1968).

25\. I. F. Ginzburg, Uspekhi Fiz. Nauk (Moscow) 179, 525 (2009) (in Ruassian).