Current fluctuations in a partially asymmetric simple exclusion process with a defect particle

We study an exclusion process on a ring comprising a free defect particle in a bath of normal particles. The model is one of the few integrable cases in which the bath particles are partially asymmetric. The presence of the free defect creates localized or shock phases according to parameter values. We use a functional approach to Bethe equations resulting from a nested Bethe ansatz to calculate exactly the mean currents and diffusion constants. The results agree very well with Monte-Carlo simulations and reveal the main modes of fluctuation in the different phases of the steady state.


I. INTRODUCTION
Models of interacting driven diffusive particles are a natural effective description for many systems found in nature, particularly in biology.Cases to which such models have been applied include the motion of RNA polymerase during DNA translation [1] and ribosome dynamics in mRNA translation [2], traffic flow on a busy street [3,4], and driven colloids in a narrow channel [5][6][7].Moreover, these models have been shown to have links to many other problems in statistical physics, including disordered polymers in random media [8], surface growth models [9] (notably, some of the models are known to lie within the KPZ universality class [10]), diffusion in strongly anisotropic materials [11], equations in fluid dynamics, such as the Burgers equation [12], and certain combinatorial problems [13].
Both in the mathematical and physics literature, there has been focus on one-dimensional systems, for which there are powerful exact methods.In particular, the asymmetric simple exclusion process (ASEP) has become the prototype of driven diffusive systems.The simplicity of the ASEP has allowed for many exact results for its stochastic dynamics to be derived.See, for instance [14][15][16][17][18] for reviews.
A good understanding of fluctuations in such minimal models is important for several reasons.As the systems to which these models are applied to, such as traffic on a highway, typically contain many fewer particles than conventional equilibrium systems, fluctuations may be important to account for finite-size effects.
In particular, analyzing the system at a microscopic level allows one to derive the properties of fluctuations without postulating their form, as would need to be done if starting from a hydrodynamic picture.This type of microscopic analysis is especially relevant as these models are far from equilibrium, which means that standard tools for describing fluctuations, such as Onsager-Machlup theory, are not applicable.Conversely, one can use exact results from microscopic dynamics to see whether general results can be derived for nonequilibrium systems from first principles.
The two main approaches that have emerged in the literature for calculating current fluctuations in ASEPs exactly are matrix product states and the Bethe ansatz.The matrix product approach was initially used to describe the spatial structure and mean current in the steady state [19].It has been extended to fluctuations in some cases [20,21], but this generalization has proved to be quite difficult.
On the other hand, the Bethe ansatz approach allows for direct calculations of the full current statistics.To be precise, it is used to calculate the scaled cumulant generating function of particle displacement, which can be done to all orders in some simple cases [22,23].
In this paper, we investigate a partially asymmetric simple exclusion process (PASEP) with a defect particle that has priority in the dynamics.Recently, the steady state of this model was solved using a matrix product ansatz and the mean current at long times was calculated [24].Moreover, it has also recently been shown that it can be solved using a coordinate Bethe ansatz [25].Building on these results, we use a functional Bethe ansatz to rederive the mean current and calculate exactly the current fluctuations around the mean, which can be related to the diffusion constant.We emphasize that this is the first time the Bethe ansatz has been used to calculate current statistics in a partially asymmetric process with a phase transition.As such, the calculations presented here involve a combination of the techniques developed for the homogeneous partially asymmetric case and totally asymmetric case with a defect.
The importance of the Bethe ansatz to the ASEP was first appreciated when it was realized that its transi-tion rate matrix (Markov operator) has a very similar structure to that of a quantum spin chain Hamiltonian [26,27].Indeed, both of them are naturally expressed as sums of tensor products of Pauli matrices.Consequently, the well-known Bethe ansatz techniques from the quantum spin chain case could be carried over to study the spectral properties of the asymmetric exclusion process [27][28][29][30].Furthermore, by using a nested Bethe ansatz, Bethe equations have been derived for a PASEP with particles of different sizes [31] and multi-species hierarchical PASEPs [32,33].However, to our knowledge, the calculation of current statistics beyond the mean in these latter models remains an open problem.
A major technical advancement was a modification of the time evolution problem that allowed current statistics to be calculated easily [22].It was shown that a deformation of the transition rate matrix that corresponds to conditioning on a large current gives the time evolution of the total particle displacement.As this satisfies a large deviation principle in the long time limit, the cumulant generating function of the current in the steady state could thus be calculated using the Bethe ansatz.
This method was used to calculate the large deviation functions of the current in the TASEP [22], the TASEP with a defect particle [34] and the PASEP [35].For the PASEP, it proved to be useful to reformulate the problem as a functional Bethe ansatz [23,35,36].This simplified calculations to the extent that allowed the cumulants to all orders to be formally expressed in terms of combinatorial objects [23].The case considered in this paper is an extension of these results to the PASEP with a defect particle, which requires a generalization of the methods for those earlier cases.
The remainder of this paper is structured as follows.In Sec.II, we define the model, review the known results obtained using matrix product states and the coordinate Bethe ansatz for this model, and state the new results for the diffusion constant.In Sec.III, we reformulate the Bethe equations in a functional form, which allows the cumulants of the current to be calculated directly.We perform these calculations in Sec.IV and Sec.V.In Sec.IV, we show that the Bethe ansatz solution reproduces biased diffusion statistics for the defect particle, as expected.In Sec.V, we calculate the cumulants of the hopping of normal particles to second order.In Sec.VI, we make some final remarks.We fill in the algebraic detail of some of the lengthier calculations in the appendices.

