Non-equilibrium thermodynamics and Phase transition of Ehrenfest urns with interactions

Ehrenfest urns with interaction that are connected in a ring is considered as a paradigm model for non-equilibrium thermodynamics and is shown to exhibit two distinct non-equilibrium steady states (NESS) of uniform and non-uniform particle distributions. As the inter-particle attraction varies, a first order non-equilibrium phase transition occurs between these two NESSs characterized by a coexistence regime. The phase boundaries, the NESS particle distributions near saddle points and the associated particle fluxes, average urn population fractions, and the relaxational dynamics to the NESSs are obtained analytically and verified numerically. A generalized non-equilibrium thermodynamics law is also obtained, which explicitly identifies the heat, work, energy and entropy of the system.


I. INTRODUCTION
The classic Ehrenfest model [1] was introduced to solve the reversal and Poincare recurrence paradox [2,3]. It describes a total of N particles in two urns that can randomly jump from one urn to the other with equal probability. The system has a Poincare cycle time of 2 N [4], providing a fundamental link between reversible microscopic dynamics and irreversible thermodynamics.
Later on, a directional jumping rate between urns was introduced [5,6], and extensions to multi-urns [7][8][9] were made. Although there are various modifications [10][11][12][13][14][15][16][17][18][19][20][21] of the classic Ehrenfest model, or even extensions by incorporating non-linear contribution [22][23][24][25], particles do not interact or the interaction is merely phenomenological. In fact, pairwise particle interaction in the same urn has been considered recently in the twourn model [26]. The interacting two-urn Ehrenfest model can exhibit phase transitions by varying the interaction strength and directional jumping rate from which the relaxation time and Poincaré cycle can be derived. The multi-urn system with interaction and unbiased directional jumping will evolve to equilibrium and has been shown to exhibit different levels of non-uniformity emerge with the coexistence of uniform and non-uniform phases [27]. Contrary to the better understood equilibrium cases, non-equilibrium statistical physics remains challenging, partly due to the lack of well-characterised states. Even for non-equilibrium steady states (NESS), it is difficult to describe non-equilibrium phase transitions between different NESS and their relationship to some microscopic models. For example, selection rules, such as maximal or minimal entropy production principles [28][29][30] have been proposed to determine the non-equilibrium states. Yet, universal guiding principles are still lacking.
In this manuscript, we consider urns with intra-urn interactions connected in a one-dimensional ring with di-rectional jumping rates. We will show that the system has non-equilibrium steady states in uniform and nonuniform phases that can coexist. For high directional jumping rates with appropriate interaction strengths, the steady states become unstable. The phase diagram will be obtained analytically, and the relaxation dynamics to the NESS will be studied. The relationship between nonequilibrium thermodynamical variables, such as the internal entropy production rate, the rate of work done applied to the system, will be shown to obey a generalized thermodynamic law. Finally, we will demonstrate that, in the coexistence region, the internal entropy production rate fails to select the favorable steady state.

II. EHRENFEST URNS IN A RING
We consider three urns as illustrated in Fig.1 as this already captures the non-equilibrium physics. The state of the system is labeled by n = (n 1 , n 2 , n 3 ) where n i is the number of particles in the i-th urn with n 1 +n 2 +n 3 = N , the fixed total particle number. Similar to previous models [26,27], we include a pairwise attractive (repulsive) interaction with negative (positive) energy J for particles in the same urn. Particles in different urns do not interact. A particle in the i-th urn (initial state n) jumps to the j-th urn (final state m) with corresponding transition probability where m i = n i − 1 and m j = n j + 1. g ≡ N Jβ where β is the inverse of effective temperature. Without interaction (g = 0), we have T m, n = 1 2 . Next, a jumping rate is introduced such that the probability of anticlockwise (clockwise) direction is p (q). For the sake of convenience, p+q = 1 is imposed for which only changes the time scale. Schematic diagram of the model. Three urns with particle numbers n1, n2, and n3 are connected in a ring. The direct jumping rate in anticlockwise (clockwise) direction is p (q). Ki→j represents net particle flow rate from the i-th to the j-th urn.
After s steps from the initial state, the probability of the state n is denoted by ρ( n, s). The master equation from the s-th to the (s + 1)-th time step can be written as holds if the particle jumps from the i-th to the j-th urn is anticlockwise, i.e., (i, j) = (1, 2), (2, 3), (3,1), and if (i, j) = (2, 1), (1,3), (3,2), which represents clockwise jumps. Finally, the net particle flow rate from the i-th to the j-th urn (see Fig.1 for illustration) is given by where m i = n i − 1, m j = n j + 1, and m ′ i = n i + 1, m ′ j = n j − 1.

