Asymmetric Exclusion Processes with Disorder: Effect of Correlations
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b Asymmetric Exclusion Processes with Disorder: Effect of Correlations
M. Ebrahim Foulaadvand ∗ Department of Nano-Science, Institute for Studies in Theoretical Physics and Mathematics (IPM),P.O. Box 19395- 5531,Tehran, Iran and Department of Physics,Zanjan University, P.O. Box 45195-313, Zanjan, Iran.
Anatoly B. Kolomeisky
Department of Chemistry, Rice University, Houston, TX 77005 USA.
Hamid Teymouri
Department of Physics, Zanjan University, P.O. Box 45195-313, Zanjan, Iran. (Dated: October 31, 2018)Multi-particle dynamics in one-dimensional asymmetric exclusion processes with disorder is inves-tigated theoretically by computational and analytical methods. It is argued that the general phasediagram consists of three non-equilibrium phases that are determined by the dynamic behavior atthe entrance, at the exit and at the slowest defect bond in the bulk of the system. Specifically, weconsider dynamics of asymmetric exclusion process with two identical defect bonds as a functionof distance between them. Two approximate theoretical methods, that treat the system as a se-quence of segments with exact description of dynamics inside the segments and neglect correlationsbetween them, are presented. In addition, a numerical iterative procedure for calculating dynamicproperties of asymmetric exclusion systems is developed. Our theoretical predictions are comparedwith extensive Monte Carlo computer simulations. It is shown that correlations play an importantrole in the particle dynamics. When two defect bonds are far away from each other the strongestcorrelations are found at these bonds. However, bringing defect bonds closer leads to the shift ofcorrelations to the region between them. Our analysis indicates that it is possible to develop a suc-cessful theoretical description of asymmetric exclusion processes with disorder by properly takinginto account the correlations.
I. INTRODUCTION
In recent years a significant attention has been devotedto investigation of low-dimensional asymmetric simpleexclusion processes (ASEPs) [1, 2, 3, 4, 5]. They play acritical role for understanding fundamental properties ofnon-equilibrium phenomena in Chemistry, Physics andBiology. ASEPs have been widely utilized for descrip-tion of traffic phenomena [4], kinetics of biopolymeriza-tion [6], protein synthesis [7, 8, 9, 10], and biologicaltransport of motor proteins [11, 12]. The advantage ofusing asymmetric exclusion processes for studying mech-anisms of non-equilibrium phenomena is due to the factthat some homogeneous versions of ASEPs can be solvedexactly via matrix-product approach and related meth-ods [1, 5, 13]. In addition, understanding of processes inASEPs can be achieved by utilizing a phenomenologicaldomain wall approach [14]. In order to have a more real-istic description of different non-equilibrium phenomenaASEPs with inhomogeneous distribution of rates are re-quired. However, there is a limited number of studiesdealing with ASEPs with disorder in the transition ratesat sites (static impurities) [8, 9, 15, 16, 17, 18, 19, 20,21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33] and withdisorder associated to particles’s hopping rates (moving ∗ Corresponding author: [email protected] impurities) [34, 35, 36, 37, 38]. In this case exact solu-tions are not obtained, and extensive Monte Carlo com-puter simulations and approximate theories are utilizedin order to understand particle dynamics. Disorder hasa strong effect on the behavior of ASEPs. Even a singledefect bond far away from the boundaries lead to dra-matic effects in the stationary properties both in closed[15, 18] and open boundary conditions [19]. It was shownrecently that the dynamics of ASEPs is also influenced byseveral defects that are close to each other [8, 10, 32], al-though the mechanism of this phenomenon is not well un-derstood. This interaction between defects is importantfor understanding several biological transport phenom-ena [8, 10]. Recently the particular case of two defectshas been extensively investigated by Monte Carlo simu-lations [10]. It has been shown that the system currentexhibits a notable dependence on the distance betweendefects with equal hopping rates. Moreover, it was foundthat the density profile is linear between defects whichmarks the existence of wandering shock between defects[10]. The case of two defective sites with equal rate hasbeen generalized to include extended objects [31]. Theo-retical efforts to analyze ASEPs with disorder have beenmostly directed to the cases with a single or few defects[8, 19, 32]. In Ref. [19] ASEP with open boundaries andwith a local inhomogeneity in the bulk has been investi-gated by arguing that the defect bond divides the systeminto two coupled homogeneous ASEPs. This theoreticalapproach can be called a defect mean-field (DMF) be-cause the mean-field assumptions are made only at theposition of local inhomogeneity. Although a good agree-ment with computer simulations has been found, therewere significant deviations in statistical properties of thephase with the maximal current that was attributed tothe neglect of correlations at the defect bond in the pro-posed theory [19]. A related approach called interact-ing subsystem approximation (ISA) has been proposed inRef. [32] for ASEPs with a single defect or several con-secutive defects (bottleneck). Here it was suggested thatdue to the defect bonds there are 3 segments in the sys-tem: two homogeneous ASEPs are coupled by a segmentthat includes all sites that surround defect bonds. Ex-plicit results have been used inside the segments, andmean-field assumptions have been utilized for particledynamics between the segments. A better agreementwith Monte Carlo computer simulations has been found,and the method was also successfully applied to describeinteractions of defects with boundaries. It was arguedthat ISA can be used for analyzing properties of generalASEPs with disorder [32, 33]. However, ISA has not beenapplied for the systems with 2 defects at finite distancesfrom each other, and because of this observation it isdifficult to apply ISA for understanding mechanisms ofmore complex inhomogeneous asymmetric exclusion pro-cesses. A slightly different method of calculations hasbeen proposed by Chou and Lakatos [8], who applied a finite segment mean-field theory (FSMFT). According tothis approach, the segment of finite length n that cov-ers the defect and surrounding sites is considered, andits dynamics is fully described by solving explicitly foreigenvectors of the corresponding transition rate matrix.The segment is then coupled in the mean-field fashion tothe rest of the system. However, this approach becomesnumerically quite involved for cluster sizes larger than ≈
20, and it also limits its applicability. Different stud-ies of asymmetric exclusion processes with disorder pointout to importance of correlations in the system. It is rea-sonable to expect that correlations are stronger near theslow defect sites. However, it is not clear how far fromthe local inhomogeneity and how fast these correlationsdecay. In addition, it is also unclear how correlationsfrom two close defects affect each other. The goal of thispaper is to investigate the role of correlations in dynamicsof ASEPs with disorder. By analyzing several analyticalapproaches in combination with extensive Monte Carlocomputer simulations it will be shown that a successfuldescription of disordered driven diffusive systems can beachieved by properly accounting for correlations near thedefect bonds.
II. MODEL AND THEORETICALDESCRIPTION
We investigate a totally asymmetric simple exclusionprocesses with disorder. In the one-dimensional latticethe particle at the site i can jump forward with the rate p i if the next site i + 1 is unoccupied, otherwise it stays at the same place. The particle can enter the systemwith the rate α if the site is empty, and it can also exitthe lattice with the rate β . When all p i = 1 we havea homogeneous ASEP for which dynamic properties areknown explicitly from exact solutions [1, 3]. ASEPs withdisorder correspond to the situation when there is inho-mogeneities in the transition rates, and p i are drawn fromarbitrary distributions. Numerous theoretical and com-putational studies indicate that in the limit of large timesthe dynamics in the system can be determined by com-paring entrance rate, exit rates and the transition rateat the slowest defect bond [19, 33]. This observation hasa significant consequence for properties of ASEPs withdisorder, yielding a generic phase diagram with 3 phases.When the entrance is a rate-limiting process the systemcan be found in the low-density phase, while for slow ex-iting the high-density phase governs the system. If therate-limiting process is the transition via the slowest de-fect bond the system is in the maximal-current phase. Inthis maximal-current phase, a segregation of density pro-file into macroscopic high and low regions occurs at thelocation of slowest defect bond. Other defects only per-turb the density profile on a local scale. However, whenthe number slowest defect bonds exceeds two or more theabove picture needs modification. Furthermore, the pre-vious studies on disordered ASEPs lack investigations oncorrelation effects induced by defects. To address thesequestions, we analyze the simplest model with 2 identi-cal defects in the bulk of the system far away from theboundaries. It was shown earlier [32] that positioningof the slow defects close to the boundaries leads only torescaling of the effective entrance and/or exit rates, andwe will not consider this possibility in this paper. Notethat in this paper we are using terms of defect bonds anddefect sites. To clarify, we define the defect site as thesite i from which the particle hopes to the site i + 1 withthe rate q <
1. Correspondingly, the bond connectingsites i and i + 1 is a defect one. A. Defect Mean-Field Theory
Consider a totally asymmetric exclusion processes withopen boundaries and with 2 slow defective sites at i = d and i = d at a distance d with d − d = d (separated by d − q <
1, in all othersites the hopping rate is equal to one. It can be seen thattwo defects divide the system into three segments.The particle dynamics inside each segment can be cal-culated exactly, however, it is assumed that there are nocorrelations between the segments. If entrance to the lat-tice is the slowest process then the system can be foundin low-density (LD) phase with stationary current andbulk densities given by J = α (1 − α ) , ρ bulk = α. (1) α β q qd d d FIG. 1: Fig.1: Schematic of ASEP with two defective sitesat i = d and i = d separated by d − d − d = d . The reduced hopping rates at each defect isequal to q . In normal sites the hoping rates are one. Similarly, when the exit becomes a bottleneck processthe system is in high-density (HD) phase with J = β (1 − β ) , ρ bulk = 1 − β. (2)Note that in both phases particle densities near the de-fect bonds will deviate from the bulk values. The moreinteresting case is when the dynamics in the system isgoverned by transitions via local inhomogeneities. Inthis phase, which has the maximal current, we expectto have density phase segregation similar to the case asingle defect ASEPs. We emphasize that all our inves-tigation in this paper is on this maximal-current phase.First consider a lattice segment after the second defect.The dynamics in this part of the system is controlledby the entrance of particle via the defect, then it has alow-density profile with unknown bulk density ρ ∗ < / J = ρ ∗ (1 − ρ ∗ ), the bulk density in this seg-ment is equal to 1 − ρ ∗ . The region between two defectscan be viewed as asymmetric exclusion process on a fi-nite lattice with d sites. The effective entrance and exitrates to this segment can be easily evaluated using ourmean-field assumptions, α eff = β eff = q (1 − ρ ∗ ) . (3)The stationary properties of the lattice segment with d sites between the defects can be evaluated explicitly byutilizing exact results for finite-size ASEPs [5]. Specifi-cally, the particle flux is given by J ( α, β, d ) = R d − (1 /β ) − R d − (1 /α ) R d (1 /β ) − R d (1 /α ) , (4)where the function R d ( x ) is defined as R d ( x ) = d +1 X p =2 ( p − d − p )! d !( d + 1 − p )! x p . (5) To understand the density profile in the segment wecan use a domain-wall picture of asymmetric exclusionprocesses [14]. Since the entrance and exit rates are thesame [see Eq. (3)], the domain wall that separates high-density and low-density blocks performs an unbiased ran-dom, leading to a linear density profile with a positiveslope. Explicit expressions for particle densities can alsobe found in Ref. [5]. The full description of dynamicsin ASEPs with two defects is obtained by solving for theunknown parameter ρ ∗ . It can be done by applying thecondition of stationarity in the particle flux, J = ρ ∗ (1 − ρ ∗ ) = J ( α eff , β eff , d ) . (6)This equation can always be solved analytically or nu-merically exactly for any number of sites between localinhomogeneities, leading to stationary particle currentsand density profiles. it is important to note that there isa particle-hole symmetry in the system because defectsare far away from the boundaries. To illustrate our ap-proach let us calculate dynamic properties of ASEPs withtwo defects for several values of the parameter d . First,let us analyze the simplest case of d = 1 with consecutivedefects in the bulk. It can be shown that for this system J ( α, β, d = 1) = αβα + β . (7)Then Eq. (6) can be written as ρ ∗ (1 − ρ ∗ ) = q (1 − ρ ∗ ) / , (8)which produces simple expressions for the bulk densityand the particle current, ρ ∗ = q/ , J = q (2 − q ) / . (9)The density l at the site between the defects can alsobe found from the condition that the flux via this site, J = ql (1 − ρ ∗ ), should be equal to the flux through othersegments, and this yields l = ρ ∗ /q = 1 / . (10)This result could also be obtained from the particle-hole symmetry arguments. Note that for q = 1 we ob-tain ρ ∗ = 1 / J = 1 / d = 2 thereare two lattice sites between the defects, and stationaryproperties of this system can also be obtained analyti-cally. From Eq. (5) one can easily derive R ( x ) = x , R ( x ) = x + x , (11)which produces the following expression for the currentin the segment between the defects J ( α, β, d = 2) = α + β α + β + α + β + αβ . (12)Using the expression for the effective entrance and exitrates for the segment between the inhomogeneities [seeEq. (3)], the condition for the stationary current leadsto ρ ∗ (1 − ρ ∗ ) = 2 q (1 − ρ ∗ )3 + 2 q (1 − ρ ∗ ) . (13)This quadratic equation can be solved, and taking thephysically reasonable root we obtain ρ ∗ = 2 q + 3 − p q − q q ; (14) J = 8 q − q − p q − q q . (15)It can be checked that for q = 1 these equations reduceto expected relations ρ ∗ = 1 / J = 1 /
4. We can alsocalculate the densities l and l at the sites between thedefects. Because of the particle-hole symmetry one canargue that l = 1 − l , (16)and the density at the first site can be found by ana-lyzing the current via the first defect, J = q (1 − ρ ∗ )(1 − l ) = ρ ∗ (1 − ρ ∗ ) . (17)Then we have l = 1 − ρ ∗ q = 4 q − q − p q − q q . (18)We have solved equation (6) for the case d = 3. In thiscase we have: R ( x ) = 2 x + 2 x + x , (19)After some lengthy but straightforward algebra we ar-rive at the following cubic equation for ρ ∗ :4 q ( ρ ∗ ) − (8 q +6 q )( ρ ∗ ) +(6 q +6 q +4) ρ ∗ − q (3+2 q ) = 0 . (20)For brevity we avoid writing the answer explicitly. An-alytical results for ASEP with 2 defects can also be ob-tained in the limit of very large distances between the q DMF, d = 1DMF, d = 2DMF, d = 3MC, d = 1MC, d = 2MC, d = 3MF, large d J FIG. 2: Fig.2: (Colour online) J vs q for d = 1 , , inhomogeneities ( d ≫ α eff = β eff ). This correspondsto a linear density profile for the segment between thedefects. Then it leads to the following expression for thecurrent ρ ∗ (1 − ρ ∗ ) = q (1 − ρ ∗ ) [1 − q (1 − ρ ∗ )] , (21)and finally we obtain ρ ∗ = q q , J = q (1 + q ) . (22)These results are identical to stationary properties ofASEP with only one local inhomogeneity far away fromthe boundaries obtained using DMF approximation [19],suggesting that 2 defects at large distances do not affecteach other [8]. For a general d equation (6) leads to apolynomial equation of order d for the unknown ρ ∗ . For d > J as a function of q for d = 1 , J is an increasingfunction of both q and d . DMF notably underestimatesthe current in comparison to the MC simulation espe-cially in the intermediate values of q . B. Interacting Subsystem Approximation
Interacting subsystem approximation (ISA) is anothermethod of calculating stationary properties of ASEPswith a single defect or a single bottleneck developed byGreulich and Schadschneider [32]. It can be easily ex-tended to the case of asymmetric exclusion processes with2 defects separated by d lattice sites. Similarly to DMF α β d l l q ql l segment 3 segment 5segment 1 segment 2 segment 4 FIG. 3: Fig.2: (Colour online) Fig.3: (Colour online) inter-acting subsystems connected via mean-field assumption. Ford > this method divides the lattice into several segments.Particle dynamics inside the segments is treated exactly,while between the segments mean-field assumptions aremade. ISA differs from DMF in the defining of segments.In DMF the position of defects separates different parts,and there is always 3 segments in the system. In ISAthe sites that are connected by the defect bond are puttogether in one segment, as shown in Fig. 3.For d = 1 there are also 3 parts in the lattice, andthe middle segment has 3 sites. For d = 2 there are4 segments and 2 middle segments (with 2 lattice siteseach) border each other. For any larger distance betweenlocal inhomogeneities ISA assumes 5 segments: see Fig.3. Note that the size of the middle segment is equal to d − d > − ρ ∗ and ρ ∗ correspondingly. Let us define l and l as the probabilities to find the particles at thecorresponding sites of the segment around the first de-fect. Similarly, l and l describe densities in the segmentaround the second defect bond. For the middle segmentwith d − x i for i = 1 , · · · , d − i -th site of this segment. As forDMF approach, the existing particle-hole symmetry sim-plifies calculations significantly. Specifically, it suggeststhat l = 1 − l , l = 1 − l , x i = 1 − x d +1 − i . (23)The overall particle current in the system can be writ-ten as J = ρ ∗ (1 − ρ ∗ ) , (24)while due to the mean-field assumptions the currentbetween the first and the second segments is equal to J = (1 − l )(1 − ρ ∗ ) = α (1 − ρ ∗ ) , (25)where α is the effective rate to enter the second seg-ment. At large times we expect to find the system in thestationary state, i.e., J = J , yielding1 − α = l = (1 − ρ ∗ ) . (26)The current between segments 2 and 3 can be pre-sented in the several ways, J = l (1 − x ) = β l = α (1 − x ) , (27)with β being the effective exit rate from the segment2, while α is the effective rate to enter the segment 3.When the system reaches stationary phase, J = J , andwe obtain β = 1 − x , α = ρ ∗ (1 − ρ ∗ )1 − x . (28)Because of the particle-hole symmetry the effective en-trance and exit rates from the segment 3 are the same, α = β . The particle current via the ASEP segmentwith N sites and with entrance and exit rates α and β ,respectively, J ( α, β, N ), can be calculated explicitly [5].Then to obtain stationary properties of ASEP with 2 de-fects in the maximal-current phase the following systemof equations should be solved, ρ ∗ (1 − ρ ∗ ) = qJ ( 1 − ρ ∗ q , − x q , ρ ∗ (1 − ρ ∗ ) = J ( ρ ∗ (1 − ρ ∗ )1 − x , ρ ∗ (1 − ρ ∗ )1 − x , d − . (29)where ρ ∗ and x are 2 unknown variables. The expres-sion on the right side of the first equation describes thecurrent inside the segment 2 and 4. Because the hoppingrate is q <
1, the effective entrance and exit rates mustbe rescaled by the same factor. The right side of the sec-ond equation gives the current inside the segment 3. Theapplication of ISA for d = 1 and d = 2 cases is different.In the case of 2 consecutive defect bonds the system isdivided only in 3 segments. The middle segment has 3sites that surround defect bonds. In the HD/LD phasethe effective entrance rate is α = 1 − ρ ∗ , and the sta-tionary properties can be obtained by solving only oneequation ρ ∗ (1 − ρ ∗ ) = qJ ( 1 − ρ ∗ q , − x q , . (30)Using Eqs. (4) and (5) for the middle segment withequal entrance and exit rates gives us ρ ∗ (1 − ρ ∗ ) = q (1 − ρ ∗ ) [2(1 − ρ ∗ ) + 3 q ]2 [2(1 − ρ ∗ ) + 3 q (1 − ρ ∗ ) + 2 q ] , (31)which can be simplified into the following expression,4( ρ ∗ ) − q +4)( ρ ∗ ) +4( q +1) ρ ∗ − q (3 q +2) = 0 . (32)This cubic equation can be solved explicitly, yielding ρ ∗ = (cid:2) q + 4 + (4 − q ) /D + D (cid:3) / , (33)where D = h − q − q + 3 √ p q − q − q + 28 q i / . (34)ISA also works differently in the case of d = 2. Thereare 4 segments in the system, and because of the neglectof correlations between segments 2 and 3 we have l (1 − l ) = l = ρ ∗ (1 − ρ ∗ ) . (35)Then the effective entrance rate to the segment 2 is α = 1 − ρ ∗ , and the effective exit rate is equal to β = l = p ρ ∗ (1 − ρ ∗ ). The unknown parameter ρ ∗ isdetermined from the equation for the stationary current, ρ ∗ (1 − ρ ∗ ) = qJ ( α q , β q , . (36)Substituting the values of the effective boundary ratesand utilizing Eq. (12) we obtain ρ ∗ p − ρ ∗ + ( ρ ∗ ) / = q p − ρ ∗ . (37)which can be recast in the form of a cubic equation:2( ρ ∗ ) − (2 q + 1)( ρ ∗ ) + ( q + 2 q ) ρ ∗ − q = 0 , (38)This equation can be solved analytically but for brevitywe do not write the solution. It can also be shown that inthe limit of d ≫ ρ ∗ = h q + 2 − p q − q + 4 i / ,J = h q − q + 3 q p q − q + 4 i / . (39)Let us now exhibit the dependence of J on q in the ISAmethod. In figure (4) we have drawn J vs q for d = 1 , d we have considered. q ISA, d = 1ISA, d = 2ISA, d = largeMC, d = 1MC, d = 2MF, d = large J FIG. 4: Fig.4: (Color online) Current vs q for d = 1 , d obtained within ISA method and MC simulation. Insimulations we have taken α = β = 0 . III. CORRELATIONS NEAR DEFECTA. Monte Carlo Simulations
In this section we aim to investigate correlations inthe vicinity of defects. We restrict ourselves to adjacenttwo-point correlations and will present our results for thegeneral two-point and multi-point correlations in a futurework. Let us now introduce the normalized connectedtwo-point correlation function C i between the neighbour-ing sites i and i + 1. This quantity is defined as follows: C i = h τ i τ i +1 i − h τ i ih τ i +1 i p h τ i i − h τ i i q h τ i +1 i − h τ i +1 i i = 1 , · · · , L − . (40)The function C i measures the correlation and it liesbetween − d = 10 and 100 each for threevalues of q . The system size is L = 500 and we havetaken α = β = 0 . T Monte Carlo steps. Each step consists of L moves. Ineach move, we randomly choose a site and update itsstatus according to ASEP rules described above. We dis-card the first T steps to ensure reaching steady state, andwe have accumulated data separated by 10 MC steps toavoid any possible temporal correlations. The value of T is taken 10 in our simulations.Two defects are symmetrically placed with respect tochain mid point. We observed that correlations are largein sites between the defects. There is a rather stronganti-correlation in the sites immediately after the firstdefect and before the second defect. The correlations are site number
240 250 260-0.5-0.4-0.3-0.2-0.100.10.20.30.40.50.60.7 q = 0.1q = 0.3q = 0.5 d = 10 C i site number
150 200 250 300 350-0.6-0.5-0.4-0.3-0.2-0.100.10.20.30.40.50.60.70.8 q = 0.1q = 0.3q = 0.5 d = 100 C i FIG. 5: Fig.5: (Color online) Profile of normalized correlationfunction at d = 10 (top) and d = 100 (bottom) for q =0 . , . , . growing up for middle sites where the maximum valueis achieved. It can be seen that correlations are greaterwhen d is increased. This is unexpected and counterintu-itive because increasing the distance between the defectsreduces their interaction. It has been observed via MCsimulations that when d is increased the current reachesasymptotically to its mean-field value J MF = q (1+ q ) [10, 31]. Therefore, one expects the correlations to ex-hibit a reducing behavior with respect to distance d butthis is not observed in our simulations. To have a deeperunderstanding, we have sketched the behavior of correla-tion profile upon varying the distance d for q = 0 . q = 0 . q , increasing the distance d be-tween the defects gives rise to enhancement of corre-lations/antocorrelations. For instance, the correlationvalue in the middle point rises from roughly 0.5 at small d ∼
10 to 0.65 for d ∼ d larger than 100. Moreover, the correlations are al-ways greater than anti-correlations. By increasing q , thecorrelations/anti-correlations are notably reduced in val-ues. This is expected since in the limit of homogeneousASEP where q → site number
150 200 250 300 350-0.6-0.5-0.4-0.3-0.2-0.100.10.20.30.40.50.60.70.8 d= 10d= 30d= 70d= 100d= 130 q = 0.1 C i site number
150 200 250 300 350-0.4-0.3-0.2-0.100.10.20.30.4 d = 10d = 30d = 70d = 100d = 130 q = 0.3 C i FIG. 6: Fig.6: (Color online) Profile of normalized correlationfunctions at various values of d for q = 0 . q = 0 . small. Here we wish to make a pause and have a discus-sion on correlations in normal ASEPs. In fact the middlesegment between two defects can be regarded as an ASEPchain with length d with equal entrance and exit rates.To the best of our knowledge, correlations in ASEP withrandom sequential update, has only been analytically dis-cussed by Derrida and Evans who obtained exact analyt-ical expression for a general two-point function and madea conjecture to generalize their findings to n -point func-tion [39]. Their study was restricted to the special case α = β = 1 and they found that long range correlationspersist in the bulk which was attributed as a boundaryeffect. In order to see if the large value of the connectedtwo-point function survives in the normal ASEP withequal entrance and exit rates, we performed MC simula-tions. Our results show that when α = β , the profile of C i reaches a small constant (almost zero) in the bulk.The correlations become large near boundaries. Thisboundary behaviour depends on whether α = β < . α = β > .
5. Figure (7) illustrates this aspect.We recall that correlations in other types of updatesuch as parallel updating has been discussed in [16, 17].It is worthwhile to examine the behavior of density pro-file between defects. Dong et al have recently shown via site number
100 200 300 400 500-0.1-0.0500.050.10.15 C i α = β = 0.3α = β = 0.8 FIG. 7: Fig.7: (Color online) Profile of normalized correlationfunctions in a normal ASEP chain with α = β . site number den s i t y
150 200 250 300 35000.10.20.30.40.50.60.70.80.91 d = 10, q = 0.1d = 30, q = 0.1d = 70, q = 0.1d = 100, q = 0.1d = 130, q = 0.1
FIG. 8: Fig.8: (Color online) Profile of density at variousvalues of d for q = 0 .
1. The profile exhibits a linear behaviourwith positive slope. This behaviour is associated to the exis-tence of wandering shock in the region between left and rightdefects. extensive MC simulations that the density profile takes alinear shape between defects [10]. This behavior remainsunchanged in ASEP with extended objects [31]. For thesake of completeness, we show some typical density pro-files in Fig. 8.The interesting point is the absence of boundary layerin this phase-segregated regime. It would be illustrativeto look at the dependence of two-point correlation func-tions at some particular sites on values of q and d . Theseresults are sketched in Fig. 9 where correlations at thefirst defect site ( d ), its rightmost site ( d + 1) and in themiddle site of the chain are considered.Note that all correlations/anti-correlations approachzero when q tends to one. Moreover, increasing defect’sseparation d increases the correlations. The dependenceon d is more interesting. While the values of correlationfunctions reach an asymptotic value at large d, the behav-ior is not monotonous. Especially for C d +1 correlationincreases up to a maximum and then begin to decrease q C o rr e l a t o r s d= 10, C [d1]d= 100, C [d1]d = 10, C[d1+1]d = 100, C [d1+1]d =10, C [mid]d= 100, C [mid] d C o rr e l a t o r
10 20 30-0.5-0.4-0.3-0.2-0.100.10.20.30.40.50.60.7
C[d1],q =0.1C[d1+1],q =0.1C[mid],q =0.1C[d1],q =0.3C[d1+1],q =0.3C[mid],q =0.3
FIG. 9: Fig.9: (Color online) Dependence of normalized cor-relation functions at three selected sites on q (top) and on d (bottom). smoothly towards its asymptotic value. The value of d where C d +1 is maximized does not show a significantdependence on q . In can be concluded that varying d can dramatically affect the system characteristics as faras correlations are considered. B. Analytical theory
Our simulation findings in the preceding section con-firms that in between the defects the correlations are no-tably higher than other sites. In this section we try todevelop a theoretical framework to capture this feature.Suppose we have two slow defects located in the bulk atsites k and l ( k < l ) respectively both with rate q . Themean occupation at site i in the steady state is denotedby n i = h τ i i , in which τ i = 0 , i . We assume that a simple MF assumption, i.e., h τ i τ i +1 i = h τ i ih τ i +1 i = n i n i +1 holds for all sites except i = k, · · · , l i.e.; defective sites themselves and all thesites between them. At these sites the correlations arestrong enough to violate the simple mean-field assump-tion. Let us introduce two-point correlation functions m i in the following way, m i = h τ i τ i +1 i . (41)There are L + l − k + 2 unknowns, namely, n , n , · · · , n L , m k , m k +1 , · · · , m l and lastly the current J . In the stationary state there exists L + 1 equationsamong these unknowns. Let us label them by A to A L .These equations can be obtained by expressing the cur-rent J in terms of the mean site densities and two pointcorrelators. The first equation, A , is J = α (1 − n ). Theequations A i for i = 1 , · · · , k − i = l + 1 , · · · , L − J = n i (1 − n i +1 ) . (42)Equation A k is: J = q h τ k (1 − τ k +1 ) i = q ( n k − m k ) . (43)Similarly, equation A l is given below J = q h τ l (1 − τ l +1 ) i = q ( n l − m l ) . (44)Equations A i ( i = k + 1 , · · · , l −
1) have the followingform J = h τ i (1 − τ i +1 ) i = ( n i − m i +1 ) . (45)lastly equation A L is as follows, J = β h τ L i = βn L . (46)We do not intend to add more equations. Unfortu-nately the above equations are nonlinear and it wouldbe a formidable task to solve them analytically. Alterna-tively, we shall utilize a numerical approach to solve thesystem of equations by exploiting their recursive struc-ture. This approach was originally introduced in [29] inthe context of disordered ASEP and was later applied tothe problem of two intersecting ASEP chains [30]. Ac-cording to this numerical scheme, we assign a value to J .Having J , it is possible to iterate equations (in forwarddirection) and obtain n up to n k . Then we proceed tofind m k by incorporating the relation J = q ( n k − m k ).At this stage it is not possible to proceed further becauseboth n k +1 and m k +1 are unknown and we have only onerelation between them : J = n k +1 − m k +1 . In order toproceed, we approximate n k +1 in the following way. Con-sider the rate equation for the 2-point function h τ k τ k +1 i which is governed by the following master equation: d h τ k τ k +1 i dt = h τ k − (1 − τ k ) τ k +1 i − h τ k τ k +1 (1 − τ k +2 ) i . (47)In the steady state, the left hand side becomes zero,and therefore two terms on the right hand side will be equal. To proceed further we have to approximate 3-point functions. This is achieved by utilizing the clustermean-field assumption [40]. According to this assump-tion we replace any three point function by the productof 2-point functions as follows, h n i n j n k i = h n i n j ih n j n k ih n j i . (48)We then replace all 2-point functions by the productof 1-point functions except m k = h τ k τ k +1 i . Then it ispossible to express 1 − n k +2 in terms of n k − , m and n k +1 . After some algebra a quadratic equation for n k +1 is obtained: n k − n k +1 − m k n k − n k +1 − m k J = 0 . (49)The physically reasonable solution for this equation isgiven by n k +1 = m k + q m k + m k Jn k − . (50)Now it is possible to find m k +1 via equation J = n k +1 − m k +1 . Analogous to the above procedure we can find n k +2 as follows: n k +2 = qm k m k +1 + q ( qm k m k +1 ) + 4 qJn k n k +1 m k +1 qn k n k +1 . (51)After having n k +2 we simply obtain m k +2 via equation J = n k +2 − m k +2 . Now it would be possible to proceediteratively after taking into account some approximation.To this end we recall the equality d h τ i τ i +1 i dt = h τ i − (1 − τ i ) τ i +1 i − h τ i τ i +1 (1 − τ i +2 ) i . (52)Putting the left hand side equal to zero in the steadystate, utilizing cluster mean-field in three point functionsand finally substituting m i +1 by m i +1 = n i +1 − J wearrive at the following equation: n i +1 = m i − m i + p ( m i − m i ) + 4 Jn i − n i m i n i − n i . (53)Note that we have approximated h τ i − τ i +1 i by themean-field relation h τ i − ih τ i +1 i . We can iterate equation(53) together with m i +1 = n i +1 − J from i = k + 2 to l − n i and m i up to i = l − i = l needs to be treated separately. Followingthe same strategy we easily find: n l = m l − m l − + q ( m l − m l − ) + 4 Jn l − n l − m l − n l − n l − . (54)0 q MF, d = 2DMF, d = 2ISA, d = 2NR, d = 2MC, d = 2 J FIG. 10: Fig.10: (Color online) Current vs q for d = 2obtained by various analytical and numerical methods andMC simulation. J approaches to 0.25 when q tends to one inaccordance to normal ASEP. From which one can compute m l . In a similar fashion,we can obtain n l +1 . We only should shift up all the sub-scripts in equation (54) by one. Now it is possible againto proceed iteratively to the end of the chain and evaluate n L which gives us the output current J out . If the givenvalue of input J were correct, the output current J out ,which is βn L , should be the same, up to a given precision,as the input value of J . By systematically increasing theinput J in an small amount δJ , we can determine the cor-rect J and correspondingly the mean densities n , · · · , n L together with correlators m k , · · · , m l . In Fig. (10) thedependence J on q obtained by the numeric scheme de-vised above is sketched. For the sake of comparison, wehave augmented the figure with the analogous graphs ob-tained by MF, DMF, ISA and MC methods.The result of the numerical scheme is almost identicalto ISA method. They both are in very good agreementwith Monte Carlo simulations. However, the numericalscheme has an advantage over ISA method in the sensethat it can easily be implemented for any d , whereas find-ing the solution of the ISA nonlinear set of equations,i.e., Eq. (29) is not an easy task even by employingadvanced numerical methods. Finally in figure (11) wehave sketched the dependence of J on d obtained fromMC simulation and the numeric scheme.The results of the numeric method are in rather goodagreement with those obtained by MC simulations. J isan increasing function of d and becomes saturated aftersome short q -dependent distance. The results confirmsthe earlier finding in [10]. Note that the length scale onwhich J recovers its single-defect value is of the same or-der of magnitude of the correlation length in the densityprofile near boundaries. Finally we would like to addthat our numerical scheme is not capable of reproducingthe profile of correlators obtained via MC simulations.Figure (12) depicts the profile of unnormalised adjacenttwo-point correlation function h τ i +1 τ i i−h τ i +1 ih τ i i for twomethods of simulation and numerical scheme. d
10 20 300.060.080.10.120.140.160.180.20.220.24 q = 0.1, NRq = 0.1, MCq = 0.3, NRq = 0.3, MCq = 0.5, NRq = 0.5, MC J FIG. 11: Fig.11: (Color online) Current vs d for q = 0 . , . . site number
150 200 250 300 350-0.07-0.06-0.05-0.04-0.03-0.02-0.0100.010.020.030.040.050.06 q= 0.3 , simulationq = 0.3 , numeric C i FIG. 12: Fig.12: (Color online) Profile of corrrelators for d = 100 and q = 0 . We see that within the numerical framework both thevalue and the extension of correlators are small in com-parison to the simulation results. The reason lies in theimplementation of a series of approximation in severalplaces in this numerical algorithm.
IV. SUMMARY AND CONCLUSIONS
An open ASEP chain with two defective sites with re-duced hopping rates q < q less than 0 . q approachesone. However, the dependence of two point correlators for a fixed q as a function of distance between defectsis not monotonous. Despite reaching to an asymptoticvalue for large distances, one observes a peak at shortdistances for the two point correlators near defects. Ourtheoretical analysis indicates that correlations are criti-cally important for dynamics of particles in disorderedASEPs. It also shows that it is possible to devise an ap-proximate method that can take into account these corre-lations, providing a satisfactory description of stationaryproperties. Acknowledgments
ABK acknowledges the support from the Welch Foun-dation (under Grant No. C-1559), and from the USNational Science Foundation (grants CHE-0237105 andNIRT ECCS-0708765). MEF expresses his gratitude toM. F. miri for their useful help. [1] B. Derrida, Phys. Rep.
65 (1998).[2] B. Schmittmann and R.K.P. Zia in:
Phase transitionsand Crtitical Phenomena , vol 17 , ed. C. Domb and L.Lebowitz (London: Academic) 1995.[3] G. M. Sch¨utz, Integrable stochastic many-body systems in Phase Transitions and Critical Phenomena edited by C.Domb and J. Lebowitz (Academic, London, 2000), Vol.19.[4] D. Chowdhury, L. Santen and A. Schadschneider, Phys.Rep.
199 (2000).[5] B. Derrida, M. R. Evans, V. Hakim and V. Pasquier, J.Phys. A: Math. Gen.
21 (2007).[11] R. Lipowsky, S. Klumpp and T. M. Nienwenhuizen, Phys.Rev. Lett. R333 (2007).[14] A.B.Kolomeisky, G.M.Sch¨utz, E.B.Kolomeisky, andJ.P.Straley, J. Phys. A: Math. Gen. , 471 (1993).[17] G.M. Sch¨utz, Phys. Rev. E , 4265 (1993).[18] S. A. Janowsky and J. L. Lebowitz, J. Stat. Phys.
83, 2004.[26] R. Juhasz, L. Santen and F. Igl´i, Phys. Rev. Lett.
22, 2006.[28] P.Pierobon, M. Mobilia, R. Kouyos and E. Frey, Phys.Rev. E,
813 (1993).[35] K. Mallick, J. Phys. A: Math. Gen , 311(1993).[40] D. Chowdhury and J.S. Wang, Phys. Rev. E65