Abstract

We present a variational study of the ground-state energies of the four-electron members of the beryllium isoelectronic sequence (Be, B⁺, C²⁺) using explicitly correlated Hylleraas-type trial wave functions formulated within a complex-scaling (complex-rotation) framework. The trial function incorporates all six interelectronic distances of the four-electron problem together with a restricted set of nonlinear variational parameters, and the nonrelativistic Hamiltonian is evaluated after the coordinate transformation $r\to r{\text{e}}^{\left\{\text{i}\theta \right\}}$ . Analytic matrix elements were derived and evaluated numerically with the computer-algebra system Maxima. For the (1s²2s²) ¹S ground state, the resulting energies are compared with Hartree-Fock and variational Monte Carlo values, and with high-accuracy Hylleraas-configuration-interaction (Hy-CI) benchmark energies available for the Be-like sequence. The present energy for neutral beryllium agrees with the Hy-CI reference to within 8 × 10⁻⁵ a.u., and the B⁺ energy agrees to within 1.1 × 10⁻³ a.u. For C²⁺, however, a deviation of about 0.41 a.u. (≈11 eV) is found, revealing a clear limitation of the present implementation for the more highly charged member of the sequence. The evolution of the optimized nonlinear parameters along the sequence is examined and is found to coincide with this breakdown, without our data allowing any causal conclusion to be drawn. The complex-scaling formalism used to construct the trial function is also discussed; the nonzero imaginary parts obtained for the nominally bound ground states are not physically defensible as resonance widths and are therefore not analyzed here as a physical result. The results show close numerical agreement with benchmark values for the lower members of the sequence, while also revealing implementation limitations that prevent the present calculations from being regarded as fully validated variational ground-state energies. These findings identify the numerical and formal improvements required before the present framework can be extended reliably to genuine four-electron autoionizing resonances.

Keywords

Electron Correlation, Hylleraas Wave Function, Complex Scaling, Four-Electron Atoms, Beryllium Isoelectronic Sequence, Variational Method

Share and Cite:

1. Introduction

Atoms and ions containing more than one electron cannot, in general, be treated exactly: the interelectronic repulsion term couples the motion of every pair of electrons and prevents an analytic separation of the Schrödinger equation for N > 1, except in a few special cases. Among multi-electron systems, those with four electrons: the beryllium atom and its isoelectronic sequence B⁺, C²⁺, N³⁺, and so on, occupy a position of particular interest. They are complex enough to exhibit genuinely non-trivial electron correlation, notably through the coupling of a compact 1s² core to a more diffuse 2s² valence shell, yet remain small enough to be treated with explicitly correlated variational methods that become computationally prohibitive for larger systems. For this reason the beryllium isoelectronic sequence has long served as a benchmark for testing correlated wave-function methods and has attracted sustained theoretical attention, from early configuration-interaction studies to present-day sub-microhartree calculations [1]-[3].

Beyond their role as a theoretical testing ground, four-electron ions are also of direct spectroscopic and astrophysical relevance: doubly and singly excited configurations of Be-like ions decay by autoionization and radiative channels that are observed in laboratory plasmas and in astrophysical sources, and accurate energies and correlation parameters for the parent ground and low-lying states are a prerequisite for interpreting such spectra [1]. This provides additional motivation, beyond benchmarking, for examining how well a compact, explicitly correlated variational description performs across the sequence as the nuclear charge increases.

Historically, approximate mean-field treatments such as the Hartree [4] and Hartree-Fock [5] methods, and Slater’s semi-empirical screening-constant approach [6], provided the first computationally tractable descriptions of multi-electron atoms; all three neglect, to varying degrees, the instantaneous correlation between electrons that becomes essential for a quantitatively accurate description of systems such as the Be-like sequence.

Two complementary theoretical ingredients are relevant to the present work. The first is the explicit treatment of electron correlation through interelectronic coordinates, introduced by Hylleraas for the helium atom [7] and subsequently generalized to more elaborate trial functions incorporating several interelectronic distances simultaneously. For two- and three-electron systems, Hylleraas-type and Hylleraas-configuration-interaction (Hy-CI) expansions routinely reach spectroscopic accuracy. For four-electron systems such as beryllium and its isoelectronic ions, explicitly correlated approaches remain computationally demanding but have nonetheless produced benchmark nonrelativistic energies of very high precision, notably through the systematic Hy-CI study of the Be-like ¹S ground states by Sims and Hagstrom [2] [3] and its extension to excited ¹S states [8] [9], through compact Hylleraas expansions for the Li⁻ to Ne⁶⁺ sequence by King et al. [10], through explicitly correlated Gaussian (ECG) calculations [11], and through the high-precision relativistic and QED studies of the beryllium spectrum by Puchalski, Komasa and Pachucki [12]-[14].

