Maximally localized dynamical quantum embedding for solving many-body correlated systems

Quantum computing opens new avenues for modeling correlated materials, which are notoriously challenging to solve due to the presence of large electronic correlations. Quantum embedding approaches, such as dynamical mean-field theory, provide corrections to first-principles calculations for strongly correlated materials, which are poorly described at lower levels of theory. Such embedding approaches are computationally demanding on classical computing architectures and hence remain restricted to small systems, limiting the scope of their applicability. Hitherto, implementations on quantum computers have been limited by hardware constraints. Here, we derive a compact representation, where the number of quantum states is reduced for a given system while retaining a high level of accuracy. We benchmark our method for archetypal quantum states of matter that emerge due to electronic correlations, such as Kondo and Mott physics, both at equilibrium and for quenched systems. We implement this approach on a quantum emulator, demonstrating a reduction of the required number of qubits. Near-term quantum computers hold many promises but remain limited to a moderate number of qubits. This Article presents a pathway for modeling correlated materials with a reduced number of qubits, bringing quantum computing to materials modeling.


INTRODUCTION
The understanding of materials with strongly correlated electrons is one of the main challenges of modern solid state physics.Triggered by the discovery of high-temperature superconductivity in copper-oxides, the study of doped Mott insulators has grown in the last decades, building on the development of theoretical tools designed to solve accurately models of strongly correlated electrons.
Despite that the exact solution of simple correlated theoretical models in two or three dimensions is lacking, accurate predictions for the properties of strongly correlated solids are obtained by using approximations [1].A central role has been played by the dynamical mean field theory (DMFT) [2], a non perturbative method which allowed for the first complete description of the Mott-Hubbard transition.This method has been extended to a variety of correlated methods and combined with density functional theory, leading to remarkable agreement with the properties of many correlated materials.DMFT more generally falls within the larger group of quantum embedding theories [3,4] (for a review see [5] and references therein), which have been widely successful at describing transition metal and f -elements into both solids and molecular systems.
The central idea of DMFT is to self-consistently map the infinite bulk system onto a so-called Anderson impurity model (AIM) with only a few interacting impurity sites embedded in an infinite non-interacting bath.The latter Anderson impurity models can be solved using high-level many-body methods, with a breadth of approaches for all have their own limitations.Indeed, the recent development of CTQMC has generated a strong activity in the field.For single-site DMFT, CTQMC yields an exact solution of the AIM problem within the statistical error bars in imaginary time.The main limitation of the approach is that the evaluation of real-frequency spectra requires a poorly conditioned analytical continuation, based on the maximum entropy method [6,7] (or some alternative strategy), which strongly limits the possibility to study fine details of the spectra.For multi-orbital or cluster extensions of DMFT, CTQMC suffers from the fermionic sign-problem as long as inter-orbital hybridizations are present, and it is therefore limited to finite temperatures.
Numerical renormalization group (NRG) [8] provides an alternative for real axis calculations and access to Kondo physics, but remains challenging to extend for multi-orbital systems.
Finally, exact diagonalization (ED) solvers are instead based on a finite discretization of the AIM, through the representation of the effective bath in terms of a small number of bath-sites.In practical implementations, the bath size (N b ) is severely limited because of the growth of the Hilbert space: its dimension scales exponentially with the total number of sites N s (bath sites and impurity orbitals).Nonetheless, the use of Lanczos-based algorithms allows to deal with large Hilbert spaces, and the discretization at low temperature [9,10] is fine enough to accurately compute the thermodynamic and static observables.In particular, the finite size effects affect the spectral functions, which are slowly converging to the continuous features of the exact DMFT solution.Notwithstanding reasonably accurate static phase diagrams obtained with the Lanczos solver for three [11,12] and five [13] orbitals system, the limitation in the bath size becomes also particularly relevant for multi-orbital AIM models, which are necessary to realistically describe transition metal oxides, as the impurity's orbitals (five for a d-manyfold) contribute to the enlargement of the Hilbert space.
Scaling the precision of the latter approach on classical computers is not feasible due to the exponential increase of the Hilbert space with number of bath sites.In the view of the recent progresses in the implementation of ED on quantum computers [14], the latter approach has gained interest due to the possibility of scaling linearly on quantum computers.However, due to current hardware limitations (noise and decoherence issues), such algorithm are limited to small number of qubits and short quantum circuits on the currently available noisy intermediate scale quantum (NISQ) computers.Hence it is of great importance to develop an ED based solver that can obtain precise results with a reduced number of bath sites, since the required number of qubits is proportional to this.
In this light, we present here a new quantum embedding methodology to resolve the Anderson impurity model (AIM) for DMFT, based on an extended ED method.Our method provides a maximally localized quantum impurity model, where the non-local self-energy component of the correlation remains minimal, and hence the AIM minimally breaks locality (DMFT is a purely local theory).As reported in this work, this comes at remarkable benefit, as the environment used in the quantum embedding approach is described by propagating correlated electrons (instead of free electrons in DMFT), and hence offers an exponentially increasing numbers of degrees of freedom for the embedding mapping, in contrast for traditional free electron representation where the scaling is linear.This is reminiscent from the representation of correlated electrons by a Green's function embedding approach, where correlations are described by hidden fictitious additional fermionic degrees of freedom [15,16].
This representation has hence the potential to improve dramatically the scope of applicability of the quantum embedding approach, whilst limiting the small number of bath sites.We report that quantum impurity models with as few as 3 bath sites can reproduce both the Kondo regime and the Mott transition, and obtain excellent agreement for dynamical magnetic susceptibilities, poising this approach as a candidate to describe 2-particle excitations such as excitons in correlated systems, such as high-Tc superconductor [17,18].Our approach aligns with recent progresses in quantum computing, where a realistic number of qubits would achieve a fine description of correlations in materials.

