Few-body systems capture many-body physics: tensor network approach

TL;DR

This paper introduces a tensor network-based method that constructs small, manageable models mimicking infinite many-body systems, enabling accurate analysis of ground states and phase transitions in higher dimensions.

Contribution

A novel scheme using an entanglement bath to accurately simulate infinite many-body systems with few-body models, facilitating experimental and theoretical studies.

Findings

Accurately captures ground-state properties and phase transitions.

Efficient and sign-problem-free approach.

Applicable to higher-dimensional lattices.

Abstract

Due to the presence of strong correlations, theoretical or experimental investigations of quantum many-body systems belong to the most challenging tasks in modern physics. Stimulated by tensor networks, we propose a scheme of constructing the few-body models that can be easily accessed by theoretical or experimental means, to accurately capture the ground-state properties of infinite many-body systems in higher dimensions. The general idea is to embed a small bulk of the infinite model in an "entanglement bath" so that the many-body effects can be faithfully mimicked. The approach we propose is efficient, simple, flexible, sign-problem-free, and it directly accesses the thermodynamic limit. The numerical results of the spin models on honeycomb and simple cubic lattices show that the ground-state properties including quantum phase transitions and the critical behaviors are accurately captured by only $\mathcal{O}(10)$ physical and bath sites. Moreover, since the few-body Hamiltonian only contains local interactions among a handful of sites, our work provides new ways of studying the many-body phenomena in the infinite strongly-correlated systems by mimicking them in the few-body experiments using cold atoms/ions, or developing novel quantum devices by utilizing the many-body features.

Tables

Table 1. Table 1: The values of the critical field h c subscript ℎ 𝑐 h_{c} with the perturbation expansions (PE) Pfeuty and Elliott ( 1971 ) , cluster quantum Monte Carlo (cQMC) Blöte and Deng ( 2002 ) , linked-cluster expansions (LCE) Coester et al. ( 2016 ) , mean field theory (MFT) Suzuki et al. ( 2013 ) , and our AOP simulations in Stages one and two.

	PE	cQMC	LCE	MFT	Stage one	Stage two
$h_c$	2.60	2.58	2.65	3	2.8	2.66

Equations

$$\hat{H}=\sum _{⟨i,j⟩}\hat{H}(i,j).$$

$$M_s\propto (h_c-h_x{)}^{{\beta }^{\ast }},$$

$$\xi \propto (h_c-h_x{)}^{-\sigma },$$

$${\hat{F}}^{\partial }(i,j)=\sum _a{\hat{F}}^L(i{)}_a\otimes {\hat{F}}^R(j{)}_a.$$

$$\hat{\mathcal{F}}(i,j{)}_{a_1{a}_2{a}_3{a}_4}={\hat{F}}^L(i{)}_{a_1}{\hat{F}}^L(i{)}_{a_2}{\hat{F}}^R(j{)}_{a_3}{\hat{F}}^R(j{)}_{a_4}{\tilde{H}}^B(i,j),$$

$$\hat{\mathcal{H}}(i,j{)}_{{\mu }_1^{\prime }{\mu }_2^{\prime }{\mu }_3^{\prime }{\mu }_4^{\prime },{\mu }_1{\mu }_2{\mu }_3{\mu }_4}=\sum _{a_1{a}_2{a}_3{a}_4}(V^{[1]}{|}_{a_1{\mu }_1{\mu }_1^{\prime }}(V^{[2]}{|}_{a_2{\mu }_2{\mu }_2^{\prime }}\hat{\mathcal{F}}(i,j{)}_{a_1{a}_2{a}_3{a}_4}|V^{[3]}{)}_{a_3{\mu }_3{\mu }_3^{\prime }}|V^{[4]}{)}_{a_4{\mu }_4{\mu }_4^{\prime }},$$

$$M_{a_1{\mu }_1{\mu }_1^{\prime },a_3{\mu }_3{\mu }_3^{\prime }}^{[1]}=\sum _{a_2{a}_4{\mu }_2{\mu }_2^{\prime }{\mu }_4{\mu }_4^{\prime }}(V^{[2]}{|}_{a_2{\mu }_2{\mu }_2^{\prime }}⟨{\tilde{A}}^{[1]}(i,j){|}_{{\mu }_1^{\prime }{\mu }_2^{\prime }{\mu }_3^{\prime }{\mu }_4^{\prime }}\hat{\mathcal{F}}(i,j{)}_{a_1{a}_2{a}_3{a}_4}|{\tilde{A}}^{[1]}(i,j){⟩}_{{\mu }_1{\mu }_2{\mu }_3{\mu }_4}|V^{[4]}{)}_{a_4{\mu }_4{\mu }_4^{\prime }},$$

$$M_{a_2{\mu }_2{\mu }_2^{\prime },a_4{\mu }_4{\mu }_4^{\prime }}^{[2]}=\sum _{a_1{a}_3{\mu }_1{\mu }_1^{\prime }{\mu }_3{\mu }_3^{\prime }}(V^{[1]}{|}_{a_1{\mu }_1{\mu }_1^{\prime }}⟨{\tilde{A}}^{[2]}(i,j){|}_{{\mu }_1^{\prime }{\mu }_2^{\prime }{\mu }_3^{\prime }{\mu }_4^{\prime }}\hat{\mathcal{F}}(i,j{)}_{a_1{a}_2{a}_3{a}_4}|{\tilde{A}}^{[2]}(i,j){⟩}_{{\mu }_1{\mu }_2{\mu }_3{\mu }_4}|V^{[3]}{)}_{a_3{\mu }_3{\mu }_3^{\prime }},$$

$$M_{a_1{\mu }_1{\mu }_1^{\prime },a_3{\mu }_3{\mu }_3^{\prime }}^{[3]}=\sum _{a_2{a}_4{\mu }_2{\mu }_2^{\prime }{\mu }_4{\mu }_4^{\prime }}(V^{[2]}{|}_{a_2{\mu }_2{\mu }_2^{\prime }}⟨{\tilde{A}}^{[3]}(i,j){|}_{{\mu }_1^{\prime }{\mu }_2^{\prime }{\mu }_3^{\prime }{\mu }_4^{\prime }}\hat{\mathcal{F}}(i,j{)}_{a_1{a}_2{a}_3{a}_4}|{\tilde{A}}^{[3]}(i,j){⟩}_{{\mu }_1{\mu }_2{\mu }_3{\mu }_4}|V^{[4]}{)}_{a_4{\mu }_4{\mu }_4^{\prime }},$$

$$M_{a_2{\mu }_2{\mu }_2^{\prime },a_4{\mu }_4{\mu }_4^{\prime }}^{[4]}=\sum _{a_1{a}_3{\mu }_1{\mu }_1^{\prime }{\mu }_3{\mu }_3^{\prime }}(V^{[1]}{|}_{a_1{\mu }_1{\mu }_1^{\prime }}⟨{\tilde{A}}^{[4]}(i,j){|}_{{\mu }_1^{\prime }{\mu }_2^{\prime }{\mu }_3^{\prime }{\mu }_4^{\prime }}\hat{\mathcal{F}}(i,j{)}_{a_1{a}_2{a}_3{a}_4}|{\tilde{A}}^{[4]}(i,j){⟩}_{{\mu }_1{\mu }_2{\mu }_3{\mu }_4}|V^{[3]}{)}_{a_3{\mu }_3{\mu }_3^{\prime }}.$$

$$|A(i,j){⟩}_{{\mu }_1{\mu }_2{\mu }_3{\mu }_4}=\sum _{\nu }|{\tilde{A}}^{[3]}(i,j){⟩}_{{\mu }_1{\mu }_2\nu {\mu }_4}{R}_{\nu {\mu }_3}^{[3]}.$$

$$\sum _{{\mu }_1{\mu }_2{\mu }_4}⟨{\tilde{A}}^{[3]}(i,j){|}_{{\mu }_1{\mu }_2{\mu }_3{\mu }_4}|{\tilde{A}}^{[3]}(i,j){⟩}_{{\mu }_1{\mu }_2{\mu }_3^{\prime }{\mu }_4}=I_{{\mu }_3{\mu }_3^{\prime }}.$$

$$\hat{\rho }(i,j)\simeq \sum _{{\mu }_1{\mu }_2{\mu }_3{\mu }_4}|A(i,j){⟩}_{{\mu }_1{\mu }_2{\mu }_3{\mu }_4}⟨A(i,j){|}_{{\mu }_1{\mu }_2{\mu }_3{\mu }_4}.$$

$${\hat{\mathcal{H}}}^{\partial }(n,x{)}_{{\mu }_x{\mu }_x^{\prime }}=\sum _{a_x}{\hat{F}}^{L(R)}(n{)}_{a_x}|V^{[x]}{)}_{a_x{\mu }_x{\mu }_x^{\prime }}.$$

$${\hat{\mathcal{H}}}^{FB}=\prod _{⟨x\in cluster,n\in bath⟩}{\hat{\mathcal{H}}}^{\partial }(n,x)\sum _{⟨i,j⟩\in cluster}[I-\tau \hat{H}(i,j)].$$

$$\tilde{\rho }(i,j,\cdots )=\sum _{\{\mu \}}|\tilde{\Phi }(i,j,\cdots ){⟩}_{\{\mu \}}⟨\tilde{\Phi }(i,j,\cdots ){|}_{\{\mu \}}.$$

$${\hat{\mathcal{H}}}^{\partial }(n,x)=I-\tau {\hat{H}}^{\partial }(n,x)+O({\tau }^2).$$