III. EQUILIBRIUM STATES
If lim s→∞ ρ( n, s) exists, it defines the steady state ρ ss ( n). Taking the limit s → ∞ in Eq.(2) would lead to an equation in which ρ ss ( n) satisfies. No closed form exists in general. But for the case of p = q = 1 2 , an analytic expression for ρ ss can be obtained up to a normalization constant. This steady state is also the equilibrium state, because it satisfies the condition of detailed balance W m, n ρ ss ( n) = W n, m ρ ss ( m) (7) which can be verified by direct substitution. Results for the general case of M urns at equilibrium can be found in Ref. [27]. For the three urns case here with the constraint n 1 + n 2 + n 3 = N , one can define the population fraction with x i ≡ n i /N and rewrite ρ ss in terms of The saddle point approximation [31] gives asymptotic form where x sp is the saddle point(s) satisfying The solutions at different coupling constant g are shown in the data of p = 0.5 in Fig.2(a). The steady net particle flow at equilibrium from Eq.(5) reads which gives K ss 1→2 = K ss 2→3 = K ss 3→1 = 0 from Eq.(10), in large N limit. At equilibrium, there is no net particle flow between any two urns. In addition, it is also easy to see from Eqs.(5) and (7) that for p = q = 1 2 , all fluxes K ss i→j = 0 at equilibrium. On the other hand, or the nonequilibrium case of p = q, there can be non-vanishing circulating fluxes as in other general NESS systems [32][33][34].

