Emergence of stationary uphill currents in 2D Ising models: the role of reservoirs and boundary conditions

TL;DR

This study demonstrates that in a 2D Ising model with non-equilibrium reservoirs, stationary uphill currents emerge below a critical temperature, especially under non-local equilibrium conditions, revealing novel transport phenomena.

Contribution

The paper provides numerical evidence of uphill currents in a 2D Ising model with non-equilibrium reservoirs, highlighting the impact of boundary conditions and temperature on these phenomena.

Findings

Uphill currents occur below the critical temperature in the ISF model.

Local equilibrium violation leads to unexpected transport phenomena.

Different boundary conditions significantly affect the emergence of uphill currents.

Abstract

We investigate the dynamics of a 2D Ising model on a square lattice with conservative Kawasaki dynamics in the bulk, coupled with two external reservoirs that pull the dynamics out of equilibrium. Two different mechanisms for the action of the reservoirs are considered. In the first, called ISF, the condition of local equilibrium between reservoir and the lattice is not satisfied. The second mechanism, called DB, implements a detailed balance condition, thus satisfying the local equilibrium property. We provide numerical evidence that, for a suitable choice of the temperature (i.e. below the critical temperature of the equilibrium 2D Ising model) and the reservoir magnetizations, in the long time limit the ISF model undergoes a ferromagnetic phase transition and gives rise to stationary uphill currents, namely positive spins diffuse from the reservoir with lower magnetization to the reservoir with higher magnetization. The same phenomenon does not occur for DB dynamics with properly chosen boundary conditions. Our analysis extends the results reported in Colangeli et al. (Phys. Rev. E, 97:030103(R), 2018), shedding also light on the effect of temperature and the role of different boundary conditions for this model. These issues may be relevant in a variety of situations (e.g. mesoscopic systems) in which the violation of the local equilibrium condition produces unexpected phenomena that seem to contradict the standard laws of transport.

Tables

Table 1. Table 1: The sorted list of ( m + , β ​ h + ) subscript 𝑚 𝛽 subscript ℎ (m_{+},\beta h_{+}) pairs contains all values ’selected’ in this manuscript.

$m_{+}$	$\beta h_{+}$	$m_{+}$	$\beta h_{+}$	$m_{+}$	$\beta h_{+}$	$m_{+}$	$\beta h_{+}$
0.000000	0.00	0.716298	0.90	0.874053	1.35	0.998778	3.70
0.148885	0.15	0.781806	1.05	0.970452	2.10	0.998894	3.75
0.291313	0.30	0.817754	1.15	0.978026	2.25	0.999329	4.00
0.537050	0.60	0.833655	1.20	0.998650	3.65	0.999955	5.35
0.604368	0.70	0.848284	1.25	0.999273	3.96	1.000000	$\infty$

Table 2. Table 2: Notation for magnetizations considered in this work. Whenever a symbol h ℎ h is used, it is related to the corresponding magnetization by m = tanh ⁡ ( β ​ h ) 𝑚 𝛽 ℎ m=\tanh(\beta h) , where β 𝛽 \beta denotes inverse temperature, c.f. Tab.  1 .

Symbol	Definition	Description
$m_{+}$	Parameter	Magnetization of reservoir ${ℛ}_{+}$
$m_{-}$	$=-m_{+}$	Magnetization of reservoir ${ℛ}_{-}$
$m_{\text{res}}$	res $\in \{-1,+1\}$	stands for $m_{+}$ or $m_{-}$
$m_{\beta }$	Eq. (2)	Onsager’s equilibrium magnetization
${\overline{m}}_x$	Eq. (17)	Mean magnetization at position $x\in \{1,..,L\}$
${\overline{m}}_L$	$={\overline{m}}_{x=L}$	Mean magnetization at the right boundary
${\overline{m}}_L^{\text{X}}$		${\overline{m}}_L$ obtained using model X
$m_{\text{eq}}$	$={\overline{m}}_L$	Mean absolute magnetization at the right boundary for $ℰ$-model
$m_{\text{eq}}^{\text{X}}$	$={\overline{m}}_L^{\text{X}}$	$m_{\text{eq}}$ obtained using model X $\in ℰ$
$m_{\text{crit}}$		Largest $m_{+}$ value for which $J=0$
$m_{\text{crit}}^{\text{X}}$		$m_{\text{crit}}$ obtained using model X $\in 𝒩$
$m_{\text{bump}}$		maximum ${\overline{m}}_x$ value close to boundary of the bump profile
$m_o$	$={\overline{m}}_L^{\text{DB-L4}}$	for the case of $m_{+}=0$
$\overline{m}$	$=L^{-1}{\sum }_{x=1}^L{\overline{m}}_x$	Mean bulk magnetization
$m\left(r\right)$		Magnetization at reduced coordinate $r\in [0,1]$
$m_u$		stability region for bump $m\in [m_u,m_{\text{crit}}]$
$m_{+}^{\text{eff}}$	Eq. (18)	effective magnetization, ${\overline{m}}_L^{\text{ISF-L4}}(\beta ,m_{+}^{\text{eff}})={\overline{m}}_L^{\text{DB-L4}}(\beta ,m_{+})$

Equations

$${\beta }_c=\frac{\ln (1+\sqrt{2})}{2}\approx 0.440686.$$

$$m_{\beta }={\left[1-\frac{1}{{\sinh }^4\left(2\beta \right)}\right]}^{1/8}.$$

$$H(\sigma |{\sigma }^o)=-\frac{1}{2}\sum _{\underset{i,j\in \Lambda }{⟨i,j⟩}}{\sigma }_i{\sigma }_j-\sum _{\underset{i\in \Lambda ,j\in {\Lambda }^o}{⟨i,j⟩}}{\sigma }_i{\sigma }_j^o,$$

$$m(\sigma )=\frac{1}{L^2}\sum _{i\in \Lambda }{\sigma }_i.$$

$$\beta h_{\text{res}}={\tanh }^{-1}(m_{\text{res}})=\frac{1}{2}\ln \left(\frac{1+m_{\text{res}}}{1-m_{\text{res}}}\right),{\text{m}}_{\text{res}}\in \{m_{-},m_{+}\}$$

$$c_{i,j}^{\text{bulk}}=\text{min}(1,e^{-\beta \Delta H}),\Delta H=H({\sigma }^{ij}|{\sigma }^o)-H(\sigma |{\sigma }^o),$$

$$[L4]{\sigma }_{(\bar{x},y)}^o={\sigma }_{(x_b,y-L/4)},\forall \bar{x}\in \{x_b\pm 1\},\forall y\in \{1,..,L\}.$$

$$[\mathbf{FB}]{\sigma }_i^o=0,\forall i\in {\Lambda }^o.$$

$$[\mathbf{SB}]\text{prob}({\sigma }_i^o)=\frac{1+{\sigma }_i^o{m}_{\text{res}}}{2},\forall i\in {\Lambda }^o.$$

$$c^{\text{ISF}}({\sigma }_i^{\prime })=\frac{1+{\sigma }_i^{\prime }{m}_{\text{res}}}{2}$$

$$c^{\text{NSF}}=0.$$