$${\hat{H}}^{FB}=\sum _{⟨i,j⟩\in cluster}\hat{H}(i,j)+\sum _{⟨x\in cluster,n\in bath⟩}{\hat{H}}^{\partial }(n,x).$$

$${\hat{H}}^{\partial }(n,x)=\sum _{\alpha {\alpha }^{\prime }}{J}_{\alpha {\alpha }^{\prime }}(n,x){\hat{\mathcal{S}}}^{{\alpha }^{\prime }}(n){\hat{S}}^{\alpha }(x).$$

$${\hat{H}}_{Inf}=\sum _n{\hat{H}}_{n,n+1},$$

$${\hat{F}}^{\partial }(i,j)=\sum _a{\hat{F}}^L(i{)}_a\otimes {\hat{F}}^R(j{)}_a.$$

$$\hat{\mathcal{F}}(i,j{)}_{a_1{a}_2}={\hat{F}}^R(i{)}_{a_1}{\hat{F}}^L(j{)}_{a_2}{\tilde{H}}^B(i,j),$$

$$\hat{\mathcal{H}}(i,j{)}_{{\mu }_1^{\prime }{\mu }_2,{\mu }_1{\mu }_2}=\sum _{a_1{a}_2}(V^{[1]}{|}_{a_1{\mu }_1{\mu }_1^{\prime }}\hat{\mathcal{F}}(i,j{)}_{a_1{a}_2}|V^{[2]}{)}_{a_2{\mu }_2{\mu }_2^{\prime }},$$

$$M_{a_1{\mu }_1{\mu }_1^{\prime },a_2{\mu }_2{\mu }_2^{\prime }}^{[1]}=⟨{\tilde{A}}^{[1]}(i,j){|}_{{\mu }_1^{\prime }{\mu }_2^{\prime }}\hat{\mathcal{F}}(i,j{)}_{a_1{a}_2}|{\tilde{A}}^{[1]}(i,j){⟩}_{{\mu }_1{\mu }_2},$$

$$M_{a_1{\mu }_1{\mu }_1^{\prime },a_2{\mu }_2{\mu }_2^{\prime }}^{[2]}=⟨{\tilde{A}}^{[2]}(i,j){|}_{{\mu }_1^{\prime }{\mu }_2^{\prime }}\hat{\mathcal{F}}(i,j{)}_{a_1{a}_2}|{\tilde{A}}^{[2]}(i,j){⟩}_{{\mu }_1{\mu }_2}.$$

$$|A(i,j){⟩}_{{\mu }_1{\mu }_2}=\sum _{\nu }|{\tilde{A}}^{[1]}(i,j){⟩}_{\nu {\mu }_2}{R}_{\nu {\mu }_1}^{[1]}.$$

$$\sum _{{\mu }_2}⟨{\tilde{A}}^{[1]}(i,j){|}_{{\mu }_1{\mu }_2}{\tilde{A}}^{[1]}(i,j){⟩}_{{\mu }_1^{\prime }{\mu }_2}=I_{{\mu }_1{\mu }_1^{\prime }}.$$

$${\hat{\mathcal{H}}}^{\partial }(n,x{)}_{{\mu }_x{\mu }_x^{\prime }}=\sum _{a_x}{\hat{F}}^{L(R)}(n{)}_{a_x}|V^{[x]}{)}_{a_x{\mu }_x{\mu }_x^{\prime }}.$$

$${\hat{H}}^{FB}=\sum _{⟨i,j⟩\in cluster}\hat{H}(i,j)+\sum _{⟨x\in cluster,n\in bath⟩}{\hat{H}}^{\partial }(n,x).$$

$${\hat{\mathcal{H}}}^{FB}=\prod _{⟨x,n⟩}{\hat{\mathcal{H}}}^{\partial }(n,x).$$

$${\hat{\mathcal{H}}}^{FB}=\sum _{\{\alpha \}}\prod _{⟨x,n⟩}{\hat{\mathcal{H}}}^{\partial }(n,x{)}^{{\alpha }_{x,n}}.$$

$$\xi \sim {\mathcal{D}}^{\kappa },$$

$$\xi =\frac{\tau }{\log {\Lambda }_0-\log {\Lambda }_1},$$

Full text

I Introduction

Quantum many-body systems in one, two and three dimensions. Investigating the ground states and low-lying states of strongly correlated quantum many-body systems is one of the most important challenges in modern physics. It lies in the centre of interest of condensed matter physics Auerbach (1994); Wen (2007), atomic, molecular and optic physics Lewenstein et al. (2012), and high energy physics Montvay and Münster (1997); Wiese (2013). Fundamentally, these systems may exhibit exotic states and phenomena, such as spin liquids Balents (2010) and topological phases Qi and Zhang (2011); Bansil et al. (2016). On the other hand, these systems have important applications in contemporary electronics, superconductivity Lee (2008), spintronics Žutić et al. (2004) and, more recently, quantum information roa . The future quantum technologies, i.e. quantum computers, quantum simulators and annealers, quantum metrology and sensing rely essentially exclusively on the use of strongly-correlated quantum lattice systems. However, the high complexity rising from strong correlations makes the exact solutions/diagonalizations impossible or inefficient in most cases. Numerical methods, benefiting from fast development of the computer technology, become nowadays the most frequently used tools capable of reliably studying such systems. Developments of new more efficient methods, with lower cost and higher accuracy, are therefore highly demanded.

Technically, one-dimensional (1D) models are the simplest, although quantum fluctuations in 1D are particularly large Giamarchi (2006). The 1D systems play important roles in electronics and spintronics, as they provide specific possibilities in controlling transport and reveal exotic excitations such as Majorana fermions Nadj-Perge et al. (2010); Suominen et al. . They can be naturally viewed as the edges of two-dimensional (2D) systems, and may correspond to edge states of these 2D systems. Wen (2007); Qi and Zhang (2011). 2D systems are obviously more demanding numerically and experimentally, whereas from a physical point of view they can be taken as playgrounds for novel concepts and exotic states such as anyonic excitations Wen (2007); Yoshioka (2002), frustrated antiferromagnetism Misguich et al. (2005); Sachdev (2008), spin liquids Savary and Balents (2017), topological order and topological phase transitions, and graphene-like systems Neto et al. (2009), etc.

In principle, the three-dimensional (3D) models are even more interesting, as they are much closer to reality of our daily experience. Because of their extreme complexity, the adoption of various approximations to treat them is totally unavoidable. The three dimension is closer to the upper critical dimension, and one may expect that the mean-field theories would work well for them. A paradigmatic example is the Bose-Hubbard model, which can be nicely explored by bosonic dynamic mean-field theory (DMFT) Anders et al. (2010). Such few, but well-controlled systems can serve as validation, calibration and benchmark for various numerical and analytical methods. Still, there are also 3D models that are extremely demanding to be understood, such as, among others, the spin ice Nisoli et al. (2013a) with pyrochlore lattice Gingras and McClarty (2014), that is a highly frustrated magnet; the Fermi–Hubbard model, which is usually invoked to describe high temperature superconductivity of cuprates that consist of strongly correlated 2D planes weakly coupled in the transverse direction (Lee (2008), see also Imriška et al. (2014) for a quantum simulation with ultracold atoms). The (3+1)D lattice gauge theories at high densities and temperature are also beyond the possibilities of the existing codes and machines. Generally, different approximate analytical methods might generate converse results, leading to unnecessary controversies in many cases. It turns out that finding reliable and efficient numerical methods to solve 3D quantum many-body problems becomes indeed imperative.

Tensor networks — state-of-the-art. The density matrix renormalization group (DMRG) White (1992) is widely recognized as a major breakthrough in the calculations of the ground states in 1D systems. Originally proposed as a mere numerical tool, the reformulation of the DMRG as a variational algorithm in terms of matrix product states (MPS) Schollwöck (2011) leads to the proposal of more general formalism, based on tensor networks (TN’s) Cirac and Verstraete (2009); Orús (2014); Verstraete et al. (2008). TNs provide a very general ansatz for the wave functions: the quantities of interest may be expressed as results of the contraction of a network of local tensors. It has rapidly evolved into a promising powerful tool to study large or even infinite size systems in two dimensions. In fact, TNs overcome most of the limitations of the standard numerical algorithms: for instance, in contrast to quantum Monte Carlo (QMC), TNs do not suffer from the notorious “negative-sign” problem White et al. (1989) and allow for an accurate access to frustrated spin systems and fermionic models away from half-filling.

To what extent the TN is feasible depends on the amount of entanglement of the states to be simulated. The efficiency (computational memory and time) of the TN approaches is also determined by the capability of the current computers. In the standard formulations, TN works for low-entangled states such as the ground states of local and gapped Hamiltonians. For these states, an area law for the entanglement entropy holds, i.e. the entanglement entropy of a subsystem (consisting of a large, but finite block) scales with the block’s boundary Eisert et al. (2010). This fact explains the efficiency of the MPS-based algorithms in 1D. For the same reason, the MPS-based algorithms (e.g. DMRG) work well for small 2D systems, but are strongly limited when the size grows White and Scalapino (1998); Stoudenmire and White (2012); Verstraete et al. (2008). By acknowledging this, many different competing approaches have been developed. Among others, a purely 2D ground-state TN ansatz, termed projected entangled pair state (PEPS), was proposed as a natural extension of MPS. PEPS fulfils the 2D area law of entanglement entropy Verstraete et al. (2008); Verstraete and Cirac ; Jordan et al. (2008); Cirac and Verstraete (2009); Orús (2014); Verstraete et al. (2006), while the multiscale entanglement renormalization ansatz (MERA) Vidal (2007a) bears particular advantages for studying critical models.

