Natural occupation numbers: When do they vanish?

The non-vanishing of the natural orbital occupation numbers of the one-particle density matrix of many-body systems has important consequences for the existence of a density matrix-potential mapping for nonlocal potentials in reduced density matrix functional theory and for the validity of the extended Koopmans' Theorem. On the basis of Weyl's theorem we give a connection between the differentiability properties of the ground state wave function and the rate at which the natural occupations approach zero when ordered as a descending series. We show, in particular, that the presence of a Coulomb cusp in the wave function leads, in general, to a power law decay of the natural occupations, whereas infinitely differentiable wave-functions typically have natural occupations that decay exponentially. We analyze for a number of explicit examples of two-particle systems that in case the wave function is non-analytic at its spatial diagonal (for instance, due to the presence of a Coulomb cusp) the natural orbital occupations are non-vanishing. We further derive a more general criterium for the non-vanishing of NO occupations for two-particle wave functions with a certain separability structure. On the basis of this criterium we show that for a two-particle system of harmonically confined electrons with a Coulombic interaction (the so-called Hookium) the natural orbital occupations never vanish.


I. INTRODUCTION
The fractional occupation numbers n j of the correlated one-body reduced density matrix (1RDM) have intrigued many scientists in the past decades. They are defined by the eigenvalue equation dx γ(x, x )φ j (x ) = n j φ j (x) (1) where the 1RDM itself is defined in terms of the usual creation and annihilation field operators as γ(x, x ) := Ψ|ψ † (x )ψ(x)|Ψ for a state |Ψ where x := rσ is space-spin coordinate. The one-particle orbitals φ j (x) in Eq. (1) are denoted as the natural orbitals (NO) whereas the eigenvalues n j are called the NO occupation numbers. As an integral kernel the 1RDM is a bounded linear Hermitian operator with an infinite but countable eigenvalue spectrum and the set of all NOs form a basis in the set of quadratically integrable functions. If the state |Ψ is fermionic it is not difficult to prove that 0 ≤ n k ≤ 1 [1]. In the following we will restrict ourselves to electronic systems such that this property holds. The fact that the occupation numbers can also have non-integer values between zero and one is one of the most distinct features of interacting systems compared to non-interacting systems which can only have integer occupation numbers typically. Therefore the occupation numbers reflect strongly the electronic correlations present in the system under consideration. A system is considered weakly correlating when the occupation numbers differ only slightly from zero or one, in which case the full many-electron wavefunction can well be approximated by a single Slater determinant (non-interacting wavefunction). If a system is strongly correlated, the occupation numbers deviate strongly from integer values and multiple determinants are required to obtain a sufficiently accurate approximation to the many-body wavefunction which captures the physics of the system. The ability of the 1RDM occupation numbers to signal strong correlation has encouraged people to develop 1RDM functional theory as an alternative to traditional density function theory (DFT) to handle strongly correlated systems such as dissociating molecules [2][3][4], Mott insulators [5] and quantum Hall systems [6], for which the current approximate density functionals fail miserably. The sum of the occupation numbers equals the number of electrons in the system. Therefore, if we order the occupation numbers, n k , from the highest to the lowest one, their values need to decay to zero sufficiently fast for k → ∞, i.e. lim k→∞ n k = 0, or even become zero after some point k max . The question whether they actually do become zero or only approach zero for k → ∞ is not only an academic question, but is also of practical interest for methods that try to build an accurate approximation to the wavefunction by making an expansion in terms of Slater determinants, e.g. configuration interactions (CI). This question has recently been addressed for the dissociating hydrogen molecule [7]. One would expect that an optimal set of orbitals exists which leads to the fastest convergence of the expansion of the wavefunction in terms of Slater determinants [1]. One can prove that if all determinants are taken into account (full CI), that the highest occupied NOs are the orbitals which give the fastest convergence towards the exact one in the L 2 -norm [8]. The NOs become even more interesting if the occupation numbers become all zero for k sufficiently large, since this would imply that only a finite set of NOs would already be sufficient to expand the full many-electron wavefunction.
The question if zero occupation numbers exist in Coulomb systems is maybe even more important for 1RDM functional theory. Basic theorems in 1RDM functional theory [9] follow similar arguments as the famous Hohenberg-Kohn theorem [10] of density functional which establishes a one-to-one correspondence between densities, potentials and non-degenerate ground states. The main difference between 1RDM functional theory and density functional theory is that the natural conjugate variable to the 1RDM is a non-local external potential of the form rather than the local potential of density-functional theory. It therefore immediately follows that the energy contribution of the nonlocal external field to the total energy is given by With this expression the Hohenberg-Kohn proof can be followed exactly as in density-functional theory and Gilbert [9] in fact did this to establish that there is a oneto-one correspondence between non-degenerate ground states |Ψ and their corresponding 1RDM γ. This is already sufficient to establish 1RDM functional theory, since the ground state energy can be written as a functional of the 1RDM E[γ]. In density-functional theory one can further prove that two different (up to a gauge) potentials can not have the same non-degenerate ground state. The analogous proof fails in 1RDM theory since there can exist nonlocal potentialsV with the property thatV for a given ground state |Ψ of some HamiltonianĤ. Such a potential can therefore always be added to this Hamiltonian without affecting the ground state (it could be that |Ψ is now an excited state but by multiplyingV by a small enough number we can ensure that |Ψ is still the ground state). Let us now see how Eq. (3) can come about. Let us first define the annihilation operator which annihilates the NO φ s from any many-body quantum state. Suppose now thatâ s |Ψ = 0 for some of the labels s, which means that the orbital φ s does not appear in any Slater determinant of a CI expansion of |Ψ . This implies that n s = Ψ|â † sâs |Ψ = 0 such that the corresponding NO occupation number vanishes. We can then construct the following one-body po-tentialV where v rs = v * sr is an arbitrary Hermitian matrix and where we sum only over the labels for which n r = n s = 0 for the state |Ψ . It is clear that this potential exactly has the propertyV |Ψ = 0. In real space this corresponds to a nonlocal spatial potential of the form of Eq. (2) where We therefore see that Eq. (3) can be satisfied whenever the state |Ψ has vanishing NO occupations. The nonvanishing of the NO occupation numbers for electronic ground states is therefore a necessary condition for the existence of a one-to-one mapping between nonlocal potentials and 1RDMs. To the best of our knowledge the answer to the question whether the necessary condition is also a sufficient one is unknown. The one-to-one mapping between non-local potentials and 1RDMs would be relevant for the foundations of linear response 1RDM functional theory and also its time-dependent extension would greatly benefit from the resulting simplifications. Other consequences of vanishing occupation numbers arise in the extended Koopmans' theorem [11][12][13][14]. The extended Koopmans' theorem is an extension to arbitrary wavefunctions of the well known theorem by Koopmans that the occupied Hartree-Fock orbital energies provide approximations to the ionization energies [15]. If the exact wavefunction is used in the extended Koopmans' procedure, even the exact ionization energies should result, provided the set of partially occupied NOs is complete, i.e. none of the occupation numbers vanishes. A less restrictive condition has been derived by Pernal and Cioslowski [16], though in practice it simply implies that none of occupation numbers should vanish. For systems with Coulombic interactions the extended Koopmans' theorem is found to hold to very high numerical accuracy [17,18] although this does not prove its validity.
We have therefore seen that the possible vanishing of NO occupation numbers has important consequences for CI expansions as well as for the validity of fundamental theorems in many-body theory. This then immediately raises the question in which cases the NO occupation numbers vanish. If none of the occupation numbers vanishes, then every NO is needed in an expansion of the ground state wave function. One general observation that one can make is that infinite expansions are typically required when expanding non-smooth functions in terms of smooth ones. In the case of electronic ground states the Coulomb interaction requires the wavefunction to have a cusp at the positions where the electrons come together of the form Ψ(r 12 → 0) = Ψ(r 12 = 0) 1 + 1 2 where r 12 := |r 1 − r 2 |. This cusp gives an infinite kinetic energy which exactly compensates the infinity from the Coulomb interaction between the electrons [19][20][21][22]. Due to this non-analytic behavior of the wavefunction, a full expansion of the wavefunction in one-electron functions requires in general all functions to be present. Hence, one may expect in the particular case of an expansion in NOs, none of the NOs should have an occupation number equal to zero, since that would imply that the NO is not required in the expansion. Although this argument sounds very reasonable, it is certainly not a proof that zero occupations do not occur in Coulomb systems. Though infinitely many occupation numbers are required to be non-zero, it might be that some of them are still zero in some special situations. In the case of the homogeneous electron gas (HEG), however, this argument can be turned into a proof. Since the NOs of the HEG are simply plane waves, the occupation numbers are then given by the momentum distribution, n(k). Kimball has shown that the momentum distribution is required to decay as 1/k 8 due to the interelectronic cusp condition [23], so the occupation numbers never become exactly zero.
In the case of the HEG we were in the fortunate situation that the NOs are plane waves and, so that their occupation numbers are simply given by the momentum distribution. For general systems we are not in such a convenient position, because a straightforward expansion in a finite basis set effectively smoothens the electronelectron cusp (4) and the argument does not apply anymore. For two-electron systems we are in a more fortunate situation, however, since for singlet two-electron systems there is a strong connection between the NOs and the wavefunction. The spatial part of the singlet two-electron wavefunction is symmetric and can therefore be diagonalized By calculating the corresponding spin-integrated 1RDM, one readily finds that the eigenfunctions are NOs and that the coefficients are related to the occupation numbers as n k = c 2 k . Though we are not in such a good position as the for the HEG, this connection is quite useful, since it allows us to connect the behavior of the occupation numbers directly to the analytic properties of the wavefunction, instead of going via the 1RDM in which much of the analytic properties are integrated out. Therefore, we will focus our attention in this paper mainly to singlet two-electron systems to demonstrate how the form of the interaction determines the analytic properties of the wavefunction, which in turn dictates the asymptotic decay of the occupation numbers.