$$c^{\text{DB}}({\sigma }_i^{\prime })=\text{min}\left[1,\left(\frac{1+m_{\text{res}}{\sigma }_i^{\prime }}{1+m_{\text{res}}{\sigma }_i}\right)e^{-\beta \Delta H}\right],\Delta H=H({\sigma }^{\prime }|{\sigma }^o)-H(\sigma |{\sigma }^o),$$

$$\frac{1+m_{\text{res}}{\sigma }_i^{\prime }}{1+m_{\text{res}}{\sigma }_i}=e^{\beta h_{\text{res}}({\sigma }_i^{\prime }-{\sigma }_i)},\forall {\sigma }_i,{\sigma }_i^{\prime }\in \{-1,+1\}$$

$$c^{\text{DB}}({\sigma }_i^{\prime })=\text{min}\left(1,e^{-\beta \Delta \mathcal{H}}\right),\Delta \mathcal{H}=\mathcal{H}({\sigma }^{\prime }|{\sigma }^o)-\mathcal{H}(\sigma |{\sigma }^o)$$

$$\mathcal{H}(\sigma |{\sigma }^o)=H(\sigma |{\sigma }^o)-\sum _{i\in (x_b,y)}{h}_{\text{res}}{\sigma }_i$$

$$m_x(t)=\frac{1}{L}\sum _{y=1}^L{\sigma }_{(x,y)}(t),\forall x\in \{1,..,L\}$$

$${\overline{m}}_x=\frac{1}{T}\sum _{t=1}^T{m}_x(t),\forall x\in \{1,..,L\}$$

$${\overline{m}}_L^{\text{ISF-L4}}(\beta ,m_{+}^{\text{eff}})={\overline{m}}_L^{\text{DB-L4}}(\beta ,m_{+})$$

$${\overline{m}}_x=\frac{2x-1-L}{L+1}{m}_{+},{\overline{m}}_L=-{\overline{m}}_1=\frac{L-1}{L+1}{m}_{+}(\beta =0).$$

Full text

Fick’s law; Uphill diffusion; 2D Ising model; Nonequilibrium steady states; Monte Carlo methods.

1. I Introduction

2. II The models

3. II.1 Kawasaki dynamics in the bulk

4. II.2 Boundary conditions

5. II.3 Spin-updating mechanisms

6. II.4 Observables

7. II.5 Implementation details

8. III Results and discussion

9. III.1 Behavior of $J$ in external magnetization-temperature space

10. III.2 Critical value $m_{\textrm{crit}}$, NSF equilibrium magnetization $m_{\textrm{eq}}$ and DB equilibrium magnetization $m_{o}$

11. III.3 Average magnetization at the boundaries

12. III.4 High temperature regime

13. III.5 Low temperature regime

14. IV Conclusions

I Introduction

The presence of uphill diffusion in particle or spin models is a remarkable problem of statistical mechanics, that was first envisaged by Darken in his pioneering work dating back to 1949 Darken49 . While Fick’s law of diffusion dictates that particles diffuse against the concentration gradient (downhill diffusion), in presence of a multi-component system or an external field particles may be found migrating up the gradient CDMP17bis ; CC17 , giving thus rise to the so-called uphill currents, confirmed both experimentally and theoretically in numerous studies Lesher1994 ; Frink2000 ; Ferri2006 ; Kuhl2007 ; Yu2007 ; Sundman2009 ; Vielzeuf2011 ; Boudin2012 ; Rougier2013 ; Laeuerer2015 . Remarkably, in CDMP16 ; CDMP17 some numerical and theoretical evidence supported the conclusion that uphill diffusion may also occur in single-component systems in presence of a phase transition. In these works, the authors consider a 1D lattice gas model coupled to external particle reservoirs, in which particles that hop on the lattice are subject to an exclusion principle and are equipped with a long range, Kac-like interaction Presutti09 , which gives rise, at sufficiently low temperatures, to a phase transition. By adopting spin variables, the authors show that when the absolute value of the magnetization at the boundaries is larger than the equilibrium mean field magnetization, a sharp interface (called the “instanton” therein) appears at the center of the lattice accompanied by a downhill current. The novelty comes when the absolute magnetization at the boundary is lowered below the equilibrium mean field value: the interface, referred to as the “bump”, is then observed moving towards one of the boundaries, and the current changes sign, namely positive spins move from the reservoir with higher magnetization towards the one with lower magnetization. One might consider, therefore, an experimental set-up in which two finite reservoirs, equipped with metastable values of magnetization, are connected by two channels: one crossed by an uphill current and another one (in which the long-range Kac interaction is absent) featuring a Fickian current. The result would indeed be a closed circuit in which, in presence of a phase transition, current spontaneously flows across the ring. This phenomenon was indeed observed and reproduced in CDMP17 by performing extensive Monte Carlo simulations of the model (note that no violation of the thermodynamic principles occurs therein, as the total energy of the system is not a conserved quantity). A step forward was then taken in CGGV18 , in which a similar phenomenology was also observed in a 2D Ising model on a square lattice coupled to two external magnetization reservoirs attached at the right and left boundaries (whose magnetizations are equal, respectively, to $m_{+}>0$ and $m_{-}=-m_{+}$). The presence of a stationary uphill current is induced by the reservoirs, whose updating mechanism at the boundaries breaks the condition of local detailed balance, cf. also (ELS1990, , Eq. 1.6). In particular, the authors show that, in presence of a ferromagnetic phase transition, a certain critical value $m_{\textrm{crit}}$ marks the transition between two different regimes (referred to, in CGGV18 , as the stable and metastable one), and characterized, respectively by downhill and uphill currents, when $m_{+}$ is, respectively, larger or smaller than $m_{\textrm{crit}}$.

We recall that in the absence of reservoirs and in the infinite volume limit, the equilibrium 2D Ising model undergoes a phase transition at the inverse critical temperature Onsager1944

[TABLE]

For inverse temperatures $\beta>\beta_{c}$ the 2D Ising model (with vanishing external magnetic field) exhibits a spontaneous magnetization given by Yang1952

[TABLE]

In CGGV18 it is claimed that the critical value $m_{\textrm{crit}}$, evaluated at $\beta=1$, can also be estimated by measuring the magnetization value $m_{\textrm{eq}}$ (that approaches, in the large volume limit, the critical value $m_{\beta}$ in Eq. (2) Onsager1944 ; Yang1952 ) evaluated at the rightmost column of an Ising model in equilibrium conditions (i.e. in the absence of reservoirs and external magnetic fields) and characterized by a conservative dynamics. The claim above thus indicates the possibility that a “nonequilibrium” observable, such as $m_{\textrm{crit}}$, characterizing the dynamics of a boundary-driven Ising model, may be estimated from the analysis of an equilibrium Ising model. We shall tackle carefully this question in Section III and will show that the two quantities $m_{\textrm{crit}}$ and $m_{\textrm{eq}}$, evaluated for some fixed $L$ and $\beta$, are generally different from one another, their deviation becoming negligible only for large values of $\beta$ and by choosing suitable boundary conditions (b.c.s). At $\beta=1$ and with the b.c. considered in CGGV18 , the two observables are, indeed, almost coinciding.

In this work, we continue along the path traced in CGGV18 and aim to investigate in more detail the effect of different b.c.s for the considered 2D Ising model, cf. also Jan1983 , and unveil novel regimes other than the stable and metastable ones, by dropping further below the absolute value of the magnetization of the reservoirs. In particular, we shall discuss the behavior of the equation of state, namely the relation between the stationary current and the value of the magnetization at the boundaries, thus extending the preliminary results obtained in CGGV18 , cf. Fig. 2 therein.

