Order by disorder: saving collective motion from topological defects in a conservative model

TL;DR

This paper investigates a 2D Hamiltonian model of interacting spins and velocities, revealing how topological defects influence collective motion and how thermal fluctuations can induce order in the system.

Contribution

It provides a detailed analysis of finite size ground states, continuum limit predictions, and the role of topological defects in a conservative spin-velocity model, including finite temperature effects.

Findings

Identification of ground states with and without topological defects.

Prediction of phase domains based on control parameters.

Observation of order-by-disorder transition at finite temperature.

Abstract

Using analytic and numerical methods, we study a $2d$ Hamiltonian model of interacting particles carrying ferro-magnetically coupled continuous spins which are also locally coupled to their own velocities. This model has been characterised at the mean field level in a parent paper. Here, we first obtain its finite size ground states, as a function of the spin-velocity 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 magnetisation presents a maximum at a finite temperature. This peculiar behaviour, akin to an order-by-disorder 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.

Figures

Equations

$$\mathcal{L}$$

$$\mathcal{H}$$

$$\mathbf{m}=\frac{1}{N}\sum _{i=1}^N{\mathbf{s}}_i,\text{mbox}and{\mathbf{v}}_G=\frac{1}{N}\sum _{i=1}^N{\dot{\mathbf{r}}}_i,$$

$$J(r)$$

$$U(r)$$

$${\mathcal{L}}_{kin}$$

$${\mathcal{L}}_{sv}$$

$${\mathcal{L}}_{int}$$

$${\mathcal{L}}_{int}\approx -\frac{1}{2}\sum _{k\ne i}U(r_{ik})+\frac{1}{2}\sum _{k\ne i}J(r_{ik})-\frac{1}{4}\sum _{k\ne i}J(r_{ik}){\theta }_{ik}^2,$$

$$\mathcal{V}\equiv \sum _{k\ne i}J(r_{ik}){\theta }_{ik}^2=\sum _{\mathbf{r}}\sum _{{\mathbf{r}}^{\prime }}\sum _{k\ne i}J(\mathbf{r},{\mathbf{r}}^{\prime })(\theta ({\mathbf{r}}^{\prime })-\theta (\mathbf{r}){)}^2\delta (\mathbf{r}-{\mathbf{r}}_{\mathbf{i}})\delta ({\mathbf{r}}^{\prime }-{\mathbf{r}}_{\mathbf{k}}).$$

$$\mathcal{V}\approx \sum _{\mathbf{r}}\sum _{\mathbf{r}+\mathbf{d}}\sum _{k\ne i}J(d)(\nabla \theta (\mathbf{r})\cdot \mathbf{d}{)}^2\delta (\mathbf{r}-{\mathbf{r}}_{\mathbf{i}})\delta (\mathbf{r}+\mathbf{d}-{\mathbf{r}}_{\mathbf{k}}).$$

$$\mathcal{V}\stackrel{a\to 0}{\to }\int d^2\mathbf{r}\rho (\mathbf{r})\int d^2\mathbf{x}\rho (\mathbf{r}+\mathbf{x})J(x)(\nabla \theta (\mathbf{r})\cdot \mathbf{x}{)}^2.$$

$$\mathcal{V}\approx \int d^2\mathbf{r}{\rho }_0^2\pi J_0\frac{{\lambda }^4}{4}(\nabla \theta (\mathbf{r}){)}^2.$$

$${\mathcal{L}}_c$$

$$\mathbf{v}$$

$$\left(\frac{{\partial }^2}{\partial t^2}-c^2△\right)\theta$$

$${\mathcal{L}}_c$$

$$\vartheta (\xi )$$

$${\vartheta }_{\Lambda }(x,y)$$

$${\vartheta }_{\pm }(x,y)=\pm \arctan \left(\frac{y-y_0}{x-x_0}\right)+\frac{R}{2}\log \left[(y-y_0{)}^2+(x-x_0{)}^2\right],$$

$$E_{pol}=\frac{1}{2}{K}^2-\frac{1}{2}z\bar{J}+E_U$$

$$E_{sol}\approx -\frac{z\bar{J}}{2}\left[\frac{1}{2}\cos \left(\frac{2\pi \bar{r}}{L}\right)+\frac{1}{2}\right]+E_U.$$

$$K_c^{pol-sol}=\sqrt{\frac{z\bar{J}}{2}\left(1-\cos \frac{2\pi \bar{r}}{L_c}\right)}=\sqrt{\frac{z\bar{J}}{2}\left(1-\cos \left[2\pi \bar{r}\sqrt{\frac{\rho }{N_c}}\right]\right)},$$

$$\Delta E_{4D-pol}^{2dXY}=E_{4D}^{2dXY}-E_{pol}^{2dXY}\approx 8\pi \bar{J}\ln \frac{d}{a},$$

$$N_c^{sol-vor}$$

$$\Delta E_{4D-pol}^h\approx E_{4D}^{2dXY}(h)-E_{pol}(h)=\Delta E_{4D-pol}^{2dXY}+8{\pi }^2\sqrt{h}\left(\sqrt{\frac{d}{a}}-1\right).$$

$$K_c^{4D-pol}=\frac{2\varphi z-\ln \sqrt{N_c}}{\pi \left(\sqrt{N_c}-1\right)}.$$

$$T$$

$$2T$$

$$X_{|m|}=N\left(⟨m^2⟩-⟨|m|{⟩}^2\right).$$

$$S_{sol}=\frac{1}{N}\ln \left(\frac{L}{\bar{r}}\right).$$

$$E_{pol}$$

$$m(N;T)$$

Full text

Collective motion, Topological defects, Frustrated magnetism, Order by Disorder

††preprint: APS/123-QED

I Introduction

Collective motion, the macroscopic ordering of velocities in a many-body system, is a key feature of assemblies of self-propelled 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 self-propelled 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 Vicsek-like models moving individuals are represented as polar self-propelled 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 Vicsek-type 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 be 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 self-propelled 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 motorised 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 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 two-dimensional 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 centre 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 take place in the mean-field (fully-connected) limit, in spite of momentum conservation Bore et al. (2016).

In the present paper, using molecular dynamics simulations and analytic calculations, we characterise the $2d$ behaviour of this model in states such that the system is homogeneous, focusing on the effects of the size, dimensionality, temperature, and spin-velocity coupling intensity. We report a finite-size magnetisation crossover that is accompanied by a collective ordering of the velocities at low enough temperatures and low enough values of the spin-velocity coupling constant, hereafter defined as $K$. This ordering breaks down at a finite value of $K$, that depends on the size of the system, as accurately predicted using simple energetic arguments. The resulting non-magnetised ground states feature topological defects that, by suppressing the magnetisation, enable the system to avoid a high kinetic energy cost. The crossover between magnetised and non-magnetised states resembles a first-order phase transition, as both families of states are metastable close to the critical value of the spin-velocity coupling. When heated up, systems prepared in a non-magnetised ground state display a maximum of the modulus of their magnetisation 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 magnetisation upon heating, and can be explained by the examination of the free energies of the various states of the system. Finally, the magnetised states are characterised 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 spin-velocity 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 low-temperature equilibrium states and predict their location in parameter space with simple energetic arguments, as well as their nature, resorting to a coarse-grained 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 first-order 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 spin-velocity 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 $N$ interacting particles that carry continuous planar spins of norm equal to $1$. The particles are confined to move in $2d$, in a periodic square box of linear size $L$. Their motion is described by the Lagrangian

