Doctor of Philosophy California Institute of Technology figures/caltech.png
2021 December 1, 2020

0000-0002-8299-9094

All rights reserved

# Finite Temperature Simulations of Strongly Correlated Systems

Chong Sun  Pasadena, California

###### Abstract

This thesis describes several topics related to finite temperature studies of
strongly correlated systems: finite temperature density matrix embedding
theory (FT-DMET), finite temperature metal-insulator transition, and quantum
algorithms including quantum imaginary time evolution (QITE), quantum Lanczos
(QLanczos), and quantum minimally entangled typical thermal states (QMETTS)
algorithms.

While the absolute zero temperature is not reachable, studies of physical and
chemical problems at finite temperatures, especially at low temperature, is
essential for understanding the quantum behaviors of materials in realistic
conditions. Here we define low temperature as the temperature regime where the
quantum effect is not largely dissipated due to thermal fluctuation. Treatment
of systems at low temperature is specially difficult compared to both high
temperature - where classical approximation can be applied - and zero
temperature where only the ground state is required to describe the system of
interest. FT-DMET is a wavefunction-based embedding scheme which can handle
finite temperature simulations of a variety of strongly correlated problems.
The "high-level in low-level" framework enables FT-DMET to tackle large bulk
sizes and capture the majority of the entanglement at the same time. FT-DMET
formulations and implementation details for both model systems and ab initio
problems are provided in Chapter 2 and Chapter 3.

Metal-insulator transition is a common but important phase transition in many
strongly correlated materials. The widely accepted scheme to distinguish an
insulator from a metal is band structure theory based on a single-particle
picture. However, insulating phases caused by disorder or strong correlation
cannot be explained merely with the band structure. In Chapter 4, we
demonstrate that electron locality/mobility is a more general criteria to
detect metal-insulator transition. We further introduce complex polarization
as the order parameter to reflect the electron locality/mobility and provide a
formalism based on thermofield theory to evaluate the complex polarization at
finite temperature.

Quantum algorithms are designed to perform simulations on a quantum device.
The infrastructure of a quantum processing unit (QPU) utilizes the
superposition property of quantum bits (qubits), and thus can potentially
outplay the classical simulations in computational scaling for certain
problems. In Chapter 5, we introduce the QITE algorithm, which can be applied
to quantum simulations of both ground state and finite temperature problems.
We further introduce a subspace method, QLanczos algorithm, and a a finite
temperature quantum algorithm, QMETTS, where QITE is used as a building block
for the two algorithms. We demonstrate above quantum algorithms with
simulations on both classical computers and quantum computers.

[logo]

###### Acknowledgements.

I am deeply blessed as a member of the Caltech community. Caltech provided me
opportunities to participate in advanced research projects via collaborations
with excellent researchers. I am truly thankful to the institute and every
member of this big family. My advisor, Garnet Kin-Lic Chan, has provided me
invaluable guidance and support throughout my graduate studies. From him, I
learned to always check my hypothesis carefully against data. Whenever the
data not seem reasonable, I should question my code first before questioning
the theory or algorithm. He also provided plenty of opportunities for me to
attend conferences and communicate with researchers in the field. Garnet is
and will be the role model as a scientist to me for the rest of my life. I
would like to thank my committee, Professor Mitchio Okumura, Professor Thomas
Miller and Professor Austin Minnich. They provided many helpful advices during
my graduate career and tried to bring the best out of me. I am also thankful
to Professor Lu Wei who was really patient and helpful with my many questions
about stimulated Raman spectroscopy. Thank you, CCE administrative staff, in
particular Alison and Elizabeth. Without the help from you, I would not have
been able to fully focus on research without worrying about many tough
errands. I am grateful to be a member of the Chan group and work with so many
awesome colleagues. Everyone in the Chan group is nice and always willing to
help. I have been working closely with Zhihao, Mario, Ushnish, and Reza, from
whom learned useful knowledge and skills. I am still close friends with
previous group members such as Boxiao, Zhendong, and Mario, who constantly
provide valuable suggestions to me when I need help. I also enjoyed group
activities. Before the pandemic, we hung out monthly and tried many good or
mediocre restaurants. We also had trips to Yosemite, Sequoia, and Universal
Studios. My graduate school life is full of fun because of the Chan group
villagers. Lastly, I would like to thank my family. My parents are the best
parents I could ever dream of. They did not have opportunities for good
education, but they value education for my sister and I, and fully support my
career as a scientist. I was lucky to have my young sister as my close friend
since childhood when most of my friends are only children in their families.
Having a smart and aggressive sibling was overall helpful to push me to work
harder. My husband, James has always been there to make me laugh when I was
unhappy with my research progress. Theoretically, my cat Jujube should thank
me for providing her a home and food, but I know she does not think in that
way. I am thankful for her company and for not making loud noises when I have
Zoom meetings.

[0.7] _To my beloved parents_

[iknowwhattodo]

## Contents

refsectionchapter*.1refsectionchapter*.2refsectionchapter*.3refsectionsection*.4refsectionsection*.5refsectionsection*.6refsectionchapter.1refsectionsection.1.1refsectionequation.1.1.3refsectionItem.13refsectionsection.1.2refsectionchapter.2refsectionsection.2.1refsectionsection.2.2refsectionsection.2.3refsectionsection.2.3refsectionequation.2.3.1refsectionequation.2.3.7refsectionequation.2.3.10refsectionequation.2.3.14refsectionequation.2.3.15refsectionequation.2.3.15refsectionequation.2.3.21refsectionsection.2.4refsectionsection.2.4refsectionsection.2.4refsectionfigure.caption.15refsectionsection.2.5refsectionchapter.3refsectionsection.3.1refsectionsection.3.2refsectionsection.3.3refsectionsection.3.3refsectionequation.3.3.1refsectionfigure.caption.20refsectionfigure.caption.21refsectionfigure.caption.21refsectionsection.3.4refsectionsection.3.5refsectionchapter.4refsectionsection.4.1refsectionsection.4.2refsectionsection.4.3refsectionequation.4.3.5refsectionequation.4.3.17refsectionsection.4.4refsectionsection.4.5refsectionsection.4.6refsectionsection.4.7refsectionchapter.5refsectionsection.5.1refsectionsection.5.2refsectionsection.5.3refsectionsection.5.4refsectionsection.5.5refsectionsection.5.6refsectionequation.5.6.22refsectionsection.5.7refsectionchapter.Arefsectionsection.A.1refsectionsection.A.2refsectionsection.A.3refsectionchapter.Brefsectionsection.B.1refsectionsection.B.2refsectionsection.B.3refsectionsection.B.4refsectionsection*.47

###### List of Figures

  1. 1 Structure of the matrix product states used in the purification approach of the finite temperature density matrix renormalization group algorithm.
  2. 2 Evaluating the total energy with density matrix embedding theory. (a) A hydrogen ring composed of $10$ atoms obeying the periodic boundary condition, and the impurity (supercell) is two adjacent atoms. The total energy equals the energy of the supercell times the number of supercells. (b) A single ligand heme molecule is divided into $5$ non-overlapping fragments, and the DMET energy of each fragment is calculated. The total energy is the sum of energies from all fragments.
  3. (a) $E_{\text{tot}}=5E_{0}$
  4. (b) $E_{\text{tot}}=\sum_{i=0}^{4}E_{i}$
  5. 3 Weight of the entanglement with the impurity on environmental sites. The weights are evaluated as the square norm of projections of environmental sites to the bath space. Due to periodicity of the system, only half of the environmental sites are shown in the figure.
  6. 4 Error in energy per site (units of $t$) of FT-DMET for the 1D Hubbard model at $U=2$ and $U=4$ (2 impurity sites and half-filling) with bath orbitals generated via the density matrix $\gamma$ (Eq. (45)) (blue lines) or lattice Hamiltonian $h$ (Eq. (43)) (orange lines) as a function of inverse temperature $\beta$. The numbers in parentheses denote the number of bath orbitals. The grey area denotes the ground state error with 2 impurity orbitals.
  7. 5 Percentage error of the FT-DMET (with 2 impurity sites and 2 bath orbitals) energy per site vs. ED (4 sites) on a non-embedded cluster with PBC and APBC boundary conditions for the 1D Hubbard model at various $U$ and $\beta$.
  8. 6 Absolute error of the FT-DMET energy per site of the 1D Hubbard model at half-filling as a function of impurity and bath size. $InBm$ denotes $n$ impurity sites and $m$ bath orbitals. Increasing impurity (blue lines); increasing bath (orange lines). The grey band depicts the ground state error with $2$ impurity sites and $2$ bath orbitals.
  9. 7 Absolute error of the FT-DMET entropy per site of the 1D Hubbard model at half-filling as a function of the number of bath sites. The right panels show the absolute entropy.
  10. 8 Energy per site (units of $t$) of a 16-site Hubbard chain with periodic boundary conditions at $U=4$ as a function of the chemical potential $\mu$ at various $\beta$ values. The difference between the DMRG and DMET ($I2B6$) energies per site is $0.01-1.4\%$. Solid lines: DMRG energies; dashed lines: DMET energies; pentagons: Bethe ansatz.
  11. 9 Energy per site versus U (units of $t$) of the 2D Hubbard model at half-filling with FT-DMET ($2\times 2$ cluster with 4 and 8 bath orbitals), DCA (34, 72 and $2\times 2$ site clusters).
  12. 10 Néel transition for the 2D Hubbard model within quantum impurity simulations. (a) Antiferromagnetic moment $m$ as a function of $T$ with various $U$ values (units of $t$); (b) Néel temperature $T_{N}$ calculated with FT-DMET, single-site DMFT and DCA. DMFT data is taken from Ref. [1], DCA/NCA data for $U=4$ is taken from Ref. [2], DCA/QMC data for $U=6$ is taken from Ref. [3], and DCA/QMC data for $U=8$ is taken from Ref. [4].
  13. (a)
  14. (b)
  15. 11 2D Hubbard antiferromagnetic moment (color scale) as a function of $T$ and $U$ (units of $t$) in FT-DMET ($2\times 2$ impurity, 4 bath orbitals.)
  16. 12 Bath orbitals from singular value decomposition (SVD) of the off-diagonal block of the mean-field density matrix. The whole square represents the mean-field density matrix of size $N\times N$, with $N$ being the total number of orbitals, and the first $n_{\text{imp}}$ columns/rows of the matrix are orbitals in the impurity. (a) Standard DMET routine computes the bath orbitals by the SVD of the off-diagonal block (orange blocks on the left and top, the block on the left corresponds to Eq. (49)). (b) Ab initio DMET computes the bath orbitals by the SVD of only the valence columns in the off-diagonal blocks (green block).
  17. (a)
  18. (b)
  19. 13 Fermi-Dirac distribution of electrons on Hartree-Fock molecular orbitals for H30 chain with STO-6G basis. $T$ is in unit Hartree.
  20. 14 Accuracy test on FT-DMRG and LT-DMRG solvers against exact diagonalization. The label "FT($x$)" stands for FT-DMRG solver with $\tau=x$, and the label "LT($x,y$)" stands for LT-DMRG solver with $x$ Davidson roots and $y$ electron deviations from half-filling for both spins.
  21. (a) $R=1.5$ Bohr
  22. (b) $R=3.0$ Bohr
  23. 15 Magnetic moment of a $22$-atom chain at ground state (left panel) and $T=0.02$ Hartree (right panel) with STO-6G, 6-31G, and CC-PVDZ basis sets.
  24. (a)
  25. (b)
  26. 16 Double occupancy of a $22$-atom chain at ground state (left panel) and $T=0.02$ Hartree (right panel) with STO-6G, 6-31G, and CC-PVDZ basis sets. The insets show a sudden change of the gradient of $D$ as a function of $R$, indicating metal to insulator transition.
  27. (a)
  28. (b)
  29. 17 Dissociation curve of hydrogen chain at $T=0.05$ Hartree compared to AFQMC. The AFQMC data is extracted from Ref. [5].
  30. 18 Staggered magnetic moment of hydrogen chain with periodic boundary condition at ground state, $T=0.02$, $T=0.05$, and $T=0.1$. The unit of $T$ is Hartree.
  31. 19 Dispersion relation and energy levels of the half-filled tight binding model for (a) $L=8$ and (b) $L=10$. Solid black dots are occupied orbitals, blank circles are unoccupied orbitals, and circles with stripes are partially occupied orbitals due to degeneracy.
  32. (a) $L=8$
  33. (b) $L=10$
  34. 20 Complex polarization of the tight binding model ($L=42$) with staggered potential $u$ at ground state (GS), $T=0.2$, $T=0.5$, and $T=1.0$, respectively.
  35. 21 Phase diagram of the tight binding model ($L=42$) with the staggered potential $u$. The blue area corresponds to $Z_{N}>0$ (insulator) and the white area corresponds to $Z_{N}=0$ (metal). The 2D plot is smoothed by Bessel interpolation. Grid: $20$ points in the $x$-axis and $10$ points in the $y$-axis.
  36. 22 Energy gap of the half-filled staggered tight binding model against the staggered potential $u$. Inset: energy gap for $u\in[0t,0.1t]$.
  37. 23 Complex polarization, magnetic moment, population on 2s orbital and energy gap of hydrogen chain with unrestricted Hartree-Fock method. Note that the complex polarization at $T=0.01$ is not presented here due to overflow.
  38. 24 Complex polarization, magnetic moment, population on 2s orbital, and energy gap of hydrogen chain from DFT with PBE functional.
  39. 25 Complex polarization, magnetic moment, population on 2s orbital, and energy gap of hydrogen chain from DFT with B3LYP functional.
  40. 26 Quantum imaginary time evolution algorithm and correlation length. (a) Schematic of the QITE algorithm. Top: imaginary-time evolution under a geometric $k$-local operator $\hat{h}[m]$ can be reproduced by a unitary operation acting on a group of $D>k$ qubits. Bottom: exact imaginary-time evolution starting from a product state requires unitaries acting on a domain $D$ that grows with $\beta$. (b,c) Left: mutual information $I(i,j)$ between qubits $i$, $j$ as a function of distance $d(i,j)$ and imaginary time $\beta$, for a 1D (b) and a 2D (c) FM transverse-field Ising model, with $h=1.25$ (1D) and $h=3.5$ (2D). The mutual information is seen to saturate at longer times. Right: relative error in the energy $\Delta E$ and fidelity $F=|\langle\Phi(\beta)|\Psi\rangle|^{2}$ between the finite-time state $\Phi(\beta)$ and infinite-time state $\Psi$ as a function of imaginary time. The noise in the 2D fidelity error at large $\beta$ arises from the approximate nature of the algorithm used.
  41. 27 Energy calculations with QITE and QLanczos algorithms. Top: QITE energy $E(\beta)$ (a) and fidelity (b) between finite-time state $\Phi(\beta)$ and exact ground state $\Psi$ as function of imaginary time $\beta$, for a 1D 10-site Heisenberg model, showing the convergence with increasing unitary domains of $D=2-8$ qubits. Bottom: QITE (dashed red, dot-dashed green lines) and QLanczos (solid red, solid green lines) energies as function of imaginary time $\beta$, for a 1D Heisenberg model with $N=20$ qubits, using domains of $D=2$ (c) and $4$ qubits (d), showing improved convergence of QLanczos over QITE. Black line is the exact ground-state energy/fidelity.
  42. 28 QITE energy evaluations. (a) QITE energy $E(\beta)$ as a function of imaginary time $\beta$ for a 6-site 1D long-range Heisenberg model, for unitary domains $D=2-6$; (b) a 4-site 1D Hubbard model with $U/t=1$, for unitary domains $D=2,4$; (d) the H2 molecule in the STO-6G basis. (c) Probability of MAXCUT detection, $P(C=C_{max})$ as a function of imaginary time $\beta$, for the $6$-site graph in the panel. Black line is the exact ground-state energy/probability of detection.
  43. 29 QITE, QLanczos, and QMETTS energies $E(\beta)$ as a function of imaginary time $\beta$ for 1-qubit field model using the QVM and QPU (qubit 14 on Aspen-1) and 2-qubit AFM transverse field Ising model using the QVM and QPU (qubit 14, 15 on Aspen-1). (a) Ground state energies for 1-qubit field model using the QVM and QPU (qubit 14 on Aspen-1); (b) ground state energies for 2-qubit AFM transverse field Ising model using the QVM and QPU (qubit 14, 15 on Aspen-1); (c) finite temperature energies for 1-qubit field model using the QVM and QPU (qubit 14 on Aspen-1) ; and (d) finite temperature energies for 2-qubit AFM transverse field Ising model using the QVM. Black lines are the exact solutions.
  44. 30 Thermal (Gibbs) average $\langle E\rangle$ at temperature $\beta$ from QMETTS for a 1D 6-site Heisenberg model (exact emulation).