The manuscript is organized as follows. The geometry and the bulk dynamics of our model are introduced in Section II, together with the boundary conditions in Section II.2, and details of the spin-updating mechanisms (due to the coupling with the reservoirs) in Section II.3. The role of the chosen spin-updating mechanism, and the comparison with a different mechanism based on detailed balance, are clarified in the same section. Definitions of observables and implementation details are provided in Sections II.4 and II.5. The results of our simulations are presented and discussed in Section III. Section IV is devoted to conclusions.

II The models

We consider non-equilibrium Monte Carlo (MC) dynamics of the nearest-neighbor ferromagnetic Ising model on a finite square 2D binary $L\times L$ lattice $\Lambda$ of linear size $L$, that is periodic only in one dimension (the vertical or synonymously, the $y$–direction). We shall denote by $\sigma_{i}\in\{-1,+1\}$ the “bulk” spin state at the lattice coordinate $i=(x,y)\in\Lambda$.

The system exhibits conservative exchange dynamics in the bulk (Section II.1), and is coupled to two infinite and interaction-free magnetization reservoirs on the horizontal direction (the $x$–direction, see Fig. 1 for a description of the set-up), denoted as $\mathcal{R}_{+}$ and $\mathcal{R}_{-}$, located, respectively, at the right and the left ends of the lattice. To this end, we define a bounded domain $\Lambda^{o}$, constituted by two vertical stripes, each made of $L$ lattice sites, located at $x=\bar{x}\in\{0,L+1\}$. We shall then denote by $\sigma^{o}_{i}$, $i=(\bar{x},y)\in\Lambda^{o}$ the ghost spins, picking up values in a set which depends on the chosen boundary model, as described in Section II.2. We denote by $N_{b}=2L^{2}+L$ the total number of bonds in the system, consisting of $L^{2}+L(L-1)$ bulk bonds between adjacent spins in $\Lambda$ and $2L$ horizontal bonds connecting the spins at the boundaries of $\Lambda$ (hereafter called boundary spins, located at $x=x_{b}\in\{1,L\}$) with the $2L$ so-called ghost spins in $\Lambda^{o}$. Note that the Ising model is equivalent to a lattice gas model via the standard mapping between $L^{2}$ spin variables $\sigma_{i}$ and occupation variables $\eta_{i}=(1+\sigma_{i})/2\in\{0,1\}$ with $\eta_{i}=1$ (resp. $\eta_{i}=0$) denoting the presence (resp. absence) of a particle.

Calling $\sigma=\{\sigma_{i}\}$, $i\in\Lambda$ and $\sigma^{o}=\{\sigma^{o}_{i}\}$, $i\in\Lambda^{o}$, respectively, the spin configuration and the boundary condition, we define the Hamiltonian

[TABLE]

where $\langle i,j\rangle$ denotes a nearest neighbor pair, and the magnetization of spin configuration $\sigma$ is given by

[TABLE]

Both the properties of the ghost spins at $\bar{x}$ (Section II.2) and the dynamics of the boundary spins at $x_{b}$ (Section II.3) are eventually affected by the two reservoirs $\mathcal{R}_{+}$ and $\mathcal{R}_{-}$. These reservoirs are fully characterized by their constant magnetizations $m_{\textrm{res}}$, respectively $m_{+}\in[0,1]$ and $m_{-}=-m_{+}$, respectively. In the sequel, when convenient, we will use

[TABLE]

in which $h_{\textrm{res}}$ can be regarded as the strength of a fictitious, dimensionless magnetic field, taking values $h_{+}\geq 0$ and $h_{-}=-h_{+}$ at the left and the right reservoirs, respectively. It will become soon clear, in Section II.3, the reason of this notation: when the mechanism of interaction with the reservoirs is based on detailed balance, $\tanh^{-1}(m_{+})$ attains indeed the structure of a dimensionless magnetic field multiplied by inverse temperature $\beta$. Table 1 provides the translation from $m_{+}$ values to $\beta h_{+}$ values considered in most plots of the following sections.

II.1 Kawasaki dynamics in the bulk

We consider the Ising system with reservoirs at inverse temperature $\beta$ and let the spins evolve following a continuous-time stochastic dynamics with two contributions: a conservative exchange dynamics in the bulk and and spin flip mechanisms at the opposing vertical boundaries. More precisely, in the bulk, the spins follow a Kawasaki dynamics, i.e. the two spins of a selected bond $\langle i,j\rangle$, with $i,j\in\Lambda$, exchange (“bond flip”) their values with rate

[TABLE]

where $\sigma^{ij}$ denotes the configuration obtained from $\sigma$ by exchanging the spins at sites $i$ and $j$. In the sequel we shall investigate the stationary dynamics corresponding to two different spin-updating mechanisms at the two vertical boundaries due the interaction with the reservoirs $\mathcal{R}_{-}$ and $\mathcal{R}_{+}$, and three different b.c.s.

II.2 Boundary conditions

Periodic b.c.s hold along the vertical direction of the considered model. The focus, hereafter, will be the investigation of different b.c.s imposed along the horizontal direction. These will indeed be found to affect the stationary magnetization profiles $\overline{m}_{x}$ (as function of grid coordinate $x$) as well as the stationary, $x$-independent current $J$, which is defined as the time-averaged rate of change of a boundary spin value.

The “L4” b.c., inspired by BP2003 and recently implemented in CGGV18 , is defined as follows

[TABLE]

The distance $L/4$ was chosen so as to make the two spins $\sigma_{(x_{b},y)}$ and $\sigma_{(x_{b},y-L/4)}$ sufficiently uncorrelated. The second b.c. we consider is the so-called “free boundary” (FB) condition, in which one takes

[TABLE]

Note that in both the L4 and FB b.c.s, the properties of the reservoirs do not enter directly. The third b.c. we shall consider will be referred to, hereafter, as the “stochastic boundary” (SB) condition. The SB rule dictates that the probability that a ghost spin takes a certain value is ruled by the magnetization $m_{\textrm{res}}$ of the reservoir next to $\sigma_{i}^{o}\in\{-1,+1\}$, by assuming that the average value of $\sigma_{i}^{o}$ is $m_{\textrm{res}}$, and reads

[TABLE]

where $m_{\textrm{res}}\in\{m_{-},m_{+}\}$ is the magnetization of the reservoir close to $x=x_{b}$.

II.3 Spin-updating mechanisms

We shall now consider two different spin-updating mechanism at the vertical boundaries of the lattice, that mimic the action of the two reservoirs on the system. As shown in Fig. 1, the system is coupled, to its horizontal boundaries, to the two reservoirs $\mathcal{R}_{+}$ and $\mathcal{R}_{-}$, having magnetization equal, respectively, to $m_{+}\in[0,1]$ and $m_{-}=-m_{+}$.

The first mechanism we shall consider in this work is the one considered in CGGV18 , called independent spin flip (ISF). According to the ISF mechanism, a selected boundary spin $\sigma_{i}$, $i=(x_{b},y)$ is updated to its new value $\sigma^{\prime}_{i}\in\{-1,1\}$ (which may coincide with $\sigma_{i}$) with a rate $c^{\textrm{ISF}}$ which is independent from the current state $\sigma_{i}$ and only dictated by the magnetization of the adjacent reservoir, i.e.:

