Do equidistant energy levels necessitate a harmonic potential?

Experimental results from literature show equidistant energy levels in thin Bi films on surfaces, suggesting a harmonic oscillator description. Yet this conclusion is by no means imperative, especially considering that any measurement only yields energy levels in a finite range and with a nonzero uncertainty. Within this study we review isospectral potentials from the literature and investigate the applicability of the harmonic oscillator hypothesis to recent measurements. First, we describe experimental results from literature by a harmonic oscillator model, obtaining a realistic size and depth of the resulting quantum well. Second, we use the shift-operator approach to calculate anharmonic non-polynomial potentials producing (partly) equidistant spectra. We discuss different potential types and interpret the possible modeling applications. Finally, by applying $n$th order perturbation theory we show that \textbf{exactly} equidistant eigenenergies cannot be achieved by polynomial potentials, except by the harmonic oscillator potential. In summary, we aim to give an overview over which conclusions may be drawn from the experimental determination of energy levels and which may not.


Introduction
Thin films on surfaces with thicknesses of some ten layers act as a quantum well state perpendicular to the film surface. The in-plane directions are bulk-like and can be separated within the 3D Schrödinger equation. The variation over the different layers can then be described by a 1D model, most commonly by a square quantum well, where the eigenenergies increase quadratically and scale with 1/d 2 , d being the well width. In many cases, however, a deviation from this standard case is observed. For epitaxially grown SnO 2 , InSb, and Pb films with low charge densities, measurements with exponents about 1.5 have been reported [1][2][3]. Furthermore, Kröger et al. [4] and Hirahara et al. [5] investigated Bi films [6][7][8][9][10] and demonstrated that these exhibit equidistant energy levels with a 1/d thickness dependence [4,5]. For this, Kröger et al. prepared epitaxially grown Bi(111) films on Si(111) samples. The thicknesses varied between 20 and 100 bilayers (i.e. 8 nm to 40 nm). The conductance G was measured as a function of temperature T . A general dependence was fitted to this data and an effective bandgap E g was extracted by describing the temperature dependence of the conductance by thermal excitation, G ∝ exp(E g /k B T ), where k B is the Boltzmann constant. Hirahara et al. prepared Bi(001) films grown on Si(111) samples. The thicknesses varied from 7 to 40 bilayers (i.e. 2.8 nm to 16 nm). The energy levels were measured with angle-resolved photoemission spectroscopy (ARPES) and were found to be equidistant. This motivates a description by the well-known harmonic oscillator potential in contrast to the usually found square well potential, which was also stated by Kröger et al. An explanation given by Hirahara et al. is the phase shift accumulation model, originating from the Bohr-Sommerfeld quantization condition, which yields the 1/d-dependence if the dispersion near the Fermi energy is linear [5]. Nevertheless, the question must be asked whether equidistant energy levels actually necessitate a harmonic potential or whether alternative potential shapes with equidistant levels exist, especially considering that only a limited energy range is examined and a finite measurement uncertainty is present. The general answer is "no" [11], which was demonstrated by previous analytic derivations resulting in different potentials by generalizing the creation/annihilation operator. For this purpose, the factorization method [12,13] or the more general shift-operator approach [14,15] was used.
First principle studies of any material are expected to result in a potential which oscillates with the atomic positions of the real structure which is put in. As the harmonic oscillator is a valid model, an expectation could be that the potential of thin Bi films would on average and at a scale larger than the atomic structure be quadratic with additional spatial variations at smaller scale. Within this manuscript we focus on the model point of view and calculate model potentials satisfying the previous expectation. We hope that this publication will serve to interpret experimental results more carefully and clarify that even when the energy levels are known, a wide range of potentials must still be considered. In the following, we assume equidistant energy levels and investigate what follows for the potential. In Sec. 2 the experimental results are modeled by a truncated harmonic oscillator, yielding good agreement of the energy levels with the analytical harmonic oscillator case up to a certain energy and providing a more physical description including the specific layer thickness. Sec. 3 treats the 1D Schrödinger equation as an inverse problem by calculating the potentials from the given equidistant energy spectrum. This is done using the shift-operator approach. We present some previous analytical results and some new numerical calculations and possible interpretations concerning real
structures. In Sec. 4 we show by perturbation theory that polynomial potentials apart from the quadratic harmonic oscillator potential cannot exhibit an infinite number of exactly equidistant energy levels.

