Hysteresis and multistability in networks of bistable stochastic elements with global interaction

We demonstrate the existence of the first-order phase transitions and hysteresis in a network of bistable stochastic elements with global interaction subject to additive white noise. Using the Fokker–Planck equation approach, we present a method which allows one to use a continuation technique (AUTO) to follow the stationary one-particle distribution density in the space of system parameters. In addition, the Gaussian approximation is employed to compute the loci of the bifurcation points.

Networks of interacting bistable stochastic elements are of great importance for studying critical phenomena [1,2], phase transitions in spin systems [3], neural networks [4,5] and decision-making processes in the groups of individuals [6].
The state of the ith element is characterized by a scalar state variable x i (t) which changes continuously in time and has two stable states given by two minima of an external bistable potential U (x). In the simplest case, the coupling between the elements is realized through the mean field, i.e. through the ensemble average of the state variable x, x = (1/N ) i x i , where N denotes the total number of elements in the network [4,5,7,8]. The global stationary state of such a network can be characterized by the one-particle distribution density P 1 (x) which is an even function in the symmetric phase and a non-symmetric function with one well-pronounced maximum in the non-symmetric phase. If x changes continuously (discontinuously) at the point of transition, the transition is known to be of the second (first) order. Phase transitions of the first and second orders in the presence of noise can be also found in networks of interacting phase oscillators. A prominent example of the latter is the noisy Kuramoto model (for a review, see e.g. [9]). It should be noticed that the state variables x i , in the Kuramoto model is confined to the interval [0; π], whereas the state variable x of a bistable element can take arbitrary (positive or negative) values. 2 The case of mean-field coupling was studied in detail in [8,10]. It was demonstrated that the order of the phase transition depends crucially on the stochastic force. According to [8], the first-order phase transition was observed for the multiplicative noise term only. Additive noise and mean-field coupling do not induce multistability.
Here, we show that additive white noise can lead to multistability and hysteresis in a network, where interaction between the elements is global and depends on the difference between their state variables. Such interaction is often called global pair interaction. We consider a network of N bistable elements interacting globally through a pair potential W (x i − x k ). The force acting between the elements i and k is given by −∂ W (x i − x k )/∂ x i . Time evolution of the state variable x i obeys the Langevin equation in the overdamped limiṫ where ξ i (t) is a Gaussian white noise with the correlation function given by . D is the strength of the noise term and U (x) = x 4 /4 − x 2 /2 is a bistable potential. It should be stressed that the above model can also be used to describe systems of interacting particles in the overdamped limit, such as polymer and colloidal solutions [11]. The Fokker-Planck equation for the one-particle distribution function P 1 (x, t) is coupled to the two-particle distribution density P 2 (x, x , t). When the number of particles N is large, the average distance between the particles is much smaller than the characteristic length of the interparticle interaction and also smaller than the distance between the two wells of the potential energy U (x). As it was demonstrated in [12,13], in this limit (N → ∞) one can use the molecular chaos approximation and set P 2 (x, x , t) = P 1 (x, t)P 1 (x , t) which leads to a closed equation for P 1 (x, t). In the density functional theory, one refers to it as a random phase approximation [14,15]. Equation (2) must be supplemented by natural boundary conditions for P 1 (x, t), i.e. lim (x→±∞) P 1 (x, t) = 0. The presence of the non-local integral term in equation (2) makes it difficult to treat. Various approximations have been suggested to deal with the non-locality. In the approximation of strictly local coupling [16,17] when any single particle experiences influence only from its nearest neighbours but not from the rest of the network, the force −∂ W (x − x )/∂ x which acts between two elements in the states x and x , becomes strictly zero if the difference (|x − x |) exceeds a certain fixed value. In that case one can approximately reduce the integral in equation (2) The Fokker-Planck equation in the local approximation then reads This approximation has been used to study transport properties of rocking ratchets [16,17].
It is worthwhile noticing that in the local approximation [16,17], the stationary solution P loc s (x) of the Fokker-Planck equation, equation (3), with a symmetric potential U (x) is always an even function of x and, therefore, there can be no phase transitions in it. To show this, 3 we solve the stationary Fokker-Planck equation, equation (3), using zero probability current condition J = 0.
From equation (4) P loc s (x) is given by a solution of the transcendental equation where C 0 is a normalization constant which ensures that is a symmetric function and g is negative which corresponds to an attractive potential W. For any sufficiently large x, there are two solutions of equation (5). However, only one of them satisfies the natural boundary condition. This shows that the stationary density P loc s (x) is also a symmetric function of x. Therefore, the phase transitions, i.e. the transitions from a symmetric P loc s (x) to a non-symmetric P loc s (x) are not captured by the local approximation used in [16,17].
Here we show that global pair interaction does allow for the phase transitions in the network. Moreover, for a certain range of parameters, the phase transition can be of the first order, resulting in the coexistence of stable symmetric and stable non-symmetric stationary distributions P s (x). This gives rise to a hysteresis, as discussed in detail below.
Further we assume short-range interaction potential decreasing exponentially with its where α is the strength of the interaction and (1/λ) represents a characteristic 'length' of interaction. The sign in front of the interaction strength α was chosen in such a way that positive values of α correspond to an attractive potential.
Note that the case of the repulsive potential W (negative α) is rather trivial, as it does not lead to phase transitions. For negative α the stationary one-particle distribution density is broadened due to repulsion, remaining however, symmetric (not shown).
The stationary distribution density P s yields zero probability current ∂ J/∂ x = 0, i.e.
Our aim is to study the bifurcations of the stationary distribution density P s in the parameter plane (α, λ). A direct way to find P s (x) at fixed values of λ and α is to numerically solve the time-dependent equation (2) in the molecular chaos approximation. However, this method is very time-consuming as one has to wait until the initial distribution has relaxed to the stationary one. The relaxation time can be extremely large close to the bifurcation point rendering the method inappropriate. Moreover, it does not allow one to find unstable stationary solutions which play an important role in the observed phenomena. Instead, we solve equation (6) using the continuation technique. Namely, any solution of equation (6) which is known numerically or analytically at a certain set of parameters can be followed ('continued') in the space of system parameters using a collocation method and subsequent evaluation of the solution using Newton's iteration scheme. On each continuation step adaptive mesh is used to discretize the solution. This technique is implemented within the freely available software package AUTO [18], which we used here.

