# Prediction of Magnetic State of UO2 within Hubbard-corrected Density-
Functional Theory: A self-consistent approach

Mahmoud Payami<EMAIL_ADDRESS>School of Physics & Accelerators, Nuclear
Science and Technology Research Institute,
AEOI, P. O. Box 14395-836, Tehran, Iran

###### Abstract

The magnetic state of UO2 was determined experimentally to be anti-
ferromagnetic. Starting from this experimental fact, researchers have
calculated other properties within the Hubbard-corrected density-functional
theory, DFT+U. Up to now, the Hubbard parameters for UO2 were usually so
chosen that the calculations give good results for some experimental data.
Also, to our knowledge there exists no valid theoretical research report on
the energetically stable magnetic state of this system. In present work,
employing the new method which is based on density-functional perturbation
theory, we have determined self-consistently the Hubbard parameters and
ground-state energies for UO2 crystal in both ferromagnetic and anti-
ferromagnetic configurations, and the calculated results show that UO2 crystal
energetically favors an anti-ferromagnetic state with a small energy
difference. In all the calculations the PBE-sol approximation was used for the
exchange-correlation energy functional.

Uranium dioxide; Ferromagnetism; Anti-ferromagnetism; Density-Functional
Theory; Hubbard Model; Mott Insulator; DFT+U.

##  I Introduction

The popular local-density approximation (LDA) [1, 2] and generalized gradient
approximation (GGA) [3] approximations for the exchange-correlation (XC)
energy functional in the density-functional theory (DFT) [4, 1] suffer from
self-interaction errors, which are significant in systems containing atoms
with localized d and f orbitals, and by over-delocalization of their
corresponding wave-functions lead to incorrect prediction of metallic behavior
for Mott insulators. A simple and low-cost workaround is using the Hubbard-
corrected DFT model to correct the correlation energy of localized orbitals in
DFT energy-functional, DFT+U, in which only on-site corrections are added to
the DFT energy functional [5]:

$E_{\rm DFT+U}=E_{\rm DFT}[n({\bf r})]+E_{\rm Hub}[{n_{m}^{I\sigma}}]-E_{\rm
dc}[n^{I\sigma}],$ (1)

where $n({\bf r})$ is electron density, $n_{m}^{I\sigma}$ are orbital
occupation numbers of atom at lattice site ${\bf R}_{I}$, and
$n^{I\sigma}=\sum_{m}n_{m}^{I\sigma}$. The last term in right hand side of Eq.
(1) is needed to avoid double counting of interactions contained in the first
and second terms. The rotationally invariant form [6] of the correction is
given by [5]:

$E_{\rm U}[{n_{mm^{\prime}}^{I\sigma}}]\equiv E_{\rm Hub}-E_{\rm
dc}=\sum_{I,\sigma}\frac{U^{I}}{2}{\rm Tr[{\bf n}^{\it I\sigma}({\bf 1}-{\bf
n}^{\it I\sigma}})],$ (2)

in which ${\bf n}^{I\sigma}$ is the atomic occupation matrix. This on-site
correction, significantly improves the over-delocalization and lead to correct
insulating properties. The coefficients $U^{I}$ are called Hubbard on-site
parameters, and for a known material, these $U^{I}$ values may be so adjusted
that the calculations results well agrees with some experimental data [7].
However, in the case of designing new materials there is no experimental data
to be used for parameters fitting and on the other hand, it is very important
for the theory describing a material to be a parameter-free one. The first
attempts in this way was using linear-response constrained-DFT (LR-cDFT)
within super-cell method [5, 8] which was somewhat inconvenient and
computationally demanding. Another scheme to estimate the parameters was named
as constrained random-phase approximation (cRPA) [9] which was recently used
by others [10] to estimate the U parameter for uranium dioxide. However, a new
method based on density-functional perturbation theory (DFPT) was recently
introduced which instead of using previous super-cells, focuses on the unit-
cell, which is more convenient and relatively fast [11, 12, 13]. Using this
new method, we have determined self-consistently the Hubbard parameters for
UO2 crystal in both ferromagnetic (FM) and ant-ferromagnetic (AFM)
configurations of uranium atoms and showed that energetically UO2 crystal
favors an AFM magnetic state with a small energy difference, which is in
agreement with experimental findings. In our recent work [14] we have
demonstrated that choosing PBE-sol approximation [15] for the XC functional
leads to excellent results from self-consistent DFT+U calculations, and
therefore we have employed it in all computations throughout this work.

The organization of this paper is as follows: Section II is dedicated to the
computational details; in Section III we present and discuss the calculated
results; and finally section IV summarizes and concludes this research.

##  II Computational details

The crystal structure of uranium dioxide is described by a simple tetragonal
lattice with a six-atoms basis, shown in Fig. 1. In FM configuration, all
spins of U atoms have the same direction along $z$-axis, while to setup AFM
structure for U atoms, we use the simple model in which the planes of U atoms
alternate their spins along $z$ direction, i.e., we assume a 1-dimensional
AFM.

Figure 1: UO2 crystal structure as simple tetragonal with six atoms basis.
Left and right figures schematically represent FM and AFM configurations,
respectively. Large grey and small red balls represent uranium and oxygen
atoms, respectively. The up-spin and down-spin atoms are shown with yellow and
green colors, respectively.

