Totally asymmetric exclusion process fed by using a non-Poissonian clock

In this article we consider the one-dimensional totally asymmetric open-boundary exclusion process fed by a process with power-law-distributed waiting times. More speciﬁcally, we use a modiﬁed Pareto distribution to deﬁne the jump rate for jumps into the system. We then characterize the propagation of ﬂuctuations through the system by kinetic Monte Carlo simulations and by numerical evaluation of the steady-state partition function.

One generalization of the TASEP studied recently is the model "with finite resources" where several TASEP chains compete for a fluctuating number of incoming particles available because of restricting the total number of particles in the system [22]. This is motivated by a biological application where the number of ribosomes available fluctuates, see also Ref. [14]. Another modification of input rates is periodic driving [23]. Also the fully non-Poissonian exclusion process has been investigated [24]. In this article we consider a single TASEP where the chain is fed with the jump rate of incoming particles being non-Poissonian with a wide distribution of waiting times for a new particle to come in. Below we call the TASEP with this modification the fed TASEP (FTASEP). In practical applications, the reason for such a waiting time distribution might be fluctuations of the number and density of particles available or the chemical potential due to other processes in the environment.

II. THE SIMULATION MODEL
In the standard one-dimensional open-boundary TASEP model there are two rate parameters, the input rate α and the output rate β (standard notations) with the rate for jumps in the "bulk" taken to be unity. The site occupation numbers are denoted by n i , where i = 1, . . . ,L can have values n i = 0,1 for empty and occupied sites, respectively, and L is the size of the system.
In the FTASEP the constant input rate α into site i = 1 of the standard TASEP is replaced by where the random variable ξ is drawn from a modified Pareto distribution, where the decay of the distribution is determined by the parameter a. The rate α(t) (and the related expected waiting time) is defined when the lattice site p = 1 becomes empty and it keeps its value until the site gets occupied. The other parameters A, γ , and δ are constants. Note that for the decay parameter we use here notation a instead of notation α, which is reserved for standard use in the context of exclusion processes.
In this paper we choose γ = δ = 1. This leads for A = 1 to input rate statistics, where and The corresponding effective waiting times are determined by the inverse quantities, valid and finite for a > 1, and finite for a > 2. It turns out that for small values of a because of the diverging waiting times, there is no reasonable way to choose the constant A so that the input rate would correspond to a given constant rate for all values of a. Since arbitrarily long waiting times are possible, there are time epochs when the lattice, no matter how large, is empty and the observed fluctuation statistics will depend on the length of the simulation. However, by choosing A = (1 + a)/a the mean input rate is constant and For increasing a well above 2, the dynamics of the FTASEP approaches that of a TASEP with a constant input rate. On the other hand, for small a the average waiting time diverges.
Since the average values of the density and current in the standard TASEP are ρ = 1/2 and J = 1/4 in the maximal current phase, long waiting times in the FTASEP for small a lead to a situation, where ρ = 1/2 and J = 1/4 are beyond reach regardless of the choice of the parameter A as the TASEP dynamics of the bulk restricts the current for times during time epochs when there are particles in the lattice.
In the simulations of this model the dynamics is generated by using a rejection-free continuous-time kinetic Monte Carlo method [25,26] by which we can handle arbitrarily long (and short) waiting times. For the TASEP simulations we choose input and output rates α = β = 1 for which there exist analytical results in the literature, see, e.g., Refs. [3,9]. For the FTASEP we choose the Pareto parameter values γ = δ = 1 varying the decay parameter a with A = (1 + a)/a and β = 1 to allow comparison with the TASEP.

