Detection and characterization of symmetry-broken long-range orders in the spin-$\frac{1}{2}$ triangular Heisenberg model

TL;DR

This paper introduces advanced numerical tools to analyze symmetry-breaking phases in the triangular Heisenberg model, revealing various magnetic orders and spin liquids, and demonstrating the robustness of topological phases.

Contribution

It develops new methods using MPS and DMRG to detect symmetry-breaking and characterize phases in frustrated quantum magnets, including identifying spin liquids and their stability.

Findings

Identified 120° magnetic order, columnar order, and algebraic spin liquid phases.

Demonstrated robustness of topological spin liquid against chiral perturbations.

Developed criteria for detecting symmetry-breaking transitions via entanglement spectrum and cumulants.

Abstract

We present new numerical tools to analyze symmetry-broken phases in the context of $SU(2)$-symmetric translation-invariant matrix product states (MPS) and density-matrix renormalization-group (DMRG) methods for infinite cylinders, and determine the phase diagram of the geometrically-frustrated triangular Heisenberg model with nearest and next-nearest neighbor (NN and NNN) interactions. The appearance of Nambu-Goldstone modes in the excitation spectrum is characterized by "tower of states" levels in the momentum-resolved entanglement spectrum. Symmetry-breaking phase transitions are detected by a combination of the correlation lengths and second and fourth cumulants of the magnetic order parameters (which we call the Binder ratio), even though symmetry implies that the order parameter itself is strictly zero. Using this approach, we have identified $120^{\circ}$ order, a columnar order, and an algebraic spin liquid (specific to width-6 systems), alongside the previously studied topological spin liquid phase. For the latter, we also demonstrate robustness against chiral perturbations.

Equations

$$H_{J_2}=J_1\sum _{⟨i,j⟩}{\mathbf{S}}_i\cdot {\mathbf{S}}_j+J_2\sum _{⟨⟨i,j⟩⟩}{\mathbf{S}}_i\cdot {\mathbf{S}}_j,$$

$${\mathbb{T}}_X(E^a)=\sum _{s^{\prime }s}⟨s^{\prime }|\hat{X}|s⟩A^{s^{\prime }†}{E}^a{A}^s,$$

$$E^{a^{\prime }}(n)\equiv \sum _{s^{\prime },s,a}{\mathcal{A}}^{s_n^{\prime }†}{M}_{a^{\prime }a}^{s^{\prime }s}{\mathcal{A}}^{s_n}{E}^a(n-1).$$

$$E^i(L+1)={\mathbb{T}}_{M_{ii}}(E^i(L))+\sum _{j>i}{\mathbb{T}}_{M_{ij}}(E^j(L)),$$

$$M=\left(\begin{array}{cc}I& X\\ & I\end{array}\right).$$

$$M^2=\left(\begin{array}{cc}I& X\\ & I\end{array}\right)\otimes \left(\begin{array}{cc}I& X\\ & I\end{array}\right)=\left(\begin{array}{cccc}I& X& X& X^2\\ & I& 0& X\\ & & I& X\\ & & & I\end{array}\right)$$

$$\Rightarrow \left(\begin{array}{ccc}I& X& X^2\\ & I& 2X\\ & & I\end{array}\right),$$

$$E^1(L+1)={\mathbb{T}}_I(E^1(L)),$$

$$E^2(L+1)$$

$$=C_X+{\mathbb{T}}_I(E^1(L)),$$

$$E^3(L+1)={\mathbb{T}}_{X^2}(I)+2{\mathbb{T}}_X(E^2(L))+{\mathbb{T}}_I(E^3(L)).$$

$$C_X$$

$$E^2(L)$$

$$e_n^{(2)}(L+1)=c_n^{(2)}+{\eta }_n{e}_n^{(2)}(L).$$

$${\tilde{C}}_X$$

$${\tilde{E}}^2(L)$$

$$e_1^{(2)}(L+1)=e_1^{(2)}(L)+c_1^{(2)},$$

$${\tilde{E}}_{(n)}^2(L+1)={\tilde{C}}_{(n)}+{\eta }_n{\tilde{E}}_{(n)}^2(L),$$

$$(1-{\mathbb{T}}_I){\tilde{E}}^2(\infty )={\tilde{C}}_X,$$

$$X^4=(\mathbf{X}\cdot \mathbf{X}{)}^2+(\mathbf{X}\otimes \mathbf{X})\cdot (\mathbf{X}\otimes \mathbf{X}),$$

$$\text{SSF}(\mathbf{k},N=\infty )=\underset{N\to \infty }{\lim }\frac{1}{N}\sum _{i,i^{\prime }}^N⟨{\mathbf{S}}_i\cdot {\mathbf{S}}_{i^{\prime }}⟩e^{i\mathbf{k}\cdot ({\mathbf{r}}_i-{\mathbf{r}}_{i^{\prime }})},$$

$$⟨{\bar{\mathbf{M}}}^2(\infty )⟩={\bar{a}}_0+{\bar{a}}_1{L}_x^{-2}+...,$$

$$O_A^{\text{FM}}=\frac{2}{\sqrt{N_A(N_A+2)}}\sqrt{⟨{\mathbf{S}}_A^2⟩},$$

$$O^{\text{stag}}=\frac{1}{L}\sqrt{⟨{\mathbf{S}}_{\text{stag}}^2⟩}$$

$$\frac{{\xi }_i(m)}{L_u}=-\frac{1}{\ln |{\eta }_i(m)|},i=2,3,4,...,$$

$$S_{EE}={\beta }_0+{\beta }_1{L}_{cut}+\frac{N_G}{2}\log (L),$$

$$S_{EE}\simeq a_0(L_y)+a_1(L_y)\log (\xi ),$$

$$a_1(L_y)={\alpha }_0+{\alpha }_1{L}_y.$$

$$⟨M^n⟩=\sum _{i=1}^n{B}_{n,i}({\kappa }_{1}L,{\kappa }_{2}L,...,{\kappa }_{n-i+1}L),$$

$$⟨M⟩$$

$$⟨M^2⟩$$

$$⟨M^3⟩$$

$$⟨M^4⟩$$

$${\mathbf{M}}_{stag}=\sum _{i=1}^{L_y}(-1{)}^i{\mathbf{S}}_i,$$

Full text

I Introduction

The groundstate of the one-dimensional nearest-neighbor Heisenberg model (originally determined in an important work by Heisenberg Heisenberg (1928)), $H_{NN}=J\sum_{i}{\bf S}_{i}\cdot{\bf S}_{i+1}$, for $J<0$, exhibits long-range ferromagnetic (FM) order, which breaks the spins’ rotational symmetry (the $SU(2)$ group) and elementary excitations are spin-waves (also known as Nambu-Goldstone bosons or Magnons; see for example Schollwöck et al., 2004; Auerbach, 1994; Sachdev, 2011). The Bethe Ansatz Bethe (1931) can be employed (e.g. see Karbach et al., 1998) to study the antiferromagnet (AFM) spin-$\frac{1}{2}$ Heisenberg model, $J>0$, which demonstrates the absence of magnetic ordering in clear contrast to the FM case. Today, we know there exist no continuous-symmetry-broken long-range order (LRO) in any one-dimensional system. In fact, magnetism in 1D and few-leg ladders is peculiarly different to higher dimensions (where LROs exist; see below), since the magnetic ordering at zero temperature is suppressed by quantum fluctuations due to the same mechanism as described by Mermin-Wagner-Hohenberg theorem for finite-temperatures Mermin and Wagner (1966); Hohenberg (1967) (i.e. due to the low cost of creating quantum long-range fluctuations, which increases the entropy). In contrast to 1D, long range magnetic ordering is possible in 2D Heisenberg-type Hamiltonians; early examples arose from studying anisotropy (Frohlich and Lieb, 1978), the AFM Heisenberg model with $S=\frac{1}{2}$ (Anderson, 1952; Oitmaa and Betts, 1978; Huse and Elser, 1988; Huse, 1988; Kennedy et al., 1988a; Reger and Young, 1988), $S\geq\frac{3}{2}$ (Neves and Perez, 1986), large-$S$ values (Anderson, 1952) and for the antiferromagnetic XY (Oitmaa and Betts, 1978; Kennedy et al., 1988b) and XXZ (Kubo and Kishi, 1988; Okabe and Kikuchi, 1988) models for all spin magnitudes. For the majority of two-dimensional magnetic materials, if there exist no frustration, the groundstate exhibit (Lhuillier, 2005; Sachdev, 2011; Farnell et al., 2014) either ferromagnetism or antiferromagnetism (i.e. the well-known bipartite Néel order Néel (1948)). It is widely believed that the celebrated Landau symmetry-breaking theorem (Landau, 1937; Landau and Ginzburg, 1950) explains the physics behind all such conventional magnetic ordering: Hamiltonians such as $H_{NN}$ contain a set of symmetries which are absent in the groundstate, a feature known as spontaneous symmetry breaking (SSB). As a result of symmetry-breaking, a well-defined order parameter exists in the model that can be used to characterize the magnetic ordering, unambiguously. After uncovering the key mechanisms behind the conventional ordering (in particular, ferromagnetism and bipartite antiferromagnetism), the field of low-dimensional quantum magnetism enjoyed a new boost of attention aimed at understanding exotic phases of quantum matter that appear in frustrated one-dimensional (Tsvelik, 2007; Gu and Wen, 2009; Auerbach, 1994; Pollmann and Turner, 2012; Pollmann et al., 2010, 2012; Turner et al., 2011) and two-dimensional (Lhuillier, 2005; Balents, 2010; Chen et al., 2010; Powell and McKenzie, 2011; Farnell et al., 2014) systems. This happened partly due to the rise of the geometrically-frustrated antiferromagnets on non-bipartite Archimedean lattices (Lhuillier, 2005; Schollwöck et al., 2004; Diep, 2004; Lacroix et al., 2011; Farnell et al., 2014; Sachdev, 2011). Interestingly, the existence of geometrical frustration is enough by itself to often lead to the ‘melting’ of the magnetic ordering, stabilizing a family of nonmagnetic phases, collectively classified as spin liquids (also known as paramagnetic states) (Kalmeyer and Laughlin, 1987; Wen et al., 1989; Wen, 1991; Read and Sachdev, 1991; Balents, 2010; Chen et al., 2010; Wen, 2002). Such quantum liquids preserve all Hamiltonian symmetries and, consequently, their existence cannot be understood through Landau’s symmetry-breaking paradigm. The search for new, hidden order parameters has been challenging theorists for the last 20 years, and has led to the discovery of even more intriguing phases of the quantum matter. A canonical example is the discovery of the topological order Wen (2002, 2007); Chen et al. (2010), such as symmetry-protected topological (SPT) ordering (including the Haldane phase and the closely relevant Affleck-Kennedy-Lieb-Tasaki groundstate (Affleck et al., 1987, 1988)) and the intrinsic topological states (Wen, 1991; Chen et al., 2010) (including the $Z_{2}$-gauge groundstate of the toric code Kitaev (2003, 2006); Kitaev and Laumann (2009)), which can only exist in $D\geq 2$.

Magnetic ordering is often identified using the scaling behavior of a static two-point (or higher order) correlation function. For a spin system, a two-point correlator can be written in terms of the (principal) correlation length, $\xi$, as $|G_{2}(i,i^{\prime})|=\langle{\bf S}_{i}\cdot{\bf S}_{i^{\prime}}\rangle\sim C+e^{-\frac{r_{ii^{\prime}}}{\xi}}+\cdots$, where $C$ stands for a constant (which can be zero), $r_{ii^{\prime}}$ is the distance between sites, and ellipses represent faster decaying terms. Different type of ordering can be defined as follows. For magnetically disordered states, with no conventional order parameter (i.e. no broken symmetry), the correlation function decays to zero exponentially fast, $C$ is zero, $\xi$ is finite, and there is a bulk gap in the excitation spectrum. In this case, instead of symmetry breaking, we have symmetry protection, giving rise to the SPT order. Such an exponential drop is observed in the Haldane phase (as an example, see the original calculations by White and HuseWhite and Huse (1993)). For true LROs, such as Néel-type AFMs and the FM state with a conventional order parameter, the correlation function tends to a constant at large distances, $\xi\rightarrow\infty$. There exists another distinct long-range phase of the quantum matter, which is referred to as a quantum critical state (or a quasi-LRO). In such phases, the correlation function decays with a power-law with distance. Power-law decaying correlation functions can be approximated as the sum of many exponential functions, as occurs in the MPS ansatz (Rommer and Östlund, 1997; Östlund and Rommer, 1995), which again translates to having diverging $\xi$, consistent with the Bethe Ansatz’ prediction for the spin-$\frac{1}{2}$ Heisenberg chain. Critical states are common in 1D quantum magnetism and appear at a transition between two gapped disordered phases with different symmetries, when the gap necessarily closes; however, they can also stabilize in an extended region, as in the XY-phase of the anisotropic Heisenberg chain (Schollwöck et al., 2004; Sachdev, 2011).

