Driven tracer with absolute negative mobility

Take a system in an equilibrium or nonequilibrium steady state and apply a small drive to it. We expect the system to move in the direction of the drive, and increasingly so with stronger drives. However, various setups have been found where intuition is violated and the response is more surprising. Among the simplest manifestations of this behavior are negative differential mobility (NDM), where the response coefficient depends on the applied perturbation in a nonmonotonic way, and absolute negative mobility (ANM), where the sign of the response coefficient is opposite to what would intuitively be expected.

Theoretical examples of NDM include uniformly driven systems [1–6] or single driven tracers in quiescent media [7–40]. Also condensed matter experiments have been performed [1, 41–44]. The appearance of NDM mostly relies on trapping mechanisms that can be implemented through e.g. complicated potentials [7, 11, 20, 21, 43] or impurities, either present by definition [9, 10, 12–15, 18] or effectively created by a slow relaxation of the surrounding medium [5, 16, 17, 19]. A pedagogical explanation for NDM is given in [9], and a modified Green–Kubo formula that accounts for NDM has been proposed in [15].

Absolute Negative Mobility has been observed in a variety of setups. Typically, one does not expect ANM to take place when the unperturbed system is in equilibrium, since, as has been argued, this would constitute a violation of the Second Law of Thermodynamics [45–47]. Thus, previous studies demonstrating ANM considered a driving field which acts on nonequilibrium steady states. These include systems with a periodic [48–58] or a random [57, 59–62] drive, random walkers [20, 45–47], strong interactions and noise in spatially periodic potentials [63–67], and others [68–71]. A different setup where ANM has been found both experimentally and theoretically involves quantum mechanical effects such as absolute negative conductivity for semiconductors, where negative conductivity is associated with a negative effective mass of the carriers (either electrons or holes) [52, 72, 73], or interactions between light and matter [51, 74–77].

In the present work we focus on cases, where the unperturbed system is in thermal equilibrium, and demonstrate that, depending on the drive mechanism, ANM can take place in such systems. This is done in the context of a model of a tracer moving on a discrete ring populated by neutral particles, which obey a simple symmetric exclusion process (SSEP) -type dynamics. We show that, when a driving force is applied to the tracer, it moves in a direction opposite to the drive. We then consider a continuum analogue of the model by studying the Langevin-type motion of a tracer particle in a narrow channel of a gas of hard disks. Here we find that the model exhibits NDM but not ANM.

In the discrete ring model introduced in the present work, the dynamics of the tracer is such that it can move by two different processes, either by hopping towards a neighboring vacant site or by exchanging its position with close bath particles. The exchange move requires enough 'free volume' in the vicinity, which places restrictions on the dynamics and is reminiscent of kinetically-constrained models [5, 12, 13, 78]. The system is studied analytically within a mean-field approximation and numerically, and the existence of ANM is demonstrated. The model may be considered as a kind of toy model for hard disks performing Brownian motion dynamics in a narrow planar channel. Therefore, we carried out Langevin-type simulation studies of this hard disks model where NDM but no ANM has been found. But we believe that some simple variants of this model, which have not yet been tested, could exhibit ANM for selected sets of parameters.

The paper is organized as follows: in section 2 we study the lattice model analytically and numerically. In section 3 we introduce the model of hard disks moving in a narrow channel and present the results of molecular dynamics simulations. In section 4 we conclude with a discussion of the results.

2.1. Definition

Consider a set of N bath particles and a tracer particle occupying L sites of a ring of length L while satisfying the simple exclusion constraint. Time is continuous and each transition occurs with a probability $\newcommand{\dd}{{\rm d}} R \dd t$ during each infinitesimal time step $\newcommand{\dd}{{\rm d}} \dd t$, where R is the rate of the transition. The bath particles are regular SSEP particles and can hop towards the site directly to their left or to their right, each with constant rate 1, under the condition that the target site is vacant. Their average density is denoted by $\newcommand{\rhob}{\overline{\rho}} \rhob = \frac{N}{L-1}$.

The tracer is different from the bath particles in two ways. First, it hops to the right and to the left with different rates that we call p and q, respectively. Second, it can exchange its position with a bath particle two sites away under the condition that the site between the tracer and the bath particle is vacant. This process takes place with rate $\newcommand{\pp}{p'} \pp$ to the right and $\newcommand{\qp}{q'} \qp$ to the left. The condition that the intermediate site has to be empty mimics the fact that, in more realistic systems of e.g. particles moving in a narrow channel, overtakes are easier when particles have more space. The allowed transitions are summarized in figure 1. Such a system can be easily simulated using the Monte Carlo algorithm.

Zoom In Zoom Out Reset image size

At large times we expect the tracer to have a finite velocity and the bath to reach a stationary nonequilibrium state in the frame of the tracer. We define the occupation variables in the frame of the tracer by $\tau_l = 0$ or 1 for $l=1, \ldots, L-1$ and use $\langle \ldots \rangle$ for ensemble averages. The tracer occupies site l  =  0. The average velocity of the tracer is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\vtr}{{v_{\rm tr}}} \newcommand{\qp}{q'} \newcommand{\pp}{p'} \displaystyle \label{eq:vtrdef} \vtr &= p (1-\langle \tau_1 \rangle) -q (1- \langle \tau_{L-1} \rangle) \nonumber \\ &\quad+ 2 \pp \langle (1-\tau_1) \tau_2 \rangle - 2 \qp \langle (1-\tau_{L-1}) \tau_{L-2}\rangle. \nonumber \end{align} \tag{ 1 }$$

