Willy Jimenez

Amel Mammar

Ana Cavalli

Software Vulnerabilities, Prevention and Detection Methods: A Review 1

Keywords: Software vulnerability, Prevention/Detection Methods, Testing

Software is a common component of the devices or systems that form part of our actual life. These systems are usually complex and are developed by different programmers. Usually programmers make mistakes in the code which could generate software vulnerabilities. A software vulnerability is a flaw or defect in the software construction that can be exploited by an attacker in order to obtain some privileges in the system. It means the vulnerability offers a possible entry point to the system. Despite the knowledge about vulnerabilities nowadays there is still a growing tendency in the number of reported vulnerabilities, reason why software security has become an important field of research. The presence of vulnerabilities in the production of software makes necessary to have tools that can help programmers to avoid or detect them in the development of the code. Thus, in relation to our on-going research on vulnerability detection, this article presents an overview of software vulnerabilities and their prevention and detection methods.

Introduction

A software vulnerability can be seen as a flaw, weakness or even an error in the system that can be exploited by an attacker in order to alter the normal behavior of the system. Because the number of software systems increases everyday also the number of vulnerabilities. Additionally, if we consider that most of the systems are exposed to multiple users (internet) and environments (operating systems for example) then it is just a matter of time that someone can launch an attack (sequence of actions) whose consequences are unpredictable in damages and cost. Usually the goal of an attacker is to gain some privileges in the system to take control of it or to obtain valuable information for its own benefit. Then it is important for the developers and general public to know about vulnerabilities and their prevention and detection.

Under the context of the European project SHIELDS, whose main objective is to bridge the gap between security experts and software developers and thereby reduce the occurrence of security vulnerabilities; we present a review of different methods to detect and prevent software vulnerabilities as well as some well known software vulnerabilities. In most of the cases vulnerabilities are caused by improper validation of the data supplied by the user. This undesired condition is used by attackers to inject faults and malicious code into the system that allows them to run their own code and applications.

For better understanding vulnerabilities the creation of models that express the set conditions that could lead or originate them is very helpful; additionally when models are well understood they could also be used for prevention. But since it is impossible to guarantee the absence of vulnerabilities in a piece of code during its creation, then it is necessary to have methods to detect them. One possibility is security software inspections, or simply manual review of the code or related documents. Also some more automated methods for vulnerability detection can be applied, which are classified into two main categories: static, when the detection is performed without running the source code; and dynamic when the program is executed in order to detect vulnerabilities. Actually, in Telecom SudParis we are doing some research about the use of models in the detection of vulnerabilities.

The organization of this paper is as follow. In section 2 we introduce some known vulnerabilities and we mention possible consequences of their exploit. Section 3 contains methods to prevent vulnerabilities. In section 4 we study vulnerability detection techniques, classified into static and dynamic according to the execution of the source code. Section 5 describes our work in progress and finally in section 6 conclusions and perspectives of this work are presented.

Software Vulnerabilities

As we have mentioned a vulnerable software system can be exploited by attackers and the system could be compromised, the attacker might take control of the system to damage it, to launch new attacks or obtain some privileged information that he can use for his own benefit. Considering this, it is important to know the different types of vulnerabilities, their prevention and detection in order to try to avoid their presence in the final software version of the system and then reduce the possibility of attacks and costly damages.

Examples of vulnerabilities

Most of the known vulnerabilities are associated to an incorrect manner of dealing with the inputs supplied by an user of the system, if these inputs are not correctly processed before using them inside the program they can generate unexpected behavior of the system. For instance, some known and frequent vulnerabilities are [START_REF] Christey | Unforgivable Vulnerabilities[END_REF]:

 Buffer overflow: it occurs usually with fixed length buffers when some data is going to be written beyond the boundaries of the current defined capacity.

This could lead to mal functioning of the system since the new data can corrupt the data of other buffers or processes. The buffer overflow can be used also to inject malicious code, and then the execution sequence of the program could be altered in order to execute the injected code and take control of the system.  XSS or cross site scripting: usually associated to web applications, consists in the injection of code in the pages accessed by other users. If exploited an attacker can bypass access controls, perform phishing, identity theft or expose connections.  SQL injection: it consists in the injection of code with the intension of exploiting the content of a database. Usually happens because the inputs are not handled correctly, the attacker can get sensitive information from the database. However, some others common vulnerabilities that can be mentioned:

 Format string bugs: it happens when external data is given to an output function as format string argument. The output function, for instance, printf in C language, generates an output according to the specifications of the format string, some directives can write to memory locations, thus the attacker can use the printf to write malicious code and change the control flow to execute it.  Integer overflows: can be of two different types, sign conversion bugs and arithmetic overflows. The first occurs when a signed integer is converted to an unsigned integer; while in the second the result of an arithmetic operation is an integer larger than the maximum integer and it is stored in an integer variable.

Vulnerability Modeling

Most of the vulnerabilities presented in the previous section could be prevented if the software is developed more carefully, avoiding the introduction of vulnerabilities that could be exploited by attackers. One solution is in the improvement of the knowledge and understanding of software developers about: known vulnerabilities, causes, threats, attacks and counter measures. Models are in fact adequate to implement such solution.

There is for instance a vulnerability model called Vulnerability Cause Graph (VCG) [START_REF] Byers | Modeling Software Vulnerabilities with Vulnerability Cause Graphs[END_REF][START_REF] Ardi | Towards a structured unified process for software security[END_REF] which "is a directed acyclic graph that contains one exit node representing the vulnerability being modeled, and any number of cause nodes, each of which represents a condition or event during software development that might contribute to the presence of the modeled vulnerability". An example of a VCG representing a known buffer overflow in xpdf (CVE-2005-3192) taken from [START_REF] Byers | Modeling Software Vulnerabilities with Vulnerability Cause Graphs[END_REF], is shown in figure 1. In this graph we can observe the different causes and possible scenarios or sequence actions that could lead to the introduction of this kind of vulnerability. The VCG is helpful to understand what can cause the vulnerability. If causes are well understood then they could be avoided in the development process.

CVE-2005-3192

Unchecked integer multiplication

External data influences buffer size

Use of nonadaptive buffers

The return value of malloc is not checked Use of malloc/calloc/ realloc/alloca

(1) (2) (3) (4) (5) (6)
Failed to check input parameter to malloc [START_REF] Wheeler | Flawfinder[END_REF] Fig. 1. Vulnerability Cause Graph

Preventing Software Vulnerabilities

Models are a first approach to deal with vulnerabilities and their understanding. However it is necessary to count on methods or procedures to prevent any risks related to vulnerabilities. In this section two possible vulnerability prevention methods developed in the literature are presented. The purpose is to evaluate the code during the construction process to detect any security defect and correct it on time without the need of performing intensive test at the end when the whole program is finished.

Software Inspection

The software inspection process consists in reading or visually inspecting the program code or documents in order to find any defects and correct them early in the development process. When the defect is found soon the less expensive it becomes to fix. However, a good inspection depends then on the ability and expertise of the inspector, and the kind of defects he is looking for. Usually during the software inspection, it is necessary to look for any possible defects during the security inspections. In the following sections we introduce two inspection methods that intend to codify the implicit knowledge of security experts regarding how to check for correct implementation of security goals and how to search for vulnerabilities.

Security Goal Indicator Trees

Security Goal Indicator Trees (SGIT) [START_REF] Peine | Security Goal Indicator Trees: A Model of Software Features that Supports Efficient Security Inspection[END_REF] focus on positive features of the software which can be verified during the inspection process. A SGIT is then a graph where the root is a security goal and its subtree are indicators or properties that can be checked for achieving that goal. However, since not all properties can be positively expressed it is possible to have also negative indicators (something that should not occur). These indicators have Boolean relations with the goal and have to be checked in order to validate the security goal. SGIT are created by security experts. A SGIT for the goal Audit Data Generation, taken from [START_REF] Peine | Security Goal Indicator Trees: A Model of Software Features that Supports Efficient Security Inspection[END_REF], is presented in figure 2, showing some dependency relations, and positive and negative indicators. Also the small box pointing to the indicator "An audit component exists" means that a specialization tree can be deployed for this indicator.

Vulnerability Inspection Diagram

Vulnerability Inspection Diagram (VID) is also a manual inspection introduced in [START_REF]D2.1 Formalism definitions and representation schemata[END_REF], the purpose is to benefit developers from the knowledge and experience of security experts in the detection of problems in the development process. Thus a VID is a flowchart-like graph that guides developers to check the software to detect the presence of vulnerabilities based on the knowledge of experts. There is a specific VID for each vulnerability class.

Security Activity Graph

Security Activity Graphs (SAGs) [START_REF] Ardi | Towards a structured unified process for software security[END_REF][START_REF] Byers | A Cause-Based Approach to Preventing Software Vulnerabilities[END_REF] are also helpful in the prevention of vulnerabilities. SAGs are a graphical representation that is associated with causes in a VCG. SAGs indicate how a particular cause can be prevented following a combination of security activities during the development process. To illustrate this, in figure 3 there is a SAG [START_REF] Byers | A Cause-Based Approach to Preventing Software Vulnerabilities[END_REF] showing different alternatives to address the cause "Lacking design to implementation traceability". Models and inspections are useful to understand and prevent vulnerabilities; nevertheless it is also necessary to count on tools that can be used by programmers in order to detect vulnerabilities during the process of software construction. Some of these tools are based on static methods, thus it is not necessary to run the code to perform the detection. In the case of dynamic methods, the code is run inside a controlled environment to perform the detection or collect program traces that can be use for such purpose. In the next section we present some static and dynamic techniques to detect vulnerabilities.

Static Techniques

Static techniques are those applied directly to the source code without running the application, the objective is to evaluate or get specific information directly from the source code without executing it. There are different techniques to perform static analysis; here we mention some of them.

Pattern Matching

Consists in searching a "pattern" string inside the source code and give as results the number of occurrences of it. For instance if we consider C language, the pattern could be any call to possible dangerous functions (vulnerable) like "getc". Pattern matching can be implemented using a simple tool like the Unix command "grep", however this method generates much false positives because there is no analysis of the results, additionally its effectively is limited since depends on the exact writing of the strings, thus additional white spaces will limit the results.

Flawfinder also uses a more elaborated pattern matching process to find possible vulnerabilities and sort them by risk level [START_REF] Wheeler | Flawfinder[END_REF]. This risk level depends on the function and on the values of the parameters of the function.

