Hamiltonian polar particles at finite temperature:
equilibrium orderbydisorder transition to collective motion.
Abstract
Using analytic and numerical methods, we study a Hamiltonian model of interacting particles carrying ferromagnetically coupled continuous spins which are also locally coupled to their own velocities. This model has been characterized at the mean field level in a parent paper. Here, we first obtain its finite size ground states, as a function of the spinvelocity coupling intensity and system size, with numerical techniques. These ground states, namely a collectively moving polar state of aligned spins, and two non moving states embedded with topological defects, are recovered from the analysis of the continuum limit theory and simple energetic arguments that allow us to predict their domains of existence in the space of control parameters. Next, the finite temperature regime is investigated numerically. In some specific range of the control parameters, the magnetization presents a maximum at a finite temperature. This peculiar behaviour, akin to an orderbydisorder transition, is explained by the examination of the free energy of the system and the metastability of the states of minimal energy. The robustness of our results is checked against the geometry of the boundary conditions and the dimensionality of space.
I Introduction
Collective motion, the macroscopic ordering of velocities in a manybody system, is a key feature of assemblies of selfpropelled particles. This phenomenon has recently drawn a lot of attention from the statistical physics community Ramaswamy (2010); Marchetti et al. (2013), and has been reported and extensively studied in a variety of systems, both natural, like flocks of birds Cavagna et al. (2010) or bacterial colonies, Zhang et al. (2010) and artificial, like active colloids Bricard et al. (2013) or vibrated grains Deseigne et al. (2010). Interacting systems of selfpropelled particles have also been the focus of many theoretical studies, starting with the celebrated Vicsek model Csahok and Vicsek (1995), and branching into many more recent analytical and numerical works Bechinger et al. (2016).
In typical Visceklike models moving individuals are represented as polar selfpropelled particles, i.e., particles set in motion by a driving force directed along the heading vector of the particle. Assuming, that the speed is always constant leads to minimal spin models where particles move along the direction of their spin and interact with their neighbours so as to align their velocity vectors, giving rise to collective motion. An important simplification of Vicsektype models is to identify the direction of motion with that of the spin of the particles. However, as soon as the particles experience other types of interactions, as simple as rigid body repulsion, it is a priori relevant to consider the spins and velocities as distinct dynamical variables. It can also be the case that the dynamics of the spin is coupled to that of the velocity, in such a way that, not only the particle aligns its velocity on its spin, but also it aligns its spins on its velocity. This has been exemplified as the cause for collective motion in the case of selfpropelled granular disks Weber et al. (2013); Lam et al. (2015), for which the steric interaction alone does not ensure alignment of the velocities, and more recently in the effective dynamics of topological defects in active nematics Shankar et al. (2018). It was also demonstrated as being a necessary condition for describing the dynamics of a motorized polar particle in a harmonic potential Dauchot and Démery (2019).
Coupling spin to velocity is expected to induce new physics, not only in the realm of active matter, but also in the context of usual Hamiltonian dynamics. Such couplings break the Galilean invariance. Broken Galilean invariance is also a hallmark of active systems. One may thus wonder if some of the remarkable features of active matter could also occur in conservative systems. Such questions have been recently discussed in the context of a conservative twodimensional particle model Bore et al. (2016), in which particles carry a continuous and classical spin. The model, defined in Eqs. (1) and (2), is characterised by ferromagnetic interactions between spins of two different particles, and a ferromagnetic coupling between the spin and velocity of the same particle. Because of this coupling, not only Galilean invariance is broken, but also the conserved linear momentum associated to translation invariance is not proportional to the velocity of the center of mass and the dynamics are not invariant under a global rotation of the spins alone. This, in principle, leaves room for collective motion. A preliminary study of the model has shown that a transition to collective motion does indeed takes place in the meanfield (fullyconnected) limit, in spite of momentum conservation Bore et al. (2016).
In the present paper, using molecular dynamics simulations and analytic calculations, we characterize the behaviour of this model at homogeneous densities, focusing on the effects of the size, dimensionality, temperature, and spinvelocity coupling intensity. We report a finitesize magnetization crossover that is accompanied by a collective ordering of the velocities at low enough temperatures and low enough values of the spinvelocity coupling constant, hereafter defined as . This ordering breaks down at a finite value of , that depends on the size of the system, as accurately predicted using simple energetic arguments. The resulting nonmagnetized ground states feature topological defects that, by suppressing the magnetization, enable the system to avoid a high kinetic energy cost. The crossover between magnetized and nonmagnetized states resembles a firstorder phase transition, as both families of states are metastable close to the critical value of the spinvelocity coupling. When heated up, systems prepared in a nonmagnetized ground state display a maximum of the modulus of their magnetization at a finite temperature. This feature is reminiscent of an “Order by Disorder” Villain et al. (1980); Shender and Holdsworth (1996) (ObD) transition, in which a system develops a spontaneous magnetization upon heating, and can be explained by the examination of the free energies of the metastable states of the system. Finally, the magnetized states are characterized by a collective ordering of the instantaneous velocities, which is however rapidly destroyed by thermal fluctuations. The above results are rather general. On one hand the choice of the spinvelocity coupling is highly constrained by the structure of the Hamiltonian dynamics. On the other hand neither the geometry nor the dimensionality seem to play a major role in our observations.
The paper is organised as follows. In Sec. II we introduce the model, emphasise the unusual features of its conservative dynamics, and describe the numerical simulations used thereafter. In Sec. III we characterise the lowtemperature equilibrium states and predict their location in parameter space with simple energetic arguments, as well as their nature, resorting to a coarsegrained continuous theory. In Sec. IV we consider the finite temperature regimes and describe an ObD phenomenon for a subset of values of the parameters. The crossover towards order is studied, and proved to behave much like a firstorder phase transition. In Sec. V we discuss the alignment of velocities in collectively moving phases. Finally, we discuss the generality of our results, show that the choice of a spinvelocity coupling is highly constrained and present our conclusions regarding the link between this model and usual collective motion in Sec VI.
Ii The Model
The model, introduced in Ref. [Bore et al., 2016], describes the dynamics of interacting particles that carry continuous planar spins of norm equal to . The particles are confined to move in , in a periodic square box of linear size . Their motion is described by the Lagrangian
(1) 
that depends on position and spin variables. The position of the th particle is denoted and its velocity , while is the distance between the centers of particles and ; is the angle that the spin forms with a reference axis and fully parametrizes the continuous spin of unit modulus. The time derivative of the spin vector is indicated by ; is the angle between the spin of particle and the one of particle . Furthermore, is the mass of each particle and its moment of inertia. is a shortranged, isotropic and purely repulsive twobody interaction potential and is a shortranged and isotropic ferromagnetic coupling between the spins. We characterise these two potentials below. For future convenience we made explicit the typical amplitudes of the twobody potential, , and the ferromagnetic coupling, . is the parameter that controls the strength of the spinvelocity coupling, the term that yields the special properties to the model. Throughout this paper, we will only consider the case as the case can be recovered by simply flipping all the spins. In the following, we will introduce dimensionless variables according to the transformations , , , , and, in the notation hereafter, we absorb in the definition of the twobody potential .
In the Hamiltonian formalism the Lagrangian (1) transforms into the Hamiltonian
(2) 
where we defined the canonical momenta and . As already stressed in Ref. [Bore et al., 2016], this specific form of the momentum turns into a direct relation between the velocity of the center of mass and the total magnetization in the microcanonical ensemble. In particular, setting the total momentum to zero, the system spontaneously develops collective motion, whenever it magnetizes. Its center of mass velocity then satisfies where
(3) 
are the intensive magnetization and the velocity of the center of mass, respectively.
Note that although free particles would tend to align their velocity on their own spin Bore et al. (2016), the collectively moving states at low energy that satisfy antialign the particle velocities with their spins.
Also because the velocity of a magnetized state grows linearly with , low potential energy moving states become very expensive in terms of kinetic energy as increases.
As a result, for high enough values of nonmagnetized states featuring topological defects were suggested as candidate ground states Bore et al. (2016).
We study this hypothesis in detail in the body of the paper and will indeed validate it.
The soft interaction potentials are given by
(4)  
(5) 
where is a Heaviside step function, is the range of the interactions, that we
fixed to 1, and we took so that both potentials are equal at half range.
The numerical analysis is carried out using “microcanonical” Molecular Dynamics (MD) simulations of the dynamics defined by the Hamiltonian in Eq. (2), which conserves energy and momentum, and we restrict ourselves to initial conditions such that . We use random initial states with uniformly distributed and drawn from centered reduced Gaussian distributions. Such configurations are placed in a square box with periodic boundary conditions and, after giving some time for the dynamics to settle in, they are subjected to either a numerical annealing or a highrate quench. Further details on the simulation and numerical methods are provided in App. A. The phase diagram of this model has been fully characterized in the case [Casiulis et al., 2019]. Because aligned spins also attract each other, the spin fluid undergoes a liquidgas phase separation at sufficiently low temperature and density Casiulis et al. (2019). In the present paper we concentrate on the physics of homogeneous phases and therefore set the density to (or, equivalently, in terms of the packing fraction defined as ).
Iii FiniteSize Ground States
iii.1 The parameter space
In order to fully understand the conditions for the observation of moving states or topological defects, we start by focusing on equilibrium states at very low temperatures. As shown in Fig. 1, we find three kinds of ground states depending on the values of and :

