1 Introduction↩︎

The lithium atom represents the simplest multi-electron atomic system beyond helium that eludes an exact analytical solution. Unlike the hydrogen atom, for which the Schrödinger equation yields exact closed-form solutions, the presence of electron-electron repulsion in lithium introduces dynamic correlation effects that necessitate approximation methods. Consequently, techniques such as perturbation theory and the variational method serve as indispensable tools in atomic quantum mechanics.

In this study, both approaches are systematically employed and contrasted with each other, as well as with reference benchmarks. This comparative analysis serves a pedagogical purpose by illustrating the limitations of the independent-particle approximation and highlighting the quantitative and conceptual significance of electron correlations in multi-electron systems. Throughout this work, time-independent perturbation theory and the variational method are applied within the framework of the Born-Oppenheimer approximation, where the nucleus is treated as stationary. The zeroth-order perturbative baseline treats the electrons as independent particles moving within the nuclear Coulomb potential, while inter-electronic repulsion terms are subsequently introduced as perturbative corrections. The variational method on the other hand utilizes a trial function with two varied parameters, each representing the shielding effect on the potential energy in each orbital, and than minimizing the expectation value of the Hamiltonian in order to estimate the ground state energy.

The ground-state electronic configuration of the lithium atom is given by \[1s^2 2s^1\] which consists of two electrons occupying the inner \(1s\) shell and a single valence electron in the \(2s\) orbital.

2 Non-Relativistic Hamiltonian of the Lithium Atom↩︎

The time-independent Schrödinger equation is expressed as [1]: \[\hat{H}\Psi = E\Psi\]

For a lithium atom with three electrons, the exact non-relativistic Hamiltonian (\(\hat{H}\)), explicitly expanding all individual kinetic and potential energy components, is given by [1]: \[\hat{H} = \sum_{i=1}^{3} \left[ -\frac{\hbar^2}{2m}\nabla_i^2 - \frac{Ze^2}{4\pi\varepsilon_0 r_i} \right] + \sum_{i<j}^{3} \frac{e^2}{4\pi\varepsilon_0 r_{ij}}\] where \(Z = 3\) represents the atomic number of lithium, \(r_i\) denotes the distance between the nucleus and the \(i\)-th electron, and \(r_{ij} = |\mathbf{r}_i - \mathbf{r}_j|\) represents the relative inter-electronic distance between electrons \(i\) and \(j\). Expanding \(\hat{H}\) term-by-term into individual electronic coordinates yields: \[\begin{align} \hat{H} &= \left[ -\frac{\hbar^2}{2m}\nabla_1^2 - \frac{3e^2}{4\pi\varepsilon_0 r_1} \right] + \left[ -\frac{\hbar^2}{2m}\nabla_2^2 - \frac{3e^2}{4\pi\varepsilon_0 r_2} \right] \\ &\quad + \left[ -\frac{\hbar^2}{2m}\nabla_3^2 - \frac{3e^2}{4\pi\varepsilon_0 r_3} \right] + \frac{e^2}{4\pi\varepsilon_0 r_{12}} + \frac{e^2}{4\pi\varepsilon_0 r_{13}} + \frac{e^2}{4\pi\varepsilon_0 r_{23}} \end{align}\]

The coupled nature of the inter-electronic repulsion terms (\(1/r_{ij}\)) mathematically prevents the separation of spatial coordinates, rendering an exact analytical solution impossible. To systematically resolve this challenge, we begin our theoretical treatment by applying Rayleigh-Schrödinger perturbation theory, isolating these problematic repulsion terms as a structural perturbation acting upon an exactly solvable independent-particle system.

3 Perturbation Theory Approximation↩︎

To systematically evaluate the electronic structure under a perturbative framework, the total non-relativistic Hamiltonian is partitioned into an unperturbed multi-particle base operator (\(\hat{H}_0\)) and a perturbation operator (\(\hat{H}^{(1)}\)) that accounts for the mutual inter-electronic repulsive interactions [2]: \[\hat{H} = \hat{H}_0 + \hat{H}^{(1)}\]

3.1 Zeroth-Order Approximation↩︎

In the zeroth-order approximation, the electron-electron Coulomb repulsion is entirely neglected. Consequently, the unperturbed Hamiltonian (\(\hat{H}_0\)) decomposes into a direct sum of independent, single-particle hydrogen-like operators [1]: \[\hat{H}_0 = \hat{h}_1 + \hat{h}_2 + \hat{h}_3 = \sum_{i=1}^{3} \left[ -\frac{\hbar^2}{2m_e}\nabla_i^2 - \frac{Ze^2}{4\pi\varepsilon_0 r_i} \right]\] The corresponding unperturbed state satisfies the independent-particle Schrödinger equation [1]: \[\hat{H}_0 \Psi^{(0)} = E^{(0)} \Psi^{(0)}\]

3.1.1 Zeroth-Order Wave Function Construction↩︎

Because the unperturbed Hamiltonian contains no electronic cross-terms, the multi-particle spatial solution can be modeled as a simple product of independent, hydrogen-like single-particle orbitals [1]: \[\psi^{(0)}(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3) = \phi_{1s}(\mathbf{r}_1)\phi_{1s}(\mathbf{r}_2)\phi_{2s}(\mathbf{r}_3)\] These spatial functions are designated as hydrogen-like because each electron effectively experiences an unshielded nuclear potential energy of \(-Ze^2 / (4\pi\varepsilon_0 r_i)\) with \(Z=3\), rendering the individual equations mathematically analogous to a single-electron system.

In general, any atomic spatial orbital \(\phi\) can be decomposed into a product of a radial wave function (\(R_{n,l}\)) and a spherical harmonic (\(Y_l^m\)) representing the angular component [2]: \[\phi_{n,l,m}(\mathbf{r}) = R_{n,l}(r) Y_l^m(\theta, \varphi)\]

For the ground-state configuration of lithium, the orbitals of interest are the \(1s\) and \(2s\) states. For \(s\)-orbitals (where the orbital angular momentum quantum number \(l=0\)), the wave function exhibits no angular dependence, and the spherical harmonic reduces to the constant \(Y_0^0 = 1/\sqrt{4\pi}\). Multiplying this angular constant by the corresponding radial functions yields the following explicit spatial profiles evaluated in perturbation theory [2]: \[\phi_{1s}(\mathbf{r}) = \frac{1}{\sqrt{\pi}} \left( \frac{Z}{a_0} \right)^{3/2} e^{-Z r / a_0}\] \[\phi_{2s}(\mathbf{r}) = \frac{1}{4\sqrt{2\pi}} \left( \frac{Z}{a_0} \right)^{3/2} \left( 2 - \frac{Z r}{a_0} \right) e^{-Z r / (2a_0)}\]

To construct a complete, physically acceptable non-relativistic quantum state that satisfies anti-symmetrization, the spin degrees of freedom must be integrated to form individual spin-orbitals, defined as \(u_k(\mathbf{x}_i) = \phi_k(\mathbf{r}_i)\sigma_k(s_i)\). Here, \(\mathbf{x}_i = (\mathbf{r}_i, s_i)\) denotes the combined spatial and spin coordinates, and the spin functions are given by [1]: \[\sigma_k(s_i) = \begin{cases} \alpha(s_i), & \text{if } m_s = +\frac{1}{2} \\ \beta(s_i), & \text{if } m_s = -\frac{1}{2} \end{cases}\]

By designating the spin-up and spin-down projection states as \(\alpha\) and \(\beta\) respectively, the active spin-orbitals for the \(1s^2 2s^1\) ground-state configuration are formulated as [2]: \[u_1(\mathbf{x}_i) = \phi_{1s}(\mathbf{r}_i)\alpha(s_i), \quad u_2(\mathbf{x}_i) = \phi_{1s}(\mathbf{r}_i)\beta(s_i), \quad u_3(\mathbf{x}_i) = \phi_{2s}(\mathbf{r}_i)\alpha(s_i)\] According to the Pauli Exclusion Principle, no two indistinguishable fermions can occupy the same quantum state simultaneously. Therefore, the two electrons residing within the identical spatial \(1s\) orbital must possess opposing spin projections (\(m_s = +1/2\) and \(m_s = -1/2\)). The single valence electron occupying the \(2s\) orbital can formally assume either spin projection without altering the scalar energy expectation values; hence, it is chosen arbitrarily as a spin-up (\(\alpha\)) state [1]. In matrix representation, these spinors are explicitly defined as column vectors [1]: \[\alpha = \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \quad \beta = \begin{pmatrix} 0 \\ 1 \end{pmatrix}\]

For many-body perturbation integrations, these spin functions must satisfy the standard orthonormality relations. The normalization conditions for both projection states are expressed using formal spin-space integration as [2]: \[\int \alpha^*(s)\alpha(s)\,ds = 1, \quad \int \beta^*(s)\beta(s)\,ds = 1\] Furthermore, the orthogonality condition, which enforces the physical requirement that a single electron cannot simultaneously occupy opposing spin projections, is defined by: \[\int \alpha^*(s)\beta(s)\,ds = 0\]

Because electrons are identical fermions, the overall multi-particle wave function must be totally antisymmetric under the permutation of any two particle indices (\(P_{ij}\Psi^{(0)} = -\Psi^{(0)}\)). This global symmetry requirement is rigorously satisfied by constructing a \(3 \times 3\) Slater determinant framework (\(\mathcal{M}\)) [2]: \[\Psi^{(0)} = \frac{1}{\sqrt{3!}} \det(\mathcal{M}) = \frac{1}{\sqrt{6}} \begin{vmatrix} u_1(\mathbf{x}_1) & u_2(\mathbf{x}_1) & u_3(\mathbf{x}_1) \\ u_1(\mathbf{x}_2) & u_2(\mathbf{x}_2) & u_3(\mathbf{x}_2) \\ u_1(\mathbf{x}_3) & u_2(\mathbf{x}_3) & u_3(\mathbf{x}_3) \end{vmatrix}\] where the arguments \((1, 2, 3)\) are explicitly written as \(\mathbf{x}_i\) to represent the combined spatial and spin coordinates of the \(i\)-th electron. Expanding this determinant explicitly yields the complete linear combination consisting of six distinct permutation terms: \[\begin{align} \Psi^{(0)} = \frac{1}{\sqrt{6}} \Big[ & u_1(\mathbf{x}_1)u_2(\mathbf{x}_2)u_3(\mathbf{x}_3) - u_1(\mathbf{x}_1)u_3(\mathbf{x}_2)u_2(\mathbf{x}_3) \\ & + u_2(\mathbf{x}_1)u_3(\mathbf{x}_2)u_1(\mathbf{x}_3) - u_2(\mathbf{x}_1)u_1(\mathbf{x}_2)u_3(\mathbf{x}_3) \\ & + u_3(\mathbf{x}_1)u_1(\mathbf{x}_2)u_2(\mathbf{x}_3) - u_3(\mathbf{x}_1)u_2(\mathbf{x}_2)u_1(\mathbf{x}_3) \Big] \end{align}\]

3.1.2 Systematic Pairwise Grouping↩︎

