Time-dependent dynamics of the three-dimensional driven lattice Lorentz gas

We consider the motion of a single tracer particle on a three-dimensional lattice in the presence of hard, immobile obstacles at low density. Starting from an equilibrium state, a constant pulling force on the tracer particle is switched on. We elaborate a complete solution for the dynamics, exact in first order of the obstacle density. The time-dependent velocity response and the fluctuations along and perpendicular to the force are derived and compared to stochastic simulations. The non-analytic dependence of the linear-response functions on the frequency, reflecting the long-time tail of the velocity autocorrelation function in equilibrium, has a non-equilibrium analogue in the non-analytic dependence of the stationary velocity as well as the diffusion coefficients as a function of the force. We show that a divergent time scale emerges, separating a regime where linear response is qualitatively correct from a regime where nonlinear effects dominate.


Introduction
The dynamical properties of systems in thermal equilibrium are encoded in time-dependent correlation functions of fluctuating observables. By the fluctuation-dissipation theorem, the correlation functions also encode the linear response to a time-dependent external perturbation. However, often the perturbations are not small and the system is driven strongly out of equilibrium.
In systems of complex transport the equilibrium dynamics exhibits persistent correlations in the form of slow algebraic decays of the correlations functions connected to transport coefficients. In particular, the velocity-autocorrelation function, which encodes the time-dependent self-diffusion, displays a t −d/2 algebraic decay due to slow diffusion of transverse momentum in fluids in d dimensions [1][2][3][4][5].
Lattice variants of Lorentz systems, where a continuous-time random walker explores a percolating cluster of unblocked lattice sites, exhibit the same phenomena and offer the prospect for analytic progress. In particular, the velocity-autocorrelation function of the lattice Lorentz model was calculated for all times in first order of the density [16]. Thus, by the fluctuationdissipation theorem, the linear response upon driving is known completely for all times.
Beyond linear response, virtually no exact solutions are available for strongly interacting systems. Typically even the stationary state of such systems is unknown and the stationary values of observables can not be calculated. However, a well-defined experimental protocol for such a nonequilibrium setup is provided by switching on a strong perturbation driving the initial equilibrium state to a new stationary one. In particular, one can gain fundamental insight into the system's properties probing suitable observables as time evolves.
There has been recent progress on the driven lattice Lorentz system in two dimensions, where the equilibrium dynamics in first order of the density [16][17][18] has been generalized to an exact analytic solution for the case of a force pulling on the tracer particle [19][20][21]. Furthermore, the complementary case of a biased intruder in a dense environment of mobile particles on a lattice has also been considered [22][23][24].
Here, we extend the time-dependent dynamics of the tracer on the lattice to the threedimensional case for arbitrary strength of the force and for all times. In particular, we will make heavy use of the recently elaborated lattice Green functions for the three-dimensional case, which will serve as input in our calculation. The problem is solved in terms of scattering theory of the tracer with a single obstacle, which encodes the complete response in first order of the density of obstacles. We derive formal expressions for the moment-generating function for the displacements of the tracer and specialize to the time-dependent velocity response and the time-dependent diffusion coefficients parallel and perpendicular to the applied force.

