Analytical and simulation studies of driven diffusive system with asymmetric heterogeneous interactions

Totally asymmetric simple exclusion process (namely, TASEP) is one of the most vital driven diffusive systems, which depicts stochastic dynamics of self-driven particles unidirectional updating along one-dimensional discrete lattices controlled by hard-core exclusions. Different with pre-existing results, driven diffusive system composed by multiple TASEPs with asymmetric heterogeneous interactions under two-dimensional periodic boundaries is investigated. By using detailed balance principle, particle configurations are extensively studied to obtain universal laws of characteristic order parameters of such stochastic dynamic system. By performing analytical analyses and Monte-Carlo simulations, local densities are found to be monotone increase with global density and spatially homogeneous to site locations. Oppositely, local currents are found to be non-monotonically increasing against global density and proportional to forward rate. Additionally, by calculating different cases of topologies, changing transition rates are found to have greater effects on particle configurations in adjacent subsystems. By intuitively comparing with pre-existing results, the improvement of our work also shows that introducing and considering totally heterogeneous interactions can improve the total current in such multiple TASEPs and optimize the overall transport of such driven-diffusive system. Our research will be helpful to understand microscopic dynamics and non-equilibrium dynamical behaviors of interacting particle systems.


Results
Asymmetric heterogeneous model. The sketch of the model is displayed in Fig. 1. Two-dimensional periodic boundaries and random update rules are applied, which consists of K equal-sized periodic subsystems. Each system size is L. As the periodic boundary, lane + i K is equivalent to i. Additionally, > K 2 is satisfied because each subsystem has two adjacent lanes. Here, a binary parameter τ = … = …  As the update rule, boundary of each subsystem is periodic in forward direction. Thus, dynamics of lane i are spatially independent, which mean that sites of it are equivalent to each other. Moreover, straight density profile of each subsystem also reflects that dynamics of each subsystem are spatially independent. Thus, dynamics of each subsystem are mainly controlled by heterogeneous interactions among adjacent subsystems. In the global view, dynamics of the proposed driven diffusive system are mainly affected by heterogeneous interactions, since they are capable of quantitatively reflecting interactions among adjacent subsystems and will affect the system's local density. Additionally, in the special case ω = ω i u i d , such totally heterogeneous system deteriorates into the partly heterogeneous one 39 . Moreover, in another special case that all transition rates are equal to each other, such totally heterogeneous system evolves into the homogeneous one.
Detailed balance analysis. We consider particle configurations and use τ { } i,j i to depict particle configurations in subsystem i. Similarly, f i is set as the possibility of presence of particles on lattices in i. Because of homogeneity of each subsystem, possibilities can be expressed with the same f i for each site. Then, the system's partition function can be expressed as: i,j i means the probability of occurrence of the situation where those N particles distribute as τ { } i,j . As for subsystem i, each particle has the same property. Therefore, subsystem i can be treated as a homogeneous system. Weight of different τ { } i,j i should be equal with each other as for specific M i . Besides, it's also equal to weight of occurrence of the situation where M i sites are occupied in subsystem i. Thus, following restraint is satisfied: In this circumstance, the weight can be divided into C L M i species. Then, Eq. (2) can be rewritten as: As for a given configuration τ { } i,j , the probability P τ { } i,j satisfies: C C where C and C′ denote the configuration before and after transition respectively. P(C) denotes the probability of occurrence of such configuration C. → ′ C C W( ) indicates the probability of transition from C to C′. Since particles can perform both hopping in each lane or changing into the other lane, four states ′ C 1 , ″ C 1 , ′ C 2 and ″ C 2 are defined. In detail, both ′ C 1 and ″ C 1 indicate the configuration caused by lane-switching. While, both ′ C 2 and ″ C 2 reflect the one caused by hopping in the bulk. Moreover, both ′ C 1 and ′ C 2 mean the configuration after transition, which is generated from the original state C. While, ″ C 1 and ″ C 2 mean the configuration before transition, which will lead to the state C. Therefore, Eq. (5) can be rewritten as: In fact, as for given ′ C 2 and ″ C 2 , are equal to the hopping rate in the bulk. Besides, as the symmetry of topology of the proposed model, the number of ′ C 2 should be equal to the one of ″ C 2 . Moreover, as for a given ″ C 2 , = ″ P C P ( ) (C ) 2 is satisfied. Because the difference between C and ″ C 2 is just the specific location of particles while the number of particles remains unchanged. Thus, following equation can be derived: Furthermore, based on update rule and symmetry of the system, any state generated from C can also lead to C by performing lane-changing behavior, which means that each ′ C 1 satisfies ″ C 1 . Similarly, each ″ C 1 also satisfies ′ C 1 . Thus, both ′ C 1 and ″ C 1 have the same physical meaning, which means that they are equivalent. Thus, the second expression in Eq. (7) can be rewritten as: Here, we consider a specific situation. C is set to satisfy K,j K . Besides, the particle number in i is set to be more than that of any rest lanes, which means that τ = 1 i,j indicate the configuration after transition that a particle in i updating into − i 1 and + i 1, respectively. Furthermore, the intuitive description of detailed balance is displayed in Fig. 3. By performing detailed balance analysis (see Methods), we can obtain: Therefore, by solving linear equations composed by Eq. (11) (see Methods), the density weight can be derived: Analytical and simulation results of characteristic parameters. By using complex analysis like intermediate value theorem etc. (see Methods), analytical results of characteristic order parameters are obtained. In details, global density ρ can be got: are addressed, which are shown in black and red boxes.   where z means root of Eq. (13). As for a specific system, ρ can be determined since L, K and N are preset. Thus, solutions of z can be obtained from Eq. (13). Then, local density ρ i can be got: In fact, ρ is heavily dependent on ρ i . Similarly, local current J i can be given:  Additionally, expectation 〈 〉 n i of particle number can be obtained:  Furthermore, variance 〈 〉 D n i of particle number can be obtained: As for each subsystem, the probability of having particles or no particles at any two positions is the same. According to the law of large numbers, n L i converges to ρ i in probability when L is large enough, which can be expressed as: Here, Eq. (18) is satisfied for arbitrary positive number δ and ε. Finally, effects of fully heterogeneous interactions on totally current and global transport are emphasized. Based on Eq. (15), the maximum value J max of the total current J total can be derived:   The number of total time steps T is set as 10 10 to obtain steady state. Besides, the final 90% of time steps are retained to ensure occurrence of steady state. In details, Fig. 4 shows relationship among ρ, p i and J i . It reveals that J i varies for a given ρ because of heterogeneous interactions among subsystems. However, as the capacity of maximum current of TASEP, J i is not monotonously changing with ρ. Moreover, the value of ρ corresponding to peak value of J i is also varied for each subsystem as such heterogeneous interactions. The maximum current in each subsystem satisfies 0.25p i , which can also be revealed from the constraint (14) and (15). Additionally, the current in each lane increases linearly with forward rate, which can be reflected from Fig. 4(c) and the above constraint. In details, heterogeneous interactions are reflected by asymmetric rates , ω = . 0 384 , ω = . between adjacent subsystems, which can lead to various site occupation rates. Then, different cluster states of particles occur. Thus, varied particle stochastic dynamics of each specific lane are presented, including various densities, diverse flows etc. In order to highlight the effect of transition rates on local current, the scaling rate r is introduced. Here, r denotes a factor that acts on preset transition rates and depicts the extent of heterogeneity of the global system. Firstly, primitive transition rates are set. Here, primitive transition rates denote the target that the scaling rate r acts on, which are preset and randomly generated. Three complete kinds of such driven diffusive system can be constructed, namely the case of totally heterogeneous interactions, partly heterogeneous ones and homogeneous one. Then, r acts on the primitive transition rate, which makes upward and downward transition rates become r and 1/r times of original values respectively. Figure 5 shows the relationship among J i , r and p i . It can be found that different channels exhibit various properties. Local current reaches the extremum with different scaling rates and contains disparate number of peaks. Since particles in the system have two degrees of freedom, the effect of forward motion also needs to be emphasized. Compared with Fig. 4, the relationship between J i and p i in Fig. 5 is similar to that reflected in Fig. 4 in terms of quantitative and qualitative evolution rules. This is because that configuration of particles in each subsystem can be determined for preset transition rates. Then, the bulk current is heavily dependent on unidirectional self-driven motions of particles. Thus, in this circumstance, the local current is proportional to forward rate.  Besides, Figs 6 and 7 display relationship among ρ, ρ i and r. In Fig. 6, r affects each channel differently. Primitive transition rates are set as identical constants. ρ i is found to evolve in different ways as various effects of r. Particularly, ρ i is the same when = r 1, since r has no effect on transition rates in this circumstance. Moreover, as the influence of r, ρ i will be different with increasing ρ. Additionally, as conserved particle number, ρ i also increases monotonously with ρ. Especially, ρ i becomes 0 when ρ = 0, as there're no particles. Similarly, as the full capacity, ρ i reaches 1 when ρ = 1. While, the case of randomly generated primitive transition rates is shown in Fig. 7 with the same effect of r on each subsystem. Compared with Fig. 6, similar evolution law of density is presented in the form of a spindle although different values of ρ i appear. It also reveals that all analytical results match well with Monte-Carlo simulations.
Then, corresponding density profiles ρ (x) i are displayed in Fig. 8.
is independent of the location x of sites as the periodic boundary condition and update rule of the , ω = . 0 431 , ω ω = .
= . 0 285, 0359  proposed system, which means that any particle in i is equivalent to each other. In other words, the interaction among particles in i is homogeneous. Moreover, Fig. 9 shows 〈 〉 n i and i are disparate as different ω i u and ω i d . Moreover, we aim at analyzing the influence of disturbance on the system coupled with changing transition rates. Thus, for simplicity, the upward transition rate ω u 1 of lane 1 is set to range from 0 to 1 here, while other transition rates except for ω u 1 remain unchanged. Corresponding transition rates are randomly generated to confirm validity of universality of calculations. More complex topologies (namely, = K 10, 20, 30, 50) rather than just previous four TASEPs 39 are calculated. Relationship among K, ρ i and ω u 1 is displayed in Fig. 10. As for Fig. 10(a), it describes = K 10, where randomly generated transition rates except for ω 1 u range from 0 to 1. Ten colorful lines show evolvements of ρ with increasing ω 1 u . Especially, green line displays lane 10, which reveals that ρ 10 increases with ω 1 u and the gradient of change of it gradually decreases, since increasing ω 1 u directly leads to the increase of particle number in lane 10. Besides, black line shows lane 1, which shows that ρ 1 decreases with ω 1 u , since increasing ω 1 u directly leads to the decrease of particle number in lane 1. Similarly, Fig. 10(b) shows = K 20. Increasing ω 1 u directly leads to decreasing ρ 1 and increasing ρ 20 . Besides, densities of other remained lanes also become varied. Moreover, Fig. 10(c) depicts = K 30. Similarly, increasing ω 1 u directly leads to decreasing ρ 1 and increasing ρ 30 , which also causes changes of other subsystems. Furthermore, Fig. 10(d) illustrates = K 50. Similarly, increasing ω 1 u directly causes decreasing ρ 1 and increasing ρ 30 , which also leads to changes of other channels.
Thus, based on Fig. 10, it reveals that no matter how system's topology changes, the disturbance will spread and affect local densities of other channels, when transition rate of one channel changes. The first and K-th channels are the most affected since ω 1 u changes. Increasing ω 1 u leads to continuously decreasing ρ 1 and increasing ρ K . Additionally, except for lanes 1 and K, the disturbance will have greater impacts on lanes 2 and − (K 1) because of stronger interactions among them. Here, it should be addressed that as changing ω u 1 is set to illustrate the influence of disturbance, increasing ω u 1 will directly lead to decreasing ρ 1 and increasing ρ K . Thus, lanes 1 and K should be emphasized. Finally, Fig. 11 displays the relationship among the absolute value ∆J of increment of total current, K and ω u 1 . It reflects that ∆J monotonically decreases for any channel number with increasing ω u 1 . When ω u 1 degrades the system from heterogeneous one to homogeneous one, ∆J becomes 0. Moreover, total current increases monotonously with increasing channel number.
Additionally, besides in technical perspective, the improvement of our work compared with pre-existing results is intuitively presented through Figs 12 and 13. In details, beyond pre-existing results, the improvement of our work can be intuitively found in Fig. 12, through studying the relationship between ρ and J total in three complete kinds of such driven diffusive system, namely totally heterogeneous interactions (namely, our work), partly heterogeneous ones and homogeneous ones. Under any conserved global density, the optimal total current obtained from our work is generally greater than that of homogeneous one or partly heterogeneous ones including pre-existing results 39,40 and the universal case of partly heterogeneous. Besides, the value of optimal total current obtained from our model is equal to that of homogeneous when ρ = . 0 5, since the optimal current in each channel reaches to the maximum. Here, the universal case of partly heterogeneous interactions in Fig. 12 is built by equal transition rates for some randomly chosen channels (namely, equal transition rates χ i for pre-selected channel i and χ ≠ χ i j for ∀i ≠ j) and randomly generated rates for remaining ones. Thus, improving pre-existing results, our work presents a method of optimizing overall transport of such driven-diffusive systems in the global perspective under conserved global density, since the optimal total current obtained from our work introducing and considering totally heterogeneous interactions is generally greater than that of homogeneous one or partly heterogeneous ones. Here, 100 TASEPs are considered in calculations to obtain enough large and universal driven-diffusive system to confirm the validity of the improvement. Moreover, matching well with analytical results, Monte-Carlo simulations are also performed to make sure the validity of the improvement of our work comparing with pre-existing results.
Moreover, in order to intuitively describe the improvement of our work compared with pre-existing results, the effect of introducing and considering heterogeneous interactions on overall transport is further studied by calculating the relationship among J total , ρ and r shown in Fig. 13. Comparisons of our work with the universal case of partly heterogeneous are presented in Fig. 13(a,b), since totally heterogeneous case will evolve into the universal case of partly heterogeneous when = r 1. While, comparisons of our work with homogeneous case are  Fig. 13(c,d), where totally heterogeneous case can evolve into homogeneous one when = r 1. Based on Fig. 13(a,c), it can be found that J total increases at first and then decreases later with increasing ρ. However, based on Fig. 13(b,d), it reveals that the optimal current of our work is much higher than that under partially heterogeneous interactions and homogeneous one. Thus, our work also reveals that introducing and considering totally heterogeneous interactions can increase the optimal total current in overall transport of self-driven particles in such driven diffusive system modelled by multiple TASEPs, which can also reflect that totally heterogeneous interactions can improve the overall transport of such interacting multi-body particle systems.
Furthermore, based on Figs 4, 5, 6, 7, 8, 9, 10, 11, 12 and 13, it should be pointed out that heterogeneous multilane TASEPs proposed here can also be equivalent to driven diffusive system coupled with Langmuir dynamics of detachment rate 1/f i and attachment rate z 12 . This is due to following reasons. Firstly, actually, these concerned multiple TASEPs can be treated as being in contact with the thermal bath that washes out any density-density correlation whereas in the case under scrutiny the global density is conserved. Then, in the view of statistical mechanics, proposed system can be treated as a grand canonical ensemble. Because each channel of our model can exchange particles with adjacent subsystems with conserved global density, which can fit well with the essence of grand canonical ensemble, namely, each system can exchange energy and particles with other systems with conserved overall potential. Hereafter, ref. 12 reported that at high kinetic rates (namely, high values of attachment rate ω A and detachment rate ω D ) the bulk density in TASEP was structureless and equal to the Langmuir equilibrium density That's to say, when high values of the attachment and detachment rates are set 12 , bulk dynamics will dominate in the competition with boundary dynamics, which means that the model can degenerate into the periodic boundary condition. Thus, based on Langmuir equilibrium density , Eqs (12) and (14), ρ i can be derived as ρ = = . Therefore, our model can be mapped into the one of ref. 12