Polar states: These are fully magnetized states with , with collective motion taking place in the direction opposite to the magnetization as imposed by the momentum conservation, .

Solitonic states: These are nonmoving states with , characterized by a continuous rotation of spins in one direction.

Vortex states: These are nonmoving states with , characterized by 4 topological point defects, 2 vortices and 2 antivortices.
Importantly enough, neither local nor global ordering of velocities is observed in the solitonic and vortex states. The reason for this is the high kinetic energy cost of the magnetized states, which increases with and like . As these parameters are increased starting from a polar state, the system produces topological defects as a way of setting the global magnetization to zero and therefore escape the collective motion imposed by momentum conservation. The strong frustration induced by the high kinetic energy cost associated to magnetized states with a finite velocity of the center of mass has thus the effect of generating topological defects in the equilibrium ground states. Note that these defects are not thermal excitations and should not be confused with the vortices observed at the BKT transition in lattice models of XY spins. Kosterlitz and Thouless (1973); Tobochnik and Chester (1979) In fact we have shown in Ref. [Casiulis et al., 2019] that no BKT transition occurs in the limit of the model (2).
iii.2 Continuum theory
The nature of the observed ground states, can be further confirmed by considering the solutions of a coarsegrained continuum theory for the magnetic properties of the model. We define , and , the kinetic, spinvelocity and interaction contributions to the Lagrangian. They read
(6)  
(7)  
(8) 
respectively, with and . The symbol indicates a double sum over and from 1 to with no repeated indices included. Focusing first on the interaction term , we assume that the spins are locally well aligned, so that , where , in the neighbourhood defined by the range of . Therefore,
(9) 
The main effect of the first two terms is to ensure the uniformity of the density profile . Accordingly, we will discard them in the following and concentrate on the last term, that we henceforth call , and want to approximate using a continuum description. To do so, we assume that there exists a discrete mesh of allowed positions for the particles, Dean (1996) with arbitrarily small spacing, so that we may define coarsegrained fields on this mesh as:
(10) 
with being Kronecker symbols. First, rewriting the sum over after the transformation , and using the rotational and translation invariances of ; second expanding the difference of to the first order, we obtain
(11) 
We now take the continuum limit by taking the mesh size to zero and find that in the absence of singular values of :
(12) 
For , with a step function, and recalling that , we finally obtain, after integration of the second integral:
(13) 
The two other contributions to the Lagrangian, namely, and are local, so that the same procedure as above can be trivially applied. Finally, we can write a continuum Lagrangian theory of the form
(14) 
where with the number of particles in one interaction volume, and is a unit vector pointing in the direction given by . This form of the Lagrangian clearly shows that plays the role of a field, and creates a nonperturbative term even in the lowtemperature expansion which otherwise gives a free theory for the XY model at low temperatures. Kosterlitz and Thouless (1973)
We can now write the EulerLagrange equations for and to describe the lowenergy regime of this theory,
(15)  
(16) 
where the additional comes from the derivative of with respect to . Interestingly, at this level of description, Eq. (15) selects solutions such that the local momentum is zero. When injecting this zeromomentum condition into the equation on , we simply recover a wave equation on : this is the usual spinwave regime Kosterlitz and Thouless (1973) at low temperatures, that describes magnetized lowtemperature states in finite size. To describe states that are close in energy to a finitesize polar moving state, we set , and corresponding to a small perturbation around a polar state aligned along , say . The Lagrangian evaluated for such states reads,
(17) 
and contains a field term that is independent of the amplitude of the velocity perturbation. In particular, it yields a term of order one, contrary to the other terms in the Lagrangian that are vanishingly small for . Therefore, close to lowtemperature polar ground states, the effective theory describing spin alignment is essentially a sineGordon theory with a field of amplitude .
The sineGordon theory is known to feature solitonic solutions and vortex patterns, as reported in previous numerical worksBarone et al. (1971); Gouvea et al. (1990). We shall here briefly discuss both families of states. Solitons are propagative solutions of the sineGordon theory whose functional form can be shown to be Barone et al. (1971)
(18) 
where is the orientation of the local magnetization field, is the usual variable of propagative solutions, with a propagation velocity , is an offset for this variable that sets the position of the center of the soliton, is the socalled polarity of the soliton, is a Lorentz factor associated to the propagation velocity, and is the inplane magnetic field of the sineGordon theory. This function represents a localized rotation of of the magnetization field along the direction. Noticing that in simulations we only ever see one such rotation over the size of the whole box, we plot in Fig. 2 the magnetic field with a constant amplitude and an orientation defined by the following function of position ,
(19) 
where is the range of the box in which we plot the field, that is centered on . This field is a perfect reproduction of the arched states that we observe at low temperatures in the numerical simulations, as shown in Fig. 2, thus justifying the name “solitons” that we give them. As there is neither damping nor forcing in the present model, the soliton propagation velocity is a priori free to select any value. Malomed (1991) In practice, this velocity is found to be very low in our simulations, but nonzero, as demonstrated by the video shown in the Supp. Mat., which exhibits an angular velocity .
The sineGordon equation also features vortex solutions, but their exact functional shape is in general hard to derive. It is customary to consider vortex structures that are solutions of the Laplace equation, and to assume that they are little deformed by the addition of the nonlinear term in the sineGordon equation Shagalov (1992); Sinitsyn et al. (2004). Onedefect solutions can then be written as
(20) 
where are the coordinates of the center of the defect, is a parameter that tunes the twist of the spiral surrounding the defect, and the sign in front of the determines whether the defect is a vortex () or an antivortex (). An example of a 4vortex configuration obtained with this functional form is shown in Fig. 2, and compared to a 4vortex numerical ground state shown in Fig. 2. In this example, we chose the positions of the centers of the defects as well as the value of to best fit the numerical observations.
Altogether the low energy excitations of the continuous theory perfectly match the very low temperature states obtained numerically, confirming the intuition of the mechanism whereby the system escapes the kinetic energy cost of the moving polar state, when or are too large.
iii.3 Energetic arguments
One can further estimate the energy of these states and thereby obtain the range of parameters in which each state is observed. The total kinetic and magnetic energy per particle in a moving polar state is
(21) 
where the velocity is everywhere assumed to be , is the number of neighbours of each particle, and is the mean value of the ferromagnetic coupling for the interparticle distances in the considered configuration. The energy contribution due to the repulsive interaction is likely to be very similar in the three different states of interest and therefore does not need to be evaluated.
The energy of a solitonic state can be approximated as follows. Let us consider that the rotation is smooth, so that two neighbours along the direction of the rotation (on average half the neighbours) are misaligned by an angle , while spins are perfectly aligned along the other direction. Then, the total energy of a solitonic state can be approximated by
(22) 
In both cases, we set the number of magnetic neighbours to 6 (in agreement with the numerical observations at the density used in the simulations), and compute the typical value of the ferromagnetic coupling by taking . For a triangular lattice one has that .
Equating now the two energies and yields an equation for the line separating polar and soliton ground states,
(23) 
which is the lower dashed black line reported in Fig. 1. It captures remarkably well the crossover between the two kinds of states. Note that approximating the cosine by its Taylor expansion up to second order, yields , and stresses that and play a similar role in breaking the magnetic order.
To predict the curve that separates the solitonic and vortex states, we use an exact expressions for the excess energy of a vortexantivortex pair compared to a fully magnetized configuration in an onlattice XY model Kosterlitz and Thouless (1973). Taking into account only nearestneighbour pairs of defects, this excess energy can be written as
(24) 
where is the distance between the two defects and is the size of the center of the defect. Here and the vortex core diameter is of the order of the particleparticle spacing . We then equate to in its Taylorexpanded form and find a critical value for the number of particles above which four vortices states are more favorable than the solitonic ones:
(25) 
This corresponds to the vertical line in Fig. 1. While it does correspond to the size above which we mostly observe vortices, we note that also plays a role in the change between soliton and vortex states that is not explained by this crude estimate of the energy. A way to take into account the effect of on the vortex state energy is to reckon that plays the same role as an inplane field in the magnetic model analog, and to use the expression for the energy of a vortexantivortex pair in a field: Gouvea et al. (1990)
(26) 
In the magnetic analogue is the amplitude of an external inplane field, while here it is an effective field due to the spinvelocity coupling . According to the effective sineGordon theory (see Sec. III.2), . An alternative and more intuitive manner to obtain the quadratic dependence on is to recall that Eq. (26) gives the excess energy of vortexantivortex pairs with respect to a magnetized configuration under the same field. As a consequence, both energies should in principle be written in the same frame. In our case, however, we are comparing the energy of a non moving vortex configuration to that of a moving polar phase. We should therefore translate all velocities by . This change of reference frame creates, however a new dependent term in the Lagrangian , which contains an effective field on spins of magnitude . Finally, replacing the field in the expression of the infield excess energy of vortices and equating to yields the following equation with a dependence:
(27) 
The resulting line is shown, when it is higher than the zerofield one, in Fig. 1. It fits numerical observations well, even though the crossover between solitons and vortices is rather broad as we are working with finite and moderately small system sizes. Note that Eqs. (23) and (27) predict that the soliton states disappear at a finite value of : this is expected as solitons carry a nonzero topological charge (the vector field winds once around one direction of the Euclidean torus), and should not survive in the thermodynamic limit, unlike pairs of vortices and antivortices (the total winding number of which is ).
Iv Thermal Fluctuations
iv.1 Magnetization as a function of Temperature
Recalling that we perform constant energy simulations, the temperature should in principle be measured using velocity fluctuations. In the present case the coupling between the velocities and the spins modifies, however, the usual equipartition relations Bore et al. (2016):
(28)  
(29) 
The latter expression shows that the fluctuations of momenta and spins are coupled as long as , in such a manner that the temperature is not directly proportional to the translational kinetic energy.
The former, on the other hand, shows that the temperature is proportional to the rotational kinetic energy, which is thus the most natural thermometer.
We start by studying the effects of the finite temperature on the magnetic properties of the homogeneous phases of the system using the (intensive) magnetization modulus , which is the order parameter of the polar phase. Fig. 3, displays the equilibrium value of against for increasing values of at fixed (main panel) or increasing values of at fixed (inset). When varying , we observe three types of behaviours.