III. RESULTS FOR LOCAL DENSITY FLUCTUATIONS BY SIMULATIONS
First, we compare the density fluctuations for the FTASEP with those reported for the TASEP by Derrida et al. in Ref. [9], in particular the exact analytical results in their seminal article. Following their approach, which facilitates the comparison of analytical and numerical results for the TASEP, we calculate the quantities s 2 (x) and s 4 (x), essentially the variance and skewness for local density fluctuations.
For direct comparison with their results, in Fig. 1 we use their division of the lattice into intervals such that the value of parameter x corresponds to an interval L(1 − x)/2 i L(1 + x)/2, where i is the index of the lattice site. With this choice, for increasing x ≡ 2k/L with k = 1,2, . . . ,L/2, the interval expands symmetrically from the middle of the system such that for x = 1 the interval covers whole system. For a such division of the lattice, Derrida et al. [9] give for s 2 (x) and s 4 (x) in the TASEP analytical formulas [see Eq. (3.7-8) in their paper], and they test them by simulations of the TASEP [see Fig. 1 in their work]. In Fig. 1 we see how the behavior of these quantities for the FTASEP for small values of a considerably differ from those for the TASEP. (Note that the vertical scale in our plot is logarithmic.) For increasing a, the results for the FTASEP approach those for the TASEP.
Since in the choice above of intervals as in Ref. [9], the interval expands symmetrically from the middle of the system (which for symmetry reasons is optimal for the analytical study of the TASEP), it does not tell about the differences between the two halves of the lattice (relevant for the FTASEP). Therefore, we next calculate in Fig. 2   The thick black curves are the theoretical results for the TASEP from Ref. [9]. Our data for a = 3.0,4.0 are hardly discernible from their theoretical curves. In these plots, for increasing x the interval expands symmetrically from the middle of the lattice to allow comparison with the theoretical results of Ref. [9].
Evidently the anomalous fluctuations in the FTASEP are more pronounced in the first half of the lattice, the ratios in general increasing for increasing a and increasing for increasing x (for the TASEP the ratios would be one [9]).
The order of the curves in Fig. 2, i.e., the data for a = 1 being between data for a = 2 and a = 3, reflects the fact that for a = 1 the waiting time distribution becomes so wide that the averaged particle density is much lower than for a 2, cf. Eqs. (5) and (6). We will return to this point shortly.
We next turn to the question of what controls local fluctuations in the FTASEP. Obviously the local particle density is an important parameter, but we expect it not to be all there is to say. To study this question, we divide the lattice into ten consecutive intervals of equal length m = 50 for L = 500, each with index p = 1, . . . 10 corresponding to intervals (p − 1)m + 1 i pm, and show the scaled kurtosis and skewness as a function of the local particle density n p /m.
First, it is evident that the particle density in the FTASEP tends to increase for increasing a, and for large a the density in the middle of the system is close to 1/2. For small a, as explained above, there are long epochs of time such that the lattice is empty and the density is very small. This is seen in the behavior of the kurtosis s 4 /(s 2 ) 2 for low densities in Fig. 3(a), which shows that the average fluctuations are completely controlled by the density for a 1. To make this more evident, we show by the thick black curve the (trivial) single-site occupancy fluctuation statistics corresponding to the given average density.
On the other hand, for a > 1 the behavior across the system for a given value of a is more complex and, in addition to the local density n p /m, the value of a becomes relevant as seen for large densities in Fig. 3(b). For very large values of a the value of kurtosis approaches from above the normal distribution limit s 4 /(s 2 ) 2 → 3. Note that larger values of n p correspond to intervals closer to the input site i = 1 so that each curve for a given value of a in Fig. 3 (and in Fig. 4 below) should be read from right to left.
For the data in Fig. 3(b) one can attempt a data collapse. A reasonable one is obtained by scaling the x axis (the horizontal axis) as x → x + (2/5) a and the y axis (the vertical axis) as y → (y − 3) (a−1)/2 . Because of the unavoidably short span of the data, one might replace (2/5) a above by e −a . Furthermore, for a > 2 no shift is needed on the horizontal axis since all the curves in that region start to deviate from the horizontal part y ≈ 3 around x = 1/2.
For the local skewness s 3 /(s 2 ) 3/2 , the simulation data across the system shown in Fig. 4 display similar qualitative features. We note that the value n p /m = 1/2 deep in the bulk corresponds to the theoretical value for the TASEP with α = β = 1 [3,9]. That value of local density is reached in our parametrization of the FTASEP only for a 2. However, if in the simulations we multiply the input rate coefficient A = (1 + a)/a in Eq. (7) by some large factor, then it is possible to reach n p /m ≈ 1/2 for smaller values of α. Then also, for example, multiplying A by 100 the kurtosis becomes close to three and skewness close to zero for α 1.4 with L = 500 (data not shown). The locations of the apparent transitions from a simple behavior for small values of a to a complex behavior for intermediate values of a and then again to a simple behavior for a large value of a also depend on L, but the qualitative overall picture remains the same.