The first term accounts for hops of the tracer one step to the right: the transition is allowed if site l  =  1 is empty, contributing a factor $1-\tau_1$, and then occurs with rate p. The third term accounts for exchanges to the right: this transition is allowed if site l  =  1 is empty and site l  =  2 is occupied (factor $\tau_2(1-\tau_1)$). It occurs with rate $\newcommand{\pp}{p'} \pp$, and the tracer moves two steps to the right (factor 2). The second and fourth terms are hops and exchanges to the left, respectively. An ensemble average of the whole expression is taken. We also define the densities $\rho_l = \langle \tau_l \rangle$ for $l=1, \ldots, L-1$.

A configuration of the system is entirely specified by the $\{\tau_l\}_{l=1, \ldots, N-1}$, supplemented with the position of the tracer in the lab frame. In the case p  =  q and $\newcommand{\pp}{p'} \newcommand{\qp}{q'} \pp = \qp$ it is clear that the rate of every allowed transition between two states of the system is equal to the rate of the inverse transition. This implies that detailed balance is satisfied and that the stationary distribution is flat.

Interesting phenomena happen when a drive is applied to the system, namely $q \neq p$ or $\newcommand{\pp}{p'} \newcommand{\qp}{q'} \qp \neq \pp$. We now present numerical results for the tracer velocity and the density profile before showing how they can be understood analytically.

2.2. Numerical results

The system defined above can be easily simulated for any values of $p = r + \frac{\delta}{2}$, $q = r -\frac{\delta}{2}$, $\newcommand{\pp}{p'} \newcommand{\rp}{r'} \newcommand{\deltap}{\delta'} \pp = \rp + \frac{\deltap}{2}$ and $\newcommand{\qp}{q'} \newcommand{\rp}{r'} \newcommand{\deltap}{\delta'} \qp = \rp - \frac{\deltap}{2}$. In figures 2 and 3 we present numerical results of the tracer velocity for fixed values of r and $r'$ as a function of the respective biases, first $\delta \neq 0$ and then $\newcommand{\deltap}{\delta'} \deltap \neq 0$. In particular, for $\delta \neq 0$ and $\newcommand{\deltap}{\delta'} \deltap = 0$ (figure 2), the curves are monotonously increasing for small densities, but start to exhibit NDM and even ANM for larger densities. On the contrary, for $\delta = 0$ and $\newcommand{\deltap}{\delta'} \deltap \neq 0$ the curves are monotonously increasing (figure 3).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figures 4 and 5 we plot the corresponding density profiles for different values of the average density. The density is found to be flat in the bulk of the system, i.e. far from the tracer, with a meniscus appearing on one side of the tracer. It appears that the change of sign in the velocity for $\delta \neq 0$ as the density is increased is accompanied by a qualitative change in the density profile, where a meniscus appears at the front of the tracer for low densities and at its back for high densities (see figure 4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the next subsection we study the system around its equilibrium state using a mean-field approximation, and compute $\newcommand{\vtr}{{v_{\rm tr}}} \vtr$ in the linear response regime.

2.3. Tracer velocity

We start by writing mean-field equations for the densities $\{\rho_l\}_{l=1, \ldots, L-1}$,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\qp}{q'} \newcommand{\pp}{p'} \newcommand{\dd}{{\rm d}} \displaystyle \label{eq:mfdensity} &\frac{\dd \rho_1}{\dd t} = \rho_2 -\rho_1 + p (1-\rho_1) \rho_2 - q (1-\rho_{L-1}) \rho_1 +\pp (1-\rho_1) \rho_2 \rho_3 - \qp (1-\rho_{L-1}) \rho_{L-2} \rho_1, \nonumber \\ &\frac{\dd \rho_2}{\dd t} = \rho_3 -2 \rho_2 + \rho_1 + p (1-\rho_1) (\rho_3 - \rho_2) + q (1-\rho_{L-1}) (\rho_1-\rho_2) \nonumber \\ &\quad +\pp (1-\rho_1) \rho_2 (\rho_4 -1) + \qp (1-\rho_{L-1}) \rho_{L-2} (1-\rho_2), \nonumber \\ & \frac{\dd \rho_l}{\dd t} = \rho_{l+1} -2 \rho_l + \rho_{l-1} + p (1-\rho_1) (\rho_{l+1} - \rho_l) + q (1-\rho_{L-1}) (\rho_{l-1}-\rho_l) \nonumber \\ & \quad+\pp (1-\rho_1) \rho_2 (\rho_{l+2} -\rho_l) + \qp (1-\rho_{L-1}) \rho_{L-2} (\rho_{l-2}-\rho_l), \qquad \qquad l=3,\ldots, L-3, \nonumber \\ & \frac{\dd \rho_{L-2}}{\dd t} = \rho_{L-1} -2 \rho_{L-2} + \rho_{L-3} + p (1-\rho_1) (\rho_{L-1} - \rho_{L-2}) + q (1-\rho_{L-1}) (\rho_{L-3}-\rho_{L-2}), \nonumber \\ &\quad + \pp (1-\rho_1) \rho_2 (1 - \rho_{L-2}) + \qp (1-\rho_{L-1}) \rho_{L-2} (\rho_{L-4}-1), \nonumber \\ &\frac{\dd \rho_{L-1}}{\dd t} = -\rho_{L-1} + \rho_{L-2} - p (1-\rho_1) \rho_{L-1} + q (1-\rho_{L-1}) \rho_{L-2} - \pp (1-\rho_1) \rho_2 \rho_{L-1} + \qp (1-\rho_{L-1}) \rho_{L-2} \rho_{L-3}, \nonumber \end{align} \tag{ 2 }$$

where we factorized correlations $\langle \tau_{l_1} \ldots \tau_{l_k}\rangle = \rho_{l_1} \ldots \rho_{l_k}$ for distinct positions. Since correlations are factorized in equilibrium, it is reasonable to expect that this approximation will give good results at least close to equilibrium. A similar technique has also been proven accurate in closely related systems, see e.g. [19, 25–35]. The tracer velocity (1) becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\vtr}{{v_{\rm tr}}} \newcommand{\qp}{q'} \newcommand{\pp}{p'} \displaystyle \label{eq:vtrmf} \vtr = p (1-\rho_1) -q (1- \rho_{L-1}) + 2 \pp (1-\rho_1) \rho_2 - 2 \qp (1-\rho_{L-1}) \rho_{L-2}. \nonumber \end{align} \tag{ 3 }$$

We consider equations (2) in the stationary state $\newcommand{\dd}{{\rm d}} \dd \rho_l/\dd t = 0$. Because of particle number conservation, they give only L  −  2 independent conditions. The missing $(L-1){\hspace{0pt}}^{\rm th}$ condition is obtained by fixing the number of particles,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:norm} \sum_{l=1}^{L-1} \rho_l = N. \nonumber \end{align} \tag{ 4 }$$

When there is no bias, $\delta = \delta' = 0$, it is easy to confirm that a flat density profile $\newcommand{\rhob}{\overline{\rho}} \rho_l = \rhob$ solves equations (2)-(4) so that the tracer velocity (3) vanishes.

For simplicity, we now consider the case where the biases δ and $\delta'$ are small and of the same order. We expand the density,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rhob}{\overline{\rho}} \newcommand{\cO}{\mathcal{O}} \displaystyle \label{eq:defsigma} \rho_l = \rhob + \delta \sigma_l + \delta' \sigma_l' + \cO(\delta^2), \nonumber \end{align} \tag{ 5 }$$

and study the solution to linear order in δ and $\newcommand{\deltap}{\delta'} \deltap$. While the tracer velocity is rather well-described to this order, this is not the case for the density profile. In order to obtain the density profile, one has to study the equations to second order in δ and $\newcommand{\deltap}{\delta'} \deltap$. This will be done in the next subsection.

We start by solving the equations for the bulk sites $l=3, \ldots, L-3$. The terms linear in δ and $\delta'$ give equations for $\sigma_l$ and $\sigma_l'$, respectively. Both bulk equations turn out to be the same,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rhob}{\overline{\rho}} \displaystyle \label{eq:bulk} (1+r (1-\rhob)) (\sigma_{l+1}-2 \sigma_l+\sigma_{l-1}) + r' \rhob (1-\rhob) (\sigma_{l+2}-2 \sigma_l+\sigma_{l-2}) &= 0,\nonumber \end{align} \tag{ 6 }$$