###### List of Tables

  1. 1 PM-AFM transition bond length $R$ (in Bohr) at ground state and low temperature.
  2. 2 Metal-insulator transition bond length $R$ (in Bohr) at ground state and low temperature.
  3. 3 QPUs: 1-qubit QITE and QLanczos.
  4. 4 QPUs: 2-qubit QITE and QLanczos.
  5. 5 QPUs: 1-qubit METTS.
  6. 6 QVM: 2-qubit QITE and QLanczos.
  7. 7 QVM: 1-qubit METTS.
  8. 8 QVM: 2-qubit METTS.

##  Chapter 1 Introduction

We live in an era where the computational power is one of the main driving
forces for science and technology development. The hardware breakthroughs in
supercomputers, graphical processing unit (GPU) and quantum computers made
heavy computational tasks possible. The development in machine learning
algorithms and artificial intelligence changed the way people live
tremendously. Many new materials and drugs are discovered via computational
simulations, saving hundreds of laboratory hours. We believe in the
computational power to bring us new knowledge and concepts, as well as to
solve fundamental problems that remain unclear for decades. In quantum
chemistry and condensed matter physics, those hard problems include the phase
diagram of high-temperature superconductors (HTSC) [6, 7], the mechanism of
nitrogen fixation [8, 9], protein folding [10], etc. The barrier for efficient
simulations of the above problems is usually either the system size is too big
or the interaction is too complicated. The strongly correlated systems,
unfortunately, have both of the above two barriers. The hallmark of strongly
correlated systems is localized orbitals such as $d$ and $f$ orbitals, where
electrons experience strong Coulomb repulsion. For instance, transition metal
compounds usually contain strong correlations due to the localized $3d$
orbitals. Strongly correlated materials attract tremendous interest of both
experimental and theoretical researchers because they exhibit a plethora of
exotic phases or behaviors: HTSC, spintronic materials [11], Mott insulators
[12], etc. Those strongly correlated behaviors evoked novel applications such
as quantum processing units [13], superconducting magnets [14, 15], and
magnetic storage [16]. Being able to simulate strongly correlated problems and
thus understand the physics behind them has been a key task for theoretical
and computational chemists.

This thesis focuses on developing theoretical and computational approaches to
simulate strongly correlated problems at finite temperature. While ground
state simulations provide basic information on the system such as ground state
energy and band gap, finite temperature is where the real-life phase
transitions happen. The complexity of a quantum many-body problem can be
described by a term called entanglement. At ground state away from the
critical point, the entanglement is bounded by the area law [17]. However, at
finite temperature, especially low temperature where the quantum effect is not
fully dissipated by thermal fluctuation, the area law is no longer valid. One
would expect the entanglement strength to decay while the entanglement length
to grow with temperature. The interplay between the entanglement strength and
entanglement length decides the complexity of the system. Normally one would
expect more computational efforts for finite temperature calculations than
ground state calculations.

The complexity of finite temperature calculations can also be understood in
the ensemble picture. Most of the physical and chemical systems can be seen as
open systems, where the thermodynamic statistics is described by the grand
canonical ensemble. In the grand canonical ensemble picture, both energy
fluctuations and particle number fluctuations are involved. The system at
temperature $T$ is fully described by the density matrix

$\hat{\rho}(T)=e^{-(\hat{H}-\mu\hat{N})/k_{B}T},$ (1)

where $\hat{H}$ is the Hamiltonian, $\mu$ is the chemical potential, $\hat{N}$
is the number operator and $k_{B}\approx 1.38\times
10^{-23}\mathrm{J}\cdot\mathrm{K}^{-1}$ is the Boltzmann constant. The
partition function is defined as the trace of the density matrix:
$\mathcal{Z}=\text{Tr}(\hat{\rho})$. If one choose the eigenstates of the
Hamiltonian $\hat{H}$ as the basis to perform the trace summation, each
eigenstate would participate in the statistics with probability

$P(n,i)=e^{-(\varepsilon_{i}^{n}-\mu n)/k_{B}T}/\mathcal{Z},$ (2)

where $\varepsilon_{i}^{n}$ is the eigenvalue of the $i$th eigenstate in the
Fock space of $n$ particles. If $\varepsilon_{i}^{n}<\mu n$, $P(n,i)$
decreases to $1/\mathcal{N}$ as temperature rises; if $\varepsilon_{i}^{n}>\mu
n$, $P(n,i)$ increases to $1/\mathcal{N}$ as temperature rises, where
$\mathcal{N}$ is the total number of eigenstates. At $T=0$, only the ground
state is involved; as one raises the temperature, the contribution from the
ground state drops and excited states enter the ensemble. Eventually at
infinite temperature, all states are equally involved with a probability
$1/\mathcal{N}$. The inclusion of many excited states is the source of the
high complexity of finite temperature simulations. For instance, for an
electronic structure problem with $L$ orbitals, where each orbital can take
four states: $|0\rangle$, $|\uparrow\rangle$, $|\downarrow\rangle$, and
$|\uparrow\downarrow\rangle$. The total number of states is
$\mathcal{N}=4^{L}$, which scales exponentially with $L$.

Albeit the high computational cost of finite temperature simulations, there
exist a variety of finite temperature algorithms that can fulfill different
computational tasks. Section 1 presents a detailed review of current finite
temperature algorithms. We hope this review could be helpful to researchers
who are interested in learning about or using finite temperature algorithms.
Section 2 provides an outline for the rest of the chapters in this thesis.

###  1 Finite temperature algorithms

At finite temperature $T$, the grand canonical ensemble average of an operator
$\hat{O}$ is evaluated by

$\langle\hat{O}\rangle(T)=\frac{\mathrm{Tr}\left(e^{-(\hat{H}-\mu\hat{N})/k_{B}T}\hat{O}\right)}{\mathrm{Tr}\left(e^{-(\hat{H}-\mu\hat{N})/k_{B}T}\right)}.$
(3)

There are generally two approaches to design a finite temperature algorithm:
(i) directly evaluate the trace in Eq. (3) by summation over the expectation
values under an orthonormal basis; (ii) imaginary time evolution from infinite
temperature. Theoretically the two approaches are all based on Eq. (3), so one
could argue that there is no big difference between the two approaches.
Technically, however, the first approach usually involves exact or approximate
diagonalization of the Hamiltonian, while the latter approach does not. In the
following, we will discuss the two approaches with some example algorithms.

####  1.1 Direct evaluation of the trace

We first discuss the non-interacting case. For a non-interacting Hamiltonian,
only one-body terms are involved, and the Hamiltonian can be simply written as
an $L\times L$ matrix, where $L$ is the number of orbitals in the system. For
most cases, this $L\times L$ Hamiltonian matrix can be directly diagonalized,
with eigenvalues $\varepsilon_{i}$ and eigenvectors $|\phi_{i}\rangle$
(molecular orbitals, MOs). A direct implementation of Eq. (3) is to construct
Slater determinants of all possible particle numbers and evaluate the traces,
where the number of Slater determinants in the summation scales exponentially
with $L$. Luckily, for non-interacting electrons, the grand canonical density
matrix can be evaluated by Fermi-Dirac equation

$\rho=\frac{1}{1+e^{(H-\mu\mathbb{I})/k_{B}T}},$ (4)

where $\mathbb{I}$ is the identity matrix. The occupation numbers on MOs are
the diagonal terms of the density matrix:
$n_{i}=1/(1+e^{(\varepsilon_{i}-\mu)/k_{B}T})$. Thus Eq. (3) can be rewritten
as

$\langle\hat{O}\rangle_{NI}(T)=\sum_{ij}\rho_{ij}\langle\phi_{j}|\hat{O}|\phi_{i}\rangle,$
(5)

where the subscript "NI" stands for "non-interacting".

Finite temperature Hartree-Fock is an example of the above approach, with the
algorithm summarized in Algorithm 1.1.

1

Construct the Fock matrix $F$ from the Hamiltonian. Define $F^{\prime}$ as
identity.

2 while  _$F\neq F^{\prime}$_ do

3

  1. 1.

Store the Fock matrix into $F^{\prime}=F$;

  2. 2.

Diagonalized $F$ to get MO energies and coefficients;

  3. 3.

Calculate the chemical potential $\mu$ by minimizing
$(N_{\text{elec}}-\sum_{i}n_{i})^{2}$, where $N_{\text{elec}}$ is the target
electron number and $n_{i}$ is the occupation number of the $i$th MO;

  4. 4.

Calculate density matrix $\rho$ from Eq. (4) by substituting $H$ with $F$;

  5. 5.

Evaluate the new Fock matrix $F$ from the density matrix $\rho$ as in ground
state Hartree-Fock algorithm;

4 end while

* Evaluate thermal observables with converged $\rho$.

Algorithm 1.1 Finite temperature Hartree-Fock algorithm

Note that in above algorithm, the convergence criteria can also be the density
matrix or MO energies.

For the interacting case, a naive approach is exact diagonalization (ED),
where all eigenstates of the Hamiltonian $\hat{H}$ are explicitly calculated
and the thermal average of an observable $\hat{A}$ is evaluated by

$\langle\hat{O}\rangle(T)=\frac{\sum_{n,i}\langle\psi_{i}^{n}|\hat{O}e^{-(\varepsilon_{i}^{n}-\mu
n)/k_{B}T}|\psi_{i}^{n}\rangle}{\sum_{n,i}\langle\psi_{i}^{n}|e^{-(\varepsilon_{i}^{n}-\mu
n)/k_{B}T}|\psi_{i}^{n}\rangle},$ (6)

where $|\psi_{i}^{n}\rangle$ is the $i$th eigenstate in the Fock space with
$n$ particles. The algorithm of ED is described in Algorithm 1.2. The expense
of ED scales exponentially with the number of orbitals $L$, and thus is only
limited to small systems. For electronic systems with two spins, the maximum
$L$ is $\sim 8$. Therefore, nearly no meaningful calculations can be done with
ED.

$\mathcal{Z}=0,O=0$;

1

for  _$n_{a}$ in $[0,L]$_ do

2 for  _$n_{b}$ in $[0,L]$_ do

3

  1. 1.

Construct Hamiltonian $H(n_{a},n_{b})$;

  2. 2.

Diagonalize $H(n_{a},n_{b})$ to get eigenvalues
$\varepsilon_{i}^{n_{a},n_{b}}$ and eigenstates
$\\{|\psi_{i}^{n_{a},n_{b}}\rangle\\}$;

  3. 3.

Evaluate
$\mathcal{Z}^{n_{a},n_{b}}=\sum_{i}e^{-(\varepsilon_{i}^{n_{a},n_{b}}-\mu(n_{a}+n_{b}))/k_{B}T}$
and
$O^{n_{a},n_{b}}=\sum_{i}e^{-(\varepsilon_{i}^{n_{a},n_{b}}-\mu(n_{a}+n_{b}))/k_{B}T}\langle\psi_{i}^{n_{a},n_{b}}|\hat{O}|\psi_{i}^{n_{a},n_{b}}\rangle$;

  4. 4.

$\mathcal{Z}$ += $\mathcal{Z}^{n_{a},n_{b}}$ ; $O$ += $O^{n_{a},n_{b}}$;

4 end for

5

6 end for

* $\langle O\rangle(T)=O/\mathcal{Z}$

Algorithm 1.2 Finite temperature exact diagonalization

One could reduce the computational cost by only including low-lying states in
the ensemble. Davidson diagonalization [18] and Lanczos algorithm [19] are two
methods that construct a smaller subspace of the Hilbert space containing the
low-lying states. In the Lanczos algorithm, starting with a normalized vector
$|\phi_{0}\rangle$, one could generate a set of orthonormal Lanczos vectors
$\\{|\phi_{m}\rangle,m=0,...,M\\}$ to span the Krylov space
$\\{|\phi_{0}\rangle,\hat{H}|\phi_{0}\rangle,...,\hat{H}^{M}|\phi_{0}\rangle\\}$
with the following steps:

  1. 1.

Apply $\hat{H}$ to $|\phi_{0}\rangle$ and split the resulting vector into
$a_{0}|\phi_{0}\rangle$ and $b_{1}|\phi_{1}\rangle$ with
$|\phi_{1}\rangle\perp|\phi_{0}\rangle$

$\hat{H}|\phi_{0}\rangle=a_{0}|\phi_{0}\rangle+b_{1}|\phi_{1}\rangle,$ (7)

where $a_{0}=\langle\phi_{0}|\hat{H}|\phi_{0}\rangle$ and $b_{1}$ is chosen so
that $|\phi_{1}\rangle$ is normalized.

  2. 2.

Iteratively apply $\hat{H}$ to $|\phi_{i}\rangle,i=1,...,M$ to get

$|\phi_{i}\rangle=b_{i}|\phi_{i-1}\rangle+a_{i}|\phi_{i}\rangle+b_{i+1}|\phi_{i+1}\rangle,$
(8)

where the iteration stops at $i=M$ with $b_{M+1}=0$ or when $b_{i}=0$ with
$i<M$.

  3. 3.

Construct the matrix representation of the Krylov space Hamiltonian as

$H^{\prime}=\begin{bmatrix}a_{0}&b_{1}&0&\cdots&0\\\
b_{1}&a_{1}&b_{2}&\cdots&0\\\ 0&b_{2}&a_{2}&\cdots&0\\\ &&&\ddots&\\\
0&0&0&\cdots&a_{M}\end{bmatrix},$ (9)

where we choose $b_{i}$ to be real numbers.

  4. 4.

Diagonalize the Krylov Hamiltonian $H^{\prime}$ to get the eigenvalues and
eigenvectors in the basis of $\\{|\phi_{i}\rangle,i=0,...,M\\}$. Note that the
$H^{\prime}$ is a tridiagonal matrix, and the typical cost to diagonalize an
$M\times M$ symmetric tridiagonal matrix is $\mathcal{O}(M^{2})$, while the
cost of diagonalizing a random symmetric $M\times M$ matrix is
$\mathcal{O}(M^{3})$.

The quality of the Krylov space depends heavily on the initial state
$|\phi_{0}\rangle$. For instance, if $|\phi_{0}\rangle$ has zero overlap with
the ground state, then the leading part of the trace summation at low
temperature is missing and the result is inaccurate. One could sample initial
states and take the average of the sample to get a better approximation. Note
that the above routine is for a system with fixed particle numbers, so to
fulfill the grand canonical ensemble, one should also sample the Fock spaces
with all possible particle numbers. For low temperature simulation, sampling
particle numbers near the targeted electron number is usually enough. We also
provide a summary of the Davidson algorithm in Appendix 29

####  1.2 Imaginary time evolution

The imaginary time evolution operator is defined as $e^{-\beta\hat{H}}$, where
$\beta$ is called the imaginary time. This approach can be used in both ground
state search and the finite temperature calculations. In the latter case,
$\beta$ has a physical meaning: the inverse temperature $\beta=1/k_{B}T$. At
$\beta=0$ (infinite temperature), the density matrix $\hat{\rho}(\beta=0)$ is
proportional to the identity operator and the system is maximally entangled.
Differentiating $\hat{\rho}(\beta)=e^{-\beta\hat{H}}$ with respect to $\beta$
is described by the Bloch equation

$\frac{\mathrm{d}\hat{\rho}}{\mathrm{d}\beta}=-\hat{H}\hat{\rho}=-\frac{1}{2}(\hat{H}\hat{\rho}+\hat{\rho}\hat{H}),$
(10)

where the last equal sign used $[\hat{H},e^{-\beta\hat{H}}]=0$. The solution
to Eq. (10) can also be written in a symmetrized form

$\hat{\rho}(\beta)=e^{-\beta\hat{H}/2}\hat{\rho}(\beta=0)e^{-\beta\hat{H}/2}.$
(11)