[TABLE]

where $m_{\textrm{res}}\in\{m_{-},m_{+}\}$ is the magnetization of the reservoir close to $x=x_{b}$. We thus flip $\sigma_{i}$ with probability $(1-\sigma_{i}m_{\textrm{res}})/2$. Note further that each reservoir plays a double role with the ISF-SB dynamics: (i) it updates the boundary spins $\sigma_{i}\rightarrow\sigma^{\prime}_{i}$, $i=(x_{b},y)\in\Lambda$ with a rate defined by Eq. (10), (ii) it acts as a real “boundary” surrounding the lattice $\Lambda$ by fixing the stochastic b.c. $\sigma^{o}$ (i.e. the stochastic character of all the ghost spins in $\Lambda^{o}$), according to Eq. (9).

A second, so-called NSF spin-updating mechanism is obtained by considering a model that is fully decoupled from the reservoirs, namely no spin flip (NSF) mechanism takes place at the leftmost and rightmost vertical boundaries of the lattice. For this model, the bulk dynamics is still defined by Eq. (6), but Eq. (10) is replaced by

[TABLE]

Since reservoirs are not part of the description for the NSF dynamics, the SB b.c. as defined in Section II.2 is no longer relevant here. However, we can still define the SB b.c. in the NSF case by considering Eq. (9) with a given $m_{+}$ (and $m_{-}=-m_{+}$) even if these values are not associated to a reservoir. It is readily seen that the NSF dynamics conserves, at any time step, the total magnetization, that is fixed by the initial configuration. The presence of a spin-updating mechanism induced by the reservoirs, or the lack thereof, affects significantly the critical value of magnetization corresponding to a regime of vanishing current, as it will be seen in Section III.2.

Let us also shortly dwell, here, on some of the features of the so-called detailed balance (DB) dynamics CMW11 ; ELS1990 . We say that the spins located at one boundary of the lattice are (locally) in equilibrium with the nearest reservoir if the rate $c^{\textrm{DB}}$ at which a spin value $\sigma_{i}^{\prime}$ is introduced at the boundary site $i$ obeys the following DB condition:

[TABLE]

where $\sigma^{\prime}$ denotes the configuration obtained from $\sigma$ by updating the spin $\sigma_{i}$ to $\sigma_{i}^{\prime}$, and the term $(1+m_{\textrm{res}}\sigma_{i}^{\prime})/2$ corresponds to the probability of drawing at random the spin $\sigma_{i}^{\prime}$ from the reservoir with magnetization $m_{\textrm{res}}$. We may then use an identity that follows from Eq. 5,

[TABLE]

to rewrite Eq. (12) more conveniently as

[TABLE]

where

[TABLE]

corresponds to the hamiltonian of an Ising model with opposing magnetic fields $h_{\textrm{res}}$ acting on all boundary spins.

II.4 Observables

The presence of the two reservoir at the boundaries of the domain $\Lambda$ with different magnetizations induces a magnetization flow across the system. The magnetization current along a given horizontal bond is measured by counting the number of positive spins that cross this bond from left to right minus the number of those crossing this bond in the opposite direction, and dividing this number by time. The stationary current $J$ then corresponds to the long time limit of this rate CGGV18 , averaged over all bonds, or alternatively, averaged over boundary bonds only, provided a stationary limit exists.

A typical stationary state configuration can be surveyed by introducing the stationary magnetization profile $\overline{m}_{x}$, $x=1,\ldots,L$, that is obtained as follows: we measure the average magnetization along the column $x$

[TABLE]

at each MC step $t$, and then we take the average $\overline{m}_{x}$ over the set of collected values in the course of MC, with $T$ steps in total,

[TABLE]

In practise, the expensive sums are never evaluated during the MC, as each spin or bond flip causes only changes of $m_{x}(t)$. The mean boundary magnetizations are $\overline{m}_{1}$ and $\overline{m}_{L}$. For the ${\cal E}$–models $\overline{m}_{L}$ is denoted by $m_{\textrm{eq}}$. The notation used in this work is summarized in Tab. 2.

By combining reservoir mechanisms and b.c.s we have defined nine different models, a set of three equilibrium models ${\cal E}\equiv\{$NSF-L4, NSF-FB, NSF-SB$\}$, a set containing the local equilibrium models ${\cal L}\equiv\{$DB-L4, DB-FB, DB-SB$\}$, and the corresponding three nonequilibrium models ${\cal N}\equiv\{$ISF-L4, ISF-FB, ISF-SB$\}$. The investigation of the ${\cal E}$, ${\cal L}$, and ${\cal N}$ dynamics, equipped with the different b.c.s, will allow us to comment, in Section III, on the relation between the magnetizations $m_{\textrm{crit}}$ and $m_{\textrm{eq}}$, mentioned above, for different values of $\beta$.

II.5 Implementation details

All results to be presented are obtained for 2D grids with $L=40$. Each MC step corresponds to a single attempted bond or boundary spin flip at a randomly selected site. If not otherwise mentioned, (i) the initial configuration, referred to as “configuration A”, has $\sigma_{i}=-1$ for $x\in\{1,..,20\}$ and $\sigma_{i}=+1$ for $x\in\{21,..,40\}$ (the configuration shown for $h_{+}=5.35$ in Fig. 4 later below), and (ii) each simulation for given set of parameters $\beta$, $m_{+}$ runs for $\tau=10^{13}/N_{b}$ steps. Contributions to observables like the current $J$ and magnetization profiles $\overline{m}_{x}$ are computed at each MC step, as each bond of spin flip gives rise to an increment of decrement of these quantities as opposed to the total energy, which is for this reason only calculated each $N_{b}$ MC steps. To save computing time further all nine models ${\cal E}$, ${\cal L}$, and ${\cal N}$ are calculated during a single run, and share their random numbers. This way identical bonds or spins are selected for an attempted MC move simultaneously, and the cost for the calculation of acceptance rates is minimized. Exponentials $\exp(-\beta\Delta H)$ are nowhere calculated but looked up from its very limited set of possible values. Relevant neighbors for the calculation of $\Delta H$ are saved once for each bond and each model, which completely eliminates the calculation of b.c.s such as L4 or periodic b.c.s in $y$–direction. We have calculated the current $J$ independently from spins at the left and right boundary. For all stationary results to be presented, they are identical within statistical errors.

To capture the details as a function of $m_{+}$, we plot results versus the dimensionless quantity $\beta h_{+}\equiv\tanh^{-1}(m_{+})$. This representation is advantageous to blow up the region of $m_{+}$ close to 1, where interesting phenomena appear in a very narrow interval of $m_{+}$ CGGV18 . To ensure that measurement points are equidistantly spaced in this representation we choose $\beta h_{+}$ equidistantly spaced between [math] and $6.5$, at a resolution of $0.05$. In addition, we calculate the result for $m_{+}=1$ (corresponding to $\beta h_{+}=\infty$) and eventually add it to plots as a filled marker. The 2D maps are obtained at a resolution of $\Delta\beta=0.01$ and $\Delta(\beta h_{+})=0.05$. A function like $m_{\textrm{crit}}(\beta)$ is extracted from the $J$–map as contour line at $J=0$. The magnetization $m^{\textrm{eff}}_{+}$ is extracted by first inverting the $\overline{m}_{L}^{\textrm{ISF-L4}}$ map, using a continuous interpolant. The quantities $m^{\textrm{eff}}_{+}$ and $\overline{m}_{L}^{\textrm{ISF-L4}}$ are introduced in Sec.III.2, see also Tab.2.