Comparison with experiments
In this section, we concentrate on the interpretation of experimental results for Bismuth films from the literature. The experimental results of Kröger et al. [4] and Hirahara et al. [5], i.e. the energy level spacings ∆E as a function of the Bi film thickness d, are depicted in Fig. 1 (colored symbols). d BL is the thickness of a Bi bilayer, which form a stable unit in (111) direction. For the Kröger data, the measured bandgap E g is interpreted as energy spacing ∆E. The inverse thickness dependence can be clearly seen. The solid line shows a fit with a d −1 thickness dependence.
To describe the experimental data and to connect the eigenenergies to a physical situation, the model of a truncated quadratic potential is used, where d is the thickness of the material and the potential at the well edge x = ± d 2 equals the well depth V 0 . Exemplary results are shown in Fig. 2 for V 0 /m * = 8 eV/m e and d = 16 nm, where m * is the effective mass and m e is the electron mass. The energy eigenvalues E n and the energy level spacings ∆E of the truncated potential are in very good agreement with the non-truncated potential (i.e. the harmonic oscillator) for V < 0.9V 0 , yielding with equidistant energy spacing ∆E, which scales like 1/d.    A regression according to Eq. (2) for reproducing the thickness dependence of the experimental data in Fig. 1 yields ∆E = 3.34 eV · d BL /d and V 0 /m * = 2.93 eV/m e . All data are in good agreement with the regression. We conclude that the experimental data may be well fitted with a truncated harmonic oscillator potential which yields the desired equidistant energy states. Nevertheless, the question remains whether this is the only possible shape for the potential. The next section thus deals with finding alternative potential shapes with the same energy spectrum.

Anharmonic oscillators with equidistant energy levels
In this section we show how the creation/annihilation operators of the harmonic oscillator can be generalized with the aim of finding further potentials with equidistant spectra. An approach to model equidistant energy levels is using the shift operator to treat the Schrödinger equation as an inverse problem, explained in the following. We first give a general introduction and show previous analytical results by Dubov et al. [14,15]. Afterwards we present new numerical calculations and possible interpretations concerning real structures.

Method
The general idea is to define for an arbitrary potential the shift operator L, which shifts the eigenstates ψ in energy: Here, ξ = mω/hx and = E/hω are the dimensionless real space and energy. If the states are equidistant, L transforms state ψ n into state ψ n+1 , resulting in the operator equation It can be solved by applying P is the momentum operator. α k (ξ) are free parameters. The first-order case (K = 1) leads to the harmonic potential U 0 = 1 2 ξ 2 and L equals the creation operator. The derivations can be found in the Supplementary Material. The solution of the second-order shift operator gives the isotonic oscillator [16][17][18][19] with either ξ > 0 or ξ < 0, depicted in Fig. 3   eigenstates and the same eigenenergies as the half harmonic oscillator, but slightly shifted in energy. The eigenstates have a similar form to the harmonic potential: The eigenstates are also shown in Fig. 3. It can be concluded that with this potential it is possible to model asymmetric minima in contrast to the pure harmonic potential. The third-order shift operator results in the differential equation with U = W + 1 2 ξ 2 . The prime denotes the first derivative. The equation can be solved as an initial value problem or a boundary value problem with any given conditions, resulting in many different solutions for the potential. Possible solutions obtained by a Darboux transformation [14,20] are The index m numbers the solutions. They are drawn in Fig. 4. A further shift of the solutions in ξ and U is valid, but not a superposition due to the non-linearity of the constituting equations. U 2m are quadratic-like potentials with a dip at the minimum. U 0 is the harmonic potential. The corresponding energy levels are equidistant with ∆ = 1.
The ground state of potential U 2m is additionally lowered by 2m. U 2m+1 are singular potentials at ξ = 0. U 1 is one solution of the second-order shift operator. The corresponding energy levels are equidistant with ∆ = 2 (∆ = 1 can be achieved by rescaling ξ =ξ/2). The eigenstates of these anharmonic potentials are also similar to the ones of the harmonic potential. Especially the potentials U 1 and U 2n+1 allow to construct asymmetric 1D models which describe equidistant energy levels, providing a much more flexible tool than a restriction to the analytical harmonic oscillator offers.
Further examples and their derivations can be found in the Supplementary Material.

