# Hubbard model on the Bethe lattice via variational uniform tree states: Metal-insulator transition and a Fermi liquid

Peter Lunts, Antoine Georges, E. Miles Stoudenmire, Matthew Fishman

► To cite this version:

Peter Lunts, Antoine Georges, E. Miles Stoudenmire, Matthew Fishman. Hubbard model on the Bethe lattice via variational uniform tree states: Metal-insulator transition and a Fermi liquid. Physical Review Research, 2021, 3 (2), pp.023054. 10.1103/PhysRevResearch.3.023054. hal-05000888

HAL Id: hal-05000888

https://college-de-france.hal.science/hal-05000888v1

Submitted on 21 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Hubbard model on the Bethe lattice via variational uniform tree states:
Metal-insulator transition and a Fermi liquid

Peter Lunts,¹,*, Antoine Georges,¹,²,³,⁴ E. Miles Stoudenmire,¹ and Matthew Fishman¹

¹Center for Computational Quantum Physics, Flatiron Institute, 162 5th Avenue, New York, New York 10010, USA

²Collège de France, 11 place Marcelin Berthelot, 75005 Paris, France

³CPHT, CNRS, École Polytechnique, IP Paris, F-91128 Palaiseau, France

⁴DQMP, Université de Genève, 24 quai Ernest Ansermet, CH-1211 Genève, Switzerland

<figure><img src="image_2.png" /></figure>

(Received 2 November 2020; revised 1 March 2021; accepted 29 March 2021; published 19 April 2021)

We numerically solve the Hubbard model on the Bethe lattice with finite coordination number $z = 3$, and determine its zero-temperature phase diagram. For this purpose, we introduce and develop the “variational uniform tree state” (VUTS) algorithm, a tensor-network algorithm which generalizes the variational uniform matrix product state algorithm to tree tensor networks. Our results reveal an antiferromagnetic insulating phase and a paramagnetic metallic phase, separated by a first-order doping-driven metal-insulator transition. We show that the metallic state is a Fermi liquid with coherent quasiparticle excitations for all values of the interaction strength $U$, and we obtain the finite quasiparticle weight $Z$ from the single-particle occupation function of a generalized “momentum” variable. We find that $Z$ decreases with increasing $U$, ultimately saturating to a nonzero, doping-dependent value. Our work demonstrates that tensor-network calculations on tree lattices, and the VUTS algorithm in particular, are a platform for obtaining controlled results for phenomena absent in one dimension, such as Fermi liquids, while avoiding computational difficulties associated with tensor networks in two dimensions. We envision that future studies could observe non-Fermi liquids, interaction-driven metal-insulator transitions, and doped spin liquids using this platform.

DOI: 10.1103/PhysRevResearch.3.023054

# I. INTRODUCTION

The Hubbard model is a cornerstone of condensed mat- ter physics. As a paradigmatic model of strongly correlated electrons [1], it is simple to formulate yet rich in behavior. In two dimensions (relevant, e.g., to cuprate superconductors) observed behaviors include, but are not limited to, antifer- romagnetism, unconventional metallic behavior characterized by a pseudogap, and deviations from Fermi-liquid theory [2–6], as well as stripe orders [7,8] closely competing with superconducting states at low temperatures [9]. But, despite decades of effort, a comprehensive understanding of the phase diagram of the two-dimensional Hubbard model has not yet been fully reached. Therefore, any solutions of the Hubbard model, whether obtained analytically or by accurate and con- trolled numerical techniques, are of great value.

The most reliable and comprehensive solutions of the Hub- bard model obtained so far have been mainly in (quasi-)one dimension [10] and in infinite dimensions (infinite lattice co- ordination number) [11–13]. In one dimension, solutions can be obtained by integrability [14] and bosonization methods

[15], as well as numerically with matrix product state (MPS) tensor-network methods [16,17]. The latter can also reliably treat quasi-one-dimensional ladder or cylindrical geometries with a small transverse size. In the limit of infinite dimen- sions the Hubbard model is again numerically tractable due to the fact that dynamical mean-field theory (DMFT) be- comes exact [13]. Using this technique, the phase diagram of the infinite-dimensional Hubbard model can be mapped out and a detailed understanding of the interaction-driven Mott- insulator-to-metal transition has been established (for reviews, see Refs. [13,18–21]).

While they provide useful insights into the physics of the two-dimensional Hubbard model, these limiting cases also have peculiarities which limit the generality of the conclusions that can be drawn from their study. In the limit of infinite dimensions, the metallic state is found to be a Fermi liquid, with interactions affecting one-particle properties in a local (momentum-independent) manner only. Hence, in this limit, the feedback of long-wavelength collective modes or even short-range spatial correlations on quasiparticle properties is entirely absent. These effects are important, in particular close to a critical point. In two dimensions, they are, for example, responsible for the formation of a pseudogap [4,22–25]. In contrast, in one dimension the low-energy excitations consist only of bosonic collective modes, associated with charge and spin degrees of freedom. There is no notion of a Fermi liquid and the metallic behavior of the Hubbard model is a Luttinger liquid which lacks coherent quasiparticles and displays spin- charge separation [15].

<footnote>*plunts@flatironinstitute.org</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. A portion of the infinite Bethe lattice with coordination number $z = 3$. The figure depicts four “generations” of the tree structure, with a particular site chosen to be the “center” site. Note that in the actual Bethe lattice there is no special site, and generations can be counted from any of the sites.</figcaption></figure>

In this work we perform a controlled and accurate numerical study of the ground state of the Hubbard model on the Bethe lattice with a finite coordination number $z$, focusing on the case $z = 3$. This is an infinite lattice that has a tree structure, where every site is connected to the same number of other sites ($z$) but there are no loops. We show a portion of this lattice in Fig. 1. This lattice provides an intermediate case between one dimension ($z = 2$) and infinite dimensions ($z = \infty$), with the key virtue that it admits controlled solutions via tensor-network methods, including away from half-filling and in the presence of strong interactions.

Exact solutions of models on the Bethe lattice have a long history in statistical mechanics [26], starting with the pioneering article of Bethe [27]. Solutions on the finite coordination number Bethe lattice provide a better approximation to thermodynamic quantities than the mean-field approximation (corresponding to the infinite-dimensional limit) [26–28]. Models studied on the Bethe lattice include classical and quantum spin models [29–35], spin-glass systems [36–39], the Bose Hubbard model [40], and models of Anderson localization [41–45]. The fermionic Hubbard model on the finite version of the $z = 3$ Bethe lattice (known as a Cayley tree) has also been studied previously using a variant of the density matrix renormalization group (DMRG) algorithm [46], but only the case of half-filling was studied (which is a charge insulator) and only local ground-state quantities given (energy, staggered magnetization and its fluctuations, and neighboring spin correlations). Since that time, there have been significant advances in DMRG and related algorithms for infinite one-dimensional systems, which we generalize here to the Bethe lattice and use to obtain our results. Notably, there has been no previous study of metallic states on the finite connectivity Bethe lattice, to the best of our knowledge.

We determine the full phase diagram of the fermionic Hubbard model on the $z = 3$ Bethe lattice, allowing for a two-site unit cell, and find an antiferromagnetic insulator and a paramagnetic metal as the only two phases. We establish the nature of the doping-driven Mott-insulator-to-metal transition (MIT). We find that it is first order, and that for every value of interaction strength there is a region of forbidden density.

Therefore, in the interaction-density plane the model exhibits phase separation at low doping levels. We find that, for all allowed values of the density, the doped metallic ground state does not display magnetic long-range order.

The restriction of the unit-cell size to two sites has the potential to exclude lower-energy states that have a smaller “ordering wave vector” (the concept of momentum does not apply on the Bethe lattice). In the two-dimensional Hubbard model in particular, initial findings of phase separation [47] have been superseded by lower-energy states with complex orders involving large two-dimensional “supercells” [8,48]. However, a very recent work suggests that actually phase separation does indeed occur, with multiple first-order transition being present in between the insulator at half-filling and the heavily doped Fermi liquid [49]. On the other hand, in one dimension, it is known from the exact solution via Bethe Ansatz that there are only two phases [14], the same ones as found in this paper, and the transition between them is continuous [50,51]. The Bethe lattice and the one-dimensional lattice share the property of being loopless, and we believe that this feature is the crucial one when it comes to the phase diagram since the orders found in two dimensions always have a two-dimensional nature (i.e., the “supercell” has a nontrivial thickness). It is therefore reasonable to expect that no other orders are present, and a two-site unit cell is sufficient to capture the physics of this model, but of course further study allowing for larger unit cells is required to definitively answer this question. On the other hand, when it comes to the order of the transition, the Bethe lattice differs from one dimension, and might be more similar to two dimensions.

The main motivation of our work is to study the metallic state. Generally, it is difficult for tensor networks to accurately describe interacting metallic states above one dimension, although much progress has been made in this direction [48,52–54]. Our work provides an alternate route to this agenda, avoiding the computational challenges of two-dimensional tensor networks, while going beyond the restrictions of one-dimensional physics. We demonstrate that the doped metal hosts coherent quasiparticles at all studied values of the interaction strength $U$, from weak to strong coupling, and determine the behavior of the quasiparticle weight as a function of $U$. This answers in the affirmative the question of whether Fermi-liquid behavior applies as soon as the peculiar kinematic constraints of one dimension are alleviated, and also provides a concrete description of a Fermi-liquid ground state with tensor networks.

We obtain our results by generalizing a recently developed MPS method, variational uniform MPS (VUMPS) [55], to tree tensor-network (TTN) states [56–58], which we dub the variational uniform tree state algorithm (VUTS). We further introduce the fermionic version of the VUTS algorithm, using the swap-gate method of Refs. [52,59]. The VUTS algorithm works directly in the thermodynamic limit, which is important in the study of models on the Bethe lattice. The alternative is to study the finite Cayley tree and perform a finite-size scaling analysis. However, the number of boundary sites on the Cayley tree is always more than half of the total, and therefore finite-size effects are unusually strong and can even lead to conclusions that do not hold on the Bethe lattice [26,60,61]. Working directly in the infinite-size limit is therefore
important for models on trees. All previous works studying quantum models on the (infinite) Bethe lattice using tensor networks have used a variant of the infinite time-evolving block decimation (ITEBD) algorithm [62]. In the one-dimensional case, the VUMPS algorithm [55] has been found to be much more efficient than other methods that work in the thermodynamic limit, such as iTEBD or earlier infinite DMRG algorithms [63–66], and indeed we find its extension in the form of the VUTS algorithm we develop to be very efficient. Our method scales as $O(\chi^{z+1})$ where $\chi$ is the bond dimension of the tensor network being optimized. For $z=3$ this scaling is significantly better than that of the most modern and accurate projected entangled pair states (PEPS) algorithms that scale as $O(\chi^{10})$ for PEPS bond dimension $\chi$ (assuming the boundary MPS bond dimension scales as $\chi^2$) [48,67–71], but is more challenging than the $O(\chi^3)$ scaling of DMRG. However, the steeper scaling is mitigated by the fact that the typical bond dimension required to reach an accurate solution generally decreases as one goes to higher dimensions and larger coordination numbers due to the monogamy of entanglement and more mean-field-like properties of the wave function.

The accuracy of tensor-network methods is often measured by the truncation error, denoted by $\epsilon_\rho$, which measures the typical loss of fidelity incurred during the truncation step of the optimization algorithm. In this work, for a bond dimension of $\chi = 100$, we are able to achieve $\epsilon_\rho < 10^{-3}$ in the most computationally challenging part of the phase diagram (most entangled ground state), and $\epsilon_\rho < 10^{-7}$ in the best cases. At face value, $\epsilon_\rho \sim 10^{-3}$ lies in-between the accuracy standards of DMRG and PEPS. For DMRG the high-end standard is $\epsilon_\rho \sim 10^{-8}$. For (i)PEPS one needs to define a different but similar truncation error denoted by $w$ [53], which currently has an accuracy standard of $w \sim 10^{-1}$ in the doped regime of the Hubbard model (computationally more challenging than the half-filled case) [8,48]. However, the more important metric for the accuracy of results is their scaling and extrapolation with either bond dimension or truncation error. In our work, we observe that all the quantities we measure have a relatively smooth dependence on the bond dimension that is well approximated to be linear for the larger bond dimensions we use. Therefore, our extrapolations to the infinite bond dimension limit are expected to be reliable. In particular, the smooth behavior in $\chi$ of the correlation functions and their properties at all distances give a similarly smooth behavior of the occupation function near the Fermi energy, which allows us to reliably extract information about quasiparticle coherence in the metallic phase. This demonstrates that our method can be used to reliably study critical phases of matter on the Bethe lattice already at the level of accuracy we achieve here.

This work suggests promising future directions for studying the behavior of strongly correlated electrons in a controlled setup. Because Fermi-liquid behavior is a rather generic feature of metallic states, this study allows to establish a controlled platform which can be used to study how Fermi-liquid behavior can be broken by further perturbations to the model considered here, or in other fermionic models. In the concluding section of this paper, we discuss possible routes towards achieving this goal. If successful, tensor-network solutions of correlated electrons on the $z=3$ Bethe lattice

