# Dimension reduction with mode transformations: Simulating
two-dimensional

fermionic condensed matter systems
with matrix-product states

###### Abstract

Tensor network methods have progressed from variational techniques based on matrix-product states able to compute properties of one-dimensional condensed-matter lattice models into methods rooted in more elaborate states such as projected entangled pair states aimed at simulating the physics of two-dimensional models. In this work, we advocate the paradigm that for two-dimensional fermionic models, matrix-product states are still applicable to significantly higher accuracy levels than direct embeddings into one-dimensional systems allow for. To do so, we exploit schemes of fermionic mode transformations and overcome the prejudice that one-dimensional embeddings need to be local. This approach takes the insight seriously that the suitable exploitation of both the manifold of matrix-product states and the unitary manifold of mode transformations can more accurately capture the natural correlation structure. By demonstrating the residual low levels of entanglement in emerging modes, we show that matrix-product states can describe ground states strikingly well. The power of the approach is demonstrated by investigating a phase transition of spin-less fermions in lattices.

Recent years have enjoyed a flourishing development of tensor network methods, entanglement-based methods that allow to describe strongly correlated quantum many-body systems Orus (2014); Verstraete et al. (2008); Eisert (2013); Szalay et al. (2015); Eisert et al. (2010). They originate from the powerful density-matrix renormalization group (DMRG) White (1992); White and Noack (1992); Schollwöck (2011), a variational method building on matrix-product states Fannes et al. (1992); Perez-Garcia et al. (2007); Rommer and Östlund (1997) that captures the physics of one-dimensional local Hamiltonian systems provably well Hastings (2007); Arad et al. (); Verstraete et al. (2008); Schuch et al. (2008). It has been applied to countless physical systems (see the reviews Schollwöck (2005, 2011) and the comprehensive web page Nishino ()) and extended to time-evolving systems Vidal (2003); Daley et al. (2004); Haegeman et al. (2016), open systems Verstraete et al. (2004); Werner et al. (2016), and the study of excited states Khemani et al. (2016). Generalizing the variational set of matrix-product states to projected entangled pair states in two spatial dimensions, new avenues for the study of strongly correlated systems with tensor networks followed Verstraete and Cirac (); Verstraete et al. (2008); Orus (2014), including studies of fermionic models Barthel et al. (2009); Kraus et al. (2010); Pižorn and Verstraete (2010); Corboz (2016).

And yet, even if the DMRG approach has originally been devised to capture one-dimensional systems only: There are regimes in which it interestingly still performs competitively well Yan et al. (2011); Depenbrock et al. (2012) even in situations that at first seem alien to that type of approach and in which area laws for entanglement entropies are violated Eisert et al. (2010). Two-dimensional strongly correlated systems can be naturally embedded in highly non-local Hamiltonian models on a line. The high degree of entanglement that renders a variational approach based on matrix-product states challenging are partially compensated by the facts that contraction is efficient, and that very large bond dimensions are accessible. DMRG produces relevant data for strongly correlated matter even in two spatial dimensions, and for systems with fermionic degrees of freedom Stoudenmire and White (2012). The significance of this insight is even strengthened by the fact that DMRG is strictly variational, so that all ground state energies generated are upped bounds. And yet, given that the entanglement structure is not fully captured by matrix-product states, there are strong limitations of direct DMRG approaches.

In this work, we bring the idea of tackling two-dimensional strongly correlated matter with effectively one-dimensional matrix-product states to a new level. We show that the potential of an effective dimensional reduction is significantly more powerful than anticipated. We do so by systematically exploiting a degree of freedom that has not sufficiently been appreciated in the study of strongly correlated condensed-matter systems: This is the degree of freedom to adaptively define suitable modes in a strongly correlated fermionic system. Its significance is already manifest when solving problems in either real or in momentum space Xiang (1996); Nishimoto et al. (2012); Legeza and Sólyom (2003); Legeza et al. (2006); Murg et al. (2010); Ehlers et al. (2015); Motruk et al. (2016); Ehlers et al. (2017). For fermionic modes, however, there is an entire freedom that can be made use of and exploited when devising variational principles. In fact, a manifold structure emerges that originates from the tensor network and mode transformation degrees of freedom. Only the joint optimization fully exploits the potential of matrix-product state approaches in the study of strongly correlated fermionic condensed-matter system. It is this serious gap in the literature that is closed in this work: We overcome the prejudice that a one-dimensional embedding necessarily has to be an embedding in real space.

