The Three-electron Harmonium Atom: the Lowest-energy Doublet and Quadruplet States the Three-electron Harmonium Atom: the Lowest-energy Doublet and Quadruplet States

Calculations of sub-μhartree accuracy employing explicitly correlated Gaussian lobe functions produce comprehensive data on the energy E(ω), its components, and the one-electron properties of the two lowest-energy states of the three-electron harmonium atom. The energy computations at 19 values of the confinement strength ω ranging from 0.001 to 1000.0, used in conjunction with a recently proposed robust interpolation scheme, yield explicit approximants capable of estimating E(ω) and the potential energy of the harmonic confinement within a few tenths of μhartree for any ω ≥ 0.001, the respective errors for the kinetic energy and the potential energy of the electron-electron repulsion not exceeding 2 μhartrees. Thanks to the correct ω → 0 asymptotics incorporated into the approximants, comparable accuracy is expected for values of ω smaller than 0.001. Occupation numbers of the dominant natural spinorbitals and two different measures of electron correlation are also computed.


I. INTRODUCTION
Electronic structures of N-electron harmonium atoms, described by the nonrelativistic Hamiltonian, 1-5 are governed by the magnitude of the confinement strength ω that measures the extent of electron correlation.At the strong-confinement (or the weak-correlation) limit of ω → ∞, harmonium atoms can be regarded as systems of independent harmonic oscillators weakly perturbed by the electron-electron repulsion.Consequently, their energies E(ω) are given by the series expansion 3, 6-8 where the coefficients {A j } depend upon N and the electronic state under consideration.The corresponding electron densities of low-energy states decrease monotonically with the distance from the origin of the coordinate system, making harmonium atoms at the weak-correlation limit very much alike their ordinary counterparts (except for the absence of the nuclear cusp).What sets apart harmonium atoms from the fully Coulombic systems is the strong-correlation (or the weakconfinement) limit of ω → 0. While ordinary atoms undergo autoionization below critical values of nuclear charge, the Hamiltonian (1) admits bounded eigenstates for all values of ω.Therefore, upon the confinement strength becoming a) Author to whom correspondence should be addressed.Electronic mail: jerzy@wmf.univ.szczecin.pl.
vanishingly small, the harmonium atoms evolve into the socalled Wigner molecules, 4,5,9,10 which are spherical Coulomb crystals 11 weakly perturbed by the kinetic energy of their constituent electrons.In the relevant series expansion for the energy, 3, 4, 6 the first two terms are simply the electrostatic energy of the pertinent Coulomb crystal and the zero-point energy of harmonic vibrations about the equilibrium geometry determined by the balance between the harmonic confinement and the electron-electron repulsion.Whereas the dynamical electron correlation dominates at the strong-confinement limit, the Wigner molecules are perfect examples of systems involving the nondynamical correlation with the consequent electron localization and the asymptotic vanishing of occupation numbers of all the natural spinorbitals. 10t is the infinite tunability of the absolute magnitudes and relative extents of the dynamical and nondynamical electron correlations that makes harmonium atoms ideal tools for calibration, testing, and benchmarking of approximate electronic structure methods of quantum chemistry and solid-state physics.In this respect, thanks to their realistic electron densities and smooth variation of electronic properties between the ω → 0 and ω → ∞ limits, these systems hold a definite advantage over the homogeneous electron gas, in which the constant-density ground state is abruptly replaced by the Wigner crystal due to crossing of energy levels at some critical value of the electron density (which in this case serves the role analogous to that of ω).
Despite the availability of exact wavefunctions and energies for certain values of ω 1-3 that has prompted several studies assessing accuracies of diverse methods of quantum chemistry, [12][13][14][15][16][17][18][19][20] the trivial nature of electron correlation in the two-electron harmonium atom limits its usefulness, especially in research on electron-pair and density-matrix functional theories.In contrast, the three-electron species, for which only numerical and asymptotic solutions of the Schrödinger equations are possible, offers insights into more intricate interplays between the same-and opposite-spin correlation effects.
Thus far, the three-electron harmonium atom has been the subject of relatively few investigations.The study of the ω → ∞ limit has resulted in closed-form expressions for the coefficients A 0 , A 1 , and A 2 , and for the 2 P − doublet ground state, and and for the 4 P + quartet first excited state. 6The analysis of the strong-correlation limit of ω → 0 has yielded and for both of the aforementioned states. 4,6 merical computations of the ground-and excited-state energies of the three-electron harmonium atom have also been carried out, involving methods ranging from the Hartree-Fock approximation 21 and low-level electron-correlation approaches 22 to Monte Carlo calculations 23 and the full configuration interaction (FCI) employed in conjunction with extrapolation to the complete basis set (CBS). 24The FCI/CBS study, in which approximate energies of the 2 P − and 4 P + states have been computed for 12 values of ω between 0.1 and 1000.0, has uncovered both advantages and disadvantages of using uncorrelated Gaussian basis sets in electronic structure calculations on harmonium atoms.On one hand, unlike those relying upon explicitly correlated basis sets, such calculations produce not only the actual energies but also their corresponding limits for specific angular momenta that are of much value in assessment of the accuracy of approximate methods involving finite basis sets.On the other hand, they afford energies of inferior accuracy and, due to the emergence of linear dependencies among the basis functions, are not feasible for smaller values of ω.
Prompted by those observations, we have recently embarked upon large-scale calculations that aimed at obtaining highly accurate energies and other electronic properties of the lowest-energy 2 P − and 4 P + states of the three-electron harmonium atom for a wide range of the confinement strengths.The results of these definitive calculations are reported in this paper.