Within the existing TN algorithms, a lot of works were done on 2D quantum as well as 3D classical models, where the simulations consist in the contractions of 3D TN’s Vidal (2007b, a); Verstraete and Cirac ; Jordan et al. (2008); Jiang et al. (2008); Levin and Nave (2007); Gu et al. (2008); Xie et al. (2009, 2012); Yang et al. (2017); Nishino and Okunishi (1996); Orús and Vidal (2009); Orús (2012); Nishino et al. (2000); Corboz (2016); Vanderstraeten et al. (2016); Phien et al. (2015); García-Sáez and Latorre (2013). This well-known quantum-classical equivalence Suzuki (1976) becomes very explicit in the TN terminology, and is utilized frequently in the TN approaches for ground-state Vidal (2007b); Verstraete and Cirac ; Jordan et al. (2008); Jiang et al. (2008) and thermodynamic Bursill et al. (1996); Li et al. (2011); Ran et al. (2012, 2013); Czarnik et al. (2012) studies on discrete and even continuous Verstraete and Cirac (2010a) systems. However, for 3D ground-state simulations, we are essentially facing the contractions of four-dimensional TN’s Teng (2017); Braiorr-Orrs et al. (2016), which are hardly treatable even with small bond dimensions. Therefore, developing efficient 3D quantum algorithms are strongly desired, in particular for infinite quantum systems.

“Bath-stimulated” methods. For 3D quantum models, many interesting issues remain to be explored or even unsolved to a large extent Balents (2010); Guo et al. (1994); Nisoli et al. (2013b); Moessner and Ramirez (2006); Ramirez (1994). They have been the subject of intensive studies in recent years and many numerical methods were developed to handle them. Several approaches were proposed beyond the standard mean-field and renormalization group methods, such as the linked cluster expansions Oitmaa et al. (2006); Joshi et al. (2015); Coester et al. (2016), and the functional renormalization group method Metzner et al. (2012); Iqbal et al. (2016). On the other hand, the numerical simulations are extremely challenging, and the finite-size algorithms, including exact diagonalization (ED), QMC and DMRG, suffer severe finite-size effects, which are quite consuming for large systems and can access infinite systems only by utilizing finite-size scaling.

To treat the correlations in many body systems, one usually starts by evoking the ideas of “mean field”, “bath” or boundary conditions. Analytical methods such as the Hartree-Fock mean-field theory and the saddle point approximation in path integral are commonly used. In fact, for lattice models the “mean field” idea goes back to the single-site Weiss method, applied first for classical magnetic models Weiss (1907). Contemporary mean field methods for lattice models include Guztwiller ansatz for bosons and/or fermions, or pairing approaches (Bogoliubov-like for bosons, or Bardeen-Cooper-Schrieffer-like for fermions) –- for an overview of these and other methods see Ref. Lewenstein et al. (2012) and references therein. In the context of the present work it is important to mention the “cluster mean field theory” (CMFT), where the mean field à la Weiss is combined with exact diagonalization on clusters (for recent developments see Yamamoto (2009); Ren et al. (2014) and references therein). It is also worth mentioning “entanglement mean field theory” (EMFT) De and Sen (2012a, b), which for spin models is formulated on few spin clusters, demanding self-consistency of entanglement properties. Both CMFT and EMFT are close to the standard MFT that can give quite accurate description of standard (Landau-Ginsburg-like) ordered and disordered phases, but typically only far from criticality.

For thermal and open systems, one popular way is to introduce a “heat bath” to mimic the interactions between the system and the environment Zurek (1982). Regarding numerical approaches, the density functional theory (DFT), also known as ab-initio first-principle calculations Sherrill (2010), was built by extending the Thomas-Fermi approximation of homogeneous electron gas to the inhomogeneous electron system Thomas (1927). Its huge success in condensed matter physics, quantum chemistry, and materials science largely relies on the simplicity and unification, “using a popular code, a standard basis, and a standard functional approximation” Sherrill (2010).

In order to handle strong correlations, several schemes were developed in the spirit of DFT. The examples include the dynamic mean-field theory Metzner and Vollhardt (1989); Georges and Kotliar (1992); Georges et al. (1996) that maps a lattice model (such as the Hubbard model) onto a quantum impurity model subject to self-consistent conditions, and the density matrix embedding theory (DMET) Knizia and Chan (2012) that was proposed aiming at a better consideration of the entanglement, thanks to the accompanying explosive advances in both quantum information science Horodecki et al. (2009) and condensed matter physics Amico et al. (2008). However, it is difficult to use these algorithms to study long-range ordered states or phase transitions. To probe the disordered ground-states (e.g. the spin liquids in the infinite frustrated systems), it was proposed to signal the disordered nature by simulating a finite system with random boundary conditions Yuste et al. .

Our proposal: mimicking many-body systems by few-body ones. In general, as illustrated in Fig. 1, the central idea of our work is to optimally find the few-body Hamiltonian to mimic the infinite model, without any prior knowledge of the ground state. The few-body model contains the physical sites in a finite cluster and the “bath” sites around it. The few-body Hamiltonian consists of two parts: the interactions among the physical sites (blue balls) within the cluster, and those [Eq. (14)] between the boundary physical sites and the bath sites (red balls). The physical-bath interactions are represented by some local Hamiltonians, which reproduce the quantum entanglement between the cluster and the bath, in such a way that the many-body effects from the infinite environment are well captured in the few-body simulations. Then the ground-state information of the infinite system is encoded in the reduced density matrix of the few-body ground state after tracing the bath degrees of freedom.

The theoretical scheme we propose is a higher-dimensional generalization of the ab-initio optimization principle (AOP) formulated with TN Ran (2016), and originally developed for infinite 1D systems with translational invariance. The idea is to find the simplest eigenvalue equations that encodes the infinite TN contraction problem. Besides its simplicity in the implementation, the 1D AOP has proved to have several computational advantages over other established algorithms, such as iTEBD and iDMRG White (1992). For the purpose of the present work, the main advantage of the AOP is its flexibility and implications in high dimensions: without any substantial conceptual changes, the AOP can be readily extended to 2D and even 3D systems with high efficiency. Furthermore, the dynamic correlation length and the first excitation gap can be straightforwardly extracted.

Our scheme consists of two stages: (1) compute physical-bath Hamiltonian and (2) solve the few-body Hamiltonian. In the first stage, by choosing the dimension $D$ of the bath site and a supercell that obeys the translational invariance, we start from the original Hamiltonian of the system and construct a set of self-consistent eigenvalue equations. Their solution gives the Hamiltonian $\mathcal{\hat{H}}^{\partial}$ [Eq. (14)] that describes the interactions between a physical and a bath sites. Such equations in fact encode an optimal zero-loop TN approximation of the state. This approximation directly enters the thermodynamic limit with a Bethe TN state ansatz Lepetit et al. (2000); Nakatani and Chan (2013), and already gives us the first glance of the ground state with good accuracy especially for the gapped states Jiang et al. (2008); Ran et al. (2012, 2013); Braiorr-Orrs et al. (2016); Lepetit et al. (2000); Nakatani and Chan (2013).

The aim of the second stage is to construct the few-body Hamiltonian $\mathcal{\hat{H}}^{FB}$, and solve its ground state $|\tilde{\Phi}\rangle$ by e.g., DMRG (with certain dimension cut-off’s denoted by $\chi$). $\mathcal{\hat{H}}^{FB}$ is formed by all the physical interactions inside a chosen cluster and several physical-bath interactions given by $\mathcal{\hat{H}}^{\partial}$. The choice of the cluster is very flexible. The ground-state properties of the infinite system is then encoded in the ground state $|\tilde{\Phi}\rangle$ of $\mathcal{\hat{H}}^{FB}$. In other words, quantities such as energy, magnetization and entanglement of the infinite ground state are obtained from the density matrix of $|\tilde{\Phi}\rangle$ by tracing all the bath degrees of freedom [Eq. (16)].

Our numerical results show, for instance, that the Heisenberg model on infinite honeycomb lattice is accurately simulated by a $\mathcal{\hat{H}}^{FB}$ that only contains $N_{p}=18$ physical sites surrounded by $N_{b}=12$ bath sites. For the 3D Heisenberg models on infinite simple cubic lattice, the ground-state properties including the critical behaviors near the quantum phase transition point are faithfully captured with only $N_{p}=8$ physical and $N_{b}=24$ bath sites. The discrepancy (such as energy) compared with the state-of-the-art TN algorithm is around $\mathcal{O}(10^{-3})$.

The algorithm built from our scheme possesses several advantages (see Appendix G). The algorithm can directly reaches the thermodynamic limit by means of the physical-bath interactions on the boundary, thus has no conventional finite-size effects compared with the finite algorithms such as ED and DMRG. The strongly-correlated effects of the infinite models are accurately considered, and the many-body features, e.g., entanglement and criticality, can be efficiently captured, thus our scheme goes beyond the mean-field-based methods such as DFT Sherrill (2010) and DMFT Metzner and Vollhardt (1989); Georges and Kotliar (1992); Georges et al. (1996). Comparing with DMFT where the original model is approximated by an impurity model in a bath, we approximate the infinite-size system into a few-body model that contains its original interactions and the emergent physical-bath interactions. The accuracy is enhanced by fully considering all interactions in the cluster, thus outperforms the Bethe TN-based algorithms Jiang et al. (2008); Ran et al. (2012, 2013); Braiorr-Orrs et al. (2016); Lepetit et al. (2000); Nakatani and Chan (2013). In higher dimensions, the computational cost of our scheme is much lower than, e.g., the TN renormalization group algorithms Levin and Nave (2007); Vidal (2007a); Verstraete and Cirac ; Verstraete and Cirac ; Nishino and Okunishi (1996); Orús and Vidal (2009); Orús (2012); Nishino et al. (2000). It has no sign problem White et al. (1989) thus can be used to simulate frustrated and fermionic systems.