Setting. The Hilbert space of interacting fermions in modes is the fermionic Fock space originating from the basis constituted by all Slater determinants with . Let us furthermore denote with the fermionic annihilation operator of mode satisfying the canonical anti-commutation relations and and abbreviate . Matrix-product state vectors in this system then take the form

(1) |

We build upon ideas of adaptive fermionic mode transformations first presented in Refs. Krumnow et al. (2016, 2019), here brought to the level of applicability to condensed-matter lattice models in two spatial dimensions. To be specific, and to exemplify the power of our approach, the example of the spin-less interacting fermionic (spin-less Fermi-Hubbard) model

(2) |

will be in the focus of attention, where is the interaction strength, the hopping amplitude is set to 1 and denotes nearest neighbours on a 2d cubic lattice with . In addition, periodic boundary condition will be imposed along both spatial dimensions, which has been considered as a major bottleneck for MPS-based approaches. This example will show-cast that state-of-the-art energies can be reached. Having said that, in the mindset of this work would be any translationally invariant Hamiltonian of the form

(3) |

including local spin degrees of freedom. That is to say, the Hamiltonian is, regardless of its locality structure, treated as a long-ranged fermionic model on a one-dimensional line equipped with a given ordering.

Methods. We optimize the single particle basis in conjunction with the MPS tensors withing multiple successive mode transformation iterations. In our implementation, a single mode transformation iteration consists of a full forward and backward DMRG sweep without basis rotations followed by some number of additional sweeps with local mode transformations that adapt the single particle basis (compare Refs. Krumnow et al. (2016); Krumnow (2018)) that also rotate the couplings in the Hamiltonian to general couplings and . Typically, we have performed seven sweeps altogether. In the first DMRG sweep we have used the dynamically extended active space (DEAS) procedure Legeza and Sólyom (2003); Szalay et al. (2015) to achieve a good approximation of the ground state as a starting configuration for basis optimization. At the end of the seventh sweep, for the symmetric superblock configuration, we have calculated the site entropies, , the two-site mutual information, , the one-particle reduced density matrix, , and the occupation number distribution with . Here for is the von-Neumann entropy of the reduced state obtained from a partial trace of the full quantum state. Therefore, and are the one- and two-site reduced density matrices. The eigenvalues of define the natural occupation (NO) numbers, , and its eigenvectors the NO-basis. Based on we have calculated an optimized ordering using the Fiedler-vector approach Barcza et al. (2011), from a new complete active space vector for the DEAS procedure Legeza and Sólyom (2003) and from a new Hartree-Fock configuration. These together with the final rotated interaction matrices with components and are all used as inputs for the subsequent mode transformation iteration.

Therefore, each mode transformation iteration is followed by a global reordering in order to avoid being trapped in local minima Krumnow (2018). The basis optimization is carried out with fixed low bond dimension and or with a systematic increase of as will be discussed below. After convergence is reached the optimized basis is chosen from the mode transformation iteration leading to lowest energy and large scale DMRG calculations are performed with increasing bond dimension or using the dynamic block state selection (DBSS) approach with fixed truncation error threshold Legeza et al. (2003); Legeza and Sólyom (2004). In the remainder of this work, we denote these data as or , respectively. In addition, a given quantity obtained with from a calculation in the optimized basis will be indicated with a tilde on the top of it, e.g., the ground state energy at bond dimension as .

