Collective Dynamics from Stochastic Thermodynamics

From a viewpoint of stochastic thermodynamics, we derive equations that describe the collective dynamics near the order-disorder transition in the globally coupled XY model and near the synchronization-desynchronization transition in the Kuramoto model. A new way of thinking is to interpret the deterministic time evolution of a macroscopic variable as an external operation to a thermodynamic system. We then find that the irreversible work determines the equation for the collective dynamics. When analyzing the Kuramoto model, we employ a generalized concept of irreversible work which originates from a non-equilibrium identity associated with steady state thermodynamics.


Introduction
Since the discovery of the fluctuation theorem [1,2,3], non-equilibrium statistical mechanics, which aims at connecting microscopic mechanics with macroscopic properties under non-equilibrium conditions, has been intensively studied.
In particular, thermodynamic concepts such as heat, work, and entropy production are seriously reconsidered so as to have a consistent thermodynamics framework for each realization of fluctuating quantities [4,5]. This framework has been referred to as stochastic thermodynamics. Owing to much effort, nowadays, it can be said that the foundation of stochastic thermodynamics has been established, and we should consider a next challenge based on the development of stochastic thermodynamics.
In the present paper, we discuss collective dynamics in systems consisting of many elements. This topic is of course one of important problems in non-equilibrium physics, but one may wonder how this problem is related to stochastic thermodynamics. Here, the first purpose of this paper is to shed light on the connection between collective dynamics and stochastic thermodynamics. A key point is that the deterministic time evolution of a macroscopic variable is interpreted as an external operation to a thermodynamic system, and the weak irreversible work is ascribed to a macroscopic friction force for the external system. The last phrase is taken from p. 192 in Ref. [4]. The combination of the friction force and the thermodynamic force gives rise to the total force. When the total force is expressed in terms of the order parameter, a differential equation of the order parameter is determined.
In section 2, we shall explain basic notions by analyzing the globally coupled XY model subjected to thermal noise. According to equilibrium statistical mechanics, the order-disorder transition point in this model is determined by a self-consistent equation for the order parameter characterizing the phase order. We then consider the time evolution of the order parameter near the transition point. Since its characteristic time scale is much longer than other variables, we interpret the time dependence of the parameter as a nearly quasi-static operation to the system. In the quasi-static limit, the so-called adiabatic theorem holds, which claims that the work is equal to the free energy change. We find that this relation is equivalent to the self-consistent equation for determining the transition point. Then, in nearly quasi-static processes, the irreversible work, which is defined as the difference between the work and the free energy change, appears slightly. Here, the irreversible work is characterized by a macroscopic friction constant. Since the irreversible work in nearly quasi-static processes is connected to fluctuations of irreversible work in the quasi-static processes, the friction constant is determined from the time correlation of a thermodynamic force at the trivial state. By calculating the friction constant, we obtain a differential equation of the order parameter.
This method is elegant but seems applicable to only thermodynamic systems. As another example of collective dynamics, in section 3, we study the Kuramoto model which is the simplest model that describes the collective synchronization [7,8]. However, there are neither thermodynamics, equilibrium statistical mechanics, nor Hamiltonian in the Kuramoto model. The situation is rather different from the globally coupled XY model. Nevertheless, when we add a noise term to the Kuramoto model, the Langevin equation for each element is similar to that of the globally coupled XY model [9]. Only difference is that there exists a non-equilibrium driving force in the Kuramoto model. Thus, from a viewpoint of stochastic thermodynamics, the analysis of the noisy Kuramoto model requires an extension of the irreversible work and the fluctuation-dissipation relation. Here comes the steady state thermodynamics of Langevin equations [10]. We already found the generalization of the irreversible work in transitions between two steady states by extending the Jarzynski equality [11] to that valid in non-equilibrium systems. By using this extended equality, we derive a formula of the friction constant in terms of time correlation functions at the trivial state. As a result, we obtain a differential equation of the order parameter near the transition point of the noisy Kuramoto model. Furthermore, we can take the noiseless limit of the equation.
It should be noted that the collective dynamics of globally coupled XY model and the Kuramoto model were studied by the so-called bifurcation analysis using a center manifold theory [7,12,13]. That is, in this paper, we do not derive new equations of the order parameters, but we present a simpler derivation method than previously known ones. In particular, if we already know the self-consistent equation, we have only to calculate the friction constant in terms of time correlation functions. The calculation is quite elementary. Furthermore, by distinguishing "static quantities" such as the free energy and "dynamic quantities" such as the friction constant, we can gasp the problem in a clear manner. Therefore, we expect that the method will be applied to systems for which the collective dynamics are not studied yet. In the last section, we argue such future problems to be studied. Throughout this paper, the Boltzmann constant is set to unity, and β is always identified with 1/T .