The second ingredient is the complex-scaling (or complex-rotation) transformation, in which the electronic coordinates are analytically continued according to $r\to r{\text{e}}^{\left\{\text{i}\theta \right\}}$ . Its mathematical foundations were established by Aguilar and Combes [15] and by Balslev and Combes [16] for one- and many-body Schrödinger operators with dilation-analytic potentials, and its use for isolating resonance poles was formalized by Simon [17]. Complex scaling has since become a standard tool in atomic and molecular physics for the direct determination of resonance energies and widths, as reviewed by Ho [18] and by Moiseyev [19]: under the transformation, bound-state poles remain fixed on the real axis, continuum branch cuts rotate into the lower half of the complex-energy plane by an angle 2θ, and resonance poles that were hidden beneath the continuum on the real axis

become exposed as isolated complex eigenvalues, $E=E_r-\frac{\text{i}}{2}\Gamma$ .

Within our research group, explicitly correlated Hylleraas-type wave functions combined with complex scaling have previously been applied to doubly excited resonance states of two- and three-electron systems (helium- and lithium-like ions and their isoelectronic series) [20]-[25], building in part on the screening-constant-by-unit-nuclear-charge (SCUNC) and modified-atomic-orbital-theory (MTOA) formalisms developed within the same group for the semi-empirical treatment of multi-electron atomic spectra [26]-[29]. Extending this combined methodology, explicit Hylleraas-type correlation together with a complex-scaling formulation of the Hamiltonian, to four-electron beryllium-like systems is a natural continuation of this program, but it also raises the bar considerably: the wave function must now depend on six interelectronic distances instead of one, and the numerical evaluation of the resulting matrix elements is correspondingly more demanding.

The calculations reported here originate from, and substantially reanalyze, an earlier Master’s dissertation carried out within our group; the present article restructures that original numerical work around a benchmark comparison with high-accuracy literature calculations and a critical assessment of its results, rather than reproducing the dissertation’s presentation. In this work, we examine how well such a four-electron, Hylleraas-correlated, complex-scaled variational formulation reproduces the (1s²2s²) ¹S ground-state energies of Be, B⁺ and C²⁺, and how the associated nonlinear variational parameters evolve with the nuclear charge Z along the sequence. We deliberately restrict the scope of the present analysis to ground-state energies, which is the quantity for which the calculation can be meaningfully validated against independent, high-accuracy benchmarks. It is worth stating explicitly why a complex-scaled formulation is used at all for states that are, by construction, bound: the rotated Hamiltonian of Section 2.2 is adopted here not because it improves the ground-state energies themselves: a real, unscaled variational calculation would in principle suffice for that purpose, but because it defines, within a single computational framework, the same trial wave function and matrix elements that would subsequently be needed to locate genuine autoionizing resonances of these four-electron systems. We do not claim, on the basis of the present data, to provide reliable resonance parameters for these systems: as discussed in Section 4, the complex-scaling procedure applied here to a nominally bound state produces nonzero imaginary energy components that are not physically defensible as genuine resonance widths, and we treat this point with the caution it requires rather than presenting it as a result. The remainder of the paper is organized as follows. Section 2 summarizes the theoretical framework: the four-electron Hamiltonian, its complex-scaled form, and the Hylleraas-type trial wave function. Section 3 describes the computational implementation. Section 4 presents and discusses the results, including a comparison with Hartree-Fock, variational Monte Carlo and Hy-CI reference values, an analysis of the deviation observed for C²⁺, and the evolution of the optimized parameters with Z. Section 5 concludes and outlines the numerical work required before the present formulation could be reliably extended to genuine four-electron autoionizing resonances.

2. Theoretical Framework

The theoretical ingredients used in this work: the nonrelativistic four-electron Hamiltonian, its complex-scaled form, and the explicitly correlated Hylleraas-type trial wave function, are outlined below at the level of detail needed to interpret the results of Section 4. Full analytic expressions for the matrix elements, which are lengthy, are not reproduced here; they follow the standard decomposition into radial and angular integrals described for related two- and three-electron problems in Refs. [20]-[23] and are available from the authors.

2.1. Nonrelativistic Hamiltonian

For a four-electron atomic system of nuclear charge Z (an infinitely heavy, fixed nucleus is assumed), the nonrelativistic Hamiltonian in atomic units reads

$H=-\frac{1}{2}{\displaystyle \sum }_{i=1}^4{\nabla }_i^2-Z{\displaystyle \sum }_{i=1}^4\frac{1}{r_i}+{\displaystyle \sum }_{i<j}\frac{1}{r_{ij}}$ (1)