High-accuracy numerical methods, such as exact diagonalization (ED), quantum Monte Carlo (QMC), (see Sandvik, 2010 for a review), and coupled cluster (Farnell et al., 2014) methods, are often used for low-temperature frustrated magnets, modeled as strongly-interacting spin Hamiltonians exhibiting many-sublattice groundstates. In this paper, we employ and expand the functionality of the finite DMRG (fDMRG) White (1992, 1993); McCulloch (2007); Schollwöck (2011), and the state-of-the-art infinite DMRG (iDMRG) McCulloch (2008); Schollwöck (2011) methods to characterize LROs of a geometrically-frustrated system, when the many-body states are constructed through the $SU(2)$-symmetric (non-Abelian) MPS and infinite MPS (iMPS) ansätze Perez-Garcia et al. (2007); McCulloch (2007, 2008); Schollwöck (2011), respectively. The latter is a translation-invariant MPS that allows the calculation of many useful quantities directly in the thermodynamic limit via transfer matrix methods. Currently, there exist few well-established numerical tools, in the context of non-symmetric DMRG, to identify LROs. In finite-system MPS studies, SSB needs to be treated carefully because in exact calculations SSB does not occur at allKoma and Tasaki (1994); Lhuillier (2005), as finite size effects induce a gap between states that would be degenerate in the thermodynamic limit. In practice, with finite-precision arithmetic symmetry breaking can occur when the finite-size gap is smaller than the characteristic energy scale set by the accuracy of the numerics (in MPS calculations, this is set by the energy scale associated with the basis truncation). This can be difficult to control, as symmetry breaking might occur as a side-effect of the numerical algorithm or it might require an additional perturbation. Infinite-size MPS (or very large finite MPS) are better behaved in this respect, where there are a variety of techniques; one can look at the scaling of the correlation length of the groundstate against MPS number of states, $m$, which distinguishes gapped and gapless states Tagliacozzo et al. (2008); Stojevic et al. (2015); McCulloch (2008), direct measurement of local magnetization order parameters, the entanglement entropyEisert et al. (2010), and the static spin structure factor (SSF – see below). However, when the Hamiltonian symmetries are preserved explicitly, the order parameter is zero by construction and a robust set of numerical tools for characterizing magnetic ordering is not readily available. Here, we introduce and verify the accuracy of two new numerical tools, in the context of $SU(2)$-symmetric iMPS/iDMRG, to characterize and locate phase transitions incorporating LROs in the triangular Heisenberg model (THM) on infinite cylinders. New tools include study of the cumulants (cf. Stuart et al., 2009 for the definitions and relevant discussions on the non-central moments and the cumulants in the context of the probability theory) and a Binder ratio of magnetization order parameters, and further developments on tower-of-states (TOS) level patterns in the momentum-resolved entanglement spectra (ES)Li and Haldane (2008); Kolley et al. (2013); Cincio and Vidal (2013).

The triangular lattice has the highest geometrical frustration in the Archimedean crystal family with a coordination number of $Z_{c}\!=\!6$. Anderson and FazekasAnderson (1973); Fazekas and Anderson (1974) argued that the high frustration of the triangular lattice might be enough to melt the long-range magnetic ordering observed for the Heisenberg model on the square lattice (e.g. see Lhuillier, 2005; Sachdev, 2011; Farnell et al., 2014). In the first work, Anderson conjectured that the spin-$\frac{1}{2}$ THM with antiferromagnetic NN bonds should stabilize a resonating-valance-bond (RVB) groundstate (i.e. the equally-weighted superposition of all possible arrangements of the singlet dimers on the lattice; RVBs are the building blocks of the quantum liquids). The failure of robust analytical and numerical studies to find an RVB groundstate motivates the search for a minimal extension to $H_{NN}$ that increases the frustration. The obvious choice is frustration through the addition of a NNN coupling term, which frustrates a LRO $120^{\circ}$-ordered arrangements of sublattices (see below and Huse and Elser, 1988; Deutscher and Everts, 1993; Capriotti et al., 1999; Jolicoeur and Le Guillou, 1989; Jolicoeur et al., 1990; Chubukov and Jolicoeur, 1992; Bernu et al., 1994; White and Chernyshev, 2007; Kaneko et al., 2014; Farnell et al., 2014; Li et al., 2015; Zhu and White, 2015; Hu et al., 2015; Iqbal et al., 2016). This led to the introduction of the $J_{1}$-$J_{2}$ THM, for which the Hamiltonian is defined as

[TABLE]

where $\langle i,j\rangle$ ($\langle\langle i,j\rangle\rangle$) indicates that the sum goes over all NN (NNN) couplings. The $SU(2)$-symmetry of $H_{J_{2}}$ can be simply realized by noticing $[H_{J_{2}},{\bf S}]=0$ (${\bf S}$ stands for the total spin vector), which means that eigenvalues of ${\bf S}^{2}$ are good quantum numbers and can be used to label groundstate symmetry sectors. Geometrical frustration forbids the bipartite Néel order as a stable groundstate of the antiferromagnetic NN model ($J_{1}>0$ and $J_{2}=0$). Consequently, one expects the groundstate, for the majority of the phase diagram of the antiferromagnetic $H_{J_{2}}$, to be a compromise, such as a $120^{\circ}$-ordered arrangement Saadatmand et al. (2015); Lhuillier (2005); Farnell et al. (2014). By now, it is well-known that the groundstate of the nearest-neighbor THM does not exhibit an RVB, but is instead a quasi-classical LRO $120^{\circ}$ state, which is less stable (Jolicoeur and Le Guillou, 1989; Jolicoeur et al., 1990; Bernu et al., 1994; Kaneko et al., 2014; Li et al., 2015) than the Néel order on the square lattice, since the sublattice magnetization of the triangular lattice is significantly reduced compared to its classical value. Because of this reduced stability, inherent to the triangular lattice, upon perturbing the Hamiltonian one may expect to see a variety of new phases. There have been some historically important semi-classical spin-wave theory (SWT) and ED studies(Jolicoeur et al., 1990; Dagotto and Moreo, 1989; Deutscher and Everts, 1993; Hirsch and Tang, 1989; Chubukov and Jolicoeur, 1992) for the model. However, such studies did not cover the physics of the whole phase diagram and were not able to capture the detailed properties of the groundstates. Previously, we elucidatedSaadatmand et al. (2015) the complete phase diagram of the $J_{1}$-$J_{2}$ THM on three-leg finite- and infinite-length cylinders to understand the crossover of 1D and 2D physics in the model. Moreover, other precise numerical approaches (Manuel and Ceccatto, 1999; Mishmash et al., 2013; Kaneko et al., 2014; Zhu and White, 2015; Hu et al., 2015; Li et al., 2015; Iqbal et al., 2016; Wietek and Läuchli, 2016; Saadatmand and McCulloch, 2016) demonstrate the existence of a spin-liquid (SL) state that stabilizes in a region ranging from $J_{2}^{\text{low}}\approx 0.05$ (Mishmash et al., 2013) up to $J_{2}^{\text{high}}\approx 0.19$ (Manuel and Ceccatto, 1999). Some numerical studies discovered magnetic orders outside this approximate SL region (see for example Jolicoeur et al., 1990; Chubukov and Jolicoeur, 1992; Bernu et al., 1994; Farnell et al., 2014; Li et al., 2015; Wietek and Läuchli, 2016). However, we suggest the detailed properties of the magnetic groundstates are still unclear in comparison to the well-understood counterparts in the (semi-)classical THM Jolicoeur et al. (1990); Chubukov and Jolicoeur (1992) and the quantum model on the three-leg cylinders. In particular, for the finite-size lattices, the largest system sizes for which the magnetic ordering of the $J_{1}$-$J_{2}$ THM was thoroughly studied is an $18\times 18$ torus (Kaneko et al., 2014) and a $30\times 3$ cylinder (Saadatmand et al., 2015).

It is noteworthy that the $J_{1}$-$J_{2}$ THM can describe some low-temperature properties of quasi-2D organic lattices, such as $\kappa$-(BEDT-TTF)2Cu2(CN)3 and Me3EtP[Pd(dmit)2]2, and inorganic materials, such as Cs2CuCl4, Cs2CuBr4, and RbFe(MoO4)2 (see Powell and McKenzie, 2011; Svistov et al., 2003; Balents, 2010 for details).

In this work, we establish the phase diagram for the $J_{1}$-$J_{2}$ THM on infinite-length cylinders with width up to 12 sites, and show that the model contains an LRO coplanar three-sublattice $120^{\circ}$ order ($J_{2}\leq 0.105(5)$), a four-fold degenerate toric-code-type $\text{Z}_{2}$-gauge spin liquid ($0.105(5)\leq J_{2}\leq 0.140(5)$; see also Saadatmand and McCulloch, 2016), an LRO collinear two-sublattice columnar order ($0.140(5)\leq J_{2}\leq 0.170(5)$ and $0.200(5)\leq J_{2}\leq 0.5$ for width-6 cylinders, and $0.140(5)\leq J_{2}\leq 0.5$ for larger cylinders), and an algebraic spin liquidWen (2002) (ASL) groundstate ($0.170(5)\leq J_{2}\leq 0.200(5)$, only for width-6 cylinders). The new tools for the cumulants and Binder ratios of the order parameter allowed us to locate the aforementioned phase transitions accurately, while the patterns of the TOS levels in the momentum-resolved ES revealed the structure of the magnetic order (or its absence). We also consider the explicit breaking of the time-reversal symmetry in the $J_{1}$-$J_{2}$ THM and the possibility of the stabilization of a chiral LRO. We note that there were recent, indecisive discussions Hu et al. (2016); Wietek and Läuchli (2016); Saadatmand and McCulloch (2016) on the robustness of the topological spin liquids against perturbing $H_{J_{2}}$ with a chiral term. Here, we confirm the non-chiral nature of such groundstates and the existence of a continuous phase transition toward the chiral spin liquid Wen et al. (1989); Kalmeyer and Laughlin (1987); Wen (2002) (CSL) phase through the study of a scalar chiral order parameter on width-8 cylinders.

The rest of the paper is organized as follows. In Sec. II, we provide the details of the employed $SU(2)$-symmetric MPS and DMRG methods (in particular, how we construct order parameter operators) and the geometry of the cylindrical lattices (in particular, the MPS mapping). In Sec. III, we present an overview of each of the phases, a schematic phase diagram for the model, and a more quantitative diagram in the form of short-range correlations. In Sec. IV, we directly measure the magnetization order parameters on some small-width ($L_{y}=3,4,5,6$) finite-length cylinders using MPS/fDMRG algorithms, to benchmark our calculations with another algorithm. Afterward, we focus on presenting our more precise iMPS/iDMRG results on infinite cylinders (having widths up to $L^{max}_{y}=12$). In Sec. V, we investigate the scaling behaviors of the correlation lengths against $m$, deep in each phase region. In Sec. VI, in order to better understand the entanglement entropy of the symmetry-broken LROs on cylinders, we study the entropy in the columnar magnetically-ordered phase. Details of our numerical tools are presented in Sec. VII and Sec. VIII, for cumulants and Binder ratios of the magnetization order parameters, and for TOS levels in the momentum-resolved ES, respectively. In Sec. IX, we test the robustness of a topological SL groundstate against chirality perturbations of the Hamiltonian to investigate the formation of long-range chiral orders, before some concluding remarks in Sec. X.

