The Lee-Yang theory of equilibrium and nonequilibrium phase transitions

The Lee-Yang theory of equilibrium and nonequilibrium phase transitions

R.A.Blythe^I; M.R.Evans^II

^IDepartment of Physics and Astronomy, University of Manchester, Manchester M13 9PL, UK

^IISchool of Physics, University of Edinburgh, Mayfield Road, Edinburgh EH9 3JZ, UK

ABSTRACT

We present a pedagogical account of the Lee-Yang theory of equilibrium phase transitions and review recent advances in applying this theory to nonequilibrium systems. Through both general considerations and explicit studies of specific models, we show that the Lee-Yang approach can be used to locate and classify phase transitions in nonequilibrium steady states.

I Introduction

In this work we seek a mathematical understanding of phase transitions in the steady state of stochastic many-body systems. Systems at equilibrium with their environment provide examples of such steady states, and the mechanisms underlying equilibrium phase transitions are long known and understood. Experimentally, one can distinguish between two types of transition: the first-order transition at which there is phase coexistence, e.g.between a high-density solid and a low-density fluid, and the continuous transition at which fluctuations and correlations grow to such an extent as to be macroscopically observable.

From a thermodynamical perspective, one can understand first-order transitions by associating with each phase a free energy. For a given set of external parameters, the phase 'chosen' by the system is that with the lowest free energy, and so a phase transition occurs when the free energies of two (or more) phases are equal. The sudden changes in macroscopically measurable quantities that take place at first-order transitions are described mathematically as discontinuities in the first derivative of the free energy. Discontinuities in higher derivatives relate to continuous (higher-order) phase transitions.

The tools of equilibrium statistical physics allow one-in principle at least-to express the free energy solely in terms of microscopic interactions. More specifically, the free energy is given by the logarithm of a partition function, a quantity that normalises the steady-state probability distribution of microscopic configurations. Initially it was not universally accepted that this approach could faithfully describe phase transitions, in particular the first-order solid-fluid transition [1]. In order to show that the statistical mechanical approach can reproduce the correct discontinuities in the free energy at a first-order transition, Lee and Yang [2] introduced a description of phase transitions concerning zeros of the partition function when generalised to the complex plane of an intensive thermodynamic quantity. Initially, the theory was couched in terms of zeros in the complex fugacity plane which is appropriate for fluids in the grand canonical (fluctuating particle number) ensemble. However, by mapping a lattice gas onto the Ising model [3], the theory was found to hold equally well for a magnet in a fixed magnetic field. Later [4-6] it became clear that the distribution of zeros in the complex temperature plane can reveal information about phase transitions in the canonical ensemble.

In section II we present a brief, self-contained discussion of the Lee-Yang theory of equilibrium phase transitions which relates nonanalytic behaviour in the free energy to zeros of the partition function. In the remainder of this article, we move our attention to phase transitions in nonequilibrium steady states, that is, those that carry currents of mass, energy or some other quantity. We give therefore a very brief introduction to the established practice of modelling nonequilibrium physics through stochastic dynamics in section. In particular, we will identify a candidate quantity for use in place of an equilibrium partition function when performing an Lee-Yang analysis of phase transitions. Finally, in section IV we review recent work in which the Lee-Yang theory has been successfully applied to specific nonequilibrium phase transitions.

II Overview of Lee-Yang theory of equilibrium phase transitions

In this section we recount the relationship of zeros of a partition function to phase transitions in the system it describes. We present here a brief but self-contained tour of the key points. For more detail (and mathematical rigour) the reader is referred to textbooks, such as [7] or, for a more advanced presentation, [8], as well as early treatments of the subject[2-6].

For concreteness we shall describe the theory in terms of a simple model system of N spins in thermal contact with a heat reservoir. In this model system the energy of the system can take values E = n with n = 0, 1, 2, ¼, M and the number of microstates corresponding to the n^th energy level is g(n). The canonical partition function is given by

where b is the inverse temperature. To investigate the zeros of this partition function, it is useful to make the change of variable z = e-b. Then, the partition function is explicitly a polynomial in z of degree M and can be factorised as

Clearly the quantities z_n are the M zeros of the partition function Z_N(z); meanwhile, k is a constant which we can safely ignore in the following.

Since all the g(n) are positive, no zero of Z_N(z) can be real and positive: that is, the zeros z_n will generally lie in the complex plane away from the physical values of z which lie on the positive real axis. We now define for all complex z except the points z = z_n the (complex) free energy

Using (2) we rewrite this as

and note that a Taylor series expansion of h_N(z) around a point z ¹ z_n has a finite radius of convergence given by the distance to the nearest zero from z. This then implies that that h_N(z) can be differentiated infinitely many times over any region of the complex plane that is devoid of partition function zeros. Since we identify a phase transition through a discontinuity in a derivative of the free energy, we see that such a transition can only occur at a point z₀ in the complex plane if there is at least one zero of the partition function Z_N(z) within any arbitrarily small region around the point z₀.

Clearly this scenario is impossible if the number of zeros M is finite, except at the isolated points z_n where the free energy exhibits a logarithmic singularity. Since such a point cannot lie on the positive real z axis, there is no scope for a phase transition in a finite spin system, such as the simple example (1). On the other hand, if the partition function zeros accumulate towards a point z₀ on the real axis as we increase the number of spins N to infinity there is the possibility of a phase transition.

In order to deal with the thermodynamic limit (see [8] for rigorous considerations) we shall assume that the limit

exists and we may write

where r(z) is the local density of zeros in the complex-z plane in the thermodynamic limit.

