PHYSICAL REVIEW E 89, 052125 (2014)
Onset of synchronization in the disordered Hamiltonian mean-field model Juan G. Restrepo* and James D. Meiss† Department of Applied Mathematics, University of Colorado at Boulder, Boulder, Colorado 80309, USA (Received 13 January 2014; published 19 May 2014) We study the Hamiltonian mean field (HMF) model of coupled Hamiltonian rotors with a heterogeneous distribution of moments of inertia and coupling strengths. We show that when the parameters of the rotors are heterogeneous, finite-size fluctuations can greatly modify the coupling strength at which the incoherent state loses stability by inducing correlations between the momenta and parameters of the rotors. When the distribution of initial frequencies of the oscillators is sufficiently narrow, an analytical expression for the modification in critical coupling strength is obtained that confirms numerical simulations. We find that heterogeneity in the moments of inertia tends to stabilize the incoherent state, while heterogeneity in the coupling strengths tends to destabilize the incoherent state. Numerical simulations show that these effects disappear for a wide, bimodal frequency distribution. DOI: 10.1103/PhysRevE.89.052125
PACS number(s): 05.70.Ln, 05.45.Xt, 05.90.+m
I. INTRODUCTION
Hamiltonian systems with long-range interactions appear in many areas of physics, including systems with gravitational or Coulomb interactions [1,2], vortex dynamics in fluids [3–5], plasma physics [6–8], and free-electron lasers [9,10]. Many of the generic properties of emergent collective behavior in these systems can be studied with the Hamiltonian mean field (HMF) model [11–13]. In this model, N inertial rotors described by their conjugate phase, θn , and angular momentum, pn , variables interact under the Hamiltonian H =
N N pn2 1K − an am cos(θm − θn ), 2In 2 N n,m=1 n=1
(1)
where the nth rotor has moment of inertia In and coupling constant an . One can also view Eq. (1) as describing an interacting set of particles on a periodic domain (the circle) with masses In , and momenta pn , that are coupled through charges an . Despite its simplicity and analytical tractability, the HMF model exhibits some of the features present in many Hamiltonian systems with long-range interactions, including violent relaxation toward long lived quasistationary states, slow collisional relaxation, and phase transitions [1,5]. Because of these similarities, the HMF model has become an iconic testbed for the study of Hamiltonian systems with long-range interactions [14–18]. Previous studies (with the exception of Refs. [1,19,20]) have assumed that the rotors have identical moments of inertia and that they are equally coupled to all other rotors; i.e., In = an = 1. However, in many physical systems there is often some form of disorder. For example, stars in self-gravitating systems have a heterogeneous mass distribution [1], and vortices in 2D turbulence have a heterogeneous circulation distribution [21]. We study the effects of disorder in the HMF model to provide insights that may be valid more generally. The model Eq. (1) allows for the study of rotors with different masses or lengths (variable In ) and disorder in the interactions due to gravity (mass an ∝ In ) or
Coulomb forces (charge ∝ an ). We will assume that In and an are chosen from a joint distribution h(I,a), and that h vanishes when I or a is negative; the case that a > 0 can also be thought of as ferromagnetic. Since the coupling strength is represented by the parameter K, without loss of generality, we can assume that the mean coupling strength is 1: a = 1, where · denotes an average over the distribution h. We will show that even a small heterogeneity in the distribution of these parameters can drastically affect the onset of synchronization. This paper is structured as follows. In Sec. II, we present our theory and derive an expression for the critical coupling strength at the onset of synchronization as a function of the distribution of the parameters I and a as well as of the distribution of initial momenta. In Sec. III we illustrate our results with numerical experiments. We present our conclusions in Sec. IV. II. EFFECT OF DISORDER ON THE ONSET OF INSTABILITY
We are interested in the effect of disorder on the transition to synchronized behavior in the system described by Hamiltonian Eq. (1). The canonical equations corresponding to Eq. (1) are pn , − Kan R sin(θn − ψ) , (2) (θ˙n ,p˙ n ) = In where the complex order parameter is defined as Reiψ ≡
(3)
analogous to that used in the Kuramoto model [13,22]. In the continuum limit, N → ∞, this system can be formulated in terms of the density ρ(θ,p,t; I,a) of rotors with phase θ , momentum p, mass I , and charge a, at time t. In this limit, the order parameter becomes
*
[email protected] † [email protected] 1539-3755/2014/89(5)/052125(5)
N 1 an eiθn ; N n=1
Re
iψ
∞
= 0
052125-1
0
∞
∞ −∞
2π
aeiθ ρ(θ,p,t; I,a)dθ dpdI da. 0
©2014 American Physical Society
JUAN G. RESTREPO AND JAMES D. MEISS
PHYSICAL REVIEW E 89, 052125 (2014)
The evolution of the density is given by the continuity equation (in this context often referred to as a Vlasov equation [14,23]): p ∂ρ ∂ρ ∂ρ + − aKR sin(θ − ψ) = 0. ∂t I ∂θ ∂p
(4)
This equation admits incoherent equilibria corresponding to densities of the form ρ = G(p; I,a)/(2π ), for which R = 0. Since the masses and charges do not evolve in time, it is convenient isolate their distribution, h(I,a), and rewrite G(p; I,a) = gI,a (p)h(I,a); thus, gI,a (p) is the distribution of momenta conditioned on I and a. A linear stability analysis of the incoherent solution, analogous to that in Refs. [12,24,25], implies the existence of a critical coupling strength Kc and critical frequency ω for the onset of instability, determined by ∞ gI,a (p) Kc ∞ ∞ 2 1=− a I h(I,a) − dp dI da, 2 0 0 −∞ p − I ω (5) ∞
0= 0
∞
0
a 2 I h(I,a)gI,a (I ω)dI da,
= ∂gI,a (p)/∂p and - denotes the principal value integral. From this analysis, one would expect that if the initial density G is known and K is increased adiabatically, the rotors would remain incoherent for K Kc (i.e., R ≈ 0), and that the incoherent state would lose its stability for K > Kc (i.e., R would become nonzero). By contrast, numerical simulations with the simple incoherent initial distribution ρ(θ,p,0; I,a) = g0 (p)h(I,a)/(2π ), for which g0 is independent of the parameters, show that this state may remain stable for values of K much larger than the predicted Kc (e.g., as we will show later, for a narrow distribution of masses with a relative width ∼0.1, the system can remain stable up to K ∼ 10Kc ). Our goal is to understand this modification to the onset of instability. We will obtain analytical results under the assumption that the distribution of effective frequencies θ˙ = p/I [cf. (2)] has a “small enough” width around its mean. Simulations will show what happens when this assumption is violated. To begin, we note that the linear stability analysis that leads to the dispersion relation Eq. (5) also shows that perturbations to the density that are independent of θ are marginally stable. In the classical analysis of Strogatz and Mirollo for the Kuramoto model, these perturbations are forbidden in order to conserve the number of oscillators [24]. For Eq. (4), however, any incoherent perturbation ρ → G(p; I,a)/(2π ) + η(p; I,a) with zero average is allowed and undamped. Such fluctuations should be expected in the simulations for finitely many oscillators, and these may modify the predicted threshold from Eq. (5). To quantify this we assume that, instead of being zero, the order parameter Eq. (3) fluctuates around a small rotating ¯ i t + z, where R¯ 1, is a coherent term, i.e., Reiψ = Re frequency that will be determined self-consistently, and z represents time-dependent fluctuations with zero mean. For the Kuramoto model it has been shown that R¯ ∼ N −1/2 [26], and numerical simulations [18] have shown the same scaling holds for the HMF model. The existence of a dominant frequency should be reasonable if the initial distribution of frequencies is narrow enough. Defining θ˜n = θn − t and p˜ n = pn − In , where gI,a (p)
Eq. (2) becomes p˜ n , p˙˜ n = −Kan R¯ sin(θ˜n ) − Kan Im(ze−iθn ). θ˙˜n = In In the absence of the fluctuating term, the rotors are decoupled and each has constant energy E¯ n = p˜ n2 /(2In ) − Kan R¯ cos(θ˜n ). We treat the term Kan Im(ze−iθn ) as a stochastic perturbation to the Hamiltonian dynamics, with the important characteristic that this perturbation conserves the total energy Eq. (1). In the context of the HMF model, it has been shown that the stationary distribution of energies in the ensemble is numerically very close to a Boltzmann distribution [18]; i.e., the density of rotors ¯ 2 ), where σ 2 is the with energy E¯ is proportional to exp(−E/σ temperature. Thus, below the onset of synchronization, letting R¯ = 0, the stationary density of rotors with momentum p, given I and a, becomes 1 (p − I )2 , (6) gI,a (p) = √ exp − 2I σ 2 2π I σ 2 a Gaussian distribution with mean I and variance I σ 2 . This density will evolve from the initial density, g0 (p), over a time scale to be determined. To determine , we note that both the discrete Eq. (2) and continuum Eq. (4) systems preserve ∞ the average momentum p, which is initially P = −∞ pg0 (p)dp and becomes ∞∞∞ 0 0 −∞ pgI,a (p)h(I,a)dpdI da = I in the stationary state. Therefore, = P I −1 .
(7)
Similarly, conservation of energy determines the temperature σ 2 in terms of the variance σ02 of the original distribution of momenta g0 (p). Equating the initial energy and the energy in the stationary state with distribution Eq. (6) gives 1 . σ 2 = I −1 σ02 + P 2 1 − I I −1 We now proceed to calculate the critical coupling strength from Eq. (5) for the marginal distribution Eq. (6). Since gI,a (I ) = 0, then ω = . Inserting this in the first equation of Eq. (5) gives our main result: 1 I −1 2 2 . (8) Kc ≈ 2 2 σ0 + P 1 − a I I −1 To compare this result with the homogeneous HMF model, consider the case h(I,a) = δ(I − I0 )δ(a − a0 ). If the initial distribution g0 (p) were itself Gaussian, then Eq. (8) predicts Kc0 = 2
σ02 . a02 I0
First let us consider the effect of heterogeneity in the distribution of masses. Since the Cauchy-Schwarz inequality implies I I −1 1, we have Kc Kc0 . Thus, mass heterogeneity enhances the stability of the incoherent state. On the other hand, whenever there is heterogeneity in the charge distribution, the factor a 2 in Eq. (8) will be larger than a02 . Thus, charge heterogeneity reduces Kc , destabilizing the incoherent state. We also note that if the initial distribution g0 (p) is not Gaussian, Eq. (8) might predict a lower Kc than what would have been obtained using g0 (p) in Eq. (5).
052125-2
ONSET OF SYNCHRONIZATION IN THE DISORDERED . . .
PHYSICAL REVIEW E 89, 052125 (2014)
FIG. 1. (Color online) Estimated critical coupling as a function of width of the mass distribution, ε, for example (i) using average momenta as indicated (symbols), and corresponding Kc predicted from Eq. (8) (solid lines). Inset: estimated (Kc − 2σ02 )/P 2 compared to the prediction 2ε2 /3. The inset also includes simulations obtained using an uniform initial distribution of momenta (see text). III. NUMERICAL EXPERIMENTS
In this section, we will illustrate our results from Sec. II with three numerical experiments. In example (i) the initial distribution g0 (p) is a Gaussian centered at P with standard deviation σ0 = 0.35. There are N = 1 000 rotors with a uniform distribution of I on the interval [1 − ε,1 + ε], but all have a = 1. In this case, it can be shown that the coefficient of variation (standard deviation over mean) of the distribution of the effective frequencies is bounded by 2(σ0 /P )2 + ε2 for ε < 0.8, so we expect our theory to apply when ε 1 and P σ0 = 0.35. First, we test our prediction Eq. (8) for Kc . For a given value of P and ε, we simulate Eq. (2) by increasing K by 0.02 every 10 000 time units from the initial value K = 0, and use the last 5 000 time units for each K to estimate R by a time average. The instability threshold was estimated as the value of K at which the averaged R last exceeds 0.1. While this is a rough estimate intended to be used for relatively small N , it allows us to see how Kc varies as the parameters of the system are changed. In Fig. 1(a) we plot the estimated Kc as a function of the width ε for the seven values of average momentum P shown. The curves show the prediction ε 2 Kc (ε) = atanh(ε) σ02 + P 2 1 − ε atanh(ε) of Eq. (8). The inset shows that the computed value of (Kc − 2σ02 )/P 2 collapses onto the solid black curve, given from the theory (for small ε and σ0 ) by 2ε2 /3. Since our results depend only on the mean P and variance σ0 of g0 (p), we also consider a uniform distribution of initial momenta g0 (p) with the same standard deviation σ = 0.35 and various P . The rescaled values of Kc are also included in the inset of Fig. 1, and they collapse onto the same curve. To illustrate the validity of Eq. (6) in the incoherent state, we plot, in Fig. 2(a), the initial masses and momenta of N = 1 000 rotors using a uniform distribution of initial momenta g0 (p), and, in Fig. 2(b), their masses and momenta at T = 10 000 for K/Kc = 0.32. The solid line in Fig. 2(b) indicates the predicted mean I of the stationary distribution gI,a (p)
FIG. 2. (Color online) Evolution of the momentum-mass distribution for 1 000 rotors each with a = 1, given a uniform distribution of masses I with I = 1 and of momenta with p = P = 5. (a) Momentum-mass distribution at t = 0 with half-width ε = 0.15 in mass and standard deviation σ0 = 0.35 in momentum, and (b) for t = 10 000 when K = 0.2. The solid line is p = I with = P = 5; (c) CDFs of momenta p conditioned on I for I = 0.9, 1, and 1.1, for t = 0, and (d) for t = 10 000. Theoretical CDFs are also shown (thin red curves).
using Eq. (7). The empirical cumulative distribution function (CDF) of momenta p conditioned on I is shown in Fig. 2(c) at t = 0, for three values of I . These were calculated from the rotors with masses in three slices (I − 0.025, I + 0.025) of the distribution. The curves are simply the theoretical CDF for the uniform distribution. In Fig. 2(d) we show the CDFs at t = 10 000 for the same three values of I averaged over the time interval [2 500, 10 000]. The curves show that the corresponding theoretical Gaussian CDFs, calculated using Eq. (6), agree reasonably well with the numerical distributions. Our analysis assumes that the initial distribution g0 (p) relaxes to the stationary distribution gI,a (p); that is, the system must remain incoherent long enough for this relaxation to occur. To illustrate this, we estimate σ 2 from the simulations 2 as the mean of the instantaneous conditional variance, σI,a . According to Eq. (8), the effective critical coupling strength 2σ 2 I −1 /a 2 should approach Kc as the distribution relaxes to Boltzmann. This quantity is shown in Fig. 3 (squares), as a function of K when it is increased, again by 0.02 each t = 10 000 time units. Initially, the effective critical coupling strength is 2σ02 , but as K grows adiabatically, the distribution relaxes to Eq. (6), and the effective critical coupling strength approaches the value predicted by Eq. (8) (dashed line). When K becomes larger than this critical value, the incoherent state loses its stability, as can be seen in the evolution of the order parameter R (circles in the figure). If K is increased too rapidly or N is too large, the relaxation of the initial distribution to its steady state might not occur and the incoherent state could become unstable earlier. Indeed, in experiments with a uniform initial distribution (not shown in the figure) the relaxation is slower than when it is a Gaussian, and so K must be increased more slowly to allow the distribution to relax to its steady state. In agreement with previous observations [17], we observe that the relaxation time scales linearly with N . Thus, our results
052125-3
JUAN G. RESTREPO AND JAMES D. MEISS
PHYSICAL REVIEW E 89, 052125 (2014)
FIG. 3. (Color online) Order parameter R (circles) and estimated 2σ 2 I −1 /a 2 (squares) as K is increased for ε = 0.1 and P = 5 in example (i). The solid line indicates the stationary value of Kc predicted by Eq. (8).
apply in situations in which N is not too large or when the system is allowed to reach equilibrium (e.g., as when we slowly changed K). In addition, if Eq. (8) predicts a value smaller than Kˆ c predicted for the initial distribution, g0 (p), our computations show that when K ∈ [Kc ,Kˆ c ] the incoherent state is only metastable, since it will become unstable as g0 (p) relaxes to Eq. (6) [27]. In conclusion, example (i) has shown how heterogeneity in the masses combined with nonzero total momentum P results in the stabilization of the incoherent state: Kc is increased compared with the case of identical oscillators. Heterogeneity in the distribution of charges a has the opposite effect. For example (ii), we compare the case in which a = 1 for all oscillators to that in which a is uniformly distributed in [0,2]. The masses are again uniformly distributed in [1 − ε,1 + ε]. Numerical estimates of Kc as a function of ε are compared to the predictions (solid lines) in Fig. 4. Note that Kc is smaller when the charges are heterogeneous, as predicted by Eq. (8). For example (iii), we show a case in which the assumption of a narrow frequency distribution is violated and our analysis breaks down. Now the distribution of masses and charges is h(I,a) = 12 {δ[I − (1 + ε)] + δ[I − (1 − ε)]}δ(a − 1),
FIG. 5. Simulated (symbols) and theoretical (solid line) values of Kc as a function of ε for example (iii). The dashed line indicates the value 2σ02 for Kc predicted in the absence of disorder.
i.e., two peaks equidistant from I = 1, with the charges set to a = 1. The initial distribution g0 (p) is a Gaussian with mean P = 5 and σ0 = 0.35. Figure 5 shows the numerical estimates (symbols) versus the predicted Kc of Eq. (8) (solid line) as a function of ε. These agree when the separation between the peaks (2ε) is small; however, when ε becomes too large, the onset of instability suddenly drops near to the value 2σ02 (dashed line) predicted in the absence of disorder. In this second regime, the stationary distribution gI,a (p) for each value of I is observed to have mean P , instead of shifting to the predicted I of Eq. (6). A more detailed study of this transition is left for future research.
IV. CONCLUSION
In conclusion, we have investigated analytically and numerically the stability of the incoherent state in the disordered Hamiltonian mean field model. We found that finite size effects can induce correlations between the momenta, moments of inertia, and coupling constants of the rotors, modifying the onset of instability of the incoherent state. Indeed, heterogeneity in the moments of inertia tends to stabilize the incoherent state, while heterogeneity in the coupling strengths tends to destabilize the incoherent state. For sharply peaked parameter distributions, we developed an analytical formula for the modified critical coupling strength. Our analysis also qualitatively describes the behavior observed for broader distributions. Finally, we discovered a novel transition for a bimodal distribution of masses. Our results provide new insights into the factors affecting the phase transition in the HMF model, an iconic testbed for the study of long-range Hamiltonian systems, and provide a motivation to search for analogous results in more specific systems.
ACKNOWLEDGMENTS FIG. 4. (Color online) Estimated (symbols) and theoretical (solid line) values of Kc as a function of ε for example (ii). Inset: enlargement of the region 0 ε 0.3.
J.G.R. acknowledges a useful discussion with Ed Ott. J.D.M. acknowledges the support of NSF Grant No. DMS1211350.
052125-4
ONSET OF SYNCHRONIZATION IN THE DISORDERED . . .
PHYSICAL REVIEW E 89, 052125 (2014)
[1] P. Chavanis, J. Vatteville, and F. Bouchet, Euro. Phys. J. B 46, 61 (2005). [2] T. Tatekawa, F. Bouchet, T. Dauxois, and S. Ruffo, Phys. Rev. E 71, 056111 (2005). [3] H. Aref, Ann. Rev. Fluid Mech. 15, 345 (1983). [4] J. Weiss, A. Provenzale, and J. McWilliams, Phys. Fluids 10, 1929 (1998). [5] P.-H. Chavanis, in Dynamics and Thermodynamics of Systems with Long-Range Interactions (Springer, Berlin, 2002), pp. 208– 289. [6] J. Tennyson, J. Meiss, and P. Morrison, Physica D 71, 1 (1994). [7] D. del Castillo-Negrete and M. C. Firpo, Chaos 12, 496 (2002). [8] N. Balmforth, Comm. Nonlinear Sci. Numer. Simulat. 17, 1989 (2012). [9] A. Antoniazzi, Y. Elskens, D. Fanelli, and S. Ruffo, Euro. Phys. J. B 50, 603 (2006). [10] R. Bachelard, T. Manos, P. De Buyl, F. Staniscia, F. S. Cataliotti, G. De Ninno, D. Fanelli, and N. Piovella, J. Stat. Mech. (2010) P06009. [11] T. Konishi and K. Kaneko, J. Phys. A 25, 6283 (1992). [12] S. Inagaki and T. Konishi, Publ. Astron. Soc. Japan 45, 733 (1993). [13] M. Antoni and S. Ruffo, Phys. Rev. E 52, 2361 (1995). [14] J. Barr´e, F. Bouchet, T. Dauxois, S. Ruffo, and Y. Y. Yamaguchi, Physica A 365, 177 (2006).
[15] A. Campa, A. Giansanti, and G. Morelli, Phys. Rev. E 76, 041117 (2007). [16] T. Dauxois, V. Latora, A. Rapisarda, S. Ruffo, and A. Torcini, in Dynamics and Thermodynamics of Systems with Long-Range Interactions (Springer, Berlin, 2002) pp. 458–487. [17] V. Latora, A. Rapisarda, and S. Ruffo, Physica D 131, 38 (1999). [18] Y. Y. Yamaguchi, J. Barr´e, F. Bouchet, T. Dauxois, and S. Ruffo, Physica A 337, 36 (2004). [19] A. Ciani, D. Fanelli, and S. Ruffo, in Long-Range Interactions, Stochasticity, and Fractional Dynamics (Springer, Berlin, 2011), pp. 83–132. [20] S. De Nigris and X. Leoncini, Phys. Rev. E 88, 012131 (2013). [21] A. Bracco, J. McWilliams, G. Murante, A. Provenzale, and J. Weiss, Phys. Fluids 12, 2931 (2000). [22] J. A. Acebr´on, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler, Rev. Mod. Phys. 77, 137 (2005). [23] F. Bouchet and T. Dauxois, Phys. Rev. E 72, 045103 (2005). [24] S. H. Strogatz and R. E. Mirollo, J. Stat. Phys. 63, 613 (1991). [25] M. Y. Choi and J. Choi, Phys. Rev. Lett. 91, 124101 (2003). [26] E. J. Hildebrand, M. A. Buice, and C. C. Chow, Phys. Rev. Lett. 98, 054101 (2007). [27] A. Pluchino, V. Latora, and A. Rapisarda, Physica D 193, 315 (2004).
052125-5
Onset of synchronization in the disordered Hamiltonian mean-field model Juan G. Restrepo* and James D. Meiss† Department of Applied Mathematics, University of Colorado at Boulder, Boulder, Colorado 80309, USA (Received 13 January 2014; published 19 May 2014) We study the Hamiltonian mean field (HMF) model of coupled Hamiltonian rotors with a heterogeneous distribution of moments of inertia and coupling strengths. We show that when the parameters of the rotors are heterogeneous, finite-size fluctuations can greatly modify the coupling strength at which the incoherent state loses stability by inducing correlations between the momenta and parameters of the rotors. When the distribution of initial frequencies of the oscillators is sufficiently narrow, an analytical expression for the modification in critical coupling strength is obtained that confirms numerical simulations. We find that heterogeneity in the moments of inertia tends to stabilize the incoherent state, while heterogeneity in the coupling strengths tends to destabilize the incoherent state. Numerical simulations show that these effects disappear for a wide, bimodal frequency distribution. DOI: 10.1103/PhysRevE.89.052125
PACS number(s): 05.70.Ln, 05.45.Xt, 05.90.+m
I. INTRODUCTION
Hamiltonian systems with long-range interactions appear in many areas of physics, including systems with gravitational or Coulomb interactions [1,2], vortex dynamics in fluids [3–5], plasma physics [6–8], and free-electron lasers [9,10]. Many of the generic properties of emergent collective behavior in these systems can be studied with the Hamiltonian mean field (HMF) model [11–13]. In this model, N inertial rotors described by their conjugate phase, θn , and angular momentum, pn , variables interact under the Hamiltonian H =
N N pn2 1K − an am cos(θm − θn ), 2In 2 N n,m=1 n=1
(1)
where the nth rotor has moment of inertia In and coupling constant an . One can also view Eq. (1) as describing an interacting set of particles on a periodic domain (the circle) with masses In , and momenta pn , that are coupled through charges an . Despite its simplicity and analytical tractability, the HMF model exhibits some of the features present in many Hamiltonian systems with long-range interactions, including violent relaxation toward long lived quasistationary states, slow collisional relaxation, and phase transitions [1,5]. Because of these similarities, the HMF model has become an iconic testbed for the study of Hamiltonian systems with long-range interactions [14–18]. Previous studies (with the exception of Refs. [1,19,20]) have assumed that the rotors have identical moments of inertia and that they are equally coupled to all other rotors; i.e., In = an = 1. However, in many physical systems there is often some form of disorder. For example, stars in self-gravitating systems have a heterogeneous mass distribution [1], and vortices in 2D turbulence have a heterogeneous circulation distribution [21]. We study the effects of disorder in the HMF model to provide insights that may be valid more generally. The model Eq. (1) allows for the study of rotors with different masses or lengths (variable In ) and disorder in the interactions due to gravity (mass an ∝ In ) or
Coulomb forces (charge ∝ an ). We will assume that In and an are chosen from a joint distribution h(I,a), and that h vanishes when I or a is negative; the case that a > 0 can also be thought of as ferromagnetic. Since the coupling strength is represented by the parameter K, without loss of generality, we can assume that the mean coupling strength is 1: a = 1, where · denotes an average over the distribution h. We will show that even a small heterogeneity in the distribution of these parameters can drastically affect the onset of synchronization. This paper is structured as follows. In Sec. II, we present our theory and derive an expression for the critical coupling strength at the onset of synchronization as a function of the distribution of the parameters I and a as well as of the distribution of initial momenta. In Sec. III we illustrate our results with numerical experiments. We present our conclusions in Sec. IV. II. EFFECT OF DISORDER ON THE ONSET OF INSTABILITY
We are interested in the effect of disorder on the transition to synchronized behavior in the system described by Hamiltonian Eq. (1). The canonical equations corresponding to Eq. (1) are pn , − Kan R sin(θn − ψ) , (2) (θ˙n ,p˙ n ) = In where the complex order parameter is defined as Reiψ ≡
(3)
analogous to that used in the Kuramoto model [13,22]. In the continuum limit, N → ∞, this system can be formulated in terms of the density ρ(θ,p,t; I,a) of rotors with phase θ , momentum p, mass I , and charge a, at time t. In this limit, the order parameter becomes
*
[email protected] † [email protected] 1539-3755/2014/89(5)/052125(5)
N 1 an eiθn ; N n=1
Re
iψ
∞
= 0
052125-1
0
∞
∞ −∞
2π
aeiθ ρ(θ,p,t; I,a)dθ dpdI da. 0
©2014 American Physical Society
JUAN G. RESTREPO AND JAMES D. MEISS
PHYSICAL REVIEW E 89, 052125 (2014)
The evolution of the density is given by the continuity equation (in this context often referred to as a Vlasov equation [14,23]): p ∂ρ ∂ρ ∂ρ + − aKR sin(θ − ψ) = 0. ∂t I ∂θ ∂p
(4)
This equation admits incoherent equilibria corresponding to densities of the form ρ = G(p; I,a)/(2π ), for which R = 0. Since the masses and charges do not evolve in time, it is convenient isolate their distribution, h(I,a), and rewrite G(p; I,a) = gI,a (p)h(I,a); thus, gI,a (p) is the distribution of momenta conditioned on I and a. A linear stability analysis of the incoherent solution, analogous to that in Refs. [12,24,25], implies the existence of a critical coupling strength Kc and critical frequency ω for the onset of instability, determined by ∞ gI,a (p) Kc ∞ ∞ 2 1=− a I h(I,a) − dp dI da, 2 0 0 −∞ p − I ω (5) ∞
0= 0
∞
0
a 2 I h(I,a)gI,a (I ω)dI da,
= ∂gI,a (p)/∂p and - denotes the principal value integral. From this analysis, one would expect that if the initial density G is known and K is increased adiabatically, the rotors would remain incoherent for K Kc (i.e., R ≈ 0), and that the incoherent state would lose its stability for K > Kc (i.e., R would become nonzero). By contrast, numerical simulations with the simple incoherent initial distribution ρ(θ,p,0; I,a) = g0 (p)h(I,a)/(2π ), for which g0 is independent of the parameters, show that this state may remain stable for values of K much larger than the predicted Kc (e.g., as we will show later, for a narrow distribution of masses with a relative width ∼0.1, the system can remain stable up to K ∼ 10Kc ). Our goal is to understand this modification to the onset of instability. We will obtain analytical results under the assumption that the distribution of effective frequencies θ˙ = p/I [cf. (2)] has a “small enough” width around its mean. Simulations will show what happens when this assumption is violated. To begin, we note that the linear stability analysis that leads to the dispersion relation Eq. (5) also shows that perturbations to the density that are independent of θ are marginally stable. In the classical analysis of Strogatz and Mirollo for the Kuramoto model, these perturbations are forbidden in order to conserve the number of oscillators [24]. For Eq. (4), however, any incoherent perturbation ρ → G(p; I,a)/(2π ) + η(p; I,a) with zero average is allowed and undamped. Such fluctuations should be expected in the simulations for finitely many oscillators, and these may modify the predicted threshold from Eq. (5). To quantify this we assume that, instead of being zero, the order parameter Eq. (3) fluctuates around a small rotating ¯ i t + z, where R¯ 1, is a coherent term, i.e., Reiψ = Re frequency that will be determined self-consistently, and z represents time-dependent fluctuations with zero mean. For the Kuramoto model it has been shown that R¯ ∼ N −1/2 [26], and numerical simulations [18] have shown the same scaling holds for the HMF model. The existence of a dominant frequency should be reasonable if the initial distribution of frequencies is narrow enough. Defining θ˜n = θn − t and p˜ n = pn − In , where gI,a (p)
Eq. (2) becomes p˜ n , p˙˜ n = −Kan R¯ sin(θ˜n ) − Kan Im(ze−iθn ). θ˙˜n = In In the absence of the fluctuating term, the rotors are decoupled and each has constant energy E¯ n = p˜ n2 /(2In ) − Kan R¯ cos(θ˜n ). We treat the term Kan Im(ze−iθn ) as a stochastic perturbation to the Hamiltonian dynamics, with the important characteristic that this perturbation conserves the total energy Eq. (1). In the context of the HMF model, it has been shown that the stationary distribution of energies in the ensemble is numerically very close to a Boltzmann distribution [18]; i.e., the density of rotors ¯ 2 ), where σ 2 is the with energy E¯ is proportional to exp(−E/σ temperature. Thus, below the onset of synchronization, letting R¯ = 0, the stationary density of rotors with momentum p, given I and a, becomes 1 (p − I )2 , (6) gI,a (p) = √ exp − 2I σ 2 2π I σ 2 a Gaussian distribution with mean I and variance I σ 2 . This density will evolve from the initial density, g0 (p), over a time scale to be determined. To determine , we note that both the discrete Eq. (2) and continuum Eq. (4) systems preserve ∞ the average momentum p, which is initially P = −∞ pg0 (p)dp and becomes ∞∞∞ 0 0 −∞ pgI,a (p)h(I,a)dpdI da = I in the stationary state. Therefore, = P I −1 .
(7)
Similarly, conservation of energy determines the temperature σ 2 in terms of the variance σ02 of the original distribution of momenta g0 (p). Equating the initial energy and the energy in the stationary state with distribution Eq. (6) gives 1 . σ 2 = I −1 σ02 + P 2 1 − I I −1 We now proceed to calculate the critical coupling strength from Eq. (5) for the marginal distribution Eq. (6). Since gI,a (I ) = 0, then ω = . Inserting this in the first equation of Eq. (5) gives our main result: 1 I −1 2 2 . (8) Kc ≈ 2 2 σ0 + P 1 − a I I −1 To compare this result with the homogeneous HMF model, consider the case h(I,a) = δ(I − I0 )δ(a − a0 ). If the initial distribution g0 (p) were itself Gaussian, then Eq. (8) predicts Kc0 = 2
σ02 . a02 I0
First let us consider the effect of heterogeneity in the distribution of masses. Since the Cauchy-Schwarz inequality implies I I −1 1, we have Kc Kc0 . Thus, mass heterogeneity enhances the stability of the incoherent state. On the other hand, whenever there is heterogeneity in the charge distribution, the factor a 2 in Eq. (8) will be larger than a02 . Thus, charge heterogeneity reduces Kc , destabilizing the incoherent state. We also note that if the initial distribution g0 (p) is not Gaussian, Eq. (8) might predict a lower Kc than what would have been obtained using g0 (p) in Eq. (5).
052125-2
ONSET OF SYNCHRONIZATION IN THE DISORDERED . . .
PHYSICAL REVIEW E 89, 052125 (2014)
FIG. 1. (Color online) Estimated critical coupling as a function of width of the mass distribution, ε, for example (i) using average momenta as indicated (symbols), and corresponding Kc predicted from Eq. (8) (solid lines). Inset: estimated (Kc − 2σ02 )/P 2 compared to the prediction 2ε2 /3. The inset also includes simulations obtained using an uniform initial distribution of momenta (see text). III. NUMERICAL EXPERIMENTS
In this section, we will illustrate our results from Sec. II with three numerical experiments. In example (i) the initial distribution g0 (p) is a Gaussian centered at P with standard deviation σ0 = 0.35. There are N = 1 000 rotors with a uniform distribution of I on the interval [1 − ε,1 + ε], but all have a = 1. In this case, it can be shown that the coefficient of variation (standard deviation over mean) of the distribution of the effective frequencies is bounded by 2(σ0 /P )2 + ε2 for ε < 0.8, so we expect our theory to apply when ε 1 and P σ0 = 0.35. First, we test our prediction Eq. (8) for Kc . For a given value of P and ε, we simulate Eq. (2) by increasing K by 0.02 every 10 000 time units from the initial value K = 0, and use the last 5 000 time units for each K to estimate R by a time average. The instability threshold was estimated as the value of K at which the averaged R last exceeds 0.1. While this is a rough estimate intended to be used for relatively small N , it allows us to see how Kc varies as the parameters of the system are changed. In Fig. 1(a) we plot the estimated Kc as a function of the width ε for the seven values of average momentum P shown. The curves show the prediction ε 2 Kc (ε) = atanh(ε) σ02 + P 2 1 − ε atanh(ε) of Eq. (8). The inset shows that the computed value of (Kc − 2σ02 )/P 2 collapses onto the solid black curve, given from the theory (for small ε and σ0 ) by 2ε2 /3. Since our results depend only on the mean P and variance σ0 of g0 (p), we also consider a uniform distribution of initial momenta g0 (p) with the same standard deviation σ = 0.35 and various P . The rescaled values of Kc are also included in the inset of Fig. 1, and they collapse onto the same curve. To illustrate the validity of Eq. (6) in the incoherent state, we plot, in Fig. 2(a), the initial masses and momenta of N = 1 000 rotors using a uniform distribution of initial momenta g0 (p), and, in Fig. 2(b), their masses and momenta at T = 10 000 for K/Kc = 0.32. The solid line in Fig. 2(b) indicates the predicted mean I of the stationary distribution gI,a (p)
FIG. 2. (Color online) Evolution of the momentum-mass distribution for 1 000 rotors each with a = 1, given a uniform distribution of masses I with I = 1 and of momenta with p = P = 5. (a) Momentum-mass distribution at t = 0 with half-width ε = 0.15 in mass and standard deviation σ0 = 0.35 in momentum, and (b) for t = 10 000 when K = 0.2. The solid line is p = I with = P = 5; (c) CDFs of momenta p conditioned on I for I = 0.9, 1, and 1.1, for t = 0, and (d) for t = 10 000. Theoretical CDFs are also shown (thin red curves).
using Eq. (7). The empirical cumulative distribution function (CDF) of momenta p conditioned on I is shown in Fig. 2(c) at t = 0, for three values of I . These were calculated from the rotors with masses in three slices (I − 0.025, I + 0.025) of the distribution. The curves are simply the theoretical CDF for the uniform distribution. In Fig. 2(d) we show the CDFs at t = 10 000 for the same three values of I averaged over the time interval [2 500, 10 000]. The curves show that the corresponding theoretical Gaussian CDFs, calculated using Eq. (6), agree reasonably well with the numerical distributions. Our analysis assumes that the initial distribution g0 (p) relaxes to the stationary distribution gI,a (p); that is, the system must remain incoherent long enough for this relaxation to occur. To illustrate this, we estimate σ 2 from the simulations 2 as the mean of the instantaneous conditional variance, σI,a . According to Eq. (8), the effective critical coupling strength 2σ 2 I −1 /a 2 should approach Kc as the distribution relaxes to Boltzmann. This quantity is shown in Fig. 3 (squares), as a function of K when it is increased, again by 0.02 each t = 10 000 time units. Initially, the effective critical coupling strength is 2σ02 , but as K grows adiabatically, the distribution relaxes to Eq. (6), and the effective critical coupling strength approaches the value predicted by Eq. (8) (dashed line). When K becomes larger than this critical value, the incoherent state loses its stability, as can be seen in the evolution of the order parameter R (circles in the figure). If K is increased too rapidly or N is too large, the relaxation of the initial distribution to its steady state might not occur and the incoherent state could become unstable earlier. Indeed, in experiments with a uniform initial distribution (not shown in the figure) the relaxation is slower than when it is a Gaussian, and so K must be increased more slowly to allow the distribution to relax to its steady state. In agreement with previous observations [17], we observe that the relaxation time scales linearly with N . Thus, our results
052125-3
JUAN G. RESTREPO AND JAMES D. MEISS
PHYSICAL REVIEW E 89, 052125 (2014)
FIG. 3. (Color online) Order parameter R (circles) and estimated 2σ 2 I −1 /a 2 (squares) as K is increased for ε = 0.1 and P = 5 in example (i). The solid line indicates the stationary value of Kc predicted by Eq. (8).
apply in situations in which N is not too large or when the system is allowed to reach equilibrium (e.g., as when we slowly changed K). In addition, if Eq. (8) predicts a value smaller than Kˆ c predicted for the initial distribution, g0 (p), our computations show that when K ∈ [Kc ,Kˆ c ] the incoherent state is only metastable, since it will become unstable as g0 (p) relaxes to Eq. (6) [27]. In conclusion, example (i) has shown how heterogeneity in the masses combined with nonzero total momentum P results in the stabilization of the incoherent state: Kc is increased compared with the case of identical oscillators. Heterogeneity in the distribution of charges a has the opposite effect. For example (ii), we compare the case in which a = 1 for all oscillators to that in which a is uniformly distributed in [0,2]. The masses are again uniformly distributed in [1 − ε,1 + ε]. Numerical estimates of Kc as a function of ε are compared to the predictions (solid lines) in Fig. 4. Note that Kc is smaller when the charges are heterogeneous, as predicted by Eq. (8). For example (iii), we show a case in which the assumption of a narrow frequency distribution is violated and our analysis breaks down. Now the distribution of masses and charges is h(I,a) = 12 {δ[I − (1 + ε)] + δ[I − (1 − ε)]}δ(a − 1),
FIG. 5. Simulated (symbols) and theoretical (solid line) values of Kc as a function of ε for example (iii). The dashed line indicates the value 2σ02 for Kc predicted in the absence of disorder.
i.e., two peaks equidistant from I = 1, with the charges set to a = 1. The initial distribution g0 (p) is a Gaussian with mean P = 5 and σ0 = 0.35. Figure 5 shows the numerical estimates (symbols) versus the predicted Kc of Eq. (8) (solid line) as a function of ε. These agree when the separation between the peaks (2ε) is small; however, when ε becomes too large, the onset of instability suddenly drops near to the value 2σ02 (dashed line) predicted in the absence of disorder. In this second regime, the stationary distribution gI,a (p) for each value of I is observed to have mean P , instead of shifting to the predicted I of Eq. (6). A more detailed study of this transition is left for future research.
IV. CONCLUSION
In conclusion, we have investigated analytically and numerically the stability of the incoherent state in the disordered Hamiltonian mean field model. We found that finite size effects can induce correlations between the momenta, moments of inertia, and coupling constants of the rotors, modifying the onset of instability of the incoherent state. Indeed, heterogeneity in the moments of inertia tends to stabilize the incoherent state, while heterogeneity in the coupling strengths tends to destabilize the incoherent state. For sharply peaked parameter distributions, we developed an analytical formula for the modified critical coupling strength. Our analysis also qualitatively describes the behavior observed for broader distributions. Finally, we discovered a novel transition for a bimodal distribution of masses. Our results provide new insights into the factors affecting the phase transition in the HMF model, an iconic testbed for the study of long-range Hamiltonian systems, and provide a motivation to search for analogous results in more specific systems.
ACKNOWLEDGMENTS FIG. 4. (Color online) Estimated (symbols) and theoretical (solid line) values of Kc as a function of ε for example (ii). Inset: enlargement of the region 0 ε 0.3.
J.G.R. acknowledges a useful discussion with Ed Ott. J.D.M. acknowledges the support of NSF Grant No. DMS1211350.
052125-4
ONSET OF SYNCHRONIZATION IN THE DISORDERED . . .
PHYSICAL REVIEW E 89, 052125 (2014)
[1] P. Chavanis, J. Vatteville, and F. Bouchet, Euro. Phys. J. B 46, 61 (2005). [2] T. Tatekawa, F. Bouchet, T. Dauxois, and S. Ruffo, Phys. Rev. E 71, 056111 (2005). [3] H. Aref, Ann. Rev. Fluid Mech. 15, 345 (1983). [4] J. Weiss, A. Provenzale, and J. McWilliams, Phys. Fluids 10, 1929 (1998). [5] P.-H. Chavanis, in Dynamics and Thermodynamics of Systems with Long-Range Interactions (Springer, Berlin, 2002), pp. 208– 289. [6] J. Tennyson, J. Meiss, and P. Morrison, Physica D 71, 1 (1994). [7] D. del Castillo-Negrete and M. C. Firpo, Chaos 12, 496 (2002). [8] N. Balmforth, Comm. Nonlinear Sci. Numer. Simulat. 17, 1989 (2012). [9] A. Antoniazzi, Y. Elskens, D. Fanelli, and S. Ruffo, Euro. Phys. J. B 50, 603 (2006). [10] R. Bachelard, T. Manos, P. De Buyl, F. Staniscia, F. S. Cataliotti, G. De Ninno, D. Fanelli, and N. Piovella, J. Stat. Mech. (2010) P06009. [11] T. Konishi and K. Kaneko, J. Phys. A 25, 6283 (1992). [12] S. Inagaki and T. Konishi, Publ. Astron. Soc. Japan 45, 733 (1993). [13] M. Antoni and S. Ruffo, Phys. Rev. E 52, 2361 (1995). [14] J. Barr´e, F. Bouchet, T. Dauxois, S. Ruffo, and Y. Y. Yamaguchi, Physica A 365, 177 (2006).
[15] A. Campa, A. Giansanti, and G. Morelli, Phys. Rev. E 76, 041117 (2007). [16] T. Dauxois, V. Latora, A. Rapisarda, S. Ruffo, and A. Torcini, in Dynamics and Thermodynamics of Systems with Long-Range Interactions (Springer, Berlin, 2002) pp. 458–487. [17] V. Latora, A. Rapisarda, and S. Ruffo, Physica D 131, 38 (1999). [18] Y. Y. Yamaguchi, J. Barr´e, F. Bouchet, T. Dauxois, and S. Ruffo, Physica A 337, 36 (2004). [19] A. Ciani, D. Fanelli, and S. Ruffo, in Long-Range Interactions, Stochasticity, and Fractional Dynamics (Springer, Berlin, 2011), pp. 83–132. [20] S. De Nigris and X. Leoncini, Phys. Rev. E 88, 012131 (2013). [21] A. Bracco, J. McWilliams, G. Murante, A. Provenzale, and J. Weiss, Phys. Fluids 12, 2931 (2000). [22] J. A. Acebr´on, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler, Rev. Mod. Phys. 77, 137 (2005). [23] F. Bouchet and T. Dauxois, Phys. Rev. E 72, 045103 (2005). [24] S. H. Strogatz and R. E. Mirollo, J. Stat. Phys. 63, 613 (1991). [25] M. Y. Choi and J. Choi, Phys. Rev. Lett. 91, 124101 (2003). [26] E. J. Hildebrand, M. A. Buice, and C. C. Chow, Phys. Rev. Lett. 98, 054101 (2007). [27] A. Pluchino, V. Latora, and A. Rapisarda, Physica D 193, 315 (2004).
052125-5
