Uphill diffusions in single and multi-species systems

Uphill diffusions constitute an intriguing phenomenon reported in a series of numerical simulations and experiments in which particles move from lower to higher density regions, at variance with the basic tenets of transport theory. In this paper we review several examples of uphill diffusions that appear in quite different frameworks. We highlight the role of the coupling with external reservoirs in the onset of particle currents with the ‘wrong’ sign, and also put forward a statistical mechanical explanation of the phenomenon for stochastic multi-species systems as well as for single-species models undergoing a phase transition.

Uphill diffusions constitute an intriguing phenomenon reported in a series of numerical simulations and experiments in which particles move from lower to higher density regions, at variance with the basic tenets of transport theory. In this paper we review several examples of uphill diffusions that appear in quite different frameworks. We highlight the role of the coupling with external reservoirs in the onset of particle currents with the 'wrong' sign, and also put forward a statistical mechanical explanation of the phenomenon for stochastic multi-species systems as well as for single-species models undergoing a phase transition.
Keywords: uphill diffusions, phase transitions, Kac potential, Fick's Law, cellular automaton (Some figures may appear in colour only in the online journal) * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Diffusion phenomena underlie a large variety of chemical, physical and biological processes, and their mathematical characterization from a molecular theory stands as one of the great achievements of statistical mechanics [1]. The derivation of the classical diffusion law was outlined by Einstein in his seminal paper [2], which is commonly regarded as one of the first attempts to shed light on the riddle of the Brownian motion. An even earlier contribution is due to Bachelier, who formulated the mathematical theory of Brownian dynamics in a work which went rather unnoticed for several decades [3].
The rigorous formulation of the scaling limit of the random motion of interacting particles on a discrete lattice has then gradually been set on firm grounds in probability theory, and we have thus learned how the details of the microscopic dynamics determine the structure of the macroscopic diffusive behavior, see e.g. [4,5].
At a phenomenological level, the classical diffusion formula can be obtained heuristically from the law of conservation of mass, by stipulating a linear relationship between the mass flux J and the gradient of the mass density ρ, J = −D∇ρ, with D > 0 the diffusion coefficient, which is usually called the Fick's law [6][7][8][9]. 5 The minus sign in the law ensures that the mass flux goes 'downhill', i.e. from regions of higher to lower density. When the flux J points in the opposite direction (namely, it goes along the density gradient), the transport is said to be 'uphill', which seemingly breaks the common intuitive picture of diffusion phenomena.
The first experimental evidence of uphill diffusions dates back to an old paper by Hartley [10], discussing the diffusion of a molecular solute in a solvent containing a gradient of concentration of another species (a salt). This was followed by a famous paper by Darken [11], where the author reports on the diffusion of carbon atoms in pairs of steel welded at one end, containing the same amount of carbon and different alloy content. The phenomenon described by Darken provides an experimental realization of an uphill migration of particles in a multispecies system, arising as a result of the competition between the density gradients of each species. More recently, the study of uphill diffusions for multi-species systems in a thermodynamic framework was treated in the review by Krishna [12], see also [13][14][15].
Another paradigmatic instance of uphill transport concerns the spinodal decomposition, a process in which a phase separation spontaneously occurs inside the unstable region of the phase diagram, without any nucleation mechanism [16,17]. In [18] the Authors characterize the spinodal decomposition studying a one-dimensional stochastic particle model on a torus, in which particles interact via a Kac potential and migrate uphill with respect to the concentration gradient. A remarkable result, outlined in [18], is the derivation of the macroscopic scaling limit of the model, expressed by a gradient flow evolution corresponding to a non-local version of the Cahn-Hilliard equation. Further work on the same and related particle models describing the spinodal decomposition can be found e.g. in [19,20].
In the last few decades, stochastic particle dynamics have provided a useful playground to evidence further breakdowns of the Fick's law. A growing research activity has focused on the investigation of stochastic multi-species particle systems. Specifically, in [21,22] 5 Different formulations of the Fick's law are available in the literature depending on the specific scientific community. While the original law in [6] employs the concentration of a given species (or, equivalently, the color density as in [7]), in this work we may also refer to the total density of the system, as in [23].
single-file diffusions are modeled via two-component simple exclusion processes coupled to external particle reservoirs. In these works it is shown that a careful tuning of the reservoir densities of each species can produce a boundary layer inside which particles move uphill. Moreover, the hydrodynamic limit of the model is also derived, along with the time evolution of the density profiles of the two species and the corresponding stationary behavior.
Further discussion of multi-component systems, whose hydrodynamic limits are studied in the perspective of uphill diffusion, is presented in [23,24]. In [23] the particles are allowed to switch randomly between two one-dimensional lattices (i.e. layers that identify the species the particles belong to) and move on each lattice with jumping rates that depend on the layer. This model mimics a reaction-diffusion dynamics which leads, in the hydrodynamic limit, to the so-called weakly coupled reaction-diffusion equation. In the presence of external reservoirs, the model can generate boundary layers and uphill diffusion. The two-species model studied in [24] evolves on a single layer characterized by a hard-core exclusion rule and by a random switch between species and boundary reservoirs. The uphill diffusion of one species is observed with finite-size lattices, but is lost in the hydrodynamic limit.
While the interactions between different species are known to give rise to counter-intuitive behaviors in biological, chemical and even social systems [25], less studied is the case of the stationary uphill diffusion in single-species models. A preliminary characterization of this phenomenon is outlined in [26], where a stochastic cellular automaton (CA) gives rise to an uphill displacement of particles in a boundary driven system. A detailed investigation of the uphill diffusion in single-species models also appears in a series of papers dealing with onedimensional particle systems with Kac potentials [27][28][29] and with two-dimensional lattice gases related to the classical Ising model [30,31].
Another example of currents equipped with the 'wrong' sign comes from systems displaying an absolute negative mobility (ANM) [32]. This refers to a migration mechanism in which particles tend to move in a direction opposite to an external driving force [33,34], and whose origin lies in the interplay between the geometry of the system and the imposed driving. Instead, in the case of uphill diffusion in single-species systems, a key role is played by the presence of a phase transition and by the coupling mechanism with the external reservoirs. Thus, despite sharing some analogies, ANM and uphill diffusion seem to originate from different microscopic mechanisms and their connection, if any, requires further investigation.
The wealth of examples and phenomena mentioned above, witnessing traces of violation of the Fick's law, demonstrates that the realm of uphill diffusions is diverse and still partially uncharted. This review work aims, therefore, at providing a contribution to the development and systematization of a statistical mechanical theory of uphill transport in the setup of stochastic particle dynamics. Namely, our tendency here-which somewhat mirrors our own research interests-is to study uphill currents by looking at single or multi-species interacting particle systems in contact with external reservoirs. For single-species systems, the stationary flux of particles is said to move uphill if the particles migrate across the system from the reservoir with lower density to the reservoir with higher density. In the case of multi-component systems, the previous definition must be suitably adapted to account for the interaction between different species. For each model in the subsequent sections we will thus specify which quantity goes uphill.
The paper is structured as follows. In section 2 we study the appearance of an uphill transport in multi-species models and focus on the macroscopic transport equations obtained, in the scaling limit, from the analysis of interacting particle systems. In section 3 we review the onset of the uphill diffusion in a 2D Ising model driven out of equilibrium through the coupling with external reservoirs. Section 4 deals with a 1D stochastic CA, in which particles interact via a Kac potential and their uphill migration is sustained by a phase transition. It is worth pointing out that, from the mathematical standpoint, a twospecies system like the one in [22] can be interpreted and studied as a single-species system endowed with two velocities, as in the case of the CA. Nevertheless, the phase transition induced by the Kac potential requires a specific mathematical treatment, like the one carried out in [26][27][28]. In section 5 we review the presence of the uphill diffusion in a class of interacting particle systems known as zero range processes. Conclusions are drawn in section 6.