Discussion
To summarize, we propose a driven diffusive system composed of K TASEPs. Different from previous work, interactions among adjacent subsystems is asymmetric heterogeneous. Generally, the whole system is coupled with periodic boundaries in a two-dimensional periodic torus. Self-driven particles can switch into adjacent channels or unidirectionally move. System's stochastic dynamics are mainly controlled by transition rates rather than hopping rate due to spatially homogeneous local densities. Based on detailed balance principle, we check master equation. Then, restraint about weight factor f i is obtained , which proves that our model can spontaneously evolve into the system coupled with symmetric transition rates. That's to say, we have generalized previous research work 39 . Moreover, by applying complex analysis, analytical results of characteristic order parameters are obtained, including ρ i , J i , ρ, 〈 〉 n i , 〈 〉 D n i and probability of configuration n L i . Besides, universal characteristics are revealed by both analytical solutions and Monte-Carlo simulations which match well with each other. Local density is found to monotonously increase with global density. Due to totally heterogeneous interactions, specific local densities are different from each other. Based on homogeneity on spatial site locations, local density profiles are spatially homogeneous, which are reflected by horizontal lines. Moreover, due to capacity of the maximum flow in TASEP, the relationship between local current and global density doesn't monotonously change, which appears different from that of densities. While, local current is proportional to it at given transition rates. Additionally, by calculating more complex topologies (namely, K = 10, 20, 30, 50) of such driven diffusive system rather than just previous four TASEPs 39 , it reveals that no matter how system's topology changes, the disturbance of changing transition rate of one subsystem will spread and affect local densities of other subsystems. As asymmetric heterogeneous interactions among adjacent subsystems, changing interaction rate of any subsystem will directly lead to variations of densities of its adjacent channels, which also means a larger change in particle configurations in adjacent ones. Because the changing of spatially dependent local density is equivalent to changing of dynamics of particle configurations, since local density is the essentially statistical average value of occupation probability of sites in each subsystem. Moreover, as each lattice position with or without particles are independent events, expectation is equal to system size multiplied by local density. Similarly, variance is equal to system size multiplied by ρ − ρ (1 ) i i . Thus, both expectation and variance depend on global density. With increasing global density, expectation monotonously increases, while variance first increases and then decreases.
Finally, in order to detailly describe the improvement of our work, we firstly interpret the previous work in this area and then detailly compare results of our work with pre-existing physical systems. Improving pre-existing results, we propose a more universal and realistic interacting multi-body particle system which is more reliable to depict real transport phenomena. Thus, we focus on constructing such driven-diffusive system by employing multi-channel TASEPs with totally heterogeneous interactions. Besides, we obtain explicit analytical solutions and Monte-Carlo simulations of such proposed driven-diffusive system to avoid using previous approximation methods that often lead to results lack of universal laws. Moreover, different with previous work, we get universal laws of characteristic order parameters of such proposed driven diffusive system by considering complex topological structures. Furthermore, we also point out that introducing and considering totally heterogeneous interactions in such driven-diffusive system can increase the optimal total current in overall transport of the global driven diffusive system, which can also reflect that totally heterogeneous interactions can improve the overall transport of such interacting multi-body particle systems. Besides, based on our research work, it can also be revealed that introducing and considering totally and fully heterogeneous interactions is a way to optimize the total current in such driven diffusive system modelled by multiple TASEPs under the circumstance of conserved global density, which can be proved by calculating the relationship between global density and total current in three complete kinds of situations (namely, our work, partly heterogeneous interactions and homogeneous), relationships among J total , ρ and r in the intuitively three-dimensional view and the relationship between J total and r with fixed value of ρ in intuitively two-dimensional view. Additionally, we also point out totally heterogeneous interaction rates can improve the total current in such multiple TASEP system and optimize the overall transport of such driven-diffusive system. Our research will be helpful to study other homeomorphism systems like coupling with memory reservoirs, which will be helpful for better understanding of stochastic particle dynamics in driven diffusive systems constituted by multiple TASEPs. Other performing work will be reported later.