Density matrix quantum Monte Carlo (DMQMC) [20, 21] is an example of the above
approach. We introduce an energy shift $\Delta E$ to the original Hamiltonian
$\hat{H}$, and Eq. (10) turns into

$\frac{\mathrm{d}\hat{\rho}}{\mathrm{d}\beta}=-\frac{1}{2}(\hat{T}\hat{\rho}+\hat{\rho}\hat{T}),$
(12)

where $\hat{T}=\hat{H}-\Delta E\hat{\mathbb{I}}$, and $\Delta E$ is slowly
adjusted to control the population. A similar concept of $\Delta E$ is also
employed in diffusion Monte Carlo (DMC) [22, 23] and full configuration
interaction quantum Monte Carlo (FCIQMC) [24, 25].

The general form of $\hat{\rho}(\beta)$ can be written as a linear combination

$\hat{\rho}(\beta)=\sum_{ij}\rho_{ij}(\beta)|\psi_{i}\rangle\langle\psi_{j}|,$
(13)

where $\\{|\psi\rangle\\}$ forms a complete orthonormal basis of the Hilbert
space. Here we choose $\\{|\psi\rangle\\}$ to be Slater determinants.
$\\{|\psi_{i}\rangle\langle\psi_{j}|\\}$ forms a basis for operators in this
Hilbert space, denoted as $\\{X_{ij}\\}$ for simplicity. Here we introduce a
term "psips" [26, 27]: each psip resides on a particular basis operator
$X_{ij}$ or site $(i,j)$ with "charge" $p_{ij}=\pm 1$. The imaginary time
evolution is divided into $N_{\beta}$ tiny steps: $\tau=\beta/N_{\beta}$. For
each step, DMQMC loops over the sample of psips and perform the following
steps:

  1. 1.

Spawning along columns of the density matrix. Starting from a psip on site
$(i,j)$, calculate the transition probabilities $\frac{1}{2}|T_{ik}|\tau$ to
spawn onto sites $(k,j)$ with $T_{ik}\neq 0$ and $i\neq k$. If the spawning
attempt is accepted, a psip is born on site $(k,j)$ with charge
$q_{kj}=\mathrm{sign}(T_{ik})q_{ij}$.

  2. 2.

Spawning along rows of the density matrix. Repeat the above step to spawn
psips from site $(i,j)$ onto sites $(i,k)$.

  3. 3.

Psips replication and death. Evaluate the diagonal sum $d_{ij}=T_{ii}+T_{jj}$
for site $(i,j)$: if $d_{ij}>0$, a copy of the psip on site $(i,j)$ is added
to the pool with probability $p_{d}=\frac{1}{2}|d_{ij}|\tau$; if $d_{ij}<0$,
the psip on site $(i,j)$ is removed with probability $p_{d}$.

  4. 4.

Annihilation. Pairs of psips on the same site with opposite charges are
removed from the pool.

The distribution of the psips generated by repeating $N_{\beta}$ times the
above procedure provides an approximation of the unnormalized density matrix
at $\beta$. The thermal average of an observable $\hat{O}$ is then calculated
by

$\langle\hat{O}\rangle(\beta)=\frac{\sum_{ij}\bar{q}_{ij}O_{ji}}{\sum_{i}\bar{q}_{ii}},$
(14)

where $\bar{q}$ is an average of density matrices evaluated from a large
number of repeats of the above imaginary time evolution process.

The main concern of the above approach is the size of the density matrix. The
number of independent elements in the density matrix is
$\sim\mathcal{N}(\mathcal{N}+1)/2$, where $\mathcal{N}$ is the Hilbert space
size which grows exponentially with the system size. Even with heavy
parallelization, DMQMC still suffers from considerable computational cost.
Moreover, the accuracy of DMQMC becomes worse as the temperature lowers,
limiting this method to applications for intermediate or high temperature
calculations.

One could circumvent evolving a density matrix by artificially constructing an
enlarged space in which the density matrix of the original system can be
obtained by partial trace from the pure state solution of the enlarged system.
The above approach is called purification [28]. The idea of purification is
the following: suppose a system $\mathcal{S}$ can be bipartitioned into two
smaller systems $\mathcal{A}$ and $\mathcal{B}$; then a state $|\Psi\rangle$
in $\mathcal{S}$ can be written as

$|\Psi\rangle=\sum_{ij}c_{ij}|A_{i}\rangle|B_{j}\rangle,$ (15)

where $\\{|A_{i}\rangle\\}$ and $\\{|B_{i}\rangle\\}$ are orthonormal bases of
$\mathcal{A}$ and $\mathcal{B}$ respectively, and $\sum_{ij}|c_{ij}|^{2}=1$.
The density matrix of the total system is
$\hat{\rho}_{\mathcal{S}}=|\Psi\rangle\langle\Psi|$, and the density matrix of
$\mathcal{A}$ can be obtained by

$\begin{split}\hat{\rho}_{\mathcal{A}}&=\text{Tr}_{\mathcal{B}}\left(\hat{\rho}_{\mathcal{S}}\right)\\\
&=\sum_{k}\langle
B_{k}|\left(\sum_{ij}c_{ij}|A_{i}\rangle|B_{j}\rangle\right)\left(\sum_{i^{\prime}j^{\prime}}c^{*}_{i^{\prime}j^{\prime}}\langle
A_{i^{\prime}}|\langle B_{j^{\prime}}|\right)|B_{k}\rangle\\\
&=\sum_{ii^{\prime}}\left(\sum_{k}c_{ik}c^{*}_{i^{\prime}k}\right)|A_{i}\rangle\langle
A_{i^{\prime}}|\\\ &=\sum_{ii^{\prime}}w_{ii^{\prime}}|A_{i}\rangle\langle
A_{i^{\prime}}|.\end{split}$ (16)

Eq. (16) has the form of a density matrix operator, with matrix elements
$w_{ii^{\prime}}$. The matrix $\mathbf{w}$ has the following properties: (i)
Hermitian; (ii) diagonal terms $w_{ii}=\sum_{k}|c_{ik}|^{2}\geq 0$; and (iii)
$\sum_{i}w_{ii}=1$. Based on the above properties, we confirm that
$\mathbf{w}$ is a density matrix.

Given a density matrix $\hat{\rho}_{\mathcal{A}}$ and basis
$\\{|A_{i}\rangle\\}$, one could also find a set of $\\{|B_{i}\rangle\\}$ to
construct a state $|\Psi\rangle$ such that $\hat{\rho}_{\mathcal{A}}$ can be
derived from the partial trace of $|\Psi\rangle\langle\Psi|$ with
$\\{|B_{i}\rangle\\}$. The above procedure is called purification. Note that
for a system $\mathcal{A}$, there exist more than one purified state
$|\Psi\rangle$, and one could choose certain $\\{|B_{i}\rangle\\}$ and
$|\Psi\rangle$ for their convenience. At infinite temperature, the density
matrix of subspace $\mathcal{A}$ can be written as

$\hat{\rho}_{\mathcal{A}}(\beta=0)=\frac{1}{N_{\mathcal{A}}}\sum_{i}|A_{i}\rangle\langle
A_{i}|,$ (17)

where $N_{\mathcal{A}}$ is the size of $\mathcal{A}$. One could introduce a
set of ancillary orbitals $\\{|\tilde{A}_{i}\rangle\\}$ which are copies of
$\\{|A_{i}\rangle\\}$ and define the purified state as

$|\Psi(\beta=0)\rangle=\frac{1}{\sqrt{N_{\mathcal{A}}}}\sum_{i}|A_{i}\rangle|\tilde{A}_{i}\rangle.$
(18)

It is easy to prove that $\hat{\rho}_{\mathcal{A}}(\beta=0)$ can be derived as
the partial trace of $|\Psi(\beta=0)\rangle\langle\Psi(\beta=0)|$ with
$\\{|\tilde{A}_{i}\rangle\\}$.

Now one could apply imaginary time evolution onto $|\Psi(\beta=0)\rangle$
instead of $\hat{\rho}_{\mathcal{A}}(\beta=0)$,

$|\Psi(\beta)\rangle\propto
e^{-\beta(\hat{H}\otimes\hat{\mathbb{I}})}|\Psi(\beta=0)\rangle,$ (19)

where $\hat{H}$ is the original Hamiltonian on $\mathcal{A}$ and
$\hat{\mathbb{I}}$ is the identity operator on $\tilde{\mathcal{A}}$. The
thermal average of operator $\hat{O}$ in $\mathcal{A}$ is simply evaluated as

$\langle\hat{O}\rangle(\beta)=\langle\Psi(\beta)|\hat{O}\otimes\hat{\mathbb{I}}|\Psi(\beta)\rangle.$
(20)

The most time consuming step in the above procedure is applying
$e^{-\beta\hat{H}}$ onto $|\Psi\rangle$. A commonly accepted way to deal with
$e^{-\beta\hat{H}}$ is Trotter-Suzuki decomposition. Again we divide $\beta$
into $N_{\beta}$ tiny steps $\tau=\beta/N_{\beta}$, and
$e^{-\beta\hat{H}}=\left(e^{-\tau\hat{H}}\right)^{N_{\beta}}$, where we
assumed that $\hat{H}$ does not change with temperature. Suppose $\hat{H}$ can
be decomposed into $\hat{H}=\hat{H}_{1}+\hat{H}_{2}+\cdots+\hat{H}_{n}$,
according to Trotter-Suzuki approximation

$e^{-\tau\hat{H}}=e^{-\tau\hat{H}_{1}/2}e^{-\tau\hat{H}_{2}/2}\cdots
e^{-\tau\hat{H}_{2}/2}e^{-\tau\hat{H}_{1}/2}+\mathcal{O}(\tau^{3}).$ (21)

Another more accurate approach is the 4th order Runge-Kutta (RK4) algorithm,
which is based on solving the differentiation form of the imaginary time
evolution

$\frac{\mathrm{d}|\Psi\rangle}{\mathrm{d}\beta}=-\hat{H}|\Psi\rangle.$ (22)

Let $t_{m}=m\tau$, then one update step in RK4 algorithm is

$|\Psi(t_{m+1})\rangle=|\Psi(t_{m})\rangle=\frac{1}{6}\tau(k_{1}+2k_{2}+2k_{3}+k_{4}),$
(23)

with initial condition $t_{0}=0$ and
$|\Psi(t_{0})\rangle=|\Psi(\beta=0)\rangle$. $k_{i}(i=1,2,3,4)$ are defined
from the $m$th step values

$\begin{split}k_{1}&=-\hat{H}|\Psi(t_{m})\rangle,\\\
k_{2}&=-\hat{H}\left(|\Psi(t_{m})\rangle+\frac{\tau}{2}k_{1}\right),\\\
k_{3}&=-\hat{H}\left(|\Psi(t_{m})\rangle+\frac{\tau}{2}k_{2}\right),\\\
k_{4}&=-\hat{H}\left(|\Psi(t_{m})\rangle+\tau k_{3}\right).\\\ \end{split}$
(24)

The error of one RK4 iteration scales as $\mathcal{O}(\tau^{5})$, and the
accumulated error is $\mathcal{O}(\tau^{4})$.

An example which adopted the purification approach is the finite temperature
density matrix renormalization group (FT-DMRG) [29] algorithm. The matrix
product state (MPS) is defined with alternating physical and ancillary sites,
as shown in Fig. 1. The operators are arranged in the same alternating manner.
The imaginary time evolution routine then follows the same procedure as
previously developed time-evolving block decimation (TEBD) [30, 31].

Figure 1: Structure of the matrix product states used in the purification
approach of the finite temperature density matrix renormalization group
algorithm.

In addition to the examples mentioned above, there exist several other finite
temperature algorithms. Minimally entangled typical thermal states (METTS)
algorithm [32, 33] which will be mentioned in Chapter 5 is another fulfillment
of finite temperature DMRG based on importance sampling. Compared to the
purification approach, METTS requires a smaller bond dimension and the
statistical error decreases as the temperature lowers. However, METTS has only
been applied to spin systems, because the original formulation does not allow
the variation of electron numbers and thus is limited to canonical ensemble.
One could potentially adapt METTS for a grand canonical ensemble by sampling
the electron numbers or introducing a set of initial states which do not
preserve the electron numbers. Determinantal quantum Monte Carlo (DQMC) [34]
and finite temperature auxiliary field quantum Monte Carlo (FT-AFQMC) [35, 36]
are two other finite temperature algorithms based on importance sampling of
Slater determinants. Both of the two QMC methods utilizes Hubbard-Stratonovich
transformation to transform the many-body imaginary time evolution operator to
single-particle operators expressed as free fermions coupled to auxiliary
fields. AFQMC applies a constrained path to alleviate the sign problem, yet
the computational cost is still non-negligible to reach low enough
temperatures with large system sizes. The dynamical mean-field theory (DMFT)
[37, 38] is an embedding method which maps a many-body lattice problem to a
many-body local problem. Since DMFT evaluates the frequencies, it can be
naturally extended to finite temperature calculations with a finite
temperature impurity solver. As most embedding methods, DMFT results are
affected by the finite size effect, and extrapolation to thermodynamic limit
(TDL) is needed to remove the artifact from the finite impurity size. All the
above numerical algorithms have their pros and cons, and one could make their
choices based on the properties of the system and evaluate the results by
careful benchmarking.

###  2 Summary of research

(a) $E_{\text{tot}}=5E_{0}$

(b) $E_{\text{tot}}=\sum_{i=0}^{4}E_{i}$

Figure 2: Evaluating the total energy with density matrix embedding theory.
(a) A hydrogen ring composed of $10$ atoms obeying the periodic boundary
condition, and the impurity (supercell) is two adjacent atoms. The total
energy equals the energy of the supercell times the number of supercells. (b)
A single ligand heme molecule is divided into $5$ non-overlapping fragments,
and the DMET energy of each fragment is calculated. The total energy is the
sum of energies from all fragments.

This thesis provides several tools to study the finite temperature behaviors
of strongly correlated materials. First we will introduce the finite
temperature density matrix embedding theory (FT-DMET) in Chapter 2. FT-DMET,
as a thermal extension of ground state DMET (GS-DMET) [39, 40, 41, 42], maps
the many-body lattice thermal problem onto an impurity thermal problem. Same
as in GS-DMET, the system is divided into non-overlapping fragments, which are
defined by a set of local orbitals (LOs). For periodic systems, the fragments
are chosen as supercells and thus all fragments are equivalent. For systems
that do not obey periodicity, extensive observables are evaluated for each
fragment and the total value of the observable is the summation of those from
all fragments. An illustration is shown in Fig. 2 of the above two cases. Note
that in the latter case, one should be careful when evaluating intensive
properties, which should either be defined for a specific fragment or
evaluated from global extensive properties. This real-space partition ensures
that most of the entanglement is retained in the fragment for systems with
short correlation lengths. When treating one fragment, we call this fragment
the impurity and the rest of the fragments the environment.

To further capture the entanglement between the impurity and the environment,
we introduce a term called bath. Bath in DMET is a subspace of the environment
which is directly entangled with the impurity, spanned by a set of basis
called bath orbitals. In strongly correlated systems, the correlation is
highly localized, and the entanglement entropy obeys the area law. One could
imagine that the bath orbitals mostly come from sites adjacent to the
impurity. In practice, the bath orbitals are derived from the Schmidt
decomposition of the total system wavefunction, which is initialized as the
mean-field wavefunction and optimized in a bootstrap manner. A nice property
of GS-DMET is that the number of bath orbitals generated from Schmidt
decomposition is exactly equal to the number of impurity orbitals, with the
assumption that the impurity is much smaller than the environment.

The key issue going from GS-DMET to FT-DMET is that the Schmidt decomposition
no longer works since the system cannot be described by one single
wavefunction. In fact the finite temperature state is described by a density
matrix of the mixed state. Remember that the Schmidt decomposition of a
wavefunction is equivalent to the singular value decomposition (SVD) of the
corresponding density matrix. In FT-DMET algorithm, we start from the mean-
field single-particle density matrix $\rho_{0}$, and apply SVD to the
impurity-environment block to generate a set of bath orbitals, as described in
the theory part in Chapter 2. Note that since the temperature enlarges the
entanglement length, one should expect more bath orbitals to cover all
impurity-environment entanglement than in GS-DMET. To do so, we continue to
apply SVD to the impurity-environment block of powers of $\rho_{0}$ to get the
rest of the bath orbitals. The algorithm is benchmarked with one- and two-
dimensional Hubbard models, and shows systematically improved accuracy by
increasing bath or impurity size.