Numerical results
Eq. (9) or its first integrals can also be solved numerically to achieve further types of potentials. We calculated the potential using the explicit Runge-Kutta method of 4th order or higher as implemented in the software package Mathematica [21]. We quantify the numerical error of the solution using the local relative residuum, i.e. the residuum normalized to the maximum absolute additive contribution to the differential equation, evaluated as a function of ξ. We get values below 10 −4 , which is sufficiently small to consider the solution correct. The Schrödinger equation is solved by discretizing the differential operator and the potential on an equidistant grid (∆x < 0.01 nm) and solving the resulting matrix eigenvalue problem as implemented in Mathematica. For this, the calculated spatial dimension is much larger than presented in the following, i.e. large enough for the potential at the boundary to be higher than the highest considered energy level, to prevent numerical errors from limiting the spatial dimension. Three different potentials and corresponding solutions of the Schrödinger equation are depicted in Fig. 5, Fig. 6, and Fig. 7. They represent three different types of solutions, named type-1, type-2, and type-3 potentials in the following. Fig. 5 shows a potential with two singularities and a minimum in between. The singularities behave like 1/x 2 , similar to the U 2n+1 potential. The wave functions are also of similar shape as the ones of the former potentials. The energy levels are roughly equidistant in the depicted range with ∆E ≈ 1040 meV. The higher states increase quadratically in energy instead of linearly. This could be due to the finite discretization and the resulting bad energetic resolution or due to a badly converged potential near the singularities. However, the 1/x 2 dependence near the singularities acts like a square box, resulting in quadratically increasing energy levels.  The second type of potentials, see example in Fig. 6, also has a 1/x 2 singularity for positive values of x, similar to U 2n+1 . Yet for x < 0, the potential oscillates and the strength of these oscillations increases with rising absolute value of x. The states of this potential can be separated into two classes: (1) States localized in the oscillation minima (red states in Fig. 6). (2) Delocalized states oscillating between the smooth wall and the potential oscillation peaks of the corresponding energy range (blue states in Fig. 6). The depicted example exhibits nearly linearly increasing energy levels. Additional small irregular variations are due to the 2 classes of states. The insets in Fig. 6 show the energy level spacings of these distinct classes. The class-1 states (upper inset) show linearly increasing energy levels with ∆E ≈ 358 meV. The class-2 states (lower inset) vary between ∆E = 300 meV and ∆E = 450 meV. The superposition of these two classes of states with different energy level spacings leads to the variations in the combined picture and although the class-2 states are not truly equidistant, the deviations from the overall linearity of all states are small. The overall average energy level spacing is ∆E ≈ 180 meV. Fig. 7 depicts the third type of potentials with no singularities, but oscillations in both directions and also displaying two different classes of states. Class 1 are states mostly localized within the potential oscillations at the left or right (red and green in Fig. 7). There are also states similar to the depicted ones, but located at the other side. Class 2 are delocalized states oscillating between the left and right potential oscillation peaks of corresponding energy (blue in Fig. 7). Energy levels increase almost linearly; the small periodic variations are due to the two classes of states. The insets in Fig. 7 show the energy level spacings of these distinct classes. It can be clearly seen that for each class linearly increasing energy levels are present. Class 1 (upper inset in Fig. 7) has ∆E = 358 meV for states located at the left (green) as well as for states located at the right (red). The states are just shifted by a small energy due to the slightly asymmetry of the potential. Class 2 (lower inset) has ∆E ≈ 267 meV. The class-2 energy level spacings vary slightly (about 10%). But as this is a differential value, deviations from linearity are sufficiently small and in good agreement with experimental results. The superposition of these two classes of states with different energy level spacings leads to the irregular variations in the combined picture. In total, the energy levels show a linear dependence with an overall average energy level spacing of ∆E ≈ 108 meV. The additional irregular variation on the large scale are small. With regard to experimental results with some error bars these irregularities may not be visible.
The presented 1D potentials allow more interpretations concerning the comparison with real structures than a pure 1D harmonic oscillator model. They can be related to different possible cases occurring in reality. For example, a thin, free-standing slab, surrounded by vacuum, has states which are trapped within the slab in the direction perpendicular to the slab and exponentially decaying into vacuum. This could be modeled by the type-1 potential of the numerical calculations, which diverges at the slab-vacuum interface. A thin layer on a substrate could be modeled by the type-2 potential. This diverges at the layer-vacuum interface and oscillates into the bulk-substrate, representing the periodic structure. A thin sandwich-layer between two other substrates could be modeled by the type-3 potential, which oscillates into both directions to describe the penetration into the periodic bulk material. Another interpretation could be to use the type-2 or type-3 potential to describe the oscillating atomic structure within the layer itself. The experimental results from literature [4][5][6] have the structure silicon-bismuth-vacuum and fall into these categories. The ∆E = 358 meV achieved in Fig. 6 is in the range of possible experimental ∆E of Fig. 1. Using the result of the fit in Eq. (2), this corresponds to a thickness of 9.3 bismuth bilayers, i.e. 3.7 nm. Within the potential in Fig. 6 this approximately matches the width of the non-oscillating part at the right near the singularity. The oscillating part represents the transition into the bulk-like part with wave lengths between 0.2 nm (for −10 nm < x < −8 nm) and  In this section spectra have been studied which are equidistant in a limited energy range or contain a subset of equidistant states. In the next section we show what can be learned from perturbation theory. The difference to the former shift-operator approach is that, in the following, searching for equidistant eigenstates means an infinite number of exactly equidistant states. The formerly shown results are less strict and do not belong to this for the following reasons: Potentials U 2m+1 and the numerical type-1 potential have non-polynomial shapes and quadratic non-integrable singularities. Potentials U 2m consider parts of the spectrum being non-equidistant (i.e. the ground state). The numerical type-2 and type-3 potentials lead to states which split into two classes with separate energy spacing each.