[TABLE]

that depends on position and spin variables. The position of the $i$-th particle is denoted $\bm{r}_{i}$ and its velocity $\bm{v}_{i}=\dot{\bm{r}}_{i}$, while $r_{ik}=|{\bm{r}}_{i}-{\bm{r}}_{k}|$ is the distance between the centers of particles $i$ and $k$;$\theta_{i}$ is the angle that the spin forms with a reference axis and fully parametrises the continuous $2d$ spin $\bm{s}_{i}$ of unit modulus. The time derivative of the spin vector is indicated by $\dot{\bm{s}}_{i}$; $\theta_{ik}$ is the angle between the spin of particle $i$ and the one of particle $k$. Furthermore, $m$ is the mass of each particle and $I$ its moment of inertia. $U$ is a short-ranged, isotropic and purely repulsive two-body interaction potential and $J$ is a short-ranged and isotropic ferromagnetic coupling between the spins. We characterise these two potentials below. For future convenience we made explicit the typical amplitudes of the two-body potential, $U_{0}$, and the ferromagnetic coupling, $J_{0}$. $K$ is the parameter that controls the strength of the spin-velocity coupling, the term that yields the special properties to the model. Throughout this paper, we will only consider the case $K\geq 0$ as the case $K\leq 0$ can be recovered by simply flipping all the spins. In the following, we will introduce dimensionless variables according to the transformations $\bm{r}/\sqrt{I/m}\to\bm{r}$, $t/\sqrt{I/J_{0}}\to t$, $K/\sqrt{mJ_{0}}\to K$, $L/J_{0}\to L$, $U_{0}/J_{0}\to U_{0}$ and, in the notation hereafter, we absorb $U_{0}$ in the definition of the two-body potential $U$.

In the Hamiltonian formalism the Lagrangian (1) transforms into the Hamiltonian

[TABLE]

where we defined the canonical momenta $\omega_{i}=\dot{\theta_{i}}$ and $\bm{p}_{i}=\dot{\bm{r}}+K\bm{s}_{i}$. 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 centre of mass and the total magnetisation in the micro-canonical ensemble. In particular, setting the total momentum $\bm{P}=\sum_{i}\bm{p}_{i}$ to zero, the system spontaneously develops collective motion, whenever it magnetises. Its centre of mass velocity then satisfies $\bm{v}_{G}+K\bm{m}=0,$ where

[TABLE]

are the intensive magnetisation and the velocity of the centre 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 $\bm{P}=\bm{0}$ anti-align the particle velocities with their spins. Also because the velocity of a magnetised state grows linearly with $K$, low potential energy moving states become very expensive in terms of kinetic energy as $K$ increases. As a result, for high enough values of $K$ non-magnetised 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

[TABLE]

where $\Theta$ is a Heaviside step function, $\sigma$ is the range of the interactions, that we fixed to 1, and we took $U_{0}=4$ 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 $\bm{P}=\bm{0}$. We use random initial states with uniformly distributed $\left\{\bm{r}_{i},\theta_{i}\right\}_{i=1,\dots,N}$ and $\left\{\bm{p}_{i},\omega_{i}\right\}_{i=1,\dots,N}$ 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 high-rate quench. We assume that the outputs of these microcanonical simulations can be interpreted from the point of view of canonical-ensemble statistical mechanics when sampling over initial conditions. This equivalence supposes that the dynamics defined by the Hamiltonian (2) are ergodic in phase space. This ergodic hypothesis is often made in sufficiently interacting Hamiltonian systems, since the structure of their dynamics guarantees that Liouville’s theorem applies, and that a flat distribution of microstates in phase space is stable when the dynamics are run. Kardar (2007) While Liouville’s theorem is not sufficient to guarantee ergodicity, here the hypothesis is supported by the results of Ref. [Bore et al., 2016], where it is shown that mean-field canonical predictions for the magnetisation as a function of the energy and momentum qualitatively reproduce the numerical results obtained with MD simulations. Further details on the simulation and numerical methods are provided in App. A. The phase diagram of this model has been fully characterised in the case $K=0$ [Casiulis et al., 2019]. Because aligned spins also attract each other, the spin fluid undergoes a liquid-gas 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 $\rho=N/L^{2}\approx 2.81$ (or, equivalently, $\phi\approx 0.55$ in terms of the packing fraction defined as $\phi=\rho\pi\sigma^{2}/16$).

III Finite-Size 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 in the lowest attainable energy configuration. They correspond to minima of the potential energy, which we interpret as states at very low temperatures and therefore call “ground states”. As shown in Fig. 1, we find three kinds of ground states depending on the values of $K$ and $N$:

1. (i)

Polar states: These are fully magnetised states with $m=1$, with collective motion taking place in the direction opposite to the magnetisation as imposed by the momentum conservation, $\bm{v_{G}}=-K\bm{m}$. 2. (ii)

Solitonic states: These are non-moving states with $m\approx 0$, characterised by a continuous rotation of spins in one direction. 3. (iii)

Vortex states: These are non-moving states with $m\approx 0$, characterised by 4 topological point defects, 2 vortices and 2 anti-vortices.

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 magnetised states, which increases with $K$ and $N$ like $NK^{2}$. As these parameters are increased starting from a polar state, the system produces topological defects as a way of setting the global magnetisation to zero and therefore escape the collective motion imposed by momentum conservation. The strong frustration induced by the high kinetic energy cost associated to magnetised states with a finite velocity of the centre 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 $K=0$ 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 coarse-grained continuum theory for the magnetic properties of the model. We define $\mathcal{L}_{kin}$, $\mathcal{L}_{sv}$ and $\mathcal{L}_{int}$, the kinetic, spin-velocity and interaction contributions to the Lagrangian. They read

[TABLE]

respectively, with $\bm{v}_{i}=\dot{\bm{r}}_{i}$ and $\omega_{i}=\dot{\theta}_{i}$. The symbol $\sum\limits_{k\neq i}$ indicates a double sum over $i$ and $k$ from 1 to $N$ with no repeated indices included. Focusing first on the interaction term $\mathcal{L}_{int}$, we assume that the spins are locally well aligned, so that $\bm{s}_{i}\cdot\bm{s}_{k}\approx 1-\theta_{ik}^{2}/2$, where $\theta_{ik}=\theta_{i}-\theta_{k}$, in the neighbourhood defined by the range of $J(r)$. Therefore,

