Tomonari Mizoguchi
email: mizoguchi@rhodia.ph.tsukuba.ac.jp

Ludovic D C Jaubert

Roderich Moessner

Masafumi Udagawa

Magnetic clustering, half-moons, and shadow pinch points as signals of a proximate Coulomb phase in frustrated Heisenberg magnets

Keywords: numbers: 75, 10, Kt

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Geometrically frustrated magnets provide a stage to realize exotic states of matter, ranging from quantum and classical spin liquids 1-3 , unconventional magnetic ordering with topological response [START_REF] Taguchi | Spin Chirality, Berry Phase, and Anomalous Hall Effect in a Frustrated Ferromagnet[END_REF][START_REF] Onose | Observation of the Magnon Hall Effect[END_REF][START_REF] Nagaosa | Topological properties and dynamics of magnetic skyrmions[END_REF] , and states accompanied by exotic phase transitions [START_REF] Kawamura | Phase Transition of the Two-Dimensional Heisenberg Antiferromagnet on the Triangular Lattice[END_REF][START_REF] Weber | Ising Transition Driven by Frustration in a 2D Classical Model with Continuous Symmetry[END_REF][START_REF] Jaubert | Three-Dimensional Kasteleyn Transition: Spin Ice in a [100] Field[END_REF] . Among them, the disordered Coulomb phase is a canonical example for discrete and continuous spins on the three-dimensional pyrochlore lattice when frustration imposes a local divergence-free constraint [START_REF] Henley | The "Coulomb phase" in frustrated systems[END_REF] , whose exotic character has been drawing considerable interest.

The Coulomb phase consists of degenerate spin configurations in absence of any spontaneous symmetry breaking. The magnetic correlations due to the local divergence-free constraint are characterized by singularities in the static magnetic structure factor, S(q) [START_REF] Isakov | Dipolar Spin Correlations in Classical Pyrochlore Magnets[END_REF] . These singularities, called pinch point, have been observed in the canonical spin ice materials, Ho 2 Ti 2 O 7 and Dy 2 Ti 2 O 7 [START_REF] Bramwell | Spin Ice State in Frustrated Magnetic Pyrochlore Materials[END_REF][START_REF] Fennell | Magnetic Coulomb Phase in the Spin Ice Ho2Ti2O7[END_REF] . On top of exotic correlations, the Coulomb phase supports fractional excitations. These excitations are easy to visualize in spin ice, where Ising spins satisfy the so-called ice rules with two spins pointing inwards and two spins pointing outwards on every tetrahedron in the ground state. A tetrahedron in a "three-in-one-out" or "one-inthree-out" configuration, breaking this ice rule, carries a gauge charge and serves as an elementary fractional excitation. In spin ice, these gauge charges are actually effective magnetic charges [START_REF] Castelnovo | Magnetic Monopoles in Spin Ice[END_REF] . By identifying the spins with their inherent "magnetic field", one can regard the gauge charge as a source or sink of the field, and assign a magnetic charge -2(+2) for "one-in-three-out" ("threein-one-out") tetrahedron states. Magnetic charges are defined from the discrete divergence of the magnetic field, i.e. the number of inward spins minus that of outward spins.

The introduction of the concept of magnetic charges turned out to be quite illuminating. The implicit assumption for a "charge" quantity is two-fold in conventional electromagnetism. Firstly, charge should be conserved. Indeed, the above mentioned magnetic charge in spin ice satisfies a local conservation in the sense that they are always created/annihilated in pair of positive and negative charges. And secondly, opposite charges are expected to attract each others. However, this second property is non-universal. In the canonical spin-ice systems, Ho 2 Ti 2 O 7 and Dy 2 Ti 2 O 7 , opposite charges indeed interact with attractive force, attributed to the longrange dipolar interaction. However, the sign of the force actually depends on the microscopic details of the system.

Indeed, recently, the role of charge interactions are drawing interest in spin ice [START_REF] Udagawa | Out-of-equilibrium dynamics and extended textures of topological defects in spin ice[END_REF][START_REF] Rau | Spin slush in an extended spin ice model[END_REF] and its two-dimensional analog [START_REF] Möller | Magnetic multipole analysis of kagome and artificial spin-ice dipolar arrays[END_REF][START_REF] Chern | Two-Stage Ordering of Spins in Dipolar Spin Ice on the Kagome Lattice[END_REF][START_REF] Chern | Magnetic charge and ordering in kagome spin ice[END_REF][START_REF] Mizoguchi | Clustering of Topological Charges in a Kagome Classical Spin Liquid[END_REF] . If the interaction is chosen unnaturally, i.e., attractive between same-sign charges, the Coulomb phase is destabilized towards the formation of same-sign-charge hexamer clustering [START_REF] Udagawa | Out-of-equilibrium dynamics and extended textures of topological defects in spin ice[END_REF][START_REF] Rau | Spin slush in an extended spin ice model[END_REF][START_REF] Mizoguchi | Clustering of Topological Charges in a Kagome Classical Spin Liquid[END_REF] . The generic tendency to clustering can be naturally understood from the competition of the two-fold properties of charges. Same-sign charges attract each others, but they cannot pair-annihilate due to charge conservation. As a result, they form stable clus-ters. Their proliferation leads to unconventional classical spin liquids. Accompanying the clustering, the magnetic correlation displays a noticeable evolution, characterized by half-moon patterns in S(q), which replaces the pinch point singularities [START_REF] Udagawa | Out-of-equilibrium dynamics and extended textures of topological defects in spin ice[END_REF][START_REF] Rau | Spin slush in an extended spin ice model[END_REF][START_REF] Mizoguchi | Clustering of Topological Charges in a Kagome Classical Spin Liquid[END_REF][START_REF] Robert | Propagation and Ghosts in the Classical Kagome Antiferromagnet[END_REF] .

The rich physics brought by the attraction of charges of the same sign naturally motivates us to generalize its analysis to the system with continuous spins. Indeed, the magnetic charge in the Ising system can be generalized to a conserved magnetic vector in continuous spin systems. With the continuous nature of magnetic vector, one can expect a variety of stable textures beyond hexamer clustering. From this viewpoint, it is interesting to look at a class of spinel oxides with 3d magnetic ions, AB 2 O 4 (A = Zn, Hg, Mg, B = Cr, Fe) [START_REF] Lee | Emergent excitations in a geometrically frustrated magnet[END_REF][START_REF] Tomiyasu | Molecular spin-liquid state in the spin-3 2 frustrated spinel HgCr2O4[END_REF][START_REF] Tomiyasu | Molecular Spin Resonance in the Geometrically Frustrated Magnet MgCr2O4 by Inelastic Neutron Scattering[END_REF][START_REF] Tomiyasu | Emergence of Highly Degenerate Excited States in the Frustrated Magnet MgCr2O4[END_REF][START_REF] Gao | Manifolds of magnetic ordered states and excitations in the almost Heisenberg pyrochlore antiferromagnet MgCr2O4[END_REF][START_REF] Yamada | Interspin iteractions in ZnFe2O4: Theoretical analysis of neutron scattering study[END_REF][START_REF] Kamazawa | Magnetic neutron scattering measurements on a single crystal of frustrated ZnFe2O4[END_REF][START_REF] Tomiyasu | A Spin Molecule Model for Geometrically Frustrated Spinel ZnFe2O4[END_REF][START_REF] Conlon | Absent pinch points and emergent clusters: farther neighbor interactions in the pyrochlore Heisenberg antiferromagnet[END_REF] . These compounds have weak magnetic anisotropy with small spinorbit interaction of 3d ions, and the classical Heisenberg model with farther-neighbor interactions is expected to give a good starting point of analysis [START_REF] Conlon | Absent pinch points and emergent clusters: farther neighbor interactions in the pyrochlore Heisenberg antiferromagnet[END_REF] . Indeed, according to inelastic neutron scattering experiments, this family lacks "pinch points" in the dynamical structure factor, which are characteristic of the Coulomb phase. Instead, diffuse scattering patterns appear at the corners of the Brillouin zone. The diffuse scatterings are attributed to the clustering of small number of spins, coined as "spin molecules". Depending on materials, molecules take the form of hexamers [START_REF] Lee | Emergent excitations in a geometrically frustrated magnet[END_REF][START_REF] Tomiyasu | Molecular spin-liquid state in the spin-3 2 frustrated spinel HgCr2O4[END_REF][START_REF] Tomiyasu | Molecular Spin Resonance in the Geometrically Frustrated Magnet MgCr2O4 by Inelastic Neutron Scattering[END_REF][START_REF] Tomiyasu | Emergence of Highly Degenerate Excited States in the Frustrated Magnet MgCr2O4[END_REF][START_REF] Gao | Manifolds of magnetic ordered states and excitations in the almost Heisenberg pyrochlore antiferromagnet MgCr2O4[END_REF] and dodecamers [START_REF] Yamada | Interspin iteractions in ZnFe2O4: Theoretical analysis of neutron scattering study[END_REF][START_REF] Kamazawa | Magnetic neutron scattering measurements on a single crystal of frustrated ZnFe2O4[END_REF][START_REF] Tomiyasu | A Spin Molecule Model for Geometrically Frustrated Spinel ZnFe2O4[END_REF] . In the work by Conlon and Chalker [START_REF] Conlon | Absent pinch points and emergent clusters: farther neighbor interactions in the pyrochlore Heisenberg antiferromagnet[END_REF] , the lack of pinch points has been attributed to weak, generic, farther-neighbor exchange, inducing hexagonal cluster scattering as observed in experiments.

