<figure><img src="image_1.png" /></figure>

HAL
open science

# Clustering of Topological Charges in a Kagome Classical Spin Liquid

Tomonari Mizoguchi, L.D.C. Jaubert, Masafumi Udagawa

► To cite this version:

Tomonari Mizoguchi, L.D.C. Jaubert, Masafumi Udagawa. Clustering of Topological Charges in a Kagome Classical Spin Liquid. Physical Review Letters, 2017, 119, pp.077207. 10.1103/Phys-RevLett.119.077207. hal-01549582

HAL Id: hal-01549582

https://hal.science/hal-01549582

Submitted on 28 Jun 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Clustering of Topological Charges in Kagome Classical Spin Liquid

Tomonari Mizoguchi,¹ L. D. C. Jaubert,²,³ and Masafumi Udagawa⁴

¹Department of Physics, University of Tokyo, Hongo, Bunkyo-ku, Tokyo 113-0033, Japan*

²Okinawa Institute of Science and Technology Graduate University, Onna-son, Okinawa 904-0495, Japan

³CNRS, Univ. Bordeaux, LOMA, UMR 5798, F-33400 Talence, France

⁴Department of Physics, Gakushuin University, Mejiro, Toshima-ku, Tokyo 171-8588, Japan

(Dated: February 14, 2017)

Fractionalization is an ubiquitous phenomenon in topological states of matter. In this work, we study the collective behavior of fractionalized topological charges and their instabilities, through the $J_1-J_2-J_3$ Ising model on a kagome lattice, which can be mapped to a model of interacting topological charges under the constraint of Gauss' law. We find that the recombination of topological charges gives rise to a yet unexplored classical spin liquid. This spin liquid is characterized by an extensive residual entropy, as well as the formation of hexamers of same-sign topological charges. The emergence of hexamers is reflected to a half-moon signal in the magnetic structure factor, which provides us a signature of this new spin liquid in neutron-scattering experiments. To study this phase, a worm algorithm has been developed which does not require the usual divergence-free condition.

PACS numbers: 75.10.Kt

Fractionalization is a hallmark of topological states of matter. In these systems, an excitation with a unit quantum number, such as charge and spin, is fractionalized into several constituents. These excitations can then condense into exotic topological phases [1]. The nature of fractionalized excitations have been studied through a number of systems, such as quasi-one-dimensional conducting polymers [2], fractional quantum Hall systems [3, 4] and quantum spin liquids (QSLs) in one and higher dimensions [5].

Among the systems showing fractionalization, QSLs are of special interest. QSL is a long-range-entangled quantum ground state without spontaneous symmetry breaking; the ground state is expressed as a superposition of a macroscopic number of product states. The realization of QSL has been intensively sought in frustrated magnets, and a number of candidate materials have been explored actively [6–13].

Quantum spin ice [14, 15] is one of the most promising systems in this context. Its classical counterpart, spin ice, is a classical spin liquid (CSL) with macroscopically degenerate ground states. CSL often serves as a constituting source of QSL, thanks to quantum fluctuations inducing a superposition between the degenerate ground states. This is why the parent CSL phase reflects several important properties of its descendant QSL state. In the case of quantum spin ice, a fractional excitation, called monopole, can be found in its classical counterpart, in which a flipped spin from the ground state is fractionalized into two half-unit charges [16, 17]. CSL, due to its simplicity compared with QSL, enables us to study rather intractable aspects of fractional excitations.

In this work, we study the cooperative phenomena of fractionalized topological charges in the kagome CSL [18–21] of the $J_1-J_2-J_3$ Ising model. The nearest neighbor (n.n.) interactions ($J_1$) alone lead to a CSL phase, composed of the constrained configurations of topological charges, analogous to the fractional monopole excitations in spin ice. The farther-neighbor interactions ($J_2, J_3$) introduce a n.n. interaction between the topological charges. The resulting phase diagram supports a variety of (dis-)ordered phases, including what is,

to the best of our knowledge, a yet unexplored CSL. This CSL is composed of coexisting hexamer clusters with an extensive residual entropy. Its anomalous spatial structure can be detected through a half-moon signal in magnetic structure factor.

Model.- We consider the kagome lattice composed of $N_{site} = L \times L \times 3$ sites under periodic boundary conditions, and Ising spin variables, $\sigma_i^z = \pm 1$ at each site, i. Throughout this Letter, we set $k_B = 1$. We define the Hamiltonian:

$$ \mathcal{H} = J_1 \sum_{\langle i,j \rangle_{\text{n.n.}}} \sigma_i^z \sigma_j^z + J_2 \sum_{\langle i,j \rangle_{\text{2nd}}} \sigma_i^z \sigma_j^z + J_3 \sum_{\langle i,j \rangle_{\text{3rd}}} \sigma_i^z \sigma_j^z, \quad (1) $$

whose couplings are illustrated in Fig. 1(a). Here, we assume the n.n. coupling as antiferromagnetic (AF), and set its value as a unit of energy, $J_1 = 1$. We also introduce the second ($J_2$) and the third ($J_3$) neighbor terms. As opposed to a series of recent works on kagome [22–26], we shall not consider the $J_3$ term across hexagons, but rather the $J_3$ term between neighboring triangles which offers, as discussed below, an elegant representation in terms of topological charges. For general combinations of $J_2$ and $J_3$, various magnetically ordered states have been studied [20, 21]. In the present work, we focus on the case of $J_2 = J_3(\equiv J)$, where a charge representation is available [27]. This model has offered a rich out-of-equilibrium physics on the three-dimensional pyrochlore lattice [28, 29]. To introduce this representation, we classify the triangles on the kagome lattice into “upward” ($\Delta$) and “downward” ($\nabla$), according to their orientations. We then define a topological charge at each triangle as $Q_p = \eta_p \sum_{i \in p} \sigma_i^z$ with $\eta_p = +1(-1)$ for $p \in \Delta(\nabla)$ [see Fig. 1(b)], in a manner analogous to the monopole in spin ice [16]. In fact, the charge has a fractional character; flipping a single spin changes $Q_p$ by $\pm 2$ on the pair of triangles the spin belongs to, meaning that the change of spin quantum number is fractionalized into a pair of positive and negative charges.

In terms of $Q_p$, the Hamiltonian can be rewritten as follows,
up to a constant term [27]

$$ \mathcal{H} = \left(\frac{1}{2} - J\right) \sum_p Q_p^2 - J \sum_{\langle p, q \rangle} Q_p Q_q. \quad (2) $$

The first term represents the self-energy of a charge, summed over all triangles $p$, while the second term is the interaction between neighboring pairs of charges on upward and downward triangles $\langle p, q \rangle$. Opposite charges attract for $J < 0$, while they repel for $J > 0$.

Since the charges are defined on a bipartite honeycomb lattice connecting the centers of triangles, at first sight, the ground state of (2) seems to show uniform (staggered) charge order for $J > 0$ ($< 0$), trivially. However, possible configurations of $Q_p$ are constrained by the underlying spin structure. The constraint is expressed by a discrete analogue of Gauss’