[TABLE]

The main effect of the first two terms is to ensure the uniformity of the density profile $\rho({\bf r})\approx\rho_{0}$. Accordingly, we will discard them in the following and concentrate on the last term, that we henceforth call $\mathcal{V}$, 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 coarse-grained fields on this mesh as:

[TABLE]

with $\delta$ being Kronecker symbols. First, rewriting the sum over $\bm{r^{\prime}}$ after the transformation $\bm{r^{\prime}}\rightarrow\bm{r}+\bm{d}$, and using the rotational and translation invariances of $J$; second expanding the difference of $\theta$ to the first order, we obtain

[TABLE]

We now take the continuum limit by taking the mesh size to zero and find that in the absence of singular values of $\bm{\nabla}\theta$:

[TABLE]

For $J(r)=J_{0}\Theta(\lambda-r)$, with $\Theta$ a step function, and recalling that $\rho({\bf r})\approx\rho_{0}$, we finally obtain, after integration of the second integral:

[TABLE]

The two other contributions to the Lagrangian, namely, $\mathcal{L}_{kin}$ and $\mathcal{L}_{sv}$ are local, so that the same procedure as above can be trivially applied. Finally, we can write a continuum Lagrangian theory of the form

[TABLE]

where $c^{2}\equiv J_{0}N_{\lambda}\lambda^{2}/16$ with $N_{\lambda}\equiv\rho_{0}\pi\lambda^{2}$ the number of particles in one interaction volume, and $\hat{\bm{e}}(\theta(\bm{r}))$ is a unit vector pointing in the direction given by $\theta$. This form of the Lagrangian clearly shows that $K\bm{v}$ plays the role of a field, and creates a non-perturbative term even in the low-temperature expansion which otherwise gives a free theory for the XY model at low temperatures. Kosterlitz and Thouless (1973)

We can now write the Euler-Lagrange equations for $\theta$ and $\bm{v}$ to describe the low-energy regime of this theory,

[TABLE]

where the additional $\pi/2$ comes from the derivative of $\hat{\bm{e}}(\theta(\bm{r}))$ with respect to $\theta$. Interestingly, at this level of description, Eq. (15) selects solutions such that the local momentum is zero. When injecting this zero-momentum condition into the equation on $\theta$, we simply recover a wave equation on $\theta$: this is the usual spin-wave regime Kosterlitz and Thouless (1973) at low temperatures, that describes magnetised low-temperature states in finite size. To describe states that are close in energy to a finite-size polar moving state, we set $\theta(\bm{r})=\theta_{0}+\delta\theta(\bm{r})$, and $\bm{v}=-K\hat{\bm{e}}(\theta_{0})+\bm{\delta v(\bm{r})}$ corresponding to a small perturbation around a polar state aligned along $\theta_{0}$, say $\theta_{0}=0$. The Lagrangian evaluated for such states reads,

[TABLE]

and contains a field term $K^{2}\cos\delta\theta$ 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 $\delta v\sim\delta\theta\ll 1$. Therefore, close to low-temperature polar ground states, the effective theory describing spin alignment is essentially a sine-Gordon theory with a field of amplitude $K^{2}$.

The sine-Gordon 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 sine-Gordon theory whose functional form can be shown to be Barone et al. (1971)

[TABLE]

where $\vartheta$ is the orientation of the local magnetisation field, $\xi=x-ct$ is the usual variable of propagative solutions, with a propagation velocity $c$, $\xi_{0}$ is an offset for this variable that sets the position of the centre of the soliton, $\chi=\pm 1$ is the so-called polarity of the soliton, $\gamma=1/\sqrt{1-c^{2}}$ is a Lorentz factor associated to the propagation velocity, and $h_{sG}$ is the in-plane magnetic field of the sine-Gordon theory. This function represents a localised rotation of $2\pi$ of the magnetisation field along the $x$ direction. Noticing that in simulations we only ever see one such rotation over the size of the whole box, we plot in Fig. 2$(a)$ the magnetic field with a constant amplitude and an orientation defined by the following function of position $\bm{r}=(x,y)$,

[TABLE]

where $2\Lambda$ is the range of the box in which we plot the field, that is centered on $(0,0)$. 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$(b)$, 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, as demonstrated by the video shown in the Supp. Mat., we find that solitonic patterns do move in simulations, and that their propagation velocity fluctuates.

The sine-Gordon 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 non-linear term in the sine-Gordon equation Shagalov (1992); Sinitsyn et al. (2004). One-defect solutions can then be written as

[TABLE]

where $(x_{0},y_{0})$ are the coordinates of the centre of the defect, $R$ is a parameter that tunes the twist of the spiral surrounding the defect, and the sign in front of the $\arctan$ determines whether the defect is a vortex ($+$) or an antivortex ($-$). An example of a 4-vortex configuration obtained with this functional form is shown in Fig. 2$(c)$, and compared to a 4-vortex numerical ground state shown in Fig. 2$(d)$. In this example, we chose the positions of the centers of the defects as well as the value of $R$ to best fit the numerical observations.

Altogether the low energy magnetic excitations of the continuous theory perfectly match the very low temperature magnetic states obtained numerically, confirming the intuition of the mechanism whereby the system escapes the kinetic energy cost of the moving polar state, when $K$ or $N$ are too large. Note, however, that this description fails to capture the velocity configuration. Indeed, in the continuum theory, we assume that the velocities (like the spins) are essentially all aligned. While this perturbative approach enables us to predict topological magnetic excitations, that are quite different from a perfect polar state, it still predicts that the velocities should be opposed to the magnetisation field everywhere, whereas in simulations soliton and vortex patterns do not display any velocity alignment, as we further discuss in Sec. V.

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

[TABLE]

where the velocity is everywhere assumed to be $-K\bm{m}$, $z$ is the number of neighbours of each particle, and $\overline{J}$ is the mean value of the ferromagnetic coupling for the inter-particle distances $\bar{r}$ in the considered configuration. The energy contribution due to the repulsive interaction $E_{U}$ 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 $2\pi\bar{r}/L$, while spins are perfectly aligned along the other direction. Then, the total energy of a solitonic state can be approximated by

[TABLE]

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 $\bar{J}=J(\bar{r})$. For a triangular lattice one has that $\bar{r}=\sigma/2\sqrt{\phi}$.

Equating now the two energies $E_{pol}$ and $E_{sol}$ yields an equation for the line separating polar and soliton ground states,