The construction of $\mathcal{\hat{H}}^{FB}$ makes it possible to investigate the many-body effects in experiments by designing the few-body models — quantum simulators described by the predicted Hamiltonians. The many-body behaviors are expected to be observed in the bulk of the few-body model. The feasibility of realizing $\mathcal{\hat{H}}^{FB}$ in cold atom experiments is supported by several facts observed in our numerical simulations: the few-body Hamiltonian has the same interaction length as the original Hamiltonian; with a proper tolerance of error, say $\mathcal{O}(10^{-2})$, the size of the few-body model can be very small. Especially for spin-$1/2$ models on simple cubic lattice, we show that it is sufficient to use only the spin-$1/2$’s as the bath sites. The few-body Hamiltonian then is just a small spin-$1/2$ system that includes some special interactions (given by $\mathcal{\hat{H}}^{\partial}$) on the boundary.

II Numerical results

Heisenberg model on honeycomb lattice. We simulate the ground-state properties of the Heisenberg model on honeycomb lattice, which is on a gapless point and considered to be challenging to simulate. The Hamiltonian is the summation of the two-body interactions as

[TABLE]

For Heisenberg model, we have $\hat{H}(i,j)=J_{x}\hat{S}^{x}(i)\hat{S}^{x}(j)+J_{y}\hat{S}^{y}(i)\hat{S}^{y}(j)+J_{z}\hat{S}^{z}(i)\hat{S}^{z}(j)$, with $\hat{S}^{\alpha}(i)$ ($\alpha=x,y,z$) the $\alpha$ component of the spin-1/2 operators on the $i$-th site and $J_{\alpha}$ the coupling constants.

In stage one, the bath is calculated by choosing two neighboring sites as the supercell. It means $\mathcal{\hat{H}}^{FB}$ that appears in this stage contains $N_{p}=2$ physical and $N_{b}=4$ bath sites. For stage two, we choose two different clusters to construct $\mathcal{\hat{H}}^{FB}$, which contains $N_{p}=8$ physical sites surrounded by $N_{b}=8$ bath sites and $N_{p}=18$ physical sites with $N_{b}=12$ bath sites, respectively (Fig. 1). We utilize finite DMRG White (1992) to solve the ground state of $\mathcal{\hat{H}}^{FB}$.

The ground-state energy $E$ with different dimensions of the bath site $D=4$ and $8$ is shown in Fig. 2. One can see that $E$ converges rapidly by increasing the dimension cut-off of DMRG $\chi$ to $E=-0.543$ and $-0.544$ for the two clusters, respectively. With larger $D$, the bath will be able to carry more entanglement and lead to a better accuracy. The accuracy will also be improved by increasing $\chi$ since the result will approach to the exact ground state of $\mathcal{\hat{H}}^{FB}$ with no DMRG error. When $\chi$ is sufficiently large, the errors inside the cluster due to the tree approximation, Trotter discretization and truncations will vanish.

For a comparison, the ground-state energy by ED on such a cluster of 18 spins (Fig. 1) with periodic boundary condition is $E=-0.561$, which suffers severe finite-size effects. The result solely by bath calculation (tree approximation) is $E=-0.540$, and by second renormalization group (SRG) Xie et al. (2009) of TN is $E=-0.545$. SRG belongs to the state-of-the-art TN approaches for simulating 2D ground states with a high accuracy. The difference compared with our results are only $\mathcal{O}(10^{-3})$.

To further investigate the effects of the finiteness of the clusters, we calculate the nearest-neighbor bond energies $e=\langle\hat{S}_{i}\hat{S}_{j}\rangle$ and magnetizations with the cluster of $N_{p}=18$ and $N_{b}=12$ (Fig. 3). The changes of both quantities in different positions of the cluster is mostly $\mathcal{O}(10^{-2})$. By comparing with the tree results in the bath calculation and SRG, we find that the bond energies and magnetization on the boundary of the cluster are very close to the tree results, and in the middle where the “boundary effects” as well as the difference between our results and the SRG are minimal.

Our simulations show that without increasing the computational cost much, the finite-size effects are suppressed by introducing the entanglement bath, and at the same time the error from the tree approximation are reduced by choosing larger clusters.

Spin models on simple cubic lattice. We investigate the ground-state properties and quantum phase transitions in the spin models on simple cubic lattice. For bath calculations, the supercell is chosen to be two neighboring sites, giving a $\mathcal{\hat{H}}^{FB}$ with $N_{p}=2$ physical and $N_{b}=10$ bath sites. In stage two, we choose a cubic with $N_{p}=8$ physical and $N_{b}=24$ bath sites to construct $\mathcal{\hat{H}}^{FB}$ (inset of Fig. 4).

The ground-state energy $E$, staggered magnetization $M_{s}$ and dynamic correlation length $\xi$ of the antiferromagnetic Heisenberg model on simple cubic lattice are shown in Fig. 4. The energy converges to $E=-0.904$, while that from QMC Prokofev et al. (1998); Huang et al. (2016); Hua on a ($10\times 10\times 10$) lattice with periodic boundary condition is $E=-0.902$. Note the result from the tree approximation in the first stage is $E=-0.892$, which is already quite accurate. For $M_{s}$ and $\xi$, the finite-size effects are much stronger for our QMC calculations. The AOP simulations show that $M_{u}=0$ (uniform magnetization), $M_{s}=0.445$ and $\xi=0.405$. Our results are consistent with the widely accepted consensus, that its ground state is an antiferromagnetic ordered (Néel) state with a short correlation length.

We investigate the quantum phase transition of the anti-ferromagnetic Ising model in a transverse field on simple cubic lattice (Fig. 5), where the Hamiltonian reads $\hat{H}=\sum_{<i,j>}\hat{S}^{x}(i)\hat{S}^{x}(j)+h\sum_{i}\hat{S}^{z}(i)$. For a comparison, we try different dimension cut-offs with $D=2$, $\chi=10$, and $D=3$, $\chi=20$. The critical field is found to be around $h_{c}=2.66$, consistent with the results from other algorithms (Table LABEL:tab-Fields).

Our results show that from the few-body Hamiltonian, the scaling behavior in the critical region can be faithfully captured and the critical exponents are consistent with the results obtained by other methods. Meanwhile, obvious improvement in Stage two is observed compared with Stage one (Bethe approximation with a small cluster). By fitting the data in Stage two of $D=3$ and $\chi=20$ near the critical point, we find

[TABLE]

with $\beta^{\ast}=0.48$, which is close to but slightly larger than the perturbation expansions result $\beta^{\ast}=0.46$. Note that the exponent from the fitting on the data of Stage one is not reliable, which gives $\beta^{\ast}=1$.

We also calculate the dynamic correlation length $\xi$, which shows a peak at the critical point and scales as

[TABLE]

with $\sigma=0.25$ near the critical point in Stage two. In Stage one, we have $\sigma=0.23$. The exponent of the (spatial) correlation length by the perturbation expansions is $\sigma=0.5$ Pfeuty and Elliott (1971) . The discrepancy might be caused by the errors from both sides. Regarding the TN algorithms, the correlation length in the critical region will diverge with the scaling of the bath dimension $D$ as well as the DMRG dimension cut-off $\chi$ (unlike $M_{s}$ which converges to zero). Thus, it is difficult to directly extract the exponent of $\xi$ with fixed dimensions. The good thing is that the algebraic behavior of $\xi$ is clearly observed. What is open is how to get an accurate value of $\sigma$ by the scaling factors versus not only $h_{x}$ but also $\chi$ and $D$. See more discussions about the error of correlation functions with TN approaches in Appendix F.

III Method: higher-dimensional ab-initio optimization principle approach

The idea of AOP scheme Ran (2016) is, without any previous knowledge of the ground state, to transfer the infinite system to a finite one embedded in an entanglement bath. In the language of TN, the idea is to encode the contraction of an infinite TN in a simplest-possible local function that can be exactly computed, with smallest-possible number of inputs. The 1D version is briefly presented in Appendix A, where we also explain the relations and differences compared with the AOP approach in higher dimensions. To present the approach in high dimensions, we take the 2D spin model with nearest-neighbor couplings on honeycomb lattice as an example. The implementation can be easily generalized to other models on 2D and 3D lattices.

Stage One: calculate the entanglement bath. The first stage is to calculate the entanglement bath represented by a set of tensors dubbed as boundary tensors. They are obtained by solving a set of self-consistent eigenvalue equations [see Eqs. (6)-(10) below]. These equations are parametrized by the Hamiltonian as well as by the boundary tensors themselves, thus they can be solved in an alternative way: starting from an arbitrary guess, we update one boundary tensor by fixing all others as the parameters of the equations, and iterate such a procedure for every tensor until the fixed point is reached.