<figure><img src="image_3.png" /><figcaption>FIG. 1. (Color online) (a) Schematic picture of the kagome lattice and spin interactions. Red, blue and green lines are the n.n., the second neighbor, and the third-neighbor interactions, respectively. Note that the coupling across a hexagon, as shown with a dashed line, is not considered. The unit vectors of a kagome lattice, $\mathbf{a}_1$ and $\mathbf{a}_2$ are represented by the orange arrows. (b) Possible charge states on triangles. The red (blue) dots represent $\sigma^z = \uparrow$ ($\downarrow$). The colors of triangles, orange, yellow, green and cyan, represent $Q = +3, +1, -1$, and $-3$, respectively. (c)-(e) The ground state configurations for (c) $J < -1$, (d) $-1 < J < 0$, (e) $J > \frac{1}{3}$, whose magnetic unit cell is represented by a black rhombus. (f) One of the ground state configurations for $0 < J < \frac{1}{3}$. The dashed circles denote hexamers. The region enclosed by the bold line is one of the same-charge clusters, and pink (light blue) circles denote the internal (boundary) sites of the cluster.</figcaption></figure>

law;

$$ \sum_{p \in D} Q_p = \sum_{i \in \partial D} \sigma_i^z \eta_{p_D(i)}, \quad (3) $$

which follows the definition of $Q_p$. Here, $D$ is an arbitrary set of triangles and $\partial D$ means the sites on its boundary. $p_D(i)$ is a triangle in $D$ that includes the site $i$. This constraint strictly limits possible charge configurations.

Firstly, if we take $D$ as all triangles on a lattice, it directly leads to a global conservation law: $\sum_p Q_p = 0$, prohibiting, e.g., uniform charge ordering. Secondly, for arbitrary region $D$, the total amount of contained charges is bounded by at most the number of boundary sites, implying that a large-scale charge segregation, especially a phase separation, is suppressed.

Phase diagram at $T = 0$. At $J = 0$, the model is reduced to the Ising model with AF n.n. interaction, for which an exact solution is available [18]. At this limit, the ground state is disordered with a residual entropy (per site), $S_{\text{n.n.}} = 0.502$. In terms of the charge representation, the ground state manifold consists of the charge configurations with $Q = \pm 1$ on all the triangles.

For finite $J$, this degeneracy is lifted, and various ordered and disordered phases are stabilized. To see this, it is instructive to rewrite the Hamiltonian (2), as

$$ \mathcal{H} = \frac{1}{2} (1 + J) \sum_p Q_p^2 - \frac{J}{2} \sum_{\langle p, q \rangle} (Q_p + Q_q)^2. \quad (4) $$

For $J < 0$, the first and second terms can be minimized simultaneously, by setting (i) $Q_p + Q_q = 0$ for all the nearest-neighbor triangle pairs, and (ii) $|Q_p| = 1$ ($= 3$) on all the triangles for $-1 < J < 0$ ($J < -1$).

The condition (i) shows there is a staggered charge ordering. For $J < -1$, the ground state is a triple charge crystal, or a fully polarized ferromagnet (FM) [Fig. 1(c)].

Similarly, for $-1 < J < 0$, the system shows a single charge ordering. In this phase, the charge configuration is uniquely determined [Fig. 1(d)]. However, the underlying spin configurations are left disordered, even showing macroscopic degeneracy, which means a CSL ground state is realized. This CSL phase follows a simple constraint: all the triangles are occupied by, 2-up 1-down (1-up 2-down) configurations, corresponding to upward (downward) spin polarization. This CSL can be mapped to the kagome ice state observed in the [111] magnetization plateau for $Dy_2Ti_2O_7$ [30–33]. Equivalent states have been observed in artificial spin ice [34–37] and itinerant systems [38].

On the other hand, the situation for $J > 0$ is quite unusual, in that like charges attract each other, in contrast to the electrodynamics we are used to. At first sight, the representation (4) leads to a clustering of macroscopic number of like charges to lower the energy. However, it is forbidden by the Gauss’ constraint (3), which suppresses large-scale charge segregation. Under the Gauss’ constraint, how many like charges can be paired together?
To answer this question, firstly, let us limit ourselves to small $J$ ($0 < J \ll 1$), where triple charges are excluded. In this case, to maximize the contacts between same-sign charges, one may find ring-type clusters around hexagons of single charges with branches attached [Fig. 1(f)]. While rings made of more than six charges are energetically possible,

they are entropically strongly suppressed. If the system is covered with this type of clusters, each charge neighbors two charges of the same sign on average, with total energy $E_{tot} = \frac{1}{2}(1 - 3J)N_{\Delta}$. $N_{\Delta}$ is the total number of triangles of both orientations. Such configuration is macroscopically degenerate and composes the ground state for small $J$, as shown below. Accordingly, we call this novel state the hexamer CSL phase.

To see how this structure is stabilized, let us focus on the energy of a like-charge cluster, i.e., a maximal set of connected like-charge triangles. The energy of the $\alpha$-th cluster is

$$E_c^{(\alpha)} = \left(\frac{1}{2} - J\right) N_{\Delta}^{(\alpha)} - J n_i^{(\alpha)} + \frac{J}{2} n_b^{(\alpha)}. \quad (5)$$

Here, $N_{\Delta}^{(\alpha)}$ is the number of triangles in the $\alpha$-th cluster. $n_i^{(\alpha)}$ and $n_b^{(\alpha)}$ are the number of internal and boundary sites of the cluster, respectively [Fig. 1(f)], which satisfy $3N_{\Delta}^{(\alpha)} = 2n_i^{(\alpha)} + n_b^{(\alpha)}$. The term $\frac{J}{2}n_b^{(\alpha)}$ expresses the interaction between clusters, and the factor $1/2$ accounts for the double counting. This definition allows us to express the total energy of the system: $E_{tot} = \sum_{\alpha} E_c^{(\alpha)}$.

Now, we apply the Gauss' law (3) assuming a cluster of positive charges ($Q_p = +1$), without loss of generality. One obtains the inequality, $N_{\Delta}^{(\alpha)} \le n_b^{(\alpha)}$, which provides a lower bound of the system energy, $E_{tot} \ge \frac{1}{2}(1 - 3J)N_{\Delta}$. As calculated previously, this lower bound is reached for the type of configurations of Fig.1(f).

When increasing $J$, the hexamer CSL shows an instability to triple charge creation. By extending the discussion above, one can rigorously show that the phase boundary is at $J = 1/3$, and a 12-site magnetic order is stabilized for $J > 1/3$ [Fig. 1(e)]. For details, see Supplemental Material [39].

**Thermodynamic quantities.**- The hexamer CSL phase for $0 < J < 1/3$ gives characteristic features in thermodynamic quantities, especially when compared to the other disordered regime, kagome ice, for $-1 < J < 0$. To see this, we choose $J = -0.1$ and $0.1$, and perform a Monte Carlo simulation and Husimi-cactus calculations to obtain $T$ dependence of specific heat ($C$) and entropy ($S$) [Figs. 2(b)-2(e)]. The two methods match quantitatively well.

