Spin-$S$ designer hamiltonians and the square lattice $S=1$ Haldane nematic

TL;DR

This paper develops a method to construct spin-$S$ lattice models that are both symmetric and efficiently simulatable, and demonstrates a new phase called Haldane nematic in an $S=1$ square lattice model.

Contribution

It introduces a general strategy for creating Marshall positive, symmetric spin-$S$ Hamiltonians and applies it to realize a novel Haldane nematic phase in an $S=1$ model.

Findings

Discovery of a Haldane nematic phase breaking lattice rotational symmetry

First-order transition between Néel and Haldane nematic phases

Efficient simulation of spin-$S$ models with arbitrary spin using world line Monte Carlo

Abstract

We introduce a strategy to write down lattice models of spin rotational symmetric Hamiltonians with arbitrary spin-$S$ that are Marshall positive and can be simulated efficiently using world line Monte Carlo methods. As an application of our approach we consider a square lattice $S=1$ model for which we design a $3\times 3$ - spin plaquette interaction. By numerical simulations we establish that our model realizes a novel "Haldane nematic" phase that breaks lattice rotational symmetry by the spontaneous formation of Haldane chains, while preserving spin rotations, time reversal and lattice translations. By supplementing our model with a two-spin Heisenberg interaction, we present a study of the transition between N\'eel and Haldane nematic phase, which we find to be of first order.

Tables

Table 1. Table 1: Quantities measures by QMC at low temperatures (inverse temperature of β = 6 ​ L 𝛽 6 𝐿 \beta=6L ) for the Heisenberg Antiferromagnet with coupling constant J = 1 𝐽 1 J=1 compared with those determined for the ground state of the same model as found from Exact Diagonalization (ED). The energies reported ( e E ​ D subscript 𝑒 𝐸 𝐷 e_{ED} and e Q ​ M ​ C subscript 𝑒 𝑄 𝑀 𝐶 e_{QMC} ) are per site and the stiffness ( ρ s E ​ D subscript superscript 𝜌 𝐸 𝐷 𝑠 \rho^{ED}_{s} and ρ s Q ​ M ​ C subscript superscript 𝜌 𝑄 𝑀 𝐶 𝑠 \rho^{QMC}_{s} ) are as described by Eq. 21 . The energy saturates to the ground state value for low enough temperatures as can be seen in Fig. 13

Lattice	Lx	Ly	S	$e_{ED}$	$e_{QMC}$	${\rho }_s^{ED}$	${\rho }_s^{QMC}$
Square	4	4	$\frac{1}{2}$	-1.20178	-1.20186(8)	0.1855	0.1849(4)
Square	2	2	1	-5.0	-5.0011(5)	1.0	1.002(2)
Square	2	2	$\frac{3}{2}$	-10.5	-10.499(5)	2.0	2.008(3)
Chain	4	1	1	-2.5	-2.5004(2)	0.2222	0.2226(4)
Chain	6	1	$\frac{3}{2}$	-5.1488	-5.1489(3)	0.2630	0.2627(3)

Equations

$${\mathbf{S}}_i=\sum _a{\mathbf{s}}_i^a.$$

$$H_2^{ij}={\mathbf{S}}_i\cdot {\mathbf{S}}_j-1=-\sum _{a,b}(\frac{1}{4}-{\mathbf{s}}_i^a\cdot {\mathbf{s}}_j^b)$$

$$H_3^{ijk}$$

$$+$$

$$H_{3\times 3}^p=H_3^{123}{H}_3^{456}{H}_3^{789}+H_3^{147}{H}_3^{258}{H}_3^{369},$$

$$H=J\sum _{⟨ij⟩}{H}_2^{ij}+Q_{3\times 3}\sum _p{H}_{3\times 3}^p.$$

$${\vec{S}}_i=\sum _{\mu =1}^{2S}{\vec{s}}_i^a$$

$$Z=T{r}_s(e^{-\beta \tilde{H}}\mathcal{P})$$

$$\mathcal{P}=\prod _i{\mathcal{P}}_i$$

