Exact and approximate energy sums in potential wells

Sums of the N lowest energy levels for quantum particles bound by potentials are calculated, emphasising the semiclassical regime N  ≫  1. Euler-Maclaurin summation, together with a regularisation, gives a formula for these energy sums, involving only the levels N  +  1, N  +  2…. For the harmonic oscillator and the particle in a box, the formula is exact. For wells where the levels are known approximately (e.g. as a WKB series), with the higher levels being more accurate, the formula improves accuracy by avoiding the lower levels. For a linear potential, the formula gives the first Airy zero with an error of order 10−7. For the Pöschl–Teller potential, regularisation is not immediately applicable but the energy sum can be calculated exactly; its semiclassical approximation depends on how N and the well depth are linked. In more dimensions, the Euler–Maclaurin technique is applied to give an analytical formula for the energy sum for a free particle on a torus, using levels determined by the smoothed spectral staircase plus some oscillatory corrections from short periodic orbits.


Introduction
This is a study of sums of the first N energies E 1 , E 2 ... of quantum systems with discrete spectra, that is concentrating on the asymptotics for large N. This work was motivated by applications in density-functional theory [1,2]. A sufficiently accurate and robust approximation to the sum of the lowest N levels of a potential could avoid the need to solve the Kohn-Sham equations, thereby increasing the number of electrons that can be treated. An independent interest is the combination of techniques we employ. The sum S(N) is a function of the discrete index N. We will need to represent it for continuous N. One way is simply to interpolate linearly between the integers, i.e.
This function is continuous but not smooth: its slope is discontinuous at integers N. Our aim here is to find formulas for S(N) that are smooth, while still coinciding with, or being close to, the exact or approximate S(N) at integers. Figure 1 displays the difference between the integer, linearly interpolated, and smooth functions S(N). Our technique, which works when all states are bound, is explained in section 2. A regularisation procedure enables the sum of the first N energies to be replaced by the sum of levels N + 1, N + 2.... In the usual cases, where the levels are known only approximately, with the levels below N being less accurate than those above, this has the advantage of avoiding contamination of the sum by the lower energies.
Sections 3-5 describe preliminary examples: simple cases where the levels are known and S(N) can be evaluated analytically. We include the harmonic oscillator (section 3), and the particle in a 1D box (section 4), simply to illustrate the general technique. For the Pöschl-Teller potential (section 5), whose depth is finite and which binds a finite number of levels, the technique is not obviously applicable; the reason for including this additional exactly solvable case is because the potential contains a parameter (its depth), whose asymptotics interact interestingly with N.
The heart of the paper is section 6, where S(N) is calculated for the odd energies of the linear potential |x|; these are the zeros of the Airy function [3,4]. This is an explicit example where the approximate evaluation of S(N) yields extraordinary accuracy when the individual levels are given, by WKB theory, in the form of a divergent series.
In more than one dimension, the energies fluctuate pseudo-randomly about those of a smoothed spectrum. This is the case whether the classical motion is integrable (level fluctuations Poissonian [5]) or chaotic (levels distributed according to random-matrix theory [6][7][8]). Then S(N) can still be evaluated using the method of section 2, with the energies being those of the fully or partially smoothed spectrum. In section 7 we illustrate this for particles on a 2-torus.
In a complementary study [9], the statistics of the energy sum (1.1) has been calculated using random-matrix theory; we thank a referee for bringing this paper to our attention.
Although some of our approximations are semiclassical, we will save writing by setting Planck's constant = 1. In the examples we treat, can be eliminated by scaling; when this is not immediately obvious, it will be explained.