III Results and discussion

In this section we present the results of our numerical simulations of the nonequilibrium steady state (NESS) attained by the ${\cal N}$–models ISF-L4, ISF-FB, and ISF-SB introduced in Section II. Moreover, when convenient, we shall also compare the results obtained with the same b.c.s in presence of DB dynamics, i.e., ${\cal L}$–models DB-L4, DB-FB, and DB-SB. The equilibrium ${\cal E}$–models NSF-L4, NSF-FB and NSF-SB will also be considered for the purpose of comparison. We aim at elucidating the nature of the NESS and the role played by the details of the b.c.s and reservoir mechanisms. In particular, we will focus on the stationary magnetization profiles $\overline{m}_{x}$ and currents $J$, defined in Section II.4. In all the simulations to be considered below, the initial datum is the “configuration A” introduced in Section II.5, and whose total magnetization is zero.

Since our models depend on $L$, $\beta$ and $m_{+}$, the corresponding NESS will also depend on the same set of parameters in a manner which is far from obvious. It is therefore of crucial importance to choose properly these quantities in our simulations. The dependence of the NESS on the inverse temperature $\beta$ is likely to be dictated by the presence of the ferromagnetic phase transition in the equilibrium system. Thus, different behaviors are to be expected in the two temperature regimes: standard phenomena (e.g. Fick’s law) at high temperatures and a more complex and intriguing picture in the low temperature phase. As representative instances of the two different regimes, we have chosen $\beta=0.3$ and $\beta=1$, respectively lower and higher than the critical value of $\beta_{c}$ in Eq. (1).

We recall that the linear size of the lattice $\Lambda$ in our simulations, for which results are presented, is chosen as $L=40$, as in CGGV18 . Indeed, to the best of our knowledge, this lattice appears sufficiently large to suppress the most evident finite-volume effects that can be observed for $L$ ranging between 10 and 40, while being sufficiently small to allow for a detailed numerical study within an acceptable computational workload. However, we cannot exclude the existence of significant finite-size effects for systems with sizes of the order of our size $L=40$. These issues will be discussed further in what follows.

The main results of CGGV18 , that will be extended and deepened here, concern the behavior of the stationary current $J$ of the ISF-L4 model, as the magnetization $m_{+}$ ($m_{-}$) of the reservoir(s) is varied. We recall them briefly here.

As the value $m_{+}$ is decreased from its maximum value 1, the flux $J(m_{+})$ is first negative and, crossing a critical value $m_{\textrm{crit}}$, it becomes positive. More precisely, for $m_{+}>m_{\textrm{crit}}$ the current is negative, flowing from the reservoir with positive magnetization $\mathcal{R}_{+}$ to the one with negative magnetization $\mathcal{R}_{-}$ (in agreement with the Fick’s law). For $m_{+}<m_{\textrm{crit}}$ the current is reversed, flowing against the magnetization gradient, that is from $\mathcal{R}_{-}$ to $\mathcal{R}_{+}$ (in violation of the Fick’s law). In the latter case we have the phenomenon called uphill diffusion. 2. 2.

The flip of the current in passing from downhill to uphill diffusion is connected to a change in the structure of the NESS CGGV18 . This can be detected by studying the stationary magnetization profile $\overline{m}_{x}$ or by inspection of the typical spin configuration of the NESS. In CGGV18 two typical magnetization profiles have been identified: instanton and bump.

The instanton profile corresponds, in the spin configurations, to the presence of an interface separating two plus and minus equally sized phases close to $\mathcal{R}_{+}$ and $\mathcal{R}_{-}$ respectively, that is with the interface placed in the middle of the lattice (cf. also SOS16 for an interpretation of the interface based on nonequilibrium thermodynamics). A sketch of the instanton profile, in a continuous limit representation $m(r)$, where $r=(x-1)/(L-1)$, is given in Fig. 2. The figure shows that, for $r$ moving from the right reservoir $\mathcal{R}_{+}$ ($r=1$) to the center of the lattice ($r=1/2$), the profile $m(r)$ decreases from $m_{+}$ (the magnetization imposed by $\mathcal{R}_{+}$) to a value which is expected to be the equilibrium spontaneous magnetization $m_{\beta}$ given in Eq. (2). At the center of the system, $m=1/2$, the magnetization jumps down to $-m_{\beta}$ and then keeps decreasing up to the left reservoir ($r=0$), where it takes the value $m_{-}$. Obviously, at finite volumes the sharp jump is replaced by a layer, i.e. a narrow region of sites in which the profile changes rapidly. Examples of steady magnetization profiles with an instanton-like shape are shown in the panel with $\beta h_{+}=5.35$ of Fig. 3 for the ISF-L4, ISF-FS, ISF-SB models. The corresponding spin configuration, for the ISF-L4 case, is shown in Fig. 4, where it is also evident that $J$ is negative in this state.

The bump profile occurs when the typical spin configuration of the NESS presents an interface close to one of the vertical boundaries. In this case a sea of (say) positive spins emanating from $\mathcal{R}_{+}$ invades the most part of the lattice $\Lambda$, leaving a small strip adjacent to $\mathcal{R}_{-}$ for the negative spins. The profile of $\overline{m}_{x}$ for $x$ ranging from the right ($x=L$) to the left ($x=1)$ boundary, increases steadily from $m_{+}$ up to a maximum value $m_{\textrm{bump}}>m_{+}$, which is reached for $x$ close to 1. Then $\overline{m}_{x}$ jumps suddenly from $m_{\textrm{bump}}$ down to $m_{-}$ at $x=1$, forming a boundary layer close to the left boundary. A similar profile is obtained by changing in the previous description the sign of magnetization and exchanging left with right. In any case, in the presence of a bump the current flows in the ‘wrong’ direction, producing uphill diffusion. The bump appears, for instance, in the ISF-L4 model for $\beta h_{+}=3.75$, see Fig. 3. In Fig. 4 it is shown, for the same value of $\beta h_{+}$, the spin configuration with the boundary layer close to $\mathcal{R}_{+}$ and the evidence for the uphill current, i.e. $J>0$. 3. 3.

The instanton profile, that sustains a negative current, is stable for $m_{+}\in(m_{\textrm{crit}},1)$. This means that, perturbing the configuration A initial datum or even taking a random initial condition for $m_{+}>m_{\textrm{crit}}$, the dynamics leads asymptotically to the instanton profile for $\overline{m}_{x}$. This is the stable phase. The situation is quite different and much more intricate for $m_{+}\in(0,m_{\textrm{crit}})$. Here the instanton loses stability and we have numerical evidence that, at least for not too low $m_{+}$ values, the bump profile is stable. Thus, there should exist a further critical value $m_{u}<m_{\textrm{crit}}$ that defines the stability region for the bump $(m_{u},m_{\textrm{crit}})$. The nature of $m_{u}$ is not well understood, but it is conjectured that its presence should be a finite-volume effect, since in the large volume limit the interval $(m_{u},m_{\textrm{crit}})$ is supposed to shrink to the empty set. Indeed, some heuristic reasoning leads to the estimate $O(L^{-2/3})$ for its length. This phase has been termed metastable in CGGV18 . Below $m_{u}$ the picture is much more complex and less studied. In the next section we will show that decreasing further $m_{+}$ towards 0, the system undergoes an intricate sequence of different stationary states. This is the unstable phase, that was also referred to as the “weakly unstable” or “chaotic region” in CGGV18 .