For small , when the ground state is the moving polar state, we observe a finitesize crossover between a (moving) ferromagnetic state at low temperatures and a non moving paramagnetic state at high temperatures (redtoorange curves of Fig. 3). This is the same crossover as the one reported in the case between the magnetized ground state and the paramagnetic high temperature phase: the polar state does not survive in the thermodynamic limit in due to low energy spin waves excitations. Casiulis et al. (2019) The effect of is to decrease the value of the magnetization modulus at finite temperatures with respect to the case. This can be understood by analyzing the spin waves excitations at low temperature within the framework of the continuumlimit theory described in Sec. III.2. Indeed, close to these spinwaves are described by a free Gaussian theory, that comes from the Taylor expansion of the spinspin interaction terms in the Hamiltonian at second order. Kosterlitz (1974) Spin waves suppress global magnetization as the the temperature increases. Furthermore, the effective spinspin coupling decreases with temperature when taking into account higherorder terms in the Taylor expansion. In the case, an effective field is generated on top of the spinspin interaction term. The direction of this effective field is opposite to the one of global magnetization, and its amplitude grows with . As a result, the effective spinspin coupling is further decreased by the addition of , leading to stronger spinwave suppression of the magnetization as increases. In short favors the bending of the magnetization field.

For large , when the ground state has zero magnetization modulus because of its topological structure, the magnetization modulus crosses over from , to a hightemperature paramagnetic regime where .