where $r_i$ is the distance of electron i from the nucleus and $r_{ij}=\left|r_i-r_j\right|$ is the distance between electrons i and j. The first sum is the kinetic energy of the four electrons, the second is the electron-nucleus attraction, and the third is the electron-electron repulsion, which prevents an exact separation of the Schrödinger equation $H\Psi =E\Psi$ for $N>1$ .

All energies in this work are expressed in Hartree atomic units (1 a.u. = 27.2114 eV [30]), with lengths in units of the Bohr radius $r_0$ . Nuclear charges are labelled by $Z$ ($Z=4,5,6$ for Be, B⁺ and C²⁺, respectively), and the four electrons are labelled $i=1,2,3,4$ , with electrons 1 and 2 associated with the 1s-like core and electrons 3 and 4 with the 2s-like valence shell in the trial wave function of Section 2.3.

2.2. Complex-Scaling Formalism

Complex scaling proceeds by analytically continuing every electronic coordinate according to $r_i\to r_i{\text{e}}^{\left\{\text{i}\theta \right\}}$ , with $\theta$ a real rotation angle [15]-[17]. Applying this transformation to Equation (1) gives the non-Hermitian, complex-scaled Hamiltonian

$H\left(\theta \right)=-\frac{1}{2}{\text{e}}^{-2\text{i}\theta }{\displaystyle \sum }_i{\nabla }_i^2-Z{\text{e}}^{-\text{i}\theta }{\displaystyle \sum }_i\frac{1}{r_i}+{\text{e}}^{-\text{i}\theta }{\displaystyle \sum }_{i<j}\frac{1}{r_{ij}}$ (2)

As illustrated schematically in Figure 1, the spectrum of $H\left(\theta \right)$ separates the different classes of stationary states of the original problem. Bound-state poles, associated with square-integrable wave functions and real, negative energies, remain unchanged by the rotation and stay on the real axis. The continuum, associated with scattering states, is represented before rotation by a branch cut running along the positive real axis from each ionization threshold; under complex scaling this cut rotates into the lower half-plane by an angle 2θ. Resonance poles, which for $\theta =0$ lie on the unphysical (second) Riemann sheet beneath the continuum and are therefore not directly accessible, become exposed once the cut has rotated sufficiently far below them. Formally, a genuine resonance corresponds to a complex eigenvalue

$E_{res}=E_r-\frac{\text{i}}{2}\Gamma$ (3)

whose real part $E_r$ gives the resonance position and whose imaginary part is related to the natural width Γ and to the state lifetime $\tau =\frac{\hslash }{\Gamma }$ . A necessary signature of a true resonance is that its complex eigenvalue becomes stationary ($\frac{\partial E}{\partial \theta }\approx 0$ ) over a finite range of $\theta$ once the rotation angle exceeds roughly$\arg \left(E_{res}\right)/2$ and the basis is sufficiently large: the well-known θ-trajectory stabilization criterion [18] [19]. For a genuinely bound state, by contrast, Γ must vanish identically, independently of $\theta$ , once the calculation has converged. This distinction is central to the discussion of Section 4.6.

2.3. Hylleraas-Type Correlated Trial Wave Function

The trial wave function used for the (1s²2s²) ¹S ground state of the four-electron system is written, following the general architecture illustrated in Figure 2, as the product of a spatial and an angular part.

Figure 1. Principle of complex scaling: bound-state poles (blue) remain on the real axis; the continuum branch cut (red) rotates by an angle 2θ under the transformation $r\to r{\text{e}}^{\left\{\text{i}\theta \right\}}$ ; a resonance pole initially hidden beneath the continuum (left) becomes exposed once the cut has rotated below it (right). Schematic, redrawn after the general presentation of Refs. [18] [19].

${\Psi }_I\left({\overrightarrow{r}}_1,{\overrightarrow{r}}_2,{\overrightarrow{r}}_3,{\overrightarrow{r}}_4\right)={\Phi }_{npq}\left({\overrightarrow{r}}_1,{\overrightarrow{r}}_2,{\overrightarrow{r}}_3,{\overrightarrow{r}}_4\right)\times {\chi }_{ijkm}\left({\overrightarrow{r}}_1,{\overrightarrow{r}}_2,{\overrightarrow{r}}_3,{\overrightarrow{r}}_4\right)$ (4)