II Methods

To obtain the variational groundstate of the THM for a wide range of FM and AFM $J_{2}$-values in Eq. (1), we set $J_{1}=1$ as the unit of the energy, and employ the single-site DMRG algorithm (incorporating density-matrix mixingWhite (2005) with subspace expansionHubig et al. (2015) and $SU(2)$ symmetryMcCulloch and Gulácsi (2002); McCulloch (2007)). In addition, we construct the operators using the efficient formalism of matrix product operators(McCulloch, 2007, 2008; Schollwöck, 2011; Hubig et al., 2017) (MPOs), which represents an operator analogous to an ordinary MPS matrix. The MPO structure provides a formulation of any polynomial operator (with an expectation value that scales polynomially with the size of the lattice) in a Schur form (an upper- or lower-triangular matrix) for infinite systemsMcCulloch (2008); Michel and McCulloch (2010) (see below for an example and also Saadatmand, 2017 for an overview), which allows the calculation of the asymptotic limit of the expectation value per site. We keep up to $m=2,000$ number of states (approximately equivalent to 6,000 states of an Abelian $U(1)$-symmetric basis) in MPS/fDMRG, and up to $m=3,000$ number of states (approximately 9,000 $U(1)$-states) in iMPS/iDMRG calculations. Due to inherent 1D nature of the MPS, a mapping between the ansatz wavefunction and the triangular lattice is necessary. For the mapping purposes, we wrap the lattice in a way to create a long (or infinite-length) $L_{x}\!\times\!L_{y}$-site cylinder ($L_{x}$ can go to infinity; we also set $L=L_{x}\!\times\!L_{y}$) as in Fig. 1. We will employ a standard notation, previously presented in Saadatmand et al., 2015 (originally developed for single-wall carbon nano-tubesWilder et al. (1998)), to specify the wrapping vectors of the cylinders, ${\bf C}_{0}$, in terms of principal lattice directions using a notation of $(\hat{{\bf a}}_{+60^{\circ}},\hat{{\bf a}}_{-60^{\circ}})$. For the majority of the calculations, we choose the so-called YC wrapping, ${\bf C}_{0}[\text{YC}]=(L_{y},-L_{y})$ (we shall use the shorthand notations of YC$L_{y}$ and YC$L_{x}\!\times\!L_{y}$ to specify different YC lattice sizes). The YC structure is the only wrapping method with a circumference that equals to $L_{y}$ ($Y$-axis now coincides with the lattice short-direction and the $X$-axis coincides with the lattice long-direction) and is the best choice for the momentum-resolved ES (see below). However, in general, the choice of $L_{y}$ and ${\bf C}_{0}$ should respect sublattice ordering (if any) of the target state to avoid frustrating the groundstate. Consequently, depending on the desired width, the YC structure cannot be always used. Therefore, in finite-$L_{x}$ fDMRG calculations, we use a YC6 structure in all regions (allowing the stabilization of up to tripartite-symmetric groundstates), and YC3, ${\bf C}_{0}[L_{y}=4]=(4,-2)$, and ${\bf C}_{0}[L_{y}=5]=(5,-4)$ cylinders only in the $120^{\circ}$ and the SL phase regions. We also consider a YC4 structure in the columnar and the SL phase regions (occasionally, the YC3 and ${\bf C}_{0}[L_{y}=4]=(5,-4)$ systems are employed in the columnar phase region, however they are frustrating some forms of the collinear ordering – see below). For finite-length cylinders, we fix $L_{x}$ to a value that after which, an increase of the cylinder’s length would not change the average bond energy in the bulk of the system up to numerical uncertainties coming from the DMRG systematic errorsSchollwöck (2011); Saadatmand (2017). In $L_{x}\!=\!\infty$ iDMRG calculations, we use YC6 and YC12 structures in all regions, reserving YC9 only for the $120^{\circ}$ region, plus YC8 and YC10 in the columnar phase. We always set an efficient mapping for the infinite cylinders that minimizes the one-dimensional range of NN and NNN interactions, as shown in Fig. 1. Finally to calculate bipartite quantities, such as reduced density matrix, $\tilde{\rho}$, and entropy of the DMRG wavefunctionsSchollwöck (2011), we make a cut that does not cross any $Y$-direction bond and creates partitions $L$ and $R$, as shown in the figure.

We now present an overview of how to calculate higher moments of a (possibly non-local) observable in a translation-invariant infinite-size system. This is required for the measurements of the cumulants and Binder ratios of the magnetization order parameters (see below). For symmetry broken (or symmetry protected) states the Binder cumulant of the (string) order parameter can be evaluated directly in the thermodynamic limitKolley et al. (2015). However in this case, because we preserve $SU(2)$ symmetry the magnetic order parameter is strictly zero and the Binder cumulant is not well defined. However, as we show below, the moments can still be used to detect the signature of magnetic ordering. Suppose we are interested in calculating the matrix elements of the moments of an order parameter MPO, ${\bf M}^{[k]}$, of dimension $\tilde{m}$, that transforms under $SU(2)$ as a rank $k$ tensor. The explicit preservation of $SU(2)$ symmetry leads to the vanishing of the order parameter, but the even moments can be non-zero. In this case, the measurement of the expectation values of the higher-order magnetic moments, $\langle M^{n}\rangle$ (of order $n$), is of interest. These can be done using the method of the transfer operator Michel and McCulloch (2010). The generalized transfer (super-)operator, $\mathbb{T}_{X}$, associated with some operator of finite support (acting on a unit cell of an iMPS), $\hat{X}$, is defined as

[TABLE]

where $A^{s}$ are ordinary MPS matrices and $E^{a}$ denote the so called ${\bf E}$-matrices. In this context, $E^{a}$ is essentially an eigenmatrix, however, the ${\bf E}$-matrices are more familiar for their role in the expectation value of an MPO $\langle\mathcal{A}|\hat{M}|\mathcal{A}\rangle$ (see McCulloch, 2007, 2008 for full details). $E^{a}$ is in principle extensive, and on a $n$-site system, can be defined recursively as

[TABLE]

In the following example, for the sake of the simplicity we assume a one-site unit cell, although in practice for a magnetically ordered system the unit cell will be at least as large as the number of sublattices; The generalization for larger unit-cell sizes is straightforward. We give an example here for the second moment, the higher moments can be obtained recursivelyMichel and McCulloch (2010). To calculate the asymptotic limit of $\langle M^{2}\rangle$, one only needs to solve the diagrammatic fixed-point equation shown in Fig. 2(a), where the $E^{a}$-matrices are connected according to

[TABLE]

which can be solved sequentially, from $E^{1}$, $E^{2}$, $\ldots$, $E^{\tilde{m}}$. In practice, Fig. 2(a) shows the fixed point at which the addition of an extra site (or unit cell) to $E^{a}$-matrices will leave the system unchanged. The MPO form of the order parameter on a unit cell can be written(McCulloch, 2007, 2008) as a super-matrix (a matrix where elements are local operators acting on a single site or unit cell of the lattice):

[TABLE]

and we can attach $SU(2)$ quantum numbers $S=0,1$ to the rows/columns. The operator for the second moment has the form

[TABLE]

where in the last step, we have collapsed the middle rows to create a $3\times 3$ matrix with new quantum numbers of $S=0,1,0$ labeling the rows (assuming here that we want to calculate only the scalar part of $X^{2}$ – for the calculation of higher moments we need the other spin projections too). We can now write the fixed point of the last MPO in the form of recursive equations for the $E^{a}$-matrices, as shown in Fig. 2(b),(c), and (d). We note that the objects that appear on the right hand sides of the figures are nothing other than the generalized transfer operators. Translating the graphical notation into equations, for example, Fig. 2(b) can be written as

[TABLE]

which means $E^{1}(L)$ is an eigenmatrix of the transfer operator, which, for a properly orthogonalized MPS is just the identity matrix, so $E^{1}(L)=I$. As a result, equations for Fig. 2(c) and (d) can be written as

[TABLE]

where $C_{X}=\mathbb{T}_{X}(I)$ is a constant matrix, and

[TABLE]

The desired expectation value is encoded in the final matrix, i.e. $\langle M^{2}\rangle=Tr(E^{3}\tilde{\rho})$. However, importantly, most of the matrix elements of $E^{3}$ do not contribute to the expectation value of the second moment per site, we need only the component of $E^{3}$ that has non-zero overlap with $\tilde{\rho}$. Note that $\tilde{\rho}$ is the right eigenmatrix of $\mathbb{T}_{I}$ with the unity eigenvalue, hence the component of $E^{3}$ that gives the expectation value is the component in the direction of the corresponding left eigenmatrix of $\mathbb{T}_{I}$.

The calculation of the matrix elements of $E^{2}$ can be done efficiently using a linear solver. To see how this works, consider the eigenmatrix expansion for the transfer operator, $\mathbb{T}_{I}=\sum_{n=1}^{m^{2}}\eta_{n}|\eta_{n})(\eta_{n}|$, to obtain the eigenvalues $\eta_{n}$ and eigenvectors $|\eta_{n})$. If we write $C_{X}$ and $E^{2}$ matrices in this $\{|\eta_{n})\}$ basis with expansion coefficients $c^{(2)}_{n}$ and $e^{(2)}_{n}(L)$,

[TABLE]

then Eq. (8) is, for each component,

[TABLE]

Following Michel and McCulloch, 2010, we further decompose the coefficients into a component parallel and components perpendicular to the identity matrix, $I$ (i.e. the left eigenmatrix of $\mathbb{T}_{I}$, which has the largest eigenvalue of $\eta_{1}=1$ due to the MPS orthogonalization condition), and define

[TABLE]

so that $C_{X}=\tilde{C}_{X}+c^{(2)}_{1}I$ and $E^{2}=\tilde{E}^{2}+e^{(2)}_{1}I$. The reason for this is that the component in the direction of the identity $e^{(2)}_{1}$ diverges in the summation, whereas the other components that are perpendicular to the identity do not. Hence, we need to find the fixed points of these parts separately.

Solving Eq. (11) for the parallel components reveal the local expectation value of $X$ per site, which is an straightforward calculation,

[TABLE]

where $c^{(2)}_{1}$ is just the expectation value of the order parameter on one site. Hence $e_{1}^{(2)}(L+1)=\sum_{i=1}^{L}\langle X_{i}\rangle$, which is zero because of the $SU(2)$ symmetry (indeed, $c^{(2)}_{1}=0$ by construction, since it is in the wrong quantum number sector for the identity eigenvector of the transfer operator). The perpendicular components lead to

[TABLE]

where now $n\geq 2$, and the eigenvalues $|\eta_{n}|<1$. Thus, Eq. (14) is of the form of the sum of a convergent geometric series. Upon taking the limit $L\rightarrow\infty$ and writing back the projection operators as the original matrices, Eq. (14) converges to a fixed-point:

[TABLE]

which is a rather simple system of linear equations, and is numerically stable because the condition number of $1-\mathbb{T}_{I}$ is simply related to the leading correlation length, $1/(1-|\eta_{2}|)\simeq\xi$. In practice, generalized minimal residual method (GMRES) is a good choice of linear solver for Eq. (15). Upon obtaining the matrix elements of $\tilde{E}^{2}$, we can proceed to calculate the final expectation value as shown in Fig. 2(e). Note that this does not require all of the matrix elements of $E^{3}$, since we only require the overlap between $E^{3}$ and the density matrix (the right eigenvector of $\mathbb{T}_{I}$ with eigenvalue 1). This means that $\langle M^{2}\rangle=L\times\text{Tr}(\tilde{\rho}\mathbb{T}_{X^{2}}(I)+2\tilde{\rho}\mathbb{T}_{X}(E^{2}(L)))$, which is demonstrated in the MPS diagrammatic equation of Fig. 2(e). I.e. the only unknown is the $E^{2}$-matrix. This is a useful optimization and rather general – in calculating the expectation value of a triangular MPO of dimension $\tilde{m}$, only the matrix elements up to $E^{\tilde{m}-1}$ are required.