Uphill diffusion in multi-species systems
As mentioned in the Introduction, the typical experimental settings where the uphill diffusion occurs are multi-species systems. In this context such 'anomalous' behavior is to be ascribed to the interaction between the different species. Microscopic models like interacting particle systems are probably the most natural and suitable tools for modeling and studying the mechanisms that lead to uphill transport. On the other hand, real-life phenomena are observed on a macroscopic scale that is classically described by ordinary differential equations (ODEs) or partial differential equations (PDEs). In this section we review some papers in which uphill diffusion and ANM are studied under this twofold perspective. We start by considering systems described by classical linear reaction-diffusion PDEs, then we will focus on the relation between the existence of hydrodynamic limits and uphill currents.

Macroscopic reaction-diffusion equations for a two-species system
Overview of the model: In this section a two-species reaction-diffusion particle system with hard-core interaction is considered. The particles lie on a one-dimensional lattice and reach a non-equilibrium steady state (NESS) sustained by the imbalance between the densities at the right and left boundaries for both species. In the hydrodynamic limit the densities of the two species, ρ (1) and ρ (2) , satisfy a reaction-diffusion equation expressed by equation (3), with vanishing cross-diffusion terms. This fact strongly affects the chance to observe uphill transport behaviors. Indeed, the total density, i.e. ρ = ρ (1) + ρ (2) , may be associated with a stationary current J directed from the reservoir with lower total density to the one with higher total density (a phenomenon called 'global uphill diffusion' in [24]). On the other hand, the situation in which a single species, say 1, has a current J (1) that goes uphill with respect to the boundary density of the same species (called 'partial uphill diffusion' in [24]) can not be realized.
In [24] a model with two species, called 1 and 2, diffusing on the unit interval is considered. In the hydrodynamic limit, the densities of the two species assumed to be smooth. The existence of uphill diffusion stems from a linear reaction-diffusion equation of the form with Dirichlet boundary conditions ρ 0 (x). The coefficients appearing in (1) are assumed to be non-negative and the diffusivity matrix positive definite. In this model the evolution of the densities is driven by two mechanisms: 1) a crossdiffusion, in which the diffusion of a species is not only proportional to its own Laplacian but also to the Laplacian of the other species, and 2) a linear reaction mechanism which makes a species switch into the other one with intensity Υ ⩾ 0. A couple of special cases are worth being considered: • Assuming in (1) a diagonal diffusivity matrix (with σ 11 , σ 22 > 0) without reaction (Υ = 0), a decoupled system is obtained in which the two species follow a standard transport process (described by the Fick's law); in this case uphill diffusion can not be achieved. • A more interesting coupled reaction-diffusion model, known as the double diffusivity model, is obtained by switching on the reaction part (Υ > 0), while keeping the diffusions decoupled (σ 12 = σ 21 = 0): ) . (3) This model has been extensively studied in literature, and in [23] it has been obtained as the hydrodynamic limit of a switching interacting particle system in which uphill diffusion (with boundary layers) occurs, see section 2.2. for more details.
Coming back to the general model (1), Casini et al [24] introduce the stationary diffusive currents of the two species as follows: and observe that one can identify a priori two different phenomena, referred to as global uphill diffusion and partial uphill diffusion. Let ρ L = ρ L be the total density at the left boundary and ρ R = ρ (1) R + ρ (2) R the one at the right boundary, and also denote by The partial uphill diffusion for the i th species occurs when the current of species i flows along the gradient of its boundary density, that is: when In what follows, we will use the terminology global uphill and partial uphill for the sake of brevity. An example of global uphill is obtained in the double diffusivity model (3), provided that σ 11 ̸ = σ 22 [23].
As shown in [24], partial uphill diffusion is possible for the model (1) with a non-diagonal diffusivity matrix Σ. For instance, by choosing σ 11 = σ 22 = Υ = 1 and σ 21 = σ 12 = 1 2 , the diffusive currents take the form 2 (e 4 − 1) for i = 1, 2. Assuming ρ (1) R , the search for partial uphill for, say, the species 1, can be formulated by seeking for boundary densities (ρ It can be seen that there are many choices of boundary densities that allow for partial uphill diffusion of the species 1 [24]. Two examples are given in figures 1 and 2. In the first example, the current J (1) (x) is non-monotonic, while in the second one it is monotonic. In both cases the current of the species 1 attains its positive minimum at the boundary.
The fact that in the previous example σ 12 and σ 21 do not vanish is not incidental, since in [24] Casini et al state that partial uphill can occur only if cross-diffusion is present, i.e. if the matrix Σ is non-diagonal (and with a non-vanishing reaction coefficient Υ).
This issue is connected with the possibility to obtain (1) as the hydrodynamic limit of a microscopic particle system. A candidate model is studied in [24], where a two-species interacting particle system with reaction and diffusion is introduced. This boundary driven interacting particle system on Λ = {1, . . . , N} is conceived in such a way that the evolution equation of the average occupation variables ρ (i) z at the sites z ∈ Λ are the space-discretized version of (1), where the unit interval is replaced by Λ and the second order derivatives by their finite difference proxies. Thus, provided that some conditions are satisfied (see next), the resulting ODE system for the densities of species ρ with z = 1, . . . , N, ρ A necessary and sufficient condition for the existence of such interacting particle system is that the diffusion coefficients (σ k,ℓ ⩾ 0) and the reaction coefficient Υ ⩾ 0 satisfy the following relations: where m ⩾ 0 and h ⩾ 0 are free parameters, see theorem 4.1 in [24] for details 6 . The study of the hydrodynamic limit for a one-parameter family of systems satisfying (7) on Z leads to a surprising result [24]: although the microscopic dynamics has non-zero crossdiffusivity terms σ 12 and σ 21 , the hydrodynamic limit is expressed by (3), i.e. it has a diagonal diffusivity matrix Σ with σ 11 = σ 22 . This fact can be traced back to a scaling of the crossdiffusivity terms that guarantees the positiveness of the rates of the rescaled generator. Casini et al [24] also argue that this is the only possible scaling leading to a PDE in the hydrodynamic limit and that Σ has necessarily the form described above. Thus, the conclusion of this analysis is that in the hydrodynamic limit of the two-species system, described by (3) with σ 11 = σ 22 , partial and global uphill are lost, even if they are present in the finite-size system [24]. This is consistent with the result of [23], where global uphill diffusion for the double diffusivity model described by (3) is found only when the diffusivities of the two species are different 7 .
As a side remark we observe that, neglecting the issue of the hydrodynamic limit and considering (1) per se, from the constraints (7) it follows that the total density ρ = ρ (1) + ρ (2) evolves according to the heat equation ∂ t ρ = σ∂ 2 x ρ, with σ = σ 11 + σ 21 . Thus, under these conditions no global uphill diffusion can be realized.

Switching interacting particle systems
Overview of the model: In this section a further example of two-species reactiondiffusion particle system is considered. Unlike the system discussed in section 2.1, in the present model several particles of different species may occupy the same site. A NESS is generated by imposing different densities at the boundaries. In the hydrodynamic limit the densities of the two species, ρ (0) and ρ (1) satisfy a reaction-diffusion equation (3) with vanishing crossdiffusion terms. As with the model of section 2.1, the total density, i.e. ρ = ρ (0) + ρ (1) , may be associated with a stationary current J directed from the reservoir with lower total density to the reservoir with higher one (a phenomenon called 'uphill diffusion' in [23]). Since J = −∇(ρ (0) + ϵρ (1) ), with 0 ⩽ ϵ ⩽ 1 denoting the hopping rate of the species 1 (the species 0 has rate 1), the Authors say that the Fick's law is not satisfied except for the case with ϵ = 1.
As mentioned above, the double diffusivity model described by the PDE system given in (3) is the macroscopic equation for the non-equilibrium density profiles of two species, namely the  fast and the slow particles, in a switching interacting particle system [23]. The particles jump on the one-dimensional lattice Λ at rate 1 (fast particles) or at rate ϵ ∈ [0, 1] (slow particles). Moreover, particles may also switch between the two species (i.e. they may change jump rate) at rate γ ∈ (0, ∞). The process has an equivalent representation in which the particles move in a latticeΛ with two layers, i.e.Λ = Λ × {0, 1}: the layer 0 holds the fast particles, while the layer 1 holds the slow ones. Switches between layers occur at rate γ. Furthermore, the cases with Λ = Z and Λ = {1, . . . , N}, are considered [23]. In the latter case, the system is driven out of equilibrium by external particle reservoirs (containing particles of both species) located at the two boundaries of the lattice. The hydrodynamic limit that leads to (3) is obtained through a standard diffusive scaling.
In the non-equilibrium setting, the reservoir of the layer ℓ ∈ {0, 1} absorbs particles at rate 1 and injects particles (of species ℓ) at rates that depend on the left and on the right densities ρ (ℓ) L and ρ (ℓ) R , see [23] for details. The resulting particle system is dual to the original system, but with the reservoirs replaced by some absorbing sites. The duality and the knowledge of the absorbing rates allow for a detailed analysis of the system. In particular, it is possible to compute the stationary measures, the stationary microscopic density profiles, the microscopic and the macroscopic currents of the two species: are the stationary macroscopic density profiles of the fast and the slow particles. The total current is constant (as it should in a stationary state) and is given by the explicit formula [23]: while the time-dependent profiles ρ 0 (x, t) , ρ 1 (x, t) satisfy the evolution equation (3): (where Υ is the limit of the switch rate γ N scaled diffusively, i.e. γ N N 2 → Υ), with boundary conditions: As a consequence of (9), the evolution of the total density ρ( satisfies the following equation: with boundary values: For ϵ = 1, i.e. when there is only one species, equation (11) reduces to the heat equation and the diffusion is standard ( i.e. it obeys the Fick's law). On the other hand, as observed in [23], when the species are different, i.e. 0 ⩽ ϵ < 1, the total flow J, which is the sum of the flows of the fast and the slow particles, is not proportional to the gradient of the total mass, since . As a consequence, when the species are different, anomalous behaviors may arise. For instance, it can be observed that non-zero currents may exist even if the total density has the same values at the boundaries, i.e. ρ R = ρ L . An example in which this phenomenon occurs is given by the case with ρ R . In fact, since the density gradient of the slow particles is opposite to that of the fast ones, the two species flow in opposite directions and, because the velocities are different, a total non-vanishing current J(x) is produced; see equation (8) with 0 ⩽ ϵ < 1. By imposing a non-zero boundary gradient to the total density, e.g. ρ R < ρ L and still assuming different velocities, uphill diffusion can occur. Indeed, while in the regime with ϵ = 1 the total density ρ(x) is linear and interpolates between the boundary values (the current J(x) is constant and positive), when 0 ⩽ ϵ < 1 the total density is non-linear, J(x) is not proportional to the gradient ∂ x ρ(x) and uphill currents appear if and only if the densities of the reservoirs ρ L satisfy the following inequality: ( The detailed analysis reported in [23] provides a quite clear picture of a possible microscopic mechanism that produces uphill diffusion at the macroscopic level 8 . This is the combination of two conditions that are both necessary: the presence of two species with different diffusivities and the mismatch of the densities of the reservoirs, that keeps the system out of equilibrium. Indeed, (13) shows that at equilibrium, i.e. when ρ L , uphill diffusion is absent.

ANM in a system close to equilibrium
Overview of the model: In the boundary driven models of sections 2.1 and 2.2 uphill currents are generated under non-equilibrium conditions. On the contrary, the model discussed in this section is a perturbation of an equilibrium system defined on a ring, where a tracer particle moves in a sea of neutral particles. Here the driving is represented by the asymmetry (bias) between the left and right hopping and exchange rates of the tracer. The observable of interest is the velocity of the tracer which, for properly chosen values of the parameters, may have a sign opposite to the bias, that acts as a thermodynamically conjugated force. This phenomenon is called 'ANM' in [34].
Evidences of currents flowing in the direction opposite to the driving can be found in systems that may differ from those described in the previous subsection. Here, we review an example [34] with two species in which the foregoing phenomenon is produced by slightly perturbing the equilibrium state. In this case, the driving is represented by the asymmetry in the rate of jump of one of the two species: one tracer and a set of bath particles. The particles occupy the sites of a discrete ring and move satisfying the simple exclusion constraint. The bath particle follows a Symmetric Simple Exclusion Process with constant rate 1, while the tracer may hop to the right and to the left with rates p = r + δ 2 and q = r − δ 2 respectively, but may also exchange its position with a bath particle two sites away, provided that the site between the tracer and the bath particle is vacant. The exchange occurs on the right with rate The biases δ and δ ′ play the role of the drift acting on the tracer. When δ = δ ′ = 0, the system is at equilibrium and the detailed balance condition is satisfied. Interesting phenomena occur when equilibrium is perturbed by introducing an asymmetry, i.e. when the biases are non-zero. Indeed, Monte Carlo (MC) simulations show that the average velocity of the tracer, v tr , is negative if δ > 0 and δ ′ = 0 for sufficiently large values of the density of bath particlesρ. This phenomenon, called ANM in [34], is studied in the linear response regime via a mean-field approximation. Decomposing the mean velocity of the tracer as the sum v tr = v tr,H + v tr,E of the contributions coming from hops towards the empty sites (v tr,H ) and from exchanges with the bath particles (v tr,E ), the following linear response relation is obtained: ( where c > 0 is a constant. This relation shows that δ/r and δ ′ /(2r ′ ) play the role of the conjugated forces of v tr,H and v tr,E , and the linear response matrix, which is symmetric according to the Onsager relations, has off-diagonal elements that change sign atρ = 1 2 , thus allowing for ANM. As it is observed in [34], '. . .the ANM found in the linear response regime is a direct result of the fact that dynamics involves two driving mechanisms, namely hoping (p/q) and exchange (p ′ /q ′ )'.
A subtler phenomenon, called absolute differential mobility (ADM), can be observed for both small and large densities. Essentially, this occurs when the response of the quantity of interest depends in a non-monotonic way on the driving. For instance, as δ is increased from zero while keepingρ = 0.25 and δ ′ = 0 fixed, the velocity v tr first grows monotonically and then, after reaching a critical value, it decreases monotonically towards zero. At large densities, e.g.ρ = 0.7, ANM and ADM occur simultaneously: the velocity is directed uphill and its absolute value is a decreasing function of δ.

Uphill diffusion in the non-equilibrium 2D Ising model
Overview of the model: In this section we review the model considered in [30,31], where the Authors study the diffusion of positive spins on a two-dimensional lattice. The driving that keeps the system out of equilibrium, thus producing a non-zero steady current J of positive spins, corresponds to the imbalance between the magnetizations m + > m − fixed by the reservoirs at the right and left vertical boundaries, respectively. Uphill diffusion is realized when J > 0, namely when, on average, the positive spins move towards the right reservoir, which possesses a higher positive magnetization m + . The simulations in the cited papers suggest that the presence of a phase transition (in the corresponding equilibrium 2D Ising model) is a necessary condition for the occurrence of uphill diffusion. Indeed, above the critical temperature, uphill diffusion is absent.
As discussed in the previous sections, usually uphill diffusion occurs in multi-component systems as the result of the competition between the density gradients of the different components of the mixture. However, the presence of a mixture is not necessary to produce uphill currents. Indeed, in the sequel we will see that this phenomenon can also occur in singlecomponent systems that exhibit a phase transition. In section 4 we will consider a singlecomponent 1D system [26][27][28] in which a phase transition takes place in the mean field limit, while in the present section we address the diffusive behavior of a non-equilibrium 2D Ising model. The latter model, numerically studied in [30], is one of the first examples in the literature of a system undergoing a phase transition and attaining a NESS with uphill currents 9 .
In the infinite volume limit, the equilibrium 2D Ising model undergoes a phase transition at the inverse critical temperature [36] with a spontaneous magnetization given by [37] for β > β c , and m β = 0 for β < β c (in both cases with vanishing external magnetic field). In [30] it was shown that, when coupled at the boundaries with two reservoirs at different magnetizations and also endowed with a conservative exchange dynamics in the bulk, the system exhibits stationary uphill magnetization currents. Obviously, the nature of the NESS depends on the various parameters of the model (e.g. the size of the system, the inverse temperature, the magnetization of the reservoirs), as well as on the imposed boundary conditions, on the details of the microscopic dynamics in the bulk and on the coupling with the external reservoirs. Thus, many non-trivial questions may be posed regarding the appearance of stationary uphill currents in the system, mostly dealing with the identification of the microscopic mechanism producing such anomalous behavior. The insight into these problems seems to be currently out of reach of the analytical approach, while numerical simulations are affected by finite-size phenomena that can hardly be assessed, since the workload needed to properly sample the NESS is quite large even for moderate system volumes [30,31]. Therefore, the crucial question arises as to whether uphill currents are 'real' phenomena or they rather represent a manifestation of finite-volume effects, as claimed on the basis of an heuristic argument in [38]. In this respect, different boundary conditions and various coupling mechanisms with the reservoirs were considered in [31], where the role of the detailed balance condition at the interface bulk/reservoirs was also analyzed.

The model
In the present subsection we deal with the setup studied in [30], in which the nearest-neighbor ferromagnetic Ising model on the square lattice Λ = {1, . . . , L} 2 is coupled at the left and right vertical boundaries to two magnetization reservoirs, see figure 3. The microscopic state The energy of interaction between nearest-neighbor spins is given by the standard 2D Ising Hamiltonian at zero external field: At a fixed inverse temperature β > 0, the spins evolve following a continuous-time stochastic dynamics with two contributions: a conservative exchange dynamics in the bulk and independent spin flips (ISF) at the boundaries. In the bulk the spins follow a Kawasaki dynamics, namely the spins of a bond ⟨i, j⟩ (i.e. with i and j such that |i − j| = 1) exchange values at rate: where σ ij denotes the configuration obtained from σ by exchanging the spins at the sites i and j. The system is driven out of equilibrium by the presence of two infinite reservoirs adjacent to the left and right vertical boundaries, denoted by R − and R + respectively, that force a magnetization m + ∈ [0, 1] on the right column and a magnetization m − = −m + on the left column, see figure 3. The mechanism that simulates the action of the reservoirs, called ISF, is obtained by changing the spins σ i at the boundaries independently with rates: with y = 1, . . . , L. Two types of boundary conditions are considered in [30]: 4) . No significant difference in the results obtained using the two boundary conditions is detected in [30]. Other boundary conditions and reservoir mechanisms will be considered below in section 3.3.
The resulting dynamics, due to the presence of the magnetization reservoirs at the vertical boundaries, is not reversible w.r.t. the Boltzmann-Gibbs measure with Hamiltonian (16) and the NESS that sets in is characterized by a uniform magnetization current flowing in the horizontal direction. The integrated current J t (b), measured along a given horizontal bond b, can be computed counting the number of positive spins that cross b from the left to the right minus the number of those crossing the same bond in the opposite direction, and then dividing the resulting number by the elapsed time of the dynamics t. The stationary current J corresponds to the long time limit of this rate, i.e. J = lim t→∞ J t (b)/t, which turns out being independent of the bond in the stationary state. Thus, in the setup shown in figure 3, J > 0 characterizes a NESS with uphill currents.
We conclude this section by remarking that the Ising model is equivalent to a lattice gas model via the transformation η i = (1 + σ i )/2 that maps the spin variables σ i into the occupation variables η i ∈ {0, 1}, with η i = 1 (resp. η i = 0) denoting the presence (resp. the absence) of a particle. This way, one can easily translate the jargon of magnetic systems into that of lattice gases and, since the magnetization per spin can be identified with a mass density, the magnetization flux turns into a mass flux. Therefore, in the present section we will call (with some abuse of terminology) Fick's law the property according to which the magnetization flows in the direction opposite to the magnetization gradient.

Structure of the NESS
The continuous-time Markov process generated by the Kawasaki dynamics with the reservoirs was studied in [30,31] via the MC method. Simulations were performed on models with several linear sizes L ⩽ 40. We review, here, the results concerning the largest lattice (i.e. L = 40), which are not severely affected by the finite-volume effects observed in systems with smaller volume.
The nature of the NESS attained by the model (and hence the presence of uphill currents) depends in a rather complicated way on the magnetization m + of the right reservoir (recall that m − = −m + ) and on the inverse temperature β. Specifically, the presence of a phase transition for the 2D Ising model at equilibrium is expected to affect crucially also the evolution of the corresponding boundary driven model. In fact, also in the non-equilibrium setting, the regimes of high (β < β c ) and low (β > β c ) temperature lead to different scenarios characterized by standard diffusion (described by the Fick's law) in the former case, and more complex transport phenomena in the latter one. As representative instances, two different regimes have been studied in [30,31] : β = 0.3 and β = 1, respectively lower and higher than the critical value β c given in equation (14).
The main results of [30], extended and deepened in [31], concern the analysis of the stationary current J(m + ) of the model, introduced in section 3.1. The behavior of the NESS and the current J(m + ), as functions of the right reservoir magnetization m + , can be summarized as follows.
(1) Uphill diffusion. The phenomenon is highlighted in figure 4, where the stationary current J(m + ) is computed for β = 1, i.e. in the low temperature regime. As m + is decreased from its maximum value 1, the flux J(m + ) is initially negative and then, after reaching a critical value m crit where J(m crit ) = 0, it becomes positive. More explicitely: for m + > m crit the current is negative, as it flows from the reservoir with positive magnetization R + to the one with negative magnetization R − (in accordance with the Fick's law), while for m + < m crit the current changes sign and flows from R − to R + (in violation of the Fick's law). Thus, uphill diffusion appears in the system. The estimate of the critical magnetization given in [30] is m crit ∼ 0.999 31, which is quite close to the equilibrium magnetization m β=1 ≈ 0.999 27. This fact suggests the quite natural guess that m crit = m β (up to finite-size effects). In [30] it is claimed that m crit can be evaluated in a independent way by measuring the magnetization m eq of the rightmost column of Λ in the absence of reservoirs. This quantity, in the presence of L4 b.c. and in the large volume limit, is conjectured to approach m β . Some numerical evidence underpinning this claim is obtained in [30] by simulating systems with increasing volumes. On the other hand, in [31] it is shown that the two quantities m crit and m eq , evaluated for some L and β, generally do not coincide, their deviation being negligible only at high temperature and with suitable boundary conditions.
(2) Spin configurations in the NESS. The transition from downhill to uphill diffusion (we are still discussing the case β = 1) is correlated with a change in the structure of the NESS occurring at the critical magnetization m crit . In this non-equilibrium setting two plus and minus phases appear close to R + and R − respectively, and an interface between them is formed. This can be detected by studying the stationary magnetization profile m x , which is computed as where m x (t) = 1 L ∑ L y=1 σ (x,y) (t) is the mean magnetization along the column x at time t, and T is the total number of MC steps in the stationary state. In [30] two typical magnetization profiles, corresponding to two different spin configurations, have been identified: the instanton and the bump.
• The instanton profile is placed in the middle of the lattice and corresponds to the presence of two equally sized phases close to R + and R − (see also [40] for an interpretation of the interface based on non-equilibrium thermodynamics). In the continuous limit representation,  • The bump profile corresponds to typical spin configurations with an interface close to one of the vertical boundaries. In these configurations a sea of positive spins emanating from R + occupies the most part of the lattice Λ and leaves a small strip of negative spins abutting on R − . For x ranging from the right (x = L) to the left (x = 1) boundary, the stationary magnetization profile m x increases steadily from m + up to a maximum value m bump > m + , which is reached for x close to 1. Then m x suddenly drops from m bump to m − at x = 1, forming a boundary layer close to the left boundary 10 . See the right panel of figure 5 for a sketch of the bump profile. In the presence of a bump the current flows in the 'wrong' direction, producing uphill diffusion. We point out that such 'wrong' behavior of the current is consistent with the magnetization profile. These findings are obtained by adopting an instanton-like initial condition (i.e. σ (x,y) = −1 for x ∈ [1, L/2] and σ (x,y) = 1 for x ∈ (L/2, L]), called configuration A in [31]. As a matter of fact, the results do not depend on this choice [30].
(3) Stability of the instanton. The results of the extensive numerical investigations of the low temperature phase [30,31], suggest that for sufficiently large values of m + the instanton profile is stable. This means that there should exist a critical value m crit > 0 such that perturbing the configuration A initial condition, or even taking a random initial condition, when m + ∈ (m crit , 1) the dynamics leads to the instanton profile in the large time limit. This is called the stable phase. Quite different and more intricate is the scenario occurring for m + below m crit . In this regime the instanton loses stability while, at least for not too small values of m + , the bump profile is stable. This hints at the fact that there should exist another critical value m u < m crit marking the domain of the stability region of the bump, i.e. the interval (m u , m crit ). The nature of this phase is not well understood, but it is conjectured that the presence of m u should be regarded as a finite-volume effect, since some heuristic reasoning leads to the estimate |m crit − m u | ∼ O(L −2/3 ) [38], see also section 3.4. This phase has been called metastable in [30]. For m + < m u the picture is even more complex and less studied. Numerical evidence shows that in the interval (0, m u ) the system undergoes an intricate sequence of transitions leading to stationary states that can be associated neither with the instanton nor with the bump profiles. This is the unstable phase, also called 'weakly unstable' or 'chaotic region' in [30].

4) Overview of the low temperature phase.
An overall view of the phenomenology as m + is decreased from 1 towards 0 (at β = 1) is shown in figure 7, taken from [31]. Since the NESS undergoes the relevant transitions in a narrow interval of m + (notice, for instance, that and R − . Other and more complex spin configurations are realized for lower values of m + . in the unstable region of the system. (5) High temperature regime. The system at high temperature (β < β c ) was studied in [31] assuming β = 0.3. In this regime the phenomenology appears much simpler than it is observed at low temperature: no evidence of transitions through different phases is detected, while the non-equilibrium current is always negative for any m + . Therefore, magnetization flows downhill according to the Fick's law. Figure 8, taken from [31], shows that the current magnetization J(m + ) behaves as a smooth decreasing function of m + . For vanishing m + the system approaches equilibrium, since m + and m − = −m + are close to each other. Correspondingly, the flux tends to 0 and the stationary magnetization profile m x becomes flat. In this limit, the positive and negative phases are weakly separated, and they appear more and more random and intertwined: as a consequence, the sharp interface is lost and the structure of the instanton profile disappears. Figure 8 evidences that at high temperature the current is negative: the magnetization flows downhill in agreement with the Fick's law.

The role of boundary conditions and updating mechanisms at the boundaries
The existence of uphill currents, with boundary conditions and reservoir mechanisms different from those used in [30], was investigated in [31]. Beside the 'L4' boundary conditions (or 'shifted boundary conditions'), already exploited in [30], new boundary conditions were considered in [31]: the so-called 'free boundary' (FB) condition, in which one takes and the 'stochastic boundary' (SB) condition: Note that the properties of the reservoirs do not enter directly the L4 and FB boundary conditions. Differently, the SB rule imposes to the ghost spin σ o i an average value that coincides with the magnetization fixed by the nearby reservoir.
The reservoir spin-updating mechanisms considered in [31] are the ISF cited above, see equation (17), and the so-called Detailed Balance Dynamics (DB) mechanism [41,42] in which a spin value σ ′ i is introduced at the boundary site i at a rate given by: where σ ′ denotes the configuration obtained from σ by updating the spin σ i to σ ′ i , and the term (1 + m res σ ′ i )/2 corresponds to the probability of drawing at random the spin σ ′ i from the reservoir with magnetization m res . This rate corresponds to a Metropolis-Hastings algorithm for the Boltzmann-Gibbs distribution with hamiltonian H(σ ′ |σ o ), in which the proposal probability of the spin flip σ → σ ′ (involving a boundary site i) is (1 + m res σ ′ i )/2. In [31] the Authors consider also a model that is fully decoupled from the reservoirs. This dynamics, called no spin flip (NSF), obviously conserves the total magnetization that is fixed by the initial configuration.
By combining reservoir mechanisms and boundary conditions, nine different models are defined in [31]: a set of three equilibrium models E ≡ {NSF-L4, NSF-FB, NSF-SB}, a set containing DB reservoirs L ≡ {DB-L4, DB-FB, DB-SB}, and the corresponding three models with ISF reservoirs N ≡ {ISF-L4, ISF-FB, ISF-SB}. For the L and N models a detailed analysis of the current J(m + ) in the parameter plane spanned by m + and β has been performed and the critical values m crit (β) at which the transitions from downhill to uphill currents occurs have been computed. Numerical evidence shows [31] that the ISF-L4 and ISF-SB models share a similar behavior, owing to the presence of broad regions in the (m + , β)-plane in which uphill currents are visible, whereas in the ISF-FB model the appearance of such currents is hindered. Moreover, the critical values m crit (β) in the first two models coincide within numerical approximation and are appreciably larger than the corresponding value measured in the ISF-FB model. In any case, positive currents occur provided that β > β crit , where β crit is a critical value that, due to finite-size effects, might significantly differ from the critical value β c in (14).
The key role played by the reservoir dynamics in the rise of positive currents can be also evinced by comparing the behavior of the N -models with the rather different phenomenology displayed by L-models. Indeed, with the L4 and FB boundary conditions the DB dynamics give rise to downhill currents in the whole domain of the parameter plane examined in [31], while with the SB boundary conditions (and with the same reservoir magnetizations) uphill currents can instead be observed. It is worth noting that, even if the dynamics with rate (22) satisfies locally the detailed balance condition with respect to the Hamiltonian (16) with fixed boundary conditions σ o , it does not match with the bulk dynamics and with the SB boundary conditions σ o that are changing during the MC evolution. Indeed, in this setting the system is evolving under a random sequence of time-dependent Hamiltonians H(σ|σ o (t)), where σ o (t) denotes a random sequence of independent configurations with distribution (21). This prevents the detailed balance condition from being satisfied, even locally 11 .

Hydrodynamics, finite size effects
The hydrodynamics of the model, i.e. the description at large length and time scales, is discussed heuristically in [38], since a rigorous analytical approach is presently not available. In the hydrodynamic limit, the state of the system is represented by the macroscopic magnetization profile at time t > 0, M t (r 1 , r 2 ), defined on the unit square [0, 1] 2 that is the limit of the lattice Λ. Owing to the vertical symmetry induced by the periodic boundary conditions, it is assumed that the profile is a function m(r 1 , t) of the sole horizontal macroscopic coordinate, i.e. m(r 1 , t) = M t (r 1 , r 2 ). In [38] the profile m(r 1 , t) is conjectured to satisfy a PDE that is parameterized by the temperature.
(1) High temperature, i.e. 0 ⩽ β < β c . In this regime there is no phase separation and the macroscopic profile is the unique solution of a non-linear heat equation, with Dirichelet boundary conditions: where D(m) > 0 is the bulk diffusion coefficient (given by the Green Kubo formula [44]) and m 0 (r 1 ) is an initial profile. In the stationary state, the magnetization m = m(r 1 ) is expected to be the unique solution of J = −D(m) dm dr1 = const. < 0, m(0) = −m + and m(1) = m + , with a profile continuously interpolating between −m + and m + , as sketched in the left panel of figure 9. Thus, in this regime the Fick's law is satisfied and the currents go downhill, no matter if m + is larger or smaller than m β . These conclusions, whose proof should follow from [42][43][44], are in good agreement with the numerical findings of [30,31].
At infinite temperature, i.e. β = 0, the process degenerates to the symmetric exclusion process with a constant bulk diffusivity D, so that (23) becomes the linear heat equation with Dirichlet boundary conditions.
(2) Low temperature, i.e β > β c . In this regime phase-coexistence with regions (interfaces) where magnetization is not slowly varying is expected. Supposing, under proper assumptions (see [38] for details), that there is only one interface at R t ∈ [0, 1], t > 0, it can be conjectured that m (r 1 , t) and R t satisfy a free boundary problem: A crucial fact in the study of the stationary magnetization m(r 1 ) is that in the low temperature regime β > β c the bulk diffusivity satisfies [44]: This property suggests that the behavior of the system should depend critically on whether m + is larger or smaller than m β . More precisely, the conjecture is the following: • With m + > m β the stationary magnetization profile has a jump (of size 2 m β ) at r 1 = 1 2 , see the central panel of figure 9. The Fick's law is satisfied and the current J = −D(m) dm dri is negative. Thus, even at low temperature the transport is downhill, provided that the reservoir magnetization is sufficiently large. This guess is qualitatively confirmed by a large amount of numerical simulations [30,31]: compare, for instance, the left panel of figure 6 (showing the instanton profile) with the theoretical prediction illustrated in the central panel of figure 9.
• With 0 < m + < m β a naive guess for the stationary magnetization profile is sketched in the right panel of figure 9. In [38] it is argued that this profile should be unstable under small perturbations of the interface and hence can not be observed in the numerical simulations. This leaves open the problem of understanding the complex nature of the NESS for m + < m β and β > β c observed in [30,31].
As was pointed out above, uphill currents can be interpreted as a finite-size effect in the nonequilibrium 2D Ising model. This was conjectured in [38,45] basing on the structure of the canonical Gibbs distribution conditioned to a magnetization m on the torus [46]. The analysis suggests that if m ∈ (m β − cL − 2 3 , m β ), with c > 0 small enough, the Gibbs distribution with magnetization m is supported by configurations with 'small' contours (i.e. of size ⩽ log L). Thus, there is no phase separation in this case and (m β − cL − 2 3 , m β ) is the plus metastable region (analogously, (−m β , −m β + cL − 2 3 ) is the minus metastable region). In [38] the presence of the bump profile ( associated with uphill currents) is ascribed to the fact that the magnetization m + of the reservoir falls in the plus metastable region. An empirical evidence in favor of this hypothesis is obtained (via numerical simulations) by observing that, keeping m + fixed and increasing L, the bump is suppressed and replaced by a more complex profile. See, for example, the double bump profile shown in figure 7 for h + = 1.15, that is associated with a downhill current.

Uphill diffusions in cellular automata
Overview of the model: In this section we review a CA studied in [26][27][28]. The model corresponds to a one-dimensional channel coupled to external particle reservoirs, in which particles are equipped with two velocities and are injected through the boundaries with specific probabilities. The dynamics constitutes a time-discrete variant of a classical Asymmetric Simple Exclusion Process, in which the driving corresponds to the difference of the densities of the two external reservoirs. The introduction of a Kac potential favors particle jumps in the direction where the particle density (within a certain mesoscopic range) is higher. Uphill diffusions are realized by choosing metastable values of the densities of the reservoirs, so that particles migrate, on average, from the lower to the higher density reservoir.
In this section we consider a remarkable instance of uphill diffusions in a 1D singlecomponent system undergoing a phase transition. The model considered in [26][27][28] is a CA simulating a particle system coupled to external mass reservoirs, which keep fixed the densities at the boundaries. In a nutshell, the phenomenon outlined in the above references can be described as follows. One calls R 1 the reservoir, with density ρ − , adjacent to the left boundary of the system, and R 2 the reservoir, with density ρ + , adjacent to the right boundary. If ρ − falls in the metastable vapor phase, while ρ + falls in the metastable liquid phase, then the system is observed to reach a steady state in which the current flows uphill from the reservoir with smaller density to the other with larger density, i.e. from R 1 to R 2 . Instead, if the densities of the two reservoirs fall in the vapor and liquid stable phases, then the current goes downhill, in agreement with the classical Fick's law, namely from R 2 to R 1 . The CA describes the stochastic evolution of particles in a 1D lattice {1, 2, . . . , L}, with L > 1 a positive integer. Beside hopping on the lattice, particles can also leave from or enter into the system through the endpoints at 1 and L. In this case, particles are said to be absorbed or injected from the reservoirs R 1 and R 2 , respectively. The two reservoirs are considered 'infinite', in the sense that they do not keep track of the number of particles which have been previously absorbed or injected into the system. The CA dynamics amounts to a parallel updating version of a weakly asymmetric simple exclusion process, which is amenable to a straightforward numerical implementation. The classical 1D symmetric simple exclusion process is a system of random walks in which particles hop to the right and to the left with equal probability, provided that the arrival site is empty. The CA studied in [26] introduces a weak asymmetry, called bias, which proves to be essential in the onset of a stationary uphill transport. Namely, the particles preferably jump in the direction where the density is higher. Note that under periodic boundary conditions this dynamics would produce a phase separation, forming a region of higher density and another one of lower density. Instead, the system considered in [26] is open, and the coupling with the external reservoirs gives rise to non-equilibrium states characterized by a phase separation and stationary (possibly uphill) currents.

The CA
Let us shortly review, here, the definition of the CA. The phase space is where v is interpreted as the velocity of a particle. Particle configurations are functions η : S → {0, 1}, η(x, v) ∈ {0, 1} denotes the occupation variable at (x, v). Here η(x) = η(x, −1) + η(x, 1) ∈ {0, 1, 2} yields the total number of particles at x. When convenient, a suffix t is also added to underline that the occupation variables are being computed at time t. There are four basic parameters entering the definition of the CA: γ −1 ∈ N, C > 0 and ρ ± ∈ [0, 1]. Specifically, γ −1 is related to the range of interaction between particles, whereas C is related to the inverse temperature β via β = 2C (the phase transition takes place for β > 1). The simulations presented in [27] refer to the case with γ −1 = 30, C = 1.25 and L = 600. The densities ρ − and ρ + of the two reservoirs R 1 and R 2 are kept fixed during a simulation, but may be changed in different simulations in order to highlight the role of the external driving exerted by the reservoirs. In [27] the Authors introduce, also, the following quantities: where Here N +,x,γ denotes the total number of particles to the right of x within distance γ −1 from x. Note that if x is sufficiently close to the right boundary, there might not be γ −1 sites located to the right of x. If there are, say, only γ −1 − m such sites, one then adds 2m fictitious phase points (y, v), with v = ±1 and where y takes m values to be regarded as m physical sites located to the right of the endpoint at L. The occupation number η(y, v) is fixed to the value ρ + , hence the fictitious m sites contribute to N +,x,γ with a term 2ρ + m. An analogous interpretation is used to compute η (−) (y), upon replacing ρ + with ρ − . The CA operates following three distinct steps. Let η denote the configuration at time t, η ′ and η ′ ′ two consecutive updates starting from η and η ′ ′ ′ the final update giving the configuration at time t + 1. The CA is defined as follows: 1. velocity flip. At all sites x ∈ {1, . . . , L} where there is only one particle, the velocity becomes +1 with probability 1 2 + ϵ x,γ and −1 with probability 1 2 − ϵ x,γ , with ϵ x,γ = Cγ 2 [N +,x,γ − N −,x,γ ] denoting the bias. At all other sites the occupation numbers are left unchanged. One then denotes by η ′ the occupation numbers after the velocity flip. 2. advection. One first deletes the particles located at (1, −1) and (L, 1) (if present), and then each of the remaining particles is moved by one lattice step in the direction of its velocity. Let, then, η ′ ′ denote the occupation numbers after this advection step. 3. exchanges with the reservoirs. With probability ρ + one inserts a particle at (L, −1) and with probability 1 − ρ + the phase space point (L, −1) is left empty. The same operations are performed at (1, 1), using the density ρ − in place of ρ + . Let the final configuration be denoted by η ′ ′ ′ .
To proceed further, it proves useful to change variables and set where σ t (x) ∈ {−1, 0, 1} and m ± ∈ [−1, 1]. The notation in equation (27) allows to interpret the particle system as a magnetic system, with σ t (x) playing the role of a spin. To achieve non-equilibrium, one requires ρ + > ρ − , and specifically one may also set ρ + + ρ − = 1, which thus implies m + > 0 and m + + m − = 0. An important quantity measured in [27] is the current, which is defined as the net number of particle exchanged, per unit time, between one reservoir and the system. The current is denoted by j R1→ch (t) and j ch→R2 (t), which correspond to the number of particles which go from R 1 to the system minus those which move from the system to R 1 in the time step t → t + 1 and, respectively, to the number of particles which go from the system to R 2 minus those moving in the opposite direction in the time step t → t + 1. Therefore, one sets: with η ′ , η ′ ′ and η ′ ′ ′ correspond to the occupation numbers, introduced above, after the three successive updates occurring in the unit time step from t to t + 1. Since the instantaneous currents j R1→ch (t) and j ch→R2 (t) turn out being strongly fluctuating, one may conveniently compute their time averages, i.e.: Stationarity is attained in the limit T → ∞, whereas the existence of the limit is guaranteed by the Birkhoff theorem. In the simulations reported in [27] the value of T is chosen empirically so as to let j T ch→R2 ≈ j T R1→ch look independent of T. The initial condition considered in [27] corresponds to leaving all phase points empty, but it was also verified that other initial conditions do not significantly alter the stationary value of the current. Fur sufficiently large values of T, the stationarity condition j T ch→R2 (m + ) − j T R1→ch (m + ) ≈ 0 holds, where typical values of the foregoing difference are of order 10 −8 , while the currents have order 10 −5 .
The results of numerical simulations of the considered CA are shown in figure 10, where the parameter T is fixed to the values T = 3 · 10 9 for m + ∈ (0, m ′ ) and m + > m ′ ′ ′ and T = 10 10 elsewhere. The black filled disks in the figure represent the values of j T ch→R2 (m + ), whereas the continuous line, denoted by j(m + ), is a continuous interpolation of j T ch→R2 (m + ), which stands as a good approximation of simulations done with the other values of m + . In [27] it is thus argued that the stationary current j(m + ) plays the role of an equation of state in a nonequilibrium setup, which is brought about by the two reservoirs with unequal magnetizations m − and m + . A striking aspect that is highlighted in figure 10 (where m ′ = 0.500, m ′ ′ = 0.825, m ′ ′ ′ = 0.912 and m iv = 0.985) is that for m + ∈ (m iv , 1] the current j(m + ) is negative (i.e. it goes downhill from m + to m − ), whereas for m + < m iv the current is positive and goes uphill from m − to m + . Furthermore, figure 10 also shows a remarkable non-monotonic behavior of j(m + ) as a function of the parameter m + . We point out that if the interactions among particles were switched off (i.e. by setting the bias to zero, ϵ x,γ ≡ 0), there would be no phase transition and the current would then flow, for any value of m + , from R 2 to R 1 .

The mesoscopic Limit
Let us turn, next, to mesoscopic units. These are obtained by scaling space and time diffusively, i.e. x → r = γx and t → τ = γ 2 t. Correspondingly, the length of the system in mesoscopic units becomes ℓ = γL. The mesoscopic limit of the model is then obtained by letting γ → 0 [26][27][28]. Note that in this limit the system becomes the real interval [0, ℓ]. One thus lets E γ denote the expectation in the CA processes (the randomness stems from the initial datum as well as from the updating rules of the CA). In [27] the following facts are assumed: 1. the limit below exists and is smooth It is then proven, in [27], that under the above assumptions the limit magnetization m(r, t) fulfills the equation: (30) is supplemented with the following boundary conditions: The solutions of the system (30) can be interpreted within the framework of statistical mechanics, which allows to connect the results of the numerical simulations presented in [27] to the classical theory of phase transitions [47].
The evolution equation (30) in [0, ℓ] equipped with periodic boundary conditions can indeed be seen as the gradient flow relative to a non-local free energy functional F(m) [48], namely for |r − r ′ | ⩽ 1 and = 0 elsewhere .
The well-known Ginzburg-Landau functional is a local approximation of F(m), in which the non-local term is replaced by a gradient squared. The corresponding gradient flow evolution becomes, then, the Cahn-Hilliard equation, which hence stands as a local approximation of (30). The key aspect, here, is that the functional F(m) conveys the thermodynamics of the system. In fact denotes the van der Waals mean field free energy, whereas its convex envelope f * * β (m) corresponds to the thermodynamic free energy. The latter is obtained by minimizing F(m)/ℓ under the constraint´m(r) = ℓs and by then letting ℓ → ∞ [47, ch 6].
It is useful to recall, also, that the equilibrium magnetization density, in the presence of a magnetic field h, obeys the mean field equation so that the magnetization in the two intervals (m * , m β ) and in (−m β , −m * ) (denoted by the two red vertical segments in figure 11) is metastable. At h = 0 the double well is symmetric and the local minima turn to global minima: these are attained at m = ±m β , with m β denoting the positive solution of equation (34) with h = 0. Therefore, ±m β are the equilibrium magnetization at the phase transition with h = 0 and β > 1. This explains why some of the characteristic parameters of the simulations are connected to the thermodynamic formalism associated to the mesoscopic equations (32). Indeed, in figure 10 the value of m iv resembles closely the value m β ≈ 0.9856 (evaluated at β = 2.5), so that the simulations reveal that the current is negative when m + is stable, namely when m + > m β , and is instead positive (i.e. it flows uphill) when m + < m β (metastable or unstable). The stationary magnetization profiles obtained in [27] by running the CA are portrayed in figures 12 and 13 in mesoscopic units. In figure 12 the value of the magnetization of the right reservoir is fixed at m + = 1 (i.e. in the stable phase), and the particle current is negative, i.e. it goes downwhill from the right to the left reservoir, as predicted by the Fick's law. Instead, in figure 13 the magnetization m + = 0.93 falls in the metastable phase, therefore the corresponding current is positive, namely it moves (uphill) from the left to the right reservoir.   , t = 10 6 (black squares) and t = 10 8 (empty circles). The black thin line denotes the initial configuration, corresponding to a step function centered at r = 5. Reproduced from [27], with permission from Springer Nature.

Uphill diffusions in ZRP
Overview of the model: In this section we review the onset of uphill diffusions in a particle system amenable to an analytical solution, known as zero range process (ZRP) [49]. The driving is represented by the imbalance of the densities of the two external reservoirs. The model also accounts for a local asymmetry in the hopping rates on the lattice. Stationary particle currents flowing from the lower to the higher density reservoir are observed when the local hopping asymmetry overcomes the driving exerted by the reservoirs. Uphill diffusions have also been studied in the setup of a ZRP in [49]. Let us shortly summarize the main findings discussed in the above reference.
One fixes a positive integer R and defines a ZRP [50] on the 1D finite lattice Λ = {1, . . . , 2R + 1} ⊂ Z, called channel in the sequel. Let the finite configuration space be Ω R = N Λ . Given n = (n 1 , . . . , n 2R+1 ) ∈ Ω R the non-negative integer n x denotes the number of particles at the site x ∈ Λ in the configuration n. One also lets u : N → R + , a positive and nondecreasing function such that u(0) = 0, denote the intensity function. Given n ∈ Ω R such that n x > 0 for some x = 1, . . . , 2R + 1, one lets n x,x±1 be the configuration obtained by moving a particle from the site x to the site x ± 1; in particular, one understands n 1,0 and n 2R+1,2R+2 to be the configurations obtained by removing a particle from the site, respectively, 1, and 2R + 1. Similarly, one denotes by n 0,1 and n 2R+2,2R+1 the configurations obtained by adding a particle to the site 1 and 2R + 1, respectively. Given p, q,p,q, α, β, γ, δ > 0, such that p + q = p +q = p + γ = q + β = 1, one sets q 1 = γ, q x = q for x = 2, . . . , R and x = R + 2, . . . , 2R + 1, q R+1 =q, p x = p for x = 1, . . . , R and x = R + 2, . . . , 2R, p R+1 =p, and p 2R+1 = β. The ZRP considered in [49] is defined as the continuous-time Markov jump process n(t) ∈ Ω R , t ⩾ 0, with rates r(n, n 0,1 ) = α and r(n, n 2R+2,2R+1 ) = δ (36) for particles injection at the boundaries, and with rates r(n, n x,x−1 ) = q x u(n x )for x = 1, . . . , 2R + 1 for bulk leftwards jumps, and r(n, n x,x+1 ) = p x u(n x )for x = 1, . . . , 2R + 1 (38) for bulk rightwards jumps, see figure 14. Note that equations (37) and (38) for x = 1 and x = 2R + 1, respectively, account for the particle removal at the boundaries. The generator of the dynamics can be written as for any real function f on Ω R . Therefore, particles hop on the lattice to the neighboring sites with rates qu(n x ) and pu(n x ), except for the lattice site at x = R + 1. Here, indeed, different hopping rates are assumed, namelyqu(n R+1 ) andpu(n R+1 ). As with the CA analyzed in section 4, also the ZRP considered here is an open system, because particles can be absorbed from the left and right reservoirs with rates γu(n 1 ) and βu(n 2R+1 ), respectively. Moreover, particles can also be released from the reservoirs into the channel with rates, respectively, α and δ.
In [49] it is proven that the considered ZRP can exhibit stationary uphill currents. More precisely, it can be shown that, for a specific choice of the parameters, the steady state is characterized by a current flowing from the reservoir with smaller density to the one with larger density. The proof of this statement will be reviewed in detail in the next sections since, unlike with most of the models discussed above, the ZRP admits, in some cases, an explicit analytical solution. In fact, once the stationary measure is identified (section 5.1), it is possible to compute the stationary current J R (section 5.2) that, in a weakly asymmetric case, can be explicitly shown to produce uphill transport (section 5.3).

Stationary measure for the ZRP
A probability measure µ R on Ω R is stationary for the considered ZRP if and only if ∑ for any function f. A sufficient condition is provided by the balance equation for any n ∈ Ω R . Consider the positive reals s 1 , . . . , s 2R+1 , called fugacities, and the product measure on the space Ω R defined as where u(k)! = 1 if k = 0 and u(k)! = u x (1) · · · u(k) if k ⩾ 1 and Applying equation (40) to the functions f(n) = n x for any x, it can be proven that ν is stationary for the considered ZRP if the fugacities s x satisfy the following set of equations [51,52]: for x = 1, . . . , R − 1 and x = R + 2, . . . , 2R. Solutions of these equations for a specific choice of the parameters p, q,p,q, α, β, γ, δ will be discussed in section 5.3.

Stationary current and density profile for the ZRP
The main quantities of interest, here, are the stationary density profile see [53] for details, and the stationary current for x = 1, . . . , 2R. The stationary current is typically computed, in the setup of ZRP, by considering a specific bond connecting two adjacent sites on the lattice. The quantity J R,x in (46), for instance, refers to the bond joining the sites x and x + 1, and represents the difference between the average number of particles hopping from the site x to the site x + 1, and the corresponding number of particles hopping from x + 1 to x. Since the stationary current does not depend on x, in the sequel we shall denote it simply by J R . Note that no specific form of the intensity function u has been invoked in the derivation of equation (46). Nevertheless, a choice for u is required in the computation of the density profile ρ x . Before focusing on some specific cases, it is worth recalling, preliminarily, the following general result [53]: which says that, for any x ∈ Λ, ρ x is an increasing function of the local fugacity. In [49] two relevant forms of the intensity function u are considered, corresponding to the independent particle and the simple exclusion-like ZRP models, in which the intensity function is given by u(k) = k and u(k) = 1 for k ⩾ 1, respectively. In these two cases one can show that Z ip x = exp{s x } and Z se x = 1/(1 − s x ) for s x < 1, respectively. Hence, using (45), one finds for the independent particle and the simple exclusion-like models, respectively.
We also remark that the hydrodynamic limit of the independent particle model considered in this section can be studied analytically. We refer the reader to [49] for details.

Conclusions
We reviewed different mathematical models leading to the appearance of uphill diffusions. A first important class of models is that of stochastic multi-species particle systems, discussed in section 2. For these models we outlined a statistical mechanical theory which clarifies the role of the coupling between the different species in producing the 'wrong' sign of the current for a specific component of the mixture. Our conclusions are in line with the results of earlier macroscopic approaches [12]. Something unexpected is instead observed in singlecomponent boundary driven systems, when the updating mechanism at the boundaries breaks the detailed balance condition and a phase transition occurs. We first reviewed the case of the non-equilibrium 2D Ising model considered in [30], in which positive spins migrate across a stationary 'bump' in the magnetization profile. The net magnetization flux goes uphill, namely from the reservoir with lower magnetization to the one with higher magnetization. We then turned to 1D models, in which the phase transition can be induced by a Kac interaction potential. We analyzed, in section 4, the main findings related to the cellular automaton studied in [27,28], for which a 'constitutive law', linking the stationary current to the reservoir magnetization, can be numerically determined. The investigation of stationary uphill currents in zero range processes -a class of interacting particle systems that can be solved exactly -was then reviewed in section 5. While many more particle models can be constructed to produce stationary uphill diffusions, we can already draw, from our investigation, some neat contour line clarifying the origin of such noteworthy transport phenomenon. In multi-species systems, the uphill migration of particles originates from the interplay between the density gradients of each species. In single-species systems, instead, the underlying mechanism is more subtle. In fact, it is related to the details of the updating mechanism at the boundaries and to the presence of a phase transition in the corresponding equilibrium model. Connections with systems displaying an absolute negative mobility are also matter of an ongoing investigation, and may unfold new perspectives in transport theory.

Data availability statement
No new data were created or analysed in this study.