Though our method is based on the TN representation of the imaginary-time evolution with Trotter-Suzuki decomposition Suzuki and Inoue (1987) like several existing methods Vidal (2007b); Verstraete and Cirac ; Jordan et al. (2008); Nishino and Okunishi (1996); Orús and Vidal (2009); Orús (2012), the idea here is to encode the TN in the eigenvalue equations Not (a); Ran (2016) instead of contracting the TN. On the other hand, the implementation in this stage is borrowed from the generalization of DMRG on an infinite tree Lepetit et al. (2000); Nakatani and Chan (2013), which can be easily extended to 3D models with high efficiency. In the DMRG language, the (convergent) boundary tensors can be understood as the infinite environment of the tree brunches.

To begin with, one chooses a supercell that obeys the translational invariance, e.g. two sites connected by a parallel bond (see the smallest shaded circle in Fig. 1), and construct the cell tensor that parametrizes the eigenvalue equations. The bulk interaction is simply the coupling between these two spins, i.e. $\hat{H}^{B}(i,j)=\hat{H}(i,j)$. For the interactions between two neighboring supercells, we define the two-body operator $\hat{F}^{\partial}(i,j)=I-\tau\hat{H}(i,j)$ and do the singular value decomposition (Fig. 6) as

[TABLE]

We dub $a$ as the boundary index.

To obtain the TN of the imaginary-time evolution, we define the cell tensor that is the product of the (shifted) bulk Hamiltonian with $\hat{F}^{L}_{a}$ and $\hat{F}^{R}_{a}$ (Fig. 6) as

[TABLE]

with $\tilde{H}^{B}(i,j)=I-\tau\hat{H}^{B}(i,j)$. Note that both $\hat{F}^{\partial}(i,j)$ and $\tilde{H}^{B}(i,j)$ can be different from the current choice. $\hat{F}^{\partial}(i,j)$ can be different when the model has long-range interactions, and $\tilde{H}^{B}(i,j)$ can be different when choosing different subsystems to define the supercell. $\hat{\mathcal{F}}(i,j)$ Not (b) can be understood as a set of quantum operators defined in the Hilbert space of the supercell (spins $i$ and $j$) labeled by the boundary indexes $a_{1}$, $a_{2}$, $a_{3}$ and $a_{4}$. Similar to 1D AOP Ran (2016), $\tau$ is in fact the Trotter step, and $\hat{\mathcal{F}}(i,j)$ gives the TN representation of $I-\tau\hat{H}$ with an error $\mathcal{O}(\tau^{2})$.

Then, with the boundary tensors $|V^{[x]})$ (guessed or previously obtained in the last iteration) and the cell tensor $\hat{\mathcal{F}}(i,j)$, we define five eigenvalue equations as

[TABLE]

By solving the leading eigenvector of $\hat{\mathcal{H}}(i,j)$ given by Eq. (6), we obtain a tensor $|A(i,j)\rangle_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ dubbed as central tensor with $\mu_{1}$, $\mu_{2}$, $\mu_{3}$ and $\mu_{4}$ called virtual indexes according to the TN terminology. The central tensor can be considered as a state in the Hilbert space of the supercell labeled by four virtual indexes.

Meanwhile, $|V^{[x]})$ is obtained as the (left) leading eigenvector of $M^{[x]}$ [Eqs. (7)-(10)]. One can see that $M^{[x]}$ is defined by the isometries $|\tilde{A}^{[x]}(i,j)\rangle$ that is obtained by the QR decomposition of $|A(i,j)\rangle$ (referring to the $x$-th virtual bond $\mu_{x}$) of the central tensor. For example for $x=3$, we have (Fig. 7)

[TABLE]

$|\tilde{A}^{[3]}(i,j)\rangle$ is orthogonal, satisfying

[TABLE]

These isometries play the role of the renormalization group flow in the tree DMRG Lepetit et al. (2000); Nakatani and Chan (2013). Similarly, $|V^{[x]})_{a_{x}\mu_{x}\mu_{x}^{\prime}}$ can be understood as a “state” defined in the space of the boundary index $a_{x}$ labeled by $\mu_{x}$ and $\mu_{x}^{\prime}$ Not (c). The graphic representations of $\hat{\mathcal{H}}(i,j)$ and $M^{[3]}$ are given in Fig. 7 as examples.

One can see that these equations are parametrized by the solutions of others, and can be solved in an alternative way in practice. One can start with four random $|V^{[x]})$’s and calculate $|A(i,j)\rangle$ by solving the leading eigenvector of Eq. (6). Then one obtains $|\tilde{A}^{[x]}(i,j)\rangle$’s using Eq. (11) and update the $|V^{[x]})$’s according to Eqs. (7)-(10). Repeat this process until the central tensor and all boundary tensors converge.

In fact, the ground-state properties can already be well extracted by the central tensor $|A(i,j)\rangle$. For example, the reduced density matrix of the supercell $\hat{\rho}(i,j)=\text{Tr}_{/(i,j)}|\Psi\rangle\langle\Psi|$ (with $|\Psi\rangle$ denoting the ground state of the infinite model) is well approximated by the central tensor as

[TABLE]

Since each boundary tensor can be understood as the environment of an infinite tree branch, the original model is actually approximated at this stage by one defined on an infinite tree. Note that when only looking at the tree locally (from one site and its nearest neighbors), it looks the same to the original lattice. Thus, the loss of information is mainly long-range, i.e., from the destruction of loops. Though it has been shown numerically by many previous work that the tree approximation is very accurate especially for gapped systems Jiang et al. (2008); Ran et al. (2012, 2013); Braiorr-Orrs et al. (2016), we are still facing the difficulty of controlling the effects (errors) brought by such an approximation. More discussions about such a tree approximation are given in the Appendix B, starting from the state ansatz behind our approach. One can also find more details in the forth section of a recent review TNrev . To further improve the precision in a systematic way, the next stage is to embed a much larger subsystem in the entanglement bath.

Stage two: construct the few-body Hamiltonian and solve it. The second stage is to choose a finite cluster and use the obtained boundary tensors to construct a few-body Hamiltonian. All interactions inside the cluster will be fully considered to reduce the error from the tree approximation. The entanglement bath mimics the environment of the infinite tree branches, thus the algorithm directly accesses the thermodynamic limite and there is no conventional finite-size error that appears in, e.g. ED, DMRG or QMC.

The embedding is based on the generalizations of $\hat{\mathcal{H}}(i,j)$ [Eq. (6)] in stage one. From the formulation given above, one can see that $\hat{\mathcal{H}}(i,j)$ is actually the product of two parts. The first is the shifted Hamiltonian that contains all interactions inside the supercell (two neighboring sites in our example), and the second is in fact the physical-bath interactions (Fig. 8), whose Hamiltonian is written as

[TABLE]

Now we extend the supercell to a chosen larger cluster, where the few-body Hamiltonian denoted by $\mathcal{\hat{H}}^{FB}$ is written as

[TABLE]

Same as $\hat{\mathcal{H}}(i,j)$, $\mathcal{\hat{H}}^{FB}$ is also formed by two terms (Fig. 8). The first term is the product of several bath Hamiltonians that mimic the interactions between the cluster and the environment, and in the second term, the summation in Eq. (15) contains all couplings inside the cluster.

The entanglement bath only “interacts” with the physical sites nearby according to the coupling distance of the original Hamiltonian. In our example with nearest-neighbor couplings, every physical sites on the boundary interact with a bath site, and thus, the number of $\hat{\mathcal{H}}^{\partial}(n,x)$ in the product above scales with the length of the boundary of the cluster. For this reason, $\mathcal{\hat{H}}^{FB}$ is the product/summation of sparse or local matrices, and its ground state can be efficiently solved by using the finite-size approaches, such as ED or DMRG.

Note that if one takes the cluster as the supercell with two sites, Eq. (15) becomes exactly Eq. (6). The bath calculation itself can be considered as using ED to solve the $\mathcal{\hat{H}}^{FB}$ that contains only the supercell and the bath. The cluster can be arbitrarily chosen according to the computational capacity, and it does not have to obey the translational invariance of the model.

With the ground state $|\tilde{\Phi}(i,j,\cdots)\rangle_{\{\mu\}}$ of $\mathcal{\hat{H}}^{FB}$, the physical properties such as energy, magnetization, etc., can be obtained from the density operator $\tilde{\rho}$ by tracing all degrees of freedom of the bath sites as

[TABLE]

Again, this is a generalization of Eq. (13).

IV Discussions about experimental realizations

Our work provides a way of using few-body experiments to mimic many-body features of infinite systems. Since the few-body Hamiltonian only contains a handful of sites with local interactions, one could design cold-atom experiments to realize it in a lab. Specifically speaking in our examples, $\hat{\mathcal{H}}^{\partial}$ is the interaction between a physical spin and an artificial spin with $D$ (bath) degrees of freedom. Here, we assume that $\hat{\mathcal{H}}^{\partial}$ is Hermitian, which should be true due to the structure of the eigenvalue equations [Eqs. (7)-(10)] of the boundary tensors, where we have $|V^{[x]})_{a\mu\mu^{\prime}}=|V^{[x]})_{a\mu^{\prime}\mu}^{\ast}$. The task here is to get the coupling constants explicitly for implementing experiments.

To this end, let us transform $\hat{\mathcal{H}}^{\partial}$ to the standard summation form. We define $\hat{H}^{\partial}$ that satisfies

[TABLE]

It means to the first order of $\tau$, $\hat{\mathcal{H}}^{\partial}$ is the evolution operator of a Hamiltonian $\hat{H}^{\partial}$ for an infinitesimal imaginary time. This relation is true because in Eq. (14), $\hat{F}^{L(R)}$ is obtained by the decomposition of $I-\tau\hat{H}^{B}$, and the boundary tensor $|V^{[x]})$ has the similar structure since it forms an continuous MPS Tirrito et al. ; Verstraete and Cirac (2010b) in the imaginary time direction.