II. NUMERICAL METHODS
Variational energies of the lowest-energy 2 P − and 4 P + states of the three-electron harmonium atom have been computed with accuracy exceeding 1 µhartree for 19 values of ω belonging to the set {0.001, 0.002, 0.005, 0.1, . . ., 1000.0} using the ansatz where (σ 1 , σ 2 , σ 3 ) is the appropriate spin function, Â is the three-electron antisymmetrizer that ensures proper permutational symmetry, P is the spatial symmetry projector, and is the Ith explicitly correlated Gaussian lobe primitive.8][29][30][31] Their popularity stems from the fact that, despite satisfying neither the nucleus-electron nor electron-electron cusp conditions, they are flexible enough to recover most of the correlation energy even with relatively short expansions (12) while allowing for rapid computations of all the necessary integrals.
The wavefunctions of the harmonium atoms are eigenfunctions of both the L2 and Lz operators with the respective eigenvalues of L(L + 1) and m L .The primitives (13) do not correspond to any particular value of L, except for the case when all the vectors {R Ip } vanish (which yields L = 0).In principle, the multiplication of these primitives by either respective angular functions 23 or polynomials in {x p }, {y p }, and {z p } (p = 1, 2, 3) 32 would produce their counterparts pertaining to higher angular momenta.However, the lobe primitives are numerically expedient and efficient when the energies are minimized with respect to all the variational parameters {C I , R Ip , α Ip , β Ipq }. 33 The spatial symmetry projector where the index m runs over all the operations { pm } belonging to a finite symmetry point group and the coefficients {a m } depend on the chosen irreducible representation, acts upon the primitive χ I in the following manner: While lower than K h , the spatial symmetries employed in the actual calculations (A 2g of the D 4h point group for the 4 P + state and u of the D ∞h point group for the 2 P − state), conserve the correct parities. 34In order to minimize the numbers of primitives, the symmetries have been enforced with simple expedients.Thus, for the 4 P + state, the vectors {R Ip } were constrained to the xy plane, yielding the A 2 symmetry of the C 4v point group and the corresponding projector ) Similarly, for the 2 P − state, the A u symmetry of the C i point group, corresponding to the projector has been enforced by restricting the vectors {R Ip } to the z axis.However, for the three smallest values of ω, where the Wigner molecule with its off-center maximum in electron density 4,5,10 emerges, such symmetry restriction results in exceedingly long expansions (12).Lowering the symmetry to A 2u of the D 4h point group and removing all the constraints on the vectors {R Ip } rectifies this problem.
Stitching of the expansions (2) and (3) while accurately reproducing the computed energies is possible with the newly developed robust interpolation scheme. 35The appropriate [ K 1+2 ] approximants read where t is the real-valued solution of the equation The parameters C K, 0 , C K, 1 , C K, K−2 , C K, K−1 , and C K, K are entirely determined by the coefficients A 0 , A 1 , A 2 , B 0 , and B 1 [Eqs.( 4)-( 11)], and whereas the remaining parameters {C K, 2 , . . ., C K, K−3 } are obtained with the least-square fit of the computed energies.Finally, the parameter ω 0 follows from minimization of the maximum absolute difference between the computed energies and those produced by the respective approximant.
The availability of differentiable approximants E(ω) allows accurate estimation of not only the total energy for any value of ω but also its components, namely 3 the kinetic energy the potential energy of the harmonic confinement and the potential energy of the electron-electron repulsion where Natural spinorbitals and their occupation numbers {n i (ω)} have been obtained by diagonalization of the finitematrix representations { IJ } of the respective 1-matrices.These representations have been computed with accurate Gaussian quadratures of the multiple integrals where σ denotes the appropriate collective spin summations and In Eq. ( 30), H k (x) is the kth normalized Hermite polynomial and I ≡ (I x , I y , I z ).The quality of the projection (29) has been optimized by maximizing the trace of with respect to ζ .In practice, employment of the basis functions (30) with 0 ≤ I x , I y , I z ≤ 10 (for the total of 11 3 = 1321 functions) has been found to produce sufficiently accurate natural spinorbitals and their occupation numbers.][38] Two of such measures, namely the index of almost idempotency, 36,38 and the correlation entropy 37,38 S(ω) have been computed.Analysis of the occupation numbers at the weak-and strong-correlation limits 10 leads to the conclusion that the proper approximant for I(ω) has the form where c K, 0 equals the number of electrons in question and t is defined by Eq. ( 19).On the other hand, since the analogous asymptotic series for S(ω) involves logarithmic terms, they are not amenable to stitching with the aforedescribed interpolation.