Since the imaginary part of this complex free energy h(z) is multi-valued it will at times be more convenient to work with the potential j(z) defined as

An expression for the density r(z) in terms of the potential j(z) is easily obtained once one recognises that ln|z| is the Green function for the two-dimensional Laplacian. Specifically, if z = x + iy, we have

Using this expression, we can take the Laplacian of both sides of (7) to find

Such an equation is familiar from electrostatics and so j(z) is analogous to an electrostatic potential.

In analogy to electrostatics, as long as we can integrate r(x) over bounded regions containing any singularities of r, the potential will be a continuous function [9]. The significance of this statement is that we can derive a rule for locating phase boundaries given a partition function. Let us suppose that around points z₁ and z₂ in the complex plane, one has (analytic) expressions for the potential j₁(z) and j₂(z) such that j₁(z) j₂(z) (i.e.not the same function over the entirety of the complex plane). In order for the potential to be continuous at all points on the complex plane, we must have a phase boundary at those values of z for which the condition

holds. This is basically the definition of a transition mentioned above in the introduction. Since j₁ and j₂ are different functions, some derivatives of j(z) will not exist at these values of z and we expect the density of zeros at these points to be non-zero. (It is also possible for zeros to be present at other points in the complex plane, but we do not need to consider this possibility here.)

Typically, a solution of (10) describes a curve C that intersects the positive real z axis at a point z₀. Having already established that we require zeros to accumulate at the point z₀ in the thermodynamic limit for a physical phase transition to take place, we are interested in the line density of zeros per unit length of this curve. We introduce the arclength s which measures the distance along C from the transition point z₀.

To obtain an expression for the line density of zeros m(s) along the curve C, consider a short section of C of length ds enclosed by an area A that has two sides parallel to C and the other two sides perpendicular to C–see Fig.1. Integrating the density r(z) over this area we have

Meanwhile, we have from (9) and an application of the divergence theorem that

in which the functions j₁(z) and j₂(z) relate to the limiting values of j(z) as the point z on the curve C is approached from either side, and the vector is the unit vector normal to C at that point.

Recall that j(z) is the real part of the complex free energy h(z) and therefore, away from any zeros of the partition function, satisfies the Cauchy-Riemann equations

in which y(z) is the imaginary part of h(z) and z = x + iy. Then we have that

where is the unit vector tangent to the curve C at the point z. Thus we recognise the scalar product in (12) as the directional derivative of the imaginary part of the free energy h(z) along the curve C. Putting this together with (11) we find that

in which y(z) = Im [h(z)] and the subscripts indicate opposite sides of a phase boundary.

Let us now assume that in the thermodynamic limit, there is a phase transition at the point z₀. Then, either side of the transition the free energy can be written as

where = z – z₀ and f₁() is valid for Re < 0 and f₂() valid for Re > 0. Note that for the free energy to be real along the real z axis, the coefficients b and c must also be real. From the criterion that the real part of h(z) is continuous across a phase boundary, we find that the boundary lies along the curve

where and are the real and imaginary parts of respectively. We consider the conditions under which a discontinuous (first order) or continuous (second or higher order) transition appear.

First-order transition

For the case where b₁¹ b₂, and the free energy has a discontinuity in its first derivative, the curve of zeros is a hyperbola that passes smoothly through the transition point z₀. Hence the tangent to the curve of zeros is parallel to the imaginary axis at z₀ and using the rule (15) we find that

Hence we see that the density of zeros at the transition point z₀ is nonzero at a first-order phase transition (i.e.one at which the first derivative of the free energy is discontinuous).

Second-order transition

If b₁ = b₂ but c₂¹ c₁ we have that the curve C obeys the equation = ± . Since the zeros of the partition function come in complex-conjugate pairs, we find that the zeros approach the point z₀ along straight lines that meet at a right-angle. If we consider the line = = s/, we find from (15) that

This reveals that at a second-order phase transition, the density of partition function zeros decreases linearly to zero at the phase transition point.

Higher-order transition

More generally, one can consider a difference in the free energies either side of the transition point to have the leading behaviour f₂(u) - f₁(u) ~ u^a. Then the condition that Re[f₂(u)- f₁(u)] = 0 suggests that the curve of zeros approaches the real axis at an angle . Unless a = 1, the curve does not pass smoothly through the real axis. In any case, the imaginary part of the free energy difference grows as |u|^a giving rise to a density of zeros that behaves as s^a-1 for small arclength s. This means that for the density of zeros to be finite at the transition point we must have a ³ 1.

In this section we have summarised what we refer to as the Lee-Yang theory of phase transitions. This describes how partition function zeros are related to phase transitions: the accumulation of zeros at a point along the physical (real, positive) axis of the control parameter gives the critical value and the density of zeros near to the accumulation point determines the order of the phase transition.

We derived a rule (10) for locating phase boundaries and a further rule (15) for finding the density of zeros along such boundaries. At a first-order transition there is a nonzero density of zeros at the transition point whereas the density decays as a power-law to zero at the transition point when the associated phase transition is continuous.