for $l=3, \ldots, L-3$, and the very same equation for the $\{\sigma_l'\}_{l=1, \dots, L-1}$. Note that in the continuum limit this equation would reduce to a Laplace equation $\partial_l^2 \sigma = 0$, giving a linear density profile. The solution of the discrete equation is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:solbulk} \sigma_l &= \alpha + \beta l +\gamma_+ X^l + \gamma_- X^{L-l}, \nonumber \\ \sigma_l' &= \alpha' + \beta' l +\gamma_+' X^l + \gamma_-' X^{L-l}, \nonumber \end{align} \tag{ 7 }$$

for $l=1, \ldots, L-1$, where

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rhob}{\overline{\rho}} \newcommand{\rp}{r'} \displaystyle \label{eq:Xdef} X &= - \left(1+ \frac{1+r (1-\rhob)}{2 \rp \rhob (1-\rhob)}\right) + \sqrt{\left(1+ \frac{1+r (1-\rhob)}{2 \rp \rhob (1-\rhob)}\right)^2-1} \nonumber \end{align} \tag{ 8 }$$

and X⁻¹ are the roots of

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rhob}{\overline{\rho}} \newcommand{\rp}{r'} \displaystyle \label{eq:Xeq} (1+r(1-\rhob)) X +\rp \rhob (1-\rhob) (1+X)^2 = 0. \nonumber \end{align} \tag{ 9 }$$

Note that $\vert X\vert < 1$.

We now have to satisfy the boundary and normalization conditions. In the large N and L limit with $\newcommand{\rhob}{\overline{\rho}} \frac{N}{L-1} = \rhob$, the normalization (4) gives $\beta = - \frac{2 \alpha}{L}$ and $\beta' = - \frac{2 \alpha'}{L}$. Let us now take the sum of the equations for $\rho_1$ and $\rho_{L-1}$. Sorting the δ and $\delta'$ terms, we again get the same equation for $\sigma_l$ and $\sigma_l'$,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rhob}{\overline{\rho}} \newcommand{\rp}{r'} \displaystyle \label{eq:s1Lm1} (1+r(1-\rhob)) (\sigma_2 - \sigma_1 + \sigma_{L-2} - \sigma_{L-1}) + \rp \rhob (1-\rhob) (\sigma_3 - \sigma_1 + \sigma_{L-3} - \sigma_{L-1}) &= 0. \nonumber \end{align} \tag{ 10 }$$

Taking the sum of the equations for $\rho_2$ and $\rho_{L-2}$ leads to the same result. For large L, we have that $\sigma_l \simeq \alpha + \gamma_+ X^l$ for $\newcommand{\cO}{\mathcal{O}} l=\cO(1)$, and $\sigma_{L-l} \simeq -\alpha + \gamma_- X^l$ for $\newcommand{\cO}{\mathcal{O}} L-l=\cO(1)$. Inserting these forms in equation (10), we get $\gamma_- = - \gamma_+$ and, similarly, $\gamma_-' = - \gamma_+'$. The linear perturbations to the densities therefore have the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:solbulks} \sigma_l &= \alpha \left(1- 2 \frac{l}{L} \right) +\gamma_+ \left(X^l - X^{L-l} \right), \nonumber \\ \sigma_l' &= \alpha' \left(1- 2 \frac{l}{L} \right) +\gamma_+' \left(X^l - X^{L-l} \right). \nonumber \end{align} \tag{ 11 }$$

The equations for, say, $\rho_1$ and $\rho_2$ give two systems of two linear equations for α and $\gamma_+$, and for $\alpha'$ and $\gamma_+'$. The solutions to these equations are obtained using Mathematica and are given in appendix A.