Perturbed 1D harmonic oscillator
In this section, we use perturbation theory to show that an exactly equidistant energy spectrum over the whole energy range cannot be achieved by a polynomial potential other than the harmonic (quadratic) potential. To this end, the 1D harmonic oscillator is perturbed with an arbitrary polynomial perturbation potential u(ξ) = N j=0 u j ξ j , − ∂ 2 ∂ξ 2 + ξ 2 + λu(ξ) ψ n (ξ) = (0) n ψ n (ξ).
As the equations are linear, we consider each monomial ξ j individually. The perturbation potential u(ξ) = ξ j in the basis of the harmonic oscillator eigenstates gives the perturbation coefficients where H n (ξ) are the Hermite polynomials. The diagonal elements u (j) kk equal the first-order energy correction (1) k . They are equal to zero for odd values of j. For even numbers j, u (j) kk is a polynomial in k of order j/2 as shown in Fig. 8 for different monomial perturbations u(ξ) = ξ j . For the smallest values of j one gets Consequently, if the harmonic oscillator is perturbed with a quadratic potential, (1) k is linear in k, reproducing the equidistant energy levels of the harmonic oscillator. If the harmonic oscillator is perturbed with a non-quadratic polynomial potential, (1) k is a superposition of different powers of k, yielding a different spectrum. Thus, only the quadratic polynomial gives equidistant first-order energy corrections.
The non-diagonal elements u   The u (j) k,k+l are shown in Fig. 9 for different l and for different polynomials u(ξ) = ξ j . The non-zero elements for a specific l can be written as k,k−l ) 2 is a polynomial in k of order j for all l. The previous diagonal case l = 0 is included. In particular, for the smallest values of j one gets p (2) k,2 = 1 , p (4) k,2 = −1 + 2k , p With this, the second order energy correction is a polynomial in k of order j − 1, as the highest order in k cancels out with the difference (u k,k−l ) 2 . Consequently, the second order energy corrections scale linearly with k only for the quadratic potential (j = 2). Treating the higher-order perturbations similarly, in the k 3j/2 -scaling two hierarchic differences lead to a cancellation of the two highest orders of k. With this, (3) k is a polynomial in k of order 3j/2 − 2 and also linear scaling only for the quadratic potential.
Another argumentation can be done using the pth order perturbation corrections, written as a recursion: The general argumentation is as follows: For the harmonic oscillator, c k scales like k x+j/2 +k x+y . Consequently, a linear dependence of k on k can only be achieved if y = 1, thus x = 0, giving j = 2, which is the harmonic oscillator. In total, for a general perturbation potential u(ξ) = N j=0 u j ξ j , the perturbed energy levels are a sum of linear functions of k for the quadratic case and a sum of functions of k with different exponents for the non-quadratic cases, which does not sum up to a linear dependence.
Though it may seem so, Sec. 3 and Sec. 4 do not contradict each other. The previous argumentation holds only for a polynomial potential which produces an equidistant spectrum for all eigenstates and with no tolerance in the level spacing. It does not apply to power series and spectra which are only equidistant in a specific energy range or approximately equidistant spectra, as required to model experimental spectra.
The derivations of this section and some more results can be found in the Supplementary Material.

