3. Phase diagram of the SU(\(N\)) antiferromagnet of spin \(S\) on a square lattice

The results presented in this chapter are the outcome from a collaboration with Francesco Parisen Toldin and Fakher F. Assaad. These findings have been published in [17], with significant portions reproduced here verbatim. My contribution to the project comprise the quantum Monte Carlo (QMC) calculations, including the implementation of the model in code. The QMC model has been designed by F. F. A., while the group-theoretical proof for the projections has been executed by F. P. T.. The interpretation of data and written text is a combined work of all authors.

3.1. Introduction

Spin systems are ubiquitous in nature and form one of the most fundamental concept in condensed matter and statistical physics. Their complex collective behavior has spurred numerous experimental and theoretical studies, aimed at understanding their nature and properties. At the same time, modeling of spin systems represents a primary theoretical laboratory to investigate fundamental physics. Starting with the classical Ising model [33], spin systems have played an crucial role in our understanding of phase transitions [37], phases of matter, frustration and disorder [73], emergent gauge theories [74, 75] and exotic critical behavior [76]. The impact of spin models extends beyond the realm of condensed matter physics, and has found application in other areas, such as information processing [77] and quantum computing [78], where the fundamental unit of information, a qubit, is a single spin-\(1/2\) system.

Fig. 3.1 Ground-state phase diagram of the SU(\(N\))-antiferromagnet model (3.1) on the square lattice, as obtained from QMC simulations. \(S\) identifies the chosen representation of the \(\mathfrak{su}(N)\) algebra of SU(\(N\)), illustrated by the Young tableau in Fig. 3.2. Striped regions indicate the part of the phase diagram where current QMC data do not allow an unambiguous identification of the phase; in such cases we indicate between parenthesis the most likely identified order. The insets show QMC data in the highlighted dimerized phases, obtained through a pinning-field approach (see Sec. 3.4.1).

In condensed matter, spin systems are realized in Mott insulators, which arise when charge fluctuations in a given unit cell are suppressed. For instance, in undoped cuprates the copper atom is in a Cu\(^{2+}\) state and corresponds to a net spin \(S=1/2\) degree of freedom. Super-exchange leads to an \(S=1/2\), SU(2) Heisenberg spin model that has been studied numerically [79, 80] and experimentally [81] at length. Higher spin SU(2) systems arise when \(2S\) electrons are localized on a single orbital and a strong Hund’s rule favors a maximal spin state with a totally symmetric wave function. For example, in the Haldane chain realized by the CsNiCl\(_3\) compound Ni\(^{2+}\) ions carry spin 1 [82]. SU(\(N\))-invariant models, for \(N>2\), naturally arise as special cases of the Kugel-Khomski model [83, 84], where spin and orbital degrees of freedom turn out to play a very symmetric role. In particular, the observed spin-orbital liquid behavior in Ba\(_3\)CuSb\(_2\)O\(_9\) [85] has been interpreted in terms of an SU(4) quantum antiferromagnet in the defining representation [86]. Beyond the solid state physics, SU(\(N\)) spin models can be realized in the realm of cold atomic gases [87, 88].

Fig. 3.2 Young tableau corresponding to the irreducible representation of \(\mathfrak{su}(N)\) considered here; \(N\) is even and \(S\) can take half-integer and integer values.

Topology plays a decisive role in the understanding of SU(2) invariant spin systems. In fact, using a spin coherent-state path integral approach to antiferromagnetic (AFM) Heisenberg chains, one identifies a Berry phase. It corresponds to the skyrmion count of the three-components normalized order parameter in 1+1 dimensions and at angle \(\theta= 2\pi S\) [89]. This provides a topological understanding of the observed differences between half-integer and integer spin chains. In two spatial dimensions, topology enters through singular skyrmion number changing events in space time: monopoles [90]. For a square lattice with \(C_4\) symmetry, only quadrupole (double) monopole events are allowed for half-integer (odd) spin by symmetry. There is no constraint on the monopole number for even values of the spin. For the plain vanilla SU(2) Heisenberg model at arbitrary spin \(S\), the spin-wave approximation captures well the ground state and topological excitations lie high in the spectrum. In this context, the theory of deconfined quantum criticality essentially poses the question of the nature of the quantum phase transition that emanates when one decreases the energy of monopoles and ultimately condenses them [76]. For half-integer spin systems, where only quadrupole monopole insertions are allowed, one can conjecture that the Hilbert space splits into four orthogonal subspaces characterized by the number of monopoles modulo four. This provides an understanding of how the fourfold degenerate valence bond solid (VBS) state emerges for condensing topological excitations of the quantum antiferromagnet [91]. Similarly, for spin-1 (spin-2) systems, condensing monopoles should generate a twofold (zerofold) degenerate disordered state.

A crucial question is how to control the monopole energy. The seminal work of Read and Sachdev [18, 20, 21] shows that the discussion above can be carried over to SU(\(N\)) spin systems, \(N\ge 2\). Furthermore, enhancing \(N\) has the potential of lowering the monopole energy. In this work, we show that it is possible to formulate negative sign-free auxiliary-field (AF) quantum Monte Carlo (QMC) simulations [8, 9, 10, 12, 92] of the SU(\(N\)) AFM spin-\(S\) Heisenberg model, for representations given by a Young tableau with \(N/2\) rows and \(2S\) columns. This generalizes the work of Ref. [93] to generic values of \(S\). Specifically, we consider the model:

(3.1)\[ \hat{H} = J\sum_{\langle i, j \rangle, a } \hat{S}^{(a)}_i\hat{S}^{(a)}_j,\]

where the sum extends over the pair of nearest-neighbor sites \(\langle i, j \rangle\), and \(a\) runs over \(N^2 - 1\) generators of the said representation of the \(\mathfrak{su}(N)\) algebra of SU(\(N\)). The main result of this work is the rich phase diagram illustrated in Fig. 3.1. Remarkably, and for each considered value of \(S=1/2,1,3/2,2\) just above the threshold value of \(N\) above which Néel order disappears, we observe four-, two- and zerofold degenerate disordered states at half-integer, odd and even values of \(S\).

This chapter is organized as follows. In Sec. 3.2 we discuss how we construct the Hamiltonian (3.1), with the spin operators in the desired representation of Fig. 3.2. In Sec. 3.3 we illustrate its actual implementation within a fermionic representation, which can be sampled by means of QMC simulations in the AF approach. In Sec. 3.4 we present and discuss our QMC results for the phase diagram of the model. In Sec. 3.5 we summarize our findings. In Appendix B.1 we discuss a formula giving the eigenvalue of the quadratic Casimir operator of a representation in terms of its Young tableau. In Appendix B.2 we prove an upper bound on the eigenvalue of the Casimir operator of the irreducible representations emerging from a tensor product of representations discussed in Sec. 3.2. In Appendix B.3 we discuss the systematic error in the QMC formulation, arising from the Trotter discretization. In Appendix B.4 we prove a lower and upper bound for a bond observable used to diagnose the phases.