Although we have discussed the theory of partition function zeros with reference to the specific system described by the partition function (1) we should note that the ideas hold much more generally. Firstly, one is not restricted to the canonical ensemble: indeed in the original exposition of the theory [2], Lee and Yang worked in the grand canonical ensemble in which the quantity generalised to the complex plane was not a function of the temperature but rather the chemical potential. Of course, in the equilibrium theory, these intensive field-like quantities play similar mathematical roles and so there is no reason why the Lee-Yang theory shouldn't apply to all of them. For historical reasons, zeros in the complex fugacity (or chemical potential) plane are often referred to as Lee-Yang zeros [2, 3] and those in the complex temperature plane Fisher zeros [4]. Also it appears that one can just as easily generalise physical 'field-like' variables (such as temperature) or 'fugacity-like' variables (such as the quantity z considered above) to the complex plane without altering the properties of the partition-function zero densities at first-order and continuous phase transitions described above.

Despite the apparently wide generality of the Lee-Yang theory of equilibrium phase transitions, proving its validity in the general case is a difficult task (although we note recent work in this direction [10]). Therefore, most rigorous results, such as those discussed in [8], tend to rely on specific properties of a particular partition function. Perhaps the most spectacular of these is that obtained by Lee and Yang in their original work. Specifically, they found that the zeros of the partition function for an Ising ferromagnet in an external field h all lie on the circle |exp(h)| = 1. The significance of this result, which does not depend on the number of spatial dimensions, lattice structure and details of the spin-spin interactions, is that if there is a phase transition induced by varying the magnetic field h, it can only occur at the point h = 0, and then only if the partition function zeros accumulate there in the thermodynamic limit.

III Nonequilibrium steady states

A nonequilibrium steady state differs from its equilibrium counterpart in that it may admit a flow of, say, energy or mass. More generally, these states have a circulation of probability in the space of microscopic configurations. In the past few years, it has become customary to model nonequilibrium systems as stochastic processes, i.e. simple models defined by local dynamical rules. Extensive study of these has revealed a range of phenomena including nonequilibrium phase transitions [11-17]. We present below the key elements of this approach to nonequilibrium physics with a view to understanding nonequilibrium phase transitions in the framework of the Lee-Yang theory described above. In particular, we will need to propose a quantity to use in place of the equilibrium partition function (1).

Let us begin by discussing how one might realise the dynamics of the equilibrium spin system of the previous section, for example in a computer simulation. The aim is to generate a sequence of spin configurations such that the frequency with which a particular configuration is generated is proportional to the Boltzmann weight f() = exp[-bE()] where E() is the energy of configuration . Usually this is achieved by choosing the next configuration ¢ from the current configuration with a probability proportional to the transition rate W(® ¢) that satisfies the detailed balance condition

Apart from the fact that this choice of transition rates guarantees convergence to the desired equilibrium distribution of microstates [18], it also has the feature that the mean current of any quantity in the steady state is zero (as it should be at equilibrium).

Since we are interested in nonequilibrium systems that support steady-state currents we must work with transition rates that do not satisfy the condition (20) or equivalent criteria (given in e.g. [19, 14, 17]). Clearly one has a lot of freedom in the choice of transition rates, and so in practice one often devises dynamics that seem physically reasonable given the phenomena one is trying to describe. Later, in section we will give concrete examples in the form of driven diffusive systems, in which particle moves are biased in a particular direction to model the effect of an external field.

In order to find the set of steady-state weights associated with a particular set of prescribed transition rates, we impose the condition that the total flow of probability into a configuration is balanced by the corresponding outflow. That is, we must have

for every configuration . In general there might be more than one solution to this set of equations for the weights f(), each corresponding to a steady state that is reached with some probability that depends on the starting configuration. We shall assume for simplicity that the steady state is unique and therefore reached with certainty from every initial state. A sufficient condition for a unique steady state is that it is possible to reach each configuration from every other via a sequence of microscopic transitions [17].

Once the steady state weights f() are known from (21) one can obtain the corresponding probability distribution of microscopic configurations P() through a normalisation

such that P() = f()/Z. Then, one can compute averages of physical quantities and look for nonanalyticities to locate nonequilibrium phase transitions as one varies the transition rates. At equilibrium we saw that the zeros of the function (1) that normalises the Boltzmann weights encodes the phase behaviour of the system. Thus we might hope that more generally, the zeros of the normalisation (22) will provide information about nonequilibrium phase behaviour.

In section IV we review recent work that suggests that the steady-state phase behaviour of certain nonequilibrium driven diffusive and reaction-diffusion systems is correctly described by the zeros of the normalisation Z. Although we are unaware of any rigorous argument for this to be the case, some suggestive evidence is provided by an observation made in [20, 17] which we now outline.

First note that the equation (21) is linear in the weights f(). This implies that these weights, and hence the normalisation, can always be written as sums of products of the transition rates W(® ¢). Now, by numbering the microscopic configurations, one can construct the transition rate matrix whose off-diagonal elements _nm are equal to the transition rates from configuration m to n and the diagonal elements _nn are negative and express the total rate of departure from configuration n. Using elementary results from matrix theory (see [20, 17] for the details), one can show that an expression for the normalisation Z that is polynomial in the transition rates (®¢) is

in which

_i are the eigenvalues of the matrix .

With each of the eigenvalues

_i is associated an eigenvector of describing a 'mode' of the stochastic process that decays exponentially with a timescale _i = 1/|_i| [17]. We have assumed that the process described by has only one steady state, and so only one of the eigenvalues is equal to zero (since clearly the relaxation time of a steady state is infinite). Equation (23) states that the normalisation Z can be written as a product of the remaining eigenvalues. Now, at a phase transition we expect diverging timescales: physically one encounters metastable (long-lived) states near first-order transition points [21] and long correlation lengths and times at continuous transitions. The presence of long timescales _i implies small eigenvalues _i of and, from (23), that the normalisation Z approaches zero. Hence it appears that it is appropriate to consider the zeros of Z as given by (22) to locate nonequilibrium phase transitions.

