Entanglement Hamiltonian of Interacting Systems: Local Temperature Approximation and Beyond

We investigate the second quantization form of the entanglement Hamiltonian (EH) of various subregions for the ground-state of several interacting lattice fermions and spin models. The relation between the EH and the model Hamiltonian itself is an unsolved problem for the ground-state of generic local Hamiltonians. In this letter, we demonstrate that the EH is practically local and its dominant components are related to the terms present in the model Hamiltonian up to a smooth spatially varying temperature even for (a) discrete lattice systems, (b) systems with no emergent conformal or Lorentz symmetry, and (c) for subsystems with non-flat boundaries, up to relatively strong interactions. We show that the mentioned local temperature at a given point decays inversely proportional to its distance from the boundary between the subsystem and the environment. We find the subdominant terms in the EH as well and show that they are severely suppressed away from the boundaries of subsystem and are relatively small near them.

We investigate the second quantization form of the entanglement Hamiltonian (EH) of various subregions for the ground-state of several interacting lattice fermions and spin models. The relation between the EH and the model Hamiltonian itself is an unsolved problem for the ground-state of generic local Hamiltonians. In this letter, we demonstrate that the EH is practically local and its dominant components are related to the terms present in the model Hamiltonian up to a smooth spatially varying temperature even for (a) discrete lattice systems, (b) systems with no emergent conformal or Lorentz symmetry, and (c) for subsystems with non-flat boundaries, up to relatively strong interactions. We show that the mentioned local temperature at a given point decays inversely proportional to its distance from the boundary between the subsystem and the environment. We find the subdominant terms in the EH as well and show that they are severely suppressed away from the boundaries of subsystem and are relatively small near them.
Introduction.-Entanglement is a unique feature of quantum mechanics and serves as an essential tool in quantum information, quantum gravity, identification of topological order, quantum phase transition, etc [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. The entanglement Hamiltonian (EH) associated with a subregion A embedded in a manifold M = A ∪ B is defined as ρ A = e −K A . Here ρ A = Tr B ρ M denotes the reduced density matrix (RDM) of A, where ρ M represents the total density matrix. One important question that arises from this definition is the relation between K A and H A , the Hamiltonian terms with support only in region A. In fact, this problem dates back to the 19th century. A cornerstone of the classical statistical mechanics is that a subsystem A at thermal equilibrium with its environment (B) is described by a thermal ensemble with K A = H A /T 0 (k B = = c = 1) where T 0 is a uniform and position independent temperature. Furthermore, the eigenstate thermalization hypothesis conjectures that the RDM of highly excited quantum states will look thermal, again with K A = H A /T 0 , where the uniform temperature T 0 in this case is dictated by the energy density [18]. In this letter, we revisit this fundamental problem and using the densitymatrix-renormalization-group (DMRG) approach we obtain the second quantization form of K A for a number of interacting model Hamiltonians and for a variety of boundary shapes and conditions. Comparing the components of K A and H A , we demonstrate that the above mentioned statements are not quite accurate and for the ground-state of local Hamiltonians, K A is indeed well-approximated by a local and non-uniform temperature rather than a uniform one.
The theoretical form of K A is known only for a limited class of continuum models with conformal symmetry (or Lorentz symmetry at zero temperature) and only for certain geometries of A (e.g., half-space or ball geometry). It is known that under these conditions: (i) K A is local, (ii) the EH density is related to the Hamiltonian density via a smooth local temperature, namely K A = x∈A d d x K A (x) = * The two first authors have contributed equally to this work. † Corresponding author. Email address: vaezi@sharif.edu T (x) , and (iii) T (x) approaches T 0 , the equilibrium temperature of the entire system, far away from ∂A (the boundary of A) and grows as v 2πr(x) at distance r near ∂A. Here, v is the group velocity of low energy excitations [1,[19][20][21][22][23][24][25][26]. We refer to these findings as the local temperature approximation (LTA) [25,26].
The LTA can be justified using the following intuitive argument. In thermal systems, the entropy density is proportional to their temperatures. On the other hand, for groundstates, instead of the thermal entropy, we deal with the entanglement entropy which is not precisely an extensive property. Nevertheless, we can still consider and gauge the contribution of individual degrees of freedom residing inside A to the overall entanglement entropy between A and B, S A . Indeed, quantum mutual information can be one candidate to quantify such local contributions. Due to the decay of quantum mutual information with distance for the ground-state of local Hamiltonians, the degrees of freedom that live near ∂A, are more entangled with those residing at B than more distant ones. Accordingly, we can assign an effective quantum local temperature to different subregions of A proportional to their contributions to S A , which as just discussed must diminish away from ∂A.
Numerically, the EH of free fermions and free bosons can be evaluated easily [28][29][30]. However, for interacting models, it becomes highly nontrivial and challenging. Recently, several studies have analyzed the EH of quasi-one-dimensional conformal invariant or integrable models [31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47]. Nonetheless, we are still lacking a systematic derivation of the EH for larger and generic interacting systems and for various boundary geometries. In this letter, we address this problem and introduce a DMRG-based algorithm that enables us to extract the EH for a broader spectrum of problems. J 1 − J 2 Heisenberg model.-Let us first discuss the form of K A for the J 1 − J 2 Heisenberg model (J 1 = 1) with the following Hamiltonian on the square lattice:  Fig. 3 are plotted respectively.
The above Hamiltonian respects a SU(2) symmetry. Hence, K A must respect SU(2) symmetry as well and thus expanded as follows: Note that there is no restriction on i = (i x , i y ) and j = (j x , j y ) except that both must belong to A. In our study, we have dropped higher order terms since the retained terms already yield satisfactory results. The DMRG technique is based on identifying the most relevant basis states of the Hilbert space [48]. Then we truncate the Hilbert space and discard the less relevant states. The number of kept states which controls the accuracy of DMRG is called the bond dimension, χ, and its default value equals 2 10 throughout this letter. The procedure of finding the truncation operators consecutively involves the computation and diagonalization of the RDM at every step of DMRG and for different subsystem sizes. Hence, ρ A is a natural byproduct of DMRG method and is available at every step. Moreover, every operator component of K A (e.g., S i .S j ) has a matrix representation in DMRG, albeit in the truncated subspace. The remaining task is to adjust the EH's couplings, g J,ij , to bring our (simplified) guess for K A , which we will denote as K A , close enough to the matrix representation of the RDM in the truncated Hilbert space achieved via DMRG: K A = − log ρ A . To this end, we need to define an appropriate cost function as a measure of the distance between K A and K A . In our investigations, we mainly utilized the Hilbert-Schmidt distance between the Green's functions, namely is the Green's function matrix achieved by DMRG for A and G A denotes its counterpart evaluated using the trial RDM, ρ A = exp − K A [49]. This cost function yields more reliable and reasonably robust results (against changing χ) in the truncated Hilbert space than other candidates, e.g., the quantum relative entropy between ρ A and ρ A . We find the latter to overfit to numerical noises, e.g., the truncation and computer's roundoff errors (see the supplemental material). In the optimization procedure, we initialized g J,ij based on general expectations from LTA, e.g., the locality of g J,ij and its linear dependence on x ij (the minimum distance between ij = i+j 2 and ∂A). Then, we employed the gradient descent algorithm and let the cost function to decide the optimum choice for g J,ij (see the supplemental material).
We first focus on J 2 = 0 Heisenberg model which is unfrustrated and is known to host a Néel order on the square lattice [50,51]. Thus, its ground-state is a symmetry broken phase with gapless Goldstone modes and does not respect the full conformal symmetry (e.g., the translational and rotational symmetries are broken). For this model, the system is always subject to the PBC along y axis.
• As the first example, we study the manifold and subsystem A depicted in Fig. 1-a, where an OBC is imposed along x. In this case, ∂A is flat and its locus is given by x b = 6 + 1/2, the line which splits columns 6 and 7. The optimum couplings, g J,ij , which reproduce the DMRG's Green's functions (with less than 0.1 % error), are plotted in Fig. 2. x and y (more precisely, β J,x (i x + 1/2) := g J,i,i+x and β J,y (i x ) := g J,i,i+ŷ ) which are independent of i y due to the y-axis translation preserving shape of A. Indeed, β J 's are the inverse local temperature profiles. As Fig. 2a suggests, β J,x and β J,y follow the same profile, albeit if we shift the argument of β x by half of the lattice spacing. This shift is due to the fact that for g J,x , the start and end points are located at different positions along x, while for g J,y , the two points have identical x values. In the supplemental material, we demonstrate the robustness of β J,x and β J,y versus χ. In Fig. 2b, we have plotted g J,xy (i x + 1/2) := g J,i,i+x+ŷ , g J,yy (i x ) := g J,i,i+2ŷ , and g J,xx (i x + 1) := g J,i,i+2x . Their values are negligible everywhere and they all die off quickly away from ∂A. These imply the locality of K A for the Heisenberg model when ∂A is flat.
• We now consider the same conditions as above, but this time with a curved ∂A as shown in Fig. 1-b. In this case, β J,x and β J,y will depend on both i x and i y . In Fig. 3, we have plotted β J,x and β J,y for two different rows marked by orange and green lines. Interestingly, β J,x and β J,y profiles display a somewhat smooth curve satisfying our expectations from LTA. The position dependence of the inverse temperature profile is more complicated in this problem, since the distance between ij = (i + j) /2 and ∂A depends on both its x and y components. In Figs. 3b and 3d the second and third neighbor couplings are plotted for the above mentioned rows.  We see that for curved boundaries between A and B, K A remains local everywhere, except close to ∂A where we observe additional terms though relatively small and subdominant.
• Now, we consider the geometry illustrated in Fig. 1-c. Since the PBC is imposed on M along both x and y directions, we have chosen N y = 4 to ensure χ = 2 10 is sufficient for DMRG's convergence. As we see in Fig. 1-c, ∂A is described  by two surfaces, one of them separates columns 6 and 7 and the other one lies between the first and last columns. As a result, LTA predicts that β J,x and β J,y must follow a parabolic form and vanish near both boundary surfaces. Fig. 4, shows our numerical results for the nearest as well as further neighbor couplings, both consistent with LTA.
• Let us now turn to the frustrated Heisenberg model with J 2 = 0.6 whose true ground-state is not well-understood, though it is conjectured to be a spin liquid phase with no classical spin order and algebraically decaying spin-spin correlations [52]. For this model, we consider the geometry depicted in Fig. 1-d. The ground-state is expected to be more entangled when J 2 /J 1 ∼ O(1). Hence, we consider N y = 4 (again N x = 24) to ensure that the ground-state is achieved reliably via χ = 2 10 in DMRG. Since the Hamiltonian contains next nearest neighbor (NNN) couplings, we expect significant values for the NNN in g J,ij as well. In Fig. 5a, β J,x , β J,y , and also β J,xy (i x + 1/2) := 1 J2 g J,i,i+x+ŷ are plotted. Again, we see that all these β J 's follow the same curve. Furthermore, Fig. 5b verifies the locality of K A everywhere except at ∂A where g J,yy is about 29% (28%) of g J,y (g J,xy ) at that location.
Hubbard model.-Here, we discuss the second quantization LTA corrections are small (compared to the leading terms) everywhere, particularly away from ∂A. The particle-hole symmetry dictates the couplings in (b) to vanish. However, due to the finite truncation error of DMRG at χ = 2 10 , we obtain nonzero, though negligible values.
form of K A for the Hubbard model on the square lattice, whose Hamiltonian is: where n i,σ = c † i,σ c i,σ . In this letter, we consider t 1 = 1, and U = 4. The above Hamiltonian enjoys a U (1) × SU (2) symmetry for generic fillings. Accordingly, K A must be expanded as follows: where n i = n i,↑ + n i,↓ denotes the total electron number on site i, and S x,y,z the three components of the spin operator at i. Again, we have discarded higher order terms as we attain satisfactory results with the above structure. In the following, we consider both doped and undoped Hubbard models. In this section, we consider the geometry shown in Fig. 1-d. • Let us start with the half-filling case. The ground-state on the square lattice is described by a Néel anti-ferromagnetic spin order and the charge/Mott gap opens up at moderate values of U [53,54]. Fig. 6, summarizes our results for the optimum couplings of the EH. Here, motivated by LTA, we define the following inverse temperatures: β t,x (i x + 1/2) := g t,i,i+x , β t,y (i x ) := g t,i,i+ŷ , and β U (i x ) := 1 U g U,i . According to Fig. 6, they all fairly follow an identical curve, albeit by considering the previously discussed 1/2 shift in the argument of β t,x . In Figs. 6b-d, we have shown the terms beyond LTA and again the locality of couplings is confirmed. Only near ∂A, the additional terms are non-negligible. Although, g V,ij is insignificant everywhere, g J,ij has decent values near ∂A, yet inferior to those of g t, ij and g U,i .
• We now study the Hubbard model at p = 1/8 doping which is expected to have a stripe order and some tendency towards superconductivity [55][56][57][58]. The system is not expected to exhibit the Lorentz invariance or conformal symmetry for these symmetry breaking phases. At finite doping, we need to define β µ (i x ) := 1 µ g µ,i as well (for the current example: µ ≈ −0.92t 1 ). As Fig. 7 implies, various β's follow the expected trend and the locality of couplings is again verified, albeit with growing corrections close to ∂A.
Summary.-Our DMRG-based algorithm allowed us to access the second quantization form of the EH for several models and subsystem shapes. We showed that the EH is local and its dominant components are related to those of the Hamiltonian itself (more specifically the stress-energy tensor) up to a single smooth local (inverse) temperature profile and con-firmed LTA. We studied the terms beyond LTA and demonstrated they are infinitesimal far away from ∂A and relatively small near it. In the supplemental material, we have provided more evidences which further corroborate our main findings. To our knowledge, the validity of LTA for the ground-state of local Hamiltonians for generic models that do not satisfy conformal algebra or even those with conformal symmetry but non-flat ∂A is an unsolved problem despite active research. Our results suggest that LTA is perhaps a legitimate assumption and applicable to a broader class of problems.
Our findings pave the way for several applications of LTA. For instance, it can be shown that LTA can practically solve the long-standing sign problem in quantum Monte Carlo and enable us to extract the ground-state properties of some unsolved interacting models. Furthermore, LTA can be employed to enhance the performance and increase the accuracy of the DMRG technique. It can also be used to recover the entire spectrum and eigenstates of an unknown Hamiltonian by having access to its reduced density matrix (or correlation functions) associated with a rather small subregion of that system [27].

SUPPLEMENTAL MATERIAL
In this supplemental material, we will delve into the details of our algorithm and discuss the advantages and disadvantages of a number of cost function candidates along with their pairwise comparison, and present more results on the entanglement Hamiltonian (EH).

A. Entanglement Hamiltonian in the truncated Hilbert space
In general, the entanglement Hamiltonian (EH) associated with subsystem A which is defined as K A := − log ρ A , can be expanded in terms of a complete basis of operators (not necessarily local) as follows: (the tensor product of Pauli matrices, σ a , a = x, y, z with σ 0 = 1, can generate a basis for all possible operators): If we are given the reduced density matrix, ρ A , we can compute its logarithm (which requires a lot of considerations and special care when performed numerically) and achieve K A (up to computer's round-off error). Having K A available, we can easily find the expansion coefficients, g α , via the following relation: where In the exact diagonalization (ED) method M α,β ∝ δ α,β .
Therefore, g α ∝ Tr A K AÔ † α . Consequently, we do not need to consider other operators if we are interested in reading the coefficient for a specific α.
On the other hand, in the density-matrix-renormalizationgroup (DMRG) algorithm, instead of ρ A , andÔ α , we have to deal with ρ A , andÔ A which are their counterparts in the truncated Hilbert space and are defined as: and similarly for other operators, where T A denotes the truncation (a.k.a. projection) operator. As before, we define K A = − log ρ A . Because of the numerous truncations involved in DMRG, the situation is now more complicated for a few reasons: (i) The matrix M is not diagonal, neither sparse. Due to consecutive truncations inherent to the DMRG method, most of operators have non-vanishing overlaps. Therefore, we must consider all possible operators, including highly non-local ones such as string or brane operators. (ii) M can be singular and have zero eigenvalues. As a result, it might not be invertible. (iii) The above method applied to DMRG is very sensitive to various sources of numerical noises and errors, such as the truncation, as well as the round-off error. Besides the possibility of singular M , another main difficulty of applying the above algorithm to DMRG is the annoying part which requires taking all possible operators into consideration. Below, we easily demonstrate that if K A contains a few relevant and dominant terms, then K A contains exactly the same couplings and structure. To this end, recall that T A is achieved upon concatenating the dominant eigenvectors of ρ A . Therefore, Therefore, Accordingly, hence, assuming K A = α g αÔα :

B. Algorithm and cost function
In the above mentioned method, for ED, we can ignore insignificant couplings since M is diagonal. Nevertheless, when we apply this method to DMRG, we have to retain all terms, no matter how infinitesimal they are due to the complex form of M . Therefore, we must come up with a better algorithm to find g α without having to consider all irrelevant terms. For that purpose, we must consider a valid cost function. From our physical intuitions and expectations, we can think of the following three choices (as of now, we drop the overline sign and keep in mind that all operators are defined in the truncated Hilbert space): • Hilbert-Schmidt distance between the Green's functions (GFs): Here ρ A denotes the reduced density matrix (RDM) achieved via DMRG for the desired subsystem and ρ A denotes the one by combining the basis operators (O α ) with g α coefficients that are yet to determine. The basis of this method is that the RDM contains all the information about the equal time correlation functions within the subsystem. Thus, if we find a RDM which recovers all the correlation functions correctly, it must be identical to the actual one.
We would like to emphasize that in DMRG, due to finite truncation error, the RDM yields more reliable results for the expectation value of simple operators (e.g., twopoint correlation functions for short and intermediate distances) and becomes less reliable for more complex operators or at long distances. Therefore, to avoid overfitting to numerical errors, instead of considering all basis operators in the evaluation of ∆ 1 , we only consider the most physically relevant operators, i.e., simple operators motivated by symmetry considerations, etc. For example, for the Hubbard model, we consider the following components first: and similarly forG. We then evaluate η a = Tr G A,a − G A,a 2 and by combining them: The exact values of w t , w µ , · · · are not crucial as long as they all have the same order of magnitude. Nonetheless, in most computations, we choose w t = w U = w µ = w J = w V = 1 for the weights.
• Quantum relative entropy (QRE) of the two reduced density matrices: In this method (which is closely related to the next cost function), we try to tune the couplings such that ρ A 's matrix form becomes very close to ρ A 's. This cost function converges significantly fast, in both ED and DMRG method. In ED where we do not have to deal with truncation errors, QRE is the superior cost function and achieves correct results. However, for DMRG, (like ∆ 3 below), it also suffers from overfitting to numerical errors, i.e., those parts of ρ A which will change upon increasing the bond dimension of DMRG, χ (i.e., the number of retained basis states of the Hilbert space). It is these matrix elements which are responsible for the issues related to the expectation value or n−point correlation functions of complex operators explained above. When χ is large enough (e.g., when the truncation error becomes less than 10 −10 ), it yields results consistent with ∆ 1 's.
• Hilbert-Schmidt distance between the two reduced density matrices: in this method we try to tune the couplings such that ρ A 's matrix form becomes very close to ρ A 's. This cost function is slowly converging even for the ED where no truncation is involved. Moreover, for DMRG, similar to ∆ 2 , it suffers from overfitting to numerical errors

FIG. 8:
We study the Hubbard model on this ladder for U = 4, t ⊥ = 0.5, and t = 1 at half filling (µ = 0). Subsystem A, whose EH is desired, is denoted by blue sites. g U,0 versus log 2 χ. As this plot clearly indicates, ∆ 1 is the most reliable cost function for DMRG and exhibits least variations. and its results are sensitive to the bond dimension, especially for small values of χ.
In section D of this supplemental material, we compare the results achieved via all three cost functions for the two-leg ladder Heisenberg and Hubbard models for several bond dimensions. Our results suggest that for large bond dimensions, all three methods yield consistent outcomes. However, for relatively small bond dimensions, it is ∆ 1 which performs better and results in couplings which are more consistent with the results of larger bond dimensions.
C. Local temperature ansatz and the initial guess for couplings Now, let us assume we study the following Hamiltonian: FIG. 10: We study the Hubbard model on this ladder for U = 4, t ⊥ = 2, and t = 1 at half filling (µ = 0). Subsystem A, whose EH is desired, is denoted by blue sites.  Fig. 10, obtained via applying GF distance (∆ 1 cost function) and for χ = 2 11 , and χ = 2 9 . (a) Various β profiles for ∆ 1 cost function for χ = 2 11 (see the main text for their definitions). (b) Most significant corrections to LTA corrections, g J for ∆ 1 cost function for χ = 2 11 . The second and third neighbors' corrections to g t are negligible due to the particle-hole symmetry. Moreover, we found g V to be irrelevant as well, and that is why they are absent in this and the following two figures. (c-d) Same as (a-b) but for χ = 2 9 .
where, due to the locality of the Hamiltonian, only certain J α 's are nonzero. We are interested in finding the second quantization form of the EH expanded as follows: Here, due to the renormalization procedure involved in tracing the degrees of freedom outside A, g α 's can be viewed as our running coupling constants which J α has flown to. Thus, in principle, any g α consistent with symmetry considerations emerge. In practice, only a small set of them will be relevant and non-negligible. In our algorithm, we are trying to find g α numerically, assuming (a subset of relevant) correlation functions are known. In our optimization algorithm, we initialized the coupling constants of the EH, g α , using LTA's ideal form. In LTA, the EH is local and its coupling constants, g α 's, are nonzero only when the corresponding couplings of the Hamiltonian (UV theory), J α 's, are nonzero. Another task in LTA is to assign a position to each operator. For simple two-point operators (such as S i .S j in the Heisenberg model, or c † i,σ c j,σ in the Hubbard  FIG. 14: We study the Heisenberg model on this ladder for J ⊥ = 0.5, and J = 1. Subsystem A, whose EH is desired, is denoted by blue sites. model), the position is defined as the average position of its components, namely ij = i+j 2 . Next, we must compute the minimum distance (geodesics) between ij and the boundary separating A and its environment, B. Let us call this minimum distance, x ij . Finally, at zero temperature (for ground-states) and for the open boundary condition (OBC), LTA attributes the following form to g α (x ij ) [1]: where is the maximum value of x ij (i.e., the linear dimension of A normal to ∂A), and v is the group velocity of low energy excitations (quasi-particles) and is model-dependent. In our algorithm, besides v, we also treated as a variational parameter. We first optimized and tuned v, and . Then, we took the optimized form of local g α (associated with v * , and * ), and using the gradient descent algorithm we optimized our cost function. We allowed all relevant couplings, including distant neighbors and non-local terms (which were absent in the system's Hamiltonian) as well as the initialized local terms to vary and deviate from their initial point. Therefore, we have not imposed locality in our procedure, although it finally emerged naturally as the optimum solution (except at the boundary of A with B, where farther neighbors became more pronounced).
Similarly, for the periodic boundary condition (PBC) at   FIG. 16: We study the Heisenberg model on this ladder for J ⊥ = 2, and J = 1. Subsystem A, whose EH is desired, is denoted by blue sites. T = 0, LTA assigns the following value to g α (x ij ) [1]: where L is size of the entire system (M ) in the direction normal to ∂A.

D. A detailed comparison between the performance of ∆1, ∆2
and ∆3 cost functions Here, we compare the EH's coefficients obtained by utilizing all three cost functions for the ladder geometry and for the Hubbard and Heisenberg models.
We first consider the undoped Hubbard model on a ladder geometry (in which U = 4, t ⊥ = 0.5, t = 1) depicted in Fig. 8. This geometry results in a highly entangled subsystem, indeed a volume law entanglement entropy. We apply all three cost functions to this system for the following five different bond dimensions: χ = 2 5 , 2 7 , 2 9 , 2 11 . The EH for this case is translationally invariant, namely g τ,ix,jx = g τ,dxij (τ = t, U, J, V ), where dx ij := j x − i x . In this section, the translational symmetry is imposed on the couplings explicitly. We first compare the EH's couplings achieved via applying the GF distance (∆ 1 ), QRE (∆ 2 ), and the RDM distance (∆ 3 ) for χ = 2 11 (see Fig. 9). With this bond dimension, we can nearly probe the ground-state properties. The coefficients of local terms in the EH are almost consistent in these three methods. On the other hand, we know that ideally the groundstate must exhibit particle-hole symmetry. Although χ = 2 11 is still insufficient for true convergence in DMRG for such a highly entangled state (χ = 2 12 seems to be enough), ∆ 1 's results reflect the particle-hole symmetry (e.g., the (renormalized) second neighbor hoppings are infinitesimal), while for those of ∆ 2 and ∆ 3 the particle-hole symmetry is visibly violated due to overfitting issues mentioned previously. Furthermore, a previous quantum Monte Carlo based study of a similar situation [2] indicated the irrelevance of g V couplings which is consistent with ∆ 1 's estimations. Additionally, perturbative studies of the EH indicate an oscillating spin-spin couplings [3][4][5] (though subdominant to the renormalized onsite interaction) which agrees well with our results via minimizing GF distance (∆ 1 ), while those of ∆ 2 and ∆ 3 exhibit deviations in addition to their overestimation for the spin-spin couplings. In Fig. 9d, we plot the normalized onsite interac- for all four bond dimensions considered in our investigations. Again, as we see in Fig. 9d, the results of the GF distance (∆ 1 ) are more robust and less sensitive to χ, despite several orders of magnitude change in χ, while those of the QRE and RDM distance display stronger fluctuations. Now, we turn to the geometry shown in Fig. 10 (where U = 4, t ⊥ = 2, t = 1) and present our results for all three cost functions in Figs. 11, 12, and 13. Here, we have defined the following inverse local temperatures: β t,x (i x + 1/2) := g t,i,i+x , β t,y (i x ) := g t,i,i+ŷ , and β U (i x ) := 1 U g U,i . For this problem, due to the PBC imposed along x direction, we found out that even with χ = 2 11 , there is still some room for DMRG to converge to the true ground-state. As a result, we still see some minor discrepancy among the results of the three methods for χ = 2 11 for subdominant and correction terms beyond LTA (though they yield highly correlated results). Nonetheless, the local terms (i.e., dominant couplings) are reasonably consistent. We have also plotted the results of χ = 2 9 for all three methods and again, ∆ 1 's results proved to be more robust and ∆ 2 and ∆ 3 's less stable. Thus, in the presence of truncation errors, we can trust the results of the GF distance more than those of the other two candidates for the cost function.
For the sake of completeness, we have also explored the robustness and the accuracy of the above three cost function candidates for the Heisenberg model, again on a ladder geometry. To this end, we first studied the geometry shown in Fig. 14 (where J ⊥ = 0.5, J = 1), and presented its results in Fig. 15. Similar to the Hubbard model case, leads to a highly entangled ground-state. Likewise, we expect translationally invariant couplings, namely g J,ix,jx = g J,dxij . We have presented g J,dx for χ = 2 5 , 2 7 , 2 9 , 2 11 for the GF and QRE cost functions (the RDM cost function yields results similar to that of the QRE). In this case as well, the GF distance turns out to be the most stable one.
Finally, we studied the Heisenberg model on the geometry illustrated in Fig. 16 (in which J ⊥ = 0.5, J = 1). Their results are presented in Figs. 17 and 18 for the GF distance and QRE, respectively. For this problem, we indeed achieved the true ground-state using χ = 2 11 . Therefore, all cost functions must achieve the same couplings. On the other hand, for smaller bond dimensions, e.g., χ = 2 7 , ∆ 1 achieves more accurate results (relative to χ = 2 11 ) than the remaining cost functions.