3.2. General formulation of the Hamiltonian

In the Hamiltonian (3.1), the operators \(S^{(a)}_i\) form an irreducible representation of the \(\mathfrak{su}(N)\) algebra. This is uniquely specified by its maximum Dynkin weight \(\Lambda_{\alpha}\) or, alternatively, by a Young tableau, from which the components \(\Lambda_{\alpha_k}\) can be read off as [94]

(3.2)\[ \Lambda_{\alpha_k} = l_k - l_{k+1},\qquad k=1,\ldots, N-1,\]

where \(l_k\) is the length of the \(k\)-th row of the Young tableau, and one can assume \(l_N=0\) for representations of \(\mathfrak{su}(N)\).

Here we consider, on each lattice site, the representation corresponding to a Young tableau illustrated in Fig. 3.2, whose corresponding maximum Dynkin weight is

(3.3)\[ \Lambda_{\alpha_k}=2S\delta_{k,N/2}.\]

The dimension of an irreducible representation can be computed with the hook-length formula, or with Weyl’s formula [94]:

(3.4)\[ \text{dim} = \prod_{i<j}^N \frac{l_i-l_j+j-i}{j-i}.\]

For the present case, we have

(3.5)\[ \text{dim}=\prod_{j=0}^{N/2-1}\frac{\left(2S+\frac{N}{2}+j\right)! j!}{\left(2S+j\right)!\left(\frac{N}{2}+j\right)!}.\]

To realize this representation, we first introduce on each lattice site \(2S\) independent irreducible representations. Their tensor product decomposes into different irreducible representations, including, in particular, the one of Fig. 3.2. In a second step we project the Hilbert space onto that of the desired representation by maximizing the quadratic Casimir operator.

Let \(\{T_a\}\), \(a=1,\ldots, N^2-1\) be a basis of the \(\mathfrak{su}(N)\) algebra. We start by introducing on each lattice site \(i\) the antisymmetric self-adjoint representation \(T_a\rightarrow \Gamma(T_a) = \hat{T}_{a,i}\). Its maximum weight in the Dynkin representation is

(3.6)\[ \Lambda_{\alpha_k} = \delta_{k,N/2},\]

which matches Eq. (3.3) for \(S=1/2\). Equivalently, in agreement with Eq. (3.2), this representation corresponds to a Young tableau with one column and \(N/2\) boxes.

Fig. 3.3 Decomposition of the tensor product of \(2S\) antisymmetric self-adjoint representations, whose maximum Dynkin weight is given in Eq. (3.6), into irreducible ones.

Next, we consider, for each lattice site \(i\), \(2S\) independent representations \(\hat{T}_{a,i,\alpha}\), \(\alpha=1\ldots 2S\). We refer to \(\alpha\) as the flavor index. The composite generators, i.e., the generators for the tensor product of the \(2S\) representations, define the spin operators appearing in Eq. (3.1) and are given by

(3.7)\[ \hat{S}_i^{(a)} = \sum_{\alpha=1}^{2S} \hat{T}_{a,i,\alpha},\]

Using Eq. (3.7), the interaction term in Eq. (3.1) is written as [1]

(3.8)\[\begin{split} \hat{H}_J &= J \sum_{\langle i\ j \rangle} \sum_{a=1}^{N^2-1}\hat{S}^{(a)}_i \hat{S}^{(a)}_j \\ &= J \sum_{\langle i\ j \rangle} \sum_{a=1}^{N^2-1}\sum_{\alpha,\beta=1}^{2S} \hat{T}_{a,i,\alpha}\hat{T}_{a,j,\beta}.\end{split}\]

The operators \(\hat{S}_i^{(a)}\) in Eq. (3.7) form a reducible representation of \(\mathfrak{su}(N)\), which decomposes into several irreducible representations, illustrated in Fig. 3.3. As proven in App. B.2, among the resulting representations, the one of Fig. 3.2 exhibits the maximum eigenvalue of the quadratic Casimir operator. To explicitly compute it, we choose a basis \(\{T_a\}\) of \(\mathfrak{su}(N)\) such that

(3.9)\[ \text{Tr}\{T_aT_b\} = \frac{1}{2}\delta_{ab}.\]

With this choice, the structure constants of the algebra are completely antisymmetric and the chosen basis is, up to a trivial normalization, self-dual with respect to the bilinear form (3.9). Thus, given an irreducible representation \(\Gamma: \mathfrak{su}(N)\rightarrow GL(d,\mathbb{C})\), we define the quadratic Casimir operator as

(3.10)\[ \hat{C}_2 = \sum_a \Gamma(T_a)\Gamma(T_a) \equiv C\id_d,\]

where we have used the fact that \(C_2\propto \id_d\) (Schur’s Lemma) to introduce the eigenvalue of the Casimir operator \(C\) [2].

Using Eq. (3.7) in Eq. (3.10), the quadratic Casimir operator on the lattice site \(i\) is

(3.11)\[ \hat{C}_{2,\Gamma_i} = \sum_{a=1}^{N^2-1}\sum_{\alpha,\beta=1}^{2S}\hat{T}_{a,i,\alpha}\hat{T}_{a,i,\beta}\]

In order to project the Hilbert space to the subspace of the desired representation, we introduce on each site a term in the Hamiltonian which favors the states with the highest Casimir value

(3.12)\[\begin{split} \hat{H}_{\text{Casimir}} &= -J_H\sum_i \hat{C}_{2,\Gamma_i}\\ &= -J_H \sum_i\sum_{a=1}^{N^2-1}\sum_{\alpha,\beta=1}^{2S}\hat{T}_{a,i,\alpha}\hat{T}_{a,i,\beta},\end{split}\]

with \(J_H>0\). The term of Eq. (3.12) effectively introduces a ferromagnetic interaction between different flavors, with coupling strength \(J_H\).

The Hamiltonian that we will solve numerically, reads:

(3.13)\[\begin{equation} \hat{H} = \hat{H}_J + \hat{H}_{\text{Casimir}}. \end{equation}\]

Importantly, \(\left[ \hat{H}_J, \hat{H}_{\text{Casimir}} \right] = 0\), such that the projection onto the desired irreducible representation turns out to be very efficient. And since the projection is a local onsite term, we expect it to scale independent from system size.

3.3. QMC formulation

3.3.1. Fermionic representation