$C$ has two peaks both for $J = -0.1$ and $J = 0.1$. The high-temperature peak, around $T \sim 1.5$, corresponds to the entropy release due to the vanishing triple charges. Below this peak, $S$ takes a value close to the residual entropy of n.n. Ising model, $S \sim S_{n.n.} = 0.502$.

At low temperatures, however, thermodynamic behaviors are quite different between these two regions. For $J = -0.1$, $C$ diverges at $T_c \sim 0.163$ due to the long-range order of charges with ferrimagnetic ordering. While a previous study interpreted this transition to be of the three-state Potts universality class [20], recent state-of-the-art analyses in presence of long-range dipolar interactions have shown it to be of the Ising universality class [34, 35], as naively expected from the spontaneous $Z_2$ symmetry breaking of the staggered charge order.

Meanwhile, for $J = 0.1$, the specific heat never diverges.

<figure><img src="image_4.png" /><figcaption>FIG. 2. (Color online) (a) Finite-temperature phase diagram of the present model. Crosses denote the crossover between the paramagnet and the hexamer CSL. Black stars denote the exact phase boundaries at $T = 0$. (b)-(e) Temperature dependence of (b) $C$ and (c) $S$ at $J = -0.1$, and that of (d) $C$ and (e) $S$ at $J = 0.1$. Red lines and green dots are results of Husimi-cactus calculations and Monte Carlo simulations with $L = 64$, respectively. (f) Temperature dependence of $\rho_{hex}$ (the number of hexagons surrounded by same-sign charges divided by the total number of hexagons) for $J = 0.1$. (g) Autocorrelation, $R_{spin}(t) \equiv \frac{\sum_i \sigma_i^2(0) \sigma_i^2(t)}{N_{site}}$, for $J = 0.1$, $T = 0.03$. Red, green, and blue dots are for the single spin flip, the loop-update algorithm, and the new worm algorithm, respectively.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 3. (Color online) Magnetic structure factors for (J, T) = (a) (-0.1, 0.5), (b) (-0.1, 0.1), (c) (0.1, 0.5), and (d) (0.1, 0.1). Black lines denote Brillouin zones.</figcaption></figure>

The second, low-temperature, peak corresponds to a crossover where the hexamer clusters develop their structure [Fig. 2(f)]. Population of hexamer clusters introduces a problem in equilib- ration, since the cluster is locally quite stable, and a standard single-spin-flip dynamics freezes. Even a loop-update algo- rithm [34, 35, 40] effective for the kagome ice region is pow- erless here, since there is no underlying divergence-free field to support it. To overcome this difficulty, we developed a new type of worm algorithm [39], which enables nonlocal charge transportation between different clusters, and efficiently equi- librates the system in this region [Fig. 2(g)].

This non-local dynamics enables us to explore the ther- modynamic behavior below the second peak. In particular, we succeed in evaluating the *T* dependence of entropy pre- cisely, and find that it converges to an unfamiliar value of *S*<sub>0</sub><sup>HCSL</sup> ≈ 0.32, near zero temperature. This residual entropy is well reproduced by applying a variant of Pauling's argu- ment [41] to a nine-spin cluster made of four connected tri- angles. As a Pauling constraint, we impose that each triangle neighbors with two, and only two, same-sign charges, giving the entropy: *S*<sub>*p*</sub><sup>HCSL</sup> = <sup>1</sup>&frasl;<sub>6</sub> log (<sup>27</sup>&frasl;<sub>4</sub>) ≈ 0.32 [39].

*Magnetic structure factor.*- The unusual spatial structure of the CSLs can be detected by analyzing their magnetic struc- ture factors, $S(q) = \sum_{i,j} (\sigma_i^z \sigma_j^z) e^{-iq \cdot (r_i - r_j)}$.

The spatial structure of kagome ice region is characterized by the singularity of $S(q)$, called pinch point [42, 43]. In fact, upon cooling, $S(q)$ develops bow-tie structures in the centers of the second Brillouin zone, even above $T_c$ [Fig. 3(a)]. Be- low $T_c$, where the kagome ice correlation is well developed, clear kinks are observed [Fig. 3(b)], coexisting with the Bragg peaks of ferrimagnetic order [36, 37, 44–46].

In contrast, at $J = 0.1$, $S(q)$ develops a quite different pat- tern. Pinch points are missing, clearly indicating the hexamer CSL phase is not a Coulomb phase. Instead, one can see the “half-moon” patterns surrounding the zone centers [Fig. 3(d)]. These patterns give the evidence for the hexamer rings, as re- cently observed on a fine-tuned point of the $J_1-J_2-J_3$ model on pyrochlores [28, 29].

*Discussion and Summary.*- We have studied the collec-

tive behaviors of topological charges and their instabilities, through the $J_1-J_2-J_3$ Ising model on the kagome lattice. We found that the fractionalized topological charges are recom- bined into a novel CSL, as a result of a keen interplay between the topological constraint and interactions. This hexamer CSL is characterized by the hexamer clusters, whose configura- tional pattern leads to a unconventional value of residual en- tropy. This spatial pattern can be detected through the half- moon signal in the magnetic structure factor.

The hexamer CSL has a number of interesting features. Among them, the self-organization of multi-scale structure is remarkable. While the starting Hamiltonian is homoge- neous, each ground state configuration composed of the hex- amer cluster coverings, shows quite inhomogeneous spatial structures. In the presence of quantum or stochastic dynamical processes, it is anticipated that this multi-scale structure leads to rich dynamical spectra, possibly reminiscent of glassy sys- tems. The hexamer CSL also offers a promising extended re- gion ($0 < J < 1/3$) for the search of unconventional quantum spin liquids [22–26, 47–49] and quantum order-by-disorder mechanisms [50–53] beyond the traditional Heisenberg anti- ferromagnet.

It is tempting to speculate if the interaction between topo- logical charges underlies hexagonal or small-cluster forma- tion found in a broad range of frustrated magnets. For instance, cluster excitations of antiferromagnetic hexagonal loops were reported in the cubic spinel $ZnCr_2O_4$ by inelastic neutron scattering [54].

Indeed, the charge representation is useful for a broad range of spin models on frustrated lattices, and it has been used to resolve the spin liquid properties [49]. We hope the viewpoint of topological charges and their recombination will be useful for further explorations of novel states of matter in frustrated magnets, and especially their topological properties.

This work is supported by JSPS KAKENHI (Nos. 26400339, 24340076, 15H05852 and 15K13533) and the Ok- inawa Institute of Science and Technology Graduate Uni- versity. T. M. is supported by Advanced Leading Graduate Course for Photon Science (ALPS).
* mizoguchi@hosi.phys.s.u-tokyo.ac.jp

[1] A. Y. Kitaev, Phys.-Usp. **44**, 131 (2001).

[2] A. J. Heeger, S. Kivelson, J. R. Schrieffer, and W.-P. Su, Rev. Mod. Phys. **60**, 781 (1988).

[3] R. B. Laughlin, Phys. Rev. Lett. **50**, 1395 (1983).

[4] V. J. Goldman and B. Su, Science **267**, 1010 (1995).

[5] L. Balents, Nature (London) **464**, 199 (2010).