Equilibrium Statistical Mechanics
Let θ i (1 ≤ i ≤ N) be a phase variable of i-th element. We denote a collection of phases (θ i ) N i=1 by θ and define the Hamiltonian as The canonical ensemble of the system is given by We want to derive the equilibrium value of the order parameter defined by with r ≥ 0. We first notice that the Hamiltonian is expressed as Although r and ϕ depend on θ, we can assume that they take the equilibrium values (with probability one) in the limit N → ∞ owing to the law of large numbers. We fix r and ϕ to these values. We then write where with H one (θ; r, ϕ) = −Kr cos(θ − ϕ).
Z one (r) is the normalization constant given by The equilibrium value of r is then determined by By expanding (8) in r, we obtain The self-consistent equation (9) becomes This indicates that the transition inverse temperature β c for fixed K is given by Indeed, there are no other solutions than the trivial solution r = 0 for β < β c , while there is another solution for β > β c .

Collective Dynamics
Next, we consider the collective dynamics of the order parameter. We assume that the time evolution of θ i is described by the Langevin equation where ξ i is Gaussian-white noise satisfying The stationary probability density is the canonical distribution (2). The problem we want to solve is to obtain a differential equation of the order parameter re iϕ . In order to set the problem explicitly, we assume the probability density at the initial time t = 0 as for a given r 0 and ϕ 0 . The probability density p(θ, t) at time t is determined uniquely. Then, in the limit N → ∞, r(t) and ϕ(t) for each t take definite values for almost all θ with respect to p(θ, t). We fix functional forms of r(t) and ϕ(t) to those. Since we can rewrite (13) as the probability density at time t is expressed as where p one is given by the solution of the Fokker-Planck equation associated with (16): with the initial condition p one (θ, 0) = p can one (θ; r 0 , ϕ 0 ). Then, r(t) and ϕ(t) satisfy which is regarded as a self-consistent equation for r(t) and ϕ(t).
Here, we make a symmetry consideration. Suppose that ϕ(t) = ϕ 0 . Then, we can derive the solution as p one (θ, t) =p one (θ − ϕ 0 , t). This means that ϕ(t) = ϕ 0 is a solution of the self-consistent equation. Below, we set ϕ(t) = ϕ 0 = 0. Now, we focus on the collective dynamics near the transition point. Explicitly, we set βK = 2 + ǫ with |ǫ| ≪ 1 for fixed K. We then expect that the slow dynamics of r(t) are characterized by a scaling form where η → 0 and t → ∞ with ηt = τ fixed; andr is a function whose functional form is independent of η. We also expect that η is related to ǫ as The question is to derive an equation forr and to determine the values of a and b.
Mathematically, we have only to analyze p one near the transition point. One can apply a center manifold theory to this system. (See section 5.7 in Ref. [7].) Instead, we consider the problem from a viewpoint of stochastic thermodynamics. Hereafter, represents the expectation with respect to this initial distribution and the noise sequence ξ. We also denote the expectation of A(θ) with respect to p one (θ, t) by A t , and can r represents the expectation of the canonical ensemble with the Hamiltonian H one (θ; r).