$$\mathcal{P}_{i}=|\uparrow\uparrow\rangle\langle\uparrow\uparrow|+|\downarrow\downarrow\rangle\langle\downarrow\downarrow|\\ +\Big{(}\frac{|\uparrow\downarrow\rangle+|\downarrow\uparrow\rangle}{\sqrt{2}}\Big{)}\Big{(}\frac{\langle\uparrow\downarrow|+\langle\downarrow\uparrow|}{\sqrt{2}}\Big{)}$$

$$\begin{array}{rl}H& =-J\sum _{⟨ij⟩}(S^2-{\vec{S}}_i.{\vec{S}}_j)\\ & =J\sum _{⟨ij⟩}{H}_{ij}\end{array}$$

$$H_{ij}={\vec{S}}_i.{\vec{S}}_j-S^2$$

$$\tilde{H}_{ij}=-J\sum_{a,b=1}^{2S}\Big{(}\frac{\mathbbm{1}}{4}-\vec{s}_{i}^{\,a}.\vec{s}_{j}^{\,b}\Big{)}$$

$$H_{biquad}=({\vec{S}}_i.{\vec{S}}_j{)}^2-\text{mathbbm}1$$

$$\tilde{H}_{biquad}=\sum_{a,b=1}^{2}\Big{[}\Big{(}\frac{\mathbbm{1}}{4}-\vec{s}^{\,a}_{i}.\vec{s}^{\,a}_{j}\Big{)}.\Big{(}\frac{\mathbbm{1}}{4}-\vec{s}^{\,b}_{i}.\vec{s}^{\,b}_{j}\Big{)}+\\ \Big{(}\frac{\mathbbm{\mathbbm{1}}}{4}-\vec{s}^{\,a}_{i}.\vec{s}^{\,b}_{j}\Big{)}.\Big{(}\frac{\mathbbm{1}}{4}-\vec{s}^{\,b}_{i}.\vec{s}^{\,a}_{j}\Big{)}+\\ {\color[rgb]{1,0,0}\frac{1}{2}\Big{(}\vec{s}^{\,a}_{i}.\vec{s}^{\,b}_{i}-\frac{\mathbbm{1}}{4}\Big{)}+\frac{1}{2}\Big{(}\vec{s}^{\,a}_{j}.\vec{s}^{\,b}_{j}-\frac{\mathbbm{1}}{4}\Big{)}}\Big{]}$$

$$\tilde{H}^{ijk}_{3}=-\sum_{a,d,b\neq c}\Big{(}\frac{\mathbbm{1}}{4}-\vec{s}^{\,a}_{i}.\vec{s}^{\,b}_{j}\Big{)}\Big{(}\frac{\mathbbm{1}}{4}-\vec{s}^{\,c}_{j}.\vec{s}^{\,d}_{k}\Big{)}+h.c.$$

$$\begin{split}\tilde{H}^{ijk}_{3}&=-\sum_{a,d,b\neq c}\Big{(}\frac{\mathbbm{1}}{4}-\vec{s}^{\,a}_{i}.\vec{s}^{\,b}_{j}\Big{)}\Big{(}\frac{\mathbbm{1}}{4}-\vec{s}^{\,c}_{j}.\vec{s}^{\,d}_{k}\Big{)}\\ &=-\Big{[}\sum_{a,b,c,d}\Big{(}\frac{\mathbbm{1}}{4}-\vec{s}^{\,a}_{i}.\vec{s}^{\,b}_{j}\Big{)}\Big{(}\frac{\mathbbm{1}}{4}-\vec{s}^{\,c}_{j}.\vec{s}^{\,d}_{k}\Big{)}\\ &-\sum_{a,b,c}\Big{(}\frac{\mathbbm{1}}{4}-\vec{s}^{\,a}_{i}.\vec{s}^{\,b}_{j}\Big{)}\Big{(}\frac{\mathbbm{1}}{4}-\vec{s}^{\,b}_{j}.\vec{s}^{\,c}_{k}\Big{)}\Big{]}\end{split}$$