II. EXPLICIT EXAMPLES
Before we present a general treatment for simple explicitly correlated wavefunctions, we first consider some specific examples for which we can solve the NOs and coefficients explicitly or at least prove that none of the NOs have a vanishing coefficient (occupation number).

A. A simple 1D Hylleraas wavefunction
Let us first consider the simplest wavefunction with a cusp in one dimension (1D) where K is a normalization constant and α(x) is an arbitrary orbital apart from the fact that it is positive, α(x) > 0. To calculate the NOs and the coefficients in the spectral expansion of the wavefunction (5), we need to solve the following eigenvalue equation Introducing the following function ϕ k (x) := φ k (x)/α(x) the eigenvalue equation can be written as Now differentiating this equation twice with respect to x 1 , we obtain the following differential equation for ϕ k (x) where we used that |x| = 2δ(x). From this equation we see that c k = 0 if η = 0, since otherwise we would have ϕ k = 0 which is not an eigenfunction. Hence, for η = 0 we can divide by c k and write the equation as In the case that a simple Slater function is used for the orbital, α(x) = e −Z|x| , the differential equation can be cast into a Bessel's differential equation. The full construction of all the NOs and their coefficients is rather technical and has been deferred to Appendix B. The result of the calculation is that none of the occupation numbers is zero and that the occupation numbers behave asymptotically as where C is a constant. We will see later that such a power law behavior is typical for wave functions that are at most a finite number of times differentiable (in our case zero times).

