Electron density analysis of two-electron systems confined by prolate spheroids with hard walls

The electron density of two-electron systems, He and H2, was analyzed when prolate spheroids with hard walls confine these systems. For this purpose, Hartree–Fock equations were solved using Roothaan's approach with a basis set defined in prolate spheroidal coordinates imposing Dirichlet boundary conditions. Total energy, its components, and orbital energies were analyzed for several confinements, and some of these results were compared with those reported by other authors to test the performance of the proposed approach. For both systems, the electron density exhibits a maximum value out of the nuclear region for extreme confinements. The chemical bond for H2 was analyzed through the concepts of the quantum theory of atoms in molecules, concluding that the chemical bond of this molecule disappears under extreme conditions. For this system, estimations of the correlation energy indicate that this is a small contribution to the total energy, and the Hartree–Fock method contains the necessary elements to describe the chemical bond for strong confinements.


Introduction
The change of the electronic structure of atoms and molecules by confinement is an active field in physics and chemistry, where many confinement effects have been studied using model potentials to mimic spatial restrictions [1][2][3][4][5][6].The hard walls model has been used for atoms and molecules to simulate the confinement by imposing Dirichlet boundary conditions such that the wave function or electron density is canceled on the system's boundaries.The electronic configuration of atoms confined by hard walls exhibits important changes; for example, some alkali or alkaline-earth metal atoms exhibit electronic transitions that induce orbital energies crossing [7,8].Such transitions have been observed experimentally [9,10] and corroborated by hard confinement.We have mentioned one example where the electronic structure of atoms is altered by hard confinement.Several examples of these systems exist [11][12][13][14][15]. Confinement imposed by soft walls is another important model potential for studying the electronic structure of confined atoms and molecules, although this confinement is not considered in this article [16][17][18][19][20][21][22][23].
Confined one-and two-electron molecules have been studied considering different shapes and methods [6,[24][25][26][27][28][29][30][31].However, the main interest in such reports is related to energy analysis, and the chemical bond has not been the main target.In this sense, Ley-Koo and Cruz solved analytically the Schrödinger equation for the + H 2 molecule confined by prolate spheroids with hard walls [24].By using the analytic solution of this molecule, Hernández-Esparza et al analyzed the chemical bond of the + H 2 molecule through the quantum theory of atoms in molecules (QTAIM) [32], focusing on the electron density changes induced by hard confinement [33].One crucial observation obtained from this study is that QTAIM detects an attractor in the region between both nuclei.This behavior is a consequence of confinement, where the kinetic energy commands the potential energy when hard walls confine this molecule.In addition to the + H 2 molecule, some efforts have been made to study molecules with more than one electron confined by hard walls, in particular for H 2 [6,[25][26][27][28][29][30] and He 2 [34].This article aims to analyze the electron density in two-electron systems, He and H 2 , confined by prolate spheroids with hard walls.The solution of the Hartree-Fock equations is through the Roothaanʼs approach using a basis set defined within the prolate spheroidal coordinate system.The QTAIM was used to analyze the chemical bond.

Methodology
In the Born-Oppenheimer approximation, the Schrödinger equation for the H 2 molecule is given in atomic unis (a.u.) by where 2f represents the internuclear distance, R AB .The treatment of this equation will be under the Hartree-Fock (HF) approach [35].Thus, the wave-function is represented by a Slater determinant built with spinorbitals, χ(r, ω) = f(r)σ(ω), with σ = α or σ = β.In this case we treat with the ground state, χ 1 = fα and χ 2 = fβ.In general, for a closed-shell system the energy has the expression where N represents the number of electrons in the system, in our case N = 2, the Coulomb integral, J ab , is defined as and the exchange integral, K ab as ò ò The set of orbitals that minimize E HF must satisfy the HF equations Within HF theory the electron density has the expression, In our methodology, the solution of the HF equations is through the Roothaan's approach where HF orbitals are expanded in terms of a basis set functions, g μ , In this case, M represents the maximum number of functions in the basis set.In this way, the matrix representation of HF equations is r , 1 3 * and the matrix c contains the coefficients m { } ( ) c i and the diagonal part of matrix ò contains the orbital energies, {ò i }.
In this article, we use prolate spheroidal coordinates, defined as where r iM represents the distance between electron i and nucleus M. The transformation to cartesian coordinates is where 1 ξ < ∞ , −1 η 1, 0 f < 2π, and f is chosen to be equal to half the internuclear distance.Because we are interested in the electron density behavior along the internuclear axis, ρ(r) and its derivatives are evaluated in Cartesian coordinates.
In our study, we use a basis set with functions of the form where P b j is a Legendre polynomial and b j is an even integer.This basis set assures that the wave function vanishes when ξ = ξ 0 .In this way we are imposing the Dirichlet boundary conditions associated to hard walls or impenetrable walls.
One-electron and two-electron integrals are evaluated before performing each self-consistent field process (SCF).Integration over η does not have spatial restrictions.Consequently, all integrals involving this variable have the following form, which is written in terms of Wigner 3-j symbols.The integration over the f variable is trivial due to the azimutal symmetry of the molecule.All integrals over the ξ coordinate can be written in terms of the following integrals where: Integration of equations ( 22) and ( 23) are carried out via iterative integration by parts.It is worth noting that for large n values the integrand is highly oscillatory which leads to precision issues.We solve possible numerical instabilities by using arbitrary precision via the GNU MPFR library [36].For each R AB and ξ 0 , the exponent α in (20) is optimized to obtain the lowest HF energy.All matrix elements and the SCF procedure were implemented through an application designed by us and based on the C programming techniques.
For the electron density analysis, ρ(r), and its gradient, ∇ρ(r), we use Cartesian coordinates since from equations ( 14)-( 15) if x = y = 0 then the z-axis is relevant for our analysis, such that and In this way ρ(r) = ρ(z).Mathematica© v11 was used for the visual analysis of these quantities.

Hydrogen atom
Before analyzing the chemical bond in the H 2 molecule, we tested our implementation with a hydrogen atom confined by a prolate spheroidal box.The results of this system from the different approaches are presented in table 1. Naturally, exact results are the main target of any new methodology designed to solve the Schrödinger equation.From table 1, our approach is convenient for describing the hydrogen atom confined by a prolate spheroidal box using three functions in the basis set.Although our approach with three functions in the basis set gives numerical results close to the exact ones, the performance of the basis set is better for small values of ξ 0 .

Helium atom
The helium atom confined by a prolate box with hard walls is another system that can be studied using our methodology.This system's total energy and electron-electron components are listed in table 2. In this table, we also report the results published by Corella-Madueño et al [41] for this system under the same confinement.
From this table, we appreciate the good performance of our numerical approach, which gives lower values than those reported by Corella-Maduño et al [41], particularly for extreme confinements.The electron-electron interaction involved in the HF method has the Coulomb contribution, J, and exchange energy, K.Both quantities are reported in table 2. This table shows that for small values of ξ 0 , J increases and K decreases almost at the same rate; this behavior is similar to that presented by atoms confined by spheres with hard walls [22].Thus, the confinement imposed by the hard walls induced more negative exchange energy values.To the best of our knowledge, this is the first time that this observation has been confirmed for two different shapes of confinement with hard walls.The occupied orbital energy represents the ionization potential, I, under the Koopmans theorem since I = -ò for the HF approach.We observe from table 2 that in the region 3.0 < ξ 0 < 4.0, this quantity is zero.The behavior presented by I under a prolate spheroidal box with hard walls is similar to that observed when a sphere imposes confinement.The electron density, ρ(r), of this atom for several values of ξ 0 is presented in figures 1 and 2. In this case, ρ(r) is evaluated over the focal axis, where the helium nucleus is at the focus on the right side.In the early stages of the confinement, the electron density is increased at the nucleus position, showing its maximum value right at this point.For extreme confinements (small values of ξ 0 ), the electron density maximum shifts towards the box's center.From these figures, it is clear that the action of extreme confinements on ρ(r) reduces the interaction of the nucleus over the electron density.

H 2 molecule
We begin our study of the H 2 system by computing the total energy of the free molecule for R AB =1.4 a.u.From our approach, we obtain E HF = −1.32066674a.u., which differs by 0.14% from the HF-limit (E HF−limit = −1.133628674a.u.[42]).Therefore, we expect a similar performance of our approach for the confined H 2 .
Batael et al have reported results for the H 2 confined by the same conditions imposed in this article [26].However, in that report, Batael et al used a correlated wave function, and consequently, they have total energies  lower than the HF limit.In table 3 we are reporting E HF from our approach, contrasted with those obtained by Batael et al.
In table 3, we do not want to compare our HF energy (E HF ) with that reported by Batael et al (E) since our wave function is not correlated.Instead, we want to report the energy difference E − E HF since this approximates the correlation energy.Values for R AB and ξ 0 are taken from table 1 of reference [26].Table 3 shows that E − E HF is almost constant, around −0.02 a.u., for different values of R AB ξ 0 ,except for R AB ξ 0 = 13.We must mention an inconsistency in a value reported by Batael et al for R AB = 0.4493 and ξ 0 = 4.52 the product R AB ξ 0 is 2.030 836 and in table 1 of reference [26] Batael is reporting as R AB ξ 0 = 2.We think that the value reported by Batael et al E = 4.6212 a.u. is a typing error, and they must check this result.From this discussion, it is clear the the correlation energy is a small component of the total energy.Therefore, the HF method is adequate to describe the chemical bond for the confined H 2 molecule.On the other side, it is important to mention that for the He atom confined by spherical hard walls, the correlation energy is almost constant for any confinement [43].Connecting the observations for H 2 and He confined by hard walls, the correlation energy for two-electron systems is almost constant for different shapes of confinements.
The total energy, E HF , and components of the energy are reported in table 4.This table shows a similar behavior to that obtained for the confined helium atom.The Coulomb electron-electron energy increases, and exchange energy decreases when the H 2 molecule is confined.In both cases, the rate change is almost the same.From here, it is evident that the kinetic energy changes rapidly when the confinement is increased.For this reason, the kinetic energy commands the behavior of the total energy in extreme confinements.Naturally, this behavior must be mapped over the electron density.
The electron density along the internuclear axis is presented in figures 3 and 4 for several values of ξ 0 .For convenience, the internuclear axis is aligned with the z-axis and for both figures, R AB = 1.4 a.u.From these figures, there are several interesting results.For some values of ξ 0 , the electron density increases and reaches a maximum value.However, this property goes to zero for small values of ξ 0 .This behavior is enhanced in figure 5. Figure 3 shows that the maximum value of ρ(r) at the nucleus is reached when this property is constant along the internuclear axis.Another interesting result from figures 3 and 4 is related to the electron density evaluated at the middle of the internuclear axis.This quantity always is increased when ξ 0 is decreased.
According to the QTAIM, a bond critical point (BCP) is determined when ∇ρ = 0 and there is one negative eigenvalue of the Hessian of ρ(r).In our case, two hydrogen atoms form a chemical bond, and we have a BCP in the middle of the internuclear axis, as can be appreciated from figure 3.However, from figure 4, it is clear that the curvature of the electron density changes in the middle of the internuclear axis and becomes a maximum; this is no longer a BCP.
Recently, Hernández-Esparza et al demonstrated that the electron density at the BCP of the + H 2 molecule increases as ξ 0 decreases [33].Furthermore, they showed that the electron density of this one electron system starts to localize in both nuclei as ξ 0 decreases to a maximum value x max 0 which is a bifurcation point at which the BCP lowers its rank and becomes a degenerate (2, −2) critical point.For x x < max 0 0 the electron density starts shifting towards the center of the focal axis, effectively destroying the common molecular structure of the + H 2 molecule.As we have discussed the H 2 molecule above, the behavior found by Hernández-Esparza et al for + H 2 is similar to the H 2 .Table 4.Total HF energy (EHF) and its components: kinetic (T), nucleus-electron (V), Coulomb electron-electron (J), exchange electron-electron (K), for the H 2 molecule confined by a prolate spheroidal box with hard walls.In this case, the inter-nuclear distance RAB = 1.4 a.u.All quantities are in atomic units.

Conclusions
Two-electron systems, He and H 2 , confined by prolate spheroids with hard walls, were studied in this article.
When the volume that contains the atom or the molecule is reduced, the energy components change.The energy change rate is similar for the Coulomb and exchange electron-electron contributions.As a preliminary result, we found that the correlation energy for the H 2 molecule is small and almost constant for several confinements.For the He atom, we found that the electron density along the focal axis grows around the nucleus as the parameter ξ 0 decreases, but as this parameter continues to decrease, the maximum electron density shifts toward the center of the box.In the case of the H 2 molecule, the electron density along the focal axis increases at the position of the nuclei as the parameter ξ 0 diminishes, and this is a constant between both nuclei for a critical value of x max 0 .For x x < max 0 0 , the electron density exhibits a maximum at the middle of the nuclei positions, and from the quantum theory of atoms in molecules, this point does not represent a bond critical point indicating that the chemical bond disappears.Our group is working to obtain a code for many-electron diatomic molecules under the Hartree-Fock method as a starting point since we are interested in correlated methods and excited states of confined diatomic molecules.

Figure 1 .
Figure 1.Electron density along the focal axis of the He atom for ξ 0 values between 20 and 1.6.

Figure 2 .
Figure 2. Electron density along the focal axis of the He atom for ξ 0 values between 1.6 and 1.25.

Figure 3 .
Figure 3. Electron density along the focal axis of the H 2 molecule for ξ 0 values between 10 and 2.

Figure 4 .
Figure 4. Electron density along the focal axis of the H 2 molecule for ξ 0 values between 2 and 1.1.

Figure 5 .
Figure 5. Electron density evaluated at a nucleus of the H 2 molecule.

Table 1 .
Total energy, in Hartrees, for the hydrogen atom confined by a prolate spheroidal cage with hard walls.

Table 2 .
Total energy for the helium atom confined by a prolate spheroidal cage with hard walls.J, K and ò represent Coulomb, exchange and orbital energy, respectively.All quantities are in Hartrees.

Table 3 .
[13]l energy for the hydrogen molecule confined by a prolate spheroidal box with hard walls.The values obtained from this approach (EHF) are compared with values obtained by Batael et al[13](E).All quantities are in Hartrees.