IV. UNIFORM AND NON-UNIFORM NON-EQUILIBRIUM STEADY STATES
So far, the recurrence relation in Eq.(2) cannot be solved analytically even for NESS. In this section, we will transform Eq.(2) into the Fokker-Planck equation. Let the (physical) time t = τ1 N s, where τ 1 is the time scale of each single step from s to s + 1. Replace ρ( n, s + 1) − ρ( n, s) by ρ( n, t + τ1 N ) − ρ( n, t) = τ1 N ∂ρ ∂t + O(( τ1 N ) 2 ). Eq.(2) can be rewritten as which is known as the Kramers-Moyal expansion [35,36]. From now on, we take τ 1 = 1 for convenience.
If we further keep terms up to O(1/N 2 ), Eq.(13) becomes the Fokker-Planck equation where The WKB approximation [37][38][39][40][41] yields the saddle points [31] x sp = (x sp 1 , x sp 2 ) as whose solutions at different g and p are shown in Fig.2. The physical meaning of Eq.(20) is that K ss 1→2 = K ss 2→3 = K ss 3→1 = K ss , i.e., a constant non-zero cyclic flux of net particle along the ring which can be calculated from Eq.(5) as For uniform NESS (x 1 = x 2 = 1 3 ), one obtains K ss u = N 6 (p − q) and, for non-uniform NESS, K ss nu can be computed using the non-uniform saddle point from Eqs. (20). The K ss u and K ss nu NESS fluxes as a function of g at different p are shown in Fig.3 indicating that the particle flux of the uniform NESS is always significantly larger than that of the non-uniform NESS. Notice that there is a coexistence region of uniform and non-uniform saddle points.
The NESS can be further analysed by expanding around the saddle point, i.e. using where the matrix c is uniquely determined by the Lyapunov equations ac −1 + c −1 a t = 2b (See Appendix A for details). The detailed balance condition can be trans- is always a saddle point of uniform population fraction. At this uniform NESS, we have which gives Up to the cyclic permutation, we consider x sp 1 the largest fraction, and x sp 2 is the occupation fraction in the next urn along the p direction. The remaining Hence its stability requires g > −3. The stable region for the non-uniform phase can be determined analytically in a similar manner and the results agree with the analytic ones of the phase boundary in Fig.4. In fact, the non-equilibrium physics of the system can be summarized by the phase diagram in Fig.4 whose phase boundaries can be determined analytically (see Appendix C for derivations). As the particle interaction is strongly repulsive (positively large g), the particles are uniformly distributed in every urn. On the other hand, the particles "prefer" to stay in the same urn (non-uniform distri-bution) if they are strongly attractive (negatively large g). In between, for low jumping rate, 1 − p s < p < p s (p s ≃ 0.6823), there is a coexistence region where both uniform and non-uniform distribution are locally stable.
There is a first order non-equilibrium phase transition between the uniform and non-uniform NESSs whose transition value of g can be determined from the analytic result of K ss u and K ss nu (see Appendix C for details) together with the results of mean steady state flux. The latter has to be determined numerically using (2) . The first order phase transition line (dashed-dotted curve) is also shown in Fig.4. It is close to the stability line of the non-uniform NESS (for magnification, see Fig.12 in Appendix C). As the jumping rate becomes higher, i.e. p > p s (or p < 1 − p s ), the system is far from equilibrium and steady states do no longer exist. There are four regions which represent uniform NESS (particles uniformly distributed in three urns), non-uniform NESS, coexistence (both uniform and non-uniform NESSs are stable), and no steady state. The stability boundaries of the uniform NESS and non-uniform NESS are denoted by the vertical dashed line and solid curve respectively. The first order phase transition boundary in the coexistence regime is denoted by the dashed-dotted curve (see Fig.12 in Appendix C for a magnification).

V. RELAXATION TO STEADY STATES
In this section, we studied how the system evolves to its steady states. When the system is initially away from its steady state, it will relax (thermalized) towards the nearby stable NESS whose dynamics near the NESS is determined by the eigenvalues of a. From Eq. (14), and expand around the saddle point x sp , we get the evolution of the average of particle numbers in the first and second urn as d dt At uniform phase, x sp 1 = x sp 2 = x sp 3 = 1 3 , its solution is describing the decaying process with oscillation, with the relaxation time τ R = 4 g + 3 (28) and the oscillation period x 2 (t) and x 3 (t) can be obtained by making cyclic transformation to Eq. (27). Near equilibrium (|p − q| ≪ 1), τ osc ≫ τ R , then the solution is simplified as and the damped oscillation is not prominent. Fig.5 shows the relaxation towards the uniform NESS for different p, with g = −2.5 at which only the uniform NESS is stable. By increasing p from 0.5 to 1.0, it can be seen from the direct numerical calculation by Eq. (2), starting from the non-uniform initial state x(0) = (1, 0, 0), the relaxation time keep almost unchanged and the oscillation in the occupation gradually appears. When g = −3.5, the system stays at the non-uniform phase. To quantify the degree of non-uniformity, we define as the "non-uniformity" parameter [27]. It is almost zero for uniform phase and becomes larger for higher non-uniformity. The relaxation towards the non-uniform NESS is illustrated by ψ(t) with g = −3.5 in Fig.6, showing a pure relaxation behavior. Starting from the uniform initial state x(0) = ( 1 3 , 1 3 , 1 3 ), the system for different p all relax to the non-uniform NESS and saturates at a high non-uniformity. This can be understood in terms of the eigenvalues of the matrix a at the non-uniform state which are always real and negative as shown in Fig.7.