IV. THEORETICAL CONSIDERATIONS: THE PARTITION FUNCTION
Analytical calculations for systems with (part of the) processes obeying the nonexponential clock are challenging. In this section we present one theoretical machinery; an interested reader will find the numerical results obtained by using it for physical quantities in Sec. V.
As a starting point for an effective theory we take the well-known mapping between the TASEP and the Brownian excursions from Ref. [9] by Derrida et al. They show that the TASEP states can be related to random walks with successive steps of length 0, ±1 such that +1 corresponds to an occupied site, −1 corresponds to an empty site, and 0 corresponds randomly to an occupied or an empty site [27].
For this approach, to write the partition function, a few functions need to be defined: First, for the excursion statistics in the TASEP, y and z being two consecutive excursions, we take [9,27,28] H (y,x) = ye −y 2 /x .
The variable x ∈ [0,1] sets the linear scale such that the whole interval [0,1] corresponds to sites i = 1, . . . ,L in the TASEP. Second, for mapping onto densities via related local chemical potentials we also need [9] I (y,x,μ) = e −(2μ−y) 2 /x , J (y,z,x,μ) = e −(2μ+y−z) 2 /x , where μ is a rescaled local density [μ p in Eq. (15) below]; for the rescaling see (20) and Ref. [27]. For the FTASEP, instead of F (y,x) in the TASEP we use at the input end of the lattice the slowly decaying Pareto form, This functional form is chosen to mimic a wide distribution of input rates, and the factor 1/ √ x in the denominator is included for similar passage from x = 1 to x = 1/k with the functions F (y,x) for the TASEP and P (y,x,a) for the FTASEP. We also note that the derivative with respect to y of P (y,x,a) is of the same form but with larger a.
For the FTASEP with a given value of a, still suppressing dependence on x, we replace in Eq. (15) above F (y 1 ) → P (y 1 ,a) and use the notations W k → W a k and Z k → Z a k . For the FTASEP we use this for the case where Eq. (7) is valid. Physically, the choices made above lead to an effective model where the heavy-tail fluctuations at the input end of the lattice in the FTASEP get propagated through the system by the TASEP dynamics.
The integration in Eq. (15) over μ p variables means that also noninteger values of the corresponding occupation numbers n p are included in averaging, which is not a problem as long as the system size is large and the location p contains a large number of lattice sites [27].
One can calculate the integrals in Eq. (15) for a selected constant value of x and from it get the partition function and, e.g., averages μ ν p for an other value of x by a simple change in integration variables; in the actual calculations we first set x = 1. Note also that it is possible to let in all of the above x depending on p with p x p = 1 such that x p L is the length of a fraction x p of the system size L for each p.
As an example, one might set x = 1/k = 1/L, where L is the number of sites in the exclusion process, and then use the connection [9,27], to obtain with the interpretation that n p is the average occupation number of a single site. In practice one can perform a calculation with x = 1 from which In any case, with L = mk = m/x the local particle density ρ p in section p = 1,2, . . . ,k from Eq. (17) is written as and

A. The TASEP
To compare with, we first consider the model for the standard TASEP. Then for x = 1, and Z k (x) = x k−1/2 Z k (x = 1). Expectation values of specific interest (see later) are for the first site (p = 1) and last site (p = k) as follows: For x = 1/k = 1/L we would have μ p = √ x μ p x=1 and μ 2 p = x μ 2 p x=1 for all 1 p k. The formulas corresponding to those in Eqs. (23) and (24) for sites 1 < p < k become more complicated and are not shown here. In terms of lattice site occupation we get for x = 1/k = 1/L, where (δn 1 ) 2 ≡ (n 1 − n 1 ) 2 . By symmetry in the TASEP with α = β = 1, we have n k = 1 − n 1 and (δn k ) 2 = (δn 1 ) 2 . Related analytical and numerical results for a certain kind of division of the lattice into large segments for the TASEP with α = 1 = β can be found in Ref. [9]. The theory becomes exact for large L and m only [9,27], and then n p ∈ [0,1]. We note that in the interpretation L = k we would have m = 1 and for k → ∞ then from the formulas above we find n 1 → 1/2 + 1/ √ π ≈ 1.064 > 1 and n k →

B. The FTASEP
Also for the FTASEP there are a few cases with certain values of a, p, k, and some μ ν p for which the partition function can be expressed analytically in terms of, e.g., hypergeometric functions and their special cases, such as the exponential integral and the Dawson integral, the actual formula given by, e.g., Mathematica depending on the order of integrals. Most of such expressions are cumbersome and difficult to use and do not provide much physical insight, see the Appendix. It is much more useful to express some of the results with the integral over y 1 left undone. Then, for the first site p = 1 with L = k for simplicity [see Eq. (19)], we find and where In Eqs. (27)-(29) some dependence on k has been canceled, see the Appendix for detailed formulas. Note the dependence on the system size L = k in the decay parameter of the integrands. These integrals are easy to calculate accurately by numerical integration, whereas the original (k 2 −k)-dimensional integrals with the integrand proportional to factors exp[−2y 2 p ] and exp[−2y p y p+1 ] become prohibitively difficult to integrate numerically for k > 4.
From the point of view of the propagation of fluctuations, the site at the end of the system (p = k) is of interest. For the variance at p = k we need, e.g., where c k = (k − 1)(4k − 7) replaces unity in Eq. (28) and where b k = (k − 2)(k − 1). With the error function being defined as an integral, Eq. (31) is a two-dimensional integral.