For calculating the 4th moment of a magnetization order parameter using $SU(2)$ symmetry, $X^{4}$ decomposes as

[TABLE]

where the dot product ${\bf X}\cdot{\bf X}=-\sqrt{3}[{\bf X}\times{\bf X}]^{[0]}$ and outer product ${\bf X}\otimes{\bf X}=\sqrt{6/5}[{\bf X}\times{\bf X}]^{[2]}$ are proportional to the $S\!=\!0$ and $S\!=\!2$ projections of the operator product, respectively, with an additional factor arising from the $SU(2)$ coupling coefficients. In general, we would need to also include the cross-product term $({\bf X}\times{\bf X})\cdot({\bf X}\times{\bf X})$ (proportional to the spin-1 projection), however, this vanishes due to antisymmetry under time reversal.

III Overview of the Phase Diagram

In this section, we present our findings for the phase boundaries and properties of Eq. (1), for different $J_{2}/J_{1}$ with $J_{1}>0$, using iDMRG and some benchmark comparisons using fDMRG.

In Fig. 3, we show the summary of the phase diagram, with four distinct phases; two phases with symmetry-broken magnetic order, a $Z_{2}$ spin liquid, and (only for the $YC6$ geometry) an algebraic spin liquid. For each of these phases, we present below visualizations of the correlation functions, obtained from well-converged iDMRG groundstates. In these visualizations, we depict spin-spin correlations, $\langle{\bf S}_{i}\cdot{\bf S}_{0}\rangle$, with respect to a reference site ${\bf S}_{0}$, using the size and the color of some spheres, and the NN correlators are depicted using the thickness and the color of some bonds. The reference site is denoted with the gray sphere. We also present the SSF up to the second Brillouin zone. Using the discrete Fourier transform of the real-space correlations to switch to the momentum space, one can write

[TABLE]

where ${\bf r}_{i}$ denotes the position vector of a spin ${\bf S}_{i}$ in the planar map of the lattice. The momentum vector, ${\bf k}$, will sweep the extended Brillouin zones. When the momentum vector coincides with the lattice’s wave vector, ${\bf Q}$, the occurrence of the condition $\lim_{N\rightarrow\infty}\frac{\text{SSF}(N)}{N}=\text{Const.}$ guarantees the existence of a true LRO. Plotting the SSF in the $(k_{x},k_{y})$ plane will reveal occurrence of strong FM correlations as Bragg peaks. However, for a fixed-$L_{y}$ infinite cylinder, one can only estimate the sums appearing in Eq. (17) using a finite length correlation. Therefore, we consider a large enough cutoff as an upper limit for $i$, namely $N_{c}$. We note that it is possible to obtain $SSF({\bf k},N\!=\!\infty)$ directly using the same method as described above for the moments (see also Michel and McCulloch, 2010), however, this is an expensive process and for calculating the entire ${\bf k}$ space it is much faster to calculate the real-space correlations and perform a Fourier transform. Here, we truncate the real-space correlation at the first point where $|\langle{\bf S}_{0}\cdot{\bf S}_{N_{c}}\rangle|\leq 10^{-5}$ is met for the nonmagnetic short-range correlated states (i.e. spin liquids) and the condition $|\langle{\bf S}_{0}\cdot{\bf S}_{N_{c}}\rangle|\leq 10^{-3}$ is met for the symmetry-broken quasi-LROs (i.e. $120^{\circ}$ and columnar states). The obtained phases are:

$J_{2}\rightarrow-\infty$: In this limit, one can readily show that the lattice decouples into three sublattices, each of which is a NN triangular lattice with bond strength $J_{2}$. In the case of vanishing interactions between sublattices ($J_{2}/J_{1}\rightarrow-\infty$), the groundstate for each sublattice is trivially a fully-saturated ferromagnet (see also Saadatmand et al., 2015) with total spin magnetization of $S^{\text{total}}_{A,B,C}=\frac{L_{u}}{2}$ per unit cell of each sublattice ($A$, $B$, or $C$). For a width-$L_{y}$ infinite-length YC structure, $L_{u}=L_{y}/3$ and $S^{\text{total}}_{A,B,C}=L_{y}/6$. The overall state can be any arbitrary mixture of three ${\bf S}^{\text{total}}_{A}$, ${\bf S}^{\text{total}}_{B}$, and ${\bf S}^{\text{total}}_{C}$ spin vectors, where they only have to follow the angular momentum summation rules. This will cause a large degeneracy for the overall groundstate, supporting total magnetization in a range of $0\leq S^{\text{total}}\leq\frac{3L_{u}}{2}$. Perturbing the Hamiltonian with a positive $J_{1}$ would then break this degeneracy and impose a $120^{\circ}$-ordered groundstate. Similarly to the case of three-leg cylindersSaadatmand et al. (2015), we find no signs of a phase transition for any $J_{2}<0$.

$J_{2}\leq 0.105(5)$: The groundstate is a coplanar quasi-classical $120^{\circ}$ order. Our investigations on infinite YC6, YC9, and YC12 structures find a three-sublattice magnetic ordered state exhibiting SSB in the thermodynamic limit (cf. Sec. VII and Sec. VIII). By imposing $SU(2)$ symmetry, the low-lying Nambu-Goldstone modes are evident and viewing the infinite cylinder as a 1D system it appears as a 1D quantum critical gapless state (cf. Sec. V). In Fig. 4(a), we present the correlation function for a YC12 groundstate at $J_{2}=-1.0$. The appearance of $\frac{L_{y}}{3}=4$ blue (ferromagnet) spheres per ring exhibiting a roughly constant size (for short distances) and all-AFM (red) bonds (throughout the cylinder) are characteristics of the phase. In Fig. 4(b), we present the SSF for a YC6 groundstate, deep in the $120^{\circ}$ phase. The formation of six strong Bragg peaks on a slightly-distorted regular hexagon is another characteristic for the phase. Using this data, we predict a wave vector of $Q_{120^{\circ}}\approx(\pm 3.64,\pm 2.09)$, which is very close to our expected theoretical value of $Q_{120^{\circ}}^{\text{theory}}=(\pm\frac{2\pi}{\sqrt{3}},\pm\frac{2\pi}{3})\approx(\pm 3.63,\pm 2.09)$ for a $120^{\circ}$ product stateJolicoeur et al. (1990); Chubukov and Jolicoeur (1992); Saadatmand et al. (2015); Sachdev (2011). We note that the correlation functions of YC6 and YC9, and SSFs of YC9 and YC12 structures in the $120^{\circ}$ phase are essentially identical to the results of Fig. 4.

$0.105(5)\leq J_{2}\leq 0.140(5)$: The groundstate is a four-fold degenerate toric-code-type $\text{Z}_{2}$ topological spin liquid (denoted by YC$L_{y}$-$\hat{\mathrm{a}}$ for the anyonic sector $\hat{\mathrm{a}}$ \in\{\mbox{\hat{\mathrm{i}}},\mbox{\hat{\mathrm{b}}},\mbox{\hat{\mathrm{f}}},\mbox{\hat{\mathrm{v}}}\}; see Saadatmand and McCulloch, 2016 for full details). In Fig. 5(a), we present the correlation function for a YC12-$\hat{\mathrm{i}}$ groundstate at $J_{2}=0.125$. The appearance of spheres with rapidly decaying radii and relatively weak all-AFM (red) bonds throughout the cylinder, is characteristic of the SL states (such a behavior of the correlations is also observed for other topological sectors and system sizes, except there exist some weak bond anisotropiesSaadatmand and McCulloch (2016); Hu et al. (2015)). In Fig. 5(b), we present the SSF for a YC10-$\hat{\mathrm{b}}$ groundstate at $J_{2}=0.125$. The spectral function is almost homogeneous, although being noisy and containing some weak diffusive peaks (compared to the strong Bragg peaks of magnetically-ordered states) reminiscent of gradual disappearance of the $120^{\circ}$ order. We notice that this overall pattern is virtually the same for all anyonic sectors and system sizes. Furthermore, our qualitative studies demonstrate that the homogeneity of the SSF is growing with increasing $L_{y}$ (not shown in the figures). For the topological SL phase, we find the lower and upper phase boundaries of $J_{2}^{\text{low}}=0.105(5)$ and $J_{2}^{\text{high}}=0.140(5)$, respectively. Using fDMRG for rather small YC6 widths (see below) we obtain similar results, $0.101(4)\leq J_{2}\leq 0.136(4)$. These phase boundaries are fairly close, but not identical, to those found by other authors (Manuel and Ceccatto, 1999; Mishmash et al., 2013; Kaneko et al., 2014; Zhu and White, 2015; Hu et al., 2015; Li et al., 2015; Iqbal et al., 2016; Wietek and Läuchli, 2016).

$0.140(5)\leq J_{2}\leq 0.5$, but excluding a region only for YC6 of $0.170(5)\leq J_{2}\leq 0.200(5)$: The groundstate is a quasi-classical collinear columnar (striped) order. Our investigations on infinite YC6, YC8, and YC10 structures (plus few more $J_{2}$ points on YC12 structures) show that that the columnar order is two-sublattice AFM state exhibiting SSB in the thermodynamic limit (cf. Sec. VII and Sec. VIII). Again, with $SU(2)$ symmetry the state appears on an infinite cylinder as 1D quantum critical. The correlation function for a YC12 groundstate at $J_{2}=0.5$ is presented in Fig. 6(a), where the appearance of robust FM stripes in the ${\bf a}_{+60^{\circ}}$-direction is clearly recognizable. In fact, the columnar order on the triangular lattice has three possible arrangements (Lecheminant et al., 1995) of FM stripes, each aligning with one of the three principal lattice directions, which are only degenerate in the thermodynamic limit. For the THM on three-leg (trivially) and four-leg cylinders (both finite and infinite-length cases), we found that the columnar order always has FM stripes in the lattice short (Y) direction, while for wider-width finite-length YC structures, FM stripes will be in either of ${\bf a}_{+60^{\circ}}$ or ${\bf a}_{-60^{\circ}}$-directions, producing only two degenerate groundstates. We numerically confirmed that, upon choosing a suitable wavefunction unit cell, iDMRG states randomly converge to one of these two states. We present the SSF for a YC8 groundstate at $J_{2}=0.5$ (with ${\bf a}_{+60^{\circ}}$-direction FM stripes), in Fig. 6(b). The formation of four, comparatively very strong Bragg peaks on a slightly-distorted regular parallelogram (with $60^{\circ}$ angles) is a characteristic of the phase. A wave vector of $Q_{\text{striped}}\approx\pm(1.82,3.18)$ can be estimated for the SSF, which is close to our expected theoretical value of $Q_{\text{striped}}^{\text{theory}}=\pm(\frac{\pi}{\sqrt{3}},\pi)\approx\pm(1.81,3.14)$ for a columnar product stateJolicoeur et al. (1990); Chubukov and Jolicoeur (1992); Saadatmand et al. (2015). We note that the SSFs of the columnar orders on YC6, YC10, and Y12 systems are rather similar to this result, however, the wave vector changes to $Q_{\text{striped}}^{\text{theory}}=\pm(\pi,-\frac{\pi}{\sqrt{3}})$, when the direction of FM stripes are switched. Our numerical calculations extend only to $J_{2}=0.5$. However we expect that there will be some additional geometry-dependent magnetically ordered phases for larger $J_{2}$ before reaching the large $J_{2}$ limit (see below).

$0.170(5)\leq J_{2}\leq 0.200(5)$, only for YC6: The YC6 geometry appears special in that we find signatures of an algebraic spin liquid, rather distinct from any other phase that we have observed in the model. Our results suggest that this phase is a quantum critical, gapless state with power-law scaling of the correlation function (cf. Sec. V) and no magnetic order. In Fig. 7, the presented correlation function and SSF appear to be reminiscent of a columnar-like ordering, but there are subtle differences. The size of spheres, representing the two-point correlation function, Fig. 7(a), decays faster than the columnar phase. In addition, the size of SSF peaks, Fig. 7(b), are considerably smaller than the typical size of the Bragg peaks in the columnar order having the same system width. In Sec. VII and VIII, below, we show that this phase has no signatures of magnetic ordering, which indicates that there are no broken symmetries and hence some kind of algebraic spin liquid. 6. 6.