A. Model definition
We consider a ring with L+1 sites, M normal particles and one defect particle.The normal particles hop to the right and left with rates p, q respectively, and the defect particle hops to the right and left with rates αp, q/α respectively.The defect also overtakes normal particles to its right and left with rates αp, q/α.
We can summarize the dynamics as follows, where we denote empty sites, the defect and normal particles as 0, 1, 2 respectively: (1) Note that the defect (denoted by 1) does not distinguish between normal particles and empty sites, and can therefore be seen as having priority in the dynamics.Because of this, one may think of the defect as a "first-class" particle, whereas the normal particles may be thought of as "second-class".
It is convenient to introduce the normal particle density and asymmetry parameters, We will take x < 1 in this paper (but α can have any positive value).

B. Large deviation theory
In investigating the long time current statistics, the central objects are the random variables Y 1 (t), Y 2 (t), Y 12 (t), which count the number of processes of type 10 → 01, 20 → 02, 12 → 21 respectively, minus the reverse processes, up to time t.
In the long time limit, t → ∞, the joint moment generating function of these variables satisfies a large deviation principle, ⟨e γ1Y1(t)+γ2Y2(t)+γ12Y12(t) ⟩ ∼ e λ(γ1,γ2,γ12)t , where γ i are the conjugate variables of Y i (t) and λ is the rate function.
The rate function λ thus directly gives the long time limit of the scaled cumulants of the hop-counting variables, Using this general formulation, one could in principle compute the statistics of any combination of Y 1 , Y 2 , Y 12 .In practice, this makes the calculations very cumbersome.Therefore, in this work, we shall focus on two important cases: (i) the net displacement of the defect, which corresponds to γ 2 = 0, γ 1 = γ 12 ̸ = 0; and (ii) the hopping statistics of the normal particles, which corresponds to γ 1 = γ 12 = 0, γ 2 ̸ = 0 [37].Thus we will use the notation where γ is the single formal parameter conjugate to the relevant variable, and a i are constants, which we can set to a 2 = 0, a 1 = a 12 = 1 to track the defect, or a 1 = a 12 = 0, a 2 = 1 to track the normal particles.
Then, for instance, the long time limit of mean current and diffusion constant of normal particles is obtained by setting a 1 = a 12 = 0, a 2 = 1 and evaluating C. Review of previous results

Phase diagram
In [24], the steady state of this model was solved using a matrix product ansatz.It was shown, that in the limit L → ∞ with ρ fixed, the system exhibits three phases, which have distinct expressions for the density profiles and current.The phase diagram consists of two localized phases (L), in which the defect only has local effects on the normal particles; and a shock phase (S), in which the defect disrupts the normal particle current, creating two macroscopic regions with different bulk densities ρ 1 , ρ 2 .These do not depend on the mean density, ρ, but are instead given purely in terms of the system parameters, α and x, The steady-state density profiles in the reference frame of the defect are shown schematically in Figure 1.The phases are delimited by the transition lines with the shock phase lying in the region ρ 2 < ρ < ρ 1 and the localized phases outside it.The phase diagram is shown in Figure 2.