V. RESULTS OBTAINED BY NUMERICAL INTEGRATION OF THE PARTITION FUNCTION
We next turn to results obtained by computing numerically the partition function and averages for the local densities n p and their fluctuations based on Eq. (15) for the TASEP and the FTASEP. For this, we start by computing μ p x=1 and μ 2 p x=1 for the FTASEP with given values of a, where p = 1, . . . ,k with k = 3-6. For the TASEP, the integrations involved can be performed exactly as shown in Sec. IV, whereas for the FTASEP some numerical integration is needed in most cases.
To study large systems of L sites, we then set m = L/k with x = 1/k and μ p = μ p x=1 √ x from Eq. (19), so local densities are with m being the number of lattice sites in interval p.
In Figs. 5 and 6 we show the local particle density and its variance for the TASEP and the FTASEP for a convenient selection of system sizes. The dashed lines display the local particle density calculated from the exact formula for singlesite occupation numbers in the TASEP with α = β = 1 given in Eq. (48) of Ref. [3]. Comparing the TASEP result for n p /m by numerical integration in Figs. 5(a) and 6(a) with the exact result for the TASEP, we conclude that numerical integration works well for large system sizes, and we can proceed to study fluctuations.
To make the structure of the curves in Figs. 5(b) and 6(b) more visible for all values of a we have normalized the data by dividing (δn p ) 2 by the variance (δn 1 ) 2 . The data for the normalized variance show that the fluctuation profile in the FTASEP differs from that in the TASEP more than could be expected from the local density data. It is even a monotonically decreasing function of the distance from the input site for a = 1-3. For increasing a it approaches the symmetric form of the fluctuation profile of the TASEP. The absolute unnormalized values (not shown) of the variance increase fast for decreasing a as was seen in the simulation data of Fig. 1(a) in the preceding section. Note that the data in Figs. 5(b) and 6(b) are due to the normalization the same, but the absolute values of (δn p ) 2 differ by a factor of 16, being larger in the first case.
It is of interest to compare the theoretical (T ) results obtained by numerical integrations with those obtained by simulations (S) for the FTASEP. We continue to use the simulation data with L = 500. Now in the effective theoretical model there are different ways to treat the first site of the lattice. In Fig. 7 we show the result of three choices: the length of the first interval in integrations m 1 = 1,10,84 lattice units and the rest of the lattice is divided evenly such that the total number of intervals is 6, i.e., p = 1,2, . . . , 6 with an evenly divided system the agreement is perfect, i.e., n p T / n p S = 1. In simulations for a = 1 the density is so low that no reasonable comparison can be made. Note that we have not attempted to fine-tune the density of the first interval in the simulation to correspond to that of the theoretical calculation. We conclude that the theoretical and the simulation approaches agree reasonably well. However, the agreement does depend on a for the given choice of m 1 since for different values of a the correlations in the input noise die out in different length scales.

VI. CONCLUSIONS
To summarize, we used continuous-time Monte Carlo simulation and numerical evaluation of the partition function to study the effect of non-Poissonian power-law-governed feeding on the totally asymmetric exclusion process. This model we call the FTASEP. In particular we studied local density fluctuations across the system to characterize propagation of the fluctuations through the system. Depending on the parameter a defining the decay of the waiting-time distribution, we observe simple density-controlled low-density behavior for small values of a, complex behavior for intermediate values of a, and the approach of the FTASEP to the TASEP for large values of a. It would be of interest to extend the study of the combined effect of non-Poissonian feeding and finite size also to the statistics of current fluctuations.

APPENDIX: SOME ANALYTICAL RESULTS
For the FTASEP for some values of a the partition function (1 + y) a+1 dy.
Corresponding formulas can be obtained for other sites of p.
The most interesting and useful, i.e., relatively simple, result is obtained for p = k, see Sec. IV B.
of a Brownian excursion, these being independent processes. Then, integrating over a section [x p−1 ,x p ] of the system, one finds L [ρ(x) − 1/2]dx i (n i − 1/2) = μ p √ L where the summation over lattice sites i runs from Lx p−1 to Lx p , valid for large L and large Lx p . This connection is obtained in Ref. [9] by using the (nowadays standard) matrix ansatz method for the TASEP.