Computational Studies and Algorithmic Research of Strongly Correlated Materials
Zhuoran He
TL;DR
This thesis explores computational methods for strongly correlated materials, analyzing phase transitions in real materials using DFT+U and developing a real-time impurity solver with DMRG for nonequilibrium dynamical mean-field theory.
Contribution
It introduces a non-standard DMRG implementation for real-time evolution of impurity models and investigates the impact of orbital ordering on computational complexity.
Findings
DFT+U captures phase transition behaviors and metastable states.
Developed a 4-MPS DMRG method for real-time impurity simulations.
Orbital ordering significantly affects entanglement and computational cost.
Abstract
This PhD thesis conducts a focused study of strongly correlated materials with localized electron orbitals. We have studied two real materials (LuNiO and VO) and one model system, i.e., the Anderson impurity model. The thesis is divided into two parts. In the first part, we use the DFT+ method, to calculate the equilibrium phase transitions of LuNiO as an example of rare-earth nickelates under a substrate-induced strain, and the nonequilibrium phase transitions of VO as an example of narrow-gap Mott insulators under laser-pulse-induced photoexcitations. The effect of adding is manifested in both materials as the change of band structure in response to the change of orbital occupancies of electrons, i.e., the soft band effect. This effect brings about competitions of electrons between different orbitals and gives rise to multiple metastable states. In the…
| Magnetic order | Symmetry | Form of | ISPIN | LDAUTYPE |
|---|---|---|---|---|
| None | – | |||
| Collinear | ||||
| Non-collinear | General |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsTheoretical and Computational Physics · Semiconductor Quantum Structures and Devices
**Computational Studies and Algorithmic
Research of Strongly Correlated Materials **
Zhuoran He
Submitted in partial fulfillment of the
requirements for the degree of
Doctor of Philosophy
in the Graduate School of Arts and Sciences
COLUMBIA UNIVERSITY
2019
© 2019
Zhuoran He
All rights reserved
**Abstract
**
Strongly correlated materials are an important topic of research in condensed matter physics. Other than ordinary solid-state physical systems, which can be well described and analyzed by the energy band theory, the electron-electron correlation effects in strongly correlated materials are far more significant. So it is necessary to develop theories and methods that are beyond the energy band theory to describe their rich and varied behaviors. Not only are there electron-electron correlations, typically the multiple degrees of freedom in strongly correlated materials, such as the charge distribution, orbital occupancies, spin orientations, and lattice structure exhibit cooperative or competitive behaviors, giving rise to rich phase diagrams and sensitive or non-perturbative responses to changes in external parameters such as temperature, strain, electromagnetic fields, etc.
This thesis is divided into two parts. In the first part, we use the density functional theory (DFT) plus Hartree-Fock corrections, i.e., the DFT+ method, to calculate the equilibrium and nonequilibrium phase transitions of LuNiO3 and VO2. The effect of adding is manifested in both materials as the change of band structure in response to the change of orbital occupancies of electrons, i.e., the soft band effect. This effect bring about competitions of electrons between different orbitals by lowering the occupied orbitals and raising the empty orbitals in energy, giving rise to multiple metastable states. In the second part, we study the dynamic mean field theory (DMFT) as a beyond band-theory method. This is a Green’s-function-based theory for open quantum systems. By selecting one lattice site of an interacting lattice model as an open system, the other lattice sites as the environment are equivalently replaced by a set of noninteracting orbitals according to the hybridization function, so the whole system is transformed into an Anderson impurity model (AIM). We studied how we can use the density matrix renormalization group (DMRG) method to perform real-time evolutions of the Anderson impurity model to understand the nonequilibrium dynamics of a strongly correlated lattice system.
We begin in Chapter with an introduction to strongly correlated materials, density functional theory (DFT) and dynamical mean-field theory (DMFT). The Kohn-Sham density functional theory and its plus correction are discussed in detail. We also demonstrate how the dynamical mean-field theory reduces the lattice sites other than the impurity site as a set of noninteracting bath orbitals.
Then in Chapters and , we show material-related studies of LuNiO3 as an example of rare-earth nickelates under substrate strain, and VO2 as an example of a narrow-gap Mott insulator in a pump-probe experiment. These are two types of strongly correlated materials with localized 3 orbitals (for Ni and V). We use the DFT+ method to calculate their band structures and study the structural phase transitions in LuNiO3 and metal-insulator transitions in both materials. The competition between the charge-ordered and Jahn-Teller distorted phases of LuNiO3 is studied at various substrate lattice constants within DFT+. A Landau energy function is constructed based on group theory to understand the competition of various distortion modes of the NiO6 octahedra. VO2 is known for its metal-insulator transition at 68C, above which temperature it’s a metal and below which it’s an insulator with a doubled unit cell. For VO2 in a pump-probe experiment, a metastable metal phase was found to exist in the crystal structure of the equilibrium insulating phase. Our work is to understand this novel metastable phase from a soft-band picture. We also use quantum Boltzmann equation to justify the prethermalization of electrons over the lifetime of the metastable metal, so that the photoinduced transition of VO2 can be understood in a hot electron picture.
Finally, in Chapters and , we show a focused study of building a real-time solver for the Anderson impurity model out of equilibrium using the density matrix renormalization group (DMRG) method, towards the goal of building an impurity solver for nonequilibrium dynamical mean-field theory (DMFT). We study both the quenched and driven single-impurity Anderson models (SIAM) in real time, evolving the wave function written in a form with 4 matrix product states (MPS) in DMRG. For the quenched model, we find that the computational cost is polynomial time if the bath orbitals in the MPSs are ordered in energy. The same energy-ordering scheme works for the driven model in the short driving period regime in which the Floquet-Magnus expansion converges. In the long-period regime, we find that the computational time grows exponentially with the physical time, or the number of periods reached. The computational cost reduces in the long run when the bath orbitals are quasi-energy ordered, which is discussed in further detail in the thesis.
Contents
-
2 Strain control of electronic phase in rare-earth nickelates
-
3 Photoinduced phase transitions in narrow-gap Mott insulators
Acknowledgements
First and foremost, I would like to thank my advisor, Professor Andrew J. Millis in the Department of Physics of Columbia University. It is a pleasant experience to work with him. I have learned a lot from the inspiring and insightful discussions with him, not only about the knowledge, but also about the way to think about a problem and the attitudes of a good researcher. All of the above would be a precious and powerful source to encourage my efforts and enlighten my journey to the future.
Then, I would like to thank Professor Chris A. Marianetti in Department of Applied Physics. I am very grateful for his help and guidance on density functional theory, group theory and my first two projects on strongly correlated materials. I would like to thank Dante Kennes, Seyoung Park, Hyowon Park, Ara Go and Jia Chen. Discussions with them are enjoyable and fruitful. I owe special thanks to Dr. Kennes for his patient guidance in DMRG in my impurity solver project. I have also benefited a lot from discussions with Dr. Hanghui Chen and Dr. Edgardo S. Solano-Carrillo.
Next, I would like to express my thanks to Professors Igor Aleiner, Andrew Millis, Sebastian Will, Lam Hui, and David Reichman for serving on my defense committee. I would like to thank Professor Allan Blaer, Professor Igor Aleiner and Professor Boris Altshuler for their stimulating courses. I would also like to thank all the faculty and staff members in Physics department.
Finally, I would like to thank my parents Songming He and Yun Lei for their everlasting love and encouragement in my life.
(This page intentionally left blank.)
Chapter 1 Introduction
In this thesis, we study the strongly correlated materials with localized electron orbitals, which show up in the Hubbard model or Anderson impurity model as an intra-orbital Hubbard term between opposite spins on the same atomic site. A Hunds coupling is often introduced to describe anisotropies of the interactions between multiple orbitals. In this introductory chapter, we will discuss the main challenges encountered and techniques employed in the research works of this dissertation. These will include the density functional theory (DFT), Wannier orbitals, Hartree-Fock corrections of localized Wannier orbitals (often called DFT+), and dynamical mean-field theory (DMFT), which are state of the art for understanding strongly correlated materials.
1.1 Density functional theory
In most textbooks on solid state physics, the band theory is an important topic to cover because of its conceptual simplicity and computational efficiency. Yet its main limitation is assuming that electrons are noninteracting or that the electron-electron interactions can be treated on a mean-field level. The density-functional theory (DFT) [1], especially the Kohn-Sham DFT [2] is a mapping of the interacting many-electron system into an effective noninteracting system that reproduces the electron density exactly if the exact exchange-correlation functional is known. In practice, approximations of the functional are developed such as the local density approximation (LDA), generalized gradient approximation (GGA), and meta-GGA, etc., and generalizations of the theoretical framework are proposed to reproduce not only the density, but also the spin-density, one-body density matrix and even pairing amplitudes in superconducting systems with an effective noninteracting model. The density functional theory is thus a non-perturbative justification of how well one can approximate an interacting many-electron system by studying an auxiliary noninteracting system in a self-consistent loop. The challenge is, of course, that the more quantities to be reproduced by the noninteracting system, the more complicated the exchange-correlation functional becomes as it can depend on more quantities in a non-local way. Here we give a theoretical formulation of the density functional theory following the Levy-Lieb constrained search formalism proposed in [3, 4]. Consider an -electron system
[TABLE]
where we will use , , to refer to the kinetic energy, potential energy in external field due to the ions, and the electron-electron interaction energy. We notice that the part of the Hamiltonian is universal in all materials, which differ only in the potential energy due to the crystal field. We also notice that the crystal field only couples to the density via the one-body potential
[TABLE]
Therefore, the ground-state energy can be found by minimizing the average value of the total Hamiltonian with respect to the wave function , i.e.,
[TABLE]
where is a normalized -electron state satisfying . It gives the density . The idea of Levy-Lieb constrained search is to break the minimization into two steps:
[TABLE]
where the first step is to minimize over the wave functions that give the density , and then the second step is to minimize over to find the ground-state density. The big triumph of the density-functional theory is that the functional is universal, i.e., independent of the crystal field that is material specific. It only depends on , i.e., the electron kinetic energy and electron-electron interactions. Both the ground-state energy and the ground-state density can be found by minimizing the universal functional plus a linear coupling term of the crystal field with the density . The grand potential minimization formalism generalizes DFT to finite temperatures and fractional occupancies.
Eq. (1.4) is, of course, only a reformulation of the many-electron problem. The universal functional is as hard to find as solving a general many-electron problem with electron-electron interactions. It therefore requires approximations to be put to work in practice. We consider the difference of between an interacting system and a noninteracting system , subtract off the classical electrostatic potential energy, and define the remaining difference as the exchange-correlation energy , i.e.,
[TABLE]
hoping that the non-local structure of can be made short-range in to be more easily approximable. The local-density approximation (LDA), for example, assumes that the exchange-correlation functional takes the form
[TABLE]
where is the exchange-correlation energy per electron at that only depends on at the same position. The formula of is often determined by calculations or simulations of the electron gas with uniform density. The generalized gradient approximation (GGA) takes into account gradient effects of the density so one would work with and meta-GGA would allow the exchange-correlation energy density to depend on higher-order gradients of the density allowed by rotational symmetry, because is a universal functional that only depends on .
Once an LDA/GGA type exchange-correlation energy is given, the minimization of Eq. (1.4) can be done by solving an auxiliary noninteracting many-electron system, known as the Kohn-Sham system. Let us rewrite the energy functional as
[TABLE]
The energy functional contains a kinetic energy that involves minimizing a many-electron wave function subject to a given density plus other terms that are directly calculable from the density via some analytic or empirical formulas. By taking the first-order variation of the total energy , we obtain
[TABLE]
where the (noninteracting) kinetic energy functional is defined by
[TABLE]
and the Kohn-Sham effective potential is given by
[TABLE]
Minimizing the energy functional is locally equivalent to minimizing an auxiliary Kohn-Sham noninteracting system , where the Kohn-Sham potential depends on the density which must be determined self-consistently. If the exchange-correlation energy is known, the ground-state densities of and would be the same. The density functional theory then allows us to calculate the density of the interacting system by doing the self-consistent loops of the auxiliary noninteracting system .
1.2 Wannier orbitals
After doing a DFT calculation, we obtain the energy bands and Bloch waves over a k-point mesh of the first Brillouin zone (BZ). We sometimes want to build a minimum model that involves as few orbitals as possible that would reproduce the DFT band structure. The maximally localized Wannier functions (orbitals) are the basis for such a construction. The Bloch waves can be written as linear superpositions of the Wannier functions via
[TABLE]
where is cell-periodic and the matrix is unitary. The Wannier functions play the role of orthogonalized atomic orbitals. They are superposed into “molecular orbitals” by the unitary matrix and then form the Bloch wave via the sum over the lattice sites , which is similar to the method of linear combination of atomic orbitals (LCAO) in the solid-state physics textbooks. One may invert the Fourier transform and unitary matrix to obtain the Wannier functions
[TABLE]
where “BZ” is the first Brillouin zone in the reciprocal lattice basis (need not be orthogonal). In our unit system, the lattice vector is an integer vector in the Bravais lattice basis. The metric between is still identity. Since is cell-periodic, the calculated from Eq. (1.12) is indeed only a function of . This ensures that we choose a translationally invariant basis.
In the definition of the multi-band Wannier functions in Eq. (1.12), the unitary matrix is a gauge freedom, which in the single-band case becomes a single -dependent phase. The multi-band case allows more freedom, which we can use to narrow the spread of the Wannier functions. We first consider the Wannier center defined by
[TABLE]
and then the second-order moment defined by
[TABLE]
The spread of the th Wannier function can then be calculated and minimized by tuning the unitary matrices . The technical details of this part are handled by Wannier90 [5]. The spread is calculated in -space using
[TABLE]
where we have introduced the decoupled cell-periodic functions
[TABLE]
with having the meaning of the LCAO wave function formed by the th “atomic orbital” . A derivation is given in Appendix A to obtain Eqs. (1.15)–(1.16) using the mathematically elegant formula .
1.3 The DFT+ method
The density functional theory (DFT) is a highly successful method for calculating the electronic band structures of real materials. The auxiliary Kohn-Sham system in practice often not only gives a good description of the electron density, but also provides a reasonably good picture of the band structure. For some strongly correlated materials with localized orbitals (typically and orbitals), however, the DFT method suffers from the band gap problem (see Fig. 1.1),
meaning that the band gap given by the method is systematically too small. Some insulators are incorrectly calculated by DFT to be metals. The DFT+ method is a computationally cheap solution (compared with e.g. the GW method used in Fig. 1.1 or other more expensive Feynman-diagram-based method) to the band gap problem by introducing a Hartree-Fock correction term to the localized orbitals. The magnitude of the correction is controlled by the Hubbard (and sometimes also Hund’s coupling ) as adjustable parameters to fit with the experimental band structure. There are also methods that can give a reasonable estimation of the range of as a guideline, such as the constrained RPA [7, 8] and self-consistent linear response theory [9], etc. Here we will not go into the details of how to estimate the parameters and , but will mainly focus on the general idea of DFT+, with some detailed discussions on parameterizing the rotationally invariant interaction matrix elements in terms of and following [10].
1.3.1 Hartree-Fock approximation for localized orbitals
The band gap problem of Kohn-Sham DFT is due to the fact that the effects of the two-body interactions on localized orbitals are not well reproduced by the one-body Kohn-Sham potential. The locality of the orbital increases the correlation effects (repulsiveness) of electron occupancy in that once an electron occupies the orbital, it becomes very difficult to get occupied by another electron. This effect of reduced double occupancy can be reproduced by the Hartree-Fock energy of the two-body interactions on the localized orbitals. The energy functional of DFT+ is given by
[TABLE]
with the two-body Hamiltonian
[TABLE]
Here and sum over the spin directions and , and sums over localized orbitals on the same site. The Hartree-Fock energy of the on-site two-body interaction is therefore given by
[TABLE]
which depends on the on-site one-particle spin-density matrix . The last term is the double-counting energy, which subtracts off the interaction effects already considered in via the real-space local density .
Depending on the magnetic order (i.e. spin symmetry) of the system, we have types of DFT+ theories commonly used in energy band solvers for real materials calculations. For example, in the Vienna Ab-initio Simulation Package (VASP) (see http://cms.mpi.univie.ac. at/wiki/index.php/LDAUTYPE), the parameter settings are listed in the table below:
The gauge symmetry corresponds to the conservation of total number of electrons. The symmetry is the spin symmetry, which is fully preserved in paramagnetic (or diamagnetic) materials, partially spontaneously broken to in collinear spin systems (including ferromagnetic, antiferromagnetic, ferrimagnetic orders, etc), and fully broken in non-collinear spin systems (e.g. frustrated systems). The symmetry may be broken as well for attractive interactions, which would open the paring channels. Such DFT+ calculations with paring effects are not yet supported in VASP for materials calculations, but are extensively studied in model systems [11, 12].
1.3.2 Hubbard and Hund’s coupling
In materials calculations, the many interaction parameters in Eq. (1.19) are often determined by only two parameters: the Hubbard and Hund’s coupling , by considering the rotational symmetry of an isolated atom. Even though in a crystal, the symmetry is lowered due to other atoms, a rotationally invariant interaction can still be a good starting point. In an isolated atom, the on-site occupation matrix is diagonal in both spin and orbital angular momenta. Eq. (1.20) then reduces to
[TABLE]
where we have introduced the short-hand notations and , which are the direct and exchange interaction matrices. The Hartree-Fock energy between two electrons and is . The Hubbard and Hund’s coupling are defined by averaging the interaction over the orbitals, i.e.,
[TABLE]
On average, the repulsion between electrons of opposite spins is the Hubbard , while the average repulsion between electrons of the same spin is , weaker than the Hubbard by the Hund’s coupling due to the exchange effect.
1.3.3 A rotationally invariant Hamiltonian
If the interactions in Eq. (1.19) arise from a rotationally invariant two-body potential between equivalent electrons (with the same and ) on the same atomic site, the matrix elements
[TABLE]
can be parameterized by a few radial parameters due to the rotational symmetry. Let us expand the two-body potential in terms of Legendre polynomials as
[TABLE]
where denotes the th-degree Legendre polynomial, and write the orbital wave functions into the form
[TABLE]
Note that all orbitals in Eq. (1.23) have the same radial function and only differ by the angular part . If the above assumptions hold approximately true for the on-site Wannier orbitals, then we can parameterize the interactions in terms of the radial integral parameters
[TABLE]
via the universal Wigner -symbols
[TABLE]
We will give a detailed derivation in Appendix B. Only even-degree radial integrals enter into because of the parity selection rule. The conservation of angular momentum is also implied by the selection rule of the Wigner -symbols. We also show in Appendix B the sum rules of in terms of the Hubbard and Hund’s parameters in Eq. (1.22) given by
[TABLE]
So the Hubbard and Hund’s are also called the isotropic and anisotropic interactions, respectively. To parameterize a rotationally invariant interaction between electrons, we need only one parameter . To parameterize interactions between electrons, we need and . For electrons we need , , and , and so on. Empirically (typically a few eVs) is most significantly affected by screening and other renormalization effects, so it needs to be specified for every material. The anisotropies , , , of the interaction are specified proportional to one parameter via the sum rule, with the ratios of different ’s kept constant and specified empirically. A common choice for anisotropy is — eV for orbitals, with no strong dependence on materials [13].
1.3.4 The double-counting term
The double-counting correction is constructed by the same idea as Eq. (1.22). Assuming looking at only the local density cannot distinguish between different on-site orbitals, the interaction energy between electrons of opposite spins is and the interaction energy between electrons of the same spin is . Therefore, the double-counting energy to be subtracted off from is given by
[TABLE]
with is the number of electrons with spin and is the total number of electrons. This is the form of double-counting energy used in VASP called the fully localized limit (FLL). There are other forms of double-counting energy as well, such as the around mean-field (AMF) form. Some recent work to make the double-counting correction more rigorous is given in [14].
1.4 Dynamical mean-field theory
The density function theory (DFT) and DFT+ theory map an interacting electron system into an effective noninteracting system with a self-consistently determined band structure. The dynamical mean-field theory (DMFT) is a beyond-band-theory method formulated based on Green’s functions. The main idea is to choose one site of an interacting lattice model as an open system, and then based on the local Green’s function of the chosen site, we simplify the other environmental lattice sites into an equivalent noninteracting bath. The lattice model is then mapped into an Anderson impurity model with only the chosen site (the impurity) having on-site interactions (Hubbard or both and for multi-orbital impurities) and other orbitals noninteracting.
The idea can be formulated in the situation of a general open quantum system, with the total Hamiltonian of the system and the environment (bath) given by
[TABLE]
where and only act on the system and the environment respectively and acts on both. In the case of DMFT, includes the on-site orbital energy and on-site interactions of the impurity, includes the cavity lattice of all other sites, and refers to the hopping terms between the impurity and the bath. The Hilbert space is a direct product of that of the system and that of the environment . The nonequilibrium Green’s function of the system defined on the Keldysh contour is given by
[TABLE]
The Keldysh contour is a trajectory on the complex plane of time to go from on the real axis to and then back to and then down the imaginary axis to . For more details of the nonequilibrium Green’s functions, see e.g. [15]. The partition function . All operators with a time label for contour ordering are still in the Schrödinger picture. Hamiltonians are allowed to physically change with time. The operators and only act on the system . Subscripts are dropped to keep the notation simple. Let’s split the trace into partial traces over and . Since operators in Eq. (1.31) are ordered by , it is permissible to factorize the exponential and permute the operators to obtain
[TABLE]
Now we define an effective action
[TABLE]
with the partition function . The action contains at all times like a “functional” of operators. Protected by the contour-ordering , the operators at different times behave like the anticommuting Grassmann numbers (for a fermionic system ). We have omitted a lot of mathematical details to show that the exponential form exists and is well-defined over the ring of Grassmann numbers. The Green’s function of the system is then written as
[TABLE]
with defined as the partition function of the open system . All of the environmental degrees of freedom have been traced out by to give rise to an effective action of the system’s degrees of freedom.
We have formulated very conceptually the effective action theory for open quantum systems. The action contains richer physics than Hamiltonians. For example, in systems with electron-phonon coupling, in terms of the electronic degrees of freedom gives rise to a time-delayed attractive two-body (four-operator) interaction mediated by the noninteracting phonons. Similarly, the effective action produced by a noninteracting fermionic bath in the situation of DMFT is a time-delayed one-body (two-operator) hopping term
[TABLE]
governed by a hybridization function . We will give a detailed derivation of (1.34) and a specific expression for the hybridization function in terms of the bath spectrum and impurity-bath coupling strengths in Appendix C. Review papers of DMFT [16, 17, 18] show that the of an interacting cavity lattice also reduces to the form in Eq. (1.34) in the infinite dimension (or infinite coordination number) limit, justifying the approximation of DMFT in high spatial dimensions.
1.5 Summary and conclusion
We have given brief introductions to state-of-the-art techniques used for calculating the electronic states of strongly correlated systems with significant on-site interactions for localized orbitals. The density functional theory in its Kohn-Sham self-consistent field formulation has proved to be highly successful for many types of real materials. For materials with localized (typically or ) orbitals, the Coulomb repulsion of electrons on these orbitals are significant and cannot be well approximated by a local Kohn-Sham field that couples to the local density. The computationally cheap solution is to use DFT+, which includes the Hartree-Fock energy of the localized orbitals to construct a nonlocal potential that couples to the orbital occupancy, or the on-site occupation matrix. Depending on the magnetic order of the system, different types of corrections can be included. A more accurate but computationally expensive solution is to use the dynamical mean-field theory (DMFT), which keeps the full interactions on the localized orbitals treated as impurities and only attempt to map the delocalized orbitals into an effective noninteracting bath. Other interesting topics such as the self-consistency conditions of nonequilibrium DMFT, generalizations of DMFT to clusters of lattice sites, and the DFT+DMFT method for real materials calculations are not discussed in this thesis.
In the following chapters, we use the DFT+ method to study strongly correlated materials in Chaps. and , and do a focused study towards building a nonequilibrium DMFT impurity solver in Chaps. and . We study the equilibrium phase transitions in LuNiO3, out-of-equilibrium phase transitions of VO2 in a pump-probe experiment, and use the density matrix renormalization group (DMRG) method as an impurity solver for real-time DMFT with quench and periodically driven Hamiltonians. There are good review papers for the DMRG method [19, 20] and its applications to real-time evolutions [21, 22] of nonequilibrium systems. We will defer our discussion of the implementation details of the DMRG method to Chaps. and .
Chapter 2 Strain control of electronic phase in rare-earth nickelates
In this work, we study the structural phase transitions and metal-insulator transitions of LuNiO3 as an example of the rare-earth nickelates NiO3 induced by a compressive or tensile substrate strain using the DFT+ method. The rare-earth nickelates crystallize in variants of the O3 perovskite structure, with the ion on the site and Ni ion on the site. The basic structural motif is a corner-shared O6 octahedron, which can have bond-length distortions and tilts that give rise to competing electronic phases with different charge and orbital orders. We use group theory to construct a Landau energy function in terms of the distortion modes based on the calculations of DFT+, to study the competition between different electronic phases on a phenomenological level. The calculation shows that under % compressive or tensile strain, the insulating charge-ordered phase destabilizes to a metallic Jahn-Teller distorted phase. The long Ni-O bonds point out of plane under compressive strain and form an in-plane checker-board pattern under tensile strain. The two Jahn-Teller distorted phases are smoothly connected due to the octahedral tilts, while the jump from the charge-ordered phase to the Jahn-Teller distorted phase is a discontinuous first-order transition at both critcal strains. It is interesting that the magnitude of the critical strains are of the order of strains accessible by epitaxial growth on substrates. Our work in this part was published in [23].
2.1 Crystal structure of rare-earth nickelates
The rare-earth nickelates have been of substantial research interest for many years. Their chemical formula is NiO3, with standing for a rare-earth element, including Sc, Y, and the lanthanide series from La to Lu. The crystal structure of the material for Lu in its ground state is shown in Fig. 2.1. The structure is characterized by corner-shared and tilted NiO6 octahedra with Ni-O bond lengths alternating in a checkerboard pattern. This bond disproportionation is sometimes referred to as “charge ordering” based on the idea that the ionic charge of the Ni ion with longer Ni-O bond lengths should be larger than that of the Ni ions with shorter Ni-O bonds. Although the actual charge difference between the sites is very small [25, 26], for simplicity we will refer to the disproportionated state as “charge ordered”. The unit cell has four inequivalent NiO6 octahedra. In the absence of charge ordering, the octahedra differ only by rotations; the charge ordering creates two classes of octahedra with different mean Ni-O bond lengths. Figure 2.1 also shows the lattice constants. The Ni-Ni distance in the basal () plane is Å, and there is a slight rhombic distortion, so the Ni-Ni bond angles are and .
We use the DFT+ calculation as our numerical experimental apparatus to simulate the effects of placing LuNiO3 on a substrate, which will typically have a square symmetry. We therefore neglect the rhombic distortion and consider square structures with and Ni-Ni bond angles in the plane. We define the -plane lattice constant . The equilibrium lattice constant is Å at which the energy is minimum. We will be interested in the consequences of a uniform compression or expansion of the lattice in the plane with the direction free to adjust.
2.2 DFT+ calculation
Our calculations use the Vienna Ab initio Simulation Package (VASP) [27, 28]. The DFT+ algorithm we use in VASP is the rotationally invariant local spin-density approximation (LSDA) that follows [10]. The Hubbard of the Ni orbitals in LuNiO3 can be obtained with various methods, e.g., constrained local-density approximation [29, 30], self-consistent linear response [9], constrained random-phase approximation [7, 8], etc. They all give values of within eV. The Hund’s coupling is estimated to be – eV. We finally chose eV and eV, as they gave a structure in Fig. 2.1 that was closest to the experimental results. Slight changes of and within their errors were tried, and no qualitative difference was found.
We did a spin-polarized calculation using the Projector augmented-wave Perdew, Burke, and Ernzerhof (PAW-PBE) pseudopotential provided by VASP. The k-point mesh we used was , and the energy cutoff of the plane-wave basis was set to eV. We found two magnetic states in the charge-ordered structure: ferromagnetic (FM) and A-type antiferromagnetic (A-AFM) states with magnitudes of magnetic moments essentially on Ni orbitals modulated by octahedral sizes. The FM state is lower in energy than the A-AFM state at all values of lattice constant a in our DFT+ calculation. All results are obtained in the FM state.
The computational unit cell was chosen to contain four LuNiO3 formula units. Defining the basal plane as the one in which strain is applied, we take two formula units in the basal plane and two displaced vertically. To mimic the effects of a substrate, the in-plane lattice constants are fixed to preset and equal values (so any in-plane rhombic distortion is neglected). and all of the intra-unit-cell degrees of freedom are allowed to relax. We slightly modified the conjugate gradient code in VASP to do this. The minimum energy of the substrate-constrained system is obtained at Å. The structure obtained is almost identical to the free structure in Fig. 2.1, except that and are made equal (the small rhombic distortion is suppressed). We then adjust the substrate lattice constant , our control parameter, away from and see how the structure changes.
2.3 Landau energy function based on group theory
The main technical part of this work is using group theory to analyze the distortion modes observed in the DFT+ crystal structures. We begin with the Landau energy function of a single NiO6 octahedron to demonstrate how group theory works in our situation. Then we consider an array of NiO6 octahedra with no tilts (rotations) and study the bond-length distortion modes. Finally, we include the effects of octahedral tilts perturbatively and see what symmetries they break.
2.3.1 An isolated NiO6 octahedron
To define notation we begin by considering one isolated NiO6 octahedron. The unstrained structure is perfectly cubic (point symmetry ) with six mutually perpendicular Ni-O bonds, which we take to lie in the , , and directions. All six bonds have the same length, Å. The distortions of interest here preserve the inversion symmetry about the Ni ion and the orthogonality of the Ni-O bonds, so that minimally a symmetry is preserved. The distortions may be expressed in terms of three modes, defined in terms of the changes , , in the bond lengths as
[TABLE]
Here is the volume expansion mode, is the (volume-preserving) -plane square-to-rhombic distortion, and is the (volume-preserving) cubic-to-tetragonal Jahn-Teller distortion in the direction. In general, the energy function of an isolated NiO6 octahedron needs to be invariant under , which is isomorphic to the permutation group of the three directions . It should therefore be a linear combination of the permutation-symmetric polynomials
[TABLE]
where we have Taylor expanded to 3rd order. The linear terms vanish because we are expanding around the equilibrium length . In terms of the modes , the quadratic terms decouple and we obtain
[TABLE]
The mode is invariant under and can be arbitrarily coupled to other modes. We absorb it into the coefficients of the Taylor expansion of and , which together form a two-dimensional irreducible representation of . We highlight the cubic coupling in the last term with coefficient . In the lattice system, this part will give rise to an important coupling between the distortion and the staggered Jahn-Teller order , which we will define later.
2.3.2 A corner-shared NiO6 array
We next consider an infinite three-dimensional array of NiO6 octahedra, still with the symmetry in the unstrained structure at each Ni site. We must now attach a momentum label to each mode. In addition, because the octahedra are corner shared, there are constraints on the allowed momenta for each distortion. The momenta of interest are , , . Note that these momenta are defined in the unit cell of the ideal cubic structure with one octahedron per unit cell. Of primary interest in interpreting the numerical results are the two-sublattice charger-order and the in-plane staggered Jahn-Teller modes, written as and , respectively. In addition, it will be useful to consider , , and , which are the volume change, uniform Jahn-Teller, and two-sublattice Jahn-Teller modes, respectively, which describe the response to a uniform strain and its coupling to a two-sublattice charge order. Modes , , and are visualized in Fig. 2.3. The DFT+ calculation shows that there are no other modes to consider than these five.
The energy function of the five modes is, in general, very complicated. A group theoretical analysis is given in Appendix D. The variables and are controlled by the lattice constant , which induces a distortion and, via Poisson-ratio effects, a nonzero volume change of opposite sign to . Both and are coupled to the order parameters , , and , and these couplings will drive the phase transitions of interest. Based on the results of Appendix D, if we express and as smooth functions of , then the Landau energy function in terms of the non-uniform distortions , , and as order parameters is given by
[TABLE]
The smoothness assumptions and are justified by the results of DFT+ calculations plotted in Fig. 2.4. The jumps in and at the critical lattice constants are much smaller than the jumps of the non-uniform modes , , and .
A further simplification can be made by noticing in Fig. 2.4 that the order parameters and , both at the point , are always simultaneously nonzero, as in the charge-ordered structure, or simultaneously zero when the order vanishes under a large enough compressive or tensile strain. The fact that and always coexist suggests that we may combine them into one order parameter. This can be done by treating the ratio as a smooth function of . The Landau function is now further reduced to one with only two order parameters:
[TABLE]
where the coefficients
[TABLE]
are independent and smooth functions of . Equation (2.5) gives the general form of the symmetry-based Landau energy function of NiO3 without considering perovskite octahedral rotations and nonorthogonal Ni-O bond angles.
2.3.3 Including octahedral rotations
We have been ignoring octahedral tilts in the previous sections. The actual structure of the material involves a GdFeO3-type rotational distortion that may be symbolically written as . The notation means that starting from the ideal cubic perovskite structure, there is a rotation by angle about the axis and by angle about the and axes. The superscript plus sign means the rotations in neighboring octahedra along the rotational axis of (the axis) are in the same direction, while the minus sign means the rotations in neighboring octahedra along the rotational axis of ( or axis) are in opposite directions. The displacement field of the rotational pattern in the plane is shown in Fig. 2.5. Since angles and are small ( in LuNiO3), we may neglect the non-Abelian aspect of rotations and treat them as an additive displacement field.
The important feature of the octahedral rotations is a breaking of the symmetry while preserving the symmetry of Eq. (2.5). The symmetry-allowed energy function of variables , , , and is given by
[TABLE]
The omitted terms include other quartic terms that are products of the quadratic ones and higher-order terms. The leading-order term that breaks the symmetry is , which is linear in . The coefficient is of order . The derivation is using group theory similar to Appendix D. The symmetry group for the energy function at fixed lattice constant is D4h. Since all axial vectors , , of the rotations and bond-length modes , are invariant under spatial inversion , only , which contains symmetry operations, is effective in actually transforming the modes. In addition to , the translations can generate possible ways of sign change according to the points of the modes, among which and are at and , , and are at . Therefore, we have totally symmetries to satisfy. Following again the rearrangement-theorem-based algorithm in Appendix D, we get the general form of the symmetry-allowed Taylor expansion of the energy function in Eq. (2.7). The symmetry is strictly preserved order by order. Switching the sizes of the larger and smaller NiO6 octahedra of the charge-ordered structure is still a symmetry of the system even in the presence of the GdFeO3-type octahedral tilts.
We therefore add the leading-order symmetry-breaking term to the original Landau function in Eq. (2.5) as a perturbation to get the symmetry right. The new Landau function is given by
[TABLE]
The added term should be small because for small rotations and . It should therefore be ineffective unless the even-power coefficients make the state unstable or nearly unstable. Aside from octahedral rotations, nonorthogonal Ni-O bond angles can also break the symmetry if the Ni-O bond that is approximately along the direction forms different angles with the and bonds. The leading-order symmetry-breaking term should also be small and linear in and can therefore be addressed on the same footing as octahedral tilts.
2.3.4 Minimum model construction
In Eq. (2.8), the effects of octahedral tilts are considered perturbatively with only the leading order term included. To understand the phase transitions in Fig. 2.4, the even-power terms can be truncated to some highest order as well. In this section, we construct a Landau energy function with the minimum number of terms in the expansion of Eq. (2.8) and the simplest strain dependence of the expansion coefficients. Based on the observations in Fig. 2.4, the model needs to have the following features:
- In the charge-ordered phase (solid lines) with essentially , the energy as a function of has a first-order transition at both critical strains. Since can only contain even powers of due to symmetry, we have
[TABLE]
with and near the transition. Thus, has three local minima at and . A model with and that is bounded below and truncated at th degree can only exhibit second-order transitions.
- In the Jahn-Teller phase (dashed line) with charge order suppressed, the energy as a function of has an avoided second-order transition structure. We have
[TABLE]
The term comes from the symmetry-breaking octahedral tilts, which is a small perturbation and gets strongly suppressed if under compressive strain, but becomes important and allows to smoothly grow from small to large values when changes sign under tensile strain. The term has an effect similar to that of an external magnetic field on a system near a ferromagnetic transition.
- Since the charge order strongly suppresses the Jahn-Teller mode (as can be seen by the jump up of at the critical tensile strain), there is a big competition term between and that should be allowed by the cubic symmetry. The simplest form is a biquadratic term, so the full energy function is constructed as
[TABLE]
The last term with stabilizes the state when the charge order is present and vice versa, which explains the absence of coexistence of and in Fig. 2.4.
2.4 Analysis of numerical results
Based on the Landau energy model constructed in §2.3, we can now interpret and understand the phase transitions in Fig. 2.4. We have also done some corroborative calculations using DFT+ for the statements in the previous section, which are shown in this section alongside our interpretations of Fig. 2.4.
2.4.1 Structural transitions and energy difference
The most significant findings of Fig. 2.4 are the discontinuous jumps of the order parameters , , and at the critical compressive and tensile strains. The transitions being first-order are corroborated by the energy difference of the stable charge-ordered (CO) and metastable Jahn-Teller distorted (JT) phases plotted in Fig. 2.6. At zero strain Å, the charge-ordered structure is lower in energy than the Jahn-Teller structure by meV per unit cell (with 4 Ni ions). Under either a compressive strain () or a tensile strain (), the Jahn-Teller structure is favored, and is reduced. At both transition points, the curve overshoots a little bit to below zero and ends where the charge-ordered structure becomes locally unstable and relaxes to the Jahn-Teller structure. Both the overshoot and the linear relation near the transitions confirm that the transitions are first order.
The transition at compressive strain does not involve the mode. The long bonds of the Jahn-Teller phase are out of plane in the direction (as indicated by the uniform mode in Fig. 2.4). We did DFT+ calculations of a series of linearly interpolated structures between the charge-ordered and Jahn-Teller distorted phases at various lattice constants close to the critical compressive strain at around Å, to reproduce the Landau energy function that gives the first-order phase transition. Results are plotted in Fig. 2.7. We see that the energy function has two locally stable minima crossing in energy as the lattice constant is changed. When the lattice constant is way above the transition point, the Jahn-Teller phase with the charge-ordering mode suppressed is locally unstable and relaxes to the charge-ordered ground state.
The transition at the critical tensile strain ( Å) involves the dying off of the charge-ordering mode and the jump up of the in-plane staggered Jahn-Teller mode . The first-order transition of is the same story as the compressive strain case. The sudden jump up of the mode is the result of the biquadratic coupling , which reduces the quadratic coefficient of from to and triggers the instability of the state. To remove the suppressive effect of to , we did DFT+ calculations with the symmetry enforced (see dashed lines in Fig. 2.4) to study the evolution of the Jahn-Teller phase as lattice constant changes in the next section.
2.4.2 Evolution of the Jahn-Teller structure
The Jahn-Teller structure with charge-ordering mode enforced numerically is plotted as dashed lines in Fig. 2.4. Here we focus on the evolution of the mode as lattice constant changes. The nonzero is a consequence of the GdFeO3 octahedral tilts, which, as previously discussed, couple linearly to the staggered component of the Jahn-Teller distortions. We do some parameter fitting in this section to understand the avoided second-order phase transition of going from very small values to suddenly very large values as lattice constant increases under a tensile strain from the substrate.
A minimum model to understand this evolution of the Jahn-Teller structure from Eq. (2.10) with strain dependence is given by
[TABLE]
where is assumed to be constant for simplicity, and is the deviation of the lattice constant from its equilibrium value Å. Equation (2.12) is formally similar to the equation describing a ferromagnet in a magnetic field. The coefficients and are like an external magnetic field in the ferromagnetic case and arise from the breaking of symmetry due to the GdFeO3 rotations. The need to allow for a strain dependence of the coefficients is shown by the zero crossing of at Å. The dependence of on strain reflects the tendency of tensile strain to favor the staggered Jahn-Teller order . Minimizing Eq. (2.12) leads to
[TABLE]
Because we have kept dependence only to linear order, it is easy to express in terms of the equilibrium Jahn-Teller amplitude via
[TABLE]
We have fit Eq. (2.14) to the data points shown in Fig. 2.8, and from the fit parameters we extracted the critical lattice constant Å at which the hypothetical cubic structure would be unstable to staggered Jahn-Teller order in the absence of charge order or GdFeO3 rotations. We observe that while the uncertainties involved in fitting a four-parameter function to the data mean that individual coefficients cannot be determined with high accuracy, the estimated is robust. It is interesting that this value is not very much larger than the value of Å at which the charge order vanishes.
2.4.3 The competition between and
Comparison of the solid and dashed lines in Fig. 2.4 shows that the staggered charge order strongly suppresses the staggered Jahn-Teller order . In the notation of Eq. (2.11), the biquadratic term is large and repulsive. In terms of the analysis of Eq. (2.12), in the presence of the charge-ordering mode , the quadratic coefficient becomes
[TABLE]
and is so much more positive that until the charge order collapses at a first-order transition, the staggered Jahn-Teller order cannot develop. There is therefore a strong competition between the two staggered orders, and .
2.5 Insulator-to-metal transitions
The structural phase transitions of LuNiO3 from its charge-ordered phase with mode into its Jahn-Teller distorted phases with the out-of-plane and in-plane staggered modes are accompanied by the collapse of band gap, i.e., insulator-to-metal transitions. This makes the structural phase transitions very interesting to study. In Fig. 2.9, we plot the energy gap as a function of the lattice constant to show how the energy gap collapses. We have slightly varied the interaction parameters and within their reasonable ranges to test the numerical sensitivity of the band gap plot. Within DFT+, the charge-ordered phase is insulating with an energy gap of about eV while the Jahn-Teller phases under both compressive and tensile strains are found to be metallic with no gap at the Fermi level. The effects of electron-electron correlations modeled as the terms couple through lattice relaxation to the distortion modes , and , which in turn determine the electron orbital energies and whether or not a band gap opens.
2.6 Summary and conclusion
We have used DFT+ and Landau theory methods to consider the effects of strain (induced by growth on a substrate with different lattice constants) on the charge-ordered state of LuNiO3. We find that the charge-ordered state plays a primary role in controlling the physics. It is the leading instability under ambient conditions, and its presence suppresses any other instabilities. However, with sufficient applied strain (within the DFT+ method, of the order of ) the system undergoes a first-order transition a non-charge-ordered state. Interestingly, for tensile strain, the non-charge-ordered state is characterized by a staggered Jahn-Teller order. In the actual crystals, the symmetry breaking induced by the GdFeO3 rotational distortion means that the staggered Jahn-Teller order does not break any additional symmetry of the system.
The actual magnitude of the strain needed to destabilize the charge order and allow other states is an important open question. While we imagine the strain as being produced by epitaxial growth on a substrate, we have not included any quantum confinement effects in our model. Also, the DFT+ method we have used is known to overestimate the tendency to charge order [31]. The charge-order phase boundary also depends on how the double-counting correction is implemented. More refined calculations, perhaps based on DFT+DMFT methods, should be employed to obtain better estimates for the strain needed to destabilize the charge order. But it is interesting that the magnitude of strain we have found is of the order of strains accessible by epitaxial growth on substrates.
Chapter 3 Photoinduced phase transitions in narrow-gap Mott insulators
In this chapter, we study the nonequilibrium dynamics of photoexcited electrons in the narrow-gap Mott insulator VO2. The material is famous for its metal-insulator transition at C, above which temperature it is metallic in a rutile () crystal structure and below which temperature it is insulating in a monoclinic () crystal structure with a doubled unit cell [32]. In a recent pump-probe experiment [33], a metastable metal phase of VO2 is found to exist for ps within an intermediate fluence range of the pump laser, as measured by ultrafast electron diffraction (UED) to have no crystal structural transition, and measured by infrared (IR) absorption to have a complete insulator-to-metal transition, while the temperature was kept at C below the transition temperature of the equilibrium phases. As a follow up work, it is found in [34] that the metastable metal phase of VO2 could be stabilized by applying an epitaxial strain.
In our work, we build a soft-band model using DFT++ to understand the metastable metal phase of VO2 in [33]. Here the “softness” of the band structure means that the self-consistent field depends on the density and orbital occupation matrix of the electronic state, which is the crucial driving force of the photoinduced phase transition. Both the on-site interactions as in DFT+ reviewed in Chap. and intersite interactions between the V-V dimers of the crystal structure of VO2 are included on the Hartree-Fock level. The initial stages of relaxation are treated using the quantum Boltzmann equation (QBE), which reveals a rapid (fs time scale) relaxation to a pseudothermal state characterized by a few parameters that vary slowly in time ( fs). We have established a momentum-averaged QBE that significantly reduces the number of dynamical variables but still captures the time scales of the main physical processes. The long-time limit is then studied by the DFT++ phase diagram, which reveals the possibility of nonequilibrium excitation to a new metastable metal phase that is qualitatively consistent with Morrison’s experiment. The general physical picture of photoexcitation driving a correlated electronic system to a new state that is not accessible in equilibrium may be applicable in similar materials. This part of our work was published in [35].
3.1 The DFT++ method for VO2
Following [36], we construct an electronic band structure for VO2 using the density functional theory (DFT)++ method, in which the basic density functional theory is supplemented by a Hartree-Fock treatment of the on-site (“+”) and intersite (“+”) - interactions. Belozerov et al. have constructed a DFT+DMFT+ theory with very similar physics [37]. The effects of the + term are a reasonable representation of the intersite self-energy terms found in the cluster DMFT calculations of [38]. Note that in the correct orbital basis, these intersite self-energy terms have only a weak frequency dependence [39]. Let us write the Kohn-Sham Hamiltonian of the electrons in their ground state as [10, 36]
[TABLE]
where comes from a density functional band calculation, is the Hartree-Fock approximation to the electron-electron interactions involving the vanadium orbitals, and is the double-counting correction. In the phase of VO2, the unit cell contains four vanadium ions, which form two dimerized pairs. We only consider interactions within one unit cell. These may be generally written as
[TABLE]
Here labels the unit cells, run over the correlated orbitals in a unit cell, and label the spins. We consider two contributions to : the on-site intra- interactions, which we take to be the rotationally invariant form [10] including both and orbitals parameterized by the Hubbard and Hund’s coupling , and intersite interactions between the two vanadium ions in each dimer. The Hartree-Fock approximation of the electron-electron interactions takes the form
[TABLE]
where in a non-spin-polarized system (like VO2)
[TABLE]
and the occupation matrix
[TABLE]
are independent of both spin and unit cell coordinate . In Eq. (3.4), has both the on-site and intersite intradimer terms. The on-site terms are the usual ones treated in standard DFT+ calculations discussed in §1.3. The intersite terms are parameterized by a single parameter and their contributions in take the form
[TABLE]
which contains only the Fock terms of the density-density interaction . The intersite Hartree terms are assumed to be already included in and are not included again in [36]. The Fock terms are orbitally diagonal, meaning that the and sum over only orbitals of the same type (e.g., , , etc.) in the two vanadium ions in a dimer. The intersite matrix element (hybridization) between different types of orbitals is typically small. In the ground-state insulating phase, only the hybridization of orbitals makes an appreciable contribution to , but in the nonequilibrium metastable states, hybridizations of other orbitals may be also important, so we will keep the terms of all five orbitals in the Hamiltonian .
We first performed a non-spin-polarized DFT+ calculation using the Vienna Ab initio Simulation Package (VASP) with the atomic positions fixed in the experimental structure [41]. We used a -point mesh of , an energy cutoff of 600 eV, and the projector-augmented wave Perdew-Burke-Ernzerhof (PAW-PBE) pseudopotential [42] in the VASP library. The on-site interactions are parameterized by eV and eV [43]. The in Eq. (3.1) is then defined as the projection of the DFT+ Hamiltonian onto a basis obtained from a Wannier fit to the 24 O- and 20 V- orbitals using Wannier90 [44] but with the on-site contributions to and the double-counting terms removed. These on-site contributions plus the intersite Fock terms in Eq. (3.6) make up the remaining terms in Eq. (3.1).
The DFT++ band structure for VO2 is plotted in Fig. 3.2 for eV. The results are in good agreement with preexisting results obtained using the GW method [45] and cluster dynamical mean-field theory (CDMFT) [38]. The validity of modeling VO2 in a renormalized band picture is corroborated in [39]. The optical gap at the Fermi level is 0.62 eV in good agreement with experiment [46]. The indirect gap between the highest occupied and lowest unoccupied Bloch states (the HOMO-LUMO gap) is 0.45 eV. The lower gap separating the V- and O- dominant bands below the Fermi level is 0.55 eV. The bonding-antibonding splitting of the orbitals arising from the dimerization of the crystal structure and enhanced by the intersite Fock interaction is eV in agreement with optical conductivity data [47]. The optical gap and the bonding-antibonding splitting are our main experimental evidences for determining and . But since the latter measurement is less accurate, the range of parameters – eV and correspondingly – eV provide equally reasonable descriptions of the material.
3.2 Initial absorption of laser energy
Next we estimate the energy range and number of electrons photoexcited in Morrison’s pump-probe experiment [33]. The wavelength of the pump laser is nm ( eV). Solving the optics problem for the experimental geometry specified in the experiment reveals that the laser fluence of – mJ/cm2 that yielded an metal initially generates – electron-hole pairs per unit cell (4 VO2), corresponding to an energy increase per unit cell of – eV. The details of the calculation are given below.
The complex dielectric constant of VO2 to the nm laser is given in [48], which yields a complex index of refraction . The index of refraction of the Si3N4 substrate is to nm. The thicknesses of the VO2 sample and the Si3N4 substrate nm and nm are given in the Supplemental Material of [33]. These data allow us to reconstruct the experimental setup in Fig. 3.3. Since the duration of the laser pulses used in the experiment is fs, which is equivalent to over oscillation periods of the nm laser, the absorption of energy from the laser pulse can be obtained to adequate approximation by solving steady-state wave equations. Nonlinear optical effects are neglected a posteriori because the density of excited particle-hole pairs is small. We may then use the formulas given in [49], assuming normal incidence ( according to he Supplemental Material of the experiment). The formula can be derived using the matrix equation
[TABLE]
which is obtained from the boundary conditions of the continuity of E and B fields and the propagation of waves in each medium. Here is the wave number in vacuum, and and are the reflectivity and transmissivity of the complex amplitudes of the E fields. The numerical result of solving Eq. (3.7) is that of the incident fluence gets reflected, gets transmitted, and the remaining gets absorbed. Then we use the density g/cm3 of VO2 in phase to calculate the unit cell volume to obtain and per unit cell.
3.3 Fixed-band QBE dynamics
In this section, we use the quantum Boltzmann equation (QBE) to study the relaxation of electrons after the laser pulse energy is initially absorbed by them to create some electron-hole pairs in the band structure of Fig. 3.2. For simplicity, we assume that the band structure is fixed, i.e., we forget about the soft-band effect due to the dependence of the self-consistent field on electron density and orbital occupancies, to estimate the relaxation time scale. We make another simplification by constructing a momentum-averaged QBE to significantly reduce the number of dynamical variables. We find that the energy gap is the main bottleneck of the relaxation dynamics, and electrons would equilibrate to a thermal state over a time scale fs. Since the soft-band effect would close or narrow the energy gap as electron-hole pairs are created (to be discussed in the next section), we expect the real electrons to reach the thermal state even faster. This allows us to understand the metastable metal phase in a hot electron picture in the next section (§3.4).
3.3.1 Formalism of the -averaged QBE
We begin with the formalism of the quantum Boltzmann equation (QBE) [50] to study the relaxation of the photoexcited electrons before energy dissipates into other slower degrees of freedom such as phonons. The quantum Boltzmann equation is a dynamical equation for the occupancies of the Bloch states in an electronic band structure, e.g.,
[TABLE]
Here is the DFT++ Hamiltonian in Eq. (3.1), sums over k-points in the first Brillouin zone, is the band index, and labels the spin. The Kohn-Sham eigenvalues do not carry a spin index in a non-spin-polarized system like VO2. The quantum Boltzmann equation treats electron-electron interactions as in Eq. (3.2) via Fermi’s golden rule, which gives the transition rates due to the two-body Hamiltonian between different Slater-determinant eigenstates of one-body Hamiltonian . While this perturbative, golden-rule-based method fails to capture important aspects of correlated electrons, the orders of magnitude of the relaxation time scales and the qualitative features of the resulting orbital distributions should be reasonably reproduced by this simplified dynamical model. In a non-spin-polarized system, the quantum Boltzmann equation (QBE) is given by
[TABLE]
where is the total number of k-points, and the matrix element is a short-hand symbol for summed over the and cases. The occupancies are single-spin quantities. The k-variables sum over only the first Brillouin zone and the Kronecker is to be interpreted as implying equivalence up to a reciprocal lattice vector to correctly impose the conservation of crystal momentum. A direct simulation of Eq. (3.9) in a general band structure is numerically difficult. The main problem comes from the energy function, which requires . To ensure the conservation of energy in each scattering process to the needed accuracy, one has to choose a very dense k-point mesh, which then leads to too many degrees of freedom to handle in a practical simulation. In order to obtain a computationally tractable model that still captures the important physics, we construct a momentum-averaged quantum Boltzmann equation, whose key variables are the energy distributions of electrons in different bands without any k-point information. Let us begin the derivation by averaging the matrix elements over the four k-variables to introduce
[TABLE]
which are the k-averaged matrix elements that only depend on the band indices . The motivation for the k-averaging comes from the local nature of the interaction defined in Eq. (3.2). The k-dependence of comes purely from the Bloch wave functions and tends to be complicated, and effectively random in real materials, so averaging over the momentum variables is reasonable. Next, we assume that the occupation numbers of the Bloch states
[TABLE]
are only functions of band index and energy . Then defining the single-spin density of states of band
[TABLE]
and the densities of occupied and empty states
[TABLE]
we derive a k-averaged QBE
[TABLE]
The band indices are kept in full. The ab-initio rate constants are obtained from Eq. (3.9) using Monte Carlo methods on a Wannier interpolated k-point mesh of . We will give a detailed derivation of Eq. (3.15) in Appendix E.
3.3.2 Simulation and analysis of numerical results
To run the simulation using Eq. (3.15), we need to specify the initial conditions, i.e., how the initially generated electron-hole pairs as calculated in §3.2 are distributed over the energies. We assume for simplicity that the laser absorption is proportional to the product of densities of states at energy separation eV. Then at , immediately after the laser pulse, we have the distributions of holes and electrons given by
[TABLE]
where satisfies and eV, the HOMO-LUMO gap. Here the subscript “tot” means to sum over all bands . The total number of electron-hole pairs is determined by the experimental laser fluence, as discussed in §3.2. Then we assume that the initially excited electrons and holes are randomly distributed over band states, i.e., for all energy , the density of occupied states in band ,
[TABLE]
is directly proportional to the density of states in band . We then evolve the distribution according to Eq. (3.15). We find that the equilibration process comes in basically two steps: the fast prethermalization (Fig. 3.4) that establishes a pseudothermal distribution characterized by a common temperature but different chemical potentials and for the electrons and holes, and then the slow evolution of thermal parameters (Fig. 3.5) to the final thermal state.
Figure 3.4 shows the initial stages of relaxation for laser fluence = mJ/cm2, comparing the calculated distribution to the distribution expected if the electrons and holes have thermalized. In the first fs after the laser pulse, the distribution of photoexcited electrons develops a tail to both high and low energies. Then in the next 1–2 fs, the electron and hole distributions thermalize. At the same time, the number of electrons and holes begins to increase due to the inverse Auger process, in which a high-energy electron scatters to a low-energy state while creating an electron-hole pair, thereby increasing the electron and hole densities and shifting the main weight in the conduction band to lower energies (a similar effect was noted in the Hubbard model by [51]). However, as the electrons thermalize, the inverse Auger scattering rate decreases rapidly since only electrons far out in the tail of the pseudothermal distribution have enough energy to down-scatter to create an electron-hole pair while still remaining in the conduction band. By fs, the electron and hole distributions are fully thermalized and the subsequent evolution can be described by the evolution of thermal parameters. For higher laser fluence = 9 mJ/cm2 (not shown) the time evolution of electron and hole distributions is qualitatively the same as shown in Fig. 3.4 and takes roughly the same time, but produces more electron-hole pairs (Fig. 3.5).
The evolution of thermal parameters, i.e., the temperature , the chemical potential of the electrons, and of the holes, is much slower as noted above. Figure 3.5 shows the results for both the low fluence = 3.7 mJ/cm2 and the high fluence = 9 mJ/cm2. The equilibration time constant approximately scales as the inverse of the square of the number of electron-hole pairs at equilibrium, which is a signature of the three-particle Auger and inverse Auger scattering processes.
Much of what happens in the simulation are explained by the rate constants . The largest rate constants are those of the hole-hole, electron-hole, and electron-electron scattering processes that do not change . The pair creation and recombination processes that change are comparatively slow. This separation of time scales has two origins: (a) the gap, which means that the processes must involve electrons in the tail of the distribution, and (b) the different orbital characters of the top of the valence band () and the bottom of the conduction band ( and ) in Fig. 3.2, which means that changes in must come from orbital-changing interactions, i.e., the pair hopping and exchange terms , which are much smaller than the orbitally diagonal interactions .
Even though the density relaxation of is much slower than prethermalization, due to the combination of small matrix element and kinetic bottleneck, our QBE-based simulation still finds that electrons in VO2 will equilibrate in hundreds of femtoseconds. The higher the laser fluence, the more electron-hole pairs are generated, and the faster the electrons equilibrate, as is shown in Fig. 3.5. Based on the qualitative picture described in §3.4 that photoexcitation generally narrows or closes the gap, reducing the bottleneck effect of electron relaxation, we expect that the beyond-fixed-band effects will lead to even faster relaxation, and to a larger final number of excited particle-hole pairs.
3.4 Soft bands in Hartree-Fock theory
In density function theory, the electronic potential is a self-consistently determined functional of the electron density, so that changes in the electron distribution will lead to changes in the band structure. This effect is greatly enhanced in extended DFT theories such as DFT+ and DFT++ because, in particular, the relative energetics of the different orbitals depends strongly on the orbital occupation matrix. This strong dependence may lead to photoinduced phase transitions if photoexcitation changes the occupancy sufficiently.
In the specific case of VO2, since the wavelength of the pump laser is typically 800 nm ( eV), the pump laser typically changes the electron distribution among the V- orbitals (see Fig. 1), but does not change the total -count or the real-space charge density significantly. We therefore argue that we may analyze the effects of photoexcitation using Eq. (3.1) with and left unchanged, but with now determined by the nonequilibrium distribution of electrons over orbitals, i.e., the Kohn-Sham Hamiltonian becomes Hartree-Fock shifted to
[TABLE]
where is the change of due to the change of the orbital occupation matrix [see Eqs. (3.2)–(3.4)] under photoexcitation. Equation (3.18) implies that the electronic band structure becomes soft in the sense that the conduction band floats down when its occupancy increases and the valence band floats up when its occupancy decreases under photoexcitation. This general picture shows that photoexcitation has the potential of closing the Mott gap and driving an insulator-metal transition, thus giving rise to new electronic phases. The total energies of different electronic states can be compared using
[TABLE]
where the expectation value is now taken using the nonequilibrium distribution. We will later use Eq. (3.19) to construct an energy landscape for nonequilibrium VO2 that will be used to interpret the experiments of [33].
3.4.1 Nonequilibrium phase transition to a metastable metal
In §3.3, we showed that electrons in VO2 relax on a sub-picosecond time scale to a thermal state with a well-defined instantaneous temperature. Here, we investigate whether the changes in orbital occupancies due to photoexcitation can lead to significant changes in the band structure, in particular the HOMO-LUMO gap. Because the system relaxes rapidly to a thermal state, we can avoid solving a dynamical Hartree-Fock equation and consider a Hartree-Fock theory in thermal states only.
We note at the outset that obtaining an insulating state in VO2 requires two effects. First, the dimerization (enhanced by an intersite correlation effect) splits the band into bonding and antibonding portions. Second, the on-site interaction produces a level splitting between and the orbitals. The dimerization gives the possibility of having a filled band, and the level splitting ensures that the band lies far enough below the other bands that it is indeed fully occupied. The equilibrium phase transition from the insulating to the metallic state involves a change in the crystal structure, removing the dimerization. An alternative possibility is that at fixed structure a population inversion of the and the bands, driven by photoexcitation, would lead to a reversal of the energy ordering, so that the non- (weakly) dimerized bands would lie lowest, creating an metal phase.
To investigate the possibility of this metal phase, we first apply the soft-band Hartree-Fock theory at temperature by calculating the shift of the bands using Eq. (3.18). We start from an occupation matrix with a high occupancy, and find at eV, eV, and eV that our system relaxes back to the conventional insulator phase shown in Fig. 3.2 in the Hartree-Fock iterations. However, at slightly increased values of , i.e., eV and 5 eV, the iterations bring us to a new self-consistent state with a high (low ) occupancy and no gap at the Fermi level: an metal phase is found!
The projected density of states of the metal phase is plotted in Fig. 3.6. We see that the density of states at the Fermi level is nonzero, so within a band picture the state is metallic. Also, the orbitals are now substantially above the Fermi level, and the bonding-antibonding splitting of the orbitals is less, reflecting the decrease in the intersite Fock terms due to the depletion of the band.
We next construct a cut across the energy landscape in Fig. 3.7 as a function of orbital occupancies with the insulator and metal phases as its local minima. To do this, we first determine for and eV the (full basis) real-space density matrix of an intermediate state as a linear interpolation between the density matrices of the two local minima. Then we introduce k-independent Lagrange multipliers to the Kohn-Sham Hamiltonian , which are adjusted so that the band occupancies reproduce this interpolated density matrix. The states obtained are the minimum energy states subject to the constraint of a linearly interpolated real-space density matrix. The energy is then evaluated by Eq. (3.19) using without the Lagrange multipliers. The resulting curve, although not necessarily the minimum energy path between the insulator and metal phases, should give a reasonable representation of the energy barrier between them. For eV, the metal phase is a state in the ghost region of the iterative Hartree-Fock dynamics with the slowest evolution, and the energy curve is plotted following the evolution to the insulating ground state. The extrapolated states at any value of cannot be obtained by linear extrapolations of real-space density matrices, as these can have occupancy eigenvalues not between [math] and . Instead, the states are obtained by tuning the orbital energies of and with respect to using the Lagrange multipliers to further raise or lower the occupancies of the orbitals.
While Fig. 3.7 shows that the metal phase has higher energy at , we find that at the state may be favored. Figure 3.8 plots the calculated HOMO-LUMO gap as a function of the energy deposited by the pump laser into the sample for realistic range of parameter values. Because the electrons equilibrate rapidly, this is equivalent to plotting against temperature, although the temperature-energy relationship is not quite linear and depends on which phase the system is in.
Two qualitatively different behaviors are seen in Fig. 3.8. For eV, eV, there is no phase transition. The bonding band in Fig. 3.2 shifts up and the and bands shift down as temperature rises, and eventually the band gap between them is closed.
But there is always a unique stable state at every temperature or energy . Similar effects are seen for eV, eV except that the curve drops more slowly and the gap closes at a slightly higher temperature. The behavior is very different for eV, eV. When the overlap of the band with and bands (indicated by a negative gap in Fig. 3.8) exceeds a certain threshold (the small circle on the green curve), the band structure undergoes a first-order phase transition to a state with an inverted population and thus a negative HOMO-LUMO gap (metallic state) occurs. Near the discontinuity, the - curve in Fig. 3.8(a) shows a singularity, but the - curve in Fig. 3.8(b) is not singular.
The metal phase may be metastable (correspond to a local energy minimum) even if it is not thermally reachable. Figure 3.9 summarizes the situation, showing by red squares (blue diamonds) the region where a thermally driven transition to the metal phase occurs (or not), and by Roman numerals (II and III) the regions where the metal phase is locally stable and (I) where only the insulator phase is locally stable. Region II is the hysteretic range in which a thermally excited metal phase could survive but the insulator-to-metal transition would require and to reach Region III.
Compared with the input energy – eV per unit cell (4 VO2) estimated in §3.2 for the experiment in [33], the transition point in Fig. 3.8(b) corresponds to a fluence about four times larger than that at which the putative metal phase was observed. At the experimental fluence level, the theory indicates that the HOMO-LUMO gap is only slightly reduced from eV in the insulating ground state to – eV [Fig. 3.8(b), eV, eV]. This discrepancy with experiment may be due to limitations of the Hartree-Fock theory, which often does not calculate the energy of correlated electrons or locate the phase boundaries accurately.
3.5 Summary and conclusion
This chapter presents a theoretical study of photoexcited VO2 motivated by Morrison’s experiment of a long-lived metallic phase created by photoexcitation in a material with a crystal structure associated with insulating equilibrium behavior. We used a band-theory-based Hartree-Fock mean-field methodology combined with quantum Boltzmann equation treatment of excited-state kinetics. The key findings of our study were (a) very rapid (fs) relaxation of the photoexcited carriers to a pseudothermal state characterized by a common temperature but different chemical potentials for the electron and hole distributions, (b) a rapid ( fs) relaxation to a thermal state with a well-defined common temperature and chemical potential, and (c) the existence of a metallic phase which is metastable at temperature and can become favored at higher temperatures (or laser fluence levels). A recent experimental report [52] of rapid ( fs) collapse of the electronic gap is consistent with our calculations in §3.3.
The key approximations of our work are the Hartree-Fock plus Fermi’s golden rule treatment of the electron-electron interactions, and neglect of electron-phonon coupling beyond thermal energy exchange. We believe that these approximations are not crucial. The important conclusion of the quantum Boltzmann and Fermi’s golden rule studies of the dynamics is that thermalization of the excited particles proceeds much faster than experimental time scales, so that experimentally relevant issues, in particular the existence of a metastable metallic state, can be addressed using steady-state arguments. Further, the local stability of the metallic phase means that as phonons take energy out of the electronic system, the system may simply remain in this phase over a long time determined by nucleation kinetics. The conclusion seems very likely to survive the inclusion of higher order effects in the dynamics. Hartree-Fock theory is normally reliable for the identification of phases, although the estimates of the locations of phase boundaries may be inaccurate. The results presented here should be viewed as indicating the theoretical possibility of a metastable metallic phase for reasonable parameters. Further investigations of this metallic phase, including more reliable determination of the phase boundaries, investigation of the processes by which the metastable state might decay, and the study of the evolution of the lattice structure, would be of considerable interest for future research.
Chapter 4 Towards a real-time impurity solver: quench dynamics
The out-of-equilibrium simulation of a quantum many-body system can be treated semi-classically as in Chap. 3 using the quantum Boltzmann equation based on Fermi’s golden rule, or fully quantum mechanically using some many-body wave function evolved in real time. Here in Chaps. 4 and 5, we will study the methods of simulating the out-of-equilibrium dynamics of the Anderson impurity model (AIM) [53], a single spin-degenerate orbital with an intra-orbital Hubbard interaction coupled to a bath of noninteracting orbitals. This model is of fundamental importance both in its own right as a solvable [54, 55] interacting electron model and, as discussed in §1.4, as an auxiliary problem for the dynamical mean-field theory [17, 56]. While a lot of work has been done to develop imaginary-time solvers [57, 58] for the Anderson impurity model to study its equilibrium properties at finite temperatures, it has been a long-standing challenge to develop efficient real-time impurity solvers for doing out-of-equilibrium simulations. There are various candidate methods towards this goal, including wave-function-based methods such as exact diagonalization (ED) [59, 60] and density matrix renormalization group (DMRG) [61], and Green’s-function-based methods such as the quantum Monte Carlo algorithm [62].
In our work presented in Chaps. 4 and 5, we use the density matrix renormalization group (DMRG) method [20], a powerful numerical technique for solving one-dimensional electron problems. In DMRG, the wave function of the system is represented by a matrix product state (MPS). Every matrix in the MPS corresponds to a local degree of freedom in some single-electron basis. The main challenge is to find the right basis of the bath orbitals so as to slow down the growth of the entanglement entropy of the MPS. Here in Chap. 4, we use DMRG to represent the noninteracting bath orbitals in energy space to study the quenched Anderson impurity model starting from a nonequilibrium direct-product state. In Chap. 5, we generalize our method to the periodically driven Anderson impurity model. The work of this chapter has been published in [63]. The work of Chap. on the driven model has been submitted and is viewable on arXiv [64].
4.1 Theory and formalism
Our theoretical studies are focused on the single-impurity Anderson model (SIAM) with one impurity -orbital coupled to a noninteracting bath. The Hamiltonian is given by
[TABLE]
The orbital has a Hubbard in its Hamiltonian and the bath has noninteracting bath orbitals. The two systems hybridize via the one-body hopping terms in . The bath orbitals are labeled by and the two spins and of electrons are labeled by . For simplicity, we take the impurity-bath coupling amplitudes to be -independent. We define the bath density of states as and consider a semicircle DOS with a half band width . The initial state that we consider is a direct-product state
[TABLE]
where the Fermi-sea state of the bath is initially half-filled. The -orbital energies and are chosen to be symmetric about the Fermi level at [math]. The situation is depicted in Fig. 4.1. The formalism generalizes to a mixed initial state with a direct-product density matrix , where satisfies the Wick’s theorem, but we will focus on a pure initial state here.
Our numerical method requires a truncation of the bath to a finite number of bath orbitals. To choose the best truncation, we calculate the hybridization function
[TABLE]
on the Keldysh contour [15]. One may refer to Appendix C for a derivation of . Here the contour goes from to a sufficiently late time and then back to without the imaginary-time part. is the contour-ordering symbol. The mean value is taken with respect to . With our choice of the semicircular DOS and constant , the hybridization function can be analytically found for to be
[TABLE]
where , is the 1st-order Struve function and is the 1st-order Bessel function. The symbols and refer to Keldysh-contour ordering. Then we fit the hybridization function to that of a finite bath with only orbitals, i.e.,
[TABLE]
where if and [math] if . In fitting Eq. (4.10) with Eq. (4.11), all real parameters and are varied to minimize the least-square error up to a maximum time. The occupancies are chosen to be either 0 or 1 to fit the and parts independently and to make the initial state of the finite bath a Slater determinant. This is possible even if the original bath was at a nonzero temperature. Since our bath is particle-hole symmetric, we choose to be even to preserve this symmetry. The number of bath orbitals controls the maximum time up to which the exact hybridization function is reproduced with good accuracy. For example, bath orbitals are enough to reach and orbitals can reach . Adding more orbitals increases the maximal time that can be reached, but does not significantly improve the accuracy of the fit at shorter times.
4.2 Numerical method
We use DMRG/MPS methods to carry out the time evolution. We represent the wave function as an entangled state between the impurity orbital and the bath, which is linearly expanded as
[TABLE]
where sums over the impurity states , , , and . Every bath state is a normalized matrix product state (MPS). The coefficients are real and nonnegative. Eq. (4.12) is a Schmidt decomposition of between the orbital and the bath if is a simultaneous eigenstate of and , the total numbers of spin-up and spin-down electrons. This representation differs from the conventional DMRG in that it removes the orbital from the MPS, enabling analysis of the entanglement among the bath orbitals. We evolve the wave function in the interaction picture of using
[TABLE]
where is the time-ordering symbol and
[TABLE]
where is the opposite spin of . The main advantage of the interaction picture is that typically has a narrower spectral radius than (bath bandwidth large compared with impurity level width ), so one can choose bigger time steps in the simulation. We evaluate Eq. (4.13) by discretizing the time evolution into time steps . The Hamiltonian used during the time step centered at is
[TABLE]
with the coupling amplitudes
[TABLE]
The errors of both the mid-point Hamiltonian and the time-averaged Hamiltonian are . The latter choice is preferred if the bath bandwidth is large compared with the -level width, because the very high and very low-energy bath orbitals are suppressed by the sinc function. To apply the Hamiltonian to the wave function in Eq. (4.12), we work in the Jordan-Wigner transformed representation with the orbital being the first orbital ( and having no Jordan-Wigner signs). The Hamiltonian in Eq. (4.15) is rewritten as
[TABLE]
where the is the Jordan-Wigner transform of . The two operators are related by
[TABLE]
so that the operators and with commute. We do the same Jordan-Wigner transform for the two spins of the same orbital, so that and still anticommute. But this is easy to handle with a local matrix. The bath operators in Eq. (4.18) is then represented by a matrix-product operator (MPO)
[TABLE]
where the matrix is left-multiplied by to pick the second row, and the matrix is right-multiplied by to pick the first column. The MPO has a bond dimension of 2. We can similarly express in terms of . The Hamiltonian can then act on following DMRG routines [20]. The final evolution scheme is given by
[TABLE]
with the exponential Taylor expanded into a 4th-order polynomial of . The narrow spectral radius ensures good unitarity of the 4th-order truncation. Since the bath operators and are long-range, we cannot locally exponentiate the Hamiltonian as in the time-evolving block decimation (TEBD) [65] method and have to Taylor expand the exponential.
We adjust the truncation error tolerance of the singular value decomposition (SVD) in DMRG according to the MPS norm so that the higher-order terms of do not take much time to calculate. The error tolerance in our code for a norm-1 MPS is set to (in terms of probability loss) per SVD truncation. This number is multiplied by (norm)-2 for MPSs with smaller norms (coefficients). If the adjusted error tolerance becomes greater than (which happens if the MPS norm is very small), the MPS is truncated to a product state. The coefficients in in Eq. (4.12) are normalized at the end of every time step. The -occupancy produced for is found to agree in decimal places with a Slater-determinant-based noninteracting code.
We parallelize the calculations of the 4 MPSs in Eq. (4.12) on processors and also use the total numbers and of spin-up and spin-down electrons as symmetries to further speed up the calculation.
4.3 Physical results compared with analytical theories
In this section we show some results obtained for the interacting SIAM with using the method and other model parameters described in previous sections. The impurity-bath coupling . This is the parameter range of interest. The impurity level width remains smaller than the band width while the Kondo temperature [66] can change by orders of magnitudes. The Kondo temperature is a Hubbard induced energy scale that measures the spin relaxation rate in the near-equilibrium regime of the SIAM in the small (or Kondo) limit.
Fig. 4.2 (a) shows the charge relaxation dynamics, obtained by starting from an initially empty orbital and a half-filled Fermi-sea state for the bath. Our choice of particle-hole symmetric parameters ensures that always equilibrates to so long as the impurity-bath coupling is not big enough to form a bound state on the impurity. We see in agreement with previous work [61, 67] that the charge equilibration proceeds relatively rapidly. The reciprocal of the time it takes to reach is plotted in Fig. 4.3 (a). At small , is inversely proportional to the -level width . For , the rate of equilibration crosses over to approximately linear in and the equilibration process in Fig. 4.2 (a) becomes more oscillatory as we are approaching the formation of a bound state on the impurity. The variation of charge equilibration rates with can be seen in calculations performed for a noninteracting SIAM (), suggesting that the charge relaxation physics is essentially due to hybridization. The Hubbard does not change the behavior of the model qualitatively.
Fig. 4.2 (b) shows the spin relaxation dynamics obtained by starting from , a fully spin-polarized orbital and the same half-filled Fermi-sea state of the bath. The magnetization relaxes much more slowly than the charge, again in agreement with previous results [68, 69]. The asymptotic behavior of v.s. shows approximately an exponential tail, with the relaxation rate plotted in Fig. 4.3 (b). is estimated by fitting v.s. to a straight line for , where is the maximum time reached in the simulation. The solid red line is a trend line. We also show as the dashed green line the analytical result — the Kondo temperature calculated using the formula in [66] and interpreted as a relaxation rate.
The Kondo result has a similar magnitude and dependence to the calculated results. The numerical differences at large arise from relaxation processes associated with valence fluctuations not included in the Kondo limit, while the more pronounced differences at small are an intermediate asymptotics effect. For small , even at the very long times () accessible to our method, the magnetization is still substantial, so the Kondo-limit expression, which gives the linear response relaxation for small magnetization (), is not applicable. Evidently, the nonlinear response (relaxation of a finite ) is stronger than the linear response. Developing a theory of the relaxation in the small and intermediate regime is an interesting open question. For intermediate , the theoretical result is within a factor of of the numerical one with the differences likely arising from the convention used for the Kondo temperature .
4.4 Logarithmic growth of entanglement entropy
A remarkable feature of the simulations reported here is the long time scales that can be reached; these time scales are necessary to reveal, for example, the magnetization decay. As we show in this section, this is possible because the maximum entanglement entropy of the 4 bath MPSs in Eq. (4.12) grows only logarithmically during the simulation, which means the long times are not exponentially hard to reach, but are of only polynomial-time complexity.
4.4.1 Entanglement entropy growth in SIAM
In this section, we compare the maximum entanglement entropy of the interacting SIAM () with a noninteracting SIAM () with at the Fermi level. Both models start from the same initial condition with an empty orbital and a half-filled bath in Fig. 1. Results of the entanglement entropy are shown in Fig. 4. The entropy growth starting from a spin-polarized impurity is numerically found to be also logarithmic but takes smaller values.
The curves in Figs. 4.4 (a) and 4.4 (b) are obtained in slightly different ways. Fig. 4.4 (a) shows the results obtained in a noninteracting simulation (using a Slater-determinant-based code) of bath orbitals all coupled to one empty orbital at the Fermi level. We then plot the entanglement entropy between the 500 bath orbitals below the Fermi level with the other 501 orbitals up to . The Fermi level is found to be close to the maximum entropy cut of the system. The data shows a logarithmic growth of entanglement entropy at all values of the impurity-bath coupling . The slopes of the curves at long times are the same; only the offset and the transient growth depend on .
Fig. 4.4 (b) shows the results of the interacting SIAM (, ) simulated using the method of Sec. 4.2. Every point on a curve is obtained in a different simulation with a different bath size . In a simulation up to time , the hybridization function is first fitted up to with the minimal number of bath orbitals (typically between ) needed to keep the root-mean-square error (RMSE) of the fit under . Then the maximum entanglement entropy seen on all bonds of the 4 bath MPSs encountered during the simulation from [math] to is plotted v.s. . Notice that may be encountered before due to the finite bath effect. So Fig. 4.4 (b) takes into account the possibility of using the finite bath effect to limit entropy growth. But still the logarithmic growth of entropy and the independence of the steady-state slope of v.s. on the impurity-bath coupling are the same as in the noninteracting SIAM in Fig. 4.4 (a). These two properties mean , and therefore the bond dimension . Hence, the interacting SIAM can be simulated in polynomial time of .
4.4.2 Analysis of entropy growth
To understand the logarithmic growth of entropy, we consider a noninteracting chain model, as is shown in Fig. 4.5. In this model, the impurity is coupled to two semi-infinite chains. We choose a constant hopping amplitude between the bath sites in each chain. By adjusting the on-site energy difference of the two chains, we can vary the densities of states, obtaining either overlapping, gapped or just touching spectra. Our computation of the entanglement entropy across the impurity site shows that we have linear growth, logarithmic growth, and saturation, respectively. In the numerical test we did, chain a was initially empty and chain was initially full. But the conclusion is found to hold for randomized initial occupancies, too. A similar noninteracting model with two semi-infinite chains directly connected via a modified hopping amplitude was studied in [70, 71] in the formalism of conformal field theory. In our model the two chains are connected via the orbital.
The logarithmic growth of is seen at a critical , at which the density of states (DOS) of the two semi-infinite chains touch at only one energy point. When , the system is gapped and entropy growth saturates. This can be explained by the lack of energy eigenstates that are extended in both regions and , which then means that particles (or holes) that are originally in cannot go into and vice versa beyond a penetration depth determined by the gap, which then puts an upper bound on the entanglement entropy between and . This energy barrier works for a general initial occupancy. Starting from any product state, so long as the semi-chains and are gapped, the entropy must saturate.
When , there is a finite overlap of the DOS of the two semi-infinite chains and we see a linear growth of entanglement entropy in Fig. 4.5. In rare cases this does not happen. For example, for a uniform chain and , the entropy growth is logarithmic rather than linear. But this behavior depends on the initial occupancy. If the occupied sites are randomized, or if the model parameters are slightly modified to deviate from a uniform chain, the expected behavior of a linear growth of the entanglement entropy is seen between and . The energy criterion guarantees that particles do not enter the forbidden regions of a noninteracting bath. But once the energy barrier is not at work, it is difficult in general, though not impossible, to organize the migrated particles into a low entanglement entropy state to make the MPS matrices small.
The logarithmic growth of entropy in Fig. 4.4 (b) can be understood as the result of arranging the bath orbitals in the MPS in energy order, so that at any bond of the MPS, the left and right parts of the bath degrees of freedom always have touching energy spectra. This argument applies to an interacting model, too, because the bath is still noninteracting, and the Hubbard U only reduces the chance for the impurity – the only bridge via which the bath orbitals can indirectly hop to one another – to be doubly occupied, thus reducing its bridging efficiency. The bath entanglement entropy of an interacting SIAM is therefore upper bounded by that of a noninteracting SIAM from this picture.
4.4.3 Bath in chain geometry
So far we have been working in the star geometry of the bath. Bath orbitals do not hop to each other directly. They only do so via the impurity. The diagonalization of bath orbitals in energy space leads to a logarithmic growth of entanglement entropy, according to the energy criterion in the previous section. In this section, we would like to emphasize again that the energy criterion is a sufficient but not necessary condition for the entropy to grow slowly. The example to give here is the evolution of the quenched SIAM in the chain geometry of the bath. The impurity is the head of the chain, which is directly connected to only one bath orbital, which in turn is connected to another bath orbital, and so on so forth. One can go from the star geometry to the chain geometry via Lanczos tridiagonalization starting from the impurity orbital, and from the chain back to the star by diagonalizing the bath. More details of the two geometries can be found in [61].
Starting from the initial state with being the same filled Fermi-sea state as in Fig. 4.1 transformed to the chain geometry, the maximum entropy on the chain (the entanglement entropy between the left and right parts of the chain at the maximum entropy cut) is still found to grow logarithmically. Fig. 4.6 shows the result of a noninteracting calculation. The initial occupancies on the chain are spatially uniform. Every site has an occupancy of per spin except the empty impurity. The entanglement entropy between sites and on the chain are then plotted in Fig. 4.1 as a function of at equal intervals of time. On top of the logarithmic background of of the equilibrium state , an entropy peak propagates like a soliton from the impurity down the chain at a speed . The maximum entanglement entropy (height of the peak) therefore increases with time logarithmically, even though there is no separation of energy spectrum on the chain, i.e., partition of the bath into different regions with different energies like in the star geometry.
Starting from an inverted half-filled Fermi-sea state with the highest energies initially occupied, the same logarithmic growth of entropy in the chain geometry is seen due to particle-hole symmetry. But starting from a product state with random initial occupancies of the bath orbitals in the star geometry, the entanglement profile transformed to the chain geometry becomes very high () even at . Also, a linear growth of entropy is seen starting from a product state in the chain geometry with randomized initial bath occupancies (see Fig. 4.7, blue line), while in the energy-ordered star geometry, the entropy growth (green line) is still logarithmic under the same condition. These results demonstrate that the logarithmic entropy growth in Fig. 4.6 is not guaranteed by the MPS basis, but is due to the initial filled Fermi-sea state. For such a special initial state, the star geometry does not have a big advantage over the chain geometry, as they both give a logarithmic growth of maximum entanglement entropy. The benefit of the star geometry is its good behavior for more general initial states.
It is important to point out, as is shown in Fig. 4.7, that the star geometry alone does not guarantee a logarithmic entropy growth. The order of the bath orbitals in the MPS matters. The initial occupancies affect the transient growth of entropy, while the asymptotic entropy growth is determined by the ordering of the bath orbital energies. The steady-state growth of is logarithmic if the bath orbitals in the MPS are energy-ordered and linear if the bath orbital energies are randomly shuffled.
4.5 Double-impurity model
In this section, we show that the logarithmic growth of entanglement entropy is not limited to the single-impurity Anderson model by doing a noninteracting simulation of a double-impurity Anderson model. The most general noninteracting double-impurity Anderson model can be pictorially represented in Fig. 4.8. Fig. 4.8 (a) is in the basis in which the 2 impurity orbitals and all bath orbitals are diagonal, which is the double-impurity version of the star geometry. Fig. 4.8 (b) shows the double-impurity version of the chain geometry by Lanczos tridiagonalizing the star geometry in Fig. 4.8 (a) starting from the two impurities. One can also tridiagonalize the bath orbitals above and below the Fermi level separately (Fig. 4.8c) to obtain the double-impurity generalization of the two semi-infinite chain model in Fig. 4.5. Since the left and right semi-chains have touching energy spectra, a logarithmic growth of entropy is expected as a critical behavior between linear growth and saturation of entropy, as is discussed previously.
Figure 4.9 shows a sample result. We chose a half-filled bath with a semicircle DOS the same as Fig. 4.1, and put two orbitals at ( is the half band width) with hopping to mimic typical crystal field splitting. The two orbitals are equally coupled to all bath orbitals. In the basis in which the two orbitals are diagonalized, their orbital energies are and the original hopping makes the two orbitals now couple to the bath differently, which is more realistic. Then we plot the entanglement entropy between the initially occupied bath orbitals and the rest of the system.
The double-impurity model has a richer dynamics than SIAM. Since both impurities are initially empty, the one below the Fermi level leaks a hole into the bath, leading to a short-term entropy peak. The steady-state growth of is still logarithmic, but the slope of v.s. is not constant. This is because the two orbitals are not at the Fermi level (one is above and one is below). Their distances in energy to the Fermi level relative to the impurity-bath coupling determine the slope, which approaches a maximum for the case of a d-orbital at the Fermi level () as gets large.
The logarithmic growth of entropy again shows that the quenched multi-impurity model is not exponentially hard in DMRG simulations, but is of only polynomial-time complexity. Whether the conclusion still holds for interacting models needs further investigation in DMRG, especially for those multi-impurity models with non-density-density (spin flipping and pair-hopping) terms, whose entanglement entropies need not be bounded by the corresponding noninteracting models.
4.6 Summary and conclusion
We have studied the growth of entanglement entropy in quenched Anderson impurity models. It is found that the entropy growth is determined by the representation of the bath orbitals in the matrix product state (MPS). The Hubbard on the impurity orbital does not change the qualitative behavior of the steady-state growth of entanglement entropy of the bath MPSs. The crucial feature controlling the entropy growth is the overlap in energy of the density of states of the two parts of the bath at the maximum entanglement entropy cut. In the star geometry of energy-ordered bath orbitals, the touching-spectra condition is satisfied at every bond, so the maximum bond dimension is power-law in . The power is upper bounded by the case of a half-filled -orbital at the Fermi level and does not grow with the impurity-bath coupling, which allows a simulation of the long-time dynamics of the quenched impurity models in polynomial time. The conclusion is likely to generalize to quenched Anderson impurity models with multiple impurities.
The growth of entanglement entropy of an interacting quantum system and the associated computational cost has been studied previously [72, 73] in terms of the integrability of the quantum model. Our study looks at the problem from a different perspective. We focus on a special class of quantum models — the impurity models — and think of the growth of entanglement entropy among the bath orbitals. Because of the sparsity of interactions in the model, the entropy growth in the noninteracting bath is controlled by the energy partitioning of the bath and the localization of bath electrons to the energies they belong to. Since the new criterion of energy-partitioning the bath is not related in obvious ways to the integrability of the whole model (bathimpurity), hopefully this new view of entropy growth of complexity can help us find new polynomial-time solvable models, parameter ranges, and special initial conditions that are not covered by the integrability criterion.
Chapter 5 Towards a real-time impurity solver: driven dynamics
With the development of experimental technology, manipulating strongly correlated electrons using a laser-induced oscillating field is becoming possible [74, 75, 76]. This makes it interesting to study the driven dynamics of strongly correlated systems, i.e., evolution of systems with time-dependent Hamiltonians. In Chap. 4, we have studied the quenched Anderson impurity model [53] and have found that the star geometry with energy-ordered bath orbitals in the matrix product state (MPS) proves to be an efficient solver that simulates the system in polynomial time due to the logarithmic growth of maximum entanglement entropy over the MPS. In this chapter, we consider an Anderson impurity model with an oscillating -orbital energy. We let the -orbital energy oscillate in a square wave across the Fermi level of the half-filled noninteracting bath. It is found that when the driving period is short so that the Floquet-Magnus expansion converges, the energy-ordered bath MPS works as well as it does for the quenched model. But when the critical period is exceeded to make the Floquet-Magnus expansion diverge, the original algorithm of using energy-ordered MPS exhibits linear growth of entanglement entropy and therefore exponential time complexity. To overcome this problem, we tried the quasi-energy-ordering algorithm and found that the long-term entropy growth gets slowed down, but at the cost of a faster short-term entropy growth for not arranging the bath orbitals in the MPS in energy order. So there is a tradeoff between the short-term and long-term computational costs. Long driving periods would favor energy ordering while short driving periods above and comparable to the critical period would favor quasi-energy ordering. Below the critical period the two methods become identical, as there is no energy aliasing effect.
5.1 Theory and method
We begin with the general formalism of a single-impurity Anderson model (SIAM) with general time-dependent model parameters. As is pointed out in [18], this is the type of Hamiltonian that could arise in a nonequilibrium single-site dynamical mean-field theory (DMFT). The time-dependent Hamiltonian is given by
[TABLE]
where and is the spin label. We go to the interaction picture of . The part in the interaction picture becomes
[TABLE]
where is the time-ordered unitary evolution from [math] to and . Since does not couple the orbital to the bath, each bath orbital evolves independently in the interaction picture as given by
[TABLE]
with denoting the opposite spin of . Notice that does not evolve in the interaction picture of and that commutes with , which together lead to Eq. (5.6b). The 4-MPS scheme developed in Chap. 4 can be applied to the general time-dependent SIAM with only one modification: in every time step , the time-averaged Hamiltonian is
[TABLE]
with the effective hopping amplitudes given by
[TABLE]
Here we assume that in one time step , the orbital energies , and Hubbard do not change by much, so one still obtains the function after the time average. The wave function is still evolved according to
[TABLE]
with the exponential Taylor expanded to 4th order of to ensure good unitarity.
Up to now everything has been general for the single-impurity Anderson model. In this paper, we consider the evolution starting from a product state
[TABLE]
where is a half-filled Fermi-sea state of the bath with a semicircle density of states (DOS) as shown in Fig. 5.1. The bath orbitals have fixed energies and fixed equal hopping amplitudes to the impurity orbital. The Hubbard on the orbital is also fixed. The only time-dependent quantity is the -orbital energy
[TABLE]
which oscillates in a square wave every half driving period . Physically, we are interested in the local quantities on the orbital. To this end, the bath can be fitted by a finite number of bath orbitals to reproduce the hybridization function in the thermodynamic limit up to a maximum time proportional to . This fit is independent of the driving of the -orbital energy and the Hubbard . Computationally, we want to study the growth of entanglement entropy of the bath, which determines the time complexity of the problem.
5.2 Noninteracting Results
Let us first do some cheap calculations of the entanglement entropy growth of the noninteracting SIAM using a standard Slater-determinant-based method to scan the complexity diagram. The driven 4-MPS scheme developed in §5.1 will be used in §5.3 for simulating the interacting SIAM. In this section, the Hubbard and the initial state is , an empty -orbital and a half-filled Fermi-sea state in Fig. 5.1. The impurity-bath coupling is fixed. Bath size .
5.2.1 Energy-ordered bath
We use the entanglement entropy between the bath orbitals below the Fermi level and the rest of the system to estimate the maximum entanglement entropy that would be encountered in an MPS-based simulation when the bath orbitals are energy-ordered. We find that for long driving periods , the convergence radius of the Floquet-Magnus expansion (see Appendix F), a small amplitude could change the logarithmic growth of entropy to linear. This is shown in Fig. 5.2, where we did a simulation with bath orbitals and driving period . The entanglement entropy between the bath orbitals below the Fermi level and the rest of the system is plotted in Fig. 5.2 over time. The quenched model exhibits a logarithmic growth of entanglement entropy over time and while the entanglement entropy growth in the periodically driven model is linear.
The critical driving period , or over the band width, separates the logarithmic growth () and linear growth () of . In Fig. 5.3, we plot the maximum growth rate of entropy v.s. the period as an envelope of the growth rate v.s. curves at fixed driving amplitudes . Each of these curves is tangent to the envelope at some points and they all intersect with zero at the same critical period . The linear growth of entanglement entropy in a driven SIAM can be intuitively understood in an entropy pumping picture. The up and down motion of the orbital acts as an elevator that transports some electrons from the occupied bath orbitals to the unoccupied bath orbitals (and holes in the opposite direction). So if the entanglement entropy increases by a constant in every period, the linear growth rate of would then be .
However, for short periods , the linear growth of entropy cannot be maintained. To understand this critical period, we consider the Floquet Hamiltonian defined by
[TABLE]
where corresponds to respectively. The time-independent Floquet Hamiltonian reproduces the unitary evolution of the time-dependent system over full periods. It turns out is the convergence radius of the Floquet-Magnus expansion of in terms of and . Within the convergence radius and for small , we have
[TABLE]
where is the SIAM Hamiltonian with , is the adjoint representation of , and is defined via its Taylor expansion. Eq. (5.15) can be derived using the formalism given in Appendix F. The convergence radius of Eq. (5.15) is . We expect this to hold also for the interacting SIAM, because in the thermodynamic limit , the spectral radius is mainly determined by the band width of the bath DOS (unless a bound state is formed on the impurity). As a result, is equal to the half band width of the bath. Since is singular at , the series expansion of Eq. (5.15) fails to converge if , i.e. .
Once the critical period is exceeded, surprising new physics emerges. For driving periods , numerics shows that the bath orbital energies in are now aliased to , breaking the original ordering of the bath orbitals. This situation is shown in Fig. 5.4. The inter-bath-orbital hopping amplitudes remain very small. The ascending order of bath orbital energies is violated because of energy aliasing. This gives overlap of (aliased) energy between the occupied and unoccupied bath orbitals and thus a linear growth of entanglement entropy over the periods becomes possible.
5.2.2 Quasi-energy-ordered bath
What happens then, if one reorders the bath orbitals in the MPS in ascending order of quasi-energy rather than energy in case of the driving period ? The initial state in the star geometry remains a product state (an MPS with bond dimension ). We use the entanglement entropy between the bath orbitals with negative quasi-energies (within ) and the rest of the system to estimate the maximum entanglement entropy that would be encountered in an MPS-based simulation when the bath orbitals are quasi-energy-ordered. becomes the used in the previous subsection when the bath orbitals are energy-ordered.
We redo the same simulation as in Fig. 5.2 using bath orbitals ordered by their quasi-energies of . The same results of as in Fig. 5.2 (b) for the quenched and driven models are reproduced. The entropies of the energy-ordered simulation in Fig. 5.2 (a) are compared with the new results in Fig. 5.5 and the time is put on log scale. It is found that the growth of is logarithmic for both the quenched and driven models. This is because the Floquet Hamiltonian is now energy-ordered (as opposed to Fig. 5.4). But the driven model is still harder to simulate than the quenched model, because the slope of the v.s. curve is greater for the driven model (see purple line) even if the quasi-energy-ordering method is used.
For the quenched model, the steady-state slope of v.s. is unchanged when the bath orbitals are quasi-energy-ordered. By comparing the blue and green lines in Fig. 5.5, we see that the steady-state intercept is shifted up by a constant , which is found to be approximately proportional to , as is summarized in Fig. 5.6 (a). This is the price to pay for not ordering the quenched bath by energy, which is better than a randomly shuffled bath (see Fig. 4.7 in Chap. 4), with entropy growing linearly with time .
The driving amplitude changes the slope of the v.s. curve. Fig. 5.6 (b) shows how the slope increases from that of the quenched model ( at fixed is equivalent to quench) to unboundedly large values proportional to . This indicates that the leading-order term in the entropy is
[TABLE]
where depends on but is found to be bounded (see Fig. 5.7). At very large , the coefficient goes down, which is likely to come from the bound state formed on the impurity. Eq. (5.16) means that the bond dimension in an MPS-based simulation using the quasi-energy-ordered algorithm is . The time complexity of the singular value decomposition (SVD) step is then .
Since the power of for the quasi-energy-ordered algorithm is unbounded for long driving periods , the complexity is still beyond polynomial time. Another drawback of quasi-energy ordering is delocalization of maximum entanglement entropy throughout the MPS, while in energy-ordered MPSs, the maximum entanglement entropy tends to concentrate near the Fermi level. This gives the quasi-energy-ordered algorithm a prefactor of the bath size . The quasi-energy-ordering method needs to overcome these short-term drawbacks before its long-term benefit becomes competitive to the energy-ordering method.
5.3 Interacting results
In the previous section, we have been estimating what would happen in an MPS-based simulation using a noninteracting (Slater-determinant-based) code. Now let us do some real MPS-based simulations of the interacting SIAM using the 4-MPS method of §5.1. We choose a fixed Hubbard and the impurity bath coupling is the same as in §5.2. We use bath orbitals to fit the hybridization function of the continuum bath DOS in Fig. 5.1 with good accuracy up to following Chap. 4. The SVD truncation error tolerance was . Noninteracting -occupancies are reproduced with decimal places as a benchmark.
5.3.1 Physical results
The results of short periods are not significantly different from the quenched SIAM with no oscillation of -orbital energy. So we plot both Figs. 5.8 and 5.9 in the long period regime . Both the energy-ordered and quasi-energy-ordered algorithms as discussed in §5.2 give the same physical results. The Hubbard suppresses the double occupancy of the -orbital for both the quenched and driven SIAMs. In Fig. 5.8, the dashed grey lines indicate the level of double occupancy in a noninteracting SIAM (estimated from the of the quenched ). The interacting double occupancy is appreciably lower than when the driving amplitude is small. For period , both (red line in Fig. 5.8 (a)) and the double occupancy (purple line) oscillate in sinusoidal waves, even though the driving signal is a square wave. When the period increases to , the wave forms approach a relaxed oscillation (Fig. 5.8 (b)). The overshoots in every period disappear in a noninteracting simulation (, not plotted), which produces simple monotonic decays to the square wave levels. So the overshoots in the waveform of the -orbital occupancy are an interaction-induced effect.
When the driving amplitude is increased, the wave form of distorts, and the relaxation to steady-state oscillation slows down, as is shown in Fig. 5.9. Also, there is an increase of the average double occupancy . At , the double occupancy in its oscillation steady state is above almost the entire period. A possible explanation might be that the oscillating -orbital energy is like a phonon mode that induces an effective intra--orbital attraction, which becomes greater than when the oscillation amplitude is big enough (, at which ). Whether this attractive interaction can lead to superconductivity is interesting for further studies.
5.3.2 Complexity results
Obtaining results in Fig. 5.9 (a) at medium to large driving amplitudes was not easy, because in the regime, the linear growth of maximum entanglement entropy makes the maximum bond dimensions in the MPSs increase exponentially with the number of periods simulated. We used some extrapolation techniques to estimate the steady-state quantities in Fig. 5.9 (b), especially for where the relaxation is slow. In this section we mainly check whether this linear entropy growth (exponential difficulty) can be helped by reordering the bath orbitals in the MPSs in quasi-energy order.
We find that even though the entropy growth in the noninteracting SIAM changes from linear to logarithmic by quasi-energy ordering the bath orbitals, as is shown in §5.2, the entropy growth for the interacting SIAM is slightly faster than logarithmic. We increase the number of bath orbitals to to reach , and then make a comparison of the energy-ordered and quasi-energy-ordered simulations in Fig. 5.10 under , . As is shown in Fig. 5.10 (b), the quasi-energy-ordered 4-MPS simulation is slower than the energy-ordered simulation in the short run. The short-term growth of entropy, e.g. in the first few periods, is faster if the energies of the bath orbitals are not ordered. In the long run, the quasi-energy ordering is more favorable. The entropy growth only slightly curves up in the v.s. plot. The long-term growth rate of entropy and v.s. in Fig. 5.10 (b) are clearly reduced. The hardness in the regime is beyond polynomial time using either method, but is significantly reduced using quasi-energy ordering.
Figure 5.11 shows the crossing time of the maximum entanglement entropies of the energy ordered and quasi-energy ordered simulations. In a wide range of driving periods the crossing time of the entropies in Fig. 5.10 exists and is minimum at intermediate driving periods at which the linear growth rate of entropy of the energy ordered method is fastest. After the entropies cross, the quasi-energy-ordered method still needs to overcome two more short-term drawbacks: a) its maximum entanglement entropy being more widespread than the energy-ordered method with maximum entanglement entropy concentrated near the Fermi level, and b) the bigger entropy at short times, before the actual CPU-times cross.
5.4 Summary and conclusion
We have generalized the 4-MPS method in Chap. 4 to time-dependent Hamiltonians to study periodically driven SIAMs. We analyzed the computational complexity in the short period and long period regimes for both the noninteracting () and interacting () models. The model behavior in the regime is not significantly different from the quenched model. This is the regime in which the Floquet-Magnus expansion converges. Both the interacting and noninteracting models are as easy to simulate as the quenched models (polynomial time). In the regime, the entropy grows linearly in the energy-ordered algorithm, which is therefore exponentially hard to reach long times, i.e. many periods. Using quasi-energy ordering reduces the entropy growth of the noninteracting model from linear to logarithmic with a coefficient that grows unboundedly with the driving period (proportional to ). For the interacting model, it also reduces the linear growth rate of entropy and the exponential hardness of the problem in the long run. But there is a tradeoff between the short-term and long-term computational costs, as is revealed in Fig. 5.11. The quasi-energy-ordering algorithm is most favored when the driving period is greater than the critical period by about a factor of , which gives the quickest linear growth rate of entanglement entropy and the earliest crossing point.
Conclusion
In this thesis, we conduct a focused study of strongly correlated systems with localized electron orbitals. We have studied two real materials (LuNiO3 and VO2) and one model system, i.e. the Anderson impurity model, using two main theories (DFT+ and DMFT) and other techniques in the appendices. The electron-electron interactions on the localized orbitals are typically strong compared with other delocalized orbitals. The on-site Coulomb interactions parameterized by Hubbard and Hund’s coupling are included on the Hartree-Fock level in DFT+ to obtain a soft band structure that depends on the orbital occupancies of the electrons. This gives rise to competitions between multiple orbitals and complex phase diagrams as revealed in the equilibrium phase transitions of LuNiO3. In the pump-probe experiment of VO2, the soft band physics leads to the collapse of the interaction-maintained energy gap in Mott insulators as an enough number of electron-hole pairs are created via photo-excitation to induce an insulator-to-metal transition.
While a Hartree-Fock-level treatment of the interactions of the localized orbitals in Chaps. 2 and 3 is conceptually simple to understand and computationally cheap to implement, this uncontrolled approximation can be quantitatively inaccurate and miss qualitative features of the interacting system such as the quasi-particle lifetime and other beyond-band-theory effects. The dynamical mean-field theory (DMFT) attempts to keep the on-site part of the interactions on the localized orbitals exactly as impurity orbitals and only reduce the delocalized orbitals into an effective noninteracting bath, instead of reducing the whole system into an effective noninteracting soft band structure as in DFT+.
Our research works in Chaps. 4 and 5 are motivated by DMFT to build a real-time impurity solver using the density matrix renormalization group (DMRG) technique. We have built a matrix product state (MPS) based impurity solver and have achieved some preliminary complexity results of both the quenched and the driven single-impurity Anderson model. The main intuition is that for the quenched model, one can fully utilize the fact that the bath is noninteracting to control the entropy growth using energy separation, so that different regions of the MPS corresponds to different energy ranges and the noninteracting bath electrons would choose to localize themselves around the energy range they belong to. When the impurity -orbital energy is periodically driven, the impurity orbital functions as an electron elevator that transports bath electrons between different energy levels. This leads to linear entropy growth for energy-ordered bath MPSs. Using quasi-energy-ordered bath MPSs reduces the long-term growth of entanglement entropy but increases the entropy growth in the short run. The entropy crossing point of the two algorithms can come very late. How to efficiently simulate time-dependent impurity Hamiltonians remains an interesting open question. Future work on the efficient simulation of driven systems should focus on analytic properties of the Floquet Hamiltonian of the driven interacting model, so that one can come with a better guess of the “good basis” that would limit the growth of entanglement entropy among the bath orbitals.
Appendix A Maximally localized Wannier functions
In this appendix, we give a mathematical derivation of the formulas for the center and spread of the multi-band Wannier functions in -space. From Eqs. (1.12)–(1.13), we can express the Wannier center of the th Wannier function in terms of the Bloch wave functions in -space via
[TABLE]
where the gauge freedom goes into the decoupled cell-periodic functions
[TABLE]
Since the Bloch waves and are -periodic, we have
[TABLE]
This proves the relation and therefore from Eq. (A.1), we have
[TABLE]
From the cell-periodicity of , we can break the integral into an integration within one unit cell followed by a sum over the unit cells. We have
[TABLE]
The sum over then leads to a k-selection rule that collapses the integral. Since both and are restricted to the first Brillouin zone, we have and thus . So the Wannier center is simply given by
[TABLE]
To calculate the spread we also need the second-order moment
[TABLE]
Even though Eq. (A.3) does not directly generalize to higher-order moments, we can use , take the complex conjugate of Eq. (A.3) to replace one by and the other by to obtain
[TABLE]
following the same procedure from Eq. (A.4) to Eq. (A.7). The metric in need not be identity in a nonorthogonal Bravais lattice. One may in fact put in any metric and Eq. (A.10) gives the corresponding second-order moment so long as the same metric is put in . Once we have Eqs. (A.7) & (A.10), we can calculate the spread of the th Wannier function and numerically minimize the total spread by tuning the unitary matrices in Eq. (1.12).
Appendix B Rotationally invariant on-site interaction
In this appendix, we give a detailed derivation of the parameterization of rotationally invariant two-body interaction tensor in terms of the invariant radial integrals . Then we prove the sum rule of the isotropic and anisotropic integrals.
Parameterization of by spherical symmetry
We require the addition theorem of Legendre polynomials
[TABLE]
Based on the assumptions of Eqs. (1.23)–(1.25), we have from Eq. (B.1) that
[TABLE]
with the radial integrals defined in Eq. (1.26). We have used in the spherical coordinates to separate the radial and angular integrals. We may next use and its relation to the Wigner -symbols
[TABLE]
to obtain the interaction matrix elements
[TABLE]
Because of the selection rules of the Wigner -symbols, the summation of is truncated to only even numbers from [math] to , which can be rewritten as with the new variable summing from [math] to as follows:
[TABLE]
The interaction tensor is parameterized by the radial integrals as linear coefficients of the universal Wigner -symbols.
Expressing and in terms of
From the definition of Hubbard in Eq. (1.22) and Eq. (B.5), we have
[TABLE]
where we have used the selection rule of Wigner -symbols to pick out the term. Then we use the sum rule
[TABLE]
Only the term survives and we obtain
[TABLE]
Next we calculate the sum
[TABLE]
The selection rule requires , which may or may not be reached within . But the sign in Eq. (B.5) is canceled. We can use symmetries of the -symbols to show that the last two symbols are in fact equal. We have
[TABLE]
which then leads to
[TABLE]
We now use the orthonormality relation of Wigner -symbols
[TABLE]
to further reduce Eq. (B.9) into
[TABLE]
The result is independent of as a consequence of the rotational symmetry. Finally, from the definition of Hund’s coupling in Eq. (1.22) and Eqs. (B.8) and (B.11), we have
[TABLE]
Eqs. (B.8) and (B.12) constitute the results of Eq. (1.28). The Hubbard and Hund’s coupling correspond to the isotropic and anisotropic parts of the interaction, respectively.
Appendix C Hybridization function of a fermionic bath
In this appendix, we give a derivation of the hybridization function of a noninteracting fermionic bath. Consider a situation as follows:
[TABLE]
The system can be strongly interacting and correlated, but the bath it is coupled to via the one-body hopping terms in is noninteracting. From Eq. (1.32), we have
[TABLE]
Since the bath conserves particle number, we have , otherwise . Also, because the bath is noninteracting, we can use Wick’s theorem to factorize into a product of one-particle (two operator) Green’s functions. There are totally contractions corresponding to permutations of the summation indices relative to to yield equal contributions. Therefore, we have
[TABLE]
Protected by , the integrals and summations give the same factor raised to power . To make this statement more explicit, we define the bath Green’s function and hybridization function to simplify the expression. The bath Green’s function is given by
[TABLE]
and the hybridization function is defined as
[TABLE]
The eigenstate Green’s function of the bath satisfies the selection rule . The hybridization function is like a superposition-state Green’s function, which contains contributions from various modes. In terms of , we have
[TABLE]
Therefore, the effective action due to a noninteracting fermion bath is
[TABLE]
Eqs. (C.5)–(C.6) will be used for calculating the hybridization function of the bath.
Appendix D Group theory analysis of the NiO6 array
In this appendix, we give a full-length group-theoretical analysis of the energy function of a corner-shared 3D array of NiO6 octahedra without tilts. To find out the symmetry-determined form of the Landau energy as a function of modes , , , , and , we need to extend our configuration space to a minimal group-invariant subspace of dimensions
[TABLE]
This is because the Jahn-Teller distortion along the direction can be rotated to and directions by to give us the mode. Similarly, rotating to and directions gives us , and can be rotated to and . The Landau energy as a function of the modes will have to be invariant under the permutations of the , , and indices due to and the translations along , , and as well. A translation along by one nearest-neighbor Ni-Ni distance, for example, will leave all modes unchanged and will let all modes change sign. Translations in all three directions can generate, in total, ways of sign change. The Landau function will therefore have to be invariant under symmetry operations which include plus translations.
The algorithm we use for determining the symmetry-allowed form of the energy is based mainly on the rearrangement theorem of group theory. We start with a general Taylor expansion of with respect to the variables in Eq. (D.4) to some required order. The truncated expansion, which is a -variate polynomial, is then transformed by each of the symmetry operations. The average of the transformed polynomials is then guaranteed to be invariant under all symmetries according to the rearrangement theorem. Once we find the symmetry-determined function of the modes, we project back to the modes we previously started with by setting the other modes , , , to zero. The general form of is then given by
[TABLE]
The functions are Taylor expandable and have the forms
[TABLE]
where and other functions and all lowercase functions that appear in Eqs. (D.6)–(D.8) are arbitrary Taylor-expandable functions. Equation (D.5) can be thought of as some advanced version of Eq. (2.3) for a single NiO6 octahedron. We used mathematica to expand out all polynomial terms, then implemented the rearrangement projection to sift out symmetry-allowed terms, and aggregated them into the Taylor expansions of the arbitrary functions , , , etc.
Finally, we approximate the uniform distortion modes and as smooth functions of the lattice constant , to simplify the energy function into
[TABLE]
with only three order parameters , , and left. The uniform modes and are treated as control parameters smoothly determined by the lattice constant and disappear from the energy function. This approximation can be justified by the calculated structures in Fig. 2.4 to see that the jumps in and at the transitions are much smaller than those in , , and .
Appendix E K-averaged quantum Boltzmann equation
In this Appendix, we give a detailed derivation of the k-averaged quantum Boltzmann equation (QBE) in Eq. (3.15) of the main text from the standard QBE in Eq. (3.9). As one can see, when the assumption in Eq. (3.11) is satisfied, the number of degrees of freedom of the system is greatly reduced and Eq. (3.9) becomes
[TABLE]
We may insert resolutions of unity
[TABLE]
for on the right-hand side of Eq. (E.1) to get
[TABLE]
Multiplying by on both sides of Eq. (E.3) and summing over in the first Brillouin zone, the left-hand side becomes
[TABLE]
using notations defined in Eqs. (3.12)–(3.14). The right-hand side of Eq. (E.3) becomes
[TABLE]
Up to this point, the treatment has been exact. Here comes the approximation: the matrix element modulus squared
[TABLE]
is replaced by the k-averaged quantity defined in Eq. (3.10), which only depends on the band indices that typically carry orbital information. This is assuming that the k-points are randomized by the scattering processes, so they can be eliminated from the dynamical variables of the distributions of occupancies.
The randomization of k-points lets us have
[TABLE]
Since is independent of , it can be taken out of the k-sums, which then give us the product of four densities of states. Then using notations in Eqs. (3.13) and (3.14), we have
[TABLE]
which reproduces Eq. (3.15). The main assumptions are the slow manifold assumption in Eq. (3.11), which reduces the number of dynamical degrees of freedom, and the local interaction and random band approximations, which justify the k-averaging of the rate constants.
In the actual implementation of Eq. (3.10) to obtain the k-averaged rate constants , it is more convenient to first randomly generate matrix elements that satisfy both momentum and energy conservation, take the sample average, and then multiply the result by a correction factor , where the integration measure
[TABLE]
does not consider momentum conservation, while
[TABLE]
does. But it turns out that according to our actual Monte Carlo data for VO2 (see Fig. E.1). Therefore, the correction factor is insignificant. This result also partly justifies the random band approximation proposed in the main text: the fact that does not make very different from the case that equals any other value, so the conservation of momentum does not make a big difference in the integration measure, i.e., .
Appendix F Floquet theory of a small-amplitude square wave
In this appendix, we obtain the Floquet Hamiltonian of a system driven by a small-amplitude square wave. In general, the Floquet Hamiltonian of a periodically driven system with period is defined by
[TABLE]
where is the time-ordering symbol. For small amplitudes we have . We can take the derivative with respect to at to obtain
[TABLE]
Let us define an expansion for the Floquet Hamiltonian
[TABLE]
Then we have following the same derivation as Eq. (F.2) that
[TABLE]
Comparing Eqs. (F.2) and (F.4), we have from the first-order terms of that
[TABLE]
where we have multiplied on both sides by from the left. Then we use the nested commutator expansion
[TABLE]
where is the adjoint representation of , and is the -fold nested commutator of with . Using this formula on both sides of Eq. (F.5), and from the square wave model
[TABLE]
we have
[TABLE]
or in functional form
[TABLE]
All functions of are defined using their Taylor expansions in Eq. (F.8). We now apply the inverse of the function of on the right-hand side to both sides and after some algebra obtain
[TABLE]
In the eigenbasis of , the matrix elements of and are related by
[TABLE]
where and are eigenstates of with eigen-energies and . Some matrix elements of can be singular when , assuming has a continuous spectral range that is symmetric about [math].
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] P. Hohenberg and W. Kohn. Inhomogeneous electron gas. Phys. Rev. , 136:B 864–B 871, Nov 1964.
- 2[2] W. Kohn and L. J. Sham. Self-consistent equations including exchange and correlation effects. Phys. Rev. , 140:A 1133–A 1138, Nov 1965.
- 3[3] Mel Levy. Electron densities in search of hamiltonians. Phys. Rev. A , 26:1200–1208, Sep 1982.
- 4[4] Elliott H. Lieb. Density functionals for coulomb systems. International Journal of Chemistry , 24:243–277, 1983.
- 5[5] Arash A. Mostofi, Jonathan R. Yates, Giovanni Pizzi, Young-Su Lee, Ivo Souza, David Vanderbilt, and Nicola Marzari. An updated version of wannier 90: A tool for obtaining maximally-localised wannier functions. Computer Physics Communications , 185(8):2309 – 2310, 2014.
- 6[6] M. van Schilfgaarde, Takao Kotani, and S. Faleev. Quasiparticle self-consistent g w 𝑔 𝑤 gw theory. Phys. Rev. Lett. , 96:226402, Jun 2006.
- 7[7] M. Springer and F. Aryasetiawan. Frequency-dependent screened interaction in ni within the random-phase approximation. Phys. Rev. B , 57:4364–4368, Feb 1998.
- 8[8] Takao Kotani. Ab initio random-phase-approximation calculation of the frequency-dependent effective interaction between 3d electrons: Ni, fe, and mno. Journal of Physics: Condensed Matter , 12(11):2413, 2000.