III.1 Behavior of $J$ in external magnetization-temperature space

We shall discuss here the 2D maps for the current $J$ in the parameter space spanned by $\beta h_{+}=\tanh^{-1}(m_{+})$ and $\beta$ for the 3 different models ISF-4, ISF-FB and ISF-SB. The MC results are portrayed in Fig. 5 (left column). Remarkably, the ISF-L4 and ISF-SB dynamics seems to yield similar results for the current, while the FB b.c.s gives rise to a different scenario, in which the onset of uphill currents is somehow hindered. It should also be noted that, with the L4 and SB b.c.s, the Onsager’s result for $m_{\beta}$, denoted by the yellow line in Fig. 5 is close, for large enough values of $\beta$, to $m_{\textrm{crit}}$ since it marks closely the transition between regions with positive currents (colored in red, in the figure) and other regions with negative currents (colored in green). For smaller values of $\beta$, the yellow line no longer stays at the border between red and green regions, hence $m_{\beta}$ and $m_{\textrm{crit}}$ substantially differ from one another, as it will also be seen in Fig. 6.

In Section III.2 we shall focus, in particular, on the behavior of the current as a function of $\beta h_{+}$ along the vertical lines at $\beta=0.3$ and $\beta=1$. An inspection of Fig. 5 reveals that positive ’uphill’ currents $J>0$ (colored in red) are observed, with the ISF mechanism, when the parameters $\beta$ and $m_{+}$ of our model are suitably tuned. In particular, the parameter $\beta$ is required to be larger than a certain critical value $\beta_{\textrm{crit}}$. Note that the latter, due to finite size effects, might even significantly differ from the critical value $\beta_{c}$ in Eq. (1) determined by Onsager (pertaining to the infinite volume Ising model). The ISF mechanism is an essential ingredient to observe uphill currents. It is worth comparing the ISF results shown in Fig. 5 (left column) with the corresponding results obtained with the DB dynamics equipped with the L4, FB and SB b.c.’s, see Fig. 5 (right column). An inspection of Fig. 5 confirms that while the L4 and FB b.c.’s give rise to a purely Fickian behavior, the DB dynamics equipped with the SB b.c. defined in Eq. (9) yields uphill currents. The onset of positive currents even in presence of DB dynamics stems from the choice of the aforementioned boundary condition, which, as it stands, does not match properly with the DB dynamics. This aspect will be clarified in Section III.3.

III.2 Critical value $m_{\textrm{crit}}$, NSF equilibrium magnetization $m_{\textrm{eq}}$ and DB equilibrium magnetization $m_{o}$

The critical value $m_{\textrm{crit}}$ separating the stable and the metastable phases is the threshold at which the current changes sign, i.e. where the current vanishes, $J(m_{\textrm{crit}})=0$. As already said in the Introduction, in CGGV18 it has been conjectured that in the large volume limit $m_{\textrm{crit}}$ should approach the equilibrium magnetization $m_{\beta}$. However, in order to take into account the finite size effects, in CGGV18 it is also claimed that $m_{\textrm{crit}}$ can be measured by evaluating the magnetization $m_{\textrm{eq}}$ of the rightmost column of $\Lambda$ in the absence of the reservoirs, i.e. at equilibrium (NSF models). The choice of the L4 b.c.s should entail that in the large volume limit $m_{\textrm{eq}}$ approaches $m_{\beta}$. Some numerical evidence of this fact is given in CGGV18 for the ISF-L4 case, where $m_{\textrm{eq}}$ is computed for $L$ ranging from 10 up to 40 and $\beta=1$. Here, we study the dependence of $m_{\textrm{eq}}$ and $m_{\textrm{crit}}$ on the parameter $\beta$, and compare the ISF-L4 case with ISF-FB and ISF-SB. We denote by $m_{\textrm{crit}}^{\textrm{X}}$, with X $\in{\cal N}$, the critical magnetization computed according to the nonequilibirum model $X$. Analogously, $m_{\textrm{eq}}^{\textrm{X}}$, with X $\in{\cal E}$ is the quantity $\overline{m}_{L}$ defined by Eq. (17) computed in the equilibrium model X.

For later use, we introduce also a quantity, named $m_{o}$, which is defined as $\overline{m}_{L}^{\textrm{DB-L4}}$ at $h_{+}=0$. We remark that, while $m_{\textrm{eq}}$ is defined for the equilibrium bulk system in which the sole conservative Kawasaki dynamics acts, in defining $m_{o}$ we retain, in the presence of local equilibrium, the action of the two reservoirs even though with the same magnetization $m_{+}=m_{-}=0$ (corresponding to $h_{+}=0$). In this situation there is no net flux across the system induced by the gradient of magnetization at the boundaries, and a NESS is not established. Thus we argue that, at least for sufficiently large systems, the two systems (i.e. the Canonical NSF-L4 and the Grand Canonical DB-L4 with $h_{+}=0$) should be equivalent and then $m_{\textrm{eq}}$ and $m_{o}$ should agree.

In Fig. 6 we show the critical vales $m_{\textrm{crit}}^{\textrm{X}}$ for $\textrm{X}\in{\cal N}$, the spontaneous magnetization $m_{\beta}$, see Eq. (2), and $m_{\textrm{eq}}^{\textrm{X}}$, X $\in{\cal E}$ for $\beta\in(0,1)$. $m_{\textrm{crit}}^{\textrm{X}}$ is the $m_{+}$–value closest to $1$ at which the current vanishes. The figure shows that the status of the conjecture regarding $m_{\textrm{eq}}^{\textrm{NSF-L4}}$ and $m_{\textrm{crit}}^{\textrm{ISF-L4}}$ in CGGV18 is strongly dependent on $\beta$. Indeed, for $\beta$ small (say $\beta<0.6$) the two quantities are appreciably different one from the other and different from $m_{\beta}$. On the other hand, $m_{\textrm{eq}}^{\textrm{NSF-L4}}$ and $m_{\textrm{crit}}^{\textrm{ISF-L4}}$ get closer and closer to $m_{\beta}$ as $\beta$ is increased. For $\beta=1$, the value studied in CGGV18 , they almost coincide. The right panel of Fig. 6 also shows the behavior of $m_{o}$, defined for the DB-L4 dynamics, as a function of $\beta$.

The role of b.c.s is also evident: for all $\beta$, the critical value $m_{\textrm{crit}}^{\textrm{ISF-FB}}$ and $m_{\textrm{eq}}^{\textrm{NSF-FB}}$ are clearly different from $m_{\textrm{crit}}^{\textrm{ISF-L4}}$ and from $m_{\textrm{eq}}^{\textrm{NSF-L4}}$ and $m_{\beta}$. Moreover, they seem to approach each other, as $\beta$ increases. It is also evident that ISF-L4 and ISF-SB coincide in the whole range of $\beta$. This is the first instance of a situation that we will meet frequently, i.e. the equivalence of ISF-L4 and ISF-SB models. Finally, we observe that all the quantities relative to models with b.c.s (i.e. L4 and SB), converge to $m_{\beta}$ as $\beta$ increases.