$$\begin{array}{rl}{H}_3^{d,ijk}& =-(\text{mathbbm}1-{\vec{S}}_i.{\vec{S}}_j)(\text{mathbbm}1-{\vec{S}}_j.{\vec{S}}_k)+h.c.\\ H_3^{d^{\prime },ijk}& =({\vec{S}}_i.{\vec{S}}_j+{\vec{S}}_j.{\vec{S}}_k-{\vec{S}}_i.{\vec{S}}_k-\text{mathbbm}1)\end{array}$$

$$\begin{array}{rl}{H}_3^{ijk}& =H_3^{d,ijk}+H_3^{d^{\prime },ijk}\\ & =-{\vec{S}}_i.{\vec{S}}_j{\vec{S}}_j.{\vec{S}}_k+\frac{1}{2}({\vec{S}}_i.{\vec{S}}_j+{\vec{S}}_i.{\vec{S}}_k+{\vec{S}}_j.{\vec{S}}_k+2)+h.c.\end{array}$$

$$\rho_{s}=\frac{\partial^{2}E(\phi)}{\partial\phi^{2}}\bigg{|}_{\phi=0}\\$$

$${\rho }_s=\frac{⟨W^2⟩}{\beta }$$

$$B_{i}(\vec{r})=\sum_{a,b}\Big{(}\frac{1}{4}-\vec{s}^{\,a}_{\vec{r}}.\vec{s}^{\,b}_{\vec{r}+i}\Big{)}\\$$

$$H_{2\times 3}^p=H_3^{123}{H}_3^{456}+H_3^{147}{H}_3^{258}$$

$$\begin{gathered} H_{Q_3}^p=(1-{\vec{S}}_1.{\vec{S}}_2)(1-{\vec{S}}_4.{\vec{S}}_5)(1-{\vec{S}}_7.{\vec{S}}_8) \\ +(1-{\vec{S}}_1.{\vec{S}}_4)(1-{\vec{S}}_2.{\vec{S}}_5)(1-{\vec{S}}_3.{\vec{S}}_6) \end{gathered}$$

Full text

Introduction: The relationship between lattice spin models and their long distance descriptions by quantum field theories is a central topic in theoretical condensed matter physics sachdev1999:qpt; fradkin2013:book. Pioneering work on the ground state of spin chains found a striking role is played by the size of the quantum spin haldane1988:berry; affleck1989:lh: while half integer spins generically realize a gapless critical phase, integer spin chains realize a topological “Haldane phase”. In the field theoretic understanding, the value of the microscopic value of the spin enters as a co-efficient of a topological term that has a dramatic effect on the spin chain phase diagram. Given this profound result in one dimension, it is natural to ask how the value of the spin-$S$ affects the phase diagrams of two dimensional quantum spin systems?

For one dimensional systems, progress in our understanding is largely due to the availability of specialized analytic sutherland2004:bm; giamarchi2004:one and numerical methods white1992:dmrg. These methods cannot be extended as effectively to two dimensions, where consequently much less is known despite intense research. The most reliable unbiased method to study field theory and quantum criticality in two dimensions are limited to models that do not suffer from the sign problem of quantum Monte Carlo kaul2013:qmc. Although the sign-free condition is very restrictive, given their unique ability to provide unbiased insight it is of great interest to build a repertoire of sign-free spin models for arbitrary spin-$S$, as has been achieved for $S=1/2$ kaul2014:design.

In this Letter we develop a systematic method to write down a large family of sign-free bipartite spin models with arbitrary spin-$S$ and multi-spin interactions that have the Heisenberg rotational symmetry. These new models open the door to study a variety of new phases and phase transitions, many of which are of great interest to the community. As a first application of our method we design a square lattice $S=1$ interaction that realizes a long anticipated “Haldane nematic” (HN) phase affleck1988:aklt; read1990:vbs. In this phase the spin system breaks lattice rotation symmetry but preserves lattice translations due to the spontaneous formation of Haldane chains either in the $x$ or $y$ direction with an associated two-fold ground state degeneracy, Fig. 1(a). Motivated in part by the Iron superconductors the HN phase has been under intense study recently (see e.g. chen2018:s1; jiang2009:s1; gong2017:s1; bilbao2015:nem; corboz). An influential work wang2015:s1 found an exactly solvable model which realizes the HN as a ground state and provided field theoretic arguments for an exotic continuous phase transition to a Néel state described by the O(4) $\sigma$-model at $\Theta=\pi$. We establish unambiguously the existence of the HN phase in our new sign free model and provide the first unbiased numerical study of the phase transition from the HN to the Néel state. We find clear evidence that the transition is first order and discuss the implications of this finding for the field theoretic scenario.