[TABLE]

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 $K^{pol-sol}_{c}L_{c}\approx\sqrt{\pi^{2}z\bar{J}\bar{r}^{2}}$, and stresses that $K$ and $L$ 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 vortex-antivortex pair compared to a fully magnetised configuration in an on-lattice XY model Kosterlitz and Thouless (1973). Taking into account only nearest-neighbour pairs of defects, this excess energy can be written as

[TABLE]

where $d$ is the distance between the two defects and $a$ is the size of the centre of the defect. Here $d=L/2$ and the vortex core diameter is of the order of the particle-particle spacing $a\simeq\bar{r}$. We then equate $\Delta E^{2dXY}_{4D-pol}$ to $E_{sol}-E_{pol}$ in its Taylor-expanded form and find a critical value for the number of particles above which four vortices states are more favorable than the solitonic ones:

[TABLE]

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 $K$ 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 $K$ on the vortex state energy is to reckon that $K$ plays the same role as an in-plane field in the magnetic model analog, and to use the expression for the energy of a vortex-antivortex pair in a field: Gouvea et al. (1990)

[TABLE]

In the magnetic analogue $h$ is the amplitude of an external in-plane field, while here it is an effective field due to the spin-velocity coupling $K$. According to the effective sine-Gordon theory (see Sec. III.2), $h=K^{2}$. An alternative and more intuitive manner to obtain the quadratic dependence on $K$ is to recall that Eq. (26) gives the excess energy of vortex-antivortex pairs with respect to a magnetised 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 $\bm{v_{G}}=-K\bm{m}$. This change of reference frame creates, however, a new $K$-dependent term in the Lagrangian $\mathcal{L}_{Gal}=-K^{2}\bm{m}\cdot\sum_{i=1}^{N}\bm{s}_{i}$, which contains an effective field on spins of magnitude $K^{2}$. Finally, replacing the field in the expression of the in-field excess energy of vortices and equating $E_{4D}^{2dXY}$ to $E_{sol}$ yields the following equation with a $K$ dependence:

[TABLE]

The resulting line is shown, when it is higher than the zero-field 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 $N$: this is expected as solitons carry a non-zero 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 [math]). Furthermore, the number of vortices or solitons, is a priori independent of the size of the system. Indeed, while defects usually appear in equilibrium due to local frustration, such structures here develop only because the total magnetisation should be zero for the kinetic energy to be low, resulting in a global frustration between the conservation of the momentum and the one of the energy. As a result, there is no typical length scale associated to frustration smaller than the size of the system itself, and the $m=0$ patterns are scale-free. One can also reason from the point of view of energetic costs: additional vortices or solitons would only create more local curvature, leading to an increased magnetic cost, but no kinetic energy gain. Finally, one could a priori circumvent defect formation by scaling $K$ with $N$, in the same way that one usually rescales interaction constants in mean-field couplings. Here, the right choice would be $K\sim N^{-1/2}$ so that the kinetic energy of a magnetised state remains unchanged when increasing $N$. However, we choose not to discuss such scalings at length in this paper, for two different reasons. First, since the collective speed in this system scales like $K$, such a coupling would still lead to vanishingly small velocities. Second, such interactions, while they are used in fully-connected descriptions, are rather artificial, and would take this model even further from realistic microscopic models.

IV Finite Temperature States

IV.1 magnetisation as a function of Temperature

Since we are performing microcanonical simulations, statistics can only be obtained at one value of the energy by sampling initial conditions. Rather than using the energy as a control parameter, we assume ergodicity and treat the system as one of statistical mechanics, as explained in Sec. II. In particular, we introduce the statistical temperature $T$, which is usually 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)

[TABLE]

The latter expression shows that the fluctuations of momenta and spins are coupled as long as $K\neq 0$, 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. Having checked that both definitions yielded the same mean value of the temperature, we heretofore only refer to the mean temperature over independent configurations at the same energy, $\langle T\rangle$, as $T$.

We start by studying the effects of the finite temperature on the magnetic properties of the homogeneous phases of the system using the (intensive) magnetisation modulus $m=\left\langle\left|{\bf m}\right|\right\rangle$, which is the order parameter of the polar phase. Fig. 3$(a)$, displays the equilibrium value of $m$ against $T$ for increasing values of $K$ at fixed $N$ (main panel) or increasing values of $N$ at fixed $K$ (inset). When varying $K$, we observe three types of behaviours.

• (i)

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

• (ii)

For large $K$, when the ground state has zero magnetisation modulus because of its topological structure, the magnetisation modulus crosses over from $m\simeq 0$, to a high-temperature paramagnetic regime where $m\sim 1/\sqrt{N}$.

• (iii)

When $K$ takes moderately large values, although the ground state remains non polar ($m(T=0)=0$), the magnetisation 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 $N$ at a fixed $K$, as shown in the inset of Fig. 3$(a)$, the same three scenarii are observed. This is another sign of the qualitative similarity between the roles played by $K$ and $N$, 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 $m$, Tobochnik and Chester (1979); Archambault et al. (1997); Casiulis et al. (2019)

[TABLE]

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 $|m|$ that occurs at the finite-size paramagnetic-ferromagnetic crossover with a peak that grows and narrows with the system size. Archambault et al. (1997) In Fig. 3$(b)$, we plot the “reduced” susceptibility per particle, $\chi_{|m|}=X_{|m|}/N$ against the temperature $T$ defined in Eq. (28), corresponding to the magnetisation curves shown in Fig. 3$(a)$. Upon increasing $K$ at a fixed $N$ (main panel), the susceptibility peak observed in the case of a polar ground state grows with $K$ and is shifted towards lower temperatures. The same behaviour is observed upon increasing $N$ at $K=0$ Casiulis et al. (2019). A peak is also observed in cases in which the magnetisation features a local maximum [green curve in Fig. 3$(b)$], and it is even higher and shifted to lower temperatures compared to the peak observed for polar states. Finally, the peak is suppressed at $K$ large enough, when the magnetisation simply crosses over from the one of a topological ground state to that of a paramagnet. Likewise, when varying $N$ at a fixed value of $K$ (inset), the same scenario is observed: the peak of the extensive susceptibility $X_{|m|}=N\chi_{|m|}$ 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 $K=0$. Casiulis et al. (2019)

IV.2 Unusual Order-by-Disorder Scenario

