Thermal Tensor Network Simulations of the Heisenberg Model on the Bethe Lattice

TL;DR

This paper extends the tensor network method to finite temperatures to study the phase diagram and magnetic properties of the Heisenberg XXZ model on the Bethe lattice, revealing various magnetic phases and excitation spectra.

Contribution

The authors develop a finite-temperature tensor network approach for the Bethe lattice, enabling detailed phase diagram and thermodynamic property calculations for the Heisenberg model.

Findings

Identified three magnetic phases: ferromagnetic, XY, and antiferromagnetic.

Discovered continuous and first-order phase transitions at specific coupling points.

Found that magnetic excitations are gapped, with Goldstone-like modes in ordered phases.

Abstract

We have extended the canonical tree tensor network (TTN) method, which was initially introduced to simulate the zero-temperature properties of quantum lattice models on the Bethe lattice, to finite temperature simulations. By representing the thermal density matrix with a canonicalized tree tensor product operator, we optimize the TTN and accurately evaluate the thermodynamic quantities, including the internal energy, specific heat, and the spontaneous magnetization, etc, at various temperatures. By varying the anisotropic coupling constant $\Delta$, we obtain the phase diagram of the spin-1/2 Heisenberg XXZ model on the Bethe lattice, where three kinds of magnetic ordered phases, namely the ferromagnetic, XY and antiferromagnetic ordered phases, are found in low temperatures and separated from the high-$T$ paramagnetic phase by a continuous thermal phase transition at $T_c$. The XY phase is separated from the other two phases by two first-order phase transition lines at the symmetric coupling points $ \Delta=\pm 1$. We have also carried out a linear spin wave calculation on the Bethe lattice, showing that the low-energy magnetic excitations are always gapped, and find the obtained magnon gaps in very good agreement with those estimated from the TTN simulations. Despite the gapped excitation spectrum, Goldstone-like transverse fluctuation modes, as a manifestation of spontaneous continuous symmetry breaking, are observed in the ordered magnetic phases with $|\Delta|\le 1$. One remarkable feature there is that the prominent transverse correlation length reaches $\xi_c=1/\ln{(z-1)}$ for $T\leq T_c$, the maximal value allowed on a $z$-coordinated Bethe lattice.

Tables

Table 1. Table 1: Magnon gaps of the XXZ model with different anisotropies Δ Δ \Delta and external fields h ℎ h . The numerical and LSWT results are compared, where the numbers in the parentheses in the TTN data are numerical fitting errors.

Model	($\Delta$, $h$)	${\Delta }_m$ (TTN, fitted)	${\Delta }_m$ (LSWT, $z=3$)
AF	(1, 0)	0.56(3)	$0.5$
XY	(0, 0)	0.41(1)	$0.3587$
FM	(-1, 0)	0.085(3)	$0.0858$
	(-1, 0.1)	0.186(1)	$0.1858$
	(-1, 0.2)	0.286(2)	$0.2858$

Equations

$$H=\sum _{⟨i,j⟩}(S_i^x{S}_j^x+S_i^y{S}_j^y+\Delta S_i^z{S}_j^z)-h\sum _i{S}_i^z,$$

$$\sum _{x,y=1}^D\sum _{\sigma ,{\sigma }^{\prime }=1}^d{\lambda }_x^2{\lambda }_y^2(T_{\alpha }^{\ast }{)}_{xyz}^{\sigma {\sigma }^{\prime }}(T_{\alpha }{)}_{xy{z}^{\prime }}^{\sigma {\sigma }^{\prime }}={\mathbb{I}}_{zz\prime },$$

$$\rho (\beta )$$

$$H_{\zeta }$$

$$P_{\zeta }=\exp \left(-\tau H_{\zeta }\right)=\prod _i\exp \left(-\tau h_{i,i+\zeta }\right).$$

$$\rho (\beta )=\rho \left(\beta /2\right)\rho \left(\beta /2\right)={\rho }^{†}\left(\beta /2\right)\rho \left(\beta /2\right).$$

$$Z(\beta )=\mathrm{Tr}\rho (\beta )=⟨\rho (\beta /2)|\rho (\beta /2)⟩.$$

$$M_{x{x}^{\prime };y{y}^{\prime };z{z}^{\prime }}^{\alpha }=\sum _{\begin{array}{c}{\sigma }_1,{\sigma }_2\end{array}}(T_{\alpha }{)}_{xyz}^{{\sigma }_1{\sigma }_2}(T_{\alpha }^{\ast }{)}_{x^{\prime }{y}^{\prime }{z}^{\prime }}^{{\sigma }_1{\sigma }_2},$$

$$\sum _{\begin{array}{c}x{x}^{\prime }y{y}^{\prime }\end{array}}(X_A{)}_{x{x}^{\prime }}(Y_A{)}_{y{y}^{\prime }}{\Lambda }_{x{x}^{\prime }y{y}^{\prime }}{M}_{x{x}^{\prime },y{y}^{\prime },z{z}^{\prime }}^B={\eta }_B^Z(Z_B{)}_{z{z}^{\prime }},$$

$$({\rho }_A{)}_{{\sigma }_1{\sigma }_2}$$

$$⟨\hat{O}{⟩}_{\beta }=\frac{\mathrm{Tr}({\rho }_A\hat{O})}{\mathrm{Tr}{\rho }_A},$$

$$S_E=-\mathrm{Tr}\left({\lambda }_{\zeta }^2\ln {\lambda }_{\zeta }^2\right),$$

$$m=\sqrt{m_x^2+m_z^2},$$

$$m^2\propto \frac{T_c-T}{T_c}.$$

$$c_V\sim e^{-\beta {\Delta }_m}.$$

$$u(\beta )-u_0\simeq g(\beta )e^{-\beta {\Delta }_m}$$

$$\frac{d\ln [u(\beta )-u_0]}{d\beta }=\frac{1}{g}\frac{dg}{d\beta }-{\Delta }_m,$$

$${\xi }_i=2{\ln }^{-1}\left|\frac{{\eta }_0}{{\eta }_i}\right|,$$

$$C^{\zeta }(l)=⟨S_i^{\zeta }{S}_{i+l}^{\zeta }⟩-⟨S_i^{\zeta }⟩⟨S_{i+l}^{\zeta }⟩,\text{ (}\zeta =x,y,z\text{)}$$

$$H=\sum _{⟨i,j⟩}(-S_i^x{S}_j^x-S_i^y{S}_j^y+\Delta S_i^z{S}_j^z)-h\sum _i{S}_i^z.$$

$$H=(-zS\Delta +h)\sum _{l,R_l}{n}_{R_l}-S\sum _{R_l}\sum _{r_{l+1}}(a_{R_l}^{†}{a}_{R_l{r}_{l+1}}+h.c.),$$

$$a_{Q_l}=\sum _{R_l}\left(\prod _{\lambda =1}^l{U}_{q_{\lambda },r_{\lambda }}\right)a_{R_l},$$

$$U_{q_{\lambda },r_{\lambda }}=\frac{1}{\sqrt{\tilde{z}}}{e}^{2\pi \mathrm{i}\cdot q_{\lambda }{r}_{\lambda }/\tilde{z}}.$$

$$\sum _{R_l,r_{l+1}}{a}_{R_l}^{†}{a}_{R_l{r}_{l+1}}=\sqrt{\tilde{z}}\sum _{Q_l}{a}_{Q_l}^{†}{a}_{Q_l,0},$$

$$\sum _{R_l}{a}_{R_l}^{†}{a}_{R_l}=\sum _{Q_l}{a}_{Q_l}^{†}{a}_{Q_l}.$$