B. A simple 3D Hylleraas wavefunction
Let us now try to extent the approach of the previous example to a three dimensional case, so that we have a wavefunction in three dimensions (3D) that satisfies the cusp condition (4) of the form where r 12 = |r 1 −r 2 |. Since the NOs are eigenfunctions of the wavefunction (5), they satisfy the eigenvalue equation Now introducing similar function as in the 1D case, ϕ k (r) := φ k (r)/α(r), the eigenvalue equation can be written as Unfortunately, just taking the Laplacian does not work, since ∇ 2 r |r − r | = 2/|r − r |. Taking the Laplacian a second time, however, we obtain the sought after delta function and our equation becomes We can now use the same argument as in the 1D case. Since the orbital α(r) will in general not vanish on some open set, c k = 0 would imply that φ k (r) = 0 for unoccupied NOs. Since such orbitals are not normalizable, our simple explicitly correlated wavefunction does not have any NOs with occupation numbers equal to zero. Since the differential equation for the functions ϕ k (r) is now fourth order and additionally in 3D, it becomes quite hard to obtain the NOs explicitly from this equation. In any case, it is exactly the cusp behavior that allows us to conclude that none of the NO occupation numbers vanish.

C. Double Harmonium
Although we have found that the exact treatment of the cusp in the simple Hylleraas wavefunction prevented the NOs to have a zero occupation number, one can wonder if the electron-electron cusp is actually essential to have only non-zero occupation numbers. This is actually not the case as one can show explicitly for a system with harmonic interactions in one dimension. The Hamiltonian for such a system can be written aŝ where λ ≥ 0. The coordinates can be decoupled by making a transformation to center-of-mass coordinates, which gives the following expression for the Hamiltonian where s : Since this is the Hamiltonian for two independent harmonic oscillators, we can immediately write spatial part of the singlet ground state wavefunction as The structure of this solution is sufficiently simple to determine the NOs and the corresponding wavefunction coefficients explicitly. Due to the purely harmonic nature of our system, one might suspect that the NOs also have the form of harmonic oscillator solutions. Indeed, one can actually calculate the NOs to be (see Appendix C) where H k (x) are the Hermite polynomials and a k normalization constants satisfying The corresponding wavefunction coefficients are We find that the wavefunction coefficients are alternating and only vanish in the limit k → ∞. Only whenω = ω (no interactions), all the coefficients become zero, except for the first one c 0 = 1. An important conclusion is now that instead of a power law behavior we now have an exponential decay of the form n k = C a k with C and a constants. We will show below that such a behavior is typical in the limit k → ∞ for infinitely differentiable wave functions. Our examples seem to indicate that there is a connection between the differentiability properties of the wave functions and the asymptotic behavior of the NO occupations. In the next Section we make this connection more precise.