III.3 Average magnetization at the boundaries

Another relevant feature of the ISF mechanism, outlined in CGGV18 , is that the stationary magnetization (averaged over the vertical direction) on the rightmost column, denoted as $\overline{m}_{L}$, see Tab.2, attains a value close to $m_{+}$ when $\beta$ is large. Figure 7 contains the result of MC simulations showing the behavior of $\overline{m}_{L}$ as a function of $\beta h_{+}$ for different values of $\beta$, and for the various b.c.s considered above for the dynamics undergoing the ISF (top panes) and the DB (botom panels) updating mechanism. It is worth noticing that, for fixed $\beta$ and $\beta h_{+}$, the b.c.’s seem not to affect significantly the resulting value of $\overline{m}_{L}$ with the ISF updating mechanism. Another observation concerns the dependence of $\overline{m}_{L}$ on $\beta$, for a given $\beta h_{+}$. Figure 7 shows that, with the ISF mechanism (top panels) $\overline{m}_{L}$ attains a value that is closer and closer to $m_{+}$ when $\beta$ is large, e.g. at $\beta=1$. For small values of $\beta$ the two values $\overline{m}_{L}$ and $m_{+}$ start to deviate significantly, because then spin diffusion from the boundaries to the bulk of the lattice plays a major role.

It is also worth noticing that in Fig. 7 the behavior of the DB-L4 and DB-SB models differ from one another. The origin of the peculiar behavior of the DB-SB model lies in the SB b.c.s defined in Eq. (9), which equips the ghost spins with an average magnetization close to $m_{+}$. On the other hand, because of the local equilibrium between the reservoir and the rightmost column that is produced by the detailed balance, ${\cal R}_{+}$ does not impose to $\overline{m}_{L}$ its “nominal” magnetization $m_{+}$ (as with the ISF dynamics): the measured value of $\overline{m}_{L}$ is, indeed, typically larger than $m_{+}$ (cf. e.g. the plot referring to the DB-L4 dynamics). The outcome of the two competitive effects is shown in the bottom right panel of Fig. 7.

Remarkably, only Fickian diffusion is observed with DB-L4 dynamics. To see that, let us first introduce the effective magnetization $m^{\textrm{eff}}_{+}$, defined as the reservoir magnetization which, with a ISF dynamics, gives rise to the same value of boundary magnetization $\overline{m}_{L}$ obtained with a DB dynamics (equipped with the same b.c.’s of the ISF dynamics) in presence of a reservoir magnetization $m_{+}$. Precisely, considering the L4 b.c.’s, one has that $m^{\textrm{eff}}_{+}$ solves the following relation:

[TABLE]

Thus, Fig. 8 reveals that, for any $\beta$, $m^{\textrm{eff}}_{+}$ is larger than $m_{\beta}$, which, as shown in Eq. (6), is also close to $m_{\textrm{crit}}^{\textrm{ISF-L4}}$. The result is, then, the onset of a stable phase, accompanied by a downhill current, cf. Section III.

Turning to the DB-SB model, the effective magnetization obtained from the solution of an equation similar to Eq. (18) turns out being smaller, for certain values of $m_{+}$, than $m_{\textrm{crit}}^{\textrm{ISF-SB}}$, in which case uphill currents can indeed be observed, see Fig. 5. Our results seem to indicate that purely Fickian diffusion is restored with the above DB-SB dynamics if one defines differently the SB boundary condition in Eq. (9), namely by replacing $m_{+}$ by $m^{\textrm{eff}}_{+}$ obtained from the solution Eq. (18), cf. the right panel of Fig. 8 (the corresponding 2D map is not shown since it looks similar to the right upper panel of Fig. 5).

Furthermore, the following argument allows us to interpret the behavior of $\overline{m}_{L}$ as a function of $m_{+}$ at $\beta=0$ for the ISF models. At $\beta=0$ all attempted bulk bond flips are accepted, $c_{i,j}^{\textrm{bulk}}=1$, and the dynamics is equivalent to that of a simple exclusion process simpleexclusion . This implies that the magnetization profile $\overline{m}_{x}$ varies strictly linearly with $x$, as it can be seen, for instance, by using duality CGGR2013 . The two coefficients characterizing the linear profile are obtained from the following two conditions: (i) For symmetry reasons $\overline{m}_{1}=-\overline{m}_{L}$, i.e., $\overline{m}_{(L+1)/2}=0$, and (ii) the magnetization at the boundary is determined with equal weight by bond flips from the adjacent layer with magnetization $\overline{m}_{L-1}$ and by the ISF updating scheme enforcing magnetization $m_{+}$. The resulting linear magnetization profile $\overline{m}_{x}$ and its values at the boundaries are

[TABLE]

As the case of $\beta=0$ provides a lower bound, this is the $\beta$-independent result for large system sizes. The observed departures reflect the finite system size and can actually be used to measure the system size.

In the following subsections we will present the results of a set of simulations at high ($\beta=0.3$) and low ($\beta=1$) temperatures. For both cases we will describe the asymptotic magnetization profile and the corresponding current. These result are taken from a much larger bunch of simulations and are selected in order to offer an highlight of the complex dynamic presented by our models. As it could be expected, the scenario below the critical temperature is much complex than above, where no anomalies occur.

III.4 High temperature regime

We consider $\beta=0.3$ as a representative case of the high temperature regime.

In Fig. 9 we display the stationary magnetization profiles $\overline{m}_{x}$ for different values of $\beta h_{+}$. In each panel the three models ISF-L4, ISF-FB and ISF-SB are considered. As it is evident from the figure, the profiles of the three different b.c.s are qualitatively comparable in that they all are non decreasing, at least for not too small $\beta h_{+}$’s, say for $\beta h_{+}\geq 0.6$. However, the b.c.s play a role: while the L4 and SB are essentially indistinguishable in all the range of $\beta h_{+}$’s, the FB coincides with the others only at high $\beta h_{+}$, including $m_{+}=1$, but decreasing $\beta h_{+}$ the curvature of FB profile appears clearly different from the others. For $\beta h_{+}<0.6$ the L4 and SB curves develop, eventually through the appearance of slightly pronounced local maxima and minima, an almost flat region in the bulk. A similar behavior is developed also for the FB case, but for smaller $\beta h_{+}$. Such profiles sustains a negative current, as we can see in Fig. 10 where $J$ is plotted versus $\beta h_{+}$. Moreover, $J$ looks as a smooth curve that decreases monotonically as $\beta h_{+}$ increases from [math]. As $\beta h_{+}\to 0$ the systems get close to equilibrium, being $m_{+}$ and $m_{-}=-m_{+}$ close since $m_{+}\to 0$. Thus $\overline{m}_{x}$ approaches the flat equilibrium profile, as observed above, and the flux approaches 0.

