The Kronig-Penney model extended to arbitrary potentials via numerical matrix mechanics
aa r X i v : . [ c ond - m a t . o t h e r] N ov The Kronig-Penney model extended to arbitrary potentials via numerical matrixmechanics
R. L. Pavelich ∗ and F. Marsiglio † Department of Physics, University of Alberta, Edmonton, AB, Canada T6G 2G7 (Dated: November 14, 2014)The Kronig-Penney model describes what happens to electron states when a confining potentialis repeated indefinitely. This model uses a square well potential; the energies and eigenstates canbe obtained analytically for a the single well, and then Bloch’s Theorem allows one to extendthese solutions to the periodically repeating square well potential. In this work we describe howto obtain simple numerical solutions for the eigenvalues and eigenstates for any confining potentialwithin a unit cell, and then extend this procedure, with virtually no extra effort, to the case ofperiodically repeating potentials. In this way one can study the band structure effects which arisefrom differently-shaped potentials. One of these effects is the electron-hole mass asymmetry. Morerealistic unit cell potentials generally give rise to higher electron-hole mass asymmetries.
I. INTRODUCTION
The Kronig-Penney model [1] remains the paradigm forthe demonstration of energy bands, separated by gaps,in periodic solids. This in turn is used to introduce thenotion of metals vs. (band) insulators, dependent on theposition of the chemical potential (and therefore depen-dent on the number of electrons) in the solid. The modelis appealing insofar as it only requires solutions to theSchr¨odinger equation for the simplest potential (a con-stant), and the equation for determining allowed energiescan be written analytically in terms of elementary func-tions. This exercise is also enlightening as it puts intopractice what might be at first viewed as a somewhat ab-stract theorem (Bloch’s Theorem), and it also illustrateshow boundary conditions are paramount in determiningthe energy spectrum. The energies themselves are read-ily determined through an explicit equation, and can bedone with a calculator. Even more importantly, perhaps,is that in one particular limit of the model, the so-calledDirac comb, where the barriers are made to be infinitelythin while their heights are made infinitely high (i.e. a se-ries of periodic δ -functions), the existence of energy gapscan be demonstrated analytically.All of these features are beneficial and even desir-able for pedagogical reasons, but often one would like todemonstrate the same ideas for more complicated (andrealistic) potentials. In other words, it naturally occursto students at some point to inquire about the generalityof the conclusions reached through a model that considersonly the simplest possible case. Furthermore, for gradu-ate students there tends to be a quantum leap from theanalytical solution of the Kronig-Penny model to the nu-merical (black box-like) solutions inherent in electronicband structure calculations. In this paper we propose tobridge this gap somewhat with simple, albeit numerical,solutions for an arbitrary unit cell potential, though still ∗ Electronic address: [email protected] † Electronic address: [email protected] in one dimension.The basic idea comes from some recent work concern-ing very simple one-body potentials.[2, 3] We will de-scribe this work in the next section, and follow it with thealterations necessary to utilize periodic (as opposed toopen) boundary conditions. The purpose of this changeis it allows us to follow up with an implementation ofthe Bloch condition, as opposed to the periodic boundarycondition. We can then compare our numerically derivedresults with analytic ones, to establish the validity andaccuracy of the method. Note that an alternate proce-dure would be to simply define a periodic potential thatrepeats a finite number of times, say ten or twelve, andsolve the Schr¨odinger equation for this potential. Thiswas done, for example, in Ref. [4], by solving the differ-ential equation using a ‘hunt and shoot’ method. Here weprefer to use the matrix diagonalization method,[2] andwe prefer to incorporate Bloch’s Theorem analytically,for the sake of efficiency and pedagogy. In the followingsection we explore a number of fairly different lookingunit cell potentials; freedom to choose the parametersthat characterize different shaped potentials makes eachpotential able to have various bandwidths and effectivemasses. One particular feature that is shape-dependentis the asymmetry in electron vs. hole effective mass. Weconclude with a summary.
II. FORMALISMA. Infinite square well
By embedding a potential in the infinite square well, oneis able to study the bound states of that potential withminimal technical overhead. One has to have access to anumerical matrix diagonalization routine; otherwise onlya freshman knowledge of integral calculus and linear al-gebra is required. The methodology is outlined in Ref. 2so only key points will be recalled here. One first expandsthe wave function in the infinite well basis: | ψ i = ∞ X m =1 c m | ψ m i (1)where the basis states are ψ n ( x ) = q a sin (cid:0) nπa x (cid:1) < x < a, E (0) n = n π ~ ma ≡ n E (0)1 , (3)where a is the width of the well, and the (unknown) co-efficients are c n . Straightforward algebra leads to thematrix diagonalization problem, ∞ X m =1 H nm c m = Ec n (4)where H nm = h ψ n | ( H + V ) | ψ m i = δ nm E (0) n + H Vnm (5)and H Vnm = h ψ n | V ( x ) | ψ m i = 2 a Z a d x sin (cid:18) nπxa (cid:19) V ( x ) sin (cid:18) mπxa (cid:19) . (6)The problem is easily rendered into a convenient dimen-sionless form by using the infinite square well width a and ground state energy E (0)1 : ∞ X m =1 h nm c m = ec n , (7)where h nm ≡ H nm /E (0)1 and e ≡ E/E (0)1 . B. Square well with periodic boundary conditions
In this work, as a first step, we use the same strategy,but instead of a box with infinite walls, we will use abox with periodic boundary conditions. Then the wavefunction satisfies periodic boundary conditions: φ ( x + a ) = φ ( x ) , (8)with solutions that are the plane wave states, φ ( x ) ∼ e ikx , (9)where k ≡ mE/ ~ . Here, k can be either negative orpositive. Imposition of the boundary conditions in Eq. 8then requires ka = 2 nπ, (10) where n is an integer: n = ... − , − , , , , .... The eigen-values are then E n = 4 n π ~ ma ! = 4 n E (0)1 . (11)Note that these differ from the eigenvalues of an infinitesquare well: (i) they have a value of 4 × those eigenvalues,for the same integer, n , but (ii) they are doubly degen-erate, and (iii) E = 0 is possible, and also constitutesthe one exception to point (ii) ( n = 0). The orthonormalbasis states are then φ n ( x ) = r a exp (cid:20) i πna x (cid:21) , (12)where n is an integer. III. HARMONIC OSCILLATOR
The first question one should ask is, does this basis work?And if so, is this basis any better or worse than the infi-nite square well basis? To this end we recall the resultsobtained in Ref. [2], for a harmonic oscillator potentialembedded in the infinite square well.
A. Infinite square well
We placed the harmonic oscillator potential inside an in-finite square well of width a that spanned 0 < x < a . Indimensionless form the potential was written v HO = V HO E (0)1 = π ~ ωE (0)1 ! (cid:18) xa − (cid:19) . (13)Using this potential, the dimensionless Hamiltonian ma-trix components are h nm = H nm E (0)1 = δ nm n + π ~ ωE (0)1 ! (cid:18) − πn ) (cid:19) +(1 − δ nm ) ~ ωE (0)1 ! g mn , (14)where g mn = ( − n + m + 14 ! (cid:18) n − m ) − n + m ) (cid:19) . (15)A diagonalization of this matrix with a truncated basisindeed converged to the correct ground state wave func-tion, and the coefficients in Eq. (1) indeed agreed to highaccuracy with those expected from the exact analyticalsolution:[5] c n = i ( n − (cid:18) π E (0)1 ~ ω (cid:19) / exp (cid:20) − n E (0)1 ~ ω (cid:21) for n odd , n even.(16) B. Periodic boundary conditions
With periodic boundary conditions, we use − a < x < a ,and the potential is given by V ( x ) = 12 mω x for − a < x < a H Vnm = h φ n | V ( x ) | φ m i = 1 a Z a/ − a/ d x e − i πnxa (cid:18) mω x (cid:19) e i πmxa . (18)In dimensionless form this becomes h Vnm ≡ H Vnm E (0)1 = ~ ωE (0)1 ! π I, (19)where I = Z / − / u e i π ( m − n ) u d u. (20)This integral is straightforward; when combined with thekinetic energy contribution, one obtains h nm = δ nm n + π ~ ωE (0)1 ! + (1 − δ nm ) ~ ωE (0)1 ! ( − m − n ( m − n ) , (21)which is clearly different from Eq. 14.As in the infinite square well case we will have totruncate this matrix, so it makes sense to arrangethe basis states in order of increasing energy. Wewill use an ordering of quantum numbers as the se-ries { , , − , , − , , − , ... } , and then truncate at some n max . We can also anticipate our numerical results bycalculating the Fourier coefficients of the harmonic oscil-lator ground state analytically. Using the analytical fromfor the ground state wave function written in the samedimensionless units as above, ψ HO ( x ) = π a ~ ωE (0)1 ! / exp − π ~ ωE (0)1 ! (cid:18) xa (cid:19) , (22) n E n / E ( ) n + v off h ω /2 π (n−1/2)/E Infinite square wellPeriodic boundary
Fig 1: Plot of the normalized energy levels vs. quantumnumber n for numerical solutions using open (infinitesquare well) and periodic boundary conditions. We used n max = 60 and ~ ω/E (0)1 = 20, and note that ourquantum number n begins at unity (not zero). Therationale for the curve n + v off is explained in Ref. [2].Both results give the correct energy eigenvalues at lowvalues of n , and both grow as n at large values of n ,reflecting the expected behaviour due to the confiningbox. Note the clear degeneracies at large n for the casewith periodic boundary condition energies.we can take the inner product with each eigenstate fromEq. 12. We obtain c n = h φ n ( x ) ψ HO ( x ) i = π E (0)1 ~ ω ! / exp − n E (0)1 ~ ω ! (23)which again differ from those in Eq. (16). IV. NUMERICAL COMPARISON
We pause for a moment to compare the accuracy andefficiency of the two basis sets. In Fig. 1 we plot thenumerical eigenenergies obtained for the two cases with n max = 60. Both methods give very similar results.In particular, either result reproduces the expected lin-ear in n results to high accuracy for low lying states,where the presence of the box is not even noticed. Thisis the regime that is physically relevant to the harmonicoscillator potential. At higher values of n both casescases produce energies that grow as n . To evaluate theefficiency, we show the results for the coefficients on a logplot in Figure 2; both methods yield similar results, withthe smooth decay breaking down where n ≈ ~ ω/E (0)1 ,which is the condition for the basis state energy equal to −10 −8 −6 −4 −2 n | c n | Infinite square wellPeriodic boundary
Fig 2: The numerically-derived absolute value of theFourier coefficients for the infinite square well and theperiodic boundary condition embedding potentials,plotted on a semilog graph. The even cases for theinfinite square well have finite nonzero results on theorder of less than 10 − and are not shown. Here, thetruncated matrix dimension is n max = 60 and( ~ ω/E (0)1 ) = 20. Analytical results cannot bedistinguished from the numerical results.the crossover (from harmonic oscillator to square well)energy for the potential. Exact analytical results agreeto 7 digits in either case. Note that agreement can besystematically improved indefinitely in both cases by en-larging the width of the embedding square well. V. BLOCH’S THEOREM
Bloch’s Theorem[7] states that the electronic wave func-tion in a periodic potential can be written as a productof a plane wave and a function which has the same peri-odicity as the potential. Alternatively one can write ψ ( x + a ) = e iKa ψ ( x ) (24)where a is the unit cell length and K is a wave vectorsuch that − π < Ka < π . A. Analytical solution to the Kronig-Penney model
In practice, the analytical treatment of the Kronig-Penney model[1] allows one to utilize Bloch’s Theoremas a boundary condition. In that problem, one uses apotential that repeats a square well of width b followed by a square barrier of width a − b : V ( x ) = V ∞ X n = −∞ θ (cid:2) x − ( na + b ) (cid:3) θ (cid:2) ( n + 1) a − x (cid:3) , (25)where θ [ x ] is the Heaviside step function. Here, the over-all repeat distance is a . For states whose energies arebelow the barrier height one can solve the Schr¨odingerequation for each region in turn: a linear combination oftwo plane waves in the ‘well’ regions, followed by a linearcombination of an exponentially decaying and an expo-nentially growing solution in the ‘barrier regions. Match-ing the wave functions and their derivatives at the inter-face determines two of the unknown coefficients; at thenext interface, a similar procedure determines two morecoefficients, but this still leaves two unknown coefficients,and on this would go ad infinitum. Bloch’s Theorem al-lows the second matching process to terminate the proce-dure, since now the two coefficients are written in termsof the original two, and we are left with four homogeneousequations with four unknowns. Since the determinant ofthe coefficients of these four equations must be zero forthere to be a solution, this results in the expressioncos ( Ka ) = cos ( k b ) cosh [ κ ( a − b )]+ κ − k k κ sin ( k b ) sinh [ κ ( a − b )] (26)where k = p mE/ ~ and κ = p m ( V − E ) / ~ .Equation (26) is thus an implicit equation for E ( K ). Inpractice, one selects a value of E ; the absolute value ofthe right-hand-side of Eq. (26) is evaluated — it is eithergreater or less than unity. If greater, then there is clearlyno solution possible, while if less than unity, then takingthe inverse cosine of this quantity gives the value of Ka for which this energy is the solution. Plots will be shownlater when comparisons are made with the numerical re-sults. B. Matrix method for the Kronig-Penney model
It should be clear from Eq. (24) why we altered theoriginal matrix method to include the case of periodicboundary conditions; now we simply replace Eq. (8) withEq. (24). How does this alter the procedure discussed inthe previous section? Quite simply, Eq. (10) is replacedby the one required by Bloch’s Theorem: ka = 2 πn → ka − Ka = 2 πn (27)where, as before, k = 2 mE/ ~ . Hence E (0) n = ~ π ma (cid:18) n + Kaπ (cid:19) = E (0)1 (cid:18) n + Kaπ (cid:19) . (28)Crucially, Eq. 28 applies only to the diagonal elements;the off-diagonal elements are unaffected, because the ba-sis states are modified just as φ n ( x ) = 1 √ a e i πna x → √ a e i πn + Kaa x = 1 √ a e i πna x e + iKx = e + iKx φ n ( x ) . (29)But this means that φ ∗ n → e − iKx φ ∗ n . Thus the extraexponentials arising from the Bloch condition cancel oneanother when we calculate H Vnm and so the off-diagonalterms remain purely potential dependent.To generate bands of energy, different values of Ka inthe region ( − π, π ) are passed to matrices suitably mod-ified by Eq. 28. The remaining matrix elements are de-termined as before. For the case of a well of width b centred between two barriers of height V , each of width( a − b ) / ρ ≡ b/a and v ≡ V /E (0)1 , we have h nm = δ nm "(cid:18) n + Kaπ (cid:19) + v (1 − ρ ) (30)+(1 − δ nm ) v ( − m − n +1 π sin [ π ( m − n ) ρ ] m − n . We repeatedly diagonalize matrices of this form for var-ious values of Ka ∈ ( − π, π ) to form the band solutions. x/a V ( x ) / E ( ) Fig 3: Central square well with v = 1 and ρ = 0 . v = V /E (0)1 = 0) we see the foldedparabolas, corresponding to plane wave energy solutions.As we turn the periodic potential on, e.g. with v = 10,we see bandgaps emerge, with larger gaps in the lowerenergy bands, and more curvature (i.e. bandwidth) inthe upper bands. Analytical solutions found by solvingEq. 26 are overlaid with the numerical ones, and are in −1 −0.5 0 0.5 10510152025 Ka/ π E / E ( ) (a)v = 0 −1 −0.5 0 0.5 1051015202530 Ka/ π E / E ( ) (b)v = 10 Fig 4: Solutions to the Kronig-Penney potential with(a) no well/barrier, and (b) a repeated well/barriersequence with barrier height v = 10. In (b) we havealso indicated with open circles some of the analyticalsolutions to Eq. (26). A schematic representation of thepotential is shown to indicate the boundedness of theenergy bands. The matrix solutions utilize n max = 60,and we have used ρ = b/a = 0 . C. General potential shapes
It should be clear that we now have freedom to studyperiodic arrays of potentials with any shape, simply byperforming many calculations (corresponding to differentvalues of Ka ) for the wave function in one unit cell . See (a) Kronig−Penney, ρ = 0.5(b) Kronig−Penney, ρ = 0.8(c) Simple harmonic oscillator V ( x ) (d) Inverted harmonic oscillator(e) Linear well Fig 5: Schematic representation of the potentials usedfor comparing energy band structures.Fig. 5 for various examples of periodic potentials, butnote that they do not necessarily have to have an analyt-ical form, nor do they necessarily require an analyticalintegration for the matrix elements. For example, wecan repeat harmonic oscillator wells. Then the requireddimensionless matrix elements are h nm = δ nm (cid:18) n + Kaπ (cid:19) + π ~ ωE (0)1 ! + (1 − δ nm ) ~ ωE (0)1 ! ( − m − n ( m − n ) . (31)For the inverted harmonic oscillator potential (seeFig. 5) we use V ( x ) = − mω h x − a i < x < a , − mω h ( x − a ) − a i a < x < a. (32)The matrix elements are then h nm = δ nm (cid:18) n + Kaπ (cid:19) + π ~ ωE (0)1 ! − (1 − δ nm ) ~ ωE (0)1 ! ( − m − n ( m − n ) . (33)Note the difference between the above equation andEq. 31 in both the off-diagonal elements (change of sign)and the diagonal elements (factor of 2 in the potentlalterm). A third example is the linear well (resembles a saw-tooth in Fig. 5), with potential defined by V ( x ) = ( A (cid:0) − xa (cid:1) < x < a , A (cid:0) xa − (cid:1) a < x < a. (34)This has matrix elements h nm = δ nm "(cid:18) n + Kaπ (cid:19) + A − (1 − δ nm ) Aπ ( m − n ) (cid:2) − ( − m − n (cid:3) . (35)As mentioned earlier, potentials can also be utilizedfor which the matrix elements have no known analyti-cal solution, or in which obtaining the analytical solu-tion would be overly cumbersome. One such potential isthe so-called “pseudo-Coulomb” potential which has theform v ( x ) = V ( x ) E (0)1 = − A p ( x − a ) + b (36)with A a positive number representing the strength, or al-ternatively, the inverted pseudo-Coulomb potential with A negative, and b is a small numerical factor introducedto prevent singularities (the true Coulomb potential isrecovered as b → D. Comparing band structures
Results for the band structures corresponding to the po-tential shapes in Fig. 5 are shown in Figs. 6, 7, 8, 9, 10.We chose parameters such that in all cases three bandswith energies less than the maximum barrier heightwould form, and in which the third band would have anenergy difference between the highest level (at Ka = ± π )and the maximum barrier potential, V max , of ∆ E = E (0)1 in each case.Note that bands with gaps form at energies above thebarrier maxima as well. However, the character of thesebands is strongly dependent on the type of periodic po-tential used. In particular, for the latter three potentialsthe gap can be quite small (not discernible, for example,on the scale of Fig. 9). Moreover for these latter three po-tentials, the minima and maxima of these higher energybands can be distinctly non-parabolic, and in fact exhibitV-shaped (or inverted V-shaped) dispersions close to theminima (or maxima).We then focused on the third band (the one closestto but lower than V max in each case and normalized thebands by setting the highest band level to E = 0 in eachcase, as shown in Fig. 11. The bandwidth varies consid-erably as shown; this can be adjusted by varying the bar-rier heights and widths, as is explicitly shown in the caseof the Kronig-Penney model (first two figures in Fig. 5);in the case of the other potentials we could include an −1 −0.5 0 0.5 10510152025303540 Ka/ π E / E ( ) Fig 6: Energy band diagram for the Kronig-Penneypotential with ρ = 0 . v = 20 . −1 −0.5 0 0.5 1051015202530 Ka/ π E / E ( ) Fig 7: Energy band diagram for the Kronig-Penneypotential with ρ = 0 . v = 10 . m ∗ elehol = 1 ~ d ǫ d K (cid:12)(cid:12)(cid:12)(cid:12) K min K max . (37)The second derivative in Eq. (37) is evaluated with a five-point fit,[9] and the ratio of the two effective masses foreach potential shape is tabulated in Table III. The two −1 −0.5 0 0.5 1051015202530 Ka/ π E / E ( ) Fig 8: Energy band diagram for the simple harmonicoscillator potential with ~ ω/E (0)1 = 4 . ǫ ′′ ele /ǫ ′′ hol ≡ m ∗ hol /m ∗ ele . Potential ǫ ′′ ele ǫ ′′ hol ǫ ′′ ele /ǫ ′′ hol K-P ( ρ = 0 .
5) 13.83 -25.35 -0.55K-P ( ρ = 0 .
8) 39.09 -70.61 -0.55Simple HO 37.84 -121.80 -0.31Inverted HO 19.83 -55.96 -0.35Linear 31.63 -102.23 -0.31
Kronig-Penney models yield essentially the same ratio,and the magnitude of this ratio can decrease consider-ably when other potential shapes are considered. Firstnote that the ratio has a magnitude that is less thanunity. This is because holes have effectively a weaker bar-rier through which to tunnel, compared with electrons.This decrease is most pronounced when either the lin-ear or simple harmonic oscillator potential is used. Thereason is that these have cusp-like barriers, so that thebarrier width is also reduced for holes compared to elec-trons; therefore the holes should have more mobility (i.e.lower effective mass) compared with electrons, over andabove the advantage already present due to the differencein effective barrier height. This asymmetry concerningelectrons and holes may be important in a new class ofsuperconducting models, where dynamic interactions aretaken into account.[10, 11] −1 −0.5 0 0.5 105101520253035404550
Ka/ π E / E ( ) Fig 9: Energy band diagram for the inverted harmonicoscillator potential with ~ ω/E (0)1 = 7 . −1 −0.5 0 0.5 10510152025303540 Ka/ π E / E ( ) Fig 10: Energy band diagram for the linear wellpotential with A = 19 . VI. SUMMARY
We hope that this paper will serve as a bridge from theclassroom to the research desk, for both undergraduateand graduate students, in the area of electronic structurecalculations. The ‘classroom’ or ‘textbook’ example forelectronic band structure has always been the Kronig-Penney model. The enlightening feature in the modelhas been the recasting of the problem of an extended (i.e.Bloch) state in terms of the problem in one unit cell, due to Bloch’s Theorem. The part that is more intimidating −1 −0.5 0 0.5 1−3−2.5−2−1.5−1−0.50
Ka/ π E / E ( ) K−P ( ρ = 0.5)Inv. HOK−P ( ρ = 0.8)Linear SHO Fig 11: The third energy band for the five cases inFig. (5), with the topmost energy set equal to zero.Note the variation in bandwidth, but, more importantly,the difference in curvatures, as tabulated in Table III.for students is the extension beyond simple square wellsto more realistic potentials; this is what this paper hastried to address.Simple procedures to obtain numerical solutions forany potential that supports bound states were illustratedin Ref. [2]. This methodology was extended to peri-odic boundary conditions for the purpose of this paper,because the next step, implementing the Bloch condi-tion, is then essentially trivial. Students can then tackletheir own favourite potential, for example the “pseudo-Coulomb” potential referred to in Eq. (36). We havegiven a number of examples already in this paper, andmentioned one particular property—the asymmetry be-tween electron and hole masses, which can have impor-tant consequences for semiconductors and possibly super-conductors. In general the more realistic potentials tendto make the barrier effectively narrower for hole thanfor electrons, with the consequence that hole masses aregenerally lower than electron masses. We also note thatmore realistic potentials can give rise to higher energybands with non-parabolic shapes.
Acknowledgments
This work was supported in part by the Natural Sciencesand Engineering Research Council of Canada (NSERC),by the Alberta iCiNano program, and by a Universityof Alberta Teaching and Learning Enhancement Fund(TLEF) grant. [1] R. de L. Kronig; W. G. Penney, “Quantum Mechanicsof Electrons in Crystal Lattices,” Proc. R. Soc. Lond. A, , 499 (1931).[2] F. Marsiglio, “The harmonic oscillator in quantum me-chanics: A third way, “Am. J. Phys. , 253 (2009).[3] V. Jelic and F. Marsiglio, “The double-well potential inquantum mechanics: a simple, numerically exact formu-lation,” Eur. J. Phys. , 600(1992).[5] See Ref. [2]. Note that this reference had a typo in Eq.(15), which has been corrected here in Eq. (16).[6] In either case, for the harmonic oscillator potential, orany potential that is even in x , we could of course in-crease the efficiency by a factor of two, by leaving outthe even integer basis functions for the infinite squarewell case (since they are zero) and using only the sumof the positive and negative coefficients for the periodicboundary condition case. In a more general case neither of these efficiencies is possible.[7] F. Bloch, “ ¨Uber die Quantenmechanik der Elektronen inKristallgittern,” Z. Physik. , 545 (1929), in German.Remarkably, he has 12 papers that are more cited thanthis one. See also F. Bloch, ‘Memories of electrons incrystals’, Proc. R. Soc. Lond. A , 24 (1980), wherehe provides some context for this Theorem.[8] D. J. Griffiths, Introduction to Quantum Mechanics , 2nded. (Pearson/Prentice Hall, Upper Saddle River, NJ,2005), pp. 224–229.[9] T. Pang,
An Introduction to Computational Physics ,2nd ed. (Cambridge University Press, Cambridge, 2006),p. 51.[10] J. E. Hirsch, “Dynamic Hubbard Model,” Phys. Rev.Lett. , 206402-1-4 (2001).[11] G. H. Bach, J. E. Hirsch and F. Marsiglio, “Two-sitedynamical mean field theory for the dynamic Hubbardmodel,” Phys. Rev. B82