Summary and conclusions
In this publication, we discuss the implications that may be drawn from the experimental observation of an equidistant energy level spacing. Though a harmonic potential comes to mind immediately, this shape of potential is sufficient, but not necessary for the equidistant spacing, especially considering that experimental results cover a finite energy range with a finite measurement uncertainty. The question is thus more involved than is visible at first sight. In this publication, we discuss different aspects of the problem.
In the first part we compare a 1D harmonic oscillator model with experimental results, looking at the analytical harmonic potential as well as on a truncated version. The energy levels measured in thin Bi films of different thicknesses and the resulting thickness-dependence fit the theoretical model with comparable energy values and dimensions. Thus, the 1D harmonic oscillator is a valid model to describe the quantum well states in thin films as shown for Bi films.
However, this still leaves the question whether it is the only possible potential to explain the measurements. Therefore, in the second part of the paper, we discuss anharmonic potentials of non-polynomial shape and apply the shift-operator approach to explicitly calculate potentials with an equidistant spectrum. We present a number of analytic and numeric examples representing different types of potentials including asymmetric minima, singularities, and oscillations. These potentials seem well-suited 1D models to describe the potential of thin slabs, thin films on substrates or thin sandwich films between two substrates. They also fit the expectation that first principle calculations of real materials should result in oscillating potentials which are in average quadratic to fit the harmonic oscillator model.
In the last part, we apply nth-order perturbation theory with a polynomial perturbation potential to the quantummechanical harmonic oscillator problem. We calculate and discuss the perturbation integrals u k,k−l and the resulting nth order energy corrections as a function of k for single monomials ξ j , showing that u 2 k,k−l is a polynomial in k of order j and arguing that every non-quadratic perturbation monomial leads to a non-quadratic k-dependence of the nth-order energy corrections, leaving only the harmonic term to produce an equidistant spectrum.
Thus we conclude that all possible anharmonic potentials producing exactly equidistant spectra must have a nonpolynomial form. Anharmonic non-polynomial potentials exist, which could describe real structures and produce spectra which mirror the experimental results.

Do equidistant energy levels necessitate a harmonic potential?
Fabian Teichert, Eduard Kuhn, Angela Thränhardt Institute of Physics, Technische Universität Chemnitz, 09107 Chemnitz, Germany E-mail address: fabian.teichert@physik.tu-chemnitz.de Abstract: In this supplementary material, the shift-operator approach [Chaos 4 (1994), 47] is derived and a few resulting potentials and states are depicted. Furthermore, The solution of the harmonic oscillator and its perturbation with different monomials is derived. A few examplary results are shown.
With this, a resulting operator equation can be set up using commutators: K is the order of the shift operator. α k (ξ) are free parameters. Inserting Eq. (7) into Eq. (6) and comparing the coefficients yields K + 2 differential equations for K + 1 unknown functions α k (ξ) and the additional unknown function U (ξ).

First-order shift operator
Inserting the first-order shift operator into Eq. (6) yields The solution for the α k yields the condition Integrating this yields the harmonic oscillator potential The corresponding shift operator is the creation resp. the annihilation operator. The corresponding solution of the Schrödinger equation is where H n (ξ) are the Hermite polynomials. The potential and the states of the first-order shift operator are depicted in Fig. 1.

Second-order shift operator
Inserting the second-order shift operator into Eq. (6) yields The solution for the α k yields the condition Multiplying with ξ and summarizing using the product rule yields Integrating and dividing by ξ 2 yields The shift operator is The solution of the Schrödinger equation can be achieved with the asymptotic behaviour and the consequential approach It follows The potential and the states of the second-order shift operator are depicted in Fig. 2.

Third-order shift operator
Inserting the thirt-order shift operator into Eq. (6) yields The solution for the α k yields the shift operator and the condition with U = W + 1 2 ξ 2 . Multiplying with ξ and summarizing using the product rule yields the first integral Multiplying with W , using the identities of Eq. (36) and summarizing using the chain rule yields another first integral Using a Darboux transformation, a possible solution can be written in the form with − 1 2 χ µ + 1 2 ξ 2 χ µ = µχ µ .
With µ = − m + 1 2 two possible classes of solutions are U 2n are quadratic potentials with a dip at the minimum. U 0 is the harmonic potential. The corresponding energy levels are equidistant with ∆ = 1. The ground state of potential U 2n is additionally lowered by 2n. U 2n+1 are singular potentials at ξ = 0. U 1 is the one of the second-order shift operator. The corresponding energy levels are equidistant with ∆ = 2 (∆ = 1 can be achieved by rescaling ξ = ξ /2). The potentials of the third-order shift operator are depicted in Fig. 3. The potentials U 2 and U 4 and corresponding lowest states are shown in Fig. 4 and Fig. 5.