To address these issues, in this paper, we consider the classical Heisenberg models on kagome and pyrochlore lattices with farther-neighbor interactions, on the highsymmetry line J 2 = J 3 = J, for arbitrary values of J > 0. We focus on cooperative paramagnetic region above magnetic ordering temperature, where the magnetic fluctuations reflect the intrinsic nature of the system, in contrast to the ordering pattern itself, which is susceptible to structural changes or other extrinsic effects.

Our main results are summarized as follows: (i) We found three distinct patterns in S(q): pinch points, halfmoons, and stars. These patterns are counterparts of the topological clusters obtained in the corresponding Ising models. (ii) The three patterns are originated in the structure of softest magnetic modes. (iii) The half-moon and star patterns can be interpreted as a shadow of pinch point, and serves as a signal of proximity to Coulomb phase. (iv) These characteristic patterns also appear in the low-energy region of dynamical structure factors, implying the possibility of experimental detection through inelastic neutron scattering.

The rest of this paper is organized as follows. In Sec. II, we first describe the model, namely the J 1 -J 2 -J 3 classical Heisenberg Hamiltonian on the kagome and pyrochlore lattices. Next, we introduce the theoretical methods; large-N calculations, Monte Carlo simulations and semiclassical Landau-Lifshitz (LL) equation. In Sec. III, the three distinct patterns of the static structure factors and their origin are discussed. Here the main arguments are based on band-structure analyses of the large-N approximation, supported by Monte Carlo simulations on the O(3) Heisenberg model. Section IV is devoted to the spin dynamics, analyzed by the LL equation. In Section V, we discuss the real-space structure of the clusters, and show that they can be understood as a continuation from the topological charge cluster obtained in the Ising counterpart. Finally, we present discussions and summary in Sec. VI. Details of the large-N approximation, Monte Carlo simulations, and the quadrupolar order parameter are described in the Appendices.

II. MODEL AND FORMALISM

A. Model

We consider a Heisenberg model on kagome and pyrochlore lattices with up to third-neighbor interactions:

H =J 1 i,j n.n. S i • S j + J 2 i,j 2nd S i • S j +J 3a i,j 3rd,a S i • S j + J 3b i,j 3rd,b S i • S j = 1 2 n,m µ,ν α=x,y,z S µ,α n H µ,ν n,m S ν,α m . (1)
J 1 , J 2 and J 3a,3b are the exchange coupling constants connecting nearest-, second-nearest-and third-nearest neighbors, as shown in Fig. 1. Notice that two different types of the third-neighbor term (J 3a and J 3b ) are distinguished in that J 3a connects the sites along edges, while J 3b connects those across a hexagon. Here, we have introduced unit cell indices, n, m, and sublattice indices, µ, ν, and expressed each sites as their combinations: i = (n, µ) and j = (m, ν). The unit cells contain respectively triangles (tetrahedra) of one orientation for kagome (pyrochlore) lattice, and the sublattice indices take µ = 1, • • • N sub , with N sub = 3 (4). S µ n = (S µ,x n , S µ,y n , S µ,z n ) are (classical) three-component vectors with unit length |S µ n | = 1. The model with general parameter set has been intensively studied on both kagome [START_REF] Chalker | Hidden order in a frustrated system: Properties of the Heisenberg Kagome antiferromagnet[END_REF][START_REF] Huse | Classical antiferromagnets on the Kagomé lattice[END_REF][START_REF] Reimers | Order by disorder in the classical Heisenberg kagomé antiferromagnet[END_REF][START_REF] Zhitomirsky | Octupolar ordering of classical kagome antiferromagnets in two and three dimensions[END_REF][START_REF] Chern | Dipolar Order by Disorder in the Classical Heisenberg Antiferromagnet on the Kagome Lattice[END_REF][START_REF] Gómez Albarracín | Degenerate and chiral states in the extended Heisenberg model on the kagome lattice[END_REF] and pyrochlore [START_REF] Conlon | Absent pinch points and emergent clusters: farther neighbor interactions in the pyrochlore Heisenberg antiferromagnet[END_REF][START_REF] Reimers | Meanfield approach to magnetic ordering in highly frustrated pyrochlores[END_REF][START_REF] Moessner | Properties of a Classical Spin Liquid: The Heisenberg Pyrochlore Antiferromagnet[END_REF][START_REF] Canals | Spin-liquid phase in the pyrochlore anti-ferromagnet[END_REF][START_REF] Chern | Partial order from disorder in a classical pyrochlore antiferromagnet[END_REF][START_REF] Okubo | Cubic and noncubic multiple-q states in the Heisenberg antiferromagnet on the pyrochlore lattice[END_REF] lattices, putting an emphasis on the low-temperature ordered states. In the present work, we focus on the parameters (J 1 , J 2 , J 3a , J 3b ) = (1, J, J, 0) with J > 0. This parameter set has recently been shown to lead to a clustering of topological charges in the corresponding Ising models for J > 0 (see Sec. V for details [START_REF] Udagawa | Out-of-equilibrium dynamics and extended textures of topological defects in spin ice[END_REF][START_REF] Rau | Spin slush in an extended spin ice model[END_REF][START_REF] Mizoguchi | Clustering of Topological Charges in a Kagome Classical Spin Liquid[END_REF] ). The value of J 1 = 1 sets the energy and temperature scales of our problem.

B. Conserved spins

It is instructive to introduce a local magnetic moment for each triangular and tetrahedral unit, n,

M n ≡ ζ n j∈n S j , (2)
where ζ n = ±1 is a sign factor distinguishing between upward (+1) and downward (-1) triangles/tetrahedra 41 . The Hamiltonian (1) can be rewritten as a function of M n :

H = 1 2 -J n |M n | 2 -J n,m M n • M m , (3)
where the summation over n is taken over both upward and downward triangles/tetrahedra, and the summation over n, m is over neighboring pairs of triangles/tetrahedra. This expression naturally accounts for the Coulomb phase at J = 0, with M n = 0 for all n, and for its stability for small J as will be discussed in details in the next sections. Equation ( 3) is a generalization of the spin to charge mapping of the corresponding Ising systems [START_REF] Udagawa | Out-of-equilibrium dynamics and extended textures of topological defects in spin ice[END_REF][START_REF] Rau | Spin slush in an extended spin ice model[END_REF][START_REF] Mizoguchi | Clustering of Topological Charges in a Kagome Classical Spin Liquid[END_REF][START_REF] Ishizuka | Spontaneous spatial inversion symmetry breaking and spin Hall effect in a spin-ice double-exchange model[END_REF] and satisfies a conservation law:

n∈D M n = j∈∂D S j , (4)
where D is a connected ensemble of triangles/tetrahedra, and ∂D is its contour. The contour ∂D is made of all spins shared between two triangles/tetrahedra, n ∈ D and m / ∈ D. This "Gauss' law" means that the internal structure of magnetic cluster is constrained by the boundary spins. Indeed, in the Ising case, the discrete variant of this Gauss' constraint strictly determines the structure of clusters, and leads to hexamer spin liquids [START_REF] Mizoguchi | Clustering of Topological Charges in a Kagome Classical Spin Liquid[END_REF] .

C. Formalism

We study the static properties of model (1) with (J 1 , J 2 , J 3a , J 3b ) = (1, J, J, 0) and J > 0, by combining classical Monte Carlo simulations and analytical large-N method. We also address the dynamics by semi-classical LL equation. In this section, we introduce the latter two methods, and all details for the classical Monte Carlo simulations will be given in App. B.

Large-N approximation

To investigate static structure factors, we employ a large-N approximation [START_REF] Isakov | Dipolar Spin Correlations in Classical Pyrochlore Magnets[END_REF][START_REF] Conlon | Absent pinch points and emergent clusters: farther neighbor interactions in the pyrochlore Heisenberg antiferromagnet[END_REF][START_REF] Garanin | Classical spin liquid: Exact solution for the infinite-component antiferromagnetic model on the kagomé lattice[END_REF][START_REF] Sen | Coulomb Phase Diagnostics as a Function of Temperature, Interaction Range, and Disorder[END_REF] . The length of classical Heisenberg spins satisfies a hard constraint |S µ n | = 1. In the large-N method, Heisenberg spins S µ n are replaced by soft-spin variables s µ n whose length is constrained on average:

(s µ n ) 2 = 1 3 . (5)
Here disordered phases are assumed with s µ n = 0. The above constraint ( 5) is enforced by introducing a Lagrange multiplier λ which satisfies

1 N site q Tr[λ Î + β Ĥ(q)] -1 = 1 3 , (6)
where the sum runs over all wavevectors q in the Brillouin zone and N site is the total number of sites. Ĥ(q) represents the Fourier transformation of the exchange matrix:

[ Ĥ(q)] µν = m H µν 0,m e iq•(Rm+rν -rµ) , (7)
where R m is the position of unit cell m with respect to the reference 0, and r µ is the position of the sublattice µ within a unit cell. The static structure factor S(q) in this formalism is given as