In Chapter 3, we further extend the FT-DMET algorithm to handle ab initio
problems. While model systems can be used to reproduce some of the behaviors
and phases in realistic lattices, being able to perform ab initio simulations
is key to achieve a complete understanding of the materials. There are two
technical differences between model systems and ab initio systems: (ii) in
most of the model systems, site basis is used which is naturally localized,
while in ab initio systems the Coulomb interaction is of long range and the
basis set used is usually not localized; (ii) in model systems, the two-body
interaction form is very simple, while the two-body interaction in an ab
initio Hamiltonian is described by a complicated rank-$4$ matrix. The above
two technical difficulties are universal for all ab initio simulations. For ab
initio FT-DMET, one also needs to deal with the large embedding space due to
the size of the supercell and the basis set, which requires necessary
truncation to the bath space. Moreover, a finite temperature impurity solver
that can handle ab initio Hamiltonian efficiently is also crucial for any
meaningful simulations. In Chapter 3, we provide solutions to the above
problems and present the ab initio FT-DMET algorithm. We further use this
algorithm to explore properties and phase transitions of hydrogen lattices.

Chapter 2 and Chapter 3 present an efficient numerical tool to simulate both
strongly-correlated model systems and ab initio systems. The next question to
answer is what order parameters we can use to capture essential thermal
properties and phase transitions at finite temperature. In Chapter 4, we will
study one of the most common but complex phase transitions: metal-insulator
transition (MIT). We argue that compared to the band structure theory which is
widely used to distinguish metal from insulator, electron locality is a more
universal criteria which can be used to detect finite temperature MIT. We
further introduce an order parameter named complex polarization to measure the
locality of electrons and provide a thermofield approach to evaluate finite
temperature complex polarization. The finite temperature complex polarization
formulation provides an easy but well-defined way to characterize MIT in any
periodic materials.

In Chapter 5, several quantum algorithms will be introduced for both ground
state and finite temperature simulations on quantum devices. With the
development of quantum computing technology, especially the hardware, it can
be foreseen that certain categories of difficult problems in classical
simulations can be solved with less effort on a quantum device. The bridge to
connect chemical problems and successful quantum simulations is efficient
quantum algorithms for noisy intermediate-scale quantum (NISQ) devices.
Several quantum algorithms have been developed to carry out quantum chemical
simulations in the past decades, including quantum phase estimation (QPE) [43,
44] and hybrid quantum-classical variational algorithms such as quantum
approximate optimization algorithm (QAOA) [45, 46, 47] and variational quantum
eigensolver (VQE) [48, 49, 50]. While the above algorithms have many
advantages as advertised, they all require quantum or classical resources that
can easily exceed the capacity of current devices. In Chapter 5, the key
quantum algorithm that will be introduced is called quantum imaginary time
evolution (QITE). As mentioned in Section 1, imaginary time evolution is an
efficient algorithm to find the ground state. If the initial state is the
identity density matrix at infinite temperature, then one could evaluate the
density matrix and thus the thermal observables at any temperature.

The conflict of implementing imaginary time evolution on a quantum device is
that the imaginary evolution operator $e^{\beta\hat{H}}$ is a non-unitary
operator, while only unitary operators are allowed on a quantum device. We
present an approach to reproduce a non-unitary operator with a rescaled
unitary operation on an enlarged domain. This approach could be flexibly
performed both exactly and approximately, depending on the computational
resources available. The result is systematically improved and converges
rapidly by increasing the size of the unitary domain. The convergence to the
ground state can be further accelerated by the quantum Lanczos algorithm
(QLanczos). QLanczos constructs a Krylov subspace with the intermediate states
in QITE simulation, and then diagonalizes the Hamiltonian in the subspace
representation to get a better approximation of the ground state. Unlike the
classical Lanczos algorithm mentioned in Section 1 where the Krylov subspace
is spanned by
$\\{|\psi_{0}\rangle,H|\psi_{0}\rangle,...,H^{m}|\psi_{0}\rangle\\}$, the
Krylov space in QLanczos is spanned by
$\\{|\psi_{0}\rangle,e^{2\tau\hat{H}}|\psi_{0}\rangle,...,e^{2m\tau\hat{H}}|\psi_{0}\rangle\\}$.
The Hamiltonian in the quantum Krylov space can be collected from the energy
measurement at each step for free and no additional measurement is needed.

The third algorithm introduced in Chapter 5 is the quantum minimally entangled
typical thermal states (QMETTS) algorithm. While the first two algorithms
(QITE and QLanczos) can be applied to both ground state and finite temperature
calculations, QMETTS is designed in particular for finite temperature
simulations. QMETTS samples a set of minimally entangled thermal states under
the thermal statistics by a repeated imaginary time evolving and then
collapsing onto the product states routine. The advantage of the QMETTS
algorithm is that the imaginary time evolution (fulfilled by QITE) always
starts from a product state, so that the entanglement will not grow too large
even at very low temperature. We present both classical and quantum
simulations on a variety of problems using the above three quantum algorithms
as examples and tests.

##  Chapter 2 Finite temperature density matrix embedding theory

###  3 Abstract

We describe a formulation of the density matrix embedding theory at finite
temperature. We present a generalization of the ground-state bath orbital
construction that embeds a mean-field finite-temperature density matrix up to
a given order in the Hamiltonian, or the Hamiltonian up to a given order in
the density matrix. We assess the performance of the finite-temperature
density matrix embedding on the 1D Hubbard model both at half-filling and away
from it, and the 2D Hubbard model at half-filling, comparing to exact data
where available, as well as results from finite-temperature density matrix
renormalization group, dynamical mean-field theory, and dynamical cluster
approximations. The accuracy of finite-temperature density matrix embedding
appears comparable to that of the ground-state theory, with at most a modest
increase in bath size, and competitive with that of cluster dynamical mean-
field theory.

###  4 Introduction

The numerical simulation of strongly correlated electrons is key to
understanding the quantum phases that derive from electron interactions,
ranging from the Mott transition [51, 52, 53, 54] to high temperature
superconductivity [55, 56, 57]. Consequently, many numerical methods have been
developed for this task. In the setting of quantum lattice models, quantum
embedding methods [58], such as dynamical mean-field theory (DMFT)[59, 60, 61,
62, 63] and density matrix embedding theory (DMET)[39, 40, 42, 41, 64, 65,
66], have proven useful in obtaining insights into complicated quantum phase
diagrams. These methods are based on an approximate mapping from the full
interacting quantum lattice to a simpler self-consistent quantum impurity
problem, consisting of a few sites of the original lattice coupled to an
explicit or implicit bath. In this way, they avoid treating an interacting
quantum many-body problem in the thermodynamic limit.

The current work is concerned with the extension of DMET to finite
temperatures. DMET so far has mainly been applied in its ground-state
formulation (GS-DMET), where it has achieved some success, particularly in
applications to quantum phases where the order is associated with large unit
cells [41, 64, 67]. The ability to treat large unit cells at relatively low
cost compared to other quantum embedding methods is due to the computational
formulation of DMET, which is based on modeling the ground-state impurity
density matrix, a time-independent quantity accessible to a wide variety of
efficient quantum many-body methods. Our formulation of finite-temperature
DMET (FT-DMET) is based on the simple structure of GS-DMET, but includes the
possibility to generalize the bath so as to better capture the finite-
temperature impurity density matrix. Bath generalizations have previously been
used to extend GS-DMET to the calculation of spectral functions and other
dynamical quantities [68, 69]. Analogously to GS-DMET, since one only needs to
compute time-independent observables, finite-temperature DMET can be paired
with the wide variety of quantum impurity solvers which can provide the
finite-temperature density matrix.

We describe the theory of FT-DMET in Section 5. In Section 6 we carry out
numerical calculations on the 1D and 2D Hubbard models, using exact
diagonalization (ED) and the finite-temperature density matrix renormalization
group (FT-DMRG) [70] as quantum impurity solvers. We benchmark our results
against those from the Bethe ansatz in 1D, and DMFT and the dynamical cluster
approximation (DCA) in 2D, and also explore the quantum impurity derived Néel
transition in the 2D Hubbard model. We finish with brief conclusions about
prospects for the method in 7.

###  5 Theory

####  5.1 Ground state DMET

In this work, we exclusively discuss DMET in lattice models (rather than for
_ab initio_ simulations [42, 40, 66, 71]). As an example of a lattice
Hamiltonian, and one that we will use in numerical simulations, we define the
Hubbard model [72, 73],

$\hat{H}=-t\sum_{\langle
i,j\rangle,\sigma}\hat{a}^{\dagger}_{i\sigma}\hat{a}_{j\sigma}-\mu\sum_{i,\sigma}\hat{a}^{\dagger}_{i\sigma}\hat{a}_{i\sigma}+U\sum_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}$
(25)

where $\hat{a}^{\dagger}_{i\sigma}$ creates an electron with spin $\sigma$ on
site $i$ and $\hat{a}_{i\sigma}$ annihilates it;
$\hat{n}_{i\sigma}=\hat{a}^{\dagger}_{i\sigma}\hat{a}_{i\sigma}$; $t$ is the
nearest-neighbour (denoted $\langle i,j\rangle$) hopping amplitude, here set
to $1$; $\mu$ is a chemical potential; and $U$ is the on-site repulsion.

The general idea behind a quantum embedding method such as DMET is to
approximately solve the interacting problem in the large lattice by dividing
the lattice into small fragments or impurities [58]. (Here we will assume that
the impurities are non-overlapping). The main question is how to treat the
coupling and entanglement between the impurities. In DMET, other fragments
around a given impurity are modeled by a set of bath orbitals. The bath
orbitals are constructed to exactly reproduce the entanglement between the
impurity and environment when the full lattice is treated at a mean-field
level (the so-called “low-level” theory). The impurity together with its bath
orbitals then constitutes a small embedded quantum impurity problem, which can
be solved with a “high-level” many-body method. The low-level lattice
wavefunction and the high-level embedded impurity wavefunction are made
approximately consistent, by enforcing self-consistency of the single-particle
density matrices of the impurities and of the lattice. This constraint is
implemented by introducing a static correlation potential on the impurity
sites into the low-level theory. The correlation potential introduced in DMET
is analogous to the DMFT self-energy. A detailed discussion of the correlation
potential including the comparison to other approaches such as density
functional theory (DFT) can be found in [58, 40].

To set the stage for the finite-temperature theory, in the following we
briefly recapitulate some details of the above steps in the GS-DMET
formulation. In particular, we discuss how to extract the bath orbitals, how
to construct the embedding Hamiltonian, and how to carry out the self-
consistency between the low-level and high-level methods. Additional details
for the GS-DMET algorithm can be found in several articles [39, 41, 42],
including the review in Ref. [42].

#####  5.1.1 DMET bath construction

Given a full lattice of $L$ sites, we define the impurity $x$ over $L_{x}$
sites, the Hilbert space of which is denoted as $\mathcal{A}^{x}$ and spanned
by a set of orthonormal basis $\\{|A^{x}_{i}\rangle\\}$. The rest of the
lattice is treated as the environment of impurity $x$, the Hilbert space of
which is denoted as $\mathcal{E}^{x}$ spanned by an orthonormal basis
$\\{|E^{x}_{i}\rangle\\}$. The Hilbert space of the entire lattice
$\mathcal{H}$ is the direct product of the two subsystem Hilbert spaces:
$\mathcal{H}=\mathcal{A}^{x}\otimes\mathcal{E}^{x}$. Any state $|\Psi\rangle$
in $\mathcal{H}$ can be written as

$|\Psi\rangle=\sum_{ij}\psi_{ij}|A^{x}_{i}\rangle|E^{x}_{j}\rangle,$ (26)

where the coefficients $\psi_{ij}$ form a $2^{n_{A}}\times 2^{n_{E}}$ matrix.
Absorbing $\psi_{ij}$ into the environment orbitals, one could rewrite Eq.
(26) as

$\begin{split}|\Psi\rangle&=\sum_{i}|A^{x}_{i}\rangle\left(\sum_{j}\psi_{ij}|E^{x}_{j}\rangle\right)\\\
&=\sum_{i}|A^{x}_{i}\rangle|B^{x}_{i}\rangle\end{split},$ (27)

where $|B^{x}_{i}\rangle=\sum_{j}\psi_{ij}|E^{x}_{j}\rangle$. Eq. (27) tells
us that the orbitals in $\mathcal{E}^{x}$ that are entangled to the impurity
$x$ are of the same size as the impurity orbitals. Note that
$\\{|B^{x}_{i}\rangle\\}$ are not orthonormal and the rest of the environment
enters as a separatable product state $|\Psi_{\text{core}}\rangle$ called
"core contribution". Let $\\{|\tilde{B}^{x}_{i}\rangle\\}$ denote the
orthonormal states derived from $\\{|B^{x}_{i}\rangle\\}$, then Eq. (27) can
be rewritten as

$|\Psi\rangle=\left(\sum_{i}\lambda_{i}|A^{x}_{i}\rangle|\tilde{B}^{x}_{i}\rangle\right)|\Psi_{\text{core}}\rangle.$
(28)

The orbitals $\\{|\tilde{B}^{x}_{i}\rangle\\}$ are directly entangled with the
impurity $x$, and thus are called bath orbitals. The space spanned by impurity
and bath is called embedding space. One can then derive the embedding state as

$|\Psi_{\text{emb}}\rangle=\sum_{i}\lambda_{i}|A^{x}_{i}\rangle|\tilde{B}^{x}_{i}\rangle.$
(29)

If $|\Psi\rangle$ is an eigenstate of the Hamiltonian $\hat{H}$ in the full
lattice, then one can prove that $|\Psi_{\text{emb}}\rangle$ is also an
eigenstate of the embedding Hamiltonian $\hat{H}_{\text{emb}}$ defined as the
projection of $\hat{H}$ onto the embedding space. The two eigenvalues are
identical. Therefore, the full lattice problem can be reduced to a smaller
embedding problem.

In practice, the exact bath orbitals are unknown since the many-body
eigenstate $|\Psi\rangle$ is the final target of the calculation. Instead, we
construct a set of approximated bath orbitals from a mean-field ("low-level")
wavefunction $|\Phi\rangle$, which is an eigenstate of a quadratic lattice
Hamiltonian $\hat{h}$. We rewrite $|\Phi\rangle$ according to Eq. (28) and Eq.
(29) in the form

$\displaystyle|\Phi\rangle=|\Phi_{\text{emb}}\rangle|\Phi_{\text{core}}\rangle.$
(30)

The single-particle density matrix $D^{\Phi}$ obtained from $|\Phi\rangle$
contains all information on the correlations in $|\Phi\rangle$. Thus the bath
orbitals can be defined from this density matrix.

We consider the impurity-environment block $D^{\Phi}_{\text{imp-env}}$
($D_{ij}$ for $i\in x,j\notin x$) of dimension $L^{x}\times(L-L^{x})$. Then
taking the thin SVD

$D^{\Phi}_{\text{imp-env}}=U\lambda B^{\dagger},$ (31)

the columns of $B$ specify the bath orbitals in the lattice basis. The bath
space is thus a function of the density matrix, denoted $B(D)$.

#####  5.1.2 Embedding Hamiltonian

After obtaining the bath orbitals, we construct the embedded Hamiltonian of
the quantum impurity problem. In GS-DMET, there are two ways to do so: the
interacting bath formulation and the non-interacting bath formulation. The
conceptually simplest approach is the interacting bath formulation. In this
case, we project the interacting lattice Hamiltonian $\hat{H}$ into the space
of the impurity plus bath orbitals, defined by the projector $\hat{P}$, i.e.
the embedded Hamiltonian is $\hat{H}_{\text{emb}}=\hat{P}\hat{H}\hat{P}$.
$\hat{H}_{\text{emb}}$ in general contains non-local interactions involving
the bath orbitals, as they are non-local orbitals in the environment. From the
embedded Hamiltonian, we compute the high-level ground-state impurity
wavefunction,

