Exactly solvable deterministic lattice model of crossover between ballistic and diffusive transport

We discuss a simple deterministic lattice gas of locally interacting charged particles, for which we show coexistence of ballistic and diffusive transport. Both, the ballistic and the diffusive transport coefficients, specifically the Drude weight and the diffusion constant, respectively, are analytically computed for particular set of generalised Gibbs states and may independently vanish for appropriate values of thermodynamic parameters. Moreover, our analysis, based on explicit construction of the matrix representation of time-automorphism in a suitable basis of the algebra of local observables, allows for an exact computation of the dynamic structure factor and closed form solution of the inhomogeneous quench problem.


Introduction
One of the main challenges of nonequilibrium statistical mechanics is a rigorous derivation, without any assumptions or approximations, of irreversible macroscopic transport laws, from the microscopic reversible equations of motion. This goal can be achieved only in certain specific interacting modes, see e.g. a very recent work [1], which can be considered as a companion to the present text. The aim of this article is to explore dynamical and transport properties of the model of a reversible cellular automaton introduced and preliminary studied in [2], consisting of locally interacting positively and negatively charged particles and freely propagating vacancies. This work is therefore connected to two active areas of research. First of all it can be viewed as a discretetime (paralel-update) and deterministic version of the exclusion processes [3][4][5], which have been widely studied [6,7]. In those models the particles obey exclusion principle, i.e. a particle can move to a neighboring site only if the site is unoccupied. Additionally, our model can be related to the gas of hard rods, which has been extensively studied [8][9][10][11][12][13][14], in a specific simple, yet nontrivial regime where the dynamics can be mapped to a spatio-temporal lattice. Due to the model's simplicity many essential questions regarding the dynamics and the transport can be answered explicitly. We hope that https://doi.org/10.1088/1742-5468/aae853 J. Stat. Mech. (2018) 123202 these results will illuminate a general understanding of dynamics in a class of similar or related interacting models.
Furthermore, our results can be put into a more modern perspective of the generalized hydrodynamics, which was developed recently [15][16][17][18][19][20] and provides an exact semi-classical description of the long time dynamics of integrable quantum systems. Semi-classical dynamics is connected to the generalization of the hard rod gas [21]. Taking into account some approximations, quantum systems can be described in terms of multi-species hard-core colliding particles [22][23][24][25][26], exactly the type of model studied here.
In the article we consider a deterministic dynamics describing the elastic scattering of charged particles, as well as a Markovian generalization, where the particles exchange their position with a certain probability upon interaction. We provide analytical expressions for dynamical (transport) coecients (specifically, the diusion constant and the Drude weight) and compare our results to the eective (hydrodynamic) description, which is typically used without rigorous justification.
In the first section we introduce the model and the underlying mathematical structure. Our model supports solitonic excitations. Characterizing the space of solitonic observables we are able to construct the set of local conserved quantities, which describe the stationary states of our model in terms of the corresponding generalized Gibbs ensemble.
The second section deals with the exact calculation of transport coecients. In the first part we introduce the linear response definitions of the charge diusion constant and the Drude weight. This is followed by the calculation of the lower bound on Drude weight by employing the Mazur inequality. The central part of the section comprises the analytical calculation of the linear response diusion constant and the Drude weight.
In the third section, a large time solution of the inhomogeneous initial value problem is obtained, showing ballistic propagation of the step function profile with diusive corrections. Depending on the density of vacancies and the imbalance of particles the system is shown to exhibit ballistic, normal or insulating behavior. The results following from the inhomogeneous initial state are compared to the hydrodynamical picture, showing perfect agreement.
In the fourth section we show how the spatio-temporal correlation functions can be obtained for particular stationary product states. The dynamics corresponds to the diusive broadening of the central peak, and the free propagation of solitary excitations.

