(Spin-)density-functional theory for open-shell systems: exact magnetization density functional for the half-filled Hubbard trimer

TL;DR

This paper derives the exact magnetization density functional for the half-filled Hubbard trimer, clarifying the relationship between DFT and spin-DFT for open-shell systems and addressing errors in common approximations.

Contribution

It provides the first explicit form of the exact magnetization functional for a model open-shell system, bridging a gap in density-functional theory.

Findings

Exact magnetization functional obtained for the Hubbard trimer

Analysis of errors in spin-restricted and unrestricted calculations

Partial correction of these errors using the exact functional

Abstract

According to the Hohenberg-Kohn theorem of density-functional theory (DFT), all observable quantities of systems of interacting electrons can be expressed as functionals of the ground-state density. This includes, in principle, the spin polarization (magnetization) of open-shell systems; the explicit form of the magnetization as a functional of the total density is, however, unknown. In practice, open-shell systems are always treated with spin-DFT, where the basic variables are the spin densities. Here, the relation between DFT and spin-DFT for open-shell systems is illustrated and the exact magnetization density functional is obtained for the half-filled Hubbard trimer. Errors arising from spin-restricted and -unrestricted exact-exchange Kohn-Sham calculations are analyzed and partially cured via the exact magnetization functional.

Equations