4
In order to be able to use AUTO [18], we need to approximate the integral in equation (6) by a certain function of x. This can be achieved by observing that the above integral is a convolution of two functions P s and ∂ W/∂ x. Both these functions are expanded in a Fourier series on the interval [−L; L] using the first M Fourier modes. With account of the fact that ∂ W/∂ x is an odd function, we obtain where ω n = 2π n/L. The length of the interval L must be chosen sufficiently large to satisfy the natural boundary conditions, i.e. P s (±L) ≈ 0. Then we can use the convolution theorem for the Fourier transforms to reduce the integral in equation (6) to Equation (6) can now be treated as a boundary value problem with a non-homogeneous term I (x) given by equation (8) and (2M + 1) integral conditions for the coefficients A 0 , A n , B n , n = 1, . . . , M To obtain the starting point for the continuation, we numerically solve the time-dependent equation (2) in the molecular chaos approximation with random initial conditions until the solution has approached a stationary distribution. To numerically solve equation (2), we use semiimplicit pseudo-spectral integration scheme as explained in detail in [19]. We take the first M = 32 modes into account and check the validity of the truncation by repeating the same calculations for the doubled number of modes M = 64. Figure 1(a) shows the stationary average x = ∞ −∞ x P s (x) dx as a function of the coupling strength α for D = 0.5 and different λ, as indicated by a number near each line. The nonsymmetric distribution corresponds to a non-vanishing x = 0, as shown by dashed and dotted lines in figure 1(b). At fixed λ, the transition from the symmetric to the non-symmetric phase is achieved by changing the interaction strength α. As we see, for λ < 5.0, the phase transition is of the second order, implying that the bifurcation is supercritical, i.e. that at the point of transition the order parameter x changes continuously. However, for larger values of λ, the bifurcation is subcritical, which results in the coexistence of stable symmetric phase x = 0 and stable non-symmetric phase x = 0 (see figure 1(b)). The unstable stationary solutions belong to the subcritical branch, i.e. to the branch between the PF point and the SN point. Interestingly the critical coupling strength α PF at the point of the (sub-or supercritical) PF bifurcation depends non-monotonically on λ. If λ is large (λ > 1) and decreases, the coupling strength α PF decreases until λ becomes small (λ < 1). After this point the further decrease of λ results in the increase of α PF , as shown by dashed lines in figure 1(a).
To understand hysteresis, we assume that λ is kept fixed at λ = 5.5 and α is being gradually increased from α = 3.4 to α = 3.8. When passing the bifurcation point at around α c = 3.65 (this point is denoted by A in figure 1(c)), the symmetric distribution becomes unstable. At this point the system undergoes a spontaneous transition to a non-symmetric phase, where the state variable x of most of the elements is in one of the two stable states. Importantly, x directly at the transition point changes abruptly as the interaction strength α increases. Conversely, if α is gradually decreased in the same limits, the system remains in the non-symmetric phase until α reaches the point of the SN bifurcation, denoted by B in figure 1(c). At this point x instantly drops to zero and the distribution function becomes symmetric. Such hysteretic behaviour is confirmed by numerical simulation of the Langevin equations, equation (1), with N = 1000 elements, as shown by symbols in figure 1(c). Notice that hysteretic behaviour can also be achieved by cooling or heating the system, i.e. by gradually changing the noise intensity D with time. This was demonstrated in [20] for a network of bistable stochastic elements with mean-field coupling. Hysteresis occurring when the parameter is being slowly changed across the bifurcation point is known in bifurcation theory and is associated with dynamical bifurcations.
A slightly larger hysteresis region obtained from the Langevin equations, equation (1) as compared to equation (6) can be explained as follows. The symmetric state becomes only weakly unstable beyond the PF bifurcation point, and, if one sets initial conditions close to this state, the characteristic growth time of the instability appears to be larger than the observation time accessible in the simulations.
To study the phase transitions analytically, we employ the Gaussian approximation and assume that the time-dependent distribution density is given by P( , where m and σ are the average and the variance of the state variable x, respectively. Using this ansatz and the time-dependent Fokker-Planck equation equation (2), we derive the time-evolution equations for m and σ With W (x − x ) = −αexp(−λ|x − x |) and P(x, t) Gaussian, the integrals in equations (10) can be solved analytically to yielḋ It should be emphasized that the bifurcation diagram obtained in the Gaussian approximation can be compared only on the qualitative level with the exact results. However, this approximation proved to be a useful and simple tool for studying phase transitions in networks of interacting elements, as has been demonstrated by many authors [4,5], [21]- [24]. From figure 1(b), it is clear that the shape of the distribution density P s (x) in the nonsymmetric phase is closer to Gaussian than it is in the symmetric phase. Therefore, we expect that the Gaussian approximation predicts more accurately the point of the SN bifurcation than that of the PF bifurcation as the former occurs at finite (m = 0).
In figure 2(a), the locus of the SN as well as the PF bifurcations is shown on the parameter plane (α, λ) at fixed noise intensity D = 0.5. Figure 2(b) demonstrates how the position of the SN line changes with the noise intensity D, which is given by a number near each line. With increasing D, the critical coupling strength at the point of the SN bifurcation α SN grows, showing that it becomes harder for the system to undergo the phase transition.
We would like to emphasize that the case of an exponential interaction potential considered here represents an intermediate situation between two limiting cases. The case of the mean-field coupling is formally obtained by setting the interaction length 1/λ to infinity (λ = 0), whereas the local approximation [16,17] is recovered for large values of λ. From figure 2(a) we can see that as λ decreases (limit of the mean-field coupling), the two bifurcation lines (SN) and (PF) collide which implies that multistability disappears and the bifurcation becomes supercritical. This is in agreement with the previous findings for the mean-field coupling [4,5,8].
The opposite case of large λ is more intriguing. When D is fixed and λ is being increased, the critical coupling strength of the SN bifurcation α SN remains finite. This shows that for a very short-range interaction (large but finite λ), the phase transition is still possible, and it is of the first order. In other words, if a sufficiently large (macroscopic) number of elements in the network will be forced into one of the two stable states, the rest of the elements will follow this trend, resulting in a phase transition to a new stable non-symmetric phase.
However, if interaction is localized to the extent that it completely vanishes at finite distances (x − x ) and is thus not described by any finite λ however large, as in [16,17], no phase transition occurs at all. The transition between intermediate global and strictly local coupling can be mathematically associated with the hysteresis in figure 1(a) becoming more pronounced until the inner and outer SN lines merge and disappear.
To conclude, we have shown that a non-local pair interaction between bistable elements driven by Gaussian white noise enables both first-and the second-order phase transitions in the network. At the point of the first-order transition, the symmetric one-particle distribution density becomes unstable leading to a spontaneous synchronous transition of the majority of the elements to a state with nonzero average x . Multistability is observed for small but finite values of the interaction length 1/λ. For the values of the interaction strength α between α SN and α PF , the stable symmetric phase coexists with a stable non-symmetric phase, which gives rise to a hysteresis.