S(q) = µ,ν s µ (-q)s ν (q) = µ,ν λ Î + β Ĥ(q) -1 µν = N sub η=1 µ,ν [ψ * η (q)] µ [ψ η (q)] ν λ + βε η (q) , (8)
where ε η (q) and ψ η (q) are respectively eigenvalues and eigenvectors of Ĥ(q) with a band index η. Their calculation can be carried out by using the premedial lattices of kagome and pyrochlore, which are respectively the honeycomb and diamond lattices. The main idea of the analytic calculation is to regard the nearest-neighbor (n.n.) exchange interaction of the original lattice as being mediated by the sites of the premedial lattice located in between. We describe this method in Appendix A.

Landau-Lifshitz equation

To investigate the dynamical properties, we numerically solve the following LL equation [START_REF] Robert | Propagation and Ghosts in the Classical Kagome Antiferromagnet[END_REF][START_REF] Conlon | Spin Dynamics in Pyrochlore Heisenberg Antiferromagnets[END_REF][START_REF] Taillefumier | Semiclassical spin dynamics of the antiferromagnetic Heisenberg model on the kagome lattice[END_REF] :

∂S i ∂t = -S i × H eff,i , (9)
where H eff,i is an effective magnetic field at site i given as

H eff,i = ∂H ∂S i =J 1   j: i,j ∈n.n. S j   + J 2   j: i,j ∈2nd S j   +J 3a   j: i,j ∈3rd,a S j   + J 3b   j: i,j ∈3rd,b S j   . (10)
In our simulation, we first prepare the initial states, which are well-thermalized with temperature T , by using single-spin updates with Metropolis algorithm. Then we solve Eq. ( 9) by using fourth order Runge-Kutta method, and compute S(q, ω) as

S(q, ω) = 1 N t Nt l=0 n µ,ν S µ 0 (0)S ν n (lδt) init ×e i[ωlδt-q•(Rn+rν -rµ)] , (11)
where • • • init represents the sample average of independently-prepared initial states; we prepared 864 states for a kagome lattice and 432 states for a pyrochlore lattice. δt is an interval of time, set as 0.01, and N t is the number of steps for the time evolution, set as 100000. In this simulation, the system size we use is 3 × 30 2 spins for a kagome lattice, and 4 × 12 3 spins for a pyrochlore lattice.

III. FOURIER-SPACE ANALYSIS: PINCH POINTS, HALF-MOONS & "STARS"

We will primarily focus our attention to understand the magnetic correlations in the disordered regimes of cooperative paramagnetism. This is why we shall not go into the details of the low-temperature ordered phases, with the exception of the high-symmetry point J = 1/2 on pyrochlore in Section V C, whose nature is particularly enlightening. This approach presents the advantage that, in the cooperative paramagnets, the properties of the kagome and pyrochlore lattices are qualitatively very similar, allowing for a parallel analysis of the two lattices.

The evolution of the correlations in the cooperative paramagnetic regimes are closely linked to the qualitative changes in the band structure obtained by large-N analysis. The discussion in this section relies heavily on the analysis of the low-energy band structure, supported by Monte Carlo simulations at finite temperatures. The outline of the ground-state phase diagram is given in Fig. 2 while the excellent agreement between analytics and numerics is illustrated in Fig. 3.

To briefly introduce the overall structure of phase diagram [Fig. 2], the small-J region, which we call region (I), FIG. 2. Ground-state phase diagram within the large-N approximation. Region (I) represents the Coulomb phase where static structure factor shows pinch points for 0 < J < J1c = 1/5 (kagome) and 1/6 (pyrochlore). For J > J1c, the flat bands are not ground state anymore and the physics is dominated by the energy minima of the dispersive band [Figs. 4 and5], giving rise to half-moon patterns in the static structure factor [Fig. 3]. The high-symmetry point at J2c = 1 (kagome) and 1/2 (pyrochlore) separates the large-J region into two parts, with qualitatively different positions of the energy minima in Fourier space [Figs. 4 and5]. In the structure factor, the half-moons evolve continuously into "star" patterns within region (II). The boundaries have been confirmed by Monte Carlo simulations.

is characterized by the pinch points in structure factor [Fig. 3(a,d)]. As increasing J, the structure factor shows qualitative changes twice. At intermediate values of J, at the beginning of region (II), the structure factor develops characteristic pattern, which we call "half-moon" after its shape [Fig. 3(b,e)]. Further increasing J, moving continuously from region (II) to (III), S(q) shows further change into the "star" pattern. Below, we will introduce the nature of each regions, separately.

A. Flat bands of the Coulomb phase

The first noticeable outcome of the large-N theory is the persistence of the flat band(s) for all values of J [Figs. 4 and5]. These flat bands -one for kagome and two for pyrochlore -are well known from the n.n. model [START_REF] Isakov | Dipolar Spin Correlations in Classical Pyrochlore Magnets[END_REF][START_REF] Canals | Spin-liquid phase in the pyrochlore anti-ferromagnet[END_REF][START_REF] Garanin | Classical spin liquid: Exact solution for the infinite-component antiferromagnetic model on the kagomé lattice[END_REF] . They represent the Coulomb spin liquid where every unit cell (triangle and tetrahedron) bears zero magnetization, {M n = 0 | ∀n}, and appear in the static structure factor as pinch points. The persistence of the flat bands -and their double degeneracy for pyrochlore -are readily understandable from Eq. ( 3) since all configurations of the Coulomb phase with M n = 0 remain degenerate in presence of the farther-neighbor coupling J. Mathematically, this persistence takes the form of the exchange matrix Ĥ(q) being a polynomial of the n.n. exchange matrix [see Apps. A, in particular Eqs. (A2) and (A8)]. As such, the two exchange matrices share the same basis of eigenvectors and the flatness of eigenvalues is transmitted from the latter to the former.

Hence one needs to consider the evolution of the dis- persive bands to understand the qualitative changes in the magnetic correlations as a function of J.

----------------------------------------------

B. From pinch points to half-moons, near J1c

The flat bands form the ground-state manifold up to J = J 1c = 1/5 for kagome and 1/6 for pyrochlore. This The kagome n.n. Heisenberg antiferromagnet (J = 0) is well known for its Coulomb phase at intermediate temperatures, followed by a coplanar regime at lower temperatures selected by thermal order by disorder [START_REF] Chalker | Hidden order in a frustrated system: Properties of the Heisenberg Kagome antiferromagnet[END_REF][START_REF] Huse | Classical antiferromagnets on the Kagomé lattice[END_REF][START_REF] Reimers | Order by disorder in the classical Heisenberg kagomé antiferromagnet[END_REF][START_REF] Zhitomirsky | Octupolar ordering of classical kagome antiferromagnets in two and three dimensions[END_REF] . The Coulomb phase is marked by a plateau in the heat capacity and pinch points in the structure factor. When the system enters the coplanar regime, the heat capacity marks a kink, and sharp peaks of scattering at q √ 3 = (4π/3, 0) appear in the structure factor [Fig. 8]. These peaks represent the onset of the √ 3 × √ 3 longrange order as T → 0 + [34], but are not Bragg peaks since there is no dipolar long-range order at finite temperature.

At J = J 1c , the softening of the band touching between the lowest dispersive and the flat bands enhances thermal fluctuations compared to J = 0. As a consequence, the coplanar regime is pushed down to lower temperatures by an order of magnitude [Fig. 8]. Noticeably, at intermediate temperatures (T = 0.022), the pinch points visible at J = 0 have disappeared in favor of the onset of the characteristic half-moons for J = J 1c . This is why the heat capacity does not show the characteristic plateau of the Coulomb phase at J 1c .

As for the pyrochlore lattice, the value of J 1c = 1/6 obtained from large-N is confirmed by simulations down to T ∼ 10 -4 [Fig. 9]. For J = J + 1c , just above the boundary inside region (II), the system orders. But it remains disordered at the boundary J 1c . As a consequence, the pinch points of the Coulomb phase are visible up to J = J 1c , replaced by half-moons as soon as the system enters region (II) [Fig. 9.(b)]. Please note that the small thickness of the pinch points for J = J 1c is due to the proximity of the half-moon regime at finite temperature.

C. Origin of the half-moons

In region (II), the structure factor develops half-moon patterns [Fig. 3 (b,e)]. What kind of information can be read from this characteristic magnetic scattering ? The intensity of scattering at particular wave vectors q is determined by the energy of magnetic modes and the corresponding real-space structure of dominant modes. In our large-N language, they are described by the shape of the energy-minima manifold, ε min (q), and the weight of the corresponding eigenfunctions, Ξ η (q) ≡ µ,ν [ψ * η (q)] µ [ψ η (q)] ν . The weight, Ξ η (q), satisfies the sum rule:

η Ξ η (q) = N sub , (12)
due to the normalization of spin length.

Basically, the high-intensity points of the half moons follow the position of the energy minima. However, the energy of magnetic modes do not account for everything. On one hand, the energy-minima manifold, ε min (q), is extended in Fourier space. More precisely, it can be defined locally as a hypersurface; a closed line for the twodimensional kagome and a closed surface for the threedimensional pyrochlore [see Figs. 6 and7, and Apps. A 5 and A 6 for more details]. On the other hand, it is clear from Fig. 3.(b,e) that the extension of the half-moons is finite. They terminate at some point and do not form closed curves as expected from the energy-minima manifold. This vanishing intensity was coined as "ghost" excitations for the kagome Heisenberg antiferromagnet [START_REF] Robert | Propagation and Ghosts in the Classical Kagome Antiferromagnet[END_REF] .