Designer Models: While it is well known that the bipartite Heisenberg model is Marshall positive for arbitrary spin-$S$, what are the most general multi-site spin-$S$ Hamiltonian operators that are sign positive? This question has been difficult to address previously because it appears daunting directly in the language of spin-$S$ operators. Following previous work todo2001:highs; kawashima1994:spinS; kawashima2004:review we take a different route – we rewrite the spin-$S$ on each of the $N_{s}$ lattice sites as $2S$ spin-1/2 “mini-spins”,

[TABLE]

We note here that the ${\bf s}^{a}_{i}$ have both a lattice index $i$ ($1\leq i\leq N_{s}$) and a mini-spin index $a$ ($1\leq a\leq 2S$), giving a total of $2SN_{s}$ mini-spins. To faithfully simulate the original problem, we have to include a projection operator, ${\cal P}=\prod_{i}{\cal P}_{i}$, where ${\cal P}_{i}$ projects out the spin-$S$ from the ${\bf s}^{a}_{i}$ basis, $Z={\rm Tr}_{\bf S}\left[e^{-\beta H({\bf S})}\right]={\rm Tr}_{\bf s}\left[e^{-\beta H({\bf s})}{\cal P}\right]$. Since ${\cal P}$ is itself sign-problem free, in the world-line approach, any model which is sign-free in the ${\bf s}^{a}_{i}$ basis gives us a sign-free spin-$S$ model!

In this manuscript we illustrate our idea using $S=1$ spins on the square lattice, but our results can be straightforwardly extended to any bipartite lattice with arbitrary spin-$S$. Consider first in the ${\bf s}$ language the $S=1$ Heisenberg model,

[TABLE]

Diagramatically we can represent each $\frac{1}{4}-{\bf s}^{a}_{i}\cdot{\bf s}^{b}_{j}$ term in the sum in the last expression as an “${\bf s}$-bond” between mini-spins $a$ and $b$ on the two sites $i$ and $j$. A representative such term is illustrated for $S=1$ with two mini-spins per site in Fig. 1(b) (there are three other such diagrams corresponding to the sum on $a,b$). Likewise, it is easy to see that the interaction with two ${\bf s}$-bonds between $i$ and $j$ corresponds to the sign free region of the biquadratic interaction, Fig. 1(c)  harada2001:bq; spsdes2019:supmat. From these examples, we make our central observation – it is much easier to write down a sign free model in the ${\bf s}$ language than directly in the spin-$S$ basis. As a non-trivial example consider interactions between three $S=1$ spins in a row. In the ${\bf s}$-bond language the most natural interaction is with a single bond between each pair of neighbors without allowing them to touch on the middle site, Fig. 1(d). Working backwards we then find this new sign-free interaction in terms of the spin-1 operators is,

[TABLE]

For $S=1$ models the three-site interaction and its physical significance has been discussed recently michaud2013:s1; chepiga2016:s1. Here we discover that in order to study such terms in a sign free way we have to include two spin terms to balance the signs. Intuitively, the three spin interaction in Fig. 1(d) is reminiscent of the famous AKLT construction affleck1988:aklt and so we can expect it to force our system into a Haldane like phase; we confirm this below. Using the three-site interaction $H^{ijk}_{3}$, we introduce a model interaction we will study in detail below. Following the idea of the J-Q model sandvik2007:deconf we construct a $3\times 3$ plaquette interaction from $H^{ijk}_{3}$,

[TABLE]

The indexing of the sites in the plaquette by numbers 1-9 is shown in Fig. 1(e). The two terms are included to preserve square lattice symmetry 111$3\times 1$ and $3\times 2$ plaquettes interactions can also be considered, but they are found to be unsuitable for the application, i.e. they are insufficient to destroy N’eel order spsdes2019:supmat. A $S=1$ version of the $Q_{3}$ term  sandvik-kawashima is also found unable to destabilize N’eel order. .