The density profile obtained in this analysis is linear in the bulk with exponential layers on both sides of the tracer. This is different from the numerical observation of a flat profile in the bulk and an exponential layer only on one side of the tracer. This discrepancy is a result of the fact that the analysis has been carried out to linear order in δ and $\newcommand{\deltap}{\delta'} \deltap$. This will be corrected in section 2.4.

The velocity of the tracer can be obtained from the mean-field expression (3) and the solutions (11), where the coefficients are given by (A.1). We separate it into two contributions, namely the one coming from the hops of the tracer towards an empty site (terms proportional to p and q in equation (3)) and the one coming from exchanges (terms proportional to $\newcommand{\pp}{p'} \pp$ and $\newcommand{\qp}{q'} \qp$ in the same equation). We can write $\newcommand{\vtr}{{v_{\rm tr}}} \newcommand{\vtrh}{v_{\rm tr,H}} \newcommand{\vtre}{v_{\rm tr,E}} \vtr = \vtrh + \vtre$, where the subscripts ${\rm H}$ and ${\rm E}$ indicate the contributions coming from hops and exchanges, respectively. Each of these pieces has a term proportional to its corresponding bias,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\vtr}{{v_{\rm tr}}} \newcommand{\deltap}{\delta'} \displaystyle \vtr_{,A} = \mu_{A,{\rm H}} \delta + \mu_{A,{\rm E}} \deltap, \nonumber \end{align} \tag{ 12 }$$

which gives four coefficients $\mu_{A, B}$ with $A, B \in \{{\rm H}, {\rm E}\}$. Explicitly, they are

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\muee}{\mu_{{\rm E},{\rm E}}} \newcommand{\rhob}{\overline{\rho}} \newcommand{\rp}{r'} \displaystyle \label{eq:musol} \mu_{{\rm H},{\rm H}} &= (1-\rhob) - r (\sigma_1 - \sigma_{L-1}) = (1-\rhob) - 2 r (\alpha + \gamma_+ X) = \frac{\rp (1-2 \rhob)^2}{2 r \rhob^2} \muee, \nonumber \\ \mu_{{\rm H},{\rm E}} &= - r (\sigma_1' - \sigma_{L-1}') = - 2 r (\alpha' + \gamma_+' X) = \frac{1-2 \rhob}{2 \rhob} \muee, \nonumber \\ \mu_{{\rm E},{\rm H}} &= 2 \rp \left((1-\rhob) (\sigma_2 - \sigma_{L-2}) + \rhob (\sigma_{L-1} - \sigma_1)\right) = 4 \rp \left((1-\rhob) (\alpha + \gamma_+ X^2) - \rhob (\alpha + \gamma_+ X) \right) = \frac{\rp (1-2 \rhob)}{r \rhob} \muee, \nonumber \end{align} \tag{ 13 }$$

with the last coefficient

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\muee}{\mu_{{\rm E},{\rm E}}} \newcommand{\rhob}{\overline{\rho}} \newcommand{\rp}{r'} \displaystyle \label{eq:muee} \muee &= 2 \rhob (1-\rhob) + 2 \rp (1-\rhob) (\sigma_2' - \sigma_{L-2}') + 2 \rp \rhob (\sigma_{L-1}'-\sigma_1') = 2 \rhob (1-\rhob) + 4 \rp (1-\rhob) (\alpha' + \gamma_+' X^2) - 4 \rp \rhob (\alpha' + \gamma_+' X) \nonumber \\ &= \frac{2 \rhob^2 (1-\rhob) r \sqrt{1+r (1-\rhob)}}{(2 \rhob - 1) (r+\rp (2 \rhob - 1)) \sqrt{1+r (1-\rhob)} + r (1-\rhob) \sqrt{1+r(1-\rhob) + 4 \rp \rhob (1-\rhob)}}. \nonumber \end{align} \tag{ 14 }$$

In appendix C we show that $\newcommand{\muee}{\mu_{{\rm E},{\rm E}}} \muee > 0$ for all $\newcommand{\rhob}{\overline{\rho}} 0 < \rhob < 1$, r  >  0, $\newcommand{\rp}{r'} \rp > 0$. For clarity, let us group the linear response coefficients (13) and (14) into a linear response matrix,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\muee}{\mu_{{\rm E},{\rm E}}} \newcommand{\rhob}{\overline{\rho}} \newcommand{\vtre}{v_{\rm tr,E}} \newcommand{\vtrh}{v_{\rm tr,H}} \newcommand{\vtr}{{v_{\rm tr}}} \newcommand{\deltap}{\delta'} \newcommand{\rp}{r'} \displaystyle \label{eq:linresmat} \left(\begin{array}{@{}c@{}} \vtrh \nonumber \\ \vtre \end{array} \right) = \frac{\rp \muee}{2 \rhob^2} \left(\begin{array}{@{}cc@{}} (1-2\rhob)^2 & 2 \rhob (1-2\rhob) \nonumber \\ 2 \rhob (1-2\rhob) & 4 \rhob^2 \end{array} \right) \left(\begin{array}{@{}c@{}} \delta/r \nonumber \\ \deltap/(2\rp) \end{array} \right), \nonumber \end{align} \tag{ 15 }$$

where the entries of the column vector on the RHS are the thermodynamically conjugate forces. In this basis the response matrix is symmetric, as expected from the Onsager relations. In (15) we see that the diagonal coefficients of the response matrix are always positive, consistent with fluctuation-dissipation relations. Conversely, the off-diagonal coefficients need not be positive and indeed they change sign at $\newcommand{\rhob}{\overline{\rho}} \rhob = \frac{1}{2}$, which allows for ANM. Thus ANM found in this model in the linear response regime is a direct result of the fact that the dynamics involves two driving mechanisms, namely hopping (p/q) and exchange ($\newcommand{\pp}{p'} \newcommand{\qp}{q'} \pp/\qp$). Note that the columns of the linear response matrix are proportional, which shows that the response of the tracer to the two driving fields is the same. This indicates that exchange and hopping are completely coupled in the sense of [79]. The total velocity becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\muee}{\mu_{{\rm E},{\rm E}}} \newcommand{\rhob}{\overline{\rho}} \newcommand{\vtr}{{v_{\rm tr}}} \newcommand{\deltap}{\delta'} \newcommand{\rp}{r'} \displaystyle \label{eq:linresvtr} \vtr = \frac{\muee}{2 r \rhob^2} \left(\rp (1-2\rhob) \delta + r \rhob \deltap \right). \nonumber \end{align} \tag{ 16 }$$

For fixed $\newcommand{\deltap}{\delta'} \delta, \deltap > 0$, the velocity starts out positive for small $\newcommand{\rhob}{\overline{\rho}} \rhob$, and changes sign for $\newcommand{\rp}{r'} \newcommand{\deltap}{\delta'} \newcommand{\rhob}{\overline{\rho}} \rhob = \left(2 - \frac{r \deltap}{\rp \delta} \right){\hspace{0pt}}^{-1}$, which is smaller than 1 for $\newcommand{\rp}{r'} \newcommand{\deltap}{\delta'} r \deltap < \rp \delta$. The prediction (16) is compared to the results of Monte-Carlo simulations in figures 2 and 3, and the agreement is very good.

Similarly, one can predict the current $\newcommand{\Jb}{\mathcal{J}_{\rm B}} \Jb$ of bath particles in the linear regime. It can be expressed as a function of the densities in the neighborhood of the tracer (see appendix B),

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Jb}{\mathcal{J}_{\rm B}} \newcommand{\qp}{q'} \newcommand{\pp}{p'} \displaystyle \label{eq:Jb} \Jb = \frac{\rho_1 - \rho_{L-1} + 2 \qp (1-\rho_{L-1}) \rho_{L-2} - 2 \pp (1-\rho_1) \rho_2}{L}. \nonumber \end{align} \tag{ 17 }$$

Using the computed values of $\rho_1$, $\rho_2$, $\rho_{L-2}$ and $\rho_{L-1}$,, one obtains analytical predictions for $\newcommand{\Jb}{\mathcal{J}_{\rm B}} \Jb$. They are compared to Monte Carlo simulations in figures 6 and 7. The slope at the origin is predicted accurately.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Note that the analysis presented in this section yields good agreement for $\newcommand{\vtr}{{v_{\rm tr}}} \vtr$ and $\newcommand{\Jb}{\mathcal{J}_{\rm B}} \Jb$, since they are determined by the average density of the sites close to the tracer (equations (3) and ((17)). The total current can be obtained as $\newcommand{\vtr}{{v_{\rm tr}}} \newcommand{\Jb}{\mathcal{J}_{\rm B}} \frac{\vtr}{L} + \Jb$ and is therefore predicted accurately by the linear analysis as well.

While these densities are well described by equations (7) and (8), the overall density profile is not. Let us now focus on correcting the discrepancy in the density profile by extending the analysis to second order in δ and $\newcommand{\deltap}{\delta'} \deltap$.

2.4. Density profile

In order to explain the form of the density profile for small biases, one has to keep higher-order terms in the expansion. We consider a system with small biases δ and $\newcommand{\deltap}{\delta'} \deltap$ and write $\newcommand{\rhob}{\overline{\rho}} \phi_l = \rho_l - \rhob$. Expanding the bulk equation to second order in $\phi_l$ gives

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rhob}{\overline{\rho}} \newcommand{\deltap}{\delta'} \newcommand{\rp}{r'} \displaystyle \label{eq:phi} 0 &= (1+r(1-\rhob)) (\phi_{l+1} - 2 \phi_l + \phi_{l-1}) + \rp \rhob (1-\rhob) (\phi_{l+2} - 2 \phi_l + \phi_{l-2}) \nonumber \\ &\quad+ \left((1-\rhob) \frac{\delta}{2} - r \phi_1 \right) (\phi_{l+1} - \phi_l) + \left(-(1-\rhob) \frac{\delta}{2} - r \phi_{L-1} \right) (\phi_{l-1} - \phi_l) \nonumber \\ &\quad + \left((1-\rhob) \rhob \frac{\deltap}{2} - \rp \rhob \phi_{1} + \rp (1-\rhob) \phi_{2} \right) (\phi_{l+2} - \phi_l) + \left(-(1-\rhob) \rhob \frac{\deltap}{2} - \rp \rhob \phi_{L-1} + \rp (1-\rhob) \phi_{L-2} \right) (\phi_{l-2} - \phi_l) \nonumber \\ &\simeq \left[1+r(1-\rhob) + 4 \rp \rhob (1-\rhob)\right] \partial_l^2 \phi + \left[ (1-\rhob) \delta + 2 \rhob (1-\rhob) \deltap - (r+2\rp \rhob) (\phi_1 - \phi_{L-1}) + 2 \rp (1-\rhob) (\phi_2 - \phi_{L-2}) \right] \partial_l \phi, \nonumber \end{align} \tag{ 18 }$$

where we approximated discrete differences by derivatives. The coefficient of the first order derivative is exactly the part of $\newcommand{\vtr}{{v_{\rm tr}}} \vtr$ linear in ϕ, and it depends on values of ϕ close to the tracer. Let us replace them by their values from the preceding section,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\deltap}{\delta'} \displaystyle \label{eq:approxphi} \phi_n - \phi_{L-n} \simeq 2 \delta (\alpha + \gamma_+ X^n) + 2 \deltap (\alpha' + \gamma_+' X^n), \nonumber \end{align} \tag{ 19 }$$

for n  =  1,2, so that the aforementioned coefficient exactly becomes the $\newcommand{\vtr}{{v_{\rm tr}}} \vtr$ as obtained in equation (16). The solution of (18) is exponential,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ee}{{\rm e}} \displaystyle \label{eq:solphi0} \phi_l = C_1 \ee^{-\frac{l}{\xi}} + C_2, \nonumber \end{align} \tag{ 20 }$$

where C₁ and C₂ are integration constants, and the decay length is

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rhob}{\overline{\rho}} \newcommand{\vtr}{{v_{\rm tr}}} \newcommand{\rp}{r'} \displaystyle \label{eq:decaylngth} \xi = \frac{1+r(1-\rhob) + 4 \rp \rhob (1-\rhob)}{\vtr}. \nonumber \end{align} \tag{ 21 }$$

We use the definition (21), where ξ can be either positive or negative, to make the presentation simpler.

The decay length is of order $\delta^{-1}$ or $\newcommand{\deltap}{\delta'} \deltap^{-1}$, which explains why we found linear profiles when we neglected $\newcommand{\cO}{\mathcal{O}} \cO(\delta^2)$ terms. One can check that expression (21) diverges at $\newcommand{\rhob}{\overline{\rho}} \rhob=\frac{1}{2}$ for $\newcommand{\deltap}{\delta'} \deltap = 0$, and at $\newcommand{\rhob}{\overline{\rho}} \rhob=1$ for any δ, $\newcommand{\deltap}{\delta'} \deltap$. This is consistent with the respective linear profiles obtained in figure 4 and the high-density calculation of section 2.5. In the low-density limit $\newcommand{\rhob}{\overline{\rho}} \rhob \to 0$ one simply gets $\xi \sim \frac{1+r}{\delta}$, which is the diffusion coefficient of the free tracer divided by the bias.

The constants C₁ and C₂ are linked by mass conservation, $\sum_{l=1}^{L-1} \phi_l = 0$, giving

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ee}{{\rm e}} \displaystyle \label{eq:consmC1C2} C_2 = - \frac{1}{L-1} \frac{\ee^{-\xi^{-1}} - \ee^{-\xi^{-1} L}}{1-\ee^{-\xi^{-1}}} C_1. \nonumber \end{align} \tag{ 22 }$$

As $\newcommand{\ee}{{\rm e}} \newcommand{\e}{{\rm e}} C_1 = \left(\ee^{-\xi^{-1}} - \ee^{-\xi^{-1}(L-1)} \right){\hspace{0pt}}^{-1} \left(\phi_1 - \phi_{L-1}\right)$, we may employ approximation (19) again in order to obtain the value of C₁. We end up with

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\deltap}{\delta'} \newcommand{\ee}{{\rm e}} \displaystyle \label{eq:solphi} \phi_l &= \frac{2 \left[ (\alpha+\gamma_+ X) \delta + (\alpha'+\gamma_+' X) \deltap \right]}{1 - \ee^{-\xi^{-1}(L-2)}} \times \left[ \ee^{-\xi^{-1} (l-1)} - \frac{1}{L-1} \frac{1 - \ee^{-\xi^{-1}(L-1)}}{1-\ee^{-\xi^{-1}}} \right]. \nonumber \end{align} \tag{ 23 }$$

The form (23) is shown to be in good agreement with simulations in figures 4 and 5.

2.5. High density regime

Here we go beyond linear response for high densities $\newcommand{\rhob}{\overline{\rho}} \rhob \simeq 1$. In that case the holes are very sparse and can be considered independent, so that we can start by studying a system with L  −  2 bath particles and only one hole. Let $l = 1, \ldots, L-1$ denote the position of the hole with respect to the tracer. Examination shows that the probability distribution of the position of the hole, P_l(t), obeys

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\qp}{q'} \newcommand{\pp}{p'} \newcommand{\dd}{{\rm d}} \displaystyle \label{eq:meqhole} \frac{\dd P_l}{\dd t} =&\, P_{l+1} - 2 P_l + P_{l-1} + \delta_{l,1} [(1-p-\pp)P_1 + (q+\qp) P_{L-1}] \nonumber \\ & + \delta_{l,L-1} [ (1-q-\qp)P_{L-1} + (\,p+\pp) P_1 ], \nonumber \end{align} \tag{ 24 }$$

with the convention $P_0 = P_L = 0$, and the normalization $\sum_{l=1}^{L-1} P_l = 1$. When the hole is far from the tracer, $l \neq 1, L-1$, it simply diffuses as seen on the first term on the right hand side of (24). The two other terms in this equation correspond to hopping and exchange processes which take place when the hole is next to the tracer.

In the stationary state and large L limit, these equations give

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\qp}{q'} \newcommand{\pp}{p'} \displaystyle \label{eq:solPY} P_l = \frac{2}{L(\,p+\pp+q+\qp)} \left((\,p+\pp-q-\qp) \frac{l}{L} +q+\qp \right). \nonumber \end{align} \tag{ 25 }$$