Numerical results. Our systematic error and convergence analysis will be given for the two dimensional lattice, since highly accurate reference data with the real space basis can also be generated. For larger system sizes, namely for and , only final results will be discussed and related figures together with additional numerical data are presented in the supplements.

First in the left panel of Fig. 1, we show the ground state energy , for as a function of mode transformation iterations using fixed bond dimensions and . Reference energies obtained in the real space basis are indicated with dashed lines for various bond dimensions up to . It very obvious indeed that exploiting mode transformations, gets significantly below even after the fourth iteration step and is below ).

In order to investigate how faithfully information beyond the ground state energy can be reproduced and predicted in the optimized basis, we depict in the right panel of Fig. 1 the operator norm of the difference of the one particle reduced density matrix over the mode transformation iterations and the real space reference data . Using the optimized basis, we also show the result for with increasing bond dimension , using different symbols. These latter data sets are basically the same for and , thus the optimal basis has already been obtained with the lower value. The error norms obtained with the real space basis are again much larger as indicated by the dashed lines. It is worth to remark that the error norm is less meaningful for very large bond dimensions since is below rendering potentially more accurate than .

In the inset of the left panel, we show the ground state energy as an inverse of the bond dimension for the real space basis and for the optimized basis with and . In the latter case, lie on the top of each other, indicating again that the optimal basis has been found with already (see red dots inside black circles).

For larger system sizes, the improvements are even more remarkable as is shown in Fig. 4 in the supplements for the lattice for different values of and for and . Here, is already lower than . In addition, reliable extrapolation with to the truncation free limit would require even significantly larger bond dimensions for the real space basis. In contrast to this, in case of the optimized basis, this is no longer an issue since is basically a flat curve. It has to be emphasized that our very accurate results have been obtained for a torus geometry. This reduces finite size effects significantly and much smaller systems sizes could lead to a reliable extrapolation to the thermodynamic limit (see the convergence of bond energy summarized in Tab. 1). Similar rapid scaling is found for the block entropy and other quantities as well.

The remarkable superiority of the optimized basis over the real space basis is due to the dramatic reduction of the entanglement via mode transformation, so that the effective one-dimensional embedding becomes feasible. As an indication of this, we depict the block entropy , in the left panel of Fig. 2 for various selected mode transformation iterations.

Here, the maximum of reduced by a full order of magnitude, as can be seen by comparing the blue (real space basis) and the black (optimized basis) curves. In addition, artifacts of the snake-like mapping of the two-dimensional lattice in real space into the one-dimensional MPS topology apparent in the blue curve are completely diminished by the basis optimization resulting in a smooth and highly symmetric profile (additional data is available in the supplements in Fig. 5). The iterative error norm of the block entropy measured between two subsequent mode transformation iterations, converges to -10 which can also be used as a criterion when to terminate the basis optimization. For larger values, the reduction is even more pronounced, leading to a state that is close to a Slater determinant. In the right panel, the maximum of for – which typically appears near the center of the chain – is laid out for various values for the real space basis and for the optimized one. While a strong dependence for is clearly visible in the real space basis, the curves basically fall on top of each other for the optimized basis. The small peak for signals the residual entanglement that cannot be removed by basis optimization which also controls the required bond dimension and thus the computational complexity. As a benchmark we have performed DMRG calculations using the DBSS approach with minimum bond dimension and a truncation error threshold . An agreement up to four digits has been obtained compared to the real space energy reference data calculated with , but we have gained a speedup by a full order of magnitude.