General method
The lower levels in S(N) can be eliminated by writing the sum in (1.1) formally as the difference of two infinite divergent sums: (2.1) Both sums must be regularised by analytic continuation. For S (∞), which contributes a constant to S(N), we accomplish this via the partition function [10] (also called the propagator for the heat equation, i.e. heat kernel); for a Hamiltonian H, this is Usually this diverges as t → 0, but its Laurent expansion may include a term proportional to t; hence the regularisation [10], evaluated at s = −1.) For the infinite sum over N + 1 ⩽ n < ∞, we use the Euler-Maclaurin formula [3]. In this, the discrete energies E n are represented by any interpolation for continuous n, and the discrete sum is expressed as an integral over continuous n, corrected by derivatives at the summation limit. Thus our theoretical expression for the energy sum is B 2k are the Bernoulli numbers [3], expressed in terms of the Riemann zeta function ζ(s) for even integers: The formally divergent integral over n must also be interpreted by the same regularisation. This is accomplished by the continuous-variable counterpart of the partition function procedure. In all cases we encounter, E n involves monomials, whose interpretation, intuitively clear but explained in detail in appendix A, is In a few cases, described in the next two sections, the series involving k in (2.4) terminates, and Euler-Maclaurin is exact. More generally, as in the case described in section 6, the series over k is semiconvergent: it is an asymptotic series, whose terms get smaller and then increase [11]. An alternative representation of S(N) follows from transforming the sum (1.1) using the Poisson summation formula [12]. This yields a series of integrals over n that are then approximated by expanding about their endpoints n = N. But this is simply a complicated way of reproducing the result of the Euler-Maclaurin representation (2.4). Poisson summation will however play a role in section 7, not in calculating the energy sum but in calculating the energies E n involved in the sum, which will include some oscillatory contributions to the level counting function [13,14] (spectral staircase). A different alternative would be to regularise not by using S (∞) but demanding S(0) = 0; this makes no difference for the examples in sections 3-5 and negligible difference in sections 6 and 7.

Harmonic oscillator
The energies, and the sum (1.1), are (3.1) To interpret this in terms of the Euler-Maclaurin representation (2.4), we need the regularised sum over all n, in terms of the partition function (2.2) and (2.3). For this case, K(t), and its Laurent expansion, can be evaluated explicitly: In the theoretical expression (2.4), we need (2.6) to interpret the infinite integral, leading to contributions that add to reproduce the exact energy sum: An interesting minor variant is the spectrum without the zero-point energy (or, equivalently, the positive eigenvalues of an angular-momentum component), where S(N) is the sum over the positive integers: The constant is S (∞) = 1 + 2 + 3 + · · ·, discussed in recreational mathematics [15]. The regularisation (2.2) gives one interpretation of the old result in terms of the Riemann zeta function ζ(s), i.e. S (∞) = ζ (−1) = − 1 12 :

1D box
For a free particle in the space 0 ⩽ x ⩽ π between hard walls, the energies (with = 2 × mass = 1), and the energy sum, are For this case, the partition function is a theta function θ 3 [3], whose Poisson transformation shows that its expansion contains no term proportional to t, so the sum over all energies is zero: Again, Euler-Maclaurin representation reproduces the exact S(N):

Pöschl-Teller potential
This exactly solvable potential, with finite depth D supporting discrete energy levels between The energies (with = mass = 1), and the energy sum, are [16] For this potential, the number of bound states is finite, so the regularisation technique of section 2 is not immediately applicable-and regularisation would bring no advantage, because the system is exactly solvable. The Euler-Maclaurin formula can be applied in its conventional form [3], involving the lower limit n = 1 of (1.1), and is easily confirmed to reproduce the exact S(N). Temporarily reinstating , it is clear from the Schrödinger equation that it can be eliminated by the scaling D → D/ (and corresponding scalings for E n and S(N)). Thus the semiclassical regime ( small) is D ≫ 1. We are interested only in the bound states, so the total number n max (D) of energies in the sum (1.1), and the corresponding maximum value S max (D) of the energy sum, and their leading-order large D approximations, are 3) For this exactly solvable potential, large D approximations do not lead to divergent series, as for the WKB approximations for general potentials [11] (an example is explored in the following section). Instead, the series is convergent, and associated with the expansion ä + · · · There are many ways to implement the large D (semiclassical) approximation of S(N;D), depending on whether and how N is incorporated by linking it to D. With no linkage, i.e. approximating S(N;D) with N fixed, all energies below N lie in the range 0 < E < O(√D), that is, near the bottom of the well. More natural is to scale N with the total number of bound states, and S with its maximum value. Again this can be done in several ways, depending on whether the scalings are done with n max and S max or n max0 and S max0 in (5.3). The scalings n max0 and S max0 lead to simpler and faster-converging large D approx imations. Therefore we represent the level indices and the Pöschl-Teller energy sum in the form With these scalings, the ranges of the level indices and the sum are 0 < ν < 1 and 0 < S sc < 1, and the semiclassical regime corresponds to D large, ν fixed. Explicitly, This large D expansion is extraordinarily accurate, even for D = 1, as figure 2 illustrates. (The alternative scaling, using n max and S max , introduces additional terms, of order D 1/2 , D −1/2 , D −3/2 ...). The case N = 1 is interesting. Since S(1) = E 1 , the energy sum is an alternative way of approximating the ground-state energy, to be compared with the direct large D expansion The first three approximations differ by half-integer powers of D, in contrast to the energy sum (5.5), whose approximations differ by integer powers of D. Therefore, as table 1 illustrates, the Table 1. Ground state energies E 1 for Pöschl-Teller well depths D = 1, 2, 3, calculated directly from the expansion (5.6), and from the scaled energy sum (5.5).
early terms of the energy sum provide a more efficient way of approximating E 1 : S sc1 and S sc2 give the same approximation as E 12 and E 13 .