Methods
Detailed balance equation analysis. For simplicity, we use the symbol A to denote . Based on Eqs (9) and (10), detailed balance equation can be derived: is obtained, which is Eq. (11).

Solving linear equations composed by Equation (11).
In fact, Eq. (11) is satisfied for arbitrary ∈ … i K {1, 2, , }. Thus, K similar equations like Eq. (11) can be obtained. Besides, these K equations constitute the following linear equations: are satisfied. Then, the rank of the matrix W needs to be investigated in order to solve Eq. (22). Additionally, the cofactor of W can be derived: Here, W 11 denotes one of the cofactors, which corresponds to the element located in the first row and the first column of W. Thus, ≠ W 0 11 is satisfied. Moreover, W is not full rank. Therefore, the rank of W is K − 1. That's to say, the solution of Eq. (22) is one-dimensional.
Then, in order to solve Eq. (22), we suppose that following equation is satisfied for arbitrary ∈ … i K {1, 2, , }: should be proved to be the solution of Eq. (22). As the symmetry of topology and updating rule, following constraints are satisfied: , where Const means a const. Actually, according to the definition of f i , there'll be no influence on the system if all weight factors increase or reduce the same multiplier at the same time. For simplicity, Const is set as one here. Then, according to Eq. (24), density weight can be expressed as shown in Eq. (12). where generalized partition function Z s ( )

Complex analysis. As for a generalized function
Here, F(s) denotes a generalized function.  where Z s ( ) based on Eqs (12), (29) and (30), the local density ρ i can be derived: