Abstract

Total energies of three-electron atomic systems along the lithium isoelectronic sequence are computed with the variational Monte Carlo (VMC) method using compact, fully correlated Slater-Jastrow trial wave functions. The calculation addresses the (1s² 2s) ²S^e ground state, the (1s² 2p) ²P^o first excited state and the metastable doubly excited (1s2s2p) ⁴P^o quartet state, and revisits, with an independent stochastic method, the states studied previously with a hypergeometric correlation factor. The trial wave functions combine a double-zeta 1s orbital, Slater-type 2s and 2p orbitals and a Padé-Jastrow correlation factor satisfying the electron-electron cusp conditions exactly; all parameters are optimised by a mixed energy-variance criterion on correlated samples. The VMC energies reproduce the full core plus correlation (FCPC) reference values to within 0.07% or better over the whole sequence, the virial ratio remains within |η - 2| ≤ 0.023 for all twenty-four states, and the 2s → 2p resonance excitation energies agree with the reference differences to about 0.1 eV or better. For the quartet state, the computed excitation energy of the lithium atom, 57.54 eV, reproduces the experimental term value 57.41 eV to 0.13 eV, and a Z-expansion fit of the sequence yields a leading coefficient 0.75022, in excellent agreement with the exact hydrogenic limit 3/4. The variational upper-bound property of the VMC energies is used to examine critically some previously published values for the ²P^o series, and the specific difficulties raised by autoionizing doublet states such as (1sns2) ²S^e, (1s2sns) ²S^e and (1s2snp) ²P^o are discussed in detail.

1. Introduction

The importance of the electron-electron interaction in multi-excitation and multi-ionization processes has been recognised since the earliest variational treatments of two-electron atoms by Hylleraas ¹, and the spectroscopy of doubly excited states has served ever since as the privileged testing ground of correlation theories ^{2, 3, 4, 5, 6, 7}. For two-electron systems, the pioneering electron-impact and synchrotron experiments ^{3, 6, 7} revealed series of autoionizing resonances whose classification required a complete departure from the independent-particle picture ². Systems with three electrons, which realise the four-body Coulomb problem, are considerably more demanding: the lithium atom and its isoelectronic ions constitute the natural prototypes ^{8, 9}, both because they are the simplest open-shell atoms and because hollow and core-excited ionic states can be produced and observed by ionization - excitation in the K shell ^{4, 5}. The ground state of the lithium isoelectronic sequence has been the object of a large body of high-precision work, from the early configuration-interaction study of Weiss ¹⁰ and the Hylleraas-type expansion of Larsson ¹¹ to the very accurate calculations of Yan and Drake ¹², Pipin and Bishop ¹³, Tong, Jönsson and Froese Fischer ¹⁴, Pestka and Woznicki ¹⁵, Perkins ¹⁶, Ho ¹⁷, King ¹⁸, McKenzie and Drake ¹⁹ and Lüchow and Kleindienst ²⁰. Core-excited and doubly excited states have been treated with quasi-projection operators by Bhatia ²¹, with the saddle-point complex-rotation method by David and Chung ²², and with density-functional techniques by Erkoç and Jansen ²³. These benchmark values delimit today, with sub-millirydberg accuracy for the low-lying states, the exact non-relativistic energies of the sequence. In the present work we use the variational Monte Carlo (VMC) method, introduced by McMillan for liquid helium ²⁴ on the basis of the Metropolis algorithm ²⁵, and brought to atomic and molecular physics by Ceperley, Chester and Kalos ²⁶, VMC evaluates the expectation value of the Hamiltonian over an explicitly correlated trial wave function by stochastic integration in the full 3N-dimensional configuration space. Two properties make it attractive for our purpose: the complete freedom in the analytic form of the trial function, which allows the interelectronic distances r_ij to be incorporated exactly through a Jastrow factor obeying the cusp conditions of Kato ²⁷, and the strict variational upper-bound property of the resulting energy, which provides an unambiguous internal check of any published variational value. Explicitly correlated Monte Carlo wave functions have been applied with success to three- and four-electron atoms by Arias de Saavedra et al. ²⁸ and by Buendía, Gálvez and Sarsa ²⁹, and the method is the standard entry point of the more sophisticated diffusion Monte Carlo (DMC) technique ³⁰.

2. Theoretical Framework

The non-relativistic Hamiltonian of a three-electron atomic system with a fixed pointlike nucleus of charge Z reads (a.u.) :

(1)

where r_i is the distance of electron i to the nucleus, r_ij = |r_i - r_j | the interelectronic distance, T the kinetic energy, C the electron-nucleus Coulomb attraction and W the electron-electron Coulomb repulsion. For any normalisable trial wave function Ψ_T possessing the symmetry of the target state, the Rayleigh quotient :

(2)

is an upper bound of the exact energy of the lowest state of that symmetry,

(3)

Let Γ be the symmetry of the state we want (²Sᵉ, ²Pᵒ or ⁴Pᵒ). For any normalisable trial function of symmetry Γ, the Rayleigh quotient (2) cannot fall below , the lowest exact energy of that symmetry, Eq. (3). We apply this bound to each of the three symmetries separately.

The VMC method rewrites the quotient (2) as:

with (4)

where R = (r₁, r₂, r₃) denotes a point of the nine-dimensional configuration space and

(5)

is the local energy. The probability density π(R) is sampled with the Metropolis algorithm ²⁵. The statistical uncertainty of is obtained by the blocking method:

(6)

Here τ is the integrated autocorrelation time, the number of Metropolis steps the walk needs to forget where it started. The M sampled configurations are therefore worth only about M/τ independent ones. And

(7)

is the estimator of the energy over M sampled configurations. The trial wave functions are of the Slater- Jastrow form: a determinantal part fixing the spin and space symmetry, multiplied by a positive, symmetric correlation factor. For the (1s²2s) ²S^e ground state we write, with electrons 1 and 2 of spin up and electron 3 of spin down,

(8)

For the (1s²2p) ²P^o state

(9)

And for the quartet (1s2s2p)⁴P^o state, whose three electrons carry parallel spins,

(10)

Where is Slater-type functions. The Jastrow exponent contains electron-electron and electron-nucleus terms of Padé form :

(11)

The Jastrow exponent has the Padé form (11). The coefficient fixes the electron–electron cusp:

= 1/2 for antiparallel spins and 1/4 for parallel spins, so that in the ⁴Pᵒ state every pair has = 1/4. The only free parameters are g, common to all pairs, and h and k for the electron–nucleus part.

The expectation values of the three terms of Equation (1) are accumulated separately during the sampling, which gives direct access to , and hence to the virial ratio:

(12)

which must equal 2 for the exact eigenstates of the Coulomb Hamiltonian (1). As emphasised already in ³¹, the separate calculation of the kinetic, electron - nucleus and electron - electron terms is a merit of the method; in the stochastic framework this separation is obtained at no extra cost. We therefore tested, for the lithium ground state, the hybrid form:

(13)

with

(14)

u₀ is the l = 0 term and u₁, up to its sign, the l = 1 term of the Legendre expansion of r_ij; the sign of u₁ is absorbed in λ₁.

(15)

with r_< = min(r_i, r_j), r_> = max(r_i, r_j) and θ_ij the angle between r_i and r_j and D is the determinantal part of Eq. (8) and J the Jastrow exponent of Eq. (11), together with a variant in which the electron - electron part of J is removed, so that the correlation is carried by the linear factor alone. All parameters were reoptimised in each case.

3. Results and Discussions

Table 1 gives the outcome. With the linear factor alone, the energy remains acceptable but the variance of the local energy explodes, from 0.079 to 188 a.u.2: a truncated expansion of r_ij in Legendre polynomials cannot cancel the Coulomb singularity 1/r_ij at each coalescence point individually, so the cusp condition is satisfied only in the mean and the local energy develops large fluctuations. At equal computational cost, the statistical error grows by the square root of the variance ratio, here a factor of about fifty. In the full hybrid (13) the optimiser drives λ₀ to 0.003, i.e. it switches the hypergeometric l = 0 term off: once the Jastrow factor handles the coalescence region exactly, the additional linear term brings no measurable energy gain and only raises the variance.

Figure 1 displays the increase in binding energy -E (in Ry) with nuclear charge Z for three states of lithium and its isoelectronic ions. This growth is physically expected because the Coulomb attraction between the nucleus and electrons, proportional to Z², dominates and makes the systems increasingly stable. The three distinct curves represent the ground state (1s²2s) ²S^e, the first excited state (1s² 2p)²P^o, and the metastable quartet state (1s2s2p) ⁴P^o, with the energy gaps between them corresponding to excitation energies. The smoothness and regularity of the curves, with error bars smaller than the symbol sizes, validate the precision of the variational Monte Carlo (VMC) method across the entire sequence. A 1/Z-expansion fit of the quartet curve yields a leading coefficient of 0.75022, in excellent agreement with the theoretical hydrogenic limit of 0.75, confirming the physical consistency of the treatment.

Table 2 collects the present VMC energies of the ground state for 3 ≤ Z ≤ 10, together with the values of the FCPC reference values of Chung and Zhu ³², the SCUNC values of Sakho ³³ and representative Hylleraas-type results ³⁴. The VMC energies lie above the FCPC references for all ions but one, as required by the variational bound (3), with deviations that never exceed 12.3 mRy in absolute value: 10.5 (2.6) mRy for Li, 9.2 (3.2) mRy for Be⁺, decreasing steadily in relative terms from 0.070% at Z = 3 to below 0.002% at Z = 10 (Figure 2). The single apparent exception, F⁶⁺, for which the VMC value lies 1.9 mRy below the reference, is not a violation of Eq. (3): the statistical uncertainty of that point is 5.1 mRy, so that the deviation represents only 0.4 standard deviations and is fully compatible with zero.This point illustrates an important methodological difference between stochastic and deterministic variational methods. In terms of correlation, the VMC ground-state energy of neutral lithium, -7.4728(13) a.u., recovers 88% of the exact correlation energy, defined with respect to the Hartree - Fock limit EHF = -7.43273 a.u. ³⁵ and the essentially exact value -7.478060 a.u. of Yan and Drake ³⁴. This proportion, obtained with only eight variational parameters, illustrates quantitatively the efficiency of the Jastrow description of the electron-electron cusp region: it is comparable to the fraction recovered by configuration-interaction expansions of several tens of terms, in line with the analyses of ^{28, 29}.

The deviations presented in Figure 2 possess two important characteristics. Two features stand out.The VMC energies lie above the FCPC values for every ion except F⁶⁺, where the difference is within one standard deviation. The ²Pᵒ deviations are on average about 1.5 times larger than the ²Sᵉ ones (10.6 against 7.2 mRy), with a maximum of 20.2 mRy at C³⁺. These results indicate that while our simple trial function proves to be effective for the description of the ground state, its ability to take into account the angular correlation of 2p electron is not as good. An issue that can be solved by adding p-type polarization functions to the determinantal part. The non-monotonic behavior near Z=6 might be caused by the delicate balance between the nuclear attraction and the Pauli repulsion of 2p electron from the 1s core.

Table 3 presents the corresponding comparison for the (1s²2p) ²P^o state, which is the lowest state of odd parity and doublet spin and is therefore also rigorously covered by the bound (3). The VMC energies again track the FCPC references from above, with deviations going from 10.8(1.6) mRy for Li down to 1.4(5.2) mRy for Ne⁷⁺, through a maximum of 20.2(2.9) mRy at C³⁺, i.e. ratios E^VMC/E^FCPC between 0.99927 and 0.99999. The present VMC upper bounds, -46.3976(29) Ry, -163.6362(54) Ry and -204.2012(52) Ry, confirm the FCPC references independently in all three cases. Likewise the variational argument is unambiguous for boron.