We emphasize that in addition to the advantage of leading us to new non-trivial sign free interactions, the mini-spin representation also offers us a simple way to construct efficient loop update algorithms for complex interactions such as Eq. (5), since we can update the ${\bf s}$ interactions using the standard deterministic algorithm using for e.g. the stochastic series expansion sandvik2010:vietri. The update of the symmetrization operator is straightforward using the directed loop algorithm  spsdes2019:supmat; directed-loops. Clearly this program of designing sign-free interactions in terms of the ${\bf s}$-bond diagrammatic representation and then into the spin operators can be extended systematically to any value of spin-$S$ and to a wide range of multi-spin interactions. Rather than elaborate on this here, we now turn to an application.

Haldane Nematic: We consider square lattice $S=1$ antiferromagnets, which have been argued to host an exotic “Haldane nematic” (HN) state in their phase diagrams. Our goal here is to establish that the sign-free model, Eq. 5 realizes this novel phase and carry our unbiased studies of the phase transitions of the destruction of HN order.

The model we study is,

[TABLE]

The first term is the usual square lattice $S=1$ Heisenberg model. The second term is our new designer interaction with a sum on $p$, which runs over the elementary $3\times 3$ plaquettes on the square lattice. We study the phase diagram as a function of $g\equiv Q_{3\times 3}/J$ and the temperature $T=1/\beta$. We work in units in which $J^{2}+Q^{2}_{3\times 3}=1$. The phase diagram inferred from our simulations is shown in Fig. 2. At $(g,T)=(0,0)$ (labelled as H) our model is the nearest neighbor $S=1$ Heisenberg model which is Néel ordered singh1990:s1. We use the conventional order parameter $\langle m^{2}\rangle$ with $m=\sum_{\bf r}e^{i(\pi,\pi)\cdot{\bf r}}S^{z}_{\bf r}/N_{s}$ to diagnose long range magnetic order. From the finite size scaling of $\langle m^{2}\rangle$ we observe that the Néel order weakens as $g$ is increased (I). At $T=0$ the Néel order is stable until we reach a coupling $g\approx 0.17$ at which Néel order is destroyed. As is well known, the Néel order cannot survive finite-$T$ Mermin-Wagner fluctuations in two dimensions.

We now present extensive numerical evidence that at $T=0$ for $g\geq 0.17$ the system transitions into the “Haldane nematic” phase (Fig. 1(a)). We first rule out a conventional VBS pattern where pairs of $S=1$ dimerize into a columnar pattern wildeboer2018:s1, which can be studied by finite size scaling of $\langle\phi^{2}\rangle$ with $\phi=\sum_{\bf r}e^{i(\pi,0)\cdot{\bf r}}B_{x}({\bf r})/N_{s}$ [with the bond operator $B_{i}({\bf r})\equiv J{\bf S}_{\bf r}\cdot{\bf S}_{{\bf r}+{\bf e}_{i}}$]. As shown in the inset of Fig. 3 $\langle\phi^{2}\rangle$ scales to zero in the thermodynamic limit indicating that in all parts of the phase diagram under study the conventional VBS order is absent. We use an order parameter okubo2015:sun $\langle\psi^{2}\rangle$ that is sensitive to breaking of rotational symmetry without picking up signals of translational symmetry breaking. $\psi=\sum_{\bf r}(B_{x}({\bf r})-B_{y}({\bf r}))/N_{s}$. Clearly a condensation of $\psi$ indicates the breaking of lattice rotational symmetry. As shown in Fig. 3 K and J clearly have long range HN order, whereas at the other points they are absent either because of Néel order (H and I) or thermal disorder (L).