$$H_{\text{FM}}=\sum _{i=0}^{\infty }(-zS\Delta +h)b_i^{†}{b}_i-S\sqrt{\tilde{z}}(b_i^{†}{b}_{i+1}+h.c.),$$

$$\epsilon (\kappa )=-zS\Delta +h-2S\sqrt{z-1}\cos \kappa ,$$

$${\Delta }_m=-zS\Delta +h-2S\sqrt{z-1}.$$

$$H=zS\Delta \sum _{l,R_l}{n}_{R_l}+S\sum _{l,R_l}\sum _{r_{l+1}}(a_{R_l}{a}_{R_l{r}_{l+1}}+a_{R_l}^{†}{a}_{R_l{r}_{l+1}}^{†}).$$

$$\sum _{R_l,r_{l+1}}{a}_{R_l}{a}_{R_l{r}_{l+1}}=\sum _{Q_l}{a}_{Q_l}{a}_{Q_l,0}.$$

$$H_{\mathrm{AF}}=\sum _{i=0}^{\infty }S[z\Delta b_i^{†}{b}_i+\sqrt{\tilde{z}}(b_i{b}_{i+1}+b_i^{†}{b}_{i+1}^{†})].$$

Full text

I Introduction

The Bethe approximation, a renowned cluster mean-field approach proposed by Bethe in 1935 Bethe (1935), has played an important role in the studies of cooperative phenomena and phase transitions of classical statistical models Pathria and Beale (2011). This approximation can be rigorously formulated on an ideal lattice of infinite Hausdorff dimension, i.e, the Bethe lattice shown in Fig. 1(a) Ostilli (2012). The Bethe lattice is also intimately connected with the dynamical mean-field theory Eckstein et al. (2005); Petrova and Moessner (2016); Semerjian et al. (2009).

A Bethe lattice is a loop-free graph where each site is connected to $z$ neighbours, i.e., $z$ is the coordination number. Figure 1(a), as an example, shows the structure of a $z=3$ Bethe lattice, whose lattice sites are all equivalent and there exists no boundary on this infinite lattice. This Bethe lattice resembles a two-dimensional honeycomb lattice locally [as emphasized in Fig. 1(b)], but it does not contain any closed loops, e.g., hexagons. A Bethe lattice becomes a Cayley tree if the lattice size is finite, where the sites are arranged in shells around a root site. In stark contrast to the Bethe lattice, there exist boundary sites, which are as many as the bulk sites, on the $z=3$ Cayley tree.

Tensor networks provide efficient and accurate representations of quantum manybody states both at zero and finite temperatures. The simple update Jiang et al. (2008), together with many other optimization schemes Verstraete et al. (2008); Orús (2014), has been widely adopted in the tensor-network simulations of quantum lattice models Li et al. (2012, 2010); Liu et al. (2014, 2015, 2016). It has also been shown that the simple update Jiang et al. (2008), or more rigorously the canonical tree tensor network (TTN) method, is numerically exact on the Bethe lattice Li et al. (2012). In addition, the density matrix renormalization group (DMRG) has also been employed to study the magnetic orders and other physical properties of the Heisenberg model on the Cayley tree Otsuka (1996); Friedman (1997); Kumar et al. (2012); Changlani et al. (2013a, b). Besides, the Bethe approximation has also been used in investigating the Fermi-Hubbard systems, where the single-particle Green’s function as well as the density of states are calculated Brinkman and Rice (1970); Kittler and Falicov (1976); Brouers and Marconi (1982). However, most of these studies are restricted in the $T=0$ properties, and there have not much investigations on the thermodynamic properties of quantum lattice models on the Bethe lattice.

In this paper, we extend the TTN approach from zero to finite temperatures and show that it provides an efficient and accurate method to simulate the thermodynamic properties of the Bethe-lattice quantum lattice models. Through the calculations of magnetic order parameters, we obtain the finite-temperature phase diagram of the anisotropic Heisenberg model. Similar as in the two-dimensional honeycomb lattice, three different magnetic ordered phases, i.e., the ferromagnetic (FM), antiferromagnetic (AF), and planar XY phases, are found on the Bethe lattice. The planar XY phase is separated from the FM and AF phases by two first-order phase transition lines at $\Delta=-1$ and $1$, respectively, again resembling the two-dimensional case Al Hajj et al. (2004).

The correlation length, as shown in Ref. Li et al. (2012), is finite on a Bethe lattice. Here we show that a thermal phase transition can nevertheless happen when the correlation length reaches a “critical” value $\xi_{c}=1/\ln{(z-1)}$ on the Bethe lattice. As revealed by the temperature dependence of thermodynamic quantities in low temperatures, the low-energy excitations of the XXZ model are always gapped. We propose a linear spin wave theory (LSWT) on the Bethe lattice, which gives insight into the low-energy excitations of the system and provides good estimates of the magnon gaps.

The paper is organized as follows. In Sec. II, we briefly introduce the Heisenberg XXZ model and the canonical TTN method on the Bethe lattice. The results obtained with this method are presented in Sec. III. In Sec. IV, we present a spin wave analysis of the model based on the $q$-representation. Finally, in Sec. V, we summarize and discuss about further applications of the method introduced in this work.

II Interacting Spin Model and Tensor Network Method

II.1 The Heisenberg XXZ model

The Hamiltonian of Heisenberg XXZ model reads

[TABLE]

where $\Delta$ is the anisotropic coupling constant, and $\langle i,j\rangle$ denotes a pair of nearest-neighboring sites. $h$ represents an external magnetic field, which is set as zero if not mentioned explicitly. Note the model with $\Delta=-1$ is equivalent to the isotropic FM Heisenberg model upon a $\pi$-rotation around the $z$-axis on one of the two sublattices, labelled by green and orange colors in Fig. 1(a), respectively. In the discussions below, we concentrate on the XXZ model defined on a Bethe lattice of coordination number $z=3$.

II.2 Canonical TTN representation of the density matrix

As shown in Fig. 1(b), the density matrix on the $z=3$ Bethe lattice can be expressed as a TTN. A local tensor $T$ defined at a node contains two physical bonds, labelled by $\sigma$ and $\sigma^{\prime}$, representing respectively the initial and final states of the density matrix. There are three geometric bonds, labelled by $\zeta$ ($=x,y,z$), through which the network is connected. Besides this, there is also a vector or diagonal matrices $\lambda_{\zeta}$ defined on each internal bond $\zeta$, whose square (i.e., $\lambda_{\zeta}^{2}$) represents the entanglement spectra between the two blocks that are separated by this bond.

Given a density matrix, its TTN representation is not uniquely determined. By inserting a pair of invertible matrices, $U$ and its inverse $U^{-1}$, on each internal bond, we do not alter the density matrix. Hence, the TTN representations contain huge gauge redundancy. Nevertheless, the gauge degree of freedom can be fixed by converting a TTN into a canonical form, in which the local tensors $T$’s and diagonal bond matrices $\lambda$’s satisfy a set of canonical equations. Taking the $z$ bond as an example, the canonical equation is

[TABLE]

where $\alpha=A$ or $B$, representing the two sublattices of the $z=3$ Bethe lattice, $d=2$ is the dimension of the local basis states of spin-1/2, and $\mathbb{I}$ is a $D\times D$ identity matrix, with $D$ the dimension of geometric bond. The canonical equations along the $x$ and $y$ bonds can be obtained from Eq. (2) through a cyclic permutation of $x$, $y$, and $z$. This kind of canonical form has already been used in one-dimensional tensor network algorithms, including the density matrix renormalization group (DMRG) White (1992), time-evolution block decimation Vidal (2007); Orús and Vidal (2008), and linearized tensor renormalization group Li et al. (2011); Dong et al. (2017), etc.