FORMALISM
Within DMFT, the effective AIM is subject to a self-consistency condition which relates the Green's function of the impurity model G(iω n ) to the so-called Weiss-field G −1 0 (iω n ), which completely characterises the AIM.For the single band case, the Hamiltonian of the AIM reads: where d † pσ (d pσ ) creates (destroys) a particle with spin σ in the d-orbitals of the uncorrelated bath (p ∈ [1, N b ]) and f † σ (f σ ) creates (destroys) a spin σ particle on the impurity, U is the static Coulomb repulsion on the impurity and V is the tunnelling amplitude between the impurity and the bath.Typically, for a fixed set of the Hamiltonian parameters, the AIM is solved (1) by using a Lanczos algorithm [31] to converge the ground state and excited states [9,12] which contribute to the thermal average.Once the eigenstates are obtained we compute the dynamical and static observables.
In the maximally localized dynamical embedding approach (MLDE) we add an additional general two-electron interaction to the bath sites, and the interaction vertex reads: The U (2) tensor is a fictitious two-body interaction between bath electrons introduced to enlarge the degrees of freedom of the approach.Although the free electron propagator takes a simple polynomial form, e.g.
, a correlated green's function is instead described by a mapping to an infinite one dimensional chain within the Kryvlov space [19].In particular, when the bath is correlated, the Weiss field incorporates the dressed propagator of the bath electrons G (1) d : where the hybridisation to the bath is denoted as ∆ (1) .We emphasize that at this level of theory, the self energy of the impurity Σ f remains nought, and hence we have a fully local theory from the perspective of the impurity.However, as we introduce a local correlation U on the impurity, two effects occur: i) the propagator of the bath G (1) d is dependent on the correlation on the impurity, and hence changes to G 2 d and ii) a non-local part of the self-energy Σ f d between the impurity and the bath emerges: We have now a set of embedding equations, which leads to a generalized Dyson equation , where the δ∆(z) is a source field that stems from the non-local correlations.Indeed, the latter term can be rationalized following a simple argument via the Migdal energy functional which in MLDE reads Part of the correlation energy spills effectively on the bath [20], leading to a correlation leakage term Z.The DMFT equations can be recovered when the leakage potential δ∆(z) and leakage correlation energy are small.As the tuning of the bath propagator allows for a very large set of parameters, these constraints can be successfully enforced via Lagrange parameters in the fit of the DMFT hybridisation, to maximise the locality of the embedding, with a concomitant exponential improvement of the bath discretization errors.It is however known that in high dimension, or when a system is strongly correlated, the electron self-energy is well separable into a local dynamical part and a static non-local contribution [21].In this respect, we reabsorb the embedding potential into a shift of the static part of the MLDE self energy, which insures that the Migdal energy remains exact.