Hard-wall linear potential
General 1D potential wells are not exactly solvable, and semiclassical approximations, in terms of series in powers of , are given by the WKB approximation [17,18]. These explicit series take the form of a 'quantum action' I (E; ), where I(E;0) is the classical action. The energies E n that must be summed to obtain S(N) are determined implicitly by a quantum condition I (E n ; ) = n − µ, where µ (a Keller-Maslov index [19]) depends on boundary conditions. Getting the energies E n requires reversion of the semiclassical series. We illustrate the essential features of approximations to the energy sum with the levels of a linear potential, about which much is known [4]. The potential is V(x) = |x|, whose wavefunction ψ(x) is the Airy function [3], and the odd and even energy levels (with = 2 × mass = 1) are given by the zeros a n and a n of the Airy function and its derivative. We consider just the odd levels, corresponding to a linear potential on the half-line x > 0 and a hard wall at x = 0: The similar formalism for the even levels, and the totality of even and odd levels, is described at the end of appendix B. By elementary scaling, the semiclassical approximation corresponds to n ≫ 1.
For the nth zero, the approximation of order M is given [3,20] by In the literature [20], the coefficients T 0 to T 9 are listed. We needed many more; appendix B describes the way we calculated them, and lists T 0 to T 20 . The asymptotics of the coefficients is given in (B.8), indicating that the terms in the series (6.2) have the generic 'factorial/power' form [11]. Therefore the smallest term is near ). This will extend the cases considered in sections 3 and 4, where Euler-Maclaurin summation is exact because the sums over m and k terminate, to the present situation, in which the energies are known only approximately. We require the regularised sum S (∞), whose value is known [4,21] to be zero. We apply (2.4) to the sum of contributions m in (6.2), using (2.6) for the regularised integrals over n, and explicitly evaluate the elementary derivatives with respect to n, and collect terms of the same order in 1/N. Thus The first few terms of the series are   Table 2 shows the optimal errors for different N, calculated in this way. Even for relatively small values of N the accuracy is extraordinary, and the fractional error agrees well with the elementary asymptotic estimate exp(-2π(N + 3/4))/(2π(N + 3/4)) 3/2 obtained from the leading terms in (6.4), the estimate (B.8) for T m , and Stirling's formula.
Note in particular the optimal error for S 1 = E 1 in table 2. Its magnitude, of order 10 −7 , should be compared with the errors in the ground-state energy calculated directly from the series (6.2) for the Airy zeros, also terminated by adding half the smallest term. Table 3 shows the comparison, indicating that the large N series for S(N), when applied to N = 1, outperforms the direct WKB eigenvalue asymptotics (6.2). The reason is that the energy sum for N = 1, when evaluated by the method of section 2, involves the energies E 2 , E 3 , ... , whose WKB approximations are more accurate than that for E 1 .