To optimize mathematical efficiency and fully exploit orbital orthonormality during subsequent expectation value calculations, the six expanded permutation terms are clustered pairwise into three distinct symmetric brackets, denoted as \(A\), \(B\), and \(C\): \[\Psi^{(0)} = \frac{1}{\sqrt{6}} (A + B + C)\] where these sub-component functional groups are analytically partitioned as: \[\begin{align} A &= u_1(\mathbf{x}_1)u_2(\mathbf{x}_2)u_3(\mathbf{x}_3) - u_2(\mathbf{x}_1)u_1(\mathbf{x}_2)u_3(\mathbf{x}_3) = \big[u_1(\mathbf{x}_1)u_2(\mathbf{x}_2) - u_2(\mathbf{x}_1)u_1(\mathbf{x}_2)\big] u_3(\mathbf{x}_3) \\ B &= u_3(\mathbf{x}_1)u_1(\mathbf{x}_2)u_3(\mathbf{x}_3) - u_1(\mathbf{x}_1)u_3(\mathbf{x}_2)u_2(\mathbf{x}_3) = \big[u_3(\mathbf{x}_1)u_1(\mathbf{x}_2) - u_1(\mathbf{x}_1)u_3(\mathbf{x}_2)\big] u_2(\mathbf{x}_3) \\ C &= u_2(\mathbf{x}_1)u_3(\mathbf{x}_2)u_1(\mathbf{x}_3) - u_3(\mathbf{x}_1)u_2(\mathbf{x}_2)u_1(\mathbf{x}_3) = \big[u_2(\mathbf{x}_1)u_3(\mathbf{x}_2) - u_3(\mathbf{x}_1)u_2(\mathbf{x}_2)\big] u_1(\mathbf{x}_3) \end{align}\]

3.1.3 Zeroth-Order Energy Baseline↩︎

The unperturbed single-particle energy eigenvalues for a hydrogen-like system are determined via the analytical solution to the radial Schrödinger equation [2]: \[E_n = -\left( \frac{m_e e^4}{2(4\pi\varepsilon_0)^2 \hbar^2} \right) \frac{Z^2}{n^2} \approx -13.6057\,\text{eV} \cdot \frac{Z^2}{n^2}\]

For the lithium nucleus (\(Z=3\)), the constituent physical constants are defined as follows: the electron rest mass \(m_e = 0.5109989\,\mathrm{MeV}/c^2\), the elementary charge \(e = 1.602176 \times 10^{-19}\,\mathrm{C}\), the vacuum permittivity \(\varepsilon_0 = 8.854188 \times 10^{-12}\,\mathrm{F/m}\), and the reduced Planck constant \(\hbar = 6.582120 \times 10^{-16}\,\mathrm{eV\cdot s}\). The pre-factor encapsulates the fundamental Rydberg energy unit (\(-13.6057\,\text{eV}\)), yielding the individual unshielded orbital energy levels directly as: \[E_{1s} = -13.6057 \left( \frac{3^2}{1^2} \right) \approx -122.45\,\text{eV}\] \[E_{2s} = -13.6057 \left( \frac{3^2}{2^2} \right) \approx -30.61\,\text{eV}\]

Since the unperturbed total Hamiltonian operator is a direct linear sum of non-interacting components, the total zeroth-order energy baseline is determined by the cumulative energies of the occupied orbital states (\(2E_{1s} + E_{2s}\)): \[E^{(0)} = 2(-122.45\,\text{eV}) + (-30.61\,\text{eV}) = -275.51\,\text{eV}\]

3.2 First-Order Perturbation Theory↩︎

The perturbation Hamiltonian accounting for the cumulative inter-electronic Coulomb repulsion is expressed as [2]: \[\hat{H}^{(1)} = \hat{H}_{12}^{(1)} + \hat{H}_{13}^{(1)} + \hat{H}_{23}^{(1)} = \frac{e^2}{4\pi\varepsilon_0 r_{12}} + \frac{e^2}{4\pi\varepsilon_0 r_{13}} + \frac{e^2}{4\pi\varepsilon_0 r_{23}}\] where the indices represent the pairwise interactions formed through the circular permutation of the three electrons, with \(r_{ij} = |\mathbf{r}_i - \mathbf{r}_j|\). The corresponding first-order energy correction is evaluated as the expectation value of this operator over the unperturbed state [1]: \[E^{(1)} = \langle \Psi^{(0)} | \hat{H}^{(1)} | \Psi^{(0)} \rangle\]

Since the unperturbed wave function \(\Psi^{(0)}\) is completely antisymmetric under coordinate exchange and the individual electrons are fundamentally indistinguishable, the inner product over each pairwise repulsion operator yields identical mathematical results [2]: \[\langle \Psi^{(0)} | \hat{H}_{12}^{(1)} | \Psi^{(0)} \rangle = \langle \Psi^{(0)} | \hat{H}_{13}^{(1)} | \Psi^{(0)} \rangle = \langle \Psi^{(0)} | \hat{H}_{23}^{(1)} | \Psi^{(0)} \rangle\]

This permutational symmetry simplifies the total expectation value down to a single interaction channel scaled by a combinatorial factor of three: \[E^{(1)} = 3 \langle \Psi^{(0)} | \hat{H}_{12}^{(1)} | \Psi^{(0)} \rangle\]

Expanding this expression using the previously established pairwise grouped brackets \(A\), \(B\), and \(C\) (Equation 18) yields: \[\begin{align} E^{(1)} &= 3 \cdot \left( \frac{1}{6} \right) \langle A + B + C | \hat{H}_{12}^{(1)} | A + B + C \rangle \\ &= \frac{1}{2} \Big[ \langle A | \hat{H}_{12}^{(1)} | A \rangle + \langle B | \hat{H}_{12}^{(1)} | B \rangle + \langle C | \hat{H}_{12}^{(1)} | C \rangle \\ &\quad + \langle A | \hat{H}_{12}^{(1)} | B \rangle + \langle B | \hat{H}_{12}^{(1)} | A \rangle + \langle A | \hat{H}_{12}^{(1)} | C \rangle \\ &\quad + \langle C | \hat{H}_{12}^{(1)} | A \rangle + \langle B | \hat{H}_{12}^{(1)} | C \rangle + \langle C | \hat{H}_{12}^{(1)} | B \rangle \Big] \end{align}\]

3.2.1 Orthogonality Reduction to Coulomb and Exchange Integrals↩︎

The matrix elements within Eq.(28) can be drastically simplified by integrating over the coordinates of the non-interacting unperturbed spectator electron. Leveraging the spatial and spin orthonormality of the single-particle spin-orbitals (\(\langle u_i | u_j \rangle = \delta_{ij}\)) [1], only the three diagonal matrix elements (\(AA\), \(BB\), and \(CC\)) provide non-vanishing contributions, which are evaluated as follows:

i. Evaluation of \(\langle A | \hat{H}_{12}^{(1)} | A \rangle\): In this configuration, electron 3 acts as the spectator particle residing in state \(u_3\). Utilizing the normalization condition \(\langle u_3(\mathbf{x}_3)|u_3(\mathbf{x}_3) \rangle = 1\), the integration reduces to the spatial and spin coordinates of the two active interacting electrons (\(\mathbf{x}_1\) and \(\mathbf{x}_2\)) [2]: \[\begin{align} \langle A|\hat{H}_{12}^{(1)}|A \rangle &= \iint [u_1(\mathbf{x}_1)u_2(\mathbf{x}_2) - u_2(\mathbf{x}_1)u_1(\mathbf{x}_2)]^* \hat{H}_{12}^{(1)} [u_1(\mathbf{x}_1)u_2(\mathbf{x}_2) - u_2(\mathbf{x}_1)u_1(\mathbf{x}_2)] \,d\mathbf{x}_1 \,d\mathbf{x}_2 \\ &= 2 \iint |\phi_{1s}(\mathbf{r}_1)|^2 \frac{e^2}{4\pi\varepsilon_0 |\mathbf{r}_1 - \mathbf{r}_2|} |\phi_{1s}(\mathbf{r}_2)|^2 \,d^3\mathbf{r}_1 \,d^3\mathbf{r}_2 \\ &= 2J_{1s,1s} \end{align}\] where \(J_{1s,1s}\) denotes the direct Coulomb integral for the two interacting electrons within the \(1s\) core shell.

ii. Evaluation of \(\langle B | \hat{H}_{12}^{(1)} | B \rangle\): Here, electron 3 occupies the spectator state \(u_2\). Performing the integration over the remaining active spin-orbital coordinates yields [2]: \[\begin{align} \langle B|\hat{H}_{12}^{(1)}|B \rangle &= \iint [u_3(\mathbf{x}_1)u_1(\mathbf{x}_2) - u_1(\mathbf{x}_1)u_3(\mathbf{x}_2)]^* \hat{H}_{12}^{(1)} [u_3(\mathbf{x}_1)u_1(\mathbf{x}_2) - u_1(\mathbf{x}_1)u_3(\mathbf{x}_2)] \,d\mathbf{x}_1 \,d\mathbf{x}_2 \\ &= 2 \iint |\phi_{1s}(\mathbf{r}_1)|^2 \frac{e^2}{4\pi\varepsilon_0|\mathbf{r}_1 - \mathbf{r}_2|} |\phi_{2s}(\mathbf{r}_2)|^2 \,d^3\mathbf{r}_1 \,d^3\mathbf{r}_2 \\ &\quad - 2 \iint \phi_{1s}^*(\mathbf{r}_1) \phi_{2s}^*(\mathbf{r}_2) \frac{e^2}{4\pi\varepsilon_0|\mathbf{r}_1 - \mathbf{r}_2|} \phi_{1s}(\mathbf{r}_2) \phi_{2s}(\mathbf{r}_1) \,d^3\mathbf{r}_1 \,d^3\mathbf{r}_2 \\ &= 2J_{1s,2s} - 2K_{1s,2s} \end{align}\] where \(J_{1s,2s}\) and \(K_{1s,2s}\) represent the direct Coulomb and quantum exchange integrals, respectively, between the \(1s\) core and \(2s\) valence states.

iii. Evaluation of \(\langle C | \hat{H}_{12}^{(1)} | C \rangle\): In this bracket, electron 3 is the spectator particle in state \(u_1\). Because \(\langle u_1(\mathbf{x}_3)|u_1(\mathbf{x}_3) \rangle = 1\), the matrix element evaluates to [2]: \[\begin{align} \langle C|\hat{H}_{12}^{(1)}|C \rangle &= \iint [u_2(\mathbf{x}_1)u_3(\mathbf{x}_2) - u_3(\mathbf{x}_1)u_2(\mathbf{x}_2)]^* \hat{H}_{12}^{(1)} [u_2(\mathbf{x}_1)u_3(\mathbf{x}_2) - u_3(\mathbf{x}_1)u_2(\mathbf{x}_2)] \,d\mathbf{x}_1 \,d\mathbf{x}_2 \\ &= 2 \iint |\phi_{1s}(\mathbf{r}_1)|^2 \frac{e^2}{4\pi\varepsilon_0|\mathbf{r}_1 - \mathbf{r}_2|} |\phi_{2s}(\mathbf{r}_2)|^2 \,d^3\mathbf{r}_1 \,d^3\mathbf{r}_2 \\ &= 2J_{1s,2s} \end{align}\] Notably, the expected exchange term in the \(C\)-bracket vanishes due to explicit spin orthogonality, as the corresponding spin integral contains the factors \(\langle \alpha|\beta\rangle = 0\) and \(\langle \beta|\alpha\rangle = 0\).

All cross-terms (such as \(\langle A | \hat{H}_{12}^{(1)} | B \rangle\)) vanish entirely because the spectator electron maps to mutually orthogonal spin-orbitals (e.g., \(\langle u_3 | u_2 \rangle = 0\) or \(\langle u_3 | u_1 \rangle = 0\)). Collecting the surviving non-zero components and substituting them back into the primary expression yields: \[E^{(1)} = \frac{1}{2} \Big[ 2J_{1s,1s} + \big(2J_{1s,2s} - 2K_{1s,2s}\big) + 2J_{1s,2s} \Big] = J_{1s,1s} + 2J_{1s,2s} - K_{1s,2s}\]

3.2.2 Numerical Substitution and First-Order Energy Results↩︎