Given an arbitrary TTN representation which generically does not satisfy these canonical equations, we can nevertheless gauge the TTN into the canonical form through a so-called canonicalization procedure elaborated in App. A.2.

II.3 Imaginary-time evolution

The tensor network representation of thermal density matrix $\rho(\beta)$ Li et al. (2011); Ran et al. (2012); Chen et al. (2017); Dong et al. (2017); Chen et al. (2018); Czarnik et al. (2012); Czarnik and Dziarmaga (2015); Czarnik et al. (2016, 2017); Kshetrimayum et al. (2019); Czarnik et al. (2019); Czarnik and Corboz can be determined by taking an imaginary-time evolution, starting from an infinitely high temperature at which $\rho(\beta)$ is represented by the identity operator. This is achieved by taking a Trotter-Suzuki decomposition for the density matrix

[TABLE]

where $\beta\equiv 1/T$ is the inverse temperature and $\tau=\beta/N$. $h_{i,i+\eta}$ is the interacting Hamiltonian between (nearest-neighboring) sites $i$ and $i+\zeta$ along $\zeta$(=$x,y,z$) directions. In practical calculations, the Trotter step $\tau$ is set as a small value, e.g., $\tau=0.01$, to control the Trotter error. For the model we study, all the local terms in each $H_{\eta}$ commute with each other, i.e., $[h_{i,i+\zeta},h_{i^{\prime},i^{\prime}+\zeta}]=0$, thus we can further decompose $\exp\left(-\tau H_{\zeta}\right)$ into a product of local evolution gates, i.e.,

[TABLE]

A detailed introduction to the update scheme of local tensors in the imaginary-time evolution is given in App. A.1.

We dub the simple update equipped with the canonicalization procedure (App. A.2) as canonical TTN scheme on the Bethe lattice. The canonical TTN approach adopted in the present work indeed improves the accuracy and stability of the calculations, especially for the case with a thermal phase transition, as benchmarked in App. A.3.

After obtaining the tensor network representation, to evaluate thermodynamic quantities, we notice that the density matrix opertor $\rho(\beta)$ is a product of two half density matrix operators $\rho(\beta/2)$

[TABLE]

This TTN operator $\rho(\beta/2)$ can be also regarded as a “supervector” $|\rho(\beta/2)\rangle$, defined in the (enlarged) product space of the “ket” and “bra” spaces. In this supervector representation, the partition function then becomes an inner product of $|\rho(\beta/2)\rangle$ and its vector conjugate, namely

[TABLE]

Thus the partition function can be obtained by contracting a bilayer TTN Dong et al. (2017).

On the Bethe lattice, contracting the above bilayer supervectors is equivalent to solving the dominant eigen-problem of the transfer tensor $M^{\alpha}$ defined by

[TABLE]

whose schematic representation is shown Fig. 2(a). Supposing that $X_{\alpha}$, $Y_{\alpha}$ and $Z_{\alpha}$ are the three dominant eigenvectors to find, $M^{\alpha}$ ($\alpha=A$ or $B$) simply transfers two of these eigenvectors in one of the sublattice to an eigenvector in the other sublattice, see Figs. 2(b-d). For example, $M^{B}$ transfers $X_{A},Y_{A}$ to $Z_{B}$ through the equation [Fig. 2(b)],

[TABLE]

where $\Lambda_{xx^{\prime}yy^{\prime}}=\lambda_{x}\lambda_{x^{\prime}}\lambda_{y}\lambda_{y^{\prime}}$. This defines a generalized eigenvalue problem associated with the transfer tensor $M^{\alpha}$ and can be solved iteratively, i.e., we start from six random vectors and perform the contractions until all vectors converge.

In stark contrast to the conventional linear eigenvalue problem, there is a gauge flexibility in the definition of eigenvalues $\eta_{\alpha}$ here. Once the normalization condition of eigenvectors $X_{\alpha}$, $Y_{\alpha}$ and $Z_{\alpha}$ is varied, the corresponding eigenvalues $\eta_{\alpha}^{X,Y,Z}$ also change, i.e, they are not uniquely defined. In the following, we fix the gauge by normalizing all the dominating eigenvectors to 1, e.g., $\lVert X_{\alpha}\rVert=\sqrt{\mathrm{Tr}(X_{\alpha}^{\dagger}X_{\alpha})}=1$. This is a convenient choice, and note that once the density matrix $\rho\left(\beta/2\right)$ is in the canonical form, c.f., Eq. (2), the above iterative contraction procedures can be skipped, since each dominant eigenvector is just identity matrices if represented as a $D\times D$ matrix.

Given the eigenvectors $X_{\alpha}$, $Y_{\alpha}$ and $Z_{\alpha}$, we can evaluate the thermal expectation values of local operators, such as the local magnetization. For example, to evalueate the expectation value of an operator $\hat{O}$ on the $A$-sublattice, we first construct the single-site reduced density matrix

[TABLE]

The expectation value is then given by

[TABLE]

which no longer depends on how the generalized dominant eigenvectors are normalized, since $\rho_{\alpha}$ appears in both the numerator and the denominator. Similarly, we can evaluate the bond energies from the two-site reduced density matrix, etc.

From the canonical TTN, we can also calculate the bipartite-entanglement entropy $S_{E}$ using the normalized entanglement spectrum $\lambda_{\zeta}$, as

[TABLE]

which reflects both the quantum entanglement and classical correlations in a thermal equilibrium state. For gapless systems, the entanglement $S_{E}$ might exhibit a universal logarithmic scaling as a function of temperatures at low $T$, e.g., in the one-dimensional quantum critical points and two-dimensional Heisenberg models Prosen and Pižorn (2007); Barthel ; Dubail ; Chen et al. (2018, 2019).

Therefore, it is of great interest to explore the scaling behaviors of $S_{E}$, particularly near the phase transition temperatures, for the XXZ model on the Bethe lattice. Note, $S_{E}$ provides a quantitive measure of the bond dimension that is needed for an accurate representation of the thermal density matrix, especially in low temperatures. Generically, the bond dimension $D$ scales exponentially with $S_{E}$, which would be saturated at low temperatures in a gapped system. On the contrary, in a gapless system $D$ would scale algebraically with inverse temperature $\beta$, given that $S_{E}\sim\ln\beta$.

III Numerical results

Here we present the thermodynamic results calculated with the canonical TTN method. In our calculations, up to $D=80$ states are retained in the geometric bonds of local tensors to ensure that the results are converged down to $T=0.05$.

To benchmark the method, we have evaluated the thermodynamic quantities with the canonical TTN in the classical limit $\Delta=\infty$ [see Eq. (1)], at which the XXZ model is reduced to the exactly soluble Ising model. In this limit, the density matrix has an exact TTN representation with a bond dimension $D=2$. We find that our numerical results agree excellently with the exact values. One can refer to App. B for details. Below, we show our TTN results of the quantum XXZ model at finite $T$.

III.1 Phase diagram

It has been well-established that on a two-dimensional bipartite lattice, say, the honeycomb or square lattice, the XXZ model with $-1\leq\Delta\leq 1$ breaks the continuous symmetry and possesses a long-range order in the ground state. However, at finite temperatures, the continuous symmetry is restored and the long-range magnetic order is melted by the low-lying excitations according to the Mermin-Wagner theorem Mermin and Wagner (1966).

Figure 3 shows the temperature dependence of the magnetic order parameter

[TABLE]