IV Application of Lee-Yang theory to nonequilibrium phase transitions

We shall now review recent progress in applying the Lee-Yang theory to nonequilibrium phase transitions. We consider first in section IV.1 driven diffusive systems, where most work has so far been focussed [21-23]. An appealing feature, discussed in more detail below, is that some one-dimensional cases have been solved exactly, the normalisation (23) calculated, and nonequilibrium phase transitions analysed. (The existence of one-dimensional phase transitions contrasts with the case of one-dimensional equilibrium models which do not admit phase transitions if the interactions are short-ranged.) Then in section IV.2 we move onto reaction-diffusion systems and directed percolation. We do not endeavour to describe all possible nonequilibrium dynamics in the following - for example, self-organised criticality has also been studied using a Lee-Yang approach [25].

IV.1 Driven diffusive systems

In their original papers on partition function zeros, Lee and Yang [2, 3] made use of the mapping between the Ising model of a magnet and the lattice gas. Essentially one associates up-spins with particles and down-spins with vacancies so that a positive interaction strength J in the Ising model results in a particle-particle attraction in the lattice gas, and a negative interaction strength gives rise to repulsive gas particles. Note that there is implicitly at most one particle per lattice site, so there is a hard-core exclusion in the lattice gas.

As discussed in section III, one can realise the dynamics of the lattice gas through a set of transition rates that satisfy the detailed balance condition (20) with respect to the Boltzmann distribution. Thirty years after Lee and Yang's work, Katz, Lebowitz and Spohn (KLS) [20, 27] introduced a driven lattice gas model in which the rate at which particles hop in the direction of an external field is enhanced and the hop-rate against the field is suppressed. This model is well-studied and many results are discussed in [28]. The principal effect of the driving field is to introduce anisotropy in the phase-separation that occurs below the critical temperature associated with the spontaneous magnetisation of the Ising ferromagnet (J > 0) in two or more spatial dimensions and in the critical exponents that characterise this transition.

As yet, the KLS model remains unsolved for general interaction strength, although in one dimension the steady state is known for some parameters [29]. The particular case of one dimension and zero interaction strength J = 0, known as the asymmetric simple exclusion process (ASEP), had already been studied-at least at a mean-field level-by biophysicists interested in the kinetics of biopolymerisation [30]. The mean-field approach predicts phase transitions in the steady state as parameters controlling the rate of insertion and extraction of particles at the boundaries are varied [31]. The existence of these phase transitions is confirmed through an exact solution of the ASEP [29-31], achieved using a powerful matrix product approach [32, 34] which has subsequently been used to solve many generalisations of the ASEP. The details of the matrix product method are not necessary for the following-suffice to say that one ends up calculating a normalisation proportional to (23) through a product of matrices, often of infinite dimension. In this way one obtains some explicit formulas for the normalisation (22) which we shall use below ¹.

The asymmetric exclusion process with open boundaries is perhaps the simplest exactly solved nonequilibrium model that exhibits both a first-order and continuous phase transition in its steady state. Therefore it is an ideal candidate for testing the hypothesis outlined in section 3 that zeros of the normalisation should accumulate towards the positive real axis in the complex plane of transition rates. Before outlining the results of this analysis (the details of which are presented in [23]) we recall the definition of the ASEP with open boundaries.

In this system, a particle on an N-site lattice can hop one site to the right at unit rate, as long as the receiving site is empty. Meanwhile, particles are inserted onto the leftmost lattice site (if empty) at a rate a and removed from the rightmost site (if occupied) at a rate b - see Fig. 2. Along the line a = b < there is a coexistence between a high- and a low-density phase, with a shock separating the two. This is indicative of a first-order transition. Meanwhile, along the lines a > , b = and b > , a = there is a continuous transition (i.e. one accompanied by diverging lengthscales) to a phase in which the particle current is independent of a and b. The phase diagram for the model is shown in Fig. 3.

The steady-state normalisation (22) for the ASEP with N sites has been calculated [32] as

It is a simple matter to use a computer algebra package to solve this equation for its zeros in the complex-a plane at fixed N and b. (Equivalently, one could look at the complex-b zeros at fixed a since Z_N(a,b) = Z_N(b,a).) In Fig.4 plots of the zeros are shown for b = 1, where a continuous transition occurs at a = , and for b = , where a first order transition occurs at a = . We immediately notice that the curve of zeros seems to intersect the real positive a axis at the correct transition point. Furthermore, the density of zeros near the first-order transition point (a = b = ) seems to be uniform and nonzero, whereas the density of zeros near the second-order transition point (a = , b = 1) seems to decrease to zero. Both of these observations are in accord with the results known for equilibrium partition function zeros discussed in section II.

As shown in [23], the distribution of zeros in the thermodynamic limit can be calculated once one knows that, for large N, the normalisation behaves as Z_N ~ A J^-NN^g. In this expression, A, J and g depend on a and b and the quantity J is the current of particles across the lattice. Thus, the complex free energy (3) in the limit N ® ¥ is

Furthermore for a smaller than the transition value a_c, J = a(1 - a) whereas for a > a_c, J = a_c(1-a_c).