As discussed in Sec. 3.2, the Hamiltonian is constructed using as basic building blocks antisymmetric self-adjoint representations, defined by the maximum weight of Eq. (3.6) or, equivalently, by a Young tableau with one column and \(N/2\) boxes. The corresponding operators \(\hat{T}_{a,i,\alpha}\) entering in Eqs. (3.8) and (3.12) can be realized by introducing, for every lattice site \(i\) and for every flavor index \(\alpha\), \(N\) nonrelativistic fermions, with creation and annihilation operators \(\cdag{i,\alpha,\sigma}\), \(\cd{i,\alpha,\sigma}\), \(\sigma=1\ldots N\), and fixing the total charge (i.e., the number of fermions) to half-filling, i.e., to \(N/2\). For every \(i\) and \(\alpha\), a basis of this Hilbert space is generated by the states

(3.14)\[\begin{equation} \left( \cdag{i,\alpha,1}\right)^{n_{i,\alpha,1}} \cdots \left(\cdag{i,\alpha,N}\right)^{n_{i,\alpha,N}} | 0 \rangle,\qquad n_{i,\alpha,\sigma} = 0,1, \end{equation}\]

with the constraint

(3.15)\[ \sum_{\sigma=1}^{N} n_{i,\alpha,\sigma} = \frac{N}{2}, \qquad \forall i, \alpha.\]

In this space, the (representation of the) \(\mathfrak{su}(N)\) generators are