[6] Z. Hiroi, M. Hanawa, N. Kobayashi, M. Nohara, H. Takagi, Y. Kato, and M. Takigawa, J. Phys. Soc. Jpn. **70**, 3377 (2001).

[7] Y. Shimizu, K. Miyagawa, K. Kanoda, M. Maesato, and G. Saito, Phys. Rev. Lett. **91**, 107001 (2003).

[8] M. P. Shores, E. A. Nytko, B. M. Bartlett, and D. G. Nocera, J. Am. Chem. Soc. **127**, 13462 (2005).

[9] Y. Okamoto, M. Nohara, H. Aruga-Katori, and H. Takagi, Phys. Rev. Lett. **99**, 137207 (2007).

[10] Y. Okamoto, H. Yoshida, and Z. Hiroi, J. Phys. Soc. Jpn. **78**, 033701 (2009).

[11] Y. Singh and P. Gegenwart, Phys. Rev. B **82**, 064412 (2010).

[12] B. Fåk, E. Kermarrec, L. Messio, B. Bernu, C. Lhuillier, F. Bert, P. Mendels, B. Koteswararao, F. Bouquet, J. Ollivier, A. D. Hillier, A. Amato, R. H. Colman, and A. S. Wills, Phys. Rev. Lett. **109**, 037208 (2012).

[13] T. Takayama, A. Kato, R. Dinnebier, J. Nuss, H. Kono, L. S. I. Veiga, G. Fabbris, D. Haskel, and H. Takagi, Phys. Rev. Lett. **114**, 077202 (2015).

[14] M. J. P. Gingras and P. A. McClarty, Rep. Prog. Phys. **77**, 056501 (2014).

[15] L. Savary and L. Balents, Rep. Prog. Phys. **80**, 016502 (2017).

[16] C. Castelnovo, R. Moessner, and S. L. Sondhi, Nature (London) **451**, 42 (2008).

[17] L. D. C. Jaubert and P. C. W. Holdsworth, Nat. Phys. **5**, 258 (2009).

[18] K. Kano and S. Naya, Prog. Theor. Phys. **10**, 158 (1953).

[19] P. Azaria, H. T. Diep, and H. Giacomini, Phys. Rev. Lett. **59**, 1629 (1987).

[20] M. Wolf and K. D. Schotte, J. Phys. A: Math. Gen. **21**, 2195 (1988).

[21] T. Takagi and M. Mekata, J. Phys. Soc. Jpn. **62**, 3943 (1993).

[22] L. Balents, M. P. A. Fisher, and S. M. Girvin, Phys. Rev. B **65**, 224412 (2002).

[23] S. V. Isakov, Y. B. Kim, and A. Paramekanti, Phys. Rev. Lett. **97**, 207204 (2006).

[24] Y.-C. He, D. N. Sheng, and Y. Chen, Phys. Rev. Lett. **112**, 137202 (2014).

[25] S. Bieri, L. Messio, B. Bernu, and C. Lhuillier, Phys. Rev. B **92**, 060407(R) (2015).

[26] Y. Iqbal, H. O. Jeschke, J. Reuther, R. Valenti, I. I. Mazin, M. Greiter, and R. Thomale, Phys. Rev. B **92**, 220404(R) (2015).

[27] H. Ishizuka and Y. Motome, Phys. Rev. B **88**, 100402(R) (2013).

[28] M. Udagawa, L. D. C. Jaubert, C. Castelnovo, and R. Moessner, Phys. Rev. B **94**, 104416 (2016).

[29] J. G. Rau and M. J. P. Gingras, Nat. Communs. **7**, 12234 (2016).

[30] K. Matsuhira, Z. Hiroi, T. Tayama, S. Takagi, and T. Sakakibara, J. Phys.: Condens. Matter **14**, L559 (2002).

[31] M. Udagawa, M. Ogata, and Z. Hiroi, J. Phys. Soc. Jpn. **71**, 2365 (2002).

[32] S. V. Isakov, K. S. Raman, R. Moessner and S. L. Sondhi, Phys. Rev. B **70**, 104418 (2004).

[33] H. Aoki, T. Sakakibara, K. Matsuhira, and Z. Hiroi, J. Phys. Soc. Jpn. **73**, 2851 (2004).

[34] G. Möller and R. Moessner, Phys. Rev. B **80**, 140409(R) (2009).

[35] G.-W. Chern, P. Mellado, and O. Tchernyshyov, Phys. Rev. Lett. **106**, 207202 (2011); G.-W. Chern and O. Tchernyshyov, Phil. Trans. R. Soc. A **370**, 5718 (2012).

[36] B. Canals, I.-A. Chioar, V.-D. Nguyen, M. Hehn, D. Lacour, F. Montaigne, A. Locatelli, T. O. Menteş, B. S. Burgos, and N. Rougemaille, Nat. Commun. **7**, 11446 (2016).

[37] O. Sendetskyi, L. Anghinolfi, V. Scagnoli, G. Möller, N. Leo, A. Alberca, J. Kohlbrecher, J. Lüning, U. Staub, and L. J. Heyderman, Phys. Rev. B **93** 224413 (2016).

[38] H. Ishizuka and Y. Motome, Phys. Rev. B **88**, 081105(R) (2013).

[39] See Supplemental Material for detailed discussions on a phase boundary in $J > 0$, new algorithms for global updates, and Pauling entropy in hexamer CSL phase.

[40] G. T. Barkema and M. E. Newman, Phys. Rev. E **57**, 1155 (1998).

[41] L. Pauling, J. Am. Chem. Soc. **57**, 2680-2684 (1935).

[42] R. W. Youngblood and J. D. Axe, Phys. Rev. B **23**, 232 (1981).

[43] S. V. Isakov, K. Gregor, R. Moessner, and S. L. Sondhi, Phys. Rev. Lett. **93**, 167204 (2004).

[44] M. E. Brooks-Bartlett, S. T. Banks, L. D. C. Jaubert, A. Harman-Clarke, and P. C. W. Holdsworth, Phys. Rev. X, **4**, 011007 (2014).

[45] S. Petit, E. Lhotel, B. Canals, M. C. Hatnean, J. Ollivier, H. Mutka, E. Ressouche, A. R. Wildes, M. R. Lees, and G. Balakrishnan, Nat. Phys. **12**, 746 (2016).

[46] J. A. M. Paddison, H. S. Ong, J. O. Hamp, P. Mukherjee, X. Bai, M. G. Tucker, N. P. Butch, C. Castelnovo, M. Mourigal, and S. E. Dutton, arXiv: 1605. 01423 (2016).

[47] B. Bauer, L. Cincio, B.P. Keller, M. Dolfi, G. Vidal, S. Trebst, and A.W.W. Ludwig, Nature Commun. **5**, 5137 (2014).

[48] K. Essafi, O. Benton, and L. D. C. Jaubert, Nat. Communs. **7**, 10297 (2016).

[49] F. L. Buessen and S. Trebst, Phys. Rev. B **94**, 235138 (2016).

[50] A. L. Chernyshev and M. E. Zhitomirsky, Phys. Rev. Lett. **113**, 237202 (2014).