We now turn to a study of the phase transition at which HN order is destroyed. We begin by simulating the model at $g=0.5$ and tuning $T$ along the vertical dashed line in Fig. 2. From Fig 3, as we move from L (no HN order) to K (HN order) to J (stronger HN order) we have clear evidence for a phase transition. If the pattern of symmetry breaking is of the form Fig. 1(a) thermal criticality is expected to be of the Ising universality class. In Fig. 4 we present a study of the histograms of the order parameter. We see that just above the critical $T$, $P(\psi)$ shows one peak at zero. As $T$ is lowered, the zero-peak splits into two symmetric peaks corresponding to spontaneous symmetry breaking, just as one expects for the Ising model. There is no evidence for a peak at zero co-existing with the non-zero peaks, which one would expect at a first order transition. A study of the scaling behavior of the $T$-dependence of the order parameter at $g=0.5$ (right panel of Fig. 4) shows conclusive evidence that the HN order parameter undergoes a continuous thermal Ising phase transition, as expected for its order parameter manifold. This provides our final piece of evidence that the broken symmetry is indeed of the Haldane nematic form illustrated in Fig. 1(a).

A final interesting question we address is the nature of the quantum phase transition between Néel-HN, labeled by a star in Fig. 2. The field theory for this phase transition has been argued to be the O(4) $\sigma$-model at topological angle $\pi$ wang2015:s1, building on previous work for $S=1/2$ tanaka2005:nvbs; senthil2006:topnlsm. Very little is known about this field theory, but a consistent scenario for a continuous transition with emergent $O(4)$ symmetry at the critical point would require only one relevant $O(3)\times Z_{2}$ anisotropy that appears as the tuning parameter $g$ in the lattice model. This delicate question has not yet been accessed in unbiased simulations. To approach this point we study the nature of the phase transition as we move down the thermal phase transition line to lower temperatures. From Fig. 4, we have seen at high-$T$ the transition is continuous and of the Ising type. In Fig. 5 we study data at $T=1$ (which is very low-T in the units in which we are working) while tuning $g$ (horizontal dashed line in Fig. 2). The histogram data shows clear evidence that the transition has become first order for the HN order parameter, with a co-existence of a peak at zero (for non-HN phase) and the finite symmetry related peaks for the HN phase. While there is no thermal phase transition for the Néel order it also shows double peaks that are incipient behavior of the first order quantum phase transition it undergoes at $g\approx 0.17$. We thus reach the conclusion that along the phase boundary line (solid curve in Fig. 2) the phase transition changes from being Ising and continuous at high-$T$ to becoming first order at low-$T$ and remains first order at the quantum phase transition, marked with a star. The change from continuous Ising to first order is expected to happen at a multi-critical point somewhere along the solid line in Fig. 2 between the two limiting cases we have studied and is expected to be described by the tricritical Ising field theory cardy1996:ren. We have not made an effort to locate this point precisely in our phase diagram in this work.

Our finding of a first order quantum transition can be interpreted in two different ways for the $O(4)$ sigma-model at $\theta=\pi$. The first is simply that the field theory itself does not have a non-trivial critical fixed point, the other is that such a fixed point exists but it has more than one relevant $O(3)\times Z_{2}$ anisotropy and thus requires more than one tuning parameter to be reached. We note that our finding is consistent with previous studies of the $S=1/2$ Néel-VBS deconfined critical point on a rectangular lattice which is expected to be described by the same field theory and anisotropies as the $S=1$ Néel-HN studied here haldane1988:berry; read1990:vbs; senthil2004:science and was also found to be first order block2013:fate.

Conclusions: We have introduced a scheme to design general multi-spin interactions for spin-$S$ models without the sign problem. Our scheme opens up the possibility to simulate a wide range of models and address the role of $S$ on quantum phase transitions in two and higher dimensions. Higher spin can introduce new phases not present for $S=1/2$, including multi-polar ordered phases and new paramagnetic phases, like the unconventional valence bond ordering we found here and quantum spin liquids. The theory of phase transitions between these new phases is largely unexplored. All of these are exciting avenues for future work.

Acknowledgments: We gratefully acknowledge useful discussion with S. Pujari and partial support from NSF DMR-1611161 and Keith B. MacAdam Graduate Excellence Fellowship. The numerical results were produced on SDSC comet cluster through the NSF supported XSEDE award TG-DMR140061 as well as the DLX cluster at UK.