Trimer classical spin liquid from interacting fractional charges

TL;DR

This paper investigates a classical spin liquid phase arising from interacting fractional charges in a checkerboard lattice, revealing a degenerate state with potential for quantum spin liquid behavior and multiple magnetization plateaus.

Contribution

It introduces a novel classical spin liquid phase formed by fractional charges mapped to trimer coverings, expanding understanding of spin liquids in frustrated lattice systems.

Findings

Identification of a classical spin liquid phase with macroscopic entropy

Mapping of spin configurations to trimer coverings on a dual lattice

Observation of multiple magnetization plateaus indicating complex screening processes

Abstract

We study a problem of interacting fractional charges with $J_1$-$J_2$-$J_3$ Ising model on a checkerboard lattice under magnetic field. As a result of the interplay between repulsive interactions and particle density tuning by magnetic field, the fractional charges form a classical spin liquid (CSL) phase. The CSL phase is composed of degenerate spin configurations, which can be mapped to the trimer covering of dual square lattice. The CSL state shows macroscopic ground-state entropy, implying the emergence of novel quantum spin liquid phase when quantum fluctuations are turned on. In addition to the CSL phase, the system exhibits multiple magnetization plateaus, reflecting the fertile screening processes of dimer-monomer mixtures.

Tables

Table 1. Table 1: The width L 𝐿 L dependence of residual entropy per site, 𝒮 0 subscript 𝒮 0 \mathcal{S}_{0} .

Width: $L$	Residual entropy
2	0.14931329
3	0.15296086
4	0.13473800
5	0.13851784
6	0.13915249
7	0.13313547

Table 2. Table 2: Each number represents type of plaquettes in Fig.  14 (a). Combinations of a row are possible configuration for p , p ′ 𝑝 superscript 𝑝 ′ p,p^{\prime} and q 𝑞 q .

p	p’	q
(0,1,3)	(1,3,5)	9
	(6,8,9)	9
(2,4,5)	(0,2,4)	8
	(1,3,5)	(6,7)
	(6,8,9)	(6,7)
6	7	5
	(6,8,9)	(3,4)
(7,8,9)	(6,8,9)	(0,6,7)
	7	(1,2)
	(0,2,4)	8
	(1,3,5)	(6,7)

Equations

$$\mathcal{H}=$$

$$+J_3\sum _{⟨i,j{⟩}_{3rd}}{\sigma }_i^z{\sigma }_j^z-h\sum _i{\sigma }_i^z,$$

$$J_2=2J,J_3=J.$$

$$Q_p={\eta }_p{S}_p,$$

$$S_p=\sum _{i\in p}{\sigma }_i^z,$$

