Anomalous Structure in Critical Screening Parameters of the ECSC Potential

1. Introduction

1.1. Background

Critical binding in Screened Coulomb potentials such as Yukawa/Debye, Hulthén, and ECSC potentials has remained a phenomenon of intense and sustained interest for nearly a century. The fundamental question concerns the precise values of the screening (Debye) length $\mathcal{D}$ at which quantum states of various angular momenta just become bound or unbound, which has both theoretical and practical ramifications. This question has stimulated a very large body of scientific work by a wide spectrum of physicists and mathematicians, for both bound states and resonances at positive energy. In our current paper we confine our attention to critical binding parameters – the same methods are suitable for calculating eigenvalues however.

Our recent paper [1] was the stimulus for exploring the first mystery addressed in this paper, and the landmark work by Rogers et al [2] and Lam and Varshni [3] are the stimulus for exploring the second mystery. Some familiarity with our first paper is strongly recommended for contextualizing the contents and methods used in this second paper.

Recent progress has been reviewed in [4,5,6,7], with applications to mapping the parameter space of the generalized exponential cosine screened coulomb potential, and attention to circular Rydberg states. A good collection of earlier citations on ECSC potential can be found in Dutt [8].

Recently (2026) we have made a first contribution [1] to this literature which can be viewed as a prelude to the current paper – the two papers naturally go together. From our very accurate computed critical binding parameters we discovered anomalous structure in the ECSC critical screening values, with a strikingly regular periodicity in the ${\mathcal{D}}_c$ vs n plots. This is most clearly visible when the generic $n^2$ trend is divided out, as shown in Figure 1 and Figure 2. The contrast between the Yukawa and ECSC curves is dramatic. Hulthén values are smooth, like those of the Yukawa.

The ECSC potential, $-{\displaystyle \frac{e^{-\mu r}}{r}}cos\left(\mu r\right)$ (where $\mu \equiv 1/\mathcal{D}$ – for simplicity we omit writing the centrifugal potential), is identical to the Yukawa/Debye potential $-{\displaystyle \frac{e^{-\mu r}}{r}}$, except (crucially) for the $cos\left(\mu r\right)$ modulating factor. This factor gives the potential an infinite sequence of positive energy barriers, and associated negative energy wells, the magnitudes of which are damped by the exponential decay and $1/r$ factors. Evidently both Yukawa and ECSC potentials reduce to the Coulomb potential $-1/r$ as $\mu \to 0$, i.e. $\mathcal{D}\to \infty$.