$J_{2}\rightarrow+\infty$: Following the arguments presented for the $J_{2}\rightarrow-\infty$ case, in the limit of $|J_{2}|\gg 1$, the physical lattice will transform to three decoupled sublattices with antiferromagnetic NN bonds of the strength $J_{2}$. For $J_{2}\rightarrow+\infty$, the groundstate on each new sublattice is the same as the overall groundstate for $J_{2}=0$, i.e., the $120^{\circ}$ order. However for few-leg ladder systems, other symmetry broken phases could appear due to the restricted geometry. As an example, for three-leg finite cylindersSaadatmand et al. (2015) in $J_{2}\rightarrow\infty$, we found that the groundstate is three weakly-coupled copies of a NNN Majumdar-Ghosh state. Interestingly, we found a similar dual Majumdar-Ghosh phase for four-leg finite cylindersSaadatmand (2017)). Consistent with the expected 2D limit, we did not observe any signature of such Majumdar-Ghosh-type phases for $L_{y}>4$ ladders. In addition, semi-classical SWT studies (Jolicoeur et al., 1990) confirms that the “order from disorder” mechanism would choose three-fold degenerate and decoupled states for $J_{2}\gg 1$, which are energetically favorable to arrange according to $120^{\circ}$ ordering. Hence, we expect that such exotic ordered phases are particular features of narrow cylinders.

To get a better quantitative insight on the phase diagram of the THM, we study the short-range (NN and NNN) spin-spin correlations, $\langle{\bf S}_{i}\cdot{\bf S}_{j}\rangle$, Fig. 8. Short-distance correlators in a crystalline phase have a repetitive pattern reflecting the bulk properties of the groundstate. In Fig. 8, we choose six reference bonds, including three NN and three NNN correlators, to build up a picture of the real-space correlations for different system widths. In the $120^{\circ}$ phase region, correlators are very nearly isotropic, where NN (NNN) bonds are all AFM (FM). On the other hand, topological spin liquids on finite-width systems contain strong anisotropiesHu et al. (2015); Zhu and White (2015), which is clearly seen in Fig. 8. As we showed previouslySaadatmand and McCulloch (2016), in the thermodynamic limit anyonic sectors $\hat{\mathrm{b}}$ and $\hat{\mathrm{f}}$ are anisotropic on finite cylinders, while $\hat{\mathrm{i}}$ and $\hat{\mathrm{v}}$ are isotropicSlagle and Xu (2014); Lu (2016). The behavior of the correlation functions is distinct in the columnar phase, where there are always two FM bonds (one is a NN and another one a NNN correlator) and four AFM bonds (two are NN and other two NNN correlators) out of the six reference bonds. The FM stripes of the columnar order can, of course, choose either of ${\bf a}_{+60^{\circ}}$ or ${\bf a}_{-60^{\circ}}$ directions, so such data-points in this region are exchangeable. Furthermore, curiously for YC6, in the ASL phase region ($0.170(5)\leq J_{2}\leq 0.200(5)$), the system temporarily restore all symmetries, again, by crossing the ${\bf a}_{\pm 60^{\circ}}$-direction bonds.

IV Direct measurement of the order parameters on finite-length cylinders

To provide a verification of the phase boundaries for comparison against our iDMRG results, we calculated two magnetization order parameters on $L_{y}\leq 6$ finite-length cylinders (small compare to the largest width of our infinite-length YC systems) using an approach originally suggested by White and ChernyshevWhite and Chernyshev (2007). Consider the arbitrary magnetization vector order parameter of ${\bf M}(m)$ for a wavefunction with $m$ number of states (the preservation of the SU(2)-symmetry causes the structural vanishing of all projection components). Upon a suitable choice of the system size and careful extrapolation toward the thermodynamic limit, non-zero values for the second moment of ${\bf M}$ (which is directly proportional to the spin susceptibility) can be derived. In White and Chernyshev’s method, one first extrapolates the order parameter linearly with the DMRG truncation errors, $\varepsilon_{m}$, toward the thermodynamic limit of $m\rightarrow\infty$ ($\varepsilon_{m}\rightarrow 0$) to calculate $\langle{\bf M}^{2}(\infty)\rangle$. Then, using only fixed aspect-ratio ($\frac{L_{y}}{L_{x}}=\text{Const.}$) system sizes, $L_{x}$ and $L_{y}$ should be simultaneously extrapolated toward the thermodynamic limit of $L\rightarrow\infty$. By employing a similar approach, plus some simple dimensional analyses and numerical examination of the magnetic moments, we suggest in the MPS constructions of the $SU(2)$ $S\!=\!0$-sector groundstates on fixed aspect-ratio cylinders ($L_{x}>L_{y}$), the normalized order parameter, ${\bf M}^{2}(\infty)$ per site, scales as

[TABLE]

where eclipses represents higher order terms in $\frac{1}{L_{x}}$ (note that Eq. (18) is only a heuristic fit; see Neuberger and Ziman, 1989 for theoretical predictions). One should note that any independent growth of $L_{x}$ and $L_{y}$ toward the $L\rightarrow\infty$ limit, can be interpreted as the existence of an infinitely long cylinder at some stage. This will collapse the system, essentially, to an inherently 1D state, for which the behavior of the magnetic moments is essentially different (see below).

The magnetic order parameters that we selected to study the phase diagram on finite-size cylinders include: the FM sublattice magnetization, defined arbitrarily on sublattice $A$,

[TABLE]

where ${\bf S}_{A}=\sum_{i\in A}{\bf S}_{i}$ is summed over all sites in sublattice $A$, and $2/\sqrt{N_{A}(N_{A}+2)}$ is a normalization factor ($N_{A}$ is the total number of sublattice-$A$ spins on the finite lattice). $O^{\text{FM}}_{A}$ is a well-defined order parameter for the $120^{\circ}$ phase. The classical $120^{\circ}$ order will result in the maximum possible value for the order parameter in the limit of $L\rightarrow\infty$, i.e. $O^{\text{FM}}_{A}[\text{classical},L\rightarrow\infty]=1$ Saadatmand (2017). The next order parameter is the staggered magnetization, ${\bf M}_{\text{stag}}$, for which the second moment is a well-defined order parameter for the columnar phase,

[TABLE]

where ${\bf S}_{\text{stag}}={\bf S}_{A}-{\bf S}_{B}$ is the staggered magnetization for sublattices $A$ and $B$. The classical columnar order will result in the maximum possible value for the order parameter in the limit of $L\rightarrow\infty$, i.e. $O_{\text{stag}}[\text{classical},L\rightarrow\infty]=1$Saadatmand (2017).

Our results for $O^{\text{FM}}$ and $O^{\text{stag}}$, in the thermodynamic limit of $L\rightarrow\infty$, are presented in Fig. 9. Individual error-bars are relatively large, but the overall behavior of the magnetization curves follow the expected pattern: there exists a small region for $J_{2}$, where both $O^{\text{FM}}(L\rightarrow\infty)$ and $O^{\text{stag}}(L\rightarrow\infty)$ are touching the zero axis (considering uncertainties), which provides SL region boundaries, $J^{\text{low}}_{2}=0.101(4)$ and $J^{\text{high}}_{2}=0.136(4)$. Next to the SL phase region, on the left, only $O^{\text{stag}}(L\rightarrow\infty)$ is touching the zero axis, while $O^{\text{FM}}(L\rightarrow\infty)$ are increasing for $J_{2}\rightarrow-\infty$ (confirming the stabilization of $120^{\circ}$ order in this region). On the other hand, next to the SL phase region on the right, $O^{\text{stag}}(L\rightarrow\infty)$ increases rapidly, indicating columnar order. Interestingly, the value of $O^{\text{FM}}(L\rightarrow\infty)$ ($O^{\text{stag}}(L\rightarrow\infty)$) is increasing (decreasing) again for large $J_{2}$. This is consistent with the existence of a multi-component $120^{\circ}$ order (three copies of a conventional $120^{\circ}$ order placed on sublattices; see Sec. III) in the $J_{2}\rightarrow\infty$ limit.

It is worth noting the magnitude of the sublattice magnetization at $J_{2}=0$ (NN model). Measurement of variants of a $120^{\circ}$ order parameter for the NN model has been in the center of attention (Jolicoeur and Le Guillou, 1989; Bernu et al., 1994; Li et al., 2015; Kaneko et al., 2014; Saadatmand et al., 2015) to understand the degree of magnetization reduction (in comparison to their classical counterparts) in such a frustrated model. As shown in Fig. 9, we predict $O^{\text{FM}}[J_{2}=0]/O^{\text{FM}}[\text{classical}]=25(8)\%$, which is considerably smaller than approximate results of $50\%$ by SWT (Jolicoeur and Le Guillou, 1989), $48\%$ by ED (Bernu et al., 1994), $40\%$ by CCM (Li et al., 2015), and $50\%$ by variational QMC (Kaneko et al., 2014).

V Correlation Lengths

For infinite cylinders, the gapped or gapless nature of the groundstate can be understood through the study of the (principal) correlation length, $\xi$, since the behavior of the magnetic ordering and the scaling behavior of the static correlation functions are connected (cf. Sec. I). Indeed, the Hastings-Oshikawa-Lieb-Schultz-Mattis theoremHastings (2004a, b) relates the size of the energy gap, $\Delta_{e}$, to $\xi$, for local, translation-invariant Hamiltonians on even-width cylinders as $\xi\leq\frac{\mathrm{Const.}}{\Delta_{e}}$ (i.e. $\xi^{-1}$ serves as an upper boundary for the gap size). For the (inherently one-dimensional) MPS ansatz, the connection between entanglement scaling and the correlation length is well-understood (Tagliacozzo et al., 2008; Stojevic et al., 2015; McCulloch, 2008). In a critical phase, the correlation length diverges with a signature power-law scaling with the number of states as $\xi(m)=\tilde{\kappa}_{c}m^{\tilde{\kappa}}$. Furthermore, in such states, the entanglement entropy diverges with a scaling of $S_{EE}\sim\log\xi$. On the other hand, for short-range gapped states $\xi$ saturates to a finite value as $m$ is increased, which in the topological spin-liquid state of the THM is short; of the order of a few lattice spacings. Interestingly, as we see below, the correlation length scaling for magnetic ordering in $SU(2)$-symmetric MPS on infinite cylinders appears rather differently than the full 2D limit. Such cylindrical magnets exhibit some signatures of true LRO (e.g. in the ES – see below), however, due to the explicit preservation of SU(2) and the dominating 1D physics of the MPS ansatz, the groundstates emerge as quasi-LRO critical states (note that the correlation length can still diverge with respect to the cylinder circumference). Nevertheless, in the iMPS representation of the wavefunction, the correlation lengths (per unit-cell) can be conveniently read from the eigenspectrum of the transfer operator, $\mathbb{T}_{I}$ (cf. Sec. II):

[TABLE]

where $\eta_{i}$ are eigenvalues of $\mathbb{T}_{I}$ (arranged as $\{|\eta_{1}|>|\eta_{2}|>|\eta_{3}|>...\}$). $\eta_{i}$ depends on the number of states, and are also labeled by an $SU(2)$ spin sector, which is the symmetry sector of the (block diagonal) transfer operator, and corresponds to the symmetry of the associated correlation function. We have discarded $i=1$, as the largest eigenvalue of $\mathbb{T}_{I}$ in an orthonormalized basis always corresponds to $\eta_{1}=1$ (belonging to the identity eigenmatrix) and the principal correlation length is the second largest eigenvalue, $\xi_{2,S}\equiv\xi_{S}$. For a phase with magnetic ordering, such as $120^{\circ}$ and columnar order, the principal correlation length is expected to belong to the $S\!=\!1$ sector, indicating that the slowest decaying correlations are in the spin-spin form. For the topological and algebraic spin liquid phases, we find that the principal correlation length is in the $S\!=\!0$ sector, indicating that the slowest decaying correlation is some kind of singlet-singlet correlator (we have not determined the exact form). An undesirable effect of the variational convergence of the groundstate using the iDMRG approach emerges from the constraint of $SU(2)$ symmetry, whereby spurious symmetry effects make the wavefunction non-injective (the spectrum of $\mathbb{T}_{I}$ contains multiple identity eigenvalues in each $S$-sector). We have removed such wavefunctions everywhere except in the immediate vicinity of the $J_{2}=0.105(5)$ transition (expectation values are still reliable), where the non-injectivity was difficult to avoid (this is likely due to the closeness of this point to the topological SL region).

