Microscopic models for uphill diffusion

TL;DR

This paper investigates microscopic particle models exhibiting uphill diffusion phenomena under various external and attractive potentials, revealing complex behaviors including phase transition-like effects through numerical analysis.

Contribution

It introduces a particle system model that captures uphill diffusion and phase transition phenomena, providing numerical evidence and theoretical insights.

Findings

Uphill diffusion observed in the model with specific potentials.

Final stationary states align qualitatively with experimental data.

Phase transition-like behavior occurs under certain conditions.

Abstract

We study a system of particles which jump on the sites of the interval $[1,L]$ of $\mathbb Z$. The density at the boundaries is kept fixed to simulate the action of mass reservoirs. The evolution depends on two parameters $\lambda'\ge 0$ and $\lambda"\ge 0$ which are the strength of an external potential and respectively of an attractive potential among the particles. When $\lambda'=\lambda"= 0$ the system behaves diffusively and the density profile of the final stationary state is linear, Fick's law is satisfied. When $\lambda'> 0$ and $\lambda"= 0$ the system models the diffusion of carbon in the presence of silicon as in the Darken experiment: the final state of the system is in qualitative agreement with the experimental one and uphill diffusion is present at the weld. Finally if $\lambda'=0$ and $\lambda">0$ is suitably large, the system simulates a vapor-liquid phase transition and we have a surprising phenomenon. Namely when the densities in the reservoirs correspond respectively to metastable vapor and metastable liquid we find a final stationary current which goes uphill from the reservoir with smaller density (vapor) to that with larger density (liquid). Our results are mainly numerical, we have convincing theoretical explanations yet we miss a complete mathematical proof.

Figures

Equations

$${\sigma }_t(x)={\eta }_t(x)-1,m_{\pm }=2{\rho }_{\pm }-1$$

$${\sigma }^{t_0,T}(x)=\frac{1}{T}\sum _{t=t_0+1}^{t_0+T}{\sigma }_t(x)$$

$$j_{+}(t)={\eta }_t(L,1)-{\eta }_{t+1}(L,-1)$$

$$j_{-}(t)={\eta }_{t+1}(1,1)-{\eta }_t(1,-1)$$

$$j_{\pm }^{t_0,T}=\frac{1}{T}\sum _{t=t_0+1}^{t_0+T}{j}_{\pm }(t)$$

$$j_{\pm }^{t_0,T}\approx \frac{1}{2}\frac{m_{+}-m_{-}}{L},L=600,L=1200$$

$$U_{\gamma }(x)=\lambda \sum _y{J}_{\gamma }(y,x){\rho }_{\mathrm{si}}(y),J_{\gamma }(y,x)=\gamma J(\gamma |x-y|),\lambda >0$$

$$J(r)=(1-r)1_{0\le r\le 1}$$

$$f_{\gamma}(x)=-\gamma^{2}\lambda\Big{(}\sum_{y=x}^{x+\gamma^{-1}}\mathbf{1}_{y\leq\frac{L}{2}}-\sum_{x-\gamma^{-1}}^{x}\mathbf{1}_{y\leq\frac{L}{2}}\Big{)}$$

