Symmetry Breaking and the Generation of Spin Ordered Magnetic States in Density Functional Theory Due to Dirac Exchange for a Hydrogen Molecule

Abstract

We study symmetry breaking in the mean field solutions to the electronic structure problem for the 2 electron hydrogen molecule within the Kohn Sham (KS) local spin density functional theory with Dirac exchange (the XLDA model). This simplified model shows behavior related to that of the (KS) spin density functional theory (SDFT) predictions in condensed matter and molecular systems. The KS solutions to the constrained SDFT variation problem undergo spontaneous symmetry breaking leading to the formation of spin ordered states as the relative strength of the non-convex exchange term increases. Numerically, we observe that with increases in the internuclear bond length, the molecular ground state changes from a paramagnetic state (spin delocalized) to an antiferromagnetic (spin localized) ground state and a symmetric delocalized (spin delocalized) excited state. We further characterize the limiting behavior of the minimizer when the strength of the exchange term goes to infinity both analytically and numerically. This leads to further bifurcations and highly localized states with varying character. Finite element numerical results provide support for the formal conjectures. Various solution classes are found to be numerically stable. However, for changes in the R parameter, numerical Hessian calculations demonstrate that these are stationary but not stable solutions.

A Appendix

We recall here the basics of the finite element methods we use in this work to numerically find the critical points of the XLDA Lagrangian. The numerical algorithms are implemented using a python implementation of the FENICS finite element package, (Dupont et al. 2003). Many of the tools we use here are discussed in more detail in the references (Bylaska et al. 2009; Hu 2014). For complete discussions of finite element methods, see the books of Axelsson and Alan Barker (2001), Bank and Dupont (1981), Braess (2001) and Brenner and Ridgway (2008). The method we develop here takes advantage of the sparsity of the FEM representation of the eigenvalue problem leading to an algorithm that scales linearly with number of basis functions. For resources on large-scale computing in computational chemistry, see Kendall et al. (2000) and Valiev et al. (2010).

1.1 I. The Finite Element Setup

Fig. 6

Finite element tetrahedron defining FEM elements and nodes. The L tetrahedral elements are identified by $e_l$. The global node are indexes by (m). For each element $e_l$ a, there is global node at each corner. Generally, local nodes belonging to individual tetrahedra are also defined. See Axelsson and Alan Barker (2001), Braess (2001), Brenner and Ridgway (2008) and Bylaska et al. (2009). These are not shown and are managed transparently by the Dolfin software (Logg et al. 2010)

We assume that the solution, $\psi _{\pm }$, exists in a bounded domain $\Omega \in {\mathscr {R}}^3$ that can be divided into a set of L non-overlapping tetrahedral elements, $\{e_l\}_{l=1}^L$ (Axelsson and Alan Barker 2001; Braess 2001; Brenner and Ridgway 2008; Bylaska et al. 2009), see Figs. 6 and 1. For the electronic structure problems, we are concerned with the atomic potentials represented by V(x) in the Hamiltonian below, (74), are singular. This leads to rapid variation in the solution to the eigenvalue problem in this region. In order to obtain accuracy, the FEM grid in this region must have a finer resolution as illustrated in Fig. 1 and discussed in Bylaska et al. (2009), Bylaska et al. (1995) and Kohn et al. (1997).

To construct the grid used in the calculation, we

1. Use BoxMesh (Logg et al. 2012) to generate a coarse mesh in a $50\times 50\times 50$ domain. The initial number of cells in each direction is 2. So the total number of tetrahedra will be 48 and the total number of vertices will be 27 in the coarse mesh.

2. Find the closest mesh grids to the nuclei and set the parameter cell_marker (Logg et al. 2012) true that tells the code to refine the mesh. If cell_marker = false, it means this mesh will NOT be refined.

3. Refine the whole grid for three cycles.