III. A LOWER BOUND ON THE DECAY RATE OF THE OCCUPATION NUMBERS
The eigenvalue equation for the NOs (7) is a Fredholm integral equation, where the two-body wavefunction Ψ is the kernel and c k are the eigenvalues. (Often one considers the characteristic values, 1/c k , instead of the eigenvalues.) In the theory of Fredholm integral equations lower bounds on the decay rate of the eigenvalues have been established on the differentiability and analyticity of the integral kernel [24][25][26]. Although only lower bounds have been found, these bounds are a nice illustration of the strong link between the differentiability of the kernel and the decay rate of the kernel.
The first result of interest is by Hille and Tamarkin, who showed that the eigenvalues of operators with analytic kernels decay exponentially. More precisely where λ k are the eigenvalues and the constant R is related to the size of the region where the kernel is analytic [25,26]. Indeed, the kernel of the double harmomium is analytic and the corresponding coefficients decay exponentially (10).
A result for finitely differentiable integral kernels by Weyl [24] is of particular interest for wavefunctions with a cusp. Weyl showed that the eigenvalues of a finitely differentiable kernel only need to decay polynomially. More precisely, if the partial derivatives of a symmetric kernel are continuous up to order p, the decay rate of its eigenvalues is bounded by where d is the dimension of the integration variable. A derivation of Weyl's result can be found in Appendix A. Let us compare Weyl's theorem with the results for the simple 1D Hylleraas atom (Sec. II A). The first derivative of the wavefunction is already discontinuous in this case, so p = 0. This gives us a rather modest lower bound on the decay rate of the coefficients, comparing to the actual quadratic decay-rate k −2 of the coefficients (see Appendix B). The bound by Weyl can be tightened by using more regularity properties of the discontinuous derivatives of the wavefunction [25,26]. Since the discontinuity is bound, the two-electron wavefunction can be put in Lip 1 (1, 2) [25,26] which gives the lower bound k −3/2 on the decay-rate of the coefficients.
Though Weyl's theorem does not put a very stringent constraint on the decay of the occupation numbers, it clearly demonstrates that wavefunctions with an interelectronic cusp will typically have only a polynomial decaying occupation number spectrum, whereas the occupation numbers of a wavefunction without a cusp will decay exponentially. Since in a finite basis set the nonanalyticity of the cusp can not be fully represented and thereby effectively smoothened, the calculated occupation numbers will always decay exponentially or faster, i.e. too fast compared to the typical polynomial decay. This is actually not surprising, since if only a finite number of orbitals is included in the calculations, all the orbitals in the complement are automatically NOs with zero occupation number. Additionally, since a finite basis set representation effectively removes the cusp from the wavefunction (4), the original argument that NOs with vanishing occupation numbers do not exist in many-body Coulomb systems, does not even apply anymore.

