Hartel R D Flick

Numerical model of melting, liquid migration and deformation during an ice cream meltdown test

A numerical model was developed which describes the combined phenomena occurring during an ice cream meltdown test which is a test used by ice cream producers to estimate the quality of ice cream. It is based on the measurement of the drip loss and the deformation during melting of an ice cream block placed onto a grid at ambient temperature. This information is related to the ice cream structure and reflects the quality of the product: even partially melted, the ice cream structure should retain the liquid phase and show low deformation. The proposed model takes into account melting, liquid migration inside the ice cream structure (which is considered as an unsaturated porous medium) and deformation. The objective in developing the model is to better understand the mechanisms of the ice cream structure evolution when it is submitted to ambient conditions.

To describe the multiphysics phenomena first, energy equation is solved to describe temperature evolution of the sample. Ice crystal melting was computed by using the enthalpic approach. After partial melting, water movement inside ice cream structure is predicted by Richard's equation. The hydraulic conductivity, the water retention characteristic and the differential moisture capacity are the parameters that control water movement inside the ice cream. Deformation of the partially melted ice cream was described by an elasticity model, Youngs' modulus being function of the fraction of melted ice. The set of equations was solved using the CFD finite element method implemented in Comsol software. The numerical model was validated with experimental results of temperature, drip loss and deformation. The model allows determining the main parameters, related to product structure, controlling the drip loss.

INTRODUCTION

During ice cream process, destabilisation of the fat emulsion occurs as a result of process conditions and mix formula. The whipping and freezing process transform some of the fat in the mix (originally very small <1µm) into three dimensional fat aggregates giving a certain microstructure; this phenomena is called fat destabilisation. According to Goff D. & Hartel (2013), fat aggregates play a very important role on the slowness of melting and shape retention during melting. In some cases fat destabilisation is needed, when the product should have some particular structural characteristics, for example in the case of extruded ice cream products. In other words, the process needs to be optimised as a function of the end product (not enough destabilisation would lead to poor shape retention whereas too much destabilisation would lead to the formation of too visible fat aggregates). The ice cream does not collapse although ice melts because water is retained in the microstructure. Fat destabilisation is promoted by several parameters that should be controlled: the total amount of fat and composition of the mix (type of fat, emulsifiers, etc) but also the process conditions such as the speed of the dasher that generates the shearing forces and whipping conditions that generate fat destabilisation.

Ice cream meltdown properties are important for two main reasons (Goff D. & Hartel 2013): melting in the mouth during consumption and meltdown related to storage at warm temperatures. Little work has been done to characterise melting of ice cream or frozen desserts and the relation with the sensory attributes during consumption. The meltdown test or drip-through test has been used by industrial to estimate the microstructure of their products. Normally this test is performed at room temperature. As heat penetrates into the ice cream and ice crystals melt, the water from melted ice dilute the unfrozen mix that drains down through the structure formed by bubbles, fat agglomerates and globules; then the mix diluted leaves the sample through a screen. Some ice cream meltdown collapse and flow totally through the screen. Other ice cream collapse only very slightly, leaving a nearly structure remaining on the screen. The internal structure and interconnectivity inside the sample influences the ability of the serum to flow inside the structure [START_REF] Hartel | Effects of Structural Attributes on Hardness and Melting Rate of Ice Cream[END_REF], [START_REF] Marshall | Effects of Overrun on Structural and Physical Properties of Ice Cream[END_REF]).

Very few work in literature deals with modelling multi-physical phenomena during ice cream meltdown. Several interconnected governing phenomena occurs at the same time: heat transfer, phase change, fluid flow through the structure and mechanical phenomena related to the elasticity of the structure formed during freezing. Modelling meltdown coupled multiphysic would allow a better understanding of the microstructure and the process influence on it.

MATERIALS AND METHODS

Mix

Mix of ice cream was prepared under typical conditions, mixing of the ingredients, pasteurization, homogenization and cooling. Then ice cream mix was stored at 5°C for 24 hrs allowing maturation process. The ingredients of the mix are: milk, dry milk, milk fat, sweetening agents, stabiliser and emulsifier. With a Brix of 26°, the milk mix density was about 1130 kg/m 3 and its initial freezing point was -2.5°C. Two formulas were considered, containing 6% or 12% of fat.

Experimental platform

The milk mix was aerated and frozen in a continuous SSHE (WCB® Model MF 50), a schematic representation of which is shown in Fig. 1.

Prior to use, the sorbet mix was stored in a refrigerated storage tank at 5 °C (1) . During the freezing process, the sorbet mix was supplied into the SSHE by means of a volumetric piston pump (2). Once the sorbet was frozen in the SSHE (3), it flowed down through the outlet pipe of the SSHE (4). A Pt100 probe measured the temperature of the product at the outlet. The freezer has a variable capacity from 25 to 75 kg/h. The inner diameter of the heat exchange tube of the freezer is 5cm and its length is 40cm (total volume of cylinder 0.785 cm 3 ). Inside the cylinder (figure 1), there is a dasher provided with two scraper blades mounted in a floating manner on the rotor (dasher volume = 0.351cm 3 ). The factory settings allowed the variation of the dasher speed from 200 to 1000 rpm. Ice cream samples with 75% overrun (volume of gas/ volume of mix) were prepared at two dasher speed: 200 and 1000 rpm. Samples were obtained from the two formula (fat content 6% or 12%). 500ml pins were filled with ice cream and placed in a freezer at -30°C and stored for the duration of experiments.

Meltdown Test