This discrepancy, symbolized by a missing arc that should connect the half moons by pair, can be attributed to the spatial character of the magnetic mode. In Fig. 10, we show the intensity maps of Ξ η (q) for a kagome lattice in the second Brillouin zone, surrounding the wave vector 0, 4π that the high-intensity regions of flat mode and lowerdispersive mode complement each other. This complementarity is originated in the sum rule, Eq. ( 12). Around the pinch point, the contribution from the highest-energy band is small, and the sum rule is satisfied only between the flat mode and the lower dispersive mode. The missing arc is attributed to the vanishing weight of lower dispersive band in the bow-tie region, where the flat band contribution is dominant. In this sense, the missing arc of the half-moon can be considered as a shadow of pinch point, explaining the 10-year old open question about the nature of the "ghost" excitations in the kagome Heisenberg antiferromagnet [START_REF] Robert | Propagation and Ghosts in the Classical Kagome Antiferromagnet[END_REF] . This missing arc signals the proximate presence of pinch point, and serves as evidence that the system is in the vicinity of a Coulomb phase.

The same scenario holds for the pyrochlore lattice: the maps are shown in Figs. 11(a)-(c) (centered at [002]) and 11(d)-(f) (centered at [111]). Again, the combination of the energy-minima surface and the intensity map Ξ η (q) gives rise to the half-moon patterns.

This explanation remains valid throughout the phase diagram for J > J 1c , and in particular as the half-moons continuously deform into star patterns.

D. From half moons to star patterns for J1c < J

In the previous subsection, we have seen how the shape of the half-moon is linked to the position of the energy minima in Fourier space. These energy minima continuously move as a function of J, and the shape of the half-moon evolves with them, as illustrated in Figs. 12 and13.

As increasing J, the radius 48 of the half-moons in- For convenience, we shall use the name of "star" for both lattices, which shall be understood as the patterns formed by connected half-moons.

The increase of the half-moon radius in Fourier space, R K,P , can be calculated analytically thanks to the large-N method

R K = 4 √ 3 arccos 1 8 1 + J 2J 2 -1 , (13)
along the ΓM direction for kagome, and

R P = 2 arccos 4J + 1 -28J 2 32J 2 , ( 14
)
along the ΓX direction for pyrochlore. These formulae are in excellent agreement with results obtained from Monte Carlo simulations in the regime of collective paramagnetism, and above any potential transition temperature [Fig. 14].

IV. SEMI-CLASSICAL DYNAMICS

In experiments, signals of anomalous magnetic correlation are sometimes observed in finite-frequency regions, through, e.g., inelastic neutron scattering. For example, it is at finite energy of the kagome Heisenberg antiferromagnet that half-moons were first observed [START_REF] Robert | Propagation and Ghosts in the Classical Kagome Antiferromagnet[END_REF] , before being stabilised as signature of the ground state, at low energy, via farther-neighbor interactions [START_REF] Udagawa | Out-of-equilibrium dynamics and extended textures of topological defects in spin ice[END_REF][START_REF] Rau | Spin slush in an extended spin ice model[END_REF][START_REF] Mizoguchi | Clustering of Topological Charges in a Kagome Classical Spin Liquid[END_REF] . Accordingly, in order to find the half-moons and stars in a realistic experimental setting, it is desirable to estimate the energy scale of corresponding magnetic excitations. In the context of large-N analysis, these magnetic patterns are associated with energy bands [Fig. 4,5], but this band energy cannot be interpreted as the frequency of experimental probes in itself. A relation connecting them was proposed under the assumption of relaxation dynamics [START_REF] Conlon | Absent pinch points and emergent clusters: farther neighbor interactions in the pyrochlore Heisenberg antiferromagnet[END_REF][START_REF] Conlon | Spin Dynamics in Pyrochlore Heisenberg Antiferromagnets[END_REF] . However, it is not a priori obvious if this assumption holds in our system. To this end, in this section, we address the dynamics of the system by solving the semiclassical LL equation in Eq. ( 9) and calculating the dynamical structure factor S(q, ω). Here we focus on regions (II) and (III) for both kagome and pyrochlore lattices.

Let us first see the results for a kagome lattice. Cuts for several frequencies are shown in Figs. 15. As expected, the characteristic patterns observed for S(q) are obtained in the low-energy sectors in both regions. In the region (II), as clearly seen in Fig. 15(a), the half-moon signal appears in the quasi-elastic regime: ω = 0, showing that this pattern dominates the long-time behavior of magnetic correlation in this region. The pattern sustains with small ω dependence in the low-frequency region. As going to intermediate energy scale comparable to n.n. 13) and ( 14)]. The error bars come from the discretisation of the Fourier space in a system of finite size (L = 30). For kagome, data were taken at T = 0.32. For pyrochlore, data were taken for a range of temperature above the transition temperature, T ∈ [0.08; 0.24].

coupling, ω ∼ 1, the signal smears out. Similarly, in the region (III), the star shape appears in the low-frequency part of the dynamical structure factor [Fig. 15(c)], and it gradually blurs towards higher energy. The remnant of star pattern remains observable up to higher energy, compared with the vanishment of half-moons in the region (II), probably attributed to the growing energy scale of magnetic modes as J, as implied by the larger band width obtained in the large-N analysis [Fig. 4 and5].

The same trend is also seen for a pyrochlore lattice, as shown in Fig. 16. The half-moon and the star patterns are clearly visible in each region. The results of both lattices show that the shadow of pinch points can be observed through the excitations in finite-frequency range, i.e., the proximity to Coulomb phase can be captured through the inelastic neutron scattering experiment.

V. REAL SPACE PICTURE: MAGNETIC CLUSTERING

In this section, we will address the real space picture, accompanying the characteristic patterns in the structure factors. We will show that half-moons and stars re- flect the formation of magnetic clusters. These magnetic clusters can be associated with the cluster of topological charges obtained in the Ising systems, through the analogy of conserved spin introduced in Eq. ( 2), with the topological charge defined in the Ising system [START_REF] Udagawa | Out-of-equilibrium dynamics and extended textures of topological defects in spin ice[END_REF][START_REF] Rau | Spin slush in an extended spin ice model[END_REF][START_REF] Mizoguchi | Clustering of Topological Charges in a Kagome Classical Spin Liquid[END_REF] .

A. Comparison with Ising systems

The half-moons and stars in S(q) are also seen in the corresponding Ising model [START_REF] Udagawa | Out-of-equilibrium dynamics and extended textures of topological defects in spin ice[END_REF][START_REF] Rau | Spin slush in an extended spin ice model[END_REF][START_REF] Mizoguchi | Clustering of Topological Charges in a Kagome Classical Spin Liquid[END_REF] :

H = i,j n.n. σ z i σ z j + J i,j 2nd σ z i σ z j + J i,j 3rd,a σ z i σ z j . (15)
For Ising degrees of freedom σ z i , the Hamiltonian can be re-written in terms of local topological charges on each triangle (for kagome) or tetrahedron (for pyrochlore),

Q n = ζ n i∈n σ z i , (16)
with ζ n = +(-)1 for a upward (downward) triangle/tetrahedron. The possible values of charges are Q n = {±3, ±1} for a kagome lattice; Q n = {±4, ±2, 0} for a pyrochlore lattice. The Hamiltonian then becomes [START_REF] Udagawa | Out-of-equilibrium dynamics and extended textures of topological defects in spin ice[END_REF][START_REF] Rau | Spin slush in an extended spin ice model[END_REF][START_REF] Mizoguchi | Clustering of Topological Charges in a Kagome Classical Spin Liquid[END_REF][START_REF] Ishizuka | Spontaneous spatial inversion symmetry breaking and spin Hall effect in a spin-ice double-exchange model[END_REF]

H = 1 2 -J n Q 2 n -J n,m Q n Q m + (const.). ( 17
)
The vector field defined in Eq. ( 2) and the Hamiltonian form of Eq. ( 3) were natural extensions of these discrete topological charges to continuous degrees of freedom. We see in Eq. ( 17) that J couples the n.n. charges. J > 0 means that same-sign charges attract each other.

The static structure factors of this Ising model show similar features to those of the Heisenberg model. The half-moons appear for 0 < J < 1 3 on kagome [START_REF] Mizoguchi | Clustering of Topological Charges in a Kagome Classical Spin Liquid[END_REF] , and for J ∼ 1 4 on pyrochlore [START_REF] Udagawa | Out-of-equilibrium dynamics and extended textures of topological defects in spin ice[END_REF][START_REF] Rau | Spin slush in an extended spin ice model[END_REF] . In both lattices, the origin of the half-moons is to a large extent due to the formation of "hexamers" [Figs. 17(e) and 17(f)], which correspond to closed loops made of at least six charges of the same sign. Branches of same-sign charges are then attached to these central hexamers, forming disordered spin-liquid phases made of large clusters of topological charges.

Here, we show that the star patterns discussed in this paper also appear in the large-J region of the Ising models: J > 1 3 for a kagome lattice [Fig. 17(a)] and J > 1 4 for a pyrochlore lattice [Fig. 17(b)]. At lower temperatures, the system orders into phases paved by small clusters of charges, with a maximal charge at the centre, surrounded by smaller charges of the same sign [Fig. 17(c,d)].