Towards more dimensons: quantum particle on a torus
The Euler-Maclaurin representation (2.4) for the energy sum requires the energies E n to be represented analytically in terms of a monotonic function of continuous n. Achieving this in d > 1 dimensions is problematic, because the successive energies fluctuate pseudo-randomly. For separable d-dimensional Hamiltonians, the energies can be represented (exactly or approximately) in terms of a set of d quantum numbers, but the successive energies usually involve very different numbers in the set, and the statistics of the fluctuations is Poissonian [5]. For general nonseparable Hamiltonians (e.g. when the classical trajectories are chaotic), there are no quantum numbers, no explicit expressions are available for the individual energies E n in increasing order, and the statistics are those of random-matrix theory [6,7]. However, it is always possible to represent the spectrum in terms of a level counting function (spectral staircase) N (E), that is the sum of two contributions [14]. These are determined uniquely within semiclassical approximation schemes: a series of terms that are smooth functions of energy (monomials in ), arising from gross features of the classical phase space; and a series of terms that are oscillatory in , arising from classical periodic orbits [13,22] and representing fluctuations about the smoothed staircase. Thus, with Θ denoting the unit step, the staircase is Key to using the representation (7.1) to calculate the energy sum is the observation that the two contributions depend differently on : in which [14] Ω(E) is the phase space volume enclosed by the energy surface E, S j (E) is the action of the periodic orbit labelled j , the amplitude A j (E) depends on the stability of the orbit, and the exponent µ is a rational fraction between 0 (for chaotic dynamics) and (d − 1)/2 (for completely integrable dynamics). The quantum condition (7.2) must be a monotonic function of E, so the derivative ∂N /∂E (i.e. the level density) must be positive. The contribution N sm (E) is always positive and O −d , and this must dominate the contributions to N osc (E), which are O T j (E) −(µ+1) , in which T j is the period of the orbit j . Since µ + 1 < d, the desired monotonicity will hold for orbits whose period is not too long, that is, whose periods are smaller than O d−µ−1 . Therefore in calculating the energy sum using (2.4) we can use approximate energy levels incorporating these short orbits in the semiclassical representation Table 3. Errors in the approximations (6.2) and (6.6) for the lowest odd energy (smallest Airy zero) for the potential |x|; smallest errors bold. To illustrate this, we employ the simplest nontrivial example: a quantum particle moving freely on a 2-torus, that is a rectangle of side lengths 2π and 2πα with opposite sides identified, or, equivalently, the periodised plane with the rectangle as unit cell. The exact energies are These are easy to compute and arrange in order, to give the sequence E n . This case is particularly challenging for semiclassical approximations because there are many degeneracies (e.g. ±l, ±m, exchanging l and m), so the spectral fluctuations are large. Its simplicity has the advantage of avoiding complications, irrelevant for our purpose here, from perimeters, curvatures, and corners [23], that occur for 'quantum billiards' with boundaries.
For this system, the spectral staircase function, in the form (7.1) can be calculated exactly using the Poisson summation formula: in which N sm (E) = παE comes from the term j = k = 0 in the sum, and in N osc (E), involving the Bessel functions, is the length of the periodic orbit labelled j ,k. (These orbits form continuous families, each consisting of paths in the periodised plane from points x, y in the rectangle to points x + 2πj, y + 2παk.) Figure 4 shows the emergence of the steps as increasing numbers of periodic orbits are included, with the concomitant loss of monotonicity, arising from the fact that although the strengths of the contributions from the longer orbits decrease, they also oscillate faster.
To determine the ordered energies E n , it is necessary to invert the quantization condition (7.2), incorporating a number of periodic orbits small enough to ensure that the staircase is monotonic. We accomplish this approximately, using the fact that the oscillatory contributions are smaller than the smooth contributions. We start with the zero-order energies E 0,n , i.e.   and include the oscillatory contributions, with a finite number of periodic orbits, by solving (7.2) iteratively. The first iteration suffices, yielding the approximate energies These are the energies to be included in the Euler-Maclaurin formula (2.4). We also require S (∞); this is zero, because the partition function is the product of two factors of the 1D form (4.2). We have confirmed numerically that at the present level of approximation, it is not necessary to include the derivative terms in (2.4). Thus the energy sum, including periodic orbits with |j |,|k| ⩽ K, is (7.9) Figure 5 shows the exact and approximate energy sums for N ⩽ 15. Even for these small values of N, the curves are hard to distinguish, even for K = 0, i.e. when no periodic orbits are included. Figure 6 shows the errors over a much larger range; as described below, these increase with N, because of the fluctuations, but are much smaller when periodic orbits are included. The errors are larger than the mean spacing between adjacent energies, which is constant (~1/πα = 0.379), but are always very small compared with the value of S(N) ~ N 2 /2πα: e.g. S(5000)~4.73 × 10 6 . Figure 7 shows magnified regions of figures 6, illustrating how the scale of the fluctuations gets smaller as more orbits are included. These calculations are for side ratio α 2 = 1/√ 2. For different α, other calculations (not shown) give similar results.
We can estimate the increase of the fluctuations with N by calculating the r.m.s. error σ K (N), where the variance σ K (N) 2 is calculated by averaging over the oscillations in S K (N). Thus, using asymptotics of the Bessel functions, using (7.7) for E 0,n , and noting that a similar calculation shows that the second sum in (7.9) grows more slowly with N, where the sum is over all j , k excluded by the square of side K in (7.9); this sum converges, and is O((K + 1) −3 ). Thus S K (N) grows as N 3/4 , consistent with the calculations illustrated in figures 6. From the argument of the Bessel functions, the oscillation frequency as N increases scales as (K + 1) 2 . As a check of (7.9), we also calculated the sum energy directly from the individual approximate levels given by (7.8). The results are indistinguishable in every detail from figures 6 and 7, probably because any improvement is masked by the fluctuations. The close similarity indicates the effectiveness of the analytical sum (7.9), generated from the Euler-Maclaurin procedure.
The numerical values in table 4 show the energy sum, and the errors for increasing periodic-orbit truncation K, for small values of N. All errors (except one) are small in comparison with the mean level spacing. For such such small N it is hard to see a consistent picture, though there is a tendency for the errors to decrease with increasing K.

