Nonequilibrium phase transitions in finite arrays of globally coupled Stratonovich models: strong coupling limit

A finite array of N globally coupled Stratonovich models exhibits a continuous nonequilibrium phase transition. In the limit of strong coupling, there is a clear separation of timescales of centre of mass and relative coordinates. The latter relax very fast to zero and the array behaves as a single entity described by the centre of mass coordinate. We compute analytically the stationary probability distribution and the moments of the centre of mass coordinate. The scaling behaviour of the moments near the critical value of the control parameter ac(N) is determined. We identify a crossover from linear to square root scaling with increasing distance from ac. The crossover point approaches ac in the limit N→∞ which reproduces previous results for infinite arrays. Our results are obtained in both the Fokker–Planck and the Langevin approach and are corroborated by numerical simulations. For a general class of models we show that the transition manifold in the parameter space depends on N and is determined by the scaling behaviour near a fixed point of the stochastic flow.


Introduction
Arrays of stochastically driven nonlinear dynamical systems may exhibit nonequilibrium phase transitions of continuous or discontinuous type, for a recent review see [1], cf also [2,3]. Concepts developed to describe equilibrium phase transitions such as symmetry or ergodicity breaking, order parameter, critical behaviour, critical exponents, etc, have been successfully transferred to noise induced nonequilibrium phase transitions. The structure of the theory will be generically of mean field type, if infinite globally coupled arrays are studied. This allows for a number of analytical results.
Remarkably, essential characteristics of phase transitions can already be found in the case of a single Stratonovich model. This is mainly due to the multiplicative nature of the noise. Models driven by additive noise do not show this peculiar property. The Langevin equation for the single-site Stratonovich model [4]- [6] reads where a is a control parameter, σ denotes the strength of the noise and W (t) is a Wiener process with autocorrelation W (t)W (s) = min(t, s). Equation (1) is interpreted in the Stratonovich sense as indicated by the symbol •. The Stratonovich model describes, e.g., the overdamped motion in a biquadratic potential U (x) = − a 2 x 2 + 1 4 x 4 where the control parameter is stochastically modulated, a → a + ξ t , with a Gaussian white noise ξ t .
The associated Fokker-Planck equation (FPE) describing the evolution of the probability density P(x, t) is Equation (2) has a weak stationary solution, a Dirac distribution δ(x) located at the common zero x = 0 of drift and diffusion coefficient, which is also a zero of the stochastic flow in equation (1). If the system is initially at x = 0, it will always stay there. Furthermore, there exist spatially extended strong stationary solutions determined up to a constant factor, P s (x) ∝ |x| 2a/σ 2 −1 exp{−(x/σ ) 2 }. P s (x) lives on S + = [0, ∞) if the initial distribution lives on S + \0, and on S − = (−∞, 0], if the initial distribution lives on S − \0. The 3 constant is determined such that the solution is normalized via integrating over the support and can be interpreted as a probability density, i.e. provided 2a/σ 2 > 0. For 2a/σ 2 0 the normalization Z diverges since the integrand in (4) scales like |x| 2a/σ 2 −1 as x → 0. In this case it can be shown (see appendix A) that a weakly normalized version converges to the known weak solution δ(x) and that x = 0 is an absorbing fixed point of the system.
If fractions of the initial distribution of given weights live on S − , on S + , and on 0, all will keep their weight and evolve to the stationary probability densities living on their respective support as guaranteed by a H -theorem [7].
The Stratonovich model exhibits a strong ergodicity breaking [8] depending on the control parameter a, since the state space decomposes into regions where the system cannot reach one region if it is initiated in a different one. For a 0 the only stationary solution is δ(x), i.e. the fixed point x = 0 of the stochastic dynamics is absorbing. Additionally, for a > 0 we have the spatially extended solution (3) living on S ± depending on the initial distribution. This is reflected by the mean value Obviously, x ± can serve as an order parameter and shows critical behaviour x ± ∼ ± √ π σ (a − a c (1)) β as a → a c (1) = 0 with β = 1.
Note that also the location of the maximum of the spatially extended density undergoes a bifurcation, x max ± = 0 for 0 < a a max c (1) = σ 2 /2 and x max ± = ±(a − a max c ) 1/2 for a a max c (1). The critical behaviour of an array of infinitely many globally coupled Stratonovich models has been thoroughly investigated in [9]. The scaling of higher moments was considered in [10], see also [11]. The stationary probability density is the solution of a nonlinear FPE which depends on the order parameter. The scaling behaviour of the order parameter is analytically determined, x ± ∼ ±(a − a c (∞)) β as a → a c (∞) with a c (∞) = −σ 2 /2 and β = sup{1/2, σ 2 /(2D)}, where D is the strength of the harmonic coupling between the systems [9]. The strong coupling limit, D → ∞, of an infinite array of globally coupled systems was analytically treated already in the pioneering paper [12], cf also [13,14]. For nearest neighbour coupling, the coupling term can be considered as a discretization of the Laplacian. A linear stability analysis of the corresponding Ginzburg-Landau and Swift-Hohenberg equations in the case of strong coupling is performed in [15].
In this paper, we investigate nonequilibrium phase transitions in finite arrays of globally coupled Stratonovich models in the strong coupling limit. We introduce centre of mass and relative coordinates and exploit that for strong coupling there is a clear separation of timescales. The relative coordinates relax very quickly to zero and the system behaves as a single entity described by the centre of mass coordinate R t . Thus, we can adiabatically eliminate the relative coordinates. The stationary probability density of the centre of mass coordinate p s (R) is 4 analytically calculated for a class of nonlinear systems and a scheme to determine the transition manifold in the parameter space is developed. For finite arrays of linearly coupled Stratonovich models the mean value R of the centre of mass coordinate is computed analytically. Near a critical value of the control parameter a the stochastic system shows scaling behaviour similar to the order parameter of the single Stratonovich model with the same critical exponent β = 1 but with a different a c (N ) which is also given analytically. Keeping a finite small distance to a c (N ) we recover for N → ∞ the known result of the self-consistent theory [9] with critical exponent β = 1/2, see above. For finite N we identify a crossover value of the control parameter a (N ). For a c (N ) < a a (N ) we have a linear scaling as for N = 1 whereas for a a (N ) a square root behaviour as for N → ∞ is observed. Our analytical results are corroborrated by numerical simulations.
Recently, finite arrays of (non) linear stochastic systems have been investigated also by Muñoz et al [10], and by Hasegawa [16].
Muñoz et al tried to obtain for multiplicative noise characteristics of the probability density of the mean field for finite N . They argued that the Langevin equation for the mean field variable is of similar form as the Langevin equation for a single system. Assuming that the multiplicative driving noise and the local field variable are uncorrelated, they inferred the scaling behaviour of the variance of an effective multiplicative noise with N , and of the critical value a c (N ) of the control parameter. They also predicted a crossover from a critical exponent β = 1 near a c (N ) to the critical exponents for N → ∞ for larger distances to a c (N ). Note that in [10] the Langevin equation was treated in the Ito-sense which leads to a shift of the critical control parameter compared to the same equation in the Stratonovich-sense.
Hasegawa considered finite systems with additive and multiplicative noise using his augmented moment method which is applicable for small noise strength. He emphasized that multiplicative noise and the local field variable are not uncorrelated in contrast to the assumption in [10] and demonstrated some consequences of such a simplification.
Our approach, though similar in spirit to [10], is controllable, valid in leading order for strong coupling D, and provides explicit analytic results which are confirmed by independent numerical simulations. It may serve as a starting point to calculate next order corrections ∼ 1/D. The paper is organized as follows. In the next section we consider two harmonically coupled Stratonovich models and show that for strong coupling the centre of mass coordinate R is the relevant degree of freedom. The mean value of R shows a critical behaviour which is analytically characterized. Section 3 deals with a class of N globally coupled systems of general kind. For strong coupling we compute analytically the stationary probability distribution p s (R) after eliminating the relative coordinates. Further, we determine the transition manifold in the parameter space where p s (R) undergoes a transition from a delta-distribution to a spatially extended solution. In section 4, we specialize to the case of N globally coupled Stratonovich models and determine the critical behaviour of the order parameter and of higher moments of R for strong coupling. Conclusions are drawn and a summary is given in section 5. In appendix A, we introduce the concept of weak normalization for the case that a spatially extended solution of the stationary FPE cannot be normalized in the naive sense. Appendix B shows that the Langevin approach both in Stratonovich-and in Ito-interpretation leads to the same results as the Fokker-Planck approach used in the main part of the paper.