Differences of total energies give access to the term differences ∆E(2s → 2p), which correspond to the strongest resonance transitions of the sequence: the red doublet of neutral lithium near 670.8 nm, the 313.1 nm doublet of Be⁺, and the ultraviolet doublets of C³⁺ (λλ 1548, 1551 Å), N⁴⁺ (λλ 1239, 1243 Å) and O⁵⁺ (λλ 1032, 1038 Å) that serve as standard diagnostics of hot astrophysical and laser-produced plasmas ³⁶. Table 6 and Figure 4 compare the VMC excitation energies with the differences of the FCPC references. The largest discrepancy along the whole sequence is 0.11 eV (C³⁺), and for five of the eight ions the VMC and reference values agree within one combined standard deviation. For neutral lithium the VMC value, 1.85(4) eV, reproduces the experimental term difference 1.848 eV to 0.2 %. This accuracy on differences is not automatic in a stochastic method: it results from the balanced quality of the trial wave functions attested independently by their virial ratios and variances (Table 4), and it degrades only slowly with Z because the missing correlation energy of the two states evolves in parallel.

Table 4 lists, for the twenty-four computed states, the expectation values ⟨T⟩, ⟨C⟩ and ⟨W⟩ of the three terms of the Hamiltonian, the virial ratio and the variance of the local energy. Several trends emerge:

The virial ratio remains confined to the interval 1.990 ≤ η ≤ 2.023 over the whole table (Figure 3), without any degradation as Z grows. Since η = 2 exactly for eigenstates, while η itself is not stationary with respect to variations of the wave function, the deviation |η - 2| provides a first-order indicator of the residual imbalance between kinetic and potential contents of each trial function. Both the kinetic energy ⟨T⟩ and the electron-electron repulsion ⟨W⟩ increase monotonically with the nuclear charge at fixed state, as observed in Ref. ⁸ for the expectation values ⟨1/r_ij ⟩: the contraction of the charge cloud with increasing Z brings the electrons closer to the nucleus and to each other. In the quartet state the repulsion ⟨W⟩ is systematically smaller than in the doublet states of the same ion for Li, 1.121 a.u. against 2.217 a.u. in the ground state. The variance σ²(E_L) grows approximately as Z² while the total energy grows as Z² as well : the relative spread of the local energy is thus essentially constant along the sequence, which explains why the relative statistical accuracy of the method improves with Z at fixed sample size. The smallest variances of the table are those of the quartet states (e.g. 0.025 a.u. for Li), reflecting again the Fermi-hole suppression of the coalescence regions where the local energy varies most rapidly; this observation quantifies, through the zero-variance principle, the common statement that quartet states of three-electron systems are `easier' to describe than their doublet counterparts ³⁷.

There are two comments to be made about the virial ratios presented in Figure 3. On the one hand, the quartet values (squares) do not sit on one side of 2. Four of them exceed 2.01 (Z = 3, 7, 9, 10), two lie just above 2 (Z = 4, 8) and two just below (Z = 5, 6). The doublets spread over the same kind of range, 1.991 ≤ η ≤ 2.023. In fact the largest deviation in Table 4 is not a quartet at all but the ²Pᵒ state of Be⁺. For the quartet states, η > 2 means that the trial functions selected by the mixed energy-variance criterion are slightly too diffuse. Under a uniform scaling R → sR one has dE/ds = ⟨T⟩(2 − η) at s = 1, so a small contraction, s ≈ η/2, lowers the energy by about

⟨T⟩(η − 2)²/4. For Ne⁷⁺ this gain is 2.9 mHa. On the other hand, in spite of the increase of the variance of the local energy with Z (Table 4), the virial ratio remains remarkably constant at 1.2 % from its true value for all the twenty-four states.

Table 5 presents the VMC total energies of the (1s2s2p) ⁴P^o state and the corresponding excitation energies with respect to the ground state. As recalled in the Introduction, this state is the lowest of its symmetry: the variational bound (3) applies to it exactly as to the ground state. It was first observed in lithium through its optical decay channels and its metastability was elucidated by Feldman and Novick ³⁸; precise term values were subsequently established by beam-foil and beam-laser spectroscopy ³⁷. Two recent studies used this same state as a benchmark for compact correlated wave functions: the analytic ultra-compact functions ⁴⁰ and their VMC optimisation in Ref. ³⁹, which reaches at least 99.9 % of the exact energy with ten parameters or fewer. Our energies sit at a comparable level; for lithium, the value below recovers about 99.8 % of the reference energy deduced from the experimental term value ³⁷, with a determinant-Jastrow form that was not tailored to this state.

For neutral lithium the VMC energy, -5.3581(5) a.u., corresponds to an excitation energy of 57.54 eV above our VMC ground state, to be compared with the experimental term value 57.41 eV ³⁷: the residual 0.13 eV is consistent with the combined missing correlation energies of the two states discussed above. Along the sequence, the excitation energy grows from 57.5 eV (Li) to 895.5 eV (Ne⁷⁺), moving from the extreme ultraviolet into the soft-X-ray range where the corresponding transitions of highly charged ions are observed in tokamak and laser-produced plasmas.

The excitation energies ΔE(2s→2p) displayed in Figure 4 deserve a physical interpretation beyond the numerical agreement. The VMC values follow the reference FCPC curve with a mean relative deviation of 0.6 %, but the largest discrepancy (0.11 eV at Z=6, C³⁺) is not accidental. For this intermediate -Z ion, the 2s and 2p orbitals are nearly degenerate in the independent-particle picture, and the correlation energy difference between the two states is particularly sensitive to the angular part of the Jastrow factor. The near-perfect agreement at Z=3 (Li) and Z=10 (Ne⁷⁺) suggests that our trial function is well adapted to the weakly and strongly correlated regimes, respectively, while the maximum deviation around Z=6 marks the transition zone where the 2s-2p degeneracy amplifies the residual fixed-form bias. This interpretation is consistent with the known difficulty of describing near-degenerate manifolds with compact wave functions, as discussed in Ref. ²⁹.

4. Conclusion

Compact Slater-Jastrow wave functions with at most eight variational parameters, sampled by the Metropolis algorithm and optimised by a mixed energy - variance criterion, have been shown to describe the (1s ²2s) ²S^e, (1s ²2p) ²P^o and (1s2s2p) ⁴P^o states of the entire lithium isoelectronic sequence up to Ne⁷⁺ with a uniform accuracy of 0.001 - 0.07 % on total energies, statistical uncertainties of a few millirydbergs, virial ratios within 1.2 % of the exact value, and resonance excitation energies accurate to a few hundredths of an electronvolt. The stochastic evaluation of the energy, which requires no analytic integration whatever the structure of the correlation factor, appears as the natural complement of the compact correlated wave functions developed in our previous work.

Our energies agree with the independent VMC calculations where they overlap; what this study adds is the uniform three-symmetry protocol and the systematic use of the bound (3) as an audit tool. The hybrid test adds a quantitative reading of this complementarity: the correlation factor and the Jastrow factor describe the same physics, but only the exponential form keeps the variance of the local energy finite under stochastic sampling.

The study of the metastable quartet state extends the approach into the doubly excited spectrum and demonstrates, through the Z-expansion fit whose leading coefficient reproduces the hydrogenic limit 3/4 to 0.03 %, the global physical consistency of the sequence of independently optimised calculations. The autoionizing doublet series, on the other hand, have been argued to require state-specific or complex-scaled extensions of the stochastic method; their treatment, together with the removal of the residual fixed-form bias by diffusion Monte Carlo and the inclusion of relativistic corrections perturbatively for the heavier members of the sequence, constitutes the programme of future work. Beyond the specific results, the present paper illustrates a methodological point that we believe to be of general interest for compact-wave-function studies of few-electron atoms.

References