[51] O. Götze and J. Richter, Phys. Rev. B **91**, 104402 (2015).

[52] A. L. Chernyshev and M. E. Zhitomirsky, Phys. Rev. B **92**, 144415 (2015).

[53] J. Oitmaa and R. R. P. Singh, Phys. Rev. B **93**, 014424 (2016).

[54] S.-H. Lee, C. Broholm, W. Ratcliff, G. Gasparovic, Q. Huang, T. H. Kim, and S.-W. Cheong, Nature (London) **418**, 856 (2002).
Supplemental Material for “Clustering of Topological Charges in Kagome Classical Spin Liquid”

Tomonari Mizoguchi,¹ L. D. C. Jaubert,²,³ and Masafumi Udagawa⁴

¹Department of Physics, University of Tokyo, Hongo, Bunkyo-ku, Tokyo 113-0033, Japan*

²Okinawa Institute of Science and Technology Graduate University, Onna-son, Okinawa 904-0495, Japan

³CNRS, Univ. Bordeaux, LOMA, UMR 5798, F-33400 Talence, France

⁴Department of Physics, Gakushuin University, Mejiro, Toshima-ku, Tokyo 171-8588, Japan

(Dated: February 14, 2017)

# I. DETERMINATION OF THE PHASE BOUNDARY IN $J > 0$

In this section, we present detailed discussions on the ground-state spin configurations and phase boundaries for $J > 0$. In the main text, we briefly discuss the ground-state structure for $0 < J \ll 1$, by estimating the energy of same-charge clusters [Eq. (5) in the main text]. Here, we give the details of the derivation, and extend the argument to larger $J$, where triple charges are also allowed to exist. In Sec. IA and Sec. IB, we summarize the basic equations and the analysis of ground-state energy for single-charge cluster, in a slightly different notation from the main text, for later convenience. In Sec. IC, we describe the derivation of ground-state spin configuration and phase boundary for general $J(> 0)$.

## A. Basic equations

We base our arguments on the Gauss’ law as introduced in Eq. (3) in the main text. This equation results from the definition of charge defined on a triangle, $p$:

$$Q_p = \eta_p \sum_{i \in p} \sigma_i^z, \quad (1)$$

where $\eta_p = +1(-1)$ for $p \in \Delta(\nabla)$. If we consider a set of triangles, $D$, and sum Eq. (1) over $D$, we obtain the Gauss’ law:

$$\sum_{p \in D} Q_p = \sum_{i \in \partial D} \eta_{p_D(i)} \sigma_i^z. \quad (2)$$

As noted in the main text, $\partial D$ means the sites on its boundary, and $p_{D(i)}$ is a triangle in $D$ that includes the site $i$. Here, we have used the fact that each site, $i$, belongs to one upward and one downward triangles. Accordingly, the summation of $\eta_{p_D(i)} \sigma_i^z$ vanishes for internal sites of $D$, leaving only the boundary contribution.

If $D$ is composed of $N_{\Delta}^D$ triangles (of both orientations), the number of internal sites, $n_i^D$, and the boundary sites, $n_b^D$ satisfy the relation:

$$3N_{\Delta}^{D} = 2n_{i}^{D} + n_{b}^{D}. \quad (3)$$

The number of boundary sites, $n_b$, limits the amount of charge in $D$. By taking the absolute value of Eq. (2), we have

$$|\sum_{p \in D} Q_p| \le \sum_{i \in \partial D} |\eta_{p_D(i)}| |\sigma_i^z| = \sum_{i \in \partial D} 1 = n_b^D. \quad (4)$$

On the basis of Eqs. (2), (3) and (4), we can obtain various information regarding the ground states.

## B. Single-charge cluster

In this subsection, we summarize the analysis of region: $0 < J \ll 1$, where only the single charges are allowed to exist. Let us define a same-charge cluster as a maximal set of connected same-charge triangles. We take the $\alpha$-th cluster composed of $N_{\Delta}^{(\alpha)}$ triangles, with $n_i^{(\alpha)}$ internal and $n_b^{(\alpha)}$ boundary sites. If we apply Eq. (4) to this cluster, we have

$$N_{\Delta}^{(\alpha)} \le n_b^{(\alpha)}, \quad (5)$$

since we only consider single charges.

Next, we compute the energy of the $\alpha$-th cluster,

$$E_c^{(\alpha)} = \left(\frac{1}{2} - J\right) N_{\Delta}^{(\alpha)} - J n_i^{(\alpha)} + \frac{J}{2} n_b^{(\alpha)}, \quad (6)$$

as in Eq. (5) in the main text. This definition allows us to express the total energy of the system, $E_{\text{tot}}$, as the summation of $E_c^{(\alpha)}$ over the clusters: $E_{\text{tot}} = \sum_{\alpha} E_c^{(\alpha)}$.

By combining Eqs. (3), (5) and (6), we obtain

$$E_c^{(\alpha)} \ge \frac{1}{2}(1 - 3J)N_{\Delta}^{(\alpha)}. \quad (7)$$

By summing up all the cluster contributions, we obtain the lower bound of the total energy:

$$E_{\text{tot}} \ge \frac{1}{2}(1 - 3J)N_{\Delta} \equiv E_{\text{tot}}^{s}, \quad (8)$$

with the total number of triangles, $N_{\Delta}$. $E_{\text{tot}}^s$ gives the lowest energy of the system, given that only single charges exist. $E_{\text{tot}}^s$, in fact, gives the ground state energy for small $J$.

<footnote>* mizoguchi@hosi.phys.s.u-tokyo.ac.jp</footnote>
### C. Triple-charge cluster

#### 1. Structure of triple-charge cluster

Next, we consider a general value of $J(>0)$, where triple charges are also allowed to exist. Firstly, we prove that there exists at most one triple charge in a same-charge cluster. To show this, let us assume there are more than two triple charges in a given cluster. Then, as a subset of cluster, there must be a string of same-charge triangles with two triple charges on both ends, which we call $\mathcal{D}$. This set of triangles satisfies $n_b^D = N_\Delta^D + 2$. Accordingly, we obtain

$$ \left| \sum_{p \in \mathcal{D}} Q_p \right| = N_{\Delta}^{\mathcal{D}} + 4 > n_b^{\mathcal{D}}. \quad (9) $$

This contradicts Eq. (4), and consequently, the hypothesis of more than two triple charges in a cluster is proven to be wrong. Accordingly, each cluster contains at most one triple charge. Hereafter, we call a cluster without triple charge, a s-cluster (sc), and a cluster with one triple charge, a t-cluster (tc). The α-th s-cluster composed of $N_\Delta^{(\alpha)}$ triangles, with $n_i^{(\alpha)}$ internal and $n_b^{(\alpha)}$ boundary sites, satisfies the relation, Eq. (5). Meanwhile, the β-th t-cluster composed of $N_\Delta^{(\beta)}$ triangles, with $n_i^{(\beta)}$ internal and $n_b^{(\beta)}$ boundary sites, satisfies

$$ N_{\Delta}^{(\beta)} + 2 \le n_b^{(\beta)}, \quad (10) $$