The spatial double integrals defining these static Coulomb and exchange interactions are explicitly written as [2]: \[J_{1s,1s} = \iint \phi_{1s}^*(\mathbf{r}_1)\phi_{1s}^*(\mathbf{r}_2) \left(\frac{e^2}{4\pi\varepsilon_0 r_{12}}\right) \phi_{1s}(\mathbf{r}_1)\phi_{1s}(\mathbf{r}_2) \,d^3\mathbf{r}_1 \,d^3\mathbf{r}_2\] \[J_{1s,2s} = \iint \phi_{1s}^*(\mathbf{r}_1)\phi_{2s}^*(\mathbf{r}_2) \left(\frac{e^2}{4\pi\varepsilon_0 r_{12}}\right) \phi_{1s}(\mathbf{r}_1)\phi_{2s}(\mathbf{r}_2) \,d^3\mathbf{r}_1 \,d^3\mathbf{r}_2\] \[K_{1s,2s} = \iint \phi_{1s}^*(\mathbf{r}_1)\phi_{2s}^*(\mathbf{r}_2) \left(\frac{e^2}{4\pi\varepsilon_0 r_{12}}\right) \phi_{2s}(\mathbf{r}_1)\phi_{1s}(\mathbf{r}_2) \,d^3\mathbf{r}_1 \,d^3\mathbf{r}_2\]

The exact analytical evaluations of these multi-center integrals for the lithium nucleus (\(Z=3\)) yield the following expressions (the detailed integral derivations are provided in Appendix A): \[J_{1s,1s} = \frac{5}{8} \left(\frac{Z e^2}{4\pi\varepsilon_0 a_0}\right), \quad J_{1s,2s} = \frac{17}{81} \left(\frac{Z e^2}{4\pi\varepsilon_0 a_0}\right), \quad K_{1s,2s} = \frac{16}{729} \left(\frac{Z e^2}{4\pi\varepsilon_0 a_0}\right)\]

Factoring out the core unit energy group \(\frac{e^2}{4\pi\varepsilon_0 a_0} = 2\,\text{Ry} \approx 27.2114\,\text{eV}\), the total first-order energy correction for \(Z=3\) evaluates quantitatively to: \[E^{(1)} = 3\left[\left( \frac{5}{8} \right) + 2\left(\frac{17}{81}\right) - \left( \frac{16}{729} \right)\right] \times 27.2114\,\text{eV} \approx 83.50\,\text{eV}\] which is in precise agreement with standard literature references [2].

Consequently, the total cumulative ground-state energy corrected to first order is given by: \[E_{\text{total}}^{(1)} = E^{(0)} + E^{(1)} \approx -275.51\,\text{eV} + 83.50\,\text{eV} = -192.01\,\text{eV}\] While this first-order correction significantly resolves the unphysical baseline energy, it still deviates from the reference non-relativistic reference ground-state energy of lithium (\(\sim -203.5\;\text{eV}\)) by approximately \(5.65\%\). This residual error underscores the necessity of computing second-order perturbation corrections to account for dynamic electron correlation effects. The reference value is adopted from [3] and the total ionization energy is calculated via adding up the three values (\(-5.391714996\;\text{eV} -75.6400970\;\text{eV} -122.45435913\;\text{eV} = -203.486119\,\text{eV}\)) shown in the reference, where the first ionization energy is measured experimentally, whereas the other two values are calculated theoretically. The final result is rounded herein \(-203.5\;\text{eV}\).

3.3 Second-Order Perturbation Theory↩︎

While first-order perturbation theory evaluates electron-electron repulsion strictly within an averaged-field framework, incorporating dynamical, instantaneous electronic correlation and virtual state mixing necessitates higher-order corrections [2] [4]. The second-order energy correction can be formally derived from the time-independent Schrödinger equation \(\hat{H}|\Psi\rangle = E|\Psi\rangle\) by expanding the operators and states in terms of a continuous ordering parameter \(\lambda\) [2] [5]: \[\begin{align} \hat{H} &= \hat{H}_0 + \lambda \hat{H}^{(1)} \nonumber \\ |\Psi\rangle &= |\Psi_0^{(0)}\rangle + \lambda |\Psi^{(1)}\rangle + \lambda^2 |\Psi^{(2)}\rangle + \dots \nonumber \\ E &= E_0^{(0)} + \lambda E^{(1)} + \lambda^2 E^{(2)} + \dots \end{align}\]

Substituting these series expansions into the Schrödinger equation and isolating the second-order (\(\lambda^2\)) components establishes the fundamental relation [1]: \[\hat{H}_0 |\Psi^{(2)}\rangle + \hat{H}^{(1)} |\Psi^{(1)}\rangle = E_0^{(0)} |\Psi^{(2)}\rangle + E^{(1)} |\Psi^{(1)}\rangle + E^{(2)} |\Psi_0^{(0)}\rangle\]

Projecting this relation onto the unperturbed ground-state bra \(\langle\Psi_0^{(0)}|\) and utilizing the Hermiticity of the unperturbed Hamiltonian (\(\langle\Psi_0^{(0)}|\hat{H}_0 = E_0^{(0)} \langle\Psi_0^{(0)}|\)), the unknown second-order wave function components cancel out identically. Assuming an intermediate-normalized zeroth-order ground state (\(\langle\Psi_0^{(0)}|\Psi_0^{(0)}\rangle = 1\)), the expression reduces to the compact bracket formulation [1]: \[E^{(2)} = \langle\Psi_0^{(0)}| \hat{H}^{(1)} - E^{(1)} |\Psi^{(1)}\rangle\]

To resolve this expression into an explicitly calculable sum over static unperturbed states, the first-order wave function correction \(|\Psi^{(1)}\rangle\)—which represents the explicit multi-body polarization of the doublet state under electron correlation—is expanded as a linear combination of the unperturbed excited eigenstates \(|\Psi_m^{(0)}\rangle\) [1]: \[|\Psi^{(1)}\rangle = \sum_{m\ne0} c_m |\Psi_{m}^{(0)}\rangle \quad \text{where} \quad c_m = \frac{\langle\Psi_{m}^{(0)}|\hat{H}^{(1)}|\Psi_{0}^{(0)}\rangle}{E_{0}^{(0)} - E_{m}^{(0)}}\]

Substituting this expansion back into the bracket relation yields the standard Rayleigh-Schrödinger second-order energy correction formula: \[E^{(2)} = \sum_{m \neq 0} \frac{\left| \langle \Psi_m^{(0)} | \hat{H}^{(1)} | \Psi_0^{(0)} \rangle \right|^2}{E_0^{(0)} - E_m^{(0)}}\]

At this stage, it is crucial to address the intrinsic spin degeneracy of the ground state. The unperturbed lithium configuration (\(1s^2 2s^1\)) is two-fold degenerate because the valence electron in the \(2s\) orbital can assume either a spin-up (\(m_s = +1/2\)) or a spin-down (\(m_s = -1/2\)) projection, both yielding the identical unperturbed energy \(E_0^{(0)}\). The presence of degeneracy formally requires the framework of Degenerate Perturbation Theory (DPT), which dictates the construction and diagonalization of the perturbation matrix within this degenerate subspace [2]: \[\mathbf{H}^{(1)} = \begin{pmatrix} \langle \alpha | \hat{H}^{(1)} | \alpha \rangle & \langle \alpha | \hat{H}^{(1)} | \beta \rangle \\ \langle \beta | \hat{H}^{(1)} | \alpha \rangle & \langle \beta | \hat{H}^{(1)} | \beta \rangle \end{pmatrix}\]

However, the physical perturbation operator—the inter-electronic Coulomb repulsion \(\hat{H}^{(1)} = \sum \frac{e^2}{4\pi\varepsilon_0 r_{ij}}\)—is purely spatial and contains no spin-dependent components. Because this operator commutes with the total spin operators, the quantum mechanical integrals separate into independent spatial and spin products. For the off-diagonal elements [in eq. 44], the inner product of the spin states strictly enforces the spin orthogonality condition: \[\langle \alpha | \hat{H}^{(1)} | \beta \rangle \propto \int \alpha^*(s) \beta(s) \,ds = 0\]

Consequently, the cross-subspace coupling terms vanish entirely, demonstrating that the perturbation matrix is inherently diagonal. Furthermore, since the spatial integrals for the diagonal elements are physically identical, no lifting of the degeneracy occurs. Because the perturbation does not couple or split these degenerate states, the DPT matrix equation mathematically collapses directly into the standard Non-Degenerate Perturbation Theory (NDPT) framework [2]. This spatial symmetry justifies the direct application of the non-degenerate second-order summation formula across individual orbital channels.

3.3.1 Excited Slater Matrix and Integral Reduction↩︎

To evaluate the matrix elements for the excited states, we insert an unpopulated higher virtual single-particle orbital, such as \(u_4 = \phi_{3s}\alpha\), into the system configuration space to construct the excited Slater matrix (\(\mathcal{M}_m\)): \[\Psi_m^{(0)} = \frac{1}{\sqrt{6}} \det(\mathcal{M}_m) = \frac{1}{\sqrt{6}} \begin{vmatrix} u_1(\mathbf{x}_1) & u_2(\mathbf{x}_1) & u_4(\mathbf{x}_1) \\ u_1(\mathbf{x}_2) & u_2(\mathbf{x}_2) & u_4(\mathbf{x}_2) \\ u_1(\mathbf{x}_3) & u_2(\mathbf{x}_3) & u_4(\mathbf{x}_3) \end{vmatrix}\]

Expanding this determinant explicitly along the first row yields the linear combination of six distinct permutation terms: \[\begin{align} \Psi_{m}^{(0)} = \frac{1}{\sqrt{6}} \Big[ &u_1(\mathbf{x}_1)u_2(\mathbf{x}_2)u_4(\mathbf{x}_3) - u_1(\mathbf{x}_1)u_4(\mathbf{x}_2)u_2(\mathbf{x}_3) \\ &- u_2(\mathbf{x}_1)u_1(\mathbf{x}_2)u_4(\mathbf{x}_3) + u_2(\mathbf{x}_1)u_4(\mathbf{x}_2)u_1(\mathbf{x}_3) \\ &+ u_4(\mathbf{x}_1)u_1(\mathbf{x}_2)u_2(\mathbf{x}_3) - u_4(\mathbf{x}_1)u_2(\mathbf{x}_2)u_1(\mathbf{x}_3) \Big] \end{align}\]

By systematically rearranging these terms, the active spatial and spin pairs can be isolated. These six expanded permutation terms are clustered pairwise into three distinct grouped brackets denoted as \(A_m\), \(B_m\), and \(C_m\): \[\Psi_m^{(0)} = \frac{1}{\sqrt{6}} (A_m + B_m + C_m)\] where the excited multi-particle sub-components are defined as: \[\begin{align} A_m &= \big[u_1(\mathbf{x}_1)u_2(\mathbf{x}_2) - u_2(\mathbf{x}_1)u_1(\mathbf{x}_2)\big] u_4(\mathbf{x}_3) \\ B_m &= \big[u_4(\mathbf{x}_1)u_1(\mathbf{x}_2) - u_1(\mathbf{x}_1)u_4(\mathbf{x}_2)\big] u_2(\mathbf{x}_3) \\ C_m &= \big[u_2(\mathbf{x}_1)u_4(\mathbf{x}_2) - u_4(\mathbf{x}_1)u_2(\mathbf{x}_2)\big] u_1(\mathbf{x}_3) \end{align}\]