The model
The model is defined on the chain (one-dimensional periodic lattice) of even size n, where each lattice site can be occupied by three types of particles: positively charged particles (+), negatively charged particles (−) or vacancies (∅). The configuration of particles at time t is denoted by s t , with s t = (s t 1 , s t 2 , ..., s t n ) and s t x ∈ {∅, +, −}. The dynamics of particles is described by the propagation rule https://doi.org/10.1088/1742-5468/aae853 that correspond to the following local mapping These rules describe elastic scattering of particles, with scattering in the first half-time step occurring between the particles on odd-even sites and in the second half-time step between particles on even-odd sites. The schematic representation of the dynamics can be seen in figure 1.
The density of charge q t x is defined as a sum of net charge on two neighboring sites, is occupied by a positive charge, −1 if a negative charge occupies the site and 0 in the case of an empty site. One should note that the total charge, is a constant of motion. To study the dynamics of charges q x , we introduce the corresponding current j x j t+1/2 (2.6) The local current j x and the charge q x satisfy the continuity equation (see appendix A for the details with appropriate precise notation).

Algebra of observables
To describe the statistical properties of the model we introduce a multiplicative commutative algebra of observables Ã R 3 , i.e. functions over (2.7) The algebra can be extended to a local algebra of functions over the lattice configuration space (Z 3 ) ⊗n by defining the local basis https://doi.org/10.1088/1742-5468/aae853 We introduce the compact notation for the local basis elements of the tensor product algebra A =Ã ⊗n , The time propagation of observables can be represented in terms of a linear map U ∈ End(A), a t (s) = a(s t ) ≡ U t a(s), (2.10) which is again composed of the local two site propagators, The expectation value (2.12) can be represented as a p . For later convenience we introduce the following two bases for Ã , depending on the density parameter ρ, , , ( 2.14) which are dual w.r.t. the maximum entropy state, The observable [0] corresponds to the identity in Ã . Let ρ = (ρ 1 + ρ 2 )/2 be an average value of the density parameter on sites 1 and 2 of a unit cell, comprising two consecutive sites, while ∆ = (ρ 1 − ρ 2 )/2 is one half of their dierence. With respect to this parametrization, the local propagator takes the following form

16)
where ρ = 1 − ρ is a density of vacancies. Note that the propagator is expressed w.r.t. the basis and the dual basis that have dierent values of parameter Δ, corresponding to the exchange of ρ 1 and ρ 2 . To be more precise, the elements of the propagator are defined by

17)
with l and k expressed in bases (2.14) with density ρ 1 , and l , k in bases with density ρ 2 . In these bases the total charge (Q) and the total current (J ) observables are expressed as

Solitons and local conservation laws
Local and pseudolocal conservation laws play a key role in the physics of local correlation functions [2,[27][28][29][30]. Here we show how the ballistic propagation of vacancies gives rise to exponentially many local conservation laws, which should be contrasted to expected behaviour in generic integrable systems, where the number of conservation laws scales linearly with the system size. First, notice that the local propagation of [20] and The solitons ensure the existence of exponentially many conserved quantities, which we now construct. Let us define the even solitonic subspace as and the odd solitonic subspace as For the elements of these subspaces, the following holds,

22)
where k = (k 1 , k 2 , . . . , k l ) denotes the places where the observables [2] appear. Therefore any solitonic observable s k ∈ S e/o , is a local density of a conserved quantity S k = x η 2x (s k ), If the length of the chain is n, the number of distinct solitonic conserved quantities S k is 2 n/2 , i.e. square-root of the number of all linearly independent observables. As already mentioned, the total charge Q is conserved as well. Therefore the number of all local integrals of motion is 2 n/2 + 1.