as a result of applying Eq. (2).

#### 2. Cluster energy

Next, we define a cluster energy in the presence of triple charge, so that the total energy is expressed as the sum of cluster energies:

$$ E_{\text{tot}} = \sum_{\alpha \in \text{sc}} E_c^{s,(\alpha)} + \sum_{\beta \in \text{tc}} E_c^{t,(\beta)}. \quad (11) $$

Here, $E_c^{s,(\alpha)}$ ($E_c^{t,(\beta)}$) is the contribution from a cluster without (with) triple charge. In the presence of surrounding triple charges, $E_c^{s,(\alpha)}$ is modified from the formula (6). Even though the cluster does not contain any triple charge itself, it can still interact with neighboring triple charges.

$$ E_c^{s,(\alpha)} = \left(\frac{1}{2} - J\right) N_{\Delta}^{(\alpha)} - J n_i^{(\alpha)} + \frac{J}{2} n_b^{ss(\alpha)} + \frac{3J}{2} n_b^{st(\alpha)}, \quad (12) $$

where $n_b^{ss(\alpha)}$ ($n_b^{st(\alpha)}$) means the number of boundary sites between single charges (single and triple charges). They satisfy

$$ n_b^{ss(\alpha)} + n_b^{st(\alpha)} = n_b^{(\alpha)}. \quad (13) $$

Accordingly, Eq. (12) can be rewritten as

$$ E_c^{s,(\alpha)} = \left(\frac{1}{2} - J\right) N_{\Delta}^{(\alpha)} - J n_i^{(\alpha)} + \frac{J}{2} n_b^{(\alpha)} + J n_b^{st(\alpha)}. \quad (14) $$

A cluster with a triple charge composed of $N_\Delta^{(\beta)}$ triangles have $N_\Delta^{(\beta)} - 1$ single, and one triple charges. We assume the triple charge neighbors with $n^{(\beta)}(=0, 1, 2, 3)$ single charges within the cluster, and $3 - n^{(\beta)}$ charges of other clusters. Accordingly, we write the cluster energy as

$$ E_c^{t,(\beta)} = \left(\frac{1}{2} - J\right)(3^2 + (N_\Delta^{(\beta)} - 1)) - J[3n^{(\beta)} + (n_i^{(\beta)} - n^{(\beta)})] \\ + \frac{J}{2}n_b^{ss(\beta)} + \frac{3J}{2}n_b^{st(\beta)} + \frac{3J}{2}n_b^{ts(\beta)} + \frac{9J}{2}n_b^{tt(\beta)}. \quad (15) $$

Here, the first and the second terms come from the self-energy of charges and the interaction within the cluster, respectively. The other terms denote four possible types of inter-cluster interactions. $n_b^{ss(\beta)}$ and $n_b^{tt(\beta)}$ denote the number of sites between single charges and triple charges, respectively; $n_b^{st(\beta)}$ ($n_b^{ts(\beta)}$) denotes the number of sites between single (triple) charges in the β-th cluster and triple (single) charges in other clusters. These numbers satisfy

$$ n_b^{ss(\beta)} + n_b^{st(\beta)} = n_b^{(\beta)} - (3 - n^{(\beta)}), \quad (16) $$

and

$$ n_b^{ts(\beta)} + n_b^{tt(\beta)} = (3 - n^{(\beta)}). \quad (17) $$

Accordingly, Eq. (15) can be rewritten as

$$ E_c^{t,(\beta)} = \left(\frac{1}{2} - J\right)(3^2 + (N_\Delta^{(\beta)} - 1)) - J[3n^{(\beta)} + (n_i^{(\beta)} - n^{(\beta)})] \\ + \frac{J}{2}[n_b^{(\beta)} - (3 - n^{(\beta)})] + \frac{3J}{2}(3 - n^{(\beta)}) + Jn_b^{st(\beta)} + 3Jn_b^{tt(\beta)}. \quad (18) $$

With these definitions of $E_c^{s,(\alpha)}$ and $E_c^{t,(\beta)}$, the total energy can be written as in Eq. (11).

Let us evaluate the total energy. By combining eqs. (11), (14) and (18), we have

$$ E_{\text{tot}} = \sum_{\alpha \in \text{sc}} \left( \frac{1}{2} - J \right) N_{\Delta}^{(\alpha)} - J n_{i}^{(\alpha)} + \frac{J}{2} n_{b}^{(\alpha)} \\ + \sum_{\beta \in \text{tc}} \left( \frac{1}{2} - J \right) (3^2 + (N_{\Delta}^{(\beta)} - 1)) - J[3n^{(\beta)} + (n_{i}^{(\beta)} - n^{(\beta)})] \\ + \frac{J}{2}[n_{b}^{(\beta)} - (3 - n^{(\beta)})] + \frac{5J}{2}(3 - n^{(\beta)}) + 2Jn_{b}^{tt(\beta)}. \quad (19) $$

Here, we have used

$$ \sum_{\alpha \in \text{sc}} n_b^{st(\alpha)} + \sum_{\beta \in \text{tc}} (n_b^{st(\beta)} + n_b^{tt(\beta)}) = \sum_{\beta \in \text{tc}} (3 - n^{(\beta)}). \quad (20) $$

Now, we apply the relation (3), and the inequalities due to Gauss' constraint, (5) and (10) to Eq. (19), and obtain the inequality:

$$ E_{\text{tot}} \ge \frac{1}{2}(1 - 3J)N_{\Delta} + 4 \sum_{\beta \in \text{tc}} (1 - Jn^{(\beta)}). \quad (21) $$
The equality holds if (i) $N_{\Delta}^{(\alpha)} = n_{b}^{(\alpha)}$, (ii) $N_{\Delta}^{(\beta)} + 2 = n_{b}^{(\beta)}$ for $\alpha \in \text{sc}$ and $\beta \in \text{tc}$, and (iii) there are no triple charges next to each other, $n_{b}^{\text{tt}(\beta)} = 0$. Since $E_{\text{tot}}^{s} = \frac{1}{2}(1-3J)N_{\Delta}$ gives the lowest energy given only s-clusters, as defined in Eq. (8), the correction of the total energy due to t-clusters can be evaluated as

$$E_{\text{tot}} - E_{\text{tot}}^{s} \geq 4 \sum_{\beta \in \text{tc}} (1 - Jn^{(\beta)}). \quad (22)$$

### 3. Ground-state energy

Finally, we discuss the ground-state energy of the system, on the basis of Eq. (22). Since $n^{(\beta)} \le 3$, Eq. (22) implies that triple charges appear for $J > 1/3$. This is indeed the case; the hexamer classical spin liquid (CSL) is destabilized at $J = 1/3$, and 12-site cluster ordering with triple charges is stabilized for $J > 1/3$, as we show below.

To prove this, we firstly estimate the lower bound of the righthand side of Eq. (22). Then, we look for the spin configuration that gives this minimum energy. To this aim, we introduce $M_n$, the number of t-clusters with $n^{(\beta)} = n$. Since a t-cluster with $n^{(\beta)} = n$ is composed of at least $n+1$ triangles, $M_n$ satisfies