$$n_{\sigma }(\mathbf{r})=⟨\Psi |{\hat{\psi }}_{\sigma }^{†}(\mathbf{r}){\hat{\psi }}_{\sigma }(\mathbf{r}|\Psi ⟩,$$

$$\left[-\frac{{\nabla }^2}{2}+V_{\sigma }(\mathbf{r})+\int \frac{n({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }+V_{\mathrm{xc},\sigma }[n_{\uparrow },n_{\downarrow }](\mathbf{r})\right]{\phi }_{j\sigma }^u(\mathbf{r})$$

$$n_{\sigma }(\mathbf{r})=\sum _{j=1}^{N_{\sigma }}|{\phi }_{j\sigma }^u(\mathbf{r}){|}^2,$$

$$n(\mathbf{r})=n_{\uparrow }(\mathbf{r})+n_{\downarrow }(\mathbf{r}),$$

$$m_z(\mathbf{r})=n_{\uparrow }(\mathbf{r})-n_{\downarrow }(\mathbf{r}).$$

$$\left[-\frac{{\nabla }^2}{2}+V(\mathbf{r})+\int \frac{n({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d{\mathbf{r}}^{\prime }+V_{\mathrm{xc}}[n](\mathbf{r})\right]{\phi }_{j\sigma }^r(\mathbf{r})$$

$$n(\mathbf{r})=2\sum _{j=1}^{N/2}|{\phi }_j(\mathbf{r}){|}^2.$$

$$n(\mathbf{r})=2\sum _{j=1}^{(N-1)/2}|{\phi }_j(\mathbf{r}){|}^2+|{\phi }_{(N+1)/2}(\mathbf{r}){|}^2.$$

$$m_z^r(\mathbf{r})=\pm |{\phi }_{(N+1)/2}(\mathbf{r}){|}^2,$$

$$m_z^r(\mathbf{r})\ne m_z(\mathbf{r}).$$

$$\hat{H}$$

$${\hat{H}}^{\mathrm{KS}}=-t\sum _{⟨i,j⟩}\sum _{\sigma }({\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma }+{\hat{c}}_{j\sigma }^{†}{\hat{c}}_{i\sigma })+\sum _j\sum _{\sigma }{\hat{V}}_{j\sigma }^{\mathrm{KS}}{c}_{j\sigma }^{†}{\hat{c}}_{j\sigma },$$

$$V_j^{\mathrm{H}}=U{n}_j,$$

$$V_{j\sigma }^{\mathrm{XX}}=-U{n}_{j\sigma }.$$

$$V_j^{\mathrm{XXr}}=-U{n}_j/2.$$

$$\left\{\sum _{i=1}^3[T(i)+V(i)+\mathbf{B}(i)\cdot \mathbf{\sigma }(i)]+\sum _{i<j}^{3}W(i,j)\right\}{\Psi }_j(1,2,3)$$

$$=\hat{H}(1,2,3){\Psi }_j(1,2,3)$$

$${\Psi }_j(1,2,3)=\sum _n{c}_{nj}{\psi }_n(1,2,3).$$

$$\mathbf{H}{\vec{c}}_j=E_j{\vec{c}}_j,$$

$$\varphi_{1}=\left(\begin{array}[]{c}1\\ 0\\ 0\end{array}\right)\;,\hskip 14.22636pt\varphi_{2}=\left(\begin{array}[]{c}0\\ 1\\ 0\end{array}\right)\;,\hskip 14.22636pt\varphi_{3}=\left(\begin{array}[]{c}0\\ 0\\ 1\end{array}\right)\;.$$

$$\psi_{1}(1,2,3)=\left|\begin{array}[]{ccc}\varphi_{1}(1)\alpha(1)&\varphi_{1}(2)\alpha(2)&\varphi_{1}(3)\alpha(3)\\ \varphi_{2}(1)\beta(1)&\varphi_{2}(2)\beta(2)&\varphi_{2}(3)\beta(3)\\ \varphi_{3}(1)\beta(1)&\varphi_{3}(2)\beta(2)&\varphi_{3}(3)\beta(3)\end{array}\right|,$$

$$V(i)=\left(\begin{array}[]{ccc}V_{1}&0&0\\ 0&V_{2}&0\\ 0&0&V_{3}\end{array}\right)_{i}\sigma_{0}(i)\;,$$

$$⟨\hat{V}{⟩}_{1,1}=\dots =⟨\hat{V}{⟩}_{6,6}$$

$$⟨\hat{V}{⟩}_{7,7}=⟨\hat{V}{⟩}_{13,13}$$

$$⟨\hat{V}{⟩}_{8,8}=⟨\hat{V}{⟩}_{14,14}$$

$$⟨\hat{V}{⟩}_{9,9}=⟨\hat{V}{⟩}_{15,15}$$

$$⟨\hat{V}{⟩}_{10,10}=⟨\hat{V}{⟩}_{16,16}$$

$$⟨\hat{V}{⟩}_{11,11}=⟨\hat{V}{⟩}_{17,17}$$

$$⟨\hat{V}{⟩}_{12,12}=⟨\hat{V}{⟩}_{18,18}$$

$$⟨\hat{V}{⟩}_{19,19}=⟨\hat{V}{⟩}_{20,20}$$

$$T(i)=\left(\begin{array}[]{ccc}0&-t&-t_{p}\\ -t&0&t\\ -t_{p}&-t&0\end{array}\right)_{i}\sigma_{0}(i),$$

$$⟨\hat{T}{⟩}_{1,7}=⟨\hat{T}{⟩}_{1,12}=t_p⟨\hat{T}{⟩}_{1,8}=⟨\hat{T}{⟩}_{1,10}=t$$

$$⟨\hat{T}{⟩}_{2,8}=⟨\hat{T}{⟩}_{2,10}=⟨\hat{T}{⟩}_{2,11}=-t⟨\hat{T}{⟩}_{2,9}=t$$

$$⟨\hat{T}{⟩}_{3,7}=⟨\hat{T}{⟩}_{3,12}=-t_p⟨\hat{T}{⟩}_{3,9}=-t⟨\hat{T}{⟩}_{3,11}=t$$

$$⟨\hat{T}{⟩}_{4,13}=⟨\hat{T}{⟩}_{4,18}=-t_p⟨\hat{T}{⟩}_{4,14}=⟨\hat{T}{⟩}_{4,16}=-t$$

$$⟨\hat{T}{⟩}_{5,14}=⟨\hat{T}{⟩}_{5,16}=⟨\hat{T}{⟩}_{5,17}=t⟨\hat{T}{⟩}_{5,15}=-t$$

Full text

I Introduction: DFT versus SDFT

Spin-density-functional theory (SDFT) Parr and Yang (1989); Dreizler and Gross (1990); Engel and Dreizler (2011) is concerned with interacting $N$-electron systems in the presence of static magnetic fields. The magnetic fields are assumed to act only on the spins of the electrons (and not on the orbital currents), giving rise to a Zeeman-like term in the many-body Hamiltonian. The majority of applications of SDFT are for situations of collinear magnetism; spin then becomes a good quantum number associated with a fixed spin quantization axis, which is usually chosen to be the $z$-axis. In this case, the basic variables of SDFT are the spin densities,

[TABLE]

where $\Psi$ is the many-body wave function, and $\hat{\psi}^{\dagger}_{\sigma}({\bf r})$ and $\hat{\psi}_{\sigma}({\bf r})$ are creation and annihilation operators for Fermions of spin $\sigma=\uparrow,\downarrow$.

The basic theorem of SDFT establishes that in an $N$-electron system under the influence of a scalar potential $V({\bf r})$ and a magnetic field along the $z$ direction, $B_{z}({\bf r})$, the many-body Hamiltonian and all quantities that follow from it are functionals of the ground-state spin densities $n_{\sigma}({\bf r})$ von Barth and Hedin (1972); Gunnarsson and Lundqvist (1976); Gidopoulos (2007). In practice, the spin densities are obtained by solving the following Kohn-Sham equation:

[TABLE]

using atomic units with $e=m=\hbar=1$. Here, the spin-dependent external potential is defined as $V_{\uparrow,\downarrow}({\bf r})=V({\bf r})\pm\mu_{B}B_{z}({\bf r})$, where $\mu_{B}$ is the Bohr magneton, and the exchange-correlation (xc) potential $V_{{\rm xc},\sigma}[n_{\uparrow},n_{\downarrow}]$ is a functional of the spin densities.

Equation (2) is known as spin-unrestricted Kohn-Sham equation Jacob and Reiher (2012) since, in general, the $\varphi^{u}_{j\uparrow}({\bf r})$ and $\varphi^{u}_{j\downarrow}({\bf r})$ can be different. From the Kohn-Sham spin orbitals one obtains, in principle, the exact spin densities,

[TABLE]

where $N_{\sigma}$ is the number of electrons with spin $\sigma$, and $N_{\uparrow}+N_{\downarrow}=N$. The total particle density then follows as

[TABLE]

and the magnetization density along the $z$-direction is given by

[TABLE]

Notice that the definition $m_{z}({\bf r})=-\mu_{B}[n_{\uparrow}({\bf r})-n_{\downarrow}({\bf r})]$ is more frequently used in the literature. Here, however, we find it convenient, following Engel and Dreizler Engel and Dreizler (2011), to omit the minus sign and the Bohr magneton $\mu_{B}$ in the definition of $m_{z}$.

The SDFT formalism outlined above is suitable for systems of electrons in magnetic fields, but, more often, it is applied to situations where external magnetic fields are absent and the system is spontaneously magnetic. This happens, for instance, in open-shell atoms and molecules with an odd number of electrons. In that case, $V_{\uparrow}=V_{\downarrow}=V$, but in general $V_{\rm xc,\uparrow}\neq V_{\rm xc,\downarrow}$. Most modern approximate xc functionals are therefore constructed in terms of the spin densities and/or the (unrestricted) spin orbitals Perdew and Kurth (2003); Cohen et al. (2012); Becke (2014).

From a strictly formal point of view, using SDFT in systems without external magnetic field is unnecessary. According to the Hohenberg-Kohn theorem of DFT Hohenberg and Kohn (1964); Kohn and Sham (1965), the Hamiltonian and all observables in such an $N$-electron systems are functionals of the ground-state density $n({\bf r})$ alone, which can be obtained from the spin-restricted Kohn-Sham equation

[TABLE]

where $\varphi^{r}_{j\uparrow}({\bf r})=\varphi^{r}_{j\downarrow}({\bf r})\equiv\varphi_{j}({\bf r})$ and $\varepsilon^{r}_{j\uparrow}=\varepsilon^{r}_{j\downarrow}\equiv\varepsilon_{j}$.

Let us first consider the straightforward case where the system is finite with an even number of electrons, which is the case for closed-shell configurations. Since then $N_{\uparrow}=N_{\downarrow}=N/2$, the ground-state density is given by a sum over doubly occupied orbitals:

[TABLE]

If, however, the number of electrons $N$ is odd, then the density must be computed as

[TABLE]

In other words, the first $N-1$ electrons are in doubly occupied orbitals, and the last ($N^{\rm th}$) electron is in a singly occupied orbital. The spin of individual Kohn-Sham electrons is not explicitly referenced; all electrons experience the same Kohn-Sham effective potential. If the exact xc potential $V_{\rm xc}[n]({\bf r})$ is used in Eq. (6), then the density $n({\bf r})$ obtained from Eqs. (7) (for even $N$) or (8) (for odd $N$) will be exact, and equal to the density obtained via SDFT by using Eq. (4).

Within spin-restricted Kohn-Sham, and for odd $N$, the magnetization along $z$ is given by

[TABLE]

where the sign indicates whether the highest singly occupied orbital is spin-up $(+)$ or spin-down $(-)$. In general, the restricted Kohn-Sham magnetization is not equal to the exact magnetization, even if the exact $V_{\rm xc}[n]({\bf r})$ is used in Eq. (6):

[TABLE]

On the other hand, since in DFT the many-body wave function $\Psi[n]$ is a functional of the density, the exact magnetization can also be formally expressed as a functional of the density, $m_{z}[n]({\bf r})$. This means that, in principle, the spin-restricted Kohn-Sham formalism of DFT is sufficient to calculate the magnetization exactly (again, no external magnetic field is present). However, the form of this density functional $m_{z}[n]({\bf r})$ is unknown. Therefore, treating open-shell, spin-polarized systems with DFT rather than SDFT is generally discouraged, unless there are specific reasons such as avoiding the so-called spin contamination problem Jacob and Reiher (2012); Perdew et al. (2009).

These differences between DFT and SDFT are well known (see also the recent discussion in Ref. Trushin and Görling (2018)). Nevertheless, there seems to be no example in the literature where the magnetization density functional $m_{z}[n]({\bf r})$ is explicitly constructed and compared with the restricted spin polarization $m_{z}^{r}({\bf r})$. In this paper, we will present such a case study, using the half-filled Hubbard trimer as model system. By comparing exact results with approximate (exchange-only) Kohn-Sham calculations, errors of the magnetization will be analyzed and shown to arise from different sources, depending on whether the Kohn-Sham scheme is restricted or unrestricted. The exact magnetization functional will be shown to provide a cure to some of the errors, but not all of them.

II The half-filled Hubbard trimer

II.1 Model

The inhomogeneous Hubbard model is defined by the following Hamiltonian Capelle and Campo Jr. (2013):

[TABLE]

where $\hat{c}_{j\sigma}^{\dagger}$ and $\hat{c}_{j\sigma}$ are creation and annihilation operators, respectively, for electrons with spin $\sigma$ on site $j$, and $\langle i,j\rangle$ indicates pairs of nearest-neighbor lattice sites. We fix the hopping parameter as $t=0.5$. The on-site interaction strength $U$ and the external potential $V_{j}$ will be treated as variable parameters.

The corresponding Kohn-Sham Hamiltonian is

[TABLE]

where the Kohn-Sham potential at lattice site $j$ is given by $V_{j\sigma}^{\rm KS}=V_{j}+V_{j}^{\rm H}+V_{j\sigma}^{\rm xc}$. The Hartree potential is simply

[TABLE]

where $n_{j}$ is the density (site occupation) on the $j^{\rm th}$ lattice point. In the spin-restricted Kohn-Sham scheme, the xc potential is constrained to satisfy $V_{j\uparrow}^{\rm xc}=V_{j\downarrow}^{\rm xc}$.

For the Hubbard model with on-site interactions, the simplest approximation for the xc potential is exact exchange (XX), given by Ullrich (2018)

[TABLE]

The spin-restricted exact exchange (XXr) potential is

[TABLE]

To include correlation effects, one could use the Bethe-ansatz local density approximation (BALDA) Capelle and Campo Jr. (2013); Lima et al. (2002, 2003); França et al. (2012). However, the BALDA runs into convergence problems on small lattices whenever any of the site occupations approaches 1 Ullrich (2018). Alternatively, correlation can be treated via orbital functionals, as discussed elsewhere Ullrich (2018); Pluhar III and Ullrich . In this paper, however, correlation effects are not included to keep things simple.

The Hubbard model, especially the half-filled Hubbard dimer, has been widely used as a test system for DFT Carrascal et al. (2015); Cohen and Mori-Sánchez (2016); Carrascal et al. (2018); França et al. (2018); Herrera et al. (2018); Ullrich (2018). In the absence of an external magnetic field, the ground state of the Hubbard dimer is a singlet, with zero total spin and no magnetization. For the purpose of comparing DFT and SDFT, the simplest suitable system is the half-filled Hubbard trimer, that is, three electrons on a three-point lattice foo . The half-filled Hubbard trimer has been discussed in the literature Shiba and Pincus (1972); Trif et al. (2010); Nossa et al. (2012); Ghassemi Tabrizi et al. (2019), but not in the context of DFT. Most notably, earlier studies of the Hubbard trimer were limited to the homogeneous case of constant potential. Here, by contrast, we consider the inhomogeneous case.

II.2 Example: from weakly to strongly correlated

II.2.1 Ground-state energies and densities

To illustrate the half-filled Hubbard trimer in the weakly and strongly correlated limits, we now consider an asymmetric example where we fix the external potential as $V_{1}=1$, $V_{2}=0$, and $V_{3}=2$ and vary the interaction strength from $U=0$ to 10. We obtain the ground state of the interacting system (in both linear and triangular configuration) through exact diagonalization; technical details are given in the Appendix. This yields the exact ground-state energy $E$, density $n_{j}$ and magnetization $m_{zj}$ for each given $U$.

We compare the exact solutions with restricted and unrestricted Kohn-Sham calculations using XX and XXr, respectively. Figure 1 shows the ground-state energies. The results are very similar for the linear and the triangular lattice: in both cases, XX agrees closely with the exact ground-state energy, approaching a final limit of $E=3$ as the electrons localize after a clear crossover from weakly to strongly correlated around $U=2$. The spin-restricted XXr Kohn-Sham calculation, on the other hand, produces a continuously increasing energy. This familiar behavior was observed earlier for the case of the half-filled (spin singlet) Hubbard dimer Carrascal et al. (2015). There, however, the spin symmetry had to be artificially broken to get better agreement with the exact results.

The ground-state densities are shown in the top panels of Fig. 2 for the linear trimer and of Fig. 3 for the triangular lattice. For XX and XXr, the density follows from Eqs. (4) and (8), respectively. The exact solution shows that the density behaves very similarly in the two cases: for small $U$, the density is much larger on the second lattice point (where the potential is $V_{2}=0$) compared to the two other lattice points (where $V_{1}=1$ and $V_{3}=2$). As $U$ increases, the electrons begin to localize; eventually, for $U\sim 10$, all lattice points are equally populated ($n_{1}=n_{2}=n_{3}=1$). The crossover between the delocalized and localized regimes is not very sharp, but occurs in a region roughly between $U=1$ and $3$.

This behavior of the density is very well reproduced by the spin-unrestricted XX calculations. On the other hand, XXr agrees less well with the exact density: on both lattices, the transition from delocalized to localized (where all $n_{j}\to 1$) is too slow.

II.2.2 Magnetizations

The results for the magnetization are shown in the middle and bottom panels of Figs. 2 and 3. For XX, the magnetization is obtained from the approximate spin densities without further approximation, using Eq. (5); in the case of XXr, the magnetization is obtained using Eq. (9), which is an approximation, even if the density and the (restricted) Kohn-Sham orbitals were exact.

The middle panels of Figs. 2 and 3 show the exact $m_{z}$, compared with unrestricted XX. Again, the linear and triangular lattices behave in a similar manner. The exact magnetization is such that, for small $U$, the first lattice point is almost fully magnetized, whereas the second and third lattice point have almost no magnetization. As $U$ increases, both lattices develop a slightly antiferromagnetic pattern: for the linear lattice, $m_{z1}$ and $m_{z3}$ approach $2/3$ and $m_{z2}$ approaches $-1/3$; for the triangular lattice the roles of $m_{z2}$ and $m_{z3}$ are reversed. It is interesting to observe that the localization of the density implies equal site occupation, but not equal magnetization, even in the case of the triangular lattice with periodic boundary conditions.

This behavior is reasonably well reproduced by the spin-unrestricted XX calculations. The agreement is excellent for both lattices up until $U\sim 2$, but for stronger interactions the antiferromagnetism is drastically overestimated: on two lattice points the magnetization approaches 1 and on the third point it approaches $-1$. This demonstrates that an SDFT calculation can produce excellent total densities and energies, but poor spin densities. This well-known phenomenon has been observed, for example, in transition-metal complexes Jacob and Reiher (2012); Conradie and Ghosh (2007); Boguslawski et al. (2011); Song et al. (2018).

The bottom panels of Figs. 2 and 3 compare the exact restricted magnetization $m_{z}^{r}$ with XXr. The exact $m_{z}^{r}$ is obtained by first calculating the exact density for a given $U$, then determining that potential which, for $U=0$, reproduces this density (this is the exact restricted $V_{j}^{\rm KS}$), and then calculating the restricted magnetization from this (see Section II.3 for further details of the inversion procedure).

Comparing $m_{z}^{r}$ with the exact $m_{z}$ reveals significant differences for moderate and high interaction strengths: mainly, the antiferromagnetic trend is not reproduced. For the linear lattice, $m_{z1}^{r}$ and $m_{z3}^{r}$ approach $1/2$ instead of $2/3$, and $m_{z2}^{r}$ goes to zero instead of $-1/3$. The XXr calculations agree quite well with the exact $m_{z}^{r}$ for all interaction strengths.

On the triangular lattice, the exact $m_{z}^{r}$ approaches a completely symmetric strongly-correlated limit at large $U$, with $m_{zj}^{r}\to 1/3$ for all lattice points. The crossover region exhibits quite dramatic changes of the behavior of the magnetization, which has an interesting underlying cause: from $U=2.7$ onwards, the second and third Kohn-Sham level become degenerate, which requires determining the correct linear combination of the degenerate orbitals to reproduce the exact density via Eq. (8). The associated magnetization, Eq. (9), is strongly affected by this choice of orbital. A similar degeneracy occurs in the XXr calculations, but only around $U=7$, and with no dramatic consequences for the magnetization. Thus, the agreement between XXr and exact $m_{z}^{r}$ is good for small $U$, but XXr approaches an unsymmetric limit at large $U$, with $m_{z1}\to 1/3$, $m_{z2}\to 2/3$, and $m_{z3}\to 0$; incidentally, this happens to be closer to the exact unrestricted $m_{z}$ than to the exact $m_{z}^{r}$.

The examples discussed here are typical for the behavior of the inhomogeneous half-filled Hubbard trimer. In general, XX, both in its spin-restricted and -unrestricted form, performs well as long as the system is weakly correlated, but fails for the magnetization at strong correlation. However, the failures of the restricted and unrestricted XX are of a qualitatively different nature, providing prime examples of the two different kinds of errors of (S)DFT, namely, those caused by approximations of the functionals, and those caused by using approximate densities (“density-driven errors”) Song et al. (2018); Kim et al. (2013); Sim et al. (2018). In the following, we will investigate these errors more closely, looking for ways in which they could possibly be mitigated.

II.3 The magnetization density functional

We now construct the functional $m_{z}[n]$ which produces the exact magnetization from a given ground-state density, for the half-filled Hubbard trimer with interaction $U$. To do this, we invert the interacting 3-body Schrödinger equation on the 3-point lattice (see Appendix), governed by the Hamiltonian (11). In other words, we construct the lattice potential $V_{j}$ that produces a given ground-state density distribution $n_{j}$ (for a given $U$), and from the solution of the Schrödinger equation with this potential we obtain the exact magnetization $m_{zj}$. This is repeated for all possible density distributions, which then defines the functional $m_{z}[n]$.

In practice, to determine the potential $V_{j}$ from the density, we minimize the mean square deviation of the associated ground-state density $n_{j}$ from a given target density $n^{t}_{j}$, using Powell’s conjugate direction method Press et al. (2007). This only requires a moderate computational effort for the systems under study.

The results are shown in Fig. 4 for the linear Hubbard trimer and in Fig. 5 for the triangle. In the graphic representation of $m_{zj}[n]$, $n_{1}$ and $n_{3}$ are chosen as the two independent values of the density; $n_{2}$ follows as $n_{2}=3-n_{1}-n_{3}$, and all densities (or site occupations) are constrained to satisfy $0\leq n_{j}\leq 2$. The domain of the functional is therefore a square area in $n_{1}-n_{3}$ space with two triangular regions in the upper right-hand and lower left-hand corners excluded.

Figures 4 and 5 illustrate the functional $m_{zj}[n]$ on the three lattice points $j=1,2,3$, and for five different interaction strengths, $U=0,1,2,5,10$. The color maps show that the range of $m_{zj}$ on the Hubbard lattices is between $-1/3$ and 1. The condition $m_{z1}+m_{z2}+m_{z3}=1$ is always satisfied. However, the detailed behavior of the magnetization is quite different for the linear chain (fixed boundary conditions) and the triangle (periodic boundary conditions).

For the noninteracting case, the magnetization lies between 0 and 1 on both lattices. As the interaction strength $U$ increases, the maps generally develop sharper features. In particular, pockets of negative magnetization begin to appear, marked by dark blue turning into black. The partial antiferromagnetism is thus a generic feature of the half-filled Hubbard trimer. For the linear lattice, only the magnetization on the central site can become negative; for the triangular lattice, it can be negative on either of the three sites.

II.4 Eliminating the functional error using the exact magnetization density functional

In Section II.2, we saw that the approximate Kohn-Sham calculations are plagued by different kinds of errors. XX generally leads to good total energies and densities, but the individual spin densities [which give rise to the magnetization via Eq. (5)] are not accurate for moderate and strong correlations. This is a clear example of a density-driven (or, more precisely, spin-density-driven) error, since the magnetization functional is exact in the unrestricted case. XXr, on the other hand, is not very accurate for total densities (unless $U$ is small), and produces magnetizations of poor quality. Hence, the errors are both density-driven and due to a poor approximation of the magnetization functional, but it is not clear a priori which is worse.

Let us now eliminate the error due to the wrong magnetization functional, and focus on the density-driven errors by themselves. Figure 6 shows the magnetizations of the linear and triangular lattices (for the same systems considered in Figs. 2 and 3) which follow from evaluating the exact functional $m_{zj}[n]$ with the self-consistent XX and XXr densities. The spin-unrestricted XX now agrees very well with the exact magnetization for all interaction strengths, reproducing the weakly antiferromagnetic behavior in both cases. The spin-density-driven error is thus clearly the dominant one in XX.

On the other hand, XXr, in spite of some improvement, is still not very close to the exact result, particularly for interaction strengths greater than $U=2$. Overall, the exact magnetization functional cannot make up for the deficiencies of the XXr total density. For better agreement, correlation effects must be included on top of XXr.

III Conclusions

In the standard SDFT approach to open-shell systems, the magnetization (or, equivalently, the spin polarization) is expressed in terms of the spin densities; the latter are obtained using a spin-unrestricted Kohn-Sham formalism, which uses the spin densities as two independent variables. However, according to the Hohenberg-Kohn theorem, all observable quantities follow in principle from the total density; an explicitly spin-polarized Kohn-Sham formalism, with two types of orbitals (spin-up and spin-down) is therefore, in principle, unnecessary. This could reduce the computational effort of such systems by half. Of course, the reason that SDFT is almost always preferred over DFT is that it is much easier to construct accurate approximations thanks to the flexibility that is afforded by having two independent variables. Nevertheless, the point of principle remains.

Here, we have constructed and analyzed the exact density functional of the magnetization, $m_{z}[n]$, for two simple model systems, the half-filled Hubbard trimer with fixed and periodic boundary conditions (chain and triangle). The functional is rather smooth for small interactions strengths $U$, but acquires more sharp features in the strongly correlated limit, where the system tends towards a partial antiferromagnetism.

Two approximations were tested: spin-unrestricted and spin-restricted exact exchange (XX and XXr). As expected, XX performs much better, giving very good total energies and densities, but fails to reproduce the correct magnetization for moderate and strong interactions. The error thus clearly lies in the individual spin densities, as demonstrated by plugging the XX total density into the exact magnetization functional. On the other hand, XXr fails as soon as the system is no longer weakly interacting, which leads to poor energies, densities and (restricted) magnetization. This is not surprising, based on what is known from other test cases such as the half-filled Hubbard dimer. The exact magnetization functional is unable to provide a cure for the shortcomings of XXr.

The case study discussed in this paper thus provides a potentially valuable proof of principle for (S)DFT. Many of today’s popular xc functionals are not optimized to produce the best spin densities, but focus more on total energies and total densities. This can be problematic if explicitly spin-dependent observables are of interest, such as electron paramagnetic resonance parameters Boguslawski et al. (2011). The construction of explicit (but approximate) magnetization functionals depending on total densities may therefore be a worthwhile effort. The example of the half-filled Hubbard trimer has shown that this can lead to quite significant improvements.

Acknowledgements.

This work was supported by Department of Energy Grant No. DE-SC0019109 and by a Research Corporation for Science Advancement Cottrell Scholar SEED Award.

Appendix A Exact diagonalization of the half-filled Hubbard trimer in a magnetic field

This Appendix gives technical details on the numerical solution of the half-filled Hubbard trimer in the presence of arbitrary scalar potentials and magnetic fields (acting as Zeeman terms on the spins only).

Consider three interacting electrons on a three-point lattice, with sites 1,2,3, arranged along a line (non-periodic boundary conditions) or as a triangle (periodic boundary conditions). The scalar potential on the three lattice sites is denoted by $V_{1,2,3}$, and the magnetic field is denoted by ${\bf B}_{1,2,3}$. The general three-particle Schrödinger equation reads:

[TABLE]

where $T(i)$ is the kinetic energy operator acting on the $i$th electron (here, nearest-neighbor hopping with amplitude $t$), $\bm{\sigma}$ is the vector of Pauli matrices, and $W(i,j)$ is the interaction between the $i$th and $j$th electron (here, on-site interaction with strength $U$).

The wave functions $\Psi_{j}(1,2,3)$ can be expanded in a basis of Slater determinants, $\{\psi_{n}(1,2,3)\}$:

[TABLE]

The expansion coefficients $c_{nj}$ and the energies $E_{j}$ follow from diagonalization of the Hamiltonian matrix $\bf H$:

[TABLE]

where the elements of the $20\times 20$ matrix $\bf H$ are given by $\langle\hat{H}\rangle_{m,n}=\langle\psi_{m}|\hat{H}|\psi_{n}\rangle$ [$\hat{H}$ is the full three-electron Hamiltonian of Eq. (A)]. The individual parts of the Hamiltonian matrix $\langle\hat{H}\rangle_{m,n}=\langle\hat{T}\rangle_{m,n}+\langle\hat{V}\rangle_{m,n}+\langle\hat{B}\rangle_{m,n}+\langle\hat{W}\rangle_{m,n}$ will now be derived.

For a three-point lattice, our basis consists of 20 Slater determinants, as schematically illustrated in Fig. 7. Each Slater determinant is an eigenfunction of $\hat{S}_{z}$ (spin along the $z$-axis), but not necessarily of $\hat{\bf S}^{2}$ (total spin): the first 18 Slater determinants, $\psi_{1}$ - $\psi_{18}$, have $S_{z}=\pm 1/2$, and the last two determinants, $\psi_{19}$ and $\psi_{20}$, have $S_{z}=\pm 3/2$. By suitable linear combination of the $\psi_{n}$, one can form 4 quartet wave functions (with $S=3/2$) and 16 doublet wave functions (with $S=1/2$). These could be used as an alternative basis to diagonalize $\hat{H}$, but the derivation of the associated Hamiltonian matrix would be more cumbersome.

The space parts of the Slater determinants are constructed using the single-particle basis functions

[TABLE]

The column vector notation refers to the three-site lattice, i.e., the components of $\varphi_{1,2,3}$ denote the value of the wave functions on lattice sites 1,2,3. The Slater determinants can then be explicitly written down; for instance,

[TABLE]

and similar for all other $\psi_{n}$, where $\alpha$ and $\beta$ are up- and down-spinors.

In this notation, the scalar potential operator acting on the $i$th electron is given by

[TABLE]

where the subscript $i$ means that the $3\times 3$ matrix acts on $\varphi_{1,2,3}(i)$, and $\sigma_{0}$ is the $2\times 2$ unit matrix (acting on the $\alpha$ and $\beta$ spinors of the $i$th electron). Here, $\hat{V}=V(1)+V(2)+V(3)$, and we define the other elements of the three-electron Hamiltonian $\hat{H}=\hat{T}+\hat{V}+\hat{\bf B}\cdot{\bm{\sigma}}+\hat{W}$ in an analogous manner.

It is easy to see that only the diagonal matrix elements of the scalar potential operator are nonzero, and we have

[TABLE]

The kinetic energy operator depends on the boundary conditions. In general, we can write

[TABLE]

where $t_{p}=0$ for the linear chain and $t_{p}=t$ for the triangle with periodic boundary conditions.

The only nonzero kinetic energy matrix elements are those that correspond to a single hopping event, without spin flip, namely:

[TABLE]

The matrix elements in the lower triangle follow from the hermiticity of the Hamiltonian matrix, i.e., $\langle\hat{T}\rangle_{m,n}=\langle\hat{T}\rangle_{n,m}$.

The only nonvanishing matrix elements of the on-site interaction are the diagonal elements

[TABLE]

The Zeeman terms in the Hamiltonian are given by

[TABLE]

We define $\hat{B}_{x}=B_{x}\sigma_{x}$ and similar for $y$ and $z$.

For the $z$-components, only the diagonal matrix elements are nonzero, and we get immediately

[TABLE]

For the $x$ and $y$-components, we collect all the spin flips. Defining $B_{\pm j}=B_{xj}\pm iB_{yj}$ for the $j$th lattice site, we get

[TABLE]

The matrix elements in the lower triangle follow directly as $\langle\hat{B}_{x}+\hat{B}_{y}\rangle_{m,n}=\langle\hat{B}_{x}+\hat{B}_{y}\rangle_{n,m}^{*}$.

From the eigenvectors of Eq. (18), we obtain the spin-density matrix from the following formula:

[TABLE]

The final result is

[TABLE]

[TABLE]

[TABLE]

From this, one obtains the density as $n=n_{\uparrow\uparrow}+n_{\uparrow\downarrow}$ and the $x-$, $y-$ and $z$-components of the magnetization as $m_{x}=n_{\uparrow\downarrow}+n_{\downarrow\uparrow}$, $m_{y}=i(n_{\uparrow\downarrow}-n_{\downarrow\uparrow})$, and $m_{z}=n_{\uparrow\uparrow}-n_{\downarrow\downarrow}$.