This gives, for one hole,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\vtre}{v_{\rm tr,E}} \newcommand{\vtrh}{v_{\rm tr,H}} \newcommand{\vtr}{{v_{\rm tr}}} \newcommand{\deltap}{\delta'} \newcommand{\rp}{r'} \newcommand{\qp}{q'} \newcommand{\pp}{p'} \displaystyle \label{eq:linresmathighole} \left(\begin{array}{@{}c@{}} \vtrh^{(1)} \nonumber \\ \vtre^{(1)} \end{array} \right) = \left(\begin{array}{@{}c@{}} p P_1 - q P_{L-1} \nonumber \\ 2 \pp P_1 - 2 \qp P_{L-1} \end{array} \right) = \frac{r \deltap - \rp \delta}{L (r + \rp)} \left(\begin{array}{@{}c@{}} -1 \nonumber \\ 2 \end{array} \right), \nonumber \end{align} \tag{ 26 }$$

to all orders in δ and $\newcommand{\deltap}{\delta'} \deltap$. For a system with not too large a number of holes $\newcommand{\rhob}{\overline{\rho}} L (1-\rhob)$, we can simply add the effect of each hole. This gives

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rhob}{\overline{\rho}} \newcommand{\vtre}{v_{\rm tr,E}} \newcommand{\vtrh}{v_{\rm tr,H}} \newcommand{\vtr}{{v_{\rm tr}}} \newcommand{\deltap}{\delta'} \newcommand{\rp}{r'} \displaystyle \label{eq:vtrhighrho} \vtr = - \vtrh = \frac{\vtre}{2} = (1-\rhob) \frac{r \deltap - \rp \delta}{r + \rp}, \nonumber \end{align} \tag{ 27 }$$