Lexical Analysis

Lexical analysis adds an additional step before applying a pattern match. In fact, the source code is transformed into a sequence of tokens, which are later compared with a vulnerability database in order to identify them. False positives number is still high because they do not consider the syntax or grammar of the program. The ITS4 [START_REF]ITS4: Software Security Tool[END_REF] tools use lexical analysis.

Parsing

Parsing is more complex than lexical analysis, thus when the source code is parsed, a representation of the program is built using a parsing tree in order to analyze the syntax and the semantics of the program. For example the parsing technique is used to detect SQL command injection attacks [START_REF] Su | The Essence of Command Injection Attacks in Web Applications[END_REF].

Type Qualifier

Type qualifiers are used to qualify types and modify the properties of variables in the programming language, as in the case of the tool Cqual [START_REF] Cqual | A tool for adding type qualifiers to C[END_REF]. Cqual is used to specify and check properties of C programs using user-defined type qualifiers, which are added to the program. The modified program is analyzed to find vulnerabilities.

Data Flow Analysis

The purpose is to determine the possible values a variable or an expression can have during the execution of the program, specially suited for buffer overflow detection. For instance data flow analysis is used in [START_REF] Kim | Applying Dataflow Analysis to Detecting Software Vulnerabilities[END_REF]. The authors take rules describing vulnerability patterns and the source code to detect locations and paths of the pattern in the program. The process is executed in three parts: pattern matching, control and data flow and flow analyzer.

Taint Analysis

It is a special case of data flow analysis where any data coming from un-trusted sources, e.g. introduced by a user, is a potential problem to the system, thus it is marked as tainted. Tainted data flow is monitored because it can not reach critical functions unless it is processed and changed to untainted.

Livshits and Lam [START_REF] Livshits | Finding Security Vulnerabilities in Java Applications with Static Analysis[END_REF] propose a static analysis framework to find vulnerabilities in Java applications. They define a Tainted Object Propagation problem class to deal with improper user input validation. Java bytecode and vulnerability specifications are employed to perform a taint object propagation and find vulnerabilities using the Eclipse platform.

Model Checking

Model Checking is a technique to automatically test if the model of a system meets its specification and it can be used to detect vulnerabilities. Usually model checking is a complex technique because the elaboration of the model is difficult, however once obtained it is easier to test the properties of the system.

A security verification framework with multi-language support was developed [START_REF] Hadjidj | Model-Checking for Software Vulnerabilities Detection with Multi-Language Support[END_REF] based on GCC compiler. Their approach uses a conventional push down system model checker for reach ability properties to verify software security properties; it is composed of three phases: security property specifications, program model extraction and property model checking, this last has as output the detected errors with execution traces.

Constraint analysis is combined with model checking [START_REF] Wang | Automated Detection of Code Vulnerabilities Based on Program Analysis and Model Checking[END_REF] in order to detect buffer overflow vulnerabilities. They trace the memory size of buffer-related variables and the code instrumented with constrains assertions before the potential vulnerable points. The vulnerability can be detected with the reach ability of the assertion using model checking. They decrease the cost of model checking slicing the program.

Dynamic Techniques

In order to dynamically detect vulnerabilities it is necessary to execute the program code, and then analyze the behavior or the answers of the system and gives a verdict. In the next part we study some of the techniques to perform dynamic detection.

Fault Injection

Fault injection is a testing technique that introduces faults in order to test the behavior of the system, some knowledge about the system is required to generate the possible faults. With fault injection is possible to find security flaws in the system [START_REF] Du | Vulnerability Testing of Software System Using Fault Injection[END_REF], in this work faults are injected into the system under test and the system behavior is observed, the failure to tolerate faults is an indicator of a potential security flaw in the system, a model is used to decide what faults to inject.

Fuzzing Testing

The idea of this test is to provide random data as input to the application in order to determine if the application can handle it correctly. Fuzzing testing is easier to implement than fault injection because the test design is simpler and previous knowledge about the system to test is not always required, additionally it is limited to the entry points of the program. Web scanners are in this tool category.

Fuzzing testing can also be improved to have a better coverage of the system. For instance recording real user inputs to fill out web forms and then utilize the collected data in the fuzz testing process to better explore web applications (reach ability) [START_REF] Mcallister | Expanding Human Interactions for In-Depth Testing of Web Applications[END_REF].

Dynamic Taint

Similar to taint analysis, however in this case the tainted data is monitored during the execution of the program to determine its proper validation before entering sensitive functions. It enables the discovering of possible input validation problems which are reported as vulnerabilities [START_REF] Chess | Dynamic Taint Propagation: Finding Vulnerabilities without Attacking[END_REF].

Sanitization

One possibility to avoid vulnerabilities due to the use of user supply data is the implementation of new incorporated functions or custom routines whose main idea is to validate or sanitize any input from the users before using it inside a program. In [START_REF] Balzarotti | Saner: Composing Static and Dynamic Analysis to Validate Sanitization in Web Applications[END_REF] they present an approach using static and dynamic analysis to detect the correctness of sanitization process in web applications that could be bypass by an attacker. They use data flow techniques to identify the flows of input values from sources to sensitive sinks or the places where the value is used. Later they apply the dynamic analysis to determine the correct sanitization process.

Our approach

Our research in Telecom SudParis evolves around the use of models for tool-based detection of vulnerabilities. In our approach a vulnerability model, the Vulnerability Cause Graph are considered as an input in order to derive a formalism called Vulnerability Detection Condition [START_REF] Mammar | Modélisation et Détection Formelles de Vulnérabilités Logicielles par le Test Passif. 4ème conférence sur la Sécurité des Architectures réseaux et des Systèmes d'information[END_REF], with the goal of automatically test the source code to detect vulnerabilities. The main idea behind this concept is to use the information provided by VCGs to point out in the code the use of a dangerous action under some particular conditions, for instance "it is dangerous to use unallocated memory". The checking for vulnerabilities is performed on execution traces of the program using the TestInv tool [START_REF] Cavalli | Two Complementary Tools for the Formal Testing of Distributed Systems with Time Constraints[END_REF]. Another tool, TestGen [START_REF] Cavalli | Two Complementary Tools for the Formal Testing of Distributed Systems with Time Constraints[END_REF], might be adapted to generate test cases for the detection of vulnerabilities.

Another possibility considered in our research is the extension of Martins et al work [START_REF] Martins | Generating Attack Scenarios for the Validation of Security Protocol Implementations[END_REF] which uses attacks trees to generate test cases to uncover protocol vulnerabilities. The objective is to extend this work toward more general programs. Additionally we will evaluate a possible integration with the modeling tool Seamonster [START_REF] Meland | SeaMonster: Providing tool support for security modeling[END_REF], it uses attack trees, in order to have a more powerful tool.

As we can see vulnerabilities are not a new topic on software field; however it is also noted that they still appear in the source code, thus it means that programmers still do not know how to deal with vulnerabilities. In order to help programmers to build better code we could use vulnerability cause graphs to teach them how the vulnerabilities are introduced thus they could learn to avoid the presence of vulnerabilities in their source code. However, in the mean time the source code should be inspected to guarantee there are no vulnerabilities, this method can be applied several times during the construction phase as advantage but requires specialists to perform the task as drawback.

Also a number of tools are available in order to detect vulnerabilities, some of them are based on static techniques, and it means the source code is analyzed without running the application while on dynamic techniques it is necessary to run it. The selection of the tools is related to the type of application to evaluate, the programming language and the type of vulnerability to detect. The static techniques cover all possible execution paths but require the source code while dynamic techniques have the difficulty of requiring the preparation of test cases and the possibility that not all paths in the program are covered, but the advantage that the problems if any, are found in the running code. Dynamic techniques have also less false positives than statics.

Finally, our current research intends to create new vulnerability detection methods based on models. In this manner we could guarantee the reusability of the tests cases and facilitate the transformation of these formal representations into the specific programming language of the tool used to perform the vulnerability detection.

Fig. 2 .

 2 Fig. 2. Security Goal Indicator Tree

Fig. 3 .

 3 Fig. 3. Security Activity Graph

The research leading to these results has received funding from the European Community's Seventh Framework Program (FP 7/2007-203) under the grant agreement number

