Benchmark calculations on the lowest-energy singlet , triplet , and quintet states of the four-electron harmonium atom

For a wide range of confinement strengths ω, explicitly-correlated calculations afford approximate energies E(ω) of the ground and low-lying excited states of the four-electron harmonium atom that are within few μhartree of the exact values, the errors in the respective energy components being only slightly higher. This level of accuracy constitutes an improvement of several orders of magnitude over the previously published data, establishing a set of benchmarks for stringent calibration and testing of approximate electronic structure methods. Its usefulness is further enhanced by the construction of differentiable approximants that allow for accurate computation of E(ω) and its components for arbitrary values of ω. The diversity of the electronic states in question, which involve both single- and multideterminantal first-order wavefunctions, and the availability of the relevant natural spinorbitals and their occupation numbers make the present results particularly useful in research on approximate density-matrix functionals. The four-electron harmonium atom is found to possess the (3)P+ triplet ground state at strong confinements and the (5)S- quintet ground state at the weak ones, the energy crossing occurring at ω ≈ 0.0240919.


I. INTRODUCTION
The ongoing quest for methods enabling accurate but computationally inexpensive estimation of correlation energies in Coulombic systems has produced a variety of quantum-chemical approaches that range from straightforward wavefunction-based formalisms to those relying on reduced density matrices or densities themselves.Despite many partial successes, the ultimate goal of uncovering a formalism capable of describing dynamical and nondynamical correlation effects with equal fidelity has not yet been attained and thus novel approaches, such as those based upon the density matrix functional theory (DMFT) 1 and breaking/restoring symmetry of the electronic wavefunctions, 2 continue to appear in chemical literature.The process of developing such approaches invariably involves not only formulation of new approximations and derivation of the respective energy expressions but also testing them against suitable benchmarks.Consequently, the availability of benchmarking and calibration tools for electronic structure methods is as important to quantum chemistry as the knowledge of diverse approximations to the solutions of the Schrödinger equation.
One of the model systems commonly employed in testing of approximate electron correlation theories is the twoelectron harmonium atom, i.e., the species described by the a) Author to whom correspondence should be addressed.Electronic mail: jerzy@wmf.univ.szczecin.plnonrelativistic Hamiltonian Ĥ = 1 2 with N = 2. Electronic properties of this model atom have been thoroughly elucidated with a combination of mathematical analysis and numerical methods. 3,4 t the weak correlation limit, i.e., for large values of the confinement strength ω, the two-electron harmonium atom closely resembles heliumlike systems.On the other hand, its strong-correlation limit of ω → 0 is characterized by purely nondynamical correlation and a complete spatial localization of electrons (i.e., Wigner crystallization).Thus, electron correlation effects with arbitrary relative strengths of the dynamical and nondynamical components can be modeled with a proper choice of ω -a fine tuning impossible to attain with ordinary atoms that, prior to the onset of Wigner crystallization, undergo spontaneous ionization upon their nuclear charges falling below certain critical values.
7][18] Unfortunately, its general usefulness is severely limited by its failure to provide any new information for those approximate electron correlation methods that are already exact for two-electron systems.
Harmonium atoms with N > 2 do not suffer from such a limitation and, due to the abundance of electronic states with diverse spin multiplicities and characters (singleor multideterminantal), 19,20 provide a much richer testing ground for approximate electronic structure formalisms.As in the case of the two-electron system, properties of these atoms possess well-defined asymptotic expansions at the ω → 0 and ω → ∞ limits.Thanks to the absence of unbound states, the coefficients that appear in the individual terms of these expansions are rigorously evaluable in closed forms.
Despite the obvious usefulness of many-electron harmonium atoms in testing, calibration, and benchmarking of approximate methods of quantum chemistry, their properties have not been investigated in sufficient detail.Thus far, only the two lowest-energy 2 P − and 4 P + states of the threeelectron species have been exhaustively studied, the values of their total, kinetic, confinement, and electron-electron repulsion energies being currently known within 1 μhartree for arbitrary values of ω. [21][22][23] The data currently available for the four-electron harmonium atom are limited to the asymptotic energy expressions valid at the ω → ∞ 20 and ω → 0 24 limits, and the results of two independent Monte Carlo studies that are in disagreement on the nature of its ground state. 19,25 rompted by this unsatisfactory state of knowledge, we report in this paper on highly accurate explicitly-correlated calculations on the lowest-energy singlet, triplet, and quintet states of this species.