IV. A PROOF BY FOURIER TRANSFORM
The example with the harmonic interaction shows that the cusp actually might not be essential to prevent occupation numbers becoming zero. It seems that any correlation that requires some r 12 behavior that can not be expressed in a finite number of simple orbital products will probably imply the absence of unoccupied NOs. We can make this idea more precise in the case of a singlet two-electron system which is limited to the form of a simple orbital product times an arbitrary correlation function depending only on r 1 − r 2 , i.e. a wavefunction of the form where α(r) > 0. Considering the situation that an NO, φ i (r), has a zero occupation number (so also c i = 0), the eigenvalue equation (7) then simplifies to where χ i (r) := α(r)φ i (r). Since this condition has the form of a convolution product, we can deconvolute it by taking the Fourier transform Provided thatf (k) = 0 almost everywhere, it follows thatχ i (k) = 0. Since α(r) = 0 almost everywhere, this implies that the NO φ i (r) = 0. Because this is not a normalizable function, we can conclude that no eigenvalues c i = 0 exist, iff (k) = 0 almost everywhere. The converse also holds. Iff (k) = 0 on some finite interval, an NO with zero occupancy exists by construction. To construct this NO, define aχ k (k) which is nonzero on this interval wheref (k) vanishes. By Fourier transforming back, we can construct the corresponding NO (φ k (r) = χ k (r)/α(r)). Since this will be a function in the null-space of the (linear) 1RDM-operator, it can not be expressed linear combination of NOs with finite occupancy. Let us show this more explicitly. Suppose that such a linear combination exists, then we write the new NO as Since the 1RDM-operator is linear, we can write its action on φ k as If we now take the inner product of this result with φ k , we find Hence, we have a contradiction, so the constructed NO with zero occupation number can not be expressed as a linear combination of the NOs with a finite occupancy. Now let us check if our proof indeed recovers the result that our previous examples do not have vanishing occupation numbers. In the case of the harmonium the correlation function is simply a Gaussian. Since the Fourier transform of a Gaussian is simply again a Gaussian, the Fourier transform is non-zero every where. Hence, no unoccupied NOs should exist, which is in agreement with our explicit construction. The situation is more complicate in the case of the simple Hylleraas wavefunction. The correlation function is not an L 2 function anymore, so we can expect distributions to appear in its Fourier transform [27]. Note that a divergent correlation function is allowed, provided the divergence of the correlation function is compensated by a stronger decay of the orbital α(r), to make the wavefunction normalizable. The Fourier transform in 1D becomes and in 3D we find The calculation of these Fourier transforms has been worked out in more detail in Appendix D. Since we find in both cases that the Fourier transform is non-zero everywhere, we recover the result that no unoccupied NOs exists as we have found before. Now let us apply the theorem to some other systems.

A. Inverse harmonic interaction
First we will consider a system with inverse harmonic interactions which has been considered before by Morrison et al. [28]. The inverse harmonic interaction is also sometimes referred to as the Calogero interaction [29][30][31]. The full Hamiltonian we will consider here is given aŝ The eigenstates of this Hamiltonian can be solved exactly by making a transformation to the centre-of-mass coordinates, which decouples the coordinates. The Hamiltonian for the centre-of-mass coordinate is simply the Hamiltonian of a harmonic oscillator, so is readily solved. The Hamiltonian for the relative coordinate is more involved, but can still be solved in terms of confluent hypergeometric functions [28,32]. Fortunately, the ground state reduces to the following particularly simple form where Γ(z) is the gamma function and α = √ 1 + 4λ − 1 /2 [33]. Since the ground state is simply the product of an orbital, α(r 1 )α(r 2 ), times a correlation function, f (r 12 ) = r α 12 , our theorem can be applied to this case. For the Fourier transform of the correlation function we find (see Appendix D for details) where δ (n) (k) denotes the n th order derivative of the delta-function. Therefore, we find that zero occupation numbers can only exist for even α, which is in agreement with the findings of Morisson et al. [28]. Actually, it is not surprising that there are only a finite number of unoccupied NOs for even α, since in that case the correlation function becomes exactly separable in r 1 and r 2 , so the wavefunction can be represented by a finite number of orbitals. For example in the simplest non-trivial case of α = 2 the wavefunction can be written as where we used χ f (r) := f (r)e − 1 2 ωr 2 as a compact notation for the various one-particle functions. Since only five orbitals are required to represent this wavefunction, only five NOs with non-zero occupation number exist. We see that the contribution from the p orbitals to the wavefunction is already diagonal, so to obtain the NOs, we only need to normalize them φ x (r) = 2ω ω π 3/2 x e − 1 2 ωr 2 and we have similar expressions for φ y (r) and φ z (r) of course [34]. Their coefficient in the spectral expansion (5) is readily obtained as c x = c y = c z = −1/ √ 15. The contribution from the other two orbitals χ 1 (r) and χ r 2 (r) is not diagonal, so has to be diagonalized. This is readily achieved by the following linear combinations where a := √ 15/4. The orbitals φ ± (r) are constructed such that they are orthonormal, so they are the required NOs. Their corresponding expansion coefficients can be calculated to be c ± = (4a ± 5)/10. This example clearly demonstrates that the cusp is actually not essential for the absence of unoccupied NOs. The absence of zero occupation numbers is caused by correlation in the full many-body wavefunction, which can not be expanded in a finite series of one-electron functions. The cusp actually causes the exact wavefunction to have such a form, hence there will be an infinite amount of non-zero occupation numbers.

