Tensor-Network Approaches to Counting Statistics for the Current in a Boundary-Driven Diffusive System

We apply tensor networks to counting statistics for the stochastic particle transport in an out-of-equilibrium diffusive system. This system is composed of a one-dimensional channel in contact with two particle reservoirs at the ends. Two tensor-network algorithms, namely, Density Matrix Renormalization Group (DMRG) and Time Evolving Block Decimation (TEBD), are respectively implemented. The cumulant generating function for the current is numerically calculated and then compared with the analytical solution. Excellent agreement is found, manifesting the validity of these approaches in such an application. Moreover, the fluctuation theorem for the current is shown to hold.


Introduction
Away from equilibrium, open systems composed of atoms and molecules manifest macroscopic currents dissipating energy and producing thermodynamic entropy [1][2][3][4]. These currents are intrinsically fluctuating with the probabilistic nature fully characterized by their large deviation properties [5][6][7]. The theoretical progress in the last three decades reveals that, for nonequilibrium systems in steady state, there exists a strikingly simple and general relation for the large deviation properties of the fluctuating currents. This relation is nowadays known as the fluctuation theorem [8][9][10][11][12][13][14][15]. At the mesoscopic level of description, the evolution of the coarse-grained state of the concerned system can be modeled as a Markov jump process. Accordingly, the probability distribution for the system states obeys a master equation. To evaluate the large deviation function for the current fluctuations, a direct method is to numerically simulate the Markov jump process with the Gillespie algorithm [16,17], which generates a sample of random trajectories. Then, we perform counting statistics for the integrated current, and the probability distribution of the current can be finally constructed from the statistics [18]. However, such direct sampling method encounters problem as it can not accurately estimate probability of rare events or large deviations given limited number of data. An alternative method is to first calculate the cumulant generating function for the currents, and then obtain the large deviation function through the Legendre-Fenchel transform. The cumulant generating function turns out to be the leading eigenvalue of the tilted generator of master equation including counting parameters [18]. Thus, it has transformed into an eigenvalue problem.
For a spatially extended system, a coarse-grained description at the mesoscopic level requires partitioning the system in space. This leads to the characterization of the states with many random variables. Accordingly, the dimension of the tilted generator is very large, growing exponentially with the number of these variables. The same issue arises in quantum many-body systems for their underlying Hilbert space. Over the past few decades, tensor networks [19][20][21][22][23][24][25][26] have emerged as one of the most powerful tools for handling such numerical complexity. Although originally developed in the context of condensed matter physics, tensor networks are presently finding applications in stochastic dynamics [27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42]. Indeed, there is a close analogy between the classical stochastic systems in nonequilibrium steady state and the quantum many-body systems in groud states. In this regard, two tensor-network approaches stand out, namely, Density Matrix Renormalization Group (DMRG) [22,23,43,44] and Time Evolving Block Decimation (TEBD) [45,46]. For quantum many-body systems, the former is a variational approach that enables the systematic search for ground states, while the latter an approach that evolves the states. If implemented Figure 1: Schematic representation of a transport channel in contact with two reservoirs.
in imaginary time, the latter can also be used to find the ground states. Correspondingly, for classical stochastic systems, these two approaches are naturally well suited for obtaining the nonequilibrium steady states.
Here, our purpose is to provide a pedagogical demonstration of how DMRG and TEBD are used to calculate the cumulant generating function for the current in a boundary-driven diffusive system. The rest of this paper is structured as follows. The stochastic description for the diffusive process is introduced in Section 2 where we discretize the system in space and then establish the master equation. The counting parameter for the particle transitions is included in the master equation, and consequently the cumulant generating function for the current is analytically obtained as the leading eigenvalue. Section 3 is devoted to the DMRG approach to counting statistics. The probability distribution is represented as a matrix product state (MPS). The generator for the master equation is formulated in terms of the local creation, annihilation, and particle number operators, and is then constructed as a matrix product operator (MPO). In Section 4, the TEBD approach to counting statistics is concerned. The small-time evolution operator for the tilted master equation is approximately decomposed into a ordered product of local ones. The numerical details and results are presented in Section 5, followed by the concluding remarks given in Section 6. The cumulant generating function is formally derived in Appendix A.