Stochastic Thermodynamics
We study the Langevin equation (16), where r is given as a function of time. We interpret the time dependent parameter as a control by an external system, without any feedback from the system. Concretely, the force Φ done by the external system is defined as By using (7), we express r(t) determined in (19) as According to equilibrium statistical mechanics, we have where F (r) is the free energy defined by F (r) = −T log Z one (r). The self-consistent equation (9) is equivalent to Φ(r(t)) t = Φ(r(t)) can r(t) . This is valid only in the limit t → ∞, and in general cases there should be the irreversible work defined by We then obtain The problem now becomes to evaluate the irreversible work in the stochastic system. The important property here is that the time scale of r is much longer than the relaxation time of the probability density for the Langevin equation (16) near the transition point. That is, the control is assumed to be performed as a nearly quasi-static process, which enables us to develop a perturbation theory. Furthermore, owing to the recent progress on the stochastic thermodynamics, we have several identities associated with thermodynamic works. By utilizing one of them, we can simplify the calculation of the irreversible work. Concretely, we start with the Jarzynski equality [11] e −βW irr = 1.
See Appendix A as for the derivation of a generalized version of (27). By combining e −x ≥ 1−x with the identity (27), we derive W irr ≥ 0, which corresponds to the second law of thermodynamics. Furthermore, from the identity (27), in the nearly quasi-static regime η → 0, we have which corresponds to the fluctuation-dissipation relation. By taking the derivative with respect to t, we obtain with The correlation time of Φ(θ, r(s)) is controlled by T . Since dr(s)/ds ≃ O(η) is much smaller than T , (29) becomes with a friction constant Essentially the same expression was obtained in Ref. [6]. Here, the expectation is taken over samples in which θ(0) is chosen obeying the canonical ensemble with r, and θ(t) is determined from the stochastic time evolution with fixed r. By combining (31) with (26), we have This determines the time evolution of r(t) uniquely. By using the expansion (10), we rewrite (33) as Recalling dr/dt = O(η) and βK = 2 + ǫ, we find that the exponents in (20) and (21) are given by a = 1 and b = 1/2. The equation forr(τ ) with τ = ηt is in the limit η → 0 and t → ∞. The positivity of γ ensures the stability of the trivial solutionr = 0 for ǫ < 0 and the non-trivial solutionr = 1 for ǫ > 0, respectively. We also note that γ > 0 implies the monotonic increment of W irr (see (31)), which is a stronger property than the second law of thermodynamics. Finally, we calculate γ(0). From the definition of γ in (32), we have Let C(t) be cos θ(t) cos θ(0) for the free Brownian motion dθ/dt = ξ which corresponds to the case r = 0 in the Langevin equation (16). We then derive where the symbol ′′ • ′′ represents the multiplication in the sense of Strotonovich. Since C(0) = 1/2, we obtain C(t) = exp(−T t)/2. The substitution of this result into (36) yields which is evaluated to be 2 at the transition point. In sum, the differential equation for r is We guess that there should be some references reporting this result, say around 1970, but we do not find them. As far as we searched, explicit calculation was presented in Ref. [14] using bifurcation analysis. Note however that the numerical coefficient of the non-linear term in Eq. (7) of Ref. [14] is not correct.