Then, the few-body Hamiltonian in Eq. (15) can be rewritten in a standard summation form as $\mathcal{\hat{H}}^{FB}=I-\tau\hat{H}^{FB}+O(\tau^{2})$ with

[TABLE]

The two summations contain the physical and physical-bath interactions, respectively, and all terms are local as discussed above. Again, $\mathcal{\hat{H}}^{FB}$ is the evolution operator of $\hat{H}^{FB}$ for an infinitesimal imaginary time to the first order of $\tau$, i.e. $\mathcal{\hat{H}}^{FB}\simeq e^{-\tau\hat{H}^{FB}}$.

The coupling constants of the physical-bath interactions can be calculated by expanding $\hat{H}^{\partial}$ as

[TABLE]

with $J_{\alpha\alpha^{\prime}}(n,x)$ the physical-bath coupling constants and $\{\hat{S}^{\alpha}\}$ and $\{\mathcal{\hat{S}}^{\alpha^{\prime}}\}$ the corresponding spin operators (including identity) that give the complete basis for the Hermitian matrices. $\{\hat{S}^{\alpha}\}$ is in fact the physical spin operators. For $\{\mathcal{\hat{S}}^{\alpha^{\prime}}\}$, one can generally choose the generators of SU(N) groups, which give a complete basis for an N-by-N Hermitian matrix. Then the bath spins should be SU(N) spins. If a symmetry McCulloch and Gulàcsi (2002); McCulloch (2007); Weichselbaum (2012) is used in the tensors, for example SU(2) symmetry for spin models, the bath spins are SU(2) spins with higher total momentum, and one will explicitly have the coefficients from the elements of $\hat{H}^{FB}$. Moreover, it is possible to translate the whole few-body Hamiltonian into the second-quantized picture, by expanding it with the bosonic or fermionic operators. The key is that the chosen operator basis have to completely expand the physical-bath Hamiltonian.

From our numerical results, we can see that the properties of the infinite model can be accurately mimicked by very small bath dimension $D$ and cluster size. Suppose we set the tolerance of the experimental error as $\mathcal{O}(10^{-2})$. In this case, the cluster can be chosen as two sites. Then we have $N_{p}=2$, with $N_{b}=4$ for honeycomb lattice and $N_{b}=10$ for simple cubic lattice. For the spin-$1/2$ models on simple cubic lattice, the dimension of the bath sites can be chosen as $D=2$. This means the bath spins are simply spin-$1/2$, same as the physical ones, which makes it easy to implement in experiments.

In short, the steps to mimic an infinite many-body system with a few-body model are as the following:

• •

Starting from the Hamiltonian of the target model [e.g. Eq. (1)], compute the physical-bath Hamiltonian $\mathcal{\hat{H}}^{\partial}$ [Eq. (14)] by our AOP algorithm.

• •

Write $\mathcal{\hat{H}}^{\partial}$ into $\hat{H}^{\partial}$ by Eq. (17), so that the total Hamiltonian of the few-body model is in a standard summation form [Eq. (18)].

• •

According to the symmetry of the system, choose a set of matrix basis to expand $\hat{H}^{\partial}$ [Eq. (19)]. The basis will determine which kind of spins will be used as the bath sites, and the expansion coefficients will be the coupling constants.

• •

Build the few-body experiment with several physical sites in the bulk and bath sites on the boundary (e.g., Fig. 1 or the inset of Fig. 4). The coupling constants in the bulk are the same as the target model, and the coupling constants on the boundary are given by the expansion coefficients of $\hat{H}^{\partial}$.

• •

Observe the properties of the bulk, which mimics the ground state of the infinite system.

V Summary and outlook

We propose an ab-initio TN approach that allows for accurate survey of the ground states of infinite many-body systems in higher dimensions by an effective few-body models embedded in an “entanglement bath”. On one hand, our scheme gives to birth to an flexible and efficient numeric algorithm for quantum lattice models. Our approach can directly access the thermodynamic limit by introducing the physical-bath interactions, which outperforms the finite-size methods such as ED and DMRG. The embedding idea allows for efficient and accurate simulations of infinite 3D quantum models, surpassing the existing TN methods. It is free from the “negative-sign” problem and can access to frustrated spin and fermionic models. It can accurately capture many-body features including entanglement, phase transitions and critical behaviors, thus it goes beyond the DFT-based approaches. It could be readily applied to other $(d\leq{3})$-dimensional systems and could be generalized to $(d\geq{4})$-dimensional models.

In practice, our numerical simulations show that with only a handful of sites, the few-body models can accurately capture the many-body features of the infinite systems. With less than 18 physical and 12 bath sites, the difference between our results and the state-of-the-art TN methods is less than $\mathcal{O}(10^{-3})$. For the spin models on simple cubic lattice, the properties of the quantum phase transitions in a magnetic field, including the phase transition point and critical exponents, are faithful captured by the few-body model containing only 8 physical and 24 bath sites.

On the other hand, the few-body Hamiltonian only contains local interactions among a handful of sites, it can be realized by, e.g., cold atoms or ions. It is possible to further improve the experiments by using the trick of synthetic gauge fields, where the higher spins, for instance, can be extended to lower spins in a synthetic dimension Celi et al. (2014). We suggest to investigate infinite many-body systems by realizing the predicted few-body Hamiltonian with cold atoms or ions. The many-body phenomenon are expected to be observed in the bulk. Furthermore, our work exhibits a new perspective of designing novel quantum devices CSW by utilizing the many-body properties that appears in the bulk of the few-body system, e.g., controlling the entanglement or quantum fluctuations by driving the system to approach or leave the critical region.

Acknowledgments

We acknowledge Leticia Tarruell, Ignacio Cirac, Emanuele Tirrito, Xi Chen, Nan Li and Jia Kong for enlightening discussions. CP appreciates ICFO (Spain) for the hospitality during her visit and is grateful to financial support from UCAS and ICFO. SJR is indebted to Nan Li for her kind help on illustrating the figures of lattices. CP thanks Yun-Hai Lan for modifying references format. This work was supported by ERC AdG OSYRIS (ERC-2013-AdG Grant No. 339106), the Spanish MINECO grants FOQUS (FIS2013-46768-P), FISICATEAMO (FIS2016-79508-P), and “Severo Ochoa” Programme (SEV-2015-0522), Catalan AGAUR SGR 874, Fundació Cellex, EU FETPRO QUIC, EQuaM (FP7/2007-2013 Grant No. 323714), and CERCA Programme / Generalitat de Catalunya. SJR acknowledges Fundació Catalunya - La Pedrera $\cdot$ Ignacio Cirac Program Chair. GS and CP were supported by the MOST of China (Grant No. 2013CB933401), the NSFC (Grant No. 11474279), and the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB07010100).

Appendix A Ab-initio optimization principle approach in one dimension

Understanding the 1D AOP Ran (2016); Tirrito et al. will be a lot of help to understand the higher-dimensional version. In this section, we present the 1D AOP scheme that does not have the Hermitian requirement Tirrito et al. . Let us take the following translationally invariant Hamiltonian as an example, which reads

[TABLE]

with $\hat{H}_{n,n+1}$ the two-body interaction.

With the Hamiltonian, the next step is to prepare the cell tensor that parameterize the self-consistent equation. This step is the same for the AOP approaches introduced in Sec. III. One firstly chooses a supercell that obeys the translational invariance, e.g. two adjacent sites. The bulk interaction is simply the coupling between these two spins, i.e. $\hat{H}^{B}(i,j)=\hat{H}(i,j)$. Then we define the two-body operator $\hat{F}^{\partial}(i,j)=I-\tau\hat{H}(i,j)$ as the shifted interaction on the boundary of the supercell and do the singular value decomposition as

[TABLE]

$\hat{F}^{L}(i)_{a}$ and $\hat{F}^{R}(j)_{a}$ are two sets of one-body operators (labeled by $a$) acting on the left and right spins of $\hat{F}^{\partial}(i,j)$, respectively.

The cell tensor is defined as the product of the (shifted) bulk Hamiltonian with $\hat{F}^{L}_{a}$ and $\hat{F}^{R}_{a}$ as

[TABLE]

with $\tilde{H}^{B}(i,j)=I-\tau\hat{H}^{B}(i,j)$.

Then, with two boundary tensors $|V^{[x]})$ (guessed or previously obtained in the last iteration) and the cell tensor $\hat{\mathcal{F}}(i,j)$, we define three eigenvalue equations (Fig. A1) as

[TABLE]

By solving the leading eigenvector of $\hat{\mathcal{H}}(i,j)$ given by Eq. (A4), we obtain the central tensor $|A(i,j)\rangle_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$. The central tensor can be considered as a state in the Hilbert space of the supercell labeled by two virtual indexes. Meanwhile, $|V^{[x]})$ is obtained as the (left) leading eigenvector of $M^{[x]}$ [Eqs. (A5) and (A6)]. $M^{[x]}$ is defined by the isometries $|\tilde{A}^{[x]}(i,j)\rangle$ that is obtained by the QR decomposition of $|A(i,j)\rangle$. For example for $x=1$, we have

[TABLE]

$|\tilde{A}^{[1]}(i,j)\rangle$ is orthogonal, satisfying

[TABLE]