Although an electostatic analogy was used in [23] to find the zero distribution, the mathematical content is the same as that used to derive the two rules (10) and (15) in section II. Applying the first rule, which demands that the real part of the free energy be continuous across a phase boundary, we find that the zeros of Z_N(a,b) should lie along the curve |a(1-a)| = a_c(1-a_c). It is therefore convenient to change variable to x = a(1-a). Then, in the complex-x plane, the zeros lie on the circle |x| = x_c = a_c(1-a_c). The density of zeros on this circle can be found by setting x = re^iq and parametrising points on the circle as s = x_cq. Then, the second rule (15) gives the density of zeros m(s) on the circle as

That is, in the complex-x plane the zeros should become evenly distributed on a circle in the thermodynamic limit. Transforming the zeros of (24) obtained at different system sizes to the complex-x plane reveals this to be the case [23].

Finally, one can show that near the intersection point between the curve of zeros in the complex-a plane and the positive real a axis, the zeros of (24) sit on the curve x = - (y² + - x_c)^1/2 where x and y are the real and imaginary parts of a. For the case b < ¹/₂, x_c < ¹/₄ and the transition is first order. One finds that the curve of zeros is smooth at the transition point a = b, and that the density of zeros is (1-2b)/[2pb(1-b)] which is nonzero. These are precisely the properties of the equilibrium partition function zeros at a first-order transition point (see section II).

At the continuous transition point (b ³ and x = ), the zeros pass through the real a axis along the line x = - |y|. Recall from section II that a transition of n^thorder has the curve of zeros meeting the positive real axis at an angle of . Here the zeros clearly approach at an angle suggesting a second-order transition. In fact, one finds the density of zeros is s at a distance s along this line, confirming that the transition is second-order.

In summary, we have found that the Lee-Yang theory of first-order and continuous phase transitions applies to the normalisation of the nonequilibrium asymmetric exclusion process just as it does to the partition function of equilibrium systems. Of course, this does not prove that the theory is generally applicable, and so there is some value in investigating other nonequilibrium steady states that exhibit phase transitions.

One such state is that initially studied by Arndt, Heinzel and Rittenberg [35, 36]. This model comprises M₊ (M_-) positively (negatively) 'charged' particles on a closed ring of N sites. When next to vacant sites, the positive particles hop to the right and negative particles to the left at unit rate, with hops in the opposite directions disallowed. Meanwhile, should two oppositely charged particles be next to one another, the following transitions can occur where the label above the arrow indicates the rate at which the transitions occur. These dynamics are illustrated in Fig.5.

Although a matrix product solution for this model is known to exist [36], it is technically difficult to work out the solution in its full generality [37]. However in the case of either a single vacancy [38] or a single negative particle, i.e.M_- = 1 [39] the steady-state normalisation is known explcitly and phase transitions have been identified. In the case M_- = 1 the zeros of the partition function in the complex-q plane have been studied [24] and, once again the continuous phase transition present in this model is correctly described by the Lee-Yang theory.

A slightly different Lee-Yang analysis, that precedes the work described above, has been used to study the nonequilbrium phase transition that occurs at a finite density of vacancies [22]. The motivation behind this study arose from computer simulations of the model [35, 36] for the case where the numbers of positive and negative particles were equal. In these simulations, three phases were initially identified: a 'solid' phase for small q, a 'fluid' phase for large q and a mixed phase comprising the background 'fluid' and a large mobile droplet in an intermediate regime 1 < q < q_cwhere q_c is some density-dependent quantity.

The steady-state normalisation Z_N,M for the case of M₊ = M_- = M positive and negative particles on the N-site ring is not known exactly (at least in the canonical ensemble), therefore Arndt [22] chose to study the generating function

in the complex-z plane. Physically this approach is equivalent to placing the ring with nonequilibrium interactions in chemical equilibrium with a particle reservoir at fixed fugacity z. One can check that as one takes the thermodynamic limit, the relative size of the fluctuations in the density r = M/N vanish, and so working at fixed fugacity z is equivalent to working at fixed density r. (In fact, this is a standard 'trick' for dealing with closed systems in which the particle number is conserved, see e.g.[40].)

The key point here, however, is that z is an equilibrium fugacity, and not a microscopic transition rate, and so the Lee-Yang theory of phase transitions described in section II ought to apply directly here, without reference to the discussion of section III. It was found [22] that when q < q_c the zeros of (27) in the complex-z plane appear to lie along ellipses, intersecting the positive real axis. Conversely, for larger q the zeros appear to describe hyperbolæ that avoid the positive real axis. After investigating numerically the density of zeros in the complex-z plane, Arndt concluded that there is a fluid-fluid phase transition at some density-dependent point q_c > 1.

Unfortunately, an exact asymptotic (i.e.large N) analysis of the generating function (27) shows that the only nonanalyticity occurs at the solid-fluid transition point q = 1 [37]. However it was noted that physical quantities vary rapidly, but continuously, at some point q > 1. This phenomenon has been explained as an abrupt increase in a correlation length to an anomalously large, but finite value [41]. We believe that resolution with the study of Arndt lies in the fact that the system sizes considered in [22] were quite small (N < 100), whereas it was suggested [37] that the distinction between this crossover behaviour and a genuine phase transition might become apparant only once N is increased above about 10⁷⁰ (an unfeasibly huge number!).

It would be interesting, therefore, to extend the numerical computation of zeros performed in [22] to much larger systems, to see how the ellipses noted above develop. For example it might well happen that instead of approaching the positive-real fugacity axis, the zeros of (27) would terminate a short distance away from it.

IV.2 Reaction-diffusion systems and directed percolation