From this point of view, the passage from the halfmoons to the stars in the static structure factor corresponds to the evolution from a disordered phase made of hexamers to the ordered phase of smaller clusters centered around a maximal charge. The similarity of halfmoons/stars between Ising and Heisenberg models suggests that short-range correlations similar to topological charge clusters also develop in the Heisenberg models, even though topological stability, with a discretised value of the topological charge, cannot be expected for the continuous spin systems. The motivation of the next sections will be to make this idea more quantitative.

B. Conserved spin correlator

To characterize the real space structure in the Heisenberg models, we focus on the conserved spin M n as a vector-field analogue of the topological charge Q n in the Ising models. The momentum-space correlator of M n is defined as

S M (q) ≡ N sub N site n,m M n • M m e -iq•(Rm-Rn) , ( 18
)
where R n,m is the coordinate at the center of the triangle/tetrahedron where M n,m is defined. Within the large-N approximation, S M (q) is represented as

S M (q) = N sub N site k,k = , ζ k ζ k N sub η=1 M η k (-q)M η k (q) = N sub η=1 µ,ν [ψ * η (q)] µ [ψ η (q)] ν [λ + βε η (q)] • F µ,ν (q). ( 19
)
Here, k = , denotes the direction of triangle/tetrahedron, r c ≡ 1 N sub µ r µ is the coordinate at the center of the upper triangle/tetrahedron,

M η k (q) = µ e iq•[ζ k (rµ-rc)] [ψ η (q)] µ , (20)
is the conserved spin of η-band, and

F µ,ν (q) = k,k = , ζ k ζ k e iq•[ζ k (rν -rc)-ζ k (rµ-rc)] , (21)
is the additional form factor. Figure 18 shows S M (q) obtained by the large-N approximation. In the region (I), S M (q) becomes very small as decreasing temperature, due to the divergence-free nature of the Coulomb phase.

In large-N sense, the flat band does not contribute to the conserved spin correlator since M µ k (q) = 0 for flat bands. For regions (II) and (III), since S M (q) is written by a linear combination of sublattice-resolved structure factors, it shows the characteristic patterns reminiscent of the static structure factor S(q). While the shadow pinch points are absent due to the additional form factor F k,k µ,ν (q), the high-intensity points of S M (q) traces the trajectory of ε min (q) in Fourier space.

The real-space correlator,

S M (R) ≡ M 0 • M R , ( 22
)
is defined on dual honeycomb (diamond) lattice for the kagome (pyrochlore) case, and given by the inverse Fourier transformation of S M (q). We show S M (R) for both cases, in Figs. 19 (kagome) and 20 (pyrochlore). The site 0 is shown with a white circle. The red (blue) dot means that the correlation M 0 • M R takes positive (negative) value, and its radius denotes the rescaled absolute value (see the captions of the figures).

As shown in Fig. 19 (a), in the region (II), positive correlation develops in n.n. sites, as can be expected from the M n -representation of Hamiltonian [Eq. ( 3)]. Moreover, a noticeable amount of correlation develops beyond n.n sites, especially surrounding hexagons, implying the clustering of spins reminiscent of the hexamer cluster in the Ising case, made of same-sign charges surrounding an hexagon. In the Heisenberg case, instead of the charge, the conserved spin shows substantial positive correlation around an hexagon.

The qualitative difference of patterns between region (II) with half-moons and region (III) with stars indicates that different types of clusters evolve in these two regions. For region (III) the positive n.n. correlation is surrounded by the negative correlations, which is reminiscent of the crystalization of double/triple-charge clusters in the Ising models [START_REF] Udagawa | Out-of-equilibrium dynamics and extended textures of topological defects in spin ice[END_REF][START_REF] Mizoguchi | Clustering of Topological Charges in a Kagome Classical Spin Liquid[END_REF] .

Despite the qualitative similarity of cluster structures between Heisenberg and Ising cases, there is one significant difference. In the Ising case, the cluster shapes are rigidly fixed in the region (II) and (III), respectively, due to the discrete spin nature of Ising spins, and do not change with J. Meanwhile in the Heisenberg case, the continuous spins allow continuous modification of clusters, and their amplitudes and cluster sizes also change continuously with J. For example, in the region (II), the cluster can be considerably long-ranged near the boundary with region (I), as implied by the small half-moon radius in the Fourier space [Fig. 14]. C. Gauge-charge ordering at the high-symmetry point J = J2c on pyrochlore

We will confirm the change of clustering patterns between regions (II) and (III) via Monte Carlo simulations, using the pyrochlore model at the boundary J = J 2c as a working example.

What happens at J 2c ? Within the large-N approximation, this is where the energy-minima manifold changes its topology, as illustrated in Figs. 6 and7. The manifold moves from enclosing the Γ point in region (II) (J 1c < J < J 2c ) to the zone corners in region (III) (J 2c < J). The model at J 2c is thus a high-symmetry point of our Hamiltonian. As we will see in this section, it confers to the J 2c boundary an advantage of simplicity particularly useful to characterise the low-temperature ordered state.

Since the energy minima cover an extended region in Fourier space, the ordering mechanism is necessarily via thermal order by disorder. For most of the phase diagram when J > J 1c , the continuous evolution of the energyminima manifold implies incommensurate order at low

M L = 1 N N i=1 S i e ıri•q L (23)
displays a clear first-order jump at the transition [Fig. 21]. Since there are eight L points for each Brillouin zone, and that each L point is shared between two adjacent Brillouin zones, the saturated value of M L is 2/8 = 1/4. Measurements of M L are especially difficult to thermalize below the transition. M L seems to converge towards its saturated value of 1/4, possibly via a second transition at very low temperature. However, further work is necessary to confirm this point. It is also possible that the order parameter M L does not saturate. Measurements of the quadrupolar order parameter M Q on the other hand thermalize quite easily to its saturated value of 2/ √ 3 at zero temperature [see App. C for the definition].

In the rest of this section, we will focus on the nature of the ground state. The saturation of the quadrupolar FIG. 22. Ground state at the boundary J = J2c on the pyrochlore lattice displayed over 32 tetrahedra. All spins are collinear, pointing either up (blue sphere) or down (red sphere). The color of the tetrahedra represent their effective topological charges as defined in Eq. ( 16): +4 (dark blue), +2 (light blue), 0 (white), -2 (orange), -4 (red). Each cubic unit cell is composed of a small cluster of five charges (including one maximal charge) and three zero charges (acting as a vacuum separating charges of opposite sign). This intervening vacuum allows for a global shift of the clusters, illustrated between the left and right panels. The spin configurations are split into two parts by an horizontal semi-transparent plane. Between the two panels, all spins above the semi-transparent plane have been shifted along the [110] direction to the next tetrahedron (at a distance of two nearest neighbours). Because the contact between charges is not modified, the energy is the same. Doing the same shift a second time gives back the initial state. From the point of view of discrete degrees of freedom, this gives a sub-extensive entropy to the ground state.

order parameter M Q [Fig. 21.(b)] indicates that all spins are collinear in the ground state. This collinearity allows to temporarily forget the continuous nature of the classical Heisenberg spins and to consider them as Ising degrees of freedom, pointing either up or down. In analogy with Sec. V A, each tetrahedron bears an effective charge Q n = {0, ±2, ±4}. At J 2c = 1/2, the "chemical potential" (i.e. the first term) of Eq. ( 17) disappears and the ground-state energy, E GS , only comes from the n.n. coupling between charges

E GS = - 1 2 n,m Q n Q m . ( 24
)
This coupling is attractive between same-sign charges. Because of the staggering value of ζ n between up and down tetrahedra in Eq. ( 16), this attraction does not give rise to ferromagnetism. The minimization of the energy is done by paving the lattice with the small clusters of Fig. 17.(d); a central charge +4 or -4 surrounded by four charges +2 or -2 respectively. This state is illustrated in Fig. 22 over four cubic unit cells (32 tetrahedra). In order to avoid contact between charges of opposite sign, each cluster is separated from the other ones by a vacuum (zero-charge tetrahedra). Furthermore, each cluster fits within the eight tetrahedra of a cubic unit cell: one maximal charge Q n = ±4, four charges Q n = ±2 and three vacuum tetrahedra Q n = 0. In order to respect the global neutrality of the system, there must be as many positive as negative clusters. Since the centres of the cubic unit cells form, by definition, a bipartite cubic lattice, the global neutrality is enforced by a staggering arrangement of the clusters, alternatively positive and negative. The magnetic unit cell of the ground state is made of 32 spins (16 tetrahedra). This structure is responsible for the peaks at (h,k,l)= (1/2, 1/2, 1/2) in the structure factor of Fig. 13 [bottom right panel]. In the actual simulations, however, the difficulty of thermalization remains in the low temperatures, because the order parameter M L does not saturate completely. We believe this is due to the sub-extensive entropy of the ground state, as illustrated in Fig. 22. We explain the origin of the degeneracy of the (1/2,1/2,1/2)state in App. D.

VI. SUMMARY AND DISCUSSIONS

We have investigated the magnetic correlations of Heisenberg models with antiferromagnetic farther-neighbor interactions on kagome and pyrochlore lattices in their cooperative paramagnetic regions. For both lattices, we found three distinct patterns of the static structure factor, S(q), namely pinch points, half-moons, and stars. We clarified the origin of these patterns by combining the band structure analysis based on a large-N approximation, and Monte Carlo simulations.