B. Hookium atom
In this section we will apply our theorem to Hooke's atom. The interaction is now the Coulomb interaction, the interaction of interest for electrons. However, the confining potential is still harmonic to allow for the separation of variables by changing the coordinates to the centre-of-mass frame. The Hamiltonian is given aŝ The ground state has the following form where N is a normalization constant. The function t(ρ) satisfies the following differential equation [35,36] for the ground state ρt − 2ρ 2 t + (˜ r − 1)ρ − λ ω/2 t = 0, where˜ r = 2 r /ω with r as the contribution to the energy from the relative coordinate (the total energy is 3 2 ω + r ). Unfortunately, this differential equation does not allow for an explicit solution, though a series solution can be constructed which even has only a finite number of terms for specific ratios λ 2 /ω. More importantly, by neglecting the second order derivative and the last term, we readily find that the solution has to behave asymptotically for large ρ as 2st (s) + s 2 + r + 3 t (s) (17) where we used that t(0) = 0. The Fourier transform of the correlation function can now directly be obtained fromt(s) as where the limits are important to include possible poles at the origin. Since the Laplace transform is a solution of a second order differential equation, it has to be an analytic function. This analyticity is carried over to the Fourier transform except for some possible irregularities located at the origin (k = 0). Since an analytic function can only be zero in an open set if it vanishes everywhere, the Fourier transform of the correlation function of the Hookium atom does not vanish in a finite region and hence, the NO coefficients (occupation numbers) do not vanish for the Hookium atom. This result is orthogonal to the claim made by Cioslowski and Pernal [37]. They have studied the behavior of the wavefunction coefficients of Hooke's atom numerically and found that the coefficients become very small and that the sign pattern of the most significant coefficients around these points change. Therefore, they concluded that the expansion coefficients have to become zero to change their sign. As a courtesy to the reader, we have repeated their calculations to obtain an accurate expression for the wavefunction and calculated the NO coefficients by diagonalizing the wavefunction directly. The results for the occupation numbers in the s-channel are shown in Fig. 1 and are identical to the ones reported in Ref. [37]. Indeed, the most significant NOs have a different sign in the large and small ω limit. However, upon closer inspection of the plot in Fig. 1, one readily sees that coefficients of different NOs actually gain in amplitude when making the transition between the weakly and strongly correlated regime, so there is actually no evidence that the NO coefficients do cross zero. The numerical results of Cioslowski and Pernal are therefore in agreement with our findings.