A large and important class of models with stochastic dynamics is provided by reaction-diffusion systems². In contrast to driven diffusive systems, where the particle-particle interactions imply conservation of particles, reaction-diffusion systems are characterised by dynamics that result in a change in particle number. Moreover there are a number of such systems that have absorbing states i.e.special configurations generally devoid of particles that once entered cannot be left. Phase transitions associated with whether the system has a finite probability of not being absorbed into such a state fall within the directed percolation and related universality classes [13]. We shall shortly discuss the second order phase transition associated with the directed percolation university class in a little detail. Meanwhile as a simple example of a reaction-diffusion system, we review a model for which the steady state that can be solved using the matrix product approach. The approach again provides us with an explicit expression for the normalisation (22) by virtue of which we can analyse its zeros in the plane of complex reaction rates.

The system in question [42] has for the dynamics at neighbouring bulk sites on a one-dimensional lattice the processes

occurring at the rates indicated, and with · representing a particle and ° an empty lattice site. It was demonstrated [42] that the matrix product scheme used for the ASEP could be generalised to cater for the steady state of the present reaction-diffusion system on a lattice of N sites with reflecting boundary conditions (i.e.which is neither periodic nor has particle input or removal at the boundaries).

In the matrix product approach, the normalisation Z_N is given by a scalar derived from C^N where C is a square matrix. If C is a finite dimensional m × m matrix and can be diagonalised, the resulting expression for Z_N has the form

in which

_n is the n^th eigenvalue of C and the coefficients a_n arise from the details of the way in which one obtains a scalar from the matrix C^N. A normalisation of this form leads to a complex free energy

in which the maximum means the eigenvalue with the largest absolute value. Then, there is the possibility of a phase boundary when the magnitude of the two largest eigenvalues of C are equal. It turns out that for the closed reaction-diffusion system introduced above, C is a 4 ×4 matrix that can be diagonalised when q²¹ 1 + k [42]. When q² = 1 +k, C cannot be diagonalised on account of the largest eigenvalue being degenerate, and a phase transition is found to occur at this point. Note that this scenario contrasts with the transfer-matrix approach to one-dimensional equilibrium systems where the partition function is also written as a product of matrices. Since all elements of the transfer matrix are positive the largest eigenvalue cannot become degenerate and therefore there can be no phase transition [12]. However there is no such restriction on the elements of C, and so eigenvalue crossing is permitted and nonequilibrium one-dimensional phase transitions can occur.

The structure of the normalisation Z_N(q,k) is even simpler in the case where particles are inserted onto and removed from the left boundary at rates a and b respectively (with the right boundary remaining reflecting), and these rates satisfy the relation a = k(q^-1 - q + b) [43]. In this case C is a 2×2 matrix and, for q²¹ 1+k, is diagonalisable with eigenvalues

We see then immediately from (29) and an application of the rule (10) that there is a phase boundary in the complex q plane along the circle |q| = . Also, since one of the eigenvalues does not depend on q, the free energy h(q,k) is a constant in one of the phases which, as the analysis of the normalisation zeros for the ASEP demonstrated, implies that the density of zeros on this circle is constant in the thermodynamic limit. In turn, this implies that the phase transition is first order, as confirmed by explicitly calculating the density profile in the two phases [43].

We finally turn our attention to a process that comprises symmetric decoagulation (that is,

® and ® taking place with equal probability in each direction) and spontaneous decay ( ® ) occurring at independent rates. Introduced as a crude model of an epidemic [44], this contact process is known to exhibit a transition from a phase in which the absorbing state (empty lattice) is reached with certainty to a phase in which there is some probability that the epidemic remains active forever in the thermodynamic (infinite system size) limit [45]. This transition occurs as the ratio between the decoagulation and decay rates is increased beyond a critical value.

Although not proven, it is widely accepted on the basis of simulation and approximate methods (such as series expansions) that the transition just described is continuous and characterised by directed percolation (DP) exponents (see [13] for an in-depth discussion of these issues). Directed percolation itself [46] is a geometric construction designed to model fluid flow through a random porous medium under the influence of gravity. Consider the rhombic lattice shown in Fig.6 and imagine 'pouring' fluid in from the uppermost site. The fluid is allowed to flow along a bond pointing diagonally downwards and then only if it is open. Each bond has a probability p of being open (and hence a probability 1-p of being blocked), and the state of each bond is independent of any other and may not change with time.

The main quantity of interest in this system is the percolation probability P_n(p) that the fluid can penetrate to a depth n in the medium (see Fig.6), averaged over all possible bond configurations. The order parameter for the model is the probability P_¥(p) that the fluid can penetrate infinitely far given a bond probability p. If p is below some critical percolation threshold p_c, the fluid only ever penetrates a finite distance-i.e.a layer becomes dry with certainty. However, for larger p, the order parameter becomes nonzero, growing as P_¥(p) ~ |p-p_c|^b with b » 0.276 near the critical point. The importance of the DP transition is that a wide range of models that have a transition into an absorbing state are expected to have that transition characterised by the DP exponents, one of which is b [13]. Despite a huge amount of interest in directed percolation, none of the models expected to belong to its universality class has been solved exactly.

Recently an attempt has been made to shed further light on the DP transition by studying the zeros of the percolation probability P_n(p) [47, 48]. At first glance, a connection between this probability and a partition function is not obvious. However, associating with each site a 'spin' state s_i, that is up if site i is connected to a point on the n^th layer and down otherwise, yields a form for P_n(p) that has the structure of a transfer-matrix representation of a partition function for an equilibrium system with three-spin interactions [49].