Through the course of basis optimization, the residual quantum correlations that have to be captured by the tensor network ansatz are significantly reduced. As a further proxy for this behaviour, one may investigate the sum of the single mode von-Neumann entropies that is reduced drastically, while pair-wise correlations reflected by get very much localized (for additional numerical data see Fig. 7). In addition, the investigation of the one-particle reduced density matrix shows that the optimized basis converges to the natural orbital basis as and tend to lie on the top of each other (Fig. 7). Therefore, here the final basis is the natural orbital basis, but the underlying basis has been systematically rotated by each mode transformation iterations. One might think that a natural step is to aim at identifying a globally optimal single particle basis could be more directly based on natural orbitals, i.e., by instead of using the local updates to the single particle basis one could rotate to the natural orbitals at the end of each mode transformation iteration. Such an approach has already been tested for quantum chemical applications Rissler et al. (2006), but a very unstable performance has been reported. In fact, we have also found that in the small- limit such an approach works acceptably, but for larger values it breaks down (see Fig. 6). The reason is that for small the optimal orbitals are Hartree-Fock like orbitals, while for large values localized orbitals seem to be more optimal. Our novel method based on fermionic mode transformation is, however, stable for all values. Importantly, it can also be used in general for interacting quantum many body systems.

Phase diagram. The power of our approach allows us to attack the physical properties such as the phase diagram of the system as well. In the limit of strong interactions, the model maps onto the anti-ferromagnetic Ising model in two dimensions and a charge density wave (CDW) phase develops. Since the hopping is restricted to nearest neighbours only, the Fermi surface takes the form of a square and perfect nesting together with Van Hove singularities providing strong arguments for an Ising transition into the CDW-ordered phase at Solyom (2011); Metzner et al. (2012). Furthermore, investigations within the Hartree-Fock approximation lead to an exponentially small order parameter in the weak coupling limit and the corrections starting from the limit, where Hartree-Fock theory becomes exact, provides only very small quantitative corrections in and even in Halvorsen et al. (1994). For this indicates a transition at and that the charge density wave order parameter is an exponential function of in the weak coupling limit. Note, however, that these simple arguments can break down as in the case of spin-less fermions in one spatial dimension, , where the model is exactly solvable because it reduces to the Heisenberg spin model and has a transition at finite Halvorsen et al. (1994). In addition, Ref. Uhrig and Vlaming (1993) has shown that there is a direct transition between the homogeneous and the CDW phases governed by phase separation, and a finite is suggested based on their obtained phase diagram. Their underlying arguments, however, have been derived for finite doping, thus an exponentially closing phase boundary between the CDW and phase separated phases together with cannot be ruled out.

In order to investigate the transition we first analyze the block entropy profiles for larger system sizes using the optimized basis and find that the peak for remains and its height increases with system size as is shown in the inset of Fig. 2. The center of the peak extracted from the spline fits ( for ) tends to shift to with while a finite value could have signaled a quantum phase transition Legeza and Sólyom (2006) at finite . In addition we compute the CDW order parameter Wang et al. (2014) as expectation value of directly, where in the real space basis and is a phase matrix with elements in a checker-board arrangement on the 2d-lattice. The real space simulations show that for large values of , takes a finite value while for it has to vanish as can be seen in Fig. 3. Unfortunately, the apparent finite size and dependencies do not allow us to conclusively decide upon the behaviour of for . Alternatively, the density-density correlation function can also be taken from the elements of the one- and two-particle reduced density matrices. The latter one has entries which can also be calculated efficiently by the DMRG method. Zgid and Nooijen (2008) Measuring these in the optimized basis and back-rotating to the real space basis, we found an agreement up to four digits between and the real space reference for and . For and the two data sets, however, began to deviate and possesses a much weaker dependence. Finite size scaling of the large scale DMRG data obtained with and is shown in the inset of Fig. 3 for various values. For large the curves scale to finite values in the thermodynamic limit, while for they show a slight downward curvature. For the extrapolated value of of the order of which might either be considered as being zero. Alternatively, after a rough extrapolation with and a spline fit on the extrapolated data (black crosses in the figure) an exponential opening of at can also be obtained. This functional form agrees to the one reported in Ref. Halvorsen et al. (1994) after some re-scaling and it is shown by a black curve in Fig. 3. Our approach hence pushes forward the capacity of the MPS based approaches to capture two dimensional strongly correlated systems significantly. Our results are in close agreement with analytic expectations (while details of the phase transition still remain an open problem numerically).