Generally, FEM nodes are located at corners, along boundaries or in the centers of tetrahedral regions (Axelsson and Alan Barker 2001; Braess 2001; Brenner and Ridgway 2008; http://hplgit.github.io). For the calculations, here, the nodes are located only at the corners of the tetrahedra. These nodes are shared by adjacent tetrahedra as in Fig. 6. Each tetrahedron l has four corner nodes. A global index identifies a node as in Fig. 6 (global node numbers in brackets). There are M global nodes in the construction. In actual calculations, a local node index identifying a corner global node with a basis function inside a particular tetrahedral is also defined in Dolfin (Axelsson and Alan Barker 2001; Braess 2001; Brenner and Ridgway 2008; Bylaska et al. 2009; http://hplgit.github.io) to identify variation associated with a node within a particular tetrahedron (Axelsson and Alan Barker 2001; Braess 2001; Brenner and Ridgway 2008; http://hplgit.github.io). The somewhat difficult book keeping problem of keeping track of the global variation in the basis functions consistent with their local behavior is taken care of nicely in the FEniCS software, see Logg et al. (2012).

For each node (with global index m and local index i) in a tetrahedral element, l, a finite element basis functions $\{\chi ^{e_l}_i\}$ is defined. In these calculations, the $\{\chi ^{e_l}_i\}$ are linear functions centered on local nodes i in element $e_l$ (Axelsson and Alan Barker 2001; Braess 2001; Brenner and Ridgway 2008; Bylaska et al. 2009; http://hplgit.github.io). For each global node m, the linear basis function $\{\chi ^{e_l}_i\}$ is 1 on global node i and zero on all other nodes contained in the tetrahedral elements containing global node m. For a particular tetrahedron, the linear basis associated with local node i of tetrahedral $e_l$, $\{\chi ^{e_l}_i\}$, has value only in tetrahedra $e_l$. Illustrations of how this works are given in [1]. The local node functions $\{\chi ^{e_l}_i\}$ can be assembled in functions centered around the global nodes with index m as the global basis functions $\eta _m(x)$.

A piecewise continuous function (here the approximated $\psi _{\pm }(x)$) can now be expanded as in Bylaska et al. (2009),

$$\begin{aligned} \psi _{\pm }(x)=\sum \limits _{m=1}^M c_{\pm ,m}\eta _{m}(x). \end{aligned}$$

(69)

Here, M is the dimension of space of global nodes and $c_{\pm ,m}$ is the coefficient of basis element $\eta _m$. The value of the $\psi _{\pm }$ on node m is $c_{\pm ,m}$.

1.2 II. The Generalized Eigenvalue Problem

With the above formulation, solving the Kohn–Sham minimization problem related to (3) leads to the generalized eigenvalue problem (in many of the following equations the ± ($+$ spin-up, − spin-down) notation has been suppressed to keep notation simple),

$$\begin{aligned} {\mathbf {H}}{\mathbf {c}}=\varvec{\epsilon }{\mathbf {S}}{\mathbf {c}}, \end{aligned}$$

(70)

or

$$\begin{aligned} \sum _n H_{mn}c_{n,k}=\epsilon _k \sum _n S_{mn}c_{n,k} \end{aligned}$$

(71)

where k identifies the kth eigenfunction,

$$\begin{aligned} S_{mn}=\int \limits _\Omega {\mathrm{d}}x \eta _m(x)\eta _n(x), \end{aligned}$$

(72)

and

$$\begin{aligned} H_{mn}=\frac{1}{2}\int \limits _\Omega {\mathrm{d}}x\nabla \eta _m(x) \nabla \eta _n(x) +\int \limits _\Omega {\mathrm{d}}x \eta _m(x) V_{\text {eff}} \eta _n(x) \end{aligned}$$

(73)

with $V_{\text {eff}}$ given by

$$\begin{aligned} {V_\mathrm{eff}=V_\mathrm{ext}(x)+V_\mathrm{ee}(\rho )+V_\mathrm{ex}(\rho ,\alpha )= V_\mathrm{ext}(x)+\int \frac{\rho (x')}{|x-x'|} {\mathrm{d}}x'+V_\mathrm{ex}(\rho _{\pm },\alpha ).}\nonumber \\ \end{aligned}$$

(74)

Here, $\rho (x)$ is the total electron density, and $\rho _\pm$ is the spin density, and $V_{\mathrm {ex}, \pm }(\rho _{\pm },\alpha )$ is given by the scaled Dirac form

$$\begin{aligned} V_{\mathrm{ex}, \pm }(\rho ,\alpha )=\alpha \rho _{\pm }^{1/3}. \end{aligned}$$

(75)

Note that, in the spin DFT, the exchange potential depends on the spin component, and thus, the effective Hamiltonian for the spin-up and spin-down orbitals is different; while the structure of the problem is the same, and hence, we keep the notation (e.g., for ${\mathbf {H}}$ and ${\mathbf {S}}$) independent of spin component. The matrix $H_{mn}$, (73), the overlap matrix $S_{mn}$, (72) and integrals over V(x) in (74) can be obtained from the FEniCS software (Logg et al. 2012). The calculation of these matrix is also carefully discussed for the electronic structure problem in Bylaska et al. (2009) and in general in Axelsson and Alan Barker (2001), Braess (2001) and Brenner and Ridgway (2008). The full potential $V_{\mathrm{eff}}$ given by (74) is a function of the density requiring that the eigenvalue problem, (70), be solved iteratively until self-consistency is achieved. We are only interested in the lowest ± energy solutions, though the methods can be modified to higher energy states as well.

1.3 III. Solution to the Generalized Eigenvalue Problem and the Associate Coulomb Problem

The objective of the calculation is the solution of the generalized eigenvalue problem, (70). However, this requires the input of a current estimate of the Classical potential $V_\mathrm{ee}$ required in $V_{\text {eff}}$, (74). This may be found as the solution to the Poisson equation

$$\begin{aligned} \Delta V_\mathrm{ee} =-4\pi \rho =-4\pi \left[ \big |\psi _{+}\big |^2+ \big |\psi _{-}\big |^2\right] \end{aligned}$$

(76)

To solve this PDE, $V_\mathrm{ee}$ is also expanded in the finite element basis as

$$\begin{aligned} V_\mathrm{ee}(x)=\sum \limits _{m=1}^Mv_m\eta _m(x). \end{aligned}$$

(77)

(76) is then represented by a system of linear equations giving the $\{v_m\}$. Given a solution to the Coulomb problem, (76), based on a current density for the iteration, the generalized eigenvalue problem, (70), is also solved in the $\{\eta _m\}_{m=1}^{M}$ finite element basis functions. Each molecular orbital is represented as

$$\begin{aligned} \psi _{\pm }=\sum _{\alpha =1}^{M} c_{\pm ,m}\eta _{m}. \end{aligned}$$

(78)

Note, for this problem, there is only one filled molecular orbital for each spin. The finite element discretization of the one-electron equation for the current iteration is given as, for $i = \pm$,

$$\begin{aligned} (T_i-\epsilon _{i})c_{i}=v_{i}. \end{aligned}$$

(79)

$c_{i}=\{c_{i1},\ldots ,c_{im}\}$ are the coefficients of molecular orbital in the expansions of finite element basis. The elements of the operator $(T_i-\epsilon _{i})$ are

$$\begin{aligned} (T_i-\epsilon _{i})_{mn}=\int _{\Omega }\left\{ \frac{1}{2}\nabla \eta _{m}\nabla \eta _{n}-\epsilon _{i}\eta _{m}\eta _{n}\right\} d x. \end{aligned}$$

(80)

The elements of $v_{i}$ are given by

$$\begin{aligned} (v_{i})_{m}=\int \eta _m(x) V_{\text {eff}, i}(x)\psi _{i}(x) {\mathrm{d}}x \end{aligned}$$

(81)

and are calculated in an iterative process in which $V_{\text {eff}, i}(x)$ is defined for step k from the results of the self-consistent solver in the prior iteration.

1.3.1 Details of the FEniCS Solver

The AMG solver is based on a V-cycle (Briggs et al. 2000; Hu 2014) with a maximum number of multi grid levels of 25. For each fine to course grid transfer, a single pre-smoothing step is taken. For each course to fine transfer, a single post-smoothing step is taken. These smoothing steps use a symmetric-SOR/Jacobi method. On the coarsest level, the course FEM equation is relaxed by Gaussian elimination. In the iteration, an energy correction step is applied to update new eigenvalues after the Helmholtz equation is solved for a set of $\epsilon ^{(k)}$ from the prior AMGCG cycle. The self-consistent solver converges when the total energy difference in two consecutive iterations is smaller than a selected tolerance. These are solved via the FEniCs code using the GMRES (Saad and Schultz 1986) and BOOMER AMG (Algebraic Multigrid (Briggs et al. 2000)) packages. The solution to this problem is of order M (Hu 2014). To improve the convergence of the solution, a preconditioning based on the algebraic multigrid method is used (Tatebe 1993; Van Emden and Yang 2002). For an introduction to multigrid methods and their application to problems in electronic structure, see Bramble (1993), Bramble and Pasciak (1987), Brandt (1986), Brandt et al. (1985), Briggs et al. (2000), Bylaska et al. (1995), Hackbusch (1985), Harrison et al. (2004), Kohn et al. (1997) and McCormick (1987).

1.4 IV. The Wavefunction and Orbital Energy Update and Iteration

1.4.1 IV.i: Some Preliminaries

The eigenvalue problem, (79), is solved using an iterative process in which for step k, the $\psi _{\pm }$ on right hand side of (79) and the orbital energies, $\epsilon _{\pm }^{(k)}$, at step k are assigned the values and functionality from the $k-1$ step. The iteration is developed with the intention of producing linear scaling in the number of FEM basis functions, M. This is achieved by developing a solver strategy that emphasizes the use of the operator $[(\nabla ^2-k^2)]^{-1}$ which is efficiently implemented in multigrid schemes.

The density functional equations leading to the values of $\epsilon _i^k$ and $\psi _i^k$ for the k values (update from $\psi _i^{(k-1)}$ to $\psi _i^{k}$) are written as

$$\begin{aligned}&{\left[ -\frac{1}{2}\nabla _{x}^2 -\epsilon _i^{(k-1)} \right] \psi _i^{k}(x)}\nonumber \\&\quad { = \left[ V_\mathrm{ext}^{(k-1)}(x) + V^{(k-1)}_\mathrm{ee}(\rho ^{(k-1)}(x)) + V_\mathrm{ex}^{(k-1)}(\rho ^{(k-1)}(x) ) \right] \psi _i^{(k-1)}(x)} \end{aligned}$$

(82)

where $V_\mathrm{ext}$ is the external potential from (1). The electron-electron Coulomb potential (calculated from FeniCS as above) is given by

$$\begin{aligned} \nabla ^2V^{(k-1)}_\mathrm{ee}(x)=-4\pi \rho ^{(k-1)}(x)=-4\pi \left[ \big |\psi ^{(k-1)}_{+}\big |^2+ \big |\psi ^{(k-1)}_{-}\big |^2\right] . \end{aligned}$$

(83)

The exchange potential is given by

$$\begin{aligned} V^{(k-1)}_\mathrm{ex}(x)=\alpha \rho ^{(k-1)}(x)^{\frac{1}{3}} . \end{aligned}$$

(84)

This is now a linear PDE of the form

$$\begin{aligned} \left[ -\frac{1}{2}\nabla ^2 -\epsilon ^{(k-1)}_i \right] \psi ^k_i(x) =f^{(k-1)}_i(x). \end{aligned}$$

(85)

Note that, all the potential terms in (82) have been collected in the function $f^{(k-1)}_i$. We calculate the solution to (85) using an efficient multigrid method. Because of the complexity of the grid, we use the AMGCG implemented in the FEniCS software (Logg et al. 2012).

1.4.2 IV.ii: Update of the Wavefunction, from $\psi ^{(k-1)}(x)$ to $\psi ^{k}(x)$

To initiate the $k^{\mathrm{th}}$ iteration ($\psi _i^{(k-1)}$ to $\psi _i^k$), we assume we have solutions $\psi _i^{(k-1)}(x)$ and $\epsilon _i^{(k-1)}$. The update of the wavefunction proceeds directly from (85) as

$$\begin{aligned} \psi ^{k}_i(x) =\left[ -\frac{1}{2}\nabla _{x}^2 -\epsilon ^{(k-1)}_i \right] ^{-1}f^{(k-1)}_i(x). \end{aligned}$$

(86)

All the functions in this equation are defined from the solution that we obtain from AMGCG.

To complete the iteration cycle, we also need an update of the orbital energy $\epsilon _i^{(k-1)}$ to $\epsilon _i^{k}$.

1.4.3 IV.iii: Update of the Orbital Energy

We assume we have $\psi _i^{(k-1)}(x)$ and $\epsilon _i^{(k-1)}$ and begin by defining two Greens functions:

The $(k-1)^{\mathrm{th}}$ Green’s function, $G^{(k-1)}_i$, with energy $\epsilon ^{(k-1)}$,

$$\begin{aligned} G^{(k-1)}_i=\biggl \{-\frac{1}{2}\nabla ^2-\epsilon _i^{(k-1)}\biggr \}^{-1} \end{aligned}$$

(87)

and a Green’s function, $G^{\mathrm{con}}_i$ with the converged DFT orbital energy (from converged solution to DFT equations), $\epsilon _i^{\mathrm{con}}$. This is given by

$$\begin{aligned} G^{\mathrm{con}}_i=\biggl \{-\frac{1}{2}\nabla ^2-\epsilon ^\mathrm{con}_i\biggr \}^{-1}. \end{aligned}$$

(88)

In the iteration, the updated $\epsilon ^{k}_i$ given by

$$\begin{aligned} \epsilon _i^{k}=\epsilon _i^{(k-1)} +\delta \epsilon _i^{k} \end{aligned}$$

(89)

is taken to be a good approximation to $\epsilon ^{\mathrm{con}}$. Using this in (88), we have

$$\begin{aligned} G^{\mathrm{con}}_i=\biggl \{-\frac{1}{2}\nabla ^2-\epsilon ^\mathrm{con}_i\biggr \}^{-1}=\biggl \{\frac{1}{2}\nabla ^2-(\epsilon _i^{(k-1)} + \delta \epsilon ^{k}_i)\biggr \}^{-1} \end{aligned}$$

(90)

the objective is to calculate an orbital energy correction from these equations using $\psi _i^{(k-1)}$.

The function $\psi ^{\mathrm{con}}_i$ satisfies the orbital PDE,

$$\begin{aligned} \psi _i^{\mathrm{con}}(x)=\biggl \{-\frac{1}{2}\nabla ^2-\epsilon ^\mathrm{con}_i\biggr \}^{-1}f^{\mathrm{con}}_i(x). \end{aligned}$$

(91)

In this equation, $\epsilon _i^{\mathrm{con}}$ is the converged orbital energy.

We assume that $\psi _i^{(k-1)}$ is a good approximation to $\psi ^{\mathrm{con}}_i$, i.e., that it approximately satisfies

$$\begin{aligned} \psi ^{(k-1)}_i(x)= & {} \biggl \{-\frac{1}{2}\nabla ^2-\epsilon ^\mathrm{con}_i\biggr \}^{-1}f_i^{(k-1)}(x)\nonumber \\= & {} \biggl \{-\frac{1}{2}\nabla ^2-(\epsilon _i^{(k-1)} + \delta \epsilon ^{k}_i) \biggr \}^{-1}f_i^{(k-1)}(x). \end{aligned}$$

(92)

Now we expand the full Green’s function (RHS) in the energy correction $\delta \epsilon ^{k}_i$ to obtain an equation that will update the orbital energy (find a correction to $\epsilon ^{(k-1)}_i$).

We use the operator identity,

$$\begin{aligned} \frac{1}{(1+a+b)}=\frac{1}{(1+a)} - \frac{1}{(1+a)}b \frac{1}{(1+a+b)}, \end{aligned}$$

(93)

to obtain

$$\begin{aligned} \begin{aligned} \biggl \{-\frac{1}{2}\nabla ^2 -\epsilon ^{(k-1)}_i -\delta \epsilon _i^{k}\biggr \} ^{-1}=&\biggl \{-\frac{1}{2}\nabla ^2 -\epsilon ^{(k-1)}_i \biggr \} ^{-1}-\biggl [\biggl \{-\frac{1}{2}\nabla ^2 -\epsilon ^{(k-1)}_i \biggr \} ^{-1} \\&\biggl \{-\delta \epsilon _i^{k}\biggr \}\biggl \{-\frac{1}{2}\nabla ^2 -\epsilon ^{(k-1)}_i -\delta \epsilon _i^{k}\biggr \} ^{-1}\biggr ]. \end{aligned} \end{aligned}$$

(94)

Iteration of this equation leads to an expression for the propagator to first order in $\delta \epsilon ^{k}_i$ as

$$\begin{aligned} \begin{aligned} \biggl \{-\frac{1}{2}\nabla ^2 -\epsilon ^{(k-1)}_i -\delta \epsilon _i^{k}\biggr \} ^{-1}=&\biggl \{-\frac{1}{2}\nabla ^2 -\epsilon ^{(k-1)}_i \biggr \} ^{-1}-\biggl [\biggl \{-\frac{1}{2}\nabla ^2 -\epsilon ^{(k-1)}_i \biggr \} ^{-1} \\&\times \biggl \{\delta \epsilon _i^{k}\biggr \}\biggl \{-\frac{1}{2}\nabla ^2 -\epsilon ^{(k-1)}_i \biggr \} ^{-1}\biggr ]. \end{aligned} \end{aligned}$$

(95)

We can use this result in (92) to give

$$\begin{aligned} \begin{aligned} \psi ^{(k-1)}_i(x)=&\biggl \{-\frac{1}{2}\nabla ^2-\epsilon _i^{(k-1)}\biggr \}^{-1}f^{(k-1)}_i(x) \\&-\biggl \{-\frac{1}{2}\nabla ^2-\epsilon _i^{(k-1)}\biggr \}^{-1}\delta \epsilon _i^{k}\biggl \{-\frac{1}{2}\nabla ^2-\epsilon _i^{(k-1)}\biggr \}^{-1}f^{(k-1)}_i(x). \end{aligned} \end{aligned}$$

(96)

This is more conveniently written in vector notation (Cohen-Tannoudji et al. 1991) as

$$\begin{aligned} \begin{aligned} \left| \psi ^{(k-1)}_i\right\rangle =&\biggl \{-\frac{1}{2}\nabla ^2-\epsilon _i^{(k-1)}\biggr \}^{-1} \left| f^{(k-1)}_i\right\rangle \\&-\biggl \{-\frac{1}{2}\nabla ^2-\epsilon _i^{(k-1)}\biggr \}^{-1} \delta \epsilon ^{k}_i\biggl \{-\frac{1}{2}\nabla ^2-\epsilon _i^{(k-1)}\biggr \}^{-1}\big |f^{(k-1)}_i\big \rangle \end{aligned} \end{aligned}$$

(97)

or

$$\begin{aligned} \begin{aligned} 0=&-\left| \psi ^{(k-1)}_i\right\rangle +\biggl \{-\frac{1}{2}\nabla ^2-\epsilon _i^{(k-1)}\biggr \}^{-1} \big |f^{(k-1)}_i\big \rangle \\&-\biggl \{-\frac{1}{2}\nabla ^2-\epsilon _i^{(k-1)}\biggr \}^{-1}\delta \epsilon ^{k}_i\biggl \{-\frac{1}{2}\nabla ^2 -\epsilon _i^{(k-1)}\biggr \}^{-1}\big |f^{(k-1)}_i\big \rangle . \end{aligned} \end{aligned}$$

(98)

Closing this equation on the left with $\left\langle f^{(k-1)}\right|$ gives a linear expression for $\delta \epsilon ^{k}_i$ which may be in terms of the $\psi ^{k}_i$, (86) , as

$$\begin{aligned} \begin{aligned} 0=&-{\bigg \langle f_i^{(k-1)}\bigg |\psi ^{(k-1)}_i\bigg \rangle }+{\bigg \langle f_i^{(k-1)}\bigg |\psi _i^{k}\bigg \rangle } \\&-\delta \epsilon ^{k}_i {\bigg \langle \psi ^{k}_i|\psi ^{k}_i\bigg \rangle }. \end{aligned} \end{aligned}$$

(99)

This may be solved for $\delta \epsilon ^{k}_i$ to obtain

$$\begin{aligned} \delta \epsilon ^{k}_i=\frac{-{\bigg \langle f^{(k-1)}_i\bigg |\psi ^{(k-1)}_i\bigg \rangle }+{\bigg \langle f_i^{(k-1)}\bigg |\psi ^{k}_i\bigg \rangle }}{{\bigg \langle \psi ^{k}_i\bigg |\psi ^{k}_i\bigg \rangle } }. \end{aligned}$$

(100)

This gives the update to $\epsilon _i^{(k-1)}$ via (89) to complete the $k^{\text {th}}$ iteration.

1.5 V. The Self-Consistent Iteration

Algorithm 1 summarizes the process followed by the self-consistent solver. An initial guess $(c_{i}^{0},\epsilon _{i}^{0}), i=1,\ldots ,n$ is given to start the self-consistent iterations. The solver stops when the total energy difference in two consecutive iterations is smaller than the tolerance TOL.

1.6 VI. Hessian Analysis

The model we investigate in this work is the local spin density approximation (LDA) (without correlation energy contributions). As above, the ground singlet state spin unrestricted density functional theory for this two-electron system defines two orbital wave functions $(\psi _+,\psi _-)$. The total energy functional is $E(\psi )$ is

$$\begin{aligned}&E_{\alpha }(\psi _+, \psi _-) \nonumber \\&\quad = \frac{1}{2} \int \big |\nabla \psi _+\big |^2 \,\mathrm {d}x + \frac{1}{2} \int \big |\nabla \psi _-\big |^2 \,\mathrm {d}x + \int V_R(x)\left( \big |\psi _+(x)\big |^2+\big |\psi _-(x)\big |^2\right) \,\mathrm {d}x \nonumber \\&\qquad + \frac{1}{2} \iint \frac{\left( \big |\psi _+(x)\big |^2+\big |\psi _-(x)\big |^2\right) \left( \big |\psi _+(y)\big |^2+\big |\psi _-(y)\big |^2\right) }{\big |x-y\big |} \,\mathrm {d}x \,\mathrm {d}y \nonumber \\&\qquad - \alpha \int \left( \big |\psi _+(x)\big |^{8/3} + \big |\psi _-(x)\big |^{8/3} \right) \,\mathrm {d}x, \end{aligned}$$

(101)

where $V_R(x)$ is the nuclear potential. The constraints on $(\psi _+,\psi _-)$ are

$$\begin{aligned} \int \big |\psi _{i}(x)\big |^2 \,\mathrm {d}x=1,i=+,- . \end{aligned}$$

(102)

We define the Lagrangian as

$$\begin{aligned}&L(\psi _+,\psi _-,\epsilon _+,\epsilon _-)\nonumber \\&\quad =\frac{1}{2} \int \big |\nabla \psi _+\big |^2 \,\mathrm {d}x + \frac{1}{2} \int \big |\nabla \psi _-\big |^2 \,\mathrm {d}x \nonumber \\&\qquad + \int V_R(x)\left( \big |\psi _+(x)\big |^2+\big |\psi _-(x)\big |^2\right) \,\mathrm {d}x \nonumber \\&\qquad + \frac{1}{2} \iint \frac{\left( \big |\psi _+(x)\big |^2+\big |\psi _-(x)\big |^2\right) \left( \big |\psi _+(y)\big |^2+\big |\psi _-(y)\big |^2\right) }{\big |x-y\big |} \,\mathrm {d}x \,\mathrm {d}y \nonumber \\&\qquad - \alpha \int \left( \big |\psi _+(x)\big |^{8/3} + \big |\psi _-(x)\big |^{8/3} \right) \,\mathrm {d}x \nonumber \\&\qquad - \epsilon _+\left( \int \big |\psi _+(x)\big |^2 \,\mathrm {d}x-1\right) -\epsilon _-\left( \int \big |\psi _-(x)\big |^2 \,\mathrm {d}x-1\right) , \end{aligned}$$

(103)

where $(\epsilon _+,\epsilon _-)$ are Lagrange multipliers.

Finding the stationary variation in (103) with respect to the functions $\psi _+$ and $\psi _-$ leads to effective one-electron eigenvalue equations.

$$\begin{aligned}&\frac{\delta L}{\delta \psi _i}=0\Rightarrow \left( -\frac{1}{2}\nabla ^2+ V_R(x)+ \int \frac{\left( \big |\psi _+(y)\big |^2+\big |\psi _-(y)\big |^2\right) }{\big |x-y\big |} \,\mathrm {d}y-\frac{4}{3}\alpha \big |\psi _{i}(x)\big |^{2/3}\right) \nonumber \\&\psi _{i}(x) = \epsilon _{i} \psi _{i}(x), i=\pm , \end{aligned}$$

(104)

where $(\psi _+,\psi _-)$ and $(\epsilon _+,\epsilon _-)$ satisfy normalization constraints.

In order to determine whether the stationary extremum of $L(\psi _+,\psi _-,\epsilon _+,\epsilon _-)$ with respect to functional variation are a maximum, a minimum or a saddle point, the second-order functional derivative (the Hessian matrix) may be analyzed (Hu 2014). In the following, the stationary solutions $(\psi _+,\psi _-)$ and their eigenvalues $(\epsilon _+,\epsilon _-)$ satisfy (104) and (102). $(\lambda _i, w_i)$ are eigenvalues and eigenvectors of the Hessian Matrix, ${\textbf {Hess}}$, defined as the solutions to the eigenvalue problem,

$$\begin{aligned} {\textbf {Hess}}\,\,w_i(y)=\int \frac{\delta ^2 L(\psi ,\epsilon )}{\delta \psi _i(x)\delta \psi _j(y)} w_i(y)\,\mathrm {d}y\bigg |_{(\psi ,\epsilon )=(\phi ,\epsilon )}=\lambda _iw_i(x), \ i=\pm ,\nonumber \\ \end{aligned}$$

(105)

where the Hessian matrix is defined as

$$\begin{aligned} {\textbf {Hess}}= \left( \begin{array}{cc} H_{11} &{} \int \frac{2\psi _-(y)}{\big |x-y\big |}\,\mathrm {d}y \psi _+(x) \\ \int \frac{2\psi _+(y)}{\big |x-y\big |}\,\mathrm {d}y \psi _-(x) &{} H_{22} \end{array} \right) \end{aligned}$$

(106)

where

$$\begin{aligned} \begin{array}{lcl} &{}&{} H_{11}=-\frac{1}{2}\nabla ^2+V_R+\int \frac{\big |\psi _-(y)\big |^2}{\big |x-y\big |}\,\mathrm {d}y-\frac{20}{9}\alpha \big |\psi _+(x)\big |^{2/3}-\epsilon _+,\\ &{}&{} H_{22}=-\frac{1}{2}\nabla ^2+V_R+\int \frac{\big |\psi _+(y)\big |^2}{\big |x-y\big |}\,\mathrm {d}y-\frac{20}{9}\alpha (\big |\psi _-(x)\big |^{2/3}-\epsilon _-. \\ \end{array} \end{aligned}$$

(107)

In addition, the eigenfunctions $w_i(x)$ satisfy the orthogonality relations

$$\begin{aligned} \begin{array}{lcl} \int \psi _i(x)w_i(x) \,\mathrm {d}x=0, i=+,-. \end{array} \end{aligned}$$

(108)

If all the eigenvalues of ${\textbf {Hess}}$ are positive, then there is no descent direction in the function space. Negative eigenvalues imply that there is a descent direction. To carry out this calculation, the integrals of the Coulomb potential required in (106) are obtained by solving the Poisson equation as in (83), and performing the numerical integrals. The calculated eigenvalues are shown in Table 1.