Model
We study the Kuramoto model [7,8] where K > 0 and ω i is a time-independent stochastic variable obeying the probability density g(ω). We assume that g(ω) = g(−ω), g(ω) ≤ g(0), and the second derivative of g(ω) at ω = 0 is positive. Note that g(ω) → 0 in ω → ∞ because dωg(ω) = 1. The collective synchronization occurs when K > 2/(πg(0)). This result was obtained by the analysis of the self-consistent equation for the order parameter (3), which corresponds to (9) in the globally coupled XY model. After that, Kuramoto and Nishikawa attempted to derive the equation that describes the collective dynamics in the Kuramoto model [15]. However, it turned out that the problem was hard to be solved. Especially, even in the linear regime around the trivial state (r = 0), the analysis was far from trivial, as pointed out in Refs. [16,17]. As one remarkable result, Ott and Antonsen derived the differential equation for the collective dynamics by noting a special solution of the non-linear equation of the distribution [18]. Note that this method relies on a special property of the model [19,20] and that it cannot be applied to general cases. Quite recently, Chiba has derived the equation of the order parameter by mathematically developing a center manifold theory with a resonance pole [13].
Since the difficulty originates from the deterministic nature of the dynamics, its noisy version has also been studied, where ξ i is Gaussian-white noise satisfying (14). This model was first proposed by Sakaguchi [9]. The self-consistent equation of the order parameter in this model was analyzed and the non-trivial solution corresponding to the synchronized state was derived [9]. Then, based on the linear stability analysis of the self-consistent solutions in the noisy Kuramoto model [16,17], bifurcation analysis was performed so as to obtain a differential equation of the order parameter near the transition point [12].
(See Ref. [21] for a story related to the development.) In this section, we study the collective dynamics near the transition point for the noisy Kuramoto model from a viewpoint of stochastic thermodynamics. We then consider the noiseless limit T → 0.

Setup of the problem
We start with re-expressing (41) by We denote by p ss one (θ i ; r, ϕ, ω i ) the stationary probability density for the Langevin equation (42) with (r, ϕ) fixed. We follow the analysis in the previous section step by step.
We assume the probability density at the initial time t = 0 as for given r 0 and ϕ 0 . The probability density p(θ, t) at time t is determined uniquely. Then, for each t, in the limit N → ∞, r(t) and ϕ(t) take definite values for almost all θ with respect to p(θ, t). We fix functional forms of r(t) and ϕ(t). Then, the probability density at time t is expressed as where p one is given by the solution of the Fokker-Planck equation associated with the Langevin equation (42): with the initial condition p one (θ, 0; ω) = p ss one (θ; r 0 , ϕ 0 , ω). Then, r(t) and ϕ(t) satisfy which is regarded as a self-consistent equation for r(t) and ϕ(t). Without loss of generality, we assume ϕ 0 = 0. Since p ss one (−θ, t; −ω) = p ss one (θ, t; ω), we find from (45)that p one (−θ, t; −ω) = p one (θ, t; ω). Then, (46) leads to ϕ(t) = 0. Hereafter, ω represents the expectation over initial conditions and noise sequences in the Langevin equation (42) with the frequency ω i = ω. We denote the expectation of A(θ) with respect to p one (θ, t; ω) by A t,ω , and ss r,ω represents the expectation with respect to p ss one (θ; r, ω). Now, let K c (T ) be the transition point of the coupling constant for the model with T fixed. We characterize the distance from the transition point by Then, in the asymptotic regime |ǫ| ≪ 1, we expect that r evolves slowly and this behavior may be characterized by a scaling form (20) with (21). The question is to derive a differential equation forr together with determining the values of a and b.

Useful Identity
Differently from the previous section, thermodynamic concepts such as the irreversible work are not established for transitions between non-equilibrium steady states. Indeed, the Jarzynski equality (27) is not available for the Langevin equation (42) due to the existence of the driving force ω i . We thus need to consider an extension of the Jarzynski equality (27). This was proposed by Hatano and Sasa [10]. By defining φ(θ; r, ω) = − log p ss one (θ; r, ω), is interpreted as a generalized irreversible work in a process starting from steady state.
(See Ref. [22] for a review of an extended framework of thermodynamics on the basis of (49).) Next, since there is no Hamiltonian, we consider a different formulation from that using the thermodynamic force (22). Concretely, by defining and by recalling (3), we have The first step of the analysis is to estimate A t,ω in the nearly quasi-static regime η → 0.
Here, as essentially the same identity as (49), we derive A ss r(t),ω = Ae (See Appendix A for the derivation.) In the limit η → 0, this identity yields A similar relation was proposed by using the identity (49) [23]. Furthermore, by noting the time scale separation, we rewrite it as with Γ(r, ω) = This is a generalized fluctuation-dissipation relation claiming that the friction constant is expressed as the time correlation function.
Finally, multiplying g(ω) with the both hand sides of (55) and integrating them in ω, we obtain where we have used (52); and γ(r) and G(r) are expressed as Since G(r) = −G(−r) (see Appendix B), G(r) can be expanded as The transition point K c (T ) is determined by which is obtained from the condition that the linear term in the right-hand side of (57) becomes zero. By substituting K = K c (1 + ǫ) and r = η 1/2r (ηt) with η = |ǫ| into (57), we obtain in the limit ǫ → 0 and t → ∞, where γ(0) and a 3 are evaluated at K = K c (T ).