We present the correlation length results for the ordered phases in Fig. 10, where we compared them against the $\xi_{S}$ from the topological spin liquidSaadatmand and McCulloch (2016). We immediately notice that the principal correlation length belongs to the $S=1$ sector for the magnetic groundstates with $120^{\circ}$ and columnar ordering, however, it switches to the $S=0$ sector for all SL states, whether they are quasi-LROs (as in ASLs) or short-range correlated (as in topological spin liquids). We can see that both the ASL and the magnetically ordered states have power-law behavior, reflecting their gapless and quantum critical natures. We emphasize that in the case of the ASL, this behavior appears to be intrinsic; however for the magnetically ordered phases the power-law correlations are a consequence of preserving $SU(2)$ symmetry. In contrast, for the topological spin liquid, the correlation length is considerably smaller in the size (order of few lattice spacings) and qualitatively begins to saturate in the large-$m$ limit, although it is surprisingly difficult to do a rigorous fit.

VI Entanglement entropy of quasi-LRO magnets

The entanglement entropy is a central quantity in the physics of the many-body systems, which provides a measure of how strongly conjunct subsystems are entangled. The entropy has proven to be a powerful numerical tool for characterizing the low-energy spectrum, detection of SSB, and topological degeneracy of the groundstate (for some examples, see Kitaev and Preskill, 2006; Levin and Wen, 2006; Eisert et al., 2010; Zhang et al., 2011; Wang et al., 2013; Metlitski and Grover, 2015). Between many different approaches to measure entropy, we employ the method of Jiang *et al.*Jiang et al. (2012) that calculates the von Neumann entropy along a bipartition cut of the cylinder, as shown in Fig. 1, since it is computationally convenient to manipulate in the context of MPS and DMRG algorithms. The bipartite von Neumann entropy is defined as $S_{EE}=-Tr(\tilde{\rho}\log\tilde{\rho})$. In terms of the eigenvalues of $\tilde{\rho}$, i.e. $\{\lambda_{i}\}$, the entropy can be written as $S_{EE}=-\sum_{i}\lambda_{i}\log\lambda_{i}$. Roughly speaking, $S_{EE}$ counts the number of entangled pairs on the bipartite boundary. $S_{EE}$ is a function of the $(D-1)$-dimensional area of the $D$-dimensional quantum system, i.e. the boundary size, $L_{cut}$ (note that $L_{cut}=L_{y}$ for the YC structure). In fact, robust theoretical studies (Hamma et al., 2005a, b; Kitaev and Preskill, 2006; Levin and Wen, 2006) proved that for interacting 2D spin systems with only local couplings and a cut size significantly larger than the correlation length, the leading term in the entropy scales with the boundary area, $S_{EE}\propto L_{cut}$, not the system volume, which is known as the area-law (the area-law was originally introduced in the context of the black holes (Bekenstein, 1973) and quantum field theory (Srednicki, 1993; Plenio et al., 2005)). However, for strictly 1D quantum critical states (in the thermodynamic limit) the condition of the boundary size being considerably larger than the correlation length cannot be met, and the $S_{EE}$ behavior is modified. In this case, the leading term in the entropy relates to the only length scale of the system, i.e. the correlation length, as $S_{EE}\sim\log(L_{eff})\sim\log(\xi)$ (Tagliacozzo et al., 2008; Stojevic et al., 2015), where $L_{eff}$ stands for the effective size of the system. For the symmetry-broken true LROs, again, the size of the cut is significantly smaller than the diverging correlation length and a logarithmic term should be added to the area-law behavior (Metlitski and Grover, 2015):

[TABLE]

where $\beta_{0}$ corresponds to a non-universal constant, which depends on the system geometry, the topological entanglement entropyKitaev and Preskill (2006); Levin and Wen (2006), spin stiffness, and the number, $N_{G}$, and the velocity of the Nambu-Goldstone excitations. In addition, $\beta_{1}$ is another non-universal constant, which depends on the short-range entanglement in the vicinity of the cut and a short-distance characteristic cutoff. For the (quasi-)LRO, $SU(2)$-symmetric, iMPS groundstates on the infinite cylinders, we find that the entropy scaling behavior is distinct. As discussed in Sec. V, MPS-ansatz symmetry broken magnets appear as quantum critical states on the cylinder. Thus, it is expected that the entropy exhibits a combination of the area-law and the critical behaviors. Our numerical measurements on an $SU(2)$-symmetric, quasi-LRO groundstate of the $J_{1}$-$J_{2}$ THM on the infinite YC structures confirms such a mixed scaling as of

[TABLE]

where

[TABLE]

The behavior of the non-universal constant of $a_{0}$ proved to be more challenging to predict, but it can only contain sub-leading corrections to the area-law term appearing in $a_{1}$ (see below).

In Fig. 11, we present our entropy measurements for the groundstates of the THM deep in the columnar phase region. Due to exponential cost of the calculations with the system width we only obtained a few wavefunctions for different $L_{y}$ in the columnar phase. However the results shown in Fig. 11 confirm the prediction of Eq. (23) and Eq. (24). In the figure, we first fit a line to the original entropy data and calculate $a_{0}$ and $a_{1}$ for each system size. Clearly, $a_{1}$-values are consistent with the area-law behavior. We measured the coefficients of $a_{1}$ as $\alpha_{0}=-0.28(1)$ and $\alpha_{1}=0.068(1)$. In contrast, there was no obvious fit possible for $a_{0}$-values, but their saturating nature for the large-$L_{y}$ limit is consistent with this term being a sub-leading correction to the mixed term containing the area-law behavior.

VII Numerical tools I: cumulants and Binder ratios of the

magnetization order parameters

In Sec. II, we constructed the theoretical framework for a method to measure the non-local moments and cumulants of the magnetic order parameters, in the context of $SU(2)$-symmetric translation-invariant MPS, where all projection components of the magnetic order parameter, ${\bf M}^{[k]}$, vanish by construction. In this case, the higher moments can play the role of the order parameter. It is convenient to connect the moments of the operators to the $i$th cumulant per site, $\kappa_{i}$, by employing

[TABLE]

where $B_{n,i}$ are partial Bell polynomials (Stuart et al., 2009) and $L$ now stands for the operator length. For some examples, we expand Eq. (25) to write the relations for the first few cumulants,

[TABLE]

The cumulants per site are obtained directly as the asymptotic large $L$ limit obtained from the summation of the tensor diagrams presented in Sec. II. For the iMPS ansatz, when the asymptotic limit is taken to derive a translation-invariant infinite-size system, one should replace the operator length with the effective system size as $L\rightarrow L_{eff}\propto\xi$ (see also Sec. VI). Below, we introduced the magnetic order parameters that are used to measure the cumulants and characterize the LROs of the THM. We first construct the MPO forms of the higher moments of a staggered magnetization (the order parameter for columnar order on cylinders with FM stripes in ${\bf a}_{+60^{\circ}}$-direction),

[TABLE]

and a tripartite magnetization (the order parameter for the $120^{\circ}$ phase),

[TABLE]

on a $L_{y}$-size unit cell. Numerical computation of the moments of such order parameters is a challenging task due to relatively large dimensions of the resulting MPOs. Nevertheless, we succeeded to calculate the second cumulant, $\kappa_{2}$, and the fourth cumulants, $\kappa_{4}$, of $M_{stag}$ and $M_{tri}$ (the odd moments vanish due to the $SU(2)$ symmetry) for a range of the groundstates. We suggest that the most useful choice of cumulants is $\kappa_{4}$, which is connected to the excess kurtosis (Stuart et al., 2009), $\gamma_{4}$, of the block distribution function associated with the operator ${\bf M}^{[k]}$:

[TABLE]

We emphasize that the above equation is only valid for the $\kappa_{1}=\kappa_{3}=0$ case. The importance of the fourth cumulant was revealed by some studies on fourth magnetic moment behavior of 2D Ising antiferromagnets (Tsai and Salinas, 1998; Malakis et al., 2007), which established $\kappa_{4}$ as an effective tool for pinpointing quantum critical points. In these studies, the scaling behavior of the fourth magnetic moment is observed to vary significantly at an Ising transition (more precisely, $\kappa_{4}$ changes sign at the critical point, and changes by many orders of magnitude nearby the critical point). Another relevant and interesting (dimensionless) quantity is the Binder cumulant (Binder, 1981a, b; Binder and Landau, 1984; Sandvik, 2010), $U_{L}=\frac{n_{H}+2}{2}(1-\frac{n_{H}}{n_{H}+2}\frac{\langle{\bf M}^{4}\rangle}{\langle{\bf M}^{2}\rangle^{2}})$, where $n_{H}$ is the number of projection spin operators used to construct the order parameter (e.g. $n_{H}=3$ for a vector magnetization). In the vicinity of a critical point, the Binder cumulant becomes independent of the system size (lower moments of the order parameter cancel out higher-order finite-size effects) and can be used to pinpoint the transition. Previously, we adoptedSaadatmand et al. (2015) $U_{L}$ of a (scalar) dimer order parameter to locate a critical point in the phase diagram of the THM on three-leg cylinders. However, until now, the scaling behavior of $U_{L}$ was less-known for the cases where the order parameter itself is strictly zero. In the limit of $L\rightarrow\infty$, as it is clear from Eq. (26), the higher-order corrections in $\langle M^{n}\rangle$ vanish and the conventional method of Binder cumulants for locating the phase transitions becomes ineffective. However, the correlation length, $\xi$, gives us a natural length scale and a rather precise process to scale a Binder-cumulant-type quantity in the vicinity of a critical point. As in the case of the entropy, Sec. VI, the key to obtaining the correct scaling of the magnetic moments of iMPS wavefunctions is to choose $L_{eff}=\tilde{s}\xi$, where $\tilde{s}$ is any fixed scaling constant. For Binder cumulant, $\tilde{s}$ has no qualitative effect except to change the value of the critical binder cumulant, similar to the role of boundary conditions for the finite-size Binder cumulant. Therefore, one can freely choose $\tilde{s}$ to obtain the most numerically stable fit. When the order parameter is zero by symmetry, so that $\kappa_{1}=\kappa_{2}=0$, the appearance of such a constant is irrelevant and only the ratio of the second and fourth cumulants plays a role. By replacing the explicit relations for $\langle{\bf M}^{2}\rangle$ and $\langle{\bf M}^{4}\rangle$ from Eq. (25) into $U_{L}$, we propose the ratio (which we call the “Binder ratio” – see also Eq. (29)):

[TABLE]

We find that numerically this combination of the moments and the correlation length removes much of the numerical noise that appears in the individual moments.

