Harmonic oscillator eigenfunction expansions, quantum dots, and effective interactions
aa r X i v : . [ phy s i c s . c o m p - ph ] M a y Harmonic oscillator eigenfunction expansions, quantum dots, and effective interactions
Simen Kvaal ∗ Centre of Mathematics for Applications, University of Oslo, N-0316 Oslo, Norway (Dated: June 6, 2018)We give a thorough analysis of the convergence properties of the configuration-interaction methodas applied to parabolic quantum dots among other systems, including a priori error estimates. Themethod converges slowly in general, and in order to overcome this, we propose to use an effective two-body interaction well-known from nuclear physics. Through numerical experiments we demonstratea significant increase in accuracy of the configuration interaction method.
I. INTRODUCTION
The last two decades, an ever-increasing amount ofresearch have been dedicated to understanding the elec-tronic structure of so-called quantum dots : semiconduc-tor structures confining from a few to several thousandselectrons in spatial regions on the nanometre scale. Insuch calculations, one typically seeks a few of the lowesteigenenergies E k of the system Hamiltonian H and theircorresponding eigenvectors ψ k , i.e., Hψ k = E k ψ k , k = 1 , · · · , k max . (1.1)One of the most popular methods is the (full) configura-tion interaction method (CI), where the many-body wavefunction is expanded in a basis of eigenfunctions of theharmonic oscillator (HO), and then necessarily truncatedto give an approximation. In fact, the so-called curse ofdimensionality implies that the number of degrees of free-dom available per particle is severely limited. It is clear,that an understanding of the properties of such basis ex-pansions is very important, as it is necessary for a priori error estimates of the calculations. Unfortunately, this isa neglected topic in the physics literature.In this article, we give a thorough analysis of the (full)configuration interaction (CI) method using HO expan-sions applied to parabolic quantum dots, and give practi-cal convergence estimates. It generalizes and refines thefindings of a recent study of one-dimensional systems, and is applicable to for example nuclear systems andquantum chemical calculations as well. We demonstratethe estimates with calculations in the d = 2 dimensionalcase for N ≤ .The main results are however somewhat discouraging.The expansion coefficients of typical eigenfunctions areshown to decay very slowly, limiting the accuracy of any practical method using HO basis functions. We thereforepropose to use an effective two-body interaction to over-come, at least partially, the slow convergence rate. Thisis routinely used in nuclear physics where the inter-particle forces are of a completely different, and basicallyunknown, nature. For electronic systems, however, theinteraction is well-known and simpler to analyze, but ef-fective interactions of the present kind have not been ap-plied, at least to the author’s knowledge. The modifiedmethod is seen to have convergence rates of at least one order of magnitude higher than the original CI method.An important point here is that the complexity of theCI calculations is not altered, as no extra non-zero ma-trix elements are introduced. All one needs is a relativelysimple one-time calculation to produce the effective in-teraction matrix elements.The HO eigenfunctions are popular for several reasons.Many quantum systems, such as the quantum dot modelconsidered here, are perturbed harmonic oscillators perse , so that the true eigenstates should be perturbationsof the HO states. Moreover, the HO has many beautifulproperties, such as complete separability of the Hamil-tonian, invariance under orthogonal coordinate changes,and thus easily computed eigenfunctions, so that comput-ing matrix elements of relevant operators becomes rela-tively simple. The HO eigenfunctions are defined on the whole of R d in which the particles live, so that trunca-tion of the domain is unnecessary. Indeed, this is one ofthe main problems with methods such as finite differenceor finite element methods. On the other hand, the HOeigenfunctions are the only basis functions with all theseproperties.The article is organized as follows. In Sec. II we discussthe harmonic oscillator and the the parabolic quantumdot model, including exact solutions for the N = 2 case.In Sec. III, we give results for the approximation proper-ties of the Hermite functions in n dimensions, and thusalso of many-body HO eigenfunctions. By approximationproperties, we mean estimates on the error k ψ − P ψ k ,where ψ is any wave function and P projects onto a fi-nite subspace of HO eigenfunctions, i.e., the model space.Here, P ψ is in fact the best approximation in the norm.The estimates will depend on analytic properties of ψ ,i.e., whether it is differentiable, and whether it falls ofsufficiently fast at infinity. To our knowledge, these re-sults are not previously published.In Sec. IV, we discuss the full configuration interactionmethod, using the results obtained in Sec. III to obtainconvergence estimates of the method as function of themodel space size. We also briefly discuss the effectiveinteraction utilized in the numerical calculations, whichare presented in Sec. V. We conclude with a discussion ofthe results, its consequences, and an outlook on furtherdirections of research in Sec. VI.We have also included an appendix with proofs of theformal propositions in Sec. III. II. THE HARMONIC OSCILLATOR ANDPARABOLIC QUANTUM DOTSA. The Harmonic Oscillator
A spinless particle of mass m in an isotropic harmonicpotential has Hamiltonian H HO = − ~ m ∇ + 12 mω k ~r k , (2.1)where ~r ∈ R d is the particle’s coordinates. By choosingproper energy and length units, i.e., ~ ω and p ~ /mω ,respectively, the Hamiltonian becomes H HO = − ∇ + 12 k ~r k . (2.2) H HO can be written as a sum over d one-dimensionalharmonic oscillators, viz, H HO = d X k =1 (cid:18) − ∂ ∂r k + 12 r k (cid:19) , (2.3)so that a complete specification of the HO eigenfunctionsis given byΦ α ,α ,...,α d ( ~r ) = φ α ( r ) φ α ( r ) · · · φ α d ( r d ) , (2.4)where φ α i ( x ), α i = 0 , , . . . are one-dimensional HOeigenfunctions, also called Hermite functions. These aredefined by φ n ( x ) = (2 n n ! π / ) − / H n ( x ) e − x / , n = 0 , , . . . , (2.5)where the Hermite polynomials H n ( x ) are given by H n ( x ) = ( − n e x ∂ n ∂x n e − x . (2.6)The Hermite polynomials also obey the recurrence for-mula H n +1 ( x ) = 2 xH n ( x ) − nH n − ( x ) , (2.7)with H ( x ) = 1 and H ( x ) = 2 x . The Hermite poly-nomial H n ( x ) has n zeroes, and the Gaussian factor in φ n ( x ) will eventually subvert the polynomial for large | x | . Thus, qualitatively, the Hermite functions can bedescribed as localized oscillations with n nodes and aGaussian “tail” as x approaches ±∞ . One can easilycompute the quantum mechanical variance(∆ x ) := Z ∞−∞ x φ n ( x ) d x = n + 12 , (2.8)showing that, loosely speaking, the width of the oscilla-tory region increases as ( n + 1 / / .The functions Φ α , ··· ,α d defined in Eqn. (2.4) are called d -dimensional Hermite functions. In the sequel, we willdefine α = ( α , · · · , α d ) ∈ I d for a tuple of non-negative integers, also called a multi-index; see Appendix A 1. Us-ing multi-indices, we may writeΦ α ( ~r ) = (cid:16) | α | α ! π d/ (cid:17) − / H α ( r ) · · · H α d ( r d ) e −k ~r k / . (2.9)The eigenvalue of φ n ( x ) is n + 1 /
2, so that the eigen-value of Φ α ( ~r ) is ǫ α = d | α | , (2.10)i.e., a zero-point energy d/ | α | the shell number of Φ α , andthe eigenspace S r ( R d ) corresponding to the eigenvalue d/ r a shell . We define the shell-truncated Hilbertspace P R ( R d ) ⊂ L ( R d ) as P R ( R d ) := span (cid:8) Φ α ( ~r ) (cid:12)(cid:12) | α | ≤ R (cid:9) = R M r =0 S r ( R d ) , (2.11)i.e., the subspace spanned by all Hermite functions withshell number less than or equal to R , or, equivalently,the direct sum of the shells up to and including R . The N -body generalization of this space, to be discussed inSection III B, is a very common model space used in CIcalculations.Since the Hermite functions constitute an orthonormalbasis for L ( R d ), P R ( R d ) → L ( R d ), in the sense that forevery ψ ∈ L ( R d ), lim R →∞ k ψ − P ψ k = 0, where P isthe orthogonal projector on P R ( R d ). Strictly speaking,we should use a symbol like P R or even P R ( R d ) for theprojector. However, R and d will always be clear from thecontext, so we are deliberately sloppy to obtain a conciseformulation. For the same reason, we will sometimessimply write P or P R for the space P R ( R d ).An important fact is that since H HO is invariant underorthogonal spatial transformations (i.e., such transforma-tion conserve energy), so is each individual shell space.Hence, each shell S r ( R d ), and also P R ( R d ), is indepen-dent of the spatial coordinates chosen.For the case d = 1 each shell r is spanned by a singleeigenfunction, namely φ r ( x ). For d = 2, each shell r hasdegeneracy r + 1, with eigenfunctionsΦ ( s,r − s ) ( ~r ) = φ s ( r ) φ r − s ( r ) , ≤ s ≤ r. (2.12)The usual HO eigenfunctions used to construct many-body wave functions are not the Hermite functionsΦ α , ··· ,α d , however, but rather those obtained by utilizingthe spherical symmetry of the HO. This gives a many-body basis diagonal in angular momentum. For d = 2we obtain the so-called Fock-Darwin orbitals given byΦ FD n,m ( r, θ ) = (cid:20) n !( n + | m | )! (cid:21) / e imθ √ π L | m | n ( r ) e − r / . (2.13)Here, n ≥ m is the azimuthal quantumnumber. The eigenvalues are ǫ n,m = 2 n + | m | + 1 . (2.14)Thus, R = 2 n + | m | is the shell number. By construc-tion, the Fock-Darwin orbitals are eigenfunctions of theangular momentum operator L z = − i∂/∂ θ with eigen-value m . Of course, we may write Φ FD n,m as a linearcombination of the Hermite functions Φ ( s,R − s ) , where0 ≤ s ≤ R = 2 n + | m | , and vice versa. The actualchoice of form of eigenfunctions is immaterial, as long aswe may identify those belonging to a given shell.The space P R =4 ( R ) is illustrated in Fig. 1 using bothHermite functions and Fock-Darwin orbitals. B. Parabolic quantum dots
We consider N electrons confined in a harmonic oscil-lator in d dimensions. This is a very common model for aquantum dot. We comment, that modelling the quantumdot geometry by a perturbed harmonic oscillator is justi-fied by self-consistent calculations, and is a widelyadopted assumption. The Hamiltonian of the quantum dot is given by H := T + U, (2.15)where T is the many-body HO Hamiltonian, given by T = N X k =1 H HO ( ~r k ) (2.16)and U is the inter-electron Coulomb interactions. In di-mensionless units the interaction is given by, U := N X i Before we discuss the approximation properties of theHermite functions, it is instructive to consider the verysimplest example of a two-electron parabolic quantum dot and the properties of the eigenfunctions, since thiscase admits analytical solutions for special values of λ andis otherwise well understood . Here, we consider d = 2 dimensions only, but the d = 3 case is similar. Wenote, that for N = 2 it is enough to study the spatialwave function, since it must be either symmetric (for thesinglet S = 0 spin state) or anti-symmetric (for the triplet S = 1 spin states). The Hamiltonian (2.15) becomes H = − 12 ( ∇ + ∇ ) + 12 ( r + r ) + λr , (2.19)where r = k ~r − ~r k and r j = k ~r j k . Introduce a setof scaled centre of mass coordinates given by ~R = ( ~r + ~r ) / √ ~r = ( ~r − ~r ) / √ 2. This coordinate changeis orthogonal and symmetric in R . This leads to theseparable Hamiltonian H = − 12 ( ∇ r + ∇ R ) + 12 ( k ~r k + k ~R k ) + λ √ k ~r k = H HO ( ~R ) + H rel ( ~r ) . A complete set of eigenfunctions of H can now be writtenon product form, viz,Ψ( ~R, ~r ) = Φ FD n ′ ,m ′ ( ~R ) ψ ( ~r ) . (2.20)The relative coordinate wave function ψ ( ~r ) is an eigen-function of the relative coordinate Hamiltonian given by H rel = − ∇ r + 12 r + λ √ r , (2.21)where r = k ~r k . This Hamiltonian can be further sepa-rated using polar coordinates, yielding eigenfunctions onthe form ψ m,n ( r, θ ) = e imθ √ π u n,m ( r ) , (2.22)where | m | ≥ u n,m is an eigenfunctionof the radial Hamiltonian given by H r = − r ∂∂r r ∂∂r + | m | r + 12 r + λ √ r . (2.23)By convention, n counts the nodes away from r = 0 of u n,m ( r ). Moreover, odd (even) m gives anti-symmetric(symmetric) wave functions Ψ( ~r , ~r ). For any given | m | ,it is quite easy to deduce that the special value λ = p | m | + 1 yields the eigenfunction u ,m = Dr | m | ( a + r ) e − r / , (2.24)where D and a are constants. The corresponding eigen-value of H r is E r = | m | + 2, and E = 2 n ′ + | m ′ | + 1 + E r .Thus, the ground state (having m = m ′ = 0, n = n ′ = 0)for λ = 1 is given byΨ ( ~R, ~r ) = D ( r + a ) e − ( r + R ) / = D √ r + √ a ) e − ( r + r ) / , n = 0 m = 0 n = 0 m = − n = 0 m = 1 n = 0 m = − n = 1 m = 0 n = 0 m = 2 n = 0 m = − n = 1 m = − n = 1 m = 1 n = 0 m = 3 n = 0 m = − n = 1 m = − n = 2 m = 0 n = 1 m = 2 n = 0 m = 4 α = (0 , α = (0 , α = (1 , α = (0 , α = (1 , α = (2 , α = (0 , α = (1 , α = (2 , α = (3 , α = (0 , α = (1 , α = (2 , α = (3 , α = (4 , o shell S Unitarily equivalent FIG. 1: Illustration of P R =4 ( R ): (Left) Fock-Darwin orbitals. (Right) Hermite functions. Basis functions with equal HOenergy are shown at same line. with D being a (new) normalization constant.Observe that this function has a cusp at r = 0, i.e.,at the origin x = y = 0 (where we have introducedCartesian coordinates ~r = ( x, y ) for the relative coordi-nate). Indeed, the partial derivatives ∂ x ψ , and ∂ y ψ , are not continuous there, and Ψ has no partial deriva-tives (in the distributional sense, see Appendix A 2) ofsecond order. The cusp stems from the famous “cuspcondition” which in simple terms states that, for a non-vanishing wave function at r = 0, the Coulomb diver-gence must be compensated by a similar divergence in theLaplacian. This is only possible if the wave functionhas a cusp.On the other hand, the non-smooth function Ψ ( ~R, ~r )is to be expanded in the HO eigenfunctions, e.g., Fock-Darwin orbitals. (Recall, that the particular represen-tation for the HO eigenfunctions are immaterial – alsowhether we use lab coordinates ~r , or centre-of-mass co-ordinates ~R and ~r , since the coordinate change is orthog-onal.) For m = 0, we haveΦ FD n, ( r ) = r π L n ( r ) e − r / , (2.25)using the fact that these are independent of θ . Thus,Ψ ( ~r ) = Φ FD0 , ( R ) u , ( r ) = Φ FD0 , ( R ) ∞ X n =0 c n Φ FD n, ( r ) , (2.26)The functions Φ FD n, ( r ) are very smooth, as is seen by not-ing that L n ( r ) = L n ( x + y ) is a polynomial in x and y , while u , ( r ) = u , ( p x + y ), so Eqn. (2.26) is ba-sically approximating a square root with a polynomial.Consider then a truncated expansion Ψ ,R ∈ P R ( R ),such as the one obtained with the CI or coupled clus-ter method. In general, this is different from P R Ψ ,which is the best approximation of the wave function in P R ( R ). In any case, this expansion, consisting of the R + 1 terms like those of Eqn. (2.26) is a very smoothfunction. Therefore, the cusp at r = 0 cannot be wellapproximated. In Section III C, we will show that the smoothnessproperties of the wave function Ψ is equivalent to a cer-tain decay rate of the coefficients c n in Eqn. (2.26) as n → ∞ . In this case, we will show that ∞ X n =0 n k | c n | < + ∞ , (2.27)so that | c n | = o ( n − ( k +1+ ǫ ) / ) . (2.28)Here, k is the number of times Ψ may be differentiatedweakly, i.e., Ψ ∈ H k ( R ), and ǫ ∈ [0 , 1) is a constant.For the function Ψ we have k = 1. This kind of esti-mate directly tells us that an approximation using onlya few HO eigenfunctions necessarily will give an errordepending directly on the smoothness k .We comment, that for higher | m | the eigenstates willstill have cusps, albeit in the higher derivatives. Indeed,we have weak derivatives of order | m | +1, as can easily bededuced by operating on ψ ,m with ∂ x and ∂ y . Moreover,recall that | m | = 1 is the S = 1 ground state, whichthen will have coefficients decaying faster than the S = 0ground state. Moreover, there will be excited states, i.e.,states with | m | > 1, that also have more quickly decayingcoefficients | c n | . This will be demonstrated numericallyin Sec. V.In fact, Hoffmann-Ostenhof et al. have shown thatnear r = 0, for arbitrary λ any local solution Ψ of( H − E )Ψ = 0 has the formΨ( ξ ) = k ξ k m P (cid:18) ξ k ξ k (cid:19) (1 + a k ξ k ) + O ( k ξ k m +1 ) , (2.29)where ξ = ( ~r , ~r ) ∈ R , and where P , deg( P ) = m ,is a hyper-spherical harmonic (on S ), and where a is a constant. This also generalizes to arbitrary N ,cf. Sec. III D. From this representation, it is manifest,that Ψ ∈ H m +1 ( R ), i.e., Ψ has weak derivatives of or-der m + 1. We discuss these results further in Sec. III D. III. APPROXIMATION PROPERTIES OFHERMITE SERIESA. Hermite functions in one dimension In this section, we consider some formal mathemati-cal propositions whose proofs are given in Appendix A 3,and discuss their importance for expansions in HO basisfunctions.The first proposition considers the one-dimensionalcase, and the second considers general, multidimensionalexpansions. The latter result has to the author’s knowl-edge not been published previously. The treatment forone-dimensional Hermite functions is similar, but notequivalent to, that given by Boyd and Hille. We stress that the results are valid for any given wave-function – not only eigenfunctions of quantum dot Hamil-tonians – assuming only that the wavefunction decays ex-ponentially as | x | → ∞ . In Appendix A 3, more generalconditions are also considered.The results are stated in terms of weak differentiabilityof the wavefunction, which is a generalization of the clas-sical notion of a derivative. The space H k ( R ) ⊂ L ( R )is roughly defined as the (square integrable) functions ψ ( x ) having k (square integrable) derivatives ∂ mx ψ ( x ),0 ≤ m ≤ k . Correspondingly, the space H k ( R n ) ⊂ L ( R n ) consists of the functions whose partial deriva-tives of total order ≤ k are square integrable. Forwavefunctions of electronic systems, it turns out that k times continuous differentiability implies k +1 times weakdifferentiability . The order k of differentiability is notalways known, but an upper or lower bound can oftenbe found through analysis. It is however important, thatthe Coulomb singularity implies that k is finite.For the one-dimensional case, we have the followingproposition: Proposition 1 (Approximation in one dimension) Let k ≥ be a given integer. Let ψ ∈ L ( R ) be expo-nentially decaying as | x | → ∞ and given by ψ ( x ) = ∞ X n =0 c n φ n ( x ) , (3.1) where φ n ( x ) is given by Eqn. (2.5). Then ψ ∈ H k ( R ) ifand only if ∞ X n =0 n k | c n | < ∞ . (3.2)We notice that the latter implies that | c n | = o ( n − ( k +1) / ) , (3.3)which shows that the more ψ ( x ) can be differentiated, thefaster the coefficients will fall off as n → ∞ . Moreover,let ψ R = P R ψ = P Rn =0 c n φ n . Then k ψ − ψ R k = ∞ X n = R +1 | c n | ! / , (3.4) which gives an estimate of how well a finite basis of Her-mite functions will approximate ψ ( x ) in the norm. Wealready notice, that for low k = 2, which is typical, thecoefficients fall off as o ( n − / ), which is rather slowly.In the general n -dimensional case, the wavefunction ψ ∈ L ( R n ) has an expansion in the n -dimensional Her-mite functions Φ α ( x ), α ∈ I n given by ψ ( x ) = X α c α Φ α ( x )= X α ··· α n c α ··· α n φ α ( x ) · · · φ α n ( x n ) . In order to obtain useful estimates on the error, we needto define the shell-weight p ( R ) by the overlap of ψ ( x )with the single shell S R , i.e., p ( R ) = k P ( S R ) ψ k = X α, | α | = R | c α | , (3.5)where P ( S R ) is the projection onto the shell. Thus, k ψ k = ∞ X R =0 p ( R ) . (3.6)For the one-dimensional case, we of course have p ( R ) = | c R | . Proposition 2 (Approximation in n dimensions) Let ψ ∈ L ( R n ) be exponentially decaying as k x k → ∞ and given by ψ ( x ) = X α c α Φ α ( x ) . (3.7) Then ψ ∈ H k ( R n ) if and only if X α | α | k | c α | = ∞ X r =0 r k p ( r ) < + ∞ . (3.8)Again, we notice that the latter implies that p ( r ) = o ( r − ( k +1) ) . (3.9)Moreover, for the shell-truncated Hilbert space P R , theapproximation error is given by k (1 − P ) ψ k = ∞ X r = R +1 p ( r ) ! / . (3.10)In applications, we often observe a decay of non-integral order, i.e., there exists an ǫ ∈ [0 , 1) such thatwe observe p ( r ) = o ( r − ( k +1+ ǫ ) ) . (3.11)This does not, of course, contradict the results. To seethis, we observe that if ψ ∈ H k ( R n ) but ψ / ∈ H k +1 ( R n ),then p ( r ) must decay at least as fast as o ( r − ( k +1) ) butnot as fast as o ( r − k +2 ). Thus, the actual decay exponentcan be anything inside the interval [ k + 1 , k + 2).Consider also the case where ψ ∈ H k ( R n ) for every k , i.e., we can differentiate it (weakly) as many timeswe like. Then p ( r ) decays faster than r − ( k +1) , for any k ≥ 0, giving so-called exponential convergence of theHermite series. Hence, functions that are best approxi-mated by Hermite series are rapidly decaying and verysmooth functions ψ . This would be the case for the quan-tum dot eigenfunctions if the inter-particle interactionswere non-singular. B. Many-body wave functions We now discuss N -body eigenfunctions of the HO in d dimensions, including spin, showing that we may iden-tify the expansion of a such with 2 N expansions in Her-mite functions in n = N d dimensions, i.e., 2 N expan-sions in HO eigenfunctions of imagined spinless particlesin n = N d dimensions. Each expansion corresponds to adifferent spin configuration.Each particle k = 1 , . . . , N has both spatial degreesof freedom ~r k ∈ R d and a spin coordinate τ k ∈ {± } ,corresponding to the z -projection S z = ± ~ of the elec-tron spin. The configuration space can thus be takenas two copies X of R d ; one for each spin value, i.e., X = R d × {± } and x k = ( ~r k , τ k ) ∈ X are the coor-dinates of particle k .For a single particle with spin, the Hilbert space is now L ( X ), with basis functions given byˆΦ i ( x ) = Φ α ( ~r ) χ σ ( τ ) , (3.12)where i = i ( α, σ ) is a new, generic index, and where χ σ is a basis function for the spinor space C .Ignoring the Pauli principle for the moment, the N -body Hilbert space is now given by H ( N ) = L ( X ) N ≡ L ( R Nd ) ⊗ ( C ) N , (3.13)i.e., each wavefunction ψ ∈ H ( N ) is equivalent to2 N spin-component functions ψ ( σ ) ∈ L ( R Nd ), σ =( σ , · · · , σ N ) ∈ {± } N . We have ψ ( x , · · · , x N ) = X σ ψ ( σ ) ( ξ ) χ σ ( τ ) , ξ ≡ ( ~r , · · · , ~r N ) , (3.14)where τ = ( τ , · · · , τ N ), and where χ σ ( τ ) = δ σ,τ are basisfunctions for the N -spinor space C N , being eigenfunc-tions for S z , i.e., corresponding to a given configurationof the N spins.The σ ’th component function ψ ( σ ) ( ξ ) ∈ L ( R Nd ) isa function of N d variables, and by considering the N d -dimensional HO as the sum of N HOs in d dimensions,it is easy to see that a basis for the L ( R Nd ) is given bythe functions Φ β ( ξ ) ≡ Φ α ( ~r ) · · · Φ α N ( ~r N ) (3.15) where ξ = ( ~r , · · · , ~r N ), and where β = ( α , · · · , α N ) isan N d -component multi-index. Correspondingly, a basisfor the complete space H ( N ) = L ( X ) N is given by thefunctionsˆΦ i ··· i N ( ξ, τ ) ≡ Φ α ( ~r ) · · · Φ α N ( ~r N ) χ σ ( τ ) , (3.16)where i k = i ( α k , σ k ). Notice, that the HO energyand hence the shell number | β | only depends on β =( α , · · · , α N ).The functions ψ ( σ ) may be expanded in the functionsΦ β , i.e., ψ ( σ ) ( ξ ) = X β c ( σ ) β Φ β ( ξ )= X α ··· α N c ( σ ) α ··· α N Φ α ( ~r ) · · · Φ α N ( ~r N ) , and we define the σ ’th shell weight p ( σ ) ( r ) as before, i.e., p ( σ ) ( r ) ≡ X | β | = r | c ( σ ) β | . (3.17)We may then apply the analysis from Section III A toeach of the spin component functions, and note that the total shell-weight is p ( r ) ≡ X i ··· i N h ˆΦ i ··· i N , ψ i δ | β | ,r = X σ p ( σ ) ( r ) (3.18)since the shell number | β | = P | α k | does not depend onthe spin configuration of the basis function.Including the Pauli principle to accommodate properwave-function symmetry does not change these consider-ations. The basis functions ˆΦ i ··· i N are anti-symmetrizedto become Slater determinants Ψ i ··· i N (see for exampleRef. 27 for details), which is equivalent to consider theprojection H AS ( N ) = P AS H ( N ) of the unsymmetrizedspace onto the antisymmetric subspace. Moreover, theprojections P R and P AS commute, so that the shell-truncated space is given by P AS ,R = span ( Ψ i , ··· ,i N : i = i < · · · < i N , X k | α k | ≤ R ) , (3.19)which is precisely the computational basis used in manyCI calculations. (See however also the discussion in Sec-tion IV.) We stress, that, P AS ,R is independent of the ac-tual one-body HO eigenfunctions used. The shell-weightof ψ ∈ H AS ( X ) is now given by p ( r ) = X i ··· i N h Ψ i ··· i N , ψ i δ | β ( i ··· i N ) | ,r , (3.20)and k P R ψ k = R X r =0 p ( r ) . (3.21)As should be clear now, studying approximation ofHermite functions in arbitrary dimensions automaticallygives the corresponding many-body HO approximationproperties, since the many-body eigenfunctions can beseen as 2 N component functions, and since the shell-truncated Hilbert space transfers to a many-body settingin a natural way. C. Two electrons revisited We return to the exact solutions of the two-electronquantum dot considered in Sec. II C. Recall, that thewave functions were on the form ψ ( r, θ ) = e imθ f ( r ) , (3.22)where f ( r ) decayed exponentially fast as r → ∞ . Assumenow, that ψ ∈ H k ( R ), i.e., that all partial derivatives of ψ of order k exists in the weak sense, viz, ∂ jx ∂ k − jy ψ ∈ L ( R ) , ≤ j ≤ k, (3.23)where x = r cos( θ ) and y = r sin( θ ). Then, by Lemma 4in the Appendix, ( a † x ) j ( a † y ) k − j ψ ∈ L ( R ) for 0 ≤ j ≤ k as well.The function ψ ( r, θ ) was expanded in Fock-Darwin or-bitals, viz, ψ ( r, θ ) = ∞ X n =0 c n Φ FD n,m ( r, θ ) . (3.24)Recall, that the shell number N for Φ FD n,m was given by N = 2 n + | m | . Thus, the shell-weight p ( N ) is in this casesimply p ( N ) = | c ( N −| m | ) / | , N ≥ | m | , (3.25)and p ( N ) = 0 otherwise. From Prop. 2, we have ∞ X N = | m | N k p ( N ) < + ∞ , (3.26)which yields | c n | = o ( n − ( k +1+ ǫ ) / ) , ≤ ǫ < , (3.27)as claimed in Sec. II C. D. Smoothness properties of many-electron wavefunctions Let us mention some results, mainly due to Hoffmann-Ostenhof et al. , concerning smoothness of many-electron wave functions. Strictly speaking, their resultsare valid only in d = 3 spatial dimensions, since theCoulomb interaction in d = 2 dimensions fails to be aKato potential, the definition of which is quite subtle and out of the scope for this article . On the other hand, it isreasonable to assume that the results will still hold true,since the analytical results of the N = 2 case is very sim-ilar in the d = 2 and d = 3 cases: The eigenfunctionsdecay exponentially with the same cusp singularities atthe origin. Consider the Schr¨odinger equation ( H − E ) ψ ( ξ ) = 0,where ξ = ( ξ , · · · , ξ Nd ) = ( ~r , · · · , ~r N ) ∈ R Nd , andwhere ψ ( ξ ) is only assumed to be a solution locally. (Aproper solution is of course also a local solution.) Recall,that ψ has 2 N spin-components ψ ( σ ) . Define a coalescepoint ξ CP as a point where at least two particles coin-cide, i.e., ~r j = ~r ℓ , j = ℓ . Away from the set of suchpoints, ψ ( σ ) ( ξ ) is real analytic, since the interaction isreal analytic there. Near a ξ CP , the wave function hasthe form ψ ( σ ) ( ξ + ξ CP ) = r k P ( ξ/r )(1 + ar ) + O ( r k +1 ) , (3.28)where r = k ξ k , P is a hyper-spherical harmonic (on thesphere S Nd − ) of degree k = k ( ξ CP ), and where a is aconstant. It is immediately clear, that ψ ( σ ) ( ξ ) is k + 1times weakly differentiable in a neighborhood of ξ CP .However, at K -electron coalesce points, i.e., at points ξ CP where K different electrons coincide, the integer k may differ. Using exponential decay of a proper eigen-function, we have ψ ( σ ) ∈ H min( k )+1 ( R Nd ). Hoffmann-Ostenhof et al. also showed, that symmetry restrictionson the spin-components due to the Pauli principle in-duces an increasing degree k of the hyper-spherical har-monic P , generating even higher order of smoothness. Ageneral feature, is that the smoothness increases with thenumber of particles.However, their results in this direction are not generalenough to ascertain the minimum of the values for k fora given wave function, although we feel rather sure thatsuch an analysis is possible. Suffice it to say, that theresults are clearly visible in the numerical calculations inSec. V.Another interesting direction of research has been un-dertaken by Yserentant, who showed that there aresome very high order mixed partial derivatives at coa-lesce points. It seems unclear, though, if this can beexploited to improve the CI calculations further. IV. THE CONFIGURATION INTERACTIONMETHODA. Convergence analysis using HO eigenfunctionbasis The basic problem is to determine a few eigenvaluesand eigenfunctions of the Hamiltonian H in Eqn. (2.15),i.e., Hψ k = E k ψ k , k = 1 , · · · , k max . (4.1)The CI method consists of approximating eigenvalues of H with those obtained by projecting the problem onto afinite-dimensional subspace H h ⊂ H ( N ). As such, it isan example of the Ritz-Galerkin variational method. We comment, that the convergence of the Ritz-Galerkinmethod is not simply a consequence of the completenessof the basis functions. We will analyze the CI methodwhen the model space is given by H h = P R H AS ( N ) = P R ( N )= span ( Ψ i , ··· ,i N : X k | α k | ≤ R ) , used in Refs. 10,32, for example, although other spacesalso are common. (We drop the subscript “AS” from nowon.) The space M R ( N ) := span (cid:26) Ψ i , ··· ,i N : max k | α k | ≤ R (cid:27) , (4.2)i.e., a cut in the single-particle shell numbers (or energy)instead of the global shell number (or energy) is alsocommon. For obvious reasons, P R ( N ) is often referredto as an “energy cut space”, while M R ( N ) is referred toas a “direct product space”.As in Sec. III A, P R is the orthogonal projector ontothe model space P R ( N ). We also define Q R = 1 − P R as the projector onto the excluded space P R ( N ) ⊥ . Thediscrete eigenvalue problem is then( P R HP R ) ψ h,k = E h,k ψ h,k , k = 1 , · · · , k max . (4.3)The CI method becomes, in principle, exact as R → ∞ .Indeed, a widely-used name for the CI method is “exactdiagonalization,” being somewhat a misnomer as only avery limited number of degrees of freedom per particle isachievable.It is clear that P R ( N ) ⊂ M R ( N ) ⊂ P NR ( N ) , (4.4)so that studying the convergence in terms of P R ( N ) issufficient. In our numerical experiments we therefore fo-cus on the energy cut model space. A comparison be-tween the convergence of the two spaces is, on the otherhand, an interesting topic for future research.Using the results in Refs. 30,33 for non-degenerateeigenvalues for simplicity, we obtain an estimate for theerror in the numerical eigenvalue E h as E h − E ≤ [1 + ν ( R )](1 + Kλ ) h ψ, Q R T ψ i , (4.5)where K is a constant, and where ν ( R ) → R → ∞ .Using T Φ β = ( N d/ | β | )Φ β and Eqn. (3.8), we obtain h ψ, Q R T ψ i = ∞ X r = R +1 (cid:18) N d r (cid:19) p ( r ) . (4.6)Assume now, that ψ ( σ ) ∈ H k ( R Nd ) for all σ , so thataccording to Proposition 2, we will have ∞ X r =0 r k p ( r ) < + ∞ (4.7) implying that rp ( r ) = o ( r − k ). We then obtain, for k > h ψ, (1 − P R ) T ψ i = o ( R − ( k − ) + o ( R − k ) . (4.8)For k = 1 (which is the worst case), we merely obtainconvergence, h ψ, (1 − P R ) T ψ i → R → ∞ . We as-sume, that R is sufficiently large, so that the o ( R − k ) termcan be neglected.Again, we may observe a slight deviation from the de-cay, and we expect to observe eigenvalue errors on theform E h − E ∼ (1 + Kλ ) R − ( k − ǫ ) , (4.9)where 0 ≤ ǫ < k ψ h − ψ k (recall that ψ h = P R ψ ), we mention that k ψ h − ψ k ≤ [1 + η ( R )] [(1 + Kλ ) h ψ, (1 − P R ) T ψ i ] / , (4.10)where η ( R ) → R → ∞ . B. Effective interaction scheme Effective interactions have a long tradition in nuclearphysics, where the bare nuclear interaction is basicallyunknown and highly singular, and where it must berenormalized and fitted to experimental data. In quan-tum chemistry and atomic physics, the Coulomb inter-action is of course well-known so there is no intrinsicneed to formulate an effective interaction. However, inlieu of the in general low order of convergence impliedby Eqn. (4.9), we believe that HO-based calculations likethe CI method in general may benefit from the use ofeffective interactions.A complete account of the effective interaction schemeoutlined here is out of scope for the present article, but werefer to Refs. 2,11,34,35,36 for details as well as numericalalgorithms.Recall, that the interaction is given by U = N X i 2) := H eff − P R T P R , (4.15)which is defined only within P R ( N = 2).The N -body effective Hamiltonian is defined by H eff := P R T P R + N X i We now present numerical results using the fullconfiguration-interaction method for N = 2–5 electronsin d = 2 dimensions. We will use both the “bare” Hamil-tonian H = T + U and the effective Hamiltonian (4.16).Since the Hamiltonian commutes with angular momen-tum L z , the latter taking on eigenvalues M ∈ Z , theHamiltonian matrix is block diagonal. (Recall, that theFock-Darwin orbitals Φ FD n,m are eigenstates of L z witheigenvalue m , so each Slater determinant has eigenvalue M = P Nk =1 m k .) Moreover, the calculations are done in abasis of joint eigenfunctions for total electron spin S andits projection S z , as opposed to the Slater determinantbasis used for convergence analysis. Such basis functionsare simply linear combinations of Slater determinantswithin the same shell, and further reduce the dimension-ality of the Hamiltonian matrix. The eigenfunctions of H are thus labeled with the total spin S = 0 , , · · · , N for even N and S = , , · · · , N for odd N , as wellas the total angular momentum M = 0 , , · · · . ( − M produce the same eigenvalues as M , by symmetry.) Wethus split P R ( N ) (or M R ( N )) into invariant subspaces P R ( N, M, S ) ( M R ( N, M, S )) and perform computationssolely within these.The calculations were carried out with a code similarto that described by Rontani et al. in Ref. 8. TableI shows comparisons of the present code with that ofTable IV of Ref. 8 for various parameters using the modelspace M R ( N, M, S ). Table I also shows the case λ = 1, N = 2, M = 0, and S = 0, whose exact lowest eigenvalueis E = 3, cf. Sec. II C. We note that there are somediscrepancies between the results in the last digits of theresults of Ref. 8. The spaces M R ( N ) were identical inthe two approaches, i.e., the number of basis functionsand the number of non-zero matrix elements producedare cross-checked and identical.We have checked that the code also reproduces theresults of Refs. 9,10,39, using the P R ( N, M, S ) spaces.Our code is described in detail elsewhere where it isalso demonstrated that it reproduces the eigenvalues ofan analytically solvable N -particle system to machineprecision. B. Experiments For the remainder, we only use the energy cut spaces P R ( N, M, S ). Figure 2 shows the development of the low-est eigenvalue E = E ( N, M, S ) for N = 4, M = 0 , , S = 0 as function of the shell truncation parame-ter R , using both Hamiltonians H and H eff . Apparently,the effective interaction eigenvalues provide estimates forthe ground state eigenvalues that are better than thebare interaction eigenvalues. This effect is attenuated0 TABLE I: Comparison of current code and Ref. 8. Figures from the latter have varying number of significant digits. We includemore digits from our own computation for reference R = 5 R = 6 R = 7 N λ M S Current Ref. 8 Current Ref. 8 Current Ref. 82 1 0 0 3 . . . . . . . . . . . . . . . . . . . . . . 046 11 . . . . . 055 11 . . . . . 691 23 . . . . . 870 23 . . . . . 15 21 . . 13 21 . . 134 0 5 29 . . 44 29 . . 31 29 . . E FIG. 2: Eigenvalues for N = 4, S = 0, λ = 2 as functionof R for H (solid) and H eff (dashed). M = 0 , with higher N , due to the fact that the two-electron ef-fective Coulomb interaction does not take into accountthree- and many-body effects which become substantialfor higher N .We take the H eff -eigenvalues as “exact” and graph therelative error in E ( N, M, S ) as function of R on a loga-rithmic scale in Fig. 3, in anticipation of the relationln( E h − E ) ≈ C + α ln R, α = − ( k − ǫ ) . (5.1)The graphs show straight lines for large R , while for small R there is a transient region of non-straight lines. For N = 5, however, λ = 2 is too large a value to reachthe linear regime for the range of R available, so in thiscase we chose to plot the corresponding error for the verysmall value λ = 0 . 2, showing clear straight lines in theerror. The slopes are more or less independent of λ , asobserved in different calculations. In Fig. 4 we show the corresponding graphs when usingthe effective Hamiltonian H eff . We estimate the relativeerror as before, leading to artifacts for the largest valuesof R due to the fact that there is a finite error in the bestestimates for the eigenvalues. However, in all cases thereare clear, linear regions, in which we estimate the slope α .In all cases, the slope can be seen to decrease by at least∆ α ≈ − R shown, indicating thegain in accuracy when using small model spaces with theeffective interaction.Notice, that for symmetry reasons only even (odd) R for even (odd) M yields increases in basis sizedim[ P R ( N, M, S )], so only these values are included inthe plots.To overcome the limitations of the two-body effectiveinteraction for higher N , an effective three-body interac-tion could be considered, and is hotly debated in the nu-clear physics community. (In nuclear physics, there arealso more complicated three-body effective forces thatneed to be included. ) However, this will lead to a hugeincrease in memory consumption due to extra nonzeromatrix elements. At the moment, there are no methodsavailable that can generate the exact three-body effectiveinteraction with sufficient precision.We stress, that the relative error decreases very slowlyin general. It is a common misconception, that if a num-ber of digits of E ( N, M, S ) is unchanged between R and R + 2, then these digits have converged. This is not thecase, as is easily seen from Fig. 3. Take for instance N = 4, M = 0 and S = 0, and λ = 2. For R = 14and R = 16 we have E = 13 . E = 13 . . × − , while the correct relative error is 1 . × − .The slopes in Fig. 3 vary greatly, showing that theeigenfunctions indeed have varying global smoothness,as predicted in Sec. III D. For ( N, M, S ) = (5 , , / α ≈ − . 2, indicating that ψ ∈ H ( R ). Itseems, that higher S gives higher k , as a rule of thumb.1 −6 −5 −4 −3 −2 −1 R R e l a t i v e e rr o r Relative error for N=2, λ = 2 M=0, S=0, α = −1.0477M=0, S=2, α = −2.1244M=3, S=0, α = −3.1749M=3, S=2, α = −4.0188 6 8 10 12 14 16 18 2010 −5 −4 −3 −2 −1 R R e l a t i v e e rr o r Relative error for N=3, λ =2 M=0, S=1, α = −1.2772M=0, S=3, α = −2.1716M=2, S=1, α = −1.3093M=2, S=3, α = −2.5417 −4 −3 −2 −1 R R e l a t i v e e rr o r Relative error for N=4, λ = 2 M=0, S=0, α = −1.4233M=0, S=4, α = −2.8023M=3, S=2, α = −1.5327M=3, S=4, α = −3.2109 6 8 10 12 14 16 18 2010 −4 −3 −2 R R e l a t i v e e rr o r Relative error for N=5, λ = 0.2 M=0, S=1, α = −1.5150M=0, S=5, α = −3.6563M=3, S=1, α = −1.8159M=3, S=5, α = −4.2117 FIG. 3: Plot of relative error using the bare interaction for various N , M and S . Clear o ( R α ) dependence in all cases. Intuitively, this is because the Pauli principle forces thewave function to be zero at coalesce points, thereby gen-erating smoothness. VI. DISCUSSION AND CONCLUSION We have studied approximation properties of Her-mite functions and harmonic oscillator eigenfunctions.This in turn allowed for a detailed convergence analy-sis of numerical methods such as the CI method for theparabolic quantum dot. Our main conclusion, is, thatfor wave functions ψ ∈ H k ( R n ) falling off exponentiallyas k x k → ∞ , the shell-weight function p ( r ) decays as p ( r ) = o ( r − k − ). Applying this to the convergence the-ory of the Ritz-Galerkin method, we obtained the esti-mate (4.5) for the error in the eigenvalues. A completecharacterization of the upper bound on the differentia- bility k , i.e., in ψ ( s ) ∈ H k , as well as a study of theconstant K in Eqn. (4.9), would complete our knowledgeof the convergence of the CI calculations.We also demonstrated numerically, that the use ofa two-body effective interaction accelerates the conver-gence by at least an order of magnitude, which showsthat such a method should be used whenever possible.On the other hand, a rigorous mathematical study of themethod is yet to come. Moreover, we have not investi-gated to what extent the increase in convergence is in-dependent of the interaction strength λ . This, togetherwith a study of the accuracy of the eigenvectors , is anobvious candidate for further investigation.The theory and ideas presented in this article shouldin principle be universally applicable. In fact, Figs. 1–3of Ref. 11 clearly indicates this, where the eigenvaluesof He as function of model space size are graphed bothfor bare and effective interactions, showing some of the2 −10 −5 R R e l a t i v e e rr o r Relative error for N=2, λ = 2 M=0, S=0, U eff M=0, S=2, U eff M=3, S=0, U eff M=3, S=2, U eff −6 −5 −4 −3 −2 −1 R R e l a t i v e e rr o r Relative error for N=3, λ = 2 M=0, S=1, U eff , α = −4.1735M=0, S=3, U eff , α = −4.7914M=3, S=1, U eff , α = −4.4535M=3, S=3, U eff , α = −4.7521 −5 −4 −3 −2 −1 R R e l a t i v e e rr o r Relative error for N=4, λ = 2 M=0, S=0, U eff , α = −4.6093M=0, S=4, U eff , α = −5.3978M=3, S=0, U eff , α = −4.7816M=3, S=4, U eff , α = −6.2366 −5 −4 −3 −2 R R e l a t i v e e rr o r Relative error for N=5, λ = 0.2 M=0, S=1, U eff , α = −3.8123M=0, S=5, U eff , α = −4.4925M=3, S=1, U eff , α = −3.6192M=3, S=5, U eff , α = −5.0544 FIG. 4: Plot of relative error using an effective interaction for various N , M and S . Clear o ( R α ) dependence in all cases, butnotice artifacts when R is large, due to errors in most correct eigenvalues. The R = 5 case does not contain enough data tocompute the slopes with sufficient accuracy. features we have discussed.Other interesting future studies would be a direct com-parison of the direct product model space M R ( N ) andour energy cut model space P R ( N ). Both techniquesare common, but may have different numerical charac-teristics. Indeed, dim[ M R ( N )] grows much quicker thandim[ P R ( N )], while we are uncertain of whether the in-creased basis size yields a corresponding increased accu-racy.We have focused on the parabolic quantum dot, firstlybecause it requires relatively small matrices to be stored,due to conservation of angular momentum, but also be-cause it is a widely studied model. Our analysis is, how-ever, general, and applicable to other systems as well,e.g., quantum dots trapped in double-wells, finite wells,and so on. Indeed, by adding a one-body potential V to the Hamiltonian H = T + U we may model othergeometries, as well as adding external fields. Acknowledgments The author wishes to thank Prof. M. Hjorth-Jensen(CMA) for useful discussions, suggestions and feedback,and also Dr. G. M. Coclite (University of Bari, Italy)for mathematical discussions. This work was financed byCMA through the Norwegian Research Council.3 APPENDIX A: MATHEMATICAL DETAILS1. Multi-indices A very handy tool for compact and unified notationwhen the dimension n of the underlying measure space R n is a parameter, are multi-indices. The set I n of multi-indices are defined as n -tuples of non-negative indices,viz, α = ( α , · · · , α n ), where α k ≥ u be a formal vector of n symbols. Moreover,let φ ( ξ ) = φ ( ξ , ξ , · · · , ξ n ) be a function. Then, define | α | ≡ α + α + · · · + α n (A.1) α ! ≡ α ! α ! · · · α n ! (A.2) α ± β ≡ ( α ± β , · · · , α n ± β n ) (A.3) u α ≡ u α u α · · · u α n n , (A.4) ∂ α φ ( ξ ) ≡ ∂ α ∂ξ α ∂ α ∂ξ α · · · ∂ α d ∂ξ α d φ ( ξ ) (A.5)In Eqn. (A.3), the result may not be a multi-index whenwe subtract two indices, but this will not be an issue forus. Notice, that Eqn. (A.5) is a mixed partial derivativeof order | α | . Moreover, we say that α < β if and only if α j < β j for all j . We define α = β similarly. We alsodefine “basis indices” e j by ( e j ) j ′ = δ j,j ′ . We comment,that we will often use the notation ∂ x ≡ ∂∂x and ∂ k ≡ ∂∂ξ k to simplify notation. Thus, ∂ α = ( ∂ α , · · · , ∂ α n n ) , (A.6)consistent with Eqn. (A.4). 2. Weak derivatives and Sobolev spaces We present a quick summary of weak derivatives andrelated concepts needed. The material is elementary andsuperficial, but probably unfamiliar to many readers, sowe include it here. Many terms will be left undefined;if needed, the reader may consult standard texts, e.g.,Ref. 45.The space L ( R n ) is defined as L ( R n ) ≡ (cid:26) ψ : R n → C : Z R n | ψ ( ξ ) | d n ξ < + ∞ (cid:27) , (A.7)where the Lebesgue integral is more general than theRiemann “limit-of-small-boxes” integral. It is importantthat we identify two functions ψ and ψ differing only ata set Z ∈ R n of measure zero. Examples of such sets arepoints if n ≥ 1, curves if n ≥ 2, and so on, and countableunions of such. For example, the rationals constitute aset of measure zero in R . Under this assumption, L ( R n )becomes a Hilbert space with the inner product h ψ , ψ i ≡ Z R n ψ ( ξ ) ∗ ψ ( ξ )d n ξ, (A.8) where the asterisk denotes complex conjugation.The classical derivative is too limited a concept for theabstract theory of partial differential equations, includingthe Schr¨odinger equation. Let C ∞ be the set of infinitelydifferentiable functions which are non-zero only in a ballof finite radius. Of course, C ∞ ⊂ L ( R n ). Let ψ ∈ L ( R n ), and let α ∈ I n be a multi-index. If there existsa v ∈ L ( R n ) such that, for all φ ∈ C ∞ , Z R n ( ∂ α φ ( ξ )) ψ ( ξ )d n ξ = ( − | α | Z R n φ ( ξ ) v ( ξ )d n ξ, (A.9)then ∂ α ψ ≡ v ∈ L is said to be a weak derivative , ordistributional derivative, of ψ . In this way, the weakderivative is defined in an average sense, using integrationby parts.The weak derivative is unique (up to redefinition on aset of measure zero), obeys the product rule, chain rule,etc.It is easily seen, that if ψ has a classical derivative v ∈ L ( R n ) it coincides with the weak derivative. Moreover,if the classical derivative is defined almost everywhere(i.e., everywhere except for a set of measure zero), then ψ has a weak derivative.The Sobolev space H k ( R n ) is defined as the subset of L ( R n ) given by H k ( R n ) ≡ (cid:8) ψ ∈ L : ∂ α ψ ∈ L , ∀ α ∈ I n , | α | ≤ k (cid:9) . (A.10)The Sobolev space is also a Hilbert space with the innerproduct ( ψ , ψ ) ≡ X α, | α |≤ k h ∂ α ψ , ∂ α ψ i , (A.11)and this is the main reason why one obtains a unifiedtheory of partial differential equations using such spaces.The space H k ( R n ) for n > H k that are unboundedon arbitrary small regions but still differentiable. (Her-mite series for such functions would still converge fasterthan, e.g., for a function with a jump discontinuity!) Forour purposes, it is enough to realize that the Sobolevspaces offer exactly the notion of derivative we need inour analysis of the Hermite function expansions. 3. Proofs of propositions We will now prove the propositions given in Sec. III,and also discuss these results on mathematical terms.Recall that the Hermite functions φ n ∈ L ( R ), where n ∈ N is a non-negative integer, are defined by φ n ( x ) = (2 n n ! √ π ) − / H n ( x ) e − x / , (A.12)where H n ( x ) is the usual Hermite polynomial.4A well-known method for finding the eigenfunctions of H HO in one dimension involves writing H HO = a † a + 12 , (A.13)where the ladder operator a is given by a ≡ √ x + ∂ x ) , (A.14)with Hermitian adjoint a † given by a † ≡ √ x − ∂ x ) . (A.15)The name “ladder operator” comes from the importantformulas aφ n ( x ) = √ nφ n − ( x ) (A.16) a † φ n ( x ) = √ n + 1 φ n +1 ( x ) , (A.17)valid for all n . This can easily be proved by using therecurrence relation (2.7). By repeatedly acting on φ with a † we generate every Hermite function, viz, φ n ( x ) = n ! − / ( a † ) n φ ( x ) . (A.18)As Hermite functions constitute a complete, orthonor-mal sequence in L ( R ), any ψ ∈ L ( R ) can be written asa series in Hermite functions, viz, ψ ( x ) = ∞ X n =0 c n φ n ( x ) , (A.19)where the coefficients c n are uniquely determined by c n = h φ n , ψ i .An interesting fact is that the Hermite functions arealso eigenfunctions of the Fourier transform with eigen-values ( − i ) n , as can easily be proved by induction byobserving firstly that the Fourier transform of φ ( x ) is φ ( k ) itself, and secondly that the Fourier transform of a † is − ia † (acting on the variable y ). It follows fromcompleteness of the Hermite functions, that the Fouriertransform defines a unitary operator on L ( R ).We now make a simple observation, namely that( a † ) k φ n ( x ) = P k ( n ) / φ n + k ( x ) , (A.20)where P k ( n ) = ( n + k )! /n ! > k for n ≥ 0. Moreover P k ( n + 1) > P k ( n ) and P k +1 ( n ) >P k ( n ).We now prove the following lemma: Lemma 1 (Hermite series in one dimension) Let ψ ∈ L ( R ) . Then1. a † ψ ∈ L ( R ) if and only if aψ ∈ L ( R ) if and onlyif P ∞ n =0 n | c n | < + ∞ , where c n = h φ n , ψ i a † ψ ∈ L ( R ) if and only if xψ, ∂ x ψ ∈ L ( R ) ( a † ) k +1 ψ ∈ L ( R ) implies ( a † ) k ψ ∈ L ( R ) ( a † ) k ψ ∈ L ( R ) if and only if ∞ X n =0 n k | c n | < + ∞ . (A.21) ( a † ) k ψ ∈ L ( R ) if and only if x j ∂ k − jx ψ ∈ L ( R ) for ≤ j ≤ k Proof: We have k a † ψ k = ∞ X n =0 ( n + 1) | c n | = k ψ k + k aψ k , (A.22)from which statement 1 follows. Statement 2 follows fromthe definition of a † , and that a † ψ ∈ L implies aψ ∈ L (since k aψ k ≤ k a † ψ k ), which again implies xψ, ∂ x ψ ∈ L . Statement 3 follows from the monotone behaviourof P k ( n ) as function of k . Statement 4 then follows. Byiterating statement 2 and using [ ∂ x , x ] = 1 and statement3, statement 5 follows. ♦ The significance of the condition a † ψ ∈ L ( R ) is thatthe coefficients c n of ψ must decay faster than for acompletely arbitrary wave function in L ( R ). Moreover, a † ψ ∈ L ( R ) is the same as requiring ∂ x ψ ∈ L ( R ),and xψ ∈ L ( R ). Lemma 1 also generalizes this factfor ( a † ) k ψ ∈ L ( R ), giving successfully quicker decay ofthe coefficients. In all cases, the decay is expressed inan average sense, through a growing weight function in asum, as in Eqn. (A.21). Since the terms in the sum mustconverge to zero, this implies a pointwise faster decay, asstated in Eqn. (A.25) below.We comment here, that the partial derivatives must beunderstood in the weak, or distributional, sense: Eventhough ψ may not be everywhere differentiable in theordinary sense, it may have a weak derivative. For ex-ample, if the classical derivative exists everywhere exceptfor a countable set (and if it is in L ), this is the weakderivative. Moreover, if this derivative has a jump dis-continuity, there are no higher order weak derivatives.Loosely speaking, since xψ ( x ) ∈ L ⇔ ∂ y ˆ ψ ( y ) ∈ L ,where ˆ ψ ( y ) is the Fourier transform, point 2 of Lemma1 is a combined smoothness condition on ψ ( x ) and ˆ ψ ( y ).Point 5 is a generalization to higher derivatives, but isdifficult to check in general for an arbitrary ψ . On theother hand, it is well-known that the eigenfunctions ofmany Hamiltonians of interest, such as the quantum dotHamiltonian (2.15), decay exponentially fast as | x | → ∞ .For such exponentially decaying functions over R , x k ψ ∈ L ( R ) for all k ≥ 0, i.e., ˆ ψ ( y ) is infinitely differentiable.We then have the following lemma: Lemma 2 (Exponential decay in 1D) Assume that x k ψ ∈ L ( R ) for all k ≥ . Then a sufficientcriterion for ( a † ) m ψ ∈ L ( R ) is ∂ mx ψ ∈ L ( R ) , i.e., ψ ∈ H m ( R ) . In fact, x k ∂ m ′ x ψ ∈ L for all m ′ < m and all k ≥ . Proof: We prove the Proposition inductively. We notethat, since ∂ mx ψ ∈ L implies ∂ m − x ψ ∈ L , the proposi-tion holds for m − m . Moreover,it holds trivially for m = 1.Assume then, that it holds for a given m , i.e., that ψ ∈ H m implies x k ∂ m − jx ψ ∈ L for 1 ≤ j ≤ m and forall k (so that, in particular, ( a † ) m ψ ∈ L ). It remains toprove, that ψ ∈ H m +1 implies x k ∂ mx ψ ∈ L , since then( a † ) m +1 ψ ∈ L by statement 5 of Lemma 1. We computethe norm and use integration by parts, viz, k x k ∂ mx ψ k = Z R x k ∂ kx ψ ∗ ( x ) ∂ kx ψ ( x )= − k h ∂ mx ψ, x k − ∂ m − x ψ i−h ∂ m +1 x ψ, x k ∂ m − x ψ i < + ∞ . The boundary terms vanish. Therefore, x k ∂ mx ψ ∈ L forall k , and the proof is complete. ♦ The proposition states that for the subset of L ( R )consisting of exponentially decaying functions, the ap-proximation properties of the Hermite functions will only depend on the smoothness properties of ψ . Moreover, thederivatives up to the penultimate order decay exponen-tially as well. (The highest order derivative may decaymuch slower.)From Lemmas 1 and 2 we extract the following impor-tant characterization of the approximating properties ofHermite functions in d = 1 dimensions: Proposition 1 (Approximation in one dimension) Let k ≥ be a given integer. Let ψ ∈ L ( R ) be given by ψ ( x ) = ∞ X n =0 c n φ n ( x ) . (A.23) Then ψ ∈ H k ( R ) if and only if ∞ X n =0 n k | c n | < ∞ . (A.24) The latter implies that | c n | = o ( n − ( k +1) / ) . (A.25) Let ψ R = P R ψ = P Rn =0 c n φ n . Then k ψ − ψ R k = ∞ X n = R +1 | c n | ! / . (A.26)This is the central result for Hermite series approxima-tion in L ( R ). Observe that Eqn. (A.25) implies thatthe error k ψ − ψ R k can easily be estimated. See alsoProp. 2 and comments thereafter.Now a word on pointwise convergence of the Her-mite series. As the Hermite functions are uniformlybounded, viz, | φ n ( x ) | ≤ . ∀ x ∈ R , (A.27) the pointwise error in ψ R is bounded by | ψ ( x ) − ψ R ( x ) | ≤ . ∞ X n = R +1 | c n | . (A.28)Hence, if the sum on the right hand side is finite, the con-vergence is uniform. If the coefficients c n decay rapidlyenough, both errors can be estimated by the dominatingneglected coefficients.We now consider expansions of functions in L ( R n ).To stress that R n may be other than the configura-tion space of a single particle, we use the notation x =( x , · · · , x n ) ∈ R n instead of ~r ∈ R d . For N electrons in d spatial dimensions, n = N d .Recall that the Hermite functions over R n are indexedby multi-indices α = ( α , · · · , α n ) ∈ I n , and thatΦ α ( x ) ≡ φ α ( x ) · · · φ α n ( x n ) . (A.29)Now, to each spatial coordinate x k define the ladder op-erators a k ≡ ( x k + ∂ k ) / √ 2. These obey [ a j , a k ] = 0 and[ a j , a † k ] = δ j,k , as can easily be verified. Let a be a formalvector of the ladder operators, viz, a ≡ ( a , a , . . . , a n ) . (A.30)For the first Hermite function, we haveΦ ( x ) ≡ π − n/ e −k x k / . (A.31)By using Eqn. (A.18) and Eqn. (A.4), we may generateall other Hermite functions, viz,Φ α ( x ) ≡ α ! − / ( a † ) α Φ ( x ) . (A.32)Given two multi-indices α and β , we define the poly-nomial P α ( β ) by P β ( α ) ≡ ( α + β )! α ! = D Y j =1 P β j ( α j ) , (A.33)where P is defined for integers as before.Since the Hermite functions Φ α constitute a basis, any ψ ∈ L ( R n ) can be expanded as ψ ( x ) = X α c α Φ α ( x ) , k ψ k = X α | c α | , (A.34)where the sum is to be taken over all multi-indices α ∈I n . Now, let β ∈ I n be arbitrary. By using Eqn. (A.17)in each spatial direction we compute the action of ( a † ) β on ψ : ( a † ) β ψ = ( a † ) β · · · ( a † n ) β n X α c α Φ α = X α c α n Y k =1 ( α k + β k )! / α k ! / Φ α + β = X α c α P β ( α ) / Φ α + β . (A.35)6Similarly, by using Eqn. (A.16) we obtain a β ψ = a β · · · a β n n X α c α Φ α = X α c α + β n Y k =1 ( α k + β k )! / α k ! / Φ α = X α c α + β P β ( α ) / Φ α , (A.36)Computing the square norm gives k ( a † ) β ψ k = X α P β ( α ) | c α | . (A.37)and k a β ψ k = X α P β ( α ) | c α + β | . (A.38)The polynomial P β ( α ) > β, α ∈ I n , and P β ( α + α ′ ) > P β ( α ) for all non-zero multi-indices α ′ = 0. There-fore, if ( a † ) β ψ ∈ L ( R n ) then a β ψ ∈ L ( R n ). However,the converse is not true for n > c α , while Eqn. (A.37) is not. (This should be con-trasted with the one-dimensional case, where aψ ∈ L ( R )was equivalent to a † ψ ∈ L ( R ).) On the other hand, as inthe n = 1 case, the condition a † k ψ ∈ L ( R n ) is equivalentto the conditions x k ψ ∈ L ( R n ) and ∂ k ψ ∈ L ( R n ).We are in position to formulate a straightforward gen-eralization of Lemma 1. The proof is easy, so we omit it. Lemma 3 (General Hermite expansions) Let ψ ∈ L ( R n ) , with coefficients c α as in Eqn. (A.34),and let β ∈ I n be arbitrary. Assume ( a † ) β ψ ∈ L ( R n ) .Then ( a † ) β ′ ψ ∈ L ( R n ) and a β ′ ψ ∈ L ( R n ) for all β ′ ≤ β . Moreover, the following points are equivalent:1. ( a † ) β ψ ∈ L ( R n ) 2. For all multi-indices γ ≤ β , x γ ∂ β − γ ψ ∈ L ( R n ) .3. P α α β | c α | < + ∞ We observe, that as we obtained for n = 1, condition2 is a combined decay and smoothness condition on ψ ,and that this can be expressed as a decay-condition onthe coefficients of ψ in the Hermite basis by 3.Exponential decay of ψ ∈ L ( R n ) as k x k → ∞ im-plies that that x γ ψ ∈ L ( R n ) for all γ ∈ I n . We nowgeneralize Lemma 2 to the n -dimensional case. Lemma 4 (Exponentially decaying functions) Assume that ψ ∈ L ( R n ) is such that for all γ ∈ I n , x γ ψ ∈ L ( R n ) . Then, a sufficient criterion for ( a † ) β ψ ∈ L ( R n ) is ∂ β ψ ∈ L ( R n ) . Moreover, for all µ ≤ β , wehave x γ ∂ β − µ ψ ∈ L ( R n ) for all γ ∈ I n such that γ k = 0 whenever µ k = 0 , i.e., the partial derivatives of lowerorder than β decay exponentially in the directions wherethe differentiation order is lower. Proof: The proof is a straightforward application ofthe n = 1 case in an inductive proof, together with thefollowing elementary fact concerning weak derivatives: If1 ≤ j < k ≤ n , and if x j ψ ( x ) and ∂ k ψ ( x ) are in L ( R n ),then, by the product rule, ∂ k ( x j ψ ( x )) = x j ( ∂ k ψ ( x )) ∈ L ( R n ). Notice, that Lemma 2 trivially generalizes to asingle index in n dimensions, i.e., to β = β k e k , since theintegration by parts formula used is valid in R n as well.Similarly, the present Lemma is valid in n − n dimensions, as it must be valid for ¯ β =(0 , β , · · · , β n ).Assume that our statement holds for n − n dimensions. As-sume then, that ∂ β ψ ∈ L ( R n ). Let φ = ∂ ¯ β ψ ∈ L .Moreover, ∂ β φ ∈ L . Since ψ is exponentially decaying,and by the product rule, x γ φ ∈ L for all γ ≥ 0. ByLemma 2, x γ ∂ β − µ φ ∈ L for all γ and 0 < µ ≤ β .Thus, x γ ∂ β − e µ ψ ∈ L . Thus, the result holds as longas µ = e µ ; or equivalently µ = e k µ k for any k . Toapply induction, let χ = x γ ∂ β − µ ψ ∈ L . Note that ∂ ¯ β χ ∈ L and x ¯ γ χ ∈ L for all ¯ γ = (0 , γ , · · · , γ n ). Butby the induction hypothesis, x ¯ γ ∂ ¯ β − ¯ µ χ ∈ L for all ¯ µ ≤ ¯ β and all ¯ γ such that ¯ γ k = 0 if ¯ µ k = 0. This yields, usingthe product rule, that x γ ∂ β − µ ψ ∈ L for all µ ≤ β andall γ such that γ k = 0 if µ k = 0, which was the hypoth-esis for n dimensions, and the proof is complete. Notice,that we have proved that ( a † ) β ψ ∈ L as a by-product. ♦ In order to generate a simple and useful result forapproximation in n dimensions, we consider the casewhere ψ decays exponentially, and ψ ∈ H k ( R n ), i.e., ∂ β ψ ( x ) ∈ L ( R n ) for all β ∈ I n with | β | = k . In thiscase, we may also generalize Eqn. (A.25). For this, weconsider the shell-weight p ( r ) defined by p ( r ) ≡ X α ∈I n , | α | = r | c α | , (A.39)where c α = h Φ α , ψ i . Then, k ψ k = P ∞ r =0 p ( r ). More-over, if P projects onto the shell-truncated Hilbert space P R ( R n ), then k P ψ k = R X r =0 p ( r ) . (A.40) Proposition 2 (Approximation in n dimensions) Let ψ ∈ L ( R n ) be exponentially decaying and given by ψ ( x ) = X α c α Φ α ( x ) . (A.41) Then ψ ∈ H k ( R n ) , k ≥ , if and only if X α | α | k | c α | = ∞ X r =0 r k p ( r ) < + ∞ . (A.42) The latter implies that p ( r ) = o ( r − ( k +1) ) . (A.43)7 Moreover, for the shell-truncated Hilbert space P R , theapproximation error is given by k (1 − P ) ψ k = ∞ X r = R +1 p ( r ) ! / . (A.44) Proof: The only non-trivial part of the proof concernsEqn. (A.42). Since ψ is exponentially decaying and since ψ ∈ H k if and only if ∂ β ψ ∈ L for all β , | β | ≤ k , weknow that P α α β | c α | < + ∞ for all β , | β | ≤ k . Since | α | k is a polynomial of order k with terms of type a β α β , a β ≥ | β | = k , we have X α | α | k | c α | = X β, | β | = k a β X α α β | c α | < + ∞ . (A.45)On the other hand, since a β ≥ β has finitely many terms, P α | α | k | c α | < + ∞ implies P α α β | c α | < + ∞ for all β , | β | = k , and thus ψ ∈ H k since ψ was exponentially decaying. ♦ ∗ Electronic address: [email protected] S. M. Reimann and M. Manninen, Rev. Mod. Phys. ,1283 (2002). S. Kvaal, M. Hjorth-Jensen, and H. Møll Nilsen,Phys. Rev. B. , 085421 (2007). M. Hjorth-Jensen, T. Kuo, and E. Osnes, Phys. Rep. ,125 (1995). T. Helgaker, P. Jørgensen, and J. Olsen, MolecularElectronic-Structure Theory (Wiley, 2002). P. Maksym and T. Chakraborty, Phys. Rev. Lett. , 108(1990). S. M. Reimann, M. Koskinen, and M. Manninen, Phys.Rev. B , 8108 (2000). N. A. Bruce and P. A. Maksym, Phys. Rev. B , 4718(2000). M. Rontani, C. Cavazzoni, D. Bellucci, and G. Goldoni,J. Chem. Phys. , 124102 (2006). A. Wensauer, M. Korkusinski, and P. Hawrylak, SolidState Comm. , 115 (2004). S. A. Mikhailov, Phys. Rev. B , 153313 (2002). P. Navr´atil, J. Vary, and B. Barrett, Phys. Rev. C ,054311 (2000). R. Ram-Mohan, Finite Element and Boundary ElementApplications in Quantum Mechanics (Oxford UniversityPress, 2002). A. Kumar, S. Laux, and F. Stern, Phys. Rev. B , 5166(1990). M. Macucci, K. Hess, and G. J. Iafrate, Phys. Rev. B ,R4879 (1997). P. A. Maksym and N. A. Bruce, Physica E , 211 (1997). T. Ezaki, N. Mori, and C. Hamaguchi, Phys. Rev. B ,6428 (1997). P. Maksym, Physica B , 233 (1998). M. B. Tavernier, E. Anisimovas, F. M. Peeters, B. Szafran,J. Adamowski, and S. Bednarek, Phys. Rev. B , 205305(2003). M. Taut, Phys. Rev. A , 3561 (1993). M. Taut, J. Phys. A: Math. Gen. , 1045 (1994). D. Pfannkuche, V. Gudmundsson, and P. A. Maksym,Phys. Rev. B , 2244 (1993). M. Hoffmann-Ostenhof, T. Hoffmann-Ostenhof, andH. Stremnitzer, Phys. Rev. Lett. , 3857 (1992). T. Kato, Commun. Pure Appl. Math. , 151 (1957). R. J. Bartlett and M. Musia l, Rev. Mod. Phys. , 291(2007). J. Boyd, J. Comp. Phys. , 382 (1984). E. Hille, Duke Math. J. , 875 (1939). S. Raimes, Many-Electron Theory (North-Holland, 1972). M. Hoffmann-Ostenhof, T. Hoffmann-Ostenhof, andH. Stremnitzer, Commun. Math. Phys. , 185 (1994). H. Yserentant, Numer. Math. , 731 (2004). I. Babuska and J. Osborn, Handbook of Numerical Analy-sis: Finite Element Methods , vol. 2 (North-Holland, 1996). I. Babuska and J. Osborn, SIAM J. Numer. Anal. , 1249(1987). M. B. Tavernier, E. Anisimovas, and F. M. Peeters,Phys. Rev. B , 125305 (pages 9) (2006). I. Babuska and J. Osborn, Math. Comp. , 275 (1989). D. Klein, J. Chem. Phys. , 786 (1974). S. Kvaal, Phys. Rev. C , 044330 (pages 9) (2008). S. Kvaal (2008), arXiv:0810.2644v1[cond-mat.mes-hall] . J. S. Ball, SIAM J. Numer. Anal. , 2311 (2003). T. Schucan and H. Weidenm¨uller, Ann. Phys. , 483(1973). S. A. Mikhailov, Phys. Rev. B , 115312 (2002). N. F. Johnson and M. C. Payne, Phys. Rev. Lett. , 1157(1991). E. Caurier, G. Martinez-Pinedo, F. Nowacki, A. Poves, andA. P. Zuker, Rev. Mod. Phys. , 427 (pages 62) (2005). M. Helle, A. Harju, and R. M. Nieminen, Phys. Rev. B ,205329 (pages 24) (2005). M. Brask´en, M. Lindberg, D. Sundholm, and J. Olsen,Phys. Rev. B , 7652 (2000). S. Corni, M. Brask´en, M. Lindberg, J. Olsen, and D. Sund-holm, Phys. Rev. B , 045313 (2003). L. Evans, Partial Differential Equations , vol. 19 of