(3.16)\[ \hat{T}_{a,i,\alpha} = \sum_{\sigma,\sigma'}\cdag{i,\alpha,\sigma}(T_a)_{\sigma\sigma'}\cd{i,\alpha,\sigma'}.\]

It is easy to check that the maximum weight in the Dynkin representation agrees with Eq. (3.6), thus providing us the needed building block to simulate the Hamiltonian (3.1).

We study the model by means of finite-temperature AF QMC [8, 9, 10] and projective AF QMC [10, 92, 95]. In this framework, we sample respectively the grand canonical and canonical ensembles at half-filling, and charge fluctuations are generally present. Therefore, we need to additionally impose the constraint of Eq. (3.15). Notice that, unlike available techniques for canonical QMC simulations [96, 97], where the global charge of the system is fixed, here we need to impose half-filling on each lattice site. To this end, we add a repulsive Hubbard \(U\)-term on each site \(i\) and flavor \(\alpha\):

(3.17)\[ \hat{H}_U = U\sum_i \sum_{\alpha=1}^{2S}\left(\hat{n}_{i,\alpha}-\frac{N}{2}\right)^2, \quad \nop{i,\alpha} \equiv \sum_{\sigma=1}^{N} \cdag{i,\alpha,\sigma}\cd{i,\alpha,\sigma}.\]

Fig. 3.4 Sketch of the structure of the QMC Hamiltonian for \(2S=2\). \(\color{green} \hat{H}_U\): Hubbard term for freezing out charge degrees of freedom. \(\color{blue} \hat{H}_\text{Casimir}\): Term for maximizing the eigenvalue of the Casimir operator. \(\color{red} \hat{H}_J\): Antiferromagnetic interaction between elemental spins.

In summary, the Hamiltonian simulated with the AF QMC method is the sum of the interaction term given in Eq. (3.8), the Casimir term [Eq. (3.12)], and the Hubbard term [Eq. (3.17)], with the operators \(\{\hat{T}_{a,i,\alpha}\}\) given in Eq. (3.16). Eqs. (3.8) and (3.12) can be further simplified using the following summation identity [98]

(3.18)\[ \sum_{a=1}^{N^2-1} (T_a)_{\sigma\sigma'}(T_a)_{\epsilon\epsilon'} = \frac{1}{2}\left(\delta_{\sigma\epsilon'}\delta_{\sigma'\epsilon}-\frac{1}{N}\delta_{\sigma\sigma'}\delta_{\epsilon\epsilon'}\right),\]

which holds for a choice of generators that satisfies Eq. (3.9). Using Eq. (3.18) and collecting the terms in Eqs. (3.8), (3.12) and (3.17), the QMC Hamiltonian is

(3.19)\[\begin{split} \begin{split} \hat{H}_{\text{QMC}} =&\hat{H}_J + \hat{H}_{\text{Casimir}} + \hat{H}_U \\ =&- \frac{J}{4} \sum_{\langle i,j \rangle, \alpha,\beta } \left\{ \hat{D}_{(i,\alpha),(j,\beta)}, \hat{D}^{\dagger}_{(i,\alpha),(j,\beta)} \right\} \\ &+ \frac{J_H}{2} \sum_{i} \sum_{\alpha>\beta} \left\{ \hat{D}_{(i,\alpha),(i,\beta)}, \hat{D}^{\dagger}_{(i,\alpha),(i,\beta)} \right\}\\ &+ U\sum_{i,\alpha}\left(\nop{i,\alpha}-\frac{N}{2}\right)^2,\\ \end{split}\end{split}\]

where

(3.20)\[ \hat{D}_{(i,\alpha),(j,\beta)} \equiv \sum_{\sigma} \cdag{i,\alpha,\sigma}\cd{j,\beta,\sigma},\]

\(\left\{\hat{A}, \hat{B}\right\} \equiv \hat{A} \hat{B} + \hat{B} \hat{A}\), and \(\nop{i,\alpha}\) as defined in Eq. (3.17). The Hamiltonian now takes the form of the Heisenberg model considered in Ref. [93] and the proof for the absence of sign problem is similar. In Fig. 3.4 we sketch the resulting interactions for the case \(S=1\).

Before proceeding, we would like to comment on the computational cost of the AF QMC algorithm [10] for this model. The total number of orbitals is given by \(L^2 2S \) such that matrix operations required to compute, e.g., the single-particle spectral function, scales as \( (L^2 2S)^3 \beta \), where \(\beta\) is the inverse temperature. It turns out, that, in contrast to the generic Hubbard model with \(L^2 2S\) sites, this is not the leading computational cost. The number of Hubbard-Stratonovich fields per imaginary time slice scales as \(L^2 S^2\). Using fast updates, refreshing one field involves \((L^2 2S)^2\) floating point operations, such that the total cost of the updating scales as \(L^6S^4\beta\). Hence large values of \(S\) are computationally expensive. In Appendix B.3, we show that the computational cost does not explicitly scale with \(N\). We note that this estimate of the computational cost does not take into account auto-correlation times.

3.3.2. Test of projections

As discussed above, the Hamiltonian (3.19) is equivalent to Eq. (3.1) in the limit \(J_H\rightarrow\infty\), and \(U\rightarrow\infty\), under which the Hilbert space is projected to the representation of Fig. 3.2. To optimally test the projections, we use the finite-temperature AF QMC method, which evaluates \(\langle \hat{O} \rangle = \text{Tr} \left[ e^{- \beta \hat{H}} \hat{O} \right] / \text{Tr} \left[ e^{- \beta \hat{H}} \right] \), where the trace runs over the grand canonical ensemble.

The interaction term of Eq. (3.8) and the Casimir term of Eq. (3.12) manifestly conserve the charge on each lattice site \(i\). Hence, in the Gibbs density matrix \(\exp(-\beta\hat{H})\), the Hubbard term factorizes out, resulting in an effective exponential suppression of the charge fluctuations,

(3.21)\[ \Braket{\left( \hat{n}_{i,\alpha}- N/2 \right)^2} \propto e^{- \beta U},\]

independent from system size. The suppression of charge fluctuations is therefore particularly efficient. This is illustrated in Fig. 3.5, where we show \(\<\left( \hat{n}_{i,1}- N/2 \right)^2\>\) in a semilogarithmic scale for \(N=2\) and \(S=1/2\) and as a function of \(\beta U\). Besides the case of a Hamiltonian containing the Hubbard interaction [Eq. (3.17)] only, for which any observable depends only on \(\beta U\), we consider the presence of the AFM interaction Eq. (3.8) for a lattice of linear size \(L=4\). In the latter case, there is an additional dependence on the inverse temperature \(\beta\), which we illustrate by considering four values. In line with Eq. (3.21), we observe an exponential suppression of the charge fluctuations as a function of \(\beta U\). Interestingly, in the interacting case, \(\<\left( \hat{n}_{i,1}- N/2 \right)^2\>\) decreases with the temperature for any given value of \(\beta U\), even for \(U=0\). This implies that the AFM coupling itself suppresses the charge fluctuations.

Fig. 3.5 Suppression of the charge fluctuations \(\<\left( \hat{n}_{i,1}- N/2 \right)^2\>\) as a function of \(\beta U\), for \(N=2\) and \(S=1/2\). We consider a Hamiltonian containing the Hubbard term only and the case of a model with an antiferromagnetic interaction [Eq. (3.8)], with coupling constant \(J=1\) on a system size \(L=4\), and for different inverse temperatures \(\beta\). The charge fluctuations fall off asymptotically as \(\exp(-\beta U )\) [Eq. (3.21)].

In Appendix B.1, we discuss a formula that gives the value of the Casimir eigenvalue in terms of the Young tableau of the representation [99]. Employing this result, in App. E.2 we determine, for the representation of Fig. 3.2:

(3.22)\[ C(N, S) = \frac{NS(2S+N)}{4}.\]

Furthermore, in Appendix B.2, we prove that Eq. (3.22) is the maximum Casimir eigenvalue among the irreducible representations arising from the tensor product of \(2S\) self-adjoint antisymmetric representations given in Eq. (3.16), and that there is a finite gap \(O(1)\) in the eigenvalues of the quadratic Casimir operators between the maximally symmetric representation of Fig. 3.2 and the other irreducible representations arising from the tensor product. Therefore, the term of Eq. (3.12) effectively selects a single representation, and the projection is efficient.

To control the projection, we compute the expectation value of the quadratic Casimir operator from the QMC simulations and compare it with the expected result of Eq. (3.22). An example of such a projection is shown in Fig. 3.6. In Fig. 3.6(a), we plot the difference between the computed and expected Casimir eigenvalue \(C\), as a function of \(\beta J_H\), for different inverse temperatures and in a semilogarithmic scale. The deviation from the expected result is exponentially suppressed in \(\beta J_H\), underscoring the effectiveness of the projection. In Fig. 3.6(b), we show, as a function of \(N\), the sampled value of \(C\) along with the expected result, and in the inset we plot their difference, which vanishes within error bars.

Fig. 3.6 Demonstrating the effectiveness of projection onto the fully symmetric representation for \(S=1\) by comparing the Casimir eigenvalue \(C(N, S)\) [Eq. (3.22)] to the sampled one \(\<\hat{C}\>\). (a) Difference between \(C(N=2, S=1)\) of the representation \(S=1\), \(N=2\), and \(\<\hat{C}\>\), as a function of the effective interaction strength \(\beta J_H\) [Eq. (3.12)] with \(J=0\), and for four inverse temperatures. Data shown are obtained for a lattice of size \(L=4\), with a Hubbard interaction \(\beta U=6\) [Eq. (3.17)], vanishing nearest-neighbor antiferromagnetic interaction \(J=0\), and a Trotter discretization \(\Delta\tau=0.1\). (b) Casimir eigenvalue as a function of \(N\). We compare the predicted value of Eq. (3.22) with the sampled Casimir eigenvalue from QMC simulations of a lattice with size \(L=4\), with a Hubbard interaction \(U=2\) [Eq. (3.17)], antiferromagnetic coupling \(J=1\) [Eq. (3.1)], projection strength \(J_H=1\) [Eq. (3.12)], and a Trotter discretization \(\Delta\tau=0.1\). In the inset we plot the difference between the sampled and expected value.

3.4. Results

3.4.1. Order parameters and phases

We have simulated the Hamiltonian Eq. (3.19) using the ALF package [11, 12], which provides a comprehensive library to program QMC simulations of interacting models of fermions, using the AF algorithm [8, 9, 10]. In particular, we used the projective formulation of the algorithm, which projects a trial wave function \(\ket{\Psi_{\rm T}}\) onto the ground state of the system. Observables are evaluated through

(3.23)\[\braket{\hat{O}} = \frac{ \Braket{\Psi_{\rm T}| e^{-\Theta\hat{H}} \hat{O} e^{-\Theta\hat{H}} |\Psi_{\rm T}} }{ \Braket{\Psi_{\rm T}| e^{-2\Theta\hat{H}} |\Psi_{\rm T}} },\]

with \(\Theta\) the projection parameter. The algorithm employs a Hubbard-Stratonovich decomposition of the interaction terms. This results in a free fermionic system, where any observable can be computed via the Wick’s theorem from the Green’s functions. The QMC method consists in a stochastic sampling of the Hubbard-Stratonovich fields. We refer to Ref. [10] for a discussion of the AF QMC method.

As trial wave function \(\ket{\Psi_{\rm T}}\), we used the half-filled ground state of

(3.24)\[\begin{align} \hat{H}_{\rm T} = \sum_{\langle i,j \rangle} \sum_{\alpha=1}^{2S} \left( \hat{D}_{(i,\alpha),(j,\alpha)} + {\rm H.c.} \right). \end{align}\]

We scaled the projection parameter \(\Theta\) with linear system size \(L\), usually comparing the results obtained with \(\Theta = L/4\) and \(\Theta = L/2\), ensuring that they reflect ground state properties. Furthermore, we chose the parameters for suppression of charge fluctuations and projection onto the maximally symmetric representation around \(U= 4/\Theta\), \(J_H = 4/\theta\), while always checking that charge fluctuations are sufficiently suppressed and \(\braket{\hat{C}} = C(N, S)\) [cf. Eq. (3.22)].

To detect the realization of different ground states, we have sampled the spin two-point function \(S(\ve{k})\) and the correlations of the dimer operator \(D_{ij}(\ve{k})\) in momentum space, defined as

(3.25)\[S(\ve{k}) \equiv \frac{1}{(N^2 -1)N_r^2} \sum_{\ve{r},a} e^{i \ve{k} \ve{r}} \Braket{\hat{S}^{(a)}_{\ve{0}} \hat{S}^{(a)}_{\ve{r}}},\]

(3.26)\[\begin{split}D_{ij}(\ve{k}) \equiv \frac{1}{(N^2 -1)N_r^2} &\sum_{\ve{r},a,b} e^{i \ve{k} \ve{r}}\cdot \\ \Bigg[ &\Braket{\left(\hat{S}^{(a)}_{\ve{0}} \hat{S}^{(a)}_{\ve{0}+\ve{e}_i}\right) \left(\hat{S}^{(b)}_{\ve{r}} \hat{S}^{(b)}_{\ve{r}+\ve{e}_j}\right) } - \Braket{\hat{S}^{(a)}_{\ve{0}} \hat{S}^{(a)}_{\ve{0}+\ve{e}_i}} \Braket{\hat{S}^{(b)}_{\ve{r}} \hat{S}^{(b)}_{\ve{r}+\ve{e}_j}} \Bigg],\end{split}\]

where \(N_r=L^2\) is the number of sites in a lattice of linear size \(L\) and \(\ve{e}_i\) is the elementary lattice unit vector on the \(i-\)th direction. The normalization in Eqs. (3.25) and (3.26) ensure a finite thermodynamic and large-\(N\) limit. Using these observables, we can distinguish the Néel state and different dimerized ground states, to be discussed below.

The AFM Néel state exhibits long-range spin-spin correlations at momentum \(\ve{k}=(\pi, \pi)\). Thus, it can be detected by the staggered magnetization \(m\)

(3.27)\[ m^2 = S(\ve{k}=(\pi,\pi)).\]

Fig. 3.7 Sketch of possible dimerized ground states. (a)-(c) The VBS states. Each of those states is fourfold degenerate, the corresponding states can be obtained by rotations and translations. (d) A “Haldane nematic” state, an equivalent state is obtained by rotations of \(90^{\circ}\). (e) A unique ground state. States (d) and (e) are at best understood within an AKLT construction, in which, very much as done in our calculation, the spin \(S\) on each site is constructed by a totally symmetric superposition of \(2S\) states that are denoted by bullets around each site on the right-hand side of (d) and (e). These bullets correspond to an irreducible representation of \(\mathfrak{su}(N)\) with one column (\(S=1/2\)) and \(N/2\) rows. In the nematic state, each spin-\(1/2\) forms a singlet with the nearest neighbor along the axis of the broken symmetry. The AKLT state is relevant for the \(S=2\) state, where each spin-\(1/2\) on a given site can be combined into a singlet with a nearest-neighbor spin-\(1/2\) without breaking a lattice symmetry.

The valence bond state (VBS) breaks the lattice rotation and translation symmetries, realizing a fourfold degenerate pattern of strong and weak dimers. This is realized by different sets of bond configurations, illustrated in Figs. 3.7(a)-(c). Beyond the commonly identified columnar order, sketched in Fig. 3.7(a), there are two additional VBS states [Figs. 3.7(b) and (c)]. Notably, all three patterns break the lattice translation symmetry, but only columnar and ladder order break the four-fold rotation symmetry. VBS order can be detected by a suitable order parameter \(\phi\) defined in terms of the dimer correlations

(3.28)\[ \phi^2 \equiv D_{xx}(\ve{k}) + D_{yy}(\ve{k}), \qquad \ve{k}=(\pi, 0).\]

We average \(\phi\) over the two equivalent momenta \((\pi,0)\), and \((0,\pi)\) as to obtain an improved estimator.

For integer values of \(S\), we investigate the possible realization of the “Haldane nematic” Affleck-Kennedy-Lieb-Tasaki (AKLT) phase [18, 19, 20, 21]. This state is twofold degenerate and breaks the rotational symmetry but, unlike the VBS state, does not break translational symmetry. We illustrate it in Fig. 3.7(d). For such a phase we have \(\phi=0\) and a suitable order parameter can be defined as [100]

(3.29)\[ \psi^2 \equiv D_{xx}(\ve{k}) + D_{yy}(\ve{k}) - D_{xy}(\ve{k}) - D_{yx}(\ve{k}), \qquad \ve{k}=\ve{0}.\]

\(\psi\) is designed to pick up rotation symmetry breaking in the dimers and therefore does also not vanish for columnar and ladder order. Therefore, \(\psi\) distinguishes plaquette with \(C_4\) symmetry from VBS order with broken \(C_4\) symmetry.

Finally, a two-dimensional version of the AKLT phase, with singlets on all bonds as sketched in Fig. 3.7(e), is also possible for \(S=2\). This non-degenerate state does not break any symmetry, therefore all previously defined order parameters vanish. In Table 3.1 we summarize the different orders.

Table 3.1 List of considered ground states with their ordering momenta in reciprocal space and matrix of order parameters defined in Eqs. (3.27), (3.28), (3.29).

Phase	Ordering momenta	\(C_4\) Lattice symmetry preserved	\(m\)	\(\phi\)	\(\psi\)
Néel	\((\pi,\pi)\)	yes	\(\ne 0\)	\(0\)	\(0\)
Columnar	\((\pi,0)\) or \((0,\pi)\)	no	\(0\)	\(\ne 0\)	\(\ne 0\)
Plaquette	\((\pi,0)\) and \((0,\pi)\)	yes	\(0\)	\(\ne 0\)	\(0\)
Ladder	\((\pi,0)\) or \((0,\pi)\)	no	\(0\)	\(\ne 0\)	\(\ne 0\)
Nematic	\((0,0)\)	no	\(0\)	\(0\)	\(\ne 0\)
2d AKLT	\((0,0)\)	yes	\(0\)	\(0\)	\(0\)

Same as in Chapter 2, we use the correlation ratio \(R_O\) to pinpoint the order:

(3.30)\[ R_O \equiv 1 - \frac{O(\ve{p}+\delta\ve{p})}{O(\ve{p})},\]

where \(O(\ve{p})\) is the two-point function of the order parameter in Fourier space, and \(\delta\ve{p}\) is the minimum nonzero momentum on a finite lattice. On the square lattice, \(\delta\ve{p}=(2\pi/L,0)\) or \(\delta\ve{p}=(0,2\pi/L)\); as usual, one can average over the two minimum displacements to obtain an improved estimator. The correlation ratio is closely related to the second-moment finite-size correlation length \(\xi\), which on a square lattice can be defined as [101, 102]

(3.31)\[ \xi = \frac{1}{2\sin (\pi/L)}\sqrt{\frac{O(\ve{p})}{O(\ve{p}+\delta\ve{p})}-1}.\]

In a disordered phase, \(\xi\) as defined in Eq. (3.31) converges to the second-moment correlation length for \(L\rightarrow\infty\), such that \(\xi/L\rightarrow 0\) and \(R_0\rightarrow 0\). In an ordered phase, due to the lack of spontaneous symmetry breaking in any finite size, \(\xi/L\) diverges for \(L\rightarrow \infty\) and, conversely, \(R_0\rightarrow 1\). In the vicinity of a critical point, \(R_O\) and \(\xi/L\) are renormalization-group invariant quantities. Their crossing can be used to locate the onset of the phase transition, rendering them powerful quantities to diagnose the ground-state order and to study phase transitions.

An ergodic QMC simulation averages over all symmetry-breaking states. As a result, we are not able to observe the ordered state directly, but have to refer to correlation functions that do not average out to zero when averaging over all degenerate ground states. Unfortunately, such an approach does not distinguish between the different VBS states illustrated in Figs. 3.7(a)-(c). To obtain additional insights we use the method of a pinning field [103, 104]. In this approach, we explicitly break the symmetry by making one AFM interaction at the origin \(J_{\text{pin}}\sum_a \hat{S}^{(a)}_{\ve{0}} \hat{S}^{(a)}_{\ve{0}+\ve{e}_{\rm x}}\) stronger than the other interactions \(J\sum_a \hat{S}^{(a)}_{\ve{r}} \hat{S}^{(a)}_{\ve{r}+\ve{e}_i}\). The resulting Hamiltonian reads:

(3.32)\[\hat{H} = J \sum_{(\ve{r}, i) \neq (\ve{0}, {\rm x})} \sum_{a} \hat{S}^{(a)}_{\ve{r}} \hat{S}^{(a)}_{\ve{r}+\ve{e}_i} + J_{\text{pin}} \sum_{a} \hat{S}^{(a)}_{\ve{0}} \hat{S}^{(a)}_{\ve{0}+\ve{e}_{\rm x}},\]

with \(J_{\text{pin}} > J\). Therefore, we explicitly choose one of multiple degenerate ground states by pinning the bond \((\ve{0}, \ve{0}+\ve{e}_{\rm x})\) and the bond observable

(3.33)\[B_i(\ve{r}) \equiv \frac{1}{C(N,S)} \Braket{ \sum_a \hat{S}^{(a)}_{\ve{r}} \hat{S}^{(a)}_{\ve{r}+\ve{e}_i}} \qquad i = {\rm x, y}\]

does not vanish as in the unpinned case. We notice that the AFM interactions in the Hamiltonian favors a minimization of \(B_i(\ve{r})\). As proven in Appendix. B.4, \(B_i(\ve{r})\geq -1\) and approaches \(-1\) when the two spins form a singlet.

We have chosen \(J_{\text{pin}}\) such that the pinned bond satisfies

(3.34)\[\begin{equation} \Braket{\sum_a \hat{S}^{(a)}_{\ve{0}} \hat{S}^{(a)}_{\ve{0}+\ve{e}_{\rm x}}} \approx \frac{1}{2} \left( \frac{1}{2 N_r} \sum_{(\ve{r}, i), a} \Braket{\hat{S}^{(a)}_{\ve{r}} \hat{S}^{(a)}_{\ve{r}+\ve{e}_i}} - C(N,S) \right). \end{equation}\]

The right-hand side corresponds to the point halfway between the background and the minimal value of \(\Braket{\sum_a \hat{S}^{(a)}_{\ve{0}} \hat{S}^{(a)}_{\ve{0}+\ve{e}_{\rm x}}}\). This results in \(J_{\text{pin}}\) between \(1.2 J\) and \(1.5 J\).

At large distances from the pinned bond, one will either be able to explicitly observe the “selected” order through \(B_i(\ve{r})\) (if \(\phi\) or \(\psi\) are nonzero), or the order will vanish.

3.4.2. \(S=1/2\)

Fig. 3.8 Correlation functions \(S(\ve{k})\) [Eq. (3.25)] and \(D_{ij}(\ve{k})\) [Eq. (3.25)] for the representation of Fig. 3.2 with \(S=1/2\) and different values of \(N\).

We first study the representations of Fig. 3.2 with \(S=1/2\). In Fig. 3.8, we show the structure factor \(S(\ve{k})\) [Eq. (3.25)], and the two combinations of dimer correlations \(D_{ij}(\ve{k})\) appearing in the definitions of the order parameter \(\phi\) [Eq. (3.28)] and \(\psi\) [Eq. (3.29)]. By considering three values of \(N\), we analyze the evolution of the order parameters across the transition between the Néel and VBS order. For \(N=2\), we realize a standard \(S=1/2\) Heisenberg model on the square lattice, which displays a Néel order in the ground state. As expected, \(S(\ve{k})\) shows a strong peak at \((\pi,\pi)\); in comparison, the other order parameters are suppressed. Upon increasing \(N\) to \(N=4\), we observe the emergence of peaks at \((\pi,0)\) and \((0,\pi)\) for the other order parameters. As shown below using the correlation ratios, though close to the phase transition to the VBS state, the ground state is still Néel ordered. For \(N=6\) we observe a strong peak at \((\pi,0)\) and \((0,\pi)\) for the \(D_{xx}(\ve{k}) + D_{yy}(\ve{k})\) order parameters, indicating the realization of the VBS phase.

Fig. 3.9 Correlation ratios of the staggered magnetization \(m\) [Eq. (3.27)], \(\phi\) [Eq. (3.28)] and \(\psi\) [Eq. (3.29)] order parameters for \(S=1/2\), as a function of \(N\), and for lattice sizes \(L=4-16\).

To obtain a reliable determination of the ground state, we study the correlation ratios of the three order parameters discussed in Sec. 3.4.1. In Fig. 3.9, we show the three correlation ratios \(R_m\), \(R_\phi\), and \(R_\psi\) as a function of \(N\), for lattice sizes \(L=4, 8, 12, 16\). The crossing plot of \(R_m\) indicates the disappearance of Néel order at \(N\approx 4\). At the same time, the curves of \(R_\phi\) and \(R_\psi\) exhibit a crossing for values of \(N>4\); in particular, for \(N=4\) both \(R_\phi\) and \(R_\psi\) decrease with the lattice size, indicating that both VBS and nematic order are short-ranged for \(N=4\). Therefore, as anticipated above, for \(N=4\) the ground state is still a Néel order. This confirms the result of Ref. [105, 106, 107]. At \(N=6\), \(R_m\) decreases with \(L\), while \(R_\phi\) and \(R_\psi\) increases in \(L\), for \(L\ge 6\); the \(L=4\) data set is dominated by finite-size effects. Accordingly, for \(N=6\) the ground state realizes VBS order.

Fig. 3.10 Real-space value of bonds \(B_i(\ve{r})\) [Eq. (3.33)], as measured after pinning the central bond, for \(S=1/2\) and lattice size \(L=18\). Due to the observed different variations in the bond strength, in order to better highlight the patterns of bond correlations we have used different color scales for the region close and far from the pinned bond. The biggest error of the outer bonds is indicated by two red lines on the color scale. We have symmetrized the results with regards to inversions \(y \rightarrow -y\), \(x \rightarrow -x\) around the pinned bond.

These observations are confirmed by the real-space plot of the bond intensity shown in Fig. 3.10, as obtained with the pinning-field method. For \(N=4\), in the vicinity of the pinned bond we observe a pattern reminiscent of the VBS order of Fig. 3.7. However, at larger distances the modulation of bond intensity quickly decays, confirming that VBS order is actually short ranged. For \(N=6\), we observe a very clear pattern of strong and weak bonds that realize the VBS order.

3.4.3. \(S=1\)

Fig. 3.11 Same as Fig. 3.8 for \(S=1\).

In studying the representations with \(S=1\), we proceed analogously to the \(S=1/2\) case discussed in Sec. 3.4.2. In Fig. 3.11, we show the order parameters in momentum space for the case \(S=1\), and three representative values of \(N\). For \(N=8\), a clear peak of the spin structure factor at \((\pi, \pi)\), along with a comparatively smoother momentum dependence of the other order parameters, indicate the presence of Néel order. At \(N=10\), we observe instead the emergence of a clear signal of the nematic order parameter at zero momentum. The VBS order parameter exhibits a similar peak at zero momentum, whose signal predominantly arises from the nematic order parameter, while at momentum \((\pi,0)\) a subdominant peak is observed. The Néel order parameter instead does not show a predominant signal at \((\pi,\pi)\), but rather equally large values at the corners of the Brillouin zone. These behaviors suggest the onset of the two-fold degenerate nematic order. For \(N=12\) a clear signal at \((\pi,0)\) momentum appears in the VBS order parameter. Along with the sharp zero-momentum value of the nematic order parameter, these findings suggest the realization of VBS order for \(N=12\).

Fig. 3.12 Same as Fig. 3.9 for \(S=1\).

As we did for the \(S=1/2\) case, the above qualitative observations on the momentum dependence of the various order parameters can be put on firm ground by examining the correlation ratios shown in Fig. 3.12. The magnetic correlation ratio \(R_m\) displays a crossing at about \(N\approx 8\), such that for \(N>8\), \(R_m\) decreases with the lattice size. On the other hand, at \(N=8\) both \(R_\phi\) and \(R_\psi\) decrease on increasing \(L\), implying that both VBS and nematic order are short ranged. Therefore, one can conclude that for \(N=8\) the ground state is antiferromagnetically ordered. The behavior of \(R_\phi\) shows rather important finite-size corrections. In fact, while curves for \(L\le 10\) cross for \(8 < N < 10\), the crossing point quickly increases with \(L\), such that for \(L\ge 12\) a crossing is found for \(10 < N <12\). In particular, \(N=10\) has a nonmonotonic behavior, increasing in \(L\) for \(L\le 10\), and decreasing for \(L\ge 12\). This observation supports the presence of a significant, but still short-ranged, VBS order, which is responsible for important finite-size corrections. The nematic correlation ratio \(R_\psi\) shows a crossing between \(N=8\) and \(N=10\). Also, here we observe a clear drift in the crossing of \(R_\psi\), although the situation for \(N=10\) is rather clear and indicates long-range order in \(\psi\). In view of these observations, and referring to Table 3.1, we conclude that a nematic ground state is realized for \(N=10\), while for \(N\ge 12\) the ground-state is VBS ordered.

Fig. 3.13 Same as Fig. 3.10 for \(S=1\).

These conclusions are nicely confirmed by the real-space plots of the bond strength obtained with the pinning-field method and shown in Fig. 3.13. For \(N=10\) we clearly observe the formation of a twofold degenerate stripe-like structure, signalling the presence of the nematic order. For \(N=12\) instead a VBS order is found. Interestingly, while for \(S=1/2\) the VBS order found at \(N=6\) (Fig. 3.10) resembles the ladder order illustrated in Fig. 3.7, for \(S=1\) and \(N=12\) the VBS pattern shown in Fig. 3.13 rather suggests the plaquette order of Fig. 3.7.

3.4.4. \(S=3/2\)

Fig. 3.14 Same as Fig. 3.8 for \(S=3/2\).

In Fig. 3.14, we show the order parameters in momentum space, and \(N=12-18\). Analogous to the cases analyzed in the previous sections, for \(N=12\) we find a signal of Néel order. In the region \(14\le N\le 16\), QMC data do not allow us to unambiguously single out the ground state. Upon increasing \(N\), the peak at \((\pi,\pi)\) in \(S(\ve{k})\) slowly decreases in magnitude. At the same time, we observe the appearance of a maximum in the nematic order parameter at zero momentum, and in the VBS order parameter for \((\pi, 0)\) and \((0,\pi)\) momenta. Eventually, for \(N=18\) the momentum structure of the order parameters more clearly favors the realization of VBS order. The observed behavior suggests a comparatively broad critical region around \(14\le N\le 18\).

Fig. 3.15 Same as Fig. 3.9 for \(S=3/2\).

In an attempt to better understand the ground-state diagram for \(S=3/2\), as for the other values of \(S\) we have analyzed the correlation ratios, shown in Fig. 3.15. Due to the increased computational costs, we restricted the simulations for the larger lattice sizes \(L\ge 12\) to the more involved cases \(12\le N\le 16\). For smaller lattice sizes \(L\le 10\), the curves for \(R_m\) appear to cross at a value of \(N\) very close, but smaller than \(N<12\). We observe, however, some drift towards larger values of \(N\) in the crossings. Furthermore, as shown in the inset of Fig. 3.15, \(R_m\) at \(N=12\) exhibits an upwards trend for \(L>10\). Together with the observed slow decrease of \(R_\phi\) and \(R_\psi\) in \(L\) for \(N=12\), this implies Néel order for \(N=12\).

For \(N\ge 14\), the \(L-\)dependence of \(R_m\) clearly rules out Néel order. On the other hand, \(R_\phi\) slowly decreases with \(L\) for \(N=14\), and for \(N=16\) it grows slightly up to \(L=8\). In both cases, QMC data for \(L\ge 8\) are indistinguishable within error bars. A similar flattening of QMC data is found in \(R_\psi\) for \(N=14, 16\). This behavior does not allow us to draw firm conclusions on the nature of the ground state for \(N=14,16\). Since a Néel state can be ruled out, a reasonable hypothesis is the realization of a VBS state, however, with a weak order parameter.

For \(N>16\), both \(R_\phi\) and \(R_\psi\) show a crossing close to \(N=18\). Furthermore, we observe a monotonic growth of \(R_\phi\) in \(L\) for \(N=18\). This leads us to conclude a VBS order for \(N=18\).

Fig. 3.16 Same as Fig. 3.10 for \(S=3/2\) and lattice size \(L=14\) \((N=14,16)\), \(L=12\) \((N=18)\).

Finally, we have studied the pattern of bond strength in real space with the pinning-field method. The results for \(N=14,16,18\) are reported in Fig. 3.16. For \(N=14,16\), despite some signs of dimerization, we do not observe a clear VBS pattern. In line with the previous analysis, for \(N=18\) we find a bond dimerization which confirms a VBS ground state.

3.4.5. \(S=2\)

As for previous values of \(S\), we begin our investigation for \(S=2\) with momentum space plots of correlation functions shown in Fig. 3.17. At \(N=16\), the spin structure factor shows a sharp peak at \((\pi,\pi)\), indicating long-range AFM order. For bigger values of \(N\) the peak weakens and broadens, suggesting short-range AFM order. The correlations for nematic and VBS order show only very broad maxima for the range \(N \in [16, 22]\), implying the absence of both of these orders. The correlation ratios plotted in Fig. 3.18 support these qualitative observations. The inset in Fig. 3.18(a) shows that while \(R_m\) decreases at \(N=16\) from \(L=4\) to \(L=12\), the trend is reversed on bigger lattices and \(R_m\) increases from \(L=12\) to \(L=16\). This indicates that \(N=16\) has a Néel ground state which is close to a competing order. Both \(R_\phi\) and \(R_\psi\) decrease with increasing system size in the investigated range \(N \in [16, 24]\). As per Table 3.1, this leaves a two-dimensional AKLT order as a ground-state candidate for \(N \in [18, 24]\). With the AKLT construction corresponding to the right-hand side of Fig. 3.7(e), we understand that each boundary site hosts an \(\mathfrak{su}(n)\) representation corresponding to one column (\(S=1/2\)) and \(N/2\) rows. Hence the boundary defines a one-dimensional chain in the aforementioned representation. It is known that for \(N \geq 4\) this chain dimerizes [105, 106, 108, 109].

To further investigate this possibility, we simulate the model on a lattice with periodic boundary conditions in the \(x\) direction and open boundary conditions along \(y\), corresponding to a cylinder geometry. Fig. 3.19 shows the results for \(S=2,N=18\) in a pinning-field approach with pinned bonds at the edge and in the bulk, respectively. The induced dimerization pattern propagates on the edge much further than in the bulk, supporting the presence of an AKLT phase with boundary corresponding to an \(S=1/2\) SU(\(N\)) chain in the totally antisymmetric self-adjoint representation.

Fig. 3.17 Same as Fig. 3.8 for \(S=2\)

Fig. 3.18 Same as Fig. 3.9 for \(S=2\)

Fig. 3.19 Real-space value of bonds \(B_i(\ve{r})\) for \(S=2\), \(N=18\), lattice size \(L \in \{8, 10 , 12\}\) and open boundary conditions in \(y\) direction. Comparison between pinning a bond at the edge (a) and a bond in the bulk (b). (c) \(B_{\rm x}(\ve{r})\) on horizontal lines through the pinned bonds.

3.5. Summary

We have studied the ground-state phase diagram of an SU(\(N\)) AFM model on the square lattice, with irreducible representations of \(\mathfrak{su}(N)\) illustrated in Fig. 3.2 and characterized by a Young tableau consisting of \(2S\) columns and \(N/2\) rows. For even values of \(N\) we have presented negative sign free QMC data that results in the rich phase diagram of Fig. 3.1. In line with field-theoretical studies [18, 20, 21], for any value of the generalized spin \(S\), we found Néel order at small values of \(N\), and a dimerized VBS state for large \(N\). The disordered states proximate to the melting of the Néel state can be naturally understood in terms of condensation of monopoles. These states turn out to be located along the \(N=8S + 2\) in the \(S\) versus \(N\) phase diagram. At \(S=1\) we observe nematic AKLT [19] states where \(C_4\) symmetry is spontaneously reduced to \(C_2\), with the emergence of spin-1 chains along one lattice direction. At \(S=2\), the AKLT construction provides an understanding of the non-degenerate state. In fact, this construction can be generalized to any spin \(S=Z/2\) system on a lattice with coordination number \(Z\); an example for the honeycomb lattice is given in Refs. [110, 111]. In our specific case, the edge state corresponds to an SU(\(N\)) spin system in an irreducible representation specified by a Young tableau with one column (\(S=1/2\)) and \(N/2\) rows. At \(N \geq 4\), such a state is known to dimerize [106]. Our simulations with open boundary conditions support this picture. For half-integer values of \(S\), a fourfold degenerate VBS state emerges. The detailed nature of the dimerization was studied using a pinning-field approach [103].

While the resulting phase diagram reproduced in Fig. 3.1 will serve as a benchmark for future studies, our findings points to future avenues of research. It would be very interesting to study in detail the quantum phase transitions between the various states. Although the present setup allows us to consider integer values of \(N\) only, it may be possible to investigate the phase transitions with a suitably defined designer Hamiltonian containing, e.g., some interactions that favor a specific phase, so as to be able to interpolate between them. The topological arguments that lead to the observed phase diagram carry over to other representation of \(\mathfrak{su}(N)\), such that further calculations with alternative methods such as stochastic series expansion [112] are certainly desirable.

---

[1]

In Eq. (3.1) and Eq. (3.8) we have implicitly assumed a choice of the basis of \(\mathfrak{su}(N)\), such that the interaction term is SU(\(N\))-invariant. This condition is satisfied by the basis choice given below in Eq. (3.9).

[2]

We notice that the operator defined in Eq. (3.10) commutes with the algebra only for completely antisymmetric structure constants. For a general choice of the base, one needs to introduce a metric tensor \(g_{ab}\) determined by the structure constants and the Casimir operator is defined as \(\sum_{ab}g^{ab}\Gamma(T_a)\Gamma(T_b)\) [94]; for completely antisymmetric structure constants \(g_{ab}\propto\delta_{ab}\). Also, for the same reason, a normalization is implicit in the definition of Eq. (3.10), discussed in App. E.1.