We present the extrapolated results of $\kappa_{2}$ and $|\kappa_{4}|$ for ${\bf M}_{stag}$ and ${\bf M}_{tri}$, in the limit of $m\rightarrow\infty$, in Fig. 12. In the figures, each data-point is the result of a separate extrapolation of the cumulants versus $m$. Upon careful numerical examination of the scaling behaviors of numerous groundstates in the various phases, we were able to establish the scaling relation of $|\kappa_{n}|=\breve{a}_{0}+\breve{a}_{1}e^{-\breve{a}_{2}m},~{}n=2,4,$ for ordered phase regions and make sense of the cumulant results in the $m\rightarrow\infty$ limit. These results show that $\kappa_{4}$ is comparatively large and negative when there is quasi-long-range magnetic ordering. Moreover, $\kappa_{2}$ is large and positive for quasi-LROs (see $120^{\circ}$ and columnar phase regions in Fig. 12). This is in contrast to the behavior near phase transitions, and within the topological and algebraic spin liquids, where we were not able to find an appropriate analytical fit for the cumulants in the $m\rightarrow\infty$ limit, as they behave irregularly or quickly decay to numerically vanishing values. A likely reason for this is that for a magnetically-ordered, $SU(2)$ $S=0$ groundstate, the moments ${\bf M}^{[k]}$ acquire a set of equally-weighted non-zero values from the limited number of recovered (purely) TOS levels by iDMRG (see below). In such a case, the distribution function would resemble a discrete uniform distribution with very large and negative $\kappa_{4}$, and large and positive $\kappa_{2}$. However, for disordered states with no symmetry breaking in the thermodynamic limit, the distribution function is expected to resemble the normal distribution centered around zero magnetization, which has vanishing $\kappa_{4}$. For $\kappa_{2}(m\rightarrow\infty)$, in Fig. 12(a), we display in bold the boundaries where we were not able to extrapolate to $m\rightarrow\infty$. These are quite close to the phase transitions indicated by fDMRG, Fig. 9 (except for the YC6 structure, where we find an additional ASL phase), which supports the validity of the iDMRG cumulant method. The same behavior was observed for $\kappa_{4}(m\rightarrow\infty)$, indicated by the black dashed-lines in Fig. 12(b). In addition, the extremely large (negative) values of $\kappa_{4}(m\rightarrow\infty)$ are consistent with our interpretation.

Our attempts to pinpoint the phase transitions of the THM on infinite cylinders, using $U_{r}$, is presented in Fig. 13 and Fig. 14. Based on these results, we argue that $U_{r}(m)$, as the ratio between $\kappa_{4}(m)$ and $\kappa_{2}(m)$, scaled with $\xi(m)$, regularly extrapolates to a finite value in the $m\rightarrow\infty$ limit everywhere, except close to (or on) a phase transition, or when the wavefunctions are non-injective (cf. Sec. V). Careful numerical examination suggests that the Binder ratios scale with a saturating behavior similar to the cumulants,

[TABLE]

In Fig. 13, we observe that in the topological SL phase region, $U_{r}(m)$ has a comparatively small value, as expected for nonmagnetic phases with $\gamma_{4}\rightarrow 0$. In addition, when there is magnetic ordering, $U_{r}(m)$ converges to a finite, negative value, while it appears different $m$-curves tend to group together. The latter should be due to the fact that the iMPS magnetic orders are quantum critical states with an $U_{r}$ independent from $L_{eff}$. Furthermore, for the ASL phase of YC6 structures (see Fig. 14(a)), very close to the expected phase transition points from the short-range correlation data, Fig. 8, and within the entire ASL phase region, $U_{r}$ diverges with $m$ (e.g. see the inset of Fig. 14(a)), where it is impossible to extrapolate to a finite $U_{r}(\infty)$. In the immediate vicinity of the transition from the $120^{\circ}$ to topological spin liquid (cf. Fig. 14(b)), it was not possible to employ Eq. (31) due to unavoidable non-injectivity of the wavefunctions. However, we suggest that the fixed-$m$ results are rather reliable and can be used to estimate a phase transition. Overall, we locate critical points of the THM from the discontinuities of $U_{r}(m\rightarrow\infty)$-lines (i.e. where there is no extrapolation possible) or when there is a significant kink in fixed-$m$ data. Based on this approach, we estimate the phase transition points of $J_{2}=0.105(5)$ between the $120^{\circ}$ and topological spin liquid states using YC6 results of Fig. 14(b), $J_{2}=0.140(5)$ between topological SL and columnar states using YC10 results of Fig. 13(a) (YC8 results of Fig. 13(b) would estimate a transition very close to this point, so we have based the final prediction on the larger-width data), and transition points of $J_{2}=0.140(5),0.170(5),0.200(5)$ encapsulating the ASL and columnar states using YC6 results of Fig. 14(a). To further validate $U_{r}$ accuracy in estimating the transitions in case of YC6 structures, we also provide a numerical approximation for the fidelity susceptibility (You et al., 2007),

[TABLE]

in Fig. 15, where we set $\delta J_{2}=0.05$. The fidelity susceptibility is known to be well-behaved and small when away from a phase boundary, but can diverge at a transition. It is clear from the figure that the diverging peaks of $\chi_{F}^{\text{approx}}$ (considering their tendency to lean toward the right) are happening relatively close to the predicted phase transitions from the Binder ratio results of Fig. 14(a).

VIII Numerical tools II: ‘TOS columns’ in the momentum-resolved

entanglement spectrum

The entanglement between the partitions of a quantum system is encoded in the spectrum of the entanglement Hamiltonian, $H_{E}=-\log(\tilde{\rho})$, i.e. $\{-\log(\lambda_{i})\}$, which is known as the ES and commonly presented using energy-level arrangements analogous to an energy spectrum. $\{-\log(\lambda_{i})\}$ can be labeled using any global-symmetry quantum number to extract more information on the symmetry nature of the state (as long as the corresponding symmetry is preserved on the bipartite cut). $H_{E}$ maintains the symmetries of a cylindrical wavefunction, however, they may exist some symmetries that are not explicitly preserved by the ansatz due to the mapping of the 2D model onto an MPS chain. Nevertheless, one can still diagonalize such a symmetry operator in the ‘auxiliary’ basis (i.e. the basis that diagonalizes $H_{E}$) to create a new set of good quantum numbers (see Kolley et al., 2013; Saadatmand and McCulloch, 2016, 2016; Cincio and Vidal, 2013 for some examples). When the $SU(2)$-symmetry is preserved in the calculation, the obvious choice for the labels is the spin $S$ quantum number (belonging to a single partition of the system). We refer to an $H_{E}$ spectrum that is plotted against $S$ (where no other label exists) as the spin-resolved ES. Kolley et al.(Kolley et al., 2013) showed that the spin-resolved ES of a magnetically ordered state on finite-length cylinders shows signatures of symmetry-breaking at the thermodynamic limit. This emerges from a key finding: the realization (Metlitski and Grover, 2015; Kolley et al., 2013; Rademaker, 2015) that the low-energy part of the ES of magnetic orders exhibits a specific type of grouped levels, known as the entanglement-spectrum TOS (also referred to as “quasi-degenerate joint states”), closely resembling the low-lying levels in the energy spectrum, known as the Anderson TOS levels (Anderson, 1952; Lhuillier, 2005; Bernu et al., 1994) (also referred to as the “Pisa tower” structure or the “thin spectrum”), which is considered as clear-cut evidence for the existence of true LROs on finite lattices. Kolley et al. established that, similar to the energy spectrum, for a fixed $S$-sector, entanglement-spectrum TOS levels are well-separated from the denser rest of the spectrum and the lowest energy levels of the ES, immediately above the TOS levels, are spin-wave states (Nambu-Goldstone modes). In this paper, we are interested in exploiting both the $S$ quantum numbers ($SU(2)$ is explicitly preserved in the iDMRG calculations), and the momenta in the cylinder $Y$-direction, $k$, i.e. the complex phase of the eigenvalues of the reduced $T_{y}$ operator, where $T_{y}$ is the translation by one site in $Y$-direction; we can decompose the operator in the same way as the Schmidt decomposition of the wavefunction Schollwöck (2011), $T_{y}=T^{L}_{y}\otimes T^{R}_{y}$, where $T^{L}_{y}$ and $T^{R}_{y}$ are the reduced operators and maintain the unitary property of the original operator. $T_{y}$ is not preserved exactly in the calculations due to the MPS mapping on the cylinder, Fig. 1, but it can be diagonalized straightforwardlySaadatmand (2017). We refer to an $H_{E}$ spectrum that is plotted against $k$ and additionally labeled by $S$, as the momentum-resolved ES, \{-\log\big{(}\lambda_{n}[k_{n},S_{n}]\big{)}\}. For a system with PBC in $Y$-direction, dihedral symmetry implies that $T_{y}^{L_{y}}=I$. As a result, the allowed momentum spacing is as $\Delta k_{n}=\frac{2\pi n}{L_{y}}$ for $n=0,1,...,L_{y}-1$. We notice that $k_{0}$, the momentum of the lowest ES level, is not fixed due to the possibility of inserting a shift in the expectation value of $T_{y}$ (one needs to first fix $k_{0}$, then measure the rest of the momenta in respect to it; physically only $\Delta k_{n}$ matters here – see also Saadatmand, 2017; Michel and McCulloch, 2010). The study of momentum-resolved forms of the ES is now finding a place in the literature of the low-dimensional quantum magnets. Another key breakthrough was the realization of that such ES can be used to fully classify anyonic sectors of chiral (Cincio and Vidal, 2013) and $\mathbb{Z}_{2}$-gauge (Zaletel et al., 2015) topological orders on infinite cylinders (Saadatmand and McCulloch, 2016). Below, we argue that the symmetry-breaking can be recognized and characterized using the momentum-resolved ES, which shows the symmetry properties even more robustly than the spin-resolved ES.

Upon careful examination of the momentum-resolved ES of the magnetic orders in the THM on infinite cylinders and noticing the underlying symmetries of the sublattices, we find that the spectrum contains exactly $N_{s}$ (number of the groundstate sublattices) column-like structures, which are the low-lying component TOS levels, independent of the system width. We shall refer to these particular patterns as ‘TOS columns’. The appearance of TOS columns is due to that, as previously discussed, the TOS levels are clear features in the low-lying ES. These columns also have a momentum structure. Consider an ideal magnetic order that consists of $N_{s}$ fully FM sublattices, represented as $\{{\bf\tilde{S}}_{1},{\bf\tilde{S}}_{2},...,{\bf\tilde{S}}_{N_{s}}\}$ ($L_{y}=0\!\mod\!N_{s}$) in a big-S notation of the spins. The $SU(2)$-symmetric groundstate is, of course, the $S_{\text{total}}\!=\!0$-singlet, constructed by adding all spins, $||{\bf\tilde{S}}_{1},{\bf\tilde{S}}_{2},...,{\bf\tilde{S}}_{N_{s}};0\rangle$ in a reduced dimension basis notation (see for example Biedenharn et al., 2009). Importantly, this is the true groundstate of the effective Hamiltonian of $H_{\text{eff}}\propto\frac{1}{\sqrt{L}}{\bf S}_{\text{total}}^{2}$ describing purely the TOS levels Kolley et al. (2013). The only non-trivial sets of unitary symmetry operations that are allowed to act on the $S_{\text{total}}\!=\!0$-singlet and leave a Heisenberg-type Hamiltonian between the sublattices unchanged (sublattices should be still arranged on the physical lattice), can be written as the cyclic translations of sublattices, $T_{\nu}$, where $\nu$ is the number of sublattices that will be shifted (for example to the right). One can then write

[TABLE]

