Networks of Coadjoint Orbits: from Geometric to Statistical Mechanics

A class of network models with symmetry group $G$ that evolve as a Lie-Poisson system is derived from the framework of geometric mechanics, which generalises the classical Heisenberg model studied in statistical mechanics. We considered two ways of coupling the spins: one via the momentum and the other via the position and studied in details the equilibrium solutions and their corresponding nonlinear stability properties using the energy-Casimir method. We then took the example $G=SO(3)$ and saw that the momentum-coupled system reduces to the classical Heisenberg model with massive spins and the position-coupled case reduces to a new system that has a broken symmetry group $SO(3)/SO(2)$ similar to the heavy top. In the latter system, we numerically observed an interesting synchronisation-like phenomenon for a certain class of initial conditions. Adding a type of noise and dissipation that preserves the coadjoint orbit of the network model, we found that the invariant measure is given by the Gibbs measure, from which the notion of temperature is defined. We then observed a surprising `triple-humped' phase transition in the heavy top-like lattice model, where the spins switched from one equilibrium position to another before losing magnetisation as we increased the temperature. This work is only a first step towards connecting geometric mechanics with statistical mechanics and several interesting problems are open for further investigation.


Introduction
The main purpose of this paper is to provide a link between geometric mechanics and statistical mechanics following recent advances in stochastic geometric mechanics [LCO08,ACC14,Hol15,ADCH18,AHS17], where the inclusion of a structure-preserving noise in Lagrangian or Hamiltonian systems have been considered. Geometric mechanics provides a mathematical framework for describing Hamiltonian mechanical systems, but with the addition of a structure-preserving noise and dissipation, a theory of statistical geometric mechanics becomes possible, as suggested in [ADCH18]. Here, we take a step further to include lattices, or more generally, networks of these systems and study their stability and phase transitions, similar to the well-known Ising or Heisenberg models.
The Ising model is one the most well-studied systems in statistical mechanics, which is a simplified model of ferromagnetism. Its significance in statistical mechanics comes from the fact that it is the simplest model that exhibits phase transition, which is ubiquitous in the theory of matter, such as in the transition from liquids to gases or in the magnetisation of materials. See for example Chandler [Cha87] for a classic reference in statistical mechanics. Relaxing the discrete nature of the Ising model to a system of interacting unit vectors in R 3 on a lattice, we obtain the classical Heisenberg model, which, again is a simple model of ferromagnetism but now with a continuous symmetry group SO(3), instead of the discrete C 2 -symmetry of the Ising model. This model exhibits a phase transition, albeit one with a different universality class due to the difference in the symmetry group.
On the other hand, geometric mechanics is a mathematical discipline that deals with the dynamics of a 'few-body' mechanical system with symmetry and mainly concerned with issues such as symmetry reduction, integrability and stability of equilibria. For basic references in geometric mechanics, see for example Arnold [Arn89], Marsden and Ratiu [MR99] and Holm [Hol08]. Although the two fields are very different, recent progress in the theory of stochastic geometric mechanics has shed some light to bridge the gap between the two. Initiated in the work by Bismut [Bis82] for canonical systems and extended to general Hamiltonian systems in Lázaro-Camí and Ortega [LCO08], stochastic geometric mechanics have found applications in areas as vast as simple mechanical systems [ADCH16], fluid dynamics [Hol15], geometric integrators [BRO09], stability theory [AGH18] and shape analysis [AHS17]. In particular, Holm [Hol15] introduced stochastic processes at the level of the variational principle for fluid systems in such a way that the momentum map is preserved, and this introduced the idea of structure preserving noise, where the added noise does not destroy the essential geometric feature of the system. In a similar spirit, adding a dissipative mechanism in Hamiltonian systems that preserve the basic geometric structure has been considered by Bloch et al. [BKMR96], where a type of dissipative bracket that dissipates energy while preserving the coadjoint orbit is introduced. A similar type of bracket that dissipates a chosen conserved quantity while preserving another, called the selective decay mechanism, is considered in Gay-Balmaz and Holm [GBH13,GBH14].
Our work is based on Arnaudon et al. [ADCH18] where both noise and dissipation is introduced in a general finite-dimensional Hamiltonian system with configuration group G, phase space T * G and symmetry group G, in such a way that it preserves the coadjoint orbits, where the reduced dynamics take place. A remarkable consequence of adding noise and dissipation in this way is that, not only is the geometric structure of the equation preserved but also the invariant measure for the stochastic dissipative system is exactly given by the Gibbs measure restricted to the coadjoint orbit, provided that we choose an isotropic noise. The emergence of the Gibbs measure in this system is key to connecting our theory with statistical mechanics, where it arises naturally as the invariant measure in a canonical ensemble, which is a system of particles in statistical equilibrium that is coupled to an external heat bath. This strongly suggests that we can approach this system from the point of view of statistical mechanics and vice versa.
Our first goal of this paper is to construct a simple lattice model that has a general symmetry Lie group G, or more generally on a network, using symmetry reduction techniques in geometric mechanics and study their deterministic properties. By considering a G-invariant canonical Hamiltonian system on phase space T * G at each node of a given graph, we regard the momentum map J : T * G → g * as the spins in our model, which interact with other spins according to the structure of the graph, where the interaction takes place if they are connected by an edge. The coupling between the spins can be taken in two ways: (1) at the reduced space g * or (2) directly on the group G. In the former, which is the easier case, we introduce the coupling between two neighbouring spins in the reduced Hamiltonian after performing symmetry reduction. We call this approach momentum coupling since the variables that are coupled is the momentum map. We will see that this approach generalises the classical Heisenberg model to include general symmetry groups. For the latter, we consider a representation of G on a vector space V at each node and couple the neighbours in the unreduced Lagrangian that depend on vectors in V * that are acted on by G. Since the neighbours are coupled with the 'positions' of a given vector in V * that is rotated around by the G-action, we call this approach position coupling. The vector that is taken here for position coupling breaks the full G-symmetry of our system and hence we apply the semi-direct product reduction theorem given in Holm et al. [HMR98] to obtain the corresponding Lie-Poisson system on the semi-direct product Lie algebra g * V * at each node. This gives rise to a new system that finds no analogues in classical lattice models. In the special case where the Lie algebra g = Lie(G) is compact and semi-simple, we were able to obtain the equilibrium solutions for both the momentum-coupled and position-coupled systems as the eigenvectors of a generalised graph Laplacian, which we construct from the underlying graph. Furthermore, we were able to classify these stationary eigenvectors into ferromagnetic and anti-ferromagnetic states, which is consistent with known equilibrium solutions of the Ising model or the Heisenberg model and we also investigated their corresponding stability properties using the Energy-Casimir method. Relaxing the condition that each spin must have unit length, many of the equilibrium solutions that we find here are new.
Our second goal is to investigate the statistical mechanical properties of the lattice model constructed above with noise and dissipation of the type considered in Arnaudon et al. [ADCH18]. In particular, we study the phase transition exhibited in the two examples that we consider here: the rigid body lattice and the heavy top lattice which are the simplest systems that can be obtained via momentum-coupling and position-coupling respectively, by taking G = SO(3). The invariant measure is found to be the Gibbs measure restricted to the coadjoint orbit for both of these systems, which allows us to introduce the notion of temperature and hence study their respective temperature phase transition behaviours. For the rigid body network, we observe a standard second-order phase transition in both the mean-field simulation and the direct simulation as expected.
However, in the direct simulation of the heavy top network, we observe a new type of phase transition behaviour where the spins become aligned to two intermediate metastable states as we decrease the temperature before settling down to the lowest energy state, instead of jumping directly to the lowest energy state. This 'triple-humped' phase transition is not captured in our mean-field simulation, which follows a standard 'single-humped' phase transition. In this study, we only looked at numerical simulations of the phase transition but much work needs to be done in the analysis to understand this phenomenon.
1.1. Structure of the paper. In section 2, we begin with a quick review of statistical mechanics and stochastic geometric mechanics and study the statistics of a single coadjoint orbit system in section 3. Then, we introduce the notion of momentum coupling in section 4 to construct a lattice model with symmetry group G and study in detail the equilibrium solutions and their stability properties. We will then add noise and dissipation to the system and show that the invariant measure is given by the Gibbs measure. In section 5, we study an example of this system by taking G = SO(3), which we call the rigid body network. In section 6, we introduce the idea of position-coupling and derive the corresponding network Lie-Poisson model using semi-direct product reduction. Noise and dissipation are then added to the system in a similar way. This theory is then illustrated concretely in section 7 by considering the case G = SO(3), which we call the heavy top network. In section 8, we study the mean-field approximations and the phase transition behaviours of the rigid body network and the heavy top network. Finally, we give a conclusion and discuss further work in section 9.
1.2. List of symbols. This work contains notations from three topics: geometric mechanics, graph theory and stochastic analysis. For clarity, we summarize the main notations we will be using here.

G
a Lie group g an element in a Lie group g Lie algebra corresponding to G µ an element in g * ξ an element in g Statistical mechanics is a physical theory used to describe the macroscopic properties, such as temperature and entropy of a many-particle system (e.g. gas or lattice) evolving as a canonical Hamiltonian system that fluctuates around a mean state. This theory is useful for systems with a large degree of freedom where first of all, solving the full system is computationally expensive and second of all, the information that we need is independent of the motion of individual particles in the system. In equilibrium statistical mechanics, one deals with a system that is in statistical equilibrium, that is, the particles are distributed according to an invariant measure of the underlying system. The system of particles are often assumed to be in statistical equilibrium according to the following three thermodynamic ensembles (1) Micro-canonical ensemble: the system is isolated and there is no energy exchange or particle exchange with the exterior. Each state is equally probable to occur. (2) Canonical ensemble: the system is coupled to a heat bath of fixed temperature and fixed average energy and there is no particle exchange. The states are distributed according to a Gibbs measure. (3) Grand-canonical ensemble: the system is coupled to a heat bath of a fixed temperature and is subject to particle exchange.
For systems with inter-particle interactions such as in lattice models, a phenomenon called phase transition is often observed as the temperature is varied, which describes an abrupt change of state in a system. A classic example of this is the transition from liquid to gas at the boiling point. Now, consider a simple 2-dimenstional lattice model with N nodes that take the values up (+1) or down (−1) at each node. The energy of this system at a particular state is given by for constants J ij > 0 and i ∼ j means that nodes i and j are adjacent on the lattice. We denote by Ω := {1, −1} N for the sample space of all possible configurations. Assuming that the system is in statistical equilibrium in a canonical ensemble with temperature T , the probability of selecting a configuration s = (s 1 , . . . , s N ) ∈ Ω with a given energy level H 0 is given by the Gibbs distribution where the normalising constant Z β is called the partition function and β is the inverse temperature, defined by β := 1 k B T , where k B is the Boltzmann constant. It is well known that there exists some T = T c , called the critical temperature such that the average magnetisation defined by vanishes for T > T c and becomes strictly positive for T < T c . This is an example of a so-called second-order phase transition and the model that we just considered is the 2D Ising model which is the simplest known system that admits a phase transition. Now, instead of taking the values {1, −1}, one can generalise this system so that the spins s i are unit vectors in R 3 and we obtain the classical Heisenberg model which also exhibits a similar second-order phase transition only if the interaction is anisotropic, in which case it is often called the XY Z, or XXZ-model.
The ideas that we discussed here are elementary in statistical mechanics and can be found in most textbooks on the topic. Here, we will refer to Chandler [Cha87] for a more thorough physical exploration of this topic and to [Rue04] for a rigorous mathematical formulation of statistical mechanics.
2.1.1. Other developments in geometric statistical mechanics. Remarkably, in the early days of modern geometric mechanics, one of the founders of the field, Jean-Marie Souriau, had already attempted to apply ideas from geometric mechanics to statistical mechanics. His original idea was to impose natural symmetries on the Gibbs distribution and obtain a vector of inverse temperatures β ∈ g in such a way that the symmetry Lie group G corresponding to g survives in the construction of statistical quantities, such as the entropy. Some of his works in this direction can be found in [Sou66,Sou74,Sou69]. However, this theory attracted very little attention at the time until only recently, where several authors modernized the original work of Souriau. We refer to [Bar14, Bar15,Mar16] as well as to [Bar16] for a complete review on the history of Souriau's theory.
In parallel to this, there has also been recent work by Gay-Balmaz and Yoshimura [GBY17a,GBY17b] based on infinite dimensional geometric mechanics and a certain class of dynamical constraints to describe various systems in the theory of non-equilibrium thermodynamics.
2.2. Stochastic geometric mechanics. A systematic theory for introducing stochastic processes into mechanics started with the work of Bismut [Bis82], where the noise was introduced at the level of the variational principle and the corresponding stochastic Euler-Lagrange equations were derived. This theory has seen a resurgence recently and has been extended by several authors such as Lázaro-Camí and Ortega [LCO08] and Bou-Rabee and Owhadi [BRO09] based on more modern geometric mechanical techniques. Also recently, Arnaudon, Chen and Cruzeiro [ACC14] and Holm [Hol15] added noise in the variational principle that is compatible with symmetry reduction to obtain a corresponding stochastic Euler-Poincaré equation. In the former, the noise was introduced at the level of the Lie group and using symmetry reduction together with taking an expectation, they obtained a dissipative deterministic equation with applications in infinite dimensional systems, such as the Navier-Stokes equation in fluid dynamics. In the latter, which we will base our work on, the noise is introduced directly into the reconstruction relation in the variational principle and followed by symmetry reduction, a fully stochastic Euler-Poincaré equation is derived. Based on this work, Arnaudon, De Castro and Holm [AHS17] studied in detail the finite dimensional analogue of this construction together with double bracket dissipation, and the implication of noise on the nonlinear stability of relative equilibria was investigated in Arnaudon, Ganaba and Holm [AGH18]. We refer to [AHS17] and references therein for more details on other related works.
In the present text, we will use the finite dimensional equations derived in Arnaudon, De Castro and Holm [AHS17], hence we will refer to this work for more details about the derivations of the stochastic equations used here. For a dynamical system with configuration group G, phase space T G and a G-invariant Lagrangian L(g,ġ), where (g,ġ) ∈ T G, we can apply the theory of reduction by symmetry to show that the Euler-Lagrange equation on T G associated to this Lagrangian is equivalent to the Euler-Poincaré equation on g, where ξ = g −1ġ ∈ g, the reduced Lagrangian is l(ξ) := L(g,ġ) and ad * is the dual of the adjoint action. From the definition of the reduced variable ξ, it is possible to reconstruct the solution on the original configuration manifold G, using the formulȧ called the reconstruction relation. The introduction of noise appears at this level in the theory of reduction by symmetry by replacing the above relation by a stochastic process where W l t are K independent standard Wiener processes, d is the stochastic evolution operator and σ l ∈ g are given Lie algebra elements which represents the amplitude and direction of the noise. We interpret this as a Stratonovich SDE, represented by the symbol • so that we can use the standard rules of calculus.
From the above stochastic reconstruction relation, one can compute the stochastic variations and the stochastic Euler-Poincaré equation will follow from the same variational problem exactly as in the deterministic case.
As we mentioned in the introduction, the link between geometric mechanics and statistical mechanics requires a mechanism of dissipation to balance the energy input of the noise and hope for the system to reach a statistical equilibrium. As given in Arnaudon et al. [AHS17], we use the double bracket dissipation introduced by Bloch et al, [BKMR96] and extended in Gay-Balmaz and Holm [GBH13,GBH14]  where h : g * → R is the reduced Hamiltonian, C is a given Casimir function of the Lie-Poisson structure and the momentum variable µ ∈ g * is conjugate to the reduced velocity ξ ∈ g via the Legendre transform. We denote by : g → g * to be the canonical isomorphism between g and its dual given an inner product, and θ ∈ R parametrizes the strength of dissipation. This equation can be derived from a variational principle where the dissipation is inserted as a force in the Lagrange-d'Alembert framework, see [AHS17] for more details.
One can see that the dissipative term and the noise preserves the coadjoint orbit where µ = µ(0). In the following, we will denote a generic coadjoint orbit by O, discarding the foot point. In particular, the Casimir functions C : g * → R (in general, the system admits several Casimir functions, but we only select one here) is conserved whereas the energy decays due to the dissipation and fluctuate due to the noise.
The statistical information of the process can be captured via the Fokker-Planck equation, which describes the time evolution of the probability distribution of the process µ. We refer to [AHS17] for the derivation of the Fokker-Planck equation for the process µ in geometric form, given by where ·, · is the natural pairing on g, {f, h}(µ) = µ, ∂f ∂µ , ∂h ∂µ is the Lie-Poisson bracket and Φ l = σ l , µ are the stochastic potentials. We now state the result in [AHS17] that forms the basis of our present work.
Theorem 2.1. The stationary distribution of the stochastic process (2.8) with an isotropic noise, that is σ i = σe i (for e i a basis of g), is the Gibbs measure on the coadjoint orbit, that is (2.11) where β = 2θ σ 2 = 1 k B T is the inverse temperature (k B is Boltzmann's constant) and Z is the normalisation constant, or partition function (2.12) Proof. One can prove this by a direct substitution into the Fokker-Planck equation (2.10) and obtain d dt P ∞ = 0.
We need to stress an important point here. As opposed to the standard theory of statistical mechanics where the partition function is an integral over the full phase space, our partition function (2.12) is defined as an integral over the coadjoint orbit and not over the whole space g * . This is a result of the nonlinear dissipation and the multiplicative noise and makes the effective phase space compact if the Lie group is compact. It also has the advantage of being an embedding in R n , rather than a general manifold where the definition of stochastic processes is more involved. Notice the requirement of isotropic noise, necessary to obtain the simple Gibbs form of the stationary distribution. In the nonisotropic case, the distribution would be close but different from the Gibbs distribution, see [AHS17] for more details.
Hereafter, we will not use the Fokker-Planck equation as we will assume implicitly that the system is at statistical equilibrium, that is, our system is distributed according to the Gibbs measure. We can therefore apply results from the theory of equilibrium statistical mechanics to study our system. One of the fundamental results is the fluctuationdissipation theorem, which states that for a system to be at a statistical equilibrium, the diffusion coefficient must be a function of the dissipation coefficient or vice versa. In our case, this relation is given by θ = 1 2 βσ 2 , and is often called Einstein's relation. For our system, statistical equilibrium states exist for all pairs (θ, σ) but when the temperature is fixed, one of the variables depends on the other. We will not go further into this issue as it is rather simple as we are dealing with a classical system, but this theorem is important especially for quantum systems. We refer to Kubo [Kub66] for an early theoretical exposition of this theorem, and the many references therein for more details on this topic.

Statistics in a single coadjoint orbit
We illustrate our statistical analysis by first considering the statistics for a single coadjoint orbit. The key object in equilibrium statistical mechanics is the partition function (2.12) since most statistical quantities can be derived from it. Unfortunately, the partition function cannot be integrated analytically in general, except for a few simple examples.
The average energy of the system is defined as (3.1) but can be obtained directly from the partition function as In the case of compact coadjoint orbits, the average energy is bounded from above and saturtes rapidly. Indeed, as T → ∞, the Gibbs distribution P ∞ (µ) → 1 for all µ ∈ O and we have the inequality Similarly, the energy fluctuation around E is defined as As for the mean energy, it saturates rapidly as T → ∞.
We finally arrive at the definition of entropy S on the coadjoint orbit, given by the famous relation or, in terms of the partition function, If the orbit is compact, the entropy saturates as well for large T .
The entropy is an important quantity as it is maximized by the Gibbs measure, under the constraint that the mean energy E is fixed. This is, in a sense the definition of the canonical ensemble in statistical mechanics.
Theorem 3.1. The Gibbs distribution (2.11) is the distribution that maximizes the entropy (3.6) of the system.
Proof. We enforce E = E 0 with a Lagrange multiplier β (this choice will be clear later) and consider the functional (3.7) The variation with respect to P in the variational principle δS = 0 gives the relation log (P(µ)) + 1 + βh(µ) = 0 . (3.8) By normalizing the solution P(µ) to 1, we obtain the Gibbs measure (2.11) with inverse temperature β.

Network of coadjoint orbits I: Momentum coupling
The two important systems that can be described using the theory of statistical mechanics are gas and lattices. The latter, which includes the Ising model and the Heisenberg model are well studied in statistical physics and exhibit many interesting properties such as phase transition. In this section, we will construct a lattice of interacting continuous spins, similar to the Heisenberg model, that take values on a coadjoint orbit, then study in details their deterministic properties and later introduce noise and dissipation to the system that preserve the coadjoint orbits. By coupling the neighbours in the momentum, we will recover the classical Heisenberg model as the simplest example using the coadjoint orbit of SO(3), that is the sphere S 2 . We will mostly restrict our exposition to compact semi-simple Lie algebras here, but extensions to more general Lie algebras should be possible in some cases.

Deterministic equations.
Our aim here is to construct the equations of motion of spins interacting on a general undirected connected network N . At each node, we will assign a coadjoint orbit O i for i = 1, . . . , N , associated to the Lie group G, which is coupled to its neighbours according to the structure of the graph. 4.1.1. Coupling in the momentum. Given an undirected, connected graph N and a left (or right) G-invariant canonical Hamiltonian system on T * G, we introduce the notion of coupling in the momentum to construct a Lie-Poisson Hamiltonian system on (g * ) ⊕N such that at each node of N , the system evolves as a Lie-Poisson system on g * , with an additional interaction term arising from its neighbours. If we suppose for now that at each node, the system evolves independently as a canonical Hamiltonian system on T * G with left G-invariant Hamiltonians H i at node i, then according to the Lie-Poisson reduction theorem, we can collectivise the Hamiltonian as H i = h i • J R where J R : T * G → g * is the momentum map corresponding to the right cotangent lift action of G on T * G and we obtain a Lie-Poisson system on g * with reduced Hamiltonian h i at each node i. The idea of coupling in the momentum is to take the interaction between the neighbours at the level of the reduced space g * , that is, we couple the momentum map J R between neighbouring bodies. Introducing the symmetric and positive definite interaction tensor we take the momentum-coupled interaction potential energy to be where d i is the number of neighbours at node i on graph N and the factor d i d j −1 is taken to normalise the expression, which will be convenient later in our analysis. We take the total energy of the form where i ∼ j means that nodes i and j are adjacent on the graph N . Taking the (−) Lie-Poisson structure on (g * ) ⊕N which is given by the sum we obtain the Lie-Poisson equatioṅ Remark 4.1.
To be more precise, the interaction tensor J ij is a map from g * at site i to g at site j, but since the Lie algebra g is identical at each node, there is no need to distinguish between them. For the case where g = Lie(G) is a compact, semi-simple Lie algebra, we can identify g * with g using the inner product ·, · = −κ(·, ·) where κ is the Killing form, which also satisfies the associativity property a, [b, c] = [a, b], c . From this, one can check that the quadratic functions are Casimirs of the Lie-Poisson bracket (4.4) and the coadjoint orbits of G ×N on (g * ) ⊕N are contained in their level sets, i.e.
where c 1 , . . . , c N are constants, which follows from the Ad * -invariance of the Killing form.
The coadjoint orbits are preserved by the dynamics due to the equivariance of the momentum map J R .
We also consider which is the sum of all the Casimirs. This is also a Casimir for this system and will be used in our analysis later. We will also use the following shorthand notation for a system of coadjoint orbits that are interconnected by a network.
Definition 4.1. Given a graph N with N nodes, a Lie group G, and a coadjoint orbit Hereafter, we consider a simple mechanical system consisting only of a purely kinetic energy term where I i : g → g * is the moment of inertia tensor assigned to node i, which is symmetric and positive definite, and an interaction potential energy term (4.2). We take the total Hamiltonian to be (4.10) In fact, this can be expressed in a more compact form using the language of graph theory, which we will review in the next section.
In the special case where the graph is regular, i.e. d 1 = · · · = d n = d and absorbing the 1 d facor in the interaction tensors J ij , we obtain the Hamiltonian which, ignoring the kinetic energy term, resembles the Hamiltonian for the Ising model or the classical Heisenberg model. Hence, the network Lie Poisson system we obtained by momentum coupling can be viewed as a generalisation of the Heisenberg model with spins taking values on a general coadjoint orbit and with an additional kinetic energy term at each node.
Remark 4.3 (Ferromagnetic vs. anti-ferromagnetic states). From this Hamiltonian, we see that the minimum energy configuration must be an aligned state, with all the µ i pointing in the same direction, and the maximum energy state must be anti-aligned. In statistical mechanics, these two states are called ferromagnetic and anti-ferromagnetic states respectively and will be important in our discussion of equilibrium solutions later.

4.1.2.
Review of graph theory. The structure of a graph is described by the adjacency matrix, which for unweighted graphs, is given by where i ∼ j means that the node i and j are adjacent, or share an edge on the graph. This matrix is symmetric if the graph is undirected, which is the case here. Each node has d i = (A1 N ) i neighbours, where 1 N = (1, . . . 1) ∈ R N , and this is called the degree at node i. From this, we define the graph Laplacian is the degree matrix, and its normalised version which is a symmetric matrix. These Laplacians are usually used to define a random walk on a graph, whereṗ = pL for a probability vector p ∈ R n of a random walker, which corresponds to a discrete version of the diffusion equation, where L plays the role of the Laplace-de Rham operator.
4.1.3. Network of coadjoint orbits. We will now extend this theory to write our system of interacting coadjoint orbits of dimension k using the language of graph theory, in particular, using the graph Laplacian.
We extend our notion of the normalised Laplacian by weighting its components by the inertia tensor and the interaction tensor to get an extended normalised Laplacian The Hamiltonian (4.10) of our system can then be expressed compactly as (4.15) Using this notation, the network Lie-Poisson equation (4.5) with Hamiltonian (4.15) can be written in the compact formμ 4.1.4. Legendre transformation and Lagrangian formulation. We obtained equation (4.16) by coupling the neighbours at the level of the reduced space after performing Lie-Poisson reduction at each node. It is then natural to ask whether we can recover the same equation by coupling the neighbours at the level of the group and then performing symmetry reduction. The answer is yes, except for a few cases where the graph Laplacian is not invertible on the full phase space.
For L to be invertible, one has to ensure that ker(L) = ∅. For the simple case where I i = I and J ij = J, this can be reduced to showing that I and J has no eigenvalues in common. Indeed, if they share an eigenvalue, then at least one of the ferromagnetic states will have a 0 eigenvalue since the extended Laplacian reduces to the standard graph Laplacian L = D − A which has a single 0 eigenvalue. However, choosing J < I, the minimum eigenvalue of L becomes strictly positive since the reduced L corresponding to ferromagnetic states become strictly diagonally dominant.
Assuming that the moment of inertia is invertible on the full space, we obtain the following reduced Lagrangian by the inverse Legendre transform which can also be obtained through the full Lagrangian on T G ×N , given by and L g denotes left diagonal action of g. We refer to appendix 2 of [Arn89] for the reconstruction of the full Lagrangian from a reduced Lagrangian of this form. It is easy to check that the Lagrangian (4.18) is left invariant under G so we can apply Euler-Poincaré reduction on (4.18) and perform Legendre transformation to obtain our network Lie-Poisson equation (4.16). However, notice that the matrix L −1 is far from being sparse, thus the interaction between sites take a complicated nonlocal (more than neighbours interactions) form on the Lagrangian side.
If one really wishes to interpret the system from the Lagrangian side with singular L, the Lagrangian can be defined in terms of the Moore-Penrose pseudoinverse L † . However by doing so, one has to restrict the reduced Lagrangian to the subspace ker(L † ) ⊥ and recover the reduced Hamiltonian on the subspace ker(L) ⊥ by Legendre transformation. We refer to [CHHM98] for more details on how to obtain an equivalent Lagrangian system for degenerate Hamiltonians, but within the context of plasma physics. 4.2. Nonlinear stability results. In this section, we will look for equilibrium solutions of the momentum-coupled network Lie-Poisson system (4.16) on a compact, semi-simple Lie algebra and find their nonlinear stability properties. In particular, we will show that the eigenvectors of the extended graph Laplacian L correspond to equilibrium solutions of this system and they are either a ferromagnetic or an antiferromagnetic state. Furthermore, we will show that the eigenvectors with the lowest and highest eigenvalues correspond to nonlinearly stable equilibria. 4.2.1. Equilibrium solutions. Since we are on a compact semi-simple Lie algebra g, we can take ad * = −ad, so equation (4.16) becomeṡ (4.19) In general, the solutions toμ = 0 are given by all µ e ∈ g * such that Lµ e is contained in the centralizer Z(µ e ) of µ e . However, here we will only consider the case where Lµ e = λ e µ e for some λ e ∈ R, which is clearly contained in Z(µ e ). So fixing a > 0, a regular value of C : g * → R, we consider relative equilibrium configurations on the level set C −1 (a) that are given by the kN eigenvectors of L, rescaled appropriately to satisfy C(µ e ) = a. Now, pairing both sides of Lµ e = λ e µ e with 1 2 µ e , we get so the eigenvalue of L is proportional to the total energy of the equilibrium state.
In the special case where J ij = J and I i = I for all i, j = 1, . . . , N , we see that these kN eigenvectors can be categorised into two types: ferromagnetic or anti-ferromagnetic states, as given in the following proposition.
Proof. Consider the orthogonal decomposition of (g * ) ⊕N into ferromagnetic and antiferromagnetic It is easy to check that V ⊥ is the orthogonal complement of V with dim(V ) = k and dim(V ⊥ ) = (N − 1)k. We claim that V and V ⊥ are invariant subspaces under the linear operation L.
To verify the former, take µ ∈ V , so µ i = √ d i µ for all i = 1, . . . , N . Then, where µ := I ext µ ∈ g ∼ = g * and we used the fact that For the other case, take where we used the fact that i A ij = d j and the fact that µ ∈ V ⊥ in the last line. Hence, This implies that there is an appropriate change of basis such that L becomes block diagonal, where v i ∈ V are the k eigenvectors of L 1 ≡ I ext and w i ∈ V ⊥ are the (N − 1)k eigenvectors of L 2 and this proves our result.
Remark 4.5. Rescaling appropriately, all of the ferromagnetic and anti-ferromagnetic equilibrium states given above exist on a level set of the summed Casimir C, but not all of them exist if we restrict to a single coadjoint orbit.
As we will see later with the SO(3) example, most of the eigenvalues λ e of L have algebraic multiplicity n λ = mult(λ e ) > 1, so if µ e 1 and µ e 2 are two eigenvectors of L sharing the same eigenvalue λ e , then their linear combination is also an eigenvector and therefore an equilibrium solution. Hence, every point in the eigenspace E(λ e ) := ker(L − λ e 1) is an equilibrium solution with dim(E(λ e )) = n λ (since L is symmetric, the algebraic and geometric multiplicity are equivalent). It is therefore possible to construct more complicated equilibrium states that are neither ferromagnetic nor anti-ferromagnetic by taking a linear combination of a ferromagnetic eigenvector and an anti-ferromagnetic eigenvector that share the same eigenvalue. 4.2.2. Nonlinear stability. We now investigate the stability properties of the equilibrium solutions found above via the energy-Casimir method, as given in [Arn66,HMRW85]. Our main result is stated as follows.
Theorem 4.6. Fixing a level set C −1 (a) for a > 0, the equilibrium solutions µ e of the network Lie-Poisson equation (4.16) corresponding to the highest and lowest eigenvalues λ e of L are nonlinearly stable provided mult(λ e ) = 1.
Proof. Consider an equilibrium solution µ e , which is an eigenvector of L. In accordance with the energy-Casimir method, we first consider the augmented Hamiltonian where Φ is an arbitrary real-valued function such that µ e is a critical point of h Φ . That is, This equation holds for any free variations δµ ∈ g, provided Φ (a) = −λ e , where λ e is the eigenvalue of L corresponding to the eigenvector µ e .
Next, we compute the second variation of h Φ .
Remark 4.7. From (4.20), we see that this corresponds to the highest and lowest energy equilibrium states on the level set C −1 (a).
We cannot deduce further about the nonlinear stability of the other equilibrium states, but as we will see with the SO(3) example in section 5, most of them are linearly unstable and the remaining few are linearly stable at best. To investigate the linear stability of the equilibrium solutions in the general setting, we first linearize the equation (4.16) by setting µ(t) = µ e + δµ(t) and dropping terms of O( 2 ) to geṫ δµ = ad µ e ((L − λ e 1)δµ) , (4.25) where we have used Lµ e = λ e µ e . We can then investigate the linear instability of the equilibrium solution µ e by looking at whether the eigenvalues of the linear operator ad µ e ((L − λ e 1) ·) has a positive real part.
4.3. Noise and dissipation. We now perturb our system (4.16) with noise and dissipation that preserves the Casimir C, following [ADCH18]. The equation is similar to the single coadjoint orbit system with equation (2.8), that is for i = 1, . . . , N , where h = 1 2 µ, Lµ and C is a given Casimir for the Lie-Poisson bracket. The noise term has kN independent Lie algebra vectors σ i,l associated to independent Wiener processes W i,l t and the dissipation considered here is a double bracket dissipation (See [GBH13,BKMR96]). One could generalise this equation further by letting θ i be node dependent, but we will not do this here to keep our equations simple and we also choose our noise to be isotropic, that is, σ i,l = σe l , where e l is a basis vector of g, which are the same at every node.
Notice that without noise, the energy dissipates as and will tend towards the equilibrium configuration with minimum energy on the coadjoint orbit. If the noise is isotropic, the stationary distribution can be computed in the same way as we did with the single orbit case.
Theorem 4.8. The stationary distribution of (4.26) with isotropic noise is the Gibbs distribution where the inverse temperature is given by The partition function is given by is the direct product of the coadjoint orbits at each lattice site, or the total coadjoint orbit of the lattice.
Proof. We refer to [ADCH18] for the proof of this formula, which is done by direct substitution into the Fokker-Planck equation (2.10) of the stochastic process (4.26).
The Gibbs distribution (4.28) provides us with the notion of temperature in this system and as we will see in section 8, we observe a phase transition for the case G = SO(3) as we vary the temperature, similar to the classical Heisenberg model. The mean field approximation can also be derived using the partition function (4.30) and will allow us to help detect the phase transition.

Example I: Networks of rigid bodies
In this section, we study in more detail the case G = SO(3), corresponding to a network of interacting rigid bodies. The corresponding Lie algebra so(3) is compact and semi-simple so we are in the setting of the general theory studied in section 4. 5.1. Equation of motion. Consider a single free rigid body with configuration group SO(3). The reduced Hamiltonian is given by the kinetic energy h(Π) = 1 2 Π · I −1 Π where Π ∈ so * (3) ∼ = R 3 is the angular momentum vector and the corresponding equation, given by the so * (3) Lie-Poisson structure, isΠ = Π × Ω, where Ω := I −1 Π is the angular velocity vector. We will not discuss the dynamics of a single body further here, as it is standard in geometric mechanics and instead refer to chapter 15 of Marsden and Ratiu [MR99] for a more detailed exposition of the system. Extending this to a network (N , SO(3), O) of interacting rigid bodies via momentum coupling with inertia tensor I i and interaction tensor J ij , we get the Hamiltonian we obtain the corresponding network Lie-Poisson equatioṅ One can check that the Casimirs for the Lie-Poisson bracket (5.2) are given by and we denote the sum of the Casimirs by The coadjoint orbit in this special case is given as follows.
is given by the level sets of the Casimirs for some constants c 1 , . . . , c N .
Proof. Since the action of SO(3) on the two-sphere S 2 (viewed as a submanifold in R 3 ) is transitive, the coadjoint orbit is exactly given by the level sets of the Casimirs.
We observe that removing the kinetic energy term in (5.1) and taking d i = d for all i = 1, . . . , N , we recover the Hamiltonian for the Heisenberg model. Hence, the rigid body network can be viewed as a Heisenberg model with mass-carrying spins.

Equilibrium solutions.
We now investigate the relative equilibrium configurations of the deterministic rigid body network and its corresponding stability properties. Since the Lie algebra so(3) is compact and semi-simple, we can apply propositions 4.4 and 4.6 to obtain information about the equilibrium configurations of the rigid body network and its corresponding nonlinear stability properties, which we summarise in the following proposition.  where : (R 3 , ×) → so * (3) is the standard Lie algebra isomorphism that takes a vector in R 3 to a skew-symmetric matrix in so * (3).  Figure 1. In this figure, we plot the maximum real part of the eigenvalues λ S of the linearised system (5.7) corresponding to equilibrium solutions Π e (eigenvectors of L) with energy λ e (eigenvalue of L corresponding to eigenvector Π e ). The blue crosses indicate the maximum real eigenvalues that are non-zero, corresponding to unstable solutions, and the red dots indicate those that are zero, corresponding to linearly stable solutions. The purple dots indicate the algebraic multiplicity n of λ e , whose values are to be read from the right axis. As expected, the equilibria corresponding to the highest and lowest energy are stable in both cases.
In figure 1, we computed the eigenvalues λ s of the linearised system (5.7) plotted against its corresponding energy λ e for two cases: (a) I = diag(1, 1, 1), J = diag(1, 2, 3) and (b) I = diag(1, 2, 3), J = diag (1, 1, 1), where we took a 20 × 20 lattice with periodic boundary condition. As expected from proposition 5.2, we see that in both cases, the highest and lowest energy configurations are linearly stable (both have multiplicity 1), however the majority of the configurations that lie between these two states are unstable. Interestingly, we also see a few linearly stable states towards the high and low end of the energy, with significantly more of them in 1(a) than in 1(b).
We observed numerically that the linearly stable states (red dots) given in figure 1(a), are parallel to the Π 3 -axis, which tend to be characterised by a regular 'argyle' pattern of aligned spins for low energy configurations (figure 2(a)) and anti-aligned spins for high energy configurations (figure 2(c)). As expected, the lowest energy state is a ferromagnetic state that is completely aligned with the Π 3 -axis and the highest energy state is an antiferromagnetic state given by a fine checkerboard pattern. These last two configurations are furthermore nonlinearly stable by proposition 5.2.
The linearly stable states in figure 1(b) can be described as follows, from lowest to highest energy. The lowest energy configurations is a ferromagnetic state that is parallel to the Π 3 -axis as expected, the next (λ = −3.5) is a four dimensional subspace consisiting of a similar 'argyle' pattern of aligned spins seen in the other case (figure 2(d)), the next (λ = −3) is again a ferromagnetic state, along the Π 1 -axis this time, the last few (4 < λ < 5) consist of 'argyle' patterns with anti-aligned spins parallel to the Π 1 -axis 20 × 20 rigid body lattice. All the solutions are parallel to one of the axes (either Π 1 or Π 3 ) and have an 'argyle' pattern of aligned spins (2(a), 2(d)) or anti-aligned spins (2(b), 2(c)). The former pattern corresponds to configurations with low energy (negative values of λ), and the latter corresponds to configurations with high energy (positive values of λ). All four solutions given here correspond to the red dots in figure 1.
( figure 2(b)) and the highest energy state is a completely anti-ferromagnetic checkerboard patterned state as expected.

Network of coadjoint orbits II: Position coupling
In this section, we consider a Lie-Poisson network of coadjoint orbits (N , G, O) that arise from a different approach to coupling the neighbours. We will call this position coupling. In contrast to momentum coupling, where coupling occurs at the level of the reduced space g * , in position coupling, the coupling occurs at the level of the group by considering a representation of G, followed by symmetry reduction, which yields a Euler-Poincaré/Lie-Poisson equation. This introduces a semi-direct product group structure into our system, analogous to the heavy top. We saw earlier that the momentum-coupled equations arise somewhat unnaturally from the Euler-Poincaré framework, however, there is no such issue for the position coupling approach. Furthermore, the equations that we obtain by position coupling are completely new to our knowledge and further investigation is necessary to understand the new phase transition behaviour that we will describe in section 8.
where nodes i and j are adjacent on the graph N and we denote by gA to be the left action of g ∈ G on a ∈ V * . At node i, we assign a Lagrangian L i (g i ,ġ i ; A i ) depending on the parameter A i that satisfies the following properties We take our full Lagrangian to be where g = (g 1 , . . . , g N ) ∈ G ×N . We can check that this Lagrangian is invariant under the (left) diagonal action of G ×N on T G ×N × (V * ) ×N , so we obtain the reduced Lagrangian L(g,ġ, a 0 ) = L(e, g −1ġ , g −1 a 0 ) . . . , ξ N ) and a = (a 1 , . . . , a N ). This allows us to apply the semi-direct product reduction theorem given in Holm et al. [HMR98] to obtain the following network Euler-Poincaré equations.
Theorem 6.1. The following statements are equivalent.
6.1.1. Legendre transformation and Lie-Poisson equation. We find the reduced Hamiltonian by taking the partial Legendre transform in the velocity variable, where µ = (µ 1 , . . . , µ N ) and µ i = δl δξ i . The Euler-Poincaré equation (6.7) then becomes a Lie-Poisson equation which is indeed Lie-Poisson, with bracket given by the (−) Lie-Poisson bracket on the semi-direct product algebra (g * V * ) ×N , which is given by Remark 6.2. Likewise, for a right invariant Lagrangian and right representation, we recover the (+) Lie-Poisson bracket. Again, we refer the readers to Holm et al. [HMR98] for more details about left vs. right group action.
For the special case where g = Lie(G) is compact, semi-simple and given a coadjoint representation of G on V * = g * , the Lie-Poisson structure becomes with the corresponding Lie-Poisson equation (6.14) We refer to Ratiu et al. [Rat81,RM82] for more details on similar systems and in particular their complete integrability in the case of the Lagrange top (which will be lost here, due to the interaction terms).
One can check that the Casimirs for the bracket (6.13) are given by giving rise to the Lie-Poisson structure (6.12). We will also denote by as the sum of the Casimirs.
6.2. Equilibrium positions. We now seek for equilibrium solutions of (6.14) that correspond to the critical points of the Hamiltonian restricted to the level sets of the summed Casimirs C 1 and C 2 . Here, we specialise to Hamiltonians of the form for inertia tensors I i : g → g * at each node i is the kinetic energy, and is the interaction potential energy. Now consider the augmented Hamiltonian h φ,ψ (µ, a) = h(µ, a) + φ(C 1 ) + ψ(C 2 ) , (6.21) where φ, ψ are arbitrary smooth functions and take its first variations, which set to 0 gives the condition for (µ e , a e ) to be an equilibrium solution. One could also choose a more general function of the two Casimirs, but this form turns out to be general enough for our purpose. Taking the variation, we get where we defined λ 1 := φ (c 1 ) and λ 2 := 2ψ (c 2 ) . (6.23) The conditions for the equilibrium solutions are hence given by δµ : I −1 µ e + λ 1 a e = 0 (6.24) δa : −D − 1 2 AD − 1 2 a e + λ 1 µ e + λ 2 a e = 0 , (6.25) and pairing (6.24) with µ e and (6.25) with a e , we find (6.27) Now, substituting (6.24) into (6.25), we obtain the following equation L(λ 1 ) a e = −λ 2 a e , (6.28) where L(λ 1 ) := −λ 2 1 I − D − 1 2 AD − 1 2 is our new extended graph Laplacian. Hence, fixing λ 1 , the equilibrium solutions a e and µ e = −λ 1 I a e again correspond to the kN eigenvectors of the graph Laplacian L(λ 1 ) with eigenvalue −λ 2 . From this, we immediately deduce the following result. The proof is similar to that of proposition 4.4.
Remark 6.4. It is important to note that the ferromagnetic and anti-ferromagnetic equilbrium solutions given above are only found if we fix λ 1 and only one of the summed Casimirs C 1 or C 2 . If instead, we want to find the equilibrium solutions on the level sets of both C 1 and C 2 while keeping the parameter λ 1 free, then one has to solve the full nonlinear equation (6.24), (6.25), (6.26), (6.27) which in general is difficult to solve.
6.3. Nonlinear stability analysis. We saw that the eigenvectors of the graph Laplacian L(λ 1 ) correspond to equilibrium solutions of (6.14), similar to the momentum-coupled case. We will now show that the equilibrium configuration corresponding to the eigenvector of L(λ 1 ) with the lowest eigenvalue is nonlinearly stable. This differs slightly from the momentum coupled case where we were able to prove that both the lowest and highest eigenvalue configurations are stable.
Proof. We apply the energy-Casimir method to assess the stability of the lowest eigenvalue configuration. The Hessian D 2 h φ,ψ (µ e , a e ) of the augmented Hamiltonian is a symmetric matrix given by D 2 h φ,ψ (µ e , a e ) = I −1 +λ 1 a T e a e λ e 1 1 +λ 1 a T e µ e λ e 1 1 +λ 1 µ T e a e λ e 2 1 − D − 1 2 AD − 1 2 +λ 2 a T e a e +λ 1 µ T e µ e , whereλ 1 := φ (c 1 ) andλ 2 := 4ψ (c 2 ). Settingλ 1 = 0, this simplifies to It is well-known that block matrices of this form are positive definite if and only if the upper left block X and the Schur complement B := Z − Y T X −1 Y are both positive definite. Since we take I to be positive definite, it follows that X is positive definite so we only need to show that B is positive definite. Written in full, in terms of the graph Laplacian L(λ e 1 ) = −(λ e 1 ) 2 I − D − 1 2 AD − 1 2 , one can check that δa T B δa = δa T (L(λ e 1 ) + λ e 2 1) δa +λ 2 (a e · δa) 2 .
Since L(λ e 1 ) is symmetric, we can diagonalise it so that L(λ e 1 ) → diag(α 1 , . . . , α kN ) with α 1 ≥ . . . ≥ α kN . Recalling that the equilibrium solution satisfies the eigenvalue problem L(λ e 1 ) a e = −λ e 2 a e , we take −λ e 2 = α kN which is the lowest eigenvalue of L(λ e 1 ) and we assume that mult(−λ e 2 ) = 1. We then have a e = √ c 2 e kN and get where δâ i for i = 1, . . . , kN are the components of δa in this new-basis. From this, it is easy to see that the Schur complement B is positive definite if we chooseλ 2 > 0. Hence, D 2 h φ,ψ (µ e , a e ) becomes positive definite and this configuration is nonlinearly stable, by the energy-Casimir method.
6.4. Noise and dissipation. Following Arnaudon et al. [ADCH18], the general stochastic equations for this system is given by The dissipative terms parametrized by θ are of double bracket form, and preserve the structure of the coadjoint orbit. Their complicated form is due to the semi-direct product structure, see [GBH13, GBH14, ADCH18] for more details. The important difference here is in the χ(a) term which contains interactions between the neighbouring spins on the network and appears in the dissipative terms. These terms are crucial for the existence of the Gibbs distribution (4.28).

Example II: Networks of heavy tops
We now consider the case G = SO(3), and take the coadjoint representation of SO(3) on so * (3) ∼ = R 3 . The neighbours are coupled by the orientation of a fixed vector Γ 0 ∈ R 3 rotated around by the SO(3) action. This breaks the symmetry in our Lagrangian so we extend our configuration manifold to SO(3) × R 3 such that the extended Lagrangian is invariant under the diagonal SO(3) action on SO(3) × R 3 . The Lie-Poisson equation that we obtain via symmetry reduction has the Lie-Poisson structure of the semi-direct product group SE(3) ∼ = SO(3) R 3 , hence we call this system the "heavy top network" as opposed to the rigid body network obtained via momentum coupling.

Equations of motion. Starting with the full Lagrangian
(interpreted as matrix multiplication), we get the reduced Hamiltonian where Π = (Π 1 , . . . , Π N ), Γ = (Γ 1 , . . . , Γ N ), Π i = R −1 iṘ i and Γ i = R −1 i Γ i 0 . Taking the semi-direct product Lie-Poisson structure we get the Lie-Poisson equatioṅ where Ω i := I −1 i Π i is the angular velocity at node i. One can check that the Casimirs for this bracket are given by so in particular, the sums are conserved by the dynamics. Now, the coadjoint orbits of the heavy top network are given as follows.
Theorem 7.1. The coadjoint orbits O = O 1 × · · · × O N of the heavy top network are given by, if c i,2 = 0, which is a four-dimensional submanifold and if c i,2 = 0, which is a two-dimensional submanifold, unless Π i 2 = 0.
Proof. We refer to theorem 1.2 in [Rat81] for the proof in the single body case. The general multi-body case considered here is an easy extension.
7.2. Noise and dissipation. From the general equation (6.29), we use C i,1 = Π i · Γ i as the Casimir for the double bracket dissipation to obtain the stochastic equation Notice that the other Casimir C i,2 is also preserved by this dissipation, thus we did not include it here. However, it is possible that including the Casimir C i,2 would change the behaviour of the system, but we leave this investigation for future work.
7.3. Equilibrium solutions. We now study the equilibrium solutions of the heavy top network and its corresponding stability properties. The classification of equilibrium solutions into ferromagnetic and anti-ferromagnetic states and its nonlinear stability property are summarised in the following proposition.
In order to investigate the stability of the other equilibrium configurations, we move to linear stability analysis. 7.3.1. Linear stability. Linearising equations (7.4) by taking Π(t) = Π e + δΠ(t) and Γ(t) = Γ e + δΓ(t) and dropping terms of O( 2 ), we get in term of Γ e only, d dt where Π e := diag( Π e 1 , . . . , Π e N ), Γ e := diag( Γ e 1 , . . . , Γ e N ) and : (R 3 , ×) → so(3) is the standard Lie algebra isomorphism that sends a vector in R 3 to a skew symmetric matrix in so(3). Hence, we can assess the linear stability of an equilibrium solution (Π e , Γ e ) by looking at the eigenvalues of the matrix on the right hand side of (7.11).  Determining the linear stability of all the equilibrium states of this system is impossible to do analytically, so we will only discuss the numerical results here. First, we solve (6.28) to find all the equilibria of the lattice given a λ 1 , and in figure 3, we plotted the values of all possible λ 2 (i.e. minus the eigenvalues of L(λ 1 )) as a function of λ 1 . In the case I = diag(1, 2, 3) and J = diag(1, 1, 1) ( figure 3(a)), we see that as λ 1 → 0, all the solutions become degenerate with multiplicity 3 (i.e. the eigenvalues λ 2 collapse) but this does not happen in the other case I = diag(1, 1, 1) and J = diag(1, 2, 3) ( figure 3(b)). This can be explained by the fact that at λ 1 = 0, which corresponds to the zero momentum casē Π = 0, the system is isotropic in case (a) and anisotropic in case (b). . The two panels in this figure display the maximum, real part of the eigenvalues λ S of the linearised system (7.11) around the equilibrium configuration Γ e (eigenvectors of L(λ 1 )) plotted against λ 2 (minus the eigenvalue of L(λ 1 ) corresponding to Γ e ) for λ 1 = 0.5. Again, the blue crosses correspond to unstable equilibria, the red dots correspond to linearly stable equilibria and the purple dots correspond to the multiplicity of λ 2 , whose values are to be read from the right axis. As expected, the configuration corresponding to the highest λ 2 (or, the lowest eigenvalue of L(λ 1 )) is stable in both cases.
In figure 4, we display the largest, real part of the eigenvalues λ s of the linearised system (7.11) plotted against λ 2 (minus the eigenvalue of L(λ 1 )) corresponding to equilibrium solutions Γ e of (7.11) for λ 1 = 0.5. As expected from proposition 7.2, the lowest eigenvalue state (or, the states with highest λ 2 ) in both cases 4(a) and 4(b) are linearly stable. This corresponds to a ferromagnetic equilibria along the Γ 3 -axis in both cases. For I = diag(1, 2, 3) and J = diag(1, 1, 1) ( figure 4(a)), the two states near this configuration with multiplicity 1 (λ 2 ≈ 4.2 and λ 2 ≈ 4.5) correspond to the other two ferromagnetic equilibria given in proposition 7.2, which we call states A and B and these are seen to be unstable. There also exists a subspace of anti-ferromagnetic equilibria between these two states (λ 2 ≈ 4.35) with multiplicity 4 that have a very small but positive λ s . We call this state C. As we took λ 1 → 0, we saw that state C becomes linearly stable for some small λ 1 , and states A and B, while they were always found to be unstable, their corresponding values of λ s tended to 0 (i.e. approaching a linearly stable state). In the other case I = diag(1, 1, 1) and J = diag(1, 2, 3) ( figure 4(b)), the linearly stable equilibria excluding the nonlinearly stable state were found to be anti-ferromagnetic states with the 'argyle' pattern similar to that in figure 2. 7.4. Numerical experiments. Before looking at phase transitions, we will present here an interesting phenomenon that arises in the heavy top lattice when I = diag(1, 2, 3) and J = diag (1, 1, 1). In figure 5, we show several simulations of the deterministic 20×20 heavy top lattice starting from different initial conditions, with or without the double bracket dissipation. The initial conditions are taken to be nearly ferromagnetic, with spins Γ i aligned to a fixed direction, except for a small perturbation at two nodes. This way, we can numerically assess the stability of the ferromagnetic equilibria. From proposition 7.2, we know that position (0, 0, 1) (i.e. the Γ 3 -axis) is nonlinearly stable, and we also observed this in our simulations. Hence, we will not display the corresponding plots here as it is not very interesting. Instead, we will display the simulations starting close to the two other ferromagnetic equilibria (0, 1, 0) in figure 5(a) and (1, 0, 0) in figure 5(c), which were seen to be unstable with or without dissipation.
There are two interesting behaviors of the system observed from the plots. The first, which is displayed in panel 5(c), is that the solution starting close to the shortest axis (1, 0, 0) with double bracket dissipation is stuck in this position for a while, then gets stuck in the middle axis position (0, 1, 0) briefly, before relaxing to the lowest energy state (0, 0, 1). This can also be seen in panel 5(a) where the solution starting close to the middle axis (0, 1, 0) is stuck there for a while before relaxing to the lowest energy state (0, 0, 1). Although the double bracket dissipation can never stabilize an unstable equilibrium solution, as proven in Bloch et al. [BKMR94], this indicates some kind of transient stability, or metastability of the shortest and middle axis equilibrium, which will be observed again in the phase transition plots in section 8. We note that this phenomenon is not apparent in the absence of dissipation and is only observed clearly in the presence of dissipation. Apart from these special cases, for instance when we start from (1, 1, 1), which is not an equilibrium, the dissipation will always drive the system towards the lowest energy position without getting stuck (see figure 5(e)).
The second observation that we made is the partial synchronisation phenomenon seen in the non-dissipative simulations (see the bold lines in 5(a), 5(c) and 5(e)). The corresponding energy (total, kinetic and potential) plots are given by the bold lines in 5(b), 5(d) and 5(f) to demonstrate the validity of the simulations. In all of the cases considered, we observe that after some time, the spins relax to a state where on average, it oscillates closely around the lowest energy configuration (i.e. ferromagnetic state along the Γ 3 -axis), despite the absence of dissipation. We call this phenomenon "partial synchronisation". This result is rather surprising as, by Liouville's theorem, the volume of the phase space is conserved by the dynamics, so by the Poincaré recurrence theorem, any states starting near the equilibrium should return sufficiently close to it after a finite time. However, in our simulations, the trajectories after some time seem to get stuck in another part of the phase space without returning to a region close to the initial conditions. This is very counter-intuitive from the original dynamics of the rigid body, where an initial condition near the unstable (middle) axis will eventually come back near it, after traversing a long trajectory on the momentum sphere. One possible avenue of investigation of this phenomenon is to study the local dynamics of the lattice where almost periodic motions were , 5(f)) in a 20 × 20 heavy top lattice, starting from different initial conditions, with dissipation (dashed lines) or without dissipation (bold lines). The blue, orange and green lines in the left panel correspond to the Γ 1 , Γ 2 and Γ 3 -components of the spins averaged over the lattice respectively and the same colors in the right panel correspond to the total, kinetic and potential energy. We observe metastability of the Γ 1 and Γ 2 axes in panels 5(c) and 5(a) in the dissipative case, and a partial synchronisation around the Γ 3 -axis in the non-dissipative case in 5(a), 5(c) and 5(e).
observed (see the videos in http://wwwf.imperial.ac.uk/~aa10213/). This phenomenon may also be tied to the complex interaction between the Γ and Π variables and the non-compactness of the phase space. However, we will not study this further in the present work and instead leave it for future investigations.

Temperature phase transitions
In this last section, we will study phase transitions in the two examples that we constructed above, namely, the rigid body network and the heavy top network. There exist various types of phase transitions, but here we will focus on second-order phase transitions that arise from varying the temperature of the system. This is characterized by a transition from an orderly state of the spins, measured by how much they are aligned with each another, to a completely disordered state. In order to detect a phase transition, one can either apply a mean field approximation of the model and try to solve it analytically or perform direct numerical simulations of the underlying dynamics of the lattice. Many other methods are available but are out of the scope of this first investigation.
8.1. Mean field approximation. The mean field approximation relies on the assumptions that (1) the microscopic system is identical at each node and (2) the spins (momentum or position) of the rigid bodies in statistical equilibrium are close to its mean. This approximation is increasingly accurate if each site has more neighbours, which is the case for high dimensional lattices. In two dimensions, the approximation fails to properly assess the critical temperature and the corresponding critical exponents, but still, give a good indication of the presence of a phase transition. 8.1.1. Mean field approximation of the rigid body network. We assume that the system is identitcal at each node, i.e. I i = I for all i = 1, . . . , N , and that the interactions between the neighbours are identical, that is, J ij = J for all i, j = 1, . . . , N for the mean field approximation to be valid. Now, define the averaged momentum where β = 2θ σ 2 = T −1 is the inverse temperature, P ∞ (·) dΠ is the Gibbs measure (4.28) and O = O 1 × · · · × O N is the total coadjoint orbit. Since the system is assumed to be identical at each node, we have O i = O for all i = 1, . . . , N .
Since we assume that the spins Π i = π + δΠ i are close to the mean π , we linearise the interaction Hamiltonian around Π to get where the first term in the second line is neglected since it is quadratic in δΠ i and is therefore very small, and the last term is also neglected since it is constant and does not contribute to the dynamics. We therefore obtain the complete mean field Hamiltonian Notice that the kinetic energy is still exact. From this Hamiltonian, the Gibbs distribution and the partition function take a simpler form, and one can check that the average momentum Π simplifies to which is now an implicit equation for the order parameter Π . This equation is difficult to solve analytically, as it involves integrals over the momentum sphere but can be numerically estimated using Monte-Carlo integration. We will display the numerical approximation of Π in the next section, compared to the full simulations.
8.1.2. Mean field approximation of the heavy top network. In a similar fashion, we can derive the mean field approximation of the heavy top network. In this case, we take our order parameter to be the averaged position, defined by where O 1 = O 1,1 ×· · ·×O N,1 and O 2 = O 1,2 ×· · ·×O N,2 are coadjoint orbits corresponding to level sets of the Casimirs C i,1 and C i,2 respectively. Again, since we assume the system to be identical at each node, we take O i,1 = O 1 and O i,2 = O 2 for all i = 1, . . . , N .
The partition function in the mean field approximation is found to be 8.2. Numerical simulations. The simplest (but computationally expensive) way to detect phase transitions in lattices is by a direct simulation of the full stochastic equation, sampling from the Gibbs measure. This is equivalent to a classical Monte-Carlo simulation of the Ising model for example but in the more general setting of continuous spins on coadjoint orbits. The phase transitions we observe are all second-order and is detected as we increase the temperature T = σ 2 2θ . The order parameter for the rigid body and the heavy top network are given by the averaged momentum and position respectively. We summarise here the phase transition behaviours of the four systems we simulated, displayed in figures 6 and 7.
(1) Rigid body with I = diag(1, 1, 1) and J = diag(1, 2, 3) (Figure 6(a)) This case corresponds to the classical Heisenberg model with anisotropic interactions (also known as XY Z-model). Notice that some simulations have anomalous low magnetisation at low temperatures. This can be explained by the fact that in some simulations, the system got stuck in a state consisting of large regions of opposite spins. A longer simulation would be necessary to see these two domains merge into one to reach the minimum energy state.
(2) Rigid body with I = diag(1, 2, 3) and J = diag(1, 1, 1). (Figure 6(b)) This case corresponds to a massive isotropic Heisenberg model with non-uniform mass. Our simulations did not get stuck in the state consisting of opposite spins, as seen in the previous case, and the mean field approximation is closer to the direct simulations. The reason that the mean field approximation is more precise in this case can be explained by the fact the interaction term, which is approximated, is isotropic but the kinetic term, which is exact, is anisotropic.
(3) Heavy top with I = diag(1, 1, 1) and J = diag(1, 2, 3). (Figure 7(a)) This case also corresponds to the classical Heisenberg model, at least from the mean field point of view (by remark 8.1), even if the full dynamics is different. Again, we observe solutions getting stuck in the state consisting of opposite spins and hence the existence of anomalies. (4) Heavy top with I = diag(1, 2, 3) and J = diag(1, 1, 1). (Figures 7(b),7(c) and 7(d)) This last case is the most interesting as it shows a more complex phase transition, which we will describe in details below.
8.2.1. Triple-humped phase transition in the heavy top network. We now discuss the last case in more details. From the direct numerical simulations, we observed a phase transition from a strongly magnetised state (large value of Γ ) along the Γ 3 -axis to a nonmagnetized state (that is, Γ = 0) as we increased the temperature. But between these two states, we also observed two intermediate phase transitions, where magnetisation along the other two axes also occur before becoming completely disordered. We call this a 'trimple-humped' phase transition. These intermediate phase transitions indicate that these unstable ferromagnetic equilibria along the Γ 1 and Γ 2 axes can still support magnetisation, which, in statistical physics is generically called a meta-stable state. We note that this phenomenon is not captured in our mean-field simulations.
From our linear stability analysis in section 7, we observed that for small values of λ 1 , equivalent to a small ratio c 1 /c 2 , the ferromagnetic equilibria along the Γ 1 and Γ 2 axes are unstable, but are close to being linearly stable. This is compatible with the observation that these intermediate phase transitions exist only for small values of c 1 /c 2 . For example, in panel 7(d) we took c 1 /c 2 = 1.15 and saw that the third phase transition is almost negligible. Increasing the value of this ratio further will also remove the other intermediate phase transition. On the other hand, decreasing c 1 /c 2 , the intermediate phase transitions persist and furthermore, will shift their respective critical temperature towards zero. In figure 7(c), where we took c 1 /c 2 = 0.8, we see that both critical temperatures T 1 and T 2 are smaller compared to those in figure 7(b), where we took c 1 /c 2 = 1.
We will not investigate this phase transition further here, and leave its mathematical understanding as a challenging open problem.

Conclusion and outlook
In this paper, we established a link between geometric mechanics and statistical mechanics by constructing a network of interacting Lie-Poisson system on g * with noise and dissipation that preserves the coadjoint orbits, which gave us a canonical ensemble for the system in statistical equilibrium. For the construction of the system, we considered two types of coupling, one where the neighbours are coupled in the reduced space and the other where the neighbours are coupled directly on the configuration group by considering a representation of the group on a given vector space. The first approach yielded a direct generalisation of the classical Heisenberg model to include general symmetry groups with an additional kinetic energy term and the second approach gave a system that is possibly new. In the special case where g is compact and semi-simple with Hamiltonian of the form kinetic + potential energy, we were able to find the equilibrium solutions of the purely deterministic system as the eigenvectors of the underlying extended graph Laplacian and In 7(b)-7(d), we also varied the Casimirs (where c 1 = Π · Γ and c 1 = Γ 2 ) to show that the three critical temperatures depend on their values, and may even disappear for large enough c 1 /c 2 . Due to the high dimension and non-compactness of the phase space, the mean field approximation requires intensive computation to be precise, explaining the visible errors, especially at low temperatures.
found that (1) for the momentum-coupled case, the equilibrium solution corresponding to the lowest and highest eigenvalue of the graph Laplacian are nonlinear stable and (2) for the position-coupled case, the equilibrium solution corresponding to the lowest eigenvalue of the graph Laplacian is nonlinearly stable. Furthermore, we showed that in both cases, these equilibrium solutions can be classified into ferromagnetic and anti-ferromagnetic states. In our numerical simulation of the rigid body lattice and the heavy top lattice, which are the simplest examples of momentum-coupled and position-coupled systems respectively, we observed a second order phase transition, similar to that in the Ising model or in the Heisenberg model. However, in the heavy top network, we also observed a 'triple-humped' phase transition, in which the system underwent two intermediate phase transitions before settling down to the lowest energy configuration as we decreased the temperature, which is unusual for simple lattice models.
In future work, we would like to investigate further this new type of phase transition behaviour that we observed for the heavy top network, which we believe is related to the metastability of the intermediate ferromagnetic states. However, even without noise, we numerically observed unusual behaviour of the heavy top network, such as when the spins start close to an arbitrary ferromagnetic state, there is an exchange between the kinetic and potential energy that causes the spins to relax to a state that oscillates closely around the stable ferromagnetic state with lowest potential energy, despite the absence of dissipation. So studying the deterministic heavy top network further could also be interesting and may help us understand this phase transition behaviour better. This partial synchronisation result deserves a more detailed study, in particular for more general networks, where synchronisation of oscillators are an active subject of research (see for example [BP02] and the many subsequent works). Other interesting phenomena can be observed for Heisenberg models on certain types of networks, such as supra-oscillations (see for example [EdNTL17] and references therein).
Regarding phase transitions, one could also compute the critical exponents and the corresponding universality class of the phase transition seen here, or even doing a more thorough analysis using Landau theory to better understand the dynamics near the critical temperature. We can also investigate different types of phase transitions, for instance by varying the external magnetic field instead of temperature, or even extending our domain from a simple lattice to a general network. One may also guess that certain topological phase transitions could even be observed in these systems, such as the Kosterlitz-Thouless transition [KT73].
As we can see, there are many interesting questions that are open for further investigation which we hope to address in future works.