$$M_{\text{tot}} \equiv \sum_{n=0}^{3} (1+n)M_n \le N_{\Delta}. \quad (23)$$

Then, we can rewrite Eq. (22) as

$$\begin{aligned} E_{\text{tot}} - E_{\text{tot}}^{s} &\ge 4 \sum_{n=0}^{3} (1 - Jn)M_n \\ &= -4JM_{\text{tot}} + 4(1 + J) \sum_{n=0}^{3} M_n. \end{aligned} \quad (24)$$

For any given value of $M_{\text{tot}}$, the total number of clusters $\sum_{n=0}^{3} M_n$ is minimized, if $M_0 = M_1 = M_2 = 0$ and $M_3 = M_{\text{tot}}/4$, leading to

$$\begin{aligned} E_{\text{tot}} - E_{\text{tot}}^{s} &\ge (1 - 3J)M_{\text{tot}} \\ &\ge \begin{cases} 0 & (0 < J \le 1/3), \\ (1 - 3J)N_{\Delta} & (J > 1/3). \end{cases} \end{aligned} \quad (25)$$

This inequality shows that the hexamer CSL is stabilized for $0 < J \le 1/3$. Meanwhile, it implies a more stable state for $J > 1/3$. Indeed, if one can realize a spin configuration that satisfies above conditions (i), (ii) and (iii), and takes $(M_0, M_1, M_2, M_3) = (0, 0, 0, N_{\Delta}/4)$, it will be the ground state in this region, with the ground-state energy, $E_{\text{tot}} - E_{\text{tot}}^{s} = (1 - 3J)N_{\Delta}$.

One can actually realize such spin configuration by packing four-triangle clusters made of a central triple charge, surrounded by three single charges. [see Fig. 1(a)]. To pack these clusters under the constraint of total charge conservation, the magnetic unit cell has $2 \times 2$

<figure><img src="image_6.png" /><figcaption>FIG. 1. (a) 4-triangle cluster and (b) a magnetic unit cell of the 12-site ordering.</figcaption></figure>

structure with 12 distinct sites, thus we call this state “12-site ordering” [see Fig. 1(b)]. The energy of this state is given by

$$E_{\text{tot}} - E_{\text{tot}}^{s} = (1 - 3J)N_{\Delta}, \quad (26)$$

which is the minimum of energy determined in Eq. (25). This completes the proof that no other phases exist for $J > 1/3$.

## II. DETAILS OF GLOBAL UPDATES

In this section, we will present the details of three types of global update algorithms. We first explain the loop update [1–3], which accelerates the relaxation of the system in low-temperature region for $-1 < J < 0$. This algorithm is useful when the magnetic charges have a long-range order, while the spins are left disordered. In this case, local spin flip creates a high-energy state by changing the charge state, and is hardly accepted, accordingly. Meanwhile, the loop update relaxes the spin configuration effectively, without changing charge configuration. Secondly, we explain how to get rid of metastable defects by using the string update which causes the pair annihilation of the defects. Finally, we introduce a new worm algorithm, which is implemented in low-temperature region for $0 < J < \frac{1}{3}$. This algorithm transports a charge between different clusters, and relaxes both spin and charge configurations without changing the energy.

### A. Loop update

In the CSL state with staggered charge order, which is stabilized for $-1 < J < 0$, local spin updates that change the charge configuration are hardly accepted in the low-temperature regions, since such updates cost a lot of energy. Hence, we need an algorithm that relaxes spin configuration without changing the charge configuration. The loop update serves this purpose [1–3]. This algorithm is illustrated in Fig. 2 and consists of finding
a closed loop of spins aligned alternately, like ↑↓↑↓... During the construction of the loop, because of the underlying staggered charge order, the spins in the loop are successively chosen with probability 1 (no choice) and 0.5 (a choice between two spins). Once a given spin is visited twice, the process stops and the global path of spins can be divided into a closed loop and a dangling tail. The spins in the closed loop are flipped, while the ones in the dangling tail are left unchanged. Since the charge configuration is left unchanged, there is no energy cost and this loop update is always accepted in Monte Carlo simulations.

<figure><img src="image_7.png" /><figcaption>FIG. 2. Schematic picture of the loop update, formed by a closed loop and a dangling tail.</figcaption></figure>

## B. String update

In the low temperature region of $-1 < J < 0$, where the staggered charge order is stabilized, same-charge clusters are sometimes left unrelaxed, and remain as metastable defects [see the black circle in Fig. 3]. These defects are difficult to remove by single spin flip, since local spin updates conserve charge locally. This would result in triple charges which are unfavored for $-1 < J < 0$. To eliminate these defects, we need to pair-annihilate two defects of opposite sign, which exist in spatially separate places.

To this aim, we implement a string update algorithm. In this update, we first choose a triangle randomly and select one of the majority spins as a starting site. Then, we make a string on which spins are aligned alternately: ↑↓↑↓..., in the same way as in the loop update algorithm introduced in Sec. II A. However, in the string update, we do not make a closed loop, but terminate a string with open ends, at randomly-chosen length. If the string update is accepted by a Metropolis argument, then all the spins along that string are flipped. This attempt transports charge from one end to the other. Accordingly, if the chosen string has defects of opposite signs on both ends, they can be pair-annihilated [see Fig. 3]. To eliminate defects, however, we need many attempts, since the starting site and the length of the string are chosen randomly.

<figure><img src="image_8.png" /><figcaption>FIG. 3. Schematic picture of the string update.</figcaption></figure>

## C. Worm algorithm

In the hexamer CSL realized for $0 < J < \frac{1}{3}$, low-temperature configurations are locally constrained by the attraction between same-sign charges. Single-spin-flip dynamics do not suffice to explore the extensively degenerate ensemble of ground states. This ergodicity loss is further complicated by the absence of any long-range charge order, which prevents the use of traditional loop algorithms [1] as described in section II A. To restore ergodicity, a new kind of worm algorithm is thus needed, able to update both spin and charge configurations.

To introduce the main idea of this algorithm, a careful understanding of the structure of same-charge clusters is required. Triangles within a cluster can be classified into three groups, according to how many like-charges they are neighboring:

* "end-point" triangles, sitting at the extremities of a cluster, with only one neighboring triangle carrying the same topological charge;

* "middle" triangles, forming the core of the branches of the clusters, with two neighboring same-sign triangles;

* "branching triangles", surrounded by three neighboring same-sign triangles.

The main idea of the algorithm developed here is to detach an end-point charge from a given cluster, and to stick it back either somewhere else onto the same cluster, or to a different cluster of the same sign. Along the path of this end-point charge, the intervening cluster of opposite sign is also updated. This is done by an extension of the string update algorithm introduced in section IIB. As illustrated in Fig. 4, we firstly choose an end-point triangle of positive (negative) sign randomly. This triangle is by definition surrounded by two triangles of opposite sign; the string starts from one of the two corresponding boundary spins with equal probability. However, in contrast to the string update introduced in Sec. IIB, the string propagates within the neighboring cluster of opposite sign (and not on a background of staggered charge order), and is not terminated at random. Alternatively, we monitor the type of triangles at each step, and stop the update as soon as we reach a middle triangle of negative (positive) sign. At this point, we flip the
constructed open string of spins. As a result, a positive (negative) end-point charge can be moved around in the system [Fig. 4].