Conclusion. In this work, we have demonstrated that matrix-product state approaches, extending known DMRG methods, are surprisingly powerful for the simulation of two-dimensional quantum many body systems even imposing periodic boundary condition along both spatial dimensions. This is possible if only the key insight is acknowledged that one is not forced to do a local basis representation. Algorithmically, this is achieved by adaptively finding the optimal basis via fermionic mode transformation, optimizing over a larger manifold than that of MPS, which leads to a dramatic reduction of the correlations and entanglement in the system. A strongly interacting model in the real space basis thus can be converted to a weakly correlated problem in the optimized basis. Due to the torus geometry, finite size dependence is significantly reduced and intermediate system sizes make it possible to carry out more reliably extrapolations to the thermodynamic limit. In fact, for the two-dimensional translationally invariant spin-less fermion model, our results strongly suggest the presence of a quantum phase transition at , but the very small values of the charge density order parameter obtained numerically in the weak coupling limit leaves an uncertainty in our conclusion. The inclusion of a hopping between next nearest neighbours, however, would distort the square Fermi-surface and perfect nesting over an extended region of the momentum space will be destroyed. This is expected to have a have major effect, and divergencies in the susceptibilities might be removed and a finite is even more likely. This behaviour also shares features with the phase diagram of spin-less fermions on the honeycomb lattice Capponi (2017). In addition, physical properties of the transformed basis is also an important question. In general, the ground state energy cannot be written as a sum of energies of quasi-particle states except for special cases. The and large limits belong to the latter case (the ground state is a product state), but the residual block entropy for reflects the general scenario. Excitation energies, however, can be expressed through quasi-particle excitations (as forthcoming work will explore).

Our basis optimization is very robust, it can be carried out with low bond dimension, and calculations using the optimized basis can easily lead to an order of magnitude speedup in computational time. In addition, our method is stable for weakly and strongly interacting systems, in general, while standard approaches, like basis transformation based on natural orbitals, tried earlier have major limitations and drawbacks. In addition, the optimized basis for the spin-full Hubbard model does not resemble the characteristics of natural orbitals which reflects the existence of much stronger residual correlations in the system (as forthcoming work will explore). Conceptually most importantly, our work overcomes the deep misconception that lower-dimensional embeddings necessarily have to capture some kind of locality. Once this prejudice is overcome, acknowledging that fermionic mode transformations are not restricted to one-dimensional embeddings, mode reductions can be brought to a new level. Due to the polynomial scaling of the non-local DMRG White and Martin (1999), that is, , a reduction of by one or two orders of magnitude will place DMRG back into the competition for simulating higher dimensional and complex problems as well. Our approach has the potential to become a standard protocol for tensor network state methods on a daily basis.

Acknowledgement. J. E. has been supported by the Templeton Foundation, the ERC (TAQ), the DFG (EI 519/15-1, EI519/14-1, CRC 183 Project B01), and the European Unions Horizon 2020 research and innovation program under Grant Agreement No. 817482 (PASQuanS). Ö. L. has been supported by the Hungarian National Research, Development and Innovation Office (NKFIH) through Grant No. K120569, by the Hungarian Quantum Technology National Excellence Program (Project No. 2017-1.2.1-NKP-2017-00001). Ö. L. also acknowledges financial support from the Alexander von Humboldt foundation and useful discussions with Jenő Sólyom, Karlo Penc, and Florian Gebhard. L. V. has been supported by the Czech Science Foundation through Grant No. 18-18940Y.

## References

