Fock space embedding theory for strongly correlated topological phases

A many-body wave function can be factorized in Fock space into a marginal amplitude describing a set of strongly correlated orbitals and a conditional amplitude for the remaining weakly correlated part. The marginal amplitude is the solution of a Schr\"odinger equation with an effective Hamiltonian that can be viewed as embedding the marginal wave function in the environment of weakly correlated electrons. Here, the complementary equation for the conditional amplitude is replaced by a generalized Kohn-Sham equation, for which an orbital-dependent functional approximation is shown to reproduce the topological phase diagram of a multiband Hubbard model as a function of crystal field and Hubbard parameters. The roles of band filling and interband fluctuations are elucidated.

A many-body wave function can be factorized in Fock space into a marginal amplitude describing a set of strongly correlated orbitals and a conditional amplitude for the remaining weakly correlated part. The marginal amplitude is the solution of a Schrödinger equation with an effective Hamiltonian that can be viewed as embedding the marginal wave function in the environment of weakly correlated electrons. Here, the complementary equation for the conditional amplitude is replaced by a generalized Kohn-Sham equation, for which an orbital-dependent functional approximation is shown to reproduce the topological phase diagram of a multiband Hubbard model as a function of crystal field and Hubbard parameters. The roles of band filling and interband fluctuations are elucidated.
First-principles calculations of topological invariants usually rely on the Kohn-Sham band structure. This is problematic for correlated materials: the topological phase inferred from a mean-field band structure need not coincide with the actual topological phase determined from the correlated many-body wave function. Although one can argue that a topological invariant cannot change as interactions are turned on adiabatically while maintaining an energy gap and the relevant symmetries, true strongly correlated topological phases could be characterized as precisely those phases that are not adiabatically connected to the noninteracting Kohn-Sham ground state.
Embedding theories have been successful in describing electronic correlations beyond standard density functionals in extended systems [1][2][3][4][5][6][7][8][9][10][11][12][13]. In most embedding schemes, a real-space fragment such as an impurity site or a small set of atomic or molecular orbitals, is embedded in its surroundings. The use of a real-space fragment, typically with only local or short-range interactions, inherently limits the nonlocality and hence quasimomentum dependence that can be described. Since topological invariants depend on the global k-dependence of the state, either through the twisting of Bloch functions in the Brillouin zone [? ] or the behavior of the many-body wave function under twisted boundary conditions [14][15][16], it is natural to ask whether alternative embedding theories might be better suited to capturing momentumdependent correlations.
This Letter proposes a novel embedding theory rooted in the exact factorization (EF) methodology [17][18][19], a scheme for splitting the many-body wave function into marginal and conditional probability amplitudes describing different degrees of freedom. It has been applied to electrons and nuclei [18][19][20][21][22], fast and slow electrons [23], electrons and photons [24,25], and electrons and phonons [26]. We turn to the problem of strong electron-electron correlation and use an extension of the exact factorization formalism to Fock space [27] to develop a novel embed-ding theory. Here, the marginal amplitude describes the strongly correlated degrees of freedom embedded in the remaining weakly correlated degrees of freedom.
For the purpose of calculating topological invariants, the key advantage of an EF-based embedding formalism lies in the ability of the explicitly correlated marginal wave function to capture the k-dependent phase information of the strongly correlated electrons, which is partly lost in approaches based on Green's functions or reduced density matrices. The contribution of the remaining weakly correlated electron bands can be adequately described through mean-field Bloch functions. Thus, one can go beyond density functional theory, while avoiding empirical models and the infeasibility of including all degrees of freedom in a many-body calculation of the solid.
To apply the exact factorization formalism to a manyelectron wave function, we start by writing it in Fock space as a superposition of products of many-body configurations where |0 are constructed from orbitals belonging to mutually orthogonal sets S and D of weakly and strongly correlated orbitals; S = s 1 s 2 . . . s N S is a string of indices labeling the orbitals and similarly for D. The |S and |D factors may have varying particle numbers in each term of Eq. (1), but N S + N D = N is fixed. Following the EF procedure [27], the marginal amplitude is defined to be and an arbitrary phase Θ D . The conditional factor 2 then satisfies the partial normalization condition We thus arrive at the factorization c SD = χ D Φ S|D .
A criterion is needed to partition the complete set of single-particle orbitals into weakly and strongly correlated sets S and D. While different strategies are possible, here we perform the separation through a criterion involving the natural occupation number bands, i.e. the bands formed by the k-dependent eigenvalues of the onebody reduced density matrix in the Brillouin zone of the crystal. Strongly correlated orbitals are defined to be those belonging to a band whose occupation numbers satisfy f lower ≤ f nk ≤ f upper with judiciously chosen f lower and f upper , while the rest are called weakly correlated. While our formalism can be applied with any choice of S and D, it will become computationally prohibitive if too many bands are included in D. We have in mind situations where correlations are concentrated in relatively few bands with only weak residual correlations in S, so that the occupation numbers of the latter are very close to 0 and 1 and can be accurately described by standard density functional approximations. Natural occupation numbers in extended systems have only been reported for one-band systems, namely the homogeneous electron gas [28][29][30] and the one-dimensional Hubbard model [31,32]. Recent calculations of a multiband Hubbard model [33] used an unfolding procedure [34] with twisted boundary conditions to derive a continuous band structure from the discrete set of natural occupation numbers and orbitals obtained from exact diagonalization. It was found that when there is a disparity in the strength of interactions in bands of different orbital character one can have simultaneously a set of strongly correlated bands satisfying f lower ≤ f nk ≤ f upper and another set with occupation numbers very close to 0 and 1. This situation might arise, for instance, in transition metal-bearing oxides, where bands with predominantly transition metal d orbital character experience a stronger Hubbard repulsion. In general, it might be necessary carry out multiple self-consistent calculations with different partitions to find the variational minimum.
Given the above choice of partition, our theory embeds a set of natural Bloch orbital bands in an environment made up of all remaining bands. This is the crux of our approach and distinguishes it from all other embedding theories, most of which rely on a real-space partition. Preserving translational symmetry by keeping entire bands intact in the correlated subspace is the key to reliably calculating topological invariants and, in turn, topological phase diagrams, which depend on nonlocal correlations beyond those confined within a realspace fragment. To see the effect of nonlocal correlations, in Fig. 1 we compare the phase diagram of the half-filled ionic Hubbard modelĤ = iσ [−tc † iσ c i+1σ + H.c. + (−1) i ∆c † iσ c iσ ] + i Un i↑ni↓ calculated by several methods: mean-field theory (MF), a renormalization group (RG) method applied to the bosonized Hamiltonian [35], density matrix embedding theory (DMET) [8], and exact diagonalization (ED) extrapolated to the thermodynamic limit using data from periodic 8-, 10and 12-site models following the approach in Ref. 36. In our DMET implementation with a 2-site embedding fragment and spin-symmetry preserving interacting bath [37], the phase boundary calculated from the polarization of the band electrons is overestimated, i.e. the band insulator to Mott insulator transition occurs at a higher U than in the exact result. This demonstrates that nonlocal correlations can be important even when a model contains only local (Hubbard) interactions. where will be referred to as the embedding Hamiltonian. The factor H SD,S D = SD|Ĥ|S D , which can be broken down into one-body and two-body contributions, induces charge fluctuations between S and D subspaces. Despite its simple appearance, the Hamiltonian in Eq. (6) is actually quite unusual in that it couples many-body config-urations with vastly different particle number, spin, etc. The next step would be to derive the equation for the conditional factor Φ S|D , which is needed to explicitly construct H DD . However, solving the coupled equations for χ and Φ S|D in their full complexity would be tantamount to solving the original Schrödinger equation. Thus, we seek an alternative path that will determine H DD as well as the strongly and weakly correlated orbitals selfconsistently. To obtain a scheme that can be applied to real materials, we couple Eq. (5) to the following generalized Kohn-Sham (GKS) equation: +v ext +v hxc +ŵ hxc |ψ nk = nk |ψ nk , wherev hxc =v hxc [n, ψ dk , χ] denotes a scalar multiplicative potential andŵ hxc [n, ψ dk , χ] is a nonlocal operator acting only in the subspace of strongly correlated natural orbitals ψ dk ∈ D. Bothv hxc andŵ hxc are functionals of the electronic density n(r), ψ dk (r), and χ D . The Hamiltonian matrix elements H DD are similarly functionals of n(r) by virtue of the Hohenberg-Kohn theorem [38]. The GKS equation can be derived by making the energy stationary with respect to variations of n(r) and ψ dk [33]. The density is determined in a nonstandard way as with fractional strongly correlated occupation numbers determined from the marginal factor according to f dk = χ|c † dkσ c dkσ |χ . The weakly correlated orbitals ψ nk ∈ S have occupation numbers 0 and 1, with the possible exception of orbitals with energies equal to the chemical potential. The coupling to the marginal equation enters implicitly through n(r) as well as the χ-dependence of v hxc (r) = δE hxc /δn(r)| φ dk and the matrix elements where the energy has been partitioned as E = T s,e + v(r)n(r)dr+E hxc with an ensemble kinetic energy functional T s,e defined through the constrained search [39] over ensembles of Slater determinants ρ s .
So far no approximations have been made. Solving Eqs. (5) and (7) self-consistently would yield the exact n(r), ψ dk (r), and χ D . To have a practical scheme, we need to specify functional approximations for v hxc (r), w hxc and H DD . For this purpose, we introduce the following approximation, which we call the Aufbau approximation, to define the conditional amplitude Φ S|D . Namely, for each configuration |D we define the Slater determinant |S Aufbau built from the N S = N − N D lowest energy weakly correlated orbitals subject to the con-ditions that (i) theŜ z eigenvalues satisfy M S + M D = M and (ii) the quasimomentum eigenvalues satisfy K D + K S = K, where M and K are the quantum numbers of |Ψ (additional symmetries could also be imposed at this stage). Since the multi-index D uniquely determines S if the weakly correlated orbitals are nondegenerate, which, for simplicity, we assume, there exists a function S Aufbau (D). Thus, we define Since |S Aufbau is a Slater determinant of KS-like orbitals, Eq. (10) allows us to construct H DD as well as the total energy as implicit functionals of n(r). Thus, for any given choice of approximate local potentialv hxc , we have specified an approximation that can be applied to real materials without further functional development.
Although solving the many-body Schrödinger equation in Eq. (5) remains a challenging task, especially in higher dimensions, it is worth emphasizing that the EF method has simplified the original problem to a degree that established many-body techniques can be applied, while retaining the coupling to all remaining electronic degrees of freedom of the solid.
Before pursuing calculations of topological phases in real materials, it is desirable to test the theory in a case where it can be compared with benchmark calculations. To this end, we calculate the topological phase diagram of a multiband ionic Hubbard model, comprising two s bands and two d bands: To study the behavior of bands with vastly disparate interactions, we take the s electrons to be noninteracting. The hopping amplitudes depend on the sublattice displacement ξ according to t s,ii+1 = t s1 = t 0 − 2g s ξ for i = odd and t s,ii+1 = t s2 = t 0 + 2g s ξ for i = even and similarly for t d,ii+1 . The onsite potentials are staggered, i.e. s,i = (−1) i ∆ s and similarly for d,i . The s and d bands are coupled by a hopping term This model has a band insulator to Mott insulator transition at a critical value of the Hubbard parameter U , similar to the single band model [36,[40][41][42][43][44][45][46][47][48][49][50][51][52][53].
We also include a crystal field term∆ sd = ∆ sd iσ (n d,iσ − n s,iσ ) to break particle-hole symmetry.
By varying ∆ sd , we control the filling of the d band and study the effect of band-filling on the quantum phase transition. We are effectively using the weakly correlated "spectator" band to dope the strongly correlated band (reminiscent of carrier doping in cuprate superconductors [54]). Previous studies involving variable band-filling in multiband Hubbard models, e.g. in connection with the orbital-selective Mott transition [55][56][57][58][59][60] and strongly correlated superconductivity [61], have tended to treat higher symmetry scenarios with the same value of intraband U for all bands and additional interband and exchange interactions.
We start by discussing the model from the mean-field perspective. The interaction-driven transition from paramagnetic to antiferromagnetic phase upon increasing U is depicted in Fig. 2a-c. The closing of the d-orbital up-spin gap (Fig. 1b) as the critical value of U is approached from above signals the transition to the paramagnetic phase. A different scenario is found for the crystal field-driven transition ( Fig. 2d-f) with varying ∆ sd , where the gapclosing defining the critical point occurs between bands of different orbital character (see Fig. 1e). Figure 3 shows the topological phase diagram obtained in the Aufbau approximation as a function of U and ∆ sd . The Aufbau solution (solid black curve) displays transitions from a band insulator (BI) to a Mott insulator (MI) as either U or ∆ sd is increased. The phase boundaries are determined from jumps of π in the marginal geometric phase γ χ = 2π 0 i χ|∂ α χ dα [62]. The phase transitions signal a discontinuous change in the macroscopic polarization P = −(e/2π)γ, which is a topological invariant quantized to P = 0 or e 2 mod e by the parity symmetry of the model. Similar behavior is well known in singleorbital ionic Hubbard models in one and two dimensions [36,41,43]. We use a small symmetry breaking dimerization ξ = 5 × 10 −5 , and the Born-von Kármán cells used in our calculations are too small to see an intermediate bond-ordered phase [44,47,52]. The Aufbau solution agrees well with the one obtained by numerical exact diagonalization (red dots), demonstrating the viability of the approximation.
The BI-MI transition is reflected in the paramagnetic to antiferromagnetic mean-field phase boundaries (dashed gray curves), which roughly follow the BI to MI transition of the correlated solution. However, a second symmetry-restoring transition is reached when ∆ sd is further increased. For U = 0, the crystal fielddriven transition occurs at exactly ∆ sd = 1 2 (∆ s + ∆ d ) = 1.80 eV, as correctly reproduced in the mean-field solution. The intercept deviates slightly in the Aufbau approximation, which does not become exact in the limit U → 0 because it does not capture all t sd -induced interband charge fluctuations. At the symmetry-breaking transition, the geometric phase of the d-orbital valence band associated with one spin, e.g. the down-spin, γ d↓ = π/a −π/a i u dk↓ |∂ k u dk↓ dk, jumps from π to 0 as ∆ sd (or U ) is increased. At the symmetry-restoring transition, the geometric phase of the opposite spin, γ d↑ , also jumps from π to 0. No such second transition was observed in either the exact diagonalization or Aufbau solutions in the investigated range.
The change in the topological invariant from the BI to MI phase coincides with a change in the topology of the natural occupation number bands as shown in Figure 4. The strongly correlated occupation number bands develop a crossing at the zone boundary, which implies a natural occupation number "band inversion." We have described this phenomenon in the Rice-Mele-Hubbard [63], and it is similar to the purity-gap closing studied in the quench dynamics of ultracold atoms [64]. With respect to other descriptors of the transition, the natural occupation number bands have the virtues that they build in the crystal symmetries through their symmetry properties and those of the natural Bloch orbitals and can be straightforwardly extended to probe real-time dynamics.
The band-filling controlled by a chemical potential. While in the former case the phase boundaries can be detected by discontinuous π-jumps, the mean-field geometric phase of the latter is blind to the paramagnetic-antiferromagnetic transition [62]. This underscores the usefulness of an EF approach built on pure states for the detection of topological phase transitions. We expect our approach to also yield accurate energies and local observables.
In summary, a novel embedding theory based on a Fock space factorization has been found to reproduce the topological phase diagram of a strongly correlated multiband system. Calculations employing the proposed Aufbau approximation can be expected to aid the ongoing search for novel strongly correlated materials.
Note added. We recently became aware of related work by Lacombe and Maitra that also develops an exact factorization-based embedding method [65]. FIG. 4. Inversion of natural occupation number bands as ∆ sd is increased through the phase boundary from ∆ sd = 1.0 to ∆ sd = 1.1 for ts = t d = 3.5, ∆s = 1.6, ∆ d = 2.0, ξ = 5×10 −5 , t sd = 0.8 and U = 6.6 (all in eV). The zone boundary (k = π) is positioned at the center to make the crossing visible.

MANY-BODY CALCULATIONS
Both the Aufbau and exact diagonalization solutions reported in Fig. 2 of the main text were carried out for a Born-von Kármán cell with length L = 3a. With 4 orbitals (2 atoms and 2 orbitals/atom) in each primitive cell, this corresponds to 12 electrons and 12 orbitals. Brillouin zone integrations and twist-averaging were performed with between 96 and 384 k points.
The solution in the Aufbau approximation was obtained by minimizing the energy with respect to χ and all orbitals (weakly and strongly correlated). Since the s band in our model is noninteracting (see Eq. (13) of the main text), we can avoid solving the GKS equation. This minimization procedure results in fractional occupation numbers for the weakly-correlated orbitals, which can therefore no longer be identified with a noninteracting KS system with a local potential. Nevertheless, since their occupation numbers are all very close to 0 or 1, namely within 1.5 × 10 −5 for t s = 3.5, t d = 3.0, ∆ s = 0.50, ∆ d = 0.55, ∆ sd = 0, ξ = 0.0080, t sd = 0.8 and U = 7.2 (all in eV), the result is essentially identical to the result of solving the GKS equations with integeroccupied weakly-correlated orbitals.
The exact diagonalization solution was obtained by diagonalizing the Hamiltonian in a restricted Hilbert space consisting of all states of the following types: where |F S is the s electron Fermi sea, N D is the number of electrons in configuration |D , K D is the quasimomentum eigenvalue, and M D is theŜ z eigenvalue.

MEAN-FIELD PHASE DIAGRAM OF MULTIBAND VS. SINGLE BAND MODELS
The crystal field parameter ∆ sd allows us to control the relative occupation of s and d orbitals. The ground state, |χ , of the embedding Hamiltonian in Eq. (5) of the main text superposes different charge states of the strongly-correlated subspace into a single coherent pure state. This is distinctly different from reduced density matrix methods, where the strongly-correlated subspace is described by a partially incoherent mixed state and its mean occupancy is controlled by a chemical potential. To explore whether this difference has any consequences for the calculation of topological phase diagrams, in Fig. 1 we compare the mean-field phase diagrams of (i) the multiband Hamiltonian defined in Eqs. (11) and (12) of the main text and (ii) the corresponding single-band Hamiltonian defined byĤ d in Eq. (11). In case (i), the occupancy of the d orbitals is controlled by the parameter ∆ sd ; in case (ii), it is controlled by the chemical potential µ. The significant differences in the phase diagrams of the two cases suggests that doping by a "spectator" band is not equivalent to doping by an external particle reservoir. Therefore, care should be taken when using the chemical potential as a parameter in a topological phase diagram. There is another important distinction between the two cases that relates to their sensitivity to topological properties. In case (i), the paramagnetic-antiferromagnetic phase boundaries is signaled by π-jumps in the up-spin and down-spin mean-field geometric phases, in case (ii) it is not. In case (ii), the phase boundaries have instead been determined by the spin-symmetry-breaking paramagnetic to antiferromagnetic transition.
FIG. 1: Mean-field phase diagram of the multiorbital Hubbard model of the main text (dashed gray curves) as a function of (U, ∆ sd ) and the corresponding single-orbital Hubbard model (blue curves) as a function of (U, µ). Parameters are ts = t d = 3.5, ∆s = 1.6, ∆ d = 2.0, ξ = 5×10 −5 and t sd = 0.8; all in eV.

ENTANGLEMENT
With the help of the Schmidt decomposition arXiv:1909.07933v2 [cond-mat.str-el] 21 Apr 2021 the entanglement between a subsystem and its environment can be quantified by the entanglement entropy where N is the number of nonzero eigenvalues λ i of the subsystem reduced density matrix The Aufbau approximation proposed in the main text delivers a wave function in the Schmidt-decomposed form Every many-body configuration |D is paired with a state |S Aufbau (D) from the environment. The entanglement entropy in the Aufbau approximation is where N D is the number of nonzero coefficients χ D and we have introduced an ordinal index, denoted i, which labels the set of all nonzero χ D coefficients. Since many of the χ D , defined by the solution of the model Hamiltonian in Eq. (5) of the main text, differ significantly from 0 and 1 in the strongly correlated regime-a fact which follows from our assumption that the natural occupation numbers of the strongly correlated bands satisfy f lower ≤ f nk ≤ f upper -our embedding theory in the Aufbau approximation captures a significant amount of entanglement between the weakly and strongly correlated subspaces.

GEOMETRIC PHASE
Similar to many topological invariants [1], the manybody macroscopic polarization can be evaluated in terms of the ground state of a twisted Hamiltonian, i.e. with an artificial magnetic flux α [2]. The geometric phase of the Aufbau ground state is π/a −π/a i u nkσ |∂ k u nkσ dk, (6) where α is an artificial magnetic flux and u nkσ (r) is the cell-periodic part of the Bloch function ψ nkσ (r). The first term is simply the many-body geometric phase [2] of the marginal factor, while the second term is a |χ D | 2weighted sum of single-particle contributions [3]. For clarity, the formula is given for the special case of a onedimensional system.