At the end of the process, the total energy of the sys-
tem is preserved, such that the attempt is always ac-
cepted by the Metropolis rule. Furthermore, the key in
this algorithm lies in the necessary respect of detailed
balance. Indeed, while the initial triangle always change
from an “end-point” to a “middle” state, the final tri-
angle always follows the reverse change of status. And
since we stop the string update at the first middle triangle
encountered, it guarantees the necessary one-to-one cor-
respondence between the string update and its reversed
process. These conditions ensure the detailed balance of
stochastic processes.

<figure><img src="image_9.png" /><figcaption>FIG. 4. Schematic picture of the worm algorithm.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 5. Temperature dependence of the average length of worms.</figcaption></figure>

Below, we show the benchmark of the algorithm. The
average length of worms, ⟨ℓ<sub>worm</sub>⟩, is shown in Fig. 5. One
can see that the average length is almost temperature-
independent and its value is about 1.67, which indicates
that the worm algorithm changes the states locally rather
than globally.

Nevertheless, the worm algorithm improves the re-
laxation drastically; to show its quantitative evidence,
we calculate the autocorrelation. The autocorrela-

tions of spins, $R_{\text{spin}} \equiv \frac{\sum_i [\sigma_i^z(0)\sigma_i^z(t)]}{N_{\text{site}}}$, and charges,
$R_{\text{charge}} \equiv \frac{\sum_{p \in \Delta} [Q_p(0)Q_p(t)]}{N_{\Delta}}$, are shown in Fig. 6, re-
spectively. Here one Monte Carlo step consists of
$N_{\text{single}}$ single spin flip, $N_{\text{loop}}$ loop updates, and $N_{\text{worm}}$
worm updates. In Fig. 6, these numbers are cho-
sen as $(N_{\text{single}}, N_{\text{loop}}, N_{\text{worm}}) = (N_{\text{site}}, 0, 0)$ for red dots,
$(N_{\text{single}}, N_{\text{loop}}, N_{\text{worm}}) = (N_{\text{site}}, N_{\text{site}}/8, 0)$ for green dots,
and $(N_{\text{single}}, N_{\text{loop}}, N_{\text{worm}}) = (N_{\text{site}}, 0, N_{\text{site}})$ for blue
dots, respectively. One can clearly see that $R_{\text{spin}}$ and
$R_{\text{charge}}$ approach to zero very rapidly, if we introduce the
new worm algorithm. On the other hand, the loop up-
date algorithm is not efficient in hexamer CSL state. In
particular, $R_{\text{charge}}$ is unchanged by loop update, since it
does not change the charge configuration.

$N_{\text{worm}}$ dependence of the autocorrelations are shown in Fig. 7. Clearly, the larger $N_{\text{worm}}$ becomes, the faster the relaxation is achieved. Even in the smallest number of $N_{\text{worm}}(=\frac{N_{\text{site}}}{32})$, $R_{\text{spin}}$ and $R_{\text{charge}}$ go to zero within 150 Monte Carlo steps. This demonstrates that the worm algorithm is indeed a powerful technique to relax the system into equilibrium.

<figure><img src="image_11.png" /><figcaption>FIG. 6. The autocorrelation of the spin and charges. ($N_{\text{single}}$, $N_{\text{loop}}$, $N_{\text{worm}}$) for red, green, and blue dots are described in the text.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 7. The autocorrelation of the spin and charges for various values of $N_{\text{worm}}$.</figcaption></figure>

### III. PAULING'S ARGUMENT IN NEW CLASSICAL SPIN LIQUID

In this section, we discuss the residual entropy of hex-
amer CSL on the basis of Pauling’s argument. As we
<figure><img src="image_13.png" /><figcaption>FIG. 8. Charge configurations on a four-triangle cluster. Here we fix the spin configuration on the central triangle. Spin degeneracy of each configuration is shown below each snapshot.</figcaption></figure>

mentioned in the main text, the idea is that we consider a four-triangle cluster consisting of nine spins as a unit [see Fig. 8]. Suppose that the central triangle has a charge $Q = +1$, without any loss of generality. For a given fixed spin configuration on the central triangle, charge configurations of the other three triangles are classified in six categories, as shown in Figs. 8(a)-(f). They are classified as follows:

(a) Three same-sign-charge pairs.

(b) One same-sign-charge pair and two opposite-sign-charge pairs: The minority spin of the central triangle connects opposite-sign charges.

(c) One same-sign-charge pair and two opposite-sign-charge pairs: The minority spin of the central triangle connects same-sign charges.

(d) Two same-sign-charge pairs and one opposite-sign-charge pair: The minority spin of the central triangle connects opposite-sign charges.

(e) Two same-sign-charge pair and one opposite-sign-charge pair: The minority spin of the central triangle connects same-sign charges.

(f) Three opposite-sign-charge pairs.

Each charge configuration has a spin degeneracy (with respect to the six exterior spins), which is also presented in the figures.

Among these configurations, (d) and (e) have two like-charge pairs, which is equal to the average number of like-charge pairs in hexamer CSL (see the main text). Thus, the residual entropy can be estimated by considering the probability that each four-triangle cluster has a configuration of either (d) or (e). Note that we also have to consider the case that the central triangle has $Q = -1$. That probability, $P_{d+e}$, can be derived as follows:

1. A four-triangle cluster consist of nine spins, so there are $2^9$ states.

2. A central triangle has either $Q = +1$ and $Q = -1$, and the number of such spin configuration is six.

3. Once we fix the spin configuration on a central triangle, the number of states with patterns of (d) and (e) is $1 + 8 = 9$.

Then we obtain

$$P_{d+e} = \frac{6 \times 9}{2^9} = \frac{3^3}{2^8}. \qquad (27)$$

The corresponding residual entropy becomes

$$S_{P}^{\text{HCSL}} = \frac{1}{N_{\text{site}}} \log \left( 2^{N_{\text{site}}} P_{d+e}^{N_c} \right) = \frac{1}{6} \log \left( \frac{3^3}{2^2} \right) \sim 0.318, \qquad (28)$$

where $N_c$ denotes the total number of four-triangle clusters,

$$N_c = \frac{1}{4} N_{\Delta} = \frac{N_{\text{site}}}{6}, \qquad (29)$$

and $2^{N_{\text{site}}}$ is the total number of states. The obtained value of the residual entropy is very close to the one obtained by the Monte Carlo simulation.

[1] G. T. Barkema and M. E. Newman, Phys. Rev. E 57, 1155 (1998).

[2] G. Möller and R. Moessner, Phys. Rev. B 80, 140409(R) (2009).

[3] G.-W. Chern, P. Mellado, and O. Tchernyshyov, Phys. Rev. Lett. 106, 207202 (2011); G.-W. Chern and O. Tchernyshyov, Phil. Trans. R. Soc. A 370, 5718 (2012).