Among the above characteristic patterns, the pinch point serves as a direct evidence of a Coulomb phase. This singular feature of S(q) comes from the structure of flat band in the large-N analysis. Meanwhile, the latter two patterns show that the system is not in the Coulomb phase any more. The half-moon and star patterns can be interpreted as complementary of pinch point, and signals the system is in the vicinity of a Coulomb phase. The half-moon signal comes from the lower dispersive band in the large-N analysis, and it is attributed to the structure of energy-minima manifold, and the spatial structure of corresponding magnetic mode. The star pattern can be obtained from half-moon, via the topological change of energy-minima manifold in Brillouin zones. In the regions where half-moon and star patterns dominate the low-energy physics, a shadow of the pinch points remains observable at finite energy [Figs. 10 and 11].

In real-space perspective, the half-moon and star patterns signal the formation of magnetic clusters. These magnetic clusters are characterized by the short-range correlation of conserved spin, which is analogous to the topological charge defined for the corresponding Ising system. Through this analogy, the half-moon and star cluster can be associated with the hexamers and triplecharge cluster obtained in the Ising system, respectively.

The analogy to Ising system can be extended to the ordering at the high-symmetry point, J = J 2c on pyrochlore system, where the structure of low-temperature ordered phase can be well understood through the concept of topological charge. In contrast, the difference from the Ising system is found in the rigidity of cluster structure: While the clusters are rigidly fixed due to the discreteness of spins in the Ising system, in the Heisenberg system, the cluster shape is flexibly changed as varying J, due to the continuous nature of spin degrees of freedom.

We further addressed the dynamical properties of the model by solving the semiclassical LL equations. As a result, we found that the characteristic half-moon and star patterns appear in the frequency-resolved structure factors, in particular in the low-frequency regime, which means the magnetic clusters dominate the long-time behavior of the dynamics. Obtained patterns in S(q, ω) can be directly accessible through the experimental probes, such as inelastic neutron scattering.

In fact, in pyrochlore compounds, several types of magnetic clusterings have been reported. In ZnCr 2 O 4 [START_REF] Lee | Emergent excitations in a geometrically frustrated magnet[END_REF] and MgCr 2 O 4 [START_REF] Tomiyasu | Molecular Spin Resonance in the Geometrically Frustrated Magnet MgCr2O4 by Inelastic Neutron Scattering[END_REF][START_REF] Tomiyasu | Emergence of Highly Degenerate Excited States in the Frustrated Magnet MgCr2O4[END_REF][START_REF] Gao | Manifolds of magnetic ordered states and excitations in the almost Heisenberg pyrochlore antiferromagnet MgCr2O4[END_REF] , six-spin composites dominate the low-energy excitations. While the proposed spatial structure is different from the hexamers obtained in our analysis, our model clearly gives a route to similar clustering around hexagons. It is tempting to point out the possibility that our hexamers may be continuously connected to the lowenergy excitations observed for these materials. As to the half-moon signal corresponding to the hexamer type clusters, recently, it was theoretically proposed for a doublelayered kagome material [START_REF] Pohle | How many spin liquids are there in Ca10Cr7O28?[END_REF] .

It is also worthwhile to look at another pyrochlore compound, ZnFe 2 O 4 . The cluster excitation observed for this compound takes a "dodecamer" form, consisting of twelve spins [START_REF] Tomiyasu | A Spin Molecule Model for Geometrically Frustrated Spinel ZnFe2O4[END_REF] . The spatial structure of this excitation is common in the triple-charge cluster obtained in the region (III) in our model. Interestingly, for this compound, large farther-neighbor coupling, J 3a > J 1 , is expected [START_REF] Yamada | Interspin iteractions in ZnFe2O4: Theoretical analysis of neutron scattering study[END_REF][START_REF] Kamazawa | Magnetic neutron scattering measurements on a single crystal of frustrated ZnFe2O4[END_REF][START_REF] Tomiyasu | A Spin Molecule Model for Geometrically Frustrated Spinel ZnFe2O4[END_REF] . Different dodecamer structure reminiscent of the kagome hexamer was also proposed for HgCr 2 O 4 [START_REF] Tomiyasu | Molecular spin-liquid state in the spin-3 2 frustrated spinel HgCr2O4[END_REF] .

In conclusion, we found new characteristic patterns in magnetic structure factors, complementary to pinch points, which signal the proximity to Coulomb phase. These patterns signal the formation of magnetic clusters, analogous to the low-energy excitations observed for pyrochlore compounds.

Polynomial expression

We consider kagome and pyrochlore lattices on the same footing, and start with introducing an N × N incident matrix, δ(1) ≡ ĥ, where we write N = N site for brevity. Each row j, and column indices j correspond to the sites of the lattice, and the matrix element takes [ δ(1) ] jj = [ ĥ] jj = 1 if j & j are connected 0 otherwise (A1)

The Hamiltonian matrix can be expressed as Ĥ = ĥ at J = 0, supposing J 1 = 1 and J 2 = J 3 = J. Generalizing δ(1) , we introduce a matrix, δ(n) , so that the element [ δ(n) ] jj takes 1, if and only if the two sites, j and j are n Manhattan distance away, and otherwise, 0. Since the squared incident matrix ( δ(1) ) 2 connects any two sites, where one can be reached from the other in two hoppings, one can obtain

( δ(1) ) 2 = z ÎN×N + x δ(1) + δ(2) , (A2)
with z = 2(N sub -1) is a number of coordination, and x = N sub -2 is a number of paths through which, starting from a site, one reaches a n.n. site of that site by two other n.n. moves. Hereafter, Î × represents the × identity matrix. Obviously, δ(2) corresponds to the part of Hamiltonian matrix describing the second-and third-neighbor interactions, so we can express Ĥ = δ(1) + J δ(2) = + J( δ(1) ) 2 -xJ δ(1) -zJ ÎN×N = (1 -xJ) ĥ + J ĥ2 -zJ ÎN×N .

(A3)

Now the Hamiltonian matrix is expressed as a polynomial of incident matrix, ĥ, the eigenvalue problem of Ĥ is reduced to that of ĥ.

Dual lattice

To solve the eigenvalue problem of ĥ, it is convenient to introduce dual lattice. For clarity, we focus on a kagome lattice, first. We start with constructing an intermediate lattice, by placing new sites on the centers of triangles, and connecting the new sites and neighboring old sites, while erasing the original bonds of kagome lattice. Secondly, from this intermediate lattice, we erase the original sites of the kagome lattice, and obtain a honeycomb lattice as a dual lattice [Fig. 23(a)]. As a dual lattice of pyrochlore lattice, we obtain a diamond lattice in a similar way [Fig. 23(b)].

The dual lattice shares the same unit cell as the original lattice. Below, we adopt the following lattice conventions. For a kagome lattice, as lattice vectors, we choose a (K) 1 with θ G (q) = argG(q). Then, we immediately obtain the eigenvalues and eigenvectors of (A15) in the following manner. First, the eigenvalues of h K 1 (q) are identical with those of h H 1 (q) up to a constant, and thus, using Eq. (A8) and the fact that z = 4, x = 1 for a kagome lattice, the two eigenvalues of (A15) are obtained as

ε K ± (q) = J|G(q)| 2 ± (1 + J)|G(q)| + 1 -4J. (A20)
Next, to obtain the corresponding eigenvectors, we consider a rectangular matrix ĥK→H (q) as we have discussed int Sec. 3 . Then, the eigenvectors are obtained as

ψ K ± (q) = h K→H (q)ψ H ± (q) ||h K→H (q)ψ H ± (q)|| . (A22)
Note that the remaining flat mode is orthogonal to two dispersive modes, and its eigenvalue is -2(1 -J). Energy minima.-We now obtain the eigenvalues of (A15), so let us discuss the properties of obtained band structure. The energy minima of the lower dispersive band can be found by solving

∂ε K -(q) ∂q x = ∂ε K -(q) ∂q y = 0, (A23)
By using the expression (A20), one obtains

∂ε K -(q) ∂q i = [2J|G(q)| -(1 + J)] ∂|G(q)| ∂q i . (A24)
Therefore, (A23) is satisfied when

(i) |G(q)| = 1+J 2J , (ii) ∂|G(q)| ∂qx = ∂|G(q)| ∂qy = 0.
Condition (ii) is satisfied at Γ, K, and M points, but they do not become energy minima. So let us examine (i). Notice that the solution of (i) in q space forms lines, rather than a set of discrete points. The solution evolves as follows. First, when 0 ≤ J ≤ 1 5 , (i) does not have solutions: in this region the flat band has the lowest energy and the static structure factors is determined by a flat band, which gives rise to pinch points in S(q). Second, when 1 5 ≤ J ≤ 1, the solution is given by a closed path enclosing Γ point [see Fig. 6(a)]. Finally, when J ≥ 1, the solution is given by a closed path enclosing K points [see Fig. 6(b)]. As we have seen in the main text, the shape of energy minima is reflected to the characteristic shape of the static structure factor, namely half-moons and stars. Derivation of phase boundaries.-The change of topology of energy minima surface can be detected by looking at ΓM line. On ΓM line, we can parametrize the momentum as q = (Q, Q √ 3 ) with 0 ≤ Q ≤ π. Then |G(q)| is given by

|G(q)| = 8 cos 2 Q 2 + 1, (A25)
and the solution of (i) is then given by

Q = 2 cos -1 1 8 1 + J 2J 2 -1 . (A26)
Figure 24 shows u K ≡ 1 8 1+J 2J