There are obviously only $N_{s}$ distinct values that $\nu$ can take, including the identity operator. Eq. (33) already implies that the TOS levels can only acquire lattice momenta of $k^{\text{TOS}}_{\nu}=\frac{2\pi\nu}{N_{s}}$ for $\nu=0,1,...,N_{s}-1$, between the equal or greater group of general ES momenta, $k_{n}$. The only complication emerges from the distribution pattern of $n^{\prime}$ TOS-levels between $N_{s}$ momenta for a fixed $S$-sector. To clarify this, let us focus on the more general case of $n^{\prime}>N_{s}$ and choose the momentum of the lowest ES level to be $k^{\text{TOS}}_{0}[S\!=\!0]=0$, presumably, corresponding to the action of $I$ on the sublattices (chosen differently in Fig. 17). Trivially, all other $(n^{\prime}-1)$-levels should arrange symmetrically in respect to $k^{\text{TOS}}_{0}[S\!=\!0]$ (there is no relative net momentum). So, they can either, altogether, fill the zero-momentum state on top of $k^{\text{TOS}}_{0}[S=0]$ or occupy $\pm k_{\nu}$ ($\nu\neq 0$) states around it. The former is not possible, due to the fact that $T_{\nu}$ ($\nu\neq 0$) and $I$ posses a distinct set of eigenvalues and therefore produce different momenta (this can be easily observed by writing the bipartite Schmidt decomposition of the $S_{\text{total}}\!=\!0$-singlet state and switch to the basis of fixed-$S$ states for $L$ or $R$ partition to reveal distinct eigenspectra of $T_{\nu}$ and $I$). In addition, we notice that some states appearing in a TOS column are not essentially TOS levels. This is partly due to the fact that the non-TOS levels are also allowed to fill $k^{\text{TOS}}_{\nu}$ states, and partly because in an MPS representation, there is always a fixed number of states kept and consequently, only the first few TOS levels of $H_{\text{eff}}$ will be recovered. Nevertheless, such initial states (having a clear gap to the higher levels) certainly follow the TOS level counting as governed by the degree of symmetry-breaking in the thermodynamic limit. I.e. for a state that fully breaks $SU(2)$-symmetry (e.g. the $120^{\circ}$ order), there are $N_{S}^{\text{TOS}}=(2S+1)$ levels grouped together, and for a state that partially breaks the $SU(2)$-symmetry down to $U(1)$ (e.g. the columnar order), there is only $N_{S}^{\text{TOS}}=1$ level per each fixed $S$-sector (not counting the degeneracy that comes from the $SU(2)$ quantum numbers themselves; the overall degeneracy of the ES levels is always $(2S+1)N_{S}^{\text{TOS}}$ – see Saadatmand, 2017; Kolley et al., 2013 for more details). We discover another striking feature in the momentum-resolved ES of symmetry-broken phases, however, this time for the states between the TOS columns: the first few Nambu-Goldstone modes exhibit sine-like dispersion patterns (as in the energy spectrum), if $L_{y}$ chosen to be large enough.

In Fig. 16, we present the momentum-resolved ES of the $120^{\circ}$ order on different width of the YC structure (for more visibility, we have limited the display of the ES levels to $S_{max}=4$ in all ES figures of this section). The presence of three characterizing TOS columns is clear for all system widths, consistent with the theory for a $N_{s}\!=\!3$-state. The low-lying levels inside the TOS columns (purely TOS levels) have a clear gap to the higher levels, which qualitatively observed to converge to a finite value, linearly with $\frac{1}{L_{y}}$, at the thermodynamic limitSaadatmand (2017). The number of low-lying levels in the TOS columns agree with the full $SU(2)$-symmetry breaking in the thermodynamic limit. That is $N_{S}^{\text{TOS}}=(2S+1)$ for all $S=0,1,2,3,4$, as previously observed by Kolley *et al.*Kolley et al. (2013). For low-lying Nambu-Goldstone modes between the TOS columns, we suggest the triangular-shape dispersion patterns of Fig. 16(c) are signs for the formation of sine-like structures, however, due to relatively small size of $L_{y}$, the $k_{n}$-resolution does not suffice to discern more details.

In Fig. 17, we present the momentum-resolved ES of the columnar order for different widths of the YC structure. The presence of two characterizing TOS columns (note that $k^{\text{TOS}}[\lambda_{0}]=\pm\pi$-columns are the same) is clear for all system widths, as predicted by the theory for a $N_{s}\!=\!2$-state. As before, the low-lying levels inside the TOS columns have a clear gap to the higher levels and observed to converge to a finite value, linearly with $\frac{1}{L_{y}}$, at the thermodynamic limitSaadatmand (2017). The partial breaking of $SU(2)$ to $U(1)$ symmetry can be confirmed by the level counting of $N_{S}^{\text{TOS}}=1$ for low-lying $S=0,1,2,3,4$-levels in the TOS columns. A sine-like dispersion pattern for the low-lying levels between the TOS columns is apparent, at least, for the larger $L_{y}=12$ system, Fig. 17(c).

In Fig. 18, we present the momentum-resolved ES of an ASL state on a $L_{y}=6$ cylinder. Clearly, there is no signature for the presence of TOS columns, which suggests the nonmagnetic nature of the phase. In addition, we observe no non-trivial degeneracy of low-lying ES levels. So, there exist no fractionalization of symmetries to identify SPT and/or some intrinsic topological ordering with anyonic excitations (see also Saadatmand and McCulloch, 2016).

IX Time-reversal symmetry-breaking and the robustness of the

topological phase against the chirality

The existence of the time-reversal symmetry is a key feature of $H_{J_{2}}$, Eq. (1). A chiral groundstate spontaneously breaks time-reversal, $\tau$, and parity reflection, $P$, symmetry, but respects the combined $P\tau$-symmetry. After consistent numerical observations of a nonmagnetic phase in the $J_{1}$-$J_{2}$ THM phase diagram (cf. Sec. I and III), the natural question is, whether the new state stabilizes due to SSB of $\tau$, which would result in a CSL. For a scenario in which the true groundstate in the SL phase region is $\text{Z}_{2}$ topological ordered (advocated by DMRG resultsZhu and White (2015); Hu et al. (2015); Saadatmand and McCulloch (2016)), we already investigatedSaadatmand and McCulloch (2016) the chirality of anyonic sectors in detail, using direct measurement of the $\tau$-operator expectation values and calculating a scalar chiral order parameter,

[TABLE]

where $\langle i,j,k\rangle$ represent a NN triangular plaquette and the sum goes over the wavefunction unit cell. We discovered that the topological sectors are all $\tau$-symmetric as the $O_{\chi}$-values observed to be small and decreasing rapidly to numerically vanishing magnitudes at the thermodynamic limit of $m\rightarrow\infty$ (furthermore, $\hat{\mathrm{b}}$ and $\hat{\mathrm{f}}$-sector are fractionalizing time-reversal symmetry). However, Hu *et al.*Hu et al. (2015) determined the $\hat{\mathrm{i}}$-sector groundstate as strongly prone to the chirality by adding directional (${\bf a}_{\pm 60^{\circ}}$-axis) anisotropy to the Hamiltonian. This is, in part, leading another question of our interest: is the SL phase robust against perturbing $H_{J_{2}}$ with a term that explicitly breaks the $\tau$-symmetry and forms a chiral long-range order? To answer this question, one can study the $J_{1}$-$J_{2}$-$J_{\chi}$ model,

[TABLE]

where $\langle i,j,k\rangle$ indicates the sum over all NN triangular plaquettes in a Hamiltonian unit cell. The phase diagram of $H_{\chi}$ is previously studied using variational QMC (Hu et al., 2016) and ED (Wietek and Läuchli, 2016) techniques, however, no clear result has emerged on the nature of the $J_{\chi}\rightarrow 0$ limit. To shed some lights on this matter, in this section we study the response of the YC8-$\hat{\mathrm{i}}$ groundstatesSaadatmand and McCulloch (2016) to the chiral field by adiabatically adding a $J_{\chi}$-term to $H_{J_{2}}$, as in Eq. (35), and finding new groundstates using the $SU(2)$-symmetric iMPS and iDMRG methods.

We present our results for the extrapolated $O_{\chi}$ in the thermodynamic limit of $m\rightarrow\infty$ in Fig. 19. We notice that, within our resolution, upon varying $J_{\chi}$, there is at least one (significant) point exposed to nonzero chiral perturbations, but has negligible $O_{\chi}(m\rightarrow\infty)$ within the error-bars. This means that the topological SL phase is robust against chirality and one needs to provide $\tau$-symmetry-breaking terms larger than a finite-value, namely $J_{\chi}^{\text{critical}}$, to impose a chiral groundstate. To further predict this small $J_{\chi}^{\text{critical}}$, we applied a fit of $\tilde{b}_{0}+\tilde{b}_{1}J_{\chi}^{\tilde{b}_{2}}$ to the data and find $J_{\chi}^{\text{critical}}=(-\frac{b_{0}}{b_{1}})^{\frac{1}{b_{2}}}=0.0014(1)$. These results also suggest the existence of a second order phase transition toward the CSL phase. This is consistent with the suggestion from Wietek and LäuchliWietek and Läuchli (2016), and may clarify the results of Hu *et al.*Hu et al. (2016), where it is unclear if $O_{\chi}$ would be zero or not in the $J_{\chi}\rightarrow 0$ limit.

X Conclusion

We have presented comprehensive results for the phase diagram of the $J_{1}$-$J_{2}$ Heisenberg model on triangular lattices, using infinite-length YC structures. Using the Binder ratio of the magnetization order parameter, $U_{r}$, Eq. (31), and TOS columns of the momentum-resolved ES, we have obtained phase boundaries and characterized the nature of the symmetry breaking magnetic order. We found that the Binder ratio reliably detects phase boundaries between magnetically ordered states, even when using $SU(2)$ symmetry, where the order parameter itself is zero by construction. We identified the $120^{\circ}$-ordered groundstate as a three-sublattice LRO with full $SU(2)$-symmetry-breaking in the thermodynamic limit; the columnar-ordered groundstate as a two-sublattice LRO with partial $SU(2)$-symmetry-breaking at the thermodynamic limit, and confirm the nonmagnetic nature of the SL states on infinite cylinders of widths up to 12 sites. In addition, we have discovered the stabilization of a new ASL phase, with power-law correlation lengths, for width-6 infinite cylinders. We have pinpointed the phase transitions between the infinite cylinder’s groundstates of the THM, precisely, using the Binder ratios. The transitions are relatively close to the phase boundaries found from the direct measurements of the local order parameters using fDMRG on $L_{y}=3,4,5,6$-cylinders, and short-range correlations and fidelity susceptibility phase diagrams from iDMRG calculations. In addition, for the columnar order, we have numerically proved that the entropies consistently obey $S_{EE}=a_{0}(L_{y})+(\alpha_{0}+\alpha_{1}L_{y})\log(\xi)$, a mixture of the area-law and the quantum critical behavior, as expected for the magnetic phases built by the inherently one-dimensional $SU(2)$-symmetric iMPS ansatz. To the best of our knowledge, a set of numerical tools to efficiently distinguish and classify LROs were previously absent in the $SU(2)$-symmetric iDMRG literature. Considering the advantages of $SU(2)$-symmetric calculations, we suggest that the proposed methods can be applied widely to detect symmetry broken states using the iMPS.

Finally, to unravel the true nature of time-reversal symmetry-breaking in the topological SL, we have investigated the robustness of YC8-$\hat{\mathrm{i}}$ sector under perturbing $H_{J_{2}}$ with a chiral term, Eq. (35) (it was previously suggested (Hu et al., 2015) that YC8-$\hat{\mathrm{i}}$ states are prone to become chiral under applying bond anisotropies to the Hamiltonian). The results of the scalar chiral order parameter, $O_{\chi}(m\rightarrow\infty)$, versus $J_{\chi}$ can be fitted using $\tilde{b}_{0}+\tilde{b}_{1}J_{\chi}^{\tilde{b}_{2}}$ with high accuracy and shows the existence of a continuous phase transition to the CSL phase at small, but non-zero, $J_{\chi}^{\text{critical}}=0.0014(1)$. Therefore, for finite-width cylinders the topological state of the THM are time-reversal symmetric, and not a chiral topological liquid.

Acknowledgements.

The authors would like to thank Jason Pillay for useful discussions. This work has been supported by the Australian Research Council (ARC) Centre of Excellence for Engineered Quantum Systems, grant CE110001013. I.P.M. also acknowledges support from the ARC Future Fellowships scheme, FT140100625. Notes added. – After completing this work, a related paper (Gong et al., 2017) appeared in which the authors study the phase diagram of the $J_{1}$-$J_{2}$-$J_{\chi}$ model, Eq. (35), on finite-$L_{x}$ cylinders using the $SU(2)$-symmetric fDMRG algorithm. In agreement to Sec. IX results, Gong et al. find a smooth phase transition from $J_{1}$-$J_{2}$ SL to a CSL at a small but finite chiral coupling strength ($J_{\chi}^{\text{critical}}\approx 0.02$ for $J_{2}=0.1$). In addition, Fig. 19 caption should specify $J_{2}=0.125$.