Current
The current of normal particles in these phases was also calculated exactly in [24] using the matrix product ansatz.In the limit L → ∞, with ρ held fixed, it is given to leading order in L by In the localized phases, the current is essentially the same as in a homogeneous PASEP.In the shock phase, the defect throttles the current.Observe that the currents in the shock phase and localized phases are equal at the points of phase transition (i.e., there is no discontinuity).However, the current in the shock phase depends linearly on density, whereas in the localized phases, it is a concave function of density.Thus the defect evidently makes the current smaller than that of a pure system.
The current in the shock phase can be further understood as follows.Let k be the mean position of the shock front in the reference frame of the defect.Defining u = k/L, total particle number conservation means, Substituting (10) in (9b), we can write the current in the shock phase as This suggests the following interpretation.The current of normal particles is controlled by the defect.The defect sees a normal particle density ρ 1 behind itself and ρ 2 in front of itself.As the defect hops forward, with rate αp, it creates a small region of density ρ 2 behind itself.This "hole" has to travel backwards L(1 − u) sites to restore the steady-state profile.Thus the defect hopping forward generates a net current of normal particles in the forward direction of magnitude αpL(1 Similarly, as the defect hops backwards, with rate px/α, it creates a small region of density ρ 1 in front of itself.This has to travel forwards Lu sites to restore the steady-state profile, giving a net current of (q/α)Lu(ρ 1 − ρ 2 ).
These two effects combine to give (11).

D. Diffusion constant
In the present work, we use a functional Bethe ansatz approach to calculate exactly the diffusion constant of normal particles, as defined in (6b).The method to obtain an exact expression is outlined in Sec.V D.
Asymptotics to leading order in inverse system size are also extracted in Sec.V D, giving the results, The comparison of both exact finite size and asymptotic results with Monte Carlo simulations is given in Figure 3.The estimates from Monte Carlo simulations are in excellent agreement with the exact finite-size results.The agreement between finite-size and asymptotic expressions is good deep in each phase, but near the phase transitions there is some discrepancy due to strong finite-size effects.On increasing the system size, the convergence of the finite-size results to the asymptotic values can be seen, albeit slowly near the phase transitions.This is shown in Figure 4.
This result reveals the main sources of fluctuations in the two phases.In the localized phases, the diffusion constant is the same as in a homogeneous PASEP with bulk density ρ [21].This intuitively makes sense, as in the localized phases, the defect does not have a macroscopic effect on the system.
In the shock steady state, we can write the diffusion constant as, This expression can be interpreted through the same picture as the current in this phase.The fluctuations are controlled by the defect, which creates low density "holes" at rate αp, that travel backwards, and high density waves at rate px/α, that travel forwards.This creates fluctuations, of magnitude L(ρ 1 − ρ 2 )(1 − u) and L(ρ 1 − ρ 2 )u respectively, whose squares have to be added weighted by their rate of creation to obtain the total variance (13).

III. FUNCTIONAL BETHE ANSATZ CALCULATION OF THE LARGE DEVIATIONS OF THE CURRENT A. Bethe ansatz solution
In [25], it was shown that the rate function, as defined in (3), can be calculated for the model (1) using a coordinate Bethe ansatz.In summary, this is done by considering the transition rate matrix encoding the dynamics (1) and applying to it a deformation which counts the processes 10 → 01, 20 → 02 and 12 → 21.The rate function (3) can be identified with the eigenvalue of the deformed transition rate matrix with the largest real part.This eigenvalue is then found by using an eigenvector of the Bethe ansatz form, where x 0 denotes the position of the defect and x 1 , . . ., x M the positions of the normal particles; σ denotes permutations of {0, 1, . . ., M }; A σ are amplitudes; and z i are complex numbers called Bethe roots.Using this ansatz, the rate function can be written in terms of the Bethe roots, The Bethe roots are to be found by solving the Bethe equations, for i = 0, . . ., M , where and B is a constant.We remark that the structure of these equations is similar to the homogeneous PASEP case [35], but with an additional constant, B, which is also a feature in the case of a TASEP with a defect [34].
In addition to the Bethe equations, we need to enforce the conditions, The condition (18) comes from the translational invariance of the steady state.This means that applying the translation x i → x i + 1 for i = 0, 1, . . ., M to the eigenvector ( 14) should not change it, which is equivalent to (18).We shall henceforth refer to this as the periodicity condition.
The condition (19) comes from the fact that λ is the eigenvalue of the deformed transition rate matrix with the largest real part.Hence, in the undeformed limit, γ → 0, it should converge to the 0 eigenvalue of the transition rate matrix.
Finally, the constant B can be fixed by multiplying ( 16) for all i.Then together with (18), we get

Ground state solution
It is instructive at this point to consider the ground state, γ → 0. In this case, the Bethe roots, z i , all converge to 1, except one root (we may choose it to be z 0 without loss of generality), which converges to a different finite value.From (25) and the condition (19), it is evident that We note that in the TASEP with a defect particle whose hopping rate in units of the normal particles' hopping rate is α, the phenomenon of all Bethe roots converging to the 1, except for one, which converges to α −1 has also been observed [34].

Change of variable
The equations take a more felicitous form if we consider the following transformation of the Bethe roots, Such transformations have been used before in Bethe ansatz solutions for the homogeneous PASEP [23,35] and the asymmetric XXZ spin chain [29].We will refer to y i as the Bethe roots from here on.The Bethe equations ( 16) become, the B-fixing equation ( 20) becomes, and the periodicity condition (18) becomes, Furthermore, the rate function (15) can be expressed as The ground state corresponds to y i → 0 for i = 1, . . ., M and where the first equation comes from ( 21) and ( 22) and the latter equation comes from plugging (27a) into (24).The structure of the Bethe equations ( 23) is similar to the homogeneous PASEP case, which has been studied previously [23,35].The approach that has proved to be most fruitful is to transform the equations into a set of functional equations, which can be solved directly to extract the behavior of functions of the roots, like λ, without explicitly calculating the roots.
Inspired by this approach, we now reformulate (23) as a functional equation for a single-variable function.Then, we further manipulate this equation in a procedure known as going "beyond the equator", which ultimately allows direct calculation of λ.

B. One-variable function formulation
We first define the single-variable polynomial, (28) where Then from (23), it follows that all Bethe roots y i are roots of S(T ) (i.e., S(y i ) = 0).However, the degree of S(T ) is L + M + 3, whereas there are only M + 1 Bethe roots.This suggests that we should also consider the polynomial Then as the Bethe roots y i are roots of both S(T ) and Q(T ), and the degree of S(T ) is higher, we conclude that S(T ) must be divisible by Q(T ) as a polynomial.This means that there exists some polynomial R(T ) of degree L + 2, for which, This is a functional equation for Q.As Q contains the same information as the Bethe roots y i , in this functional formulation, this equation plays the same role as the Bethe equation (23).
Using Q also allows us to rewrite the B-fixing equation (24) as and the periodicity condition (25) as Moreover, the rate function (26), can be written in terms of Q as Now equations (31)(32)(33) are self-consistent and can be used to solve for Q(T ) (and therefore λ) without reference to the Bethe roots y i .

C. Bethe equation beyond the equator
We now proceed to change (31) into the so-called "beyond the equator" form [23,38], which ultimately allows us to write an equation involving only one unknown function.

Bethe equation
The Bethe equation beyond the equator is given by 2Ch(T ) = P (T /x)Q(T ) − x M ξ −1 P (T )Q(T /x), (35) where P is a polynomial of degree L − M + 1 and C is a constant defined as The derivation of ( 35) is given in Appendix A.

Ground state solution
We now consider again the ground state case, γ = 0. Recalling the remark about the Bethe roots in Sec.III A, we can write down the following form of the ground state solution, The form of Q (0) (37a) follows from the location of the Bethe roots.Using this result in (36), we obtain (37c).
Putting this back into (35), the form of P (0) (37b) can be inferred.It can be verified that with (27a-27b), this solution satisfies equations ( 32), ( 33), (35) and the condition that (34) vanishes.Crucially, we note that (37c) implies that C = O(γ).This will be important for the perturbation theory, as it turns out that it is more convenient to build an expansion in C, rather than γ.
We remark that P (0) is a polynomial of order 1, even though P is of order L − M + 1.This is because the coefficients of higher powers of T are of order γ in the perturbative expansion.Indeed, the ground state is a special case because C (0) = 0, so the degrees of P (0) , Q (0) are not strictly fixed.However, for the general case, C ̸ = 0, so we must have deg P + deg Q = deg h.

D. One function reformulation of Bethe equation
At this point, the Bethe equation ( 35) still involves two unknown functions, P and Q.We follow the approach from [23] to formulate the problem in terms of a single function, w.The Bethe equation then becomes a self-consistent equation for w, which can be solved to calculate w order by order in C.

Bethe equation
After some algebraic manipulations, which are outlined in Appendix B, (35) can be written as an equation in terms of only one function, w, which reads, Here, where y 0 is the (yet undetermined) exact value of the Bethe root that does not converge to 0, and X is an operator on power series in T , with the convention for an arbitrary constant µ.To simplify some calculations, in this paper we make the choice µ = −1, though this does not affect any physical results.We remark that (38) has the same form as the functional Bethe equation derived for the pure PASEP [23].However, there the analogue of h had the simpler form (1 − T ) L /T M .In the present case, h has additional poles at y 0 , xy 0 and contains two undetermined constants: y 0 and B. As we shall see, the additional poles imply the existence of phase transitions.The equations to fix the constants y 0 , B are given in the next section, (44), (45a), leading to a closed system of equations for w(T ), B, y 0 , C.

Additional equations
We now rewrite the remaining equations ( 32), ( 33), (34) in terms of the function w.In doing this, it is helpful to define another operator P that acts on a power series u(T ) as follows, where sgn is the signum function, which we have defined with Thus P reverses the sign of the terms with non-negative powers of T .Then, after some algebra, which is given in Appendix C, we can express the B-fixing equation (32) as, and the periodicity condition (33) as, where we have introduced the shorthand To complete these equations, we need to introduce an additional equation to fix the constant C. As C is a global constant that multiplies h, this condition can be interpreted as a normalization condition.This is given by, where the notation {T k } indicates that we take the coefficient of T k in the power series expansion of the expression that follows.The derivation of this is also given in Appendix C. Now (38) and (44-46) form a complete set of equations, which can be used to solve for w(T ), B, y 0 and C. In practice, other than for some exceptionally simple cases, one needs to proceed perturbatively.
Lastly, to complete the formulation in terms of w(T ), we can use the operator P to express the rate function as

IV. HOPPING STATISTICS OF DEFECT PARTICLE
We now proceed to explicitly calculate the cumulants of the currents at long times.We first consider the net displacement of the defect (a 1 = a 12 = 1, a 2 = 0).As the defect particle simply performs biased diffusion, we should recover Skellam statistics (i.e., the difference of two Poisson random variables).This result is obvious even without using the Bethe ansatz, but it serves as a verification of the validity of the calculation.
Setting a 1 = a 12 = 1, a 2 = 0, we immediately get from (46), Applying this to (38) and noting that {T 0 } h(T ) does not vanish, it is evident that this is consistent only if C = 0, which means w(T ) = 0 at all orders.
Then from (45b) we obtain using which, (47) yields the exact result This is precisely the cumulant generating function of a Skellam distribution with parameters α, x/α, as expected.The cumulants are given by, V. HOPPING STATISTICS OF NORMAL PARTICLES Now we examine the hopping statistics of normal particles.This is obtained with the choice a 1 = a 12 = 0, a 2 = 1.We will see the different behavior of the localized and shock phases reflected in the asymptotic limit L → ∞ with ρ = M/L held fixed.
Unlike the simpler case of Sec.IV, here w(T ) ̸ = 0, and we need to make use of the Bethe equation (38).We will do this perturbatively.

A. Perturbative expansion
To solve the system (38) together with (44-46), it turns out to be convenient to expand everything (including γ) in powers of C and then eliminate C order by order.This kind of parametric expansion has been used in Bethe ansatz solutions for other exclusion processes [22,23,34].From (37c) we recall that C = O(γ), so this expansion is justified.We use the notation Z (k) to denote the k-th order term in the expansion of some variable Z in powers of C.
To obtain derivatives of λ with respect to γ, we expand both in powers of C, Inverting (52b) gives to second order in γ, Substituting this into (52a), we can express the derivatives of λ with respect to γ in terms of the coefficients in the C expansion, λ (k) , γ (k) .For the first two derivatives, we get Our task now reduces to expanding all quantities in powers of C. For example, expanding (38) yields which we match with the expansion w = w (0) + Cw (1) + C 2 w (2) + . ... This yields Note that h(0) is known by using (39), where y (0) 0 is given by (27a) and B (0) is given by (27b).The calculation of w (2) requires knowledge of h(1) and therefore of y (1) 0 and B (1) .From equations ( 44) and (45a), we see that we need the expansions of γ and δ with respect to C.
Finally, using (47) we get the expansion of the rate function λ to second order, +(α − x/α)δ (1) .(63b) This will allow us to calculate the mean current and diffusion constant of the normal particles.

B. Integral expressions
In order to evaluate expressions explicitly, we rewrite some key quantities in terms of complex contour integrals.This is also very helpful for extracting asymptotics.In the following, Γ is a small circle around the origin and a factor dT /(2πi) is implied, Equations ( 64a) and (64b) are simple applications of the residue theorem.Equations (64d), (64e) and (64f) are derived in Appendix D.
Using these equations, and the integral expressions (64), we obtain where The normalisation Z L,M can be computed explicitly as follows.We first use the definition of h(0) , (58).Then, we expand the numerator using the binomial theorem.Crucially, we observe that the terms in the expansion with powers larger than M do not contribute to the integral as they have no pole at T = 0.The remaining terms can be evaluated using the residues of the poles outside the contour (at T = y (0) 0 and T = xy (0) 0 ).Note that due to the truncation of the binomial expansion at T M , we do not have to worry about poles at infinity.After some simplification, this gives, where we have used the definition of B (0) , (27b).The numerator in F can be computed similarly.Putting this into (66), we obtain an exact expression for the current in finite-size systems.This gives an alternative expression for the mean current, which was previously calculated in [24] using a matrix product ansatz.
There, the current was defined as the flux of normal particles across a given bond, whereas in this paper it was defined as the motion of all particles everywhere in the system, (6a).Hence, to compare the two expressions, we have multiplied the former by L. The expression from the matrix product approach, J MPA , and the one obtained in the present work, J BA , read, 0 is defined by (27a).These expressions can be checked to agree numerically.In particular, one can use a symbolic programming language like Mathematica with rational numbers or integers for all parameters, as this gives exact values, without any machine error.This has been performed for various system sizes and parameter choices and the two results have been found to agree, though a rigorous, analytic proof of identical equality is lacking.

Asymptotic expressions for J and phase diagram
The result given so far, (72), is exact and can be evaluated numerically for finite system sizes.To make sense of it physically, it is beneficial to extract the asymptotic behavior in the limit L → ∞, with ρ held fixed.
The key quantity is Z L,M , as defined in (69).We can write it as where This integral has a saddle point at the solution of the equation ϕ ′ (T ) = 0.This is found to be We can always deform the contour of integration to pass through the saddle point.In doing this, we may need to pass through the poles at T = y (0) 0 and T = xy 0 , in which case the contributions of their residues must be subtracted.
The original contour is a small circle around T = 0. Hence, the poles must be subtracted if 0 > y (0) 0 > T 0 and 0 > xy (0) 0 > T 0 respectively.Rearranging these inequalities, we see that we have 3 cases.
No poles.When α > 1 or ρ < ρ 2 , both poles are outside the contour, so the integral is dominated by the saddle point.
Two poles.When α < 1 and ρ > ρ 1 , both poles are inside the contour.However, it can be verified that the contributions from their residues cancel exactly.Therefore the integral is still dominated by the saddle point.
One pole.When α < 1 and ρ 1 > ρ > ρ 2 , the pole at xy (0) 0 is inside the contour but the pole at y (0) 0 is outside the contour.In this case the integral is dominated by the residue from the pole at xy (0) 0 .In summary, for ρ 2 < ρ < ρ 1 (i.e., in the shock phase), the integral is dominated by the pole at T = xy (0) 0 , and otherwise (in the localized phases) it is dominated by the saddle point.This implies the phase diagram presented in Section II C. Using these results, we get to leading order where f is defined in (68).Plugging this into (66), we get expressions that agree with (9a-9b).

D. Diffusion constant
For the diffusion constant, given by (55) we require (59c) and (63c).Using these expressions and the integral expressions (64a), (64e), and (64f), we can write the diffusion constant as follows, where and Z L,M is defined in (69).We note that w (2) appears in (78).The function w (2) , as determined by (57c), contains h(1) , which is obtained by expanding (39) to first order in C, We see that this expression contains the constants δ (1) , B (1) .Using (85a), (64a), and (64e), δ (1) can be expressed as and B (1) is given implicitly by (62).Hence, the explicit finite-size expression for ∆ is much more complicated than that for J and we do not give it here.Instead, we proceed immediately to the asymptotic behavior.

Asymptotic expressions for ∆
The asymptotic analysis of ∆ requires us to extract the large L behavior of h(1) , X[ h(0) ] (which enter via w (2) (57c)), and G.
The asymptotic analysis is much less straightforward than it is in the calculations of the current.However, guided by known results in related systems (for instance [20,21,34,39]) as well as Monte Carlo simulations, one expects that ∆ has the scaling L 2 in the shock phase and L 3/2 in the localized phases.This helps us to determine which terms will contribute to the asymptotic behavior.
The integrals still have no-, one-and two-pole regimes.However, if one looks at certain terms (such as G) in isolation, the residues from the two poles do not cancel exactly and some super-dominant scaling seems to emerge.These super-dominant terms ultimately cancel out when all the terms comprising ∆ are considered together.Although we do not have a rigorous analytic proof of this cancellation mechanism, our results are supported by extensive numerical analysis, as well as agreement with the exact finite-size results.
Asymptotics of G. Recall that the definition of G (79) contains g (80).Importantly, we have g(T 0 ) = 0 which implies that G does not contribute at leading order in the localized phases, when the integrals are dominated by the saddle point.In the shock phase, G is dominated by the residue of the pole at T = xy (0) 0 , therefore we have Asymptotics of X[ h(0) (T )].Splitting X[ h(0) (T )] in two parts, like (64b), we eventually see that the integral involving the kernel K never contributes at leading order.Hence, to leading order X[ h(0) ](T ) ≈ h(0) (T ). (84) Finally, h(1) , as given in (81), contains a term proportional to B (1) .However, we have verified that this term does not contribute in any phase.
Thus, in the localized phase, we can write w (2) (T ) ≈ −( h(0) (T )) 2 , whereas in the shock phase we have, w (2) (T ) ≈ h(1) (T ) − ( h(0) (T )) 2 .Ultimately, the contributing terms to the diffusion constant in each phase are a. Localized phase.In the localized phase, the integrals are to be evaluated at the saddle point.In the numerator, the first correction to the saddle point has to be computed to get the leading order contribution.This calculation, although quite cumbersome, is not as difficult as might first appear, as many of the terms present in the general saddle point correction formula cancel out due to the similarity of the two terms being subtracted.The details of this calculation are given in Appendix E. Ultimately, we obtain with f as defined in (68) and ϕ as defined in (75).Simplifying this, we obtain the expression in (12a).b.Shock phase.The integrals are dominated by the residues of the pole at T = xy (0) 0 .Recall that in this phase, F ≈ f (xy (0) 0 ).For the first term inside the integral in (85b), we have The second term inside the integral in (85b) can be evaluated after some manipulation as Combining these results and (83), we get, Simplifying this, we obtain the expression in (12b).

VI. CONCLUSION
We have used the Bethe ansatz to calculate the first two scaled cumulants of the particle hopping count in a PASEP with a defect particle that has priority in the dynamics.The first scaled cumulant (mean current) (9a-9b) agrees with the known result obtained using a matrix product ansatz [24].The second scaled cumulant (diffusion constant) (12a-12b) is a novel result, which is shown to agree with Monte Carlo simulations (see Figure 3).
The asymptotics of the scaled cumulants in the limit L → ∞ with ρ held fixed were also calculated using asymptotic analysis of the underlying integral expressions.The phase transitions were shown to correspond to a transition of the integrals being dominated by a saddle point (in the localized phases) and a pole (in the shock phases).
The asymptotic results indicate that in the localized phases, the system essentially does not feel the defect, with the current statistics being the same as those in a pure system.This makes sense intuitively, as the effects of the defect are localized, so its presence is not expected to be manifested at a macroscopic level.
In the shock phase, we have argued that the current is controlled by the defect, which creates density waves in the shock profile.These are small high density packets in the low density region and low density packets in the high density region.The results derived using the Bethe ansatz are consistent with this picture and it would be of interest to investigate whether this holds for higher order statistics.
An interesting aspect of the result for the second order cumulants is that the scaling is different in the two phases.In the localized phase, we have ∆ ∼ L 3/2 , whereas in the shock phase, ∆ ∼ L 2 .Hence there is a jump discontinuity in ∆ in the thermodynamic limit at the phase transitions, whereas the current, J, is continuous, with its first derivative being discontinuous at the phase transitions.We remark that a similar feature occurs in the open boundary TASEP, where a discontinuity of ∆ appears when crossing the shock line in the phase diagram [40].
In the shock phase the current fluctuations, given by ∆ 1/2 ∼ L, remain comparable to the mean current.This is consistent with the picture we presented in Sec.II C of the motion of the defect causing density fluctuations that must travel around the system and thus create current fluctuations of order L. As the motion of the defect is independent of the system size, the relative fluctuations do not decrease in the limit of large L. Now (34) can be written as Using (45a) and (45b), the last term can be simplified to Substituting this and (C10), we recover (47).
Appendix D: Proof of integral formulae Here, we prove the integral formulae presented in Sec.V B, specifically (64d), (64e) and (64f).We first make a comment on the operator P. The output of the operator given by the definition (42) is well-defined everywhere if the argument of the operator has a finite Laurent series expansion.For rational functions, such as h(0) , we can modify the definition slightly to avoid issues with convergence.Note that from (B9a), we see that at each order in perturbation theory w(T ) is a rational function, as it is given by a finite sum of p(T ) and q(T ), which are both rational functions.
Consider an arbitrary rational function a(T ) with a pole of order n at T = 0 and a finite number of other poles at some other locations T = a i , for i = 1, . . ., m.Using the partial fraction decomposition for rational functions, we can write a(T ) as, Then (64d) and (64e) are direct applications of this with T * = Bx and T * = 1 respectively.The required identities, w(Bx) = 0 and w(1) = 0, can be verified using (35) and the definition (B7a).
The identity (64f) can be proved by a straightforward extension of this argument.Then we also need the identity w ′ (1) = 0, which can also be verified using (35) and (B7a).
Strictly speaking, our proof only shows that the integral formulae hold up to any order in perturbation theory, as for any finite k, w (k) is always a rational function with a finite number of poles.However, this is sufficient for our purposes.where we have suppressed the arguments of A(T ), ϕ(T ), f (T ) for compactness and these functions are defined in (74), (75), (68) respectively.The terms in the numerator evidently cancel each other at first order in the saddle point.Hence, we use the general formula for the first correction to the saddle point (see any standard textbook, for instance [41, chapter 6]).For instance, we have (remembering the suppressed factor of 2π), . (E2) The overall factor of −1 arises because the saddle point T 0 is on the negative real axis and the original contour around the origin is anticlockwise, therefore the steepest descent contour goes from +i∞ + T 0 to −i∞ + T 0 .Applying this, we get

FIG. 1 .
FIG.1.Schematic representation of steady-state density profiles in the reference frame of the defect.Position is given as y = k/L.The reference frame is defined such that the defect is always located at k = 0 and the normal particles diffuse on sites k = 1, . . ., L. The two localized phases, which have structure only near the defect, are labelled L L/R .The shock phase is labelled S.

FIG. 2 .
FIG. 2. Phase diagram for x = 0.01 (top left), x = 0.2 (top right), x = 0.8 (bottom left) and x = 0.99 (bottom right).The two localized phases are labelled L L/R and the shock phase is labelled S.

( 1 −
x) 2 y 0 (1 − y 0 )(1 − xy 0 ) = αe δ − 1 − x + (x/α)e −δ .(C12) a(T ) = b 1 (T ) + b 2 (T ) m i=1 (T − a i ) + b 3 (T ) T n ,(D1)where b 1 , b 2 , b 3 are (finite) polynomials in T , with deg b 3 < n.Note that the first two terms are both analytic at T = 0. Then we define the action of the operator P on a(T ) as,P[a](T ) = −b 1 (T ) − b 2 (T ) m i=1 (T − a i ) + b 3 (T ) T n .(D2)It is easy to check that this alternative definition agrees with (42) within the latter's radius of convergence.Now we wish to evaluate a(T ) at some location T * , where a(T * ) = 0.This condition, together with (D1) and (D2), implies thatP[a](T * ) = 2 b 3 (T * ) T n * .(D3)At the same time, consider an integral of the type used in Sec.V B, namely, Γ a(T )/(T * − T ), where Γ is a small circle around the origin.Out of the terms in the expansion (D1), only the last one has a pole at T = 0, so the other two vanish.Then the integral of the last term can be evaluated using the sum of the residues of the poles outside the contour.As b 3 (T ) is a finite polynomial with deg b 3 < n, the only pole outside the contour is the simple pole at T = T * .Thus we obtain overall, Appendix E: Asymptotics of ∆ in the localized phasesFirst we rewrite the expression (85a) in a form that is conducive to saddle point expansion.The diffusion constant