2 -1 as a function of J. In order that Q is between 0 and π, u K has to be between 0 and 1, which is represented by the shade. One can see that the lower bound is J = 1 5 , and the upper bound is J = 1, which correspond to J 1c and J 2c , respectively.

Nearly isotropic nature of half-moon.-The energy dispersion of the lower-dispersive band around Γ point is isotropic. Indeed, ε K -(q) can be expanded around Γ point as

ε K -(q) ∼ -2(1 -J) + 1 4 - 5J 4 q 2 + - 1 192 + 17J 192 q 4
+O(q 5 i ), (A27) with q = q 2 x + q 2 y . Equation (A27) shows that the energy dispersion is isotropic up to the order of q 4 , and this leads to nearly circular shape of the energy-minima surface in region (II) [see Fig. 6(a)].

Pyrochlore lattice

The same method can be applied to a pyrochlore lattice, so here we outline the calculations. We consider the In this section, we explain the origin of the subextensive entropy of (1/2,1/2,1/2)-state at J = J 2c , illustrated in Fig. 22. For this state, which consists of the double-charge clusters, the energy of the ground state comes from the contact between charges [Eq. (24)]. These contacts only take place within each cluster. Let us recall there is one double-charge cluster per cubic unit cell, with one double charge, four single charges, and three vacuum tetrahedra. It means that any change which does not affect the integrity and connectivity of the clusters is iso-energetic. Thanks to the vacuum surrounding every cluster, such changes are possible by shifting an entire plane of clusters. An example is given in Fig. 22. The plane of clusters in the upper part of the figure is shifted in the [110] direction between the left and right panels, while the bottom part remains fixed -the thickness of a plane is exactly one cubic unit cell. More precisely all spins along the [110] lines are shifted by a distance of two nearest neighbors; the spins along the [1 10] lines are left unchanged. Thanks to the intervening layers of vacuum tetrahedra (the one just below the transparent plane, and the one at the top of the figure), this shift does not affect any n.n. pair of charges. The resulting state is thus also a ground state. A second shift in the [110] direction gives back the initial state. The same is also true if one shifts the spins in the upper plane along the [1 10] lines, leaving the [110] lines unchanged. The addition of both shifts makes a fourth possibility. This gives 4 L ground states for a system of L planes in the [001] direction.

The same reasoning applies for planes orthogonal to the [010] and [100] directions. However, it is not possible to do successively a shift in a (100) plane followed by a shift in a (001) plane. This is because the conservation of the energy depends on the intervening vacuum layers. Visually, a plane of clusters can glide at no energy cost as long as there is a layer of vacuum to isolate it from the two planes above and below. But the shift of a (100) plane intersects the layers of vacuum orthogonal to the [001] direction; the shift of a (001) plane is now forbidden in the ground state. As a result, the entropy of the ground state is sub-extensive, of the order of L ln 4.

The configurations of Fig. 22, as well as the shift of entire planes, have been observed in snapshots of Monte Carlo simulations at low temperatures, up to fluctuations away from collinearity. The presence of the shifting planes favors a given cubic axis but not a given direction. This spontaneously breaks rotational symmetry as measured by the quadrupolar order parameter M Q . But in this picture, the correlations along the preferred cubic axis should be "paramagnetic" and the order parameter M L should vanish like 1/L in the thermodynamic limit. This is not what is observed in simulations. The reason is because the system is not made of Ising spins, but of continuous Heisenberg spins. The ground state is selected via order by disorder because of soft modes around the L points of the Brillouin zone, which favor the longrange ordered states where none of the planes are shifted.

However, the sub-extensive entropy is probably responsible for the difficulty of the simulations to thermalize at very low temperatures.

FIG. 1 .

 1 FIG. 1. Schematic picture of the model in Eq. (1) for (a) a kagome lattice and (b) a pyrochlore lattice. Red, blue, green and purple lines denote, respectively, J1, J2, J3a and J 3b . In this paper, we consider J1 = 1, J2 = J3a = J and J 3b = 0.

FIG. 3 .

 3 FIG. 3. Structure factor for kagome (top) and pyrochlore (bottom) lattices calculated by (a)-(c) large-N approximation and (d)-(f) Monte Carlo simulation. White lines denote the Brillouin zones. The model being antiferromagnetic, all the characteristic features of the scattering appear in the secondary Brillouin zones. For region (I), the pinch points are clear signatures of the divergence-free condition of the Coulomb phase. Their absence in regions (II) and (III) indicates that the system is out of the Coulomb phase. The complementary patterns of the pinch points are the half-moons (II), which adiabatically evolve into "star" shapes. The star patterns appear in region (II) and persist in (III). See Fig. 2 for the boundaries of the three regions.

FIG. 4 .

 4 FIG.4. Kagome: Band structures of H(q) on high-symmetry lines (the first row), and on a two-dimensional momentum space (the second row) for several values of J. In the second row the highest band is omitted for clarity. The position of energy minima are shown with black arrows.

FIG. 8 .

 8 FIG. 8. Kagome: Top: Heat capacity C h showing the lowtemperature bump/kink into the coplanar regime at T ≈ 0.005 for J = 0 (black dots) and at T ≈ 0.0004 for J1c = 1/5 (colored triangles). Bottom: Temperature evolution of the static structure factor S(q) for J = 0 (left) and J1c = 1/5 (right), obtained by the classical Monte Carlo simulation.

FIG. 9 .

 9 FIG. 9. Pyrochlore: (a) Heat capacity C h showing the lowtemperature ordering when J is just above J1c, but not at J1c = 1/6, down to T ∼ 10 -4 . (b) Static structure factor S(q) in the [hhl] plane at T = 0.005 showing the pinch points at J1c replaced by very small half moons above J1c.

√ 3 .

 3 FIG. 10. Kagome: The weights of eigenvectors, Ξη(q), in the second Brillouin zone centered at 0, 4π √ 3 , for (a) a dispersive band with lower energy, (b) a dispersive band with higher energy, and (c) a flat band.

FIG. 12 .

 12 FIG.12. Kagome: Evolution from half-moon to star patterns in the static structure factor S(q) for J1c < J < J2c, obtained by Monte Carlo simulation. The simulation temperature is T = 0.32.

FIG. 13 .

 13 FIG. 13. Pyrochlore: Evolution from half-moon to star patterns in the static structure factor S(q) for J1c < J < J2c, obtained by Monte Carlo simulation. The bottom right panel is below the transition temperature at the high-symmetry point J = J2c with Bragg peaks at qL = 1 2 1 2 1 2 .

FIG. 14 .

 14 FIG. 14. The evolution of the radius of the half moons in Fourier space (normalized by 2π) for pyrochlore (a) and kagome (b) agrees quantitatively between numerics (red circles) from classical Monte Carlo simulations and analytics (black line) derived from large-N calculations [Eqs. (13) and (14)]. The error bars come from the discretisation of the Fourier space in a system of finite size (L = 30). For kagome, data were taken at T = 0.32. For pyrochlore, data were taken for a range of temperature above the transition temperature, T ∈ [0.08; 0.24].

FIG. 15 .

 15 FIG. 15. Kagome: ω-cuts of S(q, ω) for (a,b) J = 0.3, T = 0.05 and (c,d) J = 1.1, T = 0.275.

FIG. 17 .

 17 FIG. 17. The star patterns in the structure factor for Ising models on the (a) kagome and (b) pyrochlore lattices, obtained by Monte Carlo simulations. The parameters J and T are described in the figure. White lines denote the Brillouin zones. The corresponding clusters of topological charges are shown in (c) and (d). The hexamers for (e) kagome and (f) pyrochlore lattices are also presented. Red (blue) dots denote the spin up (down). Orange (yellow) triangles for a kagome lattice denote Q = +3(+1); dark blue (light blue) tetrahedra for a pyrochlore lattice denote Q = +4(+2).

FIG. 18 .

 18 FIG. 18. The conserved spin correlator for a kagome lattice [(a)-(c)] and for a pyrochlore lattice [(d)-(f)]. Corresponding values of J and T are given in the figure.

FIG. 19 .

 19 FIG. 19. The conserved spin correlator in real space for a kagome lattice. Black lines denote bonds of the premedial honeycomb lattice. Colors of dots correspond to the sign (see the main text), and the central white circle denotes the origin. The radius is proportional to the absolute value, and rescaled as r = C|S M (R)| with C = 3 for (a) and C = 1 for (b).

FIG. 21 .

 21 FIG. 21. First order phase transition at the high-symmetry point J = J2c on pyrochlore, as demonstrated by the discontinuity of (a) the dipolar order parameters ML at wave vector qL and (b) the quadrupolar order parameter MQ. The black dashed lines indicate the value of saturation for each order parameter.

A 3 .

 3 Its explicit form is given as ĥK→H (q) =   e iϕ1 e -iϕ1 e iϕ2 e -iϕ2 e iϕ3 e -iϕ3

FIG. 24 .

 24 FIG.24. J dependence of u K . The solution of Eq. (A25) for 0 ≤ Q ≤ π exists for the shaded area. Note that u K becomes a complex number for J ≥ 1.

ACKNOWLEDGMENTS