The spatial factor Φ is built from incomplete hydrogenic radial functions for the two 1s core electrons and the two 2s valence electrons, multiplied by an exponential factor $\exp \left[-\left({\lambda }_1{r}_1+{\lambda }_2{r}_2+{\lambda }_3{r}_3+{\lambda }_4{r}_4\right)\right]$ and by powers $r_{12}^n{r}_{13}^p{r}_{14}^p{r}_{23}^t{r}_{24}^t{r}_{34}^r$ of the six interelectronic distances that carry the explicit electron-correlation content of the wave function (labelled in Figure 2). The nonlinear scaling parameters ${\lambda }_i$ are related to the nuclear charge and to a set of variational exponentsthrough ${\lambda }_i=\frac{Z}{{\alpha }_i{n}_i{r}_0}$ , with $r_0$ the Bohr radius. The angular factor $\chi$ combinesspherical harmonics for each electron with the appropriate coupling of the two pairs of electrons (1, 2) and (3, 4) so that the overall trial function transforms as a singlet ¹S state. This general functional form follows the family of special Hylleraas-type wave functions that has been used successfully for two- and three-electron resonance problems within our group [20]-[23]. We now give the explicit form of the spatial factor Φ and of the interelectronic correlation factors $r_{ab}^c$ that enter it.

The spatial part is written as a sum of two terms built from incomplete hydrogenic radial functions for the core (1s-like) and valence (2s-like) electrons, given here in a form condensed for readability from the original variational implementation,

$\Phi =\left(r_1{r}_2{r}_3{r}_4\right){\ell }_1{\displaystyle \sum }_{\nu }{\left[n_1^2{r}_0^2{\lambda }_1^2\left(r_1{r}_2{r}_3{r}_4\right)\right]}^{\nu }+\left(r_1{r}_2{r}_3{r}_4\right){\ell }_3{\displaystyle \sum }_{{\nu }^{\prime }}{\left[n_3^2{r}_0^2{\lambda }_3^2\left(r_1{r}_2{r}_3{r}_4\right)\right]}^{{\nu }^{\prime }}$ (5)

where $r_0$ is the Bohr radius, ${\ell }_1$ and ${\ell }_3$ are the orbital angular momenta associated with the core and valence radial functions respectively, and $n_1$ , $n_3$ are

Figure 2. Schematic architecture of the Hylleraas-type trial wave function for the four-electron system, showing the four electron-nucleus distances $r_1\text{-}{r}_4$ and the six interelectronic distances $r_{12},r_{13},r_{14},r_{23},r_{24},r_{34}$ that carry the explicit electron-correlation content of Φ.

the corresponding principal quantum numbers. This spatial factor is multiplied by the correlation factors $r_{12}^n{r}_{13}^p{r}_{14}^p{r}_{23}^t{r}_{24}^t{r}_{34}^r$ , with $\left(n,p,q\right)$ nonnegative integer exponents, and by the exponential $\exp \left[-\left({\lambda }_1{r}_1+{\lambda }_2{r}_2+{\lambda }_3{r}_3+{\lambda }_4{r}_4\right)\right]$ . Each of the six interelectronic distances entering these correlation factors is built, for a given pair of electrons, from the same elementary two-body geometry recalled in Figure 3: for two electrons at distances $r_a$ , $r_b$ from the nucleus and separated by an angle ${\theta }_{ab}$ , the interelectronic distance $r_{ab}$ obeys $r_{ab}^2=r_a^2+r_b^2-2r_a{r}_b\cos {\theta }_{ab}$ , the relation on which the Hylleraas coordinates $\left(s,t,u\right)=\left(r_a+r_b,r_b-r_a,r_{ab}\right)$ originally introduced for the two-electron problem are based [7].

Each interelectronic correlation factor $r_{ab}^c$ (for $r_a\le r_b$ ) is expanded in Legendre polynomials and hypergeometric functions of the individual electron coordinates,

$r_{ab}^c=4\pi {\displaystyle \sum }_{l=0}^{\infty }\left[\frac{1}{2\ell +1}\right]a_l^c\left(r_a,r_b\right)Y_l\left({\Omega }_a\right)Y_l\left({\Omega }_b\right)$ (6)

with $a_l^c\left(r_a,r_b\right)$ expressed through the Gauss hypergeometric function ₂F₁, following the standard construction used for two- and three-electron Hylleraas-type resonance wave functions in Refs. [20]-[23] from which the present four-electron form was adapted. The angular part $\chi$ combines products of s- and higher-ℓ spherical harmonics for each electron, coupled pairwise for electrons (1, 2) and (3, 4) with intermediate angular momenta L₁₂ and L₃₄ and spins S₁₂ and S₃₄, and antisymmetrized between the two pairs, so that the complete trial function $\Psi =\Phi \times \chi$ transforms overall as a singlet ¹S state appropriate to the (1s²2s²) ground-state configuration. The full analytic expressions for all matrix elements of the normalization, kinetic-energy, electron-nucleus and electron-electron operators built from this wave function are lengthy and are not reproduced here; they are available from the authors upon request.

Figure 3. Hylleraas coordinates for a single electron pair: radial distances $r_a$ , $r_b$ from the nucleus, interelectronic distance $r_{ab}$ , and the angle ${\theta }_{ab}$ between the two position vectors, related by $r_{ab}^2=r_a^2+r_b^2-2r_a{r}_b\cos {\theta }_{ab}$ [7]. This elementary two-body geometry underlies each of the six correlation factors entering the four-electron trial wave function of Figure 2 and Equation (6).