$\displaystyle\hat{H}_{\text{emb}}|\Psi\rangle=E|\Psi\rangle.$ (32)

If $\hat{H}$ were itself the quadratic lattice Hamiltonian $\hat{h}$, then
then $\Psi=\Phi$ and

$\hat{P}\hat{h}\hat{P}|\Phi\rangle=E|\Phi\rangle.$ (33)

Another way to write Eq. (33) for a mean-field state is

$\displaystyle[{P}{h}{P},{P}{D}^{\Phi}{P}]=0,$ (34)

where $h$ denotes the single-particle Hamiltonian matrix and $P$ is the
single-particle projector into the impurity and bath orbitals. These
conditions imply that the lattice Hamiltonian and the embedded Hamiltonian
$\hat{H}_{\text{emb}}$ share the same ground-state at the mean-field level,
which is the basic approximation in GS-DMET.

In the alternative non-interacting bath formulation, interactions on the bath
are approximated by a quadratic correlation potential (discussed below). This
formulation retains the same exact embedding property as the interacting bath
formulation for a quadratic Hamiltonian. In practice, both formulations give
similar results in the Hubbard model [65, 74], and the choice between the two
depends on the available impurity solvers; the interacting bath formulation
generates non-local two-particle interactions in the bath that not all
numerical implementations can handle. In this work, we use the interacting
bath formulation in the 1D Hubbard model where an ED solver is used. In the 2D
Hubbard model, we use the non-interacting bath formulation, where both ED and
FT-DMRG solvers are used. This latter choice is because the cost of treating
non-local interactions in FT-DMRG is relatively high (and we make the same
choice with ED solvers to keep the results strictly comparable).

#####  5.1.3 Self-consistency

To maintain self-consistency between the ground-state of the lattice mean-
field $|\Phi\rangle$, and that of the interacting embedded Hamiltonian
$|\Psi\rangle$, we introduce a quadratic correlation potential $\hat{u}$ into
$h$, i.e.

$\displaystyle\hat{h}\to\hat{h}+\hat{u},$ (35)

where $\hat{u}$ is constrained to act on sites in the impurities, i.e.
$\hat{u}=\sum_{x}\hat{u}^{x}$. To study magnetic order, we choose the form

$\displaystyle\hat{u}^{x}=\sum_{ij\in
x,\sigma\in\\{\uparrow,\downarrow\\}}u^{x}_{ij\sigma}a^{\dagger}_{i\sigma}a_{j\sigma}.$
(36)

The coefficients $u^{x}_{ij\sigma}$ are adjusted to match the density matrices
on the impurity that are evaluated from the low-level wavefunction
$|\Phi\rangle$ and from the high-level embedded wavefunction $|\Psi\rangle$.
In this work, we only match the single-particle density matrix elements of the
impurity (impurity-only matching [42]) by minimizing the cost function:

$f(u)=\sum_{i,j\in\text{imp}}(D_{ij}^{\text{low}}-D_{ij}^{\text{high}})^{2},$
(37)

where $D^{\text{low}}$ and $D^{\text{high}}$ are single-particle density
matrices of low-level and high-level solutions, respectively. For each
minimization iteration, we assume that the high-level single-particle density
matrix is fixed, and the gradient of Eq. (37) is

$\frac{\mathrm{d}f}{\mathrm{d}u}=\sum_{i,j\in\text{imp}}2(D_{ij}^{\text{low}}-D_{ij}^{\text{high}})\frac{\mathrm{d}D_{ij}^{\text{low}}}{\mathrm{d}u}.$
(38)

For ground state, Ref. [42] provided an analytical approach to evaluate
$\frac{\mathrm{d}D_{ij}^{\text{low}}}{\mathrm{d}u}$ using the first order
perturbation theory. At finite temperature, one could also evaluate the
gradient analytically as shown in the Appendix.

Note also that we will only be considering translationally invariant systems,
and thus $\hat{u}^{x}$ is the same for all impurities.

####  5.2 Ground-state expectation values

Ground-state expectation values are evaluated from the density matrices of
each high-level impurity wavefunctions $|\Psi^{x}\rangle$. Since there are
multiple impurities (in a translationally invariant system, these are
constrained to be identical), an expectation value is typically assembled from
the multiple impurity wavefunctions using a democratic partitioning [42]. For
example, given two sites $i$, $j$, where $i$ is part of impurity $x$ and $j$
is part of impurity $y$, the expectation value of $a^{\dagger}_{i}a_{j}$ is

$\displaystyle\langle
a^{\dagger}_{i}a_{j}\rangle=\frac{1}{2}[\langle\Psi^{x}|a^{\dagger}_{i}a_{j}|\Psi^{x}\rangle+\langle\Psi^{y}|a^{\dagger}_{i}a_{j}|\Psi^{y}\rangle].$
(39)

Note that the pure bath components of the high-level wavefunctions, e.g.
$\langle\Psi^{x}|a^{\dagger}_{i}a_{j}|\Psi^{x}\rangle$ for $i,j\notin x$ are
not used in defining the DMET expectation values. Instead, the democratic
partitioning is arranged such that an individual impurity embedding
contributes the correct amount to a global expectation value so long as the
impurity wavefunction produces correct expectation values for operators that
act on the impurity alone, or the impurity and bath together.

####  5.3 Finite temperature DMET

Our formulation of FT-DMET follows the same rubric as the ground-state theory:
a low-level (mean-field-like) finite-temperature density matrix is defined for
the lattice; this is used to obtain a set of bath orbitals to define the
impurity problem; a high-level finite-temperature density matrix is calculated
for the embedded impurity; and self-consistency is carried out between the two
via a correlation potential. The primary difference lies in the properties of
the bath, which we focus on below, as well as in the appearance of quantities
such as the entropy, which are formally defined from many-particle expectation
values.

#####  5.3.1 Finite temperature bath construction

In GS-DMET bath construction, the bath orbitals are directly defined from
Schmidt decomposition of the full lattice ground state wavefunction as in Eq.
(28). However, at finite temperature, the state of an open quantum system
(grand canonical ensemble) is described by a mixed state: the density matrix
is described by a linear combination of pure state density matrices. As a
consequence, the Schmidt decomposition can no longer be used to define bath
orbitals. In fact, with non-zero temperature, the entanglement becomes more
delocalized. To capture the entanglement between the impurity and environment,
a larger bath space is needed compared to that of ground state. In Fig. 3, we
plotted the weight of entanglement with the impurity as a function of distance
(in sites) from the impurity for a $100$-site tight binding model
($\hat{H}=\sum_{i}\hat{a}_{i}^{\dagger}\hat{a}_{i+1}+\text{h.c.}$). One could
see that as the temperature rises, more and more farther sites are entangled
with the impurity, and eventually all sites are uniformly and maximumly
entangled. At ground state ($T=0$), the weight decreased with distance with an
oscillating manner, with wavelength = $2$ sites; at $T=0.15$, the wavelength
increased to $6$ sites due to the smearing effect of finite temperature. The
increase of oscillating wavelength is another example of the increase of
correlation length with temperature.

Figure 3: Weight of the entanglement with the impurity on environmental sites.
The weights are evaluated as the square norm of projections of environmental
sites to the bath space. Due to periodicity of the system, only half of the
environmental sites are shown in the figure.

The difficulty of finite temperature bath orbital construction can also be
demonstrated by the commutation relation between the projected single-particle
density matrix and projected Hamiltonian. In GS-DMET, the bath orbital
construction is designed to be exact if all interactions are treated at the
mean-field level, giving rise to the commuting condition for the projected
single-particle density matrix and projected Hamiltonian in Eq. (34). At
finite temperature, the above commuting condition does not stand and one
should expect approximated bath orbitals even at mean-field level. In general,
we can look for a finite-temperature bath construction that preserves a
similar property. As pointed out in Sec. 5.2, the DMET embedding is still
exact for single-particle expectation values if the embedded projected single-
particle density matrix produces the correct expectation values in the
impurity and impurity-bath sectors, due to the use of the democratic
partitioning. We aim to satisfy this slightly relaxed condition.

The finite temperature single-particle density matrix of a quadratic
Hamiltonian $\hat{h}$ is given by the Fermi-Dirac function

${D}(\beta)=\frac{1}{1+e^{({h}-\mu)\beta}},$ (40)

where $\beta=1/k_{B}T$ ($k_{B}$ is the Boltzmann constant, $T$ is the
temperature). In the following, we fix $k_{B}=1$, thus $\beta=1/T$. If we
could find an embedding directly analogous to the ground-state construction,
we would obtain a projector ${P}$, such that the embedded density matrix
${P}{D}{P}$ is the Fermi-Dirac function of the embedded quadratic Hamiltonian,
i.e. ${P}h{P}$, i.e.

${P}{D}{P}=\frac{1}{1+e^{({P}{h}{P}-\mu)\beta}}.$ (41)

However, unlike in the ground-state theory, the non-linearity of the
exponential function means that Eq. (41) can only be satisfied exactly if
${P}$ projects back into the full lattice basis. Thus a bath orbital
construction at finite temperature is necessarily always approximate, even for
quadratic Hamiltonians.

Nonetheless, one can choose the bath orbitals to reduce the error between the
l.h.s. and r.h.s. in Eq. (41). First, we require that the equality is
satisfied only for the impurity-environment block of $D$, following the
relaxed requirements of the democratic partitioning. Second, we require the
equality to be satisfied only up to a finite order $n$ in $h$, i.e.

$\displaystyle[PDP]_{ij}=\left[\frac{1}{1+e^{({P}{h}{P}-\mu)\beta}}\right]_{ij}+O(h^{n})\quad
i\in x,j\notin x.$ (42)

Then there is a simple algebraic construction of the bath space as (see
Appendix for a proof)

$\displaystyle\\{B(h)\oplus B(h^{2})\oplus B(h^{3})\ldots B(h^{n})\\},$ (43)

where $B(h^{k})$ is the bath space derived from $h^{k}$, $k=1,...,n$. Note
that each order of $h$ adds $L_{x}$ bath orbitals to the total impurity plus
bath space.

We can alternatively choose the bath to preserve the inverse relationship
between the density matrix and Hamiltonian,

$\displaystyle[PhP]_{ij}=\mathrm{inverseFD}(PDP)+O(D^{n})\quad\mathrm{not}\
i,j\notin x,$ (44)

where $\mathrm{inverseFD}$ is the inverse Fermi-Dirac function, and the bath
space is then given as

$\displaystyle\\{B(D)\oplus B(D^{2})\oplus B(D^{3})\ldots B(D^{n})\\}.$ (45)

The attraction of this construction is that the lowest order corresponds to
the standard GS-DMET bath construction.

The above generalized bath constructions allow for the introduction of an
unlimited number of bath sites (so long as the total number of sites in the
embedded problem is less than the lattice size). Increasing the size of the
embedded problem by increasing the number of bath orbitals (hopefully)
increases the accuracy of the embedding, but it also increases the
computational cost. However, an alternative way to increase accuracy is simply
to increase the number of impurity sites. Which strategy is better is problem
dependent, and we will assess both in our numerical experiments.

#####  5.3.2 Thermal observables

The thermal expectation value of an observable $\hat{O}$ is defined as

$\langle\hat{O}(\beta)\rangle=\Tr\left[\hat{\rho}(\beta)\hat{O}\right].$ (46)

Once $\hat{\rho}(\beta)$ is obtained from the high-level impurity calculation,
for observables based on low-rank (e.g. one- and two-) particle reduced
density matrices, we evaluate Eq. (46) using the democratic partitioning
formula for expectation values in Sec. 5.2.

We will also, however, be interested in the entropy per site, which is a many-
particle expectation value. Rather than computing this directly as an
expectation value, we will obtain it by using the thermodynamic relation
$\frac{\dd S}{\dd E}=\beta$, and

$S(\beta_{0})=S(0)+\int_{E(0)}^{E(\beta_{0})}\beta(E)\dd E$ (47)

where $\beta_{0}$ is the desired inverse temperature, and $S(0)=\ln 4\approx
1.386$.

###  6 Results

####  6.1 Computational details

We benchmarked the performance of FT-DMET in the 1D and 2D Hubbard models as a
function of $U$ and $\beta$. For the 1D Hubbard model, we compared our FT-DMET
results to exact solutions from the thermal Bethe ansatz [75]. For the 2D
Hubbard model, the FT-DMET results were compared to DCA and DMFT results [76,
2, 1, 3, 4, 77]. We used large DMET mean-field lattices with periodic boundary
conditions (240 sites in 1D, $24\times 24$ sites in 2D). We used exact
diagonalization (ED) and finite temperature DMRG (FT-DMRG) as impurity
solvers. There are two common ways to carry out FT-DMRG calculations: the
purification method [70] and the minimally entangled typical thermal states
(METTS) method [78]. In this work, we used the purification method implemented
with the ITensor package [79] as the FT-DMRG impurity solver, as well as to
provide the finite lattice reference data in Fig.8. In the FT-DMRG solver, the
sites were ordered with the impurity sites coming first, followed by the bath
sites (an orthonormal basis for the set of bath sites of different orders was
constructed via singular value decomposition, and ordered in decreasing size
of the singular values) and the ancillae arranged in between each physical
site. In the 1D Hubbard model, we used ED exclusively and the interacting bath
formulation of DMET, while in the 2D Hubbard model, we used ED for the 4
impurity, 4 bath calculations, and FT-DMRG for the 4 impurity, 8 bath
calculations, both within the non-interacting bath formulation. FT-DMRG was
carried out using 4th order Runge-Kutta time evolution. To denote different
calculations with different numbers of impurity and bath orbitals, we use the
notation $InBm$, where $n$ denotes the number of impurity sites and $m$ the
number of bath orbitals.

####  6.2 1D Hubbard model

The 1D Hubbard model is an ideal test system for FT-DMET as its thermal
properties can be exactly computed via the thermal Bethe ansatz. We thus use
it to assess various choices within the FT-DMET formalism outlined above.

We first compare the relative performance of the two proposed bath
constructions, generated via the Hamiltonian in Eq. (43) or via the density
matrix in Eq. (45). In Fig. 4, we show the error in the energy per site
(measured from the thermal Bethe ansatz) for $U=2,4$ and half-filling for
these two choices. (The behaviour for other $U$ is similar). Using 4 bath
sites, the absolute error in the energy is comparable to that of the ground-
state calculation (which uses 2 bath sites) over the entire temperature range.
Although the Hamiltonian construction was motivated by the high temperature
expansion of the density matrix, the density matrix construction appears to
perform well at both low and high temperatures. Consequently, we use the
density matrix derived bath in the subsequent calculations.

Figure 4: Error in energy per site (units of $t$) of FT-DMET for the 1D
Hubbard model at $U=2$ and $U=4$ (2 impurity sites and half-filling) with bath
orbitals generated via the density matrix $\gamma$ (Eq. (45)) (blue lines) or
lattice Hamiltonian $h$ (Eq. (43)) (orange lines) as a function of inverse
temperature $\beta$. The numbers in parentheses denote the number of bath
orbitals. The grey area denotes the ground state error with 2 impurity
orbitals.

We next examine the effectiveness of the density matrix bath construction in
removing the finite size error of the impurity. As a first test, in Fig. 5 we
compare the energy error obtained with FT-DMET and $I2B2$ with a pure ED
calculation with 4 impurity sites ($I4$) and periodic (PBC) or antiperiodic
(APBC) boundary conditions, at various $U$ and $\beta$. For weak ($U=2$) to
moderate ($U=4$) coupling, FT-DMET shows a significant improvement over a
finite system calculation with the same number of sites, reducing the error by
a factor of $\sim 2-6$ depending on the $\beta$, thus demonstrating the
effectiveness of the bath. The maximum FT-DMET energy error is 8.1, 6.6, 3.1%
for $U=2,4,8$. At very strong couplings, the error of the finite system ED
with PBC approaches that of FT-DMET. This is because both the finite size
error and the effectiveness of the DMET bath decrease as one approaches the
atomic limit.

Figure 5: Percentage error of the FT-DMET (with 2 impurity sites and 2 bath
orbitals) energy per site vs. ED (4 sites) on a non-embedded cluster with PBC
and APBC boundary conditions for the 1D Hubbard model at various $U$ and
$\beta$.