As described in Sec. IV.1, a surprising feature of the model is that, for some finite value of $K$ or $N$, the magnetisation 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 Order-by-Disorder transitions Villain et al. (1980); Shender and Holdsworth (1996); Guruciaga et al. (2016) observed in frustrated magnets, whereby a system with a non-trivially degenerate ground states develops long-range 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 low-energy 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 magnetisation upon increasing the temperature is essentially due to the fact that the kinetic energy cost associated to collective motion in polar states decreases with $T$. The reason is that, as already discussed in Sec. III, the magnetisation 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 $K$ magnetised 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 first-order 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 $K=0$ case. Casiulis et al. (2019) This scenario is pictorially sketched in Fig. 4$(a)-(b)$ in terms of a Landau-like free energy representation. The fact that the system develops a spontaneous magnetisation and exhibits collective motion with a finite velocity of the centre of mass, upon increasing the temperature starting from a topological ground state with zero magnetisation, is thus due to the re-entrant 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 ObD-like transition described here is in fact only a smooth crossover, and is continuously suppressed as the system size is increased. Yet, at finite $N$ 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 fine-tune $K$ and $N$ to get a system very close to the polar-soliton ground-state crossover, and we choose $N$ to be rather small (here $N=128,K=0.25$). We then let the dynamics run at a finite but low temperature, $T=0.03$. As shown in Fig. 4$(c)$, switches between low- and high-magnetisation 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 solitonic 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 non-magnetised states beyond the range of parameters for which they are the equilibrium configurations is provided by the following observation. Consider a system with $K$ and $N$ such that at zero temperature the ground state is the polar one, and quench it from a high-temperature to a very low temperature. The result, shown in Fig. 4$(d)-(g)$, shows that because non-equilibrium vortices appear during such quenches (as already observed for $K=0$ Casiulis et al. (2019)), the system often ends up in a 4-vortex 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 and in spatially homogeneous states, 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 $E_{sol}$ as in Eq. (22). Then, we approximate the associated entropy per particle, $S_{sol}$, by considering that it is only associated to the number of ways in which one can draw a soliton pattern of size $L$ with $N$ particles, separated from their nearest neighbours by a typical distance $\bar{r}$,

[TABLE]

The energy of a polar moving state at finite temperature should depend on the total magnetisation of the system to reflect the momentum conservation constraint. Therefore, we modify the expression of $E_{pol}$ introduced in Sec. III to include this effect,

[TABLE]

where $m(K,N;T)$ is the intensive magnetisation modulus. The exact form of $m$ at all temperatures and for all values of $K$ below $K_{c}$ is not known exactly. We thus choose to drop the $K$ dependence and to approximate $m$ by the exact expression obtained within a spin-wave approximation Tobochnik and Chester (1979); Archambault et al. (1997),

[TABLE]

where $b$ is a constant that was evaluated from finite-size scaling Casiulis et al. (2019) in the case $K=0$, yielding the value $b\approx 5.5$. Finally, assuming that the system is well described by the spin-wave approximation in the considered domain of temperatures, we approximate the entropy per particle of magnetised states by counting how many ways there are to fit waves with a typical wave-length in an $L\times L$ box,

[TABLE]

where we took the typical wave-length to be the spin-wave correlation length $\xi_{sw}$, that itself is well approximated, at low temperatures, by Tobochnik and Chester (1979)

[TABLE]

We then have simple approximations for the free energies per particle $F_{pol}$ and $F_{sol}$ associated to polar and soliton states at low temperatures:

[TABLE]

We can also approximate the free energy per particle of paramagnetic states by assuming that their magnetic energy averages to [math], and that their entropy is simply given by the choice of random angles for all particles, so that

[TABLE]

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 $K$ and $N$. In particular, we can determine parameters such that $F_{sol}=F_{pol}$ for $K\geq K_{c}$. Indeed, at the temperature such that these two branches cross, the system goes from a low-magnetisation soliton to a magnetised state. The crossing between $F_{sol}$ and $F_{pol}$ therefore corresponds to the temperature $T_{max}$ at which a maximum of the magnetisation is reached, with height $m_{max}$. The results obtained when varying $K$ for $N=128$, as in Fig. 3, are shown in Fig. 5. As discussed in Sec. IV, the main effect of $K$ 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 $K$ larger than the value $K_{max}$ for which $F_{sol}$, $F_{pol}$ and $F_{para}$ meet at a single point (here, $K_{max}\approx 0.36$). 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. Because of the conservation of the total momentum, any state with a macroscopic magnetisation also moves with a collective velocity $\bm{v}_{G}=-K\bm{m}$. Therefore, one expects velocity vectors to align together (and to select the same speed) as the temperature is lowered. It is thus interesting to define an order parameter that reflects the alignment properties of the velocities,

[TABLE]

which reflects the alignment of unit velocities regardless of the velocity moduli, and is often called the polarity in the context of active matter Cavagna et al. (2010). We furthermore define an observable that measures the alignment of unit velocities onto the spins,

[TABLE]

The variations of $\mu$ and $\nu$ with temperature, obtained for increasing $K$ at fixed $N$ or increasing $N$ at fixed $K$ are shown in Fig. 6$(a)-(d)$.

For low values of $K$ and $N$, such that the system has a polar ground state, both $\mu$ and $\nu$ grow to $1$ at low temperatures. However, these curves do not behave like the magnetisation: $\mu$ and $\nu$ become significantly larger than zero only at very low temperatures (note the logarithmic scale of the $\mu$ and $\nu$ axis). This is even more striking when comparing the magnetisation and velocity fields at low temperatures, as illustrated by the snapshots in Fig. 6$(f)-(i)$. In particular, in the snapshots 6$(f)-(g)$, we observe that even deep in the magnetised phase, where the system moves with a collective velocity, both the orientations (colours) and moduli (lengths and thicknesses) of the velocity arrows seem to vary a lot across the system. This comparatively weak ordering of the velocities compared to that of the spins can be explained by focusing on the distribution of velocity vectors. Rewriting the equipartition equation (29) associated to momenta in terms of the velocities, one finds that

[TABLE]

In particular, deep in the magnetised phase, one can assume that all spins are essentially aligned and equal to the (unit) magnetisation vector, $\bm{s}_{i}\approx\bm{m}$. In this simple case, the spin-velocity terms in Eq.( 41) cancel out, so that the usual equipartition relation, $2T=\left\langle v_{i}^{2}\right\rangle-\left\langle v_{i}\right\rangle^{2}$, is recovered. Consequently, at very low temperatures, the velocity vectors are essentially drawn from a Gaussian distribution centered at $\bm{v}_{G}=-K\bm{m}$, and with an isotropic variance $T$. This distribution is sketched as a density plot in Fig. 6$(e)$, with values growing from gray to yellow. From this sketch, it is immediately visible that, as long as the standard deviation of the velocity distribution is larger than $K$, meaning that $T\gtrsim K^{2}$, there is a significant probability for velocities to be pointing in any direction other than that of $\bm{v}_{G}$. Their moduli, however, are skewed towards larger values. Hence the aspect of the snapshot in Fig. 6$(g)$, in which velocities display a low polar order, but are biased towards larger, red arrows. In practice, for polar order to be visible, like in snapshot 6$(i)$, one needs to reach temperatures such that $T\ll K^{2}$, thus ensuring that the Gaussian distribution of the velocities does not take significant values close to $\bm{v}=\bm{0}$. This constraint on the value of the temperature to observe polar order of the velocities is usually absent in models that feature collective motion. The reason for this is two-fold. First, in abstract models of self-propelled particles like the Vicsek model Vicsek et al. (1995), the speed of particles is often fixed to a nominal value, so that the distribution of the velocities is naturally off-centered and very sharp. In that sense, a conservative setting makes it harder to observe a polar order of the velocities, since it imposes constraints on the width of the speed distribution. Second, in order to stabilise the polar state, $K$ needs to be small enough, so that the velocity of the centre of mass is never that large and $\mu$ can only grow at very low temperatures.

