Synchronization and spin-flop transitions for a mean-field XY model in random field

We characterize the phase space for the infinite volume limit of a ferromagnetic mean-field XY model in a random field pointing in one direction with two symmetric values. We determine the stationary solutions and detect possible phase transitions in the interaction strength for fixed random field intensity. We show that at low temperature magnetic ordering appears perpendicularly to the field. The latter situation corresponds to a spin-flop transition.


Introduction
Synchronization processes are ubiquitous in nature and, in particular, are rooted in human life from the metabolic processes in our cells to the highest cognitive tasks we perform as a group of individuals. Synchronously flashing fireflies, chirping crickets, violinists playing in unison, applauding audiences, firings of neuron assemblies, pacemaker cell beats, etc. are ensembles of units able to organize spontaneously allowing order to arise starting from disordered configurations [30,36]. Synchrony has attracted much interest in the last decades due to its relevance in many different contexts as biology [3,25,12], chemistry [14], ecology [39], climatology [37], sociology [26], physics and engineering [34,38,42]. The fundamental features of all the mentioned examples are the presence of many objects -each of which is oscillating in a proper sense -and the phenomenon of mutual influence between oscillations. The attempt of modeling these complex systems therefore led to consider large families of microscopic interacting oscillators. Typically such systems are far from being (easily) tractable. The milestone, pionereed by Winfree, was to consider biological oscillators as phase oscillators, neglecting the amplitude. He discovered that a population of non-identical oscillators, coupled by all-to-all interaction terms, can exhibit a temporal analog of a phase transition producing a remarkable cooperative phenomenon [43]. The model was subsequently refined by Kuramoto who proposed a soluble coupled phase oscillator model, which has a simple form and undergoes a transition from an incoherent to a coherent scenario [24]. The original model was deterministic and its dynamics described by a system of ordinary differential equations where oscillators lie on the complete graph and are coupled through the sine of their phase difference. Since then, several variants have been extensively investigated: more general interaction function and/or network [4,11,33,41]; with periodic forcing [32]; stochastic [7,18,31] and noisy in random environment [8,17,27]; active rotator models [19,21,35]; and possible combinations of these cases. It is impossible to properly account for the literature in this field, we refer to the review article [1] for more details and further references. From an application viewpoint, an interesting modification of the Kuramoto model is represented by the XY model in magnetic systems. The latter is particularly important on the two-dimensional regular lattice since it provides for a microscopic description of certain systems in solid state physics, in particular superfluid helium films [23]. Further applications include superconducting films [5], the Coulomb gas model [16], Josephson junction arrays [6] and nematic liquid crystals [29]. In this paper we consider a noisy and disordered mean-field version of the XY model. We deal with a population of N planar spins, represented by angular variables lying in [0,2π), with meanfield interaction. Each spin evolves stochastically, driven by Brownian motion, coupled with all the other particles and subject to an external random field (quenched disorder) that accounts for effects of the environment or anisotropies in the medium where particles live. The interaction term is of the same form as Kuramoto's and therefore is of "imitative" type in the sense that favor particle configurations where the spins are aligned. In the limit N → +∞ this model is accurately described by an ordinary differential equation (McKean-Vlasov equation). We study the long-time behavior of this equation and obtain the full phase diagram of the stationary solutions. In particular we show that at the macroscopic level the system may undergo a transition from an incoherent to a coherent state whenever the interaction is sufficiently strong compared with the field. Indeed the cooperative behavior of the system is a result of the competition between the coupling strength and the intensity of the external field. The coupling tries to keep aligned all the elements of the population, whereas the anisotropy breaks the symmetry by imposing at each spin a privileged orientation. Differently from [7], where the zero-field case is considered, when supercritical we do not have a continuum of coherent stationary solutions. Due to the presence of disorder, a finite number of them is selected. Namely the directions the system possibly orders are either parallel or perpendicular to the randomness and moreover only magnetic orderings orthogonal to the field are chosen by the free energy. This transition where spins align perpendicularly to the field is called spin-flop transition and is very special in particular on Z 2 , since Mermin-Wagner theorem forbids phase transition in d = 2 for models with continuous symmetry [28]. As far as the XY model on Z 2 is concerned, it is worth to mention that the occurrence of a spin-flop transition has been shown in presence of an alternating external field in [15] and in presence of a uniaxial random field in [9]. In comparison with the treatment in [2] we consider a dichotomic quenched disorder. The results in [2] are derived for an even distribution of the random field, but the techniques allow for an approximate equation of the critical curve valid only in the case of small variance field density. This requirement excludes the possibility of bi-modal distributions which indeed cover our choice. For the XY model in i.i.d. dichotomic quenched disorder we are able to provide a complete and rigorous phase portrait, finding another example of spin-flop transition appearing.
Our paper is organized as follows. In Section 2 we define the model and introduce the notation. In the following Section 3 we present our results. Section 4 is devoted to some conclusions and future perspectives. Finally in the Section 5 we present the proofs.

Description of the microscopic model
Let [0,2π) be the one dimensional torus. Let also η = (η j ) N j=1 ∈ {−1, +1} N be a sequence of independent, identically distributed random variables, defined on some probability space (Ω, F, P), and distributed according to µ. We assume µ = 1 2 (δ −1 + δ +1 ). Given a configuration x = (x j ) N j=1 ∈ [0,2π) N and a realization of the random environment η, we can define the Hamiltonian where x j is the spin at site j and hη j , with h > 0, is the local magnetic field associated with the same site. Let θ, positive parameter, be the coupling strength. For given η, the stochastic process x(t) = (x j (t)) N j=1 , with t 0, is a N-spin system evolving as a Markov diffusion process on [0,2π) N , with infinitesimal generator L N acting on C 2 functions f : [0,2π) N −→ R as follows: Consider the complex quantity where 0 r N 1 gives information about the degree of alignment (magnetization) of the spins and Ψ N measures the average value they are pointing. We can reformulate the expression of the infinitesimal generator (2) in terms of (3): The expressions (1) and (4) describe a system of mean field ferromagnetically coupled spins, each with its own random field and subject to diffusive dynamics. The two terms in the Hamiltonian have different effects: the first one tends to align the spins, while the second one tends to make each of them point in the direction prescribed by the magnetic field. Observe that when η j = +1 the spin is pushed toward 0; whereas, when η j = −1 toward π.
For simplicity, the initial condition x(0) is such that (x j (0), η j ) N j=1 are independent and identically distributed with law ν of the form with 2π 0 q 0 (x, η) dx = 1, µ-almost surely. The quantity x j (t) represents the time evolution on [0, T ] of j-th spin; it is the trajectory of the single j-th spin in time. The space of all these paths is C[0, T ], which is the space of the continuous function from [0, T ] to [0,2π), endowed with the uniform topology.

Results
We first describe the dynamics of the process (2), in the limit as N → +∞, in a fixed time interval [0, T ]. To this effect we shall derive a law of large numbers based on a large deviation principle on the path space. Later, the possible equilibria of the limiting dynamics will be studied.
The first result we state concerns the dynamics of the flow of empirical measures. It is a special case of what is shown in [2,10,20], so the proof is omitted. We need some more notation. For a given q : [0,2π) × {−1, +1} −→ R, we introduce the linear operator L q , acting on f : [0,2π) × {−1, +1} −→ R as follows: where r q e iΨq := Given η ∈ {−1, +1} N , we denote by P η N the distribution on (C[0, T ]) N of the Markov process with generator (2) and initial distribution λ. We also denote by the joint law of the process and the environment. A large deviation principle for P N allows to characterize the unique limit of the sequence {ρ N (·)} N 1 and, in particular, makes possible to provide a Fokker-Planck equation useful to describe the time evolution of the limiting probability measure.
admits a unique solution in C 1 [0, T ], L 1 (dx ⊗ µ) , and q t (·, η) is a probability density on [0,2π), for µ-almost every η and every t > 0. Moreover, for every ε > 0 there exists C(ε) > 0 such that for N sufficiently large, where, by abuse of notations, we identify q t with the probability q t (x, η)µ(dη)dx In other words Theorem 3.1 states that under P N , for every t ∈ [0, T ], the sequence of empirical measures {ρ N (t)} N 1 converges weakly, as N → +∞, to a measure ρ t admitting density

Phase diagram for the mean field limit
Equation (9) describes the infinite-volume dynamics of the system governed by generator (2). We are interested in the detection of t-stationary solutions and possible phase transitions. We recall that being t-stationary solution for (9) means to satisfy L q q ≡ 0. Next result gives a characterization of stationary solutions of (9).
is measurable and q(·, η) is a probability on [0,2π). Then q is a stationary solution of (9) if and only if it is of the form where Z(η) is a normalizing factor and (r, Ψ) satisfies the self-consistency relation There is a one-to-one correspondence between equilibrium distributions (10) and solutions of the self-consistency equation (11). Therefore, our analysis reduces to the study of the self-consistency relation that corresponds to a bi-dimensional fixed point problem of the form The quantities r and Ψ are the macroscopic counterpart of r N and Ψ N in (3) and still are indicators of global coherence. Solutions with r = 0 are called paramagnetic, those with r > 0 are called ferromagnetic.
Remark 3.1. If r = 0 the stationary distribution (10) is given by where Z(η) is a normalizing factor.
, is solution of (11) for all values of the parameters.
The next proposition shows that ferromagnetic solutions for (11) appear only with specific values of average position Ψ. In those cases r may be implicitly characterized in terms of first kind modified Bessel functions Moreover, r + has to satisfy where I v (·) denotes the first kind modified Bessel function of order v.
We will state now our main theorem. It is concerned with the investigation of under what conditions ferromagnetic solutions for (11) may occur. In particular, it is shown that the system undergoes several phase transitions depending on the parameters.
The values r + ,r + andr + depend on the parameters θ and h; therefore, they vary according to the phase we are considering.
The rich scenario depicted in Theorem 3.2 can be qualitatively summarized in the phase portrait presented in Fig. 3 ). Each region represents a phase with as many ferromagnetic solutions of (11) as indicated by the numerical label. The valueh is implicitly defined by (14). The lower blue curve is θ 1 , whereas the upper red one corresponds to θ 2 . The dashed green separation line is θ and is obtained numerically. Indeed, the latter is defined by a tangency relation. More hints about this curve will be given in the proof of Theorem 3.2 in Section 5.

Discussion and future perspectives
The paper is concerned with the study of the phase diagram for a mean-field planar XY model with external field. We first determined the N → +∞ limiting distribution on the path space 6 by deriving an appropriate law of large numbers based on a large deviation principle. Then we analyzed the long-time behavior of the system. Indeed we were able to characterize the equilibria as solutions of a fixed point equation and in turn to detect several phase transitions. Two basic mechanisms play a relevant role. On the one hand, a ferromagnetic-type interaction tends to align the spins by decreasing their phase difference. On the other, the presence of a magnetic field tends to separate the population in two different groups, pointing in opposite orientations. Clearly there is a competition between interaction and field intensity. In particular, above a critical threshold for the interaction strength a phase transition from a paramagnetic to a ferromagnetic state occurs. More precisely we have shown that, whenever the interaction is sufficiently strong, a net magnetization appears spontaneously with average spin displacement either perpendicular or parallel to the direction fixed by the field. We then obtained 2, 4 or 6 different possible ferromagnetic solutions. The richness of the phase diagram is due to the addition of a dichotomic random environment.

Simulations, free energy and stability
If we integrate the model numerically, how does (r N (t), Ψ N (t)) evolve? For concreteness, suppose we fix the field intensity h and vary the coupling θ. Simulations show that for all θ less than the threshold θ 1 , the spins act as if they were uncoupled: their values become distributed around 0 or π, as prescribed by the magnetic field, starting from any initial condition (see Fig. 4.2).  in blue (resp. red) spins subject to positive (resp. negative) field are displayed. We can see that in the long-run the spins tend to align with the site-dependent magnetic fields. As a result, approximately half of the particles are around angle 0 and the other half around π, giving a null total magnetization; more precisely, we get r N (t fin ) = 0.001. In the last snapshot the limiting distributions for the two families of spins are superposed; the solid blue (resp. red) line is q(x, +1) (resp. q(x, −1)) defined by (10).
But when θ exceeds θ 1 , this paramagnetic state seems to lose stability and r N (t) grows, reflecting the formation of a small cluster of spins that are aligned, thereby generating a collective phenomenon. Eventually r N (t) saturates at some level 0 < r N (∞) < 1 (see Fig. 4.3). With further increases in θ, more and more spins are recruited into the "aligned cluster" and r N (∞) grows as shown in Fig. 4.4.
Simulations indicate that the paramagnetic solution is globally stable for θ < θ 1 ; whereas, it becomes unstable for θ > θ 1 , when ferromagnetic equilibria arise. In the multiplicity phase the numerics further suggest that the dynamics approach either r + , π 2 or r + , 3π 2 . All the other ferromagnetic stationary points are unstable. In other words, it seems there are only two attracting states for each value of θ > θ 1 . Indeed dynamical simulations are confirmed by the analysis of the free energy. The model is reversible and therefore the evolution is driven by a free energy F θ,h that corresponds (up to an additive constant) to the large deviation functional of the invariant  in blue (resp. red) spins subject to positive (resp. negative) field are displayed. We can see that in the long-run the spins point angles close to π 2 . As a result, a spontaneous net magnetization appears; indeed, we have r N (t fin ) = 0.8 and the mean orientation of the particles is Ψ N (t fin ) = 1.64. In the last snapshot the limiting distributions for the two families of spins are superposed; the solid blue (resp. red) line is q(x, +1) (resp. q(x, −1)) defined by (10).

(b)
The time-evolutions of r N and Ψ N corresponding to the current simulation are displayed.  measure. Let q be a stationary solution of (9); following [13] we obtain To study the stability of the various equilibria we found, we checked their relative heights on the free energy surface: we solved numerically the self-consistency relation (13) and we plugged the obtain pair(s) (r, Ψ) in (15). Moreover, to visualize the energy landscape, we plotted directly the surface (15). See Fig. 4.5 for an example. We tested several choices of θ and h for each region 8 of the parameter space. In all multiplicity phases the free energy functional has minima at the spin-flop points r + , π 2 and r + , 3π 2 . Whereas, in phases 4 and 6, the ferromagnetic solutions having Ψ ∈ {0, π} are either maxima or saddles for F θ,h and hence always unstable. This study supports the idea that a spin-flop transition occurs when increasing the coupling strength. Finally we would remark that, if we fix h and vary θ, there is no ordering between r + (θ), r + (θ) andr + (θ). On the other hand, for fixed θ and variable h, we have a monotone relation between the different coherence indicators seen as functions of h. Namely, in phase 4 we always get r + (h) r + (h) and in regime 6 it holdsr + (h) r + (h) r + (h).

Critical fluctuations
An important and interesting further step would be to understand how macroscopic observables fluctuate around their mean values when the system is put at the critical point. In this regime to obtain a limit theorem describing the fluctuations of the empirical measure process as N → +∞ we construct a process of the form for suitable α > 0. There are two notable features of this rescaling. On the one hand, we have a non-Gaussian spatial scale (N 1 4 instead of N 1 2 ); this implies that critical fluctuations are spatially larger than non-critical ones. On the other, the process must be observed in fast time N α t because of the phenomenon of "critical slowing down". It means that the fluctuations persist over long time scale. We would like to determine the exponent α such that (16) admits a meaningful limit in the sense of weak convergence. It is not clear a priori what to expect as time scale. The addition of disorder may have a drastic impact on the fluctuation process and change the time scale at which it exists. In the paper [8] the authors analyze how disorder affects the dynamics of critical fluctuations for two different types of interacting particle system: the Curie-Weiss and Kuramoto models in random environment. The interesting point is that when disorder is added, spin and rotator systems belong to different universality classes, which is not the case for their non-disordered counterparts. Hence the disorder is responsible for destroying universality. Roughly speaking in the Curie-Weiss model the fluctuations produced by the disorder always prevail in the critical regime: these fluctuations evolve in a time scale which is much shorter (α = 1 4 ) than the corresponding one for homogeneous system (α = 1 2 ). For rotators, the disorder does not modify the "standard" slowing down (α = 1 2 ). The question is why does it happen? The random Kuramoto model in [8] is not reversible 1 and moreover presents a discrepancy between the symmetry type of the state and the disorder variables (rotational vs. up/down). We wonder if this difference is due to the reversibility/irreversibility or rather symmetry issues. The general idea is therefore to consider two modifications of the random Kuramoto model, aimed at getting a reversible system with disorder having either up/down or rotational symmetry, and then make a comparison between the time scale of critical fluctuations. The XY model can be read as the variation that accounts for the reversibility plus up/down symmetry case. As a first step it would be interesting to investigate its critical fluctuations and compare them to those of the Curie-Weiss model.
Both stability properties of the steady solutions and the behavior of fluctuations appear to be difficult to determine analytically and are left unsolved by the present paper. To prove them 1 Kuramoto model in random environment [8]. Given a configuration x ∈ [0,2π) N and a realization of the random environment η ∈ {−1, +1} N , we can define the Hamiltonian H N (x, η) where x j is the position of rotator at site j; the disorder term hη j , with h > 0, can be interpreted as its own frequency and θ > 0 is the coupling strength. It is important to notice that the system (17) is not reversible unless  rigorously one should have the complete control over the spectrum of the linearization of the operator (8) that is an open and difficult problem at the moment. We feel that this analysis may deserve much more room and defer some more detailed work to future research.

Proofs of Proposition 3.2 and Proposition 3.3
Every stationary solution (r, Ψ) has to satisfy the self-consistency relation (11), which is equivalent to conditions where q(x, η) = [Z(η)] −1 · exp{2θr cos(Ψ − x) + 2hη cos x}. By standard trigonometric formulas, equations (19) and (20) can be rewritten as All the integrals involved can be rephrased in terms of Bessel functions. We make the main steps explicit for 2π 0 cos x q(x, +1) dx, the remaining integrals can be dealt with similarly. We have, where Γ (·) is the Gamma function and I v (·) denotes the first kind modified Bessel function of order v. The derivation of (21) is postponed to Appendix A.1. In the same manner we calculate 2π 0 cos x q(x, −1)dx = −2π(h − θr cos Ψ) and the normalizing constants By plugging what we obtained into equations (19 ) and (20 ), we get Focus on equation (20 ). It is equivalent either to sin Ψ = 0 or the term into square brackets vanishes. By Lemma A.2 in Appendix A.2 we know that the function g(z) = I 1 (2 √ It is easy to see that under (c) equation (19 ) is always satisfied, giving the statement of Proposition 3.2. Moreover, under (a) or (b) equation (19 ) reduces to (13) and this concludes the proof of Proposition 3.3.

Proof of Theorem 3.2
From Proposition 3.3 we infer that we have only to consider the cases when Ψ ∈ 0, π 2 , π, 3π 2 . We divide the study in a few steps.
First observe that sin y q(y, +1) dy + We look for a positive solution of the fixed point equation r = F 1 (r). Observe that by (13) we can rewrite F 1 as and, since we are interested in solutions r = 0, our problem translates in finding a positive solution for the equationF 1 (r) = 1, I 0 (2h) andF 1 is continuous in [0,1] for all the values of the parameters.
Hence, ifF 1 (0) > 1 thenF 1 (r) intersects the horizontal line 1 exactly once; on the contrary, wheneverF 1 (0) 1 there are no crosses. To conclude it is sufficient to notice that θ 1 (h) as defined in the statement of Theorem 3.2 equals h I 0 (2h) Analysis for Ψ ∈ {0, π}. We stick on the case Ψ = 0, the other being similar. We want to show that for θ > θ 2 (h), there is exactly one ferromagnetic solution and moreover, that there exists a further critical value θ (h), with θ (h) < θ 2 (h), such that if θ θ (h) there are only paramagnetic solutions, while if θ (h) < θ < θ 2 (h) there are two ferromagnetic ones. If we set Ψ = 0, then q(x, η) = [Z(η)] −1 · exp {2(θr + hη) cos x} and the self-consistency relation (11) is equivalent to the conditions We must show that (24) has a positive solution and that (25)  We look for a solution of the fixed point equation r = F 2 (r). We have 1] for all values of the parameters.
• F 2 is strictly increasing; indeed, the first derivative of F 2 with respect to r is given by which is a strictly positive quantity, since expected value of a variance. Moreover, it is readily seen that 2θ Var q (0) (x,+1) (cos X) .
• The second derivative is equal to and changes sign depending on the parameters, so that is not possible to conclude by a standard concavity argument. Nevertheless from numerics we see that there is at most one sign change (we checked the number of zeros for (26)  on a grid of mesh-size 0.5). As a consequence, F 2 changes the curvature at most once 2 .
Therefore we can argue as follows. Since as r → +∞ the function F 2 (r) approaches 1 from below, it must be concave for large r. Then, if θ θ 2 (h) and F 2 (r) 0 in a right-neighborhood of r = 0, F 2 (r) is strictly concave on [0,1] for any values of the parameters and hence there is no intersection with the diagonal.
if θ θ 2 (h) and F 2 (r) > 0 in a right-neighborhood of r = 0, F 2 (r) changes curvature either below or above the diagonal, giving rise to none or precisely two positive fix points. The boundary between these two regions is represented by the dashed green line θ (h) in Fig. 3.1. It has been obtained numerically and corresponds to the choice of parameters where there exists r > 0 such that F 2 (r) = r and F 2 (r) = 1. Note that the curves θ 2 and θ coincide for h ∈ 0,h and then separate at h =h.
if θ > θ 2 (h), no matter if either F 2 (r) 0 or F 2 (r) > 0 in a right-neighborhood of r = 0, the curve F 2 (r) crosses the diagonal at precisely one positive r.

A Appendix
A.1 Derivation of formula (21) We devote this section to compute 2π 0 cos x q(x, +1) dx. To shorten our notation, let us introduce constants A := θr cos Ψ + h and B := θr sin Ψ .
Therefore, we have Z(+1) By writing the trigonometric functions in terms of the complex exponential we can expand the powers of binomials to get

Now observe that
• the only non-zero terms in the triple sum are those for which h + − k − 1 = 0; • j must be odd for the whole sum to be real; and thus, after the index change j → 2j + 1, To continue we need the following technical lemma.