Prior to the melting test, ice cream was placed at -20°C for 24 hrs. Conical samples of following dimensions were used: bottom diameter d 1 =9cm, top diameter d 2 =8cm, height h=5cm. The samples were instrumented with a thermocouple T at the centre of the sample in order to follow temperature evolution during the test. Temperature was recorded by using a datalogger. Ice cream samples were placed in a temperature controlled room at 20°C where meltdown tests of ice cream have been performed by quantifying and determining the mass that drips from the product through a mesh screen as a function of the time. This way, meltdown curves (plots of melted product versus time) for the different samples have been generated. Pictures of the samples melting were obtained every 10 minutes to measure the deformation/collapse of the sample vs time. (Figure 2

MODEL DESCRIPTION

A numerical model was developed to describe the multiphysics phenomena occurring during an ice cream meltdown test. A number of assumptions were considered:

-Ice cream is considered as a porous structure constituted by fat, proteins and air bubbles filled with an aqueous sucrose solution which is more or less frozen. -Ice melting occurs over a temperature range from T f0 (onset melting temperature) to T f1 (onset freezing temperature) and concerns the freezable water mass fraction in the mix: ω f . The external temperature T a is constant and a uniform heat transfer coefficient h is applied. -Liquid flow through the sample (dripping) follows Richard's equation. The retention model of van Genuchten (1980) was considered for the effective saturation S e , the specific moisture capacity C m and the relative permeability k r in function of the pressure head H p =p/(ρ l g). The hydraulic permeability (κ) depends on the ice fraction which in turn depends on temperature and the liquid viscosity (µ l ) depends on the temperature and the sucrose concentration of the residual solution which depends on ice fraction. For κ/µ l , a sigmoid function varying from (κ/µ l ) 0 below T f0 to (κ/µ l ) 1 above T f1 was used. The initial internal pressure head is zero (saturated state). The external pressure head (lower and lateral side of the sample) varies from 0 below T f0 to H p1. above T f1 . -Deformation of the ice cream during melting is represented by a linear isotropic elasticity model. Young's modulus decreases with temperature from E 0 below T f0 to E 1 above T f1 . There is no stress on the upper and lateral surfaces. No displacement is assumed at the bottom. The small deformation approximation is used as a first approach even if rather large deformations are expected. Deformation is not taken into account in the heat and flow equations.

The set of equations to be solved is the following Energy equation ( )

not if C C ; T T T if T T ΔH ω C C with T k t T ρC p p.app f1 f0 f0 f1 f f p p.app p.app = < < - + = ∇ ⋅ ∇ = ∂ ∂ C C Flow equations (Richard's equation) T of function sigmoid a is μ κ where z) (H k μ κ t H g ρ C l p r l p l m         + ∇ ⋅ ∇ = ∂ ∂ C C [ ] ( ) ( ) [ ] 1 k 0; C 1; S ; θ θ : 0 H if S 1 1 S k ; S 1 )S θ (θ m 1 αm C ; αH 1 S ); θ (θ S θ θ : 0 H if r m e s p 2 m 1/m e p e r m 1/m e 1/m e r s m m m) 1/(1 p e r s e r p = = = = ≥ - - = - - - = + = - + = < - -
where θ is the volume of liquid per porous medium volume, θ r its residual value and θ s the total porosity.

Mechanical Equation

T of function sigmoid a is E and I ) X . ( 2ν 1 2 ) X ( X ν) 2(1 E σ where g ρ σ t X ρ t 2 2       ∇ - + ∇ + ∇ + = + ⋅ ∇ = ∂ ∂ C C C C C C C C C

Following table summarizes the parameter values

T f0 T f1 T a h ρ C p k ω f ρ l -10 -2.6 20 10 600 2000 0.6 0.6 1000 °C °C °C W⋅m -2 ⋅K -1 kg⋅m -3 J⋅kg -1 ⋅K -1 W⋅m -1 ⋅K -1 kg⋅m -3 (κ/µ l ) 0 (κ/µ l ) 1 m p α θ s θ r H p1 E 0 E 1 ν 10 -10 10 -9 1/2 1/2 1 0.5 0 -1.2 10 5 250 0.5 m 3 ⋅s⋅kg -1 m 3 ⋅s⋅kg -1 m m Pa Pa
These equations (completed by the adequate initial and boundary conditions) were solved with the finite element method implemented in the Comsol software.

RESULTS AND DISCUSSION

Figure 3 present the simulated temperature field at three different instants. It shows also the deformation of the sample. Globally, heat penetrates progressively toward the centre, the frozen part remains stiff whereas the partially melted zone is submitted to deformation and dripping.

Figure 3: Deformation and temperature field after 5min, one hour and two hours.

Figure 4 shown below, presents the temperature evolution at the centre of the sample. It appears, both for experimental and simulated evolutions, that temperature rises relatively quickly up to about -10°C, then it increases slowly because of the melting of freezable water. Above about -2°C, temperature again rises rapidly tending toward the ambient temperature. The experimental evolution is smoother than the predicted one. This could be explained by the fact that the sensor remained at the same height (2.5 cm) whereas the 'centre' of the sample goes down.

Fig. 4. Measured and simulated temperature evolutions at the centre of the ice cream sample during meltdown Figure 5 presents the drip loss. After 3 hours (180 min), the mass loss is about 55 g for both experimental and simulated values. But the simulated evolution does not agree with the measured one. The predicted drip loss occurs earlier and more suddenly. It has to be mentioned that the model parameters were not precisely adjusted in this preliminary work. The maximal permeability value should certainly be lowered and the temperature range of permeability increase seems to be not exactly the same as the range of ice melting.

CONCLUSION

A model coupling heat transfer, phase change, drip loss and deformation during the melting of a sample of ice cream was developed by using a number of simplification assumptions. With a first set of parameters, it was able to predict the core temperature evolution, the drip loss and the deformation after 3 hours of melting. Further parameter adjustment is necessary to better fit the drip loss and deformation evolutions. The model could also be improved by using a large deformation ALE (Arbitrary Lagrangian Eulerian) approach.

Fig. 1 .

 1 Fig. 1. Schematic representation of the experimental platform. 1. Refrigerated storage tank (200 litres capacity). 2. Volumetric piston pump. 3. Freezer. 4. Outlet pipe. 5. Pt100 probe inserted into the outlet pipe.

Figure 2 :

 2 Figure 2 : Pictures of a sample at t=5min, after 60 min and 120 min of melting.

Fig. 5 .Fig. 6

 56 Fig. 5. Measured and simulated drip loss

NOMENCLATURE


Karim Essafi
email: karim.essafi@oist.jp

L D C Jaubert

M Udagawa

Flat bands and Dirac cones in breathing lattices

In breathing pyrochlores and kagomes, couplings between neighbouring tetrahedra and triangles are free to differ. Breathing lattices thus offer the possibility to explore a different facet of the rich physics of these systems. Here we consider nearest-neighbour classical Heisenberg interactions, both ferromagnetic and antiferromagnetic, and study how the anisotropy of breathing lattices modifies the mode spectrum of pyrochlore and kagome systems. The nature and degeneracy of the flat bands are shown to be preserved for any value of the anisotropy. These flat bands can coexist with Dirac nodes at the Γ point when the model becomes particle-hole symmetric. We also derive the nature of the ground state for the breathing kagome lattice, which bears a spontaneous chirality when neighbouring triangles are alternatively ferromagnetic and antiferromagnetic.

Introduction

Pyrochlore and kagome antiferromagnets are archetypical examples of flat-band systems in frustrated magnetism [1,2]. In these classical spin liquids, the flat bands correspond to highly degenerate ground-state manifolds which support emergent Coulomb physics [3]. However, while the pyrochlore and kagome antiferromagnets have been rather extensively studied, little is known about their "breathing" counterparts. Pyrochlore and kagome lattices are made of corner-sharing units, respectively tetrahedra and triangles. The term "breathing" has been coined to describe lattices where spin couplings of neighbouring units are different [4]. Such geometry has first been used as a mathematical convenience to study the quantum spin liquid nature of the pristine spin-1/2 pyrochlore and kagome antiferromagnets [5][6][7][8][9][10][11][12][13]. But recently, breathing lattices have gained noticeable interest on their own as possible outcomes of spin-lattice coupling [14,15], and after being synthesized in the rare-earth based pyrochlore material Ba 3 Yb 2 Zn 5 O 11 [16][17][18][19][20], in the vanadium oxyfluoride kagome [21][22][23] and in the spinel oxides LiGaCr 4 O 8 and LiInCr 4 O 8 [4,[24][25][26][27], where the large magnetic moment of the Cr 3+ ions makes a classical approach sensible [15,28].

[NH 4 ] 2 [C 7 H 14 N][V 7 O 6 F 18 ] [
In this paper, our goal is to complement the understanding of classical breathing lattices by studying the evolution of the flat bands originally present in the uniform lattices. The anisotropy of breathing lattices will be used as a tunable parameter, exploring both antiferromagnetic and ferromagnetic couplings between classical Heisenberg spins. We will explain why the flat bands, and their degeneracy, are protected, and how Dirac cones appear at the Γ point when neighbouring units have couplings of exactly opposite sign. In addition to their band structures, the ground state of the breathing kagome lattice will be described, both for the Heisenberg and XXZ Hamiltonians.

Breathing lattices 2.1. Geometry

Pyrochlores and kagomes are frustrated lattices made respectively of corner-sharing tetrahedra and triangles [Fig. 1]. For convenience, tetrahedra and triangles shall be referred to as units. The centres of these units form bipartite lattices (respectively diamond and honeycomb). It is thus possible to define two types of units labeled by A and B [coloured in purple and red in Fig. 1].

Let R a be the position of the centre of unit a. Within a A-type tetrahedron, the positions of the 4 spin sublattices with respect to R a are

δ δ δ 1 A = 1 8   1 1 1   , δ δ δ 2 A = 1 8   -1 -1 1   , δ δ δ 3 A = 1 8   -1 1 -1   , δ δ δ 4 A = 1 8   1 -1 -1   . (1)
For a B-type tetrahedron, δ δ δ i B = -δ δ δ i A for i ∈ {1, 2, 3, 4}. All distances are given in units of the linear size of the cubic unit cell made of 16 pyrochlore sites. Within a A-type triangle, there are 3 spin sublattices whose positions with respect to R a are

δ δ δ 1 A = 1 4 -1 -1/ √ 3 , δ δ δ 2 A = 1 4 1 -1/ √ 3 , δ δ δ 3 A = 1 4 0 2/ √ 3 . (2)
For a B-type triangle, δ δ δ i B = -δ δ δ i A for i ∈ {1, 2, 3}. All distances are given in units of the linear size of the triangular unit cell made of 3 sites.

Hamiltonian

Let us consider a Hamiltonian with isotropic nearest-neighbour exchange

H = J A a∈A ij S i a • S j a + J B b∈B ij S i b • S j b , (3)
where S i a is a classical Heisenberg spin of unit length |S i a | = 1 on sublattice i belonging to the unit a. The first and second terms of Eq. ( 3) run over all A-and B-type of units respectively. It will be useful to introduce the breathing factor α = J B /J A [4].

3. Mode spectrum of the interaction matrix 3.1. Interaction matrix In pyrochlore and kagome lattices, every spin belongs to two different units, indexed a o ∈ A and b o ∈ B. This means that every spin can be written as S i ao or as S i bo with position

R i ao + δ δ δ i A = R i bo + δ δ δ i B = R i bo -δ δ δ i A .
Keeping this in mind, one can define the Fourier transform of the spin degrees of freedom for a given sublattice i.

Ŝi q = 1 √ N ′ a∈A S i a e -ıq.(Ra+δ δ δ i A ) = 1 √ N ′ b∈B S i b e -ıq.(R b +δ δ δ i B ) , (4)
S i a = 1 √ N ′ q Ŝi q e +ıq.(Ra+δ δ δ i A ) = 1 √ N ′ q Ŝi q e +ıq.(R b +δ δ δ i B ) = S i b , (5)
where N ′ is the total number of sites in a given spin sublattice. If N is the total number of spins in the system and s is the number of spin sublattices, then N ′ = N/s. The Fourier transform of Eq. ( 3) takes the form [28]

H = J A 2 q s i=1 s j=1 j =i Ŝi q • Ŝj -q e ıq.(δ δ δ i A -δ δ δ j A ) + J B 2 q s i=1 s j=1 j =i Ŝi q • Ŝj -q e ıq.(δ δ δ i B -δ δ δ j B ) (6) = J A 2 q s i=1 s j=1 j =i Ŝi q • Ŝj -q e ıq.(δ δ δ i A -δ δ δ j A ) + α e -ıq.(δ δ δ i A -δ δ δ j A ) . (7)
The dispersion relations are obtained by diagonalizing the q-dependent interaction matrix K q whose coefficients are K ii q = 0 and, for i = j,

K ij q = e ıq.(δ δ δ i A -δ δ δ j A ) + α e -ıq.(δ δ δ i A -δ δ δ j A ) = (1 + α) cos q.(δ δ δ i A -δ δ δ j A ) + ı(1 -α) sin q.(δ δ δ i A -δ δ δ j A ) . (8)
As opposed to the uniform pyrochlore [1,29] and kagome [2] lattices (α = 1), the coefficients of K q are intrinsically complex for breathing lattices [28]. However since K q is hermitian, the eigenvalues remain of course real. For most of the paper, the discussion will be focused on J A > 0 and -1 ≤ α ≤ 1, because the mode spectrum for J A < 0 and |α| > 1 can be easily derived by a rescaling of the coupling constants.

Breathing pyrochlore

The four eigenvalues of K q are given by

           λ 1 q = α + 1 -2 α 2 -α + α cos qx 2 (cos qy 2 + cos qz 2 ) + α cos qy 2 cos qz 2 + 1 λ 2 q = α + 1 + 2 α 2 -α + α cos qx 2 (cos qy 2 + cos qz 2 ) + α cos qy 2 cos qz 2 + 1 λ 3 q = -α -1 λ 4 q = -α -1 (9)
with two degenerate flat bands, as explained in Sec. 3.2.

Breathing kagome

The three eigenvalues of K q are given by

           λ 1 q = 1 2 1 + α -9α 2 -6α + 8α cos(q x ) + 16α cos( qx 2 ) cos √ 3 2 q y + 9 λ 2 q = 1 2 1 + α + 9α 2 -6α + 8α cos(q x ) + 16α cos( qx 2 ) cos √ 3 2 q y + 9 λ 3 q = -α -1 (10)
with one flat band, as explained in Sec. 3.2.

Origin of the flat bands

In breathing pyrochlore and kagome lattices, some of the bands are completely flat, as shown in ( 9) and (10). The origin of these flat bands can be clarified by the factorization of the Hamiltonian (3), as was done for the tight-binding model on complete graphs [30]. To this aim, we rewrite the Hamiltonian (3) as

H = 1 2 m,n Ĥmn S m • S n . ( 11
)
Here, S m is the spin defined at the site, m, and the summation with m and n is made over the entire lattice. The coefficient matrix, Ĥ, can be factorized by introducing an auxiliary hybrid lattice, as we show in Fig. 2 for the breathing kagome lattice. Firstly, we put sites at the centers of units, and secondly, we connect them to the original sites [Fig. The new auxiliary lattice consists of two classes of sites: those of the original lattice, and those at the centres of the units of the original lattice. We call the former class ℓ, and the latter, c. The number of sites in class ℓ and c are defined as N ℓ and N c , respectively. These are given as N ℓ = N and N c = 2N ′ = 2N/s for both lattices.

On this auxiliary lattice, we introduce two rectangular matrices, Ĥc←ℓ and Ĥℓ←c = ( Ĥc←ℓ ) T . Ĥc←ℓ is a N c × N ℓ incidence matrix, and its column (row) indices are related to the sites in ℓ (c). We set its (m, n) component, ( Ĥc←ℓ ) mn , to be 1, if the corresponding sites, m ∈ c and n ∈ ℓ are connected, and to be 0, otherwise. Being a N c × N ℓ rectangular matrix, the rank of Ĥc←ℓ is bounded by min(N c , N ℓ ) = N c1 . Given these definitions, the matrix Ĥ of Eq. ( 11) can be factorized as

Ĥ = J A Ĥℓ←c τ Ĥc←ℓ -J A (1 + α) 1 ℓ , (12)
where 1 ℓ is the N ℓ × N ℓ identity matrix and τ is a N c × N c diagonal matrix, τmn = τ m δ mn , with τ m = 1(α) if m belongs to sublattice A (B). This expression immediately leads to the protection of the flat bands present in the uniform lattices. Firstly, suppose there is a state annihilated by Ĥc←ℓ , then this state serves as an eigenstate of Ĥ with eigenvalue, -J A (1 + α), and is independent of α. Secondly, there are in fact at least (N ℓ -N c ) = (s -2)N ′ such states, since Ĥc←ℓ acts on a space of dimension N ℓ , while the rank of Ĥc←ℓ is bounded by N c . Consequently, Ĥ has at least (s -2)N ′ degenerate states with eigenvalue, -J A (1 + α). This argument proves that the s -2 = 2 (1) flat bands of the uniform pyrochlore (kagome) lattices see their energy shifted by the breathing anisotropy, but are not destroyed and conserve their degeneracy. In addition, the presence of new flat bands is allowed, which is precisely what happens for α = 0 (see below).

Gap opening when |J

A | = |J B |
The mode spectra for regular pyrochlore and kagome lattices are known to be gapless [1,2].

Here we find that for α > 0 excitations out of the flat bands remain gapless, but the breathing anisotropy opens a gap between the dispersive bands λ 1 q and λ 2 q when |α| = 1 [Figs. 3 &4]. From Eqs. ( 9) and ( 10), one gets

∆ pyrochlore = 4 |1 -|α|| , ∆ kagome = 3 |1 -|α|| . ( 13
)
Starting from α = 1 and decreasing its value, λ 1 q and λ 2 q become less and less dispersive [Figs. 3 &4] and flat at α = 0. For decoupled antiferromagnetic units, the ground state is made of three (resp. two) degenerate flat bands in pyrochlores (resp. kagomes). Figure 6. Mode spectrum of K q for breathing kagomes and α ∈ {-0.5, -1, -2}. The left column shows the entire (q x , q y ) plane, while the right column shows a cut along the q y = 0 line. For α < 0, the flat bands are not part of the ground state anymore. The gap closes at α = -1 where a Dirac cone appears at the Γ point. Wavevectors are in 2π units.

Particle-hole symmetry and Dirac cones for J

A = -J B
The breathing factor α can be seen as a way to tune the relative importance of the real and imaginary parts in K q [see Eq. ( 8)]. When α = -1, K q becomes purely imaginary, K q = ıΩ q with Ω q ∈ M(s, R), the ensemble of real matrices of size s × s. Since K q is hermitian, its eigenvalues are real. Let V q be an eigenvector of K q with eigenvalue λ q , K q V q = ıΩ q V q = λ q V q .

(14

)
The complex conjugate of Eq. ( 14) gives

ıΩ q V q = λ q V q ⇔ -ıΩ q V q = λ q V q ⇔ K q V q = -λ q V q . ( 15
)
It means that if the eigenvector V q is real, then λ q = 0. But if V q is complex, then -λ q is also an eigenvalue of the interaction matrix K q . Hence, breathing lattices possess particle-hole symmetry at α = -1 [Figs. 5 and6]. Furthermore K q becomes nil at the Γ point. It means that all the bands connect at q = 0 and that the flat bands are at zero energy.

An interesting outcome is that the dispersive bands form Dirac cones at the Γ point [Figs. 5 and6]. In linear order in |q|, the dispersion relation is λ 2 q = -λ 1 q = |q| for pyrochlores and λ 2 q = -λ 1 q = 3/2|q| for kagomes. Please note that the presence of a single Dirac cone per Brillouin zone is allowed despite the Nielsen-Ninomiya theorem because of their co-existence with flat bands [32].

4. Ground states for α < 0 4.1. Breathing pyrochlore In agreement with the work of Benton & Shannon [28], we find that the classical spin liquid of the pyrochlore antiferromagnet (α = +1) is robust for positive α but disappears as soon as α becomes negative [Figs. 3 and5]. For α < 0 the ground-state configurations belong to one of the dispersive bands, with lines of zero-energy modes along the [1k0] directions (and equivalent wavevectors) [Fig. 7]. An order-by-disorder transition is expected to select spin configurations with [100] ordering [28].

Breathing kagome

While the evolution of the band structure shares many similarities between breathing pyrochlore and kagome, their low temperature physics differ qualitatively. Indeed, there are no lines of zeroenergy fluctuations in kagome for α < 0 [Fig. 6]. Au contraire, the band structure shows clear minima at the corners of the Brillouin zone, which can be understood as follows.

In order to minimize the energy of an antiferromagnetic triangle (type A) with classical Heisenberg spins, the three spins must be coplanar with a 2π/3 angle between them. For a ferromagnetic triangle (type B), spins must be collinear. Any configuration which respects these constraints over the entire lattice belongs to the ground-state manifold, irrespectively of the value of α as long as α < 0. Such a configuration is depicted on Fig. 8 and is uniquely determined once the spin configuration of one A-type triangle is given. To prove this uniqueness, let A 1 in Fig. 8 be this initial triangle. Without any loss of generality, its spin configuration is chosen as

S 1 =   0 1 0   , S 2 = 1 2   - √ 3 -1 0   , S 3 = 1 2   √ 3 -1 0   , (16)
whose respective color codes are red, cyan and green in Fig. 8. Being ferromagnetic, the neighbouring triangles B 2 and B 3 are immediately fixed into configurations S 2 and S 3 respectively. On the other hand, the orientations of the four spins S i∈{̟,µ,ν,η} are a priori not fixed. The local antiferromagnetic constraints on triangles A i∈{2,3} indeed allow for a O(2) degree of freedom for each triangle. These continuous degrees of freedom are however constrained by the intervening ferromagnetic triangle B 1 which imposes S µ = S ν . From Eq. ( 16), this constraint can be rewritten as

S µ = S ν ⇔        √ 3 4 (1 + sin θ µ ) 1 4 (1 -3 sin θ µ ) √ 3 2 cos θ µ        =        √ 3 4 (-1 + sin θ ν ) 1 4 (1 + 3 sin θ ν ) √ 3 2 cos θ ν        ⇔ θ ν = -θ µ = π 2 ⇔ S µ = S ν = S 1 . ( 17
)
Step by step, this argument around a plaquette of six triangles can be applied to the entire lattice. The ground-state of the breathing kagome model for α < 0 is thus long-range ordered with a global degeneracy of O(3)×O(2). In particular, there are no weathervane modes [33,34]. If one assigns its corresponding color -red, cyan, green -to each of the ferromagnetic (Btype) triangles, then any ground state configuration can be seen as a three-coloring problem on the triangular lattice. The magnetization of the system is thus zero. On the other hand, the vector chirality defined as

χ χ χ ≡ 1 N c ∈{A,B} S 1 c × S 2 c + S 2 c × S 3 c + S 3 c × S 1 c (18)
is finite since all antiferromagnetic (A-type) triangles carry the same chirality while the ferromagnetic (B-type) ones bear none. Within the ground-state manifold, the vector chirality χ χ χ is uniformly distributed on a sphere of radius √ 3/2. Please note that the scalar chirality is nil since the configurations are coplanar.

Interestingly, if the interactions between spins are made anisotropic by favouring in-plane couplings over the out-of-plane ones (XXZ model), then the configuration of Fig. 8 remains one of the ground states. The only difference is that spins now lie in the kagome plane, and the ground-state manifold has Z 2 ×O(2) symmetry. The O(2) symmetry has the same origin as before (in-plane rotation of all spins) while the Z 2 symmetry corresponds to the sign of the vector chirality χ χ χ = ± √ 3/2 e z where e z is the out-of-plane unit vector.

Conclusion

We have calculated the mode spectrum of the interaction matrix K q for the breathing pyrochlore and breathing kagome lattices. While the nature and degeneracy of the flat bands are not affected by the breathing anisotropy, the shape of the dispersive bands dramatically changes as a function of α. This is best illustrated for α = -1 where a Dirac cone is formed at the Γ point and the system gains particle-hole symmetry. The extensive ground-state degeneracy of the antiferromagnetic models is lifted when α < 0, in favour of lines of zero modes for breathing pyrochlore [28] and a long-range ordered ground state with finite vector chirality for breathing kagome.

A promising outcome of our work is the observation of a Dirac cone in a particle-hole symmetric model for α = -1. The coexistence of Dirac nodes and flat bands has also been observed for example in a tight-binding model on the Lieb lattice with three species of spinless fermions [35] and on the Shastry-Sutherland lattice with Dzyaloshinskii-Moriya interactions in a field [36]. In both cases, it was possible to open and close the gap at the Dirac point between the flat band and both dispersive bands by tuning either the tunneling coefficients [35] or the magnetic field [36]. It would be interesting to see if perturbations of our model could open such a gap, which would establish breathing lattices as a favourable setup for topological phenomena [14].

Figure 1 .

 1 Figure 1. Breathing pyrochlore (left) and breathing kagome (right). The purple and red tetrahedra/triangles represent the A-and B-type of units with respectively J A and J B coupling constants. There are 4 spin sublattices in pryochlores and 3 spin sublattices in kagomes.

  2 (b)], and finally we remove the original bonds [Fig. 2 (c)]. The resultant lattice is the decorated diamond (honeycomb) lattice for breathing pyrochlore (kagome) lattice.

Figure 2 .

 2 Figure 2. The process of making an auxiliary hybrid lattice from a breathing kagome lattice. (a) A breathing kagome lattice. (b) Sites are added in the center of triangle units. (c) The bonds of original breathing kagome lattice is removed.

Figure 3 .

 3 Figure 3. Mode spectrum of the interaction matrix for breathing pyrochlores in the [hhk] plane at α = 1 (see Ref.[1]) and α = 0.5. The ground state is a doubly degenerate flat band with gapless excitations. A gap between the dispersive bands opens for |α| = 1. Wavevectors are in 2π units.

Figure 4 .

 4 Figure 4. Mode spectrum of the interaction matrix for breathing kagomes at α = 1 (see Ref.[2]) and α = 0.5. The ground state is a unique flat band with gapless excitations. A gap between the dispersive bands opens for |α| = 1. Wavevectors are in 2π units.

Figure 5 .

 5 Figure 5. Mode spectrum of K q for breathing pyrochlores and α ∈ {-0.5, -1, -2}. The left column shows the entire [hhk] plane, while the right column shows a cut along the [hhh] line. For α < 0, the flat bands are not part of the ground state anymore. The gap closes at α = -1 where a Dirac cone appears at the Γ point. Wavevectors are in 2π units.

Figure 7 .

 7 Figure 7. Mode spectrum of K q for breathing pyrochlores in the [hk0] plane at α = -0.5 and α = -1. The [1k0] directions (and equivalent wavevectors) form lines of zero-energy modes. Wavevectors are in 2π units.

Figure 8 .

 8 Figure 8. Example of a ground-state configuration of the breathing kagome model for α < 0. Each colored dot corresponds to a spin orientation given on the left.

It is actually possible to prove that rank Ĥc←ℓ = Nc -1 by counting the number of localized modes[31].

Acknowledgments This work is supported by KAKENHI (Nos. 26400339, 15H05852, 15K13533 and 16H04026). K. E. and L. D. C. J. are supported in part by the Okinawa Institute of Science and Technology Graduate University. M. U. is supported by MEXT-Supported Program for the Strategic Research Foundation at Private Universities, 2015-2019


Lionel Rosier
email: lionel.rosier@mines-paristech.fr.

On the Benjamin-Bona-Mahony equation with a localized damping

Keywords: AMS Subject Classification: 35Q53, 93D15 Benjamin-Bona-Mahony equation, unique continuation property, internal stabilization, boundary stabilization

We introduce several mechanisms to dissipate the energy in the Benjamin-Bona-Mahony (BBM) equation. We consider either a distributed (localized) feedback law, or a boundary feedback law. In each case, we prove the global wellposedness of the system and the convergence towards a solution of the BBM equation which is null on a band. If the Unique Continuation Property holds for the BBM equation, this implies that the origin is asymptotically stable for the damped BBM equation.

Introduction

The Benjamin-Bona-Mahony equation

v t -v xxt + v x + vv x = 0, (1.1)
was proposed in [START_REF] Benjamin | Model equations for long waves in nonlinear dispersive systems[END_REF] as an alternative to the Korteweg-de Vries (KdV) equation as a model for the propagation of one-dimensional, unidirectional, small amplitude long waves in nonlinear dispersive media. In the context of shallow water waves, v = v(x, t) stands for the displacement of the water surface (from rest position) at location x and time t. In the paper, we shall assume that either x ∈ R, or x ∈ (0, L) or x ∈ T = R/(2πZ) (the one-dimensional torus). The dispersive term -v xxt produces a strong smoothing effect for the time regularity, thanks to which the wellposedness theory of (1.1) is easier than for KdV (see [START_REF] Bona | Sharp well-posedness for the BBM equation[END_REF][START_REF] Roumégoux | A symplectic non-squeezing theorem for BBM equation[END_REF]). Solutions of (1.1) turn out to be analytic in time. On the other hand, the control theory is at his early stage for BBM (for the control properties of KdV, we refer the reader to the recent survey [START_REF] Rosier | Control and stabilization of the Korteweg-de Vries equation: recent progresses[END_REF]). S. Micu investigated in [START_REF] Micu | On the controllability of the linearized Benjamin-Bona-Mahony equation[END_REF] the boundary controllability of the linearized BBM equation, and noticed that the exact controllability fails due to the existence of a limit point in the spectrum of the adjoint equation. The author and B.-Y. Zhang introduced in [START_REF] Rosier | Unique continuation property and control for the Benjamin-Bona-Mahony equation[END_REF] a moving control and derived with such a control both the exact controllability and the exponential stability of the full BBM equation. For a distributed control with a fixed support, the exact controllability of the linearized BBM equation fails, so that the study of the controllability of the full BBM equation seems hard. However, it is reasonable to expect that some stability properties could be derived by incorporating some dissipation in a fixed subdomain or at the boundary. The aim of this paper is to propose several dissipation mechanisms leading to systems for which one has both the global existence of solutions and 1 a nonincreasing H 1 -norm. The conclusion that all the trajectories emanating from data close to the origin are indeed attracted by the origin is valid provided that the following conjecture is true:

Unique Continuation Property (UCP) Conjecture: There exists some number δ > 0 such that for any

v 0 ∈ H 1 (T) with v 0 H 1 (T) < δ, if the solution v = v(x, t) of v t -v xxt + v x + vv x = 0, x ∈ T, v(x, 0) = v 0 (x), x ∈ T (1.2) satisfies v(x, t) = 0 ∀(x, t) ∈ ω × (0, T ) (1.3)
for some nonempty open set ω ⊂ T and some time T > 0, then v 0 = 0 (and hence v ≡ 0).

To the best knowledge of the author, the UCP for BBM as stated in the above conjecture is still open. The main difficulty comes from the fact that the lines x = 0 are characteristic for BBM, so that the "information" does not propagate well in the x-direction. For some UCP for BBM (with additional assumptions) see [START_REF] Rosier | Unique continuation property and control for the Benjamin-Bona-Mahony equation[END_REF][START_REF] Yamamoto | One unique continuation for a linearized Benjamin-Bona-Mahony equation[END_REF]. See also [START_REF] Micu | Control and stabilization of a family of Boussinesq systems[END_REF][START_REF] Micu | Stabilization of a Boussinesq system with localized damping[END_REF] for control results for some Boussinesq systems of BBM-BBM type.

The following result is a conditional UCP in which it is assumed that the initial data is small in the L ∞ -norm and it has a nonnegative mean value. Its proof was based on the analyticity in time of the trajectories and on the use of some Lyapunov function.

Theorem 1 [START_REF] Rosier | Unique continuation property and control for the Benjamin-Bona-Mahony equation[END_REF] Let u 0 ∈ H 1 (T) be such that

T v 0 (x) dx ≥ 0, (1.4
)

and v 0 L ∞ (T) < 3. (1.5)
Assume that the solution v to (1.2) satisfies (1.3). Then v 0 = 0.

As it was noticed in [START_REF] Rosier | Unique continuation property and control for the Benjamin-Bona-Mahony equation[END_REF], the UCP for BBM cannot hold for any state in L 2 (T), for any initial data v 0 with values in {-2, 0} gives a trivial (stationary) solution of BBM. Thus, either a regularity assumption (v 0 ∈ H 1 (T)), or a bound on the norm of the initial data has to be imposed for the UCP to hold.

The paper is outlined as follows. In Section 2, we incorporate a simple localized damping in BBM equation and investigate the corresponding Cauchy problem. In Section 3, we consider another dissipation mechanism involving one derivative. The last section is concerned with the introduction of nonhomogeneous boundary conditions leading again to the dissipation of the H 1 -norm.

Stabilization of the BBM equation 2.1 Internal stabilization with a simple feedback law

We consider the BBM equation on T with a localized damping

u t -u xxt + u x + uu x + a(x)u = 0 x ∈ T, t ≥ 0, (2.1) u(x, 0) = u 0 (x) x ∈ T, (2.2)
where a is a smooth, nonnegative function on T with a(x) > 0 on a given open set ω ⊂ T. We write (2.1)-(2.2) in its integral form

u(t) = u 0 - t 0 (1 -∂ 2 x ) -1 [a(x)u + (u + u 2 2 ) x ](τ )dτ (2.3) where (1 -∂ 2 x ) -1 f denotes for f ∈ L 2 (T) the unique solution v ∈ H 2 (T) of (1 -∂ 2 x )v = f . Define the norm • s in H s (T) as n∈Z c n e inx 2 s = n∈Z |c n | 2 (1 + |n| 2 ) s .
We have the following result.

Theorem 2 Let s ≥ 0 be given. For any u 0 ∈ H s (T), there exist T > 0 and a unique solution u of

(2.3) in C([0, T ], H s (T)). Moreover, the correspondence u 0 ∈ H s (T) → u ∈ C([0, T ], H s (T)) is Lipschitz continuous. If s = 1
, the solution exists for every T > 0, and the energy u(t) 2 1 is nonincreasing. Finally, if v 0 denotes any weak limit in H 1 (T) of a sequence u(t n ) with t n ր +∞, then the solution v of system (1.2) satisfies (1.3). In particular, if the UCP conjecture holds, then v 0 = 0, so that u(t) → 0 weakly in H 1 (T) as t → +∞, hence strongly in H s (T) for all s < 1.

Proof. We proceed as in [START_REF] Bona | Boussinesq equations and other systems for smallamplitude long waves in nonlinear dispersive media[END_REF] using the estimate

||f g|| s ≤ C||f || s+1 ||g|| s+1 ∀s ≥ -1, ∀f, g ∈ H s+1 (T). (2.4)
The estimate (2.4) follows from a similar estimate proved in [START_REF] Bona | Boussinesq equations and other systems for smallamplitude long waves in nonlinear dispersive media[END_REF] for functions defined on R, namely

|| f g|| H s (R) ≤ C|| f || H s+1 (R) ||g|| H s+1 (R) ∀s ≥ -1, ∀f, g ∈ H s+1 (R), (2.5) by letting f (x) = ϕ(x)f (x), g(x) = ϕ(x)g(x)
, where f and g are viewed as 2π-periodic functions, and

ϕ ∈ C ∞ 0 (R) denotes a cut-off function such that ϕ(x) = 1 on [0, 2π]. Indeed, we notice that ||f || H s (T) ≤ || f || H s (R) ≤ C||f || H s (T)
for some constant C > 0. Note that for any s ≥ 0

||(1 -∂ 2 x ) -1 ∂ x (f g)|| s ≤ C||f g|| s-1 ≤ C||f || s ||g|| s .
Pick any u 0 ∈ H s (T). Let us introduce the operator

(Γu)(t) := u 0 - t 0 (1 -∂ 2 x ) -1 [a(x)u + (u + u 2 2 ) x ](τ )dτ. Then sup 0≤t≤T ||(Γu -Γv)(t)|| s ≤ C 2 T 0 [||a(u -v)|| s-2 + ||u -v|| s-1 + ||u -v|| s ||u + v|| s ]dτ ≤ C 3 T (1 + 2R)||u -v|| C([0,T ],H s (T))
if we assume that u and v are in the closed ball B R of radius R centered at 0 in C([0, T ], H s (T)). We pick

T := [2C 3 (1 + 2R)] -1 so that ||Γu -Γv|| C([0,T ],H s (T)) ≤ 1 2 ||u -v|| C([0,T ],H s (T)) •
On the other hand

||Γu|| C([0,T ],H s (T)) ≤ ||u 0 || s + ||Γu -Γ0|| C([0,T ],H s (T)) ≤ ||u 0 || s + R 2 ≤ R
for the choice R = 2||u 0 || s . It follows that the map Γ contracts in B R , hence it admits a unique fixed point u in B R which solves the integral equation (2.3). Furthermore, given any ρ > 0 and any u 0 , v 0 ∈ H s (T) with ||u 0 || s ≤ ρ, ||v 0 || s ≤ ρ, one easily sees that for

T = [2C 3 (1 + 4ρ)] -1 (2.6) one has ||u -v|| C([0,T ],H s (T)) ≤ 2||u 0 -v 0 || s . (2.7)
Finally assume that s = 1. Scaling in (2.1) by u yields

1 2 ||u(T )|| 2 1 - 1 2 ||u 0 || 2 1 + T 0 T a(x)|u(x, t)| 2 dxdt = 0 (2.8) (Note that u t = -(1 -∂ 2 x ) -1 [a(x)u + (u + u 2 2 ) x ] ∈ C([0, T ], H 2 (T)
) so that each term in (2.1) belongs to C([0, T ], L 2 (T)), and the integrations by parts are valid.) It follows that the map t → ||u(t)|| 1 is nonincreasing, hence it admits a nonnegative limit l as t → ∞, and that the solution of (2.1)-(2.2) is defined for all t ≥ 0. Let T be as in (2.6) with ρ = ||u 0 || 1 . Note that ||u(t)|| 1 ≤ ||u 0 || 1 for all t ≥ 0. Let v 0 be any weak limit of {u(t)} t≥0 in H 1 (T), i.e. we have that u(t n ) → v 0 weakly in H 1 (T) for some sequence t n → +∞. Extracting a subsequence if needed, we may assume that

t n+1 -t n ≥ T for n ≥ 0. (2.9) From 1 2 ||u(t n+1 )|| 2 1 - 1 2 ||u(t n )|| 2 1 + tn+1 tn T a(x)|u(x, t)| 2 dxdt = 0, we obtain that lim n→+∞ tn+1 tn T a(x)|u(x, t)| 2 dxdt = 0. (2.10) Since u(t n ) → v 0 in H s (T) for s < 1, and ||u(t n )|| s ≤ ||u(t n )|| 1 ≤ ρ, we have from (2.7) that for all s ∈ [0, 1[ u(t n + •) → v in C([0, T ], H s (T)) as n → +∞, (2.11)
where v = v(x, t) denotes the solution of

v t -v xxt + v x + vv x + a(x)v = 0, x ∈ T, t ≥ 0, v(x, 0) = v 0 (x),
x ∈ T.

Notice that v ∈ C([0, T ]; H 1 (T)), for v 0 ∈ H 1 (T). (2.10) combined to (2.9) and (2.11) yields

T 0 T a(x)|v| 2 dxdt = 0. (2.12) Therefore v ∈ C([0, T ]; H 1 (T)) solves v t -v xxt + v x + vv x = 0, x ∈ T, t ∈ (0, T ), v(x, 0) = v 0 (x) x ∈ T, v(x, t) = 0 x ∈ ω, t ∈ (0, T ).
If the UCP conjecture is true, we have that v ≡ 0 on R × (0, T ). It follows that v 0 ≡ 0, and that as t → +∞ u(t) → 0 weakly in H 1 (T), u(t) → 0 strongly in H s (T) for any s < 1.

Internal stabilization with one derivative in the feedback law

We now pay attention to the stabilization of the BBM equation by means of a "stronger" feedback law involving one derivative. More precisely, we consider the system

u t -u xxt + u x + uu x -(a(x)u x ) x = 0, x ∈ T, t ≥ 0, (2.13) u(x, 0) = u 0 (x), x ∈ T. (2.14)
where a = a(x) denotes again any nonnegative smooth function on T with a(x) > 0 on a given open set ω ⊂ T. Scaling in (2.13) yields (at least formally)

1 2 ||u(T )|| 2 1 - 1 2 ||u 0 || 2 1 + T 0 T a(x)|u x (x, t)| 2 dxdt = 0. (2.15)
The decay of the energy is quantified by an integral term involving the square of a localized H 1 -norm in (2.15). By contrast, the integral term in (2.8) involved the square of a localized L 2 -norm. This suggests that the damping mechanism involved in (2.13) acts in a much stronger way than in (2.1). As a matter of fact, in the trivial situation when a(x) ≥ C > 0 for all x ∈ T in (2.13), it is a simple exercise to establish the exponential stability in H 1 (T) for both the linearized equation and the nonlinear BBM equation for states with zero means.

In the general case when the function a(x) is supported in a subdomain of T, we obtain the following result.

Theorem 3 Let s ≥ 0 be given. For any u 0 ∈ H s (T), there exist T > 0 and a unique solution u of (2.13)-

(2.14) in C([0, T ], H s (T)). Moreover, the correspondence u 0 ∈ H s (T) → u ∈ C([0, T ], H s (T)) is Lipschitz continuous. If s = 1, one can pick any T > 0 and u(t) -[u 0 ] H 1 (T) is nondecreasing. If the UCP conjecture is valid by replacing (1.3) by v = C on ω × (0, T ), then u(t) → [u 0 ] = (2π) -1
T u 0 (x)dx weakly in H 1 (T ) hence strongly in H s (T) for s < 1 as t → ∞.

Proof. As the proof is very similar to those of Theorem 2, we limit ourselves to pointing out the only differences. For the wellposedness, we use the estimate valid for s ≥ 0

||(1 -∂ 2 x ) -1 (au x ) x || s ≤ C||u|| s .
For s = 1, we claim that (2.15) is justified by noticing that for u ∈ H 1 (T)

-(au x ) x , u H -1 ,H 1 = au x , u x L 2 ,L 2 .
Then the wellposedness statement is established as in Theorem 2. Let us proceed to the asymptotic behavior. Pick any u 0 ∈ H 1 (T) and any v 0 ∈ H 1 (T) which is the weak limit in H 1 (T) of a sequence u(t n ) with t n → ∞ and t n+1 -t n ≥ T . Let us still denote by v the solution of (1.2). Equation (2.12) has to be replaced by

T 0 T a(x)|v x (x, t)| 2 dxdt = 0. (2.16)
To prove (2.16), we notice first that u(t n + •) → v in C([0, T ], H s (T)) for s < 1, and that ||u(t n + •)|| L 2 (0,T ;H 1 (T)) ≤ const, so that, extracting a subsequence if needed, we have that

u(t n + •) → v weakly in L 2 (0, T ; H 1 (T)).
This yields, with (2.15),

T 0 T a(x)|v x | 2 dxdt ≤ lim inf n→∞ tn+T tn T a(x)|u nx | 2 dxdt = 0. Therefore v solves v t -v xxt + v x + vv x = 0, x ∈ T, t ∈ (0, T ), (2.17)
v x = 0, x ∈ ω, t ∈ (0, T ), (2.18) v ∈ C([0, T ]; H s (T)) for s < 1. (2.19)
(2.17) and (2.18

) yield v t = v x = 0 in ω × (0, T ), hence v(x, t) = C for (x, t) ∈ ω × (0, T ) for some constant C ∈ R. If the UCP Conjecture is still valid when v is constant on the band ω × (0, T ), then v ≡ C. As d[u(t)]/dt = 0, it follows that C = [u 0 ] and that as t → ∞, u(t) → [u 0 ] weakly in H 1 (T)
and strongly in H s (T) for any s < 1.

Remark 1 Similar results, but with convergences towards 0, hold for the system

u t -u xxt + u x + uu x -(a(x)u x ) x = 0, u(0, t) = u(2π, t) = 0, u(x, 0) = u 0 (x)
provided that a(x) > 0 for x ∈ (0, ε) ∪ (2π -ε, 0) for some ε > 0 and the UCP Conjecture holds.

Boundary stabilization of BBM

In this section, we are concerned with the boundary stabilization of the BBM equation. We consider the following system

u t -u xxt + u x + uu x = 0, x ∈ (0, L), t ≥ 0, (2.20)
u tx (0, t) = αu(0, t) + 1 3 u 2 (0, t), (2.21)

u tx (L, t) = βu(L, t) + 1 3 u 2 (L, t), (2.22) u(x, 0) = u 0 (x), (2.23)
where α and β are some real constants chosen so that

α > 1 2 and β < 1 2 •
Scaling in (2.20) by u yields (at least formally) (2.25) u(t) → 0 strongly in H s (0, L) for all s < 1.

d dt 1 2 L 0 (|u| 2 + |u x | 2 )dx = [uu tx ] L 0 - u 2 2 + u 3 3 L 0 = (β - 1 2 )|u(L, t)| 2 + ( 1 2 -α)|u(0, t)| 2 , ( 2
(2.26)

Proof. Let us begin with the well-posedness part. Pick any u 0 ∈ H s (0, L) with s > 1 2 . Let v = u t . Then v solves the elliptic problem

(1 -∂ 2 x )v = f, x ∈ (0, L), (2.27) v x (0) = a, v x (L) = b (2.28)
with f := -u x -uu x , a := αu(0, t) + 1 3 u 2 (0, t), b := βu(L, t) + 1 3 u 2 (L, t). Note that the solution v of (2.27)-(2.28) may be written as

v = w + g where g(x) = ax + b-a 2L x 2 and w = (1 -∂ 2 x ) -1 N (f -(1 -∂ 2 x )g) solves (1 -∂ 2 x )w = f -(1 -∂ 2 x )g w x (0) = w x (L) = 0. Thus u t = v = -(1 -∂ 2 x ) -1 N (u x + uu x ) + (1 -(1 -∂ 2 x ) -1 N (1 -∂ 2 x ))g. ( 2

.29)

We seek u as a fixed point of the integral equation

u(t) = Γ(u)(t) := u 0 + t 0 -(1 -∂ 2 x ) -1 N (u x + uu x )(τ ) +(1 -(1 -∂ 2 x ) -1 N (1 -∂ 2 x )) [αu(0, τ ) + 1 3 u 2 (0, τ )]x +(2L) -1 [βu(L, τ ) + 1 3 u 2 (l, τ ) -αu(0, τ ) - 1 3 u 2 (0, τ )]x 2 dτ.(2.30) Note that D((1 -∂ 2 x ) s 2
N ) = H s (0, L) for -1/2 < s < 3/2. Let R > 0 be given and let B R denote the closed ball in C([0, T ], H s (0, L)) of center 0 and radius R. For 1/2 < s < 5/2 and u ∈ B R , v ∈ B R , we have that (2.24). Thus ||u(t)|| H 1 (0,L) is nonincreasing, and T may be taken as large as desired. Let us turn our attention to the asymptotic behavior. Let (t n ) n≥0 be any sequence with t n → ∞ as n → ∞. Extracting a subsequence if needed, we can assume that t n+1 -t n ≥ T for each n and that, for some v 0 ∈ H 1 (0, L), u(t n ) → v 0 weakly in H 1 (0, L) (hence strongly in H s (0, L) for s < 1) as n → ∞. The continuity of the flow map (which follows at once from the fact that the map Γ is a contraction) yields

||Γ(u)(t) -Γ(v)(t)|| H s (0,L) ≤ C t 0 (1 + R)||u(τ ) -v(τ )|| H s (0,L) dτ ≤ CT (1 + R)||u -v|| C([0,T ],H s (0,L)) and ||Γ(u)(t)|| H s (0,L) ≤ ||u 0 || H s (0,L) + CT (1 + R)||u|| C([0,T ],H s (0,L)) • Pick R = 2||u 0 || H s (0,L) and T = (2C(1 + R)) -1 . Then Γ is a contraction in B R ,
u(t n + •) → v(•)
in C([0, T ], H s (0, L)) for s < with v(x, t) = 0 for x ∈ R \ (0, L), t ∈ (0, T ).

Since v ∈ C([0, T ], H s (R)) with 1/2 < s < 1, we infer from the UCP Conjecture (if true) that v ≡ 0, hence v 0 = 0.

  hence it admits a unique fixed point which solves (2.30), or (2.29) and (2.23). Assume now that s = 1. It follows from (2.29) that u t ∈ C([0, T ], H 2 (0, L)), hence we can scale by u in (2.20) to derive

  1, where v denotes the solution of (2.20)-(2.23) issued from v 0 at t = 0. Integrating (2.24) on [t n , t n+1 ] yields Extending v(x, t) by 0 for x ∈ R \ (0, L) and t ∈ (0, T ), we infer from (2.20), (2.21), (2.22) (for v) and (2.31) that v t -v txx + v x + vv x = 0 for x ∈ R, t ∈ (0, T )

	Letting n → ∞ yields
		(	1 2	-β)	0	T	|v(L, t)| 2 dt + (α -	1 2	)	0	T	|v(0, t)| 2 dt = 0.	(2.31)
	1 2	||u(t n+1 )|| 2 H 1 (0,L) -	1 2	||u(t n )|| 2 H 1 (0,L) +(	1 2	-β)	tn+1 tn	|u(L, t)| 2 dt+(α-	1 2	)	tn+1 tn	|u(0, t)| 2 dt.


Tomonari Mizoguchi

L D C Jaubert

Masafumi Udagawa

Clustering of Topological Charges in Kagome Classical Spin Liquid

Keywords: numbers: 75, 10, Kt

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.

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][7][8][9][10][11][12][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][19][20][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 fartherneighbor 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 × L × 3 sites under periodic boundary conditions, and Ising spin variables, σ z i = ±1 at each site, i. Throughout this Letter, we set k B = 1. We define the Hamiltonian:

H = J 1 i, j n.n. σ z i σ z j + J 2 i, j 2nd σ z i σ z j + J 3 i, j 3rd σ z i σ z j , (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][23][24][25][26], we shall not consider the J 3d 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 (≡ J), where a charge representation is available [27]. This model has offered a rich out-ofequilibrium physics on the three-dimensional pyrochlore lattice [28,29]. To introduce this representation, we classify the triangles on the kagome lattice into "upward" ( ) and "downward" ( ), according to their orientations. We then define a topological charge at each triangle as Q p = η p i∈p σ z i with η p = +1(-1) for p ∈ ( ) [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 ±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]

H = 1 2 -J p Q 2 p -J p,q Q p Q q . (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 p, q . 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'

(b) (c) (d) (e) (f)
; p∈D Q p = i∈∂D σ z i η p D (i) , (3)
which follows the definition of Q p . Here, D is an arbitrary set of triangles and ∂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: 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 largescale 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 n.n. = 0.502. In terms of the charge representation, the ground state manifold consists of the charge configurations with Q = ±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

H = 1 2 (1 + J) p Q 2 p - J 2 p,q (Q p + Q q ) 2 . (4)
For J < 0, the first and second terms can be minimized simultaneously, by setting (i) Q p + Q q = 0 for all the nearestneighbor 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 2 Ti 2 O 7 [30][31][32][33]. Equivalent states have been observed in artificial spin ice [34][35][36][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 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

J = 0.1. (g) Autocorrelation, R spin (t) ≡ i σ z i (0)σ z i (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. sible, 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 = 1 2 (1 -3J)N . N 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 α-th cluster is

E (α) c = 1 2 -J N (α) -Jn (α) i + J 2 n (α) b . (5)
Here, N (α) is the number of triangles in the α-th cluster. n (α) i and n (α) b are the number of internal and boundary sites of the cluster, respectively [Fig. 1

(f)], which satisfy 3N (α) = 2n (α) i + n (α) b . The term J 2 n (α)
b 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 = α E (α)
c . 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 (α) ≤ n (α) b , which provides a lower bound of the system energy, E tot ≥ 1 2 (1 -3J)N . 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 hightemperature peak, around T ∼ 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 ∼ 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 ∼ 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.

(a) 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 equilibration, since the cluster is locally quite stable, and a standard single-spin-flip dynamics freezes. Even a loop-update algorithm [34,35,40] effective for the kagome ice region is powerless 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 equilibrates the system in this region [Fig. 2(g)].

-2 -1 0 1 2 -2 -1 0 1 2 0.5 1 1.5
This non-local dynamics enables us to explore the thermodynamic behavior below the second peak. In particular, we succeed in evaluating the T dependence of entropy precisely, and find that it converges to an unfamiliar value of S HCSL 0 ≈ 0.32, near zero temperature. This residual entropy is well reproduced by applying a variant of Pauling's argument [41] to a nine-spin cluster made of four connected triangles. As a Pauling constraint, we impose that each triangle neighbors with two, and only two, same-sign charges, giving the entropy: S HCSL P = 1 6 log 27 4 ≈ 0.32 [39]. Magnetic structure factor.-The unusual spatial structure of the CSLs can be detected by analyzing their magnetic structure factors, S (q) ≡ i, j σ z i σ z j e -iq•(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)]. Below 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][45][46].

In contrast, at J = 0.1, S (q) develops a quite different pattern. 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 recently 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 recombined 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 configurational pattern leads to a unconventional value of residual entropy. This spatial pattern can be detected through the halfmoon 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 homogeneous, each ground state configuration composed of the hexamer 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 systems. The hexamer CSL also offers a promising extended region (0 < J < 1/3) for the search of unconventional quantum spin liquids [22][23][24][25][26][47][48][49] and quantum order-by-disorder mechanisms [50][51][52][53] beyond the traditional Heisenberg antiferromagnet.

It is tempting to speculate if the interaction between topological charges underlies hexagonal or small-cluster formation found in a broad range of frustrated magnets. For instance, cluster excitations of antiferromagnetic hexagonal loops were reported in the cubic spinel ZnCr 2 O 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 Okinawa Institute of Science and Technology Graduate University. T. M. is supported by Advanced Leading Graduate Course for Photon Science (ALPS).

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 groundstate structure for 0 < J 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. I A and Sec. I B, 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. I C, 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 = η p i∈p σ z i , (1)
where η p = +1(-1) for p ∈ ( ). If we consider a set of triangles, D, and sum Eq. ( 1) over D, we obtain the Gauss' law:

p∈D Q p = i∈∂D η p D (i) σ z i . (2)
As noted in the main text, ∂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 η p D (i) σ z i vanishes for internal sites of D, leaving only the boundary contribution.

If D is composed of N D triangles (of both orientations), the number of internal sites, n D i , and the boundary sites, n D b satisfy the relation:

3N D = 2n D i + n D b .
(3) * mizoguchi@hosi.phys.s.u-tokyo.ac.jp

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

| p∈D Q p | ≤ i∈∂D |η p D (i) ||σ z i | = i∈∂D 1 = n D b . (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 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 α-th cluster composed of N (α) triangles, with n boundary sites. If we apply Eq. (4) to this cluster, we have

N (α) ≤ n (α) b , (5)
since we only consider single charges. Next, we compute the energy of the α-th cluster,

E (α) c = 1 2 -J N (α) -Jn (α) i + J 2 n (α) b , (6)
as in Eq. ( 5) in the main text. This definition allows us to express the total energy of the system, E tot , as the summation of E (α) c over the clusters:

E tot = α E (α)
c . By combining Eqs. (3), ( 5) and ( 6), we obtain

E (α) c ≥ 1 2 (1 -3J)N (α) . (7)
By summing up all the cluster contributions, we obtain the lower bound of the total energy:

E tot ≥ 1 2 (1 -3J)N ≡ E s tot , (8)
with the total number of triangles, N . E s tot gives the lowest energy of the system, given that only single charges exist. E s tot , in fact, gives the ground state energy for small J.

C. Triple-charge cluster

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 samecharge 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 D. This set of triangles satisfies n D b = N D + 2. Accordingly, we obtain

| p∈D Q p | = N D + 4 > n D b . (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 (α) triangles, with n boundary sites, satisfies

N (β) + 2 ≤ n (β) b , (10)
as a result of applying Eq. (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 tot = α∈sc E s,(α) c + β∈tc E t,(β) c . (11)
Here, E s,(α) c

(E t,(β) c
) is the contribution from a cluster without (with) triple charge. In the presence of surrounding triple charges, E s,(α) c 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 s,(α) c = 1 2 -J N (α) -Jn (α) i + J 2 n ss(α) b + 3J 2 n st(α) b , (12)
where n

ss(α) b (n st(α) b
) means the number of boundary sites between single charges (single and triple charges). They satisfy

n ss(α) b + n st(α) b = n (α) b . (13)
Accordingly, Eq. ( 12) can be rewritten as

E s,(α) c = 1 2 -J N (α) -Jn (α) i + J 2 n (α) b + Jn st(α) b .( 14
)
A cluster with a triple charge composed of N (β) triangles have N (β) -1 single, and one triple charges. We assume the triple charge neighbors with n (β) (= 0, 1, 2, 3) single charges within the cluster, and 3n (β) charges of other clusters. Accordingly, we write the cluster energy as

E t,(β) c = 1 2 -J (3 2 + (N (β) -1)) -J[3n (β) + (n (β) i -n (β) )] + J 2 n ss(β) b + 3J 2 n st(β) b + 3J 2 n ts(β) b + 9J 2 n tt(β) b . ( 15
)
Here, the first and the second terms come from the selfenergy of charges and the interaction within the cluster, respectively. The other terms denote four possible types of inter-cluster interactions. n ) denotes the number of sites between single (triple) charges in the β-th cluster and triple (single) charges in other clusters. These numbers satisfy

n ss(β) b + n st(β) b = n (β) b -(3 -n (β) ), (16)
and

n ts(β) b + n tt(β) b = (3 -n (β) ). (17)
Accordingly, Eq. ( 15) can be rewritten as

E t,(β) c = 1 2 -J (3 2 + (N (β) -1)) -J[3n (β) + (n (β) i -n (β) )] + J 2 [n (β) b -(3 -n (β) )] + 3J 2 (3 -n (β) ) + Jn st(β) b + 3Jn tt(β) b . ( 18
)
With these definitions of E s,(α) c and E t,(β) c

, 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 tot = α∈sc 1 2 -J N (α) -Jn (α) i + J 2 n (α) b + β∈tc 1 2 -J (3 2 + (N (β) -1)) -J[3n (β) + (n (β) i -n (β) )] + J 2 [n (β) b -(3 -n (β) )] + 5J 2 (3 -n (β) ) + 2Jn tt(β) b . (19)
Here, we have used

α∈sc n st(α) b + β∈tc (n st(β) b + n tt(β) b ) = β∈tc (3 -n (β) ). ( 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 tot ≥ 1 2 (1 -3J)N + 4 β∈tc (1 -Jn (β) ). ( 21
)
The equality holds if (i)

N (α) = n (α) b , (ii) N (β) + 2 = n (β) b
for α ∈ sc and β ∈ tc, and (iii) there are no triple charges next to each other, n tt(β) b = 0. Since E s tot = 1 2 (1 -3J)N gives the lowest energy given only s-clusters, as defined in Eq. ( 8), the correction of the total energy due to tclusters can be evaluated as

E tot -E s tot ≥ 4 β∈tc (1 -Jn (β) ). ( 22
)
3. Ground-state energy Finally, we discuss the ground-state energy of the system, on the basis of Eq. ( 22). Since n (β) ≤ 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 (β) = n. Since a t-cluster with n (β) = n is composed of at least n + 1 triangles, M n satisfies

M tot ≡ 3 n=0 (1 + n)M n ≤ N . (23)
Then, we can rewrite Eq. ( 22) as

E tot -E s tot ≥ 4 3 n=0 (1 -Jn)M n = -4JM tot + 4(1 + J) 3 n=0 M n . (24)
For any given value of M tot , the total number of clusters

3 n=0 M n is minimized, if M 0 = M 1 = M 2 = 0 and M 3 = M tot /4, leading to E tot -E s tot ≥ (1 -3J)M tot ≥ 0 (0 < J ≤ 1/3), (1 -3J)N (J > 1/3). ( 25
)
This inequality shows that the hexamer CSL is stabilized for 0 < J ≤ 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 /4), it will be the ground state in this region, with the ground-state energy, E tot -E s tot = (1 -3J)N . 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 × 2 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 tot -E s tot = (1 -3J)N , (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][2][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 longrange 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 < 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][2][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.

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.

C. Worm algorithm

In the hexamer CSL realized for 0 < J < 1 3 , lowtemperature 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 II B. As illustrated in Fig. 4, we firstly choose an endpoint 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. II B, 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 system is preserved, such that the attempt is always accepted 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 triangle 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 correspondence between the string update and its reversed process. These conditions ensure the detailed balance of stochastic processes. Below, we show the benchmark of the algorithm. The average length of worms, worm , is shown in Fig. 5. One can see that the average length is almost temperatureindependent 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 relaxation drastically; to show its quantitative evidence, we calculate the autocorrelation.

The autocorrela-

tions of spins, R spin ≡ i [σ z i (0)σ z i (t)] Nsite
, and charges,

R charge ≡ p∈ [Qp(0)Qp(t)] N
, are shown in Fig. 6, respectively. Here one Monte Carlo step consists of N single single spin flip, N loop loop updates, and N worm worm updates. In Fig. 6, these numbers are chosen as (N single , N loop , N worm ) = (N site , 0, 0) for red dots, (N single , N loop , N worm ) = (N site , N site /8, 0) for green dots, and (N single , N loop , N worm ) = (N site , 0, N site ) for blue dots, respectively. One can clearly see that R spin and R charge approach to zero very rapidly, if we introduce the new worm algorithm. On the other hand, the loop update algorithm is not efficient in hexamer CSL state. In particular, R charge is unchanged by loop update, since it does not change the charge configuration.

N worm dependence of the autocorrelations are shown in Fig. 7. Clearly, the larger N worm becomes, the faster the relaxation is achieved. Even in the smallest number of N worm (= Nsite 32 ), R spin and R 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.

III. PAULING'S ARGUMENT IN NEW CLASSICAL SPIN LIQUID

In this section, we discuss the residual entropy of hexamer CSL on the basis of Pauling's argument. As we (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 likecharge 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 = 6 × 9 2 9 = 3 3 2 8 . (27)
The corresponding residual entropy becomes

N c = 1 4 N = N site 6 , (29)
and 2 Nsite 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 FIG. 1 .

 11 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, a 1 and a 2 are represented by the orange arrows. (b) Possible charge states on triangles. The red (blue) dots represent σ z =↑ (↓). 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 >1 3 , whose magnetic unit cell is represented by a black rhombus. (f) One of the ground state configurations for 0 < J <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.

  law

FIG. 2 .

 2 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 ρ 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) ≡ i σ z

FIG. 3 .

 3 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.

  the relation, Eq. (5). Meanwhile, the β-th t-cluster composed of N (β) triangles, with n

FIG. 1 .

 1 FIG. 1. (a) 4-triangle cluster and (b) a magnetic unit cell of the 12-site ordering.

FIG. 2 .

 2 FIG. 2. Schematic picture of the loop update, formed by a closed loop and a dangling tail.

FIG. 3 .

 3 FIG. 3. Schematic picture of the string update.

FIG. 4 .

 4 FIG. 4. Schematic picture of the worm algorithm.

FIG. 5 .

 5 FIG. 5. Temperature dependence of the average length of worms.

FIG. 6 .

 6 FIG. 6.The autocorrelation of the spin and charges. (N single , N loop , Nworm) for red, green, and blue dots are described in the text.

FIG. 7 .

 7 FIG. 7. The autocorrelation of the spin and charges for various values of Nworm.

FIG. 8 .

 8 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.

  where N c denotes the total number of four-triangle clusters,

arXiv:1702.03794v1 [cond-mat.str-el] 13 Feb 2017