We note that this pairwise coupling, antisymmetrized within each electron pair (1, 2) and (3, 4) and between the two pairs as a whole, does not by itself enforce antisymmetry under exchange of a core electron with a valence electron (e.g., electrons 2 and 3); a fully antisymmetric treatment would require either an explicit sum over all 4! permutations with the appropriate signs, or a configuration-state-function/Slater-determinant construction built from the same spatial factors. This restriction is a further simplification adopted for computational tractability and is noted here explicitly as a limitation of the present angular ansatz, alongside those already discussed in Section 4.5.

For the ground state considered here, the number of independent nonlinear parameters was reduced by imposing $\alpha =\beta$ for the two 1s-like core exponents and $\gamma =\delta$ for the two 2s-like valence exponents, together with a single additional parameter a controlling the interelectronic correlation factors. This restriction, adopted as a practical simplification of the variational search rather than derived from a formal symmetry argument, reduces the dimensionality of the nonlinear optimization at the cost of not treating the four electrons in a fully independent manner; its consequences are discussed further in Section 4.5.

2.4. Variational Procedure

The nonlinear parameters ($\alpha ,\gamma ,a$ ), together with the linear expansion coefficients implicit in the basis-size parameter $\Omega =i+j+k+m+n+p+q$ , were determined at $\theta =0$ , where the Hamiltonian of Equation (1) is Hermitian and the Rayleigh quotient

$E\left[\Psi \right] = \langle \Psi |H\left(\theta \right)|\Psi \rangle /\langle \Psi |\Psi \rangle$ (7)

Numerical implementation. In the Maxima implementation used for the calculations, the scalar parameter a reported in Table 1 enters the phase factors exp(−ia) and exp(−2ia); no separate numerical variable θ is introduced. Accordingly, a is treated here as the phase parameter of the implemented complex transformation rather than as an independent interelectronic-correlation exponent. The implemented Coulomb prefactors and the assignment of the interelectronic exponents also follow the conventions of the computational expressions used to generate the reported values. For this reason, the numerical results below are discussed primarily through direct comparison with independent literature values, while the strict Ritz upper-bound criterion is reserved for an explicitly Hermitian $\theta =0$ formulation.

is a rigorous upper bound to the exact ground-state energy in accordance with the Ritz variational principle [31]; this bound is tightened by minimizing $E\left[\Psi \right]$ with respect to ($\alpha ,\gamma ,a$ ). The complex-scaling transformation ($\theta \ne 0$ ) is applied afterward to the resulting trial function, as described in Section 3. For $\theta \ne 0$ the complex-scaled Hamiltonian $H\left(\theta \right)$ of Equation (2) is non-Hermitian, so the corresponding Rayleigh quotient no longer carries the upper-bound property of the real variational principle; the complex-scaled matrix elements are accordingly evaluated here as a generalized extension of the same analytic framework rather than as a conventional variational calculation. This distinction is used in Section 4.6 to interpret the complex energies obtained for $\theta \ne 0$ . All required matrix elements: normalization, kinetic energy, electron-nucleus attraction and electron-electron repulsion, were derived analytically in closed form and evaluated numerically as described in Section 3.

The present implementation evaluates the Rayleigh quotient of Equation (7) for a single basis function at fixed ($i,j,k,m,n,p,q$ ) rather than diagonalizing a set of such functions spanning a given Ω; the linear expansion coefficients are therefore not varied in the calculations reported here, and Ω labels a single dominant configuration rather than the dimension of a diagonalized basis. Extending the calculation to a genuine multi-configuration expansion, solved as a generalized eigenvalue problem, discussed in Section 5, would provide a systematic route for examining the C²⁺ discrepancy discussed in Section 4.5.

3. Computational Implementation

The analytic matrix elements described in Section 2 were implemented and evaluated using the open-source computer-algebra system Maxima, version 5.21.1 [32], which provides symbolic integration, matrix diagonalization and nonlinear equation solving within a single environment. For each system (Be, B⁺, C²⁺) the calculation proceeds in three stages. First, the real part of the Hamiltonian ($\theta =0$ ) is stabilized with respect to the nonlinear parameters α, γ and a, so as to identify a set of parameters for which the computed energy is stationary under small parameter variations. Second, the complex-scaling transformation is applied ($\theta \ne 0$ ) to the wave function obtained from the first stage. Third, the basis-size parameter Ω is increased and the first two stages are repeated to assess numerical stability with respect to the size of the variational expansion.