When takes moderately large values, although the ground state remains non polar (), the magnetization grows to a maximum at a finite temperature, before it crosses over to the paramagnetic regime at high temperatures. We shall focus on this intriguing phenomena in the next section IV.2.
When increasing at a fixed , as shown in the inset of Fig. 3, the same three scenarii are observed.
This is another sign of the qualitative similarity between the roles played by and , that was already pointed out in Sec. III.
It is also customary in the study of XY spins to introduce a “reduced” susceptibility associated to , Tobochnik and Chester (1979); Archambault et al. (1997); Casiulis et al. (2019)
(30) 
While it is not the usual response function of a magnetic system to an external field, this quantity captures the broadening of the distribution of that occurs at the finitesize paramagneticferromagnetic crossover with a peak that grows and shrinks with the system size. Archambault et al. (1997) In Fig. 3, we plot the “reduced” susceptibility per particle, against the temperature defined in Eq. (28), corresponding to the magnetization curves shown in Fig. 3. Upon increasing at a fixed (main panel), the susceptibility peak observed in the case of a polar ground state grows with and is shifted towards lower temperatures. The same behaviour is observed upon increasing at Casiulis et al. (2019). A peak is also observed in cases in which the magnetization features a local maximum [green curve in Fig. 3], and it is even higher and shifted to lower temperatures compared to the peak observed for polar states. Finally, the peak is suppressed at large enough, when the magnetization simply crosses over from the one of a topological ground state to that of a paramagnet. Likewise, when varying at a fixed value of (inset), the same scenario is observed: the peak of the extensive susceptibility grows, is shifted to lower temperatures, and is then suppressed. This behaviour is in sharp contrast with the smooth, algebraic evolution of the height of the peak reported in the case . Casiulis et al. (2019)
iv.2 Unusual OrderbyDisorder Scenario
As described in Sec. IV.1, a surprising feature of the model is that, for some finite value of or , the magnetization exhibits a local maximum at a finite temperature. In other words, the order parameter grows when thermal fluctuations are switched on. This behaviour is akin to OrderbyDisorder transitions Villain et al. (1980); Shender and Holdsworth (1996); Guruciaga et al. (2016) observed in frustrated magnets, whereby a system with a nontrivially degenerate ground states develops longrange order by the effect of classical or quantum fluctuations. In such geometrically frustrated spin models, the ObD phenomenon is due to a huge disproportion in the density of lowenergy excitations associated with a specific ground state. Here the mechanism must be of a different nature.
In fact, in the present case the growth of a spontaneous magnetization upon increasing the temperature is essentially due to the fact that the kinetic energy cost associated to collective motion in polar states decreases with . The reason is that, as already discussed in Sec. III, the magnetization modulus decreases as the temperature is increased due to spin waves fluctuations, and so does the kinetic energy cost induced by momentum conservation. Hence, at moderately large magnetized states become less and less costly compared to the soliton or vortex ones upon increasing the temperature, until a point where the free energy of the (moving) polar states crosses the one of the (non moving) soliton or vortex states and become the preferred minimum, in a way similar to a firstorder phase transition. When the temperature is further increased, the topological states become unstable, and the polar minima continuously merge into a single paramagnetic minimum, as in the case. Casiulis et al. (2019) This scenario is pictorially sketched in Fig. 4 in terms of a Landaulike free energy representation. The fact that the system develops a spontaneous magnetization and exhibits collective motion with a finite velocity of the center of mass, upon increasing the temperature starting from a topological ground state with zero magnetization, is thus due to the reentrant shape of the surface separating the polar states from the solitonic ones in the extended phase diagram of Fig. 1, when a vertical axis corresponding to temperature is added.
Note that this interpretation relies on the existence of metastable states, corresponding to local stable minima of the free energy at finite temperature. In the thermodynamic limit, metastablity is forbidden in finite dimensional systems due to the convexity of the free energy. The ObDlike transition described here is in fact only a smooth crossover, and is continuously suppressed as the system size is increased. Yet, at finite one can explicitly check in numerical simulations that both the topological and polar states are indeed locally stable at low enough temperature, as explained below. First, we finetune and to get a system very close to the polarsoliton groundstate crossover, and we choose to be rather small (here ). We then let the dynamics run at a finite but low temperature, . As shown in Fig. 4, switches between low and highmagnetization states can be observed, proving that both states are indeed local minima of the free energy. However, even for such small systems, one must wait very long times for switches to happen: this is an indication that the energy barrier between the two states is rather high compared to thermal fluctuations. One indeed expects that the energy barrier between a solitionic pattern and a polar state should scale linearly with the system size, since a collective reorganization of the spin degrees of freedom is necessary to switch between them.
Another indication of the existence of topological nonmagnetized states beyond the range of parameters for which they are the equilibrium configurations is provided by the following observation. Consider a system with and such that at zero temperature the ground state is the polar one, and quench it from a hightemperature to a very low temperature. The result, shown in Fig. 4, shows that because nonequilibrium vortices appear during such quenches (as already observed for Casiulis et al. (2019)), the system often ends up in a 4vortex state (whose centers being located at the equidistant points where red, yellow, green and blue regions meet).
iv.3 Qualitative Free Energy Expression
To go beyond the sketchy description of the free energy used in Fig. 4, we introduce simple expressions for the free energy densities of both families of states at low temperatures, in the same spirit as the discussion of the energy of the different ground states (Sec. III). First, the energy per particle of a soliton will be considered to be unaltered by temperature, and to remain of the same form as in Eq. (22). Then, we approximate the associated entropy per particle, , by considering that it is only associated to the number of ways in which one can draw a soliton pattern of size with particles, separated from their nearest neighbours by a typical distance ,
(31) 
The energy of a polar moving state at finite temperature should depend on the total magnetization of the system to reflect the momentum conservation constraint. Therefore, we modify the expression of introduced in Sec. III to include this effect,
(32) 
where is the intensive magnetization modulus. The exact form of at all temperatures and for all values of below is not known exactly. We thus choose to drop the dependence and to approximate by the exact expression obtained within a spinwave approximation Tobochnik and Chester (1979); Archambault et al. (1997),
(33) 
where is a constant that was evaluated from finitesize scaling Casiulis et al. (2019) in the case , yielding the value . Finally, assuming that the system is well described by the spinwave approximation in the considered domain of temperatures, we approximate the entropy per particle of magnetized states by counting how many ways there are to fit waves with a typical wavelength in an box,
(34) 
where we took the typical wavelength to be the spinwave correlation length , that itself is well approximated, at low temperatures, by Tobochnik and Chester (1979)
(35) 
We then have simple approximations for the free energies per particle and associated to polar and soliton states at low temperatures:
(36)  
(37) 
We can also approximate the free energy per particle of paramagnetic states by assuming that their magnetic energy averages to , and that their entropy is simply given by the choice of random angles for all particles, so that
(38) 
The expressions (36), (37) and (38) for the different branches of interest of the free energy density can be used to estimate which state is favoured at a given temperature for a given value of and . In particular, we can determine parameters such that for . Indeed, at the temperature such that these two branches cross, the system goes from a lowmagnetization soliton to a magnetized state. The crossing between and therefore corresponds to the temperature at which a maximum of the magnetization is reached, with height . The results obtained when varying for , as in Fig. 3, are shown in Fig. 5. As discussed in Sec. IV, the main effect of is to bring magnetic states higher up in energy at very low temperatures, so that they cross the solitonic branch at a finite temperature below the one at which it crosses the paramagnetic branch. From these simple arguments, one also predicts that the soliton branch directly crosses the paramagnetic branch for all values of larger than the value for which , and meet at a single point (here, ). These simple arguments thus allow us to rationalise the finite temperature results presented above.
V Ordering of Velocities
We are yet to describe how velocities order in polar moving states. In order to do so, is not the right object to consider, since it is forced to follow through the conservation of . It is however interesting to define an equivalent of the magnetization that is often called the polarity in the field of active matter Cavagna et al. (2010),
(39) 
which measures the alignment of unit velocities regardless of the velocity modulus. We furthermore define an observable that measures the alignment of unit velocities onto the spins,
(40) 
The variations of and with temperature, obtained for increasing at fixed or increasing at fixed are shown in Fig. 6. While velocities do magnetize at very low temperatures for low values of and , it is interesting to note that they do not behave like the magnetization: and are very small even for well magnetized states, and become significantly larger than zero only at very low temperature (note the logarithmic scale of the and axis).
One also sees that in the high or high regime, the velocities do not align on the spins in solitonic and vortex configurations. This is even more striking when comparing the magnetization and velocity fields at low temperatures, as illustrated by the snapshots in Fig. 6. The reason is that there is no true source of velocity alignment. At first sight, one would have expected the introduction of the selfalignment of the velocity on the spin in the Lagrangian to be responsible for a transfer of the ferromagnetic alignment of the spins to the velocities. This is actually not the case because the alignment of the velocities is actually entirely driven by the constraint imposed by the momentum conservation : . This constraint is only global and promotes local alignment between the velocities in a very indirect way. Also because needs to be small, in order to stabilize the polar state, the velocity of the center of mass is never that large and thus enforces only a weak constraint on the velocity alignment.
Vi Discussion and Conclusion
We studied a Hamiltonian model of interacting particles that carry continuous spins. The latter interact ferromagnetically with their neighbours, and are locally coupled to their velocities. The homogeneous phases exhibit a rich behaviour, including the existence of frustrationinduced topological defect configurations as ground states, ObDlike transition, and collective motion. Topological defects appear as a consequence of a classical spinorbit coupling, which is reminiscent of the socalled topological phase transitions observed in quantum models when a spinorbit coupling is added Hasan and Kane (2010), and could have deeper links with active matter, as parallels between quantum theories with a spinorbit coupling and flocking have recently been proposed Loewe et al. (2018). In Appendix C, we show that our results are unchanged when moving from to , or when changing the geometry of the boundary conditions. Our results can therefore be seen as robust.
From the point of view of magnetism, despite the fact that the model is at first sight a bit peculiar, one can think of it as an efficient dynamics to prepare topological defect configurations, for instance in dimension , in which case observing and describing topological excitations is notoriously challenging. Ackerman and Smalyukh (2017) From the point of view of collective motion, the model seems limited in terms of the amplitude of the velocities it gives access to, especially for large systems. One could in principle think of other coupling scheme between the spins and velocities, and hope for stronger effects. However, we show in Appendix B that the proposed spinvelocity coupling is, in fact, one of the only two ones that are compatible with welldefined Hamiltonian dynamics: adding coupling terms beyond the quadratic order in velocity, leads to a situation in which the Hamiltonian dynamics can not be obtained from the Legendre transform of the Lagrangian. Chi and He (2014) We have also checked that adding a quadratic coupling to the model with a linear coupling leads to the same phenomenology. The present model is therefore the only reasonable conservative model in which a spinvelocity coupling leads to collective motion. Yet, the case of a model with only a quadratic coupling but no linear coupling, , can be of interest to describe a situation in which a nematic ordering of the velocities onto the direction of the spin should be favoured, as in a recent model of Vicseklike particles with velocity reversals Mahault et al. (2018).
Acknowledgements.
We would like to warmly thank Michael Schindler for his great help on the numerical aspects of this paper. We would also like to thank Éric Bertin, Jorge Kurchan, JeanNoël Fuchs, and Pascal Viot for insightful discussions. Leticia F. Cugliandolo and Marco Tarzia are members of the Institut Universitaire de France.References
 Ramaswamy (2010) S. Ramaswamy, Annu. Rev. Condens. Matter Phys. 1, 323 (2010).
 Marchetti et al. (2013) M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Rev. Mod. Phys. 85, 1143(47) (2013), arXiv:1207.2929 .
 Cavagna et al. (2010) A. Cavagna, A. Cimarelli, I. Giardina, G. Parisi, R. Santagati, F. Stefanini, and M. Viale, Proc. Natl. Acad. Sci. 107, 11865 (2010).
 Zhang et al. (2010) H. P. Zhang, A. Be’er, E.L. Florin, and H. L. Swinney, Proc. Natl. Acad. Sci. 107, 13626 (2010).
 Bricard et al. (2013) A. Bricard, J.B. Caussin, N. Desreumaux, O. Dauchot, and D. Bartolo, Nature 503, 95 (2013).
 Deseigne et al. (2010) J. Deseigne, O. Dauchot, and H. Chaté, Phys. Rev. Lett. 105, 1 (2010).
 Csahok and Vicsek (1995) Z. Csahok and T. Vicsek, Phys. Rev. E 52 (1995).
 Bechinger et al. (2016) C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Rev. Mod. Phys. 88, 1 (2016).
 Weber et al. (2013) C. A. Weber, T. Hanke, J. Deseigne, S. Léonard, O. Dauchot, E. Frey, and H. Chaté, Phys. Rev. Lett. 110, 208001 (2013).
 Lam et al. (2015) K. D. N. T. Lam, M. Schindler, and O. Dauchot, New J. Phys. 17, 113056 (2015).
 Shankar et al. (2018) S. Shankar, S. Ramaswamy, M. C. Marchetti, and M. J. Bowick, Phys. Rev. Lett. 121, 108002 (2018).
 Dauchot and Démery (2019) O. Dauchot and V. Démery, Phys. Rev. Lett. 122, 068002 (2019).
 Bore et al. (2016) S. L. Bore, M. Schindler, K.D. N. T. Lam, E. Bertin, and O. Dauchot, J. Stat. Mech. Theory Exp. 2016, 033305 (2016).
 Villain et al. (1980) J. Villain, R. Bidaux, J.P. Carton, and R. Conte, J. Phys. 41, 1263 (1980).
 Shender and Holdsworth (1996) E. F. Shender and P. C. W. Holdsworth, in Fluctuations Order New Synth., edited by M. Millonas (Springer US, New York, NY, 1996) pp. 259–279.
 Casiulis et al. (2019) M. Casiulis, M. Tarzia, L. F. Cugliandolo, and O. Dauchot, J. Chem. Phys. 150, 154501 (2019).
 Kosterlitz and Thouless (1973) J. M. Kosterlitz and D. J. Thouless, J. Phys. C Solid State Phys. 6, 1181 (1973).
 Tobochnik and Chester (1979) J. Tobochnik and G. V. Chester, Phys. Rev. B 20, 3761 (1979).
 Dean (1996) D. S. Dean, J. Phys. A. Math. Gen. 29, L613 (1996).
 Barone et al. (1971) A. Barone, F. Esposito, C. J. Magee, and A. C. Scott, La Riv. del Nuovo Cim. 1, 227 (1971).
 Gouvea et al. (1990) M. E. Gouvea, F. G. Mertens, A. R. Bishop, and G. M. Wysin, J. Phys. Condens. Matter 2, 1853 (1990).
 Malomed (1991) B. A. Malomed, Phys. D Nonlinear Phenom. 52, 157 (1991).
 Shagalov (1992) A. G. Shagalov, Phys. Lett. A 165, 412 (1992).
 Sinitsyn et al. (2004) V. E. Sinitsyn, I. G. Bostrem, and A. S. Ovchinnikov, J. Phys. Condens. Matter 16, 3445 (2004).
 Kosterlitz (1974) J. M. Kosterlitz, J. Phys. C Solid State Phys 7, 1046 (1974).
 Archambault et al. (1997) P. Archambault, S. T. Bramwell, and P. C. W. Holdsworth, J. Phys. A. Math. Gen. 30, 8363 (1997).
 Guruciaga et al. (2016) P. C. Guruciaga, M. Tarzia, M. V. Ferreyra, L. F. Cugliandolo, S. A. Grigera, and R. A. Borzi, Phys. Rev. Lett. 117, 167203 (2016).
 Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
 Loewe et al. (2018) B. Loewe, A. Souslov, and P. M. Goldbart, New J. Phys. 20, 013020 (2018).
 Ackerman and Smalyukh (2017) P. J. Ackerman and I. I. Smalyukh, Phys. Rev. X 7, 011006 (2017).
 Chi and He (2014) H. H. Chi and H. J. He, Nucl. Phys. B 885, 448 (2014).
 Mahault et al. (2018) B. Mahault, X.c. Jiang, E. Bertin, Y.q. Ma, A. Patelli, X.q. Shi, and H. Chaté, Phys. Rev. Lett. 120, 258002 (2018).
 Diu et al. (1996) B. Diu, C. Guthmann, D. Lederer, and B. Roulet, Physique Statistique (Hermann, 1996).
 Verlet (1967) L. Verlet, Phys. Rev. 159, 183 (1967).
 Vest et al. (2014) J.P. Vest, G. Tarjus, and P. Viot, Mol. Phys. 112, 9 (2014).
 Nobre and Tsallis (2003) F. D. Nobre and C. Tsallis, Phys. Rev. E 68 (2003).
 Abramowitz and Stegun (1972) M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions (National Bureau of Standards, Applied Mathematics Series, 1972).
 Sausset and Tarjus (2007) F. Sausset and G. Tarjus, J. Phys. A Math. Theor. 40, 12873 (2007).