where $m_{x}=\langle S_{i}^{x}\rangle_{\beta}$ and $m_{z}=\langle S_{i}^{z}\rangle_{\beta}$ are the components within and perpendicular to the XY plane ($\langle S_{i}^{y}\rangle_{\beta}=0$ by default), respectively. The absolute values of $m$ are the same on the $A$ and $B$ sublattices, and there is no spontaneous symmetry breaking, i.e., $m=0$, in high temperatures. However, $m$ becomes finite when the temperature drops below a critical value $T_{c}$. Particularly, we find that $m^{2}$ varies linearly with temperature just below $T_{c}$, i.e.,

[TABLE]

We have checked three cases shown in Fig. 3(a) with $\Delta=0$, $1$ and $1.2$, which all fall into this scalings in the vicinity of $T_{c}$, thus indicating that the transition is mean-field-like.

Figure 4(a) shows the temperature dependence of the specific heat at several different $\Delta$ values. A $\lambda$-jump is observed at the critical point $T_{c}$, which confirms that the transition from the paramagnetic to magnetic-ordered phase is continuous (in a mean-field-like fashion). On the other hand, as shown in Fig. 4(b), the entanglement entropy curve vs. temperature exhibits a cusp, rather than a diverging peak, at $T_{c}$. The absence of divergent $S_{E}$ at $T_{c}$ is a natural consequence of the finite correlation length $\xi$ in the system (see discussions on $\xi$ in Sec. III.2 below). It allows us to perform accurate thermal simulations down to low temperatures, $T\leq T_{c}$, by retaining a finite number $D$ of bond states.

Figure 5 shows the $T$-$\Delta$ phase diagram of the Heisenberg XXZ model on the Bethe lattice, where three magnetic ordered phases are observed in low temperatures, by varying the anisotropic parameter $\Delta$. As depicted in Fig. 3(b), for $\Delta>1$ and $\Delta<-1$, the low-temperature states are AF and FM ordered, respectively. The system spontaneously breaks the $Z_{2}$ symmetry, as characterized by $m_{z}\neq 0$ and $m_{x}=0$. When $|\Delta|<1$, the planar U(1) symmetry is broken, with $m_{x}\neq 0$ and $m_{z}=0$ in low temperatures. This is in stark contrast to the corresponding two-dimensional lattice models with the same $\Delta$ parameter, where no long-range order exists at any finite temperature according to the Mermin-Wagner theorem.

On the two vertical phase boundaries (marked by the two dashed lines in Fig. 5), both $m_{x}$ and $m_{z}$ become finite, but their ratio is somewhat arbitrary. This indicates that the system is in a random mixture of $Z_{2}$-symmetry-breaking AF (or FM) and U(1)-symmetry-breaking XY phases. In other words, the U(1)$\otimes Z_{2}$ symmetry at $\Delta=-1$ or the SU(2) symmetry at $\Delta=1$ is broken, and the transitions from the planar-XY order to either FM- or AF-ordered phase is of the first order. Both $m_{x}$ and $m_{z}$ become discontinuous at $\Delta=\pm 1$, while $m$ remains continuous across these two high-symmetry points.

These two first-order phase transition lines can also be identified via other thermal measurements. Figure 6(a) shows the energy per bond $u_{b}$ as a function of $\Delta$ in the high-temperature paramagnetic phase with $T=0.625$, and deep in the symmetry-breaking phase with $T=0.125$. At $T=0.125$, the $u_{b}$ curve shows two turning points at $\Delta=\pm 1$, where the derivative $\partial u_{b}/\partial\Delta$ becomes discontinuous, as a consequence of the first-order phase transitions. On the other hand, at $T=0.625$, both $u_{b}$ and $\partial u_{b}/\partial\Delta$ vary smoothly, suggesting the absence of phase transitions vs. $\Delta$.

Similarly, the phase transitions among these three phases can be seen from the $\Delta$-dependence of the entanglement entropy $S_{E}$. As shown in Fig. 6(b), $S_{E}$ exhibits two peaks at $\Delta=\pm 1$ at low temperature $T=0.125$, which also becomes smooth at a high value $T=0.625$.

Lastly, the XXZ spin model can be mapped onto a hardcore boson model with nearest-neighboring (NN) interactions, by setting $S_{i}^{-}=b_{i}$ and $S_{i}^{z}=b_{i}^{\dagger}b_{i}-1/2$, where $b_{i}$ is a hardcore boson operator. In the boson language, the XY phase corresponds to a superfluid phase (SF) with an off-diagonal long-range order, and the FM and AF phases correspond to a Mott insulator (MI) and a charge density wave (CDW) phases, respectively. Along the two vertical phase boundary lines from low temperatures up to $T_{c}$, the SF phase coexists with the MI or CDW phases.

III.2 Low-lying excitations and quasi-Goldstone modes

From the inset of Fig. 4(a), it is clear that the specific heat $c_{V}$ decays exponentially with $\beta$ in low temperatures, no matter in which magnetically ordered phase, i.e.,

[TABLE]

This indicates that there is a finite excitation gap in the low-lying energy spectrum, quantified by the exponent $\Delta_{m}$ in the above equation. The values of $\Delta_{m}$ (Table 1) can be obtained by fitting the low-$T$ results of the specific heat $c_{V}$ or the internal energy $u$. For example, for the FM phase at $\Delta=-1^{-}$, the low-temperature internal energy is approximately described by the formula

[TABLE]

where $u_{0}=-0.375-h/2$ is the ground state energy per site, and $h$ is an external magnetic field, and $g(\beta)$ is some polynomial prefactor. From the derivative of this equation

[TABLE]

and through a polynomial fitting, the magnon gap $\Delta_{m}$ can be readily read out from the intercept at $T=1/\beta=0$. Moreover, for this FM system, one can further tune the magnon gap by changing the external magnetic fields $h$. The FM excitation gaps of various (small) magnetic fields $h=0,0.1,0.2$ are also obtained by fitting the internal energy curves and shown in Table 1, from which we find that $\Delta_{m}(h)-\Delta_{m}(0)\simeq h$. This suggests that the elementary excitations in the FM phase are magnons with spin $S=1$.

To further clarify the nature of low-lying excitations, we have evaluated the correlation length from the transfer matrix along the path $...\,B\overset{z}{\rightarrow}A\overset{x}{\rightarrow}B\overset{z}{\rightarrow}...$ on the Bethe lattice, meaning going through the $z$ bond of the $T_{B}$ tensor to $T_{A}$, and then through its $x$ bond to the next $T_{B}$, and so on. There are also other paths that can be used to define the transfer matrix. But all these paths are physically equivalent and the correlation lengths obtained thereof are also found exactly the same.

Suppose $\eta_{0}$ and $\eta_{i}$ are the dominant and the ($i+1$)’th largest eigenvalue of the transfer matrix, the correlation length can be determined as

[TABLE]

where $\xi_{1}$ is the largest correlation length the system can have, and $\xi_{i>1}$ is related to a shorter-ranged correlation function.

Alternatively, one can also estimate the correlation length $\xi$ directly from the real-space spin-spin correlation function

[TABLE]

where $l$ measures the distance between the two sites along the path. Given the $C^{\zeta}(l)$ data, the correlation length $\xi$ can be obtained by fitting the large-distance correlation function with an exponentially varying function.