VI. NON-EQUILIBRIUM THERMODYNAMICS
We further examine the relationship between various thermodynamical quantities, first for general nonequilibrium states and then for the NESS cases. The Boltzmann (Gibbs, Shannon) entropy of the system is given by where the multiplication factor N ! n1!n2!n3! is due to the degeneracy of ρ( n, t) [2]. Applying Eq.(2), the entropy production rate above can be written as (W n, m ρ( m, t) − W m, n ρ( n, t)) log W n, m ρ( m, t) W m, n ρ( n, t) where the first term is the internal entropy production rate [42], which is positive-definite and only vanishes when the system is at equilibrium (Eq. (7)). It is the entropy produced during the irreversible process [43]. The second term refers to the entropy production rate for the reversible process [44] into the system. In the following, we will show that where dE dt and dW dt are the rate of change of system energy and the rate of work done by the system, respectively. From the first law of thermodynamics (conservation law of energy), T d e S can be identified as dQ, the heat flow to the system from the environment. In general, during thermalization process, dS ≥ βdQ = βdE + βdW .
Using Eqs.(3)-(4), when the particle jumps from the i-th to the j-th urn, corresponding to the transition from state n to m, where ac (c) stands for anti-clockwise (clockwise) direc-tion. The first term is the rate of change of energy β dE dt and the second term is equal to the rate of work done which can be written as where µ ≡ β −1 log( p q ) is the effective chemical potential difference to actively drive the particle from one urn to another, and the natural boundary condition is assumed. Here Eq.(34) is proved, and hence a more general thermodynamic law is derived. Note that Eqs. (38) and (39) hold even for general non-equilibrium (non-steady state) processes. For p = q, the system is at equilibrium and dS dt = diS dt = dQ dt = dE dt = dW dt = 0. That is, all macroscopic thermodynamic quantities do not change with time. For p = q under NESS, since the system energy and entropy are functionals of the probability distribution and hence are time independent, thus one has dS dt = dE dt = 0. Using Eqs.(38)- (39), which is a positive constant, corresponding to the housekeeping heat production rate to maintain the NESS. All the work done (−dW ) to the system is dissipated (measured by the internal EP d i S) into heat energy (−dQ). Furthermore, the more non-uniform is the NESS (Fig.8), the less is the particle flow (Fig.3), and hence the less internal EP (Fig.9). Since the internal EP for the nonuniform phase is always lower, the maximal EP principle could not be used to select the favorable state in the coexistence region. The first order non-equilibrium phase transition between the uniform and non-uniform NESS can also be observed by examining the internal EP rate as the particle attraction varies. Fig.9 shows a sharp drop near some threshold as g decreases and signifying a first order transition from the high internal EP uniform NESS to the low internal EP non-uniform NESS. Interestingly, there is a connection between the internal EP rate and the non-uniformity in the NESS. As shown in Fig.10, when the relationship between 1 3 − ψ ss and 1 N diS dt /(p − q) log( p q ) are plotted, all data with different p are collapsed into a single curve. It implies the relation where the function Φ(ψ ss ) is some decreasing function, i.e., Φ ′ (ψ ss ) < 0. To have higher internal EP rates, the system should be more uniform (lower ψ ss ), or with a higher direct jumping rate (higher p). Fig.10 agrees with the conjecture that more uniform states have higher EP.