Benchmark
We now turn to a simple test of the MLDE equations for a typical AIM.The bath hybridisation used as a test case is obtained from a typical correlated material [32] and the impurity energy is set to f = −U/2 to stay at half-filling.We perform a benchmark of the MLDE approach with respect to both a continuous time Monte Carlo and Lanczos solver, which both provide in this case the same answer used as a reference.Both the imaginary part (Fig. 2.a) and real part ( Fig. 2.b) obtained by MLDE with as few as N b = 3 bath sites provide a remarkable agreement with the exact solution.We considered both the nearly free electron (NFE) limit (U < 10) and the atomic limit (U > 10), and in both cases the MLDE solution is consistently in agreement with the exact answer.Noteworthy, the strength of correlation and overall physical properties are also well captured by MLDE with N b = 2, whereas the ED solver with the same number of bath sites largely overestimates the strength of correlation (see Fig. Kondo physics We extend further the application of MLDE to a realistic study case of the deposition of a correlated Kondo molecule on a gold surface (see Fig. 4.a).Stable organic radical molecules exhibit a Kondo peak in the low temperature experimental conductance, which is due to the presence of a single unpaired electron in the highest occupied molecular orbital [22][23][24].In the gas phase these molecules are paramagnetic.Due to their low spin-orbit interaction and small hyperfine splitting they are expected to exhibit long spin-coherence times, and therefore have potential as building blocks of molecular spintronics applications [24].When brought in contact to a metal surface the system corresponds to a single impurity Anderson model (SIAM), and has been modelled in the past using CTQMC or NRG as impurity solvers [25,26].To demonstrate the capability of MLDE to describe Kondo physics we choose the 1,3,5-triphenyl-6-oxoverdazyl (TOV) organic radical molecule, which when deposited on a Au substrate has been shown experimentally to exhibit a Kondo temperature of about ≈ 37K [22].Compared to other radical molecules on surfaces this molecule has the advantage to have a well defined contact geometry, where the molecule lies flat on the surface.We use the same simulation setup and parameters as in Ref. [25], so that also the hybridization function of the SIAM is the same.Describing Kondo physics at low temperature is a notoriously difficult problem for quantum impurity solvers.In particular, the collapse of energy scales in the Kondo limit prevents typical Lanczos or ED solvers from capturing the Kondo resonance, as the finite discretization tends to introduce fictitious gaped states at very low temperature.
We performed calculation at T = 5K in the Kondo regime (see Fig. 4.b), the hybridisation remains constant up to the lowest frequency.We note that ∆ 1 (iω n ) = ∆ 2 (iω n ), which confirms that the embedding source field δ∆(iω n ) remains negligible.Furthermore, MLDE fares better with N b = 4 than the best possible fit obtained by the usual discretized Lanczos approach (N b = 12, triangles in Fig. 4.b), as the hybridisation vanishes at small frequency (in the molecule it remains constant).The imaginary part of the self-energy obtained by MLDE shows a Fermi liquid type behaviour at small frequency (see Fig. 4.c), whereas the self-energy obtained by ED shows an artificial Mott singularity, due to the bath discretization.We note that this low temperature was beyond the reach of our CTQMC solver (β = 35000).Below the Kondo temperature, we recover the Fermi liquid behavior of the self-energy (see Fig. 4.d).As the MLDE representation is compact, it opens a large degree of possible manipulation once the AIM is established.In particular, we extended the calculation to the time dynamics of the Kondo molecule with the Keldysh formalism [27] after a magnetic quench, where at time t = 0 the molecule is magnetically polarized along the e z axis, and the external magnetic field released for t > t 0 .The magnetic moment enters in a precession dynamics with frequency equal to the Kondo temperature (see Fig. 4.e).The dynamics can also be resolved below T K , where we observe additional harmonics, reminiscent from the Kondo Zeeman splitting at t = t 0 (see Ref. 28).
Mott transition We have so far focused at simple AIM systems in absence of mean-field corrections to the hybridisation.We now turn to the dynamical mean-field MLDE approach, applied to the Mott transition of the two-dimensional square-lattice Hubbard model (see Fig. 5.a).Using the MLDE as a solver for DMFT, we recover the well known metalinsulator transition (MIT) associated with the charge localization induced by the local Hubbard repulsion U .We recover with a simple MLDE and N b = 3 bath sites the spectral function of the Hubbard model, with the usual features (lower and upper Hubbard bands, quasi-particle peak below the transition U c ≈ 9, Mott gap for U > U c ).We note that interestingly the spectral function obtained by MLDE also shows the satellite peaks at the gap edge for U > U c .The MLDE charge gap ∆ reproduces the known trend ∆ = U − W , which provides a further test of the theory.The benchmark with the converged solution obtained by CTQMC is remarkable (see Fig. 5.b), and the MLDE solution essentially within the error bars of the CTQMC for U/t = 12.Across the Mott transition, the agreement between the Green's functions remains good (see Fig. 5

.c).
Bethe-Salpeter Equation We extended the calculations to the Bethe-Salpeter Equation (BSE) formalism, applied to the MLDE solution.In particular we calculate the local irreducible vertex Γ (see supplementary information), which enables the calculation of the non-local dynamical magnetic susceptibility χ mag (ω) within MLDE.We performed calculations for U/t = 12 at various temperatures, to explore the behavior of magnetic excitations across the Mott gap melting.Within the Mott gap phase at temperature T /t = 0.025 (Fig. 6.a), the Neel fluctuations are local in momentum at Q Neel = (π, π).As expected, the spin fluctuations are largely located at Q Neel in the Mott phase (T = 0.025t).Interestingly, the spin-flucutations survive for temperature as high as the melting temperature (Fig. 6.b).At very large temperature (Fig. 6.c), the fluctuations are fully incoherent and the spectrum uniform across the Brillouin zone.
We note that a known challenge for the BSE approach is the extraction of the local irreducible vertex Γ, which is obtained by calculating the two electron response function G (2) .Such a quantity requires traditionally computationally demanding collection of statistics by CTQMC, and few other alternatives exist.The vertex can be calculated in the Lehman representation (see SI), but requires a sampling over the whole Hilbert space, which is not possible for AIM with more than few bath sites.In MLDE this task is largely simplified as the Hilbert space remains compact.We provide a benchmark of the local vertex Γ with CTQMC for both U/t = 6 and U/t = 12 (see respectively Figs. 6.e and 6.g), in both cases the agreement is quantitative.This agreement is also obtained in the fermionic representation of the local magnetic susceptibility χ iν,iν (see Figs. 6.d and 6.f).