The analysis of the corresponding spin configurations shows a weak phase separation (without sharp interface, see the slope in the center of the curve for $\beta h_{+}=2.1$) that becomes less and less visible approaching $\beta h_{+}=0$. This indicates that, as the system approaches equilibrium, the positive and negative phases become more and more random and intertwined: the configuration loses the instanton profile structure. The confirmation of the fact that as $\beta h_{+}$ decreases, the spin configuration becomes more disordered is shown in Fig. 11, where the mean total Hamiltonian per bond $\overline{H}/N_{b}$ is represented for the ISF-L4, ISF-FB and ISF-SB models. Since the Hamiltonian can be written as $H(\sigma|\sigma^{0})=2N_{-}-N_{b}$, where $N_{-}$ is the number of broken bonds (a broken bond has antiparallel spins), this quantity, which is in between $-1$ and $1$, gives a measure of the amount of broken bonds on the total lattice bonds. Its value is -1 when no bond is broken and 1 when all the bonds are broken, while for a sharp vertical interface (with essentially $L$ broken bonds), this ratio is $(1-2L)/(1+2L)$, whose value is $\approx-0.975$ for $L=40$. Also in this respect the L4 and SB models are the same, while FB follows a different path, though qualitatively similar.

III.5 Low temperature regime

As a case study for the low temperature regime, we consider $\beta=1$, and write $h_{+}$ instead of $\beta h_{+}$ in this section, to simplify notation.

The analysis of the time averaged magnetization profiles shown in Fig. 3, offers a rich scenario as $h_{+}$ decreases: starting from the instanton profile, which is common to all three ISF models in the stable parameter region (see panel for $h_{+}=5.35$), we find a metastable phase where the instanton is replaced by the bump (see panel for $h_{+}=3.7$). Then continuing to decrease $h_{+}$, the system enters into the weakly unstable phase with double-bump profiles (see panel for $h_{+}=2.25$). In this case two (almost) flat regions with values close to $m_{+}$ and $m_{-}$ develop in the vicinity of the boundary with the opposite magnetization. Then moving towards the boundary, the profile of $\overline{m}_{x}$ jumps sharply to the value imposed by the reservoir. This behavior corresponds to a configuration with an interface in the middle, but with the plus and minus phases flipped, i.e. close respectively to ${\cal R}^{-}$ and ${\cal R}^{+}$. Because of that, two boundary layers appear. A configuration with this behavior is clearly visible, in Fig. 4 for $h_{+}=2.24$. Decreasing $h_{+}$ further, sharper maxima and minima appear, see e.g. $h_{+}=1.05$.

Observe that for a given $h_{+}$, while the behaviors of ISF-L4 and ISF-SB models are similar (maybe after the application of the left-right and spin flip symmetries, see e.g. $h_{+}=1.35$), the free b.c.s of the ISF-FB model generate quite different nonequilibrium stationary states. However, as in the high temperature case, the three curves seem to reconcile as $h_{+}\to 0$ (that is approaching the equilibrium) when in the bulk the profile tend to flatten, although with two peaks close to the reservoir that resemble those found in Rieder et al. Rieder1967 or the widely observed boundary resistance effects.

Figure 4 shows the current for ISF-L4 model: in contrast with the simplest high temperature case, in this low temperature region the current curve, reflecting the complexity of the magnetization profiles, is not smooth and not monotonic. It crosses three times the zero line and shows the presence of an hysteretic region (close to $h_{+}=4$) where the observed stationary state depends on the initial configuration chosen in the MC simulation. A blow up of this phenomenon is shown in Fig. 12. Cyan circles mark results obtained after $\tau=10^{13}/N_{b}$ MC steps starting from the default initial configuration A. The open squares are obtained by slowly increasing $h_{+}$, starting from the cyan configuration at $h_{+}=3.6$, while the filled squares are obtained by slowly decreasing $h_{+}$, starting from the ’cyan’ configuration at $h_{+}=3.7$. At these rates the hysteretic region is seen to extend over the range $h_{+}\in[3.46,3.91]$.

The spin configurations corresponding to the magnetization profiles of ISF-L4 (represented by red square lines in Fig. 3) are shown in the surrounding panels of Fig. 4. Note that the study of the mean total hamiltonian per bond $\overline{H}/N_{b}$ shows, see Fig. 13, wide intervals in which the curves are almost constant, meaning that in those regions the interfaces separating the plus and minus phase are maintained. It is worth noting the difference with the high temperature case, see Fig. 11, in which $\overline{H}/N_{b}$ varies smoothly without flat regions.

Figure 14 represents the current $J$ as a function of $\beta h_{+}$ for the two cases $\beta=0.3$ (left panel) and $\beta=1$ (right panel). Together with the current curves already shown in Figs. 10 and 4 for ISF-L4 model, there are the ones for ISF-FB and ISF-SB models. From the comparison it is confirmed that also for this quantity the model ISF-SB behaves like the model ISF-L4, while ISF-FB model differs from the first two, while maintaining qualitatively similar trends within the same temperature regime. Again, the difference between high and low temperature cases is evident. Figure 15 shows that at low temperature ($\beta=1$) there are intervals of $\beta h_{+}$ where the absolute value of the average magnetization $|\overline{m}|$ is non-zero, while for high temperature the mean magnetization vanishes for all $\beta h_{+}$.

IV Conclusions

We have investigated the behavior of the current $J$ and the boundary magnetization $\overline{m}_{L}$ for some 2D Ising models, coupled to external magnetization reservoirs and equipped with different b.c.s, called respectively L4, FB and SB.

In particular, we examined the sign of the stationary currents in the parameter space, spanned by the variables $\beta$ and $m_{+}$. Our MC simulations indicate that stationary uphill currents do appear for certain values of the parameters of the model, as a result of the ISF spin-updating mechanism, which breaks the condition of detailed balance. We also highlighted the relation between $m_{\textrm{crit}}$ and $m_{\textrm{eq}}$, the first being an observable characterizing the nonequilibrium dynamics of the 2D Ising model in contact with external reservoirs, while the latter refers to an equilibrium Ising model with conservative dynamics. Our simulations show that, if $L$ is fixed, the two quantities both tend to Onsager’s magnetization $m_{\beta}$ when $\beta$ is large.

Moreover, we studied in detail the behavior of the equation of state, namely the relation $J$ vs. $m_{+}$, for $\beta=0.3$ and $\beta=1$, as representative, respectively, of the high temperature and low temperature regimes. Lowering the parameter $m_{+}$ leads to novel stationary states, in which the current changes sign multiple times and more interfaces appear in the microscopic spin configurations. One open question, not addressed in this manuscript, concerns the presence of uphill diffusion in the infinite volume limit. More theoretical work is needed to clarify whether uphill currents persist in the above limit.

The analysis of 2D Ising models on a finite lattice, coupled to external reservoirs that break the condition of detailed balance, may be relevant in a variety of applications, e.g. in the investigation of mesoscopic systems, in which the notion of “local equilibrium” is not guaranteed. We thus expect uphill currents to play a major role for these systems, that do not follow the basic tenets of irreversible thermodynamics.

Acknowledgments.

The authors wish to thank Anna De Masi (University of L’Aquila) and Errico Presutti (Gran Sasso Science Institute) for inspiring this work and for the many useful comments and discussions. C. Giardinà is acknowledged for useful discussion. MC acknowledges financial support from FFABR 2017. CG and CV acknowledge financial supports from Fondo di Ateneo per la Ricerca 2016 and 2017 (UniMoRe).