$$\displaystyle\eta_{p}=\{\begin{array}[]{ll}1&(\mathrm{for}\;A\;\mathrm{sub.})\\ -1&(\mathrm{for}\;B\;\mathrm{sub.})\end{array}.$$

$$\mathcal{H}=\left(\frac{1}{2}-J\right)\sum _p{Q}_p^2-J\sum _{⟨p,q⟩}{Q}_p{Q}_q-\frac{h}{2}\sum _p{\eta }_p{Q}_p+C_1,$$

$$C_1=2(J-1)N_p.$$

$$\sum _{p\in D}{Q}_p=\sum _{i\in \partial D}{\eta }_{p_{D(i)}}{\sigma }_i^z.$$

$$|\sum _{p\in D}{Q}_p|\le \sum _{i\in \partial D}|{\eta }_{p_{D(i)}}{\sigma }_i^z|=N_{\partial D}.$$

$$\mathcal{H}=\left(\frac{1}{2}-J\right)\sum _p{Q}_p^2-J\sum _{⟨p,q⟩}{Q}_p{Q}_q.$$

$$\mathcal{H}=\left(\frac{1}{2}+J\right)\sum _p{Q}_p^2+\frac{|J|}{2}\sum _{⟨p,q⟩}(Q_p+Q_q{)}^2.$$

$$\mathcal{H}=\left(\frac{1}{2}-3J\right)\sum _p{Q}_p^2+\frac{J}{2}\sum _{⟨p,q⟩}(Q_p-Q_q{)}^2.$$

$$\mathcal{H}=\left(\frac{1}{2}+J\right)\sum _p(Q_p-{\eta }_p{S}_h{)}^2-\frac{J}{2}\sum _{⟨p,q⟩}{\left(Q_p+Q_q\right)}^2+C_2,$$

$$S_h\equiv \frac{h/2}{1+2J},$$

$$\displaystyle Q_{p}=\{\begin{array}[]{ll}+2&\mathrm{for}\;A\ \mathrm{sub.}\\ -2&\mathrm{for}\;B\ \mathrm{sub.}\end{array}.$$

$$\displaystyle Q_{p}=\{\begin{array}[]{ll}+4&\mathrm{for}\;A\ \mathrm{sub.}\\ -4&\mathrm{for}\;B\ \mathrm{sub.}\end{array},$$

$$\displaystyle\Delta E=4\bigl{[}(1+2J)S_{h}-(1+6J)\bigr{]},$$

$$h=2+12J,$$

$$E_{\mathrm{tCSL}}=\frac{2}{3}(3J-h-1)N_p.$$

$$\mathcal{H}$$

$$-J\sum _{⟨p,q⟩}(Q_p-{\eta }_p{S}_h)(Q_q-{\eta }_q{S}_h)+C_3,$$

$$\displaystyle C_{3}=\frac{N_{p}}{8}\bigl{[}16J-16-\frac{h^{2}}{1+2J}\bigr{]}.$$

$$\mathcal{H}$$

$$+J{S}_h(S_h-2)n_{20}+J{S}_h^2{n}_{00}+J(S_h-2{)}^2{n}_{22}.$$

$$4N_0=2n_{00}+n_{20},$$

$$4N_2=2n_{22}+n_{20}.$$

$$N_0+N_2=N_p.$$

$$n_{00}+n_{22}+n_{20}=2N_p.$$

$$\mathcal{H}$$

$$\displaystyle+\bigl{[}\frac{1}{2}(\frac{1}{2}+J)(S_{h}-1)^{2}+\frac{1}{4}-\frac{3}{2}J\bigr{]}n_{20}$$

$$=a_{00}{n}_{00}+a_{22}{n}_{22}+a_{20}{n}_{20},$$

$$\displaystyle\mathrm{cluster\>}D_{+}:Q_{p}=\{\begin{array}[]{ll}+2&(\mathrm{for}\;A\;\mathrm{sub.})\\ 0&(\mathrm{for}\;B\;\mathrm{sub.})\end{array},$$

Full text

I Introduction

Fractionalization is one of the central topics in condensed matter physics. In many-body interacting systems, low-energy excitations are usually described as quasiparticles, whose nonuniversal parameters, such as mass, have renormalized values, while their quantum numbers are preserved to be the same as those of the original particles. Fractionalization changes this canonical description considerably, by allowing the particles to split into subunits with smaller quantum numbers.

Fractional excitations are hosted in a number of systems, such as one-dimensional quantum liquidBethe (1931); Mourigal et al. (2013); Ogata and Shiba (1990); Kitaev (2001), fractional quantum Hall systemsLaughlin (1983); Tsui et al. (1982), and quantum spin liquids (QSLs) Kitaev (2006); Knolle et al. (2014); Castelnovo et al. (2012); Balents (2010). Among them, in the candidate compounds of QSLs, the role of fractional excitations has been highlighted in their thermodynamic and transport behaviorsTokiwa et al. (2016); Pan et al. (2016); Ross et al. (2011); Kimura et al. (2013); Sibille et al. (2018); Udagawa and Moessner (2019); Wan et al. (2016); Kourtis and Castelnovo (2016); Huang et al. (2018); Chen (2017); Banerjee et al. (2016, 2017); Do et al. (2017); Kasahara et al. (2018); Knolle et al. (2014). In particular, the relevant role of fractional excitations is established in classical spin liquids (CSLs). CSL corresponds to a high-temperature precursor of QSL, composed of a macroscopic number of degenerate classical states. As lowering temperatures, coherency develop among the classical states, turning CSL into QSL. One can find several essential properties of QSL already show up in CSL. Among other things, fractional excitations can be defined in CSL, usually as static objects.

While the appearance of fractionalization itself requires interaction between original particles, the assembly of fractional particles also hosts nontrivial many-body problems. The problem of interacting fractional particles is usually quite difficult to treat theoretically, however, at the level of CSL, rigorous theoretical treatments are sometimes possible. In a class of CSLs defined on frustrated magnets, we can clearly divide the role of interactions: the nearest-neighbor interaction leads to the formation of CSL with emergent fractional excitations, while the farther-neighbor interactions give rise to the interaction between the fractional particles. A typical example can be found in spin ice, a material realization of three-dimensional Coulomb phase, where a fractionalized particle, called monopole, dominates the low-energy properties of the systemCastelnovo et al. (2008).

Indeed, in dipolar spin ice, monopoles exhibit interesting many-body effects, due to their Coulombic attractive interactions proportional to $-1/r$Castelnovo et al. (2008). This Coulomb attraction gives rise to a liquid-gas-type phase transition elusive in a magnetic systemSakakibara et al. (2003). This phase transition is controlled by magnetic field as a tuning parameter, which is translated into a chemical potential of monopoles. This analogy naturally leads to the plasma description of long-distance physics, where the screening of charge plays an important role in the thermodynamic behavior of the system.

It is also interesting to turn a look at fractional particles under short-range interactions. In this case, the system is more susceptible to the hardcore constraint arising from the local charge conservation. In a variant of spin-ice-type frustrated magnet, described by $J_{1}$-$J_{2}$-$J_{3}$ Ising model, again the nearest-neighbor interaction ($J_{1}$) leads to the formation of CSL, and the second- to third-neighbor interactions ($J_{2},J_{3}$) give rise to the short-range interactions between themMizoguchi et al. (2017). In this system, the fractional excitations have magnetic charges as well-defined quantum numbers, and they satisfy a conservation law. In conventional electromagnetism, same-sign charges repel each other. Meanwhile in this emergent Coulomb phase, same-sign charges sometimes attract each other. However, the conservation law strictly forbids the same-sign charges to pair-annihilate. This keen competition between the conservation law and the interactions leads to the formation of novel classical spin liquidsMizoguchi et al. (2017); Udagawa et al. (2016); Rau and Gingras (2016), excitations Mizoguchi et al. (2018); Yan et al. (2018), and dynamics Udagawa et al. (2016).

The stream of results motivates us to look into the magnetization process of this CSL. Rich behaviors of the magnetization process have been explored so far in magnetic systems, especially with geometrical frustrationOno et al. (2003); Miyata et al. (2017); Alicea et al. (2009); Coletta et al. (2016); Momoi and Totsuka (2000); Yamamoto et al. (2014); Sakai and Nakano (2011); Shirata et al. (2012); Susuki et al. (2013); Onizuka et al. (2000); Matsuda et al. (2013); Miyata et al. (2013); Nakamura et al. (2014). In general, the magnetic field affects the monopoles as a staggered potential and facilitates the pair-creation of monopoles. The interaction between the induced monopoles may lead to a formation of novel CSL with macroscopic degeneracy. In fact, a topologically ordered state is proposed at a low-field plateau of the kagome Heisenberg magnetNishimoto et al. (2013).

In this paper, we investigate a magnetic phase diagram and magnetization processes of the $J_{1}$-$J_{2}$-$J_{3}$ Ising model on a checkerboard lattice, which is a two-dimensional analog of the pyrochlore lattice. We obtain a ground state phase diagram analytically up to small positive farther-neighbor interactions, and find that there is a series of exotic phases which are characterized by magnetic plateaus. Remarkably, at the $1/3$-magnetization plateau, we find a CSL phase, where the whole system is tiled with magnetic trimers. This state is a new-type of CSL with novel value of residual entropy. It implies that a new class of QSL’s emerges upon adding perturbations to this CSL state.

The rest of this paper is organized as follows. In Sec. II, we present our model and methods, with a special focus on Gauss’ law, which is central to the rigorous argument we base the existence of the CSL on. In Sec. III, we show the overall ground-state phase diagram and present the results at simple limiting cases. In Sec. IV, we show the existence of the trimer classical spin liquid from both the intuitive and rigorous arguments, and detail its properties. In Sec. V, we consider a full magnetization process and show a variety of commensurate phases, which are characterized by magnetic plateaus at $M=0,1/5,1/3,1/2,2/3,3/4$ and $1$. Finally, in Sec. VI, we summarize our results and provide some future perspectives.

II Model and Method

II.1 Model

We consider the $J_{1}$-$J_{2}$-$J_{3}$ Ising model on a checkerboard lattice (Fig. 1). The checkerboard lattice is a corner-sharing network of square units with diagonal bonds, which we simply call “plaquettes”. We consider the network composed of $N_{p}$ plaquettes, i.e., $2N_{p}$ spins, in a periodic boundary condition in both directions. The Hamiltonian of our model is defined as

[TABLE]

where $\sigma^{z}_{i}=\pm 1$ are Ising spins and $h$ represents an external magnetic field. Throughout this paper, we consider the case of $h\geq 0$ without loss of generality. We assume the nearest-neighbor (n.n.) coupling as antiferromagnetic and set its value as a unit energy ($J_{1}$=1). The n.n. coupling accounts for the interaction within the plaquettes. In addition to $J_{1}$, we introduce $J_{2}$ as the diagonal interaction, and $J_{3}$ as the other next-nearest-neighbor interactions, as shown in Fig. 1. Throughout the paper, we fix the ratio of $J_{2}$ and $J_{3}$ as

[TABLE]

This point is of special importance, since the model can be mapped to the Hamiltonian of interacting magnetic charges with a staggered potential. To describe the mapping, we first introduce a magnetic charge:

[TABLE]

where $p$ is a label of plaquettes and we define the total spin of the plaquette $p$, as

[TABLE]

and the sign factor of the plaquette as

[TABLE]

Here, we took the dual picture and regarded the checkerboard lattice as a square lattice of plaquettes, and made the bipartite decomposition of the plaquettes into sublattice $A$ and $B$ as shown in Fig. 1. Note that there are five possible values of $Q_{p}$: $Q_{p}=0,\pm 2,\pm 4$ (Fig. 2). Among them, we especially call the plaquette with $Q_{p}=0$, a vacuum plaquette, and that with $|Q_{p}|=4$ as a double charge. Since the neighboring plaquettes share one spin, the charges are not independent of each other. Indeed, the magnetic charge has a conserved nature, and its distribution is a subject to the lattice analog of Gauss’ law, as we will show later.

In terms of $Q_{p}$, we can rewrite the Hamiltonian (1) as

[TABLE]

where

[TABLE]

We show the detail of this transformation in Appendix A. In Eq. (8), the first term is a “self-energy” of the charges which is proportional to $Q_{p}^{2}$. The second term is the interaction between magnetic charges on nearest-neighbor plaquettes. The coupling constant of this interaction is $J$. Accordingly, same-sign charges repel (attract) each other for $J<0$ ($J>0$). The third term is a staggered potential for charges, which arises from the magnetic field.

II.2 Gauss’ law

The magnetic charge, $Q_{p}$, satisfies the conservation law. To see this, suppose that $D$ is an arbitrary set of plaquettes and $\partial D$ denotes the boundary sites of $D$ (see Fig. 3 ). Then, the total magnetic charges inside $D$ and the spins on $\partial D$ satisfy the following relation:

[TABLE]

Here, the boundary site belongs to the two plaquettes, one inside, and one outside $D$, and $p_{D(i)}$ stands for the former. Equation (10) means that the total charge in a certain region equals the sum of “flux”, $\eta_{p_{D(i)}}\sigma^{z}_{i}$, on its boundary. This is nothing but a lattice analog of Gauss’ law.

This Gauss’ law constrains the structure of charge configuration in the system. In particular, as developed in Ref. Mizoguchi et al., 2017, this law leads to the triangle inequality:

[TABLE]

Namely, the amount of total charge an arbitrary region $D$ can store is bounded by $N_{\partial D}$, the number of its boundary sites.

III Phase diagram

Before going into the details of analyses, we first show the overall phase diagram of the present model in Fig. 4. Each phase corresponds to a magnetization plateau, which we characterize by the magnetization per spin, $M\equiv\frac{1}{2N_{p}}\sum_{i}\sigma_{i}^{z}$. Here, we focus on the region up to small positive $J$.

In this section, we consider the simple limiting cases, and first discuss the stability of the Coulomb phase at $h=0$, then consider the magnetization process for $J\leq 0$. We also give a brief introduction to the ground states for $J>1/4$ and $J<-1/2$.

III.1 Zero-field states

Let us start with the ground states in the absence of magnetic field. Without magnetic field, Hamiltonian (8) is simplified to

[TABLE]

III.1.1 $J=0$: square ice

For $J=0$, where there is no interaction between charges, the model is reduced to the square-ice model only with n.n. exchange interactions whose ground state is given by the charge vacuum ($Q_{p}=0$ for every $p$), or the two-dimensional Coulomb phase [Fig. 2(b)], which is a typical CSL state with macroscopic degeneracy. Lieb’s rigorous argument gives the exact value of residual entropy as $\mathcal{S}_{0}^{\rm SI}=\frac{1}{2}\ln W=0.215566427$ per spin, with $W=\frac{4}{3}^{\frac{3}{2}}$ Lieb (1967).

III.1.2 $J<0$: staggered charge ordering

For $J<0$, where an attractive interaction acts between opposite charges, the first term in Eq. (12) favors a vacuum, whereas the second term favors a staggered charge ordering. To see this competition clearly, we transform Eq. (12) into

[TABLE]

With this form, one can minimize the first and the second terms simultaneously, by setting $Q_{p}=-Q_{q}$ for $p\in$ A sub. and $q\in$ B sub., and $|Q_{p}|=0\ (4)$ for $J>-1/2$ ($J<-1/2$). This solution means the Coulomb phase extends to $J=-1/2$, while the staggered magnetic charge ordering with $|Q_{p}|=4$ takes over for $J<-1/2$. In the spin language, the staggered charge ordering with $|Q_{p}|=4$ corresponds to the fully-polarized ferromagnetic state [Fig. 2(c)].

III.1.3 $J>0$: long-period phase

In contrast to the case of $J<0$, same-sign charges attract with each other for $J>0$. As a result, the Coulomb phase survives only up to $J=1/4$, and then the system turns to a complicated long-period spin ordering for $J>1/4$, as in the case of pyrochlore modelUdagawa et al. (2016).

To obtain an insight into the ground state for positive $J$, let us rewrite the Hamiltonian Eq. (12) as

[TABLE]

This expression shows the Coulomb phase is stable at least up to $J=1/6$ from a similar argument of $J<0$. It also implies proliferation of double charges at larger $J$. At first sight, it is preferable to cover the system with same signs of double charges for large $J$, however, according to the Gauss’ law, the same-sign-charge cluster accommodates at most one double charge Mizoguchi et al. (2017). As a result, the ordered phase contains mixed values of charges as shown in Fig. 2(d). The precise value of the critical point of $J=1/4$ can be obtained by the softening of the lowest excitations, in the same argument as in pyrochlore model Udagawa et al. (2016). The formal proof based on the Gauss’ law is rather involved. We summarize it in Appendix B.

III.2 Magnetization process at $J\leq 0$

Next, let us look at the magnetization process for $J\leq 0$. In this case, simple and rigorous arguments are available for the magnetization process, including the simplest non-interacting limit, $J=0$. In this region, it is convenient to rewrite the Hamiltonian (8) as

[TABLE]

with C_{2}=\frac{N_{p}}{8}\bigl{[}16J-16-\frac{h^{2}}{1+2J}\bigr{]}. Here, we define

[TABLE]

Then, in analogy with the argument in Sec. III.1.2, the ground state can be obtained by minimizing the two terms in Eq. (15), simultaneously. It is possible by minimizing $|Q_{p}-\eta_{p}S_{h}|$ while keeping the staggered charge alignment $Q_{p}+Q_{q}=0$ for any neighboring plaquettes, $p$ and $q$.

III.2.1 Square ice: $0<h<2+4J$

For small $h$, we find that the zero-field square ice, or the Coulomb phase, extends in the region: $S_{h}<1$ ($h<2+4J$).

III.2.2 Dimer phase: $2+4J<h<6+12J$

For $S_{h}\geq 1$, the system goes out of the Coulomb phase, and the staggered charge ordering appears. This region is divided into two cases. First, for $1<S_{h}<3$, (i.e., $2+4J<h<6+12J$), $Q_{p}$ is given by

[TABLE]

This phase corresponds to the half-magnetization plateau with $M=1/2$. While this phase has the staggered charge ordering, it still keeps macroscopic degeneracy in spin degrees of freedom, forming a CSL state. Its ground-state manifold can be mapped to the dimer covering problem, in a similar way to kagome iceMatsuhira et al. (2002); Udagawa et al. (2002); Isakov et al. (2004). To see this, let us consider a dual square lattice again, and place a dimer on each down spin. Then, each spin configuration on this half-magnetization plateau can be mapped to a dimer configuration on the dual square lattice (Fig. 5). The residual entropy can be exactly obtained to be $\mathcal{S}_{0}^{\rm hmp}=0.14578045$ per spin by counting the dimer configuration on a square lattice Fisher (1961).

III.2.3 The fully-polarized state: $6+12J<h$

Finally, for $3<S_{h}$, ($6+12J<h$), $Q_{p}$ is given by

[TABLE]

which corresponds to a fully-polarized ferromagnet.

IV Trimer CSL

At $J=0$, the system makes a direct transition from the zero-field square ice phase to the dimer phase at $h=2$. At $J>0$, between these two CSL’s, there appears another type of CSL, which we call “trimer CSL” after its structure. In this section, we will address how the dimer phase makes instability to this trimer CSL, by the two kinds of strategies: one is based on the instability analysis due to the creation of vacuum plaquette, which is more intuitive and gives a clearer picture of this trimer CSL, and the other is a rigorous argument based on Gauss’ law.

IV.1 Instability of the dimer phase

IV.1.1 Nucleation of vacuum plaquette

We start with the instability analysis of the dimer phase, as the magnetic field decreases. Here, we adopt the magnetic charge representation of the Hamiltonian, Eq. (8), and estimate the critical magnetic field. At half-magnetization plateau, we have $|Q_{p}|=2$, for all the plaquettes. As the magnetic field decreases, we expect the vacuum plaquettes to nucleate.

To examine this process, we take one upward spin and flip it downward, and create a pair of vacuum plaquettes. The nucleated plaquettes are dissociated, and then they are individually screened by the charged plaquettes [Fig. 6 (a)], to maximize the energy gain from the interaction term. From Eq. (15), we can estimate the energy increase associated with this process:

[TABLE]

which leads to the instability at $S_{h}=\frac{1+6J}{1+2J}$, or

[TABLE]

below which the dimer phase becomes unstable against the creation of vacuum plaquettes.

IV.1.2 Formation of the trimer CSL

Next, let us look into the state just below the critical field. In this state, all the vacuum plaquettes are surrounded by four charged plaquettes with $|Q_{p}|=2$. See the structure inside the magenta square in Fig. 6(a). Since this structure locally optimizes the energy, it is desirable to tile the whole lattice with as many of this local structure as possible. However, it is not straightforward to obtain the optimal tiling pattern.

To gain an insight, let us adopt a dimer representation we have introduced in Sec. III.2, and place a dimer on each downward spin [Fig. 6(b)]. From this viewpoint, a vacuum plaquette can be regarded as an overlapping part of two dimers. Then the tiling problem can be interpreted as fully packing the dual square lattice with overlapping dimers, or “trimers”, under the condition that the overlapping part (i.e., a vacuum plaquette) does not neighbor with each other (for an example of prohibited structure, see Fig. 7). The packing with trimers is in sharp contrast to the dimer covering problem that appears on the half-magnetization plateau. As we numerically estimate in the next subsection, the number of possible trimer configurations increases macroscopically as system size, suggesting that the state forms a CSL. Accordingly, we name this CSL as “trimer classical spin liquid” (tCSL) after its structure.

IV.1.3 Characterization of trimer CSL

In this section, we will look into the detailed character of tCSL. Firstly, the formation of tCSL is associated with the magnetization plateau. The picture of trimer covering immediately leads to the magnetization value $M=1/3$ per spin, since one out of the three plaquettes forming a trimer has the total spin $S_{p}=0$, while the other two have $S_{p}=2$. Comparing with the fully polarized ferromagnet, where $S_{p}=4$ for all the plaquettes, we can easily obtain $M=(0+2+2)/(4\cdot 3)=1/3$.

Secondly, what is the ground state energy of tCSL? To see this, again from the structure of a trimer, the number of vacuum plaquettes is $N_{0}=N_{p}/3$, while that of charged plaquette is $N_{2}=2N_{p}/3$. Since each vacuum plaquette contacts with four charged plaquettes, the number of bonds connecting the vacuum and charged plaquettes is $4N_{0}=\frac{4}{3}N_{p}$. Since the total number of bonds is $2N_{p}$, and no bonds connect two vacuum plaquettes, the number of bonds connecting the two charged plaquettes is $2N_{p}-\frac{4}{3}N_{p}=\frac{2}{3}N_{p}$. Inputting all these into Eq. (8), we obtain the energy of tCSL:

[TABLE]

Thirdly, how large is the ground state degeneracy? To gain an insight into the origin of degeneracy, let us look at one typical trimer configuration in Fig. 8. Looking at this configuration, one can find a system is divided into large-scale charge clusters. In a positive (negative) charge cluster, the charged plaquettes are placed only on A (B) sublattice, and the vacuum plaquettes are placed on B (A) sublattices. Inside one cluster, one finds a staggering pattern of charged and vacuum plaquettes. Only at the cluster boundaries, the charged plaquettes neighbor with each other. Meanwhile the vacuum plaquettes are never adjacent to each other.

This cluster structure implies there are two origins for the ground state degeneracy: the contribution from cluster placements, and the internal spin configuration within a cluster. The existence of these two types of configurational entropy is in common with the hexamer CSL found in the cousin system of $J_{1}$-$J_{2}$-$J_{3}$ Ising model on the Kagome lattice, where a novel value of residual entropy was found, as well as a characteristic spin correlation with iconic half-moon pattern in magnetic structure factorMizoguchi et al. (2017).

Here, in order to obtain the precise value of residual entropy, we resort to a transfer matrix method, by adopting a finite strip of checkerboard lattice with variable widths $L$, which is the number of plaquettes in the column of the strip (see Fig. 1), up to 7. The detail of transfer matrix method is summarized in Appendix C. We listed the values of residual entropy $\mathcal{S}_{0}$ for each width on Table 1. While the obtained value still fluctuates around $L\sim 7$, they are well within the window of $\mathcal{S}_{0}=0.136\pm 0.04$, giving the evidence for the existence of finite residual entropy. The value $\mathcal{S}_{0}\simeq 0.136$ is smaller compared with the zero-field value, $\mathcal{S}_{0}^{\rm SI}=0.21556643$, but comparable to that of the dimer phase, $\mathcal{S}_{0}^{\rm hmp}=0.14578045$.

Previously, the configurational entropy was obtained for a trimer covering problem in several contextsFroböse et al. (1996); Ghosh et al. (2007); Lee et al. (2017); Dong et al. (2018). In Ref. Froböse et al., 1996, only the angular trimers are considered. The obtained value is, if translated into the current context, $\mathcal{S}_{0}^{\rm ang}\sim 0.13846575$ per spin, comparable to our value of $\mathcal{S}_{0}$. At first sight, it may seem strange that the limitation to angular trimers do not make difference in larger residual entropy. However, their model allows the touching of vacuum plaquette in our context, which seems to compensate for the entropy reduction due to the limitation of the type of trimers. In Ref. Ghosh et al., 2007, on the other hand, the authors only discuss the covering by the line trimers, and they obtained $\mathcal{S}_{0}^{\rm line}\sim 0.07926$. The covering problems by all species of trimers are considered in Ref. Dong et al., 2018. They obtained $\mathcal{S}_{0}^{\rm tot}\sim 0.20597$. The difference between $\mathcal{S}_{0}^{\rm tot}$ and our $\mathcal{S}_{0}$ is attributed to the prohibition of vacuum touching.

The finite residual entropy qualifies this phase as CSL. This trimer CSL is similar to the classical spin liquids found in kagome Mizoguchi et al. (2017) and pyrochlore Udagawa et al. (2016) lattices previously in a sense that the same-charge attraction leads to these CSL’s. In these systems, one spin is fractionalized into a pair of magnetic charges, and the fractionalized charges make recombination into a new phase that cannot be easily inferred from the original spin degrees of freedom.

While the trimer CSL has some similarities to the previous CSL’s as mentioned above, it shows a unique feature in its structural rule: the trimer covering. To the knowledge of the authors, the trimer model does not belong to known integrable models, implying a possibility that this trimer CSL may serve as a precursor of a new type of QSL.

In this light, it is interesting to introduce quantum fluctuation to this system, in the form of, e.g., quantum transverse or exchange interactions. If QSL is actually realized on the basis of this trimer covering state, we expect it to show an anomalous magnetic charge correlation. In a trimer CSL, magnetic charges tend to make clusters, as mentioned above. If this tendency survives after the superposition of trimer configurations, it will present strong short-range charge correlation, quite distinct from the previously known class of QSL’s.

IV.2 Rigorous argument based on Gauss’ law

In the previous subsection, we derived the structure of tCSL in an intuitive way, from the instability analysis of half-magnetization plateau. Here, we rigorously show that the tCSL state gives the ground state, on the basis of Gauss’ law we introduced as Eq. (10). To this aim, we rewrite the Hamiltonian Eq. (8) in a form,

[TABLE]

where

[TABLE]

Here, we limit ourselves to small positive $J$, and focus on the magnetic field just below the half-magnetization plateau, i.e., $h\sim 2$. In the range of magnetic field under consideration, we can safely assume that all the plaquettes are occupied by 2-up 2-down or 3-up 1-down spin configuration, namely, $|Q_{p}|=0$ or $2$, and $Q_{p}=+2$ $(-2)$ exists only at the $A$ ($B$) sublattice. Under this assumption, let us define the number of charges with $|Q_{p}|=0\ (2)$ as $N_{0}$ ($N_{2}$), and the number of contacts between the plaquettes with $|Q_{p}|=0$ and $|Q_{q}|=2$ as $n_{20}$. Similarly, we define $n_{00}$ ($n_{22}$) as the number of contacts between two plaquettes with charge [math] ($2$). With these quantities, we can express the Hamiltonian, Eq. (26) as

[TABLE]

The variables in this Hamiltonian are subject to several geometrical constraints. Firstly, since each plaquette has four contacts with neighboring plaquettes on a dual square lattice, we have

[TABLE]

and

[TABLE]

Here, the factor 2 before $n_{00}$ and $n_{22}$ correct the double counting. Secondly, in the absence of double charges ($|Q_{p}|=4$), the sum of $N_{0}$ and $N_{2}$ is equal to the total number of plaquettes:

[TABLE]

Combining Eqs. (29)-(31), the contact numbers satisfy

[TABLE]

With Eqs. (29) and (30), we can eliminate $N_{0}$ and $N_{2}$ from the Hamiltonian (28), and obtain

[TABLE]

where $a_{00}$, $a_{22}$ and $a_{20}$ are coefficients of $n_{00}$, $n_{22}$ and $n_{20}$. Accordingly, the search for the ground state is now reduced to finding the combination of $(n_{00},n_{22},n_{20})$ to minimize $\mathcal{H}$, under the constraint of sum rule, Eq. (32).

The coefficients $a_{00}$, $a_{22}$ and $a_{20}$ are plotted in Fig. 9. For $h>2+12J$, $a_{22}$ is the smallest. This region corresponds to the half-magnetization plateau, where all the plaquettes are occupied with charges $|Q_{p}|=2$, and all the contacts are of $2$-$2$ type, accordingly.

Meanwhile, this plot shows that larger $n_{20}$ is preferable for $2-4J<h<2+12J$. To satisfy this condition, at first sight, the best strategy seems to put $Q_{p}=2$ on all the plaquettes of A sublattice, while $Q_{p}=0$ on B sublattice, to make all the contacts to be of 2-0 type. However, this charge configuration obviously violates the Gauss’ law, Eq. (10).

To find the optimal charge configuration under the Gauss’ constraint, let us define a positive (negative) charge cluster $D_{+}$ ($D_{-}$), as a maximal set of plaquettes with $Q_{p}=2$ placed on the A(B)-sublattice and those with $Q_{p}=0$ placed on the B(A)-sublattice [Fig. 8 (a)]. Namely,

[TABLE]

With these definitions, clusters $D_{+}$ and $D_{-}$ always touch with each other through the 0-0 or 2-2 contacts. To see this, suppose a cluster $D_{+}$ has a boundary plaquette that belongs to A (B) sublattice, then it must have charge $|Q_{p}|=2$ (0). This plaquette neighbors with a plaquette of $D_{-}$ on the B (A) sublattice, which has charge $|Q_{p}|=2$ (0).

Now, let us apply the Gauss’ law to a cluster, $D_{\alpha}$ of either type. We assume that the number of charged plaquettes inside $D_{\alpha}$ to be $N_{2}^{\alpha}$, and define the number of 0-0 and 2-2 contacts with neighboring clusters to be $n_{00}^{\alpha}$ and $n_{22}^{\alpha}$, respectively. The Gauss’ inequality, Eq. (11) leads to

[TABLE]

By summing over all the clusters in the system, we obtain

[TABLE]

Note that the factor 2 of the right-hand side comes from the double counting of bonds in the summation over clusters. Combining Eqs. (30) and (41), we obtain

[TABLE]

Now, considering the relative magnitudes of coefficients depicted in Fig. 9, at $h\lesssim 2+12J$, Eq. (42) results in the optimal solution to be $n_{00}=0$ and $n_{20}=2n_{22}$. If a certain configuration satisfies this condition, it gives one of the ground states. In fact, this condition is equivalent to the trimer covering we discussed in the previous subsection. By inserting this condition into Hamiltonian (28) with the constant term $C_{3}$ given by Eq. (27), we obtain

[TABLE]

which exactly corresponds to what we obtained from the trimer covering picture: Eq. (25). This means that the tCSL states give the ground state.

Conversely, it is possible to show that any member of the ground state manifold can be expressed by the trimer covering, i.e., the ground state is composed only of the tCSL states. To prove this, it is enough to show that if $n_{00}=0$ and $n_{20}=2n_{22}$ are satisfied, the corresponding charge configuration can be expressed in terms of the trimer covering.

To see this, suppose a cluster $D_{\alpha}\in D_{+}$. To satisfy the former condition, one needs $n_{00}^{\alpha}=0$, and in addition to that, the equality must hold in the Gauss’ inequality, Eq. (40):

[TABLE]

This equality requires the boundary spins to satisfy $\sigma_{i}^{z}=+1$, and the cluster $D_{\alpha}$ to neighbor with other clusters with charged plaquette. Accordingly, within the cluster $D_{\alpha}$, each charged plaquette shares its only down spin with its neighboring vacuum plaquette. Moreover, each vacuum plaquette shares its two down spins with two of its neighboring charged plaquettes. Consequently, if one places a dimer on each down spin, all the charged plaquettes are covered with one dimer, and all the vacuum plaquettes are covered with two dimers, resulting in a trimer covering. The same argument holds for a cluster $D_{\alpha}\in D_{-}$.

V Full magnetization process

In this section, we will address the rest of magnetic phase diagram shown in Fig. 4 for $J>0$. We limit ourselves to the region of small $J$, again. We show the schematic picture of spin configurations at each phase in Fig. 10, and the magnetization process in Fig. 11.

V.1 Instability of high-field phase

We start with the instability of the high-field fully-polarized phase [Fig. 10(a)], with decreasing $h$. In the fully-polarized phase, all the plaquettes have total spins $S_{p}=4$. A single spin flip makes an adjacent pair of plaquettes with $S_{p}=2$. The energy increase accompanying this process is

[TABLE]

This instability occurs for $\Delta E<0$, i.e.

[TABLE]

Below this boundary line, the system tries to maximize the number of neighboring pairs of plaquettes with $S_{p}=4$ and $S_{p}=2$. Consequently, we obtain the 4-plaquette ordering as shown in Fig. 10(b). This ordered phase corresponds to the $M=3/4$ plateau.

V.2 High-field instability of dimer phase

In the previous section, we addressed the low-field instability of half-magnetization plateau into tCSL. Here, we address the instability with an increasing magnetic field. In the half-magnetization plateau, all the plaquettes have $S_{p}=2$ uniformly. As the magnetic field increases, we expect a nucleation of double charge, $S_{p}=4$. The nucleation process takes qualitatively different form from the nucleation of $S_{p}=2$ plaquette from the fully-polarized phase.

In the fully-polarized phase, the nucleated $S_{p}=2$ plaquettes are always paired. In contrast, in the half-magnetization plateau, the nucleation of $S_{p}=4$ plaquettes occurs in pair, but they can be dissociated from each other. This is a sort of fractionalization, which reflects the deconfining nature of the dimer phase at the half-magnetization plateau.

As a result, the most economical excitation is a single plaquette of $S_{p}=4$ surrounded by $S_{p}=2$ plaquettes. This object costs energy

[TABLE]

which becomes negative, if

[TABLE]

Above this field, the spin configuration, as shown in Fig. 10(c), is stabilized.

If we stand on the dimer picture by placing a dimer on a down spin, the resultant state consists of an assembly of monomers screened by dimers. It follows naturally that the corresponding state makes a $M=2/3$ plateau. Most simply, three-plaquette ordering realizes this state [Fig. 10(c)], however, introduction of “stacking fault” does not increase the energy. Accordingly, this state shows a semi-macroscopic degeneracy of the order of $2^{L}$.

Obviously, this $M=2/3$ plateau cannot be continuously connected to the 4-plaquette ordering of $M=3/4$ plateau, just below the saturated phase. The transition between the two states occurs at $h=6+10J$, from a simple comparison of the energies.

V.3 Instability of square ice

Finally, we address the instability of square ice state, as increasing $h$ from 0. To this aim, rather than the dimer representation, we resort to Gauss’ inequality, following the argument in Sec. IV.2.

We start with Hamiltonian (33). From their coefficients as shown in Fig. 9, it follows that $n_{22}$ should be suppressed for smaller magnetic field. Meanwhile larger $n_{20}$ is preferable for $h>2-4J$.

Given that $n_{22}=0$, a set of geometrical equations, Eqs. (29), (30), and (32) lead to

[TABLE]

The optimal value of $n_{20}$ is constrained by the Gauss’ inequality, Eq. (41) as

[TABLE]

which leads to the inequality,

[TABLE]

combined with the geometrical constraint, Eq. (51).

The optimal value of $n_{20}$ corresponds to the equality of Eq. (53), which leads to

[TABLE]

and

[TABLE]

This state corresponds to the $M=1/5$ plateau, and we show the configuration in Fig. 10(f).

VI Summary

We have studied a problem of interacting fractional charges, taking the $J_{1}$-$J_{2}$-$J_{3}$ Ising model on a checkerboard lattice under the magnetic field. We focused on the case, where the Hamiltonian can be written in terms of the charge degrees of freedom. In particular, at small positive $J$, we found that the half-magnetization plateau destabilizes into a classical spin liquid state, as decreasing the magnetic field. The resultant CSL is expressed as an assembly of trimer covering of the dual square lattice, and we called this state the tCSL.

The tCSL state corresponds to the $1/3$ magnetization plateau, and has macroscopic ground state degeneracy, which is characterized by a novel value of residual entropy. In contrast to dimer covering, which is ubiquitous in a broad area of physics including statistical mechanics and condensed matter physics, the notion of trimer covering rarely appears. It is surprising that such elusive states can be obtained, by starting from a simple local Hamiltonian considered here. As turning on a quantum fluctuation, a novel quantum spin liquid state may further be stabilized based on a well-defined microscopic model.

Moreover, we showed the interactions among magnetic charges lead to a variety of magnetization plateaus in the applied magnetic field, reflecting the rich screening processes of dimer-monomer mixtures. Nontrivial magnetization processes are observed for a number of frustrated magnetic systems. In this regard, this work shows that the picture of interacting fractional charges gives a new viewpoint to the formation of magnetization plateaus.

Acknowledgements.

This work was supported by the JSPS KAKENHI (Nos. JP15H05852, JP16H04026, and JP17H06138), MEXT, Japan. Part of numerical calculations were carried out on the Supercomputer Center at Institute for Solid State Physics, University of Tokyo. K. T. was supported by the Japan Society for the Promotion of Science through the Program for Leading Graduate Schools (MERIT).

Appendix A Mapping to charge Hamiltonian

In this appendix, we derive the charge representation of Hamiltonian Eq. (8) from the original Hamiltonian Eq. (1). First, we rewrite Eq. (1) with $(J_{1},J_{2},J_{3})=(1,2J,J)$ as

[TABLE]

where we name the spins on the neighboring plaquettes, $p$ and $q$, as $p_{1}$-$p_{4}$ and $q_{1}$-$q_{4}$ respectively, as shown in Fig. 12. The symbol $[p,q]$ stands for the spin shared by $p$ and $q$. For instance, $[p,q]=p_{4}=q_{1}$ in Fig. 12. The first and the second term come from the nearest-neighbor interactions while the third term comes from the second-neighbor and the third-neighbor interactions which connect spins on neighboring plaquettes. The fourth term comes from magnetic field. We transform Eq. (56) by introducing the total spin of the plaquettes $S_{p}=\sum_{i\in p}\sigma^{z}_{i}$:

[TABLE]

and

[TABLE]

By introducing charge representation $Q_{p}=\eta_{p}S_{p}$, finally we obtain Eq. (8):

[TABLE]

where $C_{1}=2(J-1)N_{p}$.

Appendix B Ordered phase at $h=0$ for $J>1/4$.

In this appendix, we derive the phase boundary for $h=0$ and $J>0$, by evaluating the ground-state energy using the constraint of the Gauss’ law and some geometrical identities. We show the Coulomb phase survives up to $J=1/4$, and the ordered phase in Fig. 13 realizes for $J>1/4$.

The outline of the energy evaluation is as follows. First, we decompose an arbitrary charge configuration into assembly of maximally connected same-sign-charge clusters and vacuum plaquettes. Here the maximally connected same-sign-charge cluster means that it must be connected to opposite-sign-charges or vacuum plaquettes at the boundary. By doing so, the total energy is written by the sum of the intra-cluster energy and the inter-cluster energy. Then, we optimize the energy under the geometrical constraint and the Gauss’ law, to obtain the ground-state charge configuration.

To begin with, let us address an essential property of the same-sign-charge cluster (labeled as $\mathcal{D}$), which originates from the Gauss’ law. Namely, applying the triangle inequality Eq. (11) to $\mathcal{D}$, we have

[TABLE]

where $N_{2}^{(\mathcal{D})}$ ($N_{4}^{(\mathcal{D})}$) is a number of $|Q_{p}|=2$ ($|Q_{p}|=4$) plaquettes in $\mathcal{D}$, and $n_{b}^{(\mathcal{D})}$ denotes the number of boundary sites of $\mathcal{D}$.

Having this at hand, we now evaluate the ground state energy. To do this, we introduce the following variables. (1) $n_{i:(\ell.\ell^{\prime})}^{(\mathcal{D})}$: a number of inner sites of $\mathcal{D}$-clusters which is contact points of $|Q_{p}|=\ell$ and $|Q_{p}|=\ell^{\prime}$. (2) $n_{b:(\ell.\ell^{\prime})}^{(\mathcal{D})}$: a number of boundary sites of $\mathcal{D}$-clusters at which a plaquette with $|Q_{p}|=\ell$ inside the cluster confronts another plaquette with $|Q_{p}|=\ell^{\prime}$ outside the cluster. See Fig. 13 for examples of sites counted as $n_{i:(2.4)}^{(\mathcal{D})}$ and $n_{b:(2.4)}^{(\mathcal{D})}$. We note that the total number of inner sites $n_{i}^{(\mathcal{D})}$ and that of boundary sites $n_{b}^{(\mathcal{D})}$ satisfy

[TABLE]

Note that $n_{i:(4.4)}^{(\mathcal{D})}=0$, since the same-sign charges with $|Q_{p}|=4$ cannot touch with each other. Using these valuables, the total energy can be written as

[TABLE]

Here, the terms proportional to $n_{i}^{(\mathcal{D})}$ ($n_{b}^{(\mathcal{D})}$) correspond to intra-cluster (inter-cluster) energy. To simplify Eq. (63) by erasing $n_{i:(2.2)}^{(\mathcal{D})}$, we substitute the following identities,

[TABLE]

into Eq. (63), and then we obtain

[TABLE]

where $N_{p}=\sum_{\mathcal{D}}N_{4}^{(\mathcal{D})}+\sum_{\mathcal{D}}N_{2}^{(\mathcal{D})}+N^{\mathrm{vac}}$ is a total number of plaquettes and $N^{\mathrm{vac}}$ is a number of vacuum plaquettes. Applying the inequality of Eq. (61) to the right side of Eq. (66) and using Eq. (65), we obtain

[TABLE]

For $0<J<1/4$, all the coefficients of $n^{(\mathcal{D})}_{i:(\ell,\ell^{\prime})}$, $n^{(\mathcal{D})}_{b:(\ell,\ell^{\prime})}$ and $\left(N_{p}-N^{\mathrm{vac}}\right)$ in the right side of Eq. (67) are positive, thus the energy minimization is achieved when $n^{(\mathcal{D})}_{i:(\ell,\ell^{\prime})}=n^{(\mathcal{D})}_{b:(\ell,\ell^{\prime})}=\left(N_{p}-N^{\mathrm{vac}}\right)=0$, which is nothing but the Coulomb phase. For $J>1/4$, on the other hand, the coefficients of $n^{(\mathcal{D})}_{i:(2.4)}$ and $\left(N_{p}-N^{\mathrm{vac}}\right)$ become negative, indicating the instability of the Coulomb phase.

How can we derive the ground state for $J>1/4$? To do this, let us first derive the upper bound of $n_{i:(2.4)}^{(\mathcal{D})}$. To this end, we point out that the following identity holds:

[TABLE]

Then, combining Eqs. (61), (62), (65), and (68), after some of algebra, we obtain

[TABLE]

Combining Eqs. (69) and (67), we have

[TABLE]

Further, we can derive the lower bound of $N^{\mathrm{vac}}$ as

[TABLE]

where $N_{i:(0.0)}$ is the number of contacts between vacuum plaquettes. Combining Eqs. (71) and (70), we obtain

[TABLE]

Now, we see from Eq. (72) that the energy minimization is achieved when $n_{b:(2.2)}^{(\mathcal{D})}=n_{b:(2.4)}^{(\mathcal{D})}=n_{b:(4.0)}^{(\mathcal{D})}=n_{b:(4.2)}^{(\mathcal{D})}=n_{b:(4.4)}^{(\mathcal{D})}=0$, and the corresponding minimum energy is

[TABLE]

In this case, the equalities in Eqs. (61), (69) and (71) hold, i.e.,

[TABLE]

In fact, the ordered phase in Fig. 13 satisfied Eqs. (73)-(76). Therefore it is the ground state spin configurations for $J>1/4$.

Appendix C Transfer matrix method

In this appendix, we summarize the transfer matrix method we used for the estimation of residual entropy of tCSL. We consider the lattice geometry shown in Fig. 14(b). The lattice is composed of $LN$ plaquettes with $L$ rows and $N$ columns, and we impose periodic boundary conditions on both directions. In this method, we obtain the total number of trimer patterns by identifying the possible trimer configurations column by column.

Specifically, one plaquette can be covered by a trimer in ten different ways, as shown in Fig. 14(a). Let us make a numbering of all possible trimer configurations on the $i-$th column, and denote it as $\phi_{i}$. Since there are ten possible configurations for each plaquette, and each column consists of $L$ plaquettes, $\phi_{i}$ takes $0,1,\cdots 10^{L}-1$. However, the trimer configuration on $i$-th column restricts the possible configuration on the $(i+1)$-th column. As shown in Fig. 14(c), type $q$ trimer on a plaquette on the $i$-th column restricts the possible trimer types on adjacent plaquettes on the $(i+1)$-th column. Possible combinations of trimer types on a triplet of plaquettes, ($p,p^{\prime},q$) is summarized in a table 2.

We account for this restriction by introducing the transfer matrix $A^{(L)}$, whose element, $A_{\phi_{i},\phi_{i+1}}^{(L)}$, takes 1, if the configurations $\phi_{i}$ and $\phi_{i+1}$ are compatible, and 0, otherwise. In terms of this transfer matrix, the total number of trimer configurations, $W^{(L,N)}$, can be obtained as

[TABLE]

for sufficiently large $N$. Here, $\lambda_{L}$ is the largest eigenvalue of $A^{(L)}$. By taking $N\to\infty$, the residual entropy per site is evaluated as

[TABLE]

We numerically evaluate $\lambda_{L}$ up to $L=7$ and obtain the residual entropy as listed in table 1.