Generalized Gibbs ensemble
A generic integrable system equilibrates to the non-thermal equilibrium state, which is expected to be described by the generalized Gibbs ensemble [31][32][33]. Let us consider a homogeneous initial value problem with a translationally invariant initial state ρ i , Provided that the complete set of local charges is {S e k , S o k , Q}, the hypothetic equilibrium ensembles can be described by Here the chemical potentials β, β e/o k are obtained from the condition that the initial state expectation values of conserved charges match the ones in the equilibrium state ρ GGE . The expression for ρ GGE can be further simplified, where we introduced the product equilibrium distribution p, Note that here the dual basis vectors pertaining to even and odd sites correspond to ρ 1 for even and to ρ 2 for odd sites, however, we suppressed the explicit dependence to preserve the compactness of notation. The first two terms in equation (2.26) are obtained by expanding the exponents of even and odd charges, and noting that the exponents of even/odd solitonic observables are also even/odd solitonic observables. The expression (2.27) is obtained by expanding exp(−α e S e 1 − α o S o 1 − βQ), and choosing α e and α o such that the expansion of the exponent does not contain basis vectors [2] . This imposes the following restriction, For the purpose of hydrodynamic applications we are interested only in expectation values of the observables from the set J = {s k , j}, which includes densities of conserved quantities and their currents. The only current that is not a linear combination of charges is the charge current j. On the subspace of charges and currents, the state (2.26) can be equivalently represented as where we introduced symmetrized and anti-symmetrized charges S ± k = S k ± η 1 (S k ). Here, we take into account that, for every a ∈ J , the expectation value a S e k S o k p vanishes, implying a ρ GGE = a ρ GGE .
In the rest of the article, the results are obtained for the factorizable (separable) steady states ρ GGE = p (2.27), i.e. considering c ± k = 0. The overlaps of observables w.r.t. states p are easy to express due to the following relations, https://doi.org/10.1088/1742-5468/aae853

Linear response
One of the outstanding questions when regarding the transport phenomena is how to derive the phenomenological transport laws, e.g. Fick's law, stating that the current is proportional to the gradient of the external field ∇h, where the proportionality constant σ is the conductivity. The transport coecients in the linear response regime can be calculated in terms of the stationary state timecorrelation functions. The diusion constant D is related to the conductivity through the Einstein's relation,

2)
where χ is the static susceptibility. The conductivity can be expressed as the time integral of the current autocorrelation function, The Drude weight, defined as is the rate at which the current in the system increases when the system is exposed to the constant gradient external field. For precise definitions and derivation see appendix A. If the Drude weight is non-zero, the diusion constant diverges. In this case it is convenient to regularize it by subtracting the divergent part,

Mazur bound on Drude weight
In this subsection, we bound the Drude weight by local conservation laws, generalizing Mazur's argument [29] The following inequality holds for any stationary probability distribution p, Mazur's inequality then follows directly by setting a = 1 where {O k } is a set of conserved charges. Inserting ā into the expression (3.6), and maximizing the expression w.r.t. the set {α k } one obtains the following lower bound, Here we assumed that the set of the charges is orthogonal, Let us now proceed to the calculation of the Mazur bound. The only two solitonic charges that have nonzero overlap with the current are Additionally, the following linear combination of Q and S e/o (1) also contributes to the Mazur bound, Therefore we are able to utilize Mazur's inequality (3.7), and bound the current autocorrelation function, yielding the following lower bound on the Drude weight, In the following subsection we show that the lower bound saturates the exact result.

Exact results on time decay of current-current autocorrelation
The main property of the time-dependence of the autocorrelation functions is that time propagation can be restricted to a particular subspace, where the propagator has at most five-diagonal form. This allows for an explicit calculation of C J (t).
In order to calculate the autocorrelation we first observe that the number of local basis elements [1], occurring in the time propagated observables, is conserved. Furthermore, the basis elements [2] propagate ballistically, and due to the orthogonality (2.30) we can restrict the computation of current-current autocorrelation function, https://doi.org/10.1088/1742-5468/aae853 Formally, the argument can be recast as follows. We can consider only observables with a single [1], i.e. those spanned by any local basis vectors [α 1 α 2 · · · α r ] x with a single α i = 1, forming an invariant subspace A (1) , U A (1) ⊆ A (1) . Let P be a linear projector P : A (1) → A J , P 2 = P. Then for every a ∈ A (1) , J(1 − P )a p = 0, and A (1) /A J is invariant under ηU e and η −1 U e This implies that in order to compute the current autocorrelation function the dynamics can be restricted to A J . Specifically, by defining the reduced half-time step propagator U ,