Appendix A Numerical methods
In this Appendix we motivate the choice of Molecular Dynamics as an efficient simulation method to study the model on which we focus on. In particular, we give details on the precise strategy leading to the figures shown in the main text. Moreover, we describe the additional difficulties encountered when simulating the same system in but with continuous (Heisenberg) spins, as described in App. C.
a.1 Choice of Molecular Dynamics
Simulating a system with a spinvelocity coupling presents several unusual aspects, both fundamental and practical. On the fundamental side, the fact that the system is not Galilean invariant implies that the Gibbs measure should take its full, framedependent form Diu et al. (1996); Bore et al. (2016)
(41) 
and are, respectively, the inverse temperature and the velocity of the center of mass, and are both fixed through an external bath in the canonical ensemble. and are the Hamiltonian and the total momentum of the configuration , respectively. Therefore, because of the nonGalilean character of the model, in the canonical ensemble, the variables imposed by the bath are the temperature and the velocity. As a consequence, usual canonical ensemble sampling via a Monte Carlolike scheme is not expected to be the most interesting setting since the aim of the model is to demonstrate the emergence of a spontaneous as a response to a lowering of the temperature or energy. It seems more natural, instead, to impose a total momentum . A possibility would then be to impose a temperature and a total momentum, that is to say, to work in an intermediary ensemble that lies between the microcanonical (conserved energy and momentum) and the canonical (conserved temperature and momentum). However, little is known on the peculiarities of nonGalilean statistical mechanics, so that we choose the simpler setting of the microcanonical ensemble, that also enables us to see the actual dynamics of the system, for instance, the propagation of solitons.
In practice, we therefore integrate the Hamiltonian equations of motion,
(42)  
(43)  
(44)  
(45) 
where is the unit vector obtained by rotating the spin of an anticlockwise angle, and we choose a total momentum . A technical issue, however, still remains. In Molecular Dynamics, the equations of motion are usually integrated using numerical tricks like the Verlet algorithm Verlet (1967), that ensure “symplecticity”, meaning that the momentum and the energy are rigorously conserved. These tricks, however, usually rely on updating the two sets of conjugate Hamiltonian variables at different times. Here, these methods are inapplicable due to the coupling between the velocity and the spin. Instead, we simply reproduced the simulation strategy used in Ref.[ Bore et al., 2016], namely, a fourthorder RungeKutta scheme that conserves the total momentum only approximately. However, the latter always remains small and fluctuates with a typical magnitude of in simulations, and is therefore much smaller than any other observable we considered.
a.2 Simulation strategy
In the main text, we simulate the dynamics starting from random states with uniformly distributed and drawn from centered, reduced Gaussian distributions. Such initial states were placed into a square box with periodic boundary conditions and, after giving some time for the dynamics to settle in, are subjected either to a numerical annealing or to a highrate quench. These procedures are implemented as follows.
Numerical annealing is performed by multiplying all rotational velocities by every 100 time units in our adimensional variable, with an integration time step equal to in the same units. This method enables us to reach lowenergy states which, if the cooling is slow enough, should be equilibrium states.
Quenches are carried out by multiplying all rotational velocities and momentum components by once, at some initial time. This method violently takes the system away from equilibrium, thereby enabling us to study the subsequent equilibration dynamics.
Whenever curves represent observable mean values, we mean that we averaged the results over independent configurations obtained by sampling initial conditions and recording states at times that are far enough from each other that the corresponding configurations can be considered to be independent.
a.3 Variant: case
In order to simulate the system in 3 dimensions, we also used the RungeKutta Molecular Dynamics method. In this higher dimensional case the proper Hamiltonian dynamics are given by
(46)  
(47)  
(48)  
(49)  
(50)  
(51) 
where we called the interaction part of the Hamiltonian. These equations pose a major problem to their numerical integration, as divergences at the poles appear in the evolution equations for and , at each time step. This problem, which is usual for simulations of rotors on spheres, is here rather tricky to solve Vest et al. (2014), and a variety of strategies can be tested. In our case, we chose to write modified evolution equations in which we consider spins to be arbitrary Cartesian vectors, but that should (within numerical accuracy) conserve the spin normalisation. The idea, adapted from previous works on Heisenberg spins Nobre and Tsallis (2003) and dynamics in a spherical geometry Vest et al. (2014) consists in writing
where is the usual angular momentum of the spins. For unit spins, , as their time derivative should be perpendicular to them. We then write the evolution equations:
(52)  
(53) 
These equations are natural consequences of the definition of with respect to , the actual canonical momentum associated to the spins. The advantage of this way of writing the evolution is that, even though it is not strictly symplectic in 3 dimensions, it does conserve the modulus of the spins up to numerical errors.
Appendix B The spinvelocity coupling
b.1 Breakdown of the LagrangiantoHamiltonian transformation
As mentioned in the main text, it is tempting to explore various forms of the spinvelocity coupling, in the hope that they may influence the phenomenology of the ground states in different ways.
A first remark is that due to the rotational invariance of the spins and velocities, any linear spinvelocity coupling with a shape
(54) 
where is a shorthand notation for the velocity of the particle , and is a rotation by an angle , would be equivalent to the definition we used. The momentum conservation would now impose a fixed angle determined by between the velocities and spins at low energies.
Another possibility is to add an additional nonlinear spinvelocity term to the Lagrangian. By doing so, one could hope to find an additional velocity dependence in the Hamiltonian, and possibly a reduction of the cost of kinetic energy. Let us first discuss the term that is perhaps the most natural to add, a quadratic coupling. The spinvelocity interaction terms in the Lagrangian would then be
(55) 
where we called the coupling constant of the linear term and the one in front of the quadratic coupling for clarity. This coupling looks interesting, as for positive values of it decreases the kinetic cost of velocities aligned with the spins. With this coupling, the canonical linear momentum reads
(56) 
An issue with this expression is that, as it is not linear in , it is not easy to invert it and write explicitly the velocity in terms of and only. Moreover, this fact poses a problem in the Hamiltonian formalism,
(57) 
where is the interaction part of the Lagrangian. The Hamiltonian above can be easily rewritten in terms of the ’s and ’s only, and reads
(58) 
As hoped when defining the quadratic term in the Lagrangian, this term decreases the cost of the kinetic energy, even in the Hamiltonian. However, one can clearly see that problems arise for , since in this range of values the energy can be made arbitrarily low for . In fact, in this regime, the dynamics themselves are illdefined, because the Lagrangian is no longer convex with respect to the velocities, and this leads to unphysical results like infinite accelerations or velocities even for a single particle system. While there are ways to exploit such theories using path integral formulations Chi and He (2014) as well as interpretations of their peculiarities in quantum models, we will here restrict ourselves to . Note that the convexity condition imposed by would lead to similar restrictions for other nonlinear couplings, with several forbidden intervals for their corresponding coupling constants.
Another issue is to rewrite the Hamiltonian in terms of its canonical variables. It is here more complicated than usual, as the only way to do so is to project the canonical momentum onto the vectors of interest, , , and . After some algebra we obtain
(59)  
(60)  
(61) 
where we defined the constants
(62) 
Using these projections, after some more algebra, we can finally rewrite the Hamiltonian as a function of the canonical variables only
(63) 
where the Hamiltonian coupling constants and are, interestingly, defined differently from the Lagrangian ones,
(64) 
and where we discarded a constant term defined as
(65) 
It is a priori much more difficult and maybe even impossible to perform this inversion for any other nonlinear couplings, as the vector projections used here no longer yield linear equations between the terms of interest. In particular, as tempting as a coupling of the form , where looks, it yields a system of equations that do not allow to write a Hamiltonian in terms of the corresponding canonical momentum. Also note that higher order polynomials in yield, instead of Eq. (59), polynomial equations in with degree higher or equal to two, so that several roots exist.
Therefore, the only two reasonable couplings one can imagine to get unequivocal Hamiltonian dynamics with a spinvelocity coupling are the linear and quadratic couplings. An important message is that velocity terms in a conservative model are very constrained, so that it is in general necessary to take the system out of equilibrium by breaking energy or momentum conservation if one wants to act arbitrarily on velocities. Any other nonlinear term will lead either to convexity or to inversion issues.
Regarding the quadratic coupling introduced in this part, it does not seem to make the system simpler nor to lead to new physics relevant to collective motion for . Indeed, if we reproduce the logic of the first study of the model with