CONCLUSIONS
In conclusion we introduced a novel quantum many-body embedding techniques, the so-called Maximally Localized Dynamical Embedding approach (MLDE), which offers a robust and efficient methodology to describe electronic correlations in quantum materials.This approach generalizes the local dynamical mean field theory to minimally non-local Anderson impurity model, opening up remarkably the resolving power of the discretized impurity model.In particular, when the number of poles described by a canonical DMFT hybridisation scales linearly with the number of quantum impurity sites, the latter increases exponentially in MLDE.This opens new avenues in the realm of quantum computing, where with a realistic number of qubits (typically 10-20) MLDE would allow to describe transition metal systems on a quantum computer, with minor errors induced by the bath discretization.As the MLDE hamiltonian is minimal (typically 4 sites in total), we have shown that time evolutions (Keldysh) and vertex calculations become routinely possible at minimal cost.Finally, the correlated bath sites of MLDE also allows to describe Kondo physics at very low temperature, which is a known limitation of standard exact diagonalization DMFT approaches.
ACKNOWLEDGMENT CW acknowledges insightful and stimulating discussions with Andrew Mitchell.CW was supported by grant EP/R02992X/1 from the UK Engineering and Physical Sciences Research Council (EPSRC).I. R. acknowledges the support of the UK government department for Business, Energy and Industrial Strategy through the UK national quantum technologies programme (InnovateUK ISCF QUANTIFI project).C.L. and F. J. are supported by the EPSRC Centre for Doctoral Training in Cross-Disciplinary Approaches to Non-Equilibrium Systems (CANES, EP/L015854/1).F.J is supported by the Simons Many-Electron Collaboration.This work was performed using resources provided by the ARCHER UK National Supercomputing Service and the Cambridge Service for Data Driven Discovery (CSD3) operated by the University of Cambridge Research Computing Service (www.csd3.cam.ac.uk), provided by Dell EMC and Intel using Tier-2 funding from the Engineering and Physical Sciences Research Council (capital grant EP/P020259/1), and DiRAC funding from the Science and Technology Facilities Council (www.dirac.ac.uk).

SUPPLEMENTARY INFORMATION VERTEX CALCULATIONS IN MLDE
The dynamical susceptibility fermionic matrix is computed by inverting the Bethe-Salpeter Equation (BSE) where : and Γ imp is the impurity vertex obtained by inverted the BSE for the MLDE impurity problem.The impurity susceptibility reads: In our work, we compute the two-particle response function G 2 via the Lehman representation [29,30]: where we have introduced the closure relation, and define O 1 = c σ , O 2 = c † σ and O 3 = c σ .E i are the eigenvalues of the MLDE hamiltonian, Π are the triplet permutations, and the summation holds over all states of the Hilbert space.
We emphasize that it is usually difficult to provide energy cutoffs on the low energy converged states of the Hamiltonian (typically in Lanczos or other iterative approaches), as the Boltzmann statistics only truncates one of the sum (on index i), whereas all eigenstates are required for the other summation.MLDE provides an ideal candidate for such calculations, as the number of bath sites can be kept small with very little cost in accuracy.We note indeed that the previous equation can be arranged to show a computational complexity growing with the cube of the number of eigenstates n, where n itself grows exponentially with the number of bath sites.and compared with the exact vertex (dashed line) at temperature T /t = 0.025 and U/t = 6.Respectively χ (f) and Γ (g) obtained in the Mott phase for U/t = 12.The agreement is remarkable.
3.a).MLDE captures well the localization-delocalization transition (see Fig.3.b),where the MLDE with N b = 2 only differs near the transition, and with N b = 3 is exact.