(http://www.shields-project.eu/)


Mathieu Taillefumier

Julien Robert

Christopher L Henley

Roderich Moessner

Benjamin Canals

Semi-classical spin dynamics of the antiferromagnetic Heisenberg model on the kagome lattice

We investigate the dynamical properties of the classical antiferromagnetic Heisenberg model on the kagome lattice using a combination of Monte Carlo and molecular dynamics simulations. We find that frustration induces a distribution of timescales in the cooperative paramagnetic regime (i.e. far above the onset of coplanarity), as recently reported experimentally in deuterium jarosite. At lower temperature, when the coplanar correlations are well established, we show that the weathervane loop fluctuations control the system relaxation : the time distribution observed at higher temperatures splits into two distinct timescales associated with fluctuations in the plane and out of the plane of coplanarity. The temperature and wave vector dependences of these two components are qualitatively consistent with loops diffusing in the entropically dominated free energy landscape. Numerical results are discussed and compared with the O(N ) model and recent experiments for both classical and quantum realizations of the kagome magnets.

I. INTRODUCTION

In psychology, frustration is an emotional response to opposition or conflict. In the natural sciences, frustration is often associated with the impossibility of a system to optimize simultaneously all elementary interactions, whether they are single body, two body, or many body. It gives rise to many exciting phenomena in particular in magnets where magnetic states stay disordered despite the presence of strong interactions.

Although the first studies of spins models with competing interactions date back from the early 50s 1,2 , the terminology of frustration was introduced in the 70s mostly in the context of spin glasses [3][START_REF]See article of G. Toulouse[END_REF][START_REF] Villain | [END_REF][6] . Since then it has been associated to many unconventional low energy states such as quantum and classical spin liquids 7,8,[10][11][START_REF] Misguich | Quantum spin liquids and fractionalization[END_REF][START_REF] Hermele | [END_REF] , or more recently, to quantum, classical and artificial spin ices [START_REF] Gingras | Spin Ices, Chapter 3, Introduction to Frustrated Magnetism[END_REF][16][17][18] .

All these phenomena take place in diverse systems but each of them can be associated with a canonical representative, i.e. a minimal frustrated spin model which brings together most of the important features.

The 2D antiferromagnetic Heisenberg model on the kagome lattice (KHAFM) is one of the archetypes of such systems. The kagome lattice can be described as a lattice of triangles sharing one corner with each neighbor, the key property for creating these unusual and often highly degenerate ground states both in classical and quantum models 8,[17][18][19][20][21] .

In the quantum limit, the ground state of KHAFM is still unsettled, but recent results rather point towards a quantum disordered ground state [22][23][24][25] . In the classical limit, the equilibrium properties of KHAFM in the paramagnetic (T > |J|) and cooperative regimes (0.001|J| < T < |J|) regimes are now well known and understood down to low temperature 8,9,14,19,20,[26][27][28] (where |J| is the nearest neighbor interaction coupling constant). Finally a recent study finds a very weak magnetic order when the system is deep in the coplanar regime (T < 0.001|J|) 29 .

The dynamics of the classical Heisenberg model in both one, two or three dimensions is very rich at all temperatures because of the non linearity of the model. The dynamics in the paramagnetic regime was extensively studied during the nineties.

Early studies of the dynamics of frustrated magnets have shown that these systems are very different from their non-frustrated counterpart 30 at low temperature. In classical 3D frustrated magnets, such as pyrochlore antiferromagnets, most correlations in magnetic states decay exponentially at low temperature while the temperature dependence of the relaxation time follows a power law 7,31,32 . Similar results were found in the kagome antiferromagnets despite the more complex landscape around the ground state manifold 33,34 .

At low temperature, the absence of long range ordered ground states does not forbid short-lived spinwave like excitations 35 whose natural time scale is of the order of |J| -1 . The highly degenerate nature of the ground state manifold gives rise to additional processes that contribute at intermediate time scales 7,33,34 . The longest time scales (t|J| > 500) are the domain of both (i) spin diffusion, i.e. the stochastic propagation of the magnetization throughout large magnetic regions, and of (ii) magnetic relaxation, i.e. the gradual reorganization of the average spin configuration around which the spin waves are oscillating. All these time scales are present in KHAFM and can be studied with different experimental probes, ranging from thermal neutron scattering and muon relaxation for the short and intermediate time scales to ac-susceptibility and NMR measurements for the long time scales.

The aim of this work is to understand and characterize the dynamics beyond short time scales. We find that spin diffusion persists at low temperature (T /|J| < 10 -2 ) despite the presence of strong spin correlations. Below T /|J| = 10 -2 , the spin dynamics becomes anisotropic due to the entropic selection of coplanar states. We also find that the relaxation is mediated by large amplitude oscillations around small loops (also called weathervanes defects) and spinwaves, despite the absence of long range order.

This article is organized as follows. In section II, current knowledge of KHAFM thermal behavior is reviewed in order to provide the reader with a clear description of the magnetic structures that the spin dynamics is built on. Then, numerical procedures and technical details are given in section III. Finally, our numerical results are presented in sections IV and V respectively for the high and low temperature regimes, and compared with recent experimental results obtained in kagome systems (see Sec. VI).

II. EQUILIBRIUM PROPERTIES OF THE CLASSICAL KAGOME HEISENBERG ANTIFERROMAGNET -FOUNDATIONS FOR A SPIN DYNAMICS

The phase diagram of the classical KHAFM is depicted on Fig. 1. At high temperatures (T > |J|), the classical KHAFM is a conventional paramagnet with short range spin-spin correlations. When temperature becomes comparable to the exchange |J|, correlations appear and spins on each triangular plaquette of the kagome lattice approximately sum to zero and are oriented at 120°to one another. This local arrangement does not lead to large correlated domains because of the excitence of an uncountable number of configurations that form a highly degenerate and connected manifold associated with "origami" folding of the spin pattern [36][37][38] . As a result of this degeneracy, spin correlations decay algebraically with distance and can be associated with a socalled Coulomb phase 7,19,27,[39][40][41] . In such a phase, correlations are expected to decay algebraically with distance, with geometrical factors that depend on the chosen direction in the kagome lattice. This regime roughly covers the temperature range 5.10 -3 < T /J < 1.

When temperature is further reduced, i.e T /J < 5.10 -3 , the free energy of all spin foldings is no longer uniform and the spins, which are still locally constrained to stay at 120°within each triangle, now select a particular spin plane 74 , common over many triangles, around which they are fluctuating 20 .

This selection of coplanar states, also known as entropic ordering (or order out of disorder), is due to the additional soft degrees of freedom for the thermal fluctuations 7,20,36,42 available in the coplanar states. This was first identified as a coplanar ordering, i.e. the devel-opment of quadrupolar (or spin nematic) correlations 20 . This incipient order is not merely coplanar but was later recognized to imply octupolar order as well 26 .

Thus, in the model's ultra low temperature regime, spins fluctuate around one of the discrete coplanar ground states, in which every spin has one of three possible directions, which can be represented by the values (or colors) of the discrete spins in the 3-state Potts model on the same lattice. The coplanar ground states correspond 1-to-1 up to global rotations to Potts ground states, in which every triangle has three colors 19,27,38 , whose number is N c ≈ 1.13 N where N is the number of spins of the lattice (we will always consider finite lattices with periodic boundary conditions).

Consequently, there are essentially three different regimes. The generic paramagnetic regime with short range spin correlations, a cooperative paramagnetic regime or spin liquid regime, with algebraic correlations on finite area domains, whose area is controlled by a temperature dependent correlation length and a nematic-like regime, where correlations are enhanced via an order out of disorder phenomenon that stabilizes a common spin plane. At very low temperature magnetic ordering also appears [29].

The first studies of the dynamics of magnetic systems concentrated on the nature of spin fluctuations in the cooperative paramagnetic regime in comparing the spin dynamics of a strongly correlated disordered magnet with the dynamics of an ordered one. Following this perspective, it was shown that the KHAFM is a model with unusually high density of low lying excitation 30 at low temperatures. It was also shown that at sufficiently low temperatures (T /J ≤ 5.10 -3 ), coherent excitations are unexpectedly stable despite being built on a thermodynamically characterized disordered manifold 33 .

In this work, our interest is to understand how the natural high temperature (T J) signature of diffusive dynamics is found at lower temperatures (T J), how it terminates close to the nematic boundary even though local excitation are still present (6-sites loops), and how this model discriminates between in-plane and out-of-plane spin dynamics, all latter considerations being discussed in the intermediate 10 J -1 ≤ τ ≤ 10 4 J -1 time scale. In other words, it aims at resolving spin fluctuations whilst extending previous dynamical studies in order to cover a time range associated with magnetic relaxation rather than structured and propagative excitation.

10 -4 10 -3 10 -2 10 -1 T coplanar regime

• octupolar and nematic order

• algebraic spin correlations

• both localized and delocalized excitations

• Anisotropic dynam- ics.
spin liquid regime The paramagnetic regime has conventional diffusive behaviors in the high T regime. The first crossover occurs around T ≈ J and the system enters a regime where spin correlations develop at a length scale of the order of the spin correlation length ξ(T ). Around T ≈ 10 -3 J, the entropic selection favors coplanar states whose dynamics is anisotropic.

• s1 + s2 + s3 = 0 • dipolar-like S(0) • S(r) ∝ 1/r
tic sampling of the phase space and to integrate the non linear equations of motions are then detailed. Based on this technical framework, we justify our choice of temperature range and lattice sizes to ensure that most of the discussed results are free of finite size effects. We end this section with a short overview of the dynamics in the three temperature regimes that will be developed in the following sections.

A. Model

We consider the classical Heisenberg model

H = -J i,j s i • s j , (1)
where the summation is limited to nearest neighbors, J < 0 is the isotropic antiferromagnetic coupling constant and |s i | = 1 are classical spins on the unit sphere S 2 located at the kagome sites. The kagome lattice is described as a two dimensional triangular lattice with a triangular unit cell and displacement vectors a = a(1, 0) and b = a(-1/2, √ 3/2), with a the lattice constant. The unit cell contains three spins at positions r 1 = (0, 0), r 2 = a/2 and r 3 = b/2. The index i = (R i , α i ) in Eq. ( 1) is a compact notation that regroups both the position R i of the unit cell where the spin resides and α i its sublattice index. With these notations, the Brillouin Zone (BZ) is an hexagon with cor-

ners located at (Q a , Q b ) = ±(1/3, 1/3), ±(2/3, -1/3), ±(1/3, -2/3) in reciprocal space with (Q a , Q b ) = Q a a + Q b b .
It is convenient to express Eq. ( 1) as

H = - J 2 η l 2 η + E 0 , (2)
where E 0 is a constant energy shift and l η = i∈η s i is the total spin of triangle η. From this expression, it is possible to see that the ground state satisfies l η = 0 for all triangles, thus leading to a relative angle of ±2π/3 between neighboring spins in any ground state. In this article, our interest lies in the time evolution of the spin-pair correlations emerging in such a model. It is convenient to probe such dynamical correlations in reciprocal space by calculating the scattering function, also called dynamical structure factor

S(Q, t) = s -Q (0) • s Q (t) , (3)
with

s Q (t) = i,α s i,α (t)e -i(Ri+rα)•Q (4)
R i and r α are respectively the position of the unit cell and the coordinates of sublattice α.

In expression (3), the semi classical spin dynamics at T = 0 is described by the non-linear Bloch equations

ds i (t) dt = -Js i (t) ×   j s j (t)   , (5)
where sites j are the nearest neighbors of i. Note that the set described by Eq. ( 5) conserves the total energy E tot and magnetization M tot . we numerically integrate Eq. ( 5) in order to evaluate the scattering function. We combine the deterministic integration of the equations of motion with (hybrid) Monte Carlo simulations for generating samples of spin arrays at a given temperature. This numerical procedure is detailed in the next section.

B. Numerical procedures

The numerical integration of Eq. ( 5) has been performed up to 1024 J -1 (even up to 10 4 J -1 in the coplanar regime, see section V A) using an 8th-order explicit Runge-Kutta (RK) method with an adaptative step-size control, offering an excellent compromise between accuracy and computation time. The RK error parameter as well as the RK order have been fixed in order to preserve the Euclidean distance d = [ i (s RK i -s BS i ) 2 ] 1/2 , i.e. the distance between time trajectories obtained with the RK method and with the more robust but time consuming Burlisch-Stoer (BS) algorithm [START_REF] Hairer | Solving Ordinary differential equations I and II[END_REF] . As a result, trivial constants of motion, such as the total energy E tot and magnetization M tot , are conserved with a relative error smaller than 10 -9 .

The initial spin configurations used for the numerical integration are generated at each temperature by an hybrid Monte Carlo method using a single spin-flip Metropolis algorithm. The single spin flip algorithm becomes inefficient at low temperature because the number of rejected attempts increases due to the development of spin correlation as the system enters the liquid and the spin nematic regimes. To partially overcome this effect, we reduce the solid angle from which each spin flip trial is taken to ensure that the acceptance rate is above 0.4 at every temperature.

Thousand spin configurations are used at each temperature to evaluate the ensemble average in Eq. ( 3) while the number of Monte Carlo steps needed for decorrelation is adapted in such a way that the stochastic correlation between spin configurations is lower than 0.1.

In the coplanar regime, the stochastic correlation between spin configurations is relatively high because the number of accepted attempts is small. The system is trapped in the immediate surroundings of one given coplanar configuration which means that ensemble averaging is only representative of the initial conditions. To limit this effect, we use an hybrid Monte Carlo Metropolis algorithm that combines both over-relaxation [START_REF] Creutz | [END_REF] and the molecular dynamics method described earlier.

These two methods correspond to rather different ways of exploring the configuration space. An over-relaxation move, which fulfills the detailed balance, consists of rotations the selected spin by a random angle around its local exchange field so the system does not remain precisely in the same spin configuration when the single spin flip is rejected. However it does not prevent the system from being trapped into the immediate surroundings of one given coplanar configuration, so a huge number of Monte Carlo steps are still necessary for the system to decorrelate.

On the other hand, as shown in this paper, the molecular dynamics procedure is a very efficient way to probe different coloring states (or Potts states) related to each other by a spatially localized excitation. Indeed, twocolor closed spin loops of small size are easily flipped while integrating the equation of motions, even at tem-peratures as low as T /J = 0.0001. Thus, our method acts as a "natural" loop algorithm although the method is limited to small loops as the flipping time grows rapidly with loop size and temperature.

The numerical results have been obtained for different lattice sizes ranging from L = 144 (for the Q-resolved scattering function S(Q, t)) to L = 192 (for the autocorrelation function A(t)) with periodic boundary conditions, so the total number of spins does not exceed N = 3L 2 1.2 10 5 .

Finite size effects, which are negligible at high temperature, become particularly important at low temperature. Fig. 2-(a) shows that the evolution of the autocorrelation function

A(t) = d 2 Q S(Q, t) (6)
= i S i (0) • S i (t) (7)
at T /J = 0.0006 for different lattice sizes from L = 36 (red) to L = 144 (blue) becomes nearly independent of the system size when L > 108. Moreover the long time dynamics is affected by the rotation of the spins around the residual magnetization [7] 75 . The fit of the autocorrelation is represented for each lattice size on Fig. 2 (b) (for the fitting process, see section V C and Eq. ( 21)). While the short time relaxation (tJ < 60) does not depend on the lattice size, the long time relaxation, plotted versus L in Fig. 2 (c), does not seem to vary significantly for L > 120. Consequently, finite size effects will be neglected in the following for T /J 0.0001 and L ≥ 144. Microscopic quantities like the correlation lengths or more generally the Q-resolved scattering function S(Q, t) may however be affected by finite size effect at least in certain regions of reciprocal space. This problem will be discussed in section V.

C. Three dynamical regimes with blurred boundaries

Our main goal in this work is to probe the fluctuations around the ground state manifold. Before characterizing the relaxation dynamics and establishing in particular the temperature and wave-vector dependence of the lifetime of the correlated magnetic domains, we qualitatively discuss the temperature dependence of a global quantity, the autocorrelation function A(t) defined by Eq. (7). A more detailed study of each regime is given in sections IV and V.

The autocorrelation function A(t), shown in Fig. 3 (a)-(c), is respectively represented by dashed-red, solidgreen and dotted-dashed-blue lines for the paramagnetic, spin liquid and coplanar regimes. Linear-log (a,c) and log-log (b) scales are used to exhibit both exponential relaxations and diffusive behaviors. While it is manifest from panel (a) that a slowing down of the spin fluctuations with decreasing T is at work -as could be expected for any conventional magnetic system -we also notice that the overall shape of A(t) strongly depends on the temperature range.

In the paramagnetic regime, although A(t) = 1 -t 2 α for shortest times, the linear tail in a log-log scale above tJ 4 (see Fig. 3 (b)) is characteristic of spin diffusion expected to be found in the limit of high temperatures and long wave-lengths and times 47,48 . The signature of diffusive behavior is strongly reduced with decreasing temperature and is no longer visible in the two lowest temperature regimes. Nevertheless, it will be shown in section IV B that (i) spin diffusion is still present in a slightly reduced q-range with the onset of short-range correlations below T /J = 0.1, and (ii) this range tends to zero at the octupolar transition (at least, it becomes smaller than the wave-vector resolution, so that there is no apparent diffusive behavior in our simulations for this lattice size).

Below the paramagnetic/cooperative crossover occurring around T /J ∼ 0.1, the rough linear dependence of A(t) in a linear-log scale (see Fig. 3(a)) suggests an exponential decay e -t/τ T with a temperature dependent relaxation time τ T . Nonetheless, the detailed analysis of S(Q, t) given in the next section will highlight that τ T is Q-dependent as well, so that only an average appears in

A(t).
Finally, the most surprising feature in Fig. 3 is probably the intriguing behavior of A(t) in the octupolar regime, showing a kink in the A(t) behavior at around tJ ∼ 60. It is related to the presence of two relaxation processes that are different in nature (see sec. V for more details). In particular, it will be shown that the entropic selection (i) strongly affects the fluctuations of the groundstate manifold far above the transition toward coplanarity (T /J 0.05, see section IV B), and (ii) leads to different dynamical behavior for the in-plane and out-of-plane spin components below the transition (see sec.V).

IV. PARAMAGNETIC AND COOPERATIVE REGIMES

A. Models, predictions

In the absence of any order, the most basic dynamical process that may happen in a simple Heisenberg spin model is a stochastic process transferring spin density from a magnetic site to a neighboring one. By a succession of such thermally-activated random steps, the spin density arrives at a large distance r with a probability given by phenomenological spin-diffusion theory [45][46][47] . Since the total magnetization is a conserved quantity, the magnetization density m(r, t) must fulfill a local equation of continuity

∂m(r, t) ∂t + ∇.j(r, t) = 0. ( 8
)
If we assume that the local current j(r, t) is related to the magnetization by Fick's first law

j(r, t) = -D∇m(r, t), (9)
where D is the diffusion coefficient that depends on the details of the model, and after expressing Eq. ( 8) and ( 9) in Fourier space, the magnetization density obeys the diffusion equation

∂m(q, t) ∂t = -Dq 2 m(q, t) (10)
in the hydrodynamic regime, i.e. for large time t and wave vectors smaller than the inverse of the correlation length q < 1/ξ 46 . After integration over time of Eq.10,

S(q, t) = S(q, 0)e -Dq 2 t . (11)
Integrating over q gives rise to an autocorrelation with a tail that follows a power law A(t) t -d/2 where d is the dimension of the system 47 . At lower temperature T J, the spin dynamics becomes sensitive to the magnetic correlations which extend over scales of the order of the spin correlation length ξ, which diverges as 1/T according to the predictions for Ncomponent spin model 49 (ICSM). This model describes very well the apparition of structured spin pair correlations in classical Heisenberg systems 40,50 . Coupled to an appropriate Langevin dynamics, it becomes a powerful method to predict the temperature dependence of the dynamical properties 32 .

In this model 32 which we describe here for completeness, each spin component in the large-N limit has the normalized probability distribution e -βE with

βE = 1 2 i λs 2 i + 1 2 βJ α l 2 α . (12)
l α = i∈α s i is the sum of the soft spins -∞ < s i < ∞ forming the triangle α. The energy function (12) differs from Eq. ( 5) by an additional term that constrains the length of the spins. The Lagrange multiplier λ = 1 + O(T /J) in the limit T J is obtained by imposing s 2 i = 1/3 to each single component of the spin to mimic the behavior of Heisenberg spins 32 .

Then, the diffusive dynamics emerging from these static correlations can be described by the Langevin equation

ds i dt = Γ l ∆ il ∂E ∂s l + ξ i (t) (13)
for each spin component, and whose integration yields an analytic expression of the dynamical scattering function S(q, t). In this expression, ∆ ij = A ad ij -zδ ij is the lattice Laplacian, z is the coordination number of the lattice (z = 4 for the kagome lattice) and A ad ij is the adjacency matrix (see the appendix for details). Γ, which sets the energy scale of the dynamical processes, is the only free parameter of the model. This model contains two terms, a drift current that we take proportional to the difference of generalized forces ∂E/∂s j on each bond of the lattice, and a second term imposing thermal equilibrium described by a Gaussian noise contribution ξ i (t) on each site i of the lattice bonds. The noise term is correlated with an amplitude fixed by the requirement of thermal equilibrium:

ξ i (t)ξ j (t ) = 2T Γ∆ ij δ(t -t ). ( 14
)
This model was initially proposed for studying the dynamics of the pyrochlore antiferromagnet in the limit T J 32 . Similar results are found for the KHAFM: around the center of the Brillouin zone where the scattering function is described by Eq. ( 11) with a temperature independent diffusion coefficient. At larger wavevectors and away from the high symmetry directions, the decay rate τ -1 α ∝ T varies linearly with temperature and does not depend on q.

B. Numerical results and discussion

In the following subsections, we show that spin diffusion is observed in the hydrodynamic regime as predicted in the previous section with, however, a temperature dependent diffusion coefficient D T . At larger wave-vectors and away from the nodal lines [h, 0], [0, h] and [h, -h], where the dynamical properties are dominated by finiteenergy spin wave-like excitation 33,34 , an exponential relaxation is observed with a temperature and wave-vector dependent relaxation time τ T (Q) (sec. IV B 2), revealing the sizable effect of the entropic selection even at temperatures far above the transition toward coplanarity.

Hydrodynamic regime

In the hydrodynamic regime, expected to characterize only the long-wavelength-low-frequency response of the system (see sec.IV A), the scattering function S(Q, t) is expected to decrease exponentially with a relaxation rate τ -1 T (q) = D T q 2 where D T is the diffusion coefficient (see eq.( 11)). On the other hand, a short time expansion of

S(Q, t) S(Q, 0) = 1 - 1 2 ω 2 t 2 + 1 24 ω 4 t 4 + O(t 6
), (15) with

ω n = ∞ -∞ ω n S(Q, ω)dω/ ∞ -∞ S(Q, ω
)dω the nth moment of the normalized spectral weight function. As for the 1D case 45 , we find that the expansion up to the fourth order describes well the numerical simulations for 0 < tJ < 1. Above the spin velocity correlation time t vc = ω 2 -1/2 which is no more than a few J -1 at most wavevectors, the hydrodynamic regime appears and diffusion occurs.

We proceed as follows to extract D T . The relaxation time τ T (q) is obtained by fitting the scattering function S(Q, t) at some given temperature and wave-vectors using Eq. ( 11) for times t t vc : an example of such a fit is presented on Fig. 4 (a) for T /J = 0.17 at wave-vectors taken along the [h, h] direction around q = 0 (in this model the diffusive behavior is isotropic in Q-space, so other directions are not represented in the figure). Then, fitting the relaxation time τ T (q) versus q 2 for each temperature gives both the range of validity in q of the diffusive behavior and the diffusion coefficient D T (see Fig. 4 (b)).

The temperature dependence of D T is represented in Fig. 5 in both paramagnetic and cooperative regimes. At high temperature, D T asymptotically tends to a constant value D ∞ = 0.125(5) JSa 2 . This value is a little higher than the prediction 45

D T = π 2 √ 3 lim q→0 ω 2 q 2 ω 2 ω 4 1/2 (16)
which is obtained by considering a Lorentzian response for S(Q, ω) truncated at frequencies ωt vc > 1 to take into account the failure of the exponential approximation at times shorter than t vc (note that the coefficient

π 2 √
3 becomes π/2 if we consider a short time expansion instead of a rough cut-off, although the global expression remains identical 51,52 ). In the infinite temperature limit, expression (16) gives D ∞ 3/16r 2 0 (J/ ) zS(S + 1), with r 0 = a/2 the distance between two nearest neighbors and z = 4 the connectivity of the kagome lattice. Using JS(S + 1)/ → JS 2 for classical spins, we find that D ∞ = 3/32 JSa 2 0.09375 JSa 2 . This small discrepancy between numerical and theoretical results was already noticed in 1D systems, and is associated with the failure of the short time expansion which should be carried to higher orders 51 .

For T J, it becomes more difficult to obtain a simple theory since other processes appear beside spin diffusion. However, by considering the temperature dependence of the constant ratio ( ω 2 / ω 4 ) 1/2 in the whole temperature range, it is possible to rewrite Eq.( 16) as a function of the macroscopic susceptibility and internal energy 52

D T ∝ U (T )/χ(T ). ( 17
)
Although the O(N ) model does not reproduce quantitatively the simulations, it is possible to capture the global shape of the diffusion coefficient above T /J = 0.05 (see red curve in Fig. 5) using the analytic expressions of U (T ) and χ(T ) derived in ref. [49] Finally, in the very low temperature limit T J, the infinite-component spin model coupled to a Langevin dynamics (see Sec. IV A and Appendix A) predicts a temperature independent diffusion coefficient. Fig. 5 shows that D T reaches a plateau below T /J = 0.1 at around 0.37(1) JSa 2 . Moreover, from ref. [49] and using expression (17), we obtain a ratio (D T =0 /D ∞ ) O(N ) = 3 between zero and infinite temperature. This quantitatively agrees with the value 2.8(4) determined from our simulations, while extrapolating the value of the plateau down to T = 0 from the behavior observed around T /J ∼ 0.1 (see figure 5). At lower temperature T /J < 0.05, the O(N ) model rapidly fails since it does not capture the entropic selection of the coplanar states. The diffusion coefficient seems to diverge when the temperature reaches the octupolar transition, while the wave-vector range of validity of spin diffusion, which is restricted at low temperatures by the condition qξ < 1 with ξ the correlation length 51 , shrinks to very small wave-vectors (see inset of Fig. 5).

From our simulations, it is not possible to state that the long-wavelength diffusive behavior disappears in the coplanar regime in favor of propagative spin transfers or if it is simply reduced to a q-range smaller than the resolution δq = (N a) -1 , denoting that the correlation length becomes larger than the lattice size. In this latter case, bigger lattices should be considered to avoid finite-size effects. In any case, diffusive behavior may exist even in a long-range ordered AFM as long as non-linear effects such as interacting spin waves are significant. These interactions are particularly strong in frustrated magnets even at very low temperatures 53,54 , so spin diffusion may still be present in the coplanar regime, although being limited to very long wavelengths and negligible in intensity compared to propagative spin transfers.

Relaxation at generic wave-vector : lifetime of the ground states

The autocorrelation function A(t) gives useful information about the global relaxation of the system and may be an efficient way to probe the evolution of the stiffness with the development of correlations at low temperature. In a previous study 33 as well as in fig 3 of this paper, it is shown that a decaying exponential qualitatively describes the autocorrelation function in the paramagnetic and cooperative paramagnetic regimesat least in a certain time range, this range being highly reduced in the paramagnetic regime because of the 1/t diffusive tail. Thus, the wave-vector-averaged relaxation time can be extracted for each temperature, its inverse Γ r (T ) = τ r (T ) -1 being represented in Fig. 9. The relaxation rate Γ r (T ) goes from the constant value J in the paramagnetic regime, to an algebraic law AT α in the cooperative regime with α = 0.94(3) close to one, reflecting a slowing down of the spin fluctuations.

For comparison, note also that a similar result has been obtained in the pyrochlore antiferromagnet using simulations and phenomenological arguments 7,32 .

Nevertheless, A(t) = d 2 QS(Q, t) only provides qualitative information in the case of wave-vector dependent fluctuations. So it is necessary to study the Q-resolved scattering function S(Q, t) to understand the overall dynamical properties of the system. In the following, we assume an exponential decay of the scattering function at each Q

S(Q, t) = S(Q, 0)e -t/τ T (Q) ( 18
)
and use the same treatment as the one discussed in the previous section for long wavelengths. The relation (18) was checked to be a good approximation at most wavevectors. In particular Eq. ( 18) is justified in the paramagnetic and cooperative regimes because the long time dynamics is dominated by relaxation processes away from nodal lines. In that case, propagating excitation can be neglected in the first approximation 33 . So, extracting the relaxation time τ T (Q) from the numerical data allows us to distinguish the dynamical properties of shortrange correlated domains having a propagation-vector Q at temperature T . The open squares represent the relaxation time obtained after integration over wave-vectors away from nodal lines where τ is predicted to be q-independent from ICSM (see text), with a 1/T law which is represented in black.

instantaneous scattering function S(Q, t = 0) are shown in Fig. 7 along the [h, h] direction for several temperatures.

In the cooperative regime (Fig. 6 (b)), the relaxation time seems to be nearly independent of wave-vector in regions away from the nodal lines and the zone center where the dynamical properties are dominated by diffusion and spinwave like processes. This result is very similar to the predictions of the O(N ) model and is rather intuitive in the light of Ref. 35 since no particular correlations are favored in this temperature range, all locally ordered domains having roughly the same relaxation time. More generally, it is striking that the O(N ) model reproduces accurately the overall q-dependence of the relaxation time above T /J = 0.05.

The simulations show that τ T (Q) (Fig. 6 (c)) becomes more structured at lower temperatures. In particular, it is clear from Fig. 7 (a) and (b) that the longest relaxations coincide with the correlations peak of the √ 3 × √ 3 phase that are located around BZ vertices. This slowing down of the spin fluctuations at BZ vertices, obviously not reproduced by the O(N ) model which does not take into account order-by-disorder phenomena, is thus observed at higher temperatures than the onset of √ 3 × √ 3 static correlations, which occurs only when T /J 0.005 26 .

The temperature evolution of the relaxation time close to and away from the BZ vertices is represented in blue and red in the inset of figure 7 (b). The shaded region between these two curves symbolizes the relaxation time distribution (for wave-vectors contained in the blue triangle of Fig. 6 (b)). The open black squares are the mean value of this distribution. Although the O(N ) model neglects the wave-vector dependence of τ T (Q) below T /J = 0.05, the wave-vector-averaged relaxation time is roughly consistent with the law τ T (Q) ∝ T -1 obtained in sec.IV A at low temperatures.

To conclude, these results suggest that the fluctuations around the ground state manifold are strongly affected by the entropic selection far above the transition toward coplanarity which occurs only at T /J = 0.005. Contrary to the antiferromagnetic Heisenberg model on the py-rochlore lattice 32 , the O(N ) spins model describes only qualitatively both the diffusive and higher wave-vectors regimes in a restricted temperature range 0.05 < T /J 1.

V. COPLANAR REGIME A. Low temperature landscape

The goal of this section is to motivate and justify the investigation of the Q-resolved dynamical scattering function in the time range of interest.

In the coplanar (octupolar) regime 20,26 , the incoherent spin dynamics induced by thermal fluctuations is not the only channel of relaxation. In this low temperature phase, the low energy manifold can be thought as a neighborhood of all 3-colorings of the kagome lattice, which form a discrete manifold, all colorings being globally related one to another by continuous bridges, the collective rotation of spins belonging to 2-colored loops of the discrete manifold. As a consequence, this low energy manifold is connected, even though it must be noted that from the discrete point of view, only the dynamics within the discrete 3-colors manifold induced by the 2color loop move leads to a non-connected structure, the manifold being split into distinct Kempe sectors 28,73 , i.e. spin configurations which cannot be related one to another through loop moves only. In other words, it is only possible to go from one sector to another through the use of a defect, a highly unlikely event at low temperature. Consequently, in a typical time scale of tJ < 1000, one can consider that the system is trapped in a Kempe sector and does not escape it.

Whatever the sector the system is trapped in, there exists loops of different lengths p = 2 + 4n with n > 0 36,55 . Using periodic boundary conditions, loops can be divided into two categories; winding and non winding loops. One may expect different dynamics for these two families. Actually, we will now see that at low temperature, the microscopic spin model we are interested in, in the time range of interest, discriminate even more drastically within each family.

Let us consider two archetypal 3-coloring, the well known long range ordered q = 0 and √ 3× √ 3 spin configurations whose shortest weathervane modes are respectively infinite lines and small loops of 6 spins. In order to mimic a very low temperature regime, we introduce a small amount of energy in each phase, by uniformly randomizing each spin configuration with a ∆E 6.4 10 -4 J . Then, equations of motions are integrated and time evolution for each case is represented in Fig. 8, with (a), the autocorrelation of each spin in direct space and (b), the associated static structure factor S(q, ω = 0). While hexagonal loops are activated and their number increases with time for the √ 3 × √ 3 spin configuration, no flipped loop is detected for the q = 0 phase. Note also that for the √ 3 × √ 3 spin configuration, no loop of length greater than 6 occurs. loop of length L > 6 as well as infinite loops are therefore absent at this time scale and do not play any role in the dynamical properties. In reciprocal space, this results in a negligible diffuse spectral weight at ω = 0 for the q = 0 phase in opposition to the √ 3× √ 3 phase (see Fig. 8

(down)).

Because it is now well established that entropy stabilizes √ 3 × √ 3 correlations at low temperatures, one may expect that thermodynamically, spin configurations belong to the corresponding Kempe sector. On the time scale of the simulations, one may consider the phase space to be the 6-loops neighborhood of this configuration, i.e all accessible configurations starting from the pure √ 3 × √ 3 phase and applying non-overlapping 6loop move, keeping in mind that all operations are not commutative.

Therefore, while in reciprocal space √ 3 × √ 3 spin pair correlations gives rise to sharp peaks located at the Brillouin zone vertices, their width being inversely proportional to the correlation length ξ √ 3 (T ), the presence of flipped hexagonal loops yields an elastic diffuse spectral weight in reciprocal space, since the presence of those 'defects' in the parent periodic structure requires an infinite number of Fourier components. The form factor of one such 'defect' is represented in inset of Fig. 8 (b) and indeed results in broad bumps softly stretched along BZ edges. Consequently, the instantaneous structure factor is expected to be a superposition of both sharp and broad features located at different regions of reciprocal space, and whose origin are of different nature : probing the dynamics at different wave vectors will give information on different relaxation processes.

B. Models, predictions about time scales

The analytic approach described in Sec. IV A obviously fails to describe the fluctuations around such a ground state manifold, simply because it neglects the order by disorder phenomena occurring at very low temperature. However, using qualitative arguments, it is possible to roughly predict how the dynamical properties would evolve in the presence of an entropically induced potential well.

In a first approximation, let us consider the time evolution of a single loop diffusing in such a landscape, whose dynamics is described by a simple stochastic differential equation. In this approach, we also ignore the interactions between the weathervane modes and the spin waves sensitive not only to the ground state manifold (in the sense of internal energy), but also to the excitation spectrum.

This dynamics should have, at sufficiently low temperatures (T V L with V L the height of the free energy barrier), two distinct timescales, corresponding to (i) the required time to overcome the barrier and flip the loop, (ii) the weak loop fluctuations around the plane of coplanarity. While the latter time scale will mainly affect the out-of-plane component for sufficiently small fluctuations (i.e. at sufficiently low temperatures), flipping a loop will also influence the in-plane channel since such a motion induces a change of three-coloring.

Classically, the in-plane relaxation associated with loop motions should obviously undergo a reduction of the number of flips with decreasing temperature, described by the activation law

τ = τ 0 exp(-V L /T ). ( 19
)
An estimation of the energy barrier height V L has been obtained within Gaussian spin-wave theory 36,42,56 . In particular, it was shown that the π-periodic potential well induced by quantum fluctuations has the form V (φ) = ηL| sin(φ)|, with η = 0.14 36,55,56 and φ the angle between the "averaged" coplanar spin plane and the plane defined by the spins of the considered loop. Therefore, in the low temperature limit where quantum fluctuations dominate, V L = V (φ = π/2) ∝ L only depends on the loop length L. However, in the presence of substantial thermal fluctuations (classical limit), the barrier height is renormalized V L = T L log(2ηJS/T ) 36,42 . Combining this latter expression with equation ( 19) leads to the power law

τ = τ 0 2ηJS T L , (20)
whose exponent is the loop length. The behavior of the second time scale, i.e. weak fluctuations of the loops within the entropic potential well, strongly depends on the precise shape of the well and is more difficult to handle. Indeed, for small angle, the fluctuations of the neighboring loops cannot be neglected anymore. In particular, it was shown that taking into account these loops, interactions round out the well bottom such that its φ dependence becomes quadratic for φ φ 0 , with φ 0 = φ 2 1/2 = T /J the rms induced by thermal fluctuations 55 . The stochastic Langevin equation in a quadratic well V (φ) = d(T )φ 2 , with d(T ) an effective onsite planar anisotropy coefficient which possibly depends on the temperature, may be solved analytically. This leads to a relaxation time τ ⊥ ∝ 1/d(T ), pointing out that the temperature dependence of the out-of-plane relaxation time follows that of d(T ) : for instance, a temperature independent well bottom would lead to a constant out-of-plane relaxation time, leading to no freezing effects down to T = 0 (in this particular case, this would actually be the amplitude of the fluctuations which would tend to zero with decreasing temperature).

In the next section, we numerically test these ideas and try in particular to prove the presence of several characteristic time scales in the spin dynamics. We also qualitatively discuss the possible role of the interactions between weathervane and spin waves modes, which break the 120 • rule between neighboring spins and lead to incoherent spin fluctuations.

C. Numerical results and discussion

The autocorrelation function A(t) plotted on figure 3 exhibits different behaviors depending on the considered temperature range : paramagnetic, spin liquid or coplanar regimes. Although the autocorrelation in the spin liquid phase can be described by a single decaying exponential, such treatment is not valid (see Sec.IV B) far above the transition toward coplanarity (for T /J 0.05)-. Below the crossover temperature, Fig. 3 (c) shows a more complex behavior with at least two time scales separated by a crossover at around tJ 60. Since the coplanar regime is (by definition) anisotropic in spin space, different dynamical behaviors are expected for the in-plane and out-of-plane components, each one being likely associated with different relaxation processes. Separating the two contributions A (t) and A ⊥ (t) of the autocorrelation function in our simulations appears natural. The short time scale (tJ < 60) is governed by out-of-plane relaxation, while the in-plane relaxation governs the long time regime (see figure 9 (c) for T /J = 0.0004). From these considerations, A(t) can be split into two exponential contributions below T /J 0.005

A(t) = a ⊥ e -t/τ ⊥ + a e -t/τ , (21)
with τ ⊥ τ the relaxation times and a ⊥/ =

1 N i (S ⊥/ i ) 2 ∈ [0, 1]
the amplitudes of the out-ofplane/in-plane fluctuations such that a + a ⊥ = 1, gives a good agreement with the numerical data.

τ behaves like the relaxation time obtained in the spin liquid regime: it follows an algebraic law AT α with a slightly higher exponent (α 1.2(1)) denoting a slowing down of the spin fluctuations below the transition, the coplanarity inducing a stiffness in the spin texture. This exponent value remains qualitative and may slightly be influenced by finite size effect. Meanwhile, τ ⊥ 35(5) J -1 seems to be roughly temperature independent. These different dynamical behaviors come with a decrease of the out-of-plane spin component a ⊥ ∝ T 1/2 (see Fig. 9 (b)), in agreement with the equipartition theorem in the presence of out-of-plane quartic modes.

Considering the dynamical scattering function S(Q, t), which gives access to the wave vector dependence of the relaxation times, yields more insight about the underlying mechanism leading to the strongly different dynamical behaviors of the in-plane and out-of-plane spin components : by avoiding the wavevector-averaging effects, we can detect the regions of reciprocal space leading to such a behavior. As for the higher temperature results, S(Q, t) is dominated by quasistatic relaxation for wave vectors away from the nodal lines 76 . So following the in-plane/out-of-plane separation performed for the autocorrelation, the scattering function can be approximated for each Q-value by

S(Q, t) = S ⊥ (Q, 0)e -t/τ ⊥ q + S (Q, 0)e -t/τ q . ( 22
)
Some of those fits along the Q = (h, h) directions are shown in Fig. 10, first pointing out that the peculiar shape of the autocorrelation function in the coplanar regime is not induced by the wave-vector averaging, since the two characteristic time scales are also observed at each wave-vector value. Then, the relaxation times obtained by fitting the scattering function at each wave vector using expression (22) are plotted as intensity maps in reciprocal space for T /J = 0.001 on Fig. 11 (a) for τ ⊥ q and (b) for τ q .

Handling the in-plane scattering function is a bit tricky, seeing that the static spectral weight is a combina- (Q, T ) in reciprocal space for T /J = 0.0006, extracted by fitting the scattering function S(Q, t) at each wave-vector using equation (22). The blue lines bound the Brillouin Zones. The black thick lines and disks hide the regions where the quasistatic scattering function vanishes (nodal lines), so that using Eq. ( 22) is meaningless around these positions. (c) and (d), obtained around Q = (1/2, 1/2), respectively show the in-plane and out-of-plane characteristic fluctuation time of the loops.

tion of two components with very different signatures in reciprocal space as discussed in subsection V A : different relaxation processes or lifetimes are probed depending on the wavevector value. Around the BZ vertices, the static spectral weight is overwhelmed by the sharp peaks resulting from √ 3 × √ 3 correlated domains. Thus, the time evolution of the scattering function around these positions unveil the lifetime of these locally ordered states, which, from our simulations, seems to diverge with decreasing temperature. However, the static and dynamical properties around BZ vertices are strongly affected by finite size effects at such low temperatures since the correlation length ξ √ 3 (T ) may reach the lattice size. Consequently, it is not possible in the current work to quantitatively describe the temperature evolution of the relaxation of √ 3 × √ 3 correlations. On the other hand, the relaxation of the diffuse spectral weight at generic wavevectors is representative of loop dynamics. The corresponding time scales of those local motions in direct space are almost independent of the system size for a sufficiently large number of spins.

The in-plane components provide the average time to flip hexagonal loops in given √ 3 × √ 3 domains (which is different from the flipping motion by itself which has already been discussed in a previous article 33 ). Indeed, the in-plane spin correlations are at the first order not sensitive to small fluctuations of the loops, and full loop flips are naturally needed to alter three-coloring states. Figure 11 (a) points out that the averaged time to flip the loops is smaller than the lifetime of the √ 3 × √ 3 correlated domains. Each weathervane loop may be flipped several times before the √ 3 × √ 3 correlated domains to be relaxed. However, since flipping an hexagonal loop requires to overcome the free energy barrier separating the two neighboring three-colorings (the ones before and after the flip), a depletion of the flipping events with decreasing temperature is naturally expected. According to Eq. ( 20), the required time to flip a loop follows an algebraic law τ (T ) ∝ T α with α is equal to the loop length. Numerical data obtained around Q = (1/2, 1/2) rlu, shown in figure 11 (c), are in very good agreement with a power-law behaviorbut the fitted exponent α 1.5(2) is around four time smaller than the prediction α = 6 for hexagonal loops. This discrepancy could be due to the interactions between the local (loops) and non-local (spin waves) modes, which have been neglected and probably lead to significant incoherent thermal fluctuations. Note also that finite size effects although strongly reduced far away from BZ centers and corners (see sec. III) can not be totally excluded.

Away from the high symmetry directions of the Brillouin zone, the characteristic time-scale of the out-ofplane fluctuations seems roughly flat with wavevector. This result is consistent with the local spin motions in direct space (see Fig. 11 (b)) and suggests that the relaxation is mediated by the loop fluctuations at very low temperatures. τ ⊥ Q weakly depends on temperature and is around 28J -1 for T /J 0.005 (see Fig. 11 (d) at wave vector Q = (1/2, 1/2) rlu), as previously noticed for the Q-integrated scattering function A(t).

The presence of temperature independent spin fluctuations is remarkable for a classical system, whose dynamics generally slows down when the temperature decreases. It is however consistent with loops slightly fluctuating around the plane of coplanarity, if considering a tem-perature independent entropic well bottom V (φ) = dφ 2 (see sec. V B). Nevertheless, to go beyond these phenomenological considerations and confirm these numerical results, theoretical predictions of the precise temperature dependence of the entropic well are necessary in the limit of small angles φ.

In conclusion, the numerical results show that the weathervane loop fluctuations control the system relaxation. We identify two distinct timescales associated with the inplane and out-of-plane fluctuations and find that the temperature and wavevector dependences of these two components are qualitatively consistent with loops diffusing in the entropically induced potential well. However, the exact role of incoherent thermal fluctuations remains ambiguous and needs a better understanding. To go further, a thorough numerical study in direct space (which is now in progress) is required in order to separate more efficiently the dynamics of the (local) loop motions from the other contributions.

VI. COMPARISON WITH EXPERIMENTS

Experimental realizations of kagome antiferromagnet are often complicated by further neighbor and/or anisotropic interactions, single ion anisotropy, spinlattice coupling, chemical imperfections as well as lattice distortions [57][58][59][60][61] . The ground-state manifold is extremely unstable towards such perturbations, which may partially or totally lift the degeneracy 62 , so that any quantitative comparison with simple models becomes difficult. The disappearance of the nematic order parameter when the magnetic lattice contains defects, or the stabilization of a q = 0 ordered state when Dyaloshinski-Moryia interactions are included, are two major illustrations of the effect of perturbations 37,41,63 . Nevertheless, some compounds maintain a spin liquid behavior (often coexisting with spin freezing) down to the lowest temperatures, which show qualitative similarities with our present numerical work on the simple antiferromagnetic Heisenberg model.

As described in sections IV and V, fluctuations around the ground-state manifold show a complex behavior which changes when the magnetic system tends towards coplanarity. In the liquid regime, spin relaxation is the result of incoherent thermal fluctuations leading to an almost linear temperature dependence of the relaxation rate. Such behavior was recently observed by inelastic neutron scattering measurements in deuterium jarosite, an experimental realization of a kagome lattice with spins S = 5/2 64 . The static correlations of this system are very well reproduced by Monte Carlo simulations 64 , so our classical approach could be fruitful to describe its dynamical properties as well. Neutron measurements have been performed from T = 14 to 240 K, which, considering the effective coupling constant J cl = JS 2 = 244 K 64 , probes both the paramagnetic and liquid regimes (0.05 T /J 1). The relaxation rates obtained experimentally and numerically have the same order of magnitude and show qualitative agreement over all the probed temperature range: for instance, data collected at T = 240 K (resp. 15 K), giving T /J = 0.82 (resp. 0.06), provide a relaxation time τ J 3.2 (resp. 10.5), while 1.1(3) (resp. 6.2(6)) is obtained from numerical results. Interestingly, a linear dependence of the fluctuation rate is also observed in the quantum spin-1/2 kagome compound Cu(1,3-benzenedicarboxylate) by means of muon spin spectroscopy 65 , with a relaxation time that is one order of magnitude larger than predicted by the simulations. In this compound, a saturation of the relaxation rate is observed at lower temperatures. This could be due to the presence of sizable quantum fluctuations which are not taken into account in the present study.

When approaching coplanarity, a distribution of time scales, extending over approximately one order of magnitude for a given temperature (see for instance Fig. 7), is also observed numerically. This distribution is induced by the entropic selection that favors √ 3 × √ 3 correlations and leads to a longer lifetime for this type of spin configurations. Such a time scale distribution has been observed experimentally in the deuterium jarosite, for which a non lorentzian line shape of the quasi-elastic intensity of the neutron scattering data was observed at low temperatures 64 .

Below the transition, the collective motion of the hexagonal loops mostly control the spin relaxation in numerical studies. One consequence is the apparition of a temperature independent second timescale that is associated with the out-of-plane fluctuations of the hexagonal loops. Recent experiments on Gadolinium Gallium Garnet (GGG), a three dimensional generalization of the kagome lattice with Heisenberg spins, also reported the observation of distinct time-scales with very different temperature dependence. In this system, the different time scales are associated with the simultaneous development of short-range order dimerization dynamics, cooperative paramagnetism, static order, and finally fluctuating "protected" spin clusters, so that the time scale distribution extends over several order of magnitudes [66][67][68] . Dynamic magnetization measurements also reported that the protected spin clusters fluctuations are not thermally activated and do not depend on temperature. They concluded that the protected spin clusters are quantum dynamical objects 66 . Our results suggest that such a temperature independent behavior does not necessarily need quantum fluctuations and may also be observed in classical systems. However, for a more quantitative comparison, we should consider the real three dimensional crystal structure of GGG as well as the dipolar interactions, which have the same order of magnitude than the exchange and whose role in the dynamical properties is still unclear.

More generally, time-scale distributions are a feature of many frustrated compounds, often characterized by the coexistence of a fast dynamics together with a nonconventional glassy behavior. Unfortunately, the glassy behaviors can not be observed since the algorithm used for solving the dynamics does not accurately describe long time dynamics. Freezing effects may however be studied using stochastic spin dynamics method. Monte Carlo simulations applied to the q = 3 Potts model for instance show the presence of a freezing time-scale, associated with the rearrangement of the clusters with a typical length of few tens of spins 28 .

Finally, to complete this comparison, it is also necessary to discuss the fastest spin dynamics J, associated with spin wave propagation. Finite energy excitation exist in the two low temperature correlated regimes 33 , but theirs intensities are weak compared to the quasistatic spectral weight. These spin waves-like excitation can be identified in the scattering function S(Q, t) as small amplitude but high frequency modulations (of the order of ω J). However as shown on Fig. 12 (c,d) for two different directions in reciprocal space (see Fig. 12 (a)), spin waves excitation have clear signature when we consider the scattering function S(Q, ω) in the frequency domain. The detailed analysis of these excitation along the high symmetry direction awhere there is no quasi-static spectral weightreveal that they are propagative in both coplanar and cooperative regimes, although their lifetime τ SW < τ ⊥ τ is strongly sensitive to the selection of the coplanar ground state manifold 33 . It is there-fore intriguing that no evident dispersive features have been detected so far in liquid-like kagome compounds. In these systems, single crystals are often not available because of technical growing difficulties, so experiments are performed on powder samples. The absence of dispersive excitation could then arise from this powder averaging, which motivated us to calculate the excitation spectrum for powder samples.

The powder averaged intensity maps in (|Q|, ω) space are shown in Fig. 12 in the cooperative (e) and coplanar (f) regimes. It appears that the inelastic excitation spectrum is mostly dispersionless in both regimes in spite of existing propagative spin waves in the single-crystal scattering function (Fig. 12 (c,d)). Indeed, the quasistatic fluctuations (whose intensity is orders of magnitude larger than the spin wave spectral weight) overwhelm the excitation spectrum and blur any significant dispersive feature. Then, propagative effects may be very difficult to observe experimentally on powder samples.

Constant energy cuts of the powder averaged scattering function, displayed in Fig. 12 (g,h) for ω = 0 (red), 1 (blue) and 2 J (green), shows that an interval centered around the energy ω/J = 2 should maximize the experimental detection of a dispersive signal. Indeed, the Q-dependence of the scattering function further evolves while approaching the top of the dispersion. At this energy (ω/J = 2), the powder averaging gives rise to a slightly more intense flat band in |Q|. Its intensity is smoothly structured with the scattering vector, and gives broad maxima at different wave-vectors from the static scattering function (see Fig. 12 (g,h)). These results can be compared to experimental results recently obtained in the volborthite, a S = 1/2 kagome compound which shows no signs of long-range order down to 1.8 K in spite of an effective coupling of few tens of kelvin. Although its static correlations and excitation spectrum probably originates from a more complex exchange Hamiltonian 70,71 than the KHAFM, dispersive excitation as well as a flat band at finite energy, likely resulting from powder averaging, have been observed in inelastic neutron scattering on powder samples 70 .

VII. CONCLUSION

The antiferromagnetic Heisenberg model on the kagome lattice is blessed with very rich dynamics in all temperature regimes. Each regime is characterized by a different mechanism of relaxation. At high temperature, the relaxation of the magnetic phase is purely diffusive.

When temperature reaches the cooperative regime, spins are still disordered but algebraic spins correlations start to develop. They are responsible for the exponential relaxation of the magnetic states at short time scales with a relaxation time in 1/T in agreement with previous studies. At long time scale, spin diffusion remains but it is mediated by the dynamics of spins clusters rather than single spins as in the paramagnetic regime.

In the very low temperature regime, entropic selection favors coplanar states and a anisotropic dynamics. Although spin wave can propagate through the system, their contribution to the relaxation is negligible and limited to short times scales compared to the weathervane defects. They are however very important for activating the weathervane defects whose dynamics dominates the intermediate time regime. A careful analysis of the relaxation shows that it is anisotropic and depends on the direction of the fluctuations. The characteristic times have a different temperature dependence, the inplane component following an power law while the out-of-plane component weakly depends on temperature.

A more detailed study of the weathervane defects dynamics is needed to understand the origin of the different temperature dependence of the relaxation time observed in the lowest temperature regime.

FIG. 2 :

 2 FIG. 2: (Color online) Finite size effects on the autocorrelation function A(t) and relaxation time in the octupolar regime. (a) Autocorrelation function A(t) at T /J = 0.0006 for different lattice size from L = 36 (red) to L = 144 (blue). (b) Fit (blue) of the numerical data (red) for the different lattice sizes ranging from L = 36 to L = 144, assuming that A(t) = a e -t/τ + a ⊥ e -t/τ ⊥ . (c) relaxation time τ versus lattice size resulting from the fit shown in (b).

FIG. 3 :

 3 FIG. 3: (Color online) Temperature dependence of the semiclassical spin dynamics as releaved by the autocorrelation function A(t). Autocorrelation function A(t) versus time in the paramagnetic (T /J = 10 -0.1), cooperative (T /J = 0.1 -0.005) and coplanar (T /J < 0.005) regimes respectively represented by dashed-red, solid-green and dotted-dashedblue lines: (a) linear-log scale, (b) log-log scale, and (c) linearlog scale focusing on the coplanar regime.

FIG. 4 :

 4 FIG. 4: (Color online) Spin diffusion in the cooperative regime as revealed by the spin-spin correlation functions. (a) Scattering function S(q, t)/S(q, 0) versus tJ at T /J = 0.17 for several wave-vectors close to the zone center in the [h, h] direction. The fits, performed at each wave-vector using a decaying exponential law (see text), are represented by the thin black lines. (b) Inverse of the relaxation time τ d extracted from figure (a) as function of the wave-vector. The wave-vector region (|q| 0.1) where the spin diffusion law is not valid is shaded.

FIG. 5 :

 5 FIG. 5: (Color online) Diffusion coefficient and wave-vector range of validity of the diffusive approximation (inset) as a function of temperature. The red lines are the predictions obtained by the different models (see text); the value of D∞ is obtained by extrapolating the numerical data from T /J = 1 to 10; and the vertical dotted black line at T /J = 0.005 is for the transition toward coplanarity.

  the scattering function leads to the following form45

FIG. 6 :

 6 FIG. 6: (Color online) Temperature dependence of the relaxation time τ (Q, T ) in the cooperative regime. Intensity maps of the relaxation time τ (Q, T ) in reciprocal space for T /J = 1.31 (a), 0.077 (b) and 0.010 (c), extracted by fitting the scattering function S(Q, t) at each wave-vector using equation (18). The prediction of the infinite-component spins model stands at the bottom right for each temperature. The red lines bound the Brillouin Zones.

Fig. 6

 6 Fig.6 displays maps of τ T (Q) in reciprocal space for temperatures T /J = 1.31, 0.077 and 0.01; the red lines represent the BZ edges. For comparison, the same maps obtained from the O(N ) model (Sec. IV A and Appendix A) are shown in insets. Cuts of these maps as well as the

FIG. 7 :

 7 FIG. 7: (Color online) Influence of the temperature on S(q, t). (a) Relaxation time and (b) instantaneous scattering function S(Q, t = 0) versus Q = (h, h) at many temperatures in the paramagnetic and cooperative regimes. The vertical gray lines at h = n/3 are for BZ vertices (n = 1, 2, 4) and center (n = 3). (inset) Evolution of the relaxation time with temperature at some high symmetry positions in reciprocal space.The open squares represent the relaxation time obtained after integration over wave-vectors away from nodal lines where τ is predicted to be q-independent from ICSM (see text), with a 1/T law which is represented in black.

FIG. 8 :

 8 FIG. 8: (Color online) Dynamics of the weatherwane defects for the uniformly randomized spin configurations in the q = 0 and √ 3 × √ 3 states seen from the correlation functions si(t) • si(0) (∆E 6.4 10 -4 J). (a) autocorrelation in direct space at times tJ = 100, 350 and 1000 (resp. tJ = 100, 2000 and 4000) for the √ 3 × √ 3 (resp. q = 0) spin configuration. Blue (resp. red) disks are for positive (resp. negative) autocorrelation si(0)•si(t) at site i. (b) Resulting static scattering function in reciprocal space of both spin configurations after numerical integration. For comparison, the bottom left inset shows the structure factor of one weathervane mode. (c) Time evolution of the global autocorrelation function.

FIG. 9 :

 9 FIG. 9: (Color online) Algebraic temperature dependence of Relaxation time in the three different regimes. (a) Relaxation rate 1/τ J vs. temperature: red is for the isotropic high temperature regimes (paramagnetic and cooperative) while blue and gray respectively corresponds to the in-plane and out-of-plane fluctuation rates. (b) Out-of-plane amplitude of spin fluctuations in function of temperature. (c) Inplane/out-of-plane separation of the autocorrelation function for T /J = 0.0004.

FIG. 10 :

 10 FIG. 10: (Color online ) Scattering function S(Q, t) for different wavevectors at T /J = 0.0006 as a function of time tJ. The fits obtained using equation (22) at each wave vector are represented in red. Inset : Scattering function for Q = (0.25, 0) along the nodal line pointing out finite frequency features.

FIG. 11 :

 11 FIG. 11: (Color online) q dependence of the parallel and perpendicular relaxation times in the octupolar regime. Intensity maps of the parallel (a) and perpendicular (b) relaxation times τ ⊥/ q

FIG. 12 :

 12 FIG. 12: (Color online) Powder averaged scattering function in the octopular regime. (a,b) Static scattering functions at T /J = 0.015 and 0.0013 in Q-space: the blue lines bound the BZ while the red ones correspond to wavevector directions represented in panels (c,d). (c,d) Scattering function at the same temperatures as function of energy and scattering vector for different directions in reciprocal space. (e,f) Powder averaged dynamical scattering function versus energy and wave-vector modulus |Q|.(g,h) Constant energy cuts of the powder averaged scattering function versus wave-vector modulus |Q| for ω/J = 0 (red), 1 (blue), and 2 (green). Since the finite energy spectral weight is order of magnitudes weaker than the static one, the two latter energies have been multiplied by constant factors in order to superimpose all the constant energy cuts on the same plot.

  Schematic phase diagram of the classical Heisenberg model on the kagome lattice. This model undergoes two different crossovers when the temperature decreases.

	paramagnetic regime
	spin diffusion
	• both correlated and diffusive dy-namics
	FIG. 1: (color online)

Acknowledgements

It is a great pleasure to acknowledge discussions with S. Viefers, J. Chalker, O. Cepas, A. Ralko, P. Kopietz and S. Petit, as well as N. Shannon and L. Jaubert for critical reading of the manuscript. M. Taillefumier thanks the Max Planck Institute for the complex systems for hosting him during the preparation of this work. This work was performed on the Abel Cluster, owned by the University of Oslo and the Norwegian metacenter for High Performance Computing (NOTUR), and operated by the Department for Research Computing at USIT, the University of Oslo IT-department. M. Taillefumier acknowledges financial support from the Norwegian research council and the University of Frankfurt.

We describe in this appendix the derivation of the different expressions given in the main section Sec. IV A. The starting point of all calculations is the energy functional of the O(N ) model

where the index α represents the different triangles of the kagome lattice and l α is the sum of the components of the spins forming the triangles. Eq. (A1) differs from Eq. ( 2) by an additional term that is introduced to mimic the behaviors the Heisenberg spins whose Lagrange multiplier λ is fixed by the condition s 2 i = 1/3. Eq.A1 can conveniently be written as

where A adj ij is the adjacent matrix of the kagome lattice. By symmetry the adjacent matrix is diagonal in q space so it is possible to express Eq A2 in term of the collective variables s i (q) = R s R,i exp(iq(•R + r i ) where the index i is the sublattice index. The energy functional is then given by

A ad ij (q) are the matrix elements of the Fourier transform of the adjacency matrix A ad (q):

Then the eigenvalues of Eq. ( 13) can be deduced from the eigenvalues ν l of A ad (q) associated with the eigenmodes sl (q). We note P (q) the unitary operator that transforms the operator (A4) in the diagonal form. We find after some algebra that the eigenvalues of A ad (q) are given by

(A6)

As explained in the main section, we describe the spin dynamics with a Langevin equation given by Eq.13. The equation of motion of the collective variables s i (q) can be deduced by direct calculation of the Fourier transform of Eq. ( 13). We find that

where ξ i (q, t) = ξ i,r (t) exp(iq • r) is the Fourier transform of the white noise term ξ i (t) and [s(q)] is the vector formed by the collective variables s i (q). All indexes in Eq. (A5) refer to the sublattice index of the kagome net. After expressing Eq. (A8) in the diagonal basis we find that

where P αβ (q) are the matrix elements of the operator P (q). The solutions of Eq. (A9) are given by sα (q) = s0 α (q) exp -

using Eq. (A10), we find that the spins correlations functions are given by sα (q, t)|s β (q, 0) =

(A13) which combined with Eq. ( 4) gives rise to S(q, t) = α g α (q) sα (q, t)|s α (q, 0) , (A14) and

Using

we find that λ = 1 + O(T ) with s 2 i = 1/3 at low temperature. 73 B. Mohar and J. Salas, J. Stat. Mech., P05016 (2010). 74 By the Mermin Wagner theorem, the corresponding orders are cut-off at the longest lengthscales, leaving well-defined regimes separated by cross-overs. 75 This effect is relatively small for large systems and can be compensated using the properties of the octupolar phase without modifying the equation of motions. It is worth noting that no finite size effects has been observed for the short time scales. Therefore, the analysis of the spin wave spectrum determined in ref. [33], as well as the evolution of their lifetime versus temperature is justified. However, the long time dynamics qualitatively discussed in ref. [33] (see Fig. 3 (b)) is affected by finite size effects in the coplanar regime. 76 Note that this approximation is actually fully compatible with previous results pointing out the existence of propagative spin waves, as long as the intensity of the modes is small compared the static intensity S(Q, ω = 0). The spin wave excitations contribute to the scattering function S(Q, t) by a negligible modulation of high frequency (of the order of J).