Concluding remarks
We have shown that the technique of section 2, in which regularisation is combined with Euler-Maclaurin summation, is a widely applicable and effective way to evaluate the energy sum S(N) defined by (1.1). In some cases (sections 3 and 4) it is exact, and for 1D potentials where the energies are approximately known it can give the semiclassical energy sum with high accuracy, as illustrated (section 6) for the Airy zeros. In more dimensions, the method is effective for the energy sum of levels of the partially smoothed spectrum, when some periodic orbits are included, as illustrated (section 7) for the particle on a 2-torus. For finite-depth wells, unregularised Euler-Maclaurin works, at least for the exactly solvable Pöschl-Teller potential; for this case the natural semiclassical expansion parameter is a combination of N and the well depth.
Our preliminary analysis suggests several avenues for further research.
• Extending the linear-potential analysis of section 6 to general infinite potential wells, by applying the known [24] WKB series for the quantization condition analogous to (B.4), thus building on the leading-order correction already obtained [25], where some of the formulas were obtained by elementary means. This would incorporate the regularised sums S (∞) for general potentials [10,18,26,27]. • For finite-depth wells, going beyond the Pöschl-Teller example of section 5 to include WKB-generated energies for potentials that are not exactly solvable. In cases where the sum of all the bound states (the counterpart of S (∞)) is not known, this would require application of the Euler-Maclaurin formula directly to (1.1) in its familiar (i.e. not regularised) form, possibly incorporating the requirement S(0) = 0. • Where WKB is applied, further increasing the accuracy of the summations beyond that illustrated in section 5, by employing asymptotic techniques [11,[28][29][30][31] involving resummation of the divergent tail of the relevant series beyond their smallest terms. • For nonanalytic potentials, for example a smooth well truncated at finite depth, summing the energies would require modification of (6.4), incorporating techniques recently developed elsewhere [32]. • For higher-dimensional potential wells more general than the particle on a 2-torus, incorporating multidimensional semiclassical techniques [33] to get explicit formulas for the correction terms ('Weyl series') in the series for the smoothed spectrum N sm (E) and the periodic-orbit terms N osc (E).  By analogy with the partition function procedure (2.2) and (2.3), we define the integral in (2.6) as A double integration by parts eliminates the logarithm, enabling the contour integal to be evaluated by the residue at ξ = 0. Its value is ν/(ν + 1), completing the derivation of (2.6). This way of getting what seems an intuitively obvious result seems complicated; perhaps there is a simpler one. The reversion of series, and evaluation of the energy sum, are essentially the same as for the odd states. A formal combined quantum condition, including both cases, is z (E n ) = n − 1 2 1 2 π − 1 2 Im log Σ n,M Σ n,M + (−1) n log Σ n,M /Σ n,M , Σ n,M = Σ (z (E n (M))) , Σ n,M = Σ (z (E n (M))) , (B.12) in which n = 1,2,3..., with even n giving the odd levels and odd n giving the even levels.