For very large values of $K$ or $N$, such that the ground state is not polar and no ObD is observed at finite temperatures, $\mu$ remains low at all temperatures. This is expected, since no collective velocity develops in this case. Furthermore, $\nu$ also remains small, thus indicating that the velocities do not locally align on the spins either. In other words, in vortex and soliton structures, the velocities do not draw solitons or vortices themselves, but simply remain disordered at low temperatures: velocity alignment is destroyed. This is proof that there is no true microscopic source of velocity alignment. At first sight, one would have expected the introduction of the self-alignment 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, even in vortex or soliton states. This is not the case because the alignment of the velocities is actually entirely driven by the constraint imposed by the momentum conservation: $\bm{v}_{G}=-K\bm{m}=0$. This constraint is only global and promotes local alignment between the velocities in a very indirect way.

Finally, for intermediate values of $K$ or $N$, however, a macroscopic magnetisation is observed at finite temperatures, following the ObD scenario described in Sec. IV.2. As a consequence, a collective velocity is also observed at finite temperatures, and it is accompanied by a local maximum of both $\mu$ and $\nu$. In other words, the ObD phenomenon saves collective motion from a total extinction caused by topological defects at finite temperatures in a finite range of the parameters. Furthermore, in the range of $K$ and $N$ such that we observe switches between polar and soliton states, such as those described in Fig. 4, the collective velocity itself also switches between low and high values. This spontaneous intermittent motion is highly unusual in a conservative setting, as such behaviours are usually expected in systems with friction, like granulars, Rao and Nott (2008) or activity, such as E. Coli bacteria that alternate between “runs”, i.e. ballistic motion along straight lines, and “tumbles”, i.e. pauses in the motion during which they re-orient themselves. Berg and Purcell (1977)

VI Discussion and Conclusion

We studied a $2d$ Hamiltonian model of interacting particles that carry continuous spins. The latter interact ferro-magnetically with their neighbours, and are locally coupled to their velocities. The homogeneous phases exhibit a rich behaviour, including the existence of frustration-induced topological defect configurations as ground states, ObD-like transition, and collective motion. Topological defects appear as a consequence of a classical spin-orbit coupling, which is reminiscent of the so-called topological phase transitions observed in quantum models when a spin-orbit coupling is added Hasan and Kane (2010), and could have deeper links with active matter, as parallels between quantum theories with a spin-orbit coupling and flocking have recently been proposed Loewe et al. (2018). In Appendix C, we show that our results are unchanged when moving from $d=2$ to $d=3$, 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 $d\geq 3$, 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. A solution would be to choose $K$ such that it scales with the size of the system, but this rather unphysical choice would also lead to a vanishingly small collective speed in the large $N$ limit. One could also 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 spin-velocity coupling is, in fact, one of the only two ones that are compatible with well-defined 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 $K\neq 0$ leads to the same phenomenology. The present model is therefore the only reasonable conservative model in which a spin-velocity coupling leads to collective motion. Yet, the case of a model with only a quadratic coupling but no linear coupling, $K=0$, 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 Vicsek-like particles with velocity reversals. Mahault et al. (2018)

Finally, in the context of assemblies of self-propelled particles, regular arrangements of vortices Sumino et al. (2012); Nagai et al. (2015) and soliton-like structures DeCamp et al. (2015); Putzig et al. (2016); Patelli et al. (2019) have been described in experiments, simulations, and microscopic models. This could be an indication that such structures are generally less costly than homogeneous motion in models of collective motion, regardless of the presence or absence of activity.

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, Jean-Noël Fuchs, and Pascal Viot for insightful discussions. Leticia F. Cugliandolo and Marco Tarzia are members of the Institut Universitaire de France.

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. 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 $3d$ but with continuous $3d$ (Heisenberg) spins, as described in App. C.

A.1 Choice of Molecular Dynamics

Simulating a system with a spin-velocity 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, frame-dependent form Diu et al. (1996); Attard (2002); Bore et al. (2016)

[TABLE]

$\beta$ and $\bm{v}_{G}$ are, respectively, the inverse temperature and the velocity of the centre of mass, and are both fixed through an external bath in the canonical ensemble. $\mathcal{H}$ and $\bm{P}$ are the Hamiltonian and the total momentum of the configuration $\mathcal{C}$, respectively. Therefore, because of the non-Galilean 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 Carlo-like scheme is not expected to be the most interesting setting since the aim of the model is to demonstrate the emergence of a spontaneous $\bm{v}_{G}$ as a response to a lowering of the temperature or energy. It seems more natural, instead, to impose a total momentum $\bm{P}$. 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 non-Galilean 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,

[TABLE]

where $\hat{\bm{s}}_{\perp}$ is the unit vector obtained by rotating the spin of an anticlockwise $\pi/2$ angle, and we choose a total momentum $\bm{P}=\bm{0}$. 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 fourth-order Runge-Kutta scheme that conserves the total momentum only approximately. However, the latter always remains small and fluctuates with a typical magnitude of $10^{-17}-10^{-14}$ 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 $\left\{\bm{r}_{i},\theta_{i}\right\}_{i=1..N}$ and $\left\{\bm{p}_{i},\omega_{i}\right\}_{i=1..N}$ 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 high-rate quench. These procedures are implemented as follows.

Numerical annealing is performed by multiplying all rotational velocities by $\lambda_{A}=0.9999$ every 100 time units in our adimensional variable, with an integration time step equal to $\delta t=10^{-3}$ in the same units. This method enables us to reach low-energy 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 $\lambda_{Q}=0.10$ 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 $10^{2}-10^{3}$ 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: $3d$ case

In order to simulate the system in 3 dimensions, we also used the Runge-Kutta Molecular Dynamics method. In this higher dimensional case the proper Hamiltonian dynamics are given by

[TABLE]

where we called $H_{int}$ 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 $\phi_{i}$ and $p_{\theta_{i}}$, 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

[TABLE]

where $\bm{L_{i}}=\bm{s_{i}}\wedge\bm{p{s_{i}}}$ is the usual angular momentum of the spins. For unit spins, $L_{i}^{2}=\bm{p_{s_{i}}}^{2}=\bm{\dot{s}_{i}}^{2}$, as their time derivative should be perpendicular to them. We then write the evolution equations:

[TABLE]

These equations are natural consequences of the definition of $\bm{L}_{i}$ with respect to $\bm{p_{s_{i}}}$, 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 spin-velocity coupling

B.1 Breakdown of the Lagrangian-to-Hamiltonian transformation

As mentioned in the main text, it is tempting to explore various forms of the spin-velocity 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 spin-velocity coupling with a shape

[TABLE]

where $\bm{v}_{i}=\dot{\bm{r}}_{i}$ is a short-hand notation for the velocity of the particle $i$, and $\mathcal{R}_{\varphi}$ is a rotation by an angle $\varphi$, would be equivalent to the definition we used. The momentum conservation would now impose a fixed angle determined by $\varphi$ between the velocities and spins at low energies.

Another possibility is to add an additional non-linear spin-velocity 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 spin-velocity interaction terms in the Lagrangian would then be

[TABLE]

where we called $K_{1}$ the coupling constant of the linear term and $K_{2}$ the one in front of the quadratic coupling for clarity. This coupling looks interesting, as for positive values of $K_{2}$ it decreases the kinetic cost of velocities aligned with the spins. With this coupling, the canonical linear momentum reads

[TABLE]

An issue with this expression is that, as it is not linear in $\dot{\bm{r}}_{i}$, it is not easy to invert it and write explicitly the velocity in terms of $\bm{p}_{i}$ and $\bm{s_{i}}$ only. Moreover, this fact poses a problem in the Hamiltonian formalism,

[TABLE]

where $\mathcal{L}_{int}$ is the interaction part of the Lagrangian. The Hamiltonian above can be easily rewritten in terms of the $\bm{v}_{i}$’s and $\bm{s}_{i}$’s only, and reads

[TABLE]

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 $K_{2}\geq 1/2$, since in this range of values the energy can be made arbitrarily low for $v_{i}\to\infty$. In fact, in this regime, the dynamics themselves are ill-defined, 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 $K_{2}<1/2$. Note that the convexity condition imposed by $K_{2}<1/2$ would lead to similar restrictions for other non-linear 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, $\bm{v}_{i}$, $\bm{s}_{i}$, and $\bm{p}_{i}$. After some algebra we obtain

[TABLE]

where we defined the constants

[TABLE]

Using these projections, after some more algebra, we can finally rewrite the Hamiltonian as a function of the canonical variables only

[TABLE]

where the Hamiltonian coupling constants $\tilde{K}_{1}$ and $\tilde{K}_{2}$ are, interestingly, defined differently from the Lagrangian ones,

[TABLE]

and where we discarded a constant term $E_{0}$ defined as

[TABLE]

It is a priori much more difficult and maybe even impossible to perform this inversion for any other non-linear 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 $\hat{\bm{v}}_{i}\cdot\bm{s}_{i}$, where $\hat{\bm{v}}_{i}=\bm{v}_{i}/v_{i}$ 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 $\bm{v}_{i}\cdot\bm{s}_{i}$ yield, instead of Eq. (60), polynomial equations in $\bm{v}_{i}\cdot\bm{s}_{i}$ 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 spin-velocity 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 non-linear term will lead either to convexity or to $\bm{p}\leftrightarrow\bm{v}$ 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 $K_{1}>0$. Indeed, if we reproduce the logic of the first study of the model with $K_{1}$ only, and consider only cases with $\bm{P}=\bm{0}$, we expect collective motion only for magnetised states. In these states, assuming that we only consider states with $\bm{s}_{i}=\bm{m},\forall i$ the momentum conservation constraint is now $\bm{v}_{G}=-\tilde{K}_{1}\bm{m}$, so that the velocity of the centre of mass will just be rescaled in the same way as the coupling constants in the Hamiltonian. Therefore, we roughly expect to observe the same states as in the case $K_{2}=0$. Further proof that we do not expect any huge qualitative change when adding $K_{2}$ can be obtained by computing the partition function and observables in a mean-field approximation. The full calculation is reported in Sec. B.2, and shows that the results are essentially the same as those obtained with $K_{1}$ only Bore et al. (2016). In particular, the magnetisation curves remain almost unaltered by the addition of $K_{2}$.

The case $K_{1}=0,K_{2}>0$, however, might be interesting to study on its own, as it favours nematic order between the velocity and the spins. This is reminiscent of a recent model of active matter in which velocities are allowed to switch between alignment and anti-alignment with an internal spin degree of freedom. Mahault et al. (2018) It however looks like a different problem than the one of collective motion, as we do not expect any spontaneous flow of particles to arise in the case $K_{1}=0$, whatever the value of $K_{2}$. In light of these results, we are confident that the model with a single $K_{1}$ should contain all the physics linked to the addition of a conservative, and polar spin-velocity coupling.

B.2 Mean-field calculation with $K_{2}$

We report in this Section the mean-field study of the Hamiltonian written in Eq. (64),

[TABLE]

where

[TABLE]

In order to get a mean-field description, we follow the same strategy as in Ref. Bore et al. (2016) where case $K_{2}=0$ was dealt with in detail: we set $U=0$ and we take the fully-connected limit for the ferromagnetic coupling by setting

[TABLE]

with $J$ a constant. Then, discarding a constant term, the Hamiltonian can be rewritten as

[TABLE]

and, recalling that the Gibbs measure is altered by the fact that this system is not Galilean invariant, the mean-field partition function reads

[TABLE]

where we used the short-hand notations $d^{2N}\bm{R}=\prod_{i=1}^{N}d^{2}\bm{r}_{i}$, $d^{2N}\bm{P}=\prod_{i=1}^{N}d^{2}\bm{p}_{i}$, $d^{N}\bm{\Theta}=\prod_{i=1}^{N}d\theta_{i}$, and $d^{N}\bm{\Omega}=\prod_{i=1}^{N}d\omega_{i}$. The integrals over the positions and angular velocities factorise and can be easily computed, yielding

[TABLE]

Calling $s_{i_{x,y}}$ the projections of spins onto the $x,y$ direction, the exponent in the last expression can be rewritten in a more compact form by using a vector notation:

[TABLE]

with

[TABLE]

and the integration over $\bm{p}_{i}$ readily calculated to obtain

[TABLE]

Moreover, using the fact that $\bm{s_{i}}$ is a unit vector, one finds

[TABLE]

Note that in the case $K_{2}>1/2$, that was already flagged as non-physical right from the definition of the quadratic term, one finds that det$A$ becomes negative; and in the limit $K_{2}\rightarrow 1/2$, this determinant goes through an infinite value. The last results can be used to write

[TABLE]

Then, we notice a number of identities:

[TABLE]

and we define $Z_{\theta}$ through

[TABLE]

so that we can rewrite

[TABLE]

with

[TABLE]

It is now useful to apply a Hubbard-Stratonovich transformation,

[TABLE]