13)
it is possible to express the current-current autocorrelation function exactly as Note that if we consider an initial state with ρ 1 = ρ 2 , the basis vectors are position dependent, according to (2.14). Initially the basis on even sites corresponds to the density of particles ρ 1 and on odd sites to the density ρ 2 , however after every half-time step the two bases get exchanged. The reduced propagator reads so that the autocorrelation function (3.14) reduces to Applying the matrix U to the left we obtain an explicit expression for the time dependent correlation function

18)
If ∆ = 0 and µ = 0, i.e. for a translationally invariant state without the charge imbalance, the Drude weight vanishes and the transport is diusive with the following diusion constant and conductivity, Otherwise the transport is ballistic, and the Drude weight reads which is equal to the Mazur bound (3.11). Note that vanishing of Drude weight D = 0 implies that both µ = 0 and ∆ = 0, except for two additional cases, namely (i) ρ = 1, ∆ = 0, and (ii) ρ = |∆| = |µ| = 1/2, corresponding to trivial dynamics.

Stochastic generalization
Similar calculation can be repeated for a stochastic generalization, which corresponds to the tunneling of the particles, by allowing two additional processes (3.21) Let Γ and Γ denote the tunneling and scattering probabilities respectively, namely (±, ∓) maps to (∓, ±) with probabiliy Γ, and to (±, ∓) with probability Γ = 1 − Γ. In this case the local propagator (2.17) in the basis (2.14) reads Due to the change of dynamics which is no longer deterministic, the current has to be redefined in order for the continuity equation to hold (see appendix A for the details), Most of the discussion corresponding to the deterministic dynamics still applies and results in an explicit expression for the time-dependent correlation function, Similarly as in the deterministic case, the transport is ballistic if µ = 0 or ∆ = 0, with the same Drude weight, while in the regime µ = ∆ = 0, the transport is diusive. Setting µ = ∆ = 0 in (3.24) and evaluating the sum (A.13) yields the following expression for the diusion constant

25)
where ρ is the rescaled density, ρ =Γρ. In the deterministic limit Γ → 0, the diusion constant of the deterministic model is recovered, while for nonzero values of the scattering probability, D scales as Γ ρ −1 − 1 with polynomial corrections.

Inhomogeneous quench
In this section we consider an inhomogeneous quench problem, where the initial probability distribution is given by the product state with the density of particles ρ 1 on odd sites, and ρ 2 on even sites. Additionally we set a fixed expectation value of the charge q to µ L/R 1,2 , with the superscript indices denoting the left and the right half of the chain, while the indices 1 and 2 correspond to odd and even sites respectively. Note that we are considering the cases where the system size is divisible by 4, and the local propagation is initially applied to the sites (1, 2), (3, 4) The objective is to compute the steady state profile of the charge q, Since the number of local states [1] is preserved, the state p can be linearized, Similarly as before we can consider the subspace with a single local state [1] . The halfstep propagation on this subspace is