<figure><img src="image_4.png" /><figcaption>FIG. 2. A finite portion of the infinite tree tensor-network (TTN) state describing the many-body wave function on the Bethe lattice. The physical legs (green dashed lines) form the nodes of the lattice, while the virtual legs (straight black lines) form the edges. In this case, a single tensor A comprises the state, and the unit cell is just a single site.</figcaption></figure>

could provide a new platform for studying non-Fermi liquids [72,73] in a controlled and accurate manner. Other potential applications are the study of the interaction-driven MIT in frustrated fermionic systems, and the study of fermions on closely related treelike lattices, such as the Husimi cactus on which the Heisenberg model and other spin models have been shown to display spin-liquid phases [74]. We elaborate on all these directions and others at the end of the paper.

This paper is organized as follows. In Sec. II we describe the general VUTS method, applicable to generic Hamiltonians. In Sec. III we define the Hubbard model on the Bethe lattice and show the phase diagram obtained from the VUTS solution. Section IV discusses the calculation of the quasiparticle weight from the occupation function, as well as Luttinger's theorem. Finally, in Sec. V we summarize and discuss future directions.

## II. VARIATIONAL UNIFORM TREE STATE ALGORITHM

In this section we introduce the variational uniform tree state algorithm (VUTS), a generalization of the variational uniform matrix product state algorithm (VUMPS) [55], for optimizing infinite tree tensor-network (TTN) states. We start with a Bethe lattice of quantum degrees of freedom. For simplicity we focus on the algorithm for coordination number $z=3$, which is the value for the model studied in this paper, but the extension to general $z$ is straightforward. We use an infinite TTN as our *Ansatz* to approximate quantum states on the Bethe lattice. For $z=3$, the infinite TTN with a one-site unit cell consists of an order 4 tensor $A \in C^{\chi \times \chi \times \chi \times d}$, with one physical leg ($s$) which runs over the physical degrees of freedom $1, \dots, d$, and three virtual legs ($l_0, l_1, l_2$) that run over virtual degrees of freedom $1, \dots, \chi$. The virtual legs of neighboring tensors connect to each other, forming the same geometry as the Bethe lattice, as shown in the tensor-network diagram in Fig. 2.