and the definition $\bm{\gamma}(\bm{u})\equiv jN\bm{u}+\bm{h}+K_{1}\bm{v}_{G}$, to derive a more convenient expression for $\mathcal{I}$

[TABLE]

with

[TABLE]

Therefore, the full partition function reads

[TABLE]

In order to compute $J(\bm{u})$, we need to evaluate the integral

[TABLE]

With this aim, we use the identities Abramowitz and Stegun (1972):

[TABLE]

Then, integrating out odd terms that involve sines, and after some algebra, we find an exact expression for $f$:

[TABLE]

Taking advantage of the fact that $J(\bm{u})=\left[f(\beta\bm{\gamma},-\beta K_{2}v_{G}^{2})\right]^{N}$, we can now move on to the computation of the integral over $\bm{u}$:

[TABLE]

There is no easy way to compute this integral so that we have to resort to a saddle point approximation, by minimizing the expression in the exponential with respect to ${\bm{u}}$. After some replacements the expression to optimise reads

[TABLE]

Seeking zeros of the gradient of $\mathcal{F}$, first noticing that the $2d$ minimiser is attained for $\bm{u}$ aligned with $\bm{v_{G}}$, which fixes $\theta_{\gamma}=\mathrm{arccos}(\bm{\hat{\gamma}}\cdot\bm{\hat{v}}_{G})\equiv 0[\pi]$, we obtain an implicit equation for the optimum $\bm{u^{\star}}$:

[TABLE]

Under a critical temperature, this equation has two branches of minimizing solutions, corresponding to two minima of $\mathcal{F}$ that are not at the same height due to the presence of $K_{1}$. Bore et al. (2016) One can solve this equation numerically and compare the result to the case in which $K_{2}=0$. The results are essentially the same, except that $K_{2}$ makes the magnetisation more impervious to $\beta K_{1}v_{G}$, as shown in Fig. 7. This can be explained from the equation itself, as the sums over Bessel functions are naturally perturbative because of properties of Bessel functions of growing parameter Abramowitz and Stegun (1972), with larger corrections when $\beta K_{1}\alpha$ increases.

Appendix C Changing the geometry and dimensionality

We here discuss the generality of the results described in the main text when changing the geometry of the simulation box in $2d$, or increasing the dimension to $3d$.

C.1 Effect of the geometry

In the main text, every simulation was performed in a periodic square box. Topologically, this box is homotopic to a Euclidean (flat) torus Sausset and Tarjus (2007). This geometry imposes some constraints. For example, the total number of vortices has to be compensated by the total number of antivortices. Changing the topology of the simulation box could be interesting, as it would lead to different proportions of vortices and antivortices, and perhaps to whole new defects. The role of geometry in finite size is less clear: one could wonder whether the configurations of 2 vortices and 2 antivortices simply come from the square shape of the box, for instance. To address this question in a minimal way, we simulate our system in a periodic hexagon (which is also equivalent to a Euclidean torus Sausset and Tarjus (2007)) and in a periodic rectangle. Examples of a few ground states obtained under these conditions are shown in Fig. 8.

We observe the same phenomenology in hexagonal boxes as in square ones: at low $(K,N)$, the observed non-magnetic states are solitons, and when $K$ or $N$ increase, 2 vortices and 2 antivortices develop. In rectangular boxes, the more slender the box becomes, the more solitons are favoured. Indeed, using the argument given in the main text for the energy of a soliton, their cost decreases if the rotation is made smoother. Therefore, in rectangular boxes, the solitons will select the longest side and be favoured with respect to vortex configurations.

In summary, as long as the homotopy class of the simulation box does not change, and using the few examples shown in Fig. 8, we expect the same phenomenology as in square boxes. Possibly there will be a shift of the crossover lines in the $(K,N)$ plane that can be expected to depend on the precise geometry.

C.2 Effect of dimensionality

In this paper, we mostly worked with the model defined on a $2d$ space and with, consequently, $2d$ continuous spins. However, the crossover between magnetic and non-magnetic ground states, that is only due to an energetic trade-off, should survive in higher dimensionality of space and with spins of higher dimension. In order to check that dimensionality does not play an important role with respect to this qualitative feature, we now focus on the case $d=3$.

First, let us revisit the mean-field calculation in $3d$ and with Heisenberg spins. The main difference with respect to the $2d$ case with planar spins comes from the angular dependence in the momentum associated to the spins. More precisely, the Hamiltonian we want to study now is

[TABLE]

Let us focus on the angular part of the partition function: the $p_{\theta_{i}}$ contribution is simply a Gaussian integral that gives a constant. However, the integral over $p_{\phi_{i}}$ looks a bit more dangerous for future calculations because of the $\sin\theta_{i}$ factor. The first (Gaussian) integral reads:

[TABLE]

Then, using the fact that $0\leq\theta_{i}\leq\pi$, one can simply take the $\sin$ out of the square root for all values of $i$. Therefore, for each particle, we simply obtain a $\sin\theta_{i}$ factor that should be integrated over $\theta_{i}$. In fact, this is simply the factor we were missing in the partition function to recover $d\bm{s_{i}}=\sin\theta_{i}d\theta_{i}d\phi_{i}$! We thus find that the spins should be integrated over a sphere, and not simply over a cube in angular coordinates. In fact, the same trick can be reproduced for any higher dimension. In the end, we will simply recover the hyperspherical infinitesimal surface element. Therefore, in any $d\geq 3$, we get, after some usual algebra, the saddle-point solution

[TABLE]

where, as in the $2d$ case, we have to minimise a function of $\bm{u}$, an intermediate integration variable, defined through

[TABLE]

The value of the minimizing value, $\bm{u^{\star}}$, is given by an implicit equation, and $\bm{u^{\star}}$ is still the value of the mean magnetisation in the saddle-point approximation. The results in any dimensionality are essentially similar, except that in the absence of a rescaling of interactions with $d$, the ferromagnetic transition that we observe is sent to lower temperatures. We show the particular $d=3$ case in Fig. 9.

Having checked that, at the mean-field level, results for $d=3$ are qualitatively the same as in $d=2$, we now reproduce the same simulations as in the $2d$ case in $3d$, with classical Heisenberg spins instead of XY spins, and obtain ground states for a few values of $K$ and $N$ and a density such that no phase separation is observed. A few snapshots thus obtained in a periodic cubic box, and for a $3d$ packing fraction $\phi_{3d}\equiv\frac{N\pi\sigma^{3}}{48L^{3}}\approx 0.32$, are shown in Fig.10. and confirm that polar states and solitons are still observed in that case when varying $K$ and $N$. In Fig. 10$(d)$, which shows a slice of the soliton state shown in Fig. 10$(b)$, a $2d$-like vortex-antivortex structure is also seen. A priori, we could expect true $3d$ point defects to arise in $3d$ as well for larger systems and, possibly, larger values of $K$.