These isometries play the role of the renormalization group flow in the standard DMRG White (1992). Similarly, $|V^{[x]})_{a_{x}\mu_{x}\mu_{x}^{\prime}}$ can be understood as a “state” defined in the space of the boundary index $a_{x}$ labeled by $\mu_{x}$ and $\mu_{x}^{\prime}$. One can see that these equations are parametrized by the solutions of others, and can be solved in an alternative way in practice.

In the language of TN, Eqs. (A4)-(A6) encode an infinite TN that represents the imaginary time evolution for the ground state (Fig. A2). Define a local scalar function shown in the top of Fig. A2. It is easy to see that this function is maximized at the fixed point of the self-consistent equations. Then, one can reconstruct an infinite MPO multiplied with an MPS and its conjugate by repetitively replacing $|V^{[x]})$ by $M^{[x]}|V^{[x]})$. Still, the scalar function is maximized, meaning the MPS is the dominant eigenstate of the MPO. One can again iteratively replace one MPS by the product of the MPS and the MPO to reconstruct the whole infinite TN. Note such a reconstruction can be understood in an opposite order: by going from the bottom to the top of Fig. A2, it actually gives a contraction scheme of the TN.

There are two important constraints to realize such a reconstruction. In the step from the local scalar function to MPO, we have $(V^{[x]}|V^{[x]})=1$ since it should be the eigenstate of $M^{[x]}$. In the second step from the MPO to the whole TN, we require that the MPS is normalized, which is actually a non-local constraint. In the original proposal of the 1D AOP Ran (2016), this constraint is turned to be local with some tricks under the assumption that all eigenvalue problems are Hermitian. In the generalized 1D version Tirrito et al. presented above, the MPS is normalized because of the orthogonality of $|\tilde{A}^{[x]}(i,j)\rangle$.

Now we explain the few-body Hamiltonian that mimics the ground state of the infinite 1D chain. By reviewing Eq. (A3), the matrix [Eq. (A4)] whose eigenstate gives $|A(i,j)\rangle_{\mu_{1}\mu_{2}}$ can be written as the product of three parts: one bulk term and two boundary terms. The bulk term contains simply the physical interactions of the original model. Same to Eq. (14) in higher-dimensional AOP, the boundary parts are defined as

[TABLE]

Then same to Eq. (18), the few-body Hamiltonian for the 1D simulator is obtained as

[TABLE]

where $\hat{H}^{\partial}(n,x)$ satisfies $\hat{\mathcal{H}}^{\partial}(n,x)=I-\tau\hat{H}^{\partial}(n,x)+O(\tau^{2})$, same to Eq. (17).

Comparing the 1D and higher-dimensional AOP versions, we can see many connections, including the idea of defining the eigenvalue equations, the constraints for the encoding/reconstruction process, and the emergence of the few-body Hamiltonians. The differences are also crucial. For 1D quantum systems, we can directly encode the imaginary-time-evolution TN. In higher dimensions, the tree approximation is introduced. The essential reason is to satisfy the second constraint, which is the normalization of the state ansatz (see Appendix B). The normalization of a standard TN state (e.g., on square lattice) requires an extra loop of TN encoding or contraction. The AOP algorithm to directly encode the ($D+1$)-dimensional imaginary-time-evolution TN for $D$-dimensional quantum system ($D\geq 2$) is still an open issue. This will lead to a general form of the few-body Hamiltonians (see Appendix E).

Appendix B State ansatz behind our approach in higher dimensions

At the first stage, the ansatz is an infinite tree PEPS that optimally approximates the ground state in the rank-1 sense Lathauwer et al. (2000); Ran et al. (2013). This can be seen from the tensor network (TN) encoded in the self-consistent eigenvalue equations. Starting from Eq. (6), one can substitute each of the boundary tensors $|V^{[x]})$’s by the contraction of the other three $|V^{[x]})$’s, $|\tilde{A}^{[x]})$, $(\tilde{A}^{[x]}|$ and the cell tensor $\hat{\mathcal{F}}$ according to Eqs. (7)-(10). We are using the fact that $|V^{[x]})$ is the eigenvector of $M^{[x]}$. By doing so repetitively, an infinite tree PEPS formed by $|A\rangle$ and $|\tilde{A}^{[x]})$’s can be grown to reach the thermodynamic limit. At the same time, the TN that gives ($I-\tau\hat{H}_{tree}$) appears, where $\hat{H}_{tree}$ is the Hamiltonian defined on the tree. The local interactions of $\hat{H}_{tree}$ are exactly the same with the original model as long as one only looks at a loop-free subsystem, thus $\hat{H}_{tree}$ provides a reasonable approximation. Such a tree PEPS minimizes the energy of $\hat{H}_{tree}$.

For better understanding the approximation of the state on, e.g., an infinite square lattice, we could “grow” the tree in such a way that it fills the whole square lattice. Inevitably, some $|V^{[x]})$’s on the boundary of the tree will gather at the same site. The tensor product of these $|V^{[x]})$’s in fact gives the optimal rank-1 approximation Lathauwer et al. (2000) of the tensor that forms the bulk of tree TN (translational invariant). Now, if one uses the full-rank tensor to replace its rank-1 version (the tensor product of four $|V^{[x]})$’s), one will have the TN of $I-\tau\hat{H}$ (with $H$ the target Hamiltonian on square lattice) instead of $I-\tau\hat{H}_{tree}$, and the tree PEPS becomes the one defined on the square lattice. Such a picture can be understood in the opposite manner: imaging that one has the “correct” TN defined on the square lattice, what we do is to replaced certain tensors by its rank-1 approximations to destruct all the loops of the TN. In this sense, the tree PEPS defined on the original lattice (not actually a Cayley tree or Bethe lattice Cayley (1878); Bethe (1935)) in stage one provides the optimal loop-free approximation of the ground state, where the loops are destructed by the rank-1 tensors. It would be very helpful to refer to the figures and the discussions in Ref. [Ran et al. (2013)] that are given considering TN contractions.