and the density profile is linear,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rhob}{\overline{\rho}} \newcommand{\deltap}{\delta'} \newcommand{\rp}{r'} \displaystyle \label{eq:denhighrho} \rho_l = \rhob- (1-\rhob) \frac{\delta+\deltap}{r +\rp} \left(\frac{l}{L} - \frac{1}{2}\right). \nonumber \end{align} \tag{ 28 }$$

The results (27) and (28) are expected to be exact in the high-density limit. They also match the high-density limits of the mean-field predictions for the velocity (16) and the density profile (23). Indeed, the agreement between the predicted profile (28) and the numerical results can be checked to be very good for large densities.

More importantly, considering a system with one hole helps to shed light on the way ANM occurs, see figure 8 and the explanation in the caption. The sequence of transitions shown in figure 8 contains a step where the tracer hops to the right and is therefore favored by an increase of p. The net result of this sequence is an overall displacement of the system to the left. Symmetrically, an increase of q favors a sequence of transitions that results in a net displacement of the tracer to the right. Therefore, when p q'  >  q p' the tracer moves preferentially to the left and ANM is observed.

Zoom In Zoom Out Reset image size

Motivated by the lattice model of driven tracer presented in the preceding section, we move a small step in the direction of a more realistic system and consider a gas of hard disks in a planar narrow channel. Even though the fluctuation-dissipation relation forbids ANM in linear response, the gas could in principle exhibit ANM at large driving field. The channel is oriented parallel to the x axis with the origin of the coordinate system in its center. It is periodic in x with periodicity L_x. The channel width $\newcommand{\e}{{\rm e}} L_y' = 2 \sigma + \epsilon$ is chosen to allow passing and overtaking of the particles with diameter σ. Thus, $\newcommand{\e}{{\rm e}} \epsilon > 0$ is assumed. In this channel there are N  −  1 neutral particles with mass m_k, position ${\bf r}_k = (x_k, y_k)$, and velocity ${\bf v}_k = (v_{k, x}, v_{k, y})$, $k \in \{1, 2, \cdots, N-1\}$. In addition, there is a tracer particle with index k  =  0, which is driven by a homogeneous force ${\bf F} = (F, 0)$ parallel to the channel axis. The dynamics of all disks is assumed to follow an underdamped Langevin equation,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{\rm d}{\bf r}_k}{{\rm d}t} = {\bf v}_k \label{eom1} \nonumber \end{align} \tag{ 29 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle m_k \frac{{\rm d}{\bf v}_k}{{\rm d}t} &= {\bf F} \delta_{k,0} -\gamma {\bf v}_k + \sqrt{2 \gamma k_{\rm B} T} \xi_0, \label{eom2} + \,\{{\rm coll.}\}\nonumber \\ \nonumber \end{align} \tag{ 30 }$$

for $k \in \{0, 1, \cdots, N-1 \}$. As usual, ξ is a delta-correlated Gaussian white noise, $\delta_{k, 0}$ is unity for k  =  0 and zero otherwise, k_B is the Boltzmann constant, and T is the temperature of the bath.

The stochastic motion equations are solved to first order in the time step $\Delta t$ according to the updating formulas of Gillespie [80, 81]. Reduced units are used throughout, for which the mass of the tracer, m₀, the particle diameter σ and the energy $k_{\rm B} T$ are unity. In these units, the Langevin friction parameter γ, which determines the noise strength, is set to 2 and the time step $\Delta t$ to 10⁻³. In all simulations below, the following parameters are used: channel length L_x  =  300, total number of particles N  =  200, and tracer mass m₀  =  1.0. The masses $\newcommand{\e}{{\rm e}} m_k \equiv m$ of the neutral particles $k = 1, \cdots, 199$ and the driving force F are indicated where needed.

There are two kinds of collisions, namely collisions between particles and collisions of a particle with a long hard boundary of the channel. Particle–particle collisions are strictly elastic. The long channel boundaries, however, are thermal van Beijeren walls [82], which re-emit colliding particles in equilibrium with the boundary temperature T_b. The latter is taken to agree with the bath temperature, $T_{\rm b} = T$. Since in such a narrow channel the long thermal walls are of significant importance for the non-equilibrium transport, a short description of the van Beijeren walls used here is given in appendix D.

The dependence of the mean tracer speed $\langle v_{\rm tr} \rangle$ on the drive F is shown in figure 9 for a channel of width L_y'  =  2.6. The four curves are for different masses of the neutral particles as indicated by the labels. For low F the velocity increases almost linearly with F (Ohm's law). For larger fields, however, the curves bend over and reach a regime with a negative slope indicating negative differential mobility. This behavior is most prominent for 2.3  <  L_y'  <  2.8 and deteriorates quickly for broader channels. This is demonstrated by a comparison of the top and bottom panels of figures 10 for L_y'  =  2.8 and 2.9, respectively. For L_y'  =  2.8, one observes very strong negative mobility. But an increase of the channel width by a moderate amount to 2.9 gives a totally different picture. The slopes of the characteristic curves always remain positive.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The origin of this strange and unexpected behavior may be understood by inspecting a snapshot of the particle configuration in the neighborhood of the tracer. See figure 11. The (green) neutral particles accumulate in front of the tracer, which diffuses in the direction of F. The particles in front condense to crystal-like structures, which do not dissolve easily and which become longer and more stable with increasing force F. In the limit of very large F all neutral particles are included in this cluster, and the average speed of the tracer becomes minimal but still remains positive. The transition from the Ohmic regime for small F to this nearly blocked transport for large F is characterized by the negative differential mobility mentioned above.

Zoom In Zoom Out Reset image size

For a channel width larger than 2.8, the tracer and the neutral particles have a tendency to stay closer to one of the long channel boundaries than in the channel center. If this happens, the tracer has no difficulty to pass other neutral particles close to the other boundary which, consequently, enhances its speed and disallows NDM. Also the very narrow channels geometrically do not favor the nucleation of crystal-like clusters in front of the tracer, which could block its advance. Consequently, NDM is also not observed in this limit.

At this stage a few comments are in order:

• (i)  
  We have stressed above the importance of the boundary collisions for the nonlinear transport. We have verified that NDM also takes place—albeit for stronger forces F—if the thermal walls are replaced by simple elastically-reflecting boundaries. Thus, this effect appears to be very robust. For small F, however, the nature of this boundary becomes unimportant [40].
• (ii)  
  In our quest for absolute negative mobility, also modifications of the equations of motion (30) were investigated. In one attempt, an attractive harmonic force acting on the y coordinate of the tracer was added. This increases its position probability in the center of the channel. In another, also a short-range force between the tracer and bath particles was added, which was designed to enhance the passing probability of these particles when close by. None of these attempts have indicated ANM. Unfortunately, the parameter space increases enormously by such attempts and has not yet been explored to any significant extent.

In this work we exhibited two instances of negative response. The first one features ANM around an equilibrium state and occurs in a simple lattice model. We are able to predict this phenomenon analytically for small biases and shed light on the underlying mechanism in very dense systems. Keeping the same spirit, we then probed a more realistic system composed of hard disks in a channel. No ANM was found in that case, however we showed that a 'crystallization' of the neutral Brownian particles at large biases can lead to NDM.