V. CONCLUSION
The question of the existence of NOs with vanishing occupation numbers in many-body Coulomb systems is important for a number of practical applications. For the direct expansion of wavefunctions in finite orbital basis sets (CI expansions) it would be beneficial if only a finite number of NOs has a finite occupancy. However, the presence of unoccupied NOs will cause complications in the formal developments of 1RDM functional theory. For the extended Koopmans' theorem the existence of vanishing natural occupation numbers could even be catastrophic, since the ionization energies do not necessarily converge to the exact ones when the approximate wavefunction converges to the exact many-body state.
The divergence of the Coulomb interaction between the electrons requires the wavefunction to have a cusp at the coalescence points of the electrons (4). This non-analytic behavior can only be represented by including an infinite amount of orbitals (NOs) in the expansion of the wavefunction, so there are an infinite amount of NOs with a non-zero occupation number. However, this argument does not provide a proof that natural occupation numbers equal to zero do not exist, i.e. that the NOs with a non-zero occupation numbers form a complete set. However, we have been able to show that wavefunctions of the form α(r 1 )α(r 2 )f (r 1 − r 2 ) do not have any vanishing occupation number, if and only if the Fourier transform of the correlation function, f (r 12 ), does not vanish on an open set. The Fourier transform of the correlation function only seems to disappear if the wavefunction is separable, i.e. representable in a finite basis. Applying our theorem to the harmonium atom, we have shown that the occupation numbers do not vanish in contradiction with earlier assertions by Cioslowski and Pernal based on numerical calculations [37]. However, a more careful inspection of their results showed that only their interpretation was incorrect and that their results actually agree with our proof up to numerical accuracy.
Further, we have demonstrated that a discontinuity in the wavefunction is not required for the absence of unoccupied NOs. Even the perfectly smooth ground state of the double harmonium has no vanishing occupation numbers. More essential is the non-separability of the wavefunction, which is caused by the discontinuity of the cusp. The discontinuity of cusp does have an effect on the decay-rate of the occupation numbers. It causes the occupation numbers to decay merely polynomially compared to an exponential decay of the occupation numbers for wavefunctions without a cusp.
Though the proof of (11) is "äußerst einfach" according to Weyl [24], it might be worthwhile to expose its derivation. Weyl starts with an alternative derivation of a theorem by E. Schmidt [38] to prove that where λ k are the eigenvalues of the integral kernel K(s, t) ordered in descending order and k n (s, t) is a hermitian finite rank operator k n (s, t) = n i,j=1 with g i ∈ L 2 and k ij = k * ji . This theorem is easily understood by using the spectral representation of the integral kernel. The minimum value of the integral is achieved by using the largest eigenvalues of K(s, t) at the diagonal, k ij = λ i δ ij and the corresponding eigenfunctions for the functions g p . This choice exactly eliminates the largest eigenvalues of K(s, t) and only the smaller n + 1 eigenvalues will contribute to the integral. Any other choice for k n (s, t) will give a larger value of the integral.
For more rigor, consider the following proof. To cover the infinite dimensional case we use the Rayleigh quotient to define the eigenvalues. The first eigenvalue (and largest in magnitude) is defined as where · is the usual L 2 norm and the operatorK is defined be the action of the integral kernel K(s, t) on a function f asK The function that achieves this maximum, φ 1 (s), is the corresponding eigenfunction. The other eigenvalues are defined (found) by searching over a subspace where the previously found eigenfunctions have been projected out λ n+1 := min φ1,...,φn max f ⊥φ1,...,φn f =1 K f and the function that achieves this maximum, φ n+1 , is the corresponding eigenfunction. Using this definition we readily find for the eigenvalue of the the sum of two linear operator K 1 and K 2 [26] Since the rank of k n is only n, it has only n non-zero eigenvalues at maximum. Therefore, we find that Applying this inequality for all eigenvalues of K, we readily recover Schmidt's inequality (A1). Schmidt's inequality gives a lower bound for the integral. An upper bound can be obtained from Taylor's theorem. This procedure is probably most clearly explained at the end of Ref. [39]. To simplify the analysis, we assume without loss of generality that we integrate over a finite block with sides of length L, so that we can divide it in m d smaller blocks, where d is the dimension of our integration variable. In each of these regions we can make a Taylor expansion around its centre s 0 where p denotes a multi-index 1 · · · ∂s p d d and the remainder satisfies lim s→s0 h p (s) = 0.
By choosing all possible powers of s up to order p as basis functions for k n , we can create N = p+d d linearly independent functions per block, so N m d functions g i in total. Using these basis functions, we can set the kernel k n with n = N m d equal to the Taylor expansions of K in these blocks. The error of the Taylor expansion can now be approximated as where m → 0 as m → ∞. Now combining this inequality from the Taylor expansion of the kernel with Schmidt's inequality (A1) we have Now setting k = 2N m d and cleaning up the limit, Weyl's inequality for the asymptotic behavior of the eigenvalues (11) readily follows. Tighter bounds on the asymptotic decay of the eigenvalues can be found by using additional properties of the integral kernel [25,26].
Appendix B: NOs for the 1D model atom First we need to determine the boundary condition to be imposed on the solutions, which will give the required quantization of the expansion coefficients c k . To find the boundary conditions, we first rewrite the integral equation as Now considering the limit x → ∞, the last integral vanishes and we find that the solutions have to behave asymptotically as since the orbital should decay exponentially for the wavefunction to be normalizable. Because we are dealing with a symmetric orbital, the solutions can be separated in gerade and ungerade functions. For the even solutions only the first contribution survives, so we have the following boundary condition for x → ∞ for the gerade solutions Likewise, for the ungerade solutions only the last term survives and we find Note that this analysis of the boundary conditions is completely general for symmetric α 2 (x). Now we have found the boundary conditions for the various solutions, we turn our attention to the solution of the differential equation where the normalization constant is calculated to be Since we are looking for symmetry adapted solutions, we only need to solve this differential equation for x ≥ 0. If we define (for x ≥ 0) then for a function ϕ(x) = f s(x) we have and the differential equation in term of f (s) becomes This is the Bessel differential equation for the zeroth order Bessel functions and the general solution to this equation is with C 1 and C 2 constants and J 0 and Y 0 Bessel functions of the first and second kind respectively and I 0 and K 0 their modified counterparts. Now going back to the original function u(x) we find for λ < 0 and in the case of λ > 0 we find Now we construct the even and odd solutions by imposing the corresponding boundary conditions. First we impose the boundary conditions at x = 0. The odd solutions need to vanish at the origin, so we find whereλ := |λ|/Z and C is a normalization constant. For the even solutions the first order derivative at x = 0 needs to vanish, so for the even solutions we find . To obtain the proper quantization of the eigenvalue λ, we need to impose the proper boundary conditions for x → ∞, i.e. s → 0. From the asymptotic behavior of the Bessel functions for small s, we find that our ungerade solutions behave asymptotically for x → ∞ as where γ is the Euler-Mascheroni constant. We found before, however, that the odd solutions do not have a linear term (B2), so we must have Since I 0 (y) does not have any zero, only solutions for λ < 0 exist, which are related to the zeros of the the zeroth order Bessel function of the first kind, J 0 (y k ) = 0, as Using (6) the coefficients of the odd NOs are readily determined to be The zero's of the zeroth order Bessel function behave asymptotically as so asymptotically, the coefficients decay quadratically, c u,k = O k−2 . To construct the corresponding NOs, we use that u u,k (−x) = −u u,k (x), so The normalization constant of the ungerade NOs is readily calculated by using that the n th order Bessel functions satisfy where α nj is the j th zero of J n . Making the substitution ρ = e −Zx and taking n = 0, we find that the normalization constant for the odd NOs is given as . Now we will construct the even solutions. First we consider the large x behavior of the even solutions. From the asymptotic behavior of the Bessel functions for small s we find that for x → ∞ Since we know that the the even solutions have to behave for x → ∞ as given by the asymptotic relation in (B1), we must have that the the ratio between the linear and constant term must be equal to η, i.e. f ± λ = 0, where Note that in f + we could divide by I 1 (y), since this is an exponentially growing function and has no zeros apart from y = 0.
Let us now first consider the positive coefficients, λ > 0. The function f + (y) is constructed out of modified Bessel functions of the first and second kind, which are monotonically decreasing and increasing functions respectively. Since the logarithm is a monotonically increasing function, the function f + (y) is also a monotonically increasing function, so f + (y) can have only one zero at most. Because the ratio of the modified Bessel functions K 1 (y)/ I 1 (y) diverges at y → 0 as 1/y 2 , we have f + (y → 0) = −∞. Further, due to the logarithm f + (y) diverges logarithmically for y → ∞, so f + (y) has always exactly one zero (see Fig. 2). Hence we find that there is exactly one NO with a positive coefficient of gerade symmetry. The corresponding NO, φ + g (x), is readily constructed as where C + g is a normalization constant and y 0 is the zero of the function f + (y). The corresponding attains the The negative solutions, λ < 0, can be obtained in a similar manner. The only difference is that the function f − (y) has an infinite amount of zeros, z k , due to the infinite amount of oscillations of the usual Bessel functions. Following the same steps as for the positive orbital, we find that the gerade NOs with negative coefficients are where C − g,k are normalization constants and z k are the zeros of the function f − (y) defined. The corresponding (negative) coefficients are For large k, z k is a large number, so asymptotically the condition f − (z) = 0 becomes tan z − 3π 4 = 2 π Z η + γ + ln(z/2) .
For very large z the logarithm on the right-hand side will diverge and the tangent on the left-hand side is only large In this appendix we show in more detail how the Fourier transforms of the correlation functions of the various model systems were calculated. Let us first consider the correlation functions of the 1D systems. In the case of the double-harmonium the correlation function is a Gaussian, so the Fourier transform should be again a Gaussian. Indeed, when we work it out, we find where we used that the integration contour could be shifted through the complex plane, since there are no poles.
The other systems are three dimensional and have in common that their correlation function only depends on the length, r 12 , so we can already do the integration over the angles For the second integral we need to introduce the convergence factor again, so we evaluate Using this result in the limit α → 0 together with (D1), we find the full Fourier transform of the 3D Hylleraas correlation as in (15).
Finally let us consider the Fourier transform for the correlation function of the system with inverse harmonic interactions. For general α > 0 we have We find that the Fourier transform diverges for k → 0.
In the case of even α, however, it is not clear what happens, since we obtain a division of zero by zero. To asses what the real answer should be, we have to calculate these cases separately. For even α, we can work out the where δ (n) (k) denotes the n th order derivative of the delta-function and we used the previous result for the 3D Fourier transform of a constant function (D1). Combining these results, we find for general α > 0 the Fourier transform as reported (16). Note that this result for general F[r α ](k) can also be used to obtain the Fourier transform of the correlation function of the 3D Hylleraas wavefunction and give the same result.