- Orus (2014) R. Orus, Ann. Phys. 349, 117 (2014).
- Verstraete et al. (2008) F. Verstraete, J. I. Cirac, and V. Murg, Adv. Phys. 57, 143 (2008).
- Eisert (2013) J. Eisert, Mod. Sim. 3, 520 (2013).
- Szalay et al. (2015) S. Szalay, M. Pfeffer, V. Murg, G. Barcza, F. Verstraete, R. Schneider, and Ö. Legeza, Int. J. Quant. Chem. 115, 1342 (2015).
- Eisert et al. (2010) J. Eisert, M. Cramer, and M. B. Plenio, Rev. Mod. Phys. 82, 277 (2010).
- White (1992) S. R. White, Phys. Rev. Lett. 69, 2863 (1992).
- White and Noack (1992) S. R. White and R. M. Noack, Phys. Rev. Lett. 68, 3487 (1992).
- Schollwöck (2011) U. Schollwöck, Ann. Phys. 326, 96 (2011).
- Fannes et al. (1992) M. Fannes, B. Nachtergaele, and R. F. Werner, Commun. Math. Phys. 144, 443 (1992).
- Perez-Garcia et al. (2007) D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac, Quantum Info. Comput. 7, 401 (2007).
- Rommer and Östlund (1997) S. Rommer and S. Östlund, Phys. Rev. B 55, 2164 (1997).
- Hastings (2007) M. B. Hastings, J. Stat. Mech. 2007, P08024 (2007).
- (13) I. Arad, A. Kitaev, Z. Landau, and U. Vazirani, arxiv:1301.1162 .
- Schuch et al. (2008) N. Schuch, M. M. Wolf, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 100, 030504 (2008).
- Schollwöck (2005) U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
- (16) T. Nishino, “Unofficial DMRG, MPS, TPS, and MERA home page,” Http://quattro.phys.sci.kobe-u.ac.jp/dmrg/condmat.html.
- Vidal (2003) G. Vidal, Phys. Rev. Lett. 91, 147902 (2003).
- Daley et al. (2004) A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, J. Stat. Mech. 2004, P04005 (2004).
- Haegeman et al. (2016) J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Phys. Rev. B 94, 165116 (2016).
- Verstraete et al. (2004) F. Verstraete, J. J. Garcia-Ripoll, and J. I. Cirac, Phys. Rev. Lett. 93, 207204 (2004).
- Werner et al. (2016) A. H. Werner, D. Jaschke, P. Silvi, M. Kliesch, T. Calarco, J. Eisert, and S. Montangero, Phys. Rev. Lett. 116, 237201 (2016).
- Khemani et al. (2016) V. Khemani, F. Pollmann, and S. L. Sondhi, Phys. Rev. Lett. 116, 247204 (2016).
- (23) F. Verstraete and J. I. Cirac, cond-mat:0407066 .
- Barthel et al. (2009) T. Barthel, C. Pineda, and J. Eisert, Phys. Rev. A 80, 042333 (2009).
- Kraus et al. (2010) C. V. Kraus, N. Schuch, F. Verstraete, and J. I. Cirac, Phys. Rev. A 81, 052338 (2010).
- Pižorn and Verstraete (2010) I. Pižorn and F. Verstraete, Phys. Rev. B 81, 245110 (2010).
- Corboz (2016) P. Corboz, Phys. Rev. B 93, 045116 (2016).
- Yan et al. (2011) S. Yan, D. A. Huse, and S. R. White, Science 332, 1173 (2011).
- Depenbrock et al. (2012) S. Depenbrock, I. P. McCulloch, and U. Schollwöck, Phys. Rev. Lett. 109, 067201 (2012).
- Stoudenmire and White (2012) E. M. Stoudenmire and S. R. White, An. Rev. Cond. Mat. Phys. 3, 111 (2012).
- Xiang (1996) T. Xiang, Phys. Rev. B 53, R10445 (1996).
- Nishimoto et al. (2012) S. Nishimoto, E. Jeckelmann, F. Gebhard, and R. M. Noack, Phys. Rev. B. 65, 165114 (2012).
- Legeza and Sólyom (2003) Ö. Legeza and J. Sólyom, Phys. Rev. B 68, 195116 (2003).
- Legeza et al. (2006) Ö. Legeza, F. Gebhard, and J. Rissler, Phys. Rev. B 74, 195112 (2006).
- Murg et al. (2010) V. Murg, F. Verstraete, Ö. Legeza, and R. M. Noack, Phys. Rev. B 82, 205105 (2010).
- Ehlers et al. (2015) G. Ehlers, J. Sólyom, Ö. Legeza, and R. M. Noack, Phys. Rev. B 92, 235116 (2015).
- Motruk et al. (2016) J. Motruk, P. M. Zaletel, S. K. R. Mong, and F. Pollmann, Phys. Rev. B 93, 155139 (2016).
- Ehlers et al. (2017) G. Ehlers, S. R. White, and R. M. Noack, Phys. Rev. B. 95, 125125 (2017).
- Krumnow et al. (2016) C. Krumnow, L. Veis, Ö. Legeza, and J. Eisert, Phys. Rev. Lett. 117, 210402 (2016).
- Krumnow et al. (2019) C. Krumnow, J. Eisert, and Ö. Legeza, arXiv:1904.11999 (2019).
- Krumnow (2018) C. Krumnow, Detecting and understanding efficient structures in finite fermionic systems, Ph.D. thesis, Freie Universität Berlin (2018).
- Barcza et al. (2011) G. Barcza, Ö. Legeza, K. H. Marti, and M. Reiher, Phys. Rev. A 83, 012508 (2011).
- Legeza et al. (2003) Ö. Legeza, J. Röder, and B. A. Hess, Phys. Rev. B 67, 125114 (2003).
- Legeza and Sólyom (2004) Ö. Legeza and J. Sólyom, Phys. Rev. B 70, 205118 (2004).
- Rissler et al. (2006) J. Rissler, R. M. Noack, and S. R. White, Chem. Phys. 323, 519 (2006).
- Solyom (2011) J. Solyom, Fundamentals of the Physics of Solids, Vol. 3, Normal, Broken-Symmetry, and Correlated Systems (Springer, 2011).
- Metzner et al. (2012) W. Metzner, M. Salmhofer, C. Honerkamp, V. Meden, and K. Schönhammer, Rev. Mod. Phys. 84, 299 (2012).
- Halvorsen et al. (1994) E. Halvorsen, G. S. Uhrig, and G. Czycholl, Z. Phys. B. 94, 291 (1994).
- Uhrig and Vlaming (1993) G. S. Uhrig and R. Vlaming, Phys. Rev. Lett. 71, 271 (1993).
- Legeza and Sólyom (2006) Ö. Legeza and J. Sólyom, Phys. Rev. Lett. 96, 116401 (2006).
- Wang et al. (2014) L. Wang, P. Corboz, and M. Troyer, New J. Phys. 16, 103008 (2014).
- Zgid and Nooijen (2008) D. Zgid and M. Nooijen, J. Chem. Phys. 128, 144115 (2008).
- Capponi (2017) S. Capponi, J. Phys.: Condens. Matter 29, 043002 (2017).
- White and Martin (1999) S. R. White and R. L. Martin, J. Chem. Phys. 110, 4127 (1999).

## Appendix A Supplemental material: Additional data for larger systems

In this section, we present additional numerical data together with further scaling properties obtained for larger system sizes.

0.25 | 0.5 | 1 | 2 | 4 | 8 | |
---|---|---|---|---|---|---|

-0.6636 | -0.5845 | -0.4518 | -0.2878 | -0.1591 | -0.0823 | |

-0.6847 | -0.5992 | -0.4583 | -0.2911 | -0.1601 | -0.0825 | |

-0.6947 | -0.6059 | -0.4604 | -0.2912 | -0.1601 | -0.0825 | |

-0.6983 | -0.6079 | -0.4605 | -0.2912 | -0.1601 | -0.0825 |