VII. SUMMARY AND OUTLOOK
In this paper, we extended the Ehrenfest urn model with interactions and directional jumps, allowing for detailed investigations of the non-equilibrium steady states and associated thermodynamics. We showed that the model provides different kinds of equilibrium and nonequilibrium states. Albeit simple the model may serve a convenient paradigm system to investigate a variety of statistical physics phenomena, ranging from equilibrium to NESS and even far from equilibrium situations.
In some situations, Landau type free energy can be constructed for NESS near a continuous phase transition [45,46] or near the saddle-point(s) of NESS states [32]. However, it is still highly non-trivial to construct or establish the existence of NESS free energy in general, especially in our case of coexisting NESS related by first order transitions. On the other hand, because of the existence of probability density ρ ss ( x) at steady states, one may define the corresponding effective potential function . This NESS variable may reveals some NESS physical properties. Its numerical solution will be reported elsewhere.
At high direct jumping rate and moderate coupling constant, the system is far from equilibrium and cannot attain the steady state but limit cycle emerges instead. If the number of urns is more than three, chaotic behavior may be possible. Such models open the possibilities of investigating systems with different degree of non-equilibrium systematically. Detailed investigations will be reported elsewhere.
The multivariate linear Fokker-Planck equation for steady state reads where a and b are constant matrices of dimension D × D. b is symmetric. The natural boundary condition, ρ ss ( x)| | x|→∞ = 0 and ∂ i ρ ss ( x)| | x|→∞ = 0, is imposed. The steady state was already known [47]. In the following, we briefly outline the solution.
The form of the solution is Gaussian, where c is a symmetric matrix determined by a and b. Substitute this form into Eq.(A.42), we get two constraints, for any vector x. Eq.(A.44) is redundant (See below). Notice that cbc is symmetric but ca is not necessary to be. x t (ca)x = x t (a t c)x since they are both numbers. Hence Eq.(A.45) could be rewritten as then the steady state is also called the equilibrium state. W ( x; t + dt| x ′ ; t) is the transition rate from one state x ′ at time t to another state x at time t + dt [48], which can be expressed as Notice that Eq.(A.49) is also equivalent to Hence the uniform NESS state is stable for g > −3 (tr(a) < 0 and det(a) > 0). In addition, the eigenvalues of a at the uniform NESS state can be easily calculated to give and hence the relaxation to the uniform NESS state always has an oscillatory component.
The behavior of the equilibrium case of p = 1 2 is given in details in Ref. [27]. For the non-equilibrium case of p = 1 2 , uniform, non-uniform NESS and their bistable coexisting states also occur. The phase boundary p c (g) is determined in a similar way by the condition of saddlenode bifurcation of a pair of stable and unstable saddlepoints. For a given g, p c (g) is given by the solution of the following three equations for the three unknowns p c , x 1 and x 2 . The condition of saddle-node bifurcation in Eq.(A.61) can be written out explicitly as (e gx 1 +e gx 2 ) 2 The phase boundary of p c (g) for the saddle-node bifurcation together with the g = −3 line for stable uniform NESS are shown in the phase diagram ( Fig.3 in main text), classifying the dynamics of the 3-urns model into four regimes. In the region of g < −3 and p > p c (g), there is no stable NESS state with non-steady dynamics and the system is far from equilibrium. Furthermore, the first order transition line in the coexisting regime can be analytically determined as follows. The steady state distribution near the NESS can be approximated by the Gaussian form exp(N δ x t cδ x), where δ x ≡ x− x sp . In the coexisting regime, denote the relative weights of the uniform and non-uniform NESSs by f (g) and 1 − f (g), respectively, where the dependence on g is written out explicitly. Then the steady state distribution can be expressed as which can be analytically calculated as a function of g. Thus by numerically computing K ss (g) using the numerical solution of Eq.(2), g t can be obtained from the intersection of the curves of K ss (g) and Eq.(A.67). For given values of p in the coexistence regime, g t is obtained theoretically from the above manner and the result of the first order transition line is shown in Fig.12. The first order line is rather close to the stability boundary of the non-uniform NESS indicates that the non-uniform NESS dominates over the uniform NESS in the coexistence regime. This echoes with the observation in Fig.7 that the eigenvalues of a of the non-uniform NESS are much more negative than that of (the real part) the uniform NESS unless p is very close to the stability boundary, indicating that the non-uniform NESS is a strong attractor than that of the uniform NESS in most of the coexistence regime.