There are several issues we shall stress. Firstly, one will actually not do the above substitutions to reconstruct the tree PEPS. It is automatically encoded in the self-consistent equations. The “reconstruction picture” is proposed only to understand the ansatz behind the approach. Secondly, one may notice that the self-consistent equations proposed here are slightly different from those for the rank-1 decomposition of a single tensor Lathauwer et al. (2000). The reason is that in our case, the normalization of the PEPS should be considered when doing the rank-1 approximation. We here borrow the idea of iDMRG on the tree PEPS Lepetit et al. (2000); Nakatani and Chan (2013) to satisfy this constraint. The third issue is about the uniqueness of the reconstruction of the tree PEPS. Indeed, the contraction of three $|V^{[x]})$’s, $|\tilde{A}^{[x]})$, $(\tilde{A}^{[x]}|$ and $\hat{\mathcal{F}}$ to substitute $|V^{[x]})$ is not unique. However, it is unique when we require the presence of $|\tilde{A}^{[x]})$, $(\tilde{A}^{[x]}|$ and $\hat{\mathcal{F}}$, in order to recover the TN’s of $I-\tau\hat{H}$ as well as the tree PEPS. This is due to the uniqueness of the rank-1 decomposition, which is argued to be a concave problem Lathauwer et al. (2000).

Such a tree approximation is also closely related to the iPEPS algorithms called simple update Jiang et al. (2008); Ran et al. (2012); Braiorr-Orrs et al. (2016), where the infinite PEPS is updated by considering the local environment. After reaching the fixed point, the PEPS satisfies a set of self-consistent equations, which lead to a similar tree structure Ran et al. (2012). Even some long-rang effects are ignored, simple update are still quite accurate especially for gapped states.

Aimed at reducing the error of the tree approximation, the second stage of our approach is to construct the few-body Hamiltonian $\mathcal{\hat{H}}^{FB}$ on a larger cluster by reusing the bath obtained in the first stage, and then calculate the ground state of $\mathcal{\hat{H}}^{FB}$ with a finite-size algorithm. The ansatz behind can be considered as a generalized tree PEPS. In the center of the PEPS, the tensor contains all the physical sites inside the cluster, connected with several infinite tree brunches that are the same to those appearing in stage one. The bath sites carry the entanglement between the physical sites in the cluster and these infinite tree brunches.

Appendix C “Finite-loop” effects

Thanks to the infinite tree brunches in the PEPS ansatz, our algorithm does not suffer the conventional finite-size effect in the algorithms such as ED, QMC or DMRG. Thus, the effects from the finiteness of the cluster in the second stage are essentially different. In the first stage, the system size is already infinite because the bath encodes the information of an infinite tree in the eigenvalue equations. Only the loops beyond the supercell are destroyed in an optimal manner (rank-1 approximation of the tensors) Ran et al. (2013). In stage two, there will be no tree error inside the cluster since all interactions there are fully considered. If the cluster contains larger loops than the cell tensor used in stage one, the precision will be improved. On the other hand, there will be no improvement if one increases the size of the cluster without having larger loops. For this reason, the “finite-size effects” of AOP mean the errors caused by the finiteness of the considered loops.

Appendix D Computational cost

The motivation to use the tree approximation is its efficiency especially for 3D quantum models. The computational cost of the first stage is that of the generalized DMRG on an infinite tree PEPS Lepetit et al. (2000); Nakatani and Chan (2013), which roughly scales as $\mathcal{O}(d^{2N_{0}}D^{3z})$ with $d$ the dimension of the physical Hilbert space on one site, $N_{0}$ the number of physical sites in the supercell, $D$ the dimension of a virtual index and $z$ the coordination number of the lattice Not (d).

To solve the few-body Hamiltonian, the computational cost (leading term) with ED scales as $\mathcal{O}(d^{N}D^{N^{\partial}})$ ($N$ and $N^{\partial}$ the number of physical and bath sites, respectively), and that with DMRG scales as $\mathcal{O}[(N+N^{\partial})\max(d,D)^{3}\chi^{6}]$ ($\chi$ the bond dimension cut-off of DMRG). The cost is similar to solving a nearest-neighbor finite-size system that contains two kinds of sites, whose local Hilbert space is of dimension $d$ (physical) and $D$ (bath), respectively. Surely one can choose other algorithms to solve the few-body Hamiltonian in the second stage, such as QMC or finite PEPS algorithms Lubasch et al. (2014a); Zhao et al. (2016). Benefits from the fact that the few-body Hamiltonian is the product (or summation) of local couplings, the efficiency will be similar to that when applying to the standard (short-range) Hamiltonians. In addition, it is possible to update the bath simultaneously in stage two, and the computational cost would be approximately identical to the cluster update schemes of TN.

Appendix E General forms of few-body Hamiltonian

As discussed above, the dominant error comes from the destruction of the loops. As a consequence, the interactions between the bath and the physical sites are the tensor product of local terms

[TABLE]

It means that in the standard summation form, there are no bath-bath interactions (Fig. A3). The tree branches in the ground-state ansatz are not connected to each other from anywhere else than the central part.

One can adopt other TN algorithms such as the cluster or full update schemes Verstraete and Cirac ; Jordan et al. (2008); Wang and Verstraete ; Lubasch et al. (2014b); Nishino and Okunishi (1996); Orús and Vidal (2009); Orús (2012); Corboz (2016); Vanderstraeten et al. (2016); Xie et al. (2009, 2012) to obtain the physical-bath interactions. Then the Hamiltonian will not simply be the tensor product, but generally given by

[TABLE]

Then the bath-bath interactions will appear in the standard summation form. See the illustrations of three possible situations in Fig. A3. The extra summations will lead to another (similar) PEPS ansatz beyond the one with tree branches, which should better mimics the infinite environment. However, the computational cost with the currently known methods will become much more sensitive to the coordination number and the dimensionality of the model, making the 3D ground states extremely difficult to access.

Appendix F Discussions about imaginary-time evolution picture and criticality in higher dimensions

The idea of approximating an infinite Hamiltonian with a finite effective one has been proposed for the time evolution of 1D quantum systems Phien et al. (2012). An important difference in our work is that the “evolution” of the finite effective model is constructed not from a new $\hat{H}$ but with a shift ($I-\tau\hat{H}$) that is in fact the imaginary-time evolution operator. It brings several operational advantages for simulating the ground states, in particular, of higher-dimensional systems. The triangular structure of the Hamiltonian is avoided here, thus the eigenvalue equations for the boundary tensors have stable solutions and the entanglement bath is well-defined. The few-body Hamiltonians with the bath of higher-dimensional systems can be easily constructed as the summations of local terms.

Our scheme makes it possible to adopt the ($1+1$)-D scaling theories for characterizing criticality Vidal et al. (2003) to higher-dimensional models. It is known that any TN algorithms, essentially, cannot give directly a divergent correlation length at the critical point. For 1D quantum systems, it has been shown that at the critical point, any MPS with a finite bond dimension is gapped and possesses a finite correlation length $\xi$ Vidal et al. (2003) satisfying

[TABLE]

with $\mathcal{D}$ the bond dimension of the MPS and $\kappa$ its scaling exponent. One can see that with a finite $\mathcal{D}$, $\xi$ is always finite, and the information of the criticality is in hidden the algebraical scaling behavior when $\mathcal{D}$ increases. For the scaling of magnetic field $h$ near the critical point, the algebraic behavior of $\xi$ versus $h$ can still survive, however, the value of the exponent might be inaccurate.

For a 2D PEPS, one has to compute the contraction of a 2D TN (e.g., by iTEBD with MPS) to get its correlations using finite dimension cut-offs, and thus the results will still be finite Ran et al. (2017). To tackle this difficulty, it has been proposed that the divergence of the correlations can be studied by the scaling of the bond dimension of the MPS, from which the central charge of the conformal field theory to characterize the criticality can be accurately obtained Vidal et al. (2003); Ran et al. (2017).

In our approach, the dynamic correlation length of the ground state $\xi$ is given by the correlation length of an infinite MPS formed by $|V^{[x]})$ in the imaginary time direction, written as $|\tilde{\psi}\rangle=\sum_{\{\mu\}}\prod|V^{[x]})_{a_{x}\mu_{x}\mu_{x}^{\prime}}$. Such an MPS (dubbed as time MPS) is quasi-continuous (discretized up to the Trotter step $\tau\to 0$). Let us explain how to get $\xi$ in the AOP approach. In higher dimensions, the scheme is similar.

The dynamic correlation function of the ground state is defined as $\langle\Phi|\hat{S}e^{-\beta\hat{H}}\hat{S}|\Phi\rangle/e^{-\beta E}-\langle\Phi|\hat{S}|\Phi\rangle^{2}$ with $|\Phi\rangle$ the ground state and $E$ the ground-state energy. In our framework, it is the contraction of a TN, where the two operators are put in the same column. Thanks to the encoding scheme, such a contraction becomes the contraction of a tensor stripe (Fig. A4). This stripe is the product of $\mathcal{\hat{H}}^{FB}$’s [Eq. (15), also called the transfer matrix] with the two operators in between. The dynamic correlation length is defined as the exponent of the correlation length. One can see that such an exponent is obtained simply by

[TABLE]

with $\Lambda_{0}$ and $\Lambda_{1}$ the two largest eigenvalues of $\mathcal{\hat{H}}^{FB}$.

An advantage of the dynamic correlation properties is that we find much less finite-loop or finite-dimension-cutoff effects than the spatial correlations. This is also supported by a recent DMRG work DMRGscaling , where the finite-size effects are found to be much smaller for the dynamic correlations. Meanwhile, a finite dimensional matrix cannot give a critical spectrum. It means one cannot directly obtain a divergent correlation length at the critical point, and a scaling of the dimension would be necessary to identify the criticality. How to do such kind of scalings for 2D and 3D states is still an open question.

Appendix G Relations to other algorithms

By taking certain limits of the computational parameters, the relations among our approach and other algorithms are illustrated in Fig. A5. The simplest situation is to take the dimension of the bath sites $\rm dim(\mu_{x})=\rm dim(\mu_{x}^{\prime})=1$, and then $\hat{\mathcal{H}}^{\partial}$ can be written as a linear combination of spin operators (and identity). Thus in this case, $|V^{[x]})$ simply plays the role of a mean field. If one only uses the bath calculation of the first stage to obtain the ground-state properties, the algorithm will be reduced to the tree DMRG Lepetit et al. (2000); Nakatani and Chan (2013). If one takes the minimal supercell with $D=1$ in stage one, the entanglement bath will be reduced to a magnetic mean field. By choosing a large cluster, the DMRG White (1992) simulation in stage two becomes equivalent to the standard DMRG for solving the cluster in a mean field. If one uses $D=1$ and chooses a supercell of a tolerably large size in the first stage without entering stage two, or if one chooses a small cluster with $D=1$ in stage one and uses ED in stage two to solve the few-body Hamiltonian with a tolerably large cluster, our approach will become the ED on the corresponding finite system in a mean field. By taking proper supercell, cluster, algorithms and computational parameters, our approach outperforms others.

Appendix H Generalization to ($d\geq 4$) dimensions

Benefiting from its flexibility, it is possible to generalize our approach to even $(d\geq 4)$-dimensional quantum models. The main problem to be tackled is the computational cost. In the second stage by using DMRG for example, the cost increases polynomially with the size of the cluster, thus also polynomially with the dimensionality $d$. In the first stage with tree DMRG, the cost increase exponentially with $d$, which makes the simulations for higher-dimensional models extremely expensive. Luckily, the main task here is to solve ($2d+1$) number of self-consistent eigenvalue equations, say five [Eq. (6)-(10)] for 2D, seven for 3D and nine for 4D quantum systems. One way to lower the cost from exponential to polynomial expenses is to use a finite algorithm such as DMRG to solve each eigenvalue problem. It is certain that the stability and efficiency have to be tested.

Appendix I Open issues

Several following-up issues are to be further investigated. The flexibility allows for possible incorporating with other methods. For example, the TN techniques with symmetries McCulloch and Gulàcsi (2002); McCulloch (2007); Weichselbaum (2012) can be introduced to lower the computational cost so that much larger clusters can be reached in the second stage. Besides the tree DMRG Lepetit et al. (2000); Nakatani and Chan (2013), the other TN optimization schemes such as TN variational techniques Phien et al. (2015); Vanderstraeten et al. (2016); Corboz (2016) and tensor renormalization group algorithms Levin and Nave (2007); Orús and Vidal (2009); Orús (2012); Xie et al. (2009, 2012); Yang et al. (2017) can be adapted when the cost is tolerable. The finite-size scaling of the cluster should be explored. Our approach could also be readily generalized to higher-dimensional bosonic and fermionic lattice models. The entanglement embedding idea with the physical-bath Hamiltonian proposed here can be adopted to develop novel algorithms for infinite systems by hybridizing with other methods such as QMC, finite or tree PEPS algorithms Li et al. (2012); Lubasch et al. (2014a); Zhao et al. (2016), or the approaches in material sciences and quantum chemistry, such as DFT Sherrill (2010) and DMET Knizia and Chan (2012).