Figure 1 :
Figure 1: Model hybridisation.Real part and imaginary part of the model hybridisation used for the MLDE benchmark.

Figure 2 :
Figure 2: MLDE benchmark.(a) Imaginary and (b) real part of the MLDE with N b = 3 self-energy (dashed lines) obtained from on-site correlation U in range 2 − 15 for a test hybridisation function for a half-filled impurity (units are arbitrary) and compared with continuous time Monte Carlo (continuous line).The scales are in arbitrary units.

Figure 3 :
Figure 3: Strength of correlations.(a) Imaginary part of the self-energy obtained by ED for N b = 2 (dashed lines) compared to the exact solution (continuous line).(b) Quasi-particle weight obtained by MLDE with N b = 2 and N b = 3 compared to the continuous time Monte Carlo.The MLDE remains close to the exact solution with N b = 2, whereas for the same number of bath sites ED vastly overestimates the strength of correlations (see panel (a)).The scales are in arbitrary units.

Figure 4 :
Figure 4: Kondo physics.(a) TOV organic molecule deposited on a gold substrate at T = 5K.(b) Imaginary parts of the molecule' hybridisation function (AIM), represented by MLDE ∆ 1,2 (iω n ) (MLDE), and as obtained by the ED solver with N b = 12.In the Kondo regime, the hybridisation remains constant up to the lowest frequency, which is not captured in ED (c) Imaginary part of the self-energy obtained by MLDE shows a Fermi liquid type behaviour at small frequency (circles), whereas the self-energy shows an artificial Mott singularity, due to the bath discretization.We note that this low temperature was beyond the reach of our CTQMC solver (β = 35000).(d) Imaginary part of the self-energy in real axis frequencies, MLDE shows a small dip near the Fermi level, which is associated with the Kondo scale, whereas ED shows an artificial peak instead.(e) Time dynamics of the magnetically quenched system: at time t = 0 the system is polarised along the z-axis, upon relaxation we report the time dependence of the magnetic precession at a temperature T = 100K (T = 20K) above (below) the Kondo temperature (T K ≈ 30K).The Fourier transform of the dynamics at T = 100K shows a peak at ω = 0.0029eV, near the Kondo energy T K = 0.0027eV.The dynamics at T = 20K shows sub-harmonics reminiscent from the Zeeman splitting at time t = 0.

Figure 5 :
Figure 5: Mott transition.(a) Spectral function obtained for the square lattice Hubbard model solved by MLDE (N b = 3 for all calculations) with increasing values of the on-site dimensionless correlation ratio U/t.For U smaller than the critical value U c ≈ 9, the system remains metallic, with a quasi-particle weight (sharp peak at zero), Hubbard lower and upper bands, and band edge satellites that develop near U/t ≈ 6.In the Mott phase (U/t > 10) we obtain a Mott charge gap ∆ ≈ U − W , where W/t = 8 is the bandwidth for the square lattice.(b) Imaginary part of the self-energy obtained for U/t = 12 at different MLDE iteration of the self-consistent cycle.At iteration 14, the self-energy is converged, and within the error bar or the converged CTQMC solution.(c) Converged imaginary part of the MLDE Green's function for different values of U/t (symbols) compared with the DMFT CTQMC solution (dashed lines).Note the real part of the Green's function is zero in a particle-hole symmetric system.

Figure 6 :
Figure 6: Dynamical susceptibility.Momentum resolved spin susceptibility obtained by Bethe Salpeter with the vertex calculated with MLDE in the Hubbard model with U/t = 12, at increasing temperature T /t = 0.025 (a),T /t = 1 (b), T /t = 10 (c).We report that the magnetic susceptibility at (π, π) is enhanced as the system reaches the meltdown of the Mott gap (b), and then at very high temperature (c) becomes uniform.As the MLDE calculations only involves N b = 3 bath sites, the vertex is fully tractable at any temperature.d) Magnetic susceptibility χ and e) irreducible vertex Γ resolved in fermionic frequency iν obtained by MLDE (continuous line) and compared with the exact vertex (dashed line) at temperature T /t = 0.025 and U/t = 6.Respectively χ (f) and Γ (g) obtained in the Mott phase for U/t = 12.The agreement is remarkable.