Model and solution strategy
We consider a tracer particle performing a random walk on a cubic lattice. We set the lattice spacing to unity and define the accessible sites by Λ = {r = (x, y, z) ∈ Z 3 : x, y, z ∈ [−L/2, L/2 [}, where L ∈ N denotes the linear size of the lattice with N = L 3 sites in total. We employ periodic boundary conditions and anticipate the limit of large lattices L → ∞. The dynamics of the tracer is governed by successive jumps to its nearest-neighbors N = {±e x , ±e y , ±e z } and the waiting time of the tracer at every site is exponentially distributed with mean waiting time τ.

Empty lattice
To fix the notation, we recall the case where every site is accessible to the tracer [25]. Then, the random walk is described in terms of the conditional probability U 0 (r, t|r ) of the tracer to move from site r to site r in lag time t. The time evolution of the conditional probability U 0 (r, t|r ) is governed by a master equation of the form with transition probabilities W(d) and dimensionless rate Γ. We are interested in the following scenario: The system is initially in equilibrium and at time t = 0 a force on the tracer is switched on, driving the system into a new stationary state. The force is aligned along the x direction of the lattice with dimensionless strength F = force · (lattice spacing)/k B T . Local detailed balance, W(e x )/W(−e x ) = e F and W(e y )/W(−e y ) = W(e z )/W(−e z ) = 1, are satisfied by the following form for the transition probabilities: W(±e x ) = e ±F/2 /(4 + e F/2 + e −F/2 ) and W(±e y ) = W(±e z ) = 1/(4 + e F/2 + e −F/2 ), and we consider non-normalized rates by defining the dimensionless rate via Γ = [cosh(F/2) + 2]/3. This choice is in accordance with previous work on the driven lattice Lorentz gas in two dimensions [19,20], but also different choices are conceivable [24,26]. Let us emphasize that using normalized rates, while keeping the relative weights of the transition probabilites, changes only the time scale in a forcedependent manner such that all observables can be trivially rescaled. It is advantageous to employ a bra-ket notation such that the conditional probability is written as U 0 (r, t|r ) = r|Û 0 (t)|r in the complete orthonormal basis of position kets {|r : r ∈ Λ} with time-evolution operator Û 0 (t). The master equation for the time-evolution operator can then be written as with initial condition Û 0 (t = 0) = and the dynamics is driven by the unperturbed 'hamiltonian' The master equation (equation (2)) is formally solved by Û 0 (t) = exp(Ĥ 0 t) which can be explicitly calculated by considering the plane-wave basis with wave vector k = (k x , k y , k z ) ∈ Λ * := {(2πx/L, 2πy/L, 2πz/L) : (x, y, z) ∈ Λ} and scalar product k · r = k x x + k y y + k z z. Since the unperturbed hamiltonian Ĥ 0 is translationally invariant, it is diagonal in the plane-wave basis (equation (4)) via where δ(k, k ) denotes the Kronecker delta and the eigenvalue (k) can be calculated to The equilibrium state |p eq := (1/N) r∈Λ |r consists of equally likely lattice sites, and we denote averages in the absence of obstacles by · 0 . The intermediate scattering function 0 encodes all moments of the displacement of the tracer and can be evaluated to Thus, for the mean-displacement along the force in x direction, we obtain with bare velocity v 0 = sinh(F/2)/3τ . Similarly, for the mean-square displacement along the three directions, we calculate with bare diffusion coefficients D 0 x = cosh(F/2)/6τ and D 0 y = 1/6τ . The explicit form of the propagator in the basis of lattice sites |r is known and can be obtained from the conditional probability in the plane-wave basis: Then, in the limit of large lattices, L → ∞, the sums can be replaced by integrals over the first Brillouin zone BZ : Let us consider first the k x -integral that connects the points ∓π on the real line. Now we deform the contour to the path in the complex plane connecting first −π to −π − iF/2 by a straight line, then −π − iF/2 to π − iF/2, and last π − iF/2 to π. Then by exp(±iπ) = −1 one observes that the two paths parallel to the imaginary axis cancel, and we are left with (13) The three remaining integrals can be identified as integral representations of the modified Bessel function I m (·) of integer order m: Thus, the conditional probability in the site basis is obtained as with the respective displacements ∆x = x − x , ∆y = y − y , and ∆z = z − z . Formally the conditional probability in the Laplace domain can be written as a resolvent of the hamiltonian Ĥ 0 : In the basis of position kets, they are known as lattice Green functions with The effect of the external force is two-fold: First, there is an exponential prefactor accounting for the displacement of the tracer along the force. Second, the force leads to a nontrivial shift in the dimensionless rate Γ = 1. These observations allow expressing the full propagator in terms of the unbiased propagator in equilibrium, yet the force leads to a shift in the frequency s →s := s + (Γ − 1)/τ . Thus, the dynamics of the tracer is encoded in the propagator in equilibrium, ĝ 0 , depending on the shifted frequency s: Explicit expressions in term of elliptic integrals have been elaborated in a series of pioneering works by Joyce [27][28][29][30] (see appendices A and B for a summary).

Obstacle obstructed motion
Next, we consider the presence of hard, immobile obstacles on the lattice. When the tracer attempts to jump onto an obstacle site, the tracer remains at its initial site before the jump. Thus, the obstacle cancels transitions to and from the obstacle site. Formally, the presence of obstacles can be accounted for by considering the 'hamiltonian' Ĥ =Ĥ 0 +V where the 'potential' V contains information about the positions of the obstacles with the respective change of the transition rates. We denote the potential for a single impurity at site s i by v(s i ) ≡v i . Thus, for N I obstacles on the lattice, the total potential V can be written as a sum over all single-obstacle potentials v i with V = NI i=1v i . The obstacles are distributed randomly on the lattice and we perform an average over disorder realizations at the end of the calculation.
Since the tracer performs nearest-neighbor jumps, the only nonvanishing matrix elements for an obstacle located at the origin are the obstacle site e 0 = 0 and its nearest neighbors e ±3 = ±e z , e ±2 = ±e y , e ±1 = ±e x . Then, the matrix representation of the single obstacle potential e µ |v|e ν in the distinguished subspace (e −3 , e −2 , e −1 , e 0 , e 1 , e 2 , With this representation, we can follow the approach outlined in [21] for the square lattice. We express the dynamics of the tracer particle in the presence of obstacles in terms of the propagator Ĝ witĥ and employ the formalism of quantum mechanical scattering theory such that the propagator can be equivalently expressed aŝ with the scattering operator T =V +VĜ 0T . Here, the scattering operator encodes all possible sequences of collisions of the tracer with the obstacle disorder. Analytic progress is made by rearranging the terms in the multiple scattering expansion with the scattering-operators for a single obstacle at site s i : After performing a disorder average [·] av the problem again becomes translationally invariant and the propagator [G] av (k) := k|[Ĝ] av |k is diagonal in the plane-wave basis. Averaging equation (21) requires the matrix element of the average scattering operator k|[T] av |k . From the multiple scattering expansion (equation (22)), one infers that to first order in the obstacle density n = N I /N, only the forward-scattering amplitude t(k) := k|t|k of a single impurity contributes. For definiteness, we place this obstacle at the origin of the lattice. Thus we arrive at In our derivation so far, the tracers are allowed to start at an impurity site. Since we are interested in tracers that start on accessible sites, we multiply the disorder-averaged propagator with 1/(1 − n) = 1 + n + O(n 2 ) and collect again in powers of the obstacle density. Thus, we obtain the disorder-averaged propagator The disorder-averaged propagator encodes all moments of the displacement of the tracer particle in first order of the density and for arbitrarily strong driving.

Calculation of the forward-scattering amplitude
The scattering operator t for single obstacle v fulfills the relation accounting for repeated scattering with an obstacle located at the origin. Since the potential v displays non-vanishing contributions only in the distinguished subspace (e −3 , . . . , e 3 ), it is sufficient to consider equation (26) in the distinguished subspace only. One readily checks that (1, 1, 1, 1, 1, 1, 1)/ √ 7 is a left eigenvector to eigenvalue 0 for v (equation (19)) which reflects that v merely cancels the transition rates due to the presence of obstacles, thus conserving probability. From equation (26) one infers that the scattering matrix t displays the same left eigenvector also with eigenvalue 0.
For the calculation of the scattering matrix t for a single obstacle (equation (26)), only the matrix representation of the free propagator in the distinguished subspace is required: g 11 e F/2 g 11 g 10 e −F/2 g 11 g 11 g 20 g 11 g 00 e F/2 g 11 g 10 e −F/2 g 11 g 20 g 11 e −F/2 g 11 e −F/2 g 11 g 00 e −F/2 g 10 e −F g 20 e −F/2 g 11 e −F/2 g 11 g 10 g 10 e F/2 g 10 g 00 e −F/2 g 10 g 10 g 10 e F/2 g 11 e F/2 g 11 e F g 20 e F/2 g 10 g 00 e F/2 g 11 e F/2 g 11 g 11 g 20 e F/2 g 11 g 10 e −F/2 g 11 g 00 g 11 g 20 g 11 e F/2 g 11 For this representation, we used the observation of section 2.1 that the free propagator can be expressed in terms of the unbiased propagator, except for an exponential force bias and a shift in the frequency. Here, g 00 (s) := 0|ĝ 0 (s)|0 is the unbiased propagator for zero steps, In the case of no force on the tracer, the problem displays cubic symmetry, such that the calculation of the scattering matrix can be simplified by considering a transformation to symmetry adapted modes. In this work, we use the following basis transformation and we refer to the respective column vectors by the ordered set (|p z , |p y , |d y 2 −z 2 , |p x , |d x 2 , |s , |n ). The 3 p-modes are antisymmetric along each of the coordinate axis, the d-modes are source free quadrupolar modes, whereas the s-mode and n-mode are symmetric with respect to all orientations. Note that n| is left-eigenvector of v and t with eigenvalue 0, reflecting conservation of probability. Although the force breaks the cubic symmetry of the problem, is still advantageous to use the transformation to the new basis.
In the new basis the perturbation v = M T vM and the Green function In particular one infers that |p z , |p y , and |d y 2 −z 2 are still eigenmodes of the problem, reflecting the residual symmetry. In the transformed potential v the last row vanishes, reflecting the conservation of probability encoded in the 'neutral' mode |n . In the new basis, solving equation (26) reduces to inverting a 3 × 3 matrix problem in the subspace spanned by |p x , |d x 2 , and |s . With the help of the symmetry considerations, the three-dimensional lattice Lorentz problem is of the same complexity as its two-dimensional analogue [21], and can be implemented efficiently by computer algebra.
Then the forward-scattering element of the t-matrix t(k) = k|t|k can be readily evaluated in the transformed basis. In particular, for the matrix elements of the scattering t-matrix in the new basis which are necessary for the time-dependent velocity response and the timedependent diffusion coefficients, we obtain: The derivation of these expressions is along similar lines as for the two-dimensional case [21].

Moments of the time-dependent dynamics
In this section we elaborate the first and second cumulants of the time-dependent motion of the tracer particle upon switching on the external force. The calculation is based on taking derivatives of the disorder-averaged propagator with respect to the wave vector. We provide analytic expressions in the Laplace domain and derive the stationary behavior and expansions in terms of the force. Furthermore, we compare the numerical results to stochastic simulations at low but finite densities.

Time-dependent velocity
The mean-displacement along the force in the Laplace domain is obtained by taking a derivative of the disorder-averaged propagator with respect to the x component of the wave vector k: In the second line of equation (34), we used G 0 (k = 0) = 1/s, i∂G 0 /∂k x | k=0 = v 0 /s 2 with bare velocity v 0 = sinh(F/2)/3τ , and the relation for ∂Nt(k)/∂k x | k=0 given in equation (30). Thus, the frequency-dependent velocity v(s) of the tracer is obtained aŝ The behavior for small times t → 0 is obtained from the high-frequency behavior via With increasing time, the time-dependent velocity of the tracer, v(t), becomes more and more suppressed and approaches a stationary value for large times (figure 1). For the lower density n = 10 −3 , small deviations of the analytic result to the stochastic simulations become apparent only at the largest force (F = 10, where the bias in the rates is W(e x )/W(−e x ) = e 10 ≈ 2 · 10 4 ).
For the larger density n = 10 −2 , the analytic theory still captures forces up to F = 6 (bias e 6 ≈ 400). For larger forces, we still find qualitative agreement. The stationary velocity is obtained from the low-frequency behavior by taking the limit v(t → ∞) = lim s→0 sv(s) and can be written in the following way: where the free propagators g 00 and g 20 have to be evaluated at shifted frequency s = (Γ − 1)/τ . For small driving, the stationary velocity exhibits nonanalytic behavior of the form where the coefficients are given by  (38) is accurate for forces F 4, the exponential behavior for large forces is followed for F 6.
The approach to the stationary state can be inferred from the small-frequency expansion of equation (38). Since the scattering matrix t encodes the complete approach and depends on the shifted frequency s →s = s + (Γ − 1)/τ, it is natural to use an expansion in small shifted frequencies and forces: with coefficients In particular, for the stationary case s → (Γ − 1)/τ = F 2 /24τ + O(F 4 ), we recover the small-force expansion of the stationary velocity (equation (38)  shifted frequency, s → s and the first three terms in equation (42) reproduce the well-known singular low-frequency expansion of the velocity-autocorrelation function [32]. The expansion (equation (42)) reveals that the velocity response is nonanalytic in the shifted frequency s, which implies that e (Γ−1)t/τ dv(t)/dt displays a singular Laplace transform in s.
and applying the Tauberian theorems one finds from the term B v,3/2 F(sτ ) 3/2 the time-dependent approach of the velocity to its stationary state: For infinitesimal driving, O(F), we recover the algebraic tail t −3/2 as follows from the linear response framework with the velocity-autocorrelation function as input [32]. Yet, there is a second nonanalytic Direct comparison shows, that the second contribution becomes relevant on time scales t τ F = 24τ /F 2 . Therefore, there is a now divergent time scale in the problem of driven systems separating two different regimes. In particular, the limits F → 0 and t → ∞ do not commute: for t τ F linear response qualitatively applies, whereas for t τ F the nonequilibrium driving manifests itself in a qualitatively new behavior.
The leading correction due to driving is the onset of an exponential decay which is illustrated in figure 3 and confirmed by stochastic simulations. Therefore, the approach to the stationary state is exponential for any finite driving in striking contrast to the algebraic approach predicted by linear response. In particular, for small forces the velocity response follows for longer and longer times the predictions of linear response, however, for times larger than τ F , the exponential decay dominates and linear response becomes qualitatively wrong.
Our analysis reveals that the nonanalytic behavior of the stationary velocity (equation (38)) and the long-time tails of the approach to the stationary state are intimately connected. Both originate from a singular low-frequency expansion in the shifted frequency as induced by the bias, thereby introducing two different temporal regimes.

Fluctuations perpendicular to the force
Due to symmetry, the mean displacement perpendicular to the force vanishes. Therefore, the first nontrivial cumulant is given by the mean-square displacement ∆y(t) 2 = ∆z(t) 2 . Equivalently, the time-dependence can be discussed in terms of the time-dependent diffusion coefficient defined by For long times, we anticipate diffusive behavior such that D y := D y (t → ∞) becomes the diffusion coefficient in the stationary state. From the disorder-averaged propagator, the mean-square displacement perpendicular to the force in the frequency domain is obtained as with −∂ 2 G 0 /∂k 2 y | k=0 = 2D 0 y /s 2 and bare diffusion coefficient D 0 y = 1/6τ . The frequencydependent diffusion coefficient perpendicular to the force is then obtained from equation (49) by the rules for the Laplace transform of a derivative: The time-dependent behavior found in computer simulations is well described by the analytic solution ( figure 4). For all forces, the stationary diffusion coefficient D y (t → ∞) is smaller than the bare diffusion coefficient D 0 y , hence the disorder always suppresses the fluctuations perpendicular to the force. Nevertheless, the time-dependent approach to the long-time limit displays an interesting non-monotonic behavior for all nonvanishing forces. This is in striking contrast to the equilibrium case, where the approach to the stationary-state diffusion  (46)). Force F = 0 refers to the limiting case of vanishing force as inferred from the linear response prediction. is predicted to be completely monotone. Let us emphasize again that this feature does not depend on whether normalized or non-normalized rates are employed in the model definition, since the effects are measured relative to the bare lateral diffusion coefficient.
The stationary state diffusion coefficient D y (t → ∞) = lim s→0 sD y (s) can be elaborated analytically to where the free propagators g 00 and g 20 have to be evaluated at shifted frequency s = (Γ − 1)/τ .
Similarly to the stationary velocity, the stationary diffusion coefficient perpendicular to the force is nonanalytic and is given by the following expansion: For vanishing force, we recover the equilibrium diffusion coefficient D eq y := D 0 y + nA y,0 /τ , and consistent with the linear response result A v,1 = A y,0 . The leading nonanalytic contributions arise at order O(|F| 3 ). For large forces, the diffusion coefficient approaches D y (t → ∞) → D 0 y − n/8τ. Interestingly, the obstacles still suppress the perpend icular motion at very strong forces, although less than in the equilibrium case.
The analytic solution for the density-induced diffusion coefficient D y (t → ∞) − D eq y is corroborated by computer simulations as shown in figure 5. For large forces, the density-induced diffusion coefficient saturates at finite values τ [D y (t → ∞) − D eq y ]/n = −1/8 − A y,0 ≈ 0.130 19. The long-time behavior of the time-dependent diffusion coefficient is obtained from the joint expansion for small shifted frequencies and small forces: 3/2 , and an additional new coefficient Specializing equation (56) to the stationary state, yields equation (52) and reveals the additional identity A y,3 √ 24 = B y,1/2 + B y,3/2 /24. The first three terms in equation (56) represent the well-known low-frequency expansion in equilibrium, F = 0, including the leading nonanalytic contribution. By the fluctuation-dissipation theorem, it is identical (up to a factor of F) to the linear response of the time-dependent velocity (equation (42)).
Similar to the time-dependent velocity response (equation (46)), the time-dependent approach to the stationary diffusion coefficient D y (t → ∞) is a power law decorated by an exponential: which holds in the linear response regime t τ F . Yet, there is a second nonanalytic contribution B y,1/2 F 2 (sτ ) 1/2 leading to in the nonlinear regime t τ F . Interestingly, the two regimes are characterized by a different sign, which explains the nonmonotonic approach to the stationary diffusion coefficient D y (t → ∞) in figure 4.

Fluctuations parallel to the force
Parallel to the force, the next higher cumulant apart from the mean ∆x(t) is given by the variance or equivalently the time-dependent diffusion coefficient which both encode the fluctuations. The new ingredient is the mean-square displacement ∆x(t) 2 parallel to force. It is obtained in the frequency domain from the disorder-averaged propagator via . Both the mean-square displacement and the mean displacement squared in equation (60) grow ∼t 2 for long times, whereas the variance is anticipated to grow only linearly in time. To show that the leading contributions of both terms indeed cancel, we split off the bare contribution from the density-induced contribution in both quantities. For the mean-displacement, we write with the velocity response ∆v(t). Then the square of the mean displacement in first order of the density is given by where we discard contributions in second order of the density. Similarly for the mean-square displacement parallel to the force, we introduce the first-order density response ∆R x (t) and write Thus, in the frequency domain, we obtain for the variance: and correspondingly for the diffusion coefficient parallel to the force: For small frequencies, the asymptotic behavior of the response functions is provided by From the analytic solution, it can be shown that ∆R x,−3 = 4v 0 ∆v −1 such that the long-time behavior parallel to the force is diffusive for all forces with diffusion coefficient We refrain from giving an explicit expression for the stationary diffusion coefficient D x (t → ∞), since the formulas become lengthy. However, they can all be generated from taking appropriate derivatives of the free lattice Green functions. The time-dependent behavior of the diffusion coefficient parallel to the applied force is corroborated by computer simulations ( figure 6). For small forces we observe a qualitatively similar behavior as for the diffusion coefficient perpendicular to the force. In particular, the diffusion coefficient first decreases until a point of least diffusivity is reached, followed by an increase to its stationary value. For larger forces, the stationary diffusion coefficient strongly exceeds the short-time diffusion coefficient D 0 x (1 − n). In between there is a strong increase of the time-dependent diffusion coefficient which we refer to as transient superdiffusion.
Asymptotically for very strong forces, we anticipate that the window of superdiffusion becomes larger and larger and approaches a power-law increase of D x (t) ∼ t 2 . The chain of arguments can be taken literally from the two-dimensional analogue [20], where in the time regime considered, tracer particles move only along the force until they hit an obstacle, whereas the perpendicular motion can be ignored. Then, one can show by elementary methods, that the variance grows as Var x (t) ∼ t 3 by the fluctuations of the mean-free path lengths of the tracer. The data displayed in figure 6 are not in this asymptotic regime yet, although the diffusion coefficient increases already a factor of 3.5 at density n = 10 −3 , respectively by a decade for density n = 10 −2 .
For small forces, the stationary diffusion coefficient parallel to the force exhibits densityinduced nonanalytic contributions of the form with coefficients In equilibrium, diffusion along all lattice directions is identical which leads to A x,0 = A y,0 (equation (53)). However, already the leading correction, A x,2 F 2 , as well as the leading nonanalytic contribution, A x,3 |F| 3 , for the diffusion along the force differ from the corresponding expression perpendicular to the force. The long-time behavior of the time-dependent diffusion coefficient is obtained from the joint expansion for small shifted frequencies and small forces: with 3/2 , and additional new coefficient The joint low-force and frequency expansion is similar to the one for the perpendicular direction (equation (76)), except for the emergence of a term of order O(F 4 /s 1/2 ). Again, linear response is qualitatively valid for times t τ F , whereas in the opposite case t τ F the nonlinear regime sets in. Specializing equation (76) to the stationary case reveals the additional relation A x,3 = B x,3/2 /24 3/2 + B x,−1/2 /24 −1/2 + B x,1/2 /24 1/2 . Similar to the perpend icular diffusion coefficient, the nonmonotonic approach of the diffusion coefficient parallel to the force shown in figure 6 is explained by the different signs of the coefficients B x,3/2 and B x,1/2 in the expansion (equation (76)).
For large forces, the force-dependent behavior can be elaborated to Rather than discussing the stationary diffusion coefficient D x (t → ∞) we focus on the forceinduced diffusion coefficient as defined by with the stationary state diffusion coefficient in equilibrium, D eq x (t → ∞) := D x (t → ∞, F = 0). In the lattice there is already a strong exponential enhancement of the diffusion coefficient due to the force without obstacles. The force-induced diffusion coefficient relative to the empty lattice is compared to computer simulations in figure 7, revealing two trend reversals. For small forces, the obstacles lead to faster growing diffusion coefficient, while intermediate forces lead to a slower growth than the reference case. For very large forces we observe the purely density-induced exponential growth (equation (78)).

Summary and conclusion
We have solved for the dynamics of a tracer on a three-dimensional lattice in the presence of fixed, hard obstacles subject to a driving force along a lattice direction, switched on at time zero. The solution strategy has been taken over from the corresponding two-dimensional analogue essentially without change. The main insight is that the residual symmetry of the problem with the force is the same as in the two-dimensional case, thus requiring to solve for a 3 × 3 matrix problem. A necessary ingredient for the analytic solution are the unperturbed lattice Green functions which have been calculated only recently.
In particular, we have elaborated the time-dependent velocity response of the tracer particle for arbitrarily strong forces in first order of the obstacle density. Analytic expressions are provided for the stationary velocity, revealing nonanalytic dependences for small forces, which are a fingerprint of the long-time tails found in equilibrium. The approach to the stationary state is exponentially fast for any finite force, in striking contrast to the linear response prediction.
For the fluctuations of the tracer perpendicular to the applied force, the corrections induced by the obstacles remain small for all times and are of the order of the obstacle density itself. The nonanalytic behavior in the force for the stationary diffusion coefficient is a nonequilibrium analogue of the singular low-frequency expansion in equilibrium. The approach to the stationary state is nonmonotonic which is a manifestation of two distinct regimes separated by a diverging time scale. For times smaller than the crossover time, the equilibrium prediction is qualitatively correct, while for larger times, a new nonequilibrium regime is entered.
The situation is rather different for the fluctuations along the force. Here, the timedependent diffusion coefficient exhibits an intermediate superdiffusive behavior resulting in a strongly enhanced stationary diffusion coefficient for strong driving. This long-time diffusion coefficient grows even exponentially for strong forces.
The analytic solution is exact in first order of the obstacle density and for any applied force on the tracer. However, the range of validity of the theory is density-dependent, such that for increasing force, smaller densities have to be considered. In particular, the limit obstacle density n → 0 and the limit of large forces F → ∞ do not commute. Our calculation is always in the regime, that first the density goes to zero before the forces become large. In the opposite case of fixed obstacle density and increasing force, it has been demonstrated rigorously that there is always a critical force, such that the velocity of the tracer vanishes [34][35][36]. In particular, for such strong forces, the differential mobility always becomes negative as has been discussed in the physics literature [19,37,38]. This does not imply that the tracer is localized, rather it moves sublinearly in time which gives rise to aging. Here, the intuitive picture is that the tracer bumps into transiently trapping obstacle configurations which can be escaped only by moving against the field. Deeper and deeper traps are found as time proceeds, yielding a slowing down of the particle, reminiscent of a continuous-time random walk with power-law distributed waiting times [39][40][41][42].
The qualitative behavior of the three-dimensional lattice Lorentz gas agrees with its twodimensional analogue. The main difference is in the nonanalytic contributions for the longtime tails in equilibrium and the nonanalytic dependences in the force of the stationary-state behavior. Explicitly, in three dimensions, the velocity-autocorrelation function in equilibrium decays as ∼t −5/2 leading to a nonanalytic contribution of the velocity response on the force O(F 3 |F|), whereas in two dimensions already the leading correction is nonanalytic O(F 3 ln |F|) arising from the long-time tail ∼t −2 . The solution strategy can be adapted to the d-dimensional hypercubic lattice while the residual symmetry still requires a solution for the relevant 3 × 3 matrix problem. The asymptotic expansions for the respective lattice Green functions are known [43] and display small-frequency singularities due to the algebraic decay of the return probability. The intimate relation between the long-time tail ∼t −(d+2)/2 and the small-force behavior in the stationary state is anticipated to hold in any dimension d 2.
Effects of confinement can be mimicked by considering the dynamics of the tracer on a lattice with periodic boundary conditions and a finite extension such that the dynamics can be viewed as being on a torus. Such a strategy has been adapted for the case of a biased intruder in a dense bath of mobile particles [22,24]. In essence, the lattice Green functions involve sums rather than integrals in wave-number space in the compact directions of the lattice.
The phenomena discussed here are generic and are not restricted to the special properties of the lattice. In particular, the nonanalytic contributions in the equilibrium frequency-dependent diffusion coefficient and in the velocity response as a function of the force should persist for any density below the percolation threshold. Similarly, the fact that the obstacles are frozen and hard is not an essential ingredient for the behavior discussed and should hold also for mobile and/or soft obstacles.
The phenomenology takes over also to the continuum case for a colloidal suspension where a single probe particle is driven by an external force. For the equilibrium case, the low-density solution has been elaborated earlier [14,15,44]. In particular, the velocity-autocorrelation function displays a long-time tail due to repeated scattering events of the probe particle with the same collision partner. For the driven case, the stationary-state velocity as well as the diffusion coefficients as a function of the driving display nonanalytic behavior [45,46] in qualitative agreement with the results presented in our work. However, the time-dependent behavior, in particular the emergence of a divergent time scale separating the linear response from the nonlinear regime has not been elaborated yet.

Acknowledgments
This work has been supported by the DFG research unit FOR1394 'Nonlinear response to probe vitrification' and by the Austrian Science Fund (FWF) under grant I 2887-N27. TF gratefully acknowledges the hospitality of Université Pierre et Marie Curie where parts of the project have been performed.

Appendix A. Propagators for the three-dimensional case
For explicit expressions of the time-dependent dynamics of the tracer, we have to insert explicit expressions for the four propagators g 00 , g 10 , g 11 , and g 20 . The propagators are not independent and are connected by the resolvent equation 0|(s −Ĥ 0 )Ĝ 0 = 0|, (A.1) providing the two relations (sτ + 1)g 00 − g 10 = τ , 6(sτ + 1)g 10 − (4g 11 + g 00 + g 20 ) = 0, (A.2) which allows expressing g 10 and g 11 in terms of the propagators g 00 and g 20 . Explicit expressions for both propagators have been derived by Joyce in terms of complete elliptic integrals [30]. Starting with g 00 ,