Fig. 7(a-d) shows the correlation lengths obtained with the above two approaches. The correlation length $\xi^{z(x)}$ determined from the dominating correlators, $C^{z}(l)$ for $\Delta=\pm 1,1.2$ and $C^{x}(l)$ for $\Delta=0$, exhibits cusps exactly at the critical temperatures $T_{c}$. The value of the correlation length at the critical point, $\xi_{c}=1/\ln(2)\simeq 1.4427$, equals the critical upper bound of the correlation length on the $z=3$ Bethe lattice. As discussed in Ref. Li et al., 2012, the number of spins that correlate with a root spin at a given site grows exponentially with their distances. Therefore, a magnetic order-disorder phase transition with a diverging magnetic susceptibility happens when $\xi$, though being finite, reaches the critical value $1/\ln{(z-1)}$.

For the cases $\Delta=\pm 1$ shown in Fig. 7(a,c), the correlation lengths, $\xi^{x,y,z}$ along all three directions, are all equal to the dominant correlation length $\xi_{1}$ in the paramagnetic phase above $T_{c}$. Below $T_{c}$, the transverse excitation modes remain critical and do not change with temperature, i.e., the corresponding correlation lengths still equal the dominant one, $\xi^{x,y}=\xi_{1}=\xi_{c}$. Moreover, the longitudinal $\xi^{z}$ drops immediately below $T_{c}$ with decreasing temperature, due to the formation of magnetic order along the $z$-direction. Similar critical behavior has been observed in Fig. 7(b) with $\Delta=0$, where the degeneracy between $\xi^{x}$ and $\xi^{y}$ is broken, since the ordering of spins happens along the $x$-direction.

For the case $|\Delta|>1$ in Fig. 7(d), $\xi^{z}=\xi_{1}$ is still the dominant correlation length above $T_{c}$, which reaches the critical value $\xi_{c}$ and drops below $T_{c}$. On the other hand, $\xi^{x,y}$ never reaches the critical value $1/\ln{2}$, although they do not decay below $T_{c}$, and surpasses $\xi^{z}$ at a temperature below $T_{c}$.

The peculiar behaviors of the transverse correlation lengths observed in the ordered phases are quite remarkable. It indicates that although the true Goldstone modes are absent in the continuous symmetry breaking phases Laumann et al. (2009), the transverse excitations remain “critical” with a maximal correlation length $\xi_{c}=1/\ln{(2)}$ on the Bethe lattice. These quasi-critical transverse excitation modes are reminiscent of the Goldstone modes in a truly gapless continuous-symmetry-breaking system, and they can be regarded as somewhat “renormalized” Goldstone modes.

Distinct from the true gapless modes, and as the finite correlation lengths imply, these quasi-Goldstone modes are always gapped, and they become activated only above certain finite energy scales/temperatures, giving rise to the finite-$T$ phase transitions. Note that such kind of quasi-Goldstone modes are absent in the $Z_{2}$-symmetry-breaking phase, again similar to the two-dimensional lattices where Goldstone modes are absent for $\left|\Delta\right|>1$.

As a complementary to correlation length data, we show in Figs. 7(e-h) the absolute value of NN correlators $|\langle S_{i}^{\zeta}S_{i+1}^{\zeta}\rangle|$ ($\zeta=x,y,z$). In Fig. 7(e,g), we see distinct behaviors of NN correlators between the FM ($\Delta=-1$) and AF ($\Delta=1$) cases. In the AF phase, $|\langle S_{i}^{x(y)}S_{i+1}^{x(y)}\rangle|$ approaches a finite value in the zero temperature limit. However, in the FM case, $|\langle S_{i}^{x(y)}S_{i+1}^{x(y)}\rangle|$ decays exponentially as $T$ decreases and it approaches zero in the $T=0$ limit. This is consistent with the fact that the FM ground state is simply a direct-product state $|\uparrow,\uparrow,...,\uparrow\rangle$, while the AF ground state bears quantum entanglement. Moreover, as shown in Figs. 7(f,h), $|\langle S_{i}^{x}S_{i+1}^{x}\rangle|$ and $|\langle S_{i}^{z}S_{i+1}^{z}\rangle|$ are strongest NN correlators in the $\Delta=0$ and $1.2$ cases, respectively. In addition, other spin correlators in these two cases converge to finite values in the $T=0$ limit.

IV Linear Spin Wave Theory

Here we provide a linear spin wave analysis of the XXZ model on the Bethe lattice. The results reveal more information on the low-lying excitations, especially on the quasi-Goldstone mode in the continuous symmetry breaking phase.

IV.1 The ferromagnetic Heisenberg model

We first consider the excitations in the FM phase with $\Delta<-1$. By taking a U(1) rotation $\exp\left(\mathrm{i}\pi S_{n}^{z}\right)$ for all the spins on the $B$ sublattice, we can transform Eq. (1) into the form

[TABLE]

To carry out the linear spin-wave expansion, we choose an arbitrary site as a “root”, labelled as $R_{0}=0$, and then there exists an unique path that connects it to any other lattice site, rendering a convenient way to label the Bethe lattice sites. More specifically, we label the lattice according to the rule shown in Fig. 8: a site in the $l$-th layer from the root is represented by $l$ indices $R_{l}=\left(r_{1}r_{2}...r_{l}\right)$, where $r_{\lambda}=0,1,...,\tilde{z}-1$ ($\lambda=1,2,..,l$) denotes the way a path can choose at the $(\lambda-1)$-th branching point. On a $z$-coordinated Bethe lattice, $\tilde{z}=z$ for the $l=1$ layer and $\tilde{z}=z-1$ for the rest layers.

Assuming all spins are up polarized in the ground state, we exploit the Holstein-Primakoff (HP) transformation and take the leading approximation for the spin operators, $S^{+}=\sqrt{2S}a$, $S^{-}=\sqrt{2S}a^{\dagger}$ and $S^{z}=S-a^{\dagger}a$, where $a$, $a^{\dagger}$ are boson annihilation and creation operators, respectively. Under the linear spin-wave approximation, the Hamiltonian becomes

[TABLE]

where $n_{R_{l}}=a_{R_{l}}^{\dagger}a_{R_{l}}$ is the particle number operator, and a constant energy term is dropped for the sake of simplicity.

To diagonalize the Hamiltonian (18), we first take the following multidimensional discrete Fourier transformation

[TABLE]

where $Q_{l}=\left(q_{1}q_{2}...q_{l}\right)$ are the “quasi-momenta” which are dual to $R_{l}$, and $U$ is a unitary matrix

[TABLE]

Through this transformation, we generate a Bethe lattice with exactly the same geometry in the $q$-space, i.e., $q_{l}\in\{0,1,...,\tilde{z}-1\}$.

Under the above transformation, the hopping term in Eq. (18) becomes

[TABLE]

where $\left(Q_{l},0\right)$ denotes a site on the $(l+1)$-th shell, next to site $Q_{l}$ and with $q_{l+1}=0$. On the other hand, the occupation number term remains as a number operator in the $q$-space

[TABLE]

More details of the transformation can be found in App. C.

Equation (21) with its hermitian conjugate implies that a boson on the site $Q_{l}$ can only hop to site $\left(Q_{l},0\right)$, and vice versa. $\left(Q_{l},0\right)$ indicate a site on $(l+1)$-th shell, which is on the same branch of $Q_{l}$ and further has $q_{l+1}=0$. Therefore, for any given site $Q_{l}$ with $q_{l}\neq 0$, there exists a one-dimensional path, consisting of $Q_{l}$, $\left(Q_{l},0\right)$, $\left(Q_{l},0,0\right)$, $...$, etc.

The Hamiltonian on this half-infinite one-dimensional chain is tridiagonal, i.e.,

[TABLE]

where $b_{i}=a_{Q_{l}\underbrace{00...0}_{i}}$, and $b_{0}=a_{Q_{l}}$ is the starting node from which the one-dimensional path, depicted by a thick line in Fig. 8, is defined.