The Hubbard on-site U parameters were calculated self-consistently using the
HP code [13] included in the Quantum-ESPRESSO code package [16, 17] for both
FM and AFM configurations of U atoms (shown in Fig. 1) in the context of PBE-
sol [15] approximation to the XC. The calculations include results for both
”atomic” and ”ortho-atomic” types of projections onto Hubbard orbitals. The
electronic structure calculations are based on the solution of the KS
equations using the Quantum-ESPRESSO code package. For U and O atoms the
scalar-relativistic ultra-soft pseudo-potentials (USPP) were used which were
generated by the atomic code and generation inputs from the pslibrary [18], at
https://github.com/dalcorso/pslibrary. The valence configurations
U($6s^{2},\,6p^{6},\,7s^{2},\,7p^{0},\,6d^{1},\,5f^{3}$) and
O($2s^{2},\,2p^{4}$) were adopted in the USPP generation. Kinetic energy
cutoffs for the plane-wave expansions were chosen as 90 and 720 Ry for the
wave-functions and densities, respectively. The smearing method of Marzari-
Vanderbilt [19] for the occupations with a width of 0.01 Ry were used. For the
Brillouin-zone integrations in geometry optimizations, a $8\times 8\times 6$
grid were used; All geometries were fully optimized for total residual
pressures on unit cells to within 0.5 kbar, and residual forces on atoms to
within 10-3 mRy/a.u. To self-consistent determination of the Hubbard
parameters we have employed the HP code [13] following the flowchart shown in
Fig. 2.

Figure 2: Flowchart of SCF determination of Hubbard parameters. In the first-
step SCF and last-step structure-optimization, the meta-stable states were
avoided [20, 21]

To start the self-consistent procedure for determining the Hubbard parameters
according to Fig. 2, the convergent q-mesh of $4\times 4\times 3$ for linear
response calculations was adopted and we gave initial values for $U_{in}$ in
each of FM and AFM configurations; for the initial structure we chose simple
tetragonal structure with appropriate lattice constants consistent with cubic
structure of side 5.47$\AA$. To avoid meta-stable states, we determined
appropriate occupations of Hubbard orbitals $5f$ of uranium atoms [20, 21]. In
this second step, we start the DFPT calculation and obtain new values for
parameters named as $U_{out}$. In the third step, using the parameter
$U_{out}$ we obtained in the second step, the geometry of the system is
optimized taking care of meta-stable states. In each cycle the differences
between input and output parameters were monitored to see if the self-
consistency is reached within $\Delta$ value. For this system the self-
consistency was reached within 6 to 8 cycles in the flowchart with
$\Delta<10^{-4}$.

##  III Results and discussions

The calculations were performed at the level of DFT with on-site Hubbard
corrections (DFT+U) for both FM and AFM configurations, and the results showed
that both FM and AFM magnetic states are insulators with different values for
electronic band gaps. The self-consistent Hubbard parameters and geometric
equilibrium lattice constants for FM and AFM states differ slightly; but their
respective band gaps show relatively significant differences. The results are
summarized in Table 1.

Table 1: Self-consistent Hubbard parameters $U_{sc}$ in eV; total energy with respect to the lowest one corresponding to ortho-atomic AFM, per formula unit $\Delta E$ in eV; total and absolute magnetizations per formula unit in Bohr-magneton; equilibrium lattice constants in $\AA$; and electronic band-gaps in eV for different cases studied. PBE-sol approximation was used for the XC, which turns out to be the best one. | AFM |  |  FM
---|---|---|---
U-proj | $U_{sc}\;(eV)$ | $\Delta E$ | $M_{tot}$ | $M_{abs}$ | $a\;(c)\;(\AA)$ | $E_{g}\;(eV)$ |  | $U_{sc}\;(eV)$ | $\Delta E$ | $M_{tot}$ | $M_{abs}$ | $a\;(c)\;(\AA)$ | $E_{g}\;(eV)$
o-atomic | 2.9787 | 0.0000 | 0.00 | 2.150 | 5.4540(5.4631) | 2.2281 |  | 2.9936 | 0.0103 | 2.00 | 2.175 | 5.4541(5.4655) | 1.7102
atomic | 2.0862 | 0.0961 | 0.00 | 2.155 | 5.4546(5.4766) | 1.3699 |  | 2.1247 | 0.1140 | 2.00 | 2.190 | 5.4582(5.4770) | 0.9162

As is seen from Table 1, the AFM configuration with orthogonalized projections
on Hubbard orbitals is the lowest energy state, and therefore the
energetically stable state of UO2 is an AFM configuration. To our knowledge,
this result is obtained for the first time. The lattice constants differ by
$0.02\AA$ and all are in excellent agreement with experiment. The band gaps
for AFM are larger than corresponding values in FM configurations. It is
important to note that the self-consistent Hubbard parameters are different
for AFM and FM configurations which result from different responses of them to
an external perturbation. Therefore, applying the same empirical value for AFM
and FM states of a system and comparing the energies to decide the favored
state is shown to be incorrect. As is seen, the band gap of AFM with
orthogonalized projections gives the best agreement with experiment [22].

To compare the overall electronic states of FM and AFM configurations, we have
plotted the electronic total density of states (DOS) and shown in Fig. 3.
Inspecting the Fig. 3, one notices that the behaviors near the band edges are
different for two magnetic states, and the value of band gap for FM is smaller
than that of AFM.

Figure 3: Total electronic density of states (DOS) for AFM and FM
configurations of UO2 crystal. It is observed that in a given magnetic state
of UO2, the band gap is larger for orthogonalized projections. Moreover, the
band gaps of AFM are larger than corresponding values for FM.

##  IV Conclusions

In the study of strongly-correlated UO2 system, usually one takes an AFM
configuration for the uranium atoms which is borrowed from experimental
findings, and the best job afterwards is to calculate the Hubbard parameter U
self-consistently, as done by present authors in earlier work. However, it is
very much interesting to predict the magnetic state theoretically as well
without resorting to experimental facts. This becomes vital in designing new
novel materials. In present work, using DFPT, we have calculated self-
consistently the Hubbard on-site parameters and thereof the total energies for
the two AFM and FM states of UO2 crystal and have shown that the resulting
parameters are different and AFM configuration is energetically favored by a
small value of $0.01eV$/(formula unit).

## Acknowledgement

This work is part of research program in School of Physics and Accelerators,
NSTRI, AEOI.

## Data availability

The raw or processed data required to reproduce these results can be shared
with anybody interested upon sending an email to M. Payami.

## References

  * [1] W. Kohn and L. J. Sham, ”Self-consistent equations including exchange and correlation effects”, Phys. Rev. 140, A1133 (1965).
  * [2] J. P. Perdew and A. Zunger, ”Self-interaction correction to density-functional approximations for many-electron systems”, Phys. Rev. B 23, 5048 (1981).
  * [3] J. P. Perdew, K. Burke, and M. Ernzerhof, ”Generalized gradient approximation made simple”, Phys. Rev. Lett. 77, 3865 (1996); Erratum: Phys. Rev. Lett. 78, 1396 (1997).
  * [4] P. Hohenberg and W. Kohn, ”Inhomogeneous electron gas”, Phys. Rev. 136, B864 (1964).
  * [5] M. Cococcioni and S. de Gironcoli, ”Linear response approach to the calculation of the effective interaction parameters in the LDA+ U method”, Phys. Rev. B 71, 035105 (2005).
  * [6] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, ”Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+ U study”, Phys. Rev. B57, 1505 (1998).
  * [7] M. Payami, ”DFT+U study of UO2: Correct lattice parameter and electronic band-gap”, arXiv:2302.13381 [cond-mat.str-el], https://doi.org/10.48550/arXiv.2302.13381.
  * [8] V. L. Campo Jr and M. Cococcioni, ”Extended DFT+ U+ V method with on-site and inter-site electronic interactions”, J. Phys.: Condens. Matt. 22, 055602(2010).
  * [9] F. Aryasetiawan, K. Karlsson, O. Jepsen, and U. Schonberger, ”Calculations of Hubbard U from first-principles”, Phys. Rev. B 74, 125106 (2006).
  * [10] S. L. Dudarev, P. Liu, D. A. Andersson, et. al., ”Parametrization of LSDA+U for noncollinear magnetic configurations: Multipolar magnetism in UO2”, Phys. Rev. Materials 3, 083802 (2019).
  * [11] I. Timrov, N. Marzari, and M. Cococcioni, ”Hubbard parameters from density-functional perturbation theory”, Phys. Rev. B 98, 085127 (2018).
  * [12] I. Timrov, N. Marzari, and M. Cococcioni, ”Self-consistent Hubbard parameters from density-functional perturbation theory in the ultrasoft and projector-augmented wave formulations”, Phys. Rev. B 103, 045141 (2021).
  * [13] I. Timrov, N. Marzari, and M. Cococcioni, ”HP–A code for the calculation of Hubbard parameters using density-functional perturbation theory”, Comput. Phys. Commun. 279, 108455 (2022).
  * [14] M. Payami, S. Sheykhi, and M. R. Basaadat, ”Self-consistent on-site and inter-site Hubbard parameters within DFT+U+V for UO2 using density-functional perturbation theory”, arXiv:2306.06266v1 [cond-mat.mtrl-sci], https://doi.org/10.48550/arXiv.2306.06266 (2023).
  * [15] J. P. Perdew, et. al., ”Restoring the density-gradient expansion for exchange in solids and surfaces”, Phys. Rev. Lett. 100, 136406 (2008).
  * [16] P. Giannozzi, S. Baroni, N. Bonini, et. al., J. Phys.: Condensed Matt. 21, 395502 (2009).
  * [17] P. Giannozzi, O. Baseggio, P. Bonfà, et. al., J. Chem. Phys. 152, 154105 (2020).
  * [18] A. Dal Corso, ”Pseudopotentials periodic table: From H to Pu”, Comput. Mater. Sci. 95, 337 (2014).
  * [19] N. Marzari, et. al., ”Thermal contraction and disordering of the Al (110) surface”, Phys. Rev. Lett. 82, 3296 (1999).
  * [20] M. Payami, ”Spin-symmetry broken ground-state of UO2 in DFT+U approach: the SMC method”, Iranian J. Phys. Res. 22, 175 (2022), https://doi.org/10.47176/ijpr.22.3.81568; also arXiv:2108.12758v2 [cond-mat.mtrl-sci], https://doi.org/10.48550/arXiv.2108.12758 (2021).
  * [21] M. Payami, ”Comparison of SMC and OMC results in determining the ground-state and meta-stable states solutions for UO2 in DFT+U method”, arXiv:2302.04231v1 [cond-mat.mtrl-sci, https://doi.org/10.48550/arXiv.2302.04231 (2023).
  * [22] J. Schoenes, ”Optical properties and electronic structure of UO2”, J. Appl. Phys. 49, 1463 (1978).


# Electron localization in disordered quantum systems at finite temperatures

Chong Sun  Department of Chemistry, Rice University, Houston, TX 77005-1892
<EMAIL_ADDRESS>
###### Abstract

We study electron localization in disordered quantum systems, focusing on both
individual eigenstates and thermal states. We employ complex polarization as a
numerical indicator to characterize the system’s localization length.
Furthermore, we assess the efficacy of mean-field approximation in providing a
quantitative analysis of such systems. Through this study, we seek to provide
insight into the following aspects: the behavior of electron localization as a
function of interaction, disorder, and temperature, whether thermal states and
highly excited states exhibit similar properties in many-body localized
systems, and the reliability of the mean-field approximation in low-
interaction scenarios.

††preprint: APS/123-QED

Introduction: Many-body localization (MBL) in disordered quantum systems with
interactions [1, 2, 3, 4, 5] has drawn much attention with the potential of
preserving symmetry-broken orders and deferring quantum decoherence [6, 7]. In
such systems, disorder reduces crosstalking between subsystems, and thus
hinders quantum state thermalization [8, 9, 10]. Theoretical characterizations
of MBL include a Poissonian distribution in the ratio of adjacent energy gaps
[11], area-law and logarithm-growth of the entanglement entropy [12, 13, 14],
presence of quasi-local integrals of motion (LIOMs) [15], etc. Notably, MBL
emerges not only in the ground state but across all eigenstates. This suggests
that MBL can be viewed as a form of localization in the Fock space, where
transitions to proximate eigenstates are impeded.

Conflicting conclusions have been made in studies of the dc conductivity of
disordered systems with interactions at finite temperatures [16, 17], where a
nonzero conductivity persists at high temperatures regardless of the disorder
strength. However, these results are not necessarily at odds with established
understandings of MBL for two primary reasons. First, traditional MBL
arguments typically consider an isolated system, where thermalization follows
microcanonical statistics, in contrast to the canonical or grand canonical
frameworks often applied in finite-temperature simulations. Second, because
MBL entails localization within the Fock space, it is expected that individual
eigenstates exhibit localization characteristics, which may not be reflected
in the thermal average across multiple eigenstates. Nevertheless, the
exploration of disordered quantum systems at finite temperatures remains a
valuable pursuit, given the practical challenges of achieving complete
isolation in real-world systems.

In this study, we investigate real-space electron localization within
disordered quantum systems, employing the complex polarization operator, a
concept introduced in a series of works by Resta and Sorella in 1990s [18, 19,
20, 21, 22]. The recently proposed ”imaginary vector potential” adopted a
similar concept [23]. We provide formulations for evaluating the expectation
values of the complex polarization for both individual eigenstates and thermal
states. We demonstrate that full electron localization persists in any given
eigenstate of the system under strong disorder. At high temperatures, however,
electrons exhibit a large spread length even with strong disorder, aligning
with previous findings on finite-temperature conductivity [16, 17]. We also
assess the efficacy of mean-field simulations for such systems under weak
interactions. Utilizing the Hartree-Fock (HF) method as our mean-field
framework, we demonstrate that the HF approach provides reliable quantitative
descriptions of electron localization in conditions of strong disorder.

Electron localization: The electron localization length $\lambda$ in a one-
dimensional system can be described by the electron quadratic spread

$\begin{split}\lambda^{2}=\langle x^{2}\rangle-\langle
x\rangle^{2},\end{split}$ (1)

where $x$ is the position operator of an electron. The localization length
$\lambda$ has a finite value for a localized state and diverges for a
delocalized state. This divergence presents practical challenges for numerical
simulations. In addition, the electron localization length in a system with
periodic boundary conditions (PBC) becomes ambiguous. To address this, Resta
and Sorella introduced the complex polarization operator $Z$, which relates to
the electron quadratic spread by [20, 21]

$\begin{split}\lambda^{2}=-\frac{L^{2}}{2\pi^{2}}\log|Z|,\end{split}$ (2)

where the magnitude of $Z$ satisfies $0\leq|Z|\leq 1$. The electrons are fully
diffused when $|Z|=0$ and fully localized when $|Z|=1$. This relation is still
valid at finite temperatures.

Herein, we introduce the complex polarization operator for a one-dimensional
lattice with PBC as

$\begin{split}Z=\exp\left(i\frac{2\pi
X}{L}\right),\quad{X}=\sum_{j=1}^{N}{x}_{j},\end{split}$ (3)

where $x_{j}$ is the position operator of the $j$th electron, $L$ is the
lattice length, and $N$ is the total number of electrons. The expectation
value of $Z$, denoted as $\langle Z\rangle=|Z|e^{i\gamma}$, is a complex
number, where $\gamma$ is related to the single-point Berry phase [24]. Note
that while $X$ is a one-body operator, $Z$ serves as an $N$-body operator,
capturing the collective behavior of all electrons. Fig. 1 qualitatively
illustrates the relationship between $|Z|$ and the electron localization
length, where the position of an electron on the one-dimensional lattice is
mapped to an angle ranging from $0$ to $2\pi$. When an electron can occupy all
possible positions, $\langle Z\rangle=0$ due to the phase cancellation.
Conversely, complete confinement of the electron results in $\langle
Z\rangle=e^{i\gamma_{0}}$ and $|Z|=1$. Generalization to the multi-dimensional
formulation of Eq. (3) is simply the direct product of $Z$ in each dimension.

Figure 1: Connection between electron localization length and the magnitude of
the complex polarization $|Z|$. The position of an electron on the one-
dimensional lattice is mapped to an angle $\phi$ ranging from $0$ to $2\pi$.

Although $Z$ is a complex $N$-body operator, its application to a Slater
determinant produces another Slater determinant, as a result of Thouless
theorem [25, 26, 27]. This property facilitates the evaluation of $\langle
Z\rangle$ when the state is represented by a Slater determinant
$|\Phi\rangle$, reducing the problem to calculating the overlap between two
Slater determinants

$\begin{split}\langle\Phi|{Z}|\Phi\rangle=\langle\Phi|\Phi^{\prime}\rangle=\det\left[\mathbf{C}_{\mathrm{occ}}^{\dagger}\mathbf{C}^{\prime}_{\mathrm{occ}}\right],\end{split}$
(4)

where $|\Phi^{\prime}\rangle=Z|\Phi\rangle$, and the columns of
$\mathbf{C}_{\mathrm{occ}}$ and $\mathbf{C}^{\prime}_{\mathrm{occ}}$ denote
the occupied orbital coefficients on which $|\Phi\rangle$ and
$|\Phi^{\prime}\rangle$ are constructed, respectively. We have assumed that
$|\Phi\rangle$ is normalized. We write the matrix representations of operators
in the position/site basis.

We now turn to the finite-temperature formulations for non-interacting
systems, derived using the thermofield theory [28]. The thermal average of the
complex polarization $Z$ is given by[29]

$\begin{split}Z(\beta)&=\frac{\text{det}\left[\mathbf{Z}\rho(\beta)+I\right]}{\text{det}\left[\rho(\beta)+I\right]},\\\
\rho(\beta)&=\exp[-\beta(h-\mu)]\end{split}$ (5)

where $\beta=1/k_{B}T$ is the inverse temperature, $\mathbf{Z}$ is the matrix
form of $Z$, $\rho(\beta)$ is the density matrix, $h$ is the effective one-
body Hamiltonian (or Fock) matrix, and $\mu$ is the chemical potential. The
inclusion of $\mu$ ensures the correct thermal average of the electron number.
This derivation is within the grand canonical framework.

Extending Eq. (4) and Eq. (5) to account for a correlated state is
straightforward. Such a state, denoted as $|\Psi\rangle$, can be expressed as
a linear combination of Slater determinants. Consequently,
$\langle\Psi|Z|\Psi\rangle$ becomes a weighted sum of overlaps between various
Slater determinants.

In this study, our focus is on the magnitude of $Z$, i.e., $|Z|$. For
simplicity, we will use $Z$ to refer to $|Z|$ in the following discussions. We
choose the origin of the position at $x_{0}=0$ and limit the lattice to
contain $4n+2$ sites. This choice avoids dealing with a zero band gap
associated with a half-filled lattice of $4n$ sites. Additionally, in the
finite-temperature simulations, we set the Boltzmann constant $k_{B}=1$. The
simulations are performed with a homemade Python package[30] based on
PySCF[31, 32].

Figure 2:  Comparison between Hartree-Fock (HF) and full configuration
interaction (FCI) solutions for the disordered system at various disorder
strengths $W$. (a) Ground-state energy discrepancy $\Delta
E=E_{\mathrm{HF}}-E_{\mathrm{FCI}}$ as a function of $W$. (b) The complex
polarization discrepancy $\Delta Z=Z_{\rm HF}-Z_{\rm FCI}$ as a function of
$W$. (c) $\Delta Z$ as a function of the two-body interaction strength $V$ at
$W=5$. (d) The maximum overlap between a FCI eigenstate and HF eigenstates.
”Ground” refers to the FCI ground state, and ”Middle” indicates the eigenstate
at the center of the FCI energy spectrum.  Figure 3: Complex polarization
evaluated by the Hartree-Fock method. (a) Ground-state (G) and the middle-
state (M) result as a function of disorder strength $W$. Numbers in
parentheses indicate $t^{\prime}$ values. (b) Visualization of occupied
orbitals for the ground state and the middle state, respectively. (c) Finite-
temperature results of $Z$ values as a function of disorder strength $W$. (d)
$Z$ values as a function of temperature $T$.

Model: We investigate a spinless fermionic lattice system characterized by
disordered on-site potential, a setup previously explored in the context of
many-body localization (MBL). We employ periodic boundary conditions (PBC)
throughout the analysis. The Hamiltonian is given by

$\begin{split}H&=\sum_{\langle
i,j\rangle}\left[-t\left(a^{\dagger}_{i}a_{j}+\text{h.c.}\right)+Vn_{i}n_{j}\right]\\\
&-t^{\prime}\sum_{\langle\langle
i,j\rangle\rangle}\left(a^{\dagger}_{i}a_{j}+\text{h.c.}\right)+\sum_{i}w_{i}n_{i},\end{split}$
(6)

where $\langle i,j\rangle$ and $\langle\langle i,j\rangle\rangle$ denote
nearest and next-nearest neighbor sites, respectively, and
$n_{i}=a^{\dagger}_{i}a_{i}$ is the number operator. The hopping amplitudes
$t$ and $t^{\prime}$ promote electron movement among nearest and next-nearest
neighbors, respectively, while the two-body interaction term with strength
$V>0$ discourages adjacent occupations. Random on-site potentials $w_{i}$ are
uniformly distributed between $-W$ and $W$, introducing a disordered chemical
potential to the lattice, with $W$ signifying the disorder strength. This
disorder inhibits electron mobility and induces MBL. By convention, we set
$t=1$. We confine our studies to the half-filled regime, where the number of
electrons equals half the number of sites.

Is Hartree-Fock enough? In the weak-interaction regime, the Hartree-Fock (HF)
method typically yields accurate solutions. In the system studied in this
work, the accuracy of HF solutions is bolstered in scenarios of strong
disorder, where the one-body Hamiltonian becomes predominant. This insight
draws inspiration from the work of Bera et al. [33], who posited that natural
occupation numbers [34] could serve as indicators for MBL, manifesting a step-
like pattern at high disorder strengths. However, we propose a different
interpretation. We suggest that the observed step-like pattern emerges
primarily because, at high $W$, the one-body Hamiltonian prevails, rendering
the system’s eigenstates nearly identical to a single Slater determinant (HF
solution), whose natural occupation numbers are always either $0$ or $1$.
Therefore, we argue that the step-like distribution of natural occupation
numbers is not clearly related to MBL.

In Fig. 2, we compare the HF and full configuration interaction (FCI)
solutions for a $14$-site disordered chain, fixing $t^{\prime}=0$ and $V=1$ in
accordance with the settings used in Ref. [33]. The displayed results are
obtained from averaging over $1,000$ random samples. Fig. 2 (a) and (b)
present the differences in energy ($\Delta E$) and complex polarization
magnitude ($\Delta Z$), respectively, between the HF and FCI ground states.
Both $\Delta E$ and $\Delta Z$ exhibit rapid declines as the disorder strength
$W$ increases. The small values of $\Delta Z$ at low $W$ are attributed to the
near-zero values of $Z$. Fig. 2 (c) explores how the ground-state $Z$ varies
with the two-body interaction strength $V$, under the conditions of
$t^{\prime}=0$ and $W=5$. In the domain of weak interactions, the discrepancy
between the HF and FCI solutions is minimal[35].

Figure 4: Complex polarization evaluated by the full configuration interaction
(FCI). (a) $Z$ as a function of $W$ at various temperatures, with the shade
representing the error bar due to random sampling. The $Z$ values for the
middle state are also presented. (b) Heat map of $Z$ values in relation to
temperature $T$ and disorder strength $W$.

To determine the accuracy of HF solutions for excited states, we analyze the
normalized FCI vector coefficients $\mathbf{c}$. Our FCI calculations are
based on the molecular orbitals (MOs) generated by HF. Hence the elements of
$\mathbf{c}$ represent the overlaps between all possible HF states (Slater
determinants) and a specific FCI eigenstate, i.e.,
$c_{i}=\langle\Phi^{\mathrm{HF}}_{i}|\Psi^{\mathrm{FCI}}\rangle$. A large
value of $\max(|\mathbf{c}|)$ indicates a predominant single Slater
determinant within the FCI solution, and thus high HF accuracy. In Fig. 2 (d),
we focus on the ground state and an eigenstate located in the middle of the
FCI energy spectrum (hereafter referred to as the ”middle state”) to track the
behavior of $\max(|\mathbf{c}|)$ as $W$ increases. The two increasing curves
affirm the precision of HF solutions for both ground and excited states under
sufficiently high disorder strength $W$. This result also suggests further
investigation into the potential simplification of many-body localization
(MBL) to Anderson localization under strong disorder conditions, a question
that extends beyond the scope of this study.

HF and FCI results: We analyze the electron localization of a $30$-site chain
evaluated with the HF method, where the results are evaluated by averaging
$N_{\mathrm{rep}}=5,000$ random samples, shown in Fig. 3. The error bars are
indicated by the shading around the curves, calculated from the worst-case
sampling error $1/\sqrt{N_{\mathrm{rep}}}$. Fixing $V=1$, we explore scenarios
with ($t^{\prime}=1$) and without ($t^{\prime}=0$) next-nearest-neighbor
hopping. Fig. 3 (a) verifies that full localization can be attained in both
the ground and excited states, where we picked the middle state to showcase
the excited states[36]. The occupied MOs for the two states are shown as the
orange solid lines in Fig. 3 (b), where the MOs are sorted according to their
energy levels $\varepsilon$. The inclusion of next-nearest-neighbor hopping
discourages electron localization.

Proceeding to the finite-temperature results depicted in Fig. 3 (c), we
observe a slower growth of $Z$ with respect to $W$, with full electron
localization not observed even at very high $W$ values. This indicates that
thermal fluctuations diminish the impact of disorder, making the system’s
response to changes in $W$ less pronounced. A notable crossover at $T=0.2$
between the cases with ($t^{\prime}=1$) and without ($t^{\prime}=0$) next-
nearest-neighbor hopping suggests the influence of frustration. In Fig. 3 (d),
$Z$ as a function of $T$ again displays a crossover around $T=0.25$ at $W=10$,
underscoring the nuanced dynamics introduced by temperature and hopping
interactions.

Figure 5: Extrapolation of the $Z$ values at $W=5$, $V=1$ and $t^{\prime}=0$
in relation to the number of sites $L$. (a) FCI solutions. (b) HF solutions.
Dashed lines indicate the function fitting of the discrete data.

Fig. 4 displays the FCI results for the electron localization. Unlike the
finite-temperature HF simulations, where the grand-canonical statistics are
used, finite-temperature FCI adopted the canonical ensemble picture for a
lower computational cost. Compared to the HF simulations, the system under the
canonical statistics experiences less thermal effect. Fig. 4 (a) profiles $Z$
as a function of $W$ at various temperatures and for different eigenstates
(ground state and middle state). The simulations are performed on a $14$-site
chain with $N_{\mathrm{rep}}=1,000$ random samples, fixing $V=1$ and
$t^{\prime}=0$. At $T=0.2$, the curve behaves like the ground state, although
the maximum $Z$ value reached is slightly lower than the full localization
value, $1$. As the temperature increases, the curves display a linear manner,
similar to the HF curves. Contrary to the high-temperature curves, the $Z$
value of the middle state approaches near full localization at sufficiently
high $W$ values. Fig. 4 (b) presents a heatmap plot of $Z$ as a function of
both disorder $W$ and temperature $T$. The simulations are performed on a
small $10$-site chain with $V=1$ and $t^{\prime}=0$. The clear boundary in the
heat map further confirmed the absence of full electron localization at high
temperatures.

Thermodynamic limit: To mitigate finite-size effects, we extrapolate the above
simulations to the thermodynamic limit (TDL), setting $t^{\prime}=0$, $V=1$,
$W=5$ and $N_{\mathrm{rep}}=1,000$ across all simulations. We plot $Z$ against
$1/L$ in Fig. 5. The FCI simulations follow canonical statistics, whereas HF
simulations follow grand canonical statistics. The extrapolation reveals
distinct behaviors for ground-state and finite-temperature conditions. At the
ground state or very low temperatures, $Z$ linearly increases towards a higher
value as $L$ approaches infinity, fitting the data with $Z(L)=a+b/L$.
Conversely, at higher temperatures, $Z$ trends towards a lower value in a
polynomial manner, fitting the data with $Z(L)=a+b/L+c/L^{2}$. Despite
expectations of canonical and grand canonical solutions becoming equivalent at
TDL, slight differences are observed, especially at lower temperatures where
the system exhibits stronger correlations and HF simulations are less
accurate. Additionally, the limitation to only three FCI data points may
impact the precision of extrapolation. Integrating the heatmap from Fig. 4 (b)
with the extrapolation results, we anticipate a more defined boundary between
localized and delocalized states in the heatmap at TDL.

Conclusions: This study provided a comprehensive analysis of electron
localization within disordered quantum systems, exploring both individual
eigenstates and thermal states, utilizing complex polarization as the
theoretical indicator. We found that full electron localization can be
achieved in individual eigenstates, while only partial localization is
possible at high temperatures, despite strong disorder. Our assessment of the
Hartree-Fock method confirmed its accuracy for the system under study in
conditions of high disorder. However, for low-temperature states or highly
excited states within the weak to moderate disorder regime, correlated methods
or mean-field methods targeting the excited states, e.g., $\Delta$-SCF[37,
38], are necessary. The methodologies and insights gained from this work offer
valuable perspectives for characterizing other disordered quantum systems and
could inform the discovery of disordered materials for applications like
robust quantum memory.

Acknowledgments: We thank Garnet Chan, Sandro Sorella, Raffaele Resta, Soumya
Bera, Huanchen Zhai, Xing Zhang, and Zijian Zhang for their insightful
discussions.

## References

  * Gornyi  _et al._ [2005] I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Interacting electrons in disordered wires: Anderson localization and low-$t$ transport, Phys. Rev. Lett. 95, 206603 (2005).
  * Basko  _et al._ [2006a] D. Basko, I. Aleiner, and B. Altshuler, Metal–insulator transition in a weakly interacting many-electron system with localized single-particle states, Annals of Physics 321, 1126 (2006a).
  * Basko  _et al._ [2006b] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, On the problem of many-body localization (2006b), arXiv:cond-mat/0602510 [cond-mat.mes-hall] .
  * Nandkishore and Huse [2015] R. Nandkishore and D. A. Huse, Many-body localization and thermalization in quantum statistical mechanics, Annual Review of Condensed Matter Physics 6, 15 (2015).
  * Abanin and Papić [2017] D. A. Abanin and Z. Papić, Recent progress in many-body localization, Annalen der Physik 529, 1700169 (2017).
  * Huse  _et al._ [2013] D. A. Huse, R. Nandkishore, V. Oganesyan, A. Pal, and S. L. Sondhi, Localization-protected quantum order, Phys. Rev. B 88, 014206 (2013).
  * Parameswaran and Vasseur [2018] S. A. Parameswaran and R. Vasseur, Many-body localization, symmetry and topology, Reports on Progress in Physics 81, 082501 (2018).
  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Thermalization and its mechanism for generic isolated quantum systems, Nature 452, 854 (2008).
  * Kaufman  _et al._ [2016] A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, and M. Greiner, Quantum thermalization through entanglement in an isolated many-body system, Science 353, 794 (2016).
  * Deutsch [2018] J. M. Deutsch, Eigenstate thermalization hypothesis, Reports on Progress in Physics 81, 082001 (2018).
  * Oganesyan and Huse [2007] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev. B 75, 155111 (2007).
  * Bardarson  _et al._ [2012] J. H. Bardarson, F. Pollmann, and J. E. Moore, Unbounded growth of entanglement in models of many-body localization, Phys. Rev. Lett. 109, 017202 (2012).
  * Bauer and Nayak [2013] B. Bauer and C. Nayak, Area laws in a many-body localized state and its implications for topological order, Journal of Statistical Mechanics: Theory and Experiment 2013, P09005 (2013).
  * Serbyn  _et al._ [2013a] M. Serbyn, Z. Papić, and D. A. Abanin, Universal slow growth of entanglement in interacting strongly disordered systems, Phys. Rev. Lett. 110, 260601 (2013a).
  * Serbyn  _et al._ [2013b] M. Serbyn, Z. Papić, and D. A. Abanin, Local conservation laws and the structure of the many-body localized states, Phys. Rev. Lett. 111, 127201 (2013b).
  * Karahalios  _et al._ [2009] A. Karahalios, A. Metavitsiadis, X. Zotos, A. Gorczyca, and P. Prelovšek, Finite-temperature transport in disordered heisenberg chains, Phys. Rev. B 79, 024425 (2009).
  * Berkelbach and Reichman [2010] T. C. Berkelbach and D. R. Reichman, Conductivity of disordered quantum lattice models at infinite temperature: Many-body localization, Phys. Rev. B 81, 224429 (2010).
  * Resta [1992] R. Resta, Theory of the electric polarization in crystals, Ferroelectrics 136, 51 (1992).
  * Resta [1993] R. Resta, Macroscopic electric polarization as a geometric quantum phase, Europhysics Letters (EPL) 22, 133 (1993).
  * Resta [1998] R. Resta, Quantum-mechanical position operator in extended systems, Phys. Rev. Lett. 80, 1800 (1998).
  * Resta and Sorella [1999] R. Resta and S. Sorella, Electron localization in the insulating state, Phys. Rev. Lett. 82, 370 (1999).
  * Resta [2002] R. Resta, Why are insulators insulating and metals conducting?, Journal of Physics: Condensed Matter 14, R625 (2002).
  * Heußen  _et al._ [2021] S. Heußen, C. D. White, and G. Refael, Extracting many-body localization lengths with an imaginary vector potential, Phys. Rev. B 103, 064201 (2021).
  * Resta [2000] R. Resta, Manifestations of berry’s phase in molecules and condensed matter, Journal of Physics: Condensed Matter 12, R107 (2000).
  * Thouless [1960] D. Thouless, Stability conditions and nuclear rotations in the hartree-fock theory, Nuclear Physics 21, 225 (1960).
  * Rosensteel and Rowe [1981] G. Rosensteel and D. J. Rowe, Nondeterminantal hartree-fock theory, Phys. Rev. A 24, 673 (1981).
  * Ripka  _et al._ [1986] S. Ripka, J. Blaizot, and G. Ripka,  _Quantum Theory of Finite Systems_ (MIT Press, 1986).
  * Umezawa  _et al._ [1982] H. Umezawa, H. Matsumoto, and M. Tachiki,  _Thermo field dynamics and condensed states_ (North-Holland Publishing Company, 1982).
  * Sun [2021] C. Sun,  _Finite Temperature Simulations of Strongly Correlated Systems_ , Ph.D. thesis (2021).
  * Sun [2024] C. Sun, ComPol package github link (2024).
  * Sun  _et al._ [2018] Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K.-L. Chan, Pyscf: the python-based simulations of chemistry framework, WIREs Computational Molecular Science 8, e1340 (2018).
  * Sun  _et al._ [2020] Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov, and G. K.-L. Chan, Recent developments in the PySCF program package, The Journal of Chemical Physics 153, 024109 (2020).
  * Bera  _et al._ [2015] S. Bera, H. Schomerus, F. Heidrich-Meisner, and J. H. Bardarson, Many-body localization characterized from a one-particle perspective, Phys. Rev. Lett. 115, 046603 (2015).
  * [34] Natural occupation numbers are the eigenvalues of the one-particle reduced density matrix (1RDM). The corresponding eigenvectors are called natural orbitals .
  * [35] The observed reduction in $\Delta Z$ at $V=9$ is, once again, a consequence of the inherently small $Z$ values in this regime.
  * [36] The middle states of HF and FCI are not equivalent. The FCI middle state is chosen from the middle of the eigenenergy spectrum, while the HF one is constructed by choosing the middle set of molecular orbitals (MOs) to construct the Slater determinant. Nevertheless, the two states should be close to each other .
  * Jones and Gunnarsson [1989] R. O. Jones and O. Gunnarsson, The density functional formalism, its applications and prospects, Rev. Mod. Phys. 61, 689 (1989).
  * Gavnholt  _et al._ [2008] J. Gavnholt, T. Olsen, M. Engelund, and J. Schiøtz, $\Delta$ self-consistent field method to obtain potential energy surfaces of excited molecules on surfaces, Phys. Rev. B 78, 075441 (2008).