Stochastic series expansion algorithm for the $S=1∕2$ $XY$ model with four-site ring exchange

We describe a stochastic series expansion quantum Monte Carlo method for a two-dimensional $S=1∕2$ $XY$ model (or, equivalently, hard-core bosons at half filling) which in addition to the standard pair interaction $J$ includes a four-site term $K$ that flips spins on a square plaquette. The model has three ordered ground state phases; for $K∕J\lesssim 8$ it has long-range $xy$ spin order (superfluid bosons), for $K∕J\gtrsim 15$ it has staggered spin order in the $z$ direction (charge-density wave), and between these phases it is in a state with columnar order in the bond and plaquette energy densities. We discuss an implementation of directed-loop updates for the SSE simulations of this model and also introduce a “multibranch” cluster update which significantly reduces the autocorrelation times for large $K∕J$. In addition to the pure $J\text{−}K$ model, which in the $z$ basis has only off-diagonal terms, we also discuss modifications of the algorithm needed when various diagonal interactions are included.

Figure 1

(a) Labeling convention for the indices of an operator $O_{ijkl}$ acting on the corners of a plaquette. The label $p$ refers to the whole plaquette, so that $O_p\equiv O_{ijkl}$. (b) The two plaquette configurations between which the $K$ term can act; open and solid circles correspond to up and down spins, respectively.

Figure 2

Examples of vertices for the $J\text{−}K$ model. The solid and open circles correspond to up and down spins, respectively, before (below the bar) and after (above the bar) an operator has acted. The spins on the plaquette are here drawn on a line in order to allow for a more convenient graphical representation. The order of the sites correspond to the indices $i,j,k,l$ in Fig. 1 from left to right. (a) is one out of 16 diagonal $C$ vertices, (b) is one out of 32 $J$ vertices (which flip two spins), and (c) is one of the two $K$ vertices (which flip all four spins).

Figure 3

The linked-vertex representation corresponding to the matrix element $⟨\alpha \mid H_{6,2}{H}_{4,1}{H}_{4,2}\mid \alpha ⟩$ (left), where the basis state $\mid \alpha ⟩=\mid \downarrow \downarrow \downarrow \uparrow \uparrow \downarrow ⟩$ (right). The bidirectional links are represented as dashed lines. The site and plaquette numbering for the six-site lattice is shown to the right. The numbers at the vertex legs indicate links across the periodic propagation boundary, and the corresponding spins are hence the same as in the state $\mid \alpha ⟩$. The numbers here also correspond to the site numbering of the lattice shown to the right.

Figure 4

Two examples of vertex paths that are related in the directed-loop scheme. (a) shows a process, and its reverse, where a $C$ vertex is transformed into a $J$ vertex. (b) shows two related $J\leftrightarrow K$ transformations. In the loop construction the spin states at the entrance and exit legs are flipped. The spin states shown in the vertices here are those before the flips have been carried out.

Figure 5

A closed set of $C$ and $J$ vertex paths, with their corresponding weights $a_{ij}$ that have to satisfy the directed-loop equations.

Figure 6

A closed set of vertex paths with $J\leftrightarrow K$ transformations, labeled by their path weights $b_{ij}$.

Figure 7

A closed set of two vertex paths, where only the continue-straight process is allowed and a $J$ vertex is transformed into another $J$ vertex with probability 1.

Figure 8

Vertex transformations in the multibranch cluster update. The entrance leg is denoted by an arrow pointing into the vertices on the left. In the updated vertices to the right, the spins at the outgoing arrows have been flipped. The branching for all entrance legs and vertices not shown here are obtained by applying trivial symmetries to one of the cases shown in (a)–(i).

Figure 9

Multibranch cluster update in which two $C$ vertices are transformed into two $K$ vertices. The initial entrance leg is at the inward pointing arrow in the linked-vertex representation to the left. The resulting vertices with their arrows indicating legs visited are shown to the right.

Figure 10

Probability histogram of cluster sizes (measured as the number of vertex legs in a cluster) produced by multibranch cluster updates on a $16\times 16$ lattice at $K∕J=80$ and $\beta =32$. The average length of the operator list for this simulation was $⟨n⟩\approx 4.35\times {10}^5$ (i.e., approximately $3.5\times {10}^6$ vertex legs). Data in the upper figure are for the smaller bins, while data in the lower figure were rebinned to 100 legs per bin. All intermediate occupations were measured as zero.

Figure 11

Autocorrelation function for the spin stiffness (upper panel) and the plaquette-stripe order parameter (lower panel) in simulations of $L=16$ systems at $K∕J=4$ and 7, both at inverse temperature $K∕T=16$. Results of simulations both with and without the multibranch cluster update are shown.

Figure 12

Autocorrelation function for the plaquette-stripe order parameter at $K∕J=12$ and inverse temperature $K∕T=32$. Results with (solid curves) and without (dashed curves) multibranch clusters are compared for three different system sizes.

Figure 13

Autocorrelation function for the staggered order parameter at $K∕J=32$ and inverse temperature $K∕T=32$. Results with (solid curves) and without (dashed curves) multibranch clusters are compared for three different system sizes.

Figure 14

The staggered magnetization vs the propagation index $p$ (divided by the total number of operators $n$) in configurations generated for system sizes $L=8,16$, and 32, at $K∕J=32$, $K∕T=32$.

Figure 15

Distribution of the staggered magnetization at $K∕J=32$, $K∕T=32$.

Figure 16

Bin averages for the superfluid stiffness, the stripe structure factor, and the stripe susceptibility, for an $L=96$ lattice at $K∕J=7.91$ and inverse temperature $K∕T=64$. Each point represents an average over ${10}^4$ Monte Carlo steps.

Figure 17

Integrated autocorrelation times, Eq. (44), for the superfluid stiffness (bottom) and the squared plaquette structure factor (top) at $K∕J=7.91$. The solid lines show $\tau \propto L^2$. Error bars are not shown but are significant for the largest three lattice sizes.

Figure 18

A closed set of $C$ and $J$ vertex paths for which the corresponding directed-loop equations are different for $V=0$ and $V\ne 0$. A second set of vertex paths—the spin reverse of the one shown here—also gives the same directed-loop solution.

Figure 19

A closed set of $C$ and $J$ vertex paths used in solving the directed-loop equations for the uniform-field $J\text{−}K\text{−}h$ model.

Figure 20

Sublattice decoration of the two-dimensional square lattice, used in constructing the quantum Monte Carlo simulation of the $J\text{−}K$ model with staggered Zeeman field.

Figure 21

Schematic representation of the closed set of $C$ and $J$ vertex paths used in solving the directed-loop equations for the $J\text{−}K$ model. Blocks with the same shading represent equivalent vertex weights.

Figure 22

Reference elements of the closed $C$ and $J$ vertex-path diagrams. The element letters (left) and equation set designation (right) refer to the relevant vertex weight equation set in Tables IV–IX.