The Hamiltonians we will focus on here are isotropic and have an equivalence between all sites (the analog of translational invariance for hypercubic lattices). The ground state will potentially break this isotropy completely, and break the site equivalence down to a nontrivial unit cell. In this paper, we allow for the state to be fully anisotropic between
<figure><img src="image_5.png" /><figcaption>FIG. 3. Diagrammatic version of the gauge conditions of Eq. (1). The bonds labeled $k, n, m$ have link indices $l_k, l_n, l_m$, respectively (and the uncontracted $m$ bond on the ket has link index $l'_m$). The unlabeled dashed bond is the physical degree of freedom with index $s_i$.</figcaption></figure>

the different directions emanating from a given site, but for simplicity we focus exclusively on the case when the unit cell consists of no more than two sites (generalizing to arbitrary unit cells is straightforward). The infinite TTN state we study therefore has a two-site unit cell, and is parametrized by a set of 2z tensors $A_{i,m}$, where $i = 0, 1$ labels the location in the unit cell, and $m = 1, \dots, z$ labels the direction of the gauge (defined below). Again, each tensor has one physical index $s_i = 1, \dots, d$ and three virtual “link” indices $l_0, l_1, l_2$.

Because the TTN has no loops, it is straightforward to work in the canonical gauge, i.e. the gauge where the tensors are constrained to be orthonormal bases when viewed as a matrix from two link indices $l_n, l_k$ and the physical index $s_i$ to the remaining link index $l_m$. This constraint on the tensors is very useful for making the variational optimization faster and more stable, and is standard in a wide variety of tensor-network algorithms, particularly in one-dimensional (1D) algorithms like VUMPS and DMRG. The constraint on the tensors is written as

$$ \sum_{s_i, l_n, l_k} \bar{A}_{i,m}^{s_i, l'_m, l_n, l_k} A_{i,m}^{s_i, l_m, l_n, l_k} = \mathbb{1}_{i,m}^{l'_m, l_m}, \quad (1) $$

where we use $\bar{A}$ to denote the complex conjugation of A. The matrices $\mathbb{1}_{i,m}$ are identities. Diagrammatically, Eq. (1) is equivalent to Fig. 3. The arrows on the links denote the gauge of the A tensors (the outgoing link is the direction of the gauge). Any TTN can be brought into the form where the tensors obey Eq. (1) (or equivalently Fig. 3) by inserting a particular set of “gauge transformations” onto the link degrees of freedom, i.e., inserting a particular set of resolutions of the identity $XX^{-1}$ (where X is an invertible matrix) onto the links of the TTN. Note that the gauge transformation does not affect the observables of the system, and any TTN can be transformed into the canonical gauge efficiently.

Examples of the TTN state with a two-site unit cell in the canonical gauge are shown in Fig. 4. In the top diagram, the gauge center is $C_2$. Here, $C_2$ represents the projection of the infinite wave function of the system onto the finite-sized Hilbert space of that link of the network. The center matrices $C_m$ constitute invertible gauge transformations relating the A tensors to each other via $A_{i,m}C_m = A_{i,n}C_n$, where $i = 0, 1$

<figure><img src="image_6.png" /><figcaption>FIG. 4. The same portion of the Bethe lattice as in Fig. 2, but with the bonds labeled with $m = 0, 1, 2$, which have link indices $l_0, l_1, l_2$, respectively. Additionally, each tensor is labeled with subscripts $i, m$, where $i = 0, 1$ is the unit-cell index and $m$ is the direction of the gauge. In the top diagram, the gauge center $C_2$ is shown on a bond labeled by 2. The next equality shows that the $C_2$ tensor can be absorbed into the $A_{1,2}$ tensor to put the gauge center on the site tensor, creating $A_{1,C}$.</figcaption></figure>

and $n \neq m$. The center matrices $C_m$ also contain important information like the entanglement spectrum between the two infinite halves of the system split by that link. Additionally, the gauge center can be absorbed onto an A tensor, defining the center site tensors $A_{i,C} = A_{i,m}C_m$ for any $m$. This is shown in the lower diagram of Fig. 4. In this case, $A_{i,C}$ would represent the infinite wave function of the system projected onto a single site (and again, different spectra of that tensor relate to entanglement spectra of different bipartitions of the lattice).

We describe the algorithm for the case of $H = \sum_{\langle i,j \rangle} h_{i,j}$, where $h_{i,j}$ is a two-site operator that acts on nearest neighbors only. The case of longer-range operators is treatable using techniques like those described in Appendix C of Ref. [55]. The total energy is given by $E = \sum_{\langle i,j \rangle} \langle \psi | h_{i,j} | \psi \rangle$, and we want to minimize $E$, treating the tensor elements of our TTN as variational parameters. As in VUMPS (and many related tensor-network ground-state methods), VUTS proceeds in three main steps that are iterated until convergence:

(1) Compute the projected Hamiltonians (the Hamiltonian projected into the basis corresponding to the virtual degrees of freedom of the network) to turn the global optimization into a local optimization problem.

(2) Find the optimized tensors by minimizing the energy of the projected Hamiltonian.

(3) Update the tensor network with the new optimized tensors.
<figure><img src="image_7.png" /><figcaption>FIG. 5. (a) The first few terms of the series for $H_{1,2}$, the projected Hamiltonian contribution for one of the three branches of the infinite Bethe lattice connected to the center site tensor $A_{0,C}$ (on the $i=0$ sublattice). (b) Definition of the $h_{i,m}$ tensors used in (a), which are a sum of two local Hamiltonian environment tensors sitting on two branches of the Bethe lattice. Note that the tensor labeled $h$ represents the two-site operator term $h_{i,j}$, which the Hamiltonian is made up of.</figcaption></figure>

To begin, say we are interested in optimizing a one-site projected wave function $A_{i,C}$, as defined in Fig. 4. Step 1 requires computing infinite sums of local Hamiltonian terms, projected into the basis of our gauged TTN (defined by the tensors $A_{i,m}$), for each of the $z=3$ infinite subtrees connected to $A_{i,C}$. In order to perform the infinite sum, we focus on summing the energy contributions of a single subtree. An example for the series that needs to be summed for the $m=2$ direction in order to optimize the $A_{0,C}$ tensor is shown in Fig. 5. In Fig. 5 and onwards we use the notation $\bar{0}=1$, $\bar{1}=0$ for the unit-cell indices. We define the results of these summations as the matrices $H_{i,m}$. The summation can be carried out by making use of the fact that the sum is a geometric series. However, care has to be taken to project out infinite energy contributions to keep the series convergent (i.e., keep the norm of the solution $H_{i,m}$ from diverging). The procedure of performing the summation and projecting out the infinite energy contributions is a generalization of the one in Appendix D of Ref. [55], and we discuss it in more detail in Appendix A 1.

Once the environment tensors are found, we can proceed to step 2 of the algorithm, which we begin by optimizing $A_{i,C}$. This is done by finding the ground state of the Hamiltonian projected onto the sublattice site $i$, a standard procedure in VUMPS and DMRG. The eigenvalue equation for $A_{i,C}$ is shown diagrammatically in Fig. 6, and is solved iteratively (using a Hermitian eigensolver such as Lanczos). To find the TTN ground state, we obtain the eigenvector with the smallest eigenvalue. As in the VUMPS algorithm, in addition to optimizing $A_{i,C}$, we also optimize $C_m$. $A_{i,C}$ and $C_m$ are then used to solve for $A_{i,m}$, which make up the updated infinite TTN state (see next paragraph). The eigenvalue equation for $C_m$ is shown diagrammatically in Fig. 7.

Finally, once $A_{i,C}$, $C_m$ for all $i=0, 1, m=0, 1, 2$ are optimized, we can proceed to step 3 of the algorithm and solve for our new $A_{i,m}$ tensors by minimizing

$$ \epsilon_{i,m} = \min_{A_{i,m}^{\dagger}A_{i,m}=1} \|A_{i,C} - A_{i,m}C_m\|. \quad (2) $$

This minimization problem can be solved optimally using techniques described in Eqs. (18)–(22) of Ref. [55]. The new $A_{i,m}$ we obtain constitute our updated TTN, and steps 1–3 are repeated until convergence. Convergence is achieved when the largest error found in Eq. (2), $\epsilon_{prec} \equiv \max\{\epsilon_{i,m}\}$, falls below a chosen threshold (e.g., $\epsilon_{prec} < 10^{-12}$).

The VUTS algorithm with a one-site update, as we describe here, scales as $O(\chi^{z+1})$ which becomes $O(\chi^4)$ for the $z=3$ Bethe lattice and $O(\chi^3)$ for $z=2$, thus reducing to the scaling of VUMPS in the $z=2$ case. Additionally, a two-site update can be formulated, analogous to the two-site DMRG algorithm which is commonly used. This requires a slight modification of the algorithm where ground states of two-site and one-site projected Hamiltonians are computed (as opposed to one-site and zero-site projected Hamiltonians in the version of the algorithm described above). This can lead to improved convergence since a larger local Hilbert space is explored, but has a higher computational cost of $O(\chi^5)$ for $z=3$. We use this technique at lower bond dimensions in more challenging parts of the phase diagram (near the phase transition), and switch to the one-site algorithm later in the calculation to reach higher bond dimensions. To dynamically change the bond dimensions, we use a generalization of the subspace expansion procedure described in Appendix B of Ref. [55].

For fermionic models like the Hubbard model, we need to use a fermionic version of VUTS. We use the method outlined in Refs. [52,59]. Every tensor is now endowed with a fermion parity $Z_2$ quantum number and is parity preserving. When two tensor legs cross on a planar projection of a tensor diagram, a fermionic swap gate is placed at the crossing. In order to employ this method, we need to use a fixed ordering convention for the legs of the tensors $A_{i,m}$, $A_{i,C}$, $C_m$, which must be kept consistent in all of the diagrams in the calculation. We address the associated subtleties and details of this approach in Appendix A 2. Other symmetries beyond the $Z_2$ parity can also be used, such as U(1) particle number conservation (to fix the filling), U(1) spin projection symmetry in the $z$ direction, and also the spin SU(2) symmetry. The inclusion of these symmetries makes the tensors more sparse and should therefore make the tensor operations more efficient, allowing us to reach
<figure><img src="image_8.png" /><figcaption>FIG. 6. The eigenvalue equation for $A_{i,C}$. The sum over $m$ is a sum over the contributions from each leg of the $A_{i,C}$ tensor.</figcaption></figure>

larger bond dimensions. In this work, we only employ parity quantum numbers, and leave the use of additional symmetries for future work.

### III. MODEL AND PHASE DIAGRAM

The Hubbard Hamiltonian is given by

$$H = -t \sum_{\langle i,j \rangle} \sum_{\sigma=\uparrow,\downarrow} c_{i,\sigma}^{\dagger} c_{j,\sigma} + U \sum_{i} n_{i,\uparrow} n_{i,\downarrow} - \mu \sum_{i,\sigma} n_{i,\sigma}, \quad (3)$$

<figure><img src="image_9.png" /><figcaption>FIG. 7. The eigenvalue equation for $C_m$.</figcaption></figure>

where the site index $i$ now runs over the Bethe lattice and $n_{i,\sigma} = c_{i,\sigma}^{\dagger} c_{i,\sigma}$ is the onsite density for an electron of spin $\sigma$. We set $t=1$ and vary $U \ge 0$ and $\mu$. Since the model is particle-hole symmetric, we only need to consider $\delta\mu \equiv \mu - \frac{U}{2} \ge 0$. To obtain the phase diagram, we compute the ground state of the Hubbard model using the fermionic VUTS algorithm for several values of the bond dimension $\chi$. The details of the numerical calculation are given in Appendix B. We perform extrapolations in $\chi$, the results of which we describe below. The additional plots of the finite $\chi$ results and details of the extrapolations are in Appendix C.

At half-filling, $\delta\mu = 0$, the system is a charge insulator with antiferromagnetic order for all $U > 0$, as in one dimension. To illustrate this we compute the staggered magnetization of the bipartite sublattice $m_s = |\langle\bar{S}_A\rangle - \langle\bar{S}_B\rangle|/2$. The extrapolated values $m_s(\chi \to \infty)$ are shown in Fig. 8. We can see that for any $U > 0$ the magnetization is nonzero, tending to zero as $U \to 0$. From general mean-field theory considerations we expect $m_s \sim e^{-\bar{\nu}}$ for small $U$. However, we did not attempt to confirm this functional form numerically by systematic calculations at very small values of $U$.

To illustrate the charge gap at half-filling, we compute the onsite occupation $\langle n_i \rangle = \sum_{\sigma} \langle n_{i,\sigma} \rangle$ as a function of $\delta\mu$. We show this in Fig. 9 for a single value of $U=6$ and $\chi=50$. All other $U$ and $\chi$ look qualitatively similar. We see that there are two branches of VUTS solutions: insulating ($\langle n_i \rangle = 1$) and metallic ($\langle n_i \rangle > 1$) (equivalently differentiated by a zero and nonzero charge compressibility $\kappa = \partial\langle n_i \rangle/\partial\mu$, respectively). The insulating branch exists for $\delta\mu \le \delta\mu_1$, and the metallic branch for $\delta\mu \ge \delta\mu_2$: these two values $\delta\mu_1$ and $\delta\mu_2$ (with $\delta\mu_2 < \delta\mu_1$) are spinodal values limiting the metastability of the insulating and metallic solutions, respectively (see Appendix C for details). For each value of $\delta\mu$ the ground state is the branch with the lower energy. The energies of the two branches cross at a particular value, which we define to be $\delta\mu_c(\chi)$. At $\delta\mu_c(\chi)$, the ground state changes from insulating for $\delta\mu < \delta\mu_c$ to metallic for $\delta\mu > \delta\mu_c$. The occupation
<figure><img src="image_10.png" /><figcaption>FIG. 8. The extrapolated values of the staggered magnetization $m_s$ of the insulating state at half-filling, plotted as a function of $U$. Note the logarithmic scale used for $m_s$. The error bars are the discrepancy in the extrapolation with and without the last data point (for the points where they are absent, they are smaller than the data points).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 10. The extrapolated values of the jump in density $\delta n_c$ at the first-order transition occurring at $\delta\mu_c$ plotted as a function of $U$. The error bars, which are the discrepancy in the extrapolation with and without the last data point, are smaller than the data points.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 9. (a) The occupation versus $\delta\mu$ for $U = 6$ and $\chi = 50$ for the both metallic and insulating branches. The slope of these curves is the charge compressibility, which is another diagnostic of the phase. The metallic branch has a finite compressibility. All other values of $U$, $\chi$ look similar. (b) The energy per site of both branches in (a). Their crossing point $\delta\mu_c$ is indicated by a dashed line in both (a) and (b).</figcaption></figure>

undergoes a finite jump $\delta n(\delta\mu_c, \chi) = \langle n_i(\delta\mu_c, \chi) \rangle - 1$, indicating that this is a first-order metal-insulator transition. We estimate the size of the jump in the real system by the $\chi \to \infty$ extrapolated values, which are shown in Fig. 10. We can see they remain finite for all $U$, meaning that this is a true first-order transition, and not an artifact of finite bond dimension. Using a derivation based on the Maxwell construction (detailed in Appendix C), it can be shown that the total charge gap is given by $\Delta_c(\chi) = 2\delta\mu_c(\chi)$. In order to obtain the charge gap $\Delta_c$ for the real system, we extrapolate $\Delta_c(\chi)$ in $\chi$. The result is shown in Fig. 11. We can see that the extrapolated gap shows an exponential-like behavior at small $U$ similar to the one-dimensional case, though the precise behavior is difficult to extract numerically. At large $U$ the gap crosses over to a more linear dependence on $U$.

The first-order transition we observe implies that for every $U$ there is a range of forbidden density. Hence, in the $(U, n)$ plane the model exhibits phase separation. Our numerical method works in the grand canonical ensemble (we do not

<figure><img src="image_13.png" /><figcaption>FIG. 11. The extrapolated values of the charge gap $\Delta_c$ at half-filling plotted as a function of $U$. The error bars, which are the discrepancy in the extrapolation with and without the last data point, are smaller than the data points.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIG. 12. The staggered magnetization $m_s$ versus the chemical potential $\delta\mu$ for $U = 6$ and $\chi = 90$. The critical point $\delta\mu_c$ is indicated by a dashed line. We see the drop from $m_s > 0$ to $m_s = 0$, indicating the first-order transition from the antiferromagnetic insulator to the paramagnetic metal.</figcaption></figure>

fix particle number per unit cell), so we cannot observe this phase separation directly. However, VUTS, similar to other variational tensor-network methods like DMRG and VUMPS, can get “stuck” in local minima. We use this fact to find both branches of solutions near the transition point, even when they are metastable (i.e., not the lowest-energy states), by way of hysteresis in the numerical algorithm (see Appendix C for a detailed explanation). The resulting branches, shown in Fig. 9, can be continued to find the spinodal values of the first-order transition (see Appendix C).

The metallic ground state has no magnetic order for any value of density: we find that the staggered magnetization vanishes once $\delta\mu$ crosses $\delta\mu_c$. To illustrate the typical magnetization behavior we observe, in Fig. 12 we show the staggered magnetization $m_s$ as a function of $\delta\mu$ across $\delta\mu_c$ for $U = 6$ at a large but fixed $\chi = 90$. Once $\delta\mu$ becomes large enough to drive the system metallic, the staggered magnetization $m_s$ immediately drops to a very small value, which is zero within our error tolerance. All other values of $U$ and $\chi$ behave similarly. As $\chi$ increases, the small value of magnetization in the metal decreases further, although the behavior is not monotonic, as shown in Appendix C. Also, in Appendix B we describe the strategy we use to make sure we do not bias the magnetization of the metallic solution with our *Ansätze*.

It is interesting to compare our results on the $z = 3$ Bethe lattice to those established for the doping-driven MIT in the $z = \infty$ limit where DMFT becomes exact. Only a few studies [75–77] consider this transition while also taking into account phases with magnetic long-range order. As in our results, an antiferromagnetic insulator is found at half-filling for a range of chemical potentials, as well as a nonmagnetic metallic solution which can be stabilized for values of the chemical potential above a spinodal value $\delta\mu_2$. Furthermore, a magnetic metallic solution is found to exist in a narrow range of chemical potentials, which connects the magnetic insulator and the nonmagnetic metal. This may appear to differ from our findings, but it should be emphasized that all these studies consider only nonzero temperatures. As temperature is lowered, it is reported in Refs [75,76] that the magnetic metallic

<figure><img src="image_15.png" /><figcaption>FIG. 13. Our estimate for the truncation error $\epsilon_\rho$ in the metallic phase for $\langle n_i \rangle = 1.2$ (a) as a function of $\chi$ and (b) as a function of $U$ for the two largest $\chi$. Note that since (a) is a log-log plot, the linear form indicates an algebraic relationship.</figcaption></figure>

solution appears to exist only in an increasingly narrow interval of chemical potentials, and Ref. [75] suggested that at low temperature the MIT is a first-order transition between the magnetic insulator and the nonmagnetic metal, with a forbidden range of density corresponding to phase separation. Although, to the best of our knowledge, this has not yet been fully established directly at $T = 0$ for $z = \infty$, this conclusion is consistent with our findings on the finite coordination number lattice. In contrast, on fully frustrated $z = \infty$ lattices, which do not allow for long-range magnetic order (e.g., on the fully connected lattice with random hopping), it is established that the doping-driven MIT is second order at $T = 0$ and becomes first order only at finite temperature [13,78–80].

To conclude this section, we mention briefly the numerical accuracy of the data presented here. As noted in the Introduction, the numerical accuracy in tensor networks is generally measured by the truncation error, denoted by $\epsilon_\rho$. We show here the scaling of $\epsilon_\rho$ in the metallic state, which is the state with the largest entanglement and therefore the most computationally challenging for tensor networks. In Fig. 13 we plot our estimate of $\epsilon_\rho$ at a fixed density of $\langle n_i \rangle = 1.2$ as a function of $\chi$ for various $U$, and also as a function of $U$ at the largest values $\chi = 90, 100$. We can see that the decay with $\chi$ is algebraic, as expected for a gapless system. As a function of $U$, $\epsilon_\rho$ increases initially and then potentially saturates, although the large $U$ behavior is undetermined. Notably, we
can see that $\epsilon_\rho < 10^{-3}$ for all U and $\chi = 100$, which is a high level of accuracy. We discuss the behavior of $\epsilon_\rho$ at other points in the phase diagram in Appendix B.

## IV. QUASIPARTICLES

In this section we address the existence of quasiparticles in the system. We do this by computing the quasiparticle weight Z from the “momentum” distribution function of the ground state, which in turn is obtained from real-space correlation functions. One peculiarity of the Bethe lattice is that correlations between any two degrees of freedom sitting on individual nodes of the lattice have a maximal finite correlation length, even at criticality, due to the geometry of the lattice. However, algebraically decaying correlations reappear after a change of basis to single-particle states that are weighted sums of all nodes of a given generation emanating from a chosen center site. These bases of states unveil the traditional long-range criticality present in gapless states on the Bethe lattice. Below, we introduce a subset of these weighted states called the symmetric states, which we focus on in the rest of this section. From these symmetric states we define a quantum number which plays an analogous role on the Bethe lattice to quasimomentum on the hypercubic lattice, despite the absence of conventional translation invariance.

### A. Single-particle basis of symmetric states

For $U = 0$, the free-particle Hamiltonian was diagonalized in Ref. [81], using the symmetric set of single-particle states. These are given as follows. Choose any site to be labeled as the origin, with site label 0. Then consider all permutations of the nodes at each generation $l$ from the center. The symmetric states are those which are invariant under all such permutations. Their creation operators are given by

$$ \tilde{c}_{0,\sigma}^{\dagger} \equiv c_{0,\sigma}^{\dagger} \qquad (4) $$

and

$$ \tilde{c}_{l,\sigma}^{\dagger} \equiv \frac{1}{\sqrt{z(z-1)^{l-1}}} \sum_{\eta_1=0}^{z-1} \sum_{\eta_2 \neq \eta_1} \dots \sum_{\eta_l \neq \eta_{l-1}} c_{\eta_1+\eta_2+\dots+\eta_l, \sigma}^{\dagger} \quad (5) $$

for $l > 0$. The collection of $\eta_i$ denotes a unique path from the origin to the $l$th generation (this is the usual notation for nodes on the Bethe lattice). The state $\tilde{c}_{l,\sigma}^{\dagger}|\text{vacuum}\rangle$ is the symmetric combination of all the singly occupied spin-$\sigma$ states of the $l$th generation of the tree. These states form an orthonormal subset of all the states on the Bethe lattice, but for $U = 0$ they are the only relevant ones. In the symmetric state basis, the free-particle Hamiltonian maps onto fermions hopping on an infinite half-chain, with the first hopping amplitude equal to $\sqrt{z}$ and all the rest equal to $\sqrt{z-1}$ (remember $t$ has been set to one). The conjugate variable that replaces momentum is an angle $\theta \in [0, \pi]$, and the band energy is given by $\epsilon(\theta) = 2\sqrt{z-1} \cos \theta$. Note that the energy of a regular one-dimensional band is obtained by replacing $\theta$ with a momentum $k$ and $z$ with 2. The single-particle wave functions

$\psi_l(\theta)$ that diagonalize the Hamiltonian are given by

$$ \begin{aligned} \psi_0(\theta) &= \sqrt{\frac{2}{\pi}} \frac{\sqrt{z(z-1)} \sin(\theta)}{\sqrt{z^2 - 4(z-1)\cos^2(\theta)}}, \\ \psi_{l \neq 0}(\theta) &= \sqrt{\frac{2}{\pi}} \sin[l\theta + \gamma(z, \theta)], \end{aligned} \qquad (6) $$

$$ \gamma(z, \theta) = \begin{cases} \arcsin\left(\frac{z \sin(\theta)}{\sqrt{z^2 - 4(z-1)\cos^2(\theta)}}\right), & \theta \in [0, \frac{\pi}{2}) \\ \pi - \arcsin\left(\frac{z \sin(\theta)}{\sqrt{z^2 - 4(z-1)\cos^2(\theta)}}\right), & \theta \in [\frac{\pi}{2}, \pi]. \end{cases} $$

Once $U \neq 0$, the symmetric states are no longer enough to describe the system. Indeed, if we simply consider a Mott-type state with one particle sitting on each of the sites at generation $l = 1$ away from the center site, that state cannot be written using only the single-particle symmetric states associated with the same center site. In general, an arbitrary multiparticle Fock state cannot be constructed as a tensor product of the symmetric single-particle states. Therefore, to construct it one must employ states from other symmetry sectors. The interacting ground states we find numerically in this work therefore contain states from various symmetry sectors. However, excitations above the ground state can occur in any of these sectors, and we do not have to consider all of them. In order to tractably answer the question of existence of quasiparticles, we choose to focus on excitations in the symmetric sector. How quasiparticles in different symmetry sectors are related to each other is an interesting question for future work [82].

### B. $\theta$-distribution function for $U = 0$

The $\theta$-distribution function is calculated from the equal-time correlation functions of symmetric single-particle excitations $\langle\tilde{c}_{0,\sigma}^{\dagger}\tilde{c}_{l,\sigma}\rangle$, where 0 labels a chosen center site and $l$ labels the generation away from the center site. These can be computed as

$$ \begin{aligned} \langle \tilde{c}_{0,\sigma}^{\dagger} \tilde{c}_{l,\sigma} \rangle &= \int_{-\infty}^{0} d\omega \mathcal{A}_{0l}(\omega) \\ &= \int_{-2\sqrt{z-1}}^{\mu} d\epsilon \sqrt{\frac{z}{2\pi}} \frac{\psi_l\left[\arccos\left(-\frac{\epsilon}{2\sqrt{z-1}}\right)\right]}{\sqrt{z^2 - \epsilon^2}}, \end{aligned} \qquad (7) $$

where $\mathcal{A}_{0l}(\omega)$ is the probability of inserting an electron with frequency $\omega$ at the center site and observing it at generation $l$ at the same frequency. The occupation function in $\theta$ space is defined as

$$ n_{\sigma}(\theta, \theta') \equiv \langle \tilde{c}_{\theta, \sigma}^{\dagger} \tilde{c}_{\theta', \sigma} \rangle, \qquad (8) $$

where

$$ \tilde{c}_{\theta, \sigma} = \lim_{L \to \infty} \sqrt{\frac{\pi}{L+1}} \sum_{d=0}^{L} \psi_d(\theta) \tilde{c}_{d, \sigma} \qquad (9) $$

are the $\theta$ transforms of the symmetric state operators $\tilde{c}_{d,\sigma}$. The key difference with the usual calculation in hypercubic lattices is that $\langle\tilde{c}_{d,\sigma}^{\dagger}\tilde{c}_{d',\sigma}\rangle$ is not simply a function of $|d-d'|$, due to the fact that the symmetric states are defined relative to a chosen center site. This is illustrated in Fig. 14, where we show that $\langle\tilde{c}_{2,\sigma}^{\dagger}\tilde{c}_{4,\sigma}\rangle$ contains correlations of length 2, 4,
<figure><img src="image_16.png" /><figcaption>FIG. 14. Four generations of the $z = 3$ infinite Bethe lattice emanating from a given center site. The inner shell is at generation $d = 2$, and the outer shell is at generation $d' = 4$. Choosing a given site on the inner shell (colored red), there are three different groups of sites on the outer shell (colored orange, green, and purple) that contribute different-length correlations.</figcaption></figure>

and 6. Therefore, $n_σ(θ, θ')$ is not diagonal in $θ$. One way to understand this is to think of the entire symmetric sector parametrized by $θ$ as the $k = 0$ space on the hypercubic lattice, i.e., the one that is fully symmetric under translations. However, within this sector there is no additional symmetry of the Bethe lattice that requires the total $θ$ to be conserved in a scattering process, and therefore $n_σ(θ, θ')$ is not diagonal [82]. Carefully inserting Eq. (9) into Eq. (8) and rewriting the result in terms of $\langle \tilde{c}_{0,σ}^{\dagger} \tilde{c}_{l,σ} \rangle$ gives

$$n_{\sigma}(\theta, \theta') = \lim_{L \to \infty} \frac{\pi}{L+1} \sum_{d,d'=0}^{L} \psi_d(\theta) \psi_{d'}(\theta') \frac{z-2}{\sqrt{z(z-1)}} \\ \times a(z)^{\delta_{\min(d,d')}},_0 \sum_{r=0}^{\min(d,d')} b(z)^{\delta_{r,0} + \delta_{r,\min(d,d')}} \\ \times c(z)^{\delta_{d,d'} \delta_{r,0} - \delta_{d,0} \delta_{d',0}} \langle \tilde{c}_{0,\sigma}^{\dagger} \tilde{c}_{|d-d'|+2r,\sigma} \rangle, \quad (10)$$

where

$$a(z) = \frac{\sqrt{z(z-2)}}{(z-1)^{3/2}}, \quad (11)$$

$$b(z) = \frac{z-1}{z-2}, \quad (12)$$

$$c(z) = \sqrt{\frac{z}{z-1}}, \quad (13)$$

and $\delta_{ij}$ is the Kronecker delta function. Details of the derivation of Eq. (10) are in Appendix D. We plot the exact $n_σ(θ, θ')$ for $U = 0$ at half-filling in Fig. 15(a).

In this work, we focus exclusively on the diagonal component of the occupation function $n_σ(θ) = n_σ(θ, θ)$. This tells us the occupation of excitations that preserve total $θ$ when scattering with each other. We leave the detailed study of the full occupation function $n_σ(θ, θ')$ for future work. We plot the exact $n_σ(θ)$ in Fig. 15(b) for $U = 0$ and densities $\langle n_i \rangle = 1, 1.2$. We can see the expected behavior of $n_σ(θ) = \Theta(θ_F - θ)$, where $\Theta(x)$ is the Heaviside step function and $θ_F$ is the $θ$ analog

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>FIG. 15. (a) The occupation function $n_σ(θ, θ')$ for the half-filled case at $U = 0$. Aside from the expected step function along the diagonal, there is nontrivial off-diagonal structure. (b) The diagonal component $n_σ(θ)$ for $U = 0$ and densities $\langle n_i \rangle = 1, 1.2$, corresponding to values of $θ_F = π/2$ and $θ_F \approx 1.81$, respectively. Calculating the occupation from Eq. (10) requires a large distance cutoff $L$, which is chosen here to be (a) $L = 100$ and (b) $L = 200$. This introduces an artificial correlation length, causing the step function to be (slightly) smoothed out.</figcaption></figure>

of the “Fermi momentum” [the step function in Fig. 15(b) is slightly smoothed out due to a finite $L$ in Eq. (10)].

In order to compute the correlation function $\langle \tilde{c}_{0,σ}^{\dagger} \tilde{c}_{l,σ} \rangle$ from the VUTS numerical solution, we use the fact that even in the interacting ground state, $\langle c_{i,σ}^{\dagger} c_{j,σ} \rangle$ is still only a function of the distance $|i - j|$ (and $σ$). We can then compute $\langle c_{0,σ}^{\dagger} c_{l,σ} \rangle$ for an arbitrary branch to write

$$\langle \tilde{c}_{0,\sigma}^{\dagger} \tilde{c}_{l,\sigma} \rangle = \sqrt{z^{1-\delta_{l,0}}(z-1)^{l-1+\delta_{l,0}}} \langle c_{0,\sigma}^{\dagger} c_{l,\sigma} \rangle. \quad (14)$$

Note that the difference here from the previous considerations of this paragraph is that one of the reference points has been set to the center site. We plot $\langle \tilde{c}_{0,σ}^{\dagger} \tilde{c}_{l,σ} \rangle$ measured with VUTS in Fig. 16 for $U = 0$ and densities $\langle n_i \rangle = 1, 1.2$, along with the exact results. The finite-$\chi$ results are close to the exact ones, but, as expected, the correlation function at large enough distances decays exponentially with a finite correlation length $\xi(\chi)$. We can measure $\xi(\chi)$ by fitting $\langle \tilde{c}_{0,σ}^{\dagger} \tilde{c}_{l,σ} \rangle$ to an exponential at large $l$. The results are shown in Fig. 16(c), where the polynomial fit gives $\xi(\chi) \sim \chi^{0.84}$ for both densities.

When calculating $n_σ(θ)$ from Eq. (10), in practice we must choose a finite value of $L$. As long as we take $L$ large enough,
<figure><img src="image_19.png" /><figcaption>FIG. 16. The correlation function $\langle \tilde{c}_{0,\uparrow}^{\dagger} \tilde{c}_{r,\uparrow} \rangle$ ($\sigma = \downarrow$ gives the same) for densities (a) $\langle n_i \rangle = 1$ and (b) $\langle n_i \rangle = 1.2$ for the non-interacting case $U = 0$. We show the $\chi = 50, 100$ results along with the exact solution of Eq. (7). Also shown are the functions $\frac{1}{r} \psi_r(\theta_F) \psi_0(\theta_F)$ with (a) $\theta_F = \pi/2$ and (b) $\theta_F \approx 1.81$, which are excellent fits beyond a short-distance scale, showing Friedel oscillations due to the Fermi-surface singularity. The insets show that the $1/r$ decay is fit perfectly over a larger distance by the exact solution, while the finite-$\chi$ solutions display an exponential decay at large distances. (c) The correlation length extracted from the long-distance behavior of $\langle \tilde{c}_{0,\sigma}^{\dagger} \tilde{c}_{r,\sigma} \rangle$, plotted versus $1/\chi$. The slope for both densities is $\xi(\chi) \sim \chi^{0.84}$.</figcaption></figure>

the correlation length $\xi(\chi)$ will act as the long-distance cutoff and the value of $L$ will not have any effect. Our results are obtained using bond dimensions up to $\chi = 100$, for which the induced correlation length is $\xi(\chi) \lesssim 40$. We find that a value of $L \sim 400$ is large enough for all $\chi$ we study. The

<figure><img src="image_20.png" /><figcaption>FIG. 17. Plots of $n_\sigma(\theta)$ near $\theta_F$ for $U = 0$ at density $\langle n_i \rangle = 1.2$ and a range of $\chi$. The value of $L$ used here from Eq. (10) is $L = 400$.</figcaption></figure>

finite $\xi(\chi)$ smoothes out the step function in $n_\sigma(\theta)$, so we estimate $\theta_F$ from the location of the maximum of $|n'_\sigma(\theta)|$. In Fig. 17 we show $n_\sigma(\theta)$ for $\langle n_i \rangle = 1.2$ near $\theta_F$. The finite slope at $\theta_F$ diverges as a power law in $\chi$, as we show in Fig. 18. This indicates that $\xi(\chi)$ is the only low-energy scale in the problem, and the state is truly gapless in the $\chi \to \infty$ limit.

In summary, in this section we show how to compute the diagonal part of the occupation function in the noninteracting case, using VUTS and a careful extrapolation in the bond dimension, and achieve excellent agreement with the analytical result.

### C. Quasiparticles in the interacting system

Now, we turn on interactions. In Fig. 19 we plot $n_\sigma(\theta)$ near $\theta_F$ for $\langle n_i \rangle = 1.2$ and $U = 5$ at various $\chi$ as well as for various $U$ at $\chi = 70$. The occupation shows a form similar to that of the free case, albeit with a reduced size of the step at $\theta_F$. The value of the slope at $\theta_F$ diverges for all $U$, as we show in Fig. 18. The interacting system is therefore also gapless in the $\chi \to \infty$ limit, as expected.

The quasiparticle weight $Z$ of the symmetric state excitations is defined as

$$Z = \left( \lim_{\theta \to \theta_F} - \lim_{\theta \to \theta_F^\dagger} \right) n_\sigma(\theta). \quad (15)$$

<figure><img src="image_21.png" /><figcaption>FIG. 18. The value of $|n'_\sigma(\theta_F)|$ vs $1/\chi$ for various $U$ at density $\langle n_i \rangle = 1.2$. The straight lines on the log-log plot are fit by $|n'_\sigma(\theta_F)| \sim \chi^\alpha$, with $0.5 < \alpha < 1$.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>FIG. 19. Plots of $n_{\sigma}(\theta)$ near $\theta_{F}$ at density $\langle n_{i} \rangle = 1.2$, showing the dependence (a) on $\chi$ for fixed $U = 5$, and (b) on $U$ for fixed $\chi = 70$. The value of $L$ from Eq. (10) is $L = 400$.</figcaption></figure>

For a Fermi liquid $n_{\sigma}(\theta)$ has a step at $\theta_{F}$ and $Z > 0$, while for a Luttinger liquid the occupation function has a higher-order nonanalyticity that scales as $n_{\sigma}(\theta - \theta_{F}) \sim |\theta - \theta_{F}|^{\gamma}\text{sgn}(\theta - \theta_{F})$ for some $\gamma < 1$ and therefore $Z = 0$. Our goal is to find the true thermodynamic value of $Z$ to distinguish between these two scenarios. Of course, for a finite $\chi$, Eq. (15) will always give zero. However, we can define a quantity $Z(\chi)$ whose limit will give $Z$ in the $\chi \to \infty$ limit. We define this as

$$Z(\chi) = n_{\sigma} \left( \theta_{F} - \frac{\pi}{2\xi(\chi)} \right) - n_{\sigma} \left( \theta_{F} + \frac{\pi}{2\xi(\chi)} \right), \quad (16)$$

which satisfies the desired property because $\xi(\chi) \to \infty$ with increasing $\chi$. We choose a spacing of $\Delta\theta = \pi/\xi(\chi)$ around $\theta_{F}$ because that is roughly the resolution one expects from a finite correlation length, and therefore the convergence in $1/\chi$ should be fastest. We plot $Z(\chi)$ vs $\chi$ for $\langle n_{i} \rangle = 1.2$ and a range of $U$ in Fig. 20. The results show that the extrapolated $Z$ is (a) very close to the expected value of $Z = 1$ for the free theory and (b) finite for all $U$ we study. We plot $Z$ as a function of $U$ in Fig. 21(a), where we can see that it decreases as a function of $U$, but seems to saturate to a finite value. The saturation value is an increasing function of doping $\langle n_{i} \rangle - 1$, as illustrated by Fig. 21(b) where we plot the value of $Z(U = 20)$ as a function of doping.

We also address the question of Luttinger's theorem for $n_{\sigma}(\theta)$. We check that $\theta_{F}$ is independent of $U$ (the dependence on $\chi$ is negligible) for various densities in the range $\langle n_{i} \rangle \in$

<figure><img src="image_23.png" /><figcaption>FIG. 20. $Z(\chi)$ at density $\langle n_{i} \rangle = 1.2$, shown with linear fits.</figcaption></figure>

(1.1, 1.4). From this we conclude that Luttinger's theorem holds for all values of the density and interaction strength.

The decrease of $Z$ with increasing $U$ as well as the saturation at large $U$ to a value which increases with doping are both qualitatively consistent with results established in the $z = \infty$ (DMFT) limit and with slave-boson approaches

<figure><img src="image_24.png" /><figcaption>FIG. 21. (a) The extrapolated values of $Z$ at densities $\langle n_{i} \rangle = 1.2, 1.3, 1.4$, plotted as a function of $U$. The first two blue points are covered by the orange ones. The error bars, which are the discrepancy in the extrapolation with and without the last data point, are smaller than the data points. We can see that the curves are close to saturation at $U = 20$. (b) The values of $Z$ at $U = 20$, which seem close to the saturated values, as a function of $\langle n_{i} \rangle - 1$.</figcaption></figure>
(for general lattices) [13,83]. A distinctive aspect of these theories, however, is that the effective mass of quasiparticles is related to Z by $m^*/m = 1/Z$. On the technical level, this is due to the locality of the self-energy, while physically this reflects the inability of these approaches to capture the feedback of short-range order and collective modes into the physics of quasiparticles. In contrast, in the present case, since the connectivity is kept finite, we would expect this feedback to be present. It is therefore an outstanding question for future work to explore whether the dispersion of quasiparticles is renormalized in a different manner than Z itself, and in particular whether it is affected by short-range antiferromagnetic correlations at low doping level. This is left for future work since it requires an extension of our algorithm to the study of excited states.

V. DISCUSSION

In summary, we have introduced a numerical algorithm, (fermionic) VUTS, to study quantum (fermionic) models on the Bethe lattice. We apply it to the Hubbard model for coordination number $z = 3$, allowing for a two-site unit cell, obtain the $T = 0$ phase diagram, and study the doping-induced Mott transition. We find an antiferromagnetic insulating phase at half-filling and a paramagnetic metallic phase for the doped system, which are separated by a first-order insulator to metal phase transition. The model displays phase separation at low doping, with a range of forbidden densities. These conclusions were reached by allowing for a two-site unit cell. As discussed in the Introduction, we expect that this size unit cell is enough to capture all the correct phases and physics of the model. However, we cannot exclude the possibility that phases with more complex charge or magnetic ordering exist when allowing for a larger unit cell, which we leave as an open question for future work. By studying the diagonal component of the occupation function for momenta of the symmetric single-particle sector, we find that the quasiparticle weight is nonzero, consistent with the existence of a Fermi-liquid ground state for fermions on the Bethe lattice. We find that this Fermi-liquid state obeys Luttinger's theorem.

An interesting direction in which to extend this work would be to further characterize this Fermi-liquid state. One would like to know, for example, what happens near $\theta_F$ to the off-diagonal $n_\sigma(\theta, \theta')$ when interactions are turned on. It is also interesting to look at some of the other symmetry sectors, say the ones that leave each of the $z$ subtrees connected to the center site invariant, and see whether they have quasiparticles and how the quasiparticle weight depends on the sector. These questions are the Bethe lattice version of the important physical question of “momentum dependence” of quasiparticle properties on the Fermi surface of hypercubic lattices.

In the one-dimensional VUMPS algorithm, it has been shown that low-lying excitations above the ground state can be accurately computed [84]. An obvious question is how to extend these ideas to VUTS (this is another potential advantage of this method over imaginary-time evolution). Extension of the algorithm to excited states would allow to characterize the effective mass (dispersion) of the quasiparticles, paving the road for a study of how short-range (e.g., antiferromagnetic) correlations affect quasiparticle properties. This is a

very important question for the physics of strongly correlated electron systems. A breakdown of the $m^*/m = 1/Z$ relation would signal that this feedback is indeed present, in contrast to the infinite connectivity limit. Finally, studying the energy dependence of the quasiparticle lifetime, as well as the interactions between quasiparticles, would be a comprehensive study of Landau Fermi-liquid theory on the Bethe lattice.

It would also be interesting to look at the entanglement structure of the Fermi-liquid state, as compared to a Luttinger liquid. Since the entanglement spectrum is easily obtained on the one-dimensional and Bethe lattices from the singular value decomposition of a single bond tensor, this question could in principle be easily answered.

We also propose that the finite-$z$ Bethe lattice can be used as a computationally tractable platform for the study of how quasiparticles can be destroyed and Fermi-liquid behavior breaks down when considering other fermionic Hamiltonians on this lattice. One route to explore this, which connects to the feedback of long-wavelength collective modes or short-range spatial correlations on quasiparticle properties, is to study the vicinity of a quantum critical point. Another route is to study microscopic models that are tailor engineered to have incoherent excitations, such as the one of Ref. [85], or the multichannel Kondo lattice.

Studying the interplay between frustration and strong correlations is another promising direction for future research. Introducing a frustrating next-nearest-neighbor hopping term has been found with DMFT to yield an interaction-driven metal-insulator transition on the infinite-$z$ Bethe lattice [13,86]. An interesting question is whether this transition can also be found at finite $z$. The VUTS algorithm could also be extended to other lattices with a treelike structure, such as the Husimi cactus. The study of spin models on such lattices has revealed spin-liquid ground states [74,87] (see also Ref. [88] for considerations on the spin-ice model), opening the question of how these models behave upon doping.

Finally, increasing the temperature to a nonzero value is another interesting direction. Intuitively, some finite-temperature properties may be less sensitive to the differences between tree lattices and two- or three-dimensional hypercubic lattices. This could be done using the purification method [89–91] which has been formulated on the Bethe lattice in Ref. [92].

ACKNOWLEDGMENTS

We thank J. Agier, M. Claassen, M. Ferrero, G. Kotliar, C. Laumann, S.-S. Lee, R. Moessner, M. Rozenberg, S. White, and in particular M. Eckstein and R. Rossi for useful discussions. We thank R. Verresen for comments on the manuscript. All tensor-network calculations and the VUTS code were implemented with the ITensor Library (C++ version 3.1) [93]. The Flatiron Institute is a division of the Simons Foundation.

APPENDIX A: DETAILS OF THE VARIATIONAL UNIFORM TREE STATES ALGORITHM

1. Summing the Hamiltonian terms

Here we explain how we compute projected Hamiltonians for subtrees of the Bethe lattice, using infinite summations like
<figure><img src="image_25.png" /><figcaption>FIG. 22. (a) Transfer matrix $T_i^{n,m}$ from Eq. (A1). (b) Diagrammatic equations for the right eigenvectors $R_{i,m}$ of the transfer matrices, which represent the reduced density matrices for bipartitions of the Bethe lattice [and are also shown in Eq. (A3)]. As a reminder, the bonds labeled $k, m, n$ have link indices $l_k, l_m, l_n$, respectively, and we leave of the prime levels for simplicity.</figcaption></figure>

those shown in Fig. 5. To help us with the computation, we start by defining the transfer matrices, which are operators that act on the virtual degrees of freedom of the network. We build the transfer matrices out of the $A_{i,m}$ tensors, the tensors of the TTN that are in the orthonormal gauge. They are defined as

$$T_i^{n,m} = \sum_{s_i, l_k} \tilde{A}_{i,m}^{s_i, l'_m, l'_n, l_k} A_{i,m}^{s_i, l_m, l_n, l_k}, \quad (A1)$$

where $n, m = 0, 1, 2$ and $n \neq m$. The diagrammatic version of Eq. (A1) is shown in Fig. 22(a). Using the transfer matrix we can rewrite the gauge conditions of Eq. (1) (Fig. 3) as

$$(\mathbb{1}_{\tilde{i},n}) |T_i^{n,m}| = (\mathbb{1}_{i,m}) |, \quad (A2)$$

where we use the vectorized notation such that $|\mathbb{1}_{i,m}|$ is a length $\chi^2$ vector that is isomorphically equivalent to the $\chi \times \chi$ matrix $\mathbb{1}_{i,m}$, and the transfer matrix $T_i^{n,m}$ is treated as a $\chi^2 \times \chi^2$ matrix mapping from the bra and ket links of bond $n$ to the bra and ket links of bond $m$. We also use the notation $\bar{0} = 1, \bar{1} = 0$ for the unit-cell indices introduced in the main text. In addition, we make use of the right eigenvectors of the transfer matrices $R_{i,m}$. These are defined by the set of equations

$$\sum_{s_i, l_k, l_m, l'_m} \tilde{A}_{i,m}^{s_i, l'_m, l'_n, l_k} R_{i,m}^{l'_m, l_m} A_{i,m}^{s_i, l_m, l_n, l_k} = R_{i,n}^{l'_n, l_n}, \quad (A3)$$

or in vectorized notation in terms of the transfer matrices

$$T_i^{n,m} |R_{i,m}) = |R_{\tilde{i},n}). \quad (A4)$$

The matrices $R_{i,m}$ are the reduced density matrices of a bipartition of the state. The diagrammatic versions of Eqs. (A3) and (A4) is shown in Fig. 22(b).

Using vectorized notation, the environment matrices $H_{i,m}$ are calculated from the geometric series

$$(\mathbf{H}_i) = [(\mathbf{h}_i) + (\mathbf{h}_{\tilde{i}} | \mathbf{T}_i)] \sum_{k=0}^{\infty} (\mathbf{T}_{\tilde{i}} \mathbf{T}_i)^k, \quad (A5)$$

where $(\mathbf{H}_i|$ and $(\mathbf{h}_i|$ are vectors of vectorized matrices

$$(\mathbf{H}_i| = ((H_{i,0} | \quad (H_{i,1} | \quad (H_{i,2} |), \\
\times 0.0pt(\mathbf{h}_i| = ((h_{i,0} | \quad (h_{i,1} | \quad (h_{i,2} |), \quad (A6)$$

and $\mathbf{T}_i$ is a matrix of transfer matrices

$$\mathbf{T}_i = \begin{pmatrix} 0 & T_i^{0,1} & T_i^{0,2} \\ T_i^{1,0} & 0 & T_i^{1,2} \\ T_i^{2,0} & T_i^{2,1} & 0 \end{pmatrix}. \quad (A7)$$

Here, the $|h_{i,m}|$ are vectorized version of the tensors $h_{i,m}$, which are the projections of a local Hamiltonian term $h_{i,j}$ into the orthonormal basis of $A_{i,m}$, and are defined diagrammatically in Fig. 5(b). We remind readers that the transfer matrices of the $z = 3$ Bethe lattice are $T_i^{m,n}$ for $m, n = 0, 1, 2$ and $m \neq n$. We start by analytically summing Eq. (A5), using the fact that it is of the form of a geometric series [in other words, using that $\sum_{k=0}^{\infty} x^k = (1-x)^{-1}$, where $x = \mathbf{T}_{\tilde{i}}\mathbf{T}_i$ is the product of the two matrices of transfer matrices]. However, we also have to subtract away the infinite number of local energy contributions (one from each term in the series) which cause the series to diverge. This results in the following linear set of equations for $(H_{i,m})$:

$$(H_{i,m}| - \sum_{m',m''} (H_{i,m''}| T_i^{m'',m'} T_i^{m',m} | \mathbb{1} - |R_{i,m}) (\mathbb{1} - |R_{i,m}|)) \\
= \left( (h_{i,m}| + \sum_{m'} (h_{\tilde{i},m'}| T_i^{m',m} | \mathbb{1} - |R_{i,m}) (\mathbb{1} - |R_{i,m}|) \right). \quad (A8)$$

The sums of $m'$ and $m''$ run over 0,1,2 (for each direction in the $z = 3$ Bethe lattice), but we remind readers that $T_i^{m,n}$ only exists for $m \neq n$. The terms proportional to $|\mathbb{1}_{i,m}|$ on the left-hand side of Eq. (A8) act to project out the dominant eigenspace of the transfer matrix that contains the infinite energy contribution, and lead to a convergent geometric series. Congruently, the terms proportional to $|\mathbb{1}_{i,m}|$ on the right-hand side of Eq. (A8) zero out the onsite energies of the unit cell. Equation (A8) is a generalization of the result derived in Appendix D of Ref. [55] to the $z = 3$ Bethe lattice with a two-site unit cell. We now can think of the solution to Eq. (A8) as defining $(H_{i,m}|$, rather than defining it via the geometric series of Eq. (A5). Applying the vectorized density matrices $|R_{i,m})$ to Eq. (A8), we see that $(H_{i,m}|$ satisfies $(H_{i,m}|R_{i,m}) = 0$, so the infinite-energy shift can be seen as shifting the total energy of the environment to zero. We can see that the only unknowns in Eq. (A8) are the three environment matrices $(H_{i,m}|$. Also note that Eq. (A8) is in the form of a linear equation $(\mathbf{x}|\mathbf{A} = (\mathbf{b}|,
<figure><img src="image_26.png" /><figcaption>FIG. 23. (a) Fermionic ordering of the legs of the $A_{i,m}$, $A_{i,C}$, $C_m$ tensors. (b) An example of a tensor diagram containing fermionic swap gates. The fermionic swap gates are labeled by orange diamonds, and are present at the crossing of any two legs. Note that after the contraction of $A_{0,2}$, $\bar{A}_{0,2}$, $h$, and $A_{1,C}$, we have to put the dangling tensor legs back into the order of $A_{1,C}$, which requires three more swap gates.</figcaption></figure>

where **A** and (b| are known and (x| is an unknown vector of the environment matrices:

$$ ((H_{i,0}| \quad H_{i,1}| \quad H_{i,2}|) \begin{pmatrix} A_{0,0} & A_{0,1} & A_{0,2} \\ A_{1,0} & A_{1,1} & A_{1,2} \\ A_{2,0} & A_{2,1} & A_{2,2} \end{pmatrix} = ((b_0| \quad b_1| \quad b_2|), \quad (A9) $$

where the elements $A_{m,n}$ and ($b_n$| can be read off from Eq. (A8). In practice, we solve the set of equations defined by Eq. (A8) for $m = 0, 1, 2$ [or the equivalent Eq. (A9)] using an iterative solver (for this work we use GMRES, but others may be employed).

## 2. Fermions

As explained in the main text at the end of Sec. II, for fermionic models like the Hubbard model, we need to use a fermionic version of VUTS. We use the method outlined in Refs. [52,59]. Every tensor is endowed with a fermion parity $Z_2$ quantum number and is parity preserving. When two tensor legs cross on a planar projection of a tensor diagram, a fermionic swap gate is placed at the crossing. In order to employ this method, we need to use a fixed ordering convention for the legs of the tensors $A_{i,m}$, $A_{i,C}$, $C_m$, which must be kept consistent in all of the diagrams in the calculation. The convention that we choose is shown in Fig. 23(a). The legs of the tensors always appear in this order for tensor diagrams, regardless of their location in the diagram. An example of a diagram involving swap gates is shown in Fig. 23(b), where we draw a single term contributing to the action of the one-site projected Hamiltonian on the one-site wave function $A_{1,C}$ (the fermionic version of Fig. 6). Please note that care has to be taken in how the diagrams are drawn, and also that the site equivalence of the system is properly obeyed. For more details on the method, such as additional diagrammatic rules that must be followed to account for the fermion signs, we refer the reader to Refs. [52,59].

# APPENDIX B: DETAILS OF THE NUMERICAL CALCULATION

## 1. Convergence strategies

In this Appendix we outline some of the numerical ingredients used in obtaining the phase diagram. For a generic value of $U$ and $\delta\mu$, we start the VUTS algorithm with a

random Ansatz of small bond dimension $\chi = \chi_i \sim 4$. We then increase $\chi$ up to the desired $\chi_f$ using a generalization of the scheme explained in Appendix B of Ref. [55]. In most of the phase diagram, any update schedule, i.e., any way of getting from $\chi_i$ to $\chi_f$, leads to VUTS converging to the same final state, within the precision error $\epsilon_{prec}$. We choose $\epsilon_{prec} \lesssim 10^{-6}$ for all our calculations. However, near the phase transition $\delta\mu_c$, more care is needed in selecting the update schedule. Generically, we find that increasing $\chi$ in small increments and running several iterations after each small increase favors the insulating state, while using larger increments and a small number of iterations to get from $\chi_i$ to $\chi_f$ favors the metallic state. In this way, both branches near $\delta\mu_c$ can be obtained. Unfortunately, a given update schedule is never guaranteed to converge to a specific state, or to converge at all. Therefore, once the insulating/metallic state is found at a given $\delta\mu_0$, we can use that state as an Ansatz for the state at $\delta\mu \approx \delta\mu_0$, in order to favor the phase of the $\delta\mu_0$ state. In this way, both branches can be reliably obtained for all $\delta\mu \approx \delta\mu_c$ (see Appendix C for a discussion about continuing both branches into the regions of metastability).

The metallic ground state we find has vanishing magnetization, within our precision. In order to make sure we do not accidentally disfavor a magnetic metallic state by our choice of initial state, we do separate computations biased toward a magnetic metal by turning on a staggered external magnetic field $H \to H - B_s(S_{(0)}^z - S_{(1)}^z)$, where the subscripts label the location in the unit cell. Then, we use this state as an Ansatz for the metallic state close to $\delta\mu_c$ in the case when $B_s = 0$ (making sure that the density of the Ansatz state is close to what we expect in the $B_s = 0$ case). Even with this method, we do not find any energetically favorable magnetic metal.

## 2. Truncation error and run time

Here we address the important question of numerical accuracy and run time. In tensor-network states, the accuracy is measured by the truncation error or discarded weight $\epsilon_\rho$. In this work, we use a proxy for $\epsilon_\rho$, which is shown to be of the same order and used in Ref. [55], where it is referred to as $\|B_2\|^2$ (we compute $\|B_2\|^2$ using a trick similar to what is used in Ref. [32], in order for the computation to be manageable). Here we simply refer to $\|B_2\|^2$ as $\epsilon_\rho$. In Fig. 24 we show the dependence of $\epsilon_\rho$ on $\chi$ and $U$ in the insulating phase ($\epsilon_\rho$ is essentially independent of $\delta\mu$ in the insulator). We can see
<figure><img src="image_27.png" /><figcaption>FIG. 24. Our estimate for the truncation error $\epsilon_{\rho}$ in the insulating phase (a) as a function of $\chi$ for several $U$ and (b) as a function of $U$ for the largest $\chi = 100$.</figcaption></figure>

that (a) $\epsilon_{\rho}$ decays with $\chi$ roughly as an exponential for all $U$ and (b) $\epsilon_{\rho}$ decays with $U$ also roughly as an exponential.

In the metallic phase, $\epsilon_{\rho}$ has a non-negligible dependence on the density, so for the purposes of comparison we look at densities $\langle n_i \rangle = 1.2$ and the smallest density attainable for a given $U$ and $\chi$, $\langle n_i(\delta\mu_c, \chi) \rangle$. In Fig. 25 we plot $\epsilon_{\rho}$ for the smallest density as a function of $\chi$ for various $U$, and also as a function of $U$ at the largest $\chi = 100$ (the corresponding plots for $\langle n_i \rangle = 1.2$ are in Fig. 13). We can see that, as in Fig. 13, the decay with $\chi$ is algebraic, as expected for a gapless system. As a function of $U$, $\epsilon_{\rho}$ increases initially and then seems to saturate, although the large $U$ behavior is undetermined.

To conclude, we note that for the largest bond dimension we present, $\chi = 100$, $\epsilon_{\rho} < 10^{-3}$ for every point in parameter space we study. In order to improve the run time, we use multithreading over the block sparse tensors, which speeds up the calculation significantly. The largest bond dimension we present in this work is $\chi = 100$. We thread over 12, 24, and 48 cores, for which the run times for a single VUTS iteration at $\chi = 100$ are 24.2, 18.5, and 15 min, respectively. The number of iterations required to achieve a precision error of $\epsilon_{\text{prec}} \lesssim 10^{-6}$ depends significantly on the proximity of the Ansatz wave function to the true ground state, and could be as low as tens of iterations and as large as more than 1000 iterations. Therefore, it is very important to have good convergence strategies, as described in the previous section.

## APPENDIX C: SUPPORTING PLOTS FOR PHASE DIAGRAM AND DISCUSSION OF FIRST-ORDER TRANSITION

Here we provide more supporting plots of finite-$\chi$ results from Sec. III, and also discuss in more detail the first-order metal-insulator phase transition, including the calculation of the charge gap. We start by illustrating the behavior of the density as a function of $\delta\mu$. We show this in Fig. 26 for $U=2, 5$ and a range of $\chi$. All other values of $U$ and $\chi$ behave qualitatively similar. Now, in Fig. 27 we show several finite-$\chi$ plots, the extrapolations of which are in the main text. In Figs. 27(a)–27(c) we plot the charge gap, staggered magnetization at half-filling, and the jump in occupation at the phase transition, respectively, as a function of $\chi$ for a range of $U$. We fit the data with a linear fit, ignoring the first few points, after which the plots’ data look linear. The extrapolated values are shown in Figs. 11, 8, and 10, respectively. In Fig. 27(d) we plot the finite-$\chi$ values of the staggered magnetization on the metallic side of the phase transition for some values of $U$ (other values of $U$ behave similarly). We can see the (nonmonotonic) trend of $m_s(\chi)$ decreasing to arbitrarily small values as $\chi$ increases. For completeness we also provide the scaling of the correlation length $\xi(\chi)$ with $\chi$ for a range of $U$ at density $\langle n_i \rangle = 1.2$, shown in Fig. 28.

Now we turn to the question of the first-order metal-insulator transition. Near the transition point $\delta\mu_c$, both the insulating and metallic branches exist and can be stabilized for a range of $\delta\mu$, as can be seen in Fig. 9. The metastable parts of the branches are obtained by way of “hysteresis” in

<figure><img src="image_28.png" /><figcaption>FIG. 25. Our estimate of the truncation error $\epsilon_{\rho}$ in the metallic phase at density $\langle n_i(\delta\mu_c, \chi) \rangle$.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>FIG. 26. Fermion occupation $\langle n_i \rangle$ as a function of $\delta\mu$ for (a) $U = 2$ and (b) $U = 5$, and a range of $\chi$.</figcaption></figure>

the numerical algorithm. We start, e.g., with the insulating ground state at a given $\delta\mu_x < \delta\mu_c$ and use that as an *Ansatz* for $\delta\mu_y > \delta\mu_c$, then use the insulating solution thus found at $\delta\mu_y$ as an *Ansatz* for $\delta\mu_z > \delta\mu_y$ and so on. The metallic branch is obtained in reverse. This approach is appealing because of its similarity to the way actual hysteresis curves in a laboratory are obtained [another approach that can be parallelized is to use bond dimension growth schedules that favor one branch or the other (cf. Appendix B 1)]. However, this method of obtaining the spinodals $\delta\mu_1, \delta\mu_2$ (see Sec. III for definition) will very likely overshoot $\delta\mu_2$ and undershoot $\delta\mu_1$. In other words, we are not at all guaranteed to obtain both branches for the true width of the hysteresis window. Nonetheless, in Fig. 29 we show how our computed $\delta\mu_1, \delta\mu_2$ behave as a function of $\chi$ for some values of $U$, as well as the corresponding

spinodal values of density $n_1 = n(\delta\mu_1)$, $n_2 = n(\delta\mu_2)$. We can see that, although there is a general trend, $\delta\mu_i$ and $n_i$ are not always smooth as functions of $\chi$. Even though we cannot find the true spinodal values, the mere existence of a hysteresis of this type is another indication of the first-order nature of the transition [94]. Furthermore, Fig. 29 suggests that in the $\chi \to \infty$ limit the left spinodal value $n_2$ might vanish for small values of $U$ (although, as we stated, since the behavior is not monotonic we are not doing any extrapolation). If there is a way to stabilize this branch (perhaps with finite temperature or some other parameter), this would imply that the transition could become continuous below a certain value of $U$, which would be very interesting.

Before ending the discussion of the metastable states, we mention another point. As can be seen from Fig. 26, the

<figure><img src="image_30.png" /><figcaption>FIG. 27. (a) The charge gap $\Delta_c(\chi)$ as a function of increasing $\chi$, with a linear fit. (b) The staggered magnetization $m_s(\chi)$ of the insulating state at half-filling as a function of increasing $\chi$, with a linear fit. (c) The jump in occupation at the transition point $\delta n(\delta\mu_c, \chi)$ as a function of increasing $\chi$, with a linear fit. (d) The staggered magnetization $m_s$ of the metallic state at the critical point versus the inverse bond dimension for $U = 3, 6$. The horizontal dashed lines indicate the magnetization of the insulating state across the phase transition, computed at $\chi = 100$.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>FIG. 28. The correlation length scaling with $\chi$ for a range of $U$ and $\langle n_i \rangle = 1.2$. The fits shown are of the form $\xi \sim \chi^\alpha$ with $0.7 < \alpha < 1$.</figcaption></figure>

behavior of the density is quite linear, particularly near $\delta\mu_c$. If we extrapolate this linear behavior, we can define a value $\delta\mu^*$ at which the extrapolation intersects the $\delta\mu$ axis. $\delta\mu^*$ is the lowest possible value of $\delta\mu_2$. Interestingly, as a function of $U$ the value of $\delta\mu^*$ changes from zero to nonzero around a certain value $U^*$. The extrapolation in $\chi$ is difficult to do, due to nonmonotonicity, but we can estimate that $3 < U^* < 4$. This is reminiscent of the behavior found from DMFT in the $z = \infty$ limit, where the metallic spinodal terminates at $\delta\mu_2 = 0$ for $U < U_{c2}$ and at a finite value of $\delta\mu_2$ for $U > U_{c2}$ ($U_{c2}$ being a finite critical value) [13].

Finally, we turn to the question of computing the charge gap. The usual formula is $\Delta_c = E_0(N+1) + E_0(N-1) - 2E_0(N)$, where $E_0(N)$ is the canonical ground-state energy with $N$ particles, and $E_0(N \pm 1)$ are the canonical energies of the metallic ground states with $N \pm 1$ particles, respectively. In this case, it is not immediately clear what is the energy $E_0(N \pm 1)$ since there is no metallic state with density arbitrarily close to one in the grand canonical ensemble (in which our calculation is performed). However, in the canonical ensemble, we can take advantage of phase separation, and insert or remove a particle by creating a mixture of the two phases in the system. Indeed, working at finite volume $V$, at the critical point $\delta\mu_c$ we have an insulator with density $n_I = 1$ and a metal with density $n_M = 1 \pm \delta n$ (the $\pm$ for adding and removing a

particle). The fraction of the system in the insulating phase is given by $\alpha = 1 - \frac{1}{\delta nV}$ (the same for adding or removing a particle because of particle-hole symmetry). Computing the charge gap then gives $\Delta_c = \frac{2(\epsilon_M - \epsilon_I)}{\delta n}$, where $\epsilon_M, \epsilon_I$ are the canonical energy densities of the metallic and insulating states, respectively. The transition point $\delta\mu_c$ is defined such that the grand canonical energy densities of the two states are equal, which gives $\epsilon_M - \delta\mu_c(n_I + \delta n) = \epsilon_I - \delta\mu_c n_I$. This gives $\Delta_c = 2\delta\mu_c$, which is the naive expectation.

## APPENDIX D: OCCUPATION IN $\theta$ SPACE

Here we derive Eq. (10). From Eqs. (8) and (9) we have

$$n(\theta, \theta') = \lim_{L \to \infty} \frac{\pi}{L+1} \sum_{d, d'=0}^{L} \psi_d(\theta) \psi_{d'}(\theta') \langle \tilde{c}_d^\dagger \tilde{c}_d \rangle. \quad (D1)$$

We need to remember that the operator $\tilde{c}_d$ implies a reference to a center site, which we always keep the same. Now we expand $\langle \tilde{c}_{d'}^\dagger \tilde{c}_d \rangle$ in operators $c_i$, where $i$ labels a site on the original Bethe lattice. We show an example in Fig. 14. We use the fact that $\langle c_i^\dagger c_j \rangle = \langle c_j^\dagger c_i \rangle$ is a function of $|i-j|$, which is the length of the path connecting sites $i, j$. We first consider the case when $d, d' > 0, d \neq d'$. Since $d, d'$ label generations emanating from the same center site, the correlation function $\langle \tilde{c}_{d'}^\dagger \tilde{c}_d \rangle$ becomes

$$\langle \tilde{c}_{d'}^\dagger \tilde{c}_d \rangle = \frac{z(z-1)^{\min(d,d')-1}}{z\sqrt{(z-1)^{d+d'-2}}} \sum_{r=0}^{\min(d,d')} A_r \langle c_0^\dagger c_{|d-d'|+2r} \rangle. \quad (D2)$$

The denominator outside of the sum comes from the normalization factor of the symmetric states. The numerator outside of the sum is the number of sites of the inner shell. The coefficients $A_r$ are the number of sites on the outer shell that are connected to a given site on the inner shell by a distance $|d-d'| + 2r$. To find $A_r$ we start with $r=0$. The number of these is $A_0 = (z-1)^{|d-d'|}$. Now we increase $r$ by one, which means we go from the given site on the inner shell up one generation towards the center site, and then back down a different branch. Going back down gives us a factor of $z-2$ to get back to the inner shell, and then another factor of $(z-1)^{|d-d'|}$. For $r=2$, we need to go up two generations towards the center. Coming back down now gives us $(z-2)(z-1)(z-1)^{|d-d'|}$. For a general $0 < r < \min(d, d')$ this gives $A_r = (z-2)(z-1)^{r-1}(z-1)^{|d-d'|}$. The last term in the series is special since moving up one generation bring

<figure><img src="image_32.png" /><figcaption>FIG. 29. The spidonal values $\delta\mu_1, \delta\mu_2$ and the spinodal values of density $n_1, n_2$ for $U = 1, 3$.</figcaption></figure>
us to the center site, so the factor $z-2$ becomes $z-1$ and $A_{\min(d,d')} = (z-1)^{\min(d,d')} (z-1)^{|d-d'|}$. We can rewrite Eq. (D2) as

$$ \begin{aligned} \langle \tilde{c}_d^\dagger \tilde{c}_d \rangle &= \sqrt{(z-1)^{|d-d'|}} \sum_{r=0}^{\min(d,d')} \left(\frac{z-1}{z-2}\right)^{\delta_{r,0} + \delta_{r,\min(d,d')}} (z-2)(z-1)^{r-1} \langle c_0^\dagger c_{|d-d'|+2r} \rangle \\ &= \frac{z-2}{\sqrt{z(z-1)}} \sum_{r=0}^{\min(d,d')} \left(\frac{z-1}{z-2}\right)^{\delta_{r,0} + \delta_{r,\min(d,d')}} \langle \tilde{c}_0^\dagger \tilde{c}_{|d-d'|+2r} \rangle. \end{aligned} \quad (D3) $$

Now we look at the special cases. If $d = d' \neq 0$, everything goes through the same except the conversion from $\langle c_0^\dagger c_{|d-d'|+2r} \rangle$ to $\langle \tilde{c}_0^\dagger \tilde{c}_{|d-d'|+2r} \rangle$, which gives

$$ \langle \tilde{c}_d^\dagger \tilde{c}_d \rangle = \frac{z-2}{\sqrt{z(z-1)}} \sum_{r=0}^{d} \left(\frac{z-1}{z-2}\right)^{\delta_{r,0} + \delta_{r,d}} \left(\frac{z}{z-1}\right)^{\frac{\delta_{r,0}}{2}} \langle \tilde{c}_0^\dagger \tilde{c}_{2r} \rangle. $$

If $d' = 0$, $d \ge 0$ (equivalent to $d' \ge 0$, $d = 0$), the relation should just trivially give back $\langle \tilde{c}_0^\dagger \tilde{c}_d \rangle$. Incorporating these two special cases, Eq. (D3) can be written as

$$ \langle \tilde{c}_d^\dagger \tilde{c}_d \rangle = \frac{z-2}{\sqrt{z(z-1)}} \left( \frac{\sqrt{z}(z-2)}{(z-1)^{3/2}} \right)^{\delta_{\min(d,d'),0}} \sum_{r=0}^{\min(d,d')} \left( \frac{z-1}{z-2} \right)^{\delta_{r,0} + \delta_{r,\min(d,d')}} \sqrt{\frac{z}{z-1}}^{\delta_{d,d'}\delta_{r,0} - \delta_{d,0}\delta_{d',0}} \langle \tilde{c}_0^\dagger \tilde{c}_{|d-d'|+2r} \rangle. \quad (D4) $$

Important to note, this is not a function of $|d-d'|$ only. Therefore, we cannot simplify this further by doing the summation over $d+d'$. The final expression is given by inserting Eq. (D4) into Eq. (D1) and shown in Eq. (10)

[1] P. W. Anderson, The resonating valence bond state in La₂Cuo₄ and superconductivity, Science 235, 1196 (1987).

[2] A.-M. S. Tremblay, B. Kyung, and D. Sénéchal, Pseudogap and high-temperature superconductivity from weak to strong coupling. Towards a quantitative theory (Review Article), Low Temp. Phys. 32, 424 (2006).

[3] H. Alloul, What is the simplest model that captures the basic experimental facts of the physics of underdoped cuprates? C. R. Phys. 15, 519 (2014).

[4] O. Gunnarsson, T. Schäfer, J. P. F. LeBlanc, E. Gull, J. Merino, G. Sangiovanni, G. Rohringer, and A. Toschi, Fluctuation Diagnostics of the Electron Self-Energy: Origin of the Pseudogap Physics, Phys. Rev. Lett. 114, 236402 (2015).

[5] W. Wu, M. Ferrero, A. Georges, and E. Kozik, Controlling Feynman diagrammatic expansions: Physical nature of the pseudogap in the two-dimensional Hubbard model, Phys. Rev. B 96, 041105(R) (2017).

[6] T. Schäfer, N. Wentzell, F. Šimkovic IV, Y.-Y. He, C. Hille, M. Klett, C. J. Eckhardt, B. Arzhang, V. Harkov, F.-M. L. Régent, A. Kirsch, Y. Wang, A. J. Kim, E. Kozik, E. A. Stepanov, A. Kauch, S. Andergassen, P. Hansmann, D. Rohe, Y. M. Vilk, J. P. F. LeBlanc, S. Zhang, A. M. S. Tremblay, M. Ferrero, O. Parcollet, and A. Georges, Tracking the Footprints of Spin Fluctuations: A Multi-Method, Multi-Messenger Study of the Two-Dimensional Hubbard Model, Phys. Rev. X 11, 011058 (2021).

[7] B.-X. Zheng and G. K.-L. Chan, Ground-state phase diagram of the square lattice hubbard model from density matrix em- bedding theory, Phys. Rev. B 93, 035126 (2016).

[8] B.-X. Zheng, C.-M. Chung, P. Corboz, G. Ehlers, M.-P. Qin, R. M. Noack, H. Shi, S. R. White, S. Zhang, and G. K.-L. Chan,

Stripe order in the underdoped region of the two-dimensional hubbard model, Science 358, 1155 (2017).

[9] M. Qin, C.-M. Chung, H. Shi, E. Vitali, C. Hubig, U. Schollwöck, S. R. White, and S. Zhang (Simons Collaboration on the Many-Electron Problem), Absence of Superconductivity in the Pure Two-Dimensional Hubbard Model, Phys. Rev. X 10, 031016 (2020).

[10] E. H. Lieb and F. Y. Wu, Absence of Mott Transition in an Exact Solution of the Short-Range, One-Band Model in One Dimension, Phys. Rev. Lett. 20, 1445 (1968).

[11] A. Georges and G. Kotliar, Hubbard model in infinite dimensions, Phys. Rev. B 45, 6479 (1992).

[12] W. Metzner and D. Vollhardt, Correlated Lattice Fermions in $d = \infty$ Dimensions, Phys. Rev. Lett. 62, 324 (1989).

[13] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68, 13 (1996).

[14] F. H. L. Essler, H. Frahm, F. Gohmann, A. Klumper, and V. E. Korepin, The One-Dimensional Hubbard Model (Cambridge University Press, Cambridge, 2005).

[15] T. Giamarchi, Quantum Physics in One Dimension (Clarendon, Oxford, 2010).

[16] S. R. White, Density Matrix Formulation for Quantum Renormalization Groups, Phys. Rev. Lett. 69, 2863 (1992).

[17] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Ann. Phys. 326, 96 (2011).

[18] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Electronic structure calculations with dynamical mean-field theory, Rev. Mod. Phys. 78, 865 (2006).

[19] Edited by E. Pavarini, E. Koch, D. Vollhardt, and A. Lichtenstein, DMFT at 25: Infinite Dimensions, Modeling and
Simulation, Vol. 4, Autumn School on Correlated Electrons, Jülich (Germany), 15 Sep–19 Sep 2014 (Forschungszentrum Jülich Zentralbibliothek, Verlag, Jülich, 2014).

[20] A. Georges, Strongly correlated electron materials: Dynamical mean-field theory and electronic structure, in *Lectures on the Physics of Highly Correlated Electron Systems VIII*, edited by A. Avella and F. Mancini (American Institute of Physics, Melville, NY, 2004).

[21] M. Rozenberg, Dynamical mean-field theory and the mott transition, in *Many Body Methods for Real Materials*, edited by E. Pavarini E. Koch and S. Zhang, Autumn School on Correlated Electrons, Jülich (Germany), 16 Sep–20 Sep 2019 (Forschungszentrum Jülich GmbH Zentralbibliothek, Verlag, 2019).

[22] R. Preuss, W. Hanke, C. Gröber, and H. G. Evertz, Pseudogaps and Their Interplay with Magnetic Excitations in the Doped 2D Hubbard Model, Phys. Rev. Lett. 79, 1122 (1997).

[23] A. Macridin, M. Jarrell, T. Maier, P. R. C. Kent, and E. D'Azevedo, Pseudogap and Antiferromagnetic Correlations in the Hubbard Model, Phys. Rev. Lett. 97, 036401 (2006).

[24] B. Kyung, S. S. Kancharla, D. Sénéchal, A.-M. S. Tremblay, M. Civelli, and G. Kotliar, Pseudogap induced by short-range spin correlations in a doped Mott insulator, Phys. Rev. B 73, 165114 (2006).

[25] M. S. Scheurer, S. Chatterjee, W. Wu, M. Ferrero, A. Georges, and S. Sachdev, Topological order in the pseudogap metal, Proc. Natl. Acad. Sci. USA 115, E3665 (2018).

[26] R. J. Baxter, *Exactly Solved Models in Statistical Mechanics* (Academic, New York, 1982).

[27] H. A. Bethe, Statistical theory of superlattices, Proc. R. Soc. London A 150, 552 (1935).

[28] P. D. Gujrati, Bethe or Bethe-Like Lattice Calculations are More Reliable than Conventional Mean-Field Calculations, Phys. Rev. Lett. 74, 809 (1995).

[29] C. R. Laumann, S. A. Parameswaran, and S. L. Sondhi, Absence of goldstone bosons on the bethe lattice, Phys. Rev. B 80, 144415 (2009).

[30] D. Nagaj, E. Farhi, J. Goldstone, P. Shor, and I. Sylvester, Quantum transverse-field ising model on an infinite tree from matrix product states, Phys. Rev. B 77, 214431 (2008).

[31] A. Nagy, Simulating quantum systems on the bethe lattice by translationally invariant infinite-tree tensor network, Ann. Phys. 327, 542 (2012).

[32] W. Li, J. von Delft, and T. Xiang, Efficient simulation of infinite tree tensor network states on the bethe lattice, Phys. Rev. B 86, 195137 (2012).

[33] S. Depenbrock and F. Pollmann, Phase diagram of the isotropic spin-3/2 model on the z = 3 bethe lattice, Phys. Rev. B 88, 035138 (2013).

[34] T. Liu, S.-J. Ran, W. Li, X. Yan, Y. Zhao, and G. Su, Featureless quantum spin liquid, 1/3-magnetization plateau state, and exotic thermodynamic properties of the spin-1/2 frustrated Heisenberg antiferromagnet on an infinite Husimi lattice, Phys. Rev. B 89, 054426 (2014).

[35] G.-H. Liu, W. Li, X.-Y. Deng, G. Su, and G.-S. Tian, Existence of the 1/2 magnetization plateau in the s = 1 ising model with single-ion anisotropy on z = 3 bethe lattice, Solid State Commun. 180, 1 (2014).

[36] M. Mézard and G. Parisi, The bethe lattice spin glass revisited, Eur. Phys. J. B 20, 217 (2001).

[37] D. J. Thouless, Spin-Glass on a Bethe Lattice, Phys. Rev. Lett. 56, 1082 (1986).

[38] T. K. Kopec and K. D. Usadel, Short-range ±j interaction ising spin glass in a transverse field on a bethe lattice: a quantum-spherical approach, Phys. Status Solidi B 243, 502 (2006).

[39] C. Laumann, A. Scardicchio, and S. L. Sondhi, Cavity method for quantum spin glasses on the bethe lattice, Phys. Rev. B 78, 134424 (2008).

[40] G. Semerjian, M. Tarzia, and F. Zamponi, Exact solution of the bose-hubbard model on the bethe lattice, Phys. Rev. B 80, 014524 (2009).

[41] R. Abou-Chacra, D. J. Thouless, and P. W. Anderson, J. Phys. C: Solid State Phys. 6, 1734 (1973).

[42] A. D. Mirlin and Y. V. Fyodorov, Localization transition in the anderson model on the bethe lattice: Spontaneous symmetry breaking and correlation functions, Nucl. Phys. B 366, 507 (1991).

[43] B. L. Altshuler, Y. Gefen, A. Kamenev, and L. S. Levitov, Quasiparticle Lifetime in a Finite System: A Nonperturbative Approach, Phys. Rev. Lett. 78, 2803 (1997).

[44] S. Savitz, C. Peng, and G. Refael, Anderson localization on the bethe lattice using cages and the wegner flow, Phys. Rev. B 100, 094201 (2019).

[45] M. Dupont, N. Laflorencie, and G. Lemarie, Phys. Rev. B 102, 174205 (2020).

[46] M.-B. Lepetit, M. Cousy, and G. M. Pastor, Density-matrix renormalization study of the Hubbard model [4] on a Bethe lattice, Eur. Phys. J. B 13, 421 (2000).

[47] V. J. Emery, S. A. Kivelson, and H. Q. Lin, Phase Separation in the t-j Model, Phys. Rev. Lett. 64, 475 (1990).

[48] B. Ponsioen, S. S. Chung, and P. Corboz, Period 4 stripe in the extended two-dimensional Hubbard model, Phys. Rev. B 100, 195141 (2019).

[49] S. Sorella, The phase diagram of the Hubbard model by variational auxiliary field quantum Monte Carlo, arXiv:2101.07045.

[50] M. Imada, A. Fujimori, and Y. Tokura, Metal-insulator transitions, Rev. Mod. Phys. 70, 1039 (1998).

[51] M. Imada, Metal-insulator transition of correlated systems and origin of unusual metal, J. Phys. Soc. Jpn. 64, 2954 (1995).

[52] P. Corboz, R. Orús, B. Bauer, and G. Vidal, Simulation of strongly correlated fermions in two spatial dimensions with fermionic projected entangled-pair states, Phys. Rev. B 81, 165104 (2010).

[53] P. Corboz, Improved energy extrapolation with infinite projected entangled-pair states applied to the two-dimensional hubbard model, Phys. Rev. B 93, 045116 (2016).

[54] Q. Mortier, N. Schuch, F. Verstraete, and J. Haegeman, Resolving Fermi surfaces with tensor networks, arXiv:2008.11176.

[55] V. Zauner-Stauber, L. Vanderstraeten, M. T. Fishman, F. Verstraete, and J. Haegeman, Variational optimization algorithms for uniform matrix product states, Phys. Rev. B 97, 045145 (2018).

[56] Y.-Y. Shi, L.-M. Duan, and G. Vidal, Classical simulation of quantum many-body systems with a tree tensor network, Phys. Rev. A 74, 022320 (2006).

[57] V. Murg, F. Verstraete, Ö. Legeza, and R. M. Noack, Simulating strongly correlated quantum systems with tree tensor networks, Phys. Rev. B 82, 205105 (2010).
[58] N. Nakatani and G. K.-L. Chan, Efficient tree tensor network states (ttns) for quantum chemistry: Generalizations of the density matrix renormalization group algorithm, *J. Chem. Phys.* **138**, 134113 (2013).

[59] R. Orús, Advances on tensor network theory: symmetries, fermions, entanglement, and holography, *Eur. Phys. J. B* **87**, 280 (2014).

[60] H. J. Changlani, S. Ghosh, C. L. Henley, and A. M. Läuchli, Heisenberg antiferromagnet on cayley trees: Low-energy spectrum and even/odd site imbalance, *Phys. Rev. B* **87**, 085107 (2013).

[61] M. Ostilli, Cayley trees and bethe lattices: A concise analysis for mathematicians and physicists, *Phys. A (Amsterdam)* **391**, 3417 (2012).

[62] G. Vidal, Classical Simulation of Infinite-Size Quantum Lattice Systems in One Spatial Dimension, *Phys. Rev. Lett.* **98**, 070201 (2007).

[63] R. J. Bursill, T. Xiang, and G. A. Gehring, The density matrix renormalization group for a quantum spin chain at non-zero temperature, *J. Phys.: Condens. Matter* **8**, L583 (1996).

[64] X. Wang and T. Xiang, Transfer-matrix density-matrix renormalization-group theory for thermodynamics of one-dimensional quantum systems, *Phys. Rev. B* **56**, 5061 (1997).

[65] N. Shibata, Thermodynamics of the anisotropic Heisenberg chain calculated by the density matrix renormalization group method, *J. Phys. Soc. Jpn.* **66**, 2221 (1997).

[66] I. P. McCulloch, Infinite size density matrix renormalization group, revisited, arXiv:0804.2509.

[67] H. N. Phien, J. A. Bengua, H. D. Tuan, P. Corboz, and R. Orús, Infinite projected entangled pair states algorithm improved: Fast full update and gauge fixing, *Phys. Rev. B* **92**, 035142 (2015).

[68] P. Corboz, Variational optimization with infinite projected entangled-pair states, *Phys. Rev. B* **94**, 035133 (2016).

[69] L. Vanderstraeten, J. Haegeman, P. Corboz, and F. Verstraete, Gradient methods for variational optimization of projected entangled-pair states, *Phys. Rev. B* **94**, 155123 (2016).

[70] M. T. Fishman, L. Vanderstraeten, V. Zauner-Stauber, J. Haegeman, and F. Verstraete, Faster methods for contracting infinite two-dimensional tensor networks, *Phys. Rev. B* **98**, 235148 (2018).

[71] H.-J. Liao, J.-G. Liu, L. Wang, and T. Xiang, Differentiable Programming Tensor Networks, *Phys. Rev. X* **9**, 031041 (2019).

[72] G. R. Stewart, Non-Fermi-liquid behavior in *d*- and *f*-electron metals, *Rev. Mod. Phys.* **73**, 797 (2001).

[73] S.-S. Lee, Recent developments in non-Fermi liquid theory, *Annu. Rev. Condens. Matter Phys.* **9**, 227 (2018).

[74] P. Chandra and B. Doucot, Spin liquids on the husimi cactus, *J. Phys. A: Math. Gen.* **27**, 1541 (1994).

[75] A. Camjayi, R. Chitra, and M. J. Rozenberg, Electronic state of a doped mott-hubbard insulator at finite temperatures studied using the dynamical mean-field theory, *Phys. Rev. B* **73**, 041103(R) (2006).

[76] X. Wang, E. Gull, L. de' Medici, M. Capone, and A. J. Millis, Antiferromagnetism and the gap of a mott insulator: Results from analytic continuation of the self-energy, *Phys. Rev. B* **80**, 045101 (2009).

[77] L. Fratino, M. Charlebois, P. Sémon, G. Sordi, and A.-M. S. Tremblay, Effects of interaction strength, doping, and frustration on the antiferromagnetic phase of the two-dimensional hubbard model, *Phys. Rev. B* **96**, 241109(R) (2017).

[78] G. Moeller, Q. Si, G. Kotliar, M. Rozenberg, and D. S. Fisher, Critical Behavior Near the Mott Transition in the Hubbard Model, *Phys. Rev. Lett.* **74**, 2082 (1995).

[79] G. Kotliar, S. Murthy, and M. J. Rozenberg, Compressibility Divergence and the Finite Temperature Mott Transition, *Phys. Rev. Lett.* **89**, 046401 (2002).

[80] P. Werner and A. J. Millis, Doping-driven mott transition in the one-band hubbard model, *Phys. Rev. B* **75**, 085108 (2007).

[81] G. D. Mahan, Energy bands of the Bethe lattice, *Phys. Rev. B* **63**, 155110 (2001).

[82] We thank Martin Eckstein for pointing out these considerations.

[83] G. Kotliar and A. E. Ruckenstein, New Functional Integral Approach to Strongly Correlated Fermi Systems: The Gutzwiller Approximation as a Saddle Point, *Phys. Rev. Lett.* **57**, 1362 (1986).

[84] V. Zauner-Stauber, L. Vanderstraeten, J. Haegeman, I. P. McCulloch, and F. Verstraete, Topological nature of spinons and holons: Elementary excitations from matrix product states with conserved symmetries, *Phys. Rev. B* **97**, 235155 (2018).

[85] R. Nandkishore, M. A. Metlitski, and T. Senthil, Orthogonal metals: The simplest non-fermi liquids, *Phys. Rev. B* **86**, 045128 (2012).

[86] M. J. Rozenberg, G. Kotliar, and X. Y. Zhang, Mott-hubbard transition in infinite dimensions. II, *Phys. Rev. B* **49**, 10181 (1994).

[87] H. J. Liao, Z. Y. Xie, J. Chen, X. J. Han, H. D. Xie, B. Normand, and T. Xiang, Heisenberg antiferromagnet on the husimi lattice, *Phys. Rev. B* **93**, 075154 (2016).

[88] M. Udagawa and R. Moessner, Spectrum of Itinerant Fractional Excitations in Quantum Spin Ice, *Phys. Rev. Lett.* **122**, 117201 (2019).

[89] M. Zwolak and G. Vidal, Mixed-State Dynamics in One-Dimensional Quantum Lattice Systems: A Time-Dependent Superoperator Renormalization Algorithm, *Phys. Rev. Lett.* **93**, 207205 (2004).

[90] F. Verstraete, J. J. García-Ripoll, and J. I. Cirac, Matrix Product Density Operators: Simulation of Finite-Temperature and Dissipative Systems, *Phys. Rev. Lett.* **93**, 207204 (2004).

[91] A. E. Feiguin and S. R. White, Finite-temperature density matrix renormalization using an enlarged Hilbert space, *Phys. Rev. B* **72**, 220401(R) (2005).

[92] D.-W. Qu, W. Li, and T. Xiang, Thermal tensor network simulations of the heisenberg model on the Bethe lattice, *Phys. Rev. B* **100**, 125121 (2019).

[93] M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor software library for tensor network calculations, arXiv:2007.14822.

[94] S. Yan, D. A. Huse, and S. R. White, Spin-liquid ground state of the s = 1/2 kagome Heisenberg antiferromagnet, *Science* **332**, 1173 (2011).