Dynamical Transition in the Open-boundary Totally Asymmetric Exclusion Process
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] N ov Dynamical Transition in the Open-boundary TotallyAsymmetric Exclusion Process
A. Proeme, R. A. Blythe and M. R. Evans
SUPA, School of Physics and Astronomy, University of Edinburgh, Mayfield Road,Edinburgh EH9 3JZ, United KingdomE-mail: [email protected] , [email protected] , [email protected] Abstract.
We revisit the totally asymmetric simple exclusion process with openboundaries (TASEP), focussing on the recent discovery by de Gier and Essler that themodel has a dynamical transition along a nontrivial line in the phase diagram. Thisline coincides neither with any change in the steady-state properties of the TASEP,nor the corresponding line predicted by domain wall theory. We provide numericalevidence that the TASEP indeed has a dynamical transition along the de Gier–Esslerline, finding that the most convincing evidence was obtained from Density MatrixRenormalisation Group (DMRG) calculations. By contrast, we find that the dynamicaltransition is rather hard to see in direct Monte Carlo simulations of the TASEP. Wefurthermore discuss in general terms scenarios that admit a distinction between staticand dynamic phase behaviour.Date: October 26, 2010 ynamical Transition in the TASEP
1. Introduction
As a rare example of a class of exactly-solvable, nonequilibrium interacting particlesystems, asymmetric simple exclusion processes of various types have found favour withresearchers in statistical mechanics, mathematical physics and probability theory [1–4].These systems comprise hard-core particles hopping in a preferred direction on a one-dimensional lattice and have been used to describe systems as diverse as traffic flow [5],the dynamics of ribosomes [6] and molecular motors [7, 8], the flow of hydrocarbonsthrough a zeolite pore [9], the growth of a fungal filament [10], and in queueingtheory [11, 12]. The physical interest lies mainly in the rich phase behaviour that arisesas a direct consequence of being driven away from equilibrium.Of particular interest are the open boundary cases. Here the system can be thoughtof as being placed between two boundary reservoirs, generally of different densities. Thetwo reservoirs enforce a steady current across the lattice, and therewith a nonequilibriumsteady state. It was first argued by Krug [13] that as the boundary densities are variednonequilibrium phase transitions may occur in which steady state bulk quantities—such as the mean current or density—exhibit nonanalyticities. Such phase transitionswere seen explicitly in first a mean-field approximation [6, 14] and then in the exactsolution [15, 16] of the totally asymmetric exclusion process with open boundaries,hereafter abbreviated as TASEP, and of related processes [4]. Since then argumentshave been developed to predict the phase diagram of more general one-dimensionaldriven diffusive systems [17–19].Our interest in this work is in the distinction between two phase diagrams forthe TASEP: one of which characterises the steady-state behaviour, and the other thedynamics. The static phase diagram [15,16] is shown in Fig. 1(a) where the left boundarydensity is α and the right is 1 − β . There are three possible phases in the steady state:a low density (LD) phase where the bulk density is controlled by the left boundary andis equal to α ; a high density (HD) phase where the bulk density is controlled by theright boundary and is equal to 1 − β ; a maximal current (MC) phase where the bulkdensity is . The high and low density phases are further divided into subphases (e.g.LD1, LDII) according to the functional form of the spatial decay of the density profileto the bulk value near the non-controlling boundary. In these subphases the lengthscaleover which the decay is observable remains finite as the system size is increased. Thischange at the subphase boundaries in the form of the decay is thus not a true phasetransition in the thermodynamic sense.More recently, de Gier and Essler have performed an exact analysis of the ASEP’sdynamics [20–22] and derived the the longest relaxation times of the system bycalculating the gap in the spectrum using the Bethe ansatz. The analysis builds ondirectly related work on the XXZ spin chain with nondiagonal boundary fields [23, 24].The dynamical phase diagram they obtain, in which phases are associated with differentfunctional forms of the longest relaxation time, is illustrated in figure 1(b). The samephases (high density, low density and maximal current) are found as in the static phase ynamical Transition in the TASEP PSfrag replacements j = α (1 − α ) j = α (1 − α ) ρ = αρ = α ρ = 1 − βρ = 1 − β j = β (1 − β ) j = β (1 − β )LDI LDII HDIIHDI MC ρ = j = αβ (a) PSfrag replacements LDI ′ LDII ′ HDII ′ HDI ′ MC αβ (b) Figure 1. (a) Static and (b) dynamic phase diagrams for the TASEP. Solid linesindicate thermodynamic phase transitions at which the current and bulk density arenonanalytic. Dotted lines indicate subphase boundaries. In the static case, the densityprofile a finite distance from the boundary is nonanalytic. In the dynamic case, thelongest relaxation time exhibits a nonanalyticity. diagram, however, a hitherto unexpected dynamical transition line which subdivides thelow density and high density phases is now apparent. This line which we shall refer toas the de Gier–Essler (dGE) line replaces the subphase boundaries at α = 1 / , β < / β = 1 / , α < / ′ , LDII ′ .Although there is nothing to rule out the prediction of dGE of a dynamicaltransition at a distinct location to any static transition, the result came as something ofa surprise. This is because the static phase diagram had been successfully interpreted interms of an effective, dynamical theory thought to be relevant for late-time dynamics,referred to as domain wall theory (DWT). We shall review DWT more fully inSection 2.4. For the moment we note that DWT correctly predicts the static subphaseboundary and the associated change in the decay of the density profile, but does notpredict the dGE line. Therefore, it was previously thought that a dynamical transitionalso occurred at the subphase boundary and initial calculations of relaxation timesappeared to be consistent with this [25, 26].Our primary goal in this work is to establish numerically that a dynamical transitiondoes indeed occur along the dGE line rather than at the static subphase boundary. Weused two distinct approaches. First, we carried out direct Monte Carlo simulations ofthe TASEP dynamics and identified the dominant transient in three time-dependentobservables. We find that these simulations are consistent with a dynamical transitionat the dGE line but are not sufficiently accurate to rule out the scenario of the dynamicaltransition occurring at the subphase boundary. ynamical Transition in the TASEP
2. Model definition and phase diagrams
Although we have alluded to the basic properties of the totally asymmetric simpleexclusion process (TASEP) in the introduction, in the interest of a self-containedpresentation we provide in this section a precise model definition and full details ofthe static and dynamic phase diagrams.
The TASEP describes the biased diffusion of particles on a one-dimensional lattice with L sites. No more than one particle can occupy a given site, and overtaking is prohibited.The stochastic dynamics are expressed in terms of transition rates, for example a particleon a site i in the bulk hops to the right as a Poisson process with rate 1, but only if site i + 1 is unoccupied: At the left boundary only particle influx takes place, with rate α ,and at the right boundary only particle outflux takes place, with rate β , as shown inFig. 2. The corresponding reservoir densities are α at the left and 1 − β at the right. As noted in the introduction, the TASEP static phase diagram can be determinedfrom, for example, exact expressions for the current and density profile in the steadystate [15,16]. The current takes three functional forms according to whether it is limitedby a slow insertion rate (LD: α < β, α < ), by a slow exit rate (HD: β < α, β < )or by the exclusion interaction in the bulk (MC: α > , β > ). This behaviour of ynamical Transition in the TASEP Figure 2.
The dynamics of the totally asymmetric simple exclusion process (TASEP)with open boundaries. Arrows show moves that may take place, and the labels indicatethe corresponding stochastic rates. the current leads to the identification of the thermodynamic phases separated by solidlines in Fig. 1(a). Along these lines there are nonanalyticities in both the current andthe density far from either boundary. The relevant expressions for these quantities areshown on Fig. 1(a).The density near one of the boundaries is nonanalytic along the lines α = and β = in the HD and LD phases respectively. These subphases are shown dotted inFig. 1(a). Inspection of the functional form of the density profile reveals that theseare not true thermodynamic phase transitions, in the sense that the deviation fromthe bulk value extends only a finite distance into the bulk, and thus contributes onlysubextensively to the nonequilibrium analogue of a free energy (see e.g., [31] for thedefinition of such a quantity). To be more explicit, consider for example the behaviournear the right boundary in the low-density phase. Here the bulk density is ρ = α .The mean occupancy of the lattice site positioned a distance j from the right boundaryapproaches in the thermodynamic limit L → ∞ the form [4, 15, 16] ρ L − j ∼ α + c I ( β ) (cid:16) α (1 − α ) β (1 − β ) (cid:17) j α < β < (LDI) α + c II ( α, β ) [4 α (1 − α )] j j / < β (LDII) . (1)In these expressions, c I and c II are functions of the boundary parameters that we leaveunspecified here so as to focus on the lengthscale of the exponential decay from theright boundary. As β is increased from zero, the decay length increases until β = .Then the decay length is constant, and the exponential is modulated by a power law.The density profile at the left boundary exhibits the same kind of nonanalyticity in thehigh-density phase as α is increased through as a consequence of the particle-holesymmetry, ρ i − ( α, β ) = 1 − ρ L − i ( β, α ), exhibited by the model. The dynamic phase diagram is obtained by examining a quantity referred to as the gap by de Gier and Essler [20–22]. This is simply the largest nonzero eigenvalue ofthe transition matrix governing the continuous-time stochastic dynamics of the TASEP.More formally, one starts with the master equationdd t | P ( t ) i = M | P ( t ) i , (2)where the matrix M encodes the transition rates and | P ( t ) i is the vector of probabilitiesfor each configuration. Because M is a stochastic matrix it is guaranteed [32] by the ynamical Transition in the TASEP λ = 0 > Re( λ ) ≥ Re( λ ) ≥ . . . . (3)The spectrum corresponds to a set of relaxation times τ i = − (Re λ i ) − . The longest(non-infinite) relaxation time is τ , and the associated eigenvalue λ is the gap which wewill henceforth denote by the symbol ε . At long times the relaxation of any observableis expected to decay exponentially with a characteristic timescale τ , a fact we will lateruse to estimate the gap from Monte Carlo simulations.The thermodynamic phase boundaries are found by identifying lines along whichthe gap vanishes, indicating the coexistence of two stationary eigenstates (phases). Theexact calculations of de Gier and Essler [20, 21] show that, in the thermodynamic limit L → ∞ , the gap vanishes along the line α = β < that separates the HD and LDphases. The gap also vanishes in the entirety of the MC phase, reflecting the genericlong-range (power-law) correlations seen in this phase [15]. Thus at this level, the staticand dynamic phase diagrams coincide.Where they differ is in the subdivision of the high- and low-density phases ‡ , inwhich the gap remains finite in the limit L → ∞ . There is a region, marked LDI ′ andHDI ′ on Fig. 1(b), within which the gap assumes the asymptotic expression ε ( L ) = − α − β + 2( ab ) + 1 − π ( ab ) − ( ab ) − L − + O ( L − ) (4)in which a = 1 − αα and b = 1 − ββ . (5)Within the low-density phase ( α < β, α < ), this form of the gap applies for values of β < β c where β c ( α ) = " (cid:18) α − α (cid:19) − . (6)Likewise, in the high-density phase ( β < α, β < ), the region within which the gap isgiven by (4) is bounded by α < α c where α c ( β ) = β − β ! − . (7)In the remainder of the low-density phase, α < , β > β c the gap takes the form ε ( L ) = − α − β c + 2( ab c ) + 1 − π ( ab c ) − ( ab c ) − L − + O ( L − ) . (8)Finally, we have by symmetry that when β < , α > α c , ε ( L ) = − α c − β + 2( a c b ) + 1 − π ( a c b ) − ( a c b ) − L − + O ( L − ) . (9) ‡ de Gier and Essler [20–22] refer to these as “massive” phases by analogy with the quantum spinchains; we shall only use the terminology associated with the static phase diagram to avoid confusion. ynamical Transition in the TASEP Figure 3.
Exact L → ∞ gap (black, dashed) as a function of β for α = 0 .
3, describedby the HDI ′ equation (4) (red) for β < β c ( β c ≈ .
57, indicated), and by equation (8)(green) for β > β c . In these expressions, a c = 1 − α c α c and b c = 1 − β c β c . (10)The boundaries between the dynamic subphases—the de Gier–Essler (dGE) lines—areshown dotted in Fig. 1(b). Note that the coefficient of the L − term in this asymptoticexpansion is discontinuous across the dynamical transition line.In order to illustrate the behaviour of the gap along the static and dynamictransition lines, we plot it as a function of β at some α < in Fig. 3. The moststriking feature is the constancy of the gap above the nontrivial critical point β c . Laterin this work, we will use the constancy of the gap above a critical threshold as anempirical means to identify the dynamical transition point. An alternative approach to study the TASEP dynamics is domain wall theory(DWT) [33]. Although much simpler than the Bethe ansatz approach of de Gier andEssler, it correctly predicts the static phase boundaries and the analytical form of thegap in a region of the phase diagram. However, it predicts dynamic subphases thatare distinct from those found by de Gier and Essler, and that correspond to the staticsubphases. In the numerical study that follows, it will be important to be able todistinguish between these two sets of predictions, and so we summarise DWT here.The basis of this approach is to assume that collective relaxational dynamics areeffectively reducible to a single coordinate describing the position of an interface, ordomain wall. The wall separates a domain of density ρ − and current j − to the left froma domain of density ρ + and current j + = ρ + (1 − ρ + ) to the right. Each domain is takento possess the steady state characteristics imposed by the boundary reservoir on that ynamical Transition in the TASEP PSfrag replacements LDI LDII HDIIHDI MC j = αβ ρ + = ρ − = αρ + = 1 − βρ − = αρ + = 1 − βρ − = α ρ + = 1 − βρ − = Figure 4.
The phase diagram obtained from domain wall theory. ρ + and ρ − indicatethe densities of the left and right domains. In the domain wall theory, static anddynamic transitions coincide along the dashed lines. particular side of the wall, with ρ + and ρ − therefore as in Fig. 4. The effective theoryis expected to be exact along the line α = β < / α < / β < /
2. The motion of the wall is thendescribed microscopically as a random walker with left and right hopping rates D − and D + given respectively by imposing mass conservation on the fluxes into and out of thewall [25]: D − = j − ρ + − ρ − = α (1 − α )1 − α − β D + = j + ρ + − ρ − = β (1 − β )1 − α − β , (11)For α > β the random walk is biased to the left and in the stationary state the domainwall is localised at the left boundary and the bulk density is given by ρ + = 1 − β . For α < β the random walk is biased to the right and in the stationary state the domainwall is localised at the right boundary and the bulk density is given by ρ − = α . Thusthe first order transition at α = β is correctly predicted.Moreover, the gap in the resulting spectrum for large L given [25] is given by ε DW T = − D + − D − + 2 √ D + D − − π L + O ( L − ) ! . (12)Remarkably this expression is identical to (4) to order 1 /L . Thus in the region α < / , β < / β > / ρ + = 1 /
2, one can take j + in (11) to depend on the size ℓ of the right-hand domain. That is, one puts j + → j + ( ℓ ) equal to the stationary currentin a TASEP of size ℓ in the maximal current phase. This implies the large ℓ behaviour j + ≃
14 (1 + 32 ℓ ) (13) ynamical Transition in the TASEP Figure 5.
Comparison of the gap for L → ∞ from domain wall theory (dashed line)with the exact gap (solid line) at α = 0 . β c ≈ . This results in a modification of the density profile to an exponential spatial decaymodulated by a power law with power 3/2, similar to Eq. 1 (LDII).In brief, the DWT is remarkably successful, correctly predicting the static phasediagram (including subphases), and the exact thermodynamic gap function found byde Gier and Essler in the region α < / β < /
2. It differs in that the dynamicsubphases are not given by the dGE lines, but the static subphase transition lines. Thus,the region within which the gap is constant is different in the two theories, as indicatedby Fig. 5.
3. Evidence for the dynamical transition from Monte Carlo simulations
We now describe attempts to measure the gap in Monte Carlo simulations of the TASEP,with a particular interest in distinguishing between the dGE and DWT predictions. Asnoted previously, the ensemble average of any time-dependent observable O ( t ) shouldbe dominated at late times by an exponential decay e −| λ | t to its stationary value hOi .Recall that λ is the largest nonzero eigenvalue of the matrix M appearing in the masterequation (2), and that all nonstationary eigenvalues have negative real part. In principle,therefore, all one needs to do is pick an observable, and measure its late-time decay rate.In practice, this is made difficult by the fact that all decay modes may contribute to hO ( t ) i , and that by the time the higher modes have decayed away, the residual signal hO ( t ) i − hOi may be extremely small and swamped by the noise.Since we are interested in the time-dependence of an observable, we employ a continuous-time (Gillespie) algorithm [35] to simulate the TASEP dynamics. Moreprecisely, we maintain a list of events (i.e., a particle hopping to the next site, or enteringor leaving the system) that can take place given the current lattice configuration.A particular event i is then chosen with a probability proportional to its rate ω i asspecified in Section 2. A time variable is then incremented by an amount chosen from ynamical Transition in the TASEP P i ω i ) − . In this way, a member of the ensembleof all continuous-time trajectories of the TASEP dynamics from some prescribed initialcondition is generated with the appropriate probability. We consider first the decay of the total number of particles on the lattice, N ( t ), to itsstationary value as measured by the function R ( t ) = h N ( t ) i − h N ih N (0) i − h N i . (14)In this expression, angle brackets denote an average over an ensemble of initial conditionsand over stochastic trajectories. In each simulation run, an initial condition wasconstructed in which each site was independently occupied with a probability p = 1 − β (in the LD phase) or p = α (in the HD phase). In the bulk, these densities then relaxto the steady-state values displayed on Fig. 1(a).Once R ( t ) has been sampled over multiple trajectories, the task is to identify atime window over which one can fit an exponential decay and therewith estimate agap. The start and end points of this windows are both crucial. If it starts too early,then one may expect contributions from subdominant transients (i.e., the decays at rate λ , λ , . . . ) to systematically skew the estimate of the gap. If it ends too late, noise mayinstead dominate the estimate.The noise at the top end we handle by examining the behaviour of local decay rates µ i = ln R ( t i + i ) − ln R ( t i ) t i +1 − t i (15)where t i and t i +1 are successive time points at which R ( t ) was sampled. At late times,one should find µ i +1 /µ i →
1. Strong deviations from unity indicate the dominance ofnoise, and we rejected points after which the magnitude of this ratio exceeded a criticalvalue. For our data sets, we found that 5 was a suitable choice for this value: seeFig. 6(a) for an example.The bottom end of the window was chosen by maximising a goodness-of-fit measureto a fit of the exponential f ( t ) = a e − λt to data points within the window. We adoptedthe adjusted coefficient of determination [36] R = 1 − n − k P i ( R ( t i ) − f ( t i )) n − P i ( R ( t i ) − ¯ R ) . (16)as our goodness-of-fit measure, varying the start of the window until this quantity wasmaximised. In this expression ¯ R is the arithmetic mean of R ( t ) over the set of n times t i falling within the window, and k is the number of free parameters in the fit function f ( t ), i.e., k = 2. This goodness-of-fit measure trades off the increased quality of thefit obtained by discarding data points against the increasing uncertainty in the fittingparameters that comes with the noisier data at the top end of the window. We showin Fig. 6(b) an example of how the goodness-of-fit varies with the size of the window. ynamical Transition in the TASEP æææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææ - - t Μ i + (cid:144) Μ i (a) æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ æ
80 120 160 2000.999800.999850.999900.999951.00000 t s R (b) Figure 6.
Criteria for choosing the range of data over which to fit a single exponentialto the decay of the density to stationarity. (a) The ratio of gradients µ i +1 /µ i whichshould be close to unity where a single exponential is a good fit. At large times, thisratio becomes noisy; the presence of extreme values (here, a magnitude larger than 5)is used to identify the onset of noise. (b) The behaviour of R as a function of thestart of the window. We choose this such that R is maximised. In both cases, L = 50, α = 0 . β = 0 . We remark that the relative flatness of R above the optimal starting time can beascribed to the fact that only a single decay mode remains in this region. Although thisprocedure is a little cumbersome to apply to each data set, we found it preferable tofitting multiple exponentials to the data for R ( t ), which leads to large uncertainties inthe measured gap.We show results in Fig. 7, in which the gap is plotted as a function of β for fixed α = 0 . L = 10) datafrom Monte Carlo simulation are in agreement with results obtained through exactdiagonalisation of the matrix M . As the system size is increased beyond the sizes forwhich exact techniques remain tractable, the curves approach those predicted by dGEand DWT. However, it is not possible to discern from the largest system simulated, L = 150, which of these theories is more appropriate.We therefore attempt to extrapolate to the limit L → ∞ by fitting a function ofthe form ε + a L − + a L − + a L − to estimates of the gap at finite system sizes L butfixed α and β . This particular fitting function was found to have the best goodness-of-fit(taking into account the varying number of parameters) out of those that were tried,and furthermore had the most uniform goodness-of-fit as a function of β . The resultsof this extrapolation procedure are shown in Fig. 8. These data appear to exclude thefunctional form for the domain wall theory, although the uncertainties in the estimatedgaps are such that further evidence is needed in order to confidently rule it out. One possible means to reduce the uncertainty in the measured gap is to try differentobservables O ( t ) in the hope that contributions from subdominant decay modes—andin particular the second longest-lived mode that decays at rate λ —are reduced and ynamical Transition in the TASEP Figure 7.
The gap as a function of β at α = 0 . L ranging from L = 40 to L = 150 (from bottom to top) obtained from Monte Carlo simulations. The L → ∞ dGE (solid line) and DWT (dashed line) predictions are shown for comparison:even with the largest system L = 150 simulated, one cannot distinguish between themnumerically. Figure 8.
The gap extrapolated to L = ∞ from the data shown in Fig. 7. Again thedGE (solid line) and DWT (dashed line) predictions are shown for comparison: herethe data are potentially more compatible with the former than the latter. thereby allow a more accurate estimate of λ . We considered two candidates.The first was the autocorrelation function of the total occupancy, O ( t ) = N ( t ) N ( t + t ), where t is some fixed time point in the steady state. The analogue ofthe function R ( t ) is χ ( t ) = h N ( t + t ) N ( t ) i − h N i h N i − h N i (17)where here all averages are taken in the steady state. Numerically, this quantity is ynamical Transition in the TASEP Figure 9.
Interactions between first class (circles) and second class (hexagons)particles. Second class particles hop can forward into empty sites and exchange placeswith first-class particles behind them. Note that in our simulations, at most one second-class particle was present at any instant, and that the behaviour at the boundaries waschosen so as not to affect the motion of first-class particles (see text). slightly more convenient than N ( t ), since one can obtain h N ( t + t ) N ( t ) i by samplingat different t and t from a single simulation run that has reached the steady state.The other observable we investigated was the position of a second class particle[37–39]. Like a first-class particle, a second-class particle hops to the right into vacantsites. However, it may also exchange places with a particle occupying the site to its left.Both of these processes take place at unit rate, and first-class hop to the right at thesame rate, irrespective of whether the site in front of them is empty or occupied by asecond-class particle—see Fig. 9. So as not to affect the entry of first-class particles intothe system, a second-class particle on the left-boundary site is forced to exit if a first-class particle attempts to enter. It is reinserted as soon as the left-boundary site becomesvacant again. Second-class particles are prevented from exiting at the right boundary.In these simulations one can measure the position of the second class particle, x ( t ),starting from an initial condition where each site is occupied by a first-class particlewith density ρ (as described above) except for the left boundary site, which is alwaysinitially occupied by the second-class particle. Here, the analogue of R ( t ) is ξ ( t ) = h x ( t ) i − h x i x (0) − h x i . (18)Gaps obtained from these two functions χ ( t ) and ξ ( t ) are compared with thoseobtained from R ( t ) in Fig. 10 at various system sizes. We see that there are nosystematic differences between these functions at the system sizes studied, and thereforethat they are unlikely to provide improved estimates of the gap in the thermodynamiclimit L → ∞ . We therefore conclude that whilst Monte Carlo simulation data is possiblymore consistent with the de Gier–Essler prediction for the gap than the domain walltheory, the effect of the changing gap on the relaxation of a typical observable is verysmall and hard to disentangle from the noise.
4. Evidence for the dynamical transition from density matrixrenormalisation
Given the difficulties in measuring the gap from Monte Carlo simulations, we nowturn to a fundamentally different numerical approach, namely the density matrixrenormalisation group (DMRG) [27] that was briefly discussed in the introduction. ynamical Transition in the TASEP Figure 10.
The gap as a function of β as obtained by consideration of the relaxationof three different observables. Solid lines show data obtained from the relaxation of thedensity for system sizes L = 10 , , , ,
150 (bottom to top). Alongside the lowerthree curves are plotted data from the stationary density auto-correlation function (17)(dashed lines), and alongside the upper two data from the position of a single second-class particle (18) (large symbols). All three observables yield consistent results.
We begin by recapitulating the essential features of the DMRG approach beforedemonstrating that it does indeed show a dynamical transition along the de Gier–Esslerline.
The density matrix renormalisation group (DMRG) procedure builds an approximationto the transition matrix M that appears in the master equation (2). The basic idea isto repeatedly add sites to the system and renormalise the enlarged transition matrixso that its dimensionality remains constant as the system size grows. In this way,the approximated transition matrix remains sufficiently small that its spectrum (andin particular the gap) can be determined using standard numerical diagonalisationmethods. The success of this procedure relies on finding a reduced basis set ateach renormalisation step that allows the lowest-lying eigenstates to remain wellapproximated as sites are added.Since the TASEP has open boundaries, the extension of the lattice is achieved inthis instance by adding sites to the middle of the system rather than at one of theends [26]. The system is thus divided into two ‘blocks’: the left and right halves of thelattice. It is extended in each iteration by adding a site at the right-hand end of theleft block, and the left-hand end of the right block. In order to understand the resultsthat come out of the procedure, it is worth explaining in a little more detail how therenormalisation is actually achieved—full details are provided in Appendix A.We outline the procedure in terms of the transformation applied to the left block. ynamical Transition in the TASEP Figure 11.
Illustration of the DMRG procedure applied to the left block. Before thetransformation, the state of the block is specified by 1 ≤ p ≤ m and σ = 0 ,
1. Afteradding a new site to the end of the block, the original ℓ sites are renormalised so thatthey are specified by the two numbers 1 ≤ ˜ p ≤ m and ˜ σ = 0 ,
1. The right block is amirror image of the left, and treated the same way.
Suppose that at the start of the iteration, the block comprises ℓ lattice sites. The setof states the block can be in is specified by two ‘quantum’ numbers, p = 1 , . . . , m and σ = 0 ,
1. The quantity p indexes the configuration of the first ℓ − σ theoccupancy of the rightmost site. The first step of the renormalisation is to constructsome basis in the 2 m -dimensional space spanned by p and σ and, crucially, retain only m of these basis vectors. This defines a renormalised index ˜ p = 1 , , . . . , m for the entire ℓ -site block. A new site is then added to the right hand side of the lattice. Theconfiguration of this additional site is specified by the renormalised coordinate ˜ σ . SeeFig. 11. The reason for keeping the rightmost site in the ‘physical’ (site occupancy)basis is that one needs to project the transition matrix for hops from site ℓ to ℓ + 1 ontothe space spanned by ˜ p and ˜ σ . These transition rates are initially only known in thisphysical basis. Thus through this procedure, we obtain a (truncated) representation ofthe transition matrix for an ℓ + 1-site system in a space of the same dimensionality asthat used to describe the ℓ -site system.The same transformation is applied to the right block, the only difference being thatit is the leftmost site that is expressed in the physical basis. Interactions between the twoblocks enter through the construction of the renormalised indices ˜ p as we now describe.First, a transition matrix M ′ for the entire system of 2 ℓ sites is constructed by combiningthe matrices for the two halves, and by adding an interaction term that allows a particlein the left block to hop to the right block. Again, having the internal two sites specifiedin the physical basis helps here. This transition matrix is then diagonalised. However,it is not these eigenstates that are used to perform the renormalisation. Rather, a density matrix , which is a symmetric combination of the two longest-lived eigenstatesof M ′ , is constructed, and eigenvectors of this density matrix are instead used for therenormalisation. The idea is that this form of the density matrix allows the stationaryand longest-lived transient state to be accurately represented in large systems. Thespecific form of the density matrix, and the prescription for obtaining the truncatedbasis set, is given in Appendix A. For further details about the principles behind DMRGand other applications we refer to [27, 40, 41], and especially [29]. ynamical Transition in the TASEP ‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰‰ - - - - - L ¶ Figure 12.
Stable DMRG data for the gap (red crosses) and the exact thermodynamicgap (black line) at α = 0 . β = 0 . We used the DMRG algorithm outlined above to estimate the gap for a givencombination of α and β by starting with an exact diagonalisation of the L = 8 system,and keeping m = 16 eigenvectors of the density matrix in each renormalisation step. Inprinciple, one ought to be able to access arbitrarily large system sizes with this method.In practice—and as was also noted in [26]—the algorithm eventually goes unstable,which is typically manifested through the gap acquiring an imaginary part. We simplyignore data for system sizes where the instability is judged to have kicked in.Although we can access larger system sizes with the DMRG approach than waspossible with Monte Carlo (e.g., L up to about 250, as shown in Fig. 12), it is stillnecessary to extrapolate to the thermodynamic limit, L → ∞ . We follow a similarprocedure to that described in Section 3.1. That is, we specify a finite-size fittingfunction of the form f ( L ) = ε + a L − , and adjust the smallest value of L used for thefit. Again, we use R as a goodness-of-fit measure, i.e., Eq. (16) with t replaced by L .We show in Fig. 13(a) how the goodness of fit varies with the minimum value of L usedin the fit; choosing the optimal (largest) value yields an estimate of the gap that weobserve from Fig. 13(b) that appears more consistent with the de Gier–Essler predictionthan domain wall theory.We show the DMRG estimate of the gap obtained in this way as a function of β in Fig. 14. The agreement with the analytical prediction common to dGE and DWTbelow β < is very good. For larger β the data are scattered around an apparentlyconstant value, indicating the presence of inaccuracies in the DMRG algorithm and/orthe extrapolation of L → ∞ . If we regard these data as independent measurements ofthe same value, we can obtain an estimate of the uncertainty in the constant value ofthe gap above some critical point by simply calculating the standard deviation. We findthe DMRG gap to approach the constant value ε = − . α = 0 .
4. For ynamical Transition in the TASEP ææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææ L min R (a) æææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææææ - - - - - L min ¶ (b) Figure 13. (a) Coefficient of determination R and (b) extrapolated gap, as afunction of L min for DMRG data corresponding to α = 0 . β = 0 .
75. The optimalchoice of L min lies just below 200, for which the extrapolated gap very closely matchesthe de Gier-Essler result and clearly differs from the DWT gap. Figure 14.
DMRG estimates of the gap (red crosses), exact thermodynamic gap(black curve) and DWT gap (green curve), all as a function of β at α = 0 . β c ≈ . comparison, de Gier and Essler predict a constant value for the gap of − . . . . abovethe dynamical transition, while domain wall theory predicts − . . . . . The DMRGmeasurement is then clearly consistent with the de Gier–Essler prediction, whilst thedifference from the domain wall theory value is significant. Taking this result togetherwith the Monte Carlo data we conclude that a dynamical transition does indeed occurat the dGE line, and not where DWT predicts.
5. Discussion of the distinction between static and dynamic subphases
In this paper we have presented numerical evidence, the most convincing of which camefrom DMRG calculations, that a dynamical transition occurs in the TASEP at thedGE line rather than at the subphase boundary. However, this finding in turn raises ynamical Transition in the TASEP
18a number of open questions. For example, does the dGE line—and in particular itsdeparture from the static subphase boundary—have any physical significance? How isthe static subphase boundary manifested in the eigenvalue spectrum? We discuss thedistinction between the two phase diagrams first in general theoretical terms, and thenwith reference to different treatments of the TASEP.
One way to gain a general understanding of nonequilibrium phase transitions is throughthe analogue of partition function for a system of size L , Z L [31, 32]. It has been shownthat for any system governed by Markovian dynamics, it may be written as Z L = Y j =0 ( − λ j ) (19)where ( − λ j ) are the eigenvalues of the Markov transition matrix which defines thedynamics [32,42]. For a finite, irreducible configuration space there is a unique stationarystate corresponding to the largest eigenvalue λ = 0. By the Perron-Frobenioustheorem [43], λ is non-degenerate, therefore the gap is given by the second largesteigenvalue λ and the largest relaxation time is τ = − /λ .From our knowledge of equilibrium phase transitions, we expect a static phasetransition to occur when lim L →∞ L ln Z L exhibits nonanalyticities. We see from (19)that to obtain a static phase transition we require a major restructuring of the eigenvaluespectrum. To be general and explicit let us consider the nonequilibrium analogue of thefree energy density, h L = ln Z L L = 1 L X λ i =0 ln( − λ i ) . (20)Consider an arbitrary eigenvalue λ i in (20), possessing a nonanalyticity at a critical valueof some parameter. Unless a significant number of other eigenvalues in the spectrumconverge onto the same nonanalyticity in the thermodynamic limit, the effect of the λ i nonanalyticity is lost as L → ∞ and (20) and hence h L remain analytic.On the other hand a dynamical phase transition, as defined above, is only concernedwith the eigenvalue λ . Therefore a dynamical phase transition does not necessarilycoincide with a static phase transition. One very simple example would be if the twoleading subdominant eigenvalues λ and λ crossed at some values of a control parameter.Then the gap would behave in a nonanalytic way but Z would be analytic.Let us consider more explicitly the TASEP. From the exact solution of the TASEP[4, 15] the expression for Z L , as defined by (19), is Z L = ( αβ ) L +1 β − α [Φ L ( α ) − Φ L ( β )] (21)where Φ L ( x ) = L X p =1 p (2 L − p − L !( L − p )! x − ( p +1) (22) ynamical Transition in the TASEP L largeΦ L ( x ) ≃ x (1 − x )[ x (1 − x )] L +1 x < L √ πL x = x L √ π (2 x − L / x > . (23)Now consider, for example, the case α < /
2. As we increase β a first-order staticphase transition occurs on the line β = α where the dominant contribution to (21)changes from Φ L ( β ) to Φ L ( α ). Thus at the phase transition h L defined in (20) is non-analytic. The subphase boundary is at β = 1 / L ( β ) changes. As this is a subdominant contributionto ln Z L , h L remains analytic at the sub-phase boundary. At the dGE line there is nonoticeable change in the asymptotic behaviour of Z L and h L is analytic. We illustrate further the distinctions between the static and dynamics phase diagramsby collecting exact results for the TASEP together with the predictions of various levelof approximations such as mean field theory and domain wall theory.
Burgers Equation (Mean Field theory)
We begin by considering a mean fielddescription of the system that is given by the Burgers equation, ∂ t ρ ( x, t ) = − (1 − ρ ) ∂ x ρ + 12 ∂ x ρ , (24)where ρ ( x, t ) is a density field. Although the Burgers equation has been studiedextensively in the literature, we are not aware of any treatment of the case withprescribed boundary reservoir densities. In the appendix we give the solution of (24)subject to the boundary conditions ρ (0) = α and ρ ( L ) = 1 − β , and arbitrary initialcondition.The results for the Burgers gap ε B may be summarised as follows (restrictingourselves to the case α < / β < / ε B = − | α (1 − α ) − β (1 − β ) | (25)For β > / ε B = − (1 − α ) β = 1 /
2. Interestingly, it turns out that there is no change in the density profile atthis value so, in fact, the mean field theory predicts a dynamical transition instead of asubphase boundary at β = 1 /
2. Thus even at the level of mean-field theory, the TASEPexhibits the phenomenon of a dynamic transition that is not accompanied by a statictransition. ynamical Transition in the TASEP Figure 15.
Comparison of the gap for L → ∞ from the spectrum of the viscousBurgers equation (24) to the exact thermodynamic gap, for α = 0 . β c ≈ . Domain wall theory
The predictions of domain wall theory were reviewed inSection 2.4. As we pointed out there, DWT gives a subphase boundary at β = 1 / /L . Exact results
As was discussed in Section 2, the exact results for the stationary state[15, 16] and the spectrum [20, 21] reveal behaviour that is different from both mean-fieldand domain wall theory. In common with the latter, there is a static subphase boundaryat β = . Like the mean field theory, the location of the dynamic transition occurs at adistinct location to the static phase boundary. However, this location α c ( β ) is nontrivialand is not predicted by either of the simpler theories. Moreover, the behaviour of the gapat the static transition point β = α differs between the mean-field and exact theories:in the former, the gap has a cusp, whereas in the latter it varies analytically.
6. Conclusion
In this work we have provided numerical evidence to confirm the existence of thedynamical transition line predicted by de Gier and Essler. Ultimately the evidencecame from the approximate, but quantitatively reliable, DMRG technique.To conclude, we now return to the main open question resulting from this work,namely that of the physical significance of the dynamical transition that takes placealong the de Gier–Essler line. We remark again that, as was seen in Section 3,direct measurement of the gap in Monte Carlo simulations was very difficult. Thenonanalyticity in the leading eigenvalue was barely detectable in the stochasticsimulation data. We were therefore unable to identify, for example, whether the systemrelaxed to stationarity in a fundamentally different way either side of the transition line. ynamical Transition in the TASEP ′ region by simply imposing that the right domain density be ρ + = 1 − β c in this region. However DWT then no longer predicts a static subphaseboundary at β = 1 / α c ( β ) is not predicted.It may be that some knowledge of the form of the eigenvectors associated to thelow lying states will allow an understanding of what occurs at the dynamical transition.It would be of interest to extend the DMRG studies to investigate the form of theeigenfunctions. Also recent analytic work by Simon [44], in which eigenvectors areconstructed for the partially asymmetric exclusion process through the co-ordinate Betheansatz, may help to shed light on the matter. Acknowledgments
R.A.B. thanks the RCUK for the support of an Academic Fellowship. AP would like tothank Fabian Essler for hospitality and discussion, and acknowledges financial supportfrom the EPSRC.
Appendix A. DMRG algorithm for TASEP
In this appendix, we provide a step-by-step description of the DMRG algorithm usedto obtain the results presented in Section 4.In what follows we make use of three elementary transition matrices for the TASEP.With the ordering ( , • ) of single-site configurations we define h l = − α α ! and h r = β − β ! (A.1)for the entry and exit of particles at the left and right boundary sites. In the bulk, andwith the ordering ( , • , • , • • ) of two-site configurations, we define h b = − . (A.2)The DMRG algorithm now proceeds as follows.0. Initialise block transition matrices
We partition the system into two blocks, eachof size ℓ . The transition matrix for the left block can be written as M ( ℓ ) L = h l ⊗ I ⊗ ℓ − + ℓ − X k =1 I ⊗ k − ⊗ h b ⊗ I ⊗ ℓ − k − , (A.3) ynamical Transition in the TASEP I is the 2 × M ( ℓ ) R = ℓ − X k =1 I ⊗ k − ⊗ h b ⊗ I ⊗ l − k − + I ⊗ ℓ − ⊗ h r . (A.4)1. Diagonalise the transition matrix for the entire system
The transition matrix forthe entire system of size 2 ℓ , M (2 ℓ ) , is given by M (2 ℓ ) = M ( ℓ ) L ⊗ I ⊗ ℓ + I ⊗ ℓ − ⊗ h b ⊗ I ℓ − + I ℓ ⊗ M ( ℓ ) R . (A.5)Now we diagonalise M (2 ℓ ) , first using the Arnoldi algorithm [45] (as implementedin the ARPACK package [46] and accessed using Mathematica) to find the leftand right ground state eigenvectors h φ | , | ψ i . We also want to compute the gapand associated eigenvectors | ψ i and h φ | . Following [26], we find this can be donemore accurately by ‘shifting’ the transition matrix through M ′ = M (2 ℓ ) + ∆ | ψ ih φ | sufficiently far that the gap is now the largest eigenvector of M ′ .2. Form reduced density matrices
The reduced density matrices are defined as ρ L = 14 Tr R ( | ψ ih ψ | + | φ ih φ | + | ψ ih ψ | + | φ ih φ | ) (A.6) ρ R = 14 Tr L ( | ψ ih ψ | + | φ ih φ | + | ψ ih ψ | + | φ ih φ | ) . (A.7)Here Tr L and Tr R indicate a trace over the degrees of freedom in the left and rightblocks respectively.To be clear, let {| i i L } and {| j i R } be basis sets for the left and right blocksrespectively. A state vector | ψ i across both blocks can then be expanded as | ψ i = X ij c ij | i i L | j i R . (A.8)The trace operation Tr L ( | ψ ih ψ | ) is then defined asTr L ( | ψ ih ψ | ) = X jj ′ X i c ij c ij ′ ! | j i R h j ′ | R . (A.9)Note that this is an operator acting only on the right block. The correspondingoperation Tr R is defined analogously.3. Diagonalise the density matrices and form a truncated basis set
We first find alleigenvalues and eigenvectors of the symmetric matrices ρ L and ρ R using a densematrix diagonalisation routine from the LAPACK package [47], again accessed usingMathematica. We then form the matrices O L and O R that have as their columnsthe m eigenvectors with largest eigenvalues of ρ L and ρ R respectively. Then weconstruct projectors P L = O L ⊗ I and P R = I ⊗ O R (A.10)that will be used in the renormalisation step below.4. Enlarge system
The left block is enlarged by adding a site at its right end. Thetransition matrix for this enlarged block is f M ( ℓ +1) L = M ( ℓ ) L ⊗ I ⊗ + I ⊗ ℓ − ⊗ h b . (A.11) ynamical Transition in the TASEP f M ( ℓ +1) R = h b ⊗ I ℓ − + I ⊗ ⊗ M ( ℓ ) R . (A.12)5. Renormalise the transition matrices for the enlarged system
The renormalisationis performed by projection the enlarged transition matrices onto the subset ofdensity matrix eigenstates that was retained in step 3. Formally, this is achievedby computing M ( ℓ +1) L = P TL f M ( ℓ +1) L P L and H ( ℓ +1) R = P TR f M ( ℓ +1) R P R . (A.13)Note that after this transformation, the first ℓ sites of the system are no longerrepresented in the ‘physical’ basis (i.e., particle-hole configurations). However, theform of the projectors (A.10) ensures that the last site of the block is representedin the physical basis, which is essential for the expression (A.5) to be valid at eachstage of the renormalisation.6. Return to step 1 , putting ℓ → ℓ + 1. Appendix B. Solution of the Burgers equation with fixed boundarydensities
Here we solve the Burgers equation (i.e., the continuum limit of the TASEP within amean-field theory) ∂ t ρ ( x, t ) = − (1 − ρ ) ∂ x ρ + 12 ∂ x ρ . (B.1)subject to the boundary conditions ρ (0) = α and ρ ( L ) = 1 − β . We use the Cole-Hopftransformation [48], which involves the change of variable ρ ( x, t ) = 12 " ∂∂x ln u ( x, t ) . (B.2)Substituting into (B.1), and integrating with respect to x eventually yields ∂∂t u ( x, t ) = 12 ∂ ∂x u ( x, t ) . (B.3)The boundary conditions on ρ ( x, t ) transform to the conditions u ′ (0 , t ) = (2 α − u (0 , t ) (B.4) u ′ ( L, t ) = (1 − β ) u ( L, t ) . (B.5)We proceed by finding the eigenfunctions φ n ( x ) that satisfy12 d d x φ n ( x ) = λ n φ n ( x ) (B.6)along with the above boundary conditions at x = 0 and x = L . Thus given an initialcondition ρ ( x,
0) the function u ( x, t ) will be given by u ( x, t ) = X n ≥ c n e λ n t φ n ( x ) (B.7) ynamical Transition in the TASEP c n will depend on the initial condition. Inverting the Cole-Hopftransformation gives ρ ( x, t ) = 12 φ ′ ( x ) φ ( x ) + P n> c n e ( λ n − λ ) t φ ( x ) (cid:16) φ n ( x ) φ ( x ) (cid:17) ′ P n ≥ c n e ( λ n − λ ) t φ n ( x ) . (B.8)The final term, involving the ratio of sums, vanishes as t → ∞ , as all λ n with n > λ . Thus the stationary profile is ρ ( x ) = 12 " φ ′ ( x ) φ ( x ) . (B.9)We can define the Burgers gap, ε B via the asymptotic rate of the exponential decayof the density profile to its stationary form, viz, ε B = lim t →∞ dd t ln[ ρ ( x, t ) − ρ ( x )] = λ − λ . (B.10)The eigenfunctions given by (B.6) themselves take the form φ ( x ) = A e γx + B e − γx (B.11)whence λ = γ . The allowed values of γ are quantised due to the boundary conditions(B.4) and (B.5). We first look for solutions with λ >
0, i.e., real γ . Imposing theboundary conditions leads to the pair of equations γ ( A − B ) = (2 α − A + B ) (B.12) γ ( A e γL − B e − γL ) = (1 − β )( A e γx + B e − γx ) . (B.13)These equations are consistent only iftanh( γL ) = 2(1 − α − β ) γ (1 − α )(1 − β ) + γ . (B.14)In the limit L → ∞ , the left-hand side approaches ±
1. The resulting equation has asolution γ = 1 − α if α < and a solution γ = 1 − β if β < .Solutions with a negative λ have purely imaginary γ = ik ( k real), and by followingthrough the same procedure, one finds that the allowed values of k are the roots oftan( kL ) = 2( α + β − kk − (2 α − β − . (B.15)In the large- L limit, we find that k ∼ /L , and hence the associated eigenvalues vanishas 1 /L in the thermodynamic limit.We can now state the L → ∞ behaviour of the Burgers gap ε as α and β are varied.When both α and β are smaller than , there are two isolated positive eigenvalues λ = (2 α − , (2 β − . In this region, then, the gap behaves as ε B = − | (2 α − − (2 β − | − | α (1 − α ) − β (1 − β ) | . (B.16)Like the exact gap (4), this vanishes along the HD-LD coexistence line α = β < .However, unlike the exact gap, the mean-field gap has a nonanalyticity along this line. ynamical Transition in the TASEP α or β is increased above , one of the two solutions with positive λ ceases toexist, and the second-largest eigenvalue λ = 0 in the thermodynamic limit. Hence then,the size of the gap is simply equal to the value of the largest eigenvalue. That is, withinthe low-density phase α < , β > , ε B = − (1 − α ) α > , β < ε B = − (1 − β ) . (B.18)The Burgers gap thus exhibits nonanalyticities along the coexistence line α = β < ,and along the lines α < , β = and α = , β < . Thus the dynamic phase diagramfor the Burgers equation appears the same as the static phase diagram, Fig. 1(a), exceptthat the subphase boundaries have become dynamical transition lines. References [1] T. M. Liggett.
Stochastic interacting systems: contact, voter, and exclusion processes . SpringerVerlag, Berlin, 1999.[2] B. Derrida. An exactly soluble non-equilibrium system: the asymmetric simple exclusion process.
Phys. Rep. , 301:65, 1998.[3] O. Golinelli and K. Mallick. The asymmetric simple exclusion process: an integrable model fornon-equilibrium statistical mechanics.
J. Phys. A: Math. Gen. , 39:12679, 2006.[4] R. A. Blythe and M. R. Evans. Nonequilibrium steady states of matrix-product form: a solver’sguide.
J. Phys. A: Math. Theor. , 40:R333, 2007.[5] V. Popkov, L. Santen, A. Schadschneider, and G. M. Sch¨utz. Empirical evidence for a boundary-induced nonequilibrium phase transition.
J. Phys. A: Math. Gen. , 34:L45, 2001.[6] C. T. MacDonald, J. H. Gibbs, and A. C. Pipkin. Kinetics of biopolymerization on nucleic acidtemplates.
Biopolymers , 6:1, 1968.[7] Y. Aghababaie, G. I. Menon, and M. Plischke. Universal properties of interacting Brownianmotors.
Phys. Rev. E , 59:2578, 1999.[8] S. Klumpp and R. Lipowsky. Traffic of molecular motors through tube-like compartments.
J.Stat. Phys. , 113:233, 2003.[9] S. Chatterjee. Diffusion of a hydrocarbon mixture in a one-dimensional zeolite channel: Anexclusion model approach.
Micropor. Mesopor. Mat. , 125:143, 2009.[10] K. Sugden, M. R. Evans, W. C. K. Poon, and N. D. Read. Model of hyphal tip growth involvingmicrotubule-based transport.
Phys. Rev. E , 75:031909, 2007.[11] J. B. Martin and P. A. Ferrari. Stationary distributions of multi-type totally asymmetric exclusionprocesses.
Ann. Probab. , 35:807, 2007.[12] M. Evans, P. Ferrari, and K. Mallick. Matrix representation of the stationary measure for themultispecies TASEP.
J. Stat. Phys. , 135:217, 2009.[13] J. Krug. Boundary-induced phase transitions in driven diffusive systems.
Phys. Rev. Lett. , 67:1882,1991.[14] B. Derrida, E. Domany, and D. Mukamel. An exact solution of a one-dimensional asymmetricexclusion model with open boundaries.
J. Stat. Phys. , 69:667, 1992.[15] B. Derrida, M. R. Evans, V. Hakim, and V. Pasquier. Exact solution of a 1D asymmetric exclusionmodel using a matrix formulation.
J. Phys. A: Math. Gen. , 26:1493, 1993.[16] G. Sch¨utz and E. Domany. Phase transitions in an exactly soluble one-dimensional exclusionprocess.
J. Stat. Phys. , 72:277, 1993. ynamical Transition in the TASEP [17] V. Popkov and G. M. Sch¨utz. Steady-state selection in driven diffusive systems with openboundaries. Europhys. Lett. , 48:257, 1999.[18] J. S. Hager, J. Krug, V. Popkov, and G. M. Sch¨utz. Minimal current phase and universal boundarylayers in driven diffusive systems.
Phys. Rev. E , 63:056110, 2001.[19] S. Mukherji and S. M. Bhattacharjee. Nonequilibrium criticality at shock formation in steadystates.
J. Phys. A: Math. Gen. , 38:L285, 2005.[20] J. de Gier and F. H. L. Essler. Bethe ansatz solution of the asymmetric exclusion process withopen boundaries.
Phys. Rev. Lett. , 95:240601, 2005.[21] J. de Gier and F. H. L. Essler. Exact spectral gaps of the asymmetric exclusion process with openboundaries.
J. Stat. Mech.: Theor. Exp. , P12011, 2006.[22] J. de Gier and F. H. L. Essler. Slowest relaxation mode in the partially asymmetric exclusionprocess.
J. Phys. A: Math. Theor. , 41:485002, 2008.[23] R. I. Nepomechie and F. Ravanini. Completeness of the Bethe ansatz solution of the open XXZchain with nondiagonal boundary terms.
J. Phys. A: Math. Gen. , 36:11391, 2003.[24] R. I. Nepomechie. Functional relations and Bethe ansatz for the XXZ chain.
J. Stat. Phys. ,111:1363, 2003.[25] M. Dudzinski and G. M. Schutz. Relaxation spectrum of the asymmetric exclusion process withopen boundaries.
J. Phys. A: Math. Gen. , 33:8351, 2000.[26] Z Nagy, C Appert, and L Santen. Relaxation times in the ASEP model using a DMRG method.
J. Stat. Phys. , 109:623, 2002.[27] R. M. Noack and S. R. White.
Density-Matrix Renormalization: A New Numerical Method inPhysics , chapter 2. Lecture Notes in Physics. Springer, 1999.[28] S. R. White. Density matrix formulation for quantum renormalization groups.
Phys. Rev. Lett. ,69:2863, 1992.[29] U. Schollw¨ock. The density-matrix renormalization group.
Rev. Mod. Phys. , 77:259, 2005.[30] E. Carlon, H. Henkel, and U. Schollw¨ock. Density matrix renormalization group and reaction-diffusion processes.
Euro. Phys. J. B , 12:99, 1999.[31] R. A. Blythe and M. R. Evans. Lee-Yang zeros and phase transitions in nonequilibrium steadystates.
Phys. Rev. Lett. , 89:080601, 2002.[32] M. R. Evans and R. A. Blythe. Nonequilibrium dynamics in low-dimensional systems.
PhysicaA , 313:110, 2002.[33] A. B. Kolomeisky, G. M. Sch¨utz, E. B. Kolomeisky, and J. P. Straley. Phase diagram of one-dimensional driven lattice gases with open boundaries.
J. Phys. A: Math Gen. , 31:6911, 1998.[34] B. Derrida, M. R. Evans, and K. Mallick. Exact diffusion constant of a one-dimensional asymmetricexclusion model with open boundaries.
J. Stat. Phys. , 79:833, 1995.[35] D. T. Gillespie. Exact stochastic simulation of coupled chemical reactions.
J. Phys. Chem. ,81:2340, 1977.[36] T. P. Ryan.
Modern Regression Methods . Wiley, Hoboken, NJ, 1997.[37] C. Boldrighini, G. Cosimi, S. Frigio, and M. Grasso Nunes. Computer simulation of shock wavesin the completely asymmetric simple exclusion process.
J. Stat. Phys. , 55:611, 1989.[38] P. A. Ferrari. Shock fluctuations in asymmetric simple exclusion.
Probab. Theory Rel. , 91:81,1992.[39] B. Derrida, S. A. Janowsky, J. L. Lebowitz, and E. R. Speer. Exact solution of the totallyasymmetric simple exclusion process: Shock profiles.
J. Stat. Phys. , 73:813, 1993.[40] S. R. White. Strongly correlated electron systems and the density matrix renormalization group.
Phys. Rep. , 301:187, 1998.[41] R. M. Noack and S. R. Manmana. Diagonalization- and numerical renormalization-group-basedmethods for interacting quantum systems. arXiv:cond-mat/0510321v1 , 2005.[42] R. A. Blythe.
Nonequilibrium phase transitions and dynamical scaling regimes . PhD thesis,University of Edinburgh, 2001.[43] E. Senata.
Non-negative matrices and Markov Chains . Springer Series in Statistics. Springer, ynamical Transition in the TASEP New York, 1981.[44] D. Simon. Construction of a coordinate bethe ansatz for the asymmetric exclusion process withopen boundaries.
J. Stat. Mech.: Theor. Exp. , P07017, 2009.[45] G. H. Golub and C. F. van Loan.
Matrix Computations