II. NUMERICAL METHODS
At the weak-correlation limit of ω → ∞, the fourelectron harmonium atom possesses the single-determinantal 3 P + ground state, whereas the lowest-lying singlet and quintet excited states are the multideterminantal 1 D + and the singledeterminantal 5 S − , respectively, their energies bracketing that of another singlet state (a multideterminantal 1 S + ). 19,20 he asymptotic expressions for their energies read 20 where for the 5 S − excited state.At the strong-correlation limit of ω → 0, the three states share the same energy asymptotics of Ref. 24, The individual energy components are given by simple expressions involving E(ω) and its first derivative with respect to the confinement strength ω, namely, 4 and where T(ω), V (ω), and W (ω) stand for the kinetic energy, the potential energy of the harmonic confinement, and the potential energy of the electron-electron repulsion, respectively.For each of the three states under study, differentiable approximants to E(ω) are provided by the robust interpolation scheme, 26 which allows for stitching of the asymptotics (2) and ( 6) with the results of numerical calculations carried out for 19 values of ω belonging to the set {0.001, 0.002, 0.005, 0.1, ..., 1000.}.The ansatz (r 1 , r 2 , r 3 , r 4 , σ 1 , σ 2 , σ 3 , σ 4 ) employed in these calculations has several computational advantages that have been discussed in Ref. 21.In Eq. (10), is the Ith explicitly correlated Gaussian lobe primitive, (σ 1 , σ 2 , σ 3 , σ 4 ) is the appropriate spin function, Â is the fourelectron antisymmetrizer that ensures proper permutational symmetry, and P is the spatial symmetry projector pertaining to the B 1g and A 2g irreducible representations of the D 4h point group for the 1 D + and 3 P + states, respectively, and to the A 1u irreducible representation of the O h point group for the 5 S − state.While admittedly lower than the proper K h spatial symmetry, these symmetries retain transformation properties (such as parities) of the states under study.In particular, they are compatible with the wavefunctions of the 3 P + and 5 S − states that are eigenfunctions of the Lz operator with vanishing eigenvalues and with the wavefunction of the 1 D + state that is a symmetric combination of the eigenfunctions with the eigenvalues of −2 and 2.
The extraction of natural spinorbitals and their occupation numbers {n i (ω)} from the approximate wavefunctions [Eq.( 10)], and the subsequent derivation of the correlation measures, namely, the index of almost idempotency, 27,28 and the correlation entropy, 28,29 have been described previously. 21

III. RESULTS AND DISCUSSION
The computed coefficients of the [ 10 1+2 ] approximants to the total energies of the 1 D + , 3 P + , and 5 S − states of the fourelectron harmonium atom are compiled in Table I.A straightforward calculation employing these approximants leads to the conclusion that the energy crossing between the 3 P + , and 5 S − states occurs at ω ≈ 0.0240919, clarifying the ambiguities in the results of the previously published Monte Carlo studies. 19,25 he present data are somewhat less accurate than those for the three-electron species 21 as reflected in the estimated errors in the energies and their components, which are listed in Table II.At the strong-correlation regime, these estimates exceed those resulting from the classical expression where L 2 is the deviation of the computed expectation value of the L2 operator from its proper value of L(L + 1), and I is the moment of inertia of the Wigner molecule.

A. The 1 D + state
The lowest-energy singlet state of the four-electron harmonium atom possesses the 1 D + symmetry.At the weakcorrelation limit, it originates from a combination of the |ssp m p m | (m = −1, 0, 1).Slater determinants and remains an excited state for all the magnitudes of the coupling strength ω, always lying above its 3 P + triplet counterpart.Inspection of Table III reveals its energies at the confinement strengths of 0.01, 0.5, and 10 being more than 1 mhartree lower than those previously published. 19he multideterminantal nature of the 1 D + state is confirmed by the calculated occupation numbers of its dominant natural spinorbitals, and the asymptotic values of 1/2 and ln 2, respectively, for the index of almost idempotency and the correlation entropy at the ω → ∞ limit (Table IV).

B. The 3 P + state
At the weak-correlation limit, the 3 P + triplet ground state arises from a single |ssp m p m | (m = m = −1, 0, 1) Slater determinant.It remains the lowest-energy state down to ω ≈ 0.0240919, where it becomes the first excited state.As in the case of its singlet counterpart, its previously published energies turn out to be quite inaccurate, their deviations from the  present data (Table V) amounting to, respectively, at least 0.2 and 0.6 mhartree at the confinement strengths of 0.5 and 10, 19 and 0.1, 0.3, and 0.6 mhartree at the those of 0.01, 0.02, and 0.5. 25t 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 = α, β).At The number of basis functions in the expansion (10).
The single-determinantal character of the 3 P + state implies vanishing of both I(ω) and S(ω) at the limit of ω → ∞, which is indeed the case for both spin components of these  indices (Table VII).It is also consistent with the computed occupation numbers of the dominant natural spinorbitals (Table VIII).

C. The 5 S − state
The relatively high energy of the 5 S − state at strong confinements stems from its limiting |sp −1 p 0 p +1 | configuration.Consequently, for large values of ω this quintet is the third excited state that lies above the 1 D + and 1 S + pair of sin- glets.However, the weaker electron-electron repulsion [compare Eqs. ( 3)-( 5)] imparted by its higher multiplicity results in a steeper decrease of its energy upon weakening of the confinement.Therefore, it becomes the ground state of the four-electron harmonium atom for values of ω less than ca.0.0240919 (Table IX).Again, the computed energies are significantly lower (by as much as 0.6 mhartree) than their previously published counterparts. 19,25  contrast to those of the two other states under study, the occupation numbers of the dominant natural spinorbitals  of the lowest-energy quintet fall into very simple patterns as their ordering does not vary with ω (Table X).This behavior is undoubtedly related to both the totally symmetric nature of this state (which implies that all the natural spinorbitals are eigenfunctions of both the L2 and Lz operators) and the absence of the opposite-spin correlation.

IV. CONCLUSIONS
For a wide range of confinement strengths ω, the calculations reported in this paper afford approximate energies E(ω) of the ground and low-lying excited states of the fourelectron harmonium atom that are within few μhartree of the exact values, the errors in the respective energy components being only slightly higher.This level of accuracy constitutes an improvement of several orders of magnitude over the previously published data, 19,25 establishing a set of benchmarks for stringent calibration and testing of approximate electronic structure methods.Its usefulness is further enhanced by the construction of differentiable approximants that allow for accurate computation of E(ω) and its components for arbitrary values of ω.The diversity of the electronic states in question, which involve both single-and multideterminantal first-order wavefunctions, and the availability of the relevant natural spinorbitals and their occupation numbers make the present results particularly useful in research on approximate densitymatrix functionals.
This methodology, which works equally well within both the weak-and strong-correlation regimes, can be readily applied (though at a much higher computational cost) to harmonium atoms involving more than four electrons.One expects such calculations to clarify the ordering of electronic states in those species and generate further benchmarks.

TABLE I .
The parameters of the[ 10  1+2] approximants to the total energies of the 1 D + , 3 P + , and 5 S − states of the four-electron harmonium atom.a a See Refs.(21) and (26) for details.

TABLE II .
The estimated errors in the total energies and their components of the 1 D + , 3 P + , and 5 S − states of the four-electron harmonium atom.a a All values in μhartree; see Ref.(21)for the details of computation.

TABLE III .
The total energy and its components of the 1 D + state of the four-electron harmonium atom.a a All energies in hartree.bThenumber of basis functions in the expansion(10).This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 84.88.147.148On: Wed, 06 Aug 2014 18:42:09

TABLE IV .
The index of almost idempotency I(ω), the correlation entropy S(ω), and the occupation numbers of the dominant natural spinorbitals of the 1 D + state of the four-electron harmonium atom.
a All values for the individual spin components; the designations of the natural spinorbitals refer to their weak-correlation limits.

TABLE V .
The total energy and its components of the 3 P + state of the fourelectron harmonium atom.a a All energies in hartree.b

TABLE VI .
The individual spin components of the kinetic and confinement energies of the 3 P + state of the four-electron harmonium atom.a a All energies in hartree.This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 84.88.147.148On: Wed, 06 Aug 2014 18:42:09

TABLE VII .
The individual spin components of the index of almost idempotency I(ω) and the correlation entropy S(ω) of the 3 P + state of the fourelectron harmonium atom.

TABLE IX .
The total energy and its components of the 5 S − state of the four-electron harmonium atom.a a All energies in hartree.bThenumber of basis functions in the expansion(10).

TABLE VIII .
The occupation numbers of the dominant natural spinorbitals of the 3 P + state of the four-electron harmonium atom.a a The designations of the natural spinorbitals refer to their weak-correlation limits.This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 84.88.147.148On: Wed, 06 Aug 2014 18:42:09

TABLE X .
The index of almost idempotency I(ω), the correlation entropy S(ω), and the occupation numbers of the dominant natural spinorbitals of the 5 S − state of the four-electron harmonium atom.The designations of the natural spinorbitals refer to their weak-correlation limits.The occupation numbers are identical for the p 0 and p ±1 natural spinorbitals, whereas for their d counterparts they deviate (due to the O h symmetry employed in the calculations) by less than 10 −5 from the average values listed in the table. a