A primary goal of our current paper is to explain the origin of this sawtooth structure and its surprisingly regular periodicity $\Delta n\approx 40$. Another goal is to solve a second mystery that apparently has remained unsolved for more than 50 years: the strangely consistent linear relation between the total number of bound states $n^{\ast }$ vs $\mathcal{D}$ for both the Yukawa [2] and ECSC [3] potentials. We will quantitatively solve both of these mysteries, using the complementary methods employed in our previous paper: the Phase Method (PM), and our variant-semiclassical (“neoclassical", NC) method, which extends the accuracy and range of application of standard quasi/semi-classical method for calculating eigenvalues, here adapted instead for the purpose of calculating critical parameters. The rationale for NC is briefly explained below, and in [9]; several useful formulae using it also are given below.

The PM method revealed the sawtooth structure, and NC and PM methods allow us to explain it and its consistent period of $\approx 40$ , and also to predict the $\mathcal{D}$-values at which the higher order wells produce similar behavior. A simple equation is given to predict these, and such a prediction is confirmed using the PM for all screening lengths up to $\mathcal{D}<{10}^6$ au. We further use the PM to explore this “tick-tock” mechanism with modified ECSC potentials. A simple quartic model system is considered from the point of view of both methods to further illustrate the “tick-tock” mechanism responsible for the sawtooth structure. The NC method is further used to explain the approximately linear $\sqrt{{\mathcal{D}}_c}$ vs n and l in light of the scaling symmetry between $\mathcal{D}$ and coupling constant, which is characteristic of this class of potentials. Finally, the combination of this linear relation with a simple geometrical state counting procedure derived from it, solves the second mystery. Connection with the Coulomb SO(4) symmetry in the asymptotic limit also is discussed.

We predict these latter results are relevant to the whole class of short-range, linear coupling constant potentials that were studied in the seminal paper of Klaus and Simon [10].

1.2. Overview of New Results in This Paper

Our current paper contains a substantial number of new results. To assist the reader we list most of them below, preceded by the subsection number where they are introduced, and tagged by the method(s) used. Results that are carried over from our previous paper for continuity are shown in parentheses.

• PM: Tables of $D_c$ , accurate to 30 significant digits, for Yukawa, ECSC, and Hulthén potentials for all states $\mathcal{D}\le {10}^5$ au and $l=0-12$ evailable as supplementary data
• 1.3 (PM: sawtooth structure exists in ECSC, constant period 40 for $l=0$)
• 2.1 (review of Phase Method PM)
• 2.2 explanation of Neoclassical approach NC
• 2.3 NC: example: homogeneous potentials, with formulae
• 2.4 NC: example calculation of Yukawa number of bound states vs $\mathcal{D}$ to relative accuracy $\le {10}^{-3}$ for all l, n, two parameters each, with code
• 3.1 PM: modified potential: sawtooth structure depends on secondary well
• 3.1 PM: modified potential: period of sawtooth structure depends $\approx \sqrt{\mathrm{well}\mathrm{depth}}$
• 3.1 PM: modified potential: onset of sawtooth structure depends on height of barrier
• 3.2 NC: period 40 quantitatively explained by “tick-tock” mechanism: $1+{\mathcal{J}}_1/{\mathcal{J}}_2$
• 3.2 NC: sawtooth structure reproduced ab initio using neoclassical approach
• 3.3 NC: approximate formula for locations of tick-tock glitches from higher wells
• 3.4 PM: confirmation of prediction of first tick-tock glitch from tertiary well
• 3.5 NC+PM: example tick-tock calculation for asymmetric quartic double-well
• 3.6 (PM: plot Yukawa $n^{\ast }$ vs $\mathcal{D}$ plot, compare to Rogers 1970)
• 3.6 PM: plot ECSC $n^{\ast }$ vs $\mathcal{D}$, compare to Lam and Varshni 1972
• 3.6 PM: plot Hulthén $n^{\ast }$ vs $\mathcal{D}$, same linear trend as in seen for Yukawa and ECSC
• 3.6 PM: $n^{\ast }$ vs $\mathcal{D}$ plot shows ECSC states are clustered by principal quantum number $\overline{n}$
• 3.7 PM: ECSC $n^{\ast }$ vs $\mathcal{D}$ mystery is solved by arithmetic series accumulation of states
• 3.7 PM: plots of $\sqrt{\mathcal{D}}$ vs n and l show near-uniform density of states for all potentials
• 3.7 PM: geometrical state accumulation solves the second mystery for all three potentials, but begs the question re the reason for the near-linearity, addressed next
• 3.7 (PM: Yukawa contour plots showed striking near-linear relation between $\sqrt{\mathcal{D}}$, n, l)
• 3.7 PM: contour plots for Hulthén and ECSC also show same
• 3.7 PM: linear fit coefficients determined between $\sqrt{\mathcal{D}}$, n, l; residual contour plots
• 3.8 First order perturbation theory gives near-linear contour plot for Yukawa
• 3.9 NC: partial derivatives with respect to parameters under phase space integral gives numerical estimates of coefficients, that compare well to fits to PM data
• 3.9 NC: explanation for trivariate linear relation given in terms of scaling symmetry
• 3.10 NC: linear relation further explained in terms of $\beta =l(l+1)/\left(2\mathcal{D}\right)$ and scaling symmetry, ODE condition for strict linearity
• 3.10 NC: explicit approximate expressions for $\mathcal{J}(q,\beta )$ for Hulthén, Yukawa, ECSC
• NC these properties derive from linear coupling constant a la Klaus and Simon

1.3. Anomalous Sawtooth Structure in ${\mathcal{D}}_c/n^2$ vs n

In our previous paper [1] we presented calculations to 60-digit-accuracy of critical screening parameters (inverse screening (Debye) lengths up to $\mathcal{D}={10}^3$ au for Yukawa/Debye, Hulthén, Pseudo-Hulthén, and Exponential Cosine Screened Coulomb (ECSC) potentials, for angular momenta $l=0$ – 20. Initial results also were briefly presented up to $\mathcal{D}={10}^5$ for $l=0$ – 12 at 30 digits precision for Yukawa, and ECSC potentials. These tables, plus additional tables for the Hulthén potential, are key data used in the current paper, and are deposited as supplementary data for others to use.

The asymptotic dependence of the critical screening lengths ${\mathcal{D}}_c$ for the various bound states (indexed by the integer n for each l-manifold), to which all l curves converge for large $\mathcal{D}$, is ${\mathcal{D}}_c\propto n^2$. Dividing out this asymptotic trend reveals more interesting structure, as shown in Figures 1 and 2, for Yukawa and ECSC respectively. The sawtooth structure in ECSC is striking, and demands explanation.

There are several reasons this structure may have gone unreported previously. One is that it is a subtle deviation from the smoothness of the raw data, and most clearly seen when the predominant $n^2$ dependence is divided out. Another is the upper range of $\mathcal{D}$ has to be high enough for the periodicity of the jumps to be seen clearly. These data extended to $\mathcal{D}={10}^5$. A third possibility is that hints of such structure indeed might have been seen, but went unreported because there was no explanation for it, and an isolated glitch may have been thought to be a numerical error. Using the highly accurate and robust PM and extending the range made the sawtooth structure very clear. Its nature and origin are explained below, as well as other properties of all ECSC, Yukawa, and Hulthén potentials.

1.4. Nature of Yukawa, ECSC, and Hulthén Potentials

In the critical binding problem, we are concerned about state-binding at zero energy. In the following we will explicitly write the centrifugal potential term. The behavior is rather different depending whether $l=0$ or $l>0$, so we will normally treat these cases separately. In r-space, as is described at length in our first paper [1], screened Coulomb potentials can show structure on very long length scales, depending on the value of $\mathcal{D}=1/\mu$, and with very shallow potential-well depths, but such that bound states can form under the right conditions. $\mu \to 0$ gives the Coulomb limit.

For many purposes it is useful to change variables to $z=r/\mathcal{D}$. In that case, the Yukawa [13] (effective) potential becomes $-\mathcal{D}{e}^{-z}/z+l(l+1)/\left(2z^2\right)$. Evidently, for any finite value of l, as $\mathcal{D}\to \infty$ in the asymptotic large n limit, the finite value of l becomes irrelevant, and the behavior is independent of l. As shown in [1] the classical turning points (at zero energy) are $z_1=-W(0,-\beta )$, and $z_2=-W(-1,-\beta )$, where $\beta =l(l+1)/\left(2\mathcal{D}\right)$. W is the Lambert W function (ProductLog in Mathematica). The 0 or $-1$ in the first argument of W specifies the appropriate branch of the function. These are the roots of the equation $Y_{\mathrm{yukawa}}=z{e}^{-z}=\beta$. If $l=0$, $U_{\mathrm{yukawa}}$ is a monotonically increasing function, approaching zero as $r\to \infty$, with a Coulomb-like attractive $-1/r$ singularity as $r\to 0$. In this case $\beta \to 0$ and $z_1=0$ and $z_2=\infty$. If $l>0$, there is a positive $1/r^2$ repulsive singularity as $r\to 0$, which, combined with the long range attractive potential from the exponential term, under the right combination of $\mathcal{D}$ and l, can create a well at $r>0$ that is sufficient to bind a state. This is analyzed in detail below.

The Hulthén [14] potential is of similar form as the Yukawa Potential, but it replaces the $1/r$ factor of the exponential with a $1/\rho$ factor, where $\rho =(1-e^{-\mu r})/\mu$. This behaves much the same as the Yukawa potential at small r but approaches the constant $1/\mu$ for large r. We can use a similar analysis [1] to the others in z-space, in which $Y_{\mathrm{Hulthen}}\left(z\right)=z^2/(e^z-1)$. This function has a single peak similar to the Yukawa case, with a peak at $z_0=2+W(-2/e^2)\approx 1.59362$ and a value at the peak ${\beta }_{\max }=-W(-2/e^2)(2+W(-2/e^2)\approx .64761$, where, again W is the Lambert W function. Classical turning points must be found numerically by root-finding, but this is straightforward. The value of $\beta$ and ${\beta }_{\max }$ give a concise summary of the tradeoffs between l and $\mathcal{D}$ regarding the ability to bind states.

The ECSC [3] (effective) potential in z-space takes the form $-\mathcal{D}{e}^{-z}/z+l(l+1)/\left(2z^2\right)$. As in the Yukawa potential, for $l=0$ there is an attractive $-1/r$ singularity that is dominated as $r\to 0$ by a repulsive $\propto 1/r^2$ singularity for $l>0$. The ECSC potential differs from the Yukawa potential by the modulating $cos\left(z\right)$ factor which introduces an infinite sequence of positive energy barriers alternating with negative energy wells, both of which are affected by the value of $\mathcal{D}$ and l, that have anti-correlated effects. increasing $\mathcal{D}$ induces more state-binding, while increasing l does the opposite. The value of $\beta$ to some extent expresses the trade-off between them. This is put to further use in subsection 3.7.

For $l=0$, the turning points are simply $z_1=0$, $z_2=\pi /2$ for the primary well; $z_1=3\pi /2$, $z_2=5\pi /2$ for the secondary well; and $z_1=7\pi /2$, $z_2=9\pi /2$ for the tertiary well etc. These are discrete wells because of the positive energy barriers between them. For $l>0$, numerical root-finding of $Y_{\mathrm{ECSC}}\left(z\right)=z{e}^{-z}cosz=\beta$ is used to determine the turning points, depending on the value of $\beta$ as given above. The $Y_{\mathrm{ECSC}}\left(z\right)$ function has single discrete smooth peaks within each well; root-finding is straightforward, provided $\beta$ does not exceed the maximum value of $Y\left(z\right)$ at its peak ${\beta }_{\max }$, which provides an upper limit on l (for fixed $\mathcal{D}$) and a lower limit on $\mathcal{D}$ (for fixed l) at which a state might possibly become bound. Each well in the ECSC potential effectively has its own corresponding $Y\left(z\right)$ and ${\beta }_{\max }$, corresponding to their own ranges of z. The multi-well structures are analyzed in detail below in regard to phase space quantization of the wells.

2. Materials and Methods

All computations were done using code written by the author in Mathematica V14.3 [11].

In this paper, as in our previous one, we primarily use two methods in tandem: “Phase Method” as described in our first paper, a precise and robust fully quantum computational approach; and a variation on conventional semiclassical method that we call “neoclassical” (NC).

A note on terminology: the NC approach is principally based on classic foundational work of Dunham [15], and further elaborated by Bender [16,17]. It is adjacent to a vast body of semiclassical theory from which it differs significantly but only in a specific and narrow way, and it is not in conflict with it. In an attempt to prevent confusion, we refer to it as “neoclassical” (NC). This terminology is by analogy with classical music – neoclassical music (e.g. Stravinsky, Prokofiev) honored and respected the classical forms and ideas, but also reinterpreted them, and applied them to new ends. We do a similar thing here.

2.1. Phase Method

The Phase Method [1,9] is a robust, simple to use, and accurate method for calculating eigenvalues (and as needed, wave functions) by direct solution of the time independent Schrödinger equation for 1D, central, or other separable potentials, which can be smooth, discontinuous, or random; and it handles near-degeneracies very well. It is related t o conventional Shooting Methods, but differs from them in several ways, e.g. by seeking divergent solutions, instead of convergent ones – in what we call the “Shoot First” method. The PM is a significant extension beyond that. It automatically finds all of (a finite number of) states within a specified range energy range using a parallel evolutionary search process; importantly, it does not require any guessing of eigenvalues, assumptions about parity, or tweaking of initial conditions to obtain convergence. Accuracy of eigenvalues to 30–60 digits or more is routine. Fewer digits is fine too – the ultimate accuracy is “dialed-in” by the choice of number of digits of working precision in the calculation. The PM runs very well on ordinary desktop computers, even those with a minimal memory footprint, typically in minutes for most potentials and a few dozen eigenvalues. The run time scales more or less with the number of eigenvalues. It uses arbitrary precision arithmetic, so it is not currently suited for running on GPUs, but otherwise is highly parallel. It should run well on everything between Raspberry Pi computers (with its bundled Mathematica) to supercomputers. These calculations were mostly done on several 16GB M4 Mac Minis.

In our previous paper [1] we applied the PM to the calculation of critical screening parameters ${\mu }_c=1/{\mathcal{D}}_c$ for screened Coulomb potentials, to 60 digit accuracy. Adaptation to the critical screening problem required only minor changes to our simple eigenvalue code that is posted for download on our website. More information on the workings of the PM and consideration regarding its application near zero energy is described in [1].

2.2. “Neoclassical” Variant of Semiclassical Method

In our recent paper [1] and book [9] we introduced a continuous real parameter, $\nu$, in conjunction with phase space quantization, which needs further discussion. Its use represents a narrowly focussed but significant deviation from the standard semiclassical development, while remaining fully compatible with it. In fact its justification goes back to very early exact WKB theory. In using this modified approach here we restrict our attention to determination of energy eigenvalues, or equivalently, critical binding parameters. We do not need to compute asymptotic forms of wave functions by semiclassical means, which entails more detailed exposition [17,18,19,20,21]. Our neoclassical approach, in contrast, is very simple and direct, but still sufficient for our purposes.

As highlighted by Bender [16], the 1932 paper by Dunham [15] contains a foundational equation. It is equivalent to Bender’s [16] equation 10. Dunham extended the earliest WKB formulations by deriving an exact expansion in powers of ℏ, and expressed the energy quantization condition in terms of an infinite sum of contour integrals in the complex plane, where the integration contours enclose the classical turning points. This is elaborated in Bender and Oszag [17]. The various terms in the series involve successively higher derivatives of $p=\sqrt{2(\mathcal{E}-U)}$.

Crucially, single-valuedness of the wave function requires the infinite set of contour integrals to sum to an integer times $\pi$. This is the fundamental quantization condition for energy levels that we use. Instead of evaluating them term by term, however, we divide the total sum into two parts: the zeroth order term (classical Bohr-Sommerfeld-Wilson phase space integral), plus the total sum of all the rest, which we call $\nu \ast \pi$. The value of the first order term depends on the behavior at the turning points, and typically simplifies to a constant (e.g. $-\pi /2$); we incorporate that value into our parameter $\nu$. Each of the odd order terms higher than the first evaluate to zero, but the even order terms generally do not. This behavior is explained in Bender and Orszag [17]. The higher order terms typically introduce correction terms that can be expressed as a function of energy, or in other ways. We approximate $\nu$ as a constant plus a decreasing function of n, such as a Padé Approximant, an inverse power series in n, or a decreasing power law or exponential. For many purposes just setting $\nu$ to a constant is sufficient, but the value generally will differ from the canonical one, which is based on the truncated series. In principle this can be made exact using a suitable function of n.

Each successive term in the exact WKB expansion contains increasingly higher order derivatives of the potential; it can be viewed as a kind of gradient expansion. Many others e.g.  [16,17,20,21,23,24,25,26] have explained this approach very well. But all of them truncate the series at low order, in part because of the difficulty in evaluating the various contour integrals of higher order. Bender and Orszag go the farthest in evaluating these terms, and examining convergence of the series. Liverts and Barnea [22] evaluate higher order terms, to good effect. Smooth potentials can be rather accurately represented by the lowest order terms; some, e.g. harmonic oscillator and Morse potential, naturally terminate after a couple of terms.

This approach was developed in [9] for computing eigenvalues in homogeneous potentials, and others. Useful equations are given below using $\nu$ for homogeneous potentials. We use this NC method below for approximating critical binding parameters for the doubly-infinite set of n and l values for the Yukawa potential, which are calculated to $\approx {10}^{-3}$ relative precision, using a simple $\nu$ function, a decreasing power law needing only two parameters each in n, and in l. The calculation is implemented in about a dozen lines of code (shown below); it executes in less than a millisecond on a laptop computer.

The lowest order term in this expansion gives the most rudimentary Bohr-Sommerfeld-Wilson quantization, in which the phase space integral – i.e. twice the integral between the classical turning points of the momentum, expressed as a function of position – is equated to an integer multiple n of Planck’s constant $h=2\pi \hslash$. This BSW criterion is equivalent to requiring that the approximate wave function vanish at the classical turning points – i.e. an integer number of half-wavelengths resides between the classical turning points. This is what one gets when solving the classical Hamilton-Jacobi equation. But this usually gives a rather poor estimate of the true eigenvalue (except for the infinite square well, for which the criterion is exact), because tunneling actually occurs in most cases, which lowers the energy. Instead, rather than viewing the BSW value as a poor approximation to the true eigenvalue, we have argued [9] that it is more useful to view it as an inclusive upper bound on the true eigenvalue; and the $n-1$ BSW value is an exclusive lower bound. The value of $\nu$ interpolates between those limits in n-space. $\nu$ is, in part, a measure of the degree of tunneling.

Bender and Orszag [17] show that the usual value of ${\displaystyle \frac{1}{2}}$ (for systems with two smooth turning points) emerges naturally from direct evaluation of the first order term of Dunham’s series. Alternatively, arguments from topology and symplectic geometry [18,19,20] give the the same result; the topological concept “Maslov Index” is often applied in this context. This increment to n is an integer multiple of ${\displaystyle \frac{1}{4}}$, but it is only part of $\nu$. Yet, the still-neglected higher order terms can be quite limiting to the accuracy at low quantum numbers, and for rapidly varying potentials. These limitations of the standard formulation have caused many to view the use of phase space quantization for determination of eigenvalues for small n to be unacceptably inaccurate. In the NC approach, the implicit summation of all the terms in Dunham’s series is in principle exact, and applicable down to $n=1$, to hard potentials, or both.

The conventional low-order truncation of the series, although convenient, introduces other difficulties, which require workarounds. For central potentials with centrifugal potential terms, the “Langer Replacement” modifies the exact $l(l+1)$ factor in the centrifugal potential to ${(l+{\displaystyle \frac{1}{2}})}^2$, which effectively adds a positive $1/r^2$ term to it, giving accurate results despite use of the truncated WKB expansion. In this case the two approximations (truncation, and Langer l-shift) are arranged to cancel each other.

For the typical two-sided well, $\nu =0$ would imply there is no tunneling. If $\nu =1/2$ it implies moderate tunneling; this is the standard WKB value, and it is exact for harmonic oscillators, and Morse Potentials. Potentials that are softer (e.g. $\left|x\right|$, ${\left|x\right|}^{1/2}$) than harmonic oscillator have slightly larger values than $1/2$, but not much, and are still considerably less than 1.

Our development of both the NC and PM started with thinking about “mechanical similarity” [27]: the requirement of consistent dynamics in Lagrangian mechanics for similar-shaped trajectories, upon rescaling them in space and time. This analysis is quite different from ordinary dimensional analysis. The effects on the Lagrangian of rescaling of space and time implies that for homogeneous potentials of degree k, energy must scale with action S as $\mathcal{E}\propto S^{\mu }$ where $\mu \equiv {\displaystyle \frac{2k}{2+k}}$. Our thought, then, was that in the most basic form of semiclassical quantum mechanics, BSW quantization of the action integral $\mathcal{J}$ also quantizes the energy, and there must be a simple scaling formula for the energies in homogeneous potentials. This indeed is the case, and results are presented below, to illustrate the NC approach. After becoming aware of the existence of the exact Dunham formula, as elaborated by Bender et al [16], we realized such an approximate expression could be made exact by incorporating the “Dunham-shift” ($n\to n-\nu$). Furthermore, the phase-space-quantization with Dunham-shift approach evidently also applies to arbitrary (not only homogeneous) potentials. Of course this begs the question as to what is the correct value of $\nu$, but there is no disputing that $\nu$ indeed does have a value. That observation alone is sometimes sufficient; there are many alternative routes to answering the question of what’s $\nu$.

We adapt this NC approach to calculating critical screening lengths ${\mathcal{D}}_c$ values at zero energy. Instead of finding energy eigenvalues for a fixed potential/screening length, instead we find the critical screening length at a fixed energy, i.e. zero. The Phase Method implementation is essentially the same. In our previous paper we further used the NC approach to determine the asymptotic dependence of n on $\mathcal{D}$, and conditions for circular Rydberg states to form.

2.3. NC Applied to Homogeneous Potentials

As mentioned above, this approach allowed us [9] to derive simple formulae for the approximate eigenvalues of positively homogeneous potentials of the form

$$U\left(x\right)={\displaystyle \frac{ϵ}{2}}{\left|{\displaystyle \frac{x}{a}}\right|}^k,$$

(1)

where evidently $ϵ$ and a are energy and length scale parameters, and k (real-valued, not necessarily integer, or even rational) is the degree of homogeneity of the potential, and $\eta ={\displaystyle \frac{\hslash }{a\sqrt{m\left|ϵ\right|}}}$ is a dimensionless parameter; m is the (reduced) mass. Expressions similar to equations 2 and 3 (apart from the crucial introduction of the parameter $\nu$) have been independently derived by others, but their expressions are different in form, and the omission of $\nu$ implies that, as written, they fail for infinite square well, whereas our expressions are exact and smoothly interpolate in k.

Rather than being absorbed by the parameter $ϵ$, we have kept the prefactor of ${\displaystyle \frac{1}{2}}$ to ease comparison with the most basic benchmark, harmonic oscillator. Here we briefly present these expressions because they illustrate the utility of NC method, and also simply are useful, even for those (e.g. students) with just scientific calculators/apps at their disposal.

First, for such potentials, it is not difficult to show [9] that the eigenvalues scale exactly with the dimensional parameters as $\mathcal{E}\propto ϵ{\eta }^{\mu }$. This allows known eigenvalues for any k computed for a given set of potential parameters $ϵ$ and a, or values of ℏ, or mass m, to be rescaled for any other potential parameters (of the same k).

Second, we derived a simple generic formula, applicable for $k>-2$:

$$\mathcal{E}\propto ϵ{\left[\eta (n-\nu )\right]}^{\mu }.$$

(2)

We have since learned that a similar formula was independently derived by Quigg and Rosner [21] and perhaps others. Again, the key ingredient is the variable $\nu$, which allows these formulae to perform well at low n and for hard potentials, where usually WKB dare not tread.

This simple expression allows for easy determination of the dependence on n of eigenvalues in different homogeneous potentials. For example $\mathcal{E}\propto 1/n^2$ for a Coulomb potential, for which $\mu =-2$, and $\nu =0$. As $k\to \infty$, $\mathcal{E}\propto n^2$. Harmonic oscillators have $\mathcal{E}\propto (n-{\displaystyle \frac{1}{2}})$ with $\nu ={\displaystyle \frac{1}{2}}$. As $k\to -2$ the behavior enters an interesting realm of singular potentials and exponential scaling of bottomless eigenvalues, a discussion of which is outside the scope of the current paper.

The key question here is the precise value (or function) to choose for $\nu$. Fortunately the ansatz is useful without needing to know $\nu$ precisely. The choice of its value is guided to a first approximation by the ordinary (topological, negative) Maslov-index, but deviations from it are not only allowed but indeed required in most cases, because they account for the missing terms in the exact Dunham expansion.

In this paper we will express $\nu ={\nu }_0+\Delta \nu$, where ${\nu }_0$ is the canonical value (e.g. $1/2$) and $\Delta \nu$ is a decreasing function of n (or l). But even a constant approximate $\nu$ can improve accuracy for low n but still become irrelevant as $n\to \infty$, except for potentials with infinitely hard walls such as infinite square well. In such a case $\nu =0$ is appropriate because there is zero tunneling. For a general potential, the correct values of $\nu$ can be determined accurately by back-solving for the $\nu$ required to give the exact (or precise) values such as those from the Phase Method. This was mapped out for homogeneous potentials in [9], in which case $\Delta \nu$ was well-represented by a Padé Approximant for all n and k.

As an example of this approach, our approximate eigenvalue equation [9] has the form

$${\mathcal{E}}_{n,k}\approx {\displaystyle \frac{ϵ}{2}}{\left[(n-\nu )\eta \sqrt{\pi }{\displaystyle \frac{(1/\mu )!}{(1/k)!}}\right]}^{\mu },$$

(3)

where $n!\equiv \Gamma (n+1)$, and we number quantum states $n=1,2,\dots$. As before $\mu ={\displaystyle \frac{2k}{2+k}}$ and $\eta ={\displaystyle \frac{\hslash }{a\sqrt{m\left|ϵ\right|}}}$. This equation is exact for harmonic oscillators, and for the infinite square well, and it interpolates smoothly between $k=2$ and $k\to \infty$, and down to $k=0$. As above, the prefactor of ${\displaystyle \frac{1}{2}}$ is retained to facilitate comparison with harmonic oscillator. Other authors have derived equivalent formulae (albeit written differently), but with the key difference being our introduction of the variable $\nu$, which allows the formula to provide improved range and applicability, even for infinitely hard potentials and low n. One-sided potentials can be handled by considering only the even-n (odd parity) states (we number states $n=1,2,\dots$.

2.4. calculating Yukawa ${\mathcal{D}}_c$ Values Using Neoclassical Phase Space Method

A larger point we wish to make is that introducing these neoclassical corrections gives us permission to apply phase space quantization methods down to very low quantum numbers, and also to rapidly varying potentials. We will make extensive use of this in our analytical work later in this paper.

By way of illustration here we give a specific example of the Yukawa/Debye potential; the same method applies to other potentials. Here we set the energy $\mathcal{E}=0$ for the task of calculating critical screening lengths, in which we seek the values of $\mathcal{D}$ that just barely increment the number of bound states. We apply the Floor function to $\tilde{n}$ to do this quantization, but it is instructive to consider the continuously variable $\tilde{n}$ first.

Numerically calculating the phase space integral is straightforward, particularly in $z-space$. Later in this paper we give usefully-accurate closed-form approximate expressions for the values of these integrals. Here we do it by numerical integration. First, specify the potential for the given parameters $\mathcal{D}$ and l; find the classical turning points $z1$, $z2$ for those parameter values; numerically integrate $\sqrt{-2U}$ between them; and double the result. This gives $\mathcal{J}$; then $\tilde{n}=J/\left(2\pi \right)+\nu$, where $\nu$ is our continuously variable NC parameter, a deviation from the conventional semiclassical WKB/Maslov index.

As a concrete example, below is executable NC code to calculate ($\approx {10}^{-3}$ relative error) estimates of the number of bound states for the Yukawa Potential for specified l and $\mathcal{D}$, where $\mathcal{D}>e/2l(l+1)$ (lesser values of $\mathcal{D}$ can have no bound states and J will be complex-valued). This incorporates, in a primitive way, small dependencies of $\nu$ with n or l at low quantum numbers – estimates of $\tilde{n}$ would be rather poor using nominal semiclassical values. The second argument n has an effect only for $l=0$. The $l=0$ and $l>0$ cases are different in character because of the attractive singular potential for $l=0$. The parameters an, bn, al, bl are rough and certainly could be improved, as could the primitive functional form for $\nu$ itself; this code is just a simple example. The results can be compared with the extremely accurate PM tables given in [1]. For $l=0$, since $\nu$ depends on n, the correct value of $\mathcal{D}$ needed to obtain a particular n can be found by self-consistent iteration, or by root-finding, or tabulation + interpolation to obtain an inverse function. For $l>0$ no such iteration is needed, in this approximation. Similar code can be written for other potentials, for which the turning points can be determined by root-finding instead of analytically, as they are here. Another useful approximation can be found section 3.10.

nbound[l_,n_,D_]:=(

U=-(D Exp[-z]/z)+(l(l+1))/(2 z^2);

beta=(l(l+1))/(2 D);

W=ProductLog; z1=-W[0,-beta]; z2=-W[-1,-beta];

J=Quiet[2 NIntegrate[Sqrt[-2 U],{z,z1,z2}]];

If[l==0,nu0=0,nu0=1/2];

an=-0.034603; bn=-0.44314; al=-0.079887; bl=-0.83606;

deltanu[l1_,n1_]:=If[l1==0,an n1^bn,al l1^bl];

nb=J/(2 Pi)+nu0+deltanu[l,n];

Return[nb])

3. Results

3.1. Numerical Experiments Using PM Calculations on Modified ECSC Potentials

The contrast between the Yukawa and ECSC plots in Figures 1 and 2 is dramatic, and the multi-well structure that is present in the ECSC potential is an obvious hypothetical candidate to explain it. We first wish to determine if this is true, then determine a mechanism.

The Phase Method that was used to calculate critical binding parameters to 60 digit accuracy in our previous paper is sufficiently simple and robust that it easily can be (and was) used to perform numerical experiments on modified ECSC potentials. One question is whether the secondary ECSC potential well between $z=r/\mathcal{D}=(3\pi /2,5\pi /2)$ is essential for the observed sawtooth structure shown in Figure 2.

To test this, we did a simple numerical experiment in which the secondary well is filled-in (i.e. set to zero). This experiment used the same code and ranges that were used for calculating the accurate ${\mathcal{D}}_c=1/{\mu }_c$ tables in our previous paper, but with the modified potential. We see in Figure 3 that filling-in of the secondary well totally eliminates the sawtooth structure over the range $\mathcal{D}\le {10}^5$ au. Both $l=0$ and $l=1$ calculations are shown.

In another set of numerical experiments using the PM, the depth of the secondary well was multiplied by a “well scale factor” of wsf=1, $3/2$, and 2. The resulting sawtooth structure is shown in Figure 4. It is seen that the period scales approximately as $1/\sqrt{\mathrm{wsf}}$, which is explained below in terms of ratios of phase space integrals of primary and second wells.

In the third set of numerical experiments using the PM, the barrier between the primary and second wells was scaled by a “barrier scale factor” of bsf=${\displaystyle \frac{1}{4}}$, ${\displaystyle \frac{1}{2}}$, 1, 2, and 4. The results are shown in Figure 5. It is seen that the height of the barrier has little effect on the period, but reducing its height advances the position of first onset, and increasing it delays the onset.

From these numerical experiments it is clear that the second well is key to explaining the observed sawtooth structure over the indicated range. It will be shown below that tertiary and higher wells also contribute at much larger values of $\mathcal{D}$.

3.2. Neoclassical (NC) Explanation for Sawtooth Structure Observed for ECSC Potential

The previous subsection clearly shows the role of the second well in generating the sawtooth structure. We will analyze this structure from the point of view of a “neoclassical” (modified semiclassical, with variable $\nu$) approach described in subsection 2.2, in which the usually-neglected higher order terms in the exact infinite WKB series are approximated by an additive term: a constant, or a smoothly decreasing function of n or l. Here, for simplicity we will focus on the $l=0$ states, but the same kind of analysis applies to the higher l states, but with slightly greater complexity. Examples of this will be shown in subsection 3.10 in connection with explaining the linearity dependence of the total number of bound states on screening (Debye) length $\mathcal{D}$. As usual we will use atomic-like units in which ℏ and the reduced mass of the two-body system are both equal to 1.

For our purposes here, the phase space integral is the integral of momentum over position for a full cycle of classical motion in a bound state: $\mathcal{J}\left(\mathcal{E}\right)=2{\int }_{r_1\left(\mathcal{E}\right)}^{r_2\left(\mathcal{E}\right)}\sqrt{2(\mathcal{E}-U(r\left)\right)}dr$, where $r_1\left(\mathcal{E}\right)$ and $r_2\left(\mathcal{E}\right)$ are the classical turning points at energy $\mathcal{E}$, i.e. those points at which $\mathcal{E}=U\left(r\right)$, and $U\left(r\right)$ is the effective potential including the centrifugal potential if $l>0$.

We are interested in critical binding of states at $\mathcal{E}=0$, and in this subsection we confine our attention to $l=0$, which makes the turning points independent of $\mathcal{D}$ and l. Changing variables to $z=\mu r=r/\mathcal{D}$, we have for the first well of the ECSC potential: ${\mathcal{J}}_1=j_1\sqrt{\mathcal{D}}$ where $j_1=2{\int }_0^{\pi /2}\sqrt{2exp(-z)cos\left(z\right)/z}dz\approx 5.03251$; for the second well ${\mathcal{J}}_2=j_2\sqrt{\mathcal{D}}$ where $j_2=2{\int }_{3\pi /2}^{5\pi /2}\sqrt{2exp(-z)cos\left(z\right)/z}dz\approx 0.129588$. Crucially, for $l=0$, both phase space integrals scale exactly as $\sqrt{\mathcal{D}}$, so that dependence divides out. For $l>0$ the scaling is not exact, because an l-dependent term in the integrand prevents dividing out $\sqrt{\mathcal{D}}$, and because the classical turning points (limits of integration) also depend on $\mathcal{D}$ and l Yet the relation still holds approximately. In subsection 3.10 the interplay between the dependencies on $\mathcal{D}$ and l for $l>0$ is analyzed in detail.

The quantization condition for binding of a state within each (separate) well is of the form $n_1=\lfloor {\tilde{n}}_1\rfloor =\lfloor {\nu }_1+J_1/\left(2\pi \right)\rfloor$ and similarly for well 2. The floor function $\lfloor x\rfloor$ gives the largest integer that is $\le x$. Later we will it useful to consider the continuous-valued variable $\tilde{n}$, before the floor function is applied to it. The parameters ${\nu }_1$ and ${\nu }_2$ differ slightly from the well-known canonical values derived from topological considerations (yet still are based on the approximate, truncated WKB expansion), which are 0 and $1/2$ respectively for well 1 (with singularity) and well 2 (no singularity or hard walls). Here we choose values of ${\nu }_1=1/18$ and ${\nu }_2=6/10$. The choices of these values are representative, and were selected to give rough agreement with a few low-n and low l values, PM values, is not critical to our overall conclusions.

The combined state count for the two wells taken separately is then $\lfloor {\tilde{n}}_1\rfloor +\lfloor {\tilde{n}}_2\rfloor$. It might be argued that it it is the combined phase space area that should be quantized: $\lfloor {\tilde{n}}_1+{\tilde{n}}_2\rfloor$ but we find the results are more consistent with separate quantization of the wells. We note the inequality $\lfloor {\tilde{n}}_1\rfloor +\lfloor {\tilde{n}}_2\rfloor \le \lfloor {\tilde{n}}_1+{\tilde{n}}_2\rfloor$, so the scenario with two separately quantized wells has equal or fewer bound states. The combined state count is then $({\nu }_1+\sqrt{\mathcal{D}}{j}_1/\left(2\pi \right))+({\nu }_2+\sqrt{\mathcal{D}}{j}_2/\left(2\pi \right))$. As $\mathcal{D}$ increases, the main well accumulates bound states, and at a critical level, the second well start to bind a state. The relative rates at which they bind states is given by the ratio of ${\mathcal{J}}_2/{\mathcal{J}}_1=j_2/j_1\approx 38.8$. Crucially the common factor of $\sqrt{\mathcal{D}}$ divides out. In the sawtooth diagram the abscissa is the total state count $n_1$ +$n_2$. It’s analogous to two clocks, one ticking 39 times until the second clock tocks once, and the process repeats. The total number of states accumulated in a ticktock cycle is close to 40, which is that which is observed in the sawtooth diagram for $l=0$. The values of ${\nu }_1$ and ${\nu }_2$ affect the initial onset, but not the period. We call this the "tick-tock" mechanism.

Using this model a calculation of the curve corresponding to the $l=0$ curve in Figure 2 was generated, and it is shown in Figure 6. To do so, the function $\lfloor {\tilde{n}}_1\rfloor +\lfloor {\tilde{n}}_2\rfloor$ was evaluated on a grid of $\Delta \mathcal{D}=0.01$ au up to $\mathcal{D}={10}^5$. The transition points were detected in one line of code (First@GroupBy in Mathematica); the plot is given below. We see that it compares very well with corresponding $l=0$ curve in Figure 2 which is based on 30-digit accuracy PM calculations. In particular, we see several attributes are well-reproduced: the sawtooth structure; its nominal period of 40; the locations of the minima $n=17,57,97,137,176,216$; and the values of the curve themselves. The precise locations of the minima are affected somewhat by the approximate “neoclassical” $\nu$ parameters; the tick-tock period of “40” is simply a pure number with no adjustable parameters: $1+j_2/j_1\approx 39.83476128\dots$.

Calculating the curves for $l>0$ in the same manner is only slightly more complex, because of the dependence of the turning points on l and $\mathcal{D}$, but the essentials are the same. This sort of calculation is done below for comparison of Hulthén, Yukawa, and ECSC potentials for arbitrary l, for the purpose of explaining the strangely linear dependence of the total number of bound states vs $\mathcal{D}$.

At this point we consider the first mystery to be solved, but there is more to investigate about it.

3.3. Subsidiary Wells of Higher Order

This tick-tock mechanism also extends to the infinite number of progressively shallower wells at larger z that are characteristic of the ECSC potential, and should give rise it an infinite number of “tockitos” for sufficiently large $\mathcal{D}$. The question is where they occur. It is not difficult to estimate (and indeed numerically evaluate) these phase space action integrals in z-space, with $z=r/\mathcal{D}$. We number the higher wells as $m=1,2,3,..$. The $m_{\mathrm{th}}$ well is localized between $2\pi m-\pi /2\le z\le 2\pi m+\pi /2$. Replacing the integrand over that interval by its value at its average z value $2\pi m$, and multiplying by its width $\Delta z=\pi$, with a correction factor c on the order of one to account for its shape (which is quite similar for all m), we get $\mathcal{J}\left(m\right)\approx c\sqrt{2\pi \mathcal{D}}{e}^{-\pi m}/\sqrt{m}$ with $c\approx 1.17$.

Then in our usual way we find $\tilde{n}-\nu =J\left(m\right)/\left(2\pi \right)$ where $\nu \approx 1/2$, or a simple neoclassical adjusted value $0.6$. Solving for the screening length $\mathcal{D}$ required to bind $\tilde{n}$ states, we have $\mathcal{D}\approx {(\tilde{n}-\nu )}^2\left(2\pi m\right)e^{2\pi m}/c^2$. Evaluating the $\mathcal{D}$ required to bind one state and using $\nu =0.6$ we get for $m=1,2,3,4$ respectively $3.9\ast {10}^2$au, $4.2\ast {10}^5$au, $3.4\ast {10}^8$au, $2.4\ast {10}^{11}$ au, or approximately 21 nm, 22 $\mu$m, 18 mm, and 13 meter. Furthermore the period of such sawtooth oscillations for the $m=2$ well would be about 1293, compared to 40 for the $m=1$ well. As a consequence it might be expected for such a single tick-tock glitch to be observed experimentally (e.g. in spectra) for the $m=2$ well, in addition to the expected sawtooth/tick-tock structure from the $m=1$ well. The range of $\mathcal{D}$ needed to observe the second and higher sawtooth glitches is quite extensive, however, as it varies as the square of the state number. In the next subsection we address this using the PM.

3.4. Sawtooth Structure from Higher Index Wells

The spatial locations of the wells for $l=0$ in $z=r/\mathcal{D}$-space occur at $(0,\pi /2)$, $(3\pi /2,5\pi /2)$, $(7\pi /2,9\pi /2)$, …, or more concisely $\max (0,(2\pi m-\pi /2,2\pi m+\pi /2\left)\right)$ with $m=0,1,..$. The primary well corresponds to $m=0$, the secondary well $m=1$, etc. For $l=0$ and $m>0$ the wells are uniformly spaced in z and their shapes are very similar to each other, but with very large ($⪆500$-fold) reduction in depth between successive wells. This regularity makes estimating the well depths and phase space integrals straightforward by scaling successive wells, and using a common shape factor c on the order of 1 to approximate the integral. In this way we obtain a simple approximate expression for the phase space integral $\mathcal{J}$ for any $m>0$, which agrees sufficiently well with numerical integrals. We use a NC $\nu$ parameter with a value near $1/2$ (specifically $\nu =0.6$) for the wells $m>0$; and for $m=0$, a value near zero (specifically $\nu =0.06$), owing to the potential singularity. These deviations from the canonical semiclassical values are small but important corrections that are essential to obtain usefully accurate results for low n. They could be improved but these are sufficient for our purposes here.

Setting $n\approx \lfloor J/\left(2\pi \right)+\nu \rfloor$ we obtain an estimate for ${\mathcal{D}}_c$ for any well $m>0$ and any of its bound states $n_m$: ${\mathcal{D}}_c=(2\pi m{e}^{2\pi m}{(n_m-\nu )}^2/c^2)$ where $\nu \approx 0.6$ and $c\approx 1.2$. This agrees roughly with more accurate numerical phase space integrals and precise PM values. For simplicity here we neglect small $n_m$-dependent corrections to $\nu$, at the expense of some accuracy.

Evaluating this function we see that the first bound state of the tertiary well ($n=1$, $m=2$, $D\approx 4.0\ast {10}^5$ au), is predicted to be between the $n=13$ ($D\approx 3.6\ast {10}^5$ au) and the $n=14$ ($D\approx 4.2\ast {10}^5$ au) levels of the secondary $m=1$ well, generating a single“tick-tock” glitch there for that series. Specifically, evaluating our simple expression above with $n=1$ predicts the first bound state of the tertiary $m=2$ well should occur at $\mathcal{D}\approx 4.0\ast {10}^5$. This would give a downward jump in the ${\mathcal{D}}_c/n^2$ curve around that location, but a rather slowly varying one, because of the large $\mathcal{D}$-scale.

We subsequently tested this prediction by doing a PM calculation for the ECSC potentials $l=0$ up to a maximum value of $\mathcal{D}={10}^6$. The result is shown in Figure 7. Agreement with the prediction is excellent. If observed experimentally and unexpected, this feature could cause some confusion. The second “tockito” of $n=2,m=2$ ($D\approx 5.0\ast {10}^6$ au) state is expected to occur at $\mathcal{D}$ values that are an order of magnitude greater, i.e. hundreds of micrometers scale. The ${\mathcal{D}}_c$ sequences for the successive wells start to overlap in the region $n\approx 11-13$. The corresponding values of $\mathcal{D}$ in such cases quickly get into the macroscopic length range but those might be relevant for some systems. This is shown in Figure 8.

3.5. Illustrative Tick-Tock Model: Asymmetric Quartic Double Well

It is instructive to consider a simple double well potential as a model of this tick-tock behavior, and to consider it from the points of view of the fully quantum PM method, and the neoclassical NC method.

An example of such a potential is $U\left(x\right)=x^4-6x^2+4x$, which is plotted in Figure 10; the right panel also shows a phase space plot evaluated at zero energy. If placed in one of the wells, a classical particle with a total energy less than the top of the barrier, would remain confined to that well. In reality, some degree of quantum tunneling will occur between the wells, and such tunneling will be significantly enhanced for energies near the top of the barrier. If the discrete levels for the two separate wells were to line up with each other, and extensive tunneling were to occur, this could induce substantial splitting between the energy levels, which might have other consequences. However, if the energy levels weren’t extremely close (which in general is unlikely) the two-well system behaves as a combination of two separate wells, but with a combined total state count.

For each well, and at a specified energy, the phase space integrals ${\mathcal{J}}_{1,2}\left(\mathcal{E}\right)$ can be computed as twice the integral of the momentum between the classical turning points. At $\mathcal{E}=0$ it corresponds to the areas contained within the left and right regions shown in Figure 9, which correspond to the primary and seconday wells.

We define a continuous valued $\tilde{n}=\mathcal{J}\left(\mathcal{E}\right)/\left(2\pi \right)+\nu$. The value of $\nu$ is to a degree affected by the amount of quantum tunneling. When $\tilde{n}$ exceeds an integer, a new bound state can form. To compare results between PM and NC, in Table 1 this expression is evaluated for both the left (${\tilde{n}}_1$) and right (${\tilde{n}}_2$) wells, at the energy if the PM-computed eigenvalue.

Using the Phase Method we computed the energy levels and wave functions for this potential. In Table 1 these energy eigenvalues are given in column 2. We can see in Figure 10 that at lower energies, up to three bound states are confined to the deeper left well, until a new state appears in the shallow right well, giving four states altogether. This is the analog of the first tick-tock state. A second state is spread across both wells at the slightly negative fifth eigenvalue. The sixth and higher states are continuum states are spread across both wells.

The PM wavefunctions, and in particular, their squares, the corresponding probabilities, reinforce this picture. The states $n=1,2,3$ reside entirely within the left well, but the probability in the fourth level for the first time is concentrated in the right well. The fifth energy level which is only slightly bound has probability distributed across both wells, as does the sixth, which is at positive energy.

Figure 10. PM-calculated probability distributions (squared wavefunctions) for the asymmetric quartic double-well potential. The probability (squared wavefunction) is shown, superimposed on the potential to show where it is spatially distributed.

Figure 10. PM-calculated probability distributions (squared wavefunctions) for the asymmetric quartic double-well potential. The probability (squared wavefunction) is shown, superimposed on the potential to show where it is spatially distributed.

Table 1. Comparison of fully quantum PM calculations and neoclassical phase space quantization for the two wells. The first and last columns indicate the quantum number; the second the eigenvalue. ${\tilde{n}}_{1,2}$ are $\mathcal{J}/\left(2\pi \right)+\nu$ for the first and second well respectively. ${\nu }_1=1/2+0.01$, and ${\nu }_2=1/2+0.04$, small but significant deviations from canonical values, that are related to the extent of tunneling. $\lfloor x\rfloor$ represents the floor function, i.e. the largest integer $\le x$.

Table 1. Comparison of fully quantum PM calculations and neoclassical phase space quantization for the two wells. The first and last columns indicate the quantum number; the second the eigenvalue. ${\tilde{n}}_{1,2}$ are $\mathcal{J}/\left(2\pi \right)+\nu$ for the first and second well respectively. ${\nu }_1=1/2+0.01$, and ${\nu }_2=1/2+0.04$, small but significant deviations from canonical values, that are related to the extent of tunneling. $\lfloor x\rfloor$ represents the floor function, i.e. the largest integer $\le x$.

3.6. Solution to the n* Linearity Mystery for ECSC Potential

The classic paper by Rogers et al in 1970 calculated critical screening lengths for the lowest 45 bound states of the Yukawa/Debye potential, as well as eigenvalues and other things. They found two approximately linear relationships between the screening length $\mathcal{D}$ and two quantities: ${(g\ast )}^2$, and $n\ast$, where $g\ast$ is quantum number of the highest bound $l=0$ state, and $n\ast$, the total number of bound states of all angular momenta. Our own calculations of critical screening lengths to 60 digits and over a larger range of $\mathcal{D}$ confirmed and extended the results of Rogers et al. Figure 1 from our first paper [1] is reprised here as Figure 11. The rectangular region marked in the lower left corner indicates the region covered by Rogers et al [2]. Our calculations independently reproduce their results precisely over that region, and also substantially extend the range, over which the linearity persists. Both represent all of the bound states for their specified $\mathcal{D}$ ranges.

The first of the linear relations, the one involving $g\ast$, is easily explained, as others have done, and we did in our recent Atoms paper, using a simple phase space quantization argument, which indicates the slope of that linear plot is $4/\pi$, which agrees well with least squares fits to our PM-calculated critical values. Specifically this says that $\mathcal{J}=2{\int }_0^{\infty }\sqrt{2e^{-z}/z}dz=4\sqrt{\pi }$, and $g\ast =\lfloor \mathcal{J}/\left(2\pi \right)\rfloor$, and ${(g\ast )}^2\approx 4/\pi \mathcal{D}$.

Similarly, Lam and Varshni’s 1972 paper [3] found linear relationships (of course with different slopes than Yukawa/Debye) for the ECSC potential. The $g\ast$ plot refers to $l=0$ states. From phase space argument similar to the Yukawa case we find (using numerical integration) $\mathcal{J}=2\sqrt{\mathcal{D}}{\int }_0^{\pi /2}\sqrt{2e^{-z}cosz/z}dz\approx 5.0325\sqrt{\mathcal{D}}$, and $g\ast =J/\left(2\pi \right)$, giving a slope of ${(g\ast )}^2$ vs $\mathcal{D}$ of $0.641519$. Both Rogers et al, and Lam and Varshni allowed for nonzero intercept values, and we follow suit for comparison. Fits to our PM-calculated ${\mathcal{D}}_c$ values give ${(g\ast )}^2\approx 0.128296+0.641436\mathcal{D}$ with confidence intervals of $(0.118623,0.137969)$, $(0.641371,0.641501)$ respectively. We consider this very good agreement with the semiclassical phase space quantization estimate. Lam and Varshni’s fits give a slope estimate of $0.6078$, but they take pains to point out these values are underestimates of the slope and intercept, owing to the under-sampling of the data points. For this reason we consider their values and ours to be in sufficient agreement, although our contemporary values are much more accurate.

Yet what has been missing for more than a half-century, to our knowledge, is any explanation for the second linear relation involving the total number of bound states $n\ast$. The striking high degree of linearity is mysterious, because the underlying dependence is $n\propto {\mathcal{D}}^{1/2}$. In this subsection we present a solution to this mystery.

In Lam and Varshni’s $n\ast$ plot for the ECSC potential, in contrast with Rogers et al’s plot (and our Figure 12) for the Yukawa potential, the points do not lie as consistently on the best fit line for the $n\ast$ plot. Their state counts are slightly under-sampled, as pointed out by the authors. Here we plot our ECSC potential results in the same manner which as shown in Figure 1. These curves represent all the bound states with no undersampling issues. Lam and Varshni’s $\mathcal{D}$-range stops right before the sawtooth structure emerges in our plots.

In contrast to the Yukawa potential in Figure 11, In Lam and Varshni’s [3] Figure 2, the points on the $n\ast$ line are bunched-up into distinct groups, with gaps in-between. Some tendency to downward curvature of each group also is visible. Our Figure 12, based on our PM calculations, with range selected for comparison with [3]. shows all of the bound states up to the maximum $\mathcal{D}$ value, and it is quite revealing. First we note that only a single state first appears on the line, then a doublet, then a triplet, and so on. Inspecting the nature of the states (we computed them, so we know what they are) we find that they correspond to states of the same principal quantum number, and in the usual atomic ordering: (1s), (2s, 2p), (3s, 3p, 3d), (4s, 4p, 4d, 4f) etc. In the ECSC potential no level-crossings are observed over the entire range that was calculated up to $\mathcal{D}={10}^5$ au. The total number of states below a selected s-state $n\ast$ (which are easy to locate as above using $g\ast =\lfloor j/\left(2\pi \right)\sqrt{\mathcal{D}}\rfloor$) with $j\approx 5.0325$, is then just the arithmetic sumof the number of states in each complete group below it, which evaluates to a state-count $g\ast (g\ast -1)/2+1$, or $g{\ast }^2\propto \mathcal{D}$ for moderate $g\ast$. This is encouraging.

Although the dominant term varies as $\mathcal{D}$, the expression also includes a term proportional to $\sqrt{D}$, which reduces the slope slightly; yet the curve remains sufficiently close to linear. The quantum graininess adds a measure of apparent noise (which it is not). Fitting the arithmetic sum estimate derived above over the range $\mathcal{D}=(0,360)$ to a simple proportionality (with zero intercept) gives a function $0.300\mathcal{D}$, importantly, with a residual between fit and curve that does not deviate by more than 1 (bound state) unit, and so it is within the quantization graininess of the computed data points. The determined slope also is in reasonably good agreement with the slope $.3015$ of a linear least square fits to our own high precision PM-computed values, and also is consistent with Lam and Varshni’s self-proclaimed lower bound $0.2919$.

We conclude the linear n* plot for the ECSC potential is explained using basic phase space quantization and considerations of the arithmetic sum state accumulation that is appropriate in this case. This was relatively simple because the states of different principal quantum number $\overline{n}=n+l$ did not overlap – there were no level crossings, a fact that emerged from the accurate PM calculations. This is not the case for Yukawa and Hulthén potentials, however, so a more general approach is needed, which we investigate next.

Figure 11. Plot of $g{\ast }^2$ (square of maximum number of bound $l=0$ states; and $n\ast$ (total number of bound states of all l) as a function of screening length $\mathcal{D}$, for the Yukawa potential. The inset corresponds to the range covered by Rogers et al [2] Figure 5. g* is the higher bound $l=0$ state; n* is the maximum number of all bound states for any l.

Figure 11. Plot of $g{\ast }^2$ (square of maximum number of bound $l=0$ states; and $n\ast$ (total number of bound states of all l) as a function of screening length $\mathcal{D}$, for the Yukawa potential. The inset corresponds to the range covered by Rogers et al [2] Figure 5. g* is the higher bound $l=0$ state; n* is the maximum number of all bound states for any l.

Figure 12. Plot of g*² (square of maximum number of bound $l=0$ states; and n* (total number of bound states of all l) as a function of screening length $\mathcal{D}$, for the ECSC potential. The range is chosen for comparison with Lam and Varshni [3] Figure 2. The first tick-tock glitch is slightly above this range. Calculations were done by the Phase Method as described previously. Note the 1-2-3-4… state binding pattern observed along the n* line. These correspond to (1s), (2s, 2p), (3s, 3p, 3d) etc. A similar pattern is not found for the other potentials because the states of different principal quantum number overlap for them.

Figure 12. Plot of g*² (square of maximum number of bound $l=0$ states; and n* (total number of bound states of all l) as a function of screening length $\mathcal{D}$, for the ECSC potential. The range is chosen for comparison with Lam and Varshni [3] Figure 2. The first tick-tock glitch is slightly above this range. Calculations were done by the Phase Method as described previously. Note the 1-2-3-4… state binding pattern observed along the n* line. These correspond to (1s), (2s, 2p), (3s, 3p, 3d) etc. A similar pattern is not found for the other potentials because the states of different principal quantum number overlap for them.

Figure 13. Plot of g*² (square of maximum number of bound $l=0$ states; and n* (total number of bound states of all l) as a function of screening length $\mathcal{D}$, for the Hulthén potential. The range is limited so that no states for $l>12$ fall within the data range, as they are not represented in the dataset. Calculations were done by the Phase Method as described previously. The same kind of linear behavior is found as for Yukawa and ECSC potentials.

Figure 13. Plot of g*² (square of maximum number of bound $l=0$ states; and n* (total number of bound states of all l) as a function of screening length $\mathcal{D}$, for the Hulthén potential. The range is limited so that no states for $l>12$ fall within the data range, as they are not represented in the dataset. Calculations were done by the Phase Method as described previously. The same kind of linear behavior is found as for Yukawa and ECSC potentials.

3.7. Geometric Solution to the $n^{\ast }$ vs $\mathcal{D}$ Linear Relation Mystery

The ${\mathcal{D}}_c$ values for Yukawa and Hulthén potentials do not show the convenient grouping by principal quantum number that is seen for the ECSC potential. We can “unfold” the states according to both n and l by plotting the ${\mathcal{D}}_c$ values for each state vs n, or l, with all states being represented. These show a startling regularity, as can be seen in Figure 14 and Figure 15: the density of states in this representation is seen to be nearly uniform. This solves the mystery: the total number of states below a specified threshold $\sqrt{D_c}$ value is proportional to the area of the triangular region below it. This area grows in proportion to the square of the threshold, i.e. ${\left(\sqrt{{\mathcal{D}}_c}\right)}^2={\mathcal{D}}_c$. This begs the question as to why the potentials have that nearly uniform density, however, a question we address next.

In the previous section we found that, using accurate and complete tables of ${\mathcal{D}}_c$ calculated by the Phase Method, there is a nearly uniform density of states when $\sqrt{{\mathcal{D}}_c}$ is plotted vs n, or l. The reason for this mysterious result also requires an explanation, which we pursue in this subsection.

In our previous paper [1], we showed in passing that a contour plot of $\sqrt{{\mathcal{D}}_c}$ vs n and $l+1$ for the Yukawa potential exhibited a striking linearity in both variables. Here we extend this analysis to the ECSC and Hulthén potentials, using the tables from [1]. We find very similar results for all three potentials, as shown in Figure 16. The nearly parallel and equally spaced contours for all three potentials indicate a uniform linear gradient of $\sqrt{{\mathcal{D}}_c}$ vs n and l (or $l+1$), with an additive constant. It is not exact – the residuals are not zero, as shown in the right panels of those figures.

Bivariate linear least squares fits were performed, with best fit parameters and confidence intervals given in Table 2 and Table 3. The ratios of the coefficients are given in Table 4 for comparison with calculations to be given later.

We note that this is contrary to the strict one-to-one relation between n and l found for the Coulomb potential, in which the energy levels depend only on the sum $\overline{n}=n+l$, where $\overline{n}$ is known as the principal quantum number. Incrementing n and decrementing l by one give the same energy. This is a consequence of the famous dynamical SO(4) symmetry of the $1/r$ potential, which does not apply to screened Coulomb potentials. They do, however, show a broken-symmetry version of that 1:1 linear relation.

We predict that other potentials in this class (i.e those that satisfy the scaling law between screening length and coupling constant) will exhibit the same near-linear bivariate relation with n and l. The fundamental reason for this is explained below.

Table 2. Best least squares fits to PM-calculated $\sqrt{{\mathcal{D}}_c}$ linear combination of $l+1$ and n plus a constant.

Table 2. Best least squares fits to PM-calculated $\sqrt{{\mathcal{D}}_c}$ linear combination of $l+1$ and n plus a constant.

Table 3. Confidence Intervals (95%) for the fit parameters above.

Table 3. Confidence Intervals (95%) for the fit parameters above.

Table 4. Ratios between fit coefficients. See text for explanation.

Table 4. Ratios between fit coefficients. See text for explanation.

3.8. Perturbation Theory Calculation of Linear Correlation Contour Plot for the Yukawa Potential

First order perturbation theory also was used to calculate estimates of the ${\mathcal{D}}_c$ values as a function of l and n for the Yukawa potential. Exact Coulomb eigenfunctions and eigenvalues were used for the unperturbed states, and the full difference between the Yukawa potential and Coulomb potential $1-e^{\mu r})/r$ was used as a perturbation Hamiltonian – the centrifugal potential terms cancel exactly. Despite the fact that the unperturbed Coulomb states of a given principal quantum number but different l are degenerate, any nonzero screening parameter lifts that degeneracy, so the use of non-degenerate perturbation theory is still appropriate.

The first order corrections to the eigenvalues are simply given by the expectation value of the perturbation Hamiltonian above within each unperturbed eigenstate. All of the integrals were done symbolically in Mathematica. The perturbation increases the energy of the perturbed eigenstates for each $n,l$, and the value of $\mathcal{D}$ required to raise the perturbed eigenvalues to reach zero energy was solved for, and is the estimate of ${\mathcal{D}}_c$ for each state. These are the values shown in the contour plot below.

This first order perturbation theory is approximate, but it is strikingly similar to the others calculated by the far more accurate Phase Method calculations. Our conclusion is this approximately linear relation between $\sqrt{D}$, n, and l is intrinsic to the nature of this potential. The PM calculations shown above, and the NC calculations shown below indicate this is a shared feature of this class of potentials, which we will show is related to their scaling symmetry.

Figure 17. $\sqrt{{\mathcal{D}}_c}$ vs $n_{\overline{n}}$ and $l+1$ for Yukawa potential, as calculated by first order perturbation theory. The unperturbed potential is Coulomb, and the difference between the Yukawa and Coulomb potentials was used as a perturbation. The displayed ${\mathcal{D}}_c$ values are those values required for the perturbed eigenvalues to be shifted up to zero energy. Computations were done analytically in Mathematica, but are rough in comparison with the PM-calculated values in Figure 9. The key feature of interest is the striking linear relation between ${\mathcal{D}}_c^{1/2}$ and n and l, which is maintained despite the rough approximate nature of this first order PT calculation.

Figure 17. $\sqrt{{\mathcal{D}}_c}$ vs $n_{\overline{n}}$ and $l+1$ for Yukawa potential, as calculated by first order perturbation theory. The unperturbed potential is Coulomb, and the difference between the Yukawa and Coulomb potentials was used as a perturbation. The displayed ${\mathcal{D}}_c$ values are those values required for the perturbed eigenvalues to be shifted up to zero energy. Computations were done analytically in Mathematica, but are rough in comparison with the PM-calculated values in Figure 9. The key feature of interest is the striking linear relation between ${\mathcal{D}}_c^{1/2}$ and n and l, which is maintained despite the rough approximate nature of this first order PT calculation.

3.9. Explanation of Approximate Linear Dependence of q, n, l by NC Phase Space Arguments

The striking near-linearity shown in the contour plots above requires explanation, and again phase space arguments provide quantitative insight. The applicability of the neoclassical approach to small n levels extends its reach to problems for which traditional semiclassical analysis might be considered inappropriate.

We found an approximate linear relationship between all three quantities: $\sqrt{\mathcal{D}}$; the state number n; and the angular momentum l. In the asymptotic limit $\mathcal{D}\to \infty$ we also know the Yukawa Potential becomes the Coulomb potential, which has the known hidden SO(4) symmetry which explains the degeneracy of its energy levels of a given principal quantum number $\overline{n}=n+l$. It is simply the sum of n and l that determines the energy, so there is an inherent symmetry between n and l: increasing one of them by 1 is equivalent to decreasing the other by 1. The near one-to-one alignment of the contours in $n,l$ space in the contour plots presumably reflects the breaking of this symmetry that occurs at finite screening lengths.

We can compute estimates of this $n,l$ correlation for the Yukawa potential using phase space arguments. Analysis of the Hulthén and ECSC potentials is similar. The effective potential in r-space is

$$U\left(r\right)={\displaystyle \frac{-e^{-\mu r}}{r}}+{\displaystyle \frac{l(l+1)}{2r^2}}.$$

As is well-known, and also shown in our previous paper, a characteristic feature of this class of potentials is that scaling the screening length in the exponential by a factor is exactly equivalent to multiplying the exponential term by the same factor. First, we will consider the continuous-valued $\tilde{n}=\mathcal{J}/\left(2\pi \right)+\nu$, taking $\nu$ constant, and restrict our attention to $l>0$ for comparison with the contour plots. Changing variables to $z=r/D$, and defining $q^2=\mathcal{D}$ and ${\lambda }^2=l(l+1)$, our phase space integral becomes

$$\mathcal{J}=2{\int }_{z_1\left(\beta \right)}^{z_2\left(\beta \right)}p\left(z\right)dz=2\sqrt{2}{\int }_{z_1\left(\beta \right)}^{z_2\left(\beta \right)}\sqrt{{\displaystyle \frac{q^2{e}^{-z}}{z}}-{\displaystyle \frac{{\lambda }^2}{2z^2}}}dz,$$

and $z_1$, $z_2$ are the classical turning points, which are the locations in which the momentum p vanishes.

These evidently depend on a specific combination of $\mathcal{D}=q^2$ and ${\lambda }^2=l(l+1)$. Defining $\beta ={\lambda }^2/\left(2q^2\right)$, as in our previous paper, we have $z_1=-W(0,-\beta )$ and $z_2=-W(-1,-\beta )$ where W is the Lambert W function (ProductLog in Mathematica). An upper bound on the value of $\beta$ that might bind a state is ${\beta }_{\max }=1/e$.

The approximate symmetry between $q=\sqrt{\mathcal{D}}$ and $\lambda =\sqrt{l(l+1)}$ is evident by inspection, and it follows from the scaling symmetry mentioned above, because that is the origin of the coupling constant coefficient $\mathcal{D}=q^2$ of the exponential term. With the change of variables to z-space, the z-values between $z_1$ and $z_2$ are generally on the order of 1, and so the coefficients of of the $q^2$ and ${\lambda }^2$ terms are similar in magnitude. The function $Y\left(z\right)=z{e}^{-z}$ is the ratio of the two weight coefficients, and it is the roots of $Y\left(z\right)=\beta$ that give the turning point solutions. The peak of $Y\left(z\right)$ is at $z=1$, and its maximum is $Y\left(1\right)=1/e={\beta }_{\max }$.

We can make this more quantitative with an explicit computation. Our contour plots in Figure 15 and Figure 16 express q as a function of n and l (or $l+1$, a difference of little consequence). The neoclassical/semiclassical phase space integral gives us $\tilde{n}$ as a function of q and $\lambda$. We can differentiate under the integral with respect to parameters q and $\lambda$ (which we will take as continuous). By Leibnitz’s formula, the limits of the integral, which also depend on q and $\lambda$ also must be differentiated, but in the case of phase space integrals, those contributions vanish, because the momentum is zero at the turning points: $p\left(z_1\right)=p\left(z_2\right)=0$, identifying the momentum at zero energy $p\left(z\right)$ as $\sqrt{2}$ times the square root in this integrand. Therefore only differentiation with respect to the integrand itself contributes.

With this useful fact in hand, we find:

$${\displaystyle \frac{\partial \tilde{n}}{\partial q}}=1/\pi {\int }_{z_1\left(\beta \right)}^{z_2\left(\beta \right)}{\displaystyle \frac{\partial p\left(z\right)}{\partial q}}dz=2q/\pi {\int }_{z_1\left(\beta \right)}^{z_2\left(\beta \right)}{\displaystyle \frac{e^{-z}/z}{p\left(z\right)}}dz.$$

(4)

Similarly, differentiating with respect to $\lambda$ gives

$${\displaystyle \frac{\partial \tilde{n}}{\partial \lambda }}=1/\pi {\int }_{z_1\left(\beta \right)}^{z_2\left(\beta \right)}{\displaystyle \frac{\partial p\left(z\right)}{\partial \lambda }}dz=-\lambda /\pi {\int }_{z_1\left(\beta \right)}^{z_2\left(\beta \right)}{\displaystyle \frac{1/z^2}{p\left(z\right)}}dz.$$

(5)

These expressions are respectively proportional to q and $\lambda$, and the integrals are pure numbers that depend on $\beta$ via the classical turning points and $p\left(z\right)$. $\beta$ ranges from zero to $1/e$ for the Yukawa Potential. $1/p\left(z\right)$ diverges at the turning points but the integrated weight is not very large there. The functions $e^{-z}/z$ and $1/z^2$ also have similar z-dependence. The result is that both integrals are very similar in magnitude.

We want to compare with the contour plots that express q as a function of n and l. Here we have $\tilde{n}$ as a function of q and $\lambda$. To compare, we have $d\tilde{n}={\displaystyle \frac{\partial \tilde{n}}{\partial q}}dq+{\displaystyle \frac{\partial \tilde{n}}{\partial \lambda }}d\lambda$. For fixed $\tilde{n}$ we then have ${\displaystyle \frac{dq}{d\lambda }}=-{\displaystyle \frac{\partial \tilde{n}}{\partial \lambda }}/{\displaystyle \frac{\partial \tilde{n}}{\partial q}}$ and similarly ${\displaystyle \frac{dq}{dl}}=-{\displaystyle \frac{\partial \tilde{n}}{\partial l}}/{\displaystyle \frac{\partial \tilde{n}}{\partial q}}$.

Using differentiation under the integral as above, the ratio on the right side of this equation is readily computed numerically for a given q and l, or the composite variable $\beta$. We have evaluated this numerically over the (integer-valued) feasible pairs q, l i.e. those for which $0<\beta <1/e$, and examined the results from $l=1$–30 and $q=1$–60.

The mean and standard deviation of ${\displaystyle \frac{dq}{dn}}$ are $0.862\pm .017$ and ${\displaystyle \frac{dq}{dl}}$ are $1.132\pm .033$. These compare well with the coefficients in the linear fit of $0.853$ (dq/dn) and $1.153$ (dq/dl) obtained by linear bivariate fits over the narrower range $l+1=(1,15)$, $n=(1,15)$ that were fit to the PM data in the Yukawa Potential contour plot in subsection 3.7. We find this to be a satisfactory numerical confirmation of the approximate mutual linearity of n, q, l.

The feasible pairs of $q,l$ are those for which $\beta <{\beta }_{\max }=1/e$, i.e. $l(l+1)/\left(2q^2\right)<1/e$. This condition can be rewritten ${(2l+1)}^2-(8/e)q^2<1$, which is a region bounded by the indicated hyperbola in q, l, as shown in Figure 18. The apparently straight line relation between $l_{\max }$ and $q=\sqrt{D}$ is actually an asymptote of this hyperbola, with deviations from linearity at low q and l. The constant $1/4$ is small enough the relation is effectively a straight line. The slope of this line ${\displaystyle \frac{dq}{dl}}$ asymptotes to $1/\sqrt{2{\beta }_{\max }}$, which is in good agreement with the ratios of coefficients of the linear fits in Table 2. These imply that the deviations of the gradient from the (1,1) direction in the contour plots of Figure 16 are characteristic of the form of the potential, and are principally related to the value of $2{\beta }_{\max }$ for each potential.

We conclude that the linearity of the Rogers et al, and Lam and Varshni plots, as well as our extended range versions of them, and our contour plots, all are explainable in terms of a NC phase space analysis. Differentiating under the phase space integral $\mathcal{J}$, predicts near-linear dependence between $\sqrt{\mathcal{D}}$, l, and n. This is further seen to be a consequence of the approximate symmetry between $q^2$ and ${\lambda }^2=l(l+1)$ in the phase space integral $\mathcal{J}$, which is a consequence of the scaling relation between coupling constant and screening length in this entire class of potentials. The relation also reflects the fact that the Yukawa potential breaks the SO(4) symmetry of the Coulomb potential, in which n and l affect spectra in the combination $n+l$. We expect a similar analysis and conclusion holds for other potentials in this class (e.g. Woods-Saxon), with possible in multiwell potentials such ECSC and generalized ECSC.

3.10. Y-Function Analysis of q, n, l Near-Linear Relation

The q vs $\lambda$ correlation is even more intimate than suggested above, which becomes clear when written in terms of the composite variable $\beta =l(l+1)/\left(2\mathcal{D}\right)\equiv {\lambda }^2/\left(2q^2\right)$, and the function $Y\left(z\right)=z{e}^{-z}$ (for Yukawa potential – $Y\left(z\right)$ has a different form than this for its sibling potentials); as usual $z=r/\mathcal{D}$. $\beta$ is determined for a specified screening length $\mathcal{D}$ and angular momentum l. Only values between $0\le \beta \le 1/e$ are meaningful for the Yukawa potential, and the two roots of the equation $Y\left(z\right)=\beta$ determine the classical turning points which are the limits of integration in the phase space integral. In this case they are soluble analytically as $z_1\left(\beta \right)=-W(0,-\beta )$ and $z_2\left(\beta \right)=-W(-1,-\beta )$ where W is the Lambert W function, ProductLog in Mathematica.

As in [1] we can write the phase space integral $\mathcal{J}=qj\left(\beta \right)$ where $j\left(\beta \right)=2\sqrt{2}{\int }_{z_1\left(\beta \right)}^{z_2\left(\beta \right)}\sqrt{Y\left(z\right)-\beta }{\displaystyle \frac{dz}{z}}$. This tells us that $\mathcal{J}$ is proportional to $q=\sqrt{D}$ but with a residual dependence on q via the composite variable $\beta$. With $2\pi (\tilde{n}-\nu )=J$, and approximating $\nu$ as a constant, and recognizing that differentiating with respect to the limits gives no contribution as explained above, we have $2\pi {\displaystyle \frac{\partial \tilde{n}}{\partial q}}=j\left(\beta \right)-2\beta {\displaystyle \frac{dj\left(\beta \right)}{d\beta }}$. For strict linearity between $\tilde{n}$ and q to occur, this would require $j\left(\beta \right)-2\beta {\displaystyle \frac{dj\left(\beta \right)}{d\beta }}=\mathrm{constant}$. Solving this simple differential equation we get $j=j\left(0\right)-j_c{\beta }^{1/2}$. This is the form of $j\left(\beta \right)$ that would give strict proportionality between $\tilde{n}$ and q. What we obtain numerically approximates this fairly closely, but not exactly.

Evaluating the integral $j\left(\beta \right)$ numerically, tabulating the numerical integral on a fine $\beta$ grid, we find $j\left(0\right)\approx 7.08982$, which very well matches the exact asymptotic (large n) value $4\sqrt{\pi }$. Then, fitting the computed $j\left(0\right)-j\left(\beta \right)$ (using Mathematica’s NonlinearModelFit) to a power law we obtain an excellent fit: $j\left(0\right)-j\left(\beta \right)=j_c{\beta }^{j_p}$ with $j_c=12.0892\pm 0.0018$ and $j_p=0.53446\pm .00010$, with a coefficient of determination $R^2=0.9999992$. The power of 0.53446 is close to, but statistically very distinct from, the exact value of ${\displaystyle \frac{1}{2}}$ that strict proportionality between $\tilde{n}$ and q would require. This is consistent with the broken Coulomb SO(4) symmetry that was discussed earlier.

We can also relate these coefficients to other known quantities. The asymptotic value of $\mathcal{J}=qj\left(0\right)=2\pi \tilde{n}$; taking this form $j\left(0\right)-j\left(\beta \right)=j_c{\beta }^{j_p}$ as correct (hypothetically); and recognizing we must have $j\left({\beta }_{\max }\right)=0$, we conclude that $j\left(\beta \right)=j\left(0\right)(1-{(\beta /{\beta }_{\max })}^{j_p})$, where ${\beta }_{\max }$ is a precisely defined numerical characteristic of the form of the potential, as in $Y\left(z\right)$ analysis; then the only undetermined parameter is the power $j_p$, which we have seen is close to but slightly greater than ${\displaystyle \frac{1}{2}}$. If $j_p$ were exactly ${\displaystyle \frac{1}{2}}$, we would have ${\displaystyle \frac{d\tilde{n}}{dq}}=j_0/\left(2\pi \right)$, a constant, i.e. strict linearity. The actual exponent $\approx 0.5345$ gives better agreement with the PM data, and introduces slight nonlinearity between q and n.

The same procedure for Hulthén and ECSC potentials (primary well) gives similar results as shown in Table 5 below, which is consistent with the contour plots shown in Figure 16, and the bivariate linear least squares fits to them. This, and the perturbation theory results, support the conclusion the approximate linear relation between $\sqrt{\mathcal{D}}$, n, and l is a shared intrinsic property of this extensive class of potentials which are expressible in terms of a linear coupling constant, here $\mathcal{D}$.

With $\tilde{n}=qj\left(\beta \right)/\left(2\pi \right)+\nu$, the dependence on q and $\lambda$ or $\mathcal{D}$ or l can be determined directly from $\beta =l(l+1)/\left(2q^2\right))$ and the chain rule. Evaluating these derivatives $dn/dl$ and $dn/dq$ etc for $q=20$, $l=10$ gives values within about $1\%$ of those in Tables 2–4.

Table 5. Fit coefficients of numerically evaluated phase space integrals $\mathcal{J}\left(\beta \right)/\sqrt{\mathcal{D}}=j\left(\beta \right)\approx j_0-j_c{\beta }^{j_p}$ for the three potentials. Selecting l and $\mathcal{D}$ specifies $\beta$, which allows the phase space integral to be approximately evaluated for $0\le \beta \le {\beta }_{\max }$, i.e. $j_0\ge j\ge 0$. The slight deviations of $j_p$ values from ${\displaystyle \frac{1}{2}}$ shows the relation between $\tilde{n}$ and $q=\sqrt{\mathcal{D}}$ is close to, but not exactly, linear.

Table 5. Fit coefficients of numerically evaluated phase space integrals $\mathcal{J}\left(\beta \right)/\sqrt{\mathcal{D}}=j\left(\beta \right)\approx j_0-j_c{\beta }^{j_p}$ for the three potentials. Selecting l and $\mathcal{D}$ specifies $\beta$, which allows the phase space integral to be approximately evaluated for $0\le \beta \le {\beta }_{\max }$, i.e. $j_0\ge j\ge 0$. The slight deviations of $j_p$ values from ${\displaystyle \frac{1}{2}}$ shows the relation between $\tilde{n}$ and $q=\sqrt{\mathcal{D}}$ is close to, but not exactly, linear.

At this point we consider the nested mysteries contained within our second mystery to be solved. We believe the extended trilinear relation between $q=\sqrt{\mathcal{D}}$, n, and l (including a constant offset), which defines a plane in q, n, l space, and the approximation $\mathcal{J}\left(\beta \right)\approx qj\left(0\right)[1-{(\beta /{\beta }_{\max })}^{j_p}]$ to be potentially useful.

4. Discussion

In our previous paper, discovery of the unexpected sawtooth structure in the ${\mathcal{D}}_c/n^2$ vs n plots called out for an explanation. The unusual nature of the ECSC potential – an infinite periodic series of positive energy barriers and negative energy wells – evoked all kinds of subtle and exotic hypotheses (shape resonances, subtle tunneling splittings inducing binding transition, extended states modeled in terms of transmission lines...). However, further analysis, backed up by numerical experiments using the Phase Method (PM) on specifically-modified ECSC potentials, quickly revealed a very different and relatively simple (“tick-tock”) mechanism – as $\mathcal{D}$ increases, states become bound in the primary well (ticks), but a new channel opens up when sufficient phase space area has accumulated to bind a new state (tock) in the secondary well. This causes an irregularity in the expected smooth accumulation of states (ticks) vs $\mathcal{D}$, particularly when plotted vs the total state count (tick+tock). Crucially, the common $\sqrt{D}$ factor for the two phase space areas divides out exactly, leaving a dimensionless number $1+{\mathcal{J}}_0/{\mathcal{J}}_1=39.8348...$, where ${\mathcal{J}}_0$ and ${\mathcal{J}}_1$ are the phase space areas of the primary and second wells for $\mathcal{E}=0$. The additive 1 is related to the total state count being the sum of the two contributions. The analysis is more complex and nuanced for $l>0$, but it is still analytically tractable, in terms of $\beta =l(l+1)/\left(2\mathcal{D}\right)$, using methods applied above for the Yukawa potential.

We further extended this analysis to derive a simple approximate equation to predict the approximate locations of all such “tick-tock” glitches. The first such “tockito” for the tertiary well was predicted and subsequently confirmed using a PM calculation up to $\mathcal{D}={10}^6$; the second todckito is predicted to occur at $\approx 5\ast {10}^6$ au. Similar behavior is expected for other multi-well potentials of this type, such as generalized ECSC potentials.

As it is non-standard, the NC method required explanation. It is a simple but nontrivial departure from conventional semiclassical/WKB analysis. It fundamentally is based on the 1932 paper of Dunham, which has been highlighted by Bender [16,17] and reprised in many contemporary quantum mechanics books, but these treatments almost universally truncate to low order Dunham’s exact infinite-order gradient expansion in powers of ℏ, essentially assuming the higher order terms sum to zero or are negligible. In contrast we assume they sum to a finite (usually nonzero) value, the magnitude of which we have shown is related to the amount of tunneling into the potential barrier. This is a simple but useful reinterpretation of a small portion of the vast body of semiclassical theory that has been developed by eminent physicists and mathematicians. Fundamentally, however, we rely on the 1932 paper of Dunham, and the 1977 paper of Bender et al.

The NC approach frees us to tackle states of low quantum number n and/or sharply varying potentials. To illustrate its utility we give examples of its use for homogeneous (e.g. power law) potentials, with a formula for eigenvalues that is exact for both harmonic oscillator and infinite square well, and that smoothly handles potentials of intermediate hardness. The PM and the NC methods proved to be a potent combination. Our previous paper describes the PM; the current paper describes the NC method. Both are explained with many applications in our 2024 book [9].

The NC method next was used to solve the second mystery, which is a peculiar historically noted [2] (but unexplained) linear relation between the total number of bound states $n^{\ast }$ and the screening length $\mathcal{D}$ for the Yukawa potential. Lam and Varshni [3] found a similar relation for the ECSC potential, but a little less well-resolved. In our current paper, from the PM-calculated results first shown here, we demonstrate the linear relation between $\mathcal{D}$ and $n^{\ast }$ also applies to the Hulthén potential. Solving this mystery required solution of several other nested mysteries, one of which involved surprising linear relation that found in our previous paper between $\sqrt{D}$, n, and l. This relation revealed, through a geometrical argument, that the total number of states below a threshold ${\mathcal{D}}^{1/2}$ grows approximately as the square of that threshold, and hence the number of states grows linearly with $\mathcal{D}$. It begged the question underlying the basis of the approximate linear relation between ${\mathcal{D}}^{1/2}$, n, and l.

This we tackled using NC analysis. First we used the scaling relation (applicable to this class of potentials) between screening length and coupling constant (the factor multiplying the attractive exponential term in the potential). By differentiating with respect to $q=\sqrt{\mathcal{D}}$ and $\lambda =\sqrt{l(l+1)}$ under the phase space integral the near-symmetry of the competing terms in the potential indicate there is a similar but opposing effect of changing q and $\lambda =\sqrt{l(l+1)}$. Evaluating the integral numerically gives near-linear relations between $q=\sqrt{D}$, n, and l, which solved the second mystery.

The close relationship between n and q for all three potentials is further studied in the last part of the paper, where $j-4\sqrt{\pi }\approx -12.09{\beta }^{0.5345}$. Strict linear dependence between n and $\sqrt{\mathcal{D}}$ would require a dependence $j\left(0\right)-j\left(\beta \right)\propto {\beta }^{1/2}$. All three potentials have a power of $\beta$ close to, but slightly greater than $1/2$. We predict other potentials in this class of potentials to have similarly close correlations between n, l, and $q={\mathcal{D}}^{1/2}$, and to be subject to analysis by these same methods. We expect they will show similar linearity properties.

In this paper we have connected disparate observations that all reduce to a shared common origin, which also indicates its applicability to a far broader group of short range linear-coupling constant potentials, i.e. those addressed by Klaus and Simon [10].

• Based on accurate extensive PM calculations:
• – linear relationship between total number of bound states $n^{\ast }$ vs screening length $\mathcal{D}$
• – is explained by near-constant density of states in $\sqrt{{\mathcal{D}}_c}$ vs n and l observed in plots
• – is explained by near-linear relation between $\sqrt{{\mathcal{D}}_c}$ vs n and l, observed in contour plots and fits
• – is explained via NC analysis quantitatively for linear coupling constant potentials
• – near-symmetry between $\mathcal{D}$ and $l(l+1)$ in scaled (z-space) phase space integral
• – near-linearity related to dependence of phase space integral on $\beta =l(l+1)/\left(2\mathcal{D}\right)$.
• – Single-well phase space integral approximately given by: $\mathcal{J}=qj\left(\beta \right)\approx q(j\left(0\right)-j_c{\beta }^{j_p})=qj\left(0\right)(1-{(\beta /{\beta }_{\max })}^{j_p})$; where $q\equiv \sqrt{\mathcal{D}}$ and $j_p={\displaystyle \frac{1}{2}}+\delta$; $\delta \ll 1$
• – the orientation of the gradient of $\tilde{n}$ in q, l space depends on ${\beta }_{\max }$ for that potential.

5. Conclusions

Both the accurate full-quantum PM calculations and the NC/phase space integral approach were essential for our solutions to both top-level mysteries addressed in this paper. Tables of ${\mu }_c=1/{\mathcal{D}}_c$ for the Hulthén, Yukawa, and ECSC potentials up to $\mathcal{D}={10}^5$ au and $l=0$–12 are deposited as supplementary data.

The primary initial goal of this paper was to understand the origin and nature of the sawtooth structure that was observed in the plots of ${\mathcal{D}}_c/n^2$ vs n plots for the ECSC potential, that was discovered in our previous paper [1]. Those high precision values were calculated using our “Phase Method” (PM), which was explained in that paper.

To explain the sawtooth structure paper, we have used those tables, additional numerical experiments using the PM, and a variant-semiclassical method we call "neoclassical" or NC method. This was used to explain the sawtooth structure in terms of a simple“tick-tock” mechanism, which can be expected to also be relevant for other multi-well potentials at large screening lengths. In the NC formulation the observed period-40 emerges naturally in terms of the purely numerical ratio of phase space areas: $1+J_2/J_1$.

The theoretical justification of the NC method is described, and how and why it differs from conventional semiclassical method. A key benefit of the NC approach is that it can be applied without caveats or reservations to low quantum numbers and hard potentials. In its current primitve form it is limited to calculating energy eigenvalues and critical binding parameters, but this is sufficient for our purposes. Example applications are given to homogeneous potentials. The tick-tock mechanism is illustrated using both NC and PM methods with a highly asymmetric quartic double well potential.

A second goal of this paper has been to understand the strangely linear dependence of the total number of states $n^{\ast }$ vs screening length $\mathcal{D}$ that was observed more than 50 years ago in both Yukawa and ECSC potentials, but unexplained. This is strange because the underlying dependence is $n\propto \sqrt{\mathcal{D}}$, so the rate at which states accumulation with screening length has a very specific character. Solving this mystery turned out to require the solution of other nested mysteries contained within it, all ultimately being reducible to the well-known scaling symmetry between screening length and linear coupling constant to which this class of potentials conforms. Hulthén, Yukawa, and ECSC potentials were all investigated and characterized, but the methods and behavior is expected to apply much more broadly. We also confirmed this linear relationship for the Hulthén potential.

Many of the approaches used in this paper for Hulthén, Yukawa, and ECSC potentials are expected to be equally applicable to other potentials in its class. These specifically are the ones that can be written in terms of a linear coupling constant (here $\mathcal{D}$), as described by Klaus and Simon. Xu et al [5] cataloged nine of them in their appendix, specifically the three addressed here, plus venerable potentials such as Woods-Saxon, Manning-Rosen, Lennard-Jones, Pöschl-Teller, and modifications of them. Based on our NC/phase space analysis, we predict the near-linear relationship between $\sqrt{\mathcal{D}}$, n, and l, as shown in our contour plots and fits for Hulthén, Yukawa, and ECSC potentials, should also apply to these others (modulo features from tick-tock glitches in multiwell potentials as was found for ECSC).

We hope the tables, methods, insights, and predictions in this paper will prove to be of lasting value.