4)
where Γ = 1 − Γ. Note that after the half-time step the change of basis occured (i.e. ρ 1 ↔ ρ 2 ). Let us introduce the basis of the linear space spanned by the local charge densities,ê In this basis the full time step takes the form with the 2 × 2 blocks A, B, C given by while the initial state can be expressed as In what follows, we restrict the discussion to the deterministic case (Γ = 0, Γ = 1), since the results for the stochastic generalization can be obtained by rescaling ρ 1,2 → Γρ 1,2 . The matrix U can be diagonalized using the block Fourier transform, namely writing the eigenvector with eigenvalue λ(k) in the form r e ikr v(k), (4.9) where v is a two component vector depending on the momentum variable k ∈ [−π, π) (which becomes continuous in the thermodynamic limit n → ∞). The eigenvalue problem reduces to the following 2 × 2 matrix problem The eigenvalues and eigenvectors read (4.11) In order to obtain the full time evolution we express a part of the initial state µ 2l−1ê2l−1 + µ 2lê2l in terms of the eigenvectors, From this expression we can calculate the constants α 1,2 (k), yielding a complete timedependent profile f (x, t) (4.2) in an integral form, where α 1,2 (k) = α 1,2 (k) [1,1] · v 1,2 (k) . Let us now focus on the asymptotic shape of the profile, t → ∞. In this limit we can consider only the contribution of the leading eigenvalue, λ 1 (k), since |λ 2 (k)| < 1. Furthermore, since λ 1 (k) is an analytic function of k in the vicinity of k = 0 and λ 1 (0) = 1, |λ 1 (k = 0)| < 1, we should take into account only the contributions at k ≈ 0. In this region the leading eigenvalue can be approximated by which implies In the long time limit, the steady states can form on dierent space/time scales around the junction, depending on the type of the transport. In particular, if the transport is ballistic, the steady states arises on the light rays, v = x/t. However, in the case of diusive transport the dynamics is more localized, therefore one should consider the steady state formation along the space-time coordinates u = x/ √ t. In general the transport in our model is ballistic, therefore the steady state profile depends on the ballistic coordinate, The only non-vanishing contribution steams from y in the vicinity of the ballistic line, If u is kept finite as t → ∞ and v − γ 1 = 0, one obtains (4.18) where we introduced a simpler notation κ L/R = µ . Introducing a new variable h = k √ t and taking the limit t → ∞, the integration is reduced to the Gaussian integral. Noticing that α 1 (0) = 1 2π , we obtain the following expression, Approximating the sum by an integral and evaluating it yields the following result, The asymptotic charge profile on ballistic time scales, f (v), is a step function that moves with the velocity γ 1 .
Let us now consider the diusive region around the ballistic front Making similar simplifications as before, yields and after the identification of the sums with integrals we get the final result The solution of the diusion equation, therefore we can read out the diusion constant from the expression (4.23). Since the coordinate x corresponds to two lattice sites, the diusion constant should be rescaled, (4.25)

The hydrodynamical picture
Here we derive the velocity of the front γ 1 using the arguments of generalized hydrodynamics. The basic idea is the following: considering two half chains prepared in distinct stationary homogeneous states joined at the origin. The hydrodynamical approach assumes that the non-equilibrium steady state arises along the light-cone coordinates centred at the origin. The system is assumed to reach generalized equilibrium ρ GGE (v) on a given light-ray v = x/t [15,16]. Note that this is not surprising, since in the limit t → ∞, the subsystem between the light-rays v ± ε is infinitely large and is expected to equilibrate. The remaining slow modes characterizing the NESS are the conserved charges and their currents. We wish to compute the profiles of the charge q and the corresponding current j. Their expectation values in the GGE are where the constants c 0 , c + 1 , c − 1 are determined by fixing expectation values of the charges, The only coecient that depends on the light-ray is μ defined as µ 1 = µ ρ 1 , µ 2 = µ ρ 2 .
Inserting expectation values (4.26) and (4.27) into the continuity equation, yields the following equation for the chemical potential µ Taking into account the initial condition we obtain the solution for µ, which is in perfect accordance with (4.20).

Dynamic structure factor
Finally, we calculate the full spatio-temporal density-density correlation function C q (x, t) = q t 0 q x p . For simplicity, we restrict the discussion to the deterministic, translationally invariant case (i.e. ρ 1 = ρ 2 and Γ = 0). Similarly as before, only the observables with at most one local [2] contribute to the overlap q t 0 q x , therefore the time evolution of q 0 can be restricted to the infinite family of subalgebras A [z] , z ∈ Z, with the basis elements a z x defined as a e 0 The reduced full time step propagator U has the following matrix elements, The dual vectors (a d x ) are obtained from expressions (5.2) simply by replacing the canonical basis vectors by the corresponding dual (primed) vectors. In this basis, the dynamical structure factor C q (x,t) can be expressed as The , are mostly zero, due to the following property of the reduced time propagator U :