Noiseless limit
We consider the noiseless limit T → 0. We first rewrite a 1 as This immediately yields lim T →0 which leads to K c = 2/(πg(0)). Similarly, we obtain lim T →0 where the double prime represents the second derivative. Next, we evaluate the noiseless limit of γ(0). The method used in the estimation of a 1 and a 3 is not effective here. The heart of the calculation is to note an identity By using it, we rewrite (74) as We then obtain By substituting (77) and (80) into (62), the equation ofr becomes which coincides with Eq. (7.97) in Ref. [13]. The right-hand side corresponds to the self-consistent equation obtained by Kuramoto [7]. Furthermore, we remark This expression corresponds to Eq. (138) in Ref. [12]. Although the numerical coefficient of the latter is different from (82), there is no contradiction between the two, because |α| in Ref. [12] is equal to r/2π in this paper. (This unusual convention can be understood from Eq. (36) and Eq. (95) in Ref. [12].) More explicitly, we study a case that g(ω) is a Cauchy distribution By substituting it into the formulas given in (76), (77), and (80), we obtain lim T →0 and lim T →0 Since K c = 2∆ (that comes from (84)), (62) becomes Both the decay rate and the growth rate below and above the transition point are |K −K c |/2 in the original time scale, which are equal to the results of the linear stability analysis [16,17]. The stationary solution above the transition point is equal to the result by Kuramoto [7]. Finally, the order parameter equation presented in Ref. [18] becomes (87) near the critical point.