Projecting the pairwise perturbation operator \(\hat{H}_{12}^{(1)}\) across these excited configurations simplifies the multi-body channels directly into specific spatial transition integrals over the joint spatial-spin coordinates (\(d\mathbf{x}_1 \,d\mathbf{x}_2\)): \[\begin{align} \langle A_m | \hat{H}_{12}^{(1)} | A \rangle = 0 \\[0.3cm]\langle B_m | \hat{H}_{12}^{(1)} | B \rangle &= \iint [u_4(\mathbf{x}_1)u_1(\mathbf{x}_2) - u_1(\mathbf{x}_1)u_4(\mathbf{x}_2)]^* \hat{H}_{12}^{(1)} [u_3(\mathbf{x}_1)u_1(\mathbf{x}_2) - u_1(\mathbf{x}_1)u_3(\mathbf{x}_2)] \,d\mathbf{x}_1 \,d\mathbf{x}_2 \\ &= 2J_{3s,1s}^{\prime} - 2K_{3s,1s}^{\prime} \\[0.3cm] \langle C_m | \hat{H}_{12}^{(1)} | C \rangle &= \iint [u_2(\mathbf{x}_1)u_4(\mathbf{x}_2) - u_4(\mathbf{x}_1)u_2(\mathbf{x}_2)]^* \hat{H}_{12}^{(1)} [u_2(\mathbf{x}_1)u_3(\mathbf{x}_2) - u_3(\mathbf{x}_1)u_2(\mathbf{x}_2)] \,d\mathbf{x}_1 \,d\mathbf{x}_2 \\ &= 2J_{3s,1s}^{\prime}\end{align}\]

Consequently, the total aggregated transition matrix element (\(M\)) reduces to the following linear combination: \[M = \langle \Psi_m^{(0)} | \hat{H}_{12}^{(1)} | \Psi_0^{(0)} \rangle = 2J'_{3s,1s} - K'_{3s,1s}\]

The defining spatial double integrals for these virtual transition states are given by [6]: \[J'_{3s,1s} = \iint \phi_{3s}^*(\mathbf{r}_1)\phi_{1s}^*(\mathbf{r}_2) \left(\frac{e^2}{4\pi\varepsilon_0 r_{12}}\right) \phi_{2s}(\mathbf{r}_1)\phi_{1s}(\mathbf{r}_2) \,d^3\mathbf{r}_1 \,d^3\mathbf{r}_2\] \[K'_{3s,1s} = \iint \phi_{3s}^*(\mathbf{r}_1)\phi_{1s}^*(\mathbf{r}_2) \left(\frac{e^2}{4\pi\varepsilon_0 r_{12}}\right) \phi_{1s}(\mathbf{r}_1)\phi_{2s}(\mathbf{r}_2) \,d^3\mathbf{r}_1 \,d^3\mathbf{r}_2\]

To evaluate the energy denominator [in eq. 43], the unperturbed zeroth-order energy of the excited \(1s^2 3s^1\) configuration (\(E_{3s}^{(0)}\)) must be determined. Within the independent-particle model, the energy of a hydrogen-like state is dictated by: \[E_n = -13.6057\,\text{eV} \cdot \frac{Z^2}{n^2}\]

For the lithium atom (\(Z=3\)), the individual single-particle energy eigenvalues for the \(1s\) and \(3s\) states evaluate to: \[\begin{align} E_{1s} &= -13.6057 \cdot \frac{3^2}{1^2} \approx -122.45\,\text{eV} \\ E_{3s} &= -13.6057 \cdot \frac{3^2}{3^2} \approx -13.61\,\text{eV} \end{align}\]

Since the excited configuration consists of two core electrons in the \(1s\) orbital and a single valence electron promoted to the \(3s\) orbital, the total unperturbed energy is the linear sum of these occupied states: \[E_{3s}^{(0)} = 2E_{1s} + E_{3s} \approx 2(-122.45\,\text{eV}) + (-13.61\,\text{eV}) = -258.51\,\text{eV}\]

The resulting energy denominator required for the second-order perturbation expansion is given by: \[E_0^{(0)} - E_{3s}^{(0)} \approx -275.51\,\text{eV} - (-258.51\,\text{eV}) = -17.00\,\text{eV}\] Substituting the evaluated transition matrix element and the energy denominator into the second-order expansion yields the specific single-excitation correction: \[E_{3s}^{(2)} = \frac{(2J'_{3s,1s} - K'_{3s,1s})^2}{E_0^{(0)} - E_{3s}^{(0)}}\]

The analytical numerical integration of these spatial double transition integrals yields the following values: \[J'_{3s,1s} \approx 4.12\,\text{eV}, \quad K'_{3s,1s} \approx 0.92\,\text{eV}\] Substituting these transition values into the second-order energy expression results quantitatively in: \[E_{3s}^{(2)} = \frac{(2(4.12) - 0.92)^2}{-17.00} \approx -3.16\,\text{eV}\]

As demonstrated, the analytical evaluation of virtual state mixings involves highly intricate multi-center integrations that are practically restricted to spherically symmetric channels such as the \(3s\) state. While physical electron correlation is heavily driven by excitations into higher angular momentum states (such as \(p\) and \(d\) channels), individual single-electron transitions to these states are strictly forbidden by selection rules due to the spherical symmetry (\(L=0\)) of the ground state and the scalar nature of the Coulomb operator. Consequently, the remaining accessible single-excitation channels are confined strictly to higher \(s\)-orbitals (e.g., \(4s\), \(5s\), and \(6s\)).

To capture these successive single-particle correlation effects without encountering analytical intractability, the numerical integrations for higher virtual transitions were evaluated computationally using an optimized Python script. The resulting second-order energy corrections for these discrete excited configurations are tabulated in Appendix B, alongside comprehensive computational methodologies.

Summing the unperturbed reference baseline and the cumulative perturbation corrections—where the total second-order correction \(E^{(2)}\) converges numerically to approximately \(-4.35\,\text{eV}\)—yields the total energy of the system: \[E_{\text{total}} = E^{(0)} + E^{(1)} + E^{(2)}\] \[E_{\text{total}} \approx -275.51\,\text{eV} + 83.50\,\text{eV} - 4.35\,\text{eV} = -196.36\,\text{eV}\]

3.3.2 Virtual Double-Electron Excitations↩︎

Although single-electron virtual transitions are constrained by selection rules exclusively to higher \(s\)-orbitals, a significant portion of the residual correlation energy is governed by simultaneous double-electron excitations, where the conservation of total angular momentum is preserved globally without requiring individual orbital constraints. This missing energy component—excluding minor relativistic contributions—originates predominantly from coupled two-electron excitations involving the valence \(2s\) electron and one of the core \(1s\) electrons. Evaluating these double excitations demands mapping out non-spherical configurations (\(p^2, d^2, f^2, g^2\)) and isolating combined multi-particle states that satisfy a total angular momentum of \(L=0\).

To preserve physical validity, the total angular momentum of the excited two-electron subsystem must couple exclusively to a total value of 0 or 1, thereby conserving the total atomic angular momentum at \(L=0\). These coupling constraints dictate that when the unexcited spectator core electron is in a spin-up state (\(\alpha\)), the excited subsystem must couple to a singlet state (\(S=0\)), whereas a spin-down spectator (\(\beta\)) necessitates a triplet coupling (\(S=1\)). According to rigorous configurations maps established in literature [7], the total number of valid multi-particle configurations satisfying these conditions is forty-five for the lithium atom. The multi-body wave function encompassing these expansions can be structured as: \[\Psi_{\text{3E}} = K \cdot 2s^{\prime\prime} + \Phi_{1} + \Phi_{2}\] where \(2s''\) represents the Slater determinant adjusted for the coupled electrons, while \(\Phi_1\) and \(\Phi_2\) represent the multi-body interactions between the outer valence electron and the inner core electron of opposite and parallel spins, respectively, defined explicitly as: \[\begin{align} \Phi_{1} &= (1s)^{2}1s^{\prime\prime} + (1s1s^{\prime})1s^{\prime\prime} + (2p)^{2}1s^{\prime\prime} + (1s)^{2}2s + (2p^{\prime\prime})^{2}1s + (3d^{\prime})^{2}1s \\ \Phi_{2} &= (3p2p^{\prime\prime})[^3S]1s + (2p^{\prime\prime}3p^{\prime\prime})[^3S]1s + (5d3d^{\prime})[^3S]1s + (2p^{\prime\prime}3p^{\prime\prime})[^3S]2s \end{align}\]