Secondly, we solve the above chain Hamiltonian via a conventional Fourier transformation, and the spin wave energy spectrum is found to be

[TABLE]

where $\kappa\in[0,\pi]$ is the momentum along the “effective” chains on the $q$-space Bethe lattice. From Eq. (24), as well as the plot in Fig. 9, we find the magnon excitation energy gap as

[TABLE]

The single magnon state is an eigenstate of the original Heisenberg FM model, and thus the magnon gap $\Delta_{m}$ obtained in Eq. (25) constitutes an upper bound of the true excitation gap of the system. Matter of fact, as shown in Tab. 1, the energy gap obtained with this equation, $\Delta_{m}=3/2-\sqrt{2}\simeq 0.0858$, is very close to the value estimated from the TTN calculation. In addition, from Eq. (25) one can see clearly that $\Delta_{m}$ increases linearly with $h$, in excellent consistency with the TTN results (see also Tab. 1).

From the above discussion, we notice that each branch of magnon excitations is confined to a one-dimensional path which is formed by a symmetrized ($q=0$) superposition of $\tilde{z}$ real-space sites. This symmetry has already been exploited in some previous works Brinkman and Rice (1970); Lepetit et al. (1999); Mahan (2001); Petrova and Moessner (2016), however, the $r$-$q$ transformation we introduce here works in a more generic way.

IV.2 Antiferromagnetic Heisenberg model

For the AF model with $\Delta\geq 1$, the linear spin wave expansion can be similarly done by explicitly considering the spin orientations on the two sublattices. We start from the configuration that all the spins are upward aligned on the $A$-sublattice and downward aligned on the $B$-sublattice. Accordingly, we take a two-sublattice HP transformation, i.e., $S_{i}^{+}=\sqrt{2S}a_{i}$, $S_{i}^{-}=\sqrt{2S}a_{i}^{\dagger}$ and $S_{i}^{z}=S-a_{i}^{\dagger}a_{i}$, for $i\in A$, and $S_{j}^{+}=\sqrt{2S}a_{j}^{\dagger}$, $S_{j}^{-}=\sqrt{2S}a_{j}$ and $S_{j}^{z}=a_{j}^{\dagger}a_{j}-S$ for $j\in B$. Therefore, the Hamiltonian under the linear spin wave approximation can be written as

[TABLE]

A constant term is again omitted in obtaining the above expression.

To cope with the two-sublattice structure, we need to adopt different transformations on diffferent sublattices, namely taking $\tilde{U}_{q,r}=U_{q,r}$ on the $A$-sublattice and $\tilde{U}_{q,r}=U_{q,r}^{\dagger}$ on the $B$-sublattice. Under this transformation, the pair annihilation term becomes (App. C)

[TABLE]

The above Hamiltonian can then be effectively represented as a direct sum of the infinite-many one-dimensional models defined by the model

[TABLE]

Employing the standard Bogoliubov transformation, the energy spectrum of $H_{\mathrm{AFC}}$ is found to be (and plotted in Fig. 9)

[TABLE]

The energy gap is, $\Delta_{m}=S(z-2)$, for $z\geq 3$. Since the single-magnon excitation state is not an eigenstate of the original AF model, $\Delta_{m}$ may not be an upper bound of the true excitation gap. Nonetheless, from Tab. 1 we observe that $\Delta_{m}$ still provides a quite good estimate of the magnon gap.

IV.3 The XY model

The magnon bands of FM and AF cases become unstable when $|\Delta|<1$ (more precisely, when $|\Delta|<2\sqrt{z-1}/z$). To perform a linear spin-wave analysis for the XY phase, we needs to start from a classical state in which the spins are ordered on the XY plane.

Below, for the sake of simplicity, we consider only the case $\Delta=0$. The Hamiltonian can be equivalently written as

[TABLE]

Under the linear spin-wave approximation, it becomes

[TABLE]

To solve this problem, we use the following $r$-$q$ transformation matrix

[TABLE]

This particular choice ensures $U_{q,r}$ to be a real orthogonal matrix with the property

[TABLE]

From this, again we obtain an effective one-dimensional boson model

[TABLE]

Diagonalizing this Hamiltonian using the Bogoliubov transformation (App. D), the magnon excitation energy is found to be (Fig. 9)

[TABLE]

The excitation gap is $\Delta_{m}=S\sqrt{z}(\sqrt{z-1}-1)$, which is $(\sqrt{6}-\sqrt{3})/2\approx 0.3587$ for $S=1/2$ and $z=3$, close to the numerical results, $0.41(1)$, shown in Tab. 1.

V Summary

We have investigated the Heisenberg XXZ model using both the canonical TTN and linear spin-wave theory on the $z=3$ Bethe lattice. Through efficient and accurate tensor network simulations, we have obtained the finite-temperature phase diagram of the model. The system undergoes a second-order phase transition at finite temperature, and three kinds of magnetic ordered phases are uncovered in low $T$. The correlation lengths, as well as bipartite entanglements, though exhibiting their maximal values at the transition temperature $T_{c}$, are found to be always finite on the Bethe lattice. Correspondingly, the low-lying excitations are revealed, through both numerical TTN and analytical spin wave calculations, to be gapped, even in the parameter range where the system breaks spontaneously the continuous symmetries.

Therefore, although the conventional gapless Goldstone modes are absent, quasi-Goldstone transverse fluctuation modes have been observed in the Bethe lattice XXZ model. When the system spontaneously breaks the continuous symmetries, the transverse correlation lengths reach the “critical” value, $\xi_{c}=1/\ln{(z-1)}$, and remain there for $T\leq T_{c}$. Here $\xi_{c}$ is the maximal correlation length that is allowed on the Bethe lattice.

The results obtained on the Bethe-like lattices can be used to understand physical properties of quantum lattice models on the honeycomb, square, or other regular lattices. In addition, the canonical TTN method we proposed works very generally, and it can be applied to other fundamental quantum many-body models, such as the frustrated Heisenberg and Hubbard models defined on the Bethe-like lattices, etc.

VI acknowledgements

This work was supported by the National Natural Science Foundation of China (11834014, 11888101), the National R&D Program of China (2017YFA0302900). WL is indebted to Andreas Weichselbaum for helpful discussions. DWQ and WL would like to thank Jan von Delft for hospitality during a visit to LMU Munich, where part of this work was performed.

Appendix A Tensor update on the Bethe lattice

A.1 Bethe lattice update

An imaginary time evolution can be employed to cool down the TTN density operator on the Bethe lattice, where the local tensors are updated via a simple scheme. The details are illustrated in Fig. A1 and elaborated below. We take the the $x$-bond update as an example, and the successive projection procedures on the other two, i.e., $y$- and $z$-bonds, can be accomplished similarly. The three projection substeps on different bonds constitute a full Trotter step of small imaginary-time slice $\tau$.

(a) Absorb the environment matrices $\lambda_{y}$ and $\lambda_{z}$ to the $T$ tensors and construct the $W$ tensor [Fig. A1(a)]

[TABLE]

where $\alpha=A,B$.

(b) Perform a QR decomposition of $W_{\alpha}$ [Fig. A1(b)]

[TABLE]

which splits off the upper physical indices $\sigma_{\alpha}$ into the $R$ tensors.

(c) Construct a base tensor $B$ by combining $\lambda_{x}$ with the two adjacent tensors $R_{A}$ and $R_{B}$

[TABLE]

and apply the two-site imaginary-time evolution gate $P=e^{-\tau h_{x}}$ onto the base tensor [Fig. A1(c)]