It is not clear that being able to give P_n(p) the appearance of a partition function necessarily implies that the Lee-Yang theory should hold. In particular, P_n(p) has some features that set it apart from 'standard' equilibrium partition functions, in that it does not grow exponentially in the number of bonds N, at least in the range 0 < p < 1. It is also uncertain whether a free energy can meaningfully be defined for this system, especially in the region 0 < p < p_c in which lim_{n® ¥}P_n(p) = 0.

Given these observations, it is perhaps not surprising that the zeros of P_n(p) calculated numerically for n < 15 (leading to polynomials of degree up to N = 240) have a more complicated distribution than for the (exactly solved) models considered so far. As is evident from Fig.7, which was presented in [48], the zeros do not lie along a single curve but along a sequence of curves meeting at the critical point (as well as inside some region that encloses part of the real axis for p > 2).

By considering the sequences of zeros approaching the positive real axis, one can estimate the critical point p_c and the density of zeros as it is reached [48]. The latter yields a prediction for one of the DP critical exponents and one finds good agreement with the most precise estimates of both the transition point and the associated exponent obtained through other means. Hence we have further evidence for the applicability of Lee and Yang's ideas concerning partition function zeros to a much wider range of statistical distributions than equilibrium steady states.

V Summary and outlook

In this work we have revisited the Lee-Yang description of equilbrium phase transitions with a view to seeing whether the ideas apply to more general nonequilibrium transitions. Recently there have been a number of studies of zeros of partition-function-like quantities that arise in systems with nonequilibrium dynamics, and we have seen in our review of these works that the Lee-Yang theory, as described in sectionII, seems to hold quite generally.

We have argued that for dynamic models with a unique steady state, the normalisation defined as a sum over the steady-state configurational weights (22) serves as a suitable 'partition function' in the sense that its zeros, in the complex plane of any model parameter, should accumulate towards physical transition points in the thermodynamic limit. Furthermore, the density of zeros and angle of approach to the real axis indicate whether the transition is first-order (manifested physically through phase-coexistence) or continuous (i.e.characterised by divergent correlation lengths and times). Thus studying the zeros of the normalisation (22) provides an unambiguous classification of nonequilibrium phase transitions as do the Lee-Yang zeros in the equilibrium case.

The observation that backs up this scenario is embodied by equation (23) which reveals that the reciprocal of the steady-state normalisation is equal to the product of the characteristic relaxation times in the dynamics. Since near a phase transition one expects timescales to diverge, one also expects the normalisation to approach zero. However, a rigorous argument for this to be the case is still lacking. Moreover (23) implies a possible link between systems for which the steady state normalisation can be calculated and those for which eigenvalues of the transition matrix can in principle be calculated.

A different class of systems encompasses those whose steady state is not unique. The contact process is, in fact, an example of such a model, in which the absorbing state is reached with certainty below the critical decoagulation rate, whereas above it, and on an infinite system, a second steady state can also be reached with some nonzero probability. Since this additional steady state exists only when the lattice size becomes infinite, one must take that limit first, before taking time to infinity. Otherwise, on a finite system the steady state is simply the absorbing state and the steady-state normalisation is trivially equal to a constant. Nevertheless the work of [47, 48], which we reviewed in section IV.2, indicates that the Lee-Yang theory can be relevant when one considers other properties of the nonequilibrium system such as the percolation probability.

Although we have not discussed this in great detail here, it should be noted that the Lee-Yang approach gives a method for extrapolating to the thermodynamic limit from solutions for small system sizes. In the work of [48], numerical solutions for small systems were used successfully to estimate the transition point and density of zeros as it is approached. From this information one learns about the nature of the phase transition and, for example, can estimate the values of critical exponents. This is a common technique in equilibrium statistical physics (see e.g.[50]) and it may be the case that stochastic processes unyielding to analytical treatment could be understood this way.

Acknowledgments

We thank Bernard Derrida for helpful discussions. R.A.B. acknowledges financial support under EPSRC grant GR/R53197.

Received on 4 April, 2003