Boundary-Driven Diffusive System
The system we consider here is modeled as a one-dimensional channel in contact with two reservoirs at the ends, as illustrated in Figure 1. Distributed inside the channel are mobile particles with their densities expressed as a function of position, n(x). At the boundaries, because the reservoirs are supposed to be arbitrarily large, the densities of particles are constant in time with n L for the left reservoir and n R for the right reservoir. In the case that n L = n R , the system is driven out of equilibrium and a directional current flows from one reservoir to the other. The mobile particles undergo the Brownian motion [47][48][49][50], so that the time evolution of the densities is given by the diffusion equation where D is the diffusion coefficient determined phenomenally by direct observation of the random walk of a tagged Brownian particle according to the formula Here, x 0 (respectively, x t ) is the position at the initial time (respectively, time t), and · denotes the average taken over the trajectories.

Stochastic Description
At the mesoscopic level, the time evolution of the distribution of mobile particles in the channel can be described as a Markov jump process, which is formulated in terms of a master equation (see Ref. [51] and references therein). The fluctuations down to the mesoscopic scale can be fully characterized by such a stochastic description. For this purpose, the channel is spatially discretized into L cells, each of width ∆x and containing some number of mobile particles. These cells are labeled with indices i = 1, 2, · · · , L. For notational convenience, the indices i = 0 and i = L + 1 are respectively used to refer to the cells for left and right reservoirs. The number of particles in both cells for reservoirs are maintained constant in time, and are given by N 0 ≡N L and N L+1 ≡N R . The system state is specified by the numbers of particles in the cells. These numbers change in time due to the random transitions according to the kinetic network with the transition rates linearly depending on the particle number in the departure cell, W Here, k is the rate constant related to D by k = D/∆x 2 which can be obtained through the continuum limit 1 . The probability P(N, t) that the cells contain the particle numbers N for time t obeys the master equation whereL denotes the generator for the time evolution. The stationary solution of the above master equation is given by a product of Poisson distributions, with the mean number of particles corresponding to a linear profile of concentration,

Counting Statistics
To investigate the current fluctuations for the particle transport, we now perform counting statistics. For this purpose, the master equation (3) is modified to include a counting parameter λ for particle transitions between the left-reservoir cell and the first one for the discretized channel, yielding 1 The mean current is expressed as where Ω and Σ denote the volume and section area of a cell, respectively. From this expression, k = D/∆x 2 is immediately arrived.

The eigensolution of this tilted master equation has the form
where {C i } are constants depending on λ, and Q(λ) is the cumulant generating function. Substituting this solution into the tilted master equation and after some algebraic manipulations, we obtain If we define A = ln N L /N R as the affinity [52] for driving the system out of equilibrium, Then, the cumulant generating function (8) satisfies the Gallavotti-Cohen symmetry [11,53,54] which is also known under the name fluctuation theorem. The statistical cumulants can be obtained by taking successive derivatives with respect to the counting parameter λ. For example, the mean current and its diffusivity are given by For nonequilibrium systems described by master equation, the affinity can be systematically obtained through Schnakenberg graph analysis [55]. More details about the counting statistics for this diffusive process can be found in Ref. [14]. A formal derivation of the cumulant generating function (8) is provided in Appendix A. From the numerical perspective, the determination of cumulant generating function requires to calculate the leading eigenvalue of the tilted generatorL λ defined in Eq. (6). In principle, this can be done with various numerical routines for diagonalizing a matrix. However, this is often very difficult in practice. Since the dimension of the tilted generator grows exponentially with the number of the discretized cells. This motivates the need of tensor networks to tackle such a problem.

DMRG Approach to Counting Statistics
In this section, we present in detail the DMRG procedure to perform counting statistics for the stochastic process in the preceding section. The key differences from the standard usage for quantum lattice systems are highlighted.

Notations
First, it is helpful to introduce the quantum-mechanical notation for subsequent description convenience. We use |N i for the state that the i-th cell contains N i particle and |N for the system state specified by the numbers Accordingly, N i |, N|, and F † (N)| are used to denote their adjoints. Following the conventions in quantum mechanics, 2-norm and orthogonality are assumed, Now, an issue arises that the probability distribution is of 1-norm, N P(N) = 1. We can introduce the 2-norm probability distribution,P so that P † (N)|P(N) = 1. The 1-norm probability distribution can be recovered, There is a one-to-one correspondence between these two kinds of probability distributions. In practical implementation in code, we always use the 2-norm version, since many tensor-network techniques already developed for quantum-mechanical systems can be utilized.

MPS and MPO
In quantum many-body systems, the presence of the celebrated area laws [56] for the entanglement entropy enables an efficient numerical representation of the ground state as an MPS [57] (see Figure 2a for illustration). Here, MPS can also be used to represent the stationary solution for classical stochastic systems in nonequilibrium steady states. Moreover, we notice that the feature of joint product in solutions (4) and (7) imply that the bond dimensions of their MPS representations are both one. This is very remarkable that the computer memory used to store the stationary solution can be very low, scaling linearly with the number of cells. The physical indices denote the numbers {N i }.
Next, we show how to construct an MPO (see Figure 2b) for the tilted generatorL λ . This constitutes the most crucial step in applying the DMRG approach to counting statistics for the stochastic process. We adopt the formalism similar to the second quantization. The particle transitions between two neighboring cells can be associated with an annihilation operator acting on one cell and an creation operator acting the other. Following this reasoning, we introduce the local annihilation operator a − i and creation operator a + i , which are respectively defined by It should be noted that the operator definitions by Eqs. (18)- (19) only apply to the cells for the channel, 1 ≤ i ≤ L. For the cells corresponding to the reservoirs, i = 0, L + 1, these operators should be defined differently, Because the particle numbers in the reservoir cells are kept constant. To account for the probability loss in the tilted master equation, we should still introduce another operator a i defined by which can be intuitively termed as particle number operator. This definition applies to any cell. In this way, the tilted generator is expressed aŝ where ⊗ denotes the tensor product and the counting parameter λ is included for particle transitions between the left-reservoir cell and the first one in the discretized channel. This tilted generator is now split into two partsL with each constructed as an MPO, where the boundary tensors are explicitly expressed according to Eqs. (20)- (22) and contracted with their neighboring ones. So, there are L tensors in each MPO as expected. If we denoteŝ then the tilted generator is constructed as an MPO, whose bond dimension is 7. It should be noticed that the tilted generatorL λ is intrinsically non-symmetric or non-Hermitian, implying that (i) the right eigenvectors are not identical to the corresponding left eigenvectors; and (ii) the right (left) eigenvectors are not mutually orthogonal. This is the basic difference from the case of quantum Hamiltonian.   (33) with, e.g., F 2 being the canonical center.

Optimization
The aim of DMRG approach is to variationally find the eigensolution from which the leading eigenvalue can be computed. This can be formulated as an optimization problem under the condition F † (N)|F (N) = 1. The graphical representation of this normalization is shown in Figure 3 where the canonical form of MPS is presented due to the gauge freedom. The advantage of this canonical form lies in the efficient contraction of an MPS with its adjoint; we only need to contract the canonical centers. The shift of the canonical center can be realized with the singular value decomposition. The problem (31) can be further reformulated into a constrained optimization problem using a Lagrangian multiplier −Q, We now adopt the single-tensor optimization procedure, where only one tensor F † i is varied at once while keeping all others fixed, so we have Figure 4 shows the graphical representation of the consequence after this variational procedure. Because two MPSs are both in canonical form with the centers at F i and F † i , Eq. (33) simplifies tõ whereL λ the operator indicated as the red part in Figure 5. This now represents a standard eigenvalue problem. The dimension of the operatorL λ is greatly reduced, making it possible to be handled with normal numerical routines. We now need to find the right eigenvector F i ofL λ corresponding to the largest eigenvalue −Q. This can be achieved through the iterative method. Since eL λ > 0, the Perron-Frobenius theorem applies, and the leading eigenvalue −Q ofL λ corresponds to the real maximum eigenvalue e −Qt of eL λ t in magnitude (for some t > 0). Consequently, the right eigenvector can be asymptotically evaluated as which is then normalized Here, F i (0) is an arbitrary vector assumed to have the component of the desired eigenvector. The matrix exponential can be computed with many numerical routines [58]. Then, we have After completing the optimization on F i , we move the canonical center to a neighboring one and continue such procedure. This operation is performed back and forth between the first and last tensors until the leading eigenvalue converges. With different values of λ, we can numerically construct the desired cumulant generating function Q(λ).

TEBD Approach to Counting Statistics
In this section, we present another tensor-network approach to counting statistics for the stochastic current. Now that we want to obtain the eigensolution of the tilted generatorL λ corresponding to the leading eigenvalue, can we achieve this with direct iteration? The answer to this question is yes. This can be done with TEBD which is an algorithm relying on Trotter-Suzuki decomposition [59]. The desired eigensolution is asymptotically evaluated as where |F (N, 0) is an arbitrary initial distribution. At this point, the dimension of the evolution operator is exponentially large. For small-time step δ, however, eL λ δ can be decomposed into an ordered product of local ones of a much smaller dimension. Because the particles hop between neighboring cells, the tilted generatorL λ can be split into two parts,L λ =L odd +L even , withL odd = l 1,2,λ + l 3,4 + · · · , L even = l 2,3 + l 4,5 + · · · , where l i,i+1 represents the local operator acting simultaneously on the i-th and (i + 1)-th cells. The local operators in each part can be diagonalized efficiently and mutually commuting. It is noted here that the counting parameter is included in l 1,2,λ . The exact evolution operator can be decomposed to any order.
Here, we give the second-order one, as it is commonly used. If the transport channel is discretized into 5 cells, then the local operators are explicitly given by where I i stands for the identity operator acting on the i-th cell. The local evolution operators are alternatively applied to the MPS (see Figure 6 for illustration). Once we obtain the desired eigensolution |F (N) , the cumulant generating function can be computed, with δ being a very small real number.

Numerical Results
With the tensor-network approaches having already exhibited in the preceding sections, we now dive into the numerical details and present the results. First of all, we define the system with specific numbers which are listed in Table 1. We realize that an arbitrary number of particles can concentrate in a discretized cell of the transport channel. However, the probability of the event that the particle numbers {N i } far exceeds their mean value is exponentially small. This means that we can truncate the state space. The particle numbers in a cell only take the values N i = 0, 1, 2, · · · , N max − 1. Accordingly, the operators a − i , a − i , a defined by Eqs. (18)-(23) have finite dimension N max × N max , as also specified in Table 1. Choosing a proper value of N max is very crucial in numerical computation. Smaller value leads to inaccuracy of the results while larger value costs more computational resources.
We are now interested in the stationary distribution of the system in steady state. Numerically, this can be obtained with both DMRG and TEBD approaches for the generatorL. The 2-norm probability   Table 1.  Table 1.  Table 1.
that the i-th cell have N i particles is calculated asP(N i ) = N\NiP (N). This probability can be obtained by contracting |P(N) represented as an MPS with an auxiliary MPS given by where 1 j is a vector with all elements one, and V i a vector with N i -th element one and others zero. The tensor product ⊗ creates a bond of dimension 1, linking the vectors to form an MPS. Then, the 1-norm probability distribution P(N i ) can be recovered through Eq. (17). The mean value can also be calculated, N i = Ni N i P(N i ). Figures 7 and 8 respectively plot the profile of mean particle numbers and the probability distribution of particle number in a given cell. The profile is linear, in agreement with the theory (Eq. 5). The probability distribution is Poissonian, also in agreement with the theory. We notice in Figure 8 that the probability P(N max − 1) is indeed vanishingly small, manifesting that there is no over-truncation of the state space. Next, we numerically calculate the cumulant generating function Q(λ). The counting parameter λ takes several discrete values from 0 to the affinity A = ln N L /N R . The result is presented in Figure 9. The agreement between the numerical results and the theory (Eq. 8) is also striking. These agreements strongly manifest validity of both DMRG and TEBD approaches in such an application. Although DMRG and TEBD both give nice results for the system specified with the numbers in Table 1. DMRG is in general much more efficient than TEBD. Because, in the adopted version of DMRG we deals with one tensor at a time, whereas in TEBD we are simultaneously deals two tensors at a time. On the other hand, DMRG also gives more accurate result than TEBD does. Because, in TEBD there is an error associated with the Trotter-Suzuki decomposition.
The computer program for numerical computation is coded in C++ [60] with the ITensor library [61]. Readers are referred to Ref. [62] where a comprehensive and up-to-date snapshot of software for tensor computations is assembled.

Conclusion and Perspectives
In this paper, we studied the diffusive process of Brownian particles in a channel, with a particular interest in the current fluctuation. A nonequilibrium constraint is imposed by placing the channel in between two particle reservoirs with different concentrations. To describe the randomness of the particle transport, the channel is spatially discretized into cells so that a master equation can be written down for the time evolution of the system state specified by the particle numbers in the cells. The counting statistics is performed by modifying the master equation to include a counting parameter for particle transitions between the cells. This gives a tilted generator for the time evolution of the corresponding master equation. The cumulant generating function for the current can be obtained from the leading eigenvalue of the tilted generator. Because the transition rates are linear with the local particle concentration, an analytical solution for the cumulant generating function is obtained, exhibiting the Gallavotti-Cohen symmetry.
The focus of this paper is on the application of tensor networks to counting statistics. We presented two approaches, DMRG and TEBD, in details. The tilted generator is written in terms of the local annihilation, creation, particle number operators associated with random particle transitions. The annihilation and creation operators for the reservoir cells are defined differently to make the nonequilibrium constraint sustainable. We also emphasized the differences from their conventional applications for quantum many-body systems. We introduced the so-called 2-norm probability distribution in order to make use of the tensor-network techniques already developed in the quantum-mechanical setting. The conventional probability distribution can be recovered once the corresponding 2-norm probability distribution is obtained. In the DMRG approach, the key point is to construct an MPO for the tilted generator. This is done by constructing two MPOs for different parts of the tilted generator, and then performing tensor summation of these two MPOs. The resulting MPO is of bond dimension 7. We formulated the DMRG approach in such a way that a single tensor is optimized at a time. In the TEBD approach, the tilted evolution operator is approximately decomposed into an ordered product of local ones. Then these local evolution operators are alternatively applied to the an arbitrary distribution represented as an MPS. In both approaches, what we finally obtain is the eigensolution which is then used to calculate the cumulant generating function. When the counting parameter is set to zero, the leading eigensolution of the generator for the master equation represents the stationary distribution for the system in steady state. We compared the profile of mean particle numbers and the probability distribution of particles in a single cell with those given by the analytic solution. We remark here that the joint probability distribution of particles in several cells can also be obtained. We also compared the numerically obtained cumulant generating function with its analytical counterpart. Striking agreements are found for all these comparisons between the numerical results and theory.
The tensor-network approaches provides a good means to investigate the fluctuations for a vast majority of realistic systems, e.g., diffusion-reaction systems with many species of particles, systems with long-range electrostatic interactions between particles [63][64][65], or systems of complex geometries with more than two particle reservoirs. Moreover, we notice that the tensor-network approaches, especially DMRG, enable a very accurate determination of the cumulant generating function. This makes it possible to test Onsager reciprocal relations [66][67][68] and their generalizations to nonlinear regimes [13,[69][70][71][72][73] for spatially extended systems. It also makes it possible to study the mathematical properties for statistical cumulants. For example, whether the limited number of cumulants can be used to infer the affinity [74]. To conclude, an accurate determination of the cumulant generating function possibly reveals many underlying physics.

A Formal Derivation of the Cumulant Generating Function
In this appendix, we provide a formal derivation of the cumulant generating function (8) with the method introduced in Ref. [75]. We now consider the probability P(Z, N 1 , · · · , N L , t) that the cells contain given particle numbers and that the signed cumulated number Z of particles is transferred from the I-th to the (I + 1)-th cells during time interval [0, t]. The time evolution of this probability is ruled by the following master equation, We then define the moment generating function where η ≡ e −λ with λ being the counting parameter for the particle transfers Z. From the master equation (49), we can easily obtain the following first-order partial differential equation, which, in vectorial notations, can be written in the following form, where and h ≡ −kN L − kN R .
The identity matrix in P L is of dimension I × I, while the identity matrix in P R is (L − I) × (L − I).
Since the projection matrices satisfy the condition P L + P R = I, we thus have From the above related expressions, the cumulant generating function (60) can be written in the following form, with the global transition rates given by W − = Γ 0 · L 0 · P R · L −1 0 · P L · L 0 · 1.
These two global transition rates can be further developed as Inverting the matrix L 0 , we get So, the two global transition rates are calculated as finally giving