Computational procedure. The analytical expressions are evaluated in Maxima using the parameter a directly in the phase factors exp(−ia) and exp(−2ia). The reported calculations therefore correspond to this implemented complex-phase parametrization rather than to two numerically separate stages consisting of an independent θ = 0 optimization followed by rotation. This distinction is taken into account in the interpretation of Table 1, Table 2.

The angle of rotation θ and the exponents ($i,j,k,m,n,p,q$ ) defining a given basis function are initialized so as to satisfy $i,j,k,m,n,p,q\ge 0$ and $i+j+k+m+n+p+q=\Omega$ , and the nonlinear variational parameters are then adjusted to minimize the energy as described in Section 2.4. A full calculation for a single system required approximately 15 to 30 minutes of computation time within Maxima, reflecting the complexity of the twelve-dimensional analytic integrals that must be evaluated for a four-electron correlated wave function of this type; this is markedly longer than for the corresponding two- or three-electron calculations carried out previously within the group [20]-[23], and highlights a practical limitation, rather than a matter of principle, of the present symbolic-computation approach when it is extended to a larger basis or to excited-state configurations. Migrating the analytic matrix elements to a compiled numerical environment is a natural next step, discussed further in Section 5.

4. Results and Discussion

4.1. Optimized Nonlinear Parameters

Table 1 lists the optimized nonlinear parameters ($\alpha =\beta ,\gamma =\delta ,a$ ) obtained for the (1s²2s²) ¹S ground state of Be, B⁺ and C²⁺ at the basis size used in this work. The core-electron exponent α increases monotonically from 2.7800 (Be, $Z=4$ ) to 4.1900 (B⁺, $Z=5$ ) to $4.4190$ (C²⁺, $Z=6$ ), consistent with the expected contraction of the 1s-like core orbitals as the nuclear charge increases. The correlation parameter a also increases monotonically, from 0.3500 to 0.4500 to 0.5500. The valence-electron exponent γ, by contrast, is essentially unchanged between Be and B⁺ (0.5322 in both cases) but increases sharply to 1.6416 for C²⁺: a qualitatively different behavior from the smooth trends followed by α and a, and one that we return to in Section 4.5.

The tabulated quantity a is the phase parameter used in the factors exp(−ia) and exp(−2ia). Its numerical evolution along the sequence therefore characterizes the implemented complex transformation and should not, by itself, be interpreted as direct evidence for an increasing strength of electron correlation.

Table 1. Optimized nonlinear variational parameters for the (1s²2s²) ¹S ground state of the four-electron beryllium-like sequence.

4.2. Ground-State Energies

Table 2 compares the present ground-state energies with Hartree-Fock and variational Monte Carlo (VMC) values reported by Doma, Roston, Ahmed and Sen [33]: whose study focuses on beryllium and its isoelectronic ions under spherical confinement, but also reports the free-atom (unconfined, box radius $r_0\to \infty$ ) limit as a baseline for comparison, which is the quantity used here, and with the high-accuracy Hylleraas-configuration-interaction nonrelativistic benchmark energies of Sims and Hagstrom for the Be-like ¹S ground states, quoted to better than 20 nanohartree [2] [3]. For neutral beryllium, the present result (−14.667280 a.u.) agrees with the Hy-CI reference (−14.667356 a.u.) to within 7.6 × 10⁻⁵ a.u., an accuracy comparable to that of the VMC calculation and substantially better than Hartree-Fock. For B⁺, the present result (−24.350017 a.u.) agrees with the Hy-CI reference (−24.348884 a.u.) to within 1.13 × 10⁻³ a.u., again a reasonable level of agreement for a compact variational expansion of this type, though about an order of magnitude less accurate than for neutral beryllium.

For C²⁺, however, the present result (−36.949490 a.u.) differs from the Hy-CI reference energy (−36.534852 a.u. [3]) by 0.4146 a.u., i.e. approximately 11.3 eV [30], roughly two to three orders of magnitude larger than the deviations found for Be and B⁺, and far in excess of what any reasonable variational approximation would be expected to produce for a smoothly varying isoelectronic sequence. This magnitude of deviation should not be interpreted as evidence of unusually strong correlation effects for the C²⁺ ion; it instead indicates a specific numerical shortcoming of the present calculation or optimization for this member of the sequence, whose precise origin has not been established. We report the C²⁺ result together with this caveat rather than omit it, since it provides an indication of the present limitation of the implementation at higher nuclear charge (see Section 4.5).

Table 2 compares the calculated values directly with independent Hy-CI, VMC, and Hartree-Fock benchmarks. The Be value is close to the Hy-CI reference, whereas the B⁺ and especially C²⁺ values lie below the corresponding high-accuracy nonrelativistic references. Since the implemented calculation contains a nonzero complex phase through a, these values are not θ = 0 Hermitian Ritz upper bounds. They are therefore used here as comparative numerical outputs of the complex-phase formulation rather than as strict variational bounds.