We would like to thank Roderich Moessner for fruitful discussions, especially for the enlightening comments on the relation between pinch points and half-moons. We also thank Keisuke Tomiyasu and Jeffrey G. Rau for interesting discussions. This work was supported by the JSPS KAKENHI (Nos. JP15H05852, JP16H04026 and JP26400339), MEXT, Japan. LDCJ acknowledges support from the IdEx Bordeaux International Support and from the TQM Unit of the Okinawa Institute of Science and Technology Graduate University, as well as hospitality from Gakushuin University in Tokyo. Part of numerical calculations were carried out on the Supercomputer Center at Institute for Solid State Physics, University of Tokyo.

Note: While completing this work, we learned of a parallel study by Yan et al., which reports complementary results for a different model 51 .

Appendix A: Analytical formula for S(q) in large-N analysis

In this section, we describe the derivation of S(q) in large-N approximation, in detail. The J 1 -J 2 -J 3 model has two special properties at J 2 = J 3 ≡ J, which enables us simple analytical approach. One is the polynomial expression of Hamiltonian in terms of incident matrices, and the other is the line graph correspondence. With the help of graph-theoretical argument, we can construct a simple and systematic way to obtain the analytical expression of S(q) in Eq. (8).

= (1, 0) and a

2 , and as the coordinates of three sublattices, 1,2, and 3, r

= 1 2 , 0 . Accordingly, the coordinates of two sublattices, A and B of the dual honeycomb lattice are r (H) A = [START_REF] Onose | Observation of the Magnon Hall Effect[END_REF] 4 , 5 √ 3 12 and r

12 . For a pyrochlore lattice, the lattice vectors are a P 1 = 0, 1 2 , 1 2 , a P 2 = 1 2 , 0, 1 2 , a P 3 = 1 2 , 1 2 , 0 , and the positions of four sublattices are r P 1 = (0, 0, 0), r P 2 = 0, 1 4 , 1 4 , r P 3 = 1 4 , 0, 1 4 , r P 4 = 1 4 , 1 4 , 0 . For the dual diamond lattice, the coordinates of two sublattices (A and B) are

, and r D B = 7 8 , 7 8 , 7 8 .

Line-graph correspondence

Here, let us apply the idea of dual lattice to solve the eigenvalue problem of ĥ. Here, we focus on a kagome lattice, again. First, we look at the intermediate lattice we have introduced in the previous subsection. On this graph, we introduce N × N H rectangular matrix, ĥK←H , whose N rows correspond to sites on a kagome lattice, and N H columns correspond to the sites of a honeycomb lattice. We define ĥK←H as an incident matrix for the intermediate lattice, i.e., we set

And we define N H × N rectangular matrix, ĥH←K as ĥH←K = ( ĥK←H ) t .

The key step to solve the eigenvalue problem is the observation that the matrix ĥ can be written as a product of ĥH←K and ĥK←H : ĥ = ĥK←H ĥH←K -2 ÎN×N .

(A5)

The form (A5) immediately tells us significant information on the energy spectrum of ĥ. For N × M matrix  and M × L matrix B, it is known that

Applying these properties to ĥK←H and ĥH←K , we obtain rank( ĥK←H ĥH←K ) ≤ N H = (2/3)N.

(A6)

This inequality results in the existence of at least N -N H = 1 3 N zero modes, i.e. ĥ has 1 3 N eigenstates with degenerate eigenenergy, -2. Moreover, the matrix ĥK←HĥH←K and the inverse product, ĥH←KĥK←H share the common non-zero eigenvalues. Accordingly, given that the incident matrix of honeycomb lattice is given by

the eigenspectrum of ĥ consists of N H eigenvalues of ĥH + (N sub -2) 1, and N -N H -fold degenerate modes with eigenvalue, -2.

Momentum-space expression

The translational invariance of the Hamiltonian matrix allows us to block-diagonalize it with respect to the momentum q. For each q, we obtain N sub × N sub Hamiltonian matrix Ĥ(q) defined in Eq. (7). Due to the polynomial expression in Eq. (A2), we obtain

Therefore, Ĥ(q) is also expressed as a polynomial of the Fourier transformation of ĥ1 :

From Eq. (A10), we see the eigenvalue problem for Ĥ(q) is reduced to that for ĥ1 (q). To solve it, one can utilized the momentum space version of Eq. (A5), namely ĥ1 (q) = ĥK←H (q) ĥH←K (q) -2 ÎN sub ×N sub .

(A11)

Here ĥK←H (q) is N H sub × N sub matrix:

and ĥH←K (q) = [ ĥK←H (q)] † . (A13)

Using Eq. (A11) as well as the momentum-space analogue of Eq. (A7), one can show that the eigenvalues of ĥ1 (q) consists of eigenvalues of ĥH (q) (up to a constant) and (N sub -N H sub ) flat mode with the eigenvalue -2; ĥH (q) is given by

Note that the similar relations hold in the case of a pyrochlore lattice.

Kagome lattice

Utilizing the above idea, we now show the explicit forms of the eigenvalues and eigenvectors of our model. Let us first consider a kagome lattice. As discussed in the previous section, the exchange matrix on a kagome lattice is expressed as ĤK (q) = ĥK 1 (q) + J ĥK 2 (q), (A15) with

To obtain the eigenvalues of ĥH 1 (q), we first write down the exchange matrix on a dual honeycomb lattice:

with

The eigenvalues of ĥH 1 (q) are given by ε (H) ± (q) = ±|G(q)| and the corresponding eigenvectors are

±e iθ G (q)/2 e -iθ G (q)/2 , (A19) exchange matrix on a original pyrochlore lattice: ĤP (q) = ĥP 1 (q) + J ĥP 2 (q), (A28)

where

and

The polynomial form of the Hamiltonian (A28) with respect to ĥP 1 (q) is obtained as

since z = 6 and x = 2 for a pyrochlore lattice. Next, the exchange matrix for the dual diamond lattice is given as ± (q) = ±|g(q)|, and the corresponding eigenvectors are

±e iθg(q)/2 e -iθg(q)/2 , (A34) with θ g (q) = arg g(q). Then, using the argument in Sections A 3 and A 4, we obtain the eigenvalues of the Hamiltonian (A28) as

The corresponding eigenvectors are given as

where ĥP←D (q) is a rectangular matrix

and

, φ 2 = qx-qy-qz 8

, φ 3 = -qx+qy-qz 8

, and

. The rest of eigenvectors, i.e. two flat modes, are orthogonal to ψ P,± (q) and their eigenenergy is -2(1 -J). . The solution of Eq. (A39) for 0 ≤ Q ≤ π exists for the shaded area.

Energy minima.-The minima of ε P -(q) is obtained by solving

Again (ii) is satisfied at several high-symmetry points, which turn out not to be energy minima, so let us focus on (i). When 0 ≤ J ≤ 1 6 , (i) does not have solutions, the lowest-energy band in this region is the flat band. Then, for 1 6 ≤ J ≤ 1 2 , the solution is a surface enclosing Γ point [Fig. 7(a)]. Finally, for J ≥ 1 2 , a solution is a surface enclosing the zone corners [Fig. 7(b)].

Derivation of the phase boundaries.-Similar to the case of a kagome lattice, the phase boundaries for a pyrochlore lattice are determined by the presence/absence of the energy minima on ΓL line. On ΓL line, the momentum is parametrized as q = (Q, Q, Q) with 0 ≤ Q ≤ π. Then the condition for the energy minima is given by |g(q)| = 6 cos Q + 10, (A39) and its solution of (i) is

We plot u P ≡ -12J 2 +4J-1 48J 2

in Fig. 25. We again examine the condition that Q is between 0 and π (a shaded area of Fig. 25), and find that the lower (upper) bound is J = 1 6 1 2 .

Appendix B: Monte Carlo simulations

Monte Carlo simulations are performed on systems of classical O(3) spins on the kagome and pyrochlore lattices, whose system sizes are respectively 12L 2 and 16L 3 sites. To decorrelate the system, we use jointly the heatbath method, over-relaxation and parallel tempering. Thermalization is made in two steps: first a slow annealing from high temperature to the temperature of measurement T during t e Monte Carlo steps (MCs) followed by t e MCS at temperature T . After thermalization, measurements are done every 10 MCs during t m = 10t e MCs. All temperatures are given in units of J 1 = 1. The details of each simulation are • Fig. 3: L = 30 for both lattices and t m = 10 5 , 10 6

MCs for the pyrochlore and kagome lattice respectively.

• Fig. 8, top: t m = 10 6 MCs, and the error bars are coming from an average over n runs with different initial configurations, where n = 50 for L < 15 and n = 20 for L > 15.

• Fig. 8, bottom: L = 20 and t m = 10 6 MCs.

• Fig. 9: L = 8 (a) and L = 16 (b) and t m = 10 6 MCs. The error bars in (a) are coming from an average over 6 runs with different initial configurations; when not visible, they are smaller that the dots.

• Fig. 12: L = 30 and t m = 10 5 MCs.

• Fig. 13: L = 20 and t m = 10 5 MCs.

• Fig. 14: L = 30 and t m = 10 5 MCs.

• Fig. 21: L ∈ {6, 8, 10} and t m ∈ {10 7 , 2.10 7 , 10 7 } MCs respectively.

Appendix C: Rank-two tensor order parameter

The rank-two tensor order parameter is time-reversal invariant and measures the on-site quadrupolar order. For a pyrochlore lattice of N sites, the rank-two tensor is defined following Ref. 49,

where

The order parameter used in Fig. 21.(b) comes from the norm of all quadrupole moments

and is saturated when all spins are collinear, taking the value 2/ √ 3.