It is natural to ask which kind of modification would allow ANM in the hard-disks model in a way similar to the lattice model. Since the drive is homogeneous, the model as studied in this work cannot be expected to exhibit ANM in linear response around equilibrium. However, in principle it could exhibit ANM at large drive, something which we have not observed. Extending the model, e.g. by introducing a y-dependent force or an extra potential, would introduce a second driving field that could reproduce the behaviour observed in the lattice model. Emulating the exchange move in the hard-disks system is however not such an easy task. Other variants also come to mind, like changing the masses or the radii of the particles, not to mention their shapes. The parameter space being extremely large, we leave this question open for future research.

We thank M Aizenmann, H van Beijeren, B Derrida, Y Kafri, A Kundu, S N Majumdar and A Miron for interesting discussions about this problem. One of us (HAP) wants to acknowledge the hospitality and support of the Weizmann Institute of Science, where parts of the work reported here have been performed. The support of the Israel Science Foundation (ISF) is gratefully acknowledged. The molecular dynamics simulations were carried out on the Vienna Scientific Cluster (VSC). We are grateful for the generous allocation of computer resources.

Here we give the coefficients that appear in the density profile calculation of section 2.3 and a few limits. Their expressions are

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rhob}{\overline{\rho}} \newcommand{\rp}{r'} \displaystyle \label{eq:ag} \alpha &= \frac{1}{\Xi} \frac{(1-\rhob) \rhob}{2 (1-X)} \left(1+r(1-\rhob)+4 \rp \rhob (1-\rhob) + 2 \rp (1-\rhob) X (2 \rhob - 1)\right), \nonumber \\ \gamma_+ &= \frac{1}{\Xi} \frac{\rp (1-\rhob)^2 \rhob (1-2 \rhob)}{X (1-X)}, \nonumber \\ \alpha' &= \frac{1}{\Xi} \frac{(1-\rhob) \rhob }{2 (1-X)} \left(-r (1-\rhob) + (2 \rhob-1)(1+2\rp (1+X)(1-\rhob) \rhob)\right), \nonumber \\ \gamma_+' &= \frac{1}{\Xi} \frac{r (1-\rhob)^2 \rhob^2}{X (1-X)}, \nonumber \end{align} \tag{ A.1 }$$

with the shorthand

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rhob}{\overline{\rho}} \newcommand{\rp}{r'} \displaystyle \label{eq:defxi} \Xi = \rp^2 (1+X) (1-2\rhob)^2 (1-\rhob) \rhob +r \rhob +r^ 2 \rhob (1-\rhob) +\rp (1-2\rhob)^2 +r \rp (1-\rhob) (1+\rhob (-2+(3+X)\rhob)), \nonumber \end{align} \tag{ A.2 }$$

and X is defined in (8).

For low densities $\newcommand{\rhob}{\overline{\rho}} \rhob \to 0$ we get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rhob}{\overline{\rho}} \newcommand{\rp}{r'} \displaystyle \label{eq:lowrhocoeffs} X &\sim - \frac{\rp}{1+r} \rhob, \nonumber \\ \Xi &\sim \rp (1+r), \nonumber \\ \alpha \sim - \alpha' &\sim \frac{\rhob}{2 \rp}, \nonumber \\ \gamma_+ \sim \frac{1}{r \rhob} \gamma_+' &\sim - \frac{1}{\rp}. \nonumber \end{align} \tag{ A.3 }$$

In the high-density limit $\newcommand{\rhob}{\overline{\rho}} \rhob \to 1$ we have

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rhob}{\overline{\rho}} \newcommand{\rp}{r'} \displaystyle \label{eq:highrhocoeffs} X &\sim - \rp (1-\rhob), \nonumber \\ \Xi &\sim r + \rp, \nonumber \\ \alpha \sim \alpha' &\sim \frac{1-\rhob}{2 (r+\rp)}, \nonumber \\ \gamma_+ \sim - \frac{\rp}{r} \gamma_+' &\sim \frac{1-\rhob}{r+\rp}. \nonumber \end{align} \tag{ A.4 }$$

The values of α and $\alpha'$ are consistent with the high-density calculation (28) and the terms involving γ and $\gamma'$ are negligible in the high-density limit.

One can also simplify

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rhob}{\overline{\rho}} \newcommand{\rp}{r'} \displaystyle \label{eq:signC1} \alpha'+\gamma_+' X &= (2 \rhob - 1) \frac{\rhob (1-\rhob)}{2 \Xi (1-X)} \nonumber \\ &\quad \times (r(1-\rhob) + 1 + 2 \rp (1+X) (1-\rhob) \rhob), \nonumber \end{align} \tag{ A.5 }$$

which shows that the sign of the exponential layer changes when $\newcommand{\rhob}{\overline{\rho}} \rhob = \frac{1}{2}$ for $\delta = 0$.

Here we prove expression (17) for the bath particle current. In the lab frame, let us consider a given link between two sites, say 0 and 1, and examine the processes leading to a hop of a bath particle across this link.

• One possibility is that the tracer occupies site 0, that site 1 is vacant and that site 2 is occupied by a bath particle. In that case an exchange can occur and the bath particle can cross the designated link from right to left. The required configuration occurs with probability $\frac{\langle (1-\tau_1) \tau_2 \rangle}{L}$, the transition happens with a rate $\newcommand{\pp}{p'} \pp$ leads to an algebraic current  −1 across the link.
• A similar contribution is obtained from the case where the tracer occupies  −1, 0 is vacant and 1 is occupied, giving another factor $\newcommand{\pp}{p'} -\pp \frac{\langle (1-\tau_1) \tau_2 \rangle}{L}$. Two symmetric processes happen with rate $\newcommand{\qp}{q'} \qp$ instead of $\newcommand{\pp}{p'} \pp$, giving $\newcommand{\qp}{q'} 2 \qp \frac{\langle (1-\tau_{L-1}) \tau_{L-2} \rangle}{L}$.
• If 0 is occupied and 1 is vacant, the bath particle can hop with rate 1. This can happen if, initially, the bath particle occupies one of the positions $l=1, \ldots, L-2$ in the tracer frame, each one with probability $\frac{1}{L}$. For each l this gives a factor $\frac{\langle \tau_l (1-\tau_{l+1}) \rangle}{L}$.
• If 0 is vacant and 1 occupied, a symmetrical process takes place. This gives a factor $\frac{\langle \tau_l (1-\tau_{l-1}) \rangle}{L}$ for each $l=2, \ldots, L-1$.

Adding all the factors, we get

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Jb}{\mathcal{J}_{\rm B}} \newcommand{\qp}{q'} \displaystyle \label{eq:Jbproof} \Jb &= \frac{1}{L} \langle -2 \qp (1-\tau_{1}) \tau_{2} + 2 \qp (1-\tau_{L-1}) \tau_{L-2} \rangle \nonumber \\ & \quad+ \frac{1}{L} \langle \sum_{l=1}^{L-2} \tau_l (1-\tau_{l+1}) - \sum_{l=2}^{L-1} \tau_l (1-\tau_{l-1}) \rangle. \nonumber \end{align} \tag{ B.1 }$$

The sums on the second line cancel out except for the difference $\tau_1 - \tau_{L-1}$ and (B.1) eventually gives expression (17) after using the mean-field approximation.

Here we show that $\newcommand{\muee}{\mu_{{\rm E},{\rm E}}} \muee$ as obtained in (14) is always strictly positive for $\newcommand{\rp}{r'} \newcommand{\rhob}{\overline{\rho}} (\rhob, r, \rp) \in (0, 1) \times (0, \infty) \times (0, \infty)$. The numerator obviously vanishes at $\newcommand{\rhob}{\overline{\rho}} \rhob =0$, $\newcommand{\rhob}{\overline{\rho}} \rhob = 1$ and r  =  0 and is strictly positive otherwise. We express the denominator as a function of $\newcommand{\rhob}{\overline{\rho}} U = r (1-\rhob)$ and $\newcommand{\rp}{r'} \newcommand{\rhob}{\overline{\rho}} V = \rp \rhob (1-\rhob)$,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rhob}{\overline{\rho}} \displaystyle \label{eq:denommuee} D = (2 \rhob - 1) (\frac{U}{1-\rhob}+\frac{2 \rhob-1}{\rhob(1-\rhob)} V) \sqrt{1+U} + U \sqrt{1+U + 4 V}, \nonumber \end{align} \tag{ C.1 }$$

and we will show that D  >  0 for all $\newcommand{\rhob}{\overline{\rho}} (\rhob, U, V) \in (0, 1) \times (0, \infty) \times (0, \infty)$. For fixed U and V, $D \to + \infty$ for both $\newcommand{\rhob}{\overline{\rho}} \rhob \to 0$ and $\newcommand{\rhob}{\overline{\rho}} \rhob \to 1$. The denominator D reaches a minimum inside the interval, for

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rhob}{\overline{\rho}} \displaystyle \label{eq:rhomin} \rhob_{\rm min} = \frac{-V + \sqrt{U V + V^2}}{U}. \nonumber \end{align} \tag{ C.2 }$$

