Driven Lorentz model in discrete time
Dan Shafir Alessio Squarcini Stanislav Burov Thomas Franosch

APS/123-QED



dansh5d@gmail.com
Physics Department, Bar-Ilan University, Ramat Gan 5290002, Israel

alessio.squarcini@uibk.ac.at
Institut für Theoretische Physik, Universität Innsbruck, Technikerstraße 21A, A-6020 Innsbruck, Austria

stasbur@gmail.com
Physics Department, Bar-Ilan University, Ramat Gan 5290002, Israel
 
thomas.franosch@uibk.ac.at
Institut für Theoretische Physik, Universität Innsbruck, Technikerstraße 21A, A-6020 Innsbruck, Austria



§ ABSTRACT 
We consider a tracer particle performing a random walk on a two-dimensional lattice in the presence of immobile hard obstacles. Starting from equilibrium, a constant force pulling on the particle is switched on, driving the system to a new stationary state. Our study calculates displacement moments in discrete time (number of steps N) for an arbitrarily strong constant driving force, exact to first order in obstacle density. 
We find that for fixed driving force F, the approach to the terminal discrete velocity scales as ∼ N^-1exp(- N F^2 / 16) for small F, differing significantly from the ∼ N^-1 prediction of linear response. 
Besides a non-analytic dependence on the force and breakdown of Einstein's linear response, our results show that
fluctuations in the directions of the force are enhanced in the presence of obstacles. Notably, the variance grows as 
∼ N^3 (superdiffusion) for F →∞ at intermediate steps, reverting to normal diffusion (∼ N) at larger steps, a behavior previously observed in continuous time but demonstrated here in discrete steps for the first time. Unlike the exponential waiting time case, the superdiffusion regime starts immediately at N=1.
The framework presented allows considering any type of waiting-time distribution between steps and transition to continuous time using subordination methods. 
Our findings are also validated through computer simulations.


Driven Lorentz model in discrete time
    Thomas Franosch
    September 9, 2024
=====================================










§ INTRODUCTION


The transport of molecules, colloids and particulate matter in disordered media is ubiquitous in natural, industrial and technological processes. Such systems have been extensively studied using a random-walk-on-a-lattice approach. One of the best-known examples is the continuous-time random walk <cit.> (CTRW), popularized by Montrol and Scher  <cit.> to model charge transport in amorphous materials and later used very successfully for the description of transport in porous and biological mediums <cit.>. The idea behind CTRW is to build upon the classical random walk, which is a succession of random steps, by introducing a waiting time distribution between steps, i.e., the medium is a form of energetic landscape which gives rise to waiting or trapping times. It has been shown that when these waiting times follow a scale-free distribution where the mean waiting time diverges, it leads to intriguing phenomena such as anomalous diffusion, weak ergodicity breaking and aging, seen in single quantum dots <cit.>, transport of biomolecules inside the cell <cit.> and glassy systems <cit.> just to name a few.

A second prominent model is the lattice Lorentz gas <cit.>, describing obstructed transport in heterogeneous environments such as the crowded world inside biological cells. The model consists of a tracer particle performing random walk on a lattice where a fraction of the sites (density) is occupied by immobile hard obstacles. These obstacles, placed at random lattice locations, are treated with reflecting boundary conditions. 
Many works probe the characteristics of the Lorentz model system by studying the response of the particle to an external driving force <cit.>.
The emphasis is usually for a near-neighbor hopping process where the average waiting time between jumps is finite (often an exponential distribution). Already for finite average waiting time between hops, the presence of obstacles alters the dynamics of the driven system in a non-trivial manner since repeated collisions of the tracer with obstacles introduce correlations and persistent memory <cit.>. 
For example, in Ref. <cit.>, a first-order expansion in the density of immobile obstacles already presents a surprising force-dependent exponential decay towards the steady-state drift velocity; In contrast, a very dense system study reveals a surprising very short lived initial high velocity value that abruptly drops to a terminal 'low' value <cit.>.

The combination of obstacles with different types of distributions of waiting times between steps – which we refer to as temporal disorder, has received limited attention so far. One instance of this is a power law distribution resulting in diverging mean waiting times.
In our work, we find the moments of displacement in the domain of number of steps, i.e. discrete time, for arbitrarily strong constant driving force, accurate to first order in the obstacle density. This theoretical approach will establish a framework for future research that will allow to consider any type of temporal disorder and transition to continuous time by means of the method of subordination <cit.>; i.e., a summation of conditional probability on all the possible outcomes of the number of steps during total time t of the process. 
Our solution technique employs a scattering formalism borrowed from quantum mechanics <cit.> that has been successfully applied to the analytic study of the driven lattice Lorentz gas in continuous time with an exponential distribution between steps <cit.>. 

In this work, the discrete nature of working in number of steps leads to different mathematical challenges and results compared to the continuous case, since we will be dealing mainly with summation techniques (generating functions) instead of integrals (Laplace transforms). We show that for fixed driving F, the step-dependent approach (N discrete time) towards the terminal discrete velocity (average position divided by steps) is exponentially fast rather than a power law decay as predicted by linear response.
Furthermore, we provide the first and second moment to first order in the obstacle density. Showing the dependence on the force is complex and non-analytic, indicating breaking of Einstein's linear response even for small forces.
The intuitive picture is that obstacles suppress the fluctuations in the direction of the force. Our results indicate that for forces large enough, increasing disorder leads to an enhancement. We show that there is a window at intermediate values of steps N where the variance grows as ∼ N^α with a true exponent of α=3 in the limit F →∞ which then drops to normal diffusion (α=1) at large steps. Such intriguing behavior has been already documented for the lattice Lorentz gas <cit.>. We provide the first evidence that such an anomalous behavior occurs also in the discrete domain of the number of steps. Our findings are validated by high-precision stochastic simulations.

As a consequence of the approach developed in this paper, our framework will find natural application to study the impact of heavy-tailed distribution of waiting times with a diverging average in the presence of obstacles.
Another example that deserves to be considered in the presence of obstacles is the case of a temporal quenched disorder which leads to correlations and memory effects <cit.>. A quenched temporal disorder by itself is known to exhibit surprising effects, including mobility enhancement in confined geometries <cit.> and non self-averaging leading to universal fluctuations of diffusivity <cit.>.
Such interpretations which combine temporal and obstacle disorder may be a more complete picture to describe transport in a wide variety of systems and can lead to new theoretical advancements in the field.



§ THE MODEL AND SOLUTION TECHNIQUE
 
In the two-dimensional lattice Lorentz gas, a tracer particle performs a random walk on a square lattice of size L × L (where L ∈ℕ) defined by the collection of sites 𝐫∈Λ={(x, y) ∈ℕ×ℕ : (1 ≤ x , y ≤ L) }.
Here the lattice spacing a is set to unity for convenience such that the lattice sites assume only integer values and the total number of steps performed is N.
We assume that the number of sites L^2 is very large and approaching the limit L →∞, i.e., the thermodynamic limit.
At the boundaries we employ periodic boundary condition.
In every step the tracer performs a nearest-neighbor jump of size 𝐝∈𝒩={±𝐞_x, ±𝐞_y} where 𝐞_x and 𝐞_y are perpendicular unit vectors in the x and y direction respectively. 
The lattice consists of free sites accessible to the tracer as well as sites with randomly placed immobile
hard obstacles of density n (fraction of excluded sites). If the tracer attempts to jump onto an obstacles site, it remains at its initial position before the jump but still the counter for the discrete time is increased by one.
At zero steps (discrete time), a constant force acting on the tracer is switched on, and we use the thermal equilibrium state in the absence of driving F = 0 as the initial condition, i.e., the particle is equally likely to be anywhere on the L^2 sites at 'time' N=0.
The external force pulls the tracer only along the x direction of the lattice and the strength of the force is characterized by the dimensionless force F= force× a / k_B T, where k_B is Boltzmann's constant and T is the temperature. 
The transition probabilities for a single step W(𝐝) obey detailed balance, W(𝐞_x) / W(-𝐞_x)=e^F in the x direction and W(𝐞_y) / W(-𝐞_y)=1 in the y direction. 
Correspondingly, we choose the transition probabilities parallel and perpendicular to the applied force as

    W(±𝐞_x)=Γ e^± F / 2,

and

    W(±𝐞_y)=Γ,

respectively as depicted in Fig. <ref> where

    Γ = 1 / (e^F / 2+e^-F / 2+2),

is the normalization factor.




We use an approach similar to Ref. <cit.> to find the propagator (defined below) which will enable us to derive the moments of the displacement as a function of the number of steps N.
We first begin for the case of no obstacles, i.e., the free system, since it is more simple and to establish our formalism.
It is convenient to exploit the analogy of the master equation to a Schrödinger equation. 
Hence we consider the Hilbert space of lattice functions Λ→ℂ spanned by the orthonormal basis of position kets |𝐫⟩. 
We denote by  p_N(𝐫)  the probability     for the tracer to be at 𝐫 after N jumps and define an abstract ket |p_N ⟩ := ∑_𝐫∈Λ p_N(𝐫) |𝐫⟩. The probabilities can thus be obtained as p_N(𝐫) = ⟨𝐫 | p_N ⟩.
We assume that the force is switched on at step N=0 such that the thermal equilibrium state |p_0 ⟩ evolves towards a 
new stationary state | p_st⟩ = lim_N→∞ | p_N ⟩ indicated by the subscript "st".
We define the state of the system at step N=0 to be the equilibrium state of the empty lattice with no driving, ⟨𝐫 | p_0 ⟩ = 1/L^2, meaning the probability to find the tracer at any of the L^2 sites of the lattice is uniform.
We define the single jump matrix for the free system with no obstacles indicated by the subscript 0 as


    M̂_0=∑_𝐫∑_𝐝∈𝒩 W(𝐝)|𝐫⟩⟨𝐫-𝐝|,


which entails the rule to propagate probabilities in time by | p_N+1⟩ = M̂_0 | p_N⟩.
Here the matrix element ⟨𝐫 | M̂_0 | 𝐫'⟩ is the transition probability from 𝐫' to 𝐫.
Now we define the free propagator to be the z-transform (sometimes called generating function) <cit.> of the single jump matrix M̂_0 


    Ĝ_0(z) = ∑_N=0^∞ (M̂_0)^N z^N.


It is convenient to perform a spatial Fourier transform, defined by


    |𝐤⟩=1/L∑_𝐫∈Λexp (i𝐤·𝐫)|𝐫⟩ ,


where 𝐤=(k_x, k_y) ∈Λ^*={(2 π n_x / L, 2 π n_y / L): (n_x, n_y) ∈Λ} and 𝐤·𝐫=k_x x+k_y y. 
Since the matrix M̂_0 is translationally invariant, in the plane wave basis becomes diagonal,  i.e., ⟨𝐤 | M̂_0 | 𝐤' ⟩ = λ (𝐤) δ_𝐤, 𝐤' with λ(𝐤) the eigenvalues, sometimes referred to as the characteristic function <cit.>.
The jumps are independent of each other thus the free propagator Ĝ_0 is also diagonal in k-space,

    G_0(z, 𝐤)    =   ⟨𝐤 | Ĝ_0(z) | 𝐤⟩
       =   ∑_N [λ(𝐤)]^N z^N = [1-z λ(𝐤)]^-1.

The function G_0(z, 𝐤) is also called the moment generating function since (-i)^m∂^m G_0(z, 𝐤) / ∂ k_x^m  |_𝐤=0 is  the m-th moment of displacement along the x axis (parallel to the direction of the applied force) <cit.>.
This will provide the solution for the first and second moments in the obstacle-free system, a similar approach will be taken if the general propagator for the case of obstacles G(z,𝐤) is known.
Thus the solution strategy of finding the moments will be achieved by finding the general propagator G(z,𝐤).
This would yield the z -transform of the moments and later we show how we switch back to N-space.
We now turn to finding  G(z,𝐤) in the case of obstacles. 



The general case with obstacles can be obtained by relying on the scattering formalism borrowed from quantum mechanics <cit.>.
The dynamics in the presence of randomly distributed obstacles on the lattice is generated by the modified single jump propagator M̂ = M̂_0 + V̂, where V̂ cancels the transitions from and to the obstacles. 
Furthermore, in our calculations we allow the tracer to start at an obstacle site, which then remains immobile.
The potential v̂_1(𝐬_1) for a single obstacle at site 𝐬_1 cancels the transition probabilities from and into the obstacle position.
Therefore the only non vanishing elements of the matrix v̂_1(𝐬_1) involve the obstacle site and its immediate four neighbors.
In Sec. <ref> we will explicitly show the resulting matrix v̂_1(𝐬_1). 
Now for J obstacles (or impurities) the total potential will be V̂ = ∑_i=1^Jv̂_i.
We define the obstacle density as

    n = J / L^2.

Strictly speaking this definition does not properly account for obstacles that could occur as neighbors since then they will affect each others potential. 
But in the limit of large lattices L→∞ and small densities n such realizations rarely occur and therefore we neglect them in our work. 
The propagator in the presence of a fixed obstacle realization V̂ = ∑_i=1^Jv̂_i  is related to the free propagator Ĝ_0 via a Lippmann-Schwinger equation <cit.>


    Ĝ = Ĝ_0 + Ĝ_0 V̂Ĝ.


By iterating we can express Ĝ as


    Ĝ   = Ĝ_0 + Ĝ_0 V̂( Ĝ_0 + Ĝ_0 V̂Ĝ)
       = Ĝ_0 + Ĝ_0 V̂Ĝ_0 + Ĝ_0 V̂Ĝ_0 V̂Ĝ ,


and this can be further expended.
By following this line of reasoning, the propagator can be expressed as follows

    Ĝ = Ĝ_0 + Ĝ_0 T̂Ĝ_0 ,

with the scattering matrix T̂ defined by


    T̂   =   V̂+V̂Ĝ_0 V̂+V̂Ĝ_0 V̂Ĝ_0 V̂+⋯
       =   ∑_i=1^Jv̂_i+∑_j, k=1^Jv̂_j Ĝ_0 v̂_k+∑_l, m, n=1^Jv̂_l Ĝ_0 v̂_m Ĝ_0 v̂_n+⋯.


We use the single-obstacle scattering operator t̂_i which encodes all possible collisions with obstacle v̂_i defined in z-space by


    t̂_i=v̂_i+v̂_i Ĝ_0 v̂_i+v̂_i Ĝ_0 v̂_i Ĝ_0 v̂_i+⋯,


and express T̂ through a multiple scattering expansion, which reads


    T̂=∑_i=1^Jt̂_i+∑_j, k=1 
     j ≠ k^Jt̂_j Ĝ_0 t̂_k+∑_l, m, n=1 
     l ≠ m, m ≠ n^Jt̂_l Ĝ_0 t̂_m Ĝ_0 t̂_n+⋯.


We are not interested in the dynamics for individual realizations of the obstacle configurations but only in disorder-averaged properties. We use [·]_av to indicate an average over all realizations of the disorder. In particular, after disorder-averaging the system is translationally invariant. The convenient starting point for disorder averaging is Eq. (<ref>). In particular, 
[T̂]_av is translationally invariant in space and thus diagonal in the plane-wave basis.
Contributions in first order of the density can now be identified with the forward-scattering amplitude t_i(z, 𝐤)=⟨𝐤|t̂_i| 𝐤⟩ of a single obstacle placed in a random location on the lattice, therefore


    ⟨𝐤|[T̂]_av| 𝐤⟩ = n L^2 t_i(z, 𝐤)+O(n^2).


Eq. (<ref>) states that the disorder averaged scattering operator is proportional, in 𝐤-space, to the forward scattering amplitude of a single obstacle.
The Lippmann-Schwinger equation <cit.> allows us to express t̂_i as


    t̂_i=v̂_i+v̂_i Ĝ_0 t̂_i=v̂_i+t̂_i Ĝ_0 v̂_i,


therefore


    ⟨𝐫|t̂_i| 𝐫^'⟩   =   ⟨𝐫|v̂_i(1-Ĝ_0 v̂_i)^-1| 𝐫^'⟩
       =   ⟨𝐫|(1-v̂_i Ĝ_0)^-1v̂_i| 𝐫^'⟩.


Since v̂_i only has nonvanishing contributions for the obstacle site and its nearest neighbors, the calculation of  ⟨𝐫|t̂_i| 𝐫^'⟩ in Eq. (<ref>) reduces to a 5 × 5  matrix inversion problem.

The real-space matrix elements ⟨𝐫|Ĝ_0| 𝐫^'⟩ of the free propagator for sites around the obstacle are expressed in terms of complete elliptic integrals of the first and second kind (we show this in Sec.  <ref>).
The forward-scattering amplitude is then obtained by a change of basis with,


    L^2 t_i(z, 𝐤)=∑_𝐫, 𝐫^' e^i𝐤·(𝐫-𝐫^')⟨𝐫|t̂_i| 𝐫^'⟩ ,


where the sum effectively extends only over the obstacle site and its nearest neighbors.
The disorder-averaged propagator to first order in the density of obstacles n (defined in Eq. (<ref>)) is then according to Eqs. (<ref>) and (<ref>) 


    [G]_av(z, 𝐤)=G_0(z, 𝐤)+n L^2 G_0(z, 𝐤)^2 t_i(z, 𝐤)+O(n^2).


We still need to account for one more correction.
Since the starting position is random, the tracer can start at an obstacle site where it will stay immobile forever. This obstacle can be surrounded by additional obstacle sites but to first order in the density this can be ignored by assuming obstacle sites are far enough from each other.
Starting movement at an obstacle site is not physical, therefore we correct for this behavior by simply multiplying the propagator [G]_av(z, 𝐤) with 1 / (1-n) = 1 + n +O(n^2), where (1 - n) is the fraction of free lattice sites. We denote the corrected propagator by [G]_av^c(z, 𝐤) and keep only terms to first order in n


    [G]_av^c(z, 𝐤)    =    G_0(z, 𝐤)+n [ G_0(z, 𝐤) 
           + L^2 G_0(z, 𝐤)^2 t_i(z, 𝐤) ] +O(n^2).





§ THE PROPAGATOR
 

As shown in Sec. <ref>, finding the propagator in the case of obstacles in 𝐤-space will enable us to calculate the moments of displacement.
According to Eq. (<ref>), the missing piece is the single obstacle forward-scattering amplitude t_i(z, 𝐤)=⟨𝐤|t̂_i| 𝐤⟩.
Since t̂_i is provided in Eq. (<ref>) in real-space, we have to find the ingredients ⟨𝐫 | Ĝ_0 | 𝐫' ⟩ and ⟨𝐫 | v̂_i | 𝐫' ⟩ and then transition to 𝐤-space.




 §.§ The single obstacle potential matrix v̂_i
 

For a reflective obstacle placed at 𝐬_i the new transition probabilities (M̂ + v̂_i) should cancel the transition from and to the obstacle site. 
Hence v̂_i is a 5 × 5 matrix with elements ⟨𝐫 | v̂_i | 𝐫' ⟩ that correspond to the transition from 𝐫' to 𝐫 with the basis  𝐫,𝐫^'∈{𝐬_i ±𝐞_y, 𝐬_i ±𝐞_x, 𝐬_i } and we set 𝐬_i for convenience to be zero.
The column index is 𝐫' and the row index is 𝐫.
The reflective potential cancels the transition from and to the obstacle site, hence ⟨𝐬_i|v̂_i| 𝐬_i-𝐝⟩=-W(𝐝) and ⟨𝐬_i-𝐝|v̂_i| 𝐬_i⟩=-W(-𝐝).
If the tracer is at a neighboring site and attempts to jump to 𝐬_i it is scattered back to its original position, therefore ⟨𝐬_i-𝐝|v̂_i| 𝐬_i-𝐝⟩=W(𝐝).
Finally if the tracer is at the obstacle site it is stuck there forever, ⟨𝐬_i|v̂_i| 𝐬_i⟩=1.
We order the rows and columns for the matrix forms of the operators in the real-space basis via the scheme 

    [   1  ; 2 3 4;   5   ]

where the obstacle site is located at the origin 0 numbered by 3.
Consequently the order is 𝐞_y, -𝐞_x, 0, 𝐞_x, -𝐞_y.
The resulting potential matrix in our basis is now


    v̂_i =Γ([         1         0        -1         0         0;         0   e^F / 2 -e^-F / 2         0         0;        -1  -e^F / 2     1 / Γ -e^-F / 2        -1;         0         0  -e^F / 2  e^-F / 2         0;         0         0        -1         0         1 ]).




 §.§ The matrix elements of Ĝ_0(z)
 
The obstacle-free propagator is provided in Eq. (<ref>) in 𝐤-space, to find the matrix elements in 𝐫-space we perform an inversion to find

    ⟨𝐫 | Ĝ_0(z) | 𝐫'⟩ =  ∫_-π^π∫_-π^πd 𝐤/(2 π)^2exp(-i 𝐤·(𝐫-𝐫^')) /1-zλ(𝐤).


The entries of the matrix Ĝ_0 are the same as defined in Sec. <ref>, with
𝐫,𝐫^' =  𝐞_y, -𝐞_x, 0, 𝐞_x,-𝐞_y.
The eigenvalues λ(𝐤) of the single jump matrix M̂_0 follow directly from translational invariance in the 𝐤-basis


    ⟨𝐤 | M̂_0 | 𝐤⟩ =
        ∑_𝐝∈𝒩exp[i 𝐤·𝐝] W(𝐝)


with the result


    λ(𝐤)    =    2 Γ[ cos (k_x) cosh (F/2)   
           + i sin (k_x) sinh (F/2) + cos (k_y) ],


Γ was defined in Eq. (<ref>) to be the normalization of the transition probabilities.
We show in the Appendix <ref> that the matrix elements of the free propagator can also be expressed in the following form


    ⟨𝐫 | Ĝ_0(z) | 𝐫'⟩   =    e^F (x - x') / 2∫_0^∞ d s e^-s
          ×   I_|x-x'|(2 Γ s z ) I_|y-y'|(2 Γ s z) ,


where I_m(…) is a modified Bessel function of the first kind of integer
order m. 
Here, the exponential outside of the integral accounts for the asymmetry introduced by the bias and the remaining terms are the solution of the symmetrical problem but with the diffusion coefficient set to be Γ(F) instead of Γ(F=0). 
As such we denote the obstacle-free unbiased part by


    g_x-x',y-y' = ∫_0^∞ d s  e^-s I_|x-x'|(2 Γ s z ) I_|y-y'|(2 Γ s z),


and we obtain ⟨𝐫 | Ĝ_0(z) | 𝐫'⟩ = e^F (x - x') / 2 g_x-x',y-y'. By using the symmetry of the unbiased part g_x y, with respect to x, y and also the fact that g_x y=g_y x we obtain for the matrix Ĝ_0 (z) the following representation 


    Ĝ_0(z)=([          g_00  e^F / 2 g_11          g_10 e^-F / 2 g_11          g_20; e^-F / 2 g_11          g_00 e^-F / 2 g_10     e^-F g_20 e^-F / 2 g_11;          g_10  e^F / 2 g_10          g_00 e^-F / 2 g_10          g_10;  e^F / 2 g_11      e^F g_20  e^F / 2 g_10          g_00  e^F / 2 g_11;          g_20  e^F / 2 g_11          g_10 e^-F / 2 g_11          g_00 ]).


We note that by virtue of the underlying dihedral symmetry <cit.> that only four elementary propagators are needed to be explicitly calculated, they are: g_00,  g_10, g_20 and g_11.
As illustrated in Appendix <ref> these four propagators satisfy a set of linear relationship that eventyally reduce the number of independent propagators from four to two.
The resulting relations are

    g_00   = 1 + 4 Γ z g_10, 
    
            g_10   = Γ z ( g_00 + g_20 + 2 g_11 ).

We turn now to evaluating g_00 and g_11.
From Eq. (<ref>),

    g_00 = ∫_0^∞ d s e^-s I_0(2 Γ s z)^2 = 2/πK(16 Γ ^2 z^2)


    g_11   = ∫_0^∞ d s e^-s I_1(2 Γ s z)^2
       = 2/π (4 Γ z)^2{[2-(4 Γ z)^2] K[(4 Γ z)^2]
       -2 E[(4 Γ z)^2]}

where


    K(k) = ∫_0^π / 2d α/√(1-k sin ^2 α)


and

    E(k) = ∫_0^π /2√(1-k sin ^2 α)  d α

are the complete elliptic integral of the first and second kind respectively.
From Eq. (<ref>), g_10 and g_20 are easily found.



 §.§ Scattering matrix t̂_i
 
To sum up the work done so far; we have found all the ingredients required to calculate the single-obstacle scattering matrix t̂_i in z-space and in the plane-wave basis. By using 
Eq.(<ref>) together with Eq.(<ref>), we find

    L^2 t_i(z, 𝐤)   =L^2⟨𝐤|t̂_i| 𝐤⟩=
       ∑_𝐫, 𝐫^'exp[i 𝐤·(𝐫-𝐫^')]⟨𝐫|v̂_i(1-Ĝ_0 v̂_i)^-1| 𝐫^'⟩.

where v̂_i is given in Eq.(<ref>) and Ĝ_0 in Eq. (<ref>) is a function of the four obstacle-free propagators, g_00,    g_10,    g_20,    g_11. The problem turns into a 5 × 5 matrix inversion problem which we solve using computer algebra. 
Finding t_i(z, 𝐤) allows us to obtain the propagator [G]_av^c(z, 𝐤) for the system with obstacles using Eq. (<ref>). As mentioned before, this propagator is also the moment generating function. Thus, we derive the moments of displacement by taking the appropriate derivative of [G]_av^c(z, 𝐤).
The m-th moment of displacement along the x axis, i.e., in parallel to the direction of the applied force (see Fig. <ref>), is now

    ⟨x̃^m(z) ⟩    := ∑_N=0^∞⟨ x(N)^m ⟩ z^N 
       =  (-i)^m∂^m [G]_av^c(z, 𝐤) / ∂ k_x^m |_𝐤=0.

Note that in the z domain, we need to distinguish x̃^m(z) from x̃(z)^m.
We used the notation ⟨⋯⟩ in order to indicate the average over the randomness (many obstacle realizations and different trajectories).
In what follows, we investigate the behavior of these moments for any number of steps N, and in particular the convergence towards the terminal velocity v_∞ := lim_N→∞⟨ x(N) ⟩ / N and the variance [ ⟨ x(N)^2 ⟩ - ⟨ x(N) ⟩ ^2 ] as a function of N. The moments are derived for arbitrary force F correct to first order in the density of obstacles n and compared to computer simulations to find the range of validity.



§ FIRST MOMENT

Equation  (<ref>) allows computing the first moment of the displacement for any z and force F,

    ⟨x̃ (z) ⟩ =    -i ∂/∂ k_x[ (1 + n) G_0 (z, 𝐤) 
       + n L^2 G_0 (z, 𝐤)^2 t_i(z, 𝐤)] |_𝐤=0.

We remind the reader that G_0 (z, 𝐤) is the obstacle-free propagator and is provided by Eq.  (<ref>) giving G_0 (z, 𝐤=0) = 1/(1-z). 
By expanding Eq. (<ref>) we obtain 

    ⟨x̃ (z) ⟩ =   . (1+n) (-i ∂/∂ k_x G_0(z,𝐤) ) |_𝐤=0
         + n [ -2iL^2 G_0(z,𝐤) t_i(z,𝐤) ( ∂ G_0/∂ k_x) .
         +. . G_0(z,𝐤)^2 ( -iL^2 ∂ t_i(𝐤)/∂ k_x) ] |_𝐤=0.

Additionally, by symmetry of the problem and verified using computer algebra, we find t_i(z,𝐤=0)=0.
We now have an expression for ⟨x̃ (z) ⟩,


    ⟨x̃ (z) ⟩   -    (1 + n) ⟨x̃_0(z) ⟩ = 
       =    n  1/(1-z)^2.(-i L^2 ∂ t_i(z, 𝐤)/∂ k_x) |_𝐤=0,


we find it using computer algebra. The resulting expression is extremely long and provided in the appendix [Eq. (<ref>)].
Here ⟨x̃_0(z) ⟩ is the solution of the bare dynamics which we know after transitioning to N-space must be ⟨ x_0(N) ⟩ = tanh (F/4) N. This is found by averaging over the displacement of a single step (see Sec. <ref>) and multiplying by N.
As the solution in Eq. (<ref>) is in z-space, we look for the solution for any number of steps N by inverting the process.
Inverting the z-transform in Eq. (<ref>), we obtain the average displacement

    ⟨ x(N) ⟩ = 1/N!. d^N/dz^N⟨x̃ (z) ⟩|_z=0.

This provides the result for any N and F.
Successively, we will derive in Sec. <ref> the terminal velocity and we will use this formula for ⟨ x(N) ⟩ in order to test the convergence rate towards the terminal velocity in theory and simulations.



§ TERMINAL VELOCITY
 
To determine the asymptotic behavior of the discrete velocity, v_∞ := lim_N→∞( ⟨ x(N) ⟩ / N ), we make use of the Tauberian theorem <cit.> (see the appendix for the exact method).
The Tauberian Theorem allows transitioning to N-space and to find ⟨ x(N) ⟩ for large N from the behavior of its generating function ⟨x̃ (z) ⟩ =∑_N=0^∞⟨ x(N) ⟩ z^N when z→ 1.
We use in Eq. (<ref>) the asymptotic relation, .(-i L^2  ∂ t_i(z, 𝐤)/∂ k_x )|_𝐤=0 = Δ v_∞+O(1-z) obtained in Appendix  <ref> where Δ v_∞ is a function of  F only.
Consequently the 1/(1-z)^2 term in Eq. (<ref>) transforms to a linear growth in N and we obtain,


    ⟨ x(N) ⟩∼ (1+n) v_0 N + n Δ v_∞ N,


as N→∞. Here  v_0=tanh (F/4) denotes the velocity of the obstacle-free lattice. Therefore,

    v_∞   =   lim_N →∞ (⟨ x(N) ⟩) / N 
       =    v_0 + n (v_0 +Δ v_∞).


The terminal behavior is now found for arbitrarily strong driving F and small densities n, while the complete expression is very long [see Eq. (<ref>)], we can elaborate the behavior for small forces. Relying on computer algebra, we find

    v_∞ = D_x F + n/16(π/4 -1) F^3 log (F) + O(F^3),

which immediately highlights the nonanalytic behavior in F.
Here D_x = [ 1 - n (π - 1) ] / 4 corresponds to the diffusion coefficient in the x-direction in the absence of force as expected from linear response.
To make a connection with the continuous-time case, N needs to be treated as a random variable, where the duration τ of a step is taken from a distribution ψ(τ). This is essentially the renewal theorem; the m-th moment in the time domain would now be given by  ⟨ x(t)^m ⟩ = ∑_N ⟨ x(N)^m ⟩ Q_N(t),
where Q_N(t) is the probability for having performed exactly N steps until time t. For large times, this can be simplified, and the number of steps is switched with the total time passed t divided by the expected time τ of a jump (as long as it is finite). When the times between steps is exponentially distributed with a mean of one, we immediately recover the expression already found in a previous work <cit.> by just switching N with t in Eq. (<ref>), therefore validating our result.

Our result in Eq. (<ref>) holds as long as obstacles are positioned far enough from each other such that sequences of collisions involving obstacles that have been encountered previously can be ignored.
Then, in our derivations, the assumption that [T]_av(𝐤), the disorder-averaged scattering operator in 𝐤-space, is the sum of nL^2 identical single-obstacle forward-scattering amplitudes is correct [Eq. (<ref>)].
To investigate the range of validity of the low-density expansion, we plot the terminal response as a function of F in Fig. <ref>. For each curve, there is a peak indicating a transition for stronger driving where the particle is more often stuck on an obstacle instead of going around it. Thus, the lattice Lorentz model displays negative differential mobility, a phenomenon that was seen in Ref. <cit.> as well.
While the theoretical lines in Fig. <ref> agree nicely with simulations, the suppression of the stationary velocity is underestimated at strong forces indicating that contributions of higher order in the density become relevant. In other words, the range of validity of our approach for the small obstacle densities is force-dependent.







§ APPROCH  TOWARDS THE TERMINAL VELOCITY
 

Next we investigate in more detail the instantaneous velocity defined by 

    v(N) := ⟨ x(N) - x(N-1) ⟩    for  N=1,…

where we set v(0)  = 0.  In the z-domain, it is directly related to the displacement via

    ṽ(z) = (1-z) ⟨x̃(z) ⟩

The behavior of the instantaneous velocity at large step sizes is reflected in the poles of the z-transform. We have already identified a pole of order 2 in ⟨x̃(z) ⟩ which translates to a simple pole in ṽ(z) reflecting that the velocity approaches a constant v(N) → v_∞ as N→∞. Further singular behavior emerges from the single-obstacle t-matrix t̂_i (z ,𝐤) [Eq. (<ref>)]. Since the inversion of the 5× 5 matrix does not give rise to new singular behavior, all non-analytic properties are inherited by the matrix elements, i.e., by the 

    ⟨𝐫 | Ĝ_0(z) | 𝐫'⟩ =  ∫_-π^π∫_-π^πd 𝐤/(2 π)^2exp(-i 𝐤·(𝐫-𝐫^')) /1-zλ(𝐤).

For the case of F=0, there is a singularity at z → 1 is originating from λ(𝐤 = 0) = 1. 

In Sec. <ref> we have established a linear dependence of these integrals on the complete elliptic integrals of the first and second kinds, K and E, respectively.
The leading term arises from the expansion of K around z=1,

    . ⟨0 | Ĝ_0(z) | 0⟩|_F=0 = 2/πK(z^2) ∼ -1/πlog(1-z),

while E is subdominant (see Ref. <cit.>).  Upon Taylor expansion of the logarithm

    - log (1-z)=    ∑_N=1^∞z^N/N,

we read off that logarithm corresponds to an algebraic decay ∝ 1/N for N→∞. 
This corresponds to the decay rate expected from linear response when F → 0 as seen in Fig. <ref>. 
For the case of F>0, the integral in Eq. (<ref>) is no longer divergent at z=1 but the singularity is shifted. 
 Using the same argument as in Eq. (<ref>), the leading term is found to 

    ⟨0 | Ĝ_0(z) | 0⟩ = 2/πK(16 Γ ^2 z^2) ∼ -1/πlog( 1-4 Γ z ),

i.e. we anticipate a singularity at z=1/4 Γ > 1. By the scaling property, the z-transform of (4Γ)^-N [v(N)-v_∞] is ṽ(z/4Γ)- v_∞/(1- z / 4Γ )  ∝ -log(1-z) for z→ 1 , and we infer that the 
large-step behavior of the velocity in this case is 

    v(N) - v_∞∝exp[ N log (4 Γ) ] /N.

Therefore the power law tail ∼1/N no longer determines the asymptotic large-N behavior for any finite bias F; rather the decay rate is exponentially fast  which can be elaborated further for small F≪ 1 to be ∼exp (-N F^2 / 16) / N. A comparison of this analytical prediction with simulation is shown in  Fig. <ref>. 
This result is consistent with the continuous case with exponential waiting times <cit.> where it was shown that the regular algebraic decay ∼ t^-1 is elaborated to ∼ t^-1exp (-F^2 t / 16) for F → 0. 
The onset time of the exponential behavior can be found by comparing the decay rate with that of linear response, i.e., when exp[N log( 4 Γ)] / N ≪ 1/N. We find N ≫ |1 /  log(4 Γ)| = N_f, where N_f is the onset step-time. For small F, N_f ∼ 16/F^2, thus we establish that as the forces are smaller the onset time is delayed significantly as can be seen in Fig. <ref>.

We remark that if we instead examine the average velocity, ⟨ x(N) ⟩ / N, the exponential decay will become subdominant and we will observe the asymptotic behavior ∼ 1/N for N ≫ 1. 








§ SECOND MOMENT

By employing the same approach as for the first moment, the second moment of displacement in the x direction, sometimes called mean-squared displacement (MSD) is found by taking the second derivative of the propagator in Eq. (<ref>). From Eq. (<ref>) we obtain,

    ⟨x̃^2(z) ⟩ =    (-i)^2 ∂^2/∂ k_x^2[ (1+n) G_0 (z, 𝐤) 
       + n L^2 G_0 (z, 𝐤)^2 t_i(z, 𝐤)] |_𝐤=0.

The first term here is the obstacle-free system which we denote by ⟨x̃_0^2(z) ⟩. 
Using G_0 (z, 𝐤=0) = (1-z)^-1 and that t_i(z, 𝐤=0)=0, after some simplification we are left with

    ⟨x̃^2(z) ⟩   =    (1 + n) ⟨x̃_0^2(z) ⟩ - n [ 4 .( i L^2 ∂ t_i (z, 𝐤)/∂ k_x)|_𝐤=0. 
          ×tanh(F/4) z/(1-z)^3
           + .. 1/(1-z)^2( L^2 ∂^2/∂ k_x^2 t_i (z, 𝐤) )|_𝐤=0].

Making use of the inverse z-transform as in Eq. (<ref>), we find the exact solution for ⟨ x(N)^2 ⟩ for any N and F.
This result is used later to find the variance and how it behaves for small number of steps N in Sec. <ref>. We continue to find the asymptotic behavior of the second moment in the limit of N →∞.



§ SECOND MOMENT  IN THE REGIME OF LARGE STEP NUMBERS
 
We use the same methodology as in Sec. <ref> to determine the asymptotic behavior of the MSD, ⟨ x^2 (N) ⟩, in the regime of N≫ 1.
We therefore use in Eq. (<ref>) the (1-z) expansions of the derivatives of the scattering matrix t_i (z,𝐤) up to order O(1-z). 
Relying on the results derived in Appendix <ref>, we obtain .(-i L^2 ∂ t_i(z, 𝐤) / ∂ k_x )|_𝐤=0 = Δ v_∞+O(1-z) and 
. ( L^2 ∂^2 t_i(z,𝐤) / ∂^2 k_x ) |_𝐤=0=-2 (v_0)^2 z/(1-z) +c_1(F) +O(1-z). Here, c_1(F) here is very long expression that depends only on F, its explicit form can be inferred from Eq. (<ref>).
Consequently, using again the Tauberian theorem, the 1/(1-z)^2 term in Eq. (<ref>) transform to N and the z/(1-z)^3 term transforms to N^2 / 2.
Collecting results, we find asymptotically

    ⟨ x(N)^2 ⟩∼    (1+n) ⟨ x_0(N)^2 ⟩
       +n [ 2 v_0 Δ v_∞ N^2 + v_0^2 N^2 + c_1(F) N ],

for N→∞. Here the ⟨ x_0 (N)^2 ⟩ term is determined by the bare diffusion dynamics, ⟨ x_0(N)^2 ⟩ = ( 1-tanh ^2 (F/4) ) N/2 + tanh ^2 (F/4) N^2.
This result is used later to determine the behavior of the variance in the regime N≫ 1 in Sec. <ref>. For small forces, the expression in Eq. (<ref>) can be simplified further, using Δ v_∞ = -π F / 4 +(π / 4 - 1) F^3 log (F) / 16 + O(F^3) and approximating c_1=-π/2 - (0.333005 +0.106403 log (F)) F^2 + O(F^4), c_1 is given in full detail in the Appendix [Eq. (<ref>)], we obtain

    ⟨ x(N)^2 ⟩ - (1+n) ⟨ x_0(N)^2 ⟩≈
          -n {( 2 π - 1 )/16 F^2 N^2 + [ π/2 + (1/3 +1/10log (F) ) F^2 ] N}.

The force-free term [ 1-n(π -1) ] N / 2 in the second moment of the displacement indicates that when no force is applied (F=0), the obstacles still obstruct the movement of the tracer.
This can be seen in the plot [see Fig. <ref>] of the second moment of displacement in the x direction. Higher-order corrections show that there is a complex non-linear dependence on the force, specifically the logarithmic dependence on F that becomes even more relevant in the variance since the N^2 dependent terms drop, as we show in the next section.






§ THE VARIANCE
 
In order to determine the variance of the displacement [x(N)]:=⟨ x(N)^2 ⟩ - ⟨ x(N) ⟩^2, we remark that we take contributions only to first order in O(n) when taking the square   of ⟨x̃ (z) ⟩ in Eq. (<ref>). Then, the exact solution is determined by taking the inverse z-transform of ⟨x̃^2 (z) ⟩ - ⟨x̃ (z) ⟩^2 as was done in Eq. (<ref>) for the first moment.
In Fig. <ref> and Fig. <ref> we plot the variance in the direction of the force against numerical simulations for small number of steps N. 
For this purpose we use the variance of the bare dynamics [x_0(N)] = ( 1- tanh^2 (F/4) ) N / 2.
We see in Fig. <ref> that even for relatively small forces (in the range 1.0 ≲ F ≲ 1.7), the behavior is non-monotonic in N as the system is driven strongly out of equilibrium. An effect that does not appear at all for F=0 (compare Fig. <ref>).
While non-monotonic behavior and breaking of linear response is expected for large forces, we see an additional effect. The system transitions from a negative contribution to the variance in the obstacle density to a positive one when the force is increased enough, F ≳ 1.5 [Fig. <ref>]. Such a behavior has been seen in previous works for the lattice Lorentz model <cit.>, and we investigate it further in Sec. <ref>. 











§ SUPERDIFFUSION AT INTERMEDIATE STEPS
 
The intuitive picture is that obstacles suppress the fluctuations in the direction of the force. Our results indicate that for forces large enough, increasing disorder leads to an enhancement.
The step-dependent behavior of the variance can be quantified in more detail by considering the step-dependent diffusion coefficient defined by a discrete step derivative [Fig.<ref>(c)],

    D(N) := 1/2[[(x(N)] - [(x(N-1)] ],

and the local exponent α = α (N) [Fig. <ref>(b)] defined by a discrete logarithmic step derivative,

    α(N):= ln[ [x(N)] ] - ln[ [x(N-1)] ]/ln[ N ] - ln[ N - 1 ].

Thus, ordinary diffusion corresponds to α=1, whereas local subdiffusive and superdiffusive behavior is indicated by α<1 and α>1, respectively.  Transport at strong driving is dominated by a superdiffusive regime which grows with increasing strength of the driving [Fig. <ref>(b)] while the velocity keeps dropping [Fig. <ref>(a)]




The superdiffusion for stronger forces [Fig. <ref>(b)] can be rationalized as follows. While the average velocity keeps dropping with increasing force [Fig. <ref>(a)], there is a disparity between the trajectories of the tracer that have not yet hit an obstacle and follow a free path, with those that are stuck on an obstacle. Initially, this disparity is increased with the number of steps as indicated by a discrete time window of superdiffusion with a growing exponent α > 1 in Fig. <ref>(a). Once we reach discrete time-scales that are much larger than the time of the mean-free path, regular diffusion is recovered but with a significantly larger diffusion coefficient compared to the bare dynamics. This happens as long as the tracer can eventually go around the obstacle. 
For increasing forces, the number of steps it takes to go around an obstacle is increased, which in turn increases this disparity since the trajectory with the free path can cover a larger distance during this time. 

To investigate this effect, we  develop an asymptotic model similar to Ref. <cit.>. At large forces F≫ 1, the tracer's trajectory can be approximated as performing jumps only in the directions of the force in one dimension until it hits an obstacle. Once this happens, the tracer is stuck and stays there. Therefore the asymptotic model is the following: at each step the tracer has a probability to step along the direction of the force and hit an obstacle with probability p = n or move forward with probability 1-n. The probability of displacement x = j after N jumps now reads

    ℙ[x=j | N] =    q^j p+q^j(1-p) δ_j N
    
        =    p+δ_j N[1-(N+1) p]+O(p^2),
          j=0, …, N,    q=1-p.

Where we have approximated to first order in the obstacle density O(p=n) by using the approximation (1-p)^N = 1-N p + O(p^2). The first and second moments can now be calculated,

    ⟨ x(N)⟩   =   ∑_j=0^N j ℙ[x=j|N]
       =   1/2 n N(N+1)+ N[1-(N+1) n]+O(n^2),


    ⟨ x(N)^2⟩   =   ∑_j=0^N j^2 ℙ(x=j|N)
       =   1/6 n N(N+1)(2 N+1) +N^2[1-(N+1) n] 
           +O(n^2).

Hence the variance is

    [x(N)] =   1/2 (1+n) N [1-tanh ^2(F/4)] 
       + nN/6 + nN^2/2 + nN^3/3 + O(n^2),

corrected by the empty lattice (the first term) since for F≫ 1, tanh (F/4) ≈1, and it drops, but in simulations [Fig. <ref>] we still plot the behavior for finite values of F.
The diffusion coefficient from Eq. (<ref>) is now

    D(N) = 1/4(1+n)[ 1- tanh^2 ( F^2/4) ] + 1/2 n (N+1)^2.

Equation (<ref>) suggests that the true exponent of superdiffusion is α=3 as is corroborated by simulations in Fig. <ref>. The reason for the decay of α to unity  at large steps N, is the fact that eventually for any finite forces F the tracer can go around the obstacle if we wait long enough. Therefore, regular diffusion α = 1 is eventually regained for N ≫ 1 but with a considerably higher diffusion coefficient [Fig. <ref>(c)] compared to the bare one defined by D_0=[1 - tanh ^2 (F/4)] / 4.




§ THE DIFFUSION COEFFICIENT IN THE REGIME OF LARGE STEPS 
 
We now determine the diffusion coefficient for large step numbers, N ≫ 1, defined using the variance by D_∞ = lim_N→∞[x(N)]/2N.
We rewrite the MSD in this regime, which was found in Eq. (<ref>), as

    ⟨ x(N)^2 ⟩   =    (1+n) (2 D_0 N + v_0^2 N^2) 
          +n [ 2 v_0 Δ v_∞ N^2 + v_0^2 N^2 + c_1(F) N ],

where D_0 = [ 1 - tanh ^2 (F/4)] / 4 is the diffusion coefficient of the bare dynamics in the direction of F.
Expanding now ⟨ x(N) ⟩^2 from Eq. (<ref>) to first order in the obstacle density n we find,

    ⟨ x(N) ⟩^2    =    [v_0 + n (v_0 + Δ v_∞)]^2 N^2 
       =    [v_0^2 + 2 n v_0(v_0 + Δ v_∞)] N^2 +O(n^2).

The diffusion coefficient is now obtained,

    D_∞   =   [ ⟨ x(N)^2 ⟩ - ⟨ x(N) ⟩^2 ] / N 
       =    (1+n) D_0 + n   c_1(F) /2 + O(n^2).

We plot the case of F=0 in Fig. <ref> where c_1(F=0)= -π / 2. The full function c_1(F) is given in the appendix and in Fig. <ref> we plot D_∞ in the regime of large steps. Higher-order corrections show that there is a complex non-linear dependence on the force, specifically the logarithmic dependence on F via c_1 ≈ -π /2 - (0.33 + 0.11 ln (F)) F^2 + O(F^4), as was also mentioned previously in Sec. <ref>.







§ SUMMARY AND CONCLUSIONS
 
In this work we have considered the driven Lorentz model of a tracer particle hopping on a two-dimensional lattice where a fraction of the sites is inaccessible and acts as hard obstacles.
We derive a method to find the exact solution for the first and second moment of the displacement in z-space [Eq. (<ref>) and Eq. (<ref>)] to first order in the obstacle density n for any value of z (the discrete Laplace transform) and force F.
Our result is correct to first order in the obstacle densities where interactions with single obstacles at a time dominate the process. 
Meaning, obstacles are assumed to be far enough from each other so the effects of trapping by a cluster of obstacles on the tracer can be ignored. 
Considering the mean displacement as a function of steps, ⟨ x(N) ⟩, the first-order expansion in the force obeys linear response in terms of the obstacle density [Eq. (<ref>)], consistent with previous works <cit.>.
The next correction term to the terminal velocity already includes a logarithm, v_∞ = D_x F + (n/16) (π/4 -1) F^3 ln (F) + O(F^3).
Simulations [Fig. <ref>-<ref>] show a good match with theory even for low values of the number of steps.
We note that the point where the theory breaks down depends both on the obstacle density and the magnitude of the force as can be seen the simulation results.
For  larger values of F, the range of validity of our theory shrinks to smaller and smaller obstacle densities. Reversely, the role of clusters becomes more and more relevant, the larger force is.

By switching to continuous time with an exponential waiting time PDF we are able to recover all previous result from <cit.>. 
This is expected since for finite average waiting times τ, the tracer on average performs a jump every τ (arbitrary units), and the connection between N and t is N = t / τ (for large values of N).
It was predicted that for exponential waiting times the obstacles introduce anti-correlations into the system that break the power-law decay ∼ t^-1 towards the terminal velocity expected from linear response <cit.>.
In particular the decay has an exponential dependence on the force and acts as ∼ t^-1exp (-F^2 t / 16) for F → 0 <cit.>.
This is consistent with the domain of discrete time, where each step is distributed according to a delta function, we find ∼ N^-1 exp (-F^2 N / 16) with an onset step-time of N_f ∼ 16/F^2 for small F (see Sec. <ref>).



Investigating the behavior of the variance at an intermediate number of steps shows non-monotonic behavior in N even at small values of force F [Fig <ref>]. For larger forces there is a window of superdiffusion at intermediate values of steps N where the variance behaves as [x(N)] ∝ N^3 [Eq.(<ref>)] until it drops to regular diffusion with a linear dependence on N [Fig. <ref>].
In the lattice Lorentz model, the superdiffusion can be traced back to the rapid increase of the variance of the free path lengths as the tracer performs a purely directed motion along the field until it hits an obstacle.
Although increasing the applied force on the tracer can reduce its travel time between different obstacles, it will increase the time it spends trapped by the obstacles.
Such non-monotonic behavior and superdiffusivity has been found in previous works for the continuous case of exponentially distributed times between steps <cit.>. 
The superdiffusion regime starts immediately for discrete time [Fig. <ref>(b)], 
while in the exponential waiting time case there is still a subdiffusive regime at small number of steps of the order of the density n <cit.>.
Our results show that a superdiffusive regime is a generic feature occurring naturally in the lattice Lorentz gas regardless of the time statistics between steps.

The developed methodology can be used in future research to consider systems where different types of disorder exist, such as non-static obstacles or the quenched temporal disorder that is present in the quenched trap model <cit.>. 


This work was supported by the  Israel Science Foundation Grant No. 2796/20. AS acknowledges FWF Der Wissenschaftsfonds for funding through the Lise-Meitner Fellowship (Grant DOI 10.55776/M3300). TF gratefully acknowledges support by the Austrian Science Fund (FWF) (Grant DOI 10.55776/M3300).




§ THE OBSTACLE-FREE PROPAGATOR Ĝ_0(Z)
 


 §.§ The matrix elements
 
We explicitly evaluate the elements of the matrix Ĝ_0(z) using a representation that makes use of I_m(…), a modified Bessel function of the first kind of integer order m.
We use an integral representation of a fraction, ∫_0^∞ ds e^-s a = 1/a for a>0, to convert Eq. (<ref>) into a  more familiar form


    ⟨𝐫 | Ĝ_0(z) | 𝐫'⟩   =   ∫_-π^πdk_x/2 π e^-i k_x (x-x')
          ×∫_-π^πdk_y/2 π  e^-i k_y (y-y')∫_0^∞ ds e^-s[1-z λ (𝐤)] .


Substituting Eq. (<ref>) into Eq. (<ref>) we obtain 


    ⟨𝐫 | Ĝ_0(z) | 𝐫'⟩ = 
       ∫_0^∞ ds e^-s∫_-π^πdk_x/2 πexp[2sz Γ( cos (k_x) cosh (F/2) 
       + i sin (k_x) sinh (F/2) ) - i k_x ( x - x' ) ] 
       ×∫_-π^πdk_y/2 πexp [2s z Γcos (k_y) -i k_y ( y - y' ) ] .


We then use the relation (<cit.>, Eq. (64))


    ∫_-π^πd k/2 πexp [-i k m] exp [αcos (k)+ i βsin (k)]= 
          [α+β/√(α^2-β^2)]^m I_m(√(α^2-β^2)).


Notice that the positions x and y were defined to be an integer number (between 1 and L) in Sec. <ref>. 
Therefore the matrix elements of the free propagator are


    ⟨𝐫 | Ĝ_0(z) | 𝐫'⟩   =    e^F (x - x') / 2∫_0^∞ d s e^-s
          ×   I_|x-x'|(2 Γ s z ) I_|y-y'|(2 Γ s z) ,


where the exponential outside of the integral accounts for the asymmetry introduced by the bias and the remaining terms are the solution of the symmetrical problem but with the normalization factor set to be Γ(F) instead of Γ(F=0). 


 §.§ Force-free propagators relations
 
We show that the four terms of the free symmetrical propagates, g_00, g_10, g_11, g_20, appearing in the matrix elements of Ĝ_0(z) [see Eq. <ref>],

    ⟨𝐫 | Ĝ_0(z) | 𝐫'⟩ = e^F (x - x') / 2 g_x-x',y-y',

are related to each another.
Thus only the knowledge of two is required to know the other two as well.
The symmetrical case g_xy with no obstacles in Eq. (<ref>) from its definition is also represented in the form <cit.>

    g_xy = ∑_N=0^∞ z^N P_N(𝐫) ,

where 𝐫=(x,y) and just for convenience we switch to the notation that P_N(𝐫) = ⟨𝐫 | p_n ⟩ is the probability to be at position 𝐫 after N steps, given we started at zero. 
The process is translationally invariant and we are free to choose our starting condition to be at the origin, P_0(𝐫) = δ_𝐫,0 
(here δ denotes the Kronecker symbol)
the system obeys a Markovian process 

    P_N+1(𝐫) = ∑_𝐝∈𝒩 W(𝐝) P_N (𝐫 - 𝐝) ,

multiplying by z^N and summing for N from 0 to ∞

    1/z∑_N=1^∞ z^N P_N(𝐫) = ∑_𝐝∈𝒩 W(𝐝) g_𝐫-𝐝 ,

now using P_0(𝐫) = δ _𝐫, 0 we obtain

    g_𝐫 = δ _𝐫, 0 + z ∑_𝐝∈𝒩 W(𝐝) g_𝐫-𝐝.

And using 𝐫 = 0 and 𝐫 = (1,0) we find the relations 

    g_00   = 1 + 4 Γ z g_10
    
            g_10   = Γ z ( g_00 + g_20 + 2 g_11 ).




§ TAUBERIAN THEOREM

The Tauberian Theorem <cit.> allows transitioning to N-space and find ⟨ x(N) ⟩ for large N from the behavior of its generating function ⟨x̃ (z) ⟩ =∑_N=0^∞⟨ x(N) ⟩ z^N when z→ 1.
To be more precise, given that near z → 1 the function behaves as

    ⟨x̃ (z) ⟩∼1/(1-z)^γ Y(1/1-z),

where γ is some positive number and Y(u) is a slowly varying function of u, i.e.,

    lim _u →∞Y(C u)/Y(u)=1,

for any positive constant C.
If the sequence {⟨ x(N) ⟩} is monotonic (at least starting from some value of N) then

    ⟨ x(N) ⟩≅1/Γ(γ) N^γ-1 Y(N),

where Γ is the Gamma-function.
For example the 1/(1-z)^2 term in Eq. (<ref>) is transformed to N  using γ =2 and Y(u)=1 and we obtain


    ⟨ x(N) ⟩   =    (1+n)⟨ x_0(N) ⟩ - n π F/4 N 
       =   (1 - n (π -1))F/4 N 
       =    D_x F N


to first order in the force F and the regime N ≫ 1 correct to first order in the obstacle density n.



§ ASYMPTOTIC BEHAVIOR FOR Z→ 1
 
In this section we provide the exact asymptotics (z → 1) of the scattering matrix for any force F, used to determine the stationary velocity v_∞ and diffusion coefficient D_∞ appearing in the figures.

The first derivative of t_i(z, 𝐤) can be computed with Wolfram Mathematica asymptotically for z → 1 with the result



    .(-i L^2 ∂ t_i(z, 𝐤)/∂ k_x) |_𝐤=0 = 
           -4 (π -2 E) sinh(F/2) [(π -2 E) cosh ^4(F/4)+2 Ksinh ^2(F/4)]/[8 Ecosh ^4(F/4)-K(4 cosh(F/2)+cosh (F)-5)] (-2 Ecosh ^2(F/4)+2 Ksinh ^2(F/4)+π)-2 tanh(F/4) +O(1-z) 
          = Δ v_∞ + O(1-z),



where we abbreviated K = K(^4 (F / 4)) and E = E(^4 (F/4)) and the derivative of t_i(z,𝐤) was evaluated up to terms of order O(1-z).

Similarly, the second derivative of t_i(z, 𝐤) can be expanded using computer algebra for z→ 1 with the result



    .L^2 ∂^2/∂^2 k_xt_i (z, 𝐤) |_𝐤=0 = 
          sech^2(F/4) [U_1 + πK U_2 +2 E U_3 +π ^3 (8 cosh(F/2)+5 cosh (F)-5) cosh ^4(F/4)]/2 π[8 cosh ^4(F/4) E-(4 cosh(F/2)+cosh (F)-5) K] (2 sinh ^2(F/4) K-2 cosh ^2(F/4) E+π)
            -2 tanh ^2(F/4) z/1-z + O(1-z)
           = c_1(F) -2 (v_0)^2 z/1-z +O(1-z).



where


    U_1    =    128 sinh ^2(F/4) cosh ^6(F/4) KE^2,
    
        U_2    =   π(7 cosh(F/2)+cosh(3 F/2)-6 cosh (F)-2)-4 sinh ^4(F/4) (16 cosh(F/2)+3 cosh (F)-3) K ,
    
        U_3    =    6 πsinh ^4(F/2) K-16 sinh ^4(F/4) (cosh(F/2)+3) cosh ^2(F/4) K^2+π ^2 (1-9 cosh (F)) cosh ^4(F/4).


For small 0 < F≪ 1 this can be further simplified to fourth order in the force to

    .L^2 ∂^2/∂^2 k_x t_i (z, 𝐤) |_k⃗=0=c_1(F) - (F^2/8) z/1-z + O(1-z).

where c_1(F) in the regime of small F takes the form,

    
    c_1    =   F^2/64 π (π -2){ 2 [π ^2 (π -6)+16] log (F)-112 log (2)+π{16+π  [22+42 log (2)-π  (13+log (128))]}}
           -π/2 + O(F^4).





§ SIMULATIONS
 
In this section we provide details on the method of simulations.
Fluctuations in the bare dynamics (without obstacles) are much stronger than the obstacle induced ones, especially when the number of steps is small. 

Therefore to construct the simulations, it is very useful to adapt the approach from ref. <cit.> which was also successfully adapted in ref. <cit.>.
Let us denote the difference in position of a test particle for the same sequence of trial moves with and without obstacles (excluded sites) by δ x. Clearly, the average magnitude of δ x is proportional to the density of obstacles. Let us write the total position of a particle in the presence of hard obstacles as 

    x(N) = x_0(N) + δ x.

⟨δ x ⟩ can be computed directly from simulations and converges much faster than ⟨ x(N) ⟩ while x_0(N) is the position of the bare-dynamics (no obstacles) which is known. 
To further speed up simulations in the small-N regime, we use a second trick. We simulate obstacle positions only in the effective region of the trajectory. Meaning, obstacles only appear in the area the test particle can reach up to step number N. We then normalize the result to account for the shifted dynamics by using the fact that for obstacle outside the effective region, the mean position would just be that of the bare dynamics x_0(N). 
Same methodology were used to determine x(N)^2.