Two coupled Stratonovich systems
We consider a pair of particles with coordinates x 1 (t) and x 2 (t) in a biquadratic potential which are coupled harmonically with positive coupling strength D and each subjected to independent Gaussian white noise of strength σ . The system of Langevin equations reads where In contrast to equation (1) no exact solution of system (6) is known. The joint probability density P( where, adopting the notation of [17], denote drift and diffusion coefficients, respectively. One can show that the system (7) exhibits no detailed balance. Hence, there is no easy way to obtain analytically the stationary solution P s .
For strong coupling, however, a systematic analytical approach is possible. With increasing coupling strength the particles become tightly glued together and move as a single entity. Therefore it appears natural to introduce centre of mass and relative coordinates. Simulations of equation (6) show that indeed the stationary distribution of the relative coordinatep s (r ) becomes very sharp for large values of D, cf figure 1.
The stationary distribution of the centre of mass p s (R) shows behaviour which is similar to the distribution of a single Stratonovich model. For large values of a, we have a monomodal distribution which vanishes at the boundaries of the support, cf figure 1. For small values of a the distribution p s (R) diverges as R → 0 in a normalizable way, cf figure 2. For even smaller values of a all trajectories x i (t) approach zero and the distributions of both r and R are δ-distributions. Accordingly, the mean value R undergoes a continuous transition at a critical value of a. In the following, we analytically calculate p s (R) and R and its scaling characteristics in the strong coupling limit, D → ∞. We introduce the centre of mass coordinate R(t) and the relative coordinate r (t) by with the inverse transformation With the Langevin equations (6) then transform to where the transformed Wiener processes W i (t) are defined as The FPE associated with (14) and (15) governing the probability density of centre of mass and relative coordinates P(R, r ; t) reads where the Fokker-Planck operators are Note that only L r depends on D. In the strong coupling limit D → ∞ the relative coordinate vanishes, r t → 0, on a very fast timescale of the order 1/D, cf equation (15). Hence, the stationary probability density factorizes to P s (R, r ) = p s (R)δ(r ) with a Dirac distribution for the relative coordinate. In this case there is no flow related to the relative coordinate r , i.e. L r P = L r R P ≡ 0, since for any suitable function ϕ Integrating equation (17) with respect to r yields in the stationary case Similarly as for the single Stratonovich model, there is always a weak solution δ(R).
for a a c (2), There is a strong ergodicity breaking when a crosses a c (2). The mean value R ± calculated with (23) is and scales like R ± ∼ ± √ 2π σ (a − a c (2)) β with β = 1 as a → a c (2).

Adiabatic elimination of relative coordinates
In the following, we demonstrate that the strategy sketched above can be generalized for a class of N coupled systems. We consider with i = 1, . . . , N and where f and g are smooth (with no singularities) and twice differentiable chosen such that the stochastic process x(t) = {x i (t), i = 1, . . . , N } has natural boundaries at infinity [4,18]. Both f and g may depend on a d-dimensional set of control parameters a. D > 0 is the coupling strength of the harmonic attraction. Note that we have absorbed a factor σ , the strength of the noise, in the function g.
The FPE for the joint probability density P(x; t), x = {x i , i = 1, . . . , N }, associated to (26) reads Using the shorthand f i = f (x i ), g i = g(x i ) and g i = ∂ x i g i , drift coefficient and diffusion matrix are given by It is advantageous to introduce centre of mass and relative coordinates {R, r}, r = {r k , k = 2, . . . , N }, by the linear transformations 9 The inverse transformation is given by x k = R + r k , for k = 2, . . . , N .
Observing the rules for linear transformations we have Drift and diffusion coefficients in the new coordinates are given by, cf also [17], where y and z stand for the new coordinates R, r k and r l , respectively. Again, the FPE determining P(R, r; t) has the form ∂ t P = LP, with L = L R + L r + L r R where Explicitly, the new drift and diffusion coefficients are All arguments in f i , g i and g i have to be expressed by (R, r), see equations (32) and (33). Note that only D r k depends on the coupling strength D explicitly. For large times the probability density P(R, r; t) converges to a stationary probability density, cf [7], determined by LP s (R, r) which has a weak solution In the strong coupling limit all fluctuations of the relative coordinates vanish. The system is concentrated on the centre of mass and moves stochastically as a whole, combined particle. The probability density of the centre of mass p s (R) can be determined by integrating P s (R, r) over all relative coordinates. Performing this integration we obtain from the stationary FPE provided that the boundary terms associated with the relative coordinates vanish. In the strong coupling limit we have P s (R, r) ∝ δ(r), (47) holds in any case and leads to where L R is given by (37). From (40) and (42) we infer drift and diffusion for r = 0 as The spatially extended strong solution of (48) is given by provided that the normalization constant Z is finite. Whether or not this is the case depends on the scaling behaviour of the functions f (R) and g(R) near a common zero R 0 which, if existing, will build a boundary of the support. This is explained in detail in the next subsection. Equation (50) shows that in the strong coupling limit the diffusion coefficient D R,R scales like σ 2 /N , cf also [10,16]. For the infinite system and finite noise strength σ the stationary probability density of the centre of mass p s (R) is a Dirac measure located at one of the attractive zeros of the drift coefficient (49), depending on the initial conditions.
For D → ∞ all particles are strongly correlated. The variance of the coordinate x i (R, r) of an arbitrary system i calculated with P s (R, r) = p s (R)δ(r) is Due to the strong correlations, the variance of the centre of mass scales like N 0 in contrast to the case of uncorrelated systems where the central limit theorem predicts a scaling like N −1 .
Note that for small coupling, D < σ 2 /2, one expects that the relative fluctuations in finite systems become large for times of order ln N ; for larger times the mean field approximation will completely break down [19].

Determination of the transition manifold
There will be a strong ergodicity breaking if the state space decomposes into different regions with the property that one region will not be accessible if we start in a different one [8].
For multiplicative noise, zeros of the stochastic flow separate the state space into mutually non-accessible regions. If we place the system initially on such a zero, i.e. on a fixed point of the stochastic dynamics, it will stay there forever. Accordingly, the FPE has a weak solution, a δ-distribution living on that fixed point. If any trajectory in the neighbourhood 'asymptotically' reaches the fixed point (the fixed point is absorbing), there will be no spatially extended probability density in this neighbourhood. The spatially extended stationary solution of the FPE cannot be normalized in the naive sense. The weak normalization procedure leads to the weak solution.
If trajectories cannot reach the fixed point, the stationary solution of the FPE will be normalizable and we will have a spatially extended probability density living on the support bounded by the fixed point. This properties can be exploited to determine the transition manifold in the parameter space.
We suppose that f (R, a) and g(R, a) near a common zero R 0 have the following scaling behaviour where m f , m g > 0. Near R 0 we have for the stationary solution (51) of the reduced FPE (48) For m f − 2m g > −1 the integral in (55) gives a contribution ∝ ε m f −2m g +1 at the upper boundary which vanishes for ε → 0 so that in leading order p s (R 0 + ε) ∝ |ε| m g (N −2) . The exponent m g (N − 2) is negative only for N = 1. In this case, if m g 1 the singularity of p s at R = R 0 is not normalizable and we have only a weak stationary solution p s (R) = δ(R − R 0 ). Note that for N 2 coupled systems of this kind the singularity of p s does not occur.
For m f − 2m g = −1 the integral gives a logarithmic contribution (A f /A 2 g )ln|ε| which leads to If the exponent in (56) is smaller than −1 the density p s (R) will diverge for R → R 0 and will not be normalizable in a naive way. The weak normalization procedure leads to a Dirac measure located at R 0 . If the exponent is larger than −1 the density p s (R) will be normalizable and we will have a spatially extended probability density. The transition manifold A c in the control parameter space R d is determined by the condition that the exponent in (56) is −1, For In the first case p s (R) is normalizable and we have a spatially extended stationary probability density. In the latter case the weak normalization procedure yields a Dirac measure at R 0 . A change in the sign of A f induced by tuning a control parameter is associated with a change of the stability of the fixed point R 0 of the deterministic flow f (R) and leads to a significant alteration of the ergodic properties. In a vicinity of R 0 the behaviour of the stochastic system is dominated by the deterministic flow. The transition manifold does not depend on the system size and the amplitude A g in contrast to (57). If the system lives on the d − 1 dimensional transition manifold (58), that is A f = 0, the scaling of the deterministic flow is not described by (53) but by f (R 0 + ε) ∼ B f ε n f with n f > m f . The system's behaviour, now depending on B f , A g , n f , m g and N , could be classified in more detail repeating the above procedure.

N coupled Stratonovich models
Now we return to our specific example and consider N globally coupled Stratonovich models in the strong coupling limit. For drift and noise function we have The common zero of f and g is R 0 = 0 with m f − 2m g = −1, and A f = a, A g = σ . Inserting this in (57), we obtain an explicit representation of the transition curve in the two-dimensional parameter space, Given the noise strength σ we have which reproduces the results for N = 1, for N = 2 (see above), and for N → ∞. Figure 3 compares results from simulation and the asymptotic result (61) for a c (N ) and illustrates the finite size scaling a c (N ) − a c (∞) = σ 2 /(2N ) for strong coupling D 1.
For initial values x i > 0 (or x i < 0 ) ∀i the stationary distribution for the centre of mass (51) lives on S + or S − , respectively, and is given by Keeping N finite, the first moment scales as a → a c (N ) like since (z) ∼ 1/z as z → 0 [20]. Note that also the higher moments R n scale linearly with a − a c (N ). 6 Trajectories of x i (t) are numerically determined by an adapted Heun method [30] with time step 0.01 for different control parameters a a c . Values of a for which x i (t) comes very close to zero are discarded to exclude a < a c . After the transient period the steady state temporal averages are built. Supposing a power law R ∝ (a − a c ) β the data are displayed for different test values of a c in a double logarithmic plot. After a visual control we determine a c by maximizing the linear correlation coefficient as proposed in [14], β is the slope of the corresponding line. 7 The transition points a c are determined observing numerically the behaviour of R(t) for an initial state where all x i (0) = R 0 with R 3 0 R 0 for a short period of time such that the distribution of the x i (t) becomes not too broad. Then, the evolution is essentially governed by the linear part of the Langevin-equation of R(t), cf appendix B.
Trajectories of x i (t)aregeneratedbyastochasticRunge-Kuttamethodwithstep0.01(seefootnote5).Generically, for a < a c (N ) the trajectory of R(t) will relax towards zero, whereas for a > a c (N ) it will increase provided R 0 was smaller than the saturation value. Finally, the estimates of a c (N ) obtained from different trajectories are averaged. A similar procedure has been exploited in [9] for N → ∞; cf also [1].  (62). On the left we show results for different system sizes N = 2 (circles), 4 (squares) and 8 (triangles) for a = 1 > a max c . The histograms are obtained from 10 6 realizations generated by a stochastic Euler method 8 . On the right we show only results for N = 8 for a c < a = −0.42 < a max c ; here 5×10 6 realizationsweregeneratedbyastochasticRunge-Kuttaalgorithm(see footnote5).σ 2 =1.
Our results are analytically derived for the strong coupling limit in a controllable approach. We note that both the critical and the crossover value of the control parameter are in accordance with the values proposed on different grounds in [10] for weak and intermediate noise, provided the shift due to the Ito interpretation used there is taken into account.

Conclusions
In this paper, we have determined the characteristics of a continuous nonequilibrium phase transition in a finite array of globally coupled Stratonovich models in the limit of strong coupling D → ∞. In this limit there is a clear separation of the timescales governing the evolution of the centre of mass coordinate and the relative coordinates: the characteristic time of the relative coordinate scales as 1/D and thus becomes short in the strong coupling limit. The slow centre of mass coordinate enslaves the fast relative coordinates, its mean value serves as an order parameter. This allows a controllable and consistent treatment both in the Fokker-Planck and the Langevin description which is corroborated by numerical simulations. The reduction of a highdimensional problem to a problem of low dimension is of course inspired by generalizations of slaving and adiabatic elimination techniques and the concept of centre manifolds to stochastic systems developed in a different context [22]- [24], cf also [25,26].
We have analytically determined both the critical value of the control parameter a c (N ) and the scaling behaviour of the order parameter and of higher moments. With increasing distance from a c (N ) a crossover from linear to square root behaviour is found. For N → ∞ the known scaling behaviour is reproduced. The formal results, i.e. the computation of the stationary distribution of the centre of mass coordinate (up to a quadrature) and the determination of the transition manifold are given for a general class of systems.
Our approach may serve as a starting point to calculate next order corrections in 1/D. In a multiscale analysis, we have to take into account that for finite but large D the distribution of the relative coordinates is, though very sharp, of finite width.
The observation that a solution of the stationary FPE which is not normalizable in a naive way converges to the weak solution if weakly normalized is certainly of interest in a broader context.
solutions which live on a support which is bounded by zeros of the stochastic flow or by natural boundaries at infinity. Under certain conditions the spatially extended solution may diverge at a zero of the stochastic flow so strongly that it is not normalizable and therefore cannot be considered as a probability density. Here, we introduce the concept of weak normalization and show that in the latter case the weakly normalized solution converges to the Dirac distribution living on that zero.
We assume that the unnormalized solutionP s (x) lives on [x 0 , b), where x 0 is a zero of the stochastic flow and scales for x → x 0 as P s (x) ∼ const (x − x 0 ) α , α < −1. (A.1) The normalization integral diverges at the lower boundary and scales like Note that the averages here are with respect to realizations ofξ 2 . We now observe that the resulting equation for R does not depend onξ 2 , therefore R = R, and obtain from which the threshold a c (2) = −σ 2 /4 follows.
The system (14) and (15) in Stratonovich sense can be written in a compact form as dρ = f(ρ)dt + j=1,2 g ( j) (ρ) • d W j (t), where ρ = (R, r ) T . The equivalent Ito system is dρ = (f + 1/2 j ∂ ρ g ( j) g ( j) )dt + j g ( j) d W j (t), where the drift term is modified by the Ito shift; ∂ ρ g ( j) is the shorthand of a Jacobian, cf e.g. [29]. For our system we have g (1) = σ/ √ 2 ρ and g (2) = σ/ √ 2 (r, R) T . The Ito shift amounts to σ 2 /2 ρ so that the equivalent Ito version of (14) reads Again, R feels only the average of the terms associated with the fast process r , we have r 2 = 0 but now also the second average vanishes since in the Ito calculus r (t)d W 2 (t) = r (t) d W 2 (t) = 0 which results in This is indeed the Ito equivalent to equation (B.3) which can be seen observing that in the single variable case the Ito shift is simply 1/2 g g = σ 2 /4 R. For arbitrary N the same procedure leads to a c (N ) = −(σ 2 /2)(1 − 1/N ) as obtained above.