The value of this minimum is

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\rhob}{\overline{\rho}} \displaystyle \label{eq:Dmin} \min_{\rhob \in (0,1)} D = &\,U \left(\sqrt{1+U+4V} - \sqrt{1+U}\right) + 2 \sqrt{V} \left(\sqrt{U+V} - \sqrt{V}\right), \nonumber \end{align} \tag{ C.3 }$$

which is clearly positive, hence the positivity of $\newcommand{\muee}{\mu_{{\rm E},{\rm E}}} \muee$.

For a thermal wall it is required that a colliding particle exchanges energy with the wall such that it is reflected with a velocity corresponding to the equilibrium temperature $T_{\rm b}$ of the wall. In the simplest version only the normal component of the velocity, $v_n = {\bf v} \cdot {\bf n}({\bf q})$, is mapped onto v_n', and the parallel component v_p remains unchanged [82, 83]. Here, ' refers to the outgoing particle, and ${\bf n}({\bf q})$ is a unit vector at position ${\bf q}$ and normal to the wall, pointing inward. Instead, we require in the following that the particle is reflected from the wall with an outgoing angle as in a specular reflection, but with its speed $\newcommand{\e}{{\rm e}} v \equiv \vert {\bf v}\vert$ mapped according to $v' = g(v) v$. To determine $g(v)$, a detailed balance condition is imposed, which requires that the collision frequencies with impact velocities ${\bf v}$ and $- {\bf v}'$ are equal [82]. For particles in equilibrium with the wall, the number of particles reflected into ${\bf v}'$ is the same as under specular reflection. The average number of collisions per unit wall area with velocity between ${\bf v}$ and ${\bf v} + d{\bf v}$ is given by [82]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle n_c({\bf v},{\bf q}) d{\bf v} = - v_n f^{(1)}({\bf v}.{\bf q}) d{\bf v}, \nonumber \end{align} \tag{ D.1 }$$

where $f^{(1)}({\bf v}, {\bf q})$ is the one-particle distribution function, which in the canonical ensemble becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f^{(1)}({\bf v},{\bf q}) = n({\bf q}) \left(\frac{m \beta}{2 \pi} \right)^{d/2} \exp \left(- \frac{\beta m v^2}{2 } \right). \label{detailed} \nonumber \end{align} \tag{ D.2 }$$

$n({\bf q})$ is the particle density at ${\bf q}$. As usual, $\beta = 1/k_{\rm B} T_{\rm b}$, k_B is the Boltzmann constant, and d is the dimension. In the following we consider a planar system, d  =  2. Then the detailed balance condition becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle -v_n\exp\left(-\frac{\beta m {\bf v}^2}2\right){\rm d}{\bf v} = v'_n\exp\left(-\frac{\beta m({\bf v}')^2}2\right) d{\bf v}'. \label{three} \nonumber \end{align} \tag{ D.3 }$$

Using polar coordinates, ${\rm d}{\bf v} = v {\rm d}v {\rm d}\phi; \;\; v_n = v \cos{\phi},$ one finds

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & v \cos(\phi) \exp\left(-\frac{\beta m {\bf v}^2}2\right) {\rm d}v {\rm d}\phi \nonumber \\ &= v' \cos(\phi') \exp\left(-\frac{\beta m({\bf v}')^2}2\right) {\rm d}v'{\rm d}\phi'. \nonumber \end{align} \tag{ D.4 }$$

The integral over the angles on both sides of this equation cancel each other. Using the dimensionless abbreviations $X = \sqrt{(m \beta / 2)} \;v$ and $X' = \sqrt{(m \beta/2)} \;v'$, a final integration on both sides yields

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle X' e^{-X'^2} - \frac{\sqrt{\pi}}{2} {\rm erf}(X') = {\rm const} - X e^{-X^2} + \frac{\sqrt{\pi}}{2} {\rm erf}(X), \label{g} \nonumber \end{align} \tag{ D.5 }$$

where ${\rm erf}$ denotes the error function. The integration constant

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle {\rm const} = - \frac{\sqrt{\pi}}{2}\nonumber \end{align*}$$

is determined from the requirement that $v \to 0$ is mapped into $v' \to \infty$, and $v \to \infty$ into $v' \to 0$. For given X, respective v, a numerical solution of equation (D.5) yields the map $v' = g(v) v$. The new velocity components normal and parallel to the wall after the collision become

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle v'_n = g(v) v_n; \;\;\; v'_p = - g(v) v_p. \nonumber \end{align} \tag{ D.6 }$$

This map is time reversible and deterministic, $(- {\bf v}')' = {\bf v}$. For large (small) velocities the energy transfer due to a wall collision is negative (positive). Since this condition prevails whenever the instantaneous kinetic temperature of the gas is larger (smaller) than the wall temperature, $T > T_{\rm b}$ (respective $T < T_{\rm b}$), energy is transferred from the gas (wall) to the wall (gas) as required for a thermostat.

A crucial test is the distribution function of the velocity component perpendicular to the wall, which in equilibrium is expected to behave as [84]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Phi(v_n) = m \beta \vert v_n\vert \exp \left(- \frac{m \beta v_n^2}{2 } \right) \label{d2} \nonumber \end{align} \tag{ D.7 }$$

The numerically obtained distribution is in full agreement with this theoretical prediction.