$$f_{\gamma }(x)={\gamma }^2\lambda 1_{|x-\frac{L}{2}|\le {\gamma }^{-1}}\times \{\begin{array}{ll}x-(\frac{L}{2}-{\gamma }^{-1}),& x\le \frac{L}{2}\\ (\frac{L}{2}+{\gamma }^{-1})-x,& x\ge \frac{L}{2}\end{array}$$

$${ϵ}_{x,\gamma }=\frac{1}{2}\beta f_{\gamma }(x)={\gamma }^2\frac{\beta \lambda }{2}{1}_{|x-\frac{L}{2}|\le {\gamma }^{-1}}\times \{\begin{array}{ll}x-(\frac{L}{2}-{\gamma }^{-1}),& x\le \frac{L}{2}\\ (\frac{L}{2}+{\gamma }^{-1})-x,& x\ge \frac{L}{2}\end{array}$$

$$\underset{\gamma \to 0}{\lim }\underset{{\gamma }^{2}t\to \tau ;\gamma x\to r}{\lim }E[{\sigma }_t(x)]=m(r,\tau )$$

$$\frac{\partial }{\partial t}m(r,t)=-\frac{\partial }{\partial r}j(r,t),m(0,t)=m_{-},m(\ell ,t)=m_{+}$$

$$j=\frac{1}{2}\Big{(}-\frac{\partial m}{\partial r}+\beta(1-m^{2})\{\lambda\int dx\rho_{\rm si}(x)\mathbf{1}_{|r-x|\leq 1}\big{(}\mathbf{1}_{x<r}-\mathbf{1}_{x\geq r}\big{)}\}\Big{)}$$

$$j=-\frac{1}{2}\frac{\partial m}{\partial r},\text{for ∣r−2ℓ∣>1 }$$

$$\underset{\ell \to \infty }{\lim }{m}^{(\ell )}(\frac{\ell }{2}+r)=A(r),|r|\le 1$$

$$A(r)=\{\begin{array}{ll}\tanh \{\beta (\lambda /2)[(1+r{)}^2-1]\},& r\in [-1,0]\\ \tanh \{\beta (\lambda /2)[1-(1-r{)}^2]\},& r\in [0,1]\end{array}$$

$$\underset{\ell \to \infty }{\lim }{m}^{(\ell )}(r\ell )=M(r)=\{\begin{array}{ll}-2\alpha r,& r\in [0,\frac{1}{2})\\ 2\alpha (1-r),& r\in (\frac{1}{2},1]\end{array},\alpha =A(1)=\tanh \{\beta \lambda /2\}$$

$$\underset{\ell \to \infty }{\lim }\ell j^{(\ell )}=\tanh \{\beta \lambda /2\}$$

$$\left|j_{\pm }^{t_0,T}-\frac{1}{L}\tanh \{\beta \lambda /2\}\right|\le 2\cdot 10^{-5}$$

$$F(m)=\int dr\Big{(}-\frac{S(m(r))}{\beta}+m(r)U(r)\Big{)}$$

$$S(m)=-\frac{1-m}{2}\log \frac{1-m}{2}-\frac{1+m}{2}\log \frac{1+m}{2}$$

$$U(r)=\lambda {\int }_0^{\ell }d{r}^{\prime }{\rho }_{\mathrm{si}}(r^{\prime })J(|r-r^{\prime }|)={\int }_0^{\ell /2}d{r}^{\prime }J(|r-r^{\prime }|)$$

$$j=-\chi \frac{d}{dr}\frac{\delta }{\delta m(r)}F(m),\chi =\frac{\beta }{2}(1-m^2)$$

$$\int dxU(|r-x|)\{\frac{d}{dx}{\rho }_{\mathrm{si}}(x)\}=U(|r-\frac{\ell }{2}|)$$

$$N_{+,x,\gamma }=\sum _{y=x+1}^{x+{\gamma }^{-1}}{\eta }^{(+)}(y),N_{-,x,\gamma }=\sum _{y=x-{\gamma }^{-1}}^{x-1}{\eta }^{(-)}(y),x\in [1,L]$$

$$\beta =2C=2.5$$

$$\displaystyle F(m)=\int\Big{(}-\frac{m^{2}}{2}-\frac{S(m)}{\beta}\Big{)}+\frac{1}{4}\int\int J(r,r^{\prime})[m(r)-m(r^{\prime})]^{2}$$

$$f_{\beta }(m):=-\frac{m^2}{2}-\frac{S(m)}{\beta }$$

$$m_{\beta }=\tanh \{\beta m_{\beta }\},m_{\beta }>0$$

$$j(r,t)=-\frac{1}{2}\big{\{}\frac{\partial m(r,t)}{\partial r}-\beta[1-m(r,t)^{2}]\int dxJ(|r-x|)\frac{\partial m(x,t)}{\partial x}\big{\}}$$

Full text

Keywords:— Stochastic cellular automata, Kac potential, Fourier law and phase transition, Uphill diffusion.

1 Introduction

Uphill diffusion is a phenomenon which appears when the current flows along the gradient in contrast with the Fick law which states that the current is proportional to minus the gradient. We are considering the case of mass diffusion so that the current is the mass flux and the gradient is the gradient of the mass density.

There are mainly two cases where uphill diffusion appears, the first one is when the system is a mixture of two or more components (or one component but several conserved quantities), the second one when the system undergoes a phase separation. The macroscopic explanation in the first case is that the current of say component 1 has a contribution $j_{1}$ proportional to minus the gradient of its density (in agreement with the Fick law) but also contributions coming from the gradients of the other components which may be larger than $j_{1}$ and with the opposite sign. The basic reference is an old paper by Darken, [4], where he gives experimental evidence of the phenomenon [15, 16], see also [10, 11] for an updated survey and Refs. [1, 6, 7, 8, 12, 14] for related recent results. We will present here a simple particle model which reproduces qualitatively the main features of the Darken experiments.

Uphill diffusion appears also in one component systems at phase transition, say a vapor-liquid transition. In such a case there is a density interval $(\rho^{\prime},\rho^{\prime\prime})$ so that if the density$\rho$ is $\leq\rho^{\prime}$ the system is in its vapor phase while if $\rho\geq\rho^{\prime\prime}$ then it is in its liquid phase. If we put a mass $\rho|\Lambda|$ of fluid in a region $\Lambda$ with $\rho\in(\rho^{\prime},\rho^{\prime\prime})$ we observe a mass flux which gives rise to a non homogeneous final density profile with vapor at density $\rho^{\prime}$ in a subregion $\Delta$ of $\Lambda$ while in the complement the phase is liquid at density $\rho^{\prime\prime}$. Thus during the phase separation mass has flown from the lower density in $\Delta$ to the larger density in $\Lambda\setminus\Delta$: it has gone uphill.

We call this phenomenon a “transient uphill diffusion” as in the end there is no current and we distinguish it from a “steady uphill diffusion”. The latter arises when a single component fluid in contact with a left and a right mass reservoir at density $\rho_{-}<\rho_{+}$ reaches a stationary state with positive current, namely mass flows from the reservoir at lower density to the one at larger density. We have observed this phenomenon in computer simulations of a particle system which models a vapor-liquid phase transition. The system is put in contact with mass reservoirs which keep fixed the densities at the boundaries. We have seen that if the left reservoir fixes a density $\rho_{-}$ in the metastable vapor phase while the right reservoir density is $\rho_{+}$ in the metastable liquid phase then the system reaches a stationary state where the current is positive, namely flows from left to right, i.e. mass goes from the reservoir at small density to that with larger density. Instead if the reservoirs densities are in the vapor and liquid stable phases (i.e. $\rho_{-}<\rho^{\prime}$, $\rho_{+}>\rho^{\prime\prime}$), then the current is negative and goes downhill.

The final stationary state when the current is positive could be either one where in most of the space the fluid is liquid with a small region close to the left boundary where there is a sharp transition from vapor to liquid or symmetrically one where in most of the space the fluid is vapor with a small region close to the right boundary where there is a sharp transition from liquid to vapor. In both cases the density profile is decreasing except in the transition region at either one of the boundaries, thus the steady current goes downhill in most of the space and uphill at the transition. The final stationary state is determined by the initial conditions and by random fluctuations.

We are not aware that such a “steady uphill diffusion” has been observed earlier and it certainly deserves to investigate whether analogous phenomena are present in more general systems.

In Section 2 we consider a particle model which describes normal diffusion; by adding a suitable potential we obtain in Section 3 a system which simulates the Darken experiment with carbon diffusing in the presence of silicon. In Section 4 we modify the model of Section 2 by adding an attractive force among the particles which gives rise to uphill diffusion in the presence of phase transition. In Appendix A we prove a theorem stated in Section 3.

2 A microscopic model for diffusion

We want to describe fluids where the evolution is diffusive and convection is negligible. These are the main features that our models should try to catch.

• •

Particles undergo very frequent collisions, their velocities change rapidly and erratically so that their motion looks diffusive, convection being absent.

• •

There is a strong repulsive force when particles are too close to each other which makes the density bounded.

• •

The system is confined in a cylindrical vessel with a horizontal axis. The two extremal faces are in contact with mass reservoirs which keep the mass density at the boundaries fixed. We suppose a planar symmetry in the vertical planes orthogonal to the axis of the cylinder.

We will first describe our model designed for implementation on a computer, then show some computer simulations and finally discuss how well it catches the above physical requests.

The model is one dimensional, space and time are discrete. Particles are confined in the interval $\{1,...,L\}$: $L$, a positive integer, is the spatial size of the system. Particles have only velocities equal to $+1$ and $-1$. Time is discrete: $t=0,1,2,\dots$ so that at each time step particles move from one site to the next one (right or left according to their velocity), we will say later what happens at the boundaries. There is an exclusion rule which prevents two particles with same velocity to stay on a same site, thus the local density is necessarily $\leq 2$. Before moving particles change randomly their velocities, however the exclusion rule prevents changes at sites where there are two particles (as they must have opposite velocities). The precise algorithm used to update the particles configurations is as follows.

Particles configurations are described by sequences $\eta=\{\eta(x,v),x\in[1,L],v\in\{-1,1\}\}$ with $\eta(x,v)\in\{0,1\}$ the occupation variable at the phase space point $(x,v)$. We denote by $\eta(x)=\eta(x,-1)+\eta(x,1)$ the total occupation at $x$ and add a suffix $t$ when the occupation variables are computed at time $t$, $t=0,1,2\dots$. The unit time step updating is obtained as the result of three successive operations starting from a configuration $\eta$ and ending with a configuration $\eta^{\prime\prime\prime}$, we denote by $\eta^{\prime}$ and $\eta^{\prime\prime}$ the configurations at the intermediate steps.

velocity flip. At all sites $x\in[1,L]$ where there is only one particle we update its velocity to become $+1$ with probability $\frac{1}{2}$ and $-1$ with same probability $\frac{1}{2}$ (such velocity flips are independent of each other). At all other sites the occupation numbers are left unchanged. We denote by $\eta^{\prime}$ the occupation numbers after the velocity flip updating. 2. 2.

advection. After deleting the particles at $(1,-1)$ and $(L,1)$ (if present) we let each one of the remaining particles move by one lattice step in the direction of their velocity. We denote by $\eta^{\prime\prime}$ the occupation numbers after this advection step. 3. 3.

boundaries updating. Let $\rho_{\pm}\in[0,1]$ and call $2\rho_{\pm}$ the density of the right, respectively left reservoir. Then with probability $\rho_{+}$ we put a particle at $(L,-1)$ and with probability $1-\rho_{+}$ we leave $(L,-1)$ empty. We do independently the same operations at $(1,1)$ but with $\rho_{-}$ instead of $\rho_{+}$. What we get is the final configuration $\eta^{\prime\prime\prime}$.

Let us next see how the model behaves. We have run several computer simulations, we report below some of them. It is convenient here and in the sequel to change variables writing

[TABLE]

thus $\sigma_{t}(x)\in\{-1,0,1\}$ and $m_{\pm}\in[-1,1]$. The above change of variables, which simplifies some formulas below, has also a physical meaning in terms of magnetic systems with $\sigma_{t}(x)$ a spin, we refer to [3] for details. We fix an initial datum where the variables $\eta(x,v)$ are independent and take values $0,1$ with same probability. For any choice of the initial datum we run the above algorithm for a time $t_{0}+T$ and measure for each $x\in[1,L]$ the time average

[TABLE]

The current at time $t$ from the system to the right reservoir is

[TABLE]

which counts as positive the particles which leave the system from the right and as negative those which enter from the right. Analogously the current at time $t$ from the the left reservoir to the system is

[TABLE]

We have also measured the averaged currents

[TABLE]

The simulations we report here are done with $L=600$ and $L=1200$, $t_{0}=10^{9}$, $T=10^{5}$, $m_{+}=0.5$, and $m_{-}=-0.5$ in both cases. In Fig. 1 we have plotted $m_{L}(r),r\in[L^{-1},1]$ by setting $m_{L}(r)=\sigma^{t_{0},T}(Lr)$ when $L=600$ and $L=1200$, the difference between the two profiles is negligible. The currents $j_{\pm}^{t_{0},T}$ are essentially equal to each other and

[TABLE]

In conclusion the simulations show that the system obeys the Fick law with constant diffusion coefficient equal to 1. For $L$ large the stationary profile $m_{L}(r)$ is linear connecting $m_{-}$ to $m_{+}$.

Let us finally discuss how well the particle model catches the physical requests stated at the beginning of the section. As shown from the simulations it indeed describes a diffusive fluid. The flip velocity updating however is not realistic, the collisions in a real fluid are not as simple and there are correlations between successive collisions, statistical independence that we assume here is the main issue in the derivation of the Boltzmann and other kinetic equations. The assumption that the speed is 1 is also unrealistic but it is the easiest way to achieve a description on the lattice which is more easily implementable on the computer. In the model the local density $\eta_{t}(x)$ is always $\leq 2$, in real systems a bound on the density comes from strong repulsive forces at short distances, as in Lennard-Jones, in our model it is simply achieved by forbidding velocity flips when two particles are on the same site. Clearly it is not the true reason yet it does the job. The restriction to one dimension reflects the assumption of planar symmetry on vertical planes and it is therefore quite acceptable. The action of the reservoirs updating is to keep the average density at $(L,-1)$ equal to $\rho_{+}$ and at $(1,1)$ equal to $\rho_{-}$. From the simulations we find that in average $\eta_{t}(x,v)\approx\eta_{t}(x,-v)$ so that the action of the reservoirs is to fix in the average the density at $1$ and $L$ equal to $2\rho_{\mp}$, i.e. the density of the reservoirs.

3 The Darken experiment

In [4] Darken reports of experiments which show uphill diffusion of carbon, we refer in particular to the case of Fig. 2 in [4] where carbon diffuses in a welded specimen where the silicon content is concentrated on the left of the weld (and negligible on the right). We refer to [4] for the details of the experiment. We model the carbon atoms using the particle model of the previous section but we need to modify the updating rules to take into account the presence of silicon. As stated in [4] we may neglect the diffusion of silicon so that we suppose that the stationary silicon density is equal to 1 (in appropriate density units) to the left of the weld and to 0 afterwards, i.e. $\rho_{\rm si}(x)=\mathbf{1}_{x\leq\frac{L}{2}}$. Carbon does not like to stay where the silicon is, hence the carbon will feel a positive potential $U_{\gamma}(x)$:

[TABLE]

where $\gamma^{-1}$ is a positive integer and

[TABLE]

Namely $U_{\gamma}(x)$ is a sum of the contributions $J_{\gamma}(y,x)\rho_{\rm si}(y)$ exerted from all silicon atoms. The interaction strength $J_{\gamma}(x,y)$ has range which scales as $\gamma^{-1}$, $\gamma^{-1}$ is a parameter of the model which on physical grounds should be much larger than the interatomic distance (which in our model is the distance between two successive sites and thus equal to 1) but also much smaller than the macroscopic size of the specimen, $L$ in our model.

Our choice of $J(r)$ is of course quite arbitrary, it has the advantage to give a simple formula for the force (defined as usual as minus the gradient of the potential):

[TABLE]

Thus $f_{\gamma}(x)$ is simply equal to $-\gamma^{2}\lambda$ times the difference between the number of silicon particles to the right and left of $x$ in a range $\gamma^{-1}$, hence $f_{\gamma}(x)$ is directed to the right and active only in a neighborhood of the weld:

[TABLE]

To represent the force in our model we argue that $\delta v$, the average velocity change per unit time due to the force, should be proportional to $\beta f_{\gamma}(x)$, $\beta$ the inverse temperature (as thermal fluctuations dampen the effect of the force). We can put this in our model where velocities are only $\pm 1$ by changing the probability of the velocity flips. We thus modify the updating rules of the previous section only in the flip velocity step: at sites where there is only one particle its velocity is updated to be $+1$ with probability $\frac{1}{2}+\epsilon_{x,\gamma}$ and equal to $-1$ with probability $\frac{1}{2}-\epsilon_{x,\gamma}$, thus the average velocity after the flip is $2\epsilon_{x,\gamma}$ which has the desired value $\beta f_{\gamma}(x)$ if

[TABLE]

As a consequence there is a bias to the right when close to the weld. (Recall that all the other updating rules are left unchanged).

Since we want to single out the effect of the force due to the silicon we take a homogeneous initial datum where as in Section 2 the variables $\eta(x,v)$ are independent and each one has an average equal to $1/2$. Then the average density, i.e. the average of $\eta(x)$, is equal to 1. Also the reservoirs have density $1$, i.e. $m_{\pm}=0$. Having defined the model and the initial datum we can now run the simulations. We take the size of the system $L$ equal to 600 or to $1200$, $\gamma^{-1}$ equal to 30 or to 60. We call $r=\gamma x$ the space measured in mesoscopic units and $\ell=\gamma L$ the size of the system in mesoscopic units.

The simulations show that the density profile is only weakly dependent on $\gamma$, hence they suggest that it may have a limit when $\gamma\to 0$. Indeed, under suitable assumptions on the initial distribution and a propagation of chaos property we can prove, proceeding as in [3], that for any $t\geq 0$ and $r\in(0,\ell)$

[TABLE]

where $m(r,t)$ satisfies the conservation law

[TABLE]

The current $j=j(r,t)$ is equal to

[TABLE]

The curly bracket is the continuum version of the sum in (3.3). In particular

[TABLE]

To check the validity of (3.7) we compare the profiles obtained in the simulations with the stationary solution of (3.7):

*Stationary problem: find a constant $j^{(\ell)}$ and a function $m^{(\ell)}(r)$, $r\in[0,\ell]$, so that $m^{(\ell)}(0)=m^{(\ell)}(\ell)=0$ and (3.8) is satisfied with $j^{(\ell)}$ and $m^{(\ell)}$. *

In Appendix A we prove:

Theorem 1**.**

The above “Stationary problem” has for each $\ell$ a unique solution $\{j^{(\ell)},m^{(\ell)}(x),x\in[0,\ell]\}$. Such a solution has the following properties: $j^{(\ell)}>0$, $|m^{(\ell)}(x)|<1$ for all $x$ and $\frac{d}{dx}m^{(\ell)}(x)=-2j^{(\ell)}$ for all $|x-\frac{\ell}{2}|>1$. Moreover

[TABLE]

while

[TABLE]

Finally:

[TABLE]

In 3 we compare the theoretical profile $A(r)$ and the average profile $\sigma^{t_{0},T}(x)$ which indicates that we are with good approximation close to the mesoscopic and macroscopic limits when $\gamma^{-1}=30$ and $L=1200$. The measured current is $j_{\pm}^{t_{0},T}\approx 7.2\cdot 10^{-4}$ so that by (3.11) and recalling that the mesoscopic current $j^{(\ell)}$ is related to the measured current by a factor $\gamma^{-1}$ we get

[TABLE]

which again shows that when $\gamma^{-1}=30$ and $L=\gamma^{-1}\ell=1200$ the system behaves with good approximation as in the macroscopic limit.

3.1 The mesoscopic theory

The evolution equation (3.7) which describes the dynamics of our model in the limit $\gamma\to 0$ has a nice physical interpretation. In fact let $F(m)$, $m\in L^{\infty}([0,\ell],[-1,1])$, be the free energy functional

[TABLE]

where $S(m)$ is the entropy and $U(r)$ the potential generated by the silicon, namely

[TABLE]

($J(|r|)$ being defined in (3.2)). Then $j$ in (3.8) is equal to

[TABLE]

$\chi$ being the mobility (of the carbon atoms). (3.16) is the usual constitutive law which states that the current is minus the mobility times the gradient of the chemical potential (which, according to thermodynamics, is the derivative of the free energy with respect to the density, recall that the density is $m(r)+1$). In our case where there is a contribution to the free energy coming from the force exerted by the silicon, the current $j$ is not only given by minus the gradient of the density but it has an additional contribution given by the second term in (3.8). The curly bracket in (3.8) has a clear physical meaning: the silicon atoms in the interval $dx$, i.e. $dx\rho_{\rm si}(x)$ generate a force field at $r$ which is $\lambda\mathbf{1}_{|r-x|\leq 1}$ if $x<r$ (hence positive) and a negative force $\lambda\mathbf{1}_{|r-x|\leq 1}$ if $x>r$. The same expression can be rewritten as

[TABLE]

Recalling that $U(|r-x|)=\lambda J(|r-x|)$ the left hand side of (3.17) is $\lambda$ times the weighted average of the density gradient of $\rho_{\rm si}(x)$. Thus the second term in (3.8) is $\lambda$ times the mobility times the averaged density gradient of $\rho_{\rm si}$. This is what expected from thermodynamics if the averaging weight was a delta function. Such an approximation would be valid if $m(r)$ were slowly varying at the edge, but this is not the case: when $\ell$ increases $m(r)$ becomes smoother but only away from the edge!

The above proves that there is a uphill diffusion with a mass flux from the left to the right reservoirs, despite they have the same density. The work done to ensure such a flow is provided by the force exerted on carbon by the silicon atoms. More generally we may take other densities for the reservoirs, for instance $m_{+}>0$ and $m_{-}=-m_{+}$. In such a case the analogue of Theorem 1 would give a positive current (and hence an uphill diffusion) if $m_{+}<\tanh\{{\beta\lambda/2}\}$, while the current would be downhill and hence negative if $m_{+}>\tanh\{{\beta\lambda/2}\}$. see Fig. 4

3.2 The macroscopic limit

The macroscopic limit is obtained by letting $\ell\to\infty$ while expressing the mesoscopic profiles in macroscopic units. By Theorem 1 the limit profile is $M(r)$, $r\in[0,1]$, which by (3.10) has a constant negative slope $-2\tanh\{\beta\lambda/2\}$ except at $r=\frac{1}{2}$ where it jumps from $-\tanh\{\beta\lambda/2\}$ to $\tanh\{\beta\lambda/2\}$. Fick’s law is satisfied and the differential equation in the macroscopic limit splits in two equations: one in $(0,\frac{1}{2})$ and the other in $(\frac{1}{2},1)$. The boundary conditions are $m(0)=0$ and $m(\frac{1}{2})=-\tanh\{\beta\lambda/2\}$ for the first equation and $m(\frac{1}{2})=\tanh\{\beta\lambda/2\}$ and $m(1)=0$ for the second one. Thus in the macroscopic scaling the action of the force (due to silicon) is represented by boundary conditions at the discontinuity, this is a boundary layer problem which requires an analysis of the boundary layer in the stretched mesoscopic variables, as done in Theorem 1.

4 “Steady uphill diffusion” and phase transitions

In this section we discuss uphill diffusion due to phase separation referring to results obtained in [2], [3] and [5]. As already mentioned in the introduction there are several surprising facts which are not entirely understood.

The phase transitions that we consider are of liquid-vapor type and the context is the one proposed by van der Waals. As explained by van der Waals it is the presence of long range attractive forces which is responsible for the phase transition. The way to implement his ideas in particle systems was first proposed by Kac with the introduction of Kac potentials, these are potentials which scale with a parameter $\gamma>0$, the range scaling as $\gamma^{-1}$ and the strength of the potential as $\gamma^{d}$ ($d$ the dimension of the space), so that the total interaction of a point with the others stays finite as $\gamma\to 0$. In [9] and [13] it was shown that the Gibbsian statistical mechanics with Kac potentials reproduces the van der Waals theory in the limit $\gamma\to 0$.

To implement all that we simply go back to the basic model of Section 2 and add an interaction of Kac type among particles. This is just what we did in Section 3 but the force is now given by the same diffusing particles of the system and not by an external force (which in the previous section was exerted by the silicon atoms). In suitable units we suppose that the strength of the force (which in Section 3 was denoted by $\lambda$) is now equal to 1. The model we obtain is that considered by the same authors in [2] and [3]. Referring to the system of Section 3 we only have to modify the velocity flip updating which is now as follows:

Velocity flip. At all sites $x\in[1,L]$ where there is only one particle we update its velocity to become $+1$ with probability $\frac{1}{2}+\epsilon_{x,\gamma}$ and $-1$ with probability $\frac{1}{2}-\epsilon_{x,\gamma}$, $\epsilon_{x,\gamma}=C\gamma^{2}[N_{+,x,\gamma}-N_{-,x,\gamma}]$; at all other sites the occupation numbers are left unchanged. We have set

[TABLE]

where $\eta^{(+)}(y)=\eta(y)$ if $y\in[1,L]$ and $\eta^{(+)}(y)=2\rho_{+}$ if $y>L$; similarly $\eta^{(-)}(y)=\eta(y)$ if $y\in[1,L]$ and $\eta^{(-)}(y)=2\rho_{-}$ if $y<1$, recall that $2\rho_{\pm}$ is the density of the right, respectively left reservoir. We choose $C=1.25$ and $\gamma^{-1}=30$ so that the definition is well posed because $(2\gamma^{-1})C\gamma^{2}=2.5/30<\frac{1}{2}$, $(2\gamma^{-1})$ being an upper bound for $|N_{+,x,\gamma}-N_{-,x,\gamma}|$.

By (3.5) and recalling that $\lambda=1$ the above choice implies that particles are in contact with an environment which keeps the inverse temperature $\beta$ equal to

[TABLE]

In Section 5 of [3] it is shown that in the limit $\gamma\to 0$ considered in the previous section the evolution is ruled again by the conservation law (3.7) with the current $j(r,t)$ given by (3.16) where the free energy functional $F(m)$ is now given by

[TABLE]

with $m(r)=m_{\pm}$ if $r\geq\ell$ and respectively $r\leq 0$. The first term on the right hand side, namely

[TABLE]

is the van der Waals mean field free energy, which is a convex function for $\beta\leq 1$ while for $\beta>1$ becomes a double well with minima at $\pm m_{\beta}$ where

[TABLE]

In our model $\beta=2.5$ hence we are in the phase transition regime. The values $|m|\geq m_{\beta}$ define the stable phases, the interval $|m|<m_{\beta}$ is the spinodal region. Inside the spinodal region the set $|m|<m^{*}$, $m^{*}>0:\beta(1-(m^{*})^{2})=1$, is unstable while the region $m^{*}<|m|<m_{\beta}$ is metastable.

The current $j(r,t)$ given by (3.16) with $F(m)$ as in (4.3) is:

[TABLE]

If we suppose that $\frac{\partial m(x,t)}{\partial x}\approx\frac{\partial m(r,t)}{\partial r}$ in the support of $J(|r-x|)$ then

[TABLE]

Thus when $|m(r,t)|<m^{*}$, i.e. in the unstable region, the current has the same sign of the gradient of $m$ and the diffusion coefficient is negative. This is at the basis of the macroscopic explanation of the uphill diffusion.

In our context the approximation $\frac{\partial m(x,t)}{\partial x}\approx\frac{\partial m(r,t)}{\partial r}$ in the support of $J(|r-x|)$ is shaky because in the spatial region where $m$ is unstable its values vary significantly. In our simulations we set $m_{+}=-m_{-}>0$ and decrease $m_{+}$ starting from its maximal value 1.

In Fig. 5 we see a negative (downhill) current till $m_{+}>m_{\beta}$ while in Fig. 6 it becomes positive as $m_{+}<m_{\beta}$, in this latter case the current goes uphill from the reservoir with smaller to the one with larger density. The stationary profile has also a significant change, when $m_{+}>m_{\beta}$ it is smooth away from a small neighborhood of the middle point, where instead has a sharp jump going from $\approx-m_{\beta}$ to $m_{\beta}$. Instead when $m_{+}<m_{\beta}$ and metastable (i.e. $m_{+}>m^{*}$) the jump moves to one of the endpoints, in the simulation presented in Fig. 6 it goes to the left boundary where it jumps from $m_{-}$ to a value $b(m_{-})$ which is larger than $m_{+}$, it then decreases smoothly toward the value $m_{+}$ reached at the right boundary.

When $m_{+}<m^{*}$ the current is still positive (i.e. uphill) but the profile has a more complex structure, we refer to [3] for details. In [3] we give some theoretical explanation of these phenomena but a complete theory with mathematical proofs is still missing.

The canonical system in statistical mechanics to study phase transitions is the Ising model in $d\geq 2$ dimensions with ferromagnetic nearest neighbor interactions. There are preliminary results with computer simulations in the $d=2$ case (obtained by Colangeli, Giardinà, Giberti, Vernia) which show again uphill diffusion in essential agreement with what described above.

Acknowledgments:

M.C. acknowledges useful discussions with Ellen Moons and Andrea Muntean (Karlstad University).

Appendix A The stationary problem

In this section we prove Theorem 1. We will first prove in Corollary 3 below the existence of stationary solutions, namely we will prove that there are a constant $j$ and a function $m(x)$, $x\in[0,\ell]$, such that

[TABLE]

where

[TABLE]

In particular

[TABLE]

Observe that $h(x):=g(\frac{\ell}{2}+x)$, $|x|\leq 1$, is equal to

[TABLE]

and therefore independent of $\ell$.

We start by proving the following lemma:

Lemma 2**.**

Denote by $m_{(j)}(x)$, $x\geq 0$, $j\in\mathbb{R}$, the solution of

[TABLE]

Then

[TABLE]

Proof. Take $j^{\prime}>j$. Then by (A.3) $m_{(j)}(x)=-2jx>m_{(j^{\prime})}(x)=-2j^{\prime}x$ for $x\in(0,\frac{\ell}{2}-1]$. Suppose by contradiction that there is $y>\frac{\ell}{2}-1$ such that $m_{(j)}(x)>m_{(j^{\prime})}(x)$ for $x<y$ and $m_{(j)}(y)=m_{(j^{\prime})}(y)$. Then $(1-m_{(j)}(y)^{2})g(y)=(1-m_{(j^{\prime})}(y)^{2})g(y)$ and therefore

[TABLE]

which contradicts the inequality $m_{(j)}(x)>m_{(j^{\prime})}(x)$ valid for $x<y$, hence the first statement in (A.6).

To prove the second statement in (A.6) we first consider $j<0$. In such a case $\displaystyle{\frac{d}{dx}m_{(j)}(x)\geq-2j>0}$ for all $x$, hence for any $x>0$ $\lim_{j\to-\infty}m_{(j)}(x)=\infty$. When $j>0$ we define

[TABLE]

In $\mathcal{X}$ we have $\displaystyle{\frac{d}{dx}m_{(j)}(x)\leq-2j<0}$. As a consequence if $y\in\mathcal{X}$ then $x\in\mathcal{X}$ for all $x\geq y$. For $j$ large enough $\frac{\ell}{2}-1\in\mathcal{X}$ and therefore

[TABLE]

hence $\lim_{j\to\infty}m_{(j)}(\ell)=-\infty$. ∎

Corollary 3**.**

There is a unique solution $\{j^{(\ell)},m^{(\ell)}(x),x\in[0,\ell]\}$ of (A.1). Furthermore $j^{(\ell)}>0$ and $|m^{(\ell)}(x)|<1$ for all $x\in[0,\ell]$.

Proof. $m_{(j)}(\ell)$ is a continuous, strictly decreasing function of $j$ which converges to $\pm\infty$ as $j\to\mp\infty$, hence there is a unique $j^{(\ell)}$ such that $m_{(j^{(\ell)})}(\ell)=0$ and $m^{(\ell)}(x):=m_{(j^{(}\ell))}(x)$, hence it is the unique solution of (A.1).

Call for notational simplicity $m(x)$ and $j$ the unique solution of (A.1). Suppose by contradiction that $j\leq 0$, then, by (A.3), $dm(x)/dx\geq 0$ for all $x$ and $dm(x)/dx>0$ for some $x\in(\frac{\ell}{2}-1,\frac{\ell}{2}+1)$, which yields $m(\ell)>0$, while $m(\ell)=0$ by (A.1), thus $j>0$.

Suppose again by contradiction that there is $y$ such that $m(x)<1$ for all $x<y$ and $m(y)=1$. By (A.3), $dm(y)/dy=-2j<0$ hence the contradiction because it would mean that $m(x)>1$ for $x<y$ and $y-x$ small enough. Suppose again by contradiction that there is $y<\ell$ such that $m(y)=-1$, then $y\in\mathcal{X}$, see (A.7), and as argued in the proof of Lemma 2, this implies $\ell\in\mathcal{X}$, while $m(\ell)=0$ by (A.1). ∎

Lemma 4**.**

Let $\{j^{(\ell)},m^{(\ell)}\}$ be as in Corollary 3, then

[TABLE]

Proof. By (A.3) $m^{(\ell)}(\frac{\ell}{2}-1)=-2j^{(\ell)}(\frac{\ell}{2}-1)$ and $m^{(\ell)}(\frac{\ell}{2}+1)=j^{(\ell)}(\ell-[\frac{\ell}{2}+1])$ hence (A.8) having used that $|m^{(\ell)}|<1$. ∎

Let $a(x|\alpha)$, $|x|\leq 1$, $\alpha\in(-1,0)$, be the solution of

[TABLE]

where $h(x)$ has been defined in (A.4). Explicitly:

[TABLE]

Lemma 5**.**

Let $a^{(\ell)}(x)=a(x|\alpha)$ with $\alpha=m^{(\ell)}(\frac{\ell}{2}-1)$. Then

[TABLE]

Proof. Since $|m^{(\ell)}|<1$ and, by (A.10), $|a^{(\ell)}|<1$

[TABLE]

hence (A.11). ∎

Proof of (3.9). Let $\ell_{n}$ be any sequence such that $\ell_{n}\to\infty$ and such that $m^{(\ell_{n})}(\frac{\ell_{n}}{2}-1)$ has a limit, call it $-\alpha^{*}$. Then by (A.8) $\alpha^{*}>0$ and $m^{(\ell_{n})}(\frac{\ell_{n}}{2}+1)\to\alpha^{*}$. By Lemma 5 and (A.10)

[TABLE]

and $a(1|-\alpha^{*})=\alpha^{*}$. By (A.10)

[TABLE]

and this implies $\alpha^{*}=\tanh\{\beta\lambda\}$. ∎

Proof of (3.11). From (A.8) and (3.9) we have

[TABLE]

∎

Proof of (3.10). Let $r\in[0,\frac{1}{2})$ then

[TABLE]

The analogous statement holds for $r>1/2$ and (3.10) is proved. ∎