Table 2. Ground-state energies (a.u.) of the four-electron beryllium-like sequence: present work compared with Hartree-Fock and variational Monte Carlo values and with the Hylleraas-configuration-interaction reference energies of Sims and Hagstrom. ΔE = E(present) - E(reference).

4.3. Energy Trends and Deviation from the Reference Values

Figure 4 shows the ground-state energies of Table 2 as a function of nuclear charge Z, together with the Hartree-Fock, VMC and Hy-CI reference values. On this scale the four sets of values follow closely the same smooth trend, and the deviation identified for C²⁺ is not visually apparent, since it corresponds to only about 1% of the total binding energy at Z = 6. Figure 5 makes the same comparison explicit on a logarithmic scale by plotting $\left|\Delta E\right|=\left|E_{\text{present}}- E_{\text{reference}}\right|$ for each system. The deviation increases from 7.65 × 10⁻⁵ a.u. (Be) to 1.13 × 10⁻³ a.u. (B⁺) to 4.15 × 10⁻¹ a.u. (C²⁺), a jump of more than two orders of magnitude between B⁺ and C²⁺ that is not present between Be and B⁺. This pattern, rather than a gradual loss of accuracy with increasing Z, is consistent with a localized numerical issue specific to the C²⁺ calculation.

Figure 4. Ground-state energies of the Be-like sequence (Be, B⁺, C²⁺) as a function of nuclear charge Z: present work compared with variational Monte Carlo and Hartree-Fock values [33] and with the Hy-CI reference energies [2] [3] (Table 2).

Figure 5. Absolute deviation |ΔE| of the present ground-state energies from the Hy-CI reference values [2] [3], on a logarithmic scale. The deviation for C²⁺ is roughly two to three orders of magnitude larger than for Be and B⁺ (Sections. 4.3-4.5).

4.4. Evolution of the Variational Parameters along the Sequence

Figure 6 shows the optimized values of α, γ and a as a function of Z, based on the values already listed in Table 1. The core exponent α and the correlation parameter a both increase smoothly and monotonically across the sequence, a behavior consistent with the expected contraction of the electronic charge distribution as the nuclear charge increases and with a progressively more important role of the explicitly correlated term as Z increases. The valence exponent γ, however, shows a qualitatively different pattern: it is essentially frozen between Be and B+ and then increases abruptly by a factor of about three for C²⁺. We note that this marked change in γ coincides with the large deviation of the C²⁺ energy from the Hy-CI reference value discussed in Section 4.2-4.3. We stress that this is an observed coincidence in the data, not a demonstrated causal relationship: establishing whether the anomalous γ value is a symptom of an underlying optimization or basis-set problem, or an independent artefact, would require further calculations: for example, an independent optimization of C²⁺ using different starting parameters or a larger basis, that lie outside the scope of the present work.

In Figure 6, the smooth increase of a reflects the phase parameter used in the complex factors and is not assigned a purely correlation-related physical meaning. The abrupt change in γ for C²⁺ is retained as an observed numerical feature; the present results do not establish a causal relation between this change and the energy discrepancy.

Figure 6. Evolution of the optimized nonlinear variational parameters $\alpha =\beta$ , $\gamma =\delta$ and a with nuclear charge Z along the Be-like sequence (Table 1). The sharp increase in γ between B⁺ and C²⁺ coincides with the large energy deviation noted for C²⁺ (Secs. 4.4-4.5).

4.5. On the Origin of the C²⁺ Discrepancy

Several practical factors could, in principle, contribute to the discrepancy identified for C²⁺. The constraint $\alpha =\beta$ and $\gamma =\delta$ adopted for computational convenience (Section 2.3) removes some of the variational freedom that would otherwise allow the two core and two valence electrons to relax independently; if this constraint interacts unfavorably with the optimization landscape at $Z=6$ , it could produce a spuriously poor local minimum without a corresponding failure at $Z=4$ or $Z=5$ . Similarly, the basis-size parameter Ω used in this work, while adequate for Be and B⁺, was not independently re-optimized for C²⁺, and a systematic basis-size study: of the kind routinely used to establish convergence in high-accuracy Hy-CI and ECG calculations, was not carried out here. We also cannot exclude a transcription or bookkeeping error specific to the $Z=6$ case, given that the associated variational parameter γ changes so much more sharply than for the other two systems. Distinguishing between these possibilities is beyond the scope of the present study, and we deliberately refrain from asserting any specific cause.