Concluding Remarks
In this paper, we have studied collective dynamics from a viewpoint of stochastic thermodynamics. The most important achievement is that we can obtain the order parameter equation (81) quickly. The key step in the derivation is to utilize the fluctuation-dissipation relation (54) that is derived from the non-equilibrium identity (53). Owing to this identity, we have only to calculate time correlation functions for a free Brownian particle driven on a ring, in addition to the previously known selfconsistent equations [7,9]. As is understood from the derivation method, the noiseless limit T → 0 should be taken after the scaling limit ǫ → 0 and t → ∞ is considered. When both T are ǫ are finite, our theory provides a good approximation for ǫ ≪ T . On the contrary, the calculation method cannot be applied to the noiseless Kuramoto model. Nevertheless, one may expect that the behavior for the case ǫ ≪ T ≪ 1 is close to that for ǫ ≪ 1 and T = 0. This expectation is true for some cases, but not always valid. For example, it was pointed out in Ref. [17] that when g(ω) is zero expect for [−ω 0 , ω 0 ] with some positive ω 0 , the order parameter in the noiseless Kuramoto model relaxes to the trivial state in a power law form for K < K c , which is not observed for the case ǫ ≪ T ≪ 1. We need to develop a different formulation if we want to understand the behavior of the noiseless Kuramoto model correctly [24].
Although we focus on the simplest model of coupled oscillators, one can study more general cases such that the interaction includes higher harmonics e.g. sin(θ i − θ j ) + h sin 2(θ i − θ j ). See Ref. [25] for a self-consistent equation, Ref. [26] for the analysis using the center manifold theory of the noisy case, and Ref. [27] for the generalized center manifold theory for the noiseless case. According to Ref. [27], the early attempts [25,26] have some mistakes. See Ref. [28] for a recent study. It should be noted that the symmetry property that leads to G(r) = −G(−r) and ϕ = ϕ 0 is broken for the case h = 0. This makes the calculation complicated. More importantly, it was shown that the value of the critical exponent b changes discontinuously in the noiseless limit. It would be a good problem to obtain a fresh view of this phenomenon by applying the method in this paper.
So far, we have assumed that N → ∞. In finite but large N cases, we naturally expect that small Gaussian noise is added to the equation for the order parameter. We want to theoretically derive this stochastic equation. For example, one may start with the exact stochastic equation for the distribution which is referred to as Dean's equation [29]. Writing the path-integral expression for the history of p, one may combine the WKB analysis with the techniques in this paper. It is a challenging problem to complete the formulation. See Refs. [14,31,30] for arguments on finite size fluctuations.
Obviously, the exact determination of the differential equation for the order parameter relies on the mean field nature of the model. When we attempt to study models in finite dimensions, further techniques will be necessary so as to derive the time evolution of a spatially modulated order parameter. Then, a local stationary distribution for given spatial configurations of r and ϕ should be a reference state or an unperturbed state. Although it is a highly non-trivial problem to derive the equation, we should start this analysis seriously, because we have the simplest derivation of the collective dynamics in the mean-field model. The collective dynamics of coupled oscillators defined on random networks and complex networks are also worthwhile to be studied [32,33].
Finally, we briefly mention a recent work in which the Navier-Stokes equation is derived from Hamiltonian particles systems using a non-equilibrium identity [34]. This derivation method is formally correct and the most compact in existing approaches. Simplifying calculation enables us to extract the essence of the derivation problem, and thus we can now carefully review previous studies by Mori [35], Mclennan [36], Zubarev [37], and Esposito and Marra [38], which will be reported elsewhere. However, the method in Ref. [34] involves some mathematical assumptions such as convergences of time correlation functions. Now, look into the Kuramoto model again. If we set T = 0 in the integrand of (79), the friction constant γ(0) diverges. Thus, the formal calculation does not make a sense. Similarly, in the argument of the hydrodynamic equation, we should check the well-defined nature of the dissipation constants. Maybe related to this issue, we point out that arbitrarily small noise is introduced even in mathematically deriving the Euler equation [39].
Stochastic thermodynamics formalizes thermodynamic concepts of fluctuating quantities. It is obvious from this definition that the framework is useful for analyzing small machines such as molecular motors. In addition to such direct application, universal formulas found in stochastic thermodynamics may be applied to several nonequilibrium dynamics. We hope that this paper will stimulate many researchers who work on various subjects. transition matrix T (x → y; α). That is, P ss (x; α) satisfies x T (x → y; α)P ss (x; α) = P ss (y; α). (A.1) We then define the dual transition matrix T * (x → y; α) by P ss (x; α)T (x → y; α) = P ss (y; α)T * (y → x; α).
It should be noted that x T * (y → x; α) = 1. We then have a trivial identity Here, by multiplying A(x N )P ss (x N ; α N ) to the both hand sides and taking the summation over histories (x n ) N n=0 , we obtain A ss = P ss (x 1 ; α 1 ) P ss (x 1 ; α 0 ) · · · P ss (x N ; α N ) where A for a trajectory dependent quantity A represents Note that (A.4) is also valid for Markov chains on real numbers. Next, we study the following Langevin equation for a phase variable θ: dθ dt = f (θ; α) + ξ, (A. 6) where f (θ + 2π; α) = f (θ; α) and ξ is Gaussian white noise satisfying ξ(t)ξ(t ′ ) = 2T δ(t − t ′ ). We denote the stationary probability density by p ss (θ; α). We discretize the Langevin equation with a time interval ∆t. Since this defines the Markov chain on real numbers, we have the identity (A.4). Then, by taking the limit ∆t → 0, we obtain the identity (53). When we set A = 1, it becomes the identity (49).