• [1] C.Domb, The critical point: a historical introduction to the modern theory of critical phenomena (Taylor & Francis, London, 1996).
• [2] C.N. Yang and T.D. Lee, Phys. Rev. 87, 404 (1952).
• [3] T.D. Lee and C.N. Yang, Phys. Rev. 87, 410 (1952).
• [4] M.E. Fisher, The nature of critical points, in Lectures in Theoretical Physics, edited by W.E. Brittin, volume7C, p.1 (University of Colorado Press, Boulder, 1965).
• [5] S.Grossmann and W.Rosenhauer, Z. Physik 218, 437 (1969).
• [6] S.Grossmann and V.Lehmann, Z. Physik 218, 449 (1969).
• [7] L.E. Reichl, A Modern Course in Statistical Physics (University of Texas Press, Austin, 1980).
• [8] D.Ruelle, Statistical Mechanics: Rigorous Results (W A Benjamin, New York, 1969).
• [9] O.D. Kellogg, Foundations of Potential Theory (Dover, New York, 1953).
• [10] M.Biskup, C.Borgs, J.T. Chayes, L.J. Kleinwaks, and R.Kotecký, Phys. Rev. Lett. 84, 4794 (2000).
• [11] J.Marro and R.Dickman, Nonequilibrium Phase Transitions in Lattice Models (Cambridge University Press, Cambridge, 1999).
• [12] M.R. Evans, Braz. J. Phys. 30, 42 (2000).
• [13] H.Hinrichsen, Adv. Phys. 49, 815 (2000).
• [14] D.Mukamel, Phase Transitions in Nonequilibrium Systems, in Soft and Fragile Matter, edited by M.E. Cates and M.R. Evans, Vol.53 of Scottish Universities Summer School in Physics, p. 237 (Institute of Physics Publishing, Bristol, 2000).
• [15] G.M. Schütz, Exactly Solvable Models for Many-Body Systems Far From Equilibrium, in Phase Transitions and Critical Phenomena, edited by C.Domb and J.L. Lebowitz, Vol.19, p.1 (Academic Press, London, 2001).
• [16] R.Stinchcombe, Adv. Phys. 50, 431 (2001).
• [17] M.R. Evans and R.A. Blythe, Physica A 313, 110 (2002).
• [18] O.Narayan and A.P. Young, Phys. Rev. E 64, 021104 (2001).
• [19] N.G. van Kampen, Stochastic processes in physics and chemistry (North-Holland, Amsterdam, 1992) second edition.
• [20] R.A. Blythe, Nonequilibrium Phase Transitions and Dynamical Scaling Regimes Ph.D. thesis, University of Edinburgh (2001).
• [21] B. Gaveau and L. S. Schulman, J. Math. Phys. 39, 1517 (1998).
• [22] P.F. Arndt, Phys. Rev. Lett. 84, 814 (2000).
• [23] R.A. Blythe and M.R. Evans, Phys. Rev. Lett. 89, 080601 (2002).
• [24] F.H. Jafarpour, J. Phys. A: Math. Gen. 36, 7497 (2003).
• [25] B. Cessac and J. L. Meunier, Phys. Rev. E 65 036131 (2002).
• [26] S.Katz, J.L. Lebowitz, and H.Spohn, Phys. Rev. B 28, 1655 (1983).
• [27] S.Katz, J.L. Lebowitz, and H.Spohn, J. Stat. Phys 34, 497 (1984).
• [28] B.Schmittmann and R.K.P. Zia, Statistical Mechanics of Driven Diffusive Systems, Vol.17 of Phase Transitions and Critical Phenomena (Academic Press, London, 1995).
• [29] J.Krug and H.Spohn, Kinetic Roughening of Growing Surfaces, in Solids Far From Equilibrium, edited by C.Godrèche (Cambridge University Press, Cambridge, 1992).
• [30] C.T. MacDonald, J.H. Gibbs, and A.C. Pipkin, Biopolymers 6, 1 (1968).
• [31] B.Derrida, E.Domany, and D.Mukamel, J. Stat. Phys. 69, 667 (1992).
• [32] B.Derrida, M.R. Evans, V.Hakim, and V.Pasquier, J. Phys. A: Math. Gen. 26, 1493 (1993).
• [33] G.Schütz and E.Domany, J. Stat. Phys. 72, 277 (1993).
• [34] B.Derrida, Phys. Rep. 301, 65 (1998).
• [35] P.F. Arndt, T.Heinzel, and V.Rittenberg, J. Phys. A: Math. Gen. 31, L45 (1998).
• [36] P.F. Arndt, T.Heinzel, and V.Rittenberg, J. Stat. Phys. 97, 1 (1999).
• [37] N.Rajewsky, T.Sasamoto, and E.R. Speer, Physica A 279, 123 (2000).
• [38] T.Sasamoto, Phys. Rev. E 61, 4980 (2000).
• [39] F.H. Jafarpour, J. Phys. A: Math. Gen. 33, 8673 (2000).
• [40] B.Derrida, S.A. Janowsky, J.L. Lebowitz, and E.R. Speer, J. Stat. Phys. 73, 813 (1993).
• [41] Y.Kafri, E.Levine, D.Mukamel, and J.Torok, J. Phys. A: Math. Gen. 35, L459 (2002).
• [42] H.Hinrichsen, S.Sandow, and I.Peschel, J. Phys. A: Math. Gen. 29, 2643 (1996).
• [43] F.H. Jafarpour, J. Phys. A: Math. Gen. 36, 7497 (2003).
• [44] T.E. Harris, Ann. Prob. 2, 969 (1974).
• [45] D.Griffeath, Additive and Cancellative Interacting Particle Systems, Vol. 724 of Lecture Notes in Mathematics (Springer-Verlag, Berlin, 1979).
• [46] S.R. Broadbent and J.M. Hammersley, Proc. Camb. Phil. Soc. 53, 629 (1957).
• [47] P.F. Arndt, S.R. Dahmen, and H.Hinrichsen, Physica A 295, 128 (2001).
• [48] S.M. Dammer, S.R. Dahmen, and H.Hinrichsen J. Phys. A: Math. Gen. 35, 4527 (2002).
• [49] R.J. Baxter and A.J. Guttmann, J. Phys. A: Math. Gen. 21, 3193 (1988).
• [50] W.Janke and R.Kenna, Nucl. Phys. B Proc. Supp. 106, 905 (2002). See also hep-lat/0112032.
• [51] S.Redner, Scaling Theories of Diffusion-Controlled and Ballistically-Controlled Bimolecular Reactions, in Nonequilibrium Statistical Mechanics in One Dimension, edited by V.Privman, (Cambridge University Press, Cambridge, 1997), chapter 1.

1

As yet, there is no comprehensive review of applications of the matrix-product method, although one should consult [34, 13, 15, 16] and references therein.

2

As an area of ongoing research, new results are emerging continuously and there is no up-to-date review that contains them all. Nevertheless, one can consult [51, 13] for an introduction.

Publication Dates

• Publication in this collection
  11 Nov 2003
• Date of issue
  Sept 2003

History

• Accepted
  04 Apr 2003
• Received
  04 Apr 2003