The listed configurations are the subset of excited-electron arrangements that satisfy the required quantum numbers of the target Lithium state,where the orbital symbols (\(1s\) \(2p\) \(3d\) etc.) specify the occupied orbital symmetries, primes (both \('\) and \(''\) ) simply do label different correlated orbitals of the same type as the unprimed ones and do not refer to "more excitement" etc., and the symbol \(^{3}S\) after the spatial configurations of some of the pairs denotes that that pair is specifically coupled to a triplet S state. Incorporating these simultaneous double excitations up to the virtual \(5g\) shell yields an additional correlation correction of \(-0.814\,\text{eV}\) [7]. This shifts the total converged non-relativistic energy to \(-197.174\,\text{eV}\), successfully reducing the relative error against the reference baseline [3] to \(3.11\%\). This result is integrated into the scope of this work to demonstrate the physical significance of simultaneous virtual multi-particle excitations and to expand the theoretical framework. Because these terms were extracted from reference data rather than evaluated directly, this corrected value is omitted from subsequent graphical comparisons and discussion sections.

3.3.3 Comparison of Zeroth, First, and Second-Order Energy Corrections↩︎

A systematic comparison of the sequential perturbation levels reveals a monotonic convergence toward the non-relativistic reference ground-state energy of approximately \(-203.5\,\text{eV}\) [3]. As illustrated in Fig. 1, the zeroth-order approximation, which entirely neglects inter-electronic repulsion, yields a severely underestimated energy of \(-275.51\,\text{eV}\), resulting in a substantial relative error of \(35.84\%\). Introducing the static Coulomb and exchange interactions via the first-order correction drastically shifts the energy eigenvalue upward to \(-192.01\,\text{eV}\), minimizing the error to \(5.65\%\). Finally, incorporating dynamic electron correlation through second-order virtual \(s\)-orbital transitions refines the total calculated energy to \(-196.36\,\text{eV}\), driving the relative error down to \(3.5\%\).

This systematic trend quantitatively demonstrates how consecutive higher-order perturbation terms successfully diminish physical errors by capturing the missing electron correlation energy. Physically, the unperturbed baseline is unrealistically negative because the non-interacting model overestimates the binding energy between the unshielded nucleus and the individual electrons. The subsequent first and second-order corrections introduce the necessary repulsive potential energy, shifting the eigenvalues upward and systematically converging toward the exact reference benchmark.

3.3.4 Higher-Order Perturbation Corrections↩︎

Although first and second-order perturbation theories capture the bulk of static Coulomb repulsion and primary pair-correlation effects, achieving full spectroscopic accuracy necessitates the consideration of higher-order terms (\(n \geq 3\)). These terms account for increasingly complex multi-electronic virtual excitations and non-linear correlation couplings. Within the Rayleigh-Schrödinger framework under intermediate normalization (\(\langle \Psi_0^{(0)} | \Psi_0^{(k)} \rangle = \delta_{0k}\)), the general \(n\)-th order energy correction is defined recursively by projecting the perturbation operator onto the \((n-1)\)-th order wave function correction [2]: \[E^{(n)} = \langle \Psi_0^{(0)} | \hat{H}^{(1)} | \Psi_0^{(n-1)} \rangle\]

Although evaluating the perturbed wave function \(\Psi_0^{(n-1)}\) becomes computationally demanding for complex many-body atomic systems, Wigner’s \(2k+1\) theorem adopted in the perturbation theory in all fundamental quantum mechanical references that is originally formulated in Wigner’s book [8] mathematically optimizes this process by demonstrating that knowledge of the perturbed wave function up to order \(k\) is sufficient to calculate the exact energy correction up to order \(2k+1\).

4 Variational Method Using Two Effective Nuclear Charge Parameters Estimation↩︎

4.1 Theoretical Framework and Trial Wavefunction↩︎

While perturbation theory provides a systematic correction based on the unperturbed nuclear field, it struggles to fully capture the dynamic screening effects introduced by the electronic electric fields [4] [1]. To establish a more accurate upper bound for the ground state energy of the multi-electron structure, the variational method is employed in conjunction with the frozen-core mean field approximation.[4] [2]

Because electrons are indistinguishable fermions, the total wavefunction must be completely antisymmetric, which is naturally enforced using a Slater determinant. This formalism inherently accounts for standard Coulomb repulsion and introduces exchange energy (an attractive quantum mechanical correction between electrons with parallel spins). In the frozen-core mean field approximation method [6], the framework is to freeze the states of two electrons to calculate the average electrostatic potential field (\(V_{\text{eff}}\)) projected onto the third. This approximation approach is assumed for the calculation of potential energies of each electron pair, whose integrals are then solved analytically by exploiting the gamma function. Finally, the results of each of the potential calculations and the kinetic energies are summed up to construct a trial energy function which is then to be minimized in order to find a value for the constant \(\alpha\) (i.e. the effective nuclear charge).
To implement this minimization mathematically, we construct a trial wave function (Ansatz) featuring a variational parameter (\(\alpha\)). \[E(\alpha) = E_{trial} = \frac{\langle \psi | \hat{H} | \psi \rangle}{\langle \psi | \psi \rangle} \ge E_{ground}\] Physically, (\(\alpha\)) represents the effective nuclear charge (\(Z_{\text{eff}}\)) felt by the individual electrons due to the shielding effect of the core charge distribution[2] [1]. Since \(Z_{\text{eff}}\) is always smaller then the total Z, \(E_{\text{ground}}\) will always be smaller then the \(E_{\text{trial}}\). The used Hamiltonian operator is the same of that in the perturbation part which is given in Equation 3.

4.2 Trial Wavefunctions and Mathematical Identities↩︎

To proceed with the variational method, we construct hydrogenic trial wavefunctions. To maintain correct SI dimensionality, we explicitly incorporate the Bohr radius (\(a_0\)) alongside the dimensionless variational parameter \(\alpha\) where \(\alpha = Z_{\text{eff}}\) [[9]][10]: \[\begin{align} \psi_{1s}(r) &= \sqrt{\frac{\alpha^3}{\pi a_0^3}} e^{-\alpha r/a_0} \\ \psi_{2s}(r) &= \sqrt{\frac{\alpha^3}{32\pi a_0^3}} \left( 2 - \frac{\alpha r}{a_0} \right) e^{-\alpha r/(2a_0)} \end{align}\]

Because the probability wavefunctions must be normalized over all space, the spherical volume element integrates out the angular dependence: \[\int d^3{r} = 4\pi \int_{0}^{\infty} r^2 dr\]

For the subsequent radial integrations, we will frequently utilize the standard complete gamma function identity: \[\int_{0}^{\infty} r^n e^{-br} dr = \frac{n!}{b^{n+1}}\]

4.3 Kinetic Energy Expectation Values \(\langle T \rangle\)↩︎

The kinetic energy expectation value in SI units is evaluated using the operator \(\hat{T} = -\frac{\hbar^2}{2m_e}\nabla^2\), which can be alternatively expressed in its symmetric gradient form: \[\langle T \rangle = \frac{\hbar^2}{2m_e} \int |\nabla \psi|^2 d^3{r}\]

Let us evaluate this for the \(1s\) orbital. The radial gradient is: \[\nabla \psi_{1s} = \frac{\partial \psi_{1s}}{\partial r} = -\frac{\alpha}{a_0} \sqrt{\frac{\alpha^3}{\pi a_0^3}} e^{-\alpha r/a_0}\]

Substituting this into the kinetic energy integral yields: \[\begin{align} \langle T_{1s} \rangle &= \frac{\hbar^2}{2m_e} (4\pi) \int_{0}^{\infty} \left( -\frac{\alpha}{a_0} \sqrt{\frac{\alpha ^3}{\pi a_0^3}} \right)^2 e^{-2\alpha r/a_0} r^2 dr \\ &= \frac{2\hbar^2 \alpha^5}{m_e a_0^5} \int_{0}^{\infty} r^2 e^{-2\alpha r/a_0} dr \end{align}\]

Applying the gamma integral (equation 72) identity with \(n=2\) and \(b=2\alpha/a_0\): \[\langle T_{1s} \rangle = \frac{2\hbar^2 \alpha^5}{m_e a_0^5} \left( \frac{2!}{(2\alpha/a_0)^3} \right) = \frac{\hbar^2}{2m_e a_0^2} \alpha^2\]

Recognizing that the physical constant group \(\frac{\hbar^2}{2m_e a_0^2}\) represents exactly one Rydberg of energy (\(\approx 13.6057\;\text{eV}\)), we can express the kinetic energy expectations directly in electron-volts. Using the general principal quantum number dependence (\(\langle T \rangle \propto 1/n^2\)), the kinetic energies for the \(n=1\) and \(n=2\) states are: \[\langle T_{1s} \rangle = \frac{a^2}{1^2}(13.6057\;\text{eV}) = \alpha^2 (13.6057\;\text{eV})\] \[\langle T_{2s} \rangle = \frac{\alpha^2}{2^2} (13.6057\;\text{eV}) = \frac{\alpha^2}{4} (13.6057\;\text{eV})\] which align with the result in [9]. In the Hamiltonian of multi-electron atoms, there is only one kinetic energy expectation value for each orbital because the nucleus is static. However, there are two potential terms that must be physically separated because the electron-nucleus potential energy \(V_{N}\) represents the attractive interaction binding electrons to the nucleus, whereas the electron-electron repulsion \(V_{ee}\) accounts for the opposing repulsive force between the electrons themselves. Including both terms is absolutely essential to accurately calculate the total energy of the system and to properly model the shielding effect that the electrons in the lithium atom exert on one another.

4.4 Nucleus Potential Energy \(\langle V_{N} \rangle\)↩︎

The Coulomb interaction between the electrons and the nucleus (\(Z=3\)) in SI units includes the potential operator \(-\frac{Z e^2}{4\pi\varepsilon_0 r}\). Evaluating the expectation value for the \(1s\) orbital: \[\begin{align} \langle V_{N, 1s} \rangle &= 4\pi \int_{0}^{\infty} \left( \sqrt{\frac{\alpha^3}{\pi a_0^3}} e^{-\alpha r/a_0} \right) \left(-\frac{3e^2}{4\pi\varepsilon_0 r}\right) \left( \sqrt{\frac{\alpha^3}{\pi a_0^3}} e^{-\alpha r/a_0} \right) r^2 dr \\ &= -\frac{12 e^2 \alpha^3}{4\pi\varepsilon_0 a_0^3} \int_{0}^{\infty} r e^{-2\alpha r/a_0} dr \\ &= -\frac{3 e^2 \alpha^3}{\pi\varepsilon_0 a_0^3} \left( \frac{1!}{(2\alpha/a_0)^2} \right) = -\frac{3 e^2 \alpha^3}{\pi\varepsilon_0 a_0^3} \left( \frac{a_0^2}{4\alpha^2} \right) = -3\alpha \left( \frac{e^2}{4\pi\varepsilon_0 a_0} \right) \end{align}\]

The physical constant group \(\left( \frac{e^2}{4\pi\varepsilon_0 a_0} \right)\) represents exactly 2 Rydbergs (\(\approx 27.2114\;\text{eV}\)). Thus, the potential energy can be expressed directly in eV: \[\langle V_{N, 1s} \rangle = -3 \alpha(27.2114\;\text{eV}) = 2(-3 \alpha)(1\;\text{Rydberg)} = -6 \alpha (13.6057\text{eV})\]

Similarly, for the \(2s\) orbital, using the general scaling relation \(\langle V \rangle \propto -\frac{Z \cdot Z_{\text{eff}}}{n^2}\): \[\langle V_{N, 2s} \rangle = -\frac{3 \alpha}{4} \left( \frac{e^2}{4\pi\varepsilon_0 a_0} \right) = 2(-\frac{3}{4} \alpha) (1\;\text{Rydberg}) =- \frac{3}{2} \alpha(13.6057\;\text{eV})\]

4.5 Electron-Electron Repulsion \(\langle V_{ee} \rangle\)↩︎

For the inter-electronic repulsion, we utilize the multipole expansion where \(\frac{1}{|{r}_1 - {r}_2|} \rightarrow \frac{1}{r_>}\), with \(r_> \equiv \max(r_1, r_2)\). The SI operator is \(\frac{e^2}{4\pi\varepsilon_0 r_{12}}\).

\[\begin{align} \langle V_{1s, 1s} \rangle &= \frac{e^2}{4\pi\varepsilon_0} (4\pi)^2 \left( \frac{ \alpha^3}{\pi a_0^3} \right)^2 \int_{0}^{\infty} r_1^2 e^{-2 \alpha r_1/a_0} \Bigg[ \int_{0}^{r_1} \frac{1}{r_1} r_2^2 e^{-2 \alpha r_2/a_0} dr_2 \\ &\quad + \int_{r_1}^{\infty} \frac{1}{r_2} r_2^2 e^{-2 \alpha r_2/a_0} dr_2 \Bigg] dr_1 \end{align}\]

Applying integration by parts for the terms inside the brackets (which physically represents calculating the classical electrostatic potential of a spherical charge cloud), we obtain: \[\text{Inner Bracket} = \frac{a_0^3}{4 \alpha ^3 r_1} \left[ 1 - e^{-2 \alpha r_1/a_0} \left(1 + \frac{\alpha r_1}{a_0}\right) \right]\]

Substituting this back into the main outer integral yields: \[\begin{align} \langle V_{1s, 1s} \rangle &= \frac{e^2}{4\pi\varepsilon_0} \left(\frac{4 \alpha ^3}{a_0^3}\right) \int_{0}^{\infty} \left( r_1 e^{-2 \alpha r_1/a_0} - r_1 e^{-4 \alpha r_1/a_0} - \frac{\alpha}{a_0} r_1^2 e^{-4 \alpha r_1/a_0} \right) dr_1 \\ &= \frac{e^2}{4\pi\varepsilon_0} \left(\frac{4 \alpha ^3}{a_0^3}\right) \left[ \frac{1!}{(2 \alpha /a_0)^2} - \frac{1!}{(4 \alpha /a_0)^2} - \frac{\alpha}{a_0}\frac{2!}{(4 \alpha /a_0)^3} \right] \\ &= \frac{e^2}{4\pi\varepsilon_0} \left(\frac{4 \alpha ^3}{a_0^3}\right) \left[ \frac{a_0^2}{4 \alpha ^2} - \frac{a_0^2}{16 \alpha ^2} - \frac{2a_0^2}{64 \alpha ^2} \right] \\ &= \frac{e^2}{4\pi\varepsilon_0 a_0} (4\alpha ) \left( \frac{5}{32} \right) = \frac{5}{8}\alpha \left( \frac{e^2}{4\pi\varepsilon_0 a_0} \right) \\ &= 2(\frac{5}{8}\alpha ) (1\;\text{Rydberg}) = \frac{5}{4}\alpha (13.6057\;\text{eV}) \end{align}\]

By equivalent spatial integration methods, the Coulomb integral for the \(1s\) and \(2s\) interaction evaluates to: \[\langle V_{1s, 2s} \rangle = \frac{17}{81}\alpha \left( \frac{e^2}{4\pi\varepsilon_0 a_0} \right) = 2(\frac{17}{81}\alpha) (1\;\text{Rydberg}) = \frac{34}{81}\alpha (13.6057\;\text{eV})\]

4.6 Total Energy and Minimization↩︎

The total ground-state energy of the \(1s^2 2s^1\) configuration is the linear combination of its nine constituent expectation values: two kinetic terms and two nuclear potential terms for the \(1s\) core, one kinetic and one nuclear potential term for the \(2s\) valence electron, one core-core (\(1s-1s\)) repulsion term, and two core-valence (\(1s-2s\)) repulsion terms.

Factoring out the standard Rydberg energy scale, the parametric energy function \(E(\alpha)\) in SI units is (factoring out the Rydberg constant from all terms) assembled as: \[\begin{align} E(\alpha) &= 2\langle T_{1s} \rangle + \langle T_{2s} \rangle + 2\langle V_{N, 1s} \rangle + \langle V_{N, 2s} \rangle + \langle V_{1s, 1s} \rangle + 2\langle V_{1s, 2s} \rangle \\ E(\alpha) &= 2\left(\alpha^2\right) + \frac{\alpha^2}{4} + 2(-6\alpha) + \left(-\frac{3\alpha}{2}\right) + \frac{5\alpha}{4} + 2\left(\frac{34\alpha}{81}\right) \\ E(\alpha) &= \frac{9\alpha^2}{4} - \frac{3697}{324}\alpha \end{align}\]

Since there are 2 electrons in 1s orbital, the energies \(T_{1s}\) and \(V_{N,1s}\) are multiplied by a factor of 2, and the same is correct for \(V_{1s,2s}\) because each electron in the 1s orbital interact with the electron in 2s orbital.

According to the variational principle, this function represents an upper bound to the true ground-state energy. To find the optimal effective nuclear charge that minimizes the system’s energy, we differentiate \(E(\alpha)\) with respect to the variational parameter \(\alpha\) and set it to zero: \[\frac{dE}{d\alpha} = \frac{9}{2}\alpha - \frac{3697}{324} = 0 \implies \alpha \approx 2.536\]

This value indicates that the \(1s\) electrons experience an effective nuclear charge of approximately \(+2.536e\), explicitly demonstrating the shielding effect. Substituting the optimized \(\alpha\) parameter back into the energy equation yields the variational minimum energy limit, which we convert directly to electron-volts: \[\begin{align} E_{\text{min}} = \left[ \frac{9}{4}(2.536)^2 - \frac{3697}{324}(2.536) \right](1\;\text{Rydberg})\\ \approx -14.467 (13.6057\;\text{eV}) \approx {-196.83\;\text{eV}} \end{align}\]

The non-relativistic reference limit for the lithium ground state is approximately \(-203.5\;\text{eV}\) [3]. The simple single-parameter variational ansatz captures the energy to within \(\sim 3.3\%\) accuracy, confirming its effectiveness despite the absence of dynamic angular correlation.

It should be emphasized that the single-parameter variational framework presented in this study operate strictly within the non-relativistic regime. While this approach effectively models the primary electronic screening and correlation effects, achieving true spectroscopic precision necessitates the inclusion of relativistic and quantum electrodynamic (QED) corrections which were not considered in this study. For a recent and comprehensive theoretical treatment incorporating these relativistic effects into the ground-state calculation of lithium using single-parameter variational approach, the reader is referred to the detailed framework presented in [10] which obtains an error percentage of \(\sim 2.5\%\) that is closer to the reference value than the result obtained here \(\sim 3.3\%\).

4.7 Two Parameter Estimation↩︎

Although one parameter estimation yielded a result 3.3% close to the reference result, it is more proper to expand the estimation to a second parameter. That is even more reasonable because of the two orbital nature of the Lithium atom. For this, let the previously calculated shielding parameter be \(\alpha\), that is the shielding effect felt by electrons in the 1s orbital. Therefore we introduce \(\beta\) which is the shielding effect felt in the 2s orbital. [9] [11]

The calculation algorithm of a second shielding parameter follows the same steps as that of the first one. The difference is seen, however, in the definition of the hydrogenic wave functions of the orbitals. Specifically in the definition of the 2s orbital’s wave function where \(\beta\) is used instead of \(\alpha\), and is defined as \(\beta = Z-1\) where the total atomic number \(Z\) is lowered significantly by a factor of 1 to add even more shielding effect due to the 1s electrons. The wave functions are defined as:

\[\psi(1s) = \sqrt{\frac{\alpha^3}{\pi}} e^{-\alpha r}\]

where \(\psi(1s)\) represents the spatial distribution of the core electrons, with \(\alpha\) acting as the optimized effective charge under core-core screening. \[\psi(2s) = \sqrt{\frac{\beta^3}{32\pi}} (2 - \beta r) e^{-\frac{\beta r}{2}}\] where \(\psi(2s)\) denotes the wave function of the valence electron, modulated by \(\beta\) to simulate the robust shielding provided by the underlying \(1s^2\) shell.

By evaluating the corresponding Hamiltonian matrix elements (factoring out the Rydberg energy unit from all the terms) using these non-orthogonal trial states within the Slater determinant framework, the individual kinetic, potential, Coulomb, and exchange expectation values are explicitly derived as:

Core Kinetic Energy (\(T_{1s}\)):

Represents the expectation value of the kinetic energy for a single electron residing in the inner \(1s\) shell: \[\langle T_{1s} \rangle = \frac{\alpha^2}{2}\]

Valence Kinetic Energy (\(T_{2s}\)):

Describes the expectation value of the kinetic energy for the outer-shell \(2s\) electron: \[\langle T_{2s} \rangle = \frac{\beta^2}{8}\]

Core Electron-Nucleus Potential (\(V_{N,1s}\)):

Quantifies the attractive Coulomb potential energy between the nucleus of charge \(Z\) and a core \(1s\) electron: \[\langle V_{N,1s} \rangle = -Z\alpha\]

Valence Electron-Nucleus Potential (\(V_{N,2s}\)):

Represents the attractive Coulomb potential energy acting on the valence \(2s\) electron: \[\langle V_{N,2s}(\beta) \rangle = \frac{-Z\beta}{4}\]

Core-Core Coulomb Repulsion (\(V_{1s1s}\)):

Accounts for the direct electrostatic repulsion energy between the two spin-paired electrons sharing the spatial \(1s\) orbital: \[\langle V_{1s1s}(\alpha) \rangle = \frac{5}{8}\alpha\]

Direct Core-Valence Repulsion (\(V_{1s2s}\)):

Represents the classical direct Coulomb repulsion integral between a core \(1s\) electron and the valence \(2s\) electron: \[\langle V_{1s2s}(\alpha, \beta) \rangle = \alpha\beta \frac{\beta^4 + 10\alpha\beta^3 + 8\alpha^4 + 20\alpha^3\beta + 12\alpha^2\beta^2}{(2\alpha + \beta)^5}\]

Kinetic Energy Exchange Contribution (\(T_{1s2s}\)):

A non-classical kinetic term that arises because the \(1s\) and \(2s\) spatial wave functions are no longer strictly orthogonal when governed by different screening parameters (\(\alpha \neq \beta\)): \[\langle T_{1s2s}(\alpha, \beta) \rangle = -4\sqrt{2}\alpha^{\frac{5}{2}}\beta^{\frac{5}{2}} \frac{\beta - 4\alpha}{(2\alpha + \beta)^4}\] This particular term is worth noting as it is non-existing in the case of one shielding parameter. That is because when the shielding effect is reduced to one parameter only inter-orbital interaction converges to zero due to the orthogonality of the wave functions of both orbitals.

Core-Valence Exchange Repulsion (\(V_{12,12}\)):

This inter-electronic exchange integral represents the energy correction due to the Pauli exclusion principle as the electrons have parallel spin states(both electrons have a spin state of 2 as shown in the subscript): \[\langle V_{12,12}(\alpha, \beta) \rangle = 16\alpha^3\beta^3 \frac{13\beta^2 + 20\alpha^2 - 30\beta\alpha}{(\beta + 2\alpha)^7}\]

Spatial Overlap Integral (\(S_{1s,2s}\)):

is a measurment of the actual spatial non-orthogonality [1] which is due to the wave function mixing between both orbitals: \[\langle S_{1s2s}(\alpha, \beta) \rangle = 32\sqrt{2}\alpha^{\frac{3}{2}}\beta^{\frac{3}{2}} \frac{\alpha - \beta}{(2\alpha + \beta)^4}\]

The next step is to sum the separate energies to find the trial energy as a function of both \(\alpha\) and \(\beta\):

\[\langle E(\alpha, \beta) \rangle = \frac{ \begin{array}{l} 2\langle T_{1s}(\alpha) \rangle + \langle T_{2s}(\beta) \rangle - \langle T_{1s}(\alpha) \rangle S_{1s2s}(\alpha, \beta)^2 - 2\langle T_{1s2s}(\alpha, \beta) \rangle S_{1s2s}(\alpha, \beta)... \\ +2\langle V_{N1s}(\alpha) \rangle + \langle V_{N2s}(\beta) \rangle - \langle V_{N1s}(\alpha) \rangle S_{1s2s}(\alpha, \beta)^2 - 2\langle V_{N1s2s}(\alpha, \beta) \rangle S_{1s2s}(\alpha, \beta)... \\ +2\langle V_{1s2s}(\alpha, \beta) \rangle + \langle V_{1s1s}(\alpha) \rangle - 2\langle V_{1112}(\alpha, \beta) \rangle S_{1s2s}(\alpha, \beta) - \langle V_{1212}(\alpha, \beta) \rangle \end{array} }{1 - S_{1s2s}(\alpha, \beta)^2}\]

The same reasoning used in the previous parameter estimation methodology for multiplying some of the terms by a factor of 2 is also applied here, particularly to account for the presence of two electrons in the 1s orbital, and thus 2 pairs with the electron in the 2s orbital each. The term in the denominator is simply to ensure the normalization of the results of the total energy expectation value.

The final step is minimizing the trial energy expectation value, and that follows the same logic as before the substituting with the values of \(\alpha\) and \(\beta\) in the formula of \(\langle E(\alpha, \beta) \rangle\), yielding (putting back the 2 Rydberg energy units):

\[\begin{gather} \text{Given} \quad \frac{\partial}{\partial\alpha}\langle E(\alpha, \beta) \rangle = 0 \quad \frac{\partial}{\partial\beta}\langle E(\alpha, \beta) \rangle = 0 \\[10pt] \begin{pmatrix} \alpha \\ \beta \end{pmatrix} = \text{Find}(\alpha, \beta) \quad \begin{pmatrix} \alpha \\ \beta \end{pmatrix} = \begin{pmatrix} 2.6797 \\ 1.8683 \end{pmatrix} \quad \\[10pt] \langle E(\alpha, \beta) \rangle = -201.187 \text{ eV} \end{gather}\] The reference value is \(E_{ref} = -203.5\;\text{eV}\) [3]. The error percentage then can be calculated as: \[\left| \frac{\langle E(\alpha, \beta) \rangle - E_{ref}}{E_{ref}} \right| = 1.13\%\] The results all agree with those in [9].

5 Conclusion↩︎

In this study, the ground-state energy of the lithium atom (\(1s^2 2s^1\)) was systematically investigated by deploying and contrasting two foundational pillars of quantum mechanics: Rayleigh-Schrödinger perturbation theory and the variational method. By rigorously enforcing the Pauli exclusion principle via a totally antisymmetric Slater determinant framework, the non-interacting unperturbed baseline energy was established at \(-275.51\,\text{eV}\).

First-order perturbation theory successfully incorporated the averaged electrostatic inter-electronic repulsion through the analytical evaluation of static Coulomb and quantum exchange integrals, providing a massive first-order upward shift to \(-192.01\,\text{eV}\). To account for dynamic electron correlation, higher-order virtual single-electron excitation channels were computationally evaluated. Resolving the spherically symmetric virtual transitions from the \(3s\) up to the \(6s\) Rydberg shells yielded a cumulative second-order energy correction of \(-4.35\,\text{eV}\), culminating in a total perturbative ground-state energy of \(-196.36\,\text{eV}\).
Fig. 2 illustrates the radial probability density functions for the \(1s\) and \(2s\) (as given in Appendix A.1) orbitals across sequential perturbation treatments along with a curve produced using the fitted parameter obtained upon applying the variational approach for each orbital to be used as a reference point as it is the closest obtained result to the reference energy value(\(-203.5\) eV). Comparing the perturbation order correction curves of the same orbital, it is clearly seen that the zeroth and first-order profiles exhibit highly localized, unshielded hydrogen-like behaviors, while the second-order corrected curves visibly shift closer to realistic multi-body configurations where the probability is more distributed over space. Comparing the curves of both orbitals from different perturbation orders, it is notable that the valence \(2s\) orbital distribution becomes broader and exhibits a noticeable outward spatial relaxation compared to the \(1s\) curves, demonstrating that the perturbation formalism successfully captures effective physical shielding profiles and electron-electron avoidance configurations without the explicit apriori inclusion of variational shielding parameters into the calculation.

To better capture the non-linear relaxation and radial screening effects introduced by the electronic charge clouds, a non-orthogonal multi-parameter variational approach was subsequently implemented. By introducing independent variational shielding parameters for the inner and outer shells, this model explicitly map the differential screening landscapes of the atom. Minimizing the global energy functional isolated the optimal effective nuclear charges at \(\alpha = 2.6797\) for the \(1s\) core and \(\beta = 1.8683\) for the \(2s\) valence electron. This multi-parameter optimization successfully established a highly accurate upper energy bound of \(-201.187\,\text{eV}\).

A definitive comparative insight emerges when evaluating these advanced theoretical frameworks against the non-relativistic reference ground-state benchmark of approximately \(-203.5\,\text{eV}\)[3], as visually summarized in Fig. 3. The multi-parameter variational framework outpaced the second-order single-excitation perturbative treatment, shrinking the relative error margin to a mere \(1.13\%\), whereas the second-order perturbation expansion converged to a \(3.51\%\) relative error.

The superior accuracy of the variational approach is directly attributed to its capacity to dynamically adjust the underlying orbital dimensions to simulate physical shielding. Rather than experiencing the bare electrostatic force of the unshielded nucleus (\(Z=3\)), the inner core and outer valence electrons adjust to effective atomic numbers of \(Z_{\text{eff},1s} = 2.6797\) and \(Z_{\text{eff},2s} = 1.8683\), respectively.

Finally, the residual energy gap between the optimized variational limit and the reference baseline highlights the persistent frontiers of many-body atomic physics. This remaining discrepancy is predominantly driven by angular electron correlation effects—such as simultaneous multi-particle double-excitations into non-spherical (\(p^2, d^2\)) spatial manifolds—alongside minor relativistic and fine-structure corrections. Ultimately, both frameworks successfully demonstrate how systematic mathematical refinements can transform a crude independent-particle approximation into a highly predictive model of multi-electron atomic structures.

6 Detailed Evaluation of Coulomb and Exchange Integrals↩︎

6.1 Radial Functions and Corresponding Probability Distributions↩︎

The standard unperturbed hydrogenic single-particle radial functions for a general bare nuclear charge \(Z\) are defined below, where the characteristic atomic length scale is governed by the Bohr radius (\(a_0\)) and energies scale in terms of the Rydberg constant (\(1\,\text{Ry} \approx 13.6057\;\text{eV}\)) [2]: \[\begin{align} R_{1s}(r) &= 2\left(\frac{Z}{a_0}\right)^{3/2}e^{-\frac{Zr}{a_0}} \\ R_{2s}(r) &= \frac{1}{2\sqrt{2}}\left(\frac{Z}{a_0}\right)^{3/2}\left(2-\frac{Zr}{a_0}\right)e^{-\frac{Zr}{2a_0}} \end{align}\] The corresponding radial probability distributions, mapping the isotropic spatial profiles of the shells, are given by [2]: \[\begin{align} P_{1s}(r) &= r^2|R_{1s}(r)|^2 \\ P_{2s}(r) &= r^2|R_{2s}(r)|^2 \end{align}\] The python code and its detailed description provided in [12] demonstrate the algorithm used to compute the radial probability density functions and plot the results as a function of the distance from the nucleus for each perturbation correction as well as the variational approximation in order to provide further insight regarding the frameworks used in this study.

6.2 General Radial Reduction via Multipole Expansion↩︎

The inter-electronic Coulomb repulsion operator \(1/r_{12} = 1/|\mathbf{r}_1 - \mathbf{r}_2|\) can be fundamentally represented using the standard spherical multipole expansion: \[\frac{1}{|\mathbf{r}_1 - \mathbf{r}_2|} = \sum_{l=0}^{\infty} \sum_{m=-l}^{l} \frac{4\pi}{2l+1} \frac{r_{<}^{l}}{r_{>}^{l+1}} Y_{lm}^*(\Omega_1) Y_{lm}(\Omega_2)\] For multi-electron systems confined strictly to spherically symmetric \(s\)-orbitals (\(l=0, m=0\)), the angular integrations over the solid angles (\(d\Omega_1 d\Omega_2\)) isolate the monopole term exclusively. Due to the geometric orthogonality of spherical harmonics, all higher-order angular channels (\(l \ge 1\)) vanish identically. Thus, the local repulsion operator collapses strictly to its isotropic monopole bound, \(1/r_{>}\), where \(r_{>} \equiv \max(r_1, r_2)\) and \(r_{<} \equiv \min(r_1, r_2)\) [2].

Splitting the radial domain to accommodate this piecewise boundary conditions transforms the generic spatial double integral for inter-electronic interactions into a structured two-region functional layout: \[J_{ab} = \int_0^\infty r_1^2 \rho_a(r_1) \left[ \frac{1}{r_1} \int_0^{r_1} r_2^2 \rho_b(r_2)\,dr_2 + \int_{r_1}^\infty r_2 \rho_b(r_2)\,dr_2 \right] dr_1\] Physically, this represents the classical electrostatic potential energy separating two localized charge clouds. In direct harmony with Gauss’s Law, the bracketed term dictates the effective potential field generated by electron \(b\): the charge density enclosed within the interior sphere of radius \(r_1\) acts as a centralized point charge, whereas the exterior distribution acts as a uniform potential shell.

6.3 Explicit Derivation of the \(1s\text{–}1s\) Coulomb Integral (\(J_{1s,1s}\))↩︎

To evaluate the static repulsion within the core shell, we isolate the \(1s\text{–}1s\) Coulomb expectation value. Let the local effective potential generated by the \(1s\) density cloud be designated as \(U_{1s}(r_1)\): \[U_{1s}(r_1) = \mathcal{I}_1(r_1) + \mathcal{I}_2(r_1) = \frac{1}{r_1} \int_{0}^{r_1} r_2^2 \rho_{1s}(r_2) \, dr_2 + \int_{r_1}^{\infty} r_2 \rho_{1s}(r_2) \, dr_2\] To streamline the evaluation, we define the scaled scaling variable \(\alpha = 2Z/a_0\), simplifying the core density expression to \(\rho_{1s}(r) = \frac{\alpha^3}{2} e^{-\alpha r}\).

Step 1: Evaluation of the Interior Integral \(\mathcal{I}_1(r_1)\)
Applying standard integration by parts to the interior radial boundary yields: \[\mathcal{I}_1(r_1) = \frac{1}{r_1} \int_{0}^{r_1} r_2^2 \left( \frac{\alpha^3}{2} e^{-\alpha r_2} \right) dr_2 = \frac{1}{r_1} \left[ 1 - e^{-\alpha r_1} \left( \frac{\alpha^2 r_1^2}{2} + \alpha r_1 + 1 \right) \right]\]

Step 2: Evaluation of the Exterior Integral \(\mathcal{I}_2(r_1)\)
Similarly, evaluating the exterior shell component from \(r_1\) to infinity gives: \[\mathcal{I}_2(r_1) = \int_{r_1}^{\infty} r_2 \left( \frac{\alpha^3}{2} e^{-\alpha r_2} \right) dr_2 = e^{-\alpha r_1} \left( \frac{\alpha^2 r_1}{2} + \frac{\alpha}{2} \right)\]

Step 3: Construction of the Global Core Potential \(U_{1s}(r_1)\)
Summing the two sub-domains triggers an algebraic cancellation of the linear spatial terms, leaving a compact, closed-form classical screening field: \[U_{1s}(r_1) = \mathcal{I}_1(r_1) + \mathcal{I}_2(r_1) = \frac{1}{r_1} - e^{-\alpha r_1} \left( \frac{\alpha}{2} + \frac{1}{r_1} \right)\]

Step 4: Final Integration of \(J_{1s,1s}\)
Projecting this effective electrostatic potential back across the outer core distribution yields: \[\begin{align} J_{1s,1s} &= \int_{0}^{\infty} r_1^2 \left( \frac{\alpha^3}{2} e^{-\alpha r_1} \right) \left[ \frac{1}{r_1} - e^{-\alpha r_1} \left( \frac{\alpha}{2} + \frac{1}{r_1} \right) \right] dr_1 \\ &= \frac{\alpha^3}{2} \left[ \int_{0}^{\infty} r_1 e^{-\alpha r_1} \,dr_1 - \int_{0}^{\infty} e^{-2\alpha r_1} \left( \frac{\alpha r_1^2}{2} + r_1 \right) \,dr_1 \right] \end{align}\] Evaluating these standard continuous definitive channels via the Euler gamma identity (Equation 72) resolves to: \[J_{1s,1s} = \frac{\alpha^3}{2} \left[ \frac{1}{\alpha^2} - \frac{1}{8\alpha^2} - \frac{1}{4\alpha^2} \right] = \frac{\alpha^3}{2} \left( \frac{5}{8\alpha^2} \right) = \frac{5\alpha}{16}\] Restoring the physical variables (\(\alpha = 2Z/a_0\)) maps the final \(1s\text{–}1s\) Coulomb repulsion energy strictly as a function of the core nuclear boundary: \[J_{1s,1s} = \frac{5Ze^2}{32\pi \epsilon_0 a_0} \quad \xrightarrow{\text{for } Z=3} \quad J_{1s,1s} \approx 51.02\;\text{eV}\]

6.4 Explicit Derivation of the \(1s\text{–}2s\) Coulomb Integral (\(J_{1s,2s}\))↩︎

To solve the inter-shell interaction \(J_{1s,2s}\), the effective monopole potential fields are mapped over the nodally complex \(2s\) orbital density distribution: \[U_{2s}(r_1) = \frac{1}{r_1} \int_0^{r_1} r_2^2 \rho_{2s}(r_2)\,dr_2 + \int_{r_1}^\infty r_2 \rho_{2s}(r_2)\,dr_2\] Performing systematic integration by parts over the multi-termed polynomial array of the \(2s\) shell structures yields the complete spatial screening field: \[U_{2s}(r_1) = \frac{1}{r_1} - e^{-\frac{Zr_1}{a_0}} \left( \frac{1}{r_1} + \frac{3Z}{4a_0} + \frac{Z^2r_1}{4a_0^2} + \frac{Z^3 r_1^2}{8a_0^3} \right)\] Next, this comprehensive potential profile is evaluated across the inner core density \(\rho_{1s}(r_1)\): \[J_{1s,2s} = 4\left( \frac{e^2}{4 \pi \epsilon_0}\right)\left(\frac{Z}{a_0}\right)^3 \int_0^\infty \left[ r_1 e^{-\frac{2Zr_1}{a_0}} - e^{-\frac{3Zr_1}{a_0}} \left( r_1 + \frac{3Z}{4a_0}r_1^2 + \frac{Z^2}{4a_0^2}r_1^3 + \frac{Z^3}{8a_0^3}r_1^4 \right) \right] \,dr_1\] Invoking the gamma identity independently for each unique polynomial channel yields: \[\begin{align} \int_0^\infty r_1 e^{-\frac{2Zr_1}{a_0}} \,dr_1 &= \frac{a_0^2}{4Z^2} \\ \int_0^\infty r_1 e^{-\frac{3Zr_1}{a_0}} \,dr_1 &= \frac{a_0^2}{9Z^2} \\ \frac{3Z}{4a_0} \int_0^\infty r_1^2 e^{-\frac{3Zr_1}{a_0}} \,dr_1 &= \frac{3Z}{4a_0} \left( \frac{2a_0^3}{27Z^3} \right) = \frac{a_0^2}{18Z^2} \\ \frac{Z^2}{4a_0^2} \int_0^\infty r_1^3 e^{-\frac{3Zr_1}{a_0}} \,dr_1 &= \frac{Z^2}{4a_0^2} \left( \frac{6a_0^4}{81Z^4} \right) = \frac{a_0^2}{54Z^2} \\ \frac{Z^3}{8a_0^3} \int_0^\infty r_1^4 e^{-\frac{3Zr_1}{a_0}} \,dr_1 &= \frac{Z^3}{8a_0^3} \left( \frac{24a_0^5}{243Z^5} \right) = \frac{a_0^2}{81Z^2} \end{align}\] Summing the collective exterior screened contributions results in: \[\frac{a_0^2}{Z^2} \left( \frac{1}{9} + \frac{1}{18} + \frac{1}{54} + \frac{1}{81} \right) = \frac{16a_0^2}{81Z^2}\] Subtracting this grouped valuation from the localized lead channel isolates the analytical solution: \[J_{1s,2s} = 4\left(\frac{e^2}{4 \pi \epsilon_0}\right)\left(\frac{Z}{a_0}\right)^3 \left[ \frac{a_0^2}{4Z^2} - \frac{16a_0^2}{81Z^2} \right] = \frac{17Z}{81a_0}\left(\frac{e^2}{4 \pi \epsilon_0}\right)\] Evaluating this final expression for the lithium baseline (\(Z=3\)) yields: \[J_{1s,2s} = \frac{17(3)}{81a_0}\left(\frac{e^2}{4 \pi \epsilon_0}\right) \approx 17.135\;\text{eV} \implies 2J_{1s,2s} \approx 34.27\;\text{eV}\] This closed-form mathematical expression aligns perfectly with standard quantum chemistry reference literature [2].

6.5 Explicit Derivation of the Exchange Integral (\(K_{1s,2s}\))↩︎

Because the active spatial operators within the non-classical exchange channel map coordinate permutations identically across both coordinates, structural symmetry reduces the spatial double integral to: \[K_{1s,2s} = 2\left(\frac{e^2}{4 \pi \epsilon_0a_0}\right)\int_0^\infty r_1 f(r_1)\left[\int_0^{r_1} r_2^2 f(r_2)\,dr_2\right]\,dr_1\] where the core-valence cross-density distribution is defined as \(f(r) = R_{1s}(r)R_{2s}(r) = \frac{Z^3}{\sqrt{2}}(2-Zr)e^{-\alpha r}\) with a defined scaling index \(\alpha = 3Z/2\).

Evaluating the incomplete interior radial channel yields: \[\int_0^{r_1} r_2^2 f(r_2)\,dr_2 = \frac{Z^3}{\sqrt{2}} \left[ 2 \int_0^{r_1} r_2^2 e^{-\alpha r_2} \,dr_2 - Z \int_0^{r_1} r_2^3 e^{-\alpha r_2} \,dr_2 \right]\] Resolving the integration by parts arrays reveals a critical algebraic phenomenon: the constant, linear, and quadratic spatial boundaries within the exponential factors sum to zero. The solitary surviving term is the cubic envelope: \[\int_0^{r_1} r_2^2 f(r_2)\,dr_2 = \frac{Z^3}{\sqrt{2}} e^{-\alpha r_1} \left( \frac{2r_1^3}{3} \right)\] Projecting this collapsed inner profile directly into the outer integration layer simplifies the exchange value to: \[K_{1s,2s} = \left(\frac{e^2}{4 \pi \epsilon_0a_0}\right)\frac{2Z^6}{3} \int_0^\infty \left( 2r_1^4 - Z r_1^5 \right) e^{-3Zr_1} \,dr_1\] Evaluating the definitive continuous matrices via the complete gamma function maps the solution as: \[K_{1s,2s} = \left(\frac{e^2}{4 \pi \epsilon_0a_0}\right)\frac{2Z^6}{3} \left[ 2 \frac{4!}{(3Z)^5} - Z \frac{5!}{(3Z)^6} \right] = \left(\frac{e^2}{4 \pi \epsilon_0a_0}\right)\frac{16}{729}Z\] For \(Z=3\), this evaluates directly to \(K_{1s,2s} \approx 1.79\;\text{eV}\), matching reference values [2].

6.6 The Second-Order Virtual Transition Integrals (\(3s\) Channel)↩︎

To isolate the dynamic correlation contributions from the lowest virtual \(s\)-wave channel, we utilize the unperturbed hydrogenic \(3s\) single-particle radial function: \[R_{3s}(r) = \frac{2Z^{3/2}}{81\sqrt{3}} \left( 27 - 18Zr + 2Z^2r^2 \right) e^{-Zr/3}\] The second-order Coulomb transition matrix element (\(J'_{3s,1s}\)) defines the electrostatic interaction between the static \(1s\) core distribution and the non-local \(2s \rightarrow 3s\) transition density: \[J'_{3s,1s} = \iint \phi_{3s}^*(\mathbf{r}_1)\phi_{1s}^*(\mathbf{r}_2) \left( \frac{1}{r_{12}} \right) \phi_{2s}(\mathbf{r}_1)\phi_{1s}(\mathbf{r}_2) \,d^3\mathbf{r}_1 \,d^3\mathbf{r}_2\] Integrating the angular variables reduces this to a radial task utilizing the exact core potential \(U_{1s}(r_1)\): \[J'_{3s,1s} = \int_0^\infty r_1^2 \left[ R_{3s}(r_1)R_{2s}(r_1) \right] U_{1s}(r_1) \,dr_1\] Expanding the cross-density radial product \(R_{3s}(r_1)R_{2s}(r_1)\) yields the following polynomial array: \[R_{3s}(r_1)R_{2s}(r_1) = \frac{Z^3}{81\sqrt{6}} \left( 54 - 63Zr_1 + 22Z^2r_1^2 - 2Z^3r_1^3 \right) e^{-5Zr_1/6}\] Evaluating this spatial projection analytically against the core screening field isolates the transition value: \[J'_{3s,1s} \approx 4.12\;\text{eV}\]

6.7 The Transition Exchange Integral (\(K'_{3s,1s}\))↩︎

The non-classical transition exchange integral accounts for spatial coordinate exchange between the virtual excited states for parallel spin channels: \[K'_{3s,1s} = \iint \phi_{3s}^*(\mathbf{r}_1)\phi_{1s}^*(\mathbf{r}_2) \left( \frac{1}{r_{12}} \right) \phi_{1s}(\mathbf{r}_1)\phi_{2s}(\mathbf{r}_2) \,d^3\mathbf{r}_1 \,d^3\mathbf{r}_2\] Radially, this is evaluated as a mutual overlap integral between two distinct mixed transition distributions, \(f_A(r) = R_{3s}(r)R_{1s}(r)\) and \(f_B(r) = R_{2s}(r)R_{1s}(r)\): \[K'_{3s,1s} = \int_0^\infty r_1^2 f_A(r_1) \left[ \frac{1}{r_1} \int_0^{r_1} r_2^2 f_B(r_2) \,dr_2 + \int_{r_1}^\infty r_2 f_B(r_2) \,dr_2 \right] dr_1 \approx 0.92\;\text{eV}\] Combining these individual transition matrices maps the global perturbative numerator contribution: \[2J'_{3s,1s} - K'_{3s,1s} = 2(4.12\;\text{eV}) - 0.92\;\text{eV} \approx 7.33\;\text{eV}\]

---

7 Computational Evaluation of Higher-Order Virtual Transitions↩︎

To resolve the quantitative contributions of high-quantum-number virtual \(s\)-orbitals (\(n \ge 3\)) to the second-order perturbation energy (\(E^{(2)}\)), numerical integrations were executed using the Python environment. The algorithmic engine relies on the scipy.integrate and scipy.special libraries to process the rapid spatial oscillations of the excited single-particle functions.

7.1 Computational Architecture↩︎

The processing script was structured into a modular framework mapping the underlying physical mechanics:

• Radial Wave Functions (R_ns): For virtual states where \(n \ge 3\), Generalized Laguerre Polynomials (scipy.special.eval_genlaguerre) were deployed to precisely track the multiple internal nodes and high spatial extent of the excited states.

• Isotropic Simplification (multipole_term): Because the active configuration is restricted to spherically symmetric \(s\)-waves (\(L=0\)), the multi-center potential reductions collapse identically to the \(1/r_{>}\) monopole term, drastically optimizing processing efficiency.

• Numerical Integration (compute_integral): The six-dimensional spatial integrals were analytically downscaled to two-dimensional radial coordinates. SciPy’s adaptive quadrature algorithm (dblquad) was utilized to calculate the spatial profiles. To guarantee rigorous algorithmic convergence over the infinite radial domain, absolute (epsabs) and relative (epsrel) error metrics were maintained strictly at \(1\times10^{-5}\).

• Energy Conversion and Iteration: Matrix evaluations natively computed in atomic units (Hartrees) were transformed to electron-volts via the scaling index \(27.2114\;\text{eV/Hartree}\). The dynamic unperturbed energy denominators (\(\Delta E\)) and separate second-order correlation values (\(E^{(2)}\)) were evaluated iteratively across individual principal quantum shells from \(n=3\) to \(n=6\).

7.2 Calculated Energy Contributions↩︎

The variables tracks within the computational data loop are defined as follows:

• \(J^{\prime}\) and \(K^{\prime}\): The direct Coulomb and quantum exchange transition matrix elements, mapping the electrostatic repulsion and Pauli correlation channels, respectively.

• \(\Delta E\): The unperturbed energy denominator, calculated as the zero-order eigenvalue difference between the ground-state configuration and the virtual excited state (\(E_0^{(0)} - E_{ns}^{(0)}\)).

• \(E^{(2)}\) Contribution: The net second-order correlation energy injected by that specific principal quantum channel.

As the principal quantum number \(n\) advances, the energy gap \(\Delta E\) widens significantly, and the spatial overlap between the compact ground states and the diffuse virtual orbitals decays exponentially. This induces a rapid asymptotic convergence in the correlation corrections, as detailed in Table 1.

The underlying Python scripts and automation repositories are hosted publicly on GitHub [12].

7.3 Selection Rules and Angular Constraints↩︎

These virtual expansions were strictly confined to \(s\)-orbital structures (\(l=0\)). Because the lithium ground state is completely spherically symmetric (\(L=0\)) and the electrostatic Coulomb perturbation operator acts as a pure scalar, any single-electron transition into higher angular momentum configurations (such as \(p, d,\) or \(f\) blocks) yields angular integrals that vanish identically due to orthonormality constraints [6]. This systematic computational layout successfully outlines the boundary limits of \(s\)-wave perturbation theory, illustrating the quantitative necessity for multi-parameter variational methods to close the residual electronic correlation gap.

References↩︎