III. RESULTS AND DISCUSSION
Among all the three-electron Slater determinants built from one-particle wavefunctions of a harmonic oscillator with a given ω, those corresponding to the configurations ssp m (m = −1, 0, 1) possess the lowest energy.However, according to the present calculations, the triply-degenerate 2 P − doublet is the ground state of the three-electron harmonium atom for all the values of ω rather than only at the weak-correlation limit.Inspection of the computed energies (Table I) reveals that the strong-confinement expansion (2) truncated at three terms yields values of E(ω) accurate to within 1% for ω > 0.4.Conversely, the two-term weak-confinement expansion (3) produces energy estimates of such accuracy only for ω < 0.02.At large values of ω, E(ω) is dominated by the kinetic energy and the potential energy of the harmonic confinement.For these energy components, it is possible to compute the individual-spin contributions T s (ω) and V s (ω) (s = α, β, Table II).At the ω → ∞ limit, both the T α (ω) / T β (ω) and V α (ω) / V β (ω) ratios tend to (5/4)+(3/4) (3/4) = 8/3 (Fig. 1).On the other hand, at the limit of ω → 0, where a classical description becomes valid, the ratios become equal to 2, simply reflecting the totalities of particles with individual spins.The sp m p m (m = m = −1, 0, 1) configurations of the ω → ∞ limit give rise to the triply-degenerate 4 P + quartet that remains the first excited state of the three-electron harmonium atom for all the values of ω (Table III).As in the case of the ground state, the truncated weak-and strong-correlation expansions for E(ω) perform rather poorly, producing energy estimates accurate to within 1% only for ω > 0.14 and ω < 0.008, respectively.The doublet-quartet energy splitting, which is significant at strong confinements, vanishes as the confinement strength tends to zero. 4 It is instructive to compare the highly accurate energies of the present study with those published previously. 23,24 n general, excellent (sub-µhartree) agreement with the Monte Carlo results 23 is observed except for the energy of the 2 P − state at ω = 0.5, where the present value is lower by 6 µhartree (Table I).Despite judicious extrapolation employed, the FCI/CBS energies 24 fare much worse, the energy errors steadily rising with a decreasing confinement strength: 2 µhartree at ω = 1000.0,8 µhartree at ω = 10.0, 35 µhartree at ω = 0.5, and 44 µhartree at ω = 0.1 for the 2 P − state; 1 µhartree at ω = 1000.0, 1 µhartree at ω = 10.0, 8 µhartree at ω = 0.5, and 13 µhartree at ω = 0.1 for the 4 P + state.
The quality of the data compiled in Tables I and III is readily assessed by analysing the behavior of the [ K 1+2 ] approximants ( 18) and ( 19) that smoothly interpolate energies between the weak-and strong-correlation limits (Fig. 2).It is found that, after initially decreasing with K in a steady fashion, the maximum absolute deviations of the approximate energies from their "exact" computed counterparts reach temporary plateaus at K = 10, where they amount to 0.3 and 0.4 µhartree for the 2 P − and 4 P + states, respectively.This extent of numerical noise is consistent with the convergence 2. The dependence of the total energy on the confinement strength for the three-electron harmonium atom.The blue dots and lines: the 2 P − state, the red dots and lines: the 4 P + state; the dots, the solid and the dotted lines denote the computed data, the approximants, and the large-ω expansion (2), respectively, whereas the dotted green line denotes the small-ω expansion (3) that is common to both the states. of energies observed in the variational calculations that is suggestive of their accuracy in the order of 0.1 µhartree. 39he high quality of the computed energies and wavefunctions is independently confirmed by the impressive agreement between the "exact" expectation-value energy components and their counterparts calculated according to Eqs. ( 25)-( 28) from the [ 10 1+2 ] energy approximant (with the parameters listed in Table IV) and its first derivative, the maximum absolute deviations amounting to 1.4 (1.4) µhartree for T(ω), 0.4 (0.5) µhartree for V (ω), and 2.0 (1.9) µhartree for W (ω) of the 2 P − ( 4 P + ) state.
In calibration and benchmarking of approximate electron correlation methods, the knowledge of accurate occupation numbers of natural spinorbitals is as important as that of energies.The computed occupation numbers of dominant natural spinorbitals, compiled in Tables V-VII are expected to be accurate to within ca. 10 −5 .In both the cases of the strongly   and VII) natural spinorbitals, the occupation numbers vary smoothly with ω, the leading asymptotic term beyond the constant of zero or one being proportional to ω −1 at the weakcorrelation limit.The vanishing of the occupation numbers at certain values of ω, encountered previously in the case of the two-electron harmonium atom, 3 is not observed.It is not clear at present whether this vanishing is peculiar to the twoelectron species or it is not found in the present study because of the limited range of ω considered.Both the index of almost idempotency [Eq.(31)] and the correlation entropy [Eq.(32)] decrease monotonically with ω (Figs. 3 and 4).In the case of the 2 P − ground state, this behavior is also observed for individual spin components of these measures of electron correlation.The computed values of I(ω) are consistent with their ω → 0 limits that equal the respective numbers of electrons (we note in passing that the "linear entropies" of the two-electron harmonium atom displayed in Fig. 1

IV. CONCLUSIONS
The present calculations of sub-µhartree accuracy produce data on the energy E(ω), its components, and the oneelectron properties of the two lowest-energy states of the three-electron harmonium atom that are both comprehensive and definitive.When employed in conjunction with the recently proposed 35 robust interpolation scheme, the energy computations at 19 values of the confinement strength ω ranging from 0.001 to 1000.0 yield explicit approximants capable of estimating E(ω) and the potential energy of the harmonic confinement V (ω) within a few tenths of µhartree for any ω ≥ 0.001, the respective errors for the kinetic energy T(ω) and the potential energy of the electron-electron repulsion W (ω) not exceeding 2 µhartrees.Thanks to the correct ω → 0 asymptotics incorporated into the approximants, comparable accuracy is expected for values of ω smaller than 0.001.It is also worth mentioning that all the coefficients {C K, k } of the energy approximant are positive, which eliminates roundoff errors due to partial term cancellations.The same is true about the respective approximants for the index of almost idempotency.
Comparison of the computed energies with those previously published reveals some discrepancies that apparently stem from the substandard performance of the extrapolation employed in the workout of the FCI/CBS data 24 and an error of unknown origin in one instance of the results of Monte Carlo calculations. 23Overall, the excellent performance of both the explicitly correlated Gaussian lobe functions and the robust interpolation scheme in calculations on the threeelectron harmonium atom bode well for analogous studies of systems with larger numbers of particles.
FIG.2.The dependence of the total energy on the confinement strength for the three-electron harmonium atom.The blue dots and lines: the 2 P − state, the red dots and lines: the 4 P + state; the dots, the solid and the dotted lines denote the computed data, the approximants, and the large-ω expansion (2), respectively, whereas the dotted green line denotes the small-ω expansion (3) that is common to both the states.

FIG. 3 .FIG. 4 .
FIG.3.The dependence of the index of almost idempotency on the confinement strength for the three-electron harmonium atom.Blue dots and lines: the α electrons of the 2 P − state, red dots and lines: the β electrons of the 2 P − state, green dots and lines: all the electrons of the 2 P − state, orange dots and lines: all the electrons of the 4 P + state; the dots and the solid lines denote the computed data and the approximants, respectively.

TABLE I .
The total energy and its components of the 2 P − ground state of the three-electron harmonium atom.a a All energies in hartree.bThenumber of basis functions in the expansion(12).c Calculations within the D 4h point group, see the text for explanation.

TABLE II .
The individual-spin components of the kinetic and confinement energies of the 2 P − ground state of the three-electron harmonium atom.a a All energies in hartree.

TABLE III .
The total energy and its components of the 4 P + first excited state of the three-electron harmonium atom.a a All energies in hartree.bThenumber of basis functions in the expansion(12).

TABLE IV .
The parameters of the[ 10  1+2] approximants for the total energies of the 2 P − and 4 P + states of the three-electron harmonium atom.a a See Eqs.(18) and(19).

TABLE V .
The occupation numbers of the strongly occupied natural spinorbitals of the 2 P − and 4 P + states of the three-electron harmonium atom.

TABLE VII .
The occupation numbers of the dominant weakly occupied natural spinorbitals of the 4 P + first excited state of the three-electron harmonium atom.
of Ref.38violate such asymptotics and thus are obviously in error).The approximants(33)with K = 8

TABLE VI .
The occupation numbers of the dominant weakly occupied natural spinorbitals of the 2 P − ground state of the three-electron harmonium atom.

TABLE VIII .
The parameters of the approximants(33)for the index of almost idempotency of the 2 P − and 4 P + states of the three-electron harmonium atom.