{'fraction_non_alphanumeric': 0.04767744902652154, 'fraction_numerical': 0.014717154683427871, 'mean_word_length': 4.191802626343017, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
DMRG-study of current and activity fluctuations near non-equilibrium phase transitions

Cumulants of a fluctuating current can be obtained from a free energy-like generating function which for Markov processes equals the largest eigenvalue of a generalized generator. We determine this eigenvalue with the DMRG for stochastic systems. We calculate the variance of the current in the different phases, and at the phase transitions, of the totally asymmetric exclusion process. Our results can be described in the terms of a scaling ansatz that involves the dynamical exponent z. We also calculate the generating function of the activity near the absorbing state transition of the contact process. Its scaling properties can be expressed in terms of known critical exponents.

Cumulants of a fluctuating current can be obtained from a free energy-like generating function which for Markov processes equals the largest eigenvalue of a generalized generator. We determine this eigenvalue with the DMRG for stochastic systems. We calculate the variance of the current in the different phases, and at the phase transitions, of the totally asymmetric exclusion process. Our results can be described in the terms of a scaling ansatz that involves the dynamical exponent z. We also calculate the generating function of the activity near the absorbing state transition of the contact process. Its scaling properties can be expressed in terms of known critical exponents.

PACS numbers:
Physical systems that are in contact with two reservoirs at a different temperature or chemical potential, develop a heat or particle current [1]. In macroscopic systems, fluctuations of these currents can often be neglected. As is the case in equilibrium systems, one can however expect that such fluctuations become important in mesoscopic systems and in the vicinity of a nonequilibrium critical point [2].
The statistics of current fluctuations in mesoscopic conductors have received a lot of attention in the past decade [3], since they can, for example, give insight on correlated electron transport. It is nowadays possible to measure experimentally third and higher order cumulants of the current in problems of charge transport [4,5]. Theoretically, these cumulants can be obtained as derivatives of a generating function. This function has many similarities to the free energy in equilibrium systems.
In the present Letter, we focus on the scaling of the current distribution in one-dimensional (classical) stochastic models such as the (a)symmetric exclusion process. This stochastic process is a standard model of non-equilibrium statistical mechanics [6,7]. Rigorous results are known for the current distribution in this model both on a ring and for open boundaries [8,9,10,11]. Moreover, several approximate and numerical approaches to this problem have been developped: simulation techniques that sample rare events [12,13], renormalisation approaches [14] and perturbation techniques [15]. Here we apply for the first time the density matrix renormalisation group (DMRG) to the investigation of current fluctuations. We illustrate the method for the current of the totally asymmetric exclusion process, but the technique is more general. As an example we also present results on the total number of changes of configuration (a quantity that has been called activity [16]) in the contact process [2].
In the totally asymmetric exclusion process (TASEP), each site of a one-dimensional lattice of L sites can be empty or occupied by at most one particle. The dynamics of the model is a continuous time Markov process in which a particle hops to its right neighbor with unit rate provided that site is empty. At the left boundary particles enter the system with rate α, while at the right boundary they leave it with rate β. Asymptotically, the TASEP reaches a non-equilibrium steady state (NESS) in which a nonzero current flows through the system. It is by now well known that the TASEP has three distinct phases [17]: in the low density (LD) phase (α < 1/2, β > α) the average current J (per bond and in the thermodynamic limit) equals α(1 − α) while in the high density (HD) phase (β < 1/2, α > β) it is β(1 − β). Finally, in the maximal current (MC) phase, J = 1/4. Let J L (t) be the total current through all bonds up to time t during a realisation of the stochastic process. The statistical properties of this current can be obtained from its generating function where the average is taken over the realisations of the stochastic process. The distribution of J L (t) at large times can be determined from µ(s, L) by a Legendre transformation while the average current J(L), its variance ∆(L) and higher cumulants can be found as derivatives of µ(s, L): Let σ i be a spin variable which equals 1 if the site i is vacant and −1 if it is occupied. The state of the system is then characterized by the probability P (C; t) to be in a given microstate C = {σ 1 , . . . , σ L }. This probability evolves according to the master equation [18] dP where H is the generator of the stochastic process. The properties of the NESS of the stochastic process can be determined from the (right) eigenvector of H with the largest eigenvalue [18]. Similarly, it is not difficult to show that generating functions such as µ(s, L) can be obtained as the largest eigenvalue of a modified generator H(s) [6]. For the current J L (t), H(s) equals Here we have used the "quantum" notation for stochastic systems [19]. The operators s + i and s − i respectively destroy and create a particle at site i, while n i and v i count the number of particles and vacancies at that site.
Formulated this way, determining µ(s, L) is mathematically similar to finding the ground state energy of a quantum spin or fermion chain. One of the most succesful numerical techniques to study low temperature properties of quantum chains is the DMRG [20,21]. More recently, this method has been extended to stochastic systems [22,23] where the main difference is that in general the generator H is non-Hermitian. Here we apply the method for the first time to operators such as H(s) which are neither Hermitian nor stochastic. We found that with the DMRG it is possible to obtain µ(s, L) numerically exact for systems up to L = 60 with only modest computing facilities. Since there are no essential new ingredients in the method as such [24], we focus here on the results.
Firstly, in order to test the method we have calculated µ(s, L) for the symmetric exclusion process (SEP) for which this function is known for large L values [25]. In the SEP, particles can hop both to the right and left with equal rate. At its boundaries, the system is in contact with particle baths of density ρ a and ρ b . In Fig. 1, we show typical results for Lµ(s, L) for various L-values, together with the exact result (full line). As can be seen, there is a fast convergence towards the asymptotic results.
Going back to the TASEP, we calculated µ(s, L) in the various regions of the phase diagram. The cumulants J(L) and ∆(L) are then determined by numerical differentiation. As an example, we present in Fig. 2 our results for J(L)/(L + 1) and ∆(L) in the MC phase. Also shown are the exact results for J(L) per bond obtained from the matrix product ansatz [17]. The numerical data coincide with the exact ones within the accuracy. The variance of the current was so far not determined exactly. We find that ∆(L) increases as L σ . The corrections to this power law behavior are strong and cannot be neglected for the system sizes we studied. In order to get reliable exponent estimates we have used the BST-algorithm [26]. We find that in the MC phase, but also at the transition line between the MC and LD (or HD)-phase, σ = 1.50(2). In the LD (and HD) phase, σ changes to 2.01(4). Finally, along the coexistence line between HD and LD phases, we find σ = 2.03(3). These results strongly suggest that M = lim L→∞ ∆(L)/L 2 behaves as an order parameter: it is zero in the MC-phase and non-zero in both the LD and HD phases. Given the similarities between the generating function and the free energy, it is natural to ask about the scaling properties of µ near a phase transition. To focus at-tention, we consider the transition line between LD and MC phases (α = 1/2, β > 1/2). We propose that under a rescaling with a factor b the singular part of µ transforms as µ(s, ∆α, L) ∼ b −z µ(b ys s, b yα ∆α, L/b) (5) where ∆α = α − 1/2 and z is the dynamical exponent. We conjecture that z replaces the dimension d that appears in the scaling of the free energy because µ(s, L) is a quantity per "unit of time", whereas the free energy is per unit of volume. The exponent y s is a new exponent associated with current fluctuations, and y α is like a thermal exponent in equilibrium critical phenomena. From (2) and (5), it follows that J(L) ∼ L −z+ys and ∆(L) ∼ L −z+2ys at the transition. From the exact results on J(L) and our data on the variance, we find z = 1.50(2) and y s = 1.50 (2). This value of z agrees with that determined by the Bethe-ansatz [27], z = 3/2, thus providing strong support to the scaling form (5). We conjecture that also y s = 3/2. Finally, y α can be obtained from ∂J(L)/∂α. This derivative can easily be calculated from the exact results, and gives y α = 1/2 [24]. To test our scaling ansatz further, we investigate the variance of the current which should scale as with H a scaling function. To match the numerically determined behavior of the variance in the different phases, H(x) should be constant for x > 0, and linear in x for small x < 0. This implies that M goes to zero linearly as the LD-MC transition line is approached from below. In Fig. 3, we show our data for ∆(L, ∆α)L z−2ys versus L yα (∆α < 0). The scaling is well satisfied and gives exponent values close to those determined above [28]. In a similar way, we also checked that the scaling of ∂ 2 µ/∂s 2 as a function of s at the transition line is well satisfied [24]. Finally, also data on the third cumulant of the current can be well described by our scaling proposal [24]. As a second application of our approach we study the one-dimensional contact process (CP). In this model, each site of a lattice can be occupied by at most one particle. An occupied site becomes empty with rate 1, while an empty site becomes occupied with a rate ζλ/2. Here ζ is the number of occupied neighbors. When λ < λ c the process reaches an absorbing state in which all sites are empty. For λ > λ c , and in an infinite system, the model reaches a NESS with a finite density ρ of particles. The contact process [29] is a standard model for phase transitions out of an absorbing state [2]. It is known from extensive numerical investigations that its phase transition belongs to the universality class of directed percolation [30]. The scaling properties of various quantities near λ c are well characterized [2]. Here we are interested in the activity K L (t) of the model, which gives the total number of changes of configuration in a realization of the process up to time t. The generating function of K L (t) is This function can again be obtained as the largest eigenvalue of a generator which in this case equals (n 0 = n L+1 = 0). A finite system will always reach the absorbing state asymptotically. To avoid this, we allow the creation of particles at the boundary sites. Following the reasonings made for the TASEP, we expect that near the absorbing state transition, π scales as π(s, ∆λ, L) = b −z π(b yK s, b 1/ν ⊥ ∆λ, L/b) Here ∆λ = λ − λ c . The exponents z = ν /ν ⊥ = 1.5805 and ν ⊥ = 1.09684(6) are known numerically [31] while y K is a new exponent. It is possible to express y K in terms of other, known, exponents. From the dynamics of the model one can show that K L (t) [24] obeys In the NESS, the second term in (8) approaches zero, whereas the first one becomes equal to 2Lρ. The scaling of ρ is well known and therefore the average activity should scale as Here F is a scaling function and β = 0.27649 [31]. Since K(L) is also the first derivative of π we get from (7) and (9): y K = 1 + (ν − β)/ν ⊥ = 2.3284. We have calculated π(s, ∆λ, L) using the DMRG. In Fig. 4, we show our results for the variance of the activity as a function of ∆λ and L. At criticality, we find that the average activity diverges as L .746 (2) , while its variance goes as L 3.08 (2) . These exponents are close to −z + y K = 0.7479 and −z + 2y K = 3.0763 predicted by the scaling (7). Other evidence of (7) can be seen in Fig. 5 where we present a scaling plot of K(s, L) = ∂π/∂s as a function of s at ∆λ = 0. This quantity should scale as L −z+yK G(L yK s).
The numerical data again support this prediction [28]. In summary, we determined the generating function of the current in the TASEP and of the activity in the CP with the DMRG. We proposed a scaling form for these generating functions which is supported by all numerical data. We believe that this scaling is quite general and can be applied to other models as well.