[TABLE]

(d) Reshape $\tilde{B}$ into a matrix by grouping together $\sigma_{A}$ with $x_{A}^{\prime}$ as a single index, and $\sigma_{B}$ with $x_{B}^{\prime}$ into another. Then we perform a matrix SVD

[TABLE]

as shown in Fig. A1(d). Note that the matrix dimension of $\tilde{\lambda}_{x}$ is enlarged by $d^{2}$ times after the bond evolution, which needs to be truncated by retaining only the largest $D$ singular values and corresponding bond bases. After this proper truncation of bond states, we update the tensors $\tilde{R}_{a}$, $\tilde{R}_{b}$ and $\tilde{\lambda}_{x}$ accordingly.

(e) Combine $\tilde{R}$ to the corresponding $Q$ tensors, and spit off the $\lambda_{y(z)}$ matrices, that have been absorbed into $Q$ in steps (a,b), by multiplying their inverse matrices $\lambda^{-1}_{y(z)}$ to $Q$. Thus we obtain the updated $\tilde{T}$ tensors [Fig. A1(e)]

[TABLE]

where again $\alpha=A,B$ denoting two sublattices.

In the groundstate optimization Li et al. (2012), the above procedure suffices to produce accurate results. The truncation errors do not accumulate, and $\tau$ can be tuned to a sufficiently small value, say, $\tau\sim 10^{-4}$ in the end of projections. However, in the finite-temperature simulations, now the truncation errors accumulate, we therefore need to carefully optimize the truncations in every single step to improve the overall performance. Note that $e^{-\tau h_{\zeta}^{i}}$ ($\zeta=x,y,z$) is not unitary and breaks the orthogonality of bond basis, and thus TTN deviates the canonical form after each step of imaginary-time evolution. Therefore, a canonicalization procedure of the TTN is needed to restore the orthogonality of bond bases and optimize the truncations.

A.2 Canonicalization procedure of TTN

Following a similar line of standard procedure developed in the matrix product Vidal (2007); Orús and Vidal (2008); Orús (2014) as well as tensor product states Ran et al. (2012), we present below the canonicalization of TTN on the Bethe lattice.

(i) As shown in Fig. A2(a), we decompose the dominant eigenvectors of the transfer tensors as

[TABLE]

where $\alpha=A,B$. Since in practice $X_{\alpha},Y_{\alpha}$, and $Z_{\alpha}$ matrices are symmetric, the decomposition can be done via the eigenvalue or the Cholesky decomposition.

(ii) Insert pairs of reciprocal matrices, $u_{\alpha}u_{\alpha}^{-1}$, $v_{\alpha}v_{\alpha}^{-1}$ and $w_{\alpha}w_{\alpha}^{-1}$, to the three geometrical bonds, as shown in Fig. A2(b). Order of the matrix multiplication has also been specified by the arrows, e.g., we contract the first index of $u_{A}$ with $\lambda_{x}$ and the second index of $u_{A}^{-1}$ with $T_{A}$.

(iii) Now we perform the bond update by combining the bond matrices and then perform a SVD. As shown in Fig. A2(c), we take the $x$-bond as an example, i.e.,

[TABLE]

$U_{x}$ and $V_{x}$ are unitary matrices, and $\lambda^{\prime}_{x}$ is used to update the bond diagonal matrix.

(iv) As shown in Fig. A2(d), we absorb $u^{-1}$, $v^{-1}$ and $w^{-1}$ matrices, as well as the adjacent unitary matrices $U$ and $V$, into the $T_{\alpha}$ tensors, and update $T_{A}^{\prime}$ and $T_{B}^{\prime}$.

After the above procedure on the $x$ bond (and simultaneously on the $y$ and $z$ bonds), the updated tensors $T_{\alpha}^{\prime}$ satisfy the canonical conditions [see Eq. (2) of the main text], and the dominating eigenvectors $X_{\alpha},Y_{\alpha},Z_{\alpha}$ are now gauged into identities.

A.3 Simple vs. canonical schemes

Now we provide some numerical benchmarks of the simple and canonical update schemes, showing the advantage of the latter in both the accuracy and robustness. Here by canonical scheme we mean the combined procedure using techniques introduced in both Secs. A.1 and A.2; while by simple scheme we mean a poorman’s approach where the canonicalization operations in Sec. A.2 are skipped.

We consider the Heisenberg model ($\Delta=1$), and compare the entanglement entropy $S_{E}$ and the number of iterations $N_{\rm{iter}}$ required to reach a convergence in determining the dominant eigenvectors in Eq. (8) of the main text. Although the entanglement entropy $S_{E}$ defined in Eq. (11) is rigorously defined only for canonical TTN, we can nevertheless take “$S_{E}$” from the simple scheme as a measurement of entanglement for comparisons.

As seen in Fig. A3(a), the two $S_{E}$ curves almost coincide for $T>T_{c}$. However, their behaviors start to differ at the critical temperature $T_{c}$. The curve of the canonical scheme shows a cusp at the critical temperature and slowly converges to a smaller zero-temperature entanglement value, while that of the simple scheme still rises smoothly until collapsing at certain lower temperature below $T_{c}$. After this “jump”, the simple scheme curve lies almost on top of the canonical curve, and both converge to the $T=0$ entanglement value in nearly the same rate.

Correspondingly, as shown in Fig. A3(b), $N_{\rm{iter}}$ peaks at $T_{c}$ in the canonical scheme, and it peaks both at $T_{c}$ and the “jump” point at a lower $T$ in the simple scheme curve. The second peak in $N_{\rm{iter}}$ in the simple scheme belongs to an numerical “artifact”, since no phase transition really takes place there.

Apart from the “artifact” in $S_{E}$, the more severely accumulated errors in the simple scheme, as the temperature lows down, may cause other problems. In practice, sometimes the simple scheme is found to generate “wrong” metastable thermodynamic results, e.g, magnetic moments, at low temperatures $T<T_{c}$.

To conclude, the canonical scheme turns out to be more accurate and robust, and it is thus mostly adopted in our practical simulations.

Appendix B Classical Ising model on the Bethe lattice

In this appendix, we provide the rigorous solution of the classical Ising model on the Bethe lattice via transfer tensor techniques. The TTN algorithms introduced in Sec. II can also be employed to compute this classical model, and the comparisons between the numerical and rigorous results thus provide a first benchmark of the TTN algorithm.

The Hamiltonian (energy) of the classical Ising model reads $H=-J\sum_{\langle i,j\rangle}\sigma_{i}\sigma_{j}$, where $J=1$ is the energy scale, $\langle i,j\rangle$ means NN sites on the Bethe lattice, and $\sigma=\pm 1$ denotes the classical Ising variables.

This Ising model can be solved exactly through a number of essentially equivalent methods, including the self-similarity Baxter (2013); Ostilli (2012) and cavity approaches Mézard and Parisi (2001); Ostilli (2012), etc. By solving the generalized eigenvalue problem, we find the dominant eigenvectors of the transfer tensor $\mathcal{T}$, from which we further obtain the exact expression of thermal quantities.

Firstly, we rewrite the partition function as a TTN, i.e.,

[TABLE]

where $\mathcal{T}$ is a transfer tensor, and $\sigma_{1(2)}$ sits in the two central sites in the bond-centered [see, e.g., Fig. 1(b)]. The $r_{1}r_{2}...r_{l}$ labelling in $\sigma_{r_{1}r_{2}...r_{l}}$ follows the convention shown in Fig. 8, but $r_{i}$ starts from 1 instead of 0 here.