(5.7)
Furthermore, the submatrices U d d have block three-diagonal structure, with the block dimensions 1 × 1, 1 × 2 and 2 × 2 (for explicit expression see appendix B), therefore it is possible to use a similar approach as in section 4. Introducing the Fourier basis, each of the infinite submatrices U d d is reduced to a finite matrix, which depends on the parameter k (see appendix B for the details). In the basis [1] (k) g [2] (k) . . . , (5.9) the time propagator takes the following form, with the matrices A and B given by In the Fourier basis (5.4), the dynamical structure factor corresponds to The calculation can be split into a ρ-dependent diusive part and a ballistic contribution proportional to µ 2 , (5.14) The diusive contribution can be evaluated rather easily, since O (D) = Q 0 , where the matrix U (D) corresponds to the central 2 × 2 block of the reduced propagator U (5.10), In the large time limit, the contribution corresponds to the normal distribution, In order to obtain the ballistic contribution, we first observe that for any constants c e/o ∈ C, the following holds and the matrix Diagonalizing the transposed matrix U (B)T and calculating the relevant overlaps, we obtain In the large time limit, the first integral reduces to δ x,t + δ x,−t , while the second one correspond to the normal distribution (5.15). Taking into account both of the contributions, we finally obtain the asymptotic shape of the structure factor, The profile consists of the central peak that spreads diusively, and two spikes that propagate ballistically, see figure 2.

Conclusion
In the paper we explored the transport properties of a simple reversible cellular automaton modelling a hard-point interacting gas of charged particles. Depending on the state, the model exhibits diusive, ballistic and insulating behaviour. If the state is dimerized (i.e. the particle densities on odd and even sites dier) or in the case of charge imbalance, the transport is ballistic with diusive corrections. If the state is translationally invariant and the densities of both species of particles coincide, the ballistic contribution vanishes and the transport is purely diusive. When the density of particles equals 1, the dynamics is frozen and the system behaves as an insulator.
We constructed the set of conserved quantities and the corresponding set of stationary states described by the generalized Gibbs ensemble. By obtaining the explicit expressions for the time dependence of autocorrelation functions we were able to calculate the Drude weight and the diusion constant. The Drude weight was shown to match the Mazur lower bound perfectly. We also analytically solved the inhomogeneous initial state problem, which corresponds to the step function charge density profile moving with a constant velocity and the diusive corrections at large times. The velocity of the step function profile was shown to match the velocity obtained from the hydrodynamic consideration. Furthermore, we calculated the structure factor exactly. The asymptotic spatio-temporal correlation profile consists of two ballistically moving δ-spikes and a diusively broadening central peak.
Our model can be used as a rigorous benchmark of the physical properties of locally interacting lattice systems with conservation laws, showing agreement with the proposed eective descriptions such as hydrodynamics or the Mazur inequality, which are typically used without a rigorous justification. The main question that remains is whether a similar description of complete time evolution can be obtained for more complicated interacting classical or quantum systems, say for a typical integrable systems. https://doi.org/10.1088/1742-5468/aae853

J. Stat. Mech. (2018) 123202
In general, however, the connection between the odd and the even currents can be deduced directly from the relations (A.4) and (A.3), and reads The linear response quench setup is consistent with the one considered in [2]. The system in the initial stationary state p is kicked out of equilibrium at time t = 0, by a weak constant external field K(h)-a constant force h: (A.7) We note that for convenience we (only here) label our spatial lattice from x = −n to x = n, thus consisting of 2n + 1 sites and considering open boundary conditions. At time t, the current induced by the force takes the following form j LR (t) = 1 2 ( j Due to the staggered structure, the current is averaged over two lattice sites and because of the discreteness of time, it is additionally shifted for half of the time step (see (A.1)). To make the notation more compact we introduced the following convention Note that we subtracted the initial value of the current, since we are not interested in the current already present in the initial state, but rather in the current that is induced on top of the initial state current due to the perturbation. Writing out the propagator explicitly,

Appendix B. The spatio-temporal correlation function
Here we present some of the details of calculation that were omitted in section 5. All the nonvanishing submatrices U d d (5.3) are block 3-diagonal. The blocks of the submatrix U 00 are of the size 2 × 2 and the blocks in U d 0 for d = 0 have the dimension 2 × 1,