As a second test, in Fig. 6 we compare increasing the number of impurity sites
versus increasing the number of bath orbitals generated in Eq. (45) for
various $U$ and $\beta$. Although complex behaviour is seen as a function of
$\beta$, we roughly see two patterns. For certain impurity sizes, (e.g. $I4$)
it can be slightly more accurate to use a larger impurity with an equal number
of bath sites, than a smaller impurity with a larger number of bath sites.
(For example, at $U=8$, one can find a range of $\beta$ where $I4B4$ gives a
smaller error than $I2B6$). However, there are also some impurity sizes which
perform very badly; for example $I3B3$ gives a very large error, because the
(short-range) antiferromagnetic correlations do not properly tile between
adjacent impurities when the impurities are of odd size. Thus, due to these
size effects, convergence with impurity size is highly non-monotonic, but
increasing the bath size (by including more terms in Eq. (45)) is less prone
to strong odd-size effects. The ability to improve the quantum impurity by
simply increasing the number of bath sites, is expected to be particularly
relevant in higher-dimensional lattices such as the 2D Hubbard model, where
ordinarily to obtain a sequence of clusters with related shapes, it is
necessary to increase the impurity size by large steps. Nonetheless,
convergence with bath size is also not strictly monotonic, as also illustrated
in Fig. 7, where we see that the error in the $I2B4$ entropy can sometimes be
less than that of $I2B6$ for certain ranges of $\beta$. For the largest
embedded problem $I2B6$, the maximum error in the entropy is $4\times 10^{-3}$
and $2\times 10^{-2}$ for $U=4$ and $8$, respectively.

Figure 6: Absolute error of the FT-DMET energy per site of the 1D Hubbard
model at half-filling as a function of impurity and bath size. $InBm$ denotes
$n$ impurity sites and $m$ bath orbitals. Increasing impurity (blue lines);
increasing bath (orange lines). The grey band depicts the ground state error
with $2$ impurity sites and $2$ bath orbitals. Figure 7: Absolute error of the
FT-DMET entropy per site of the 1D Hubbard model at half-filling as a function
of the number of bath sites. The right panels show the absolute entropy.

The preceding calculations were all carried out at half-filling. Thus, in Fig.
8 we show FT-DMET calculations on the 1D Hubbard model away from half-filling
at $U=4$. We chose to simulate a finite Hubbard chain of 16-sites with PBC in
order to readily generate numerically exact reference data using FT-DMRG
(using a maximum bond dimension of $2000$ and an imaginary time step of
$\tau=0.025$). The agreement between the FT-DMRG energy per site and that
obtained from the thermal Bethe ansatz can be seen at half-filling,
corresponding to a chemical potential $\mu=2$. We see excellent agreement
between FT-DMET and FT-DMRG results across the full range of chemical
potentials, and different $\beta$, suggesting that FT-DMET works equally well
for doped systems as well as for undoped systems.

Figure 8: Energy per site (units of $t$) of a 16-site Hubbard chain with
periodic boundary conditions at $U=4$ as a function of the chemical potential
$\mu$ at various $\beta$ values. The difference between the DMRG and DMET
($I2B6$) energies per site is $0.01-1.4\%$. Solid lines: DMRG energies; dashed
lines: DMET energies; pentagons: Bethe ansatz.

####  6.3 2D Hubbard model

The 2D Hubbard model is an essential model of correlation physics in
materials. We first discuss the accuracy of FT-DMET for the energy of the 2D
Hubbard model at half-filling, shown in Fig. 9. The FT-DMET calculations are
performed on a $2\times 2$ impurity, with 4 bath orbitals ($I4B4$) (green
diamond markers) and 8 bath orbitals ($I4B8$) (red triangular markers). The
results are compared to DCA calculations with clusters of size $34$ (orange
circle markers), $72$ (blue square markers) [77], and $2\times 2$ (light blue
hexagon markers) (computed for this work). The DCA results with the size $72$
cluster can be considered here to represent the thermodynamic limit. The
DCA($2\times 2$) data provides an opportunity to assess the relative
contribution of the FT-DMET embedding to the finite size error; in particular,
one can compare the difference between FT-DMET and DCA(72) to the difference
between DCA($2\times 2$) and DCA(72). Overall, we see that the FT-DMET
energies with 8 bath orbitals are in good agreement with the DCA(72) energies
across the different $U$ values, and that the accuracy is slightly better on
average than that of DCA($2\times 2$). The maximum error in the $I4B8$
impurity compared to thermodynamic limit extrapolations of the DCA energy [77]
is found at $U=4$ and is in the range of 1-2%, comparable to errors observed
in ground-state DMET at this cluster size (e.g. the error in GS-DMET at $U=4$
and $U=8$ is 0.3% and 1.8%, respectively). In the $\beta=8$ case, the FT-DMET
calculations with two different bath sizes give very similar results; at low
temperature, the bath space constructed by the FT procedure is similar to that
of the ground state, and the higher order bath sites do not contribute very
relevant degrees of freedom. Thus even the smaller bath achieves good accuracy
in the low temperature FT-DMET calculations.

Figure 9: Energy per site versus U (units of $t$) of the 2D Hubbard model at
half-filling with FT-DMET ($2\times 2$ cluster with 4 and 8 bath orbitals),
DCA (34, 72 and $2\times 2$ site clusters).

(a)

(b)

Figure 10: Néel transition for the 2D Hubbard model within quantum impurity
simulations. (a) Antiferromagnetic moment $m$ as a function of $T$ with
various $U$ values (units of $t$); (b) Néel temperature $T_{N}$ calculated
with FT-DMET, single-site DMFT and DCA. DMFT data is taken from Ref. [1],
DCA/NCA data for $U=4$ is taken from Ref. [2], DCA/QMC data for $U=6$ is taken
from Ref. [3], and DCA/QMC data for $U=8$ is taken from Ref. [4].

A central phenomenon in magnetism is the finite-temperature Néel transition.
In the thermodynamic limit, the 2D Hubbard model does not exhibit a true Néel
transition, but shows a crossover [80]. However, in finite quantum impurity
calculations, the crossover appears as a phase transition at a nonzero Néel
temperature. Fig. 10(a) shows the antiferromagnetic moment $m$ calculated as
$m=\frac{1}{2L_{x}}\sum_{i}^{L_{x}}|n_{i\uparrow}-n_{i\downarrow}|$ as a
function of temperature $T$ for various $U$ values. As a guide to the eye, we
fit the data to a mean-field magnetization function
$m=a\tanh\left(bm/T\right)$, where $a$ and $b$ are parameters that depend on
$U$. The FT-DMET calculations are performed with a $2\times 2$ impurity and
$8$ bath orbitals, using a finite temperature DMRG solver with maximal bond
dimension $M=600$ and time step $\tau=0.1$. With this, the error in $m$ from
the solver is estimated to be less than 10-3. $m$ drops sharply to zero as $T$
is increased signaling a Néel transition. The Néel temperature $T_{N}$ is
taken as the point of intersection of the mean-field fitted line with the $x$
axis; assuming this form of the curve, the uncertainty in $T_{N}$ is invisible
on the scale of the plot. The plot of $T_{N}$ versus $U$ is shown in Fig.
10(b), showing that the maximal $T_{N}$ occurs at $U=6$. Similar $T_{N}$
calculations on the 2D Hubbard model with single site DMFT [1] and DCA[3, 4,
2] are also shown in Fig. 10(b) for reference. Note that the difference in the
DMFT results [1] and single-site DCA (formally equivalent to DMFT) [3, 4]
likely arise from the different solvers used. The behaviour of $T_{N}$ in our
$2\times 2$ FT-DMET calculations is quite similar to that of the 4-site DCA
cluster. In particular, we see in DCA that the $T_{N}$ values obtained from
calculations with a single-site cluster ($N_{c}=1$) are higher than the
$T_{N}$ values obtained from calculations with a 4-site cluster ($N_{c}=4$).

An alternative visualization of the Néel transition in FT-DMET is shown in
Fig. 11. The FT-DMET calculations here were performed with a $2\times 2$
impurity and $4$ bath orbitals using an ED solver. Though less quantitatively
accurate than the $8$ bath orbital simulations, these FT-DMET calculations
still capture the qualitative behavior of the Néel transition. Focusing on the
dark blue region of the phase diagram, one can estimate the maximal $T_{N}$ to
occur near $U\approx 9$, an increase over the maximal Néel temperature using
the $8$ bath orbital impurity model. This increase in the maximal $T_{N}$
appears similar to that which happens when moving from a 4-site cluster to a
1-site cluster in DCA in Fig. 10.

Figure 11: 2D Hubbard antiferromagnetic moment (color scale) as a function of
$T$ and $U$ (units of $t$) in FT-DMET ($2\times 2$ impurity, 4 bath orbitals.)

###  7 Conclusions

To summarize, we have introduced a finite temperature formulation of the
density matrix embedding theory (FT-DMET). This temperature formulation
inherits most of the basic structure of the ground-state density matrix
embedding theory, but modifies the bath construction so as to approximately
reproduce the mean-field finite-temperature density matrix. From numerical
assessments on the 1D and 2D Hubbard model, we conclude that the accuracy of
FT-DMET is comparable to that of its ground-state counterpart, with at most a
modest increase in size of the embedded problem. From the limited comparisons,
it also appears to be competitive in accuracy with the cluster dynamical mean-
field theory for the same sized cluster. Similarly to ground-state DMET, we
expect FT-DMET to be broadly applicable to a wide range of model and _ab
initio_ problems of correlated electrons at finite temperature [64, 71].

##  Chapter 3 Ab initio finite temperature density matrix embedding theory

###  8 Abstract

This work describes the framework of the finite temperature density matrix
embedding theory (FT-DMET) for ab initio simulations of solids. We introduce
the implementation details including orbital localization, density fitting
treatment to the two electron repulsion integrals, bath truncation, lattice-
to-embedding projection, and impurity solvers. We apply this method to study
the thermal properties and phases of hydrogen lattices. We provide the finite
temperature dissociation curve, paramagnetic-antiferromagnetic transition, and
metal-insulator transition of the hydrogen chain.

###  9 Introduction

The numerical study of the many-electron problem has been playing a profound
role in understanding the electronic behaviors in molecules and materials. One
big challenge for current numerical methods is the description of strong
electron correlations, which requires non-trivial treatment of the electron-
electron interaction beyond the mean-field level. A variety of numerical
algorithms have been invented in the past decades to treat strong electron
correlations, including post-Hartree-Fock quantum chemistry methods such as
CCSD [81, 82], DMRG and its multi-dimensional alternatives [83, 84, 85, 86,
87], the QMC family such as AFQMC [88, 89, 23], and embedding methods such as
DMET [90, 91]. During the past decades, noticeable progress has been made in
the study of strongly correlated models such as one-dimensional and two-
dimensional Hubbard models[92, 93, 94], while the ab initio study of strongly
correlated solids is rare. Compared to model systems where forms of the two-
electron interaction are usually simple, ab initio Hamiltonians contain much
more complicated two-electron terms with size $N^{4}$, where $N$ is the number
of orbitals. This complexity brings higher computational costs. Therefore, an
efficient method that can handle the realistic Hamiltonian accurately is
crucial for understanding the physics behind real materials.

The hydrogen lattice is believed to be the simplest chemical system with a
straightforward analog to the Hubbard model. A thorough comparison between the
hydrogen lattice and Hubbard model could provide insights of the roles of (i)
the long range correlation and (ii) the multi-band effect (with basis set
larger than STO-6G). The numerical study of hydrogen chain can be traced back
to the 70s with simple theoretical tools such as many body perturbation theory
(MBPT)[95]. The rapid development of numerical algorithms made it possible to
achieve a better accuracy and thus plausible conclusions[96, 97, 98, 99, 100,
101, 102, 5]. Motta et al. benchmarked the equation of state[101] and explored
the ground state properties[102] of the hydrogen chain with various popular
numerical methods including DMRG and AFQMC. Despite the numerous ground state
simulations, the finite temperature study of hydrogen lattices is rare, while
the finite temperature study is crucial for understanding the temperature-
related phase diagrams. Liu et. al. studied the finite temperature behaviors
of hydrogen chain with the minimal basis set (STO-6G), and identified the
signature of Pomeranchuk effect[5]. However, the minimal basis set hindered
the exploration of more interesting phenomena caused by the multi-band effect.
A more thorough study beyond the minimal basis set is needed to reach a
quantitative observation of the finite temperature behaviors of the hydrogen
lattices.

In this work, we apply ab initio finite temperature density matrix embedding
theory (FT-DMET) [103] algorithm to study metal-insulator and magnetic
crossovers in periodic one-dimensional and two-dimensional hydrogen lattices
as a function of temperature $T$ and H-H bond distance $R$. We also explore
how basis set size influences the shape of the phases by comparing the results
with STO-6G, 6-31G, and CC-PVDZ basis sets. The rest of the article is
organized as follows: in Section 10, we present the formulation and
implementation details of ab initio FT-DMET, including orbital localization,
tricks to reduce the cost due to the two electron repulsion terms, bath
truncation, impurity solver, and thermal observables. In Section 11, we
demonstrate the ab initio FT-DMET algorithm by studies of the dissociation
curves and phase transitions in a one-dimensional periodic hydrogen lattice.
We finalize this article with conclusions in Section 12.

###  10 Ab initio FT-DMET

In our previous work[103], we introduced the basic formulation of FT-DMET for
lattice models. Going from lattice models to chemical systems, there are
several practical difficulties[104]: (i) the definition of impurity relies on
the localization of the orbitals; (ii) the number of orbitals in the impurity
can be easily very large depending on the infrastructure of the supercell and
the basis set; (iii) manipulating two-electron repulsion integrals in a
realistic Hamiltonian is usually very expensive; (iv) an impurity solver which
can handle ab initio Hamiltonians efficiently at finite temperature is
required. In the rest of this section, we discuss solutions to the above
challenges and provide implementation details of ab initio FT-DMET.

####  10.1 Orbital localization

Since we are dealing with periodic lattices, the whole lattice problem is
described with Bloch (crystal) orbitals in the momentum space. Thus crystal
atomic orbitals (AOs) $\\{\phi^{\textbf{k}}_{\mu}(\textbf{r})\\}$ are a
natural choice. The definition of impurity, however, is based on real space
localized orbitals [105]. Therefore we define a two-step transformation from
Bloch orbitals to localized orbitals (LOs) $\\{w_{i}(\mathbf{r})\\}$.

$\begin{split}w_{i}^{\mathbf{R}}(\mathbf{r})=&\frac{1}{\sqrt{N_{\mathbf{k}}}}\sum_{\mathbf{k}}e^{-i\mathbf{k}\cdot\mathbf{R}}w_{i}^{\mathbf{k}}(\mathbf{r}),\\\
w_{i}^{\mathbf{k}}(\mathbf{r})=&\sum_{\mu}\phi^{\textbf{k}}_{\mu}(\textbf{r})C_{\mu
i}^{\textbf{k}},\end{split}$ (48)

where $C_{\mu i}^{\textbf{k}}$ transforms AOs in momentum space
$\\{\phi^{\textbf{k}}_{\mu}(\textbf{r})\\}$ into LOs in momentum space
$w_{i}^{\mathbf{k}}(\mathbf{r})$, and LOs in real space
$w_{i}^{\mathbf{R}}(\mathbf{r})$ are derived by a Wannier summation over the
local crystal orbitals $w_{i}^{\mathbf{k}}(\mathbf{r})$.

With the ab initio periodic system expressed in LOs, one could choose the
impurity to be spanned by the LOs in a single unit cell or supercell. In the
rest of this paper, we choose the impurity to be the supercell at the lattice
origin.

To define the localization coefficients $C_{\mu i}^{\textbf{k}}$ in Eq. (48),
we use a bottom-up strategy: transform from AO computational basis to LOs.
This strategy uses linear algebra to produce LOs, and thus avoids dealing with
complicated optimizations. There are several choices of LOs from the bottom-up
strategy: Löwdin and meta-Löwdin orbitals [106], natural AOs (NAO) [107, 108],
and intrinsic AOs (IAO) [40]. In this work, we used the $\mathbf{k}$-space
unrestricted Hartree-Fock (KUHF) function with density fitting in the quantum
chemistry package PySCF[109, 110] to generate a set of crystal MOs. Then we
applied an adapted IAO routine to generate a set of crystal IAOs from the
crystal MOs with $\mathbf{k}$-point sampling. The crystal IAOs generated from
this routine are valence orbitals that exactly span the occupied space of the
mean-field calculation. Note that the number of crystal IAOs is equal to the
size of minimal basis only. To carry out calculations beyond the minimal
basis, we construct the rest nonvalence orbitals to be projected AOs (PAOs)
[111], orthogonalized with Löwdin orthogonalization [112]. This IAO+PAO
strategy has been used in previous ground state DMET calculations [42, 104].

####  10.2 Bath truncation and finite temperature bath

In the standard DMET routine, the bath orbitals used to construct the
embedding space are obtained from the SVD of the mean-field off-diagonal
density matrix between the impurity and the remaining lattice (called
environment) $\gamma^{\mathbf{R}\neq\mathbf{0},\mathbf{0}}$

$\gamma^{\mathbf{R}\neq\mathbf{0},\mathbf{0}}=B^{\mathbf{R}\neq\mathbf{0}}\Lambda
V^{\mathbf{0}{\dagger}},$ (49)

where $B^{\mathbf{R}}$ defines the coefficients of bath orbitals in the LO
basis. Thus we can construct the projection matrix in real space

$P^{\mathbf{R}}=\begin{pmatrix}\mathbf{I}&\mathbf{0}\\\
\mathbf{0}&\mathbf{B^{R\neq 0}}\end{pmatrix}.$ (50)

The projection in momentum space can be derived from Eq. (50) with Wannier
transformation

$P^{\mathbf{k}}=\sum_{\mathbf{R}}e^{-i\mathbf{k\cdot R}}P^{\mathbf{R}}.$ (51)

Note that the projection matrices $P^{\mathbf{R}}$ and $P^{\mathbf{k}}$ are in
the basis of LOs, and to obtain the $\mathbf{k}$-space transformation matrix
in AOs, simply multiply $C^{\mathbf{k}}$ from Eq. (48) to the left of
$P^{\mathbf{k}}$.

From Eq. 49, one generates a set of bath orbitals with the same size as the
impurity. This setting is valid and efficient for model systems, and the
embedding space is purely constructed with valence bands. However, for ab
inito calculations, low-lying core and high-energy virtual impurity orbitals
will not entangle with the environment, and thus with the bath orbitals. In
practice, this results in singular values in the SVD of Eq. (49), leading to
difficulties in the convergence of the DMET self-consistency procedure. To
overcome this difficulty, we use the following strategy [42]: we identify the
impurity orbitals as core, valence, and virtual orbitals, and then only take
valence columns of the off-diagonal density matrices of the off-diagonal
density matrix to construct the bath orbitals, as illustrated in Fig. 12. With
this strategy, the size of bath orbitals is equal to the size of valence
impurity orbitals, and thus the number of embedding orbitals is reduced from
$2n_{\text{imp}}$ to $n_{\text{imp}}+n_{\text{val}}$. Note that if
pseudopotential is included in the calculation, there are no core orbitals.

(a)

(b)

Figure 12: Bath orbitals from singular value decomposition (SVD) of the off-
diagonal block of the mean-field density matrix. The whole square represents
the mean-field density matrix of size $N\times N$, with $N$ being the total
number of orbitals, and the first $n_{\text{imp}}$ columns/rows of the matrix
are orbitals in the impurity. (a) Standard DMET routine computes the bath
orbitals by the SVD of the off-diagonal block (orange blocks on the left and
top, the block on the left corresponds to Eq. (49)). (b) Ab initio DMET
computes the bath orbitals by the SVD of only the valence columns in the off-
diagonal blocks (green block).

At finite temperature, electronic occupation numbers of the energy levels are
ruled by the Fermi-Dirac distribution

$f(\varepsilon_{i})=\frac{1}{1+e^{\beta(\varepsilon_{i}-\mu)}},$ (52)

where $\varepsilon_{i}$ is the energy of the $i$th molecular orbital,
$\beta=\frac{1}{T}$ is the inverse temperature (we set the Boltzmann’s
constant $k_{B}=1$) and $\mu$ is the chemical potential or the energy of the
Fermi level at ground state. When $\beta=\infty$, Eq. (52) reproduces the
ground state electron number distribution: when $\varepsilon_{i}<\mu$, the
occupation number is $1$ (occupied orbitals), and when $\varepsilon_{i}>\mu$,
the occupation number is $0$ (virtual orbitals). However, when $\beta$ is
finite, the electronic occupation number on virtual orbitals is no longer $0$.
The extreme case is when $\beta=0$ where all energy levels are uniformly
occupied with occupation number $f=0.5$. Therefore, the ground state bath
construction described previously is no longer suitable to provide an accurate
embedding Hamiltonian. There are generally two strategies: (1) include part of
the core and virtual orbitals into the off-block for SVD; (2) obtain
additional bath orbitals from higher powers of the mean-field density matrix
[103]: take the valence columns of the off-diagonal blocks of $\gamma$,
$\gamma^{2}$, …, $\gamma^{l}$ and apply SVD to them, respectively, to get $l$
sets of bath orbitals, then put the bath orbitals together and perform
orthogonalization to produce the final bath orbitals. The disadvantage of the
first strategy is obvious: as temperature gets higher, the Fermi-Dirac curve
in Fig. 13 gets flatter, and thus more non-valence orbitals are needed.
Compared to the first strategy, the latter strategy generally requires less
number of bath orbitals. For most systems, truncating $l$ to $2$ or $3$ is
already enough for the whole temperature spectrum, therefore the number of
embedding orbitals is $n_{\text{imp}}+ln_{\text{val}}$. Since the number of
valence orbitals is much smaller than that of the non-valence orbitals, DMET
with bath orbitals derived from the second strategy is more economic and
stable. In this paper, we adopt the second strategy for our FT-DMET
calculations.

Figure 13: Fermi-Dirac distribution of electrons on Hartree-Fock molecular
orbitals for H30 chain with STO-6G basis. $T$ is in unit Hartree.

####  10.3 Embedding Hamiltonian

There are two choices of constructing the embedding Hamiltonian: (i)
interacting bath formalism and (ii) non-interacting bath formalism [42]. We
pick the interacting bath formalism to restore most of the two-body
interactions. The embedding Hamiltonian constructed from interacting bath
formalism has the form

$H_{\text{emb}}=\sum{pq}F^{\text{emb}}_{pq}c^{{\dagger}}_{p}c_{q}-\mu\sum_{p\in\text{imp}}c^{{\dagger}}_{p}c_{p}+\frac{1}{2}\sum_{pqrs}\left(pq|rs\right)c^{{\dagger}}_{p}c^{{\dagger}}_{r}c_{s}c_{q}.$
(53)

Note that we use $p,q,r,s$ to index embedding orbitals and $i,j,k,l$ to index
lattice orbitals. A chemical potential $\mu$ is added to only apply on the
impurity, making sure that the number of electrons on the impurity is correct
during the DMET self-consistency.

The embedding Fock matrix $F^{\text{emb}}$ is obtained by projecting the
lattice Fock in AOs to the embedding orbitals. Using
$\tilde{P}^{\mathbf{k}}=C^{\mathbf{k}}P^{\mathbf{k}}$ to denote the projection
operator, one computes the embedding Fock matrix by

$\tilde{F}=\frac{1}{N_{\mathbf{k}}}\sum_{\mathbf{k}}\tilde{P}^{\mathbf{k}\dagger}F^{\mathbf{k}}\tilde{P}^{\mathbf{k}}$
(54)

where $F^{\mathbf{k}}$ is the lattice Fock matrix in $\mathbf{k}$-space AO
basis. To avoid double counting, we subtract the contribution of the embedding
electron repulsion integrals (ERIs) from $\tilde{F}$

$F^{\text{emb}}_{pq}=\tilde{F}_{pq}-\left[\sum_{rs}\left(pq|rs\right)\gamma^{\text{emb}}_{sr}-\frac{1}{2}\left(pr|sq\right)\gamma^{\text{emb}}_{rs}\right]$
(55)

where $\gamma^{\text{emb}}$ is the lattice density matrix rotated to the
embedding space.

The time-consuming part is the construction and projection of the two-electron
ERIs to the embedding space. To reduce the cost, we use density fitting [113,
114] to convert the four-center ERIs to the three-center ERIs,

$\left(\mu\mathbf{k}_{\mu}\nu\mathbf{k}_{\nu}|\kappa\mathbf{k}_{\kappa}\lambda\mathbf{k}_{\lambda}\right)\approx\sum_{L\mathbf{k}_{L}}\left(\mu\mathbf{k}_{\mu}\nu\mathbf{k}_{\nu}|L\mathbf{k}_{L}\right)\left(L\mathbf{k}_{L}|\kappa\mathbf{k}_{\kappa}\lambda\mathbf{k}_{\lambda}\right)$
(56)

where $L\mathbf{k}_{L}$ is the auxiliary basis and only three $\mathbf{k}$
indices are independent due to the conservation of momentum:
$\mathbf{k}_{L}=\mathbf{k}_{\mu}-\mathbf{k}_{\nu}+n\mathbf{q}$ ($n\mathbf{q}$
is the integer multiple of reciprocal lattice vectors). The auxiliary basis
used in this work is a set of chargeless Gaussian crystal orbitals with the
divergent part of the Coulomb term treated in Fourier space [114]. Density
fitting with the above auxiliary basis is called Gaussian density fitting
(GDF). In practice, we first transform three-center ERIs from the lattice
orbitals to the embedding orbitals with cost
$\mathcal{O}\left(n_{\mathbf{k}}^{2}n_{L}n_{\text{lat}}n_{\text{emb}^{2}}+n_{\mathbf{k}}^{2}n_{L}n_{\text{lat}}^{2}n_{\text{emb}}\right)$;
then we convert the three-center ERIs back to the four-center ERIs in the
embedding space with cost $n_{\mathbf{k}}n_{L}n_{\text{emb}}^{2}$. The
computational cost is significantly reduced compared to direct transformation
with cost $\mathcal{O}\left(n_{\mathbf{k}^{3}n_{\text{lat}^{5}}}\right)$.

(a) $R=1.5$ Bohr

(b) $R=3.0$ Bohr

Figure 14: Accuracy test on FT-DMRG and LT-DMRG solvers against exact
diagonalization. The label "FT($x$)" stands for FT-DMRG solver with $\tau=x$,
and the label "LT($x,y$)" stands for LT-DMRG solver with $x$ Davidson roots
and $y$ electron deviations from half-filling for both spins.

####  10.4 Impurity solver

An accurate finite temperature algorithm is required as the impurity solver.
In this work, we use homemade finite temperature exact diagonalization (FT-ED)
and finite temperature density matrix renormalization group (FT-DMRG) for
small and large impurity problems, respectively. In particular, there are two
ways to implement FT-DMRG: (1) imaginary time evolution from an enlarged
Hilbert space, also known as the purification approach [29] (referred as FT-
DMRG) ; and (2) using Davidson diagonalization to generate a set of low-energy
levels to be used in the grand canonical statistics (referred as low
temperature DMRG, LT-DMRG). While FT-DMRG can be used for the whole
temperature spectrum, LT-DMRG is especially for low temperature calculations.
Because most of the phase transitions happen at the low temperature regime,
LT-DMRG can provide accurate enough calculations with lower cost compared to
FT-DMRG.

Since FT-DMRG is based on imaginary time evolution from inverse temperature
$\beta=0$, the entanglement grows rapidly as $\beta$ increases, and at low
temperature a bond dimension that is much larger than the ground state bond
dimension is required. Another error source of FT-DMRG is the imaginary time
step $\tau=\beta/N$, where $N$ is the number of time steps. For symmetrized
Trotter-Suzuki approximation, the local trucation error is on the order of
$\mathcal{O}(\tau^{3})$, while the total accumulated error is on the order of
$\mathcal{O}(\tau^{2})$. If the 4th order Runge-Kutta (RK4) method is used,
the local truncation error is on the order of $\mathcal{O}(\tau^{5})$ and the
total accumulated error is on the order of $\mathcal{O}(\tau^{4})$. The FT-
DMRG used in the calculations of this work used the RK4 method. Note that
since the matrix product state (MPS) truncation is applied at every time step,
having a too small $\tau$ will lead to a large accumulation of MPS truncation
errors. Therefore, one needs to choose the $\tau$ value to be not too small to
introduce large MPS truncation errors and not too big to introduce large time
evolution truncation errors. The error source of LT-DMRG is from the
truncation of the grand canonical summation and the number of roots in the
Davidson diagonalization. Generally the ground state bond dimension is enough
for the low temperature calculations with LT-DMRG.

An assessment of the accuracy of FT-DMRG and LT-DMRG solvers is shown in Fig.
14. The embedding system is composed of two impurity orbitals and two bath
orbitals, generated from a $6$-site hydrogen chain with the STO-6G basis at
$R=1.5$ and $3.0$ Bohr at $\beta=20$. Exact diagonalization (ED) is chosen as
the exact reference. In Fig. 14, we try to understand the role of imaginary
time step $\tau$ in FT-DMRG solver and the roles of the number of Davidson
roots and the size of the truncated grand canonical space in the LT-DMRG
solver. At $\beta=20$, the smaller $\tau$ (red lines) gave a smaller error
compared to $\tau=0.2$ (blue lines), and the FT-DMRG results converged at
$M\sim 300$. The errors of LT-DMRG solver do not change too much with the bond
dimension $M$, so $M=100$ is already enough for a $4$-site system. The
accuracy at $R=1.5$ Bohr is generally better than $R=3.0$ Bohr, since larger
$R$ corresponds to stronger correlation. At larger $R$, one needs to include a
larger number of Davidson roots to achieve high enough accuracy with LT-DMRG.
Generally with large enough bond dimension $M$, FT-DMRG could provide more
accurate results. However, when the embedding system is too large to use a
large $M$, one could consider the LT-DMRG method. In DMET calculations, we use
$ED$ solver for $L_{\text{emb}}<8$ embedding problems and use FT-DMRG solver
for larger problems.

####  10.5 Thermal observables

In order to identify the metal-insulator transition and the Néel transition
and explore the mechanism behind the crossings, we compute the following order
parameters: staggered magnetic moment $m$, double occupancy $D$, complex
polarization $Z$, spin-spin correlation functions $\mathcal{C}_{ss}$, and
charge-charge correlation functions $\mathcal{C}_{cc}$.

_Staggered magnetic moment._ The staggered magnetic momentum is calculated as

$m=\frac{1}{N^{\text{imp}}}\sum_{i\in\text{imp}}\frac{|n_{i,\uparrow}-n_{i,\downarrow}|}{2},$
(57)

where $N^{\text{imp}}$ is the total number of H atoms in the impurity
(supercell), and $n_{i,\uparrow}$ and $n_{i,\downarrow}$ are electron numbers
on $i$th atom with up spin and down spin, respectively. Note that if one uses
Bohr magneton $\mu_{B}=\frac{e\hbar}{2m_{e}}=\frac{1}{2}$ as the unit, then
one would drop $2$ in the denominator in Eq. (57). To evaluate
$n_{i,\uparrow}$ on the $i$th atom, we first compute the one-particle impurity
density matrix for up-spin in the IAO basis, and then sum up the diagonal
terms that belong to the $i$th atom. For example, when STO-6G basis is used,
the occupation numbers on $1s$ orbital and $2s$ orbital of atom-$i$ sum up to
the electron density on atom-$i$. $n_{i,\downarrow}$ is evaluated in the same
way from the down-spin one-particle impurity density matrix.

_Double occupancy._ The double occupancy measures the probability of two
electrons with opposite spins occupying the same hydrogen atom, calculated by

$D=\frac{1}{N^{\text{imp}}}\sum_{i\in\text{imp}}\langle\hat{n}_{i\uparrow}\hat{n}_{i,\downarrow}\rangle.$
(58)

Note that the hat on $\hat{n}_{i\uparrow}$ denotes that it is an operator, not
a number, with $n_{i\uparrow}=\langle\hat{n}_{i\uparrow}\rangle$. Since there
are multiple bands on each atom, we expand $\hat{n}_{i\uparrow}$ as

$\hat{n}_{i\uparrow}=\sum_{w}\hat{n}^{w}_{i\uparrow},$ (59)

where $w$ is the index of the bands on the $i$-th atom ( e.g., $1s$, $2s$,
$2p_{x},2p_{y},2p_{z},...$). Therefore, the precise expression of double
occupancy becomes

$D=\frac{1}{N^{\text{imp}}}\sum_{i\in\text{imp}}\sum_{w,w^{\prime}\in
i}\langle\hat{n}^{w}_{i\uparrow}\hat{n}^{w^{\prime}}_{i,\downarrow}\rangle.$
(60)

_Complex polarization._ Complex polarization measures the mobility of
electrons, and thus can be used as an indicator of metal-insulator transition.
The definition of complex polarization on $z$ direction is

$Z=\langle e^{i\frac{2\pi}{L}\hat{z}}\rangle,$ (61)

where $L$ is the chain length and $\hat{z}$ is the location operator in the
$z$-direction. When $Z=0$, electrons are delocalized and the system is
metallic; when $Z=1$, electrons are localized and the system is insulating. At
mean-field level, the ground state is a Slater determinant $|\phi\rangle$ of
occupied orbitals, so the complex polarization is evaluated by

$Z=\langle\phi|e^{i\frac{2\pi}{L}\hat{z}}|\phi\rangle,$ (62)

which is equivalent to

$Z=\text{Det}\left[C_{\text{occ}}^{{\dagger}}e^{i\frac{2\pi}{L}z}C_{\text{occ}}\right],$
(63)

where $C_{\text{occ}}$ represents the coefficients of occupied orbitals.

At finite temperature, we use a thermofield approach [115] from our recent
work (see Chapter 4). We construct the infinite temperature determinant with
an enlarged Hilbert space $\tilde{\phi}$, and thermofield operators of the
Hamiltonian $\tilde{H}$ and position operator $\tilde{z}$. Then the finite
temperature complex polarization is evaluated by

$Z(\beta)=\frac{1}{\mathcal{Z}}\langle\tilde{\phi}|e^{-\beta(\tilde{H})}e^{i\frac{2\pi}{L}\tilde{z}}|\tilde{\phi}\rangle,$
(64)

where
$\mathcal{Z}\langle\tilde{\phi}|e^{-\beta(\tilde{H})}|\tilde{\phi}\rangle$ is
the partition function.

_Spin-spin correlation and charge-charge correlation functions._ We define the
two correlation functions as follows:

$\begin{split}\mathcal{C}^{ss}_{i}=&\langle(\hat{n}_{0}^{\uparrow}-\hat{n}_{0}^{\downarrow})(\hat{n}_{i}^{\uparrow}-\hat{n}_{i}^{\downarrow})\rangle-\langle\hat{n}_{0}^{\uparrow}-\hat{n}_{0}^{\downarrow}\rangle\langle\hat{n}_{i}^{\uparrow}-\hat{n}_{i}^{\downarrow}\rangle,\\\
\mathcal{C}^{cc}_{i}=&\langle(\hat{n}_{0}^{\uparrow}+\hat{n}_{0}^{\downarrow})(n_{i}^{\uparrow}+n_{i}^{\downarrow})\rangle-\langle\hat{n}_{0}^{\uparrow}+\hat{n}_{0}^{\downarrow}\rangle\langle
n_{i}^{\uparrow}+n_{i}^{\downarrow}\rangle,\end{split}$ (65)

where $\hat{n}_{i}^{\sigma}$ is the electron density operator of spin $\sigma$
on site $i$.

###  11 Results

In this section, we show some preliminary results of FT-DMET calculations on
the hydrogen chain system with periodic boundary condition. First, we examine
the basis set effect on a $22$-atom chain with $2$ atoms in the impurity. Fig.
15 shows the magnetic moment at both ground state and $T=0.02$ Hartree
calculated by DMET with STO-6G, 6-31G, and CC-PVDZ basis sets. Paramagnetic-
antiferromagnetic (PM-AFM) transition is observed at both ground state and
$T=0.02$ Hartree. A very interesting behavior of the magnetic moment at ground
state is observed: with STO-6G, the magnetic moment drops at $R>3.0$ Bohr,
while with larger basis sets, this drop did not happen. The reason for the
above behaviors could be due to the loss of entanglement between different
sites at large $R$. Imagine at $R=\infty$, the system should behave as $22$
individual atoms, and one should expect the ground state to be paramagnetic.
With more diffused orbitals (e.g., $2s$ and $2p$ orbitals), however, the
entanglement between different sites decays slower as $R$ increases. Note that
since the impurity size is only $2$ atoms, one only needs to consider the
entanglement between adjacent sites. Once the impurity gets larger, a drop of
magnetic moment with larger basis sets should also be expected. At $T=0.02$
(left panel), the magnetic moment computed with all three basis sets dropped
as $R$ increases, as a consequence of thermal dissipation.

(a)

(b)

Figure 15: Magnetic moment of a $22$-atom chain at ground state (left panel)
and $T=0.02$ Hartree (right panel) with STO-6G, 6-31G, and CC-PVDZ basis sets.

We further show the double occupancy from the above simulation settings in
Fig. 16. A clear change of the gradient of $D$ as a function of $R$ is
observed for both ground state and $T=0.02$ Hartree, indicating a metal to
insulator transition. The transition $R$ is around $1.6\sim 1.8$ Bohr, which
agrees with the transition $R$ of PM-AFM transition, resulting in a PM metal
phase and AFM insulating phase.

(a)

(b)

Figure 16: Double occupancy of a $22$-atom chain at ground state (left panel)
and $T=0.02$ Hartree (right panel) with STO-6G, 6-31G, and CC-PVDZ basis sets.
The insets show a sudden change of the gradient of $D$ as a function of $R$,
indicating metal to insulator transition.

Next we increase the total number of atoms in the hydrogen chain to $50$ atoms
to eliminate the finite size effect of the total system size. STO-6G is used
as the basis set. The impurity is composed of two hydrogen atoms, and solved
by ED. We first present the energy calculations and dissociation curve of the
hydrogen chain, shown in Fig. 17, where the energy per electron at $T=0.05$ is
compared to FT-AFQMC [5] results. The AFQMC calculation used the STO-6G basis
set and a supercell with $10$ atoms and $5$ $k$ points. The two energy curves
predicted the same dissociating trend and equilibrium point ($\sim 1.8$ Bohr).
However, the DMET energies are all lower than the AFQMC energies, which could
be due to the finite impurity size effect.

Figure 17: Dissociation curve of hydrogen chain at $T=0.05$ Hartree compared
to AFQMC. The AFQMC data is extracted from Ref. [5].

We then examine the staggered magnetic moment $m$ as a function of inter-
atomic distance $R$ at ground state, $T=0.02$, $T=0.05$, and $T=0.1$ in Fig.
18. Compared to Fig. 15 where the paramagnetic-antiferromagnetic (PM-AFM)
transition happened around $R=1.6$ Bohr, we observed the PM-AFM transition at
$R=1.0\sim 1.5$ region for $T<0.1$, which could be an effect of the finite
total system size. Although $T=0.02$ is considered as very close to the ground
state, the magnetic moment at $T=0.02$ drops earlier than the ground state
curve as $R$ increases. This behavior is due to the thermal dissipation of the
magnetic order. As $R$ grows larger, the atoms are far apart from each other,
and thus the electron-electron correlation between different sites is weaker,
and eventually not enough to preserve the long-range AFM order, which lead to
the drop of magnetic order at large $R$ as shown in the figure. Even adding a
small temperature, the flip of the spin can happen to destroy the long-range
AFM order.

Figure 18: Staggered magnetic moment of hydrogen chain with periodic boundary
condition at ground state, $T=0.02$, $T=0.05$, and $T=0.1$. The unit of $T$ is
Hartree.

###  12 Conclusion

In this work, we generalized the previously described finite temperature
density matrix embedding theory to study ab initio solid, and employed the
method to study the hydrogen chain problem. Despite the simplicity of the
hydrogen chain lattice compared to other periodic systems, it exhibits a
variety of intriguing behaviors including paramagnetic-antiferromagnetic
transition and metal-insulator transition at both ground state and finite
temperature. At finite temperature, we observed thermal dissipation for the
magnetic order at large inter-atomic distance. We further confirmed the
stabilizing effect from multi-band basis sets. Since this work is not
completely finished, in the future we will apply this finite temperature
algorithm to a larger set of ab initio solids including the two-dimensional
and three-dimensional hydrogen lattices, transition metal oxides, and
challenging systems such as cuprate-based high temperature superconductors.

##  Chapter 4 Finite temperature complex polarization and metal-insulator
transition

###  13 Abstract

Metal-insulator transition is a fundamental yet complicated topic in condensed
matter physics and material science. Band structure theory has been widely
used in explaining why insulators are insulating and why metals are
conducting, but can fail to describe insulation phases caused by strong
correlations or disorder. Electron locality/mobility serves as a more general
indicator of metallic and insulating phases: when the system is metallic, the
electrons are delocalized and can flow freely; when the system is insulating,
the electrons are localized. The standard deviation (or second cumulant
moment) of the electron position is used as the order parameter of the
electron localization, which is directly related to the complex polarization
of the system. The complex polarization is widely accepted as a new indicator
of the metal-insulator transition at ground state: when the complex
polarization equals zero, the second cumulant moment of the position is
diverged and the system sits in the metallic phase; when the complex
polarization is nonzero, the electrons are localized and the insulating phase
wins. In this work, we present the finite temperature formulation of the
complex polarization. We also introduce a thermofield approach to compute the
complex polarization with thermal Slater determinant. We demonstrate how
finite temperature complex polarization works as an indicator of metal-
insulator transition at low temperature with a modified tight binding model
and hydrogen chain system. In the hydrogen chain case, we also compare the
metal-insulator transition with the paramagnetic-ferromagnetic transition,
electron population, and energy gap to study the origin of the insulating and
metallic phases, respectively.

###  14 Introduction

Phase transition happens when a system undergoes a macroscopic change from one
phase to another phase due to the variation of control parameters such as
temperature, magnetic field, chemical substitution, and pressure. Near the
critical point between the two phases, the physical properties of the bulk
changes dramatically with respect to even a minor perturbation in control
parameters. Metal-insulator transition (MIT) is among the most common phase
transitions, yet the microscopic cause and the physics behind the phenomena is
nontrivial. From the elementary physics textbooks, we learned that metals are
conducting when an electrical field is applied, while insulators do not allow
electrons to flow freely. However, this is a rather vague and bipartite
definition, which is not able to answer questions such as (1) what are the
microscopic driving forces for conductivity? (2) what are the causes for MIT?
and (3) how does one characterize MIT?

In the past, the microscopic featurization of insulators and metals are
generally described by the band structure theory [116, 117]. The band
structure theory describes the movement of a single electron in a periodic
solid, with the mean-field effect from the other electrons. According to band
structure theory, if the Fermi level sits in a band gap, then the system is
insulating; if the Fermi level crosses a band, the system is conducting.
However, the band structure theory is based on independent electron
approximation and is only limited to crystalline systems. The insulating
behavior caused by disorder or electron-electron correlation cannot be
captured accurately by the band structure theory [118, 119, 120, 121]. A more
general description is to use the electron localization to distinguish metal
and insulator: when the electrons are _localized_ , the system is insulating,
and when the electrons are _delocalized_ , the system is conducting. A widely
accepted approach to evaluate electron localization is based on the theory of
polarization [122, 123, 124, 125, 126, 127, 128, 129], where the macroscopic
polarization is connected to Berry phase [130]. A more straightforward
indicator of electron localization is the second cumulant moment of the
electron position operator describing the spread of electrons [129, 131, 123,
124, 125]. A value that connects to both the macroscopic polarization and the
second cumulant moment is complex polarization: the phase of the complex
polarization is the Berry phase, while the second cumulant moment can be
evaluated from the modulo of the complex polarization [132, 133, 134].
Moreover, the DC conductivity according to the ground state Kubo formula [135]
is also related to the modulo of the complex polarization.

Ground state complex polarization and its connection to macroscopic
polarization, electron localization, and DC conductivity have been thoroughly
studied and discussed in the past [129, 131, 133]. However, discussion about
finite temperature complex polarization is rare, regardless of the
significance of this parameter. In this work, we introduce the formulation of
finite temperature complex polarization and discuss its relationship with
electron localization. We also present a mean-field level implementation of
finite temperature complex polarization under thermofield theory [136, 137,
138, 115]. In Section 15, we introduce the ground state formulation of complex
polarization and electron localization, where we first use the single particle
case as a simplified example and then generalize the single particle case to
many-body mean-field formulation. In Section 16, we extend the ground state
formulation to the finite temperature version, and present a thermofield
approach to evaluate the complex polarization. In Section 17, we apply the
finite temperature formulation to a modified tight binding model both
analytically and numerically, presenting a preliminary analysis of how complex
polarization provides information of metal-insulator transition. In Section
18, we choose hydrogen chain as an example of computing finite temperature
complex polarization for ab initio systems and explore the temperature-induced
metal-insulator transition in the hydrogen chain. We finalize this article
with a summary and outlook in Section 19.

###  15 Ground state complex polarization and electron localization

The many-body complex polarization $Z_{N}^{(\alpha)}$ was first introduced as
the ground state expectation value of certain unitary many-body operators
[128]. We start by defining a general form of the unitary many-body operator

$\hat{U}(\mathbf{k})=e^{i\mathbf{k}\cdot\hat{\mathbf{R}}},$ (66)

where $\mathbf{k}$ is an arbitrary three-dimension vector and
$\hat{\mathbf{R}}$ is the three-dimensional position operator, with
$\hat{R}^{\alpha}\psi(r_{1},r_{2},r_{3})=r_{\alpha}\psi(r_{1},r_{2},r_{3}),\alpha=1,2,3$.

We introduce three $\mathbf{k}$ vectors with notation
$\mathbf{\kappa}^{(\alpha)}(\alpha=1,2,3)$, defined as

$\kappa^{(\alpha)}_{\beta}=\frac{2\pi}{L}\delta_{\alpha\beta},$ (67)

which can be explicitly written as

$\begin{split}\kappa^{(1)}&=\left(\frac{2\pi}{L},0,0\right),\\\
\kappa^{(2)}&=\left(0,\frac{2\pi}{L},0\right),\\\
\kappa^{(3)}&=\left(0,0,\frac{2\pi}{L}\right).\end{split}$ (68)

The ground state complex polarization is then defined as the expectation
values of the unitary many-body operators with the above three vectors:

$Z_{N}^{(\alpha)}=\langle\Psi_{0}|\hat{U}(\mathbf{\kappa}^{\alpha})|\Psi_{0}\rangle,$
(69)

where $N$ is the number of electrons and $|\Psi_{0}\rangle$ represents the
ground state of the system of interest. The complex polarization
$Z_{N}^{(\alpha)}$ in Eq. (69) can be explicitly written as

$Z_{N}^{(\alpha)}=|Z_{N}^{(\alpha)}|e^{i\gamma^{(\alpha)}_{N}},$ (70)

where $|Z_{N}^{(\alpha)}|\in[0,1]$ is the modulo of $Z_{N}^{(\alpha)}$ and
$\gamma^{(\alpha)}_{N}$ is the phase of $Z_{N}^{(\alpha)}$, referred as the
Berry phase. In this article, we will not discuss the macroscopic
polarization, therefore the phase of Eq. (70) will not be mentioned. In fact,
with a centrosymmetric choice of the origin, the complex polarization will
always remain real.

####  15.1 Electron localization

We start by considering a problem with one particle in a one-dimensional
potential wall. The locality of the particle can be measured by the quadratic
spread, or the second cumulant moment of the position $x$, defined as

$\langle\delta
x^{2}\rangle=\langle\psi|x^{2}|\psi\rangle-\langle\psi|x|\psi\rangle^{2},$
(71)

where $|\psi\rangle$ is the ground state of the particle in a box.
$\langle\delta x^{2}\rangle$ is finite when the state $|\psi\rangle$ is
bounded and diverges for an unbounded state. Let $n(x)=|\psi(x)|^{2}$ be the
electron density, then we can rewrite Eq. (71) as

$\langle\delta x^{2}\rangle=\int_{-\infty}^{\infty}\mathrm{d}x\