Efficient contractions of this infinite TTN can be implemented by finding the dominant (generalized) eigenvectors $V_{\sigma}$ of the transfer tensor $\mathcal{T}^{\sigma_{1}\sigma_{2}}_{\sigma_{0}}=e^{\beta\sigma_{0}(\sigma_{1}+\sigma_{2})}$, i.e.,

[TABLE]

with the (generalized) eigenvalue $\eta$. By writing it down explicitly, we have

[TABLE]

Since Eq. (A11) does not constitute a linear system of equations, one has a gauge degrees of freedom in determining the eigenvalue $\eta$ as well as the eigenvector $V_{\pm}$. Nevertheless, we can eliminate $\eta$ from the equations, and arrive at a cubic equation after some rearrangement,

[TABLE]

where $x=\sqrt{V_{+}/V_{-}}$. Note that the eigenvalue can now be expressed as $\eta=V_{+}(x/e^{\beta}+e^{\beta}/x)^{2}$, which can not be uniquely determined and depends on the normalization condition of $V_{\sigma}$ (see related discussions in Sec. II.3).

The (unnormalized) probability distribution for two neighboring spin variables is $\rho_{\sigma_{1}\sigma_{2}}=e^{\beta\sigma_{1}\sigma_{2}}V_{\sigma_{1}}V_{\sigma_{2}}$. By summing over $\sigma_{2}$, we can get the distribution of a single spin $\sigma_{1}$, from which the magnetization can be derived as $\langle\sigma\rangle=\sum_{\sigma_{1}}\sigma_{1}\rho_{\sigma_{1}}/\sum_{\sigma_{1}}\rho_{\sigma_{1}}=(x^{3}-1)/(x^{3}+1)$. Note that $x$ is just the root of Eq. (A12), and thus the local magnetization can be uniquely determined.

The order parameter, i.e., the spontaneous magnetization $m$, reads

[TABLE]

with $\omega=e^{2\beta}$ and the critical point $\beta_{c}=1/T_{c}=\ln(3)/2\simeq 0.5493$.

Apparently, the root $x=1$ (and thus $m=0$) corresponds to the paramagnetic solution, while the two roots of the remaining quadratic equation in Eq. (A12) reflects the two-fold degenerate FM states. The $\pm$ sign represents the spin-up and spin-down solutions, respectively, given that the discriminant $(1-e^{2\beta})^{2}-4>0$, i.e., $\beta>\beta_{c}=\ln(3)/2$.

The internal energy per bond $u_{b}$ can be determined from $\rho_{\sigma_{1}\sigma_{2}}$, which is

[TABLE]

It is clear that the $u_{b}$ curve exhibits a singular point at the transition temperature $\beta_{c}$.

On the other hand, the Bethe-lattice Ising model can also be solved by the TTN techniques. By performing a decomposition following Eq. (3), and then an imaginary-time evolution procedure, we obtain the thermal density matrix $\rho(\beta)=e^{\beta\sum_{\langle i,j\rangle}\sigma_{i}\sigma_{j}}$ of the classical Ising model on the Bethe lattice. On top of that, we can further calculate the thermal quantities, including the magnetization, energy, and specific heat, etc.

In Fig. A4, we compare the TTN results of the spontaneous magnetization and the specific heat to the exact solution, where excellent agreements can be observed. This can be ascribed to the absence of Trotter errors in the calculations, and also to no essential truncations in the procedure of cooling.

It is also of interest to note, in Eq. (A13), that the spontaneous magnetization $m\sim(\frac{T_{c}-T}{T_{c}})^{\kappa}$ with $\kappa=1/2$, when $T$ approaches the transition temperature $T_{c}$ (from low temperature side). In addition, the internal energy $u_{b}$ curve is continuous at $T_{c}$, while the slope $c_{V}\equiv\partial u_{b}/\partial T$ shows a discontinuity in Fig. A4(b), suggesting a mean-field-type critical exponent $\alpha=0$.

Appendix C The $r$-$q$ transformation

Here we provide more details on the $r$-$q$ transformations of the bilinear terms, including the hopping, on-site occupation number, pair-creation or annihilation terms, etc.

Firstly, we check that the real-space on-site term remains as on-site term in the $q$-space

[TABLE]

For the NN hopping term, we have

[TABLE]

where $\tilde{z}=z-1$ for $l>1$ layers and $\tilde{z}=z$ for $l=1$ one.

In the HAF model, we have pair-creation and annihilation operators, where additional care needs to be taken of, i.e.,

[TABLE]

Therefore, the pair creation and annihilation operators can be transformed into $q$-space in a well organised way, e.g.,

[TABLE]

which again falls into an effective 1D chain geometry.

Appendix D The Bogoliubov transformation in the XY model

In this appendix, we provide the details of the Bogoliubov transformation for the XY model. We start from the 1D half-chain bosonic Hamiltonian, i.e., Eq. (33) in the main text. By ignoring the “impurity site” at the center, performing a Fourier transformation on an infinite chain (without changing the energy dispersion curve), we arrive at $H_{\rm{1D}}=\sum_{\kappa\geq 0}H_{\kappa}$, where

[TABLE]

with $\gamma_{\kappa}=\sqrt{z-1}\operatorname{cos}\kappa$.

The sub-Hamiltonian $H_{\kappa}$ can be rewritten as $H_{\kappa}=\frac{1}{2}B^{\dagger}\mathcal{H}_{\kappa}B$, where $B=(b_{\kappa},b_{-\kappa},b_{\kappa}^{\dagger},b_{-\kappa}^{\dagger})^{T}$, and

[TABLE]

where we have omitted the term $-(z+\gamma_{\kappa})$ which is a constant after summing over $\kappa$.

In the Bogoliubov transformation, we find a matrix $\Theta$ so that (1) $B=\Theta\tilde{B}$, where $\tilde{B}=(\beta_{\kappa},\beta_{-\kappa},\beta_{\kappa}^{\dagger},\beta_{-\kappa}^{\dagger})^{T}$ still represents boson operators obeying the bosonic statistics, and (2) $H_{\kappa}=\frac{1}{2}\tilde{B}^{\dagger}\Theta^{\dagger}\mathcal{H}_{\kappa}\Theta\tilde{B}=\frac{1}{2}\tilde{B}^{\dagger}\mathbf{D}\tilde{B}$ is in a diagonal form.

To maintain the boson statistics in condition (1), we require $\Theta^{\dagger}\Lambda_{z}\Theta=\Lambda_{z}$, where

[TABLE]

and $\mathbf{I}_{2\times 2}$ is a $2\times 2$ identity matrix.

Combing together conditions (1) and (2), we have $(\Lambda_{z}\mathcal{H}_{\kappa})\Theta=\Theta\Lambda_{z}\Theta^{\dagger}\mathcal{H}_{\kappa}\Theta=\Theta\Lambda_{z}\mathbf{D}$, which implies $\Theta$ and $\mathbf{D}$ can be found by diagonalizing $\Lambda_{z}\mathcal{H}_{\kappa}$. After some calculations, and by observing that the Hamiltonian matrix in Eq. (A20) is block diagonal, we arrive at

[TABLE]

where the positive root with $\kappa>0$ constitutes the magnon spectrum in Eq. (34), i.e., $\epsilon(\kappa)=\tilde{\epsilon}(\kappa)$. The corresponding Bogolon (annihilation) operator turns out to be

[TABLE]

with restriction $\kappa\in[0,\pi]$, where $\theta_{\kappa}$ can be determined from $\operatorname{tanh}(\theta_{\kappa})=-\gamma_{\kappa}/(z+\gamma_{\kappa})$.