What can be stated on the basis of the present data is the following. The Hylleraas-correlated, complex-scaling variational formulation described in Section 2 agrees closely with the Hy-CI reference ground-state energy for neutral beryllium, shows reasonable agreement for B⁺, and exhibits a large, currently unexplained deviation for C²⁺. We regard this as a genuine limitation of the present implementation at the higher end of the isoelectronic sequence considered here, and we report it explicitly rather than limit the comparison to the two systems for which agreement is favorable, since the ion C²⁺ was investigated under the same conditions as Be and B⁺ and omitting it selectively would misrepresent the actual performance of the method.

4.6. Complex Scaling and the Interpretation of the Imaginary Energy Component

As summarized in Section 2.2, the complex-scaling transformation $r\to r{\text{e}}^{\left\{\text{i}\theta \right\}}$ is applied throughout this work to the same (1s²2s²) ¹S configuration that is treated, in Sections. 4.1-4.5, as the ground state of each system. As noted in Section 2.4, $H\left(\theta \right)$ is non-Hermitian for $\theta \ne 0$ , so its expectation value is not subject to the same reality and positivity constraints as a Hermitian variational calculation; within the present implementation, this expectation value $\langle \Psi |H\left(\theta \right)|\Psi \rangle /\langle \Psi |\Psi \rangle$ carries a nonzero imaginary part for all three systems. Taken at face value through Equation (3), it would correspond to a resonance half-width of several atomic units, comparable in magnitude to a sizable fraction of the total binding energy itself. Such a width is not physically defensible for a state that is, by construction, the nondegenerate ground state of a bound four-electron system, which must have $\Gamma =0$ independently of $\theta$ once convergence has been reached. We therefore do not present these quantities as physical resonance widths, and they are not included in Table 1, Table 2 or in any figure. The present data do not allow the origin of this nonzero imaginary component to be established; a systematic analysis of θ-dependence and basis-size convergence, identified in Section 5 as a priority for future work, would be required to clarify the issue. What the complex-scaling formalism does provide, and what is retained here, is the computational framework, the rotated Hamiltonian of Equation (2), within which the Hylleraas-correlated ground-state energies of Sections. 4.1-4.5 were obtained; its use for genuine resonance determination in this four-electron context remains an open problem rather than a demonstrated capability.

4.7. Computational Performance

Each system required approximately 15 to 30 minutes of computation within Maxima 5.21.1 for a single set of nonlinear parameters and basis size, reflecting the complexity of the twelve-dimensional analytic integrals involved. This cost limits the systematic basis-size and θ-convergence studies identified in Sections. 4.5 and 4.6 as necessary next steps, and motivates implementing the analytic matrix elements in a compiled numerical framework for future extensions of this work.

5. Conclusion and Perspectives

The optimized nonlinear parameters satisfy $\alpha =\beta$ and $\gamma =\delta$ throughout the calculations, reflecting the equivalent treatment of the two electrons within each coupled pair. We have examined the performance of a Hylleraas-correlated, complex-scaling variational formulation for the ground-state energies of the four-electron beryllium isoelectronic sequence, Be, B⁺ and C²⁺. The trial wave function explicitly incorporates all six interelectronic distances of the four-electron problem, and the associated matrix elements were derived analytically and evaluated numerically within Maxima. Compared against high-accuracy Hylleraas-configuration-interaction benchmark energies, the approach reproduces the ground-state energy of neutral beryllium to within 8 × 10⁻⁵ a.u. and that of B⁺ to within 1.1 × 10⁻³ a.u., while showing a substantial, currently unexplained deviation of about 0.41 a.u. for C²⁺ that coincides with an abrupt change in one of the optimized nonlinear parameters. We have reported this limitation explicitly, together with the observation, not itself resolved here, that the complex-scaling procedure, applied to what is a nominally bound ground state, currently yields nonzero imaginary energy components that cannot be defended as physical resonance widths and are accordingly excluded from the quantitative results of this work.

Taken together, these results indicate that the present formulation provides a compact description of the lower members of the Be-like sequence and also delineates the numerical studies needed to establish its range of quantitative applicability. For C²⁺, an independent basis-size analysis, several starting values of the nonlinear parameters, and relaxation of the computational constraints would help determine the origin of the larger deviation. A complementary study of the θ-dependence of the complex-scaled eigenvalues is likewise required before the formulation is applied to genuinely autoionizing four-electron configurations such as (1s²2s nℓ) states of Be and its isoelectronic ions, for which the real and imaginary parts of a stabilized complex eigenvalue would represent the resonance position and width. Implementing the analytical matrix elements in a more efficient numerical environment would facilitate these convergence and stability studies by reducing the computational cost per data point.

The numerical results are consequently interpreted with emphasis on internal consistency and comparison with independent benchmarks rather than as claims of high-accuracy θ = 0 ground-state energies. In particular, the close agreement obtained for Be is encouraging, but a definitive variational validation requires an independently defined Hermitian θ = 0 calculation satisfying the Ritz bound.

List of Symbols and Abbreviations

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References