A Chain, a Bath, a Sink, and a Wall

Abstract

We numerically investigate out-of-equilibrium stationary processes emerging in a Discrete Nonlinear Schrödinger chain in contact with a heat reservoir (a bath) at temperature $T_L$ and a pure dissipator (a sink) acting on opposite edges. Long-time molecular-dynamics simulations are performed by evolving the equations of motion within a symplectic integration scheme. Mass and energy are steadily transported through the chain from the heat bath to the sink. We observe two different regimes. For small heat-bath temperatures $T_L$ and chemical-potentials, temperature profiles across the chain display a non-monotonous shape, remain remarkably smooth and even enter the region of negative absolute temperatures. For larger temperatures $T_L$, the transport of energy is strongly inhibited by the spontaneous emergence of discrete breathers, which act as a thermal wall. A strongly intermittent energy flux is also observed, due to the irregular birth and death of breathers. The corresponding statistics exhibit the typical signature of rare events of processes with large deviations. In particular, the breather lifetime is found to be ruled by a stretched-exponential law.

1. Introduction

The study of nonequilibrium thermodynamics of systems composed of a relatively small number of particles is motivated by the need for a deeper theoretical understanding of the statistical laws leading to the possibility of manipulating small-scale systems like biomolecules, colloids, or nano-devices. In this framework, statistical fluctuations and size effects play a major role and cannot be ignored as it is customary to do in their macroscopic counterparts.

Arrays of coupled classical oscillators are representative models of such systems and have been studied intensively in this context [1,2,3]. In particular, the Discrete Nonlinear Schrödinger (DNLS) equation has been widely investigated in various domains of physics as a prototype model for the propagation of nonlinear excitations [4,5,6]. In fact, it provides an effective description of electronic transport in biomolecules [7] as well as of nonlinear waves propagation in layered photonic or phononic systems [8,9]. More recently, a renewed interest for this multi-purpose equation emerged in the physics of gases of ultra-cold atoms trapped in optical lattices (e.g., see [10] and references therein for a recent survey).

The DNLS dynamics are characterized by two conserved quantities: the energy density h and the mass density a (also termed norm—see the following section for their definitions). Therefore, a generic thermodynamic state of the DNLS equation can be seen as a point in the $(a,h)$ plane. The equilibrium thermodynamics of the DNLS model was studied in a seminal paper by Rasmussen et al. [11] within the grand-canonical formalism. Here, it was shown that an equation of state can be formally derived with the help of transfer–integral techniques. Accordingly, any thermodynamic equilibrium state $(a,h)$ can be equivalently represented in terms of a temperature T and a chemical potential $\mu$ which can be numerically determined by implementing suitable microcanonical definitions [12,13]. In Reference [11], it was found that a line of infinite–temperature equilibrium states in the $(a,h)$ plane separates standard thermodynamic states, characterized by a positive T, from those characterized by negative absolute temperatures. The existence of negative temperatures can be traced back to the properties of the two conserved quantities of the system. From the point of view of dynamics, it was realized that negative temperatures can be associated to the presence of intrinsically localized excitations, named discrete breathers (DB) (see e.g., [14,15,16]). In a series of important papers, Rumpf developed entropy-based arguments to describe the asymptotic states above this infinite-temperature line [17,18,19,20]. It has been later found that negative–temperature states can form spontaneously via the dynamics of the DNLS equation. They persist over extremely long time scales, that might grow exponentially with system size as a result of an effective mechanism of ergodicity–breaking [21].

A related question is whether the dynamics of the system influences non-equilibrium properties when the system exchanges energy and/or mass with the environment. In a series of papers [21,22,23], it has been found that when pure dissipators act at both edges of a DNLS chain (a case sometimes called boundary cooling [24,25,26,27]) the typical final state consists of an isolated static breather embedded in an almost empty background. The breather collects a sensible fraction of the initial energy and it is essentially decoupled from the rest of the chain. The spontaneous creation of localized energy spots out of fluctuations has further consequences on the relaxation to equipartition, since the interaction with the remaining part of the chain can become exponentially weak [28,29]. A similar phenomenon occurs after a quench from high to low temperatures in oscillator lattices [30]. Also, boundary driving by external forces may induce non-linear localization [31,32,33].

When, instead, the chain is put in contact with thermal baths at its edges, non-equilibrium stationary states characterized by a gradient of temperature and chemical potential emerge [13]. Local equilibrium is typically satisfied so that the overall state of the chain can be conveniently represented as a path in the $(a,h)$ plane connecting two points corresponding to the thermodynamic variables imposed by the two reservoirs. The associated transport of mass and energy is typically normal (diffusive) and can be described in terms of the Onsager formalism. However, peculiar features, such as non-monotonous temperature profiles [34] or persistent currents [35], are found as well as a signature of anomalous transport in the low-temperature regime [36,37].

In this paper we consider a set-up where one edge is in contact with a heat reservoir at temperature $T_L$ and chemical potential ${\mu }_L$, while the other interacts with a pure dissipator, i.e., a mass sink. The original motivation for studying this configuration was to better understand the role of DBs in thermodynamic conditions. At variance with standard set-ups [1,2], this is conceptually closer to a semi-infinite array in contact with a single reservoir. In fact, on the pure-dissipator side, mass can only flow out of the system.

The presence of a pure dissipator forces the corresponding path to terminate close to the point $(a=0,h=0)$, which is singular both in T and $\mu$. This is indeed the point where the infinite- and the zero-temperature lines collapse. Therefore, slight deviations may easily lead to crossing the $\beta =0$ line, where $\beta =1/T$ is the inverse temperature. This is indeed the typical scenario observed for small $T_L$, when $\beta$ smoothly changes sign twice before approaching the dissipator (see Figure 1 for a few sampled paths). The size of the negative-temperature region increases with $T_L$ and appears to be stable over long time scales. Upon further increasing $T_L$, a different stationary regime is found, characterized by strong fluctuations of mass and energy flux. In practice, the negative-temperature region first extends up to the dissipator edge and then it progressively shrinks in favour of a positive-temperature region (on the other side of the chain). In this regime, the dynamics are controlled by the spontaneous formation (birth) and disappearance (death) of discrete breathers.

In Section 2 we introduce the model and briefly recall the definition of the main observables. Section 3 is devoted to a detailed characterization of the low-temperature phase, while the far-from-equilibrium phase observed for large $T_L$ is discussed in Section 4. This is followed by the analysis of the statistical properties of the birth/death process of large-amplitude DBs, illustrated in Section 5. Finally, in Section 6 we summarize the main results and comment about possible relationships with similar phenomena previously reported in the literature.

2. Model and Methods

We consider a DNLS chain of size N and with open boundary conditions, whose bulk dynamics is ruled by the equation

$$i{\dot{z}}_n=-2{\left|z_n\right|}^2{z}_n-z_{n-1}-z_{n+1}$$

(1)

where ($n=1,\dots ,N$) and $z_n=(p_n+i{q}_n)/\sqrt{2}$ are complex variables, with $q_n$ and $p_n$ being standard conjugate canonical variables. The quantity $a_n={\left|z_n\left(t\right)\right|}^2$ can be interpreted as the number of particles, or, equivalently, the mass in the lattice site n at time t. Upon identifying the set of canonical variables $z_n$ and $i{z}_n^{\ast }$, Equation (1) can be read as the equation of motion generated by the Hamiltonian functional

$$H=\sum _{n=1}^N\left(|z_n{|}^4+z_n^{\ast }{z}_{n+1}+z_n{z}_{n+1}^{\ast }\right)$$

(2)

through the Hamilton equations ${\dot{z}}_n=-\partial H/\partial \left(i{z}_n^{\ast }\right)$. We are dealing with a dimensionless version of the DNLS equation: the nonlinear coupling constant and the hopping parameters, which are usually indicated explicitly in the Hamiltonian (2), have been set equal to unity. Accordingly, also the time variable t is expressed in arbitrary adimensional units. Without loss of generality, this formulation has the advantage of simplifying numerical simulations.

A relevant property of the DNLS dynamics is the existence of a second conserved quantity beside the total energy H, namely the total mass

$$A=\sum _{n=1}^N{\left|z_n\right|}^2.$$

(3)

As a result, an equilibrium state is specified by two parameters, the mass density $a=A/N\ge 0$ and the energy density $h=H/N$. The first reconstruction of the equilibrium phase-diagram $(a,h)$ of the DNLS equation was reported in [11]. It is reproduced in Figure 1 for our choice of parameters, where the lower solid line defines the $(T=0)$ ground-state line $h=a^2-2a$ for different values of the mass density a. The ground-state corresponds to a uniform state with constant amplitude and constant phase-differences $z_n=\sqrt{a}{e}^{i(\mu t+\pi n)}$, with $\mu =2(a-1)$. States below this curve are not physically accessible. The positive-temperature region lies above the ground-state up to the red dashed line $h=2a^2$, which corresponds to the infinite-temperature $(\beta =0)$ line. In this limit, the grand-canonical equilibrium distribution becomes proportional to $\exp \left(\beta \mu A\right)$, where the finite (negative) product $\beta \mu$ implies a diverging chemical potential. Equilibrium states at infinite temperature are therefore characterized by an exponential distribution of the amplitudes $P\left(\right|z_n{|}^2)=a^{-1}{e}^{-|z_n{|}^2/a}$ and random phases. Finally, states above the $\beta =0$ line belong to the so-called negative-temperature region [11,21].

Finite-temperature equilibrium states do not allow for straightforward analytical treatments. However, one can determine the relation $a(T,\mu )$, $h(T,\mu )$ numerically, by putting the system in interaction with an external reservoir (see below) that imposes T and $\mu$ and by measuring the corresponding equilibrium densities. This method was adopted in Reference [13] to identify the $\mu =0$ line shown in Figure 1. Upon increasing $\mu$, the curve moves to the right in the $(a,h)$ diagram and becomes more and more vertical (data not shown). Isothermal lines $T=c$ can be found analogously: they roughly follow the profile of the $T=0$ line and, upon increasing c, they span all the positive-temperature region, approaching the infinite-temperature line for $c\to +\infty$ [13].

A non-equilibrium steady state can be represented as a path in the $(a,h)$-parameter space, where $a\left(x\right)$ is the mass density, $h\left(x\right)$ the energy density, and $x=n/N$ the rescaled position along the chain. In our set-up, the first site of the chain ($n=1$) is in contact with a reservoir at temperature $T_L$ and chemical potential ${\mu }_L$. This is ensured by implementing the non-conservative Monte-Carlo dynamics described in Reference [13]. In a few words, the reservoir performs random perturbations $\delta z_1=(\delta p_1+i\delta q_1)/\sqrt{2}$ of the state variable $z_1$ that are accepted or rejected according to a grand-canonical Metropolis cost-function $\exp [-T_L(\Delta H-{\mu }_L\Delta A)]$, where $\Delta H$ and $\Delta A$ are respectively the variations of energy and mass produced by $\delta z_1$. The perturbations $\delta p_1$ and $\delta q_1$ are independent random variables extracted from a uniform distribution in the interval $[-R,R]$. The opposite site ($n=N$) interacts with a pure stochastic dissipator that sets the variable $z_N$ equal to zero, thus absorbing an amount of mass equal to $|z_N{|}^2$. Both the heat bath and the dissipator are activated at random times, whose separations are independent and identically distributed variables uniformly distributed within the interval $[t_{min},t_{max}]$. On average, this corresponds to simulating an interaction process with decay rate $\gamma \sim {\overline{t}}^{-1}$, where $\overline{t}=(t_{max}+t_{min})/2$. Notice that different prescriptions, such as for example a Poissonian distribution of times with average $\overline{t}$ or a constant pace equal to $\overline{t}$, do not introduce any relevant modification in the dynamical and statistical properties of the model (1). Finally, the Hamiltonian dynamics between successive interactions with the environment has been generated by implementing a symplectic, 4th-order Yoshida algorithm [38]. We have verified that a time step $\Delta t=2\times {10}^{-2}$ suffices to ensure suitable accuracy.

Throughout the paper we deal with measurements of temperature profiles. Since the DNLS Hamiltonian is not separable, the standard relation between temperature and local kinetic energy does not hold. Due to the existence of the second conserved quantity A, it is necessary to make use of the microcanonical definition provided in [39] and further extended in [12]. Its derivation is founded on the thermodynamic relation

$$\beta =T^{-1}={\left.\frac{\partial S}{\partial E}\right|}_{A=M}$$

(4)

for a system with total energy $H=E$, total mass $A=M$ and entropy S. The calculation amounts to derive a measure of the hyper-surface at constant energy and mass in the phase space and to compute its variation with respect to E at constant mass. The general expression of T is a thermal average of a non-local function of the dynamical variables $z_n$ and $z_n^{\ast }$ and is rather involved; we refer to [13,21] and the related bibliography for theoretical and computational details.

In what follows, we consider a situation where all parameters, other than $T_L$, are kept fixed. In particular, we have chosen ${\mu }_L=0$, $R=0.8$ and $\overline{t}=3\times {10}^{-2}$, with $t_{max}$ and $t_{min}$ of order ${10}^{-2}$. We have verified that the results obtained for this choice of the parameter values are general. A more detailed account of the dependence of the results on the thermostat properties will be reported elsewhere.

Finally, we recall the observables that are typically used to characterize a steady-state out of equilibrium: the mass flux

$$j_a=2\langle \mathrm{Im}\left(z_n^{\ast }{z}_{n+1}\right)\rangle ,$$

(5)

and the energy flux

$$j_h=2\langle \mathrm{Re}\left({\dot{z}}_n{z}_{n+1}^{\ast }\right)\rangle ,$$

(6)

where the angular brackets denote a time average.

3. Low-Temperature Regime: Coupled Transport and Negative Temperatures

In the left panel of Figure 2, we report the average profile of the inverse temperature $\beta \left(x\right)$ as a function of the rescaled site position $x=n/N$, for different values of the temperature of the thermostat. A first “anomaly” is already noticeable for relatively small $T_L$: the profile is non monotonous (see for example the curve for $T_L=1$). This feature is frequently encountered when a second quantity, besides energy, is transported [13,40,41]. In the present setup, this second thermodynamic observable is the chemical potential $\mu \left(x\right)$, set equal to zero at the left edge. Rather than plotting $\mu \left(x\right)$, in the right panel of Figure 2, we have preferred to plot the more intuitive mass density $a\left(x\right)$. There we see that the profile for $T_L=1$ deviates substantially from a straight line, indicating that the thermodynamic properties vary significantly along the chain and suggesting that the lattice might not be long enough to ensure an asymptotic behavior.

To clarify this point, we have performed simulations for different values of N. The results for $T_L=1$ are reported in Figure 3, where we plot the local temperature T as a function of x. All profiles start from $T=1$, the value imposed by the thermostat and, after an intermediate bump, eventually attain very small values. Since neither the temperature nor the chemical potential are directly imposed by the purely dissipating “thermostat”, it is not obvious to predict the asymptotic value of the temperature (and the chemical potential). The data reported in the inset suggest a sort of logarithmic growth with N, but this is not entirely consistent with the results obtained for $T_L=3$ (see below).

If transport were normal and N were large enough, the various profiles should collapse onto each other, but this is far from the case displayed in Figure 3. The main reason for the lack of convergence is the growth of the temperature bump. This is because, upon further increases of N, the system spontaneously crosses the infinite temperature line and enters the negative-temperature region. For $T_L=3$, this “transition” has already occurred for $N=4095$, as it can be seen in Figure 2. The crossings of the infinite temperature points ($\beta =0$) at the boundaries of the negative temperature region correspond to infinite (negative) values of the chemical potential $\mu$ in such a way that the product $\beta \mu$ remains finite at these turning points (data not reported). To our knowledge, this is the first example of negative-temperature states robustly obtained and maintained in nonequilibrium conditions in a chain coupled with a single reservoir at positive temperature. In order to shed light on the thermodynamic limit, we have performed further simulations for different system sizes. In Figure 4, we report the results obtained for $T_L=3$ and N ranging from 511 to 4095. In Figure 4a, we see that the negative-temperature region is already entered for $N=1023$. Furthermore, its extension grows with N, suggesting that in the thermodynamic limit it would cover the entire profile but the edges.

Since non-extensive stationary profiles have been previously observed both in a DNLS and a rotor chain (i.e., the XY-model in $d=1$) at zero temperature and in the presence of chemical potential gradients [42], it is tempting to test to what extent an anomalous scaling (say $n/\sqrt{N}$) can account for the observations. In the inset of Figure 4a, we have rescaled the position along the lattice by $\sqrt{N}$. For relatively small but increasing values of $n/\sqrt{N}$ we do see a convergence towards an asymptotic curve, which smoothly crosses the $\beta =0$ line. This suggests that close to the left edge, positive temperatures extend over a range of order $\mathcal{O}\left(\sqrt{N}\right)$, thereby covering a non-extensive fraction of the chain length. The scaling behavior in the rest of the chain is less clear, but it is possibly a standard extensive dynamics characterized by a finite temperature on the right edge. A confirmation of the anomalous scaling in the left part of the chain is obtained by plotting the profiles of h and a again as a function of $n/\sqrt{N}$ (see the panels (b) and(c) in Figure 4, respectively). Further information can be extracted from the scaling behavior of the stationary mass flux $j_a$. In Figure 5, we report the average value of $j_a$ as a function of the lattice length. There, we see that $j_a$ decreases roughly as $N^{-1/2}$. At a first glance this might be interpreted as a signature of energy super-diffusion, but it is more likely due to the presence of a pure dissipation on the right edge (in analogy to what seen in the XY-model [43]).

In stationary conditions, mass, and energy fluxes are constant along the chain. This is not necessarily true for the heat flux, as it refers only to the incoherent component of the energy transported across the chain. More precisely, the heat flux is defined as $j_q\left(x\right)=j_h-j_a\mu \left(x\right)$ [43] Since $j_a$ and $j_h$ are constant, the profile of the heat flux $j_q$ is essentially the same of the $\mu$ profile (up to a linear transformation). In Figure 6a, we report the heat flux for $T_L=1$. It is similar to the temperature profile displayed in Figure 3. It is not a surprise to discover that $j_q$ is larger where the temperature is higher. Figure 6b in the same Figure 6 refers to $T_L=3$. A very strange shape is found: the flux does not only changes sign twice but exhibits a singular behavior in correspondence of the change of sign, as if a sink and source of heat were present in these two points, where the chemical potential and the local temperature diverge (see, e.g., the red dashed line in Figure 6b representing the $\beta$ profile). The scenario looks less awkward if the entropy flux $j_s=j_q/T$ is monitored. For $T_L=1$, the bump disappears and we are in the presence of a more “natural” shape (see Figure 6d) More important is to note that the singularities displayed by $j_q$ for $T_L=3$ are almost removed since they occur where $T\to \infty$ (we are convinced that the residual peaks are due to a non perfect identification of the singularities). If one removed the singular points, the profile of the entropy flux $j_s=j_q/T$ has a similar shape for the cases $T_L=1$ and $T_L=3$. A more detailed analysis of the scenario is, however, necessary in order to provide a solid physical interpretation.

4. High Temperature Regime: DB Dominated Transport

Let us now turn our attention to the high-temperature case. As shown in Figure 2a, for sufficiently large $T_L$ values, the positive-temperature region close to the dissipator disappears (this is already true for $T_L=5$) and, at the same time, the positive-temperature region on the left grows. In other words, negative temperatures are eventually restricted to a tiny region close to the dissipator side. This stationary state is induced by the spontaneous formation and the destruction of large DBs close to the dissipator. On average, such a process gives rise to locally steep amplitude profiles that are reminiscent of barriers raised close to the right edge of the chain, see Figure 2b. As it is well known, DBs are localized nonlinear excitations typical of the DNLS chain. Their phenomenology has been widely described in a series of papers where it has been shown that they emerge when energy is dissipated from the boundaries of a DNLS chain [21,22,23]. In fact, when pure dissipators act at both chain boundaries, the final state turns out to be an isolated DB embedded in an almost empty background. In view of its localized structure and the fast rotation, the DB is essentially uncoupled from the rest of the chain and, a fortiori, from the dissipators. One cannot exclude that a large fluctuation might eventually destroy the DB, but this would be an extremely rare event.

In the setup considered in this paper, DB formation is observed in spite of one of the two dissipators being replaced by a reservoir at finite temperature. DBs are spontaneously produced close to the dissipator edge only for sufficiently high values of $T_L$. Due to its intrinsic nonlinear character, this phenomenon cannot be described in terms of standard linear-response arguments. In particular, the temperature reported in the various figures cannot be interpreted as the temperature of specific local-equilibrium state: it is at best the average over the many different macrostates visited during the simulation. Actually, it is even possible for some of theses macrostates to deviate substantially from equilibrium. Therefore, we limit ourselves to some phenomenological remarks. The spontaneous formation of small breathers close to the right edge drastically reduces the dissipation and contributed to a further concentration of energy through the merging of the DBs into fewer larger ones and, eventually, to a single DB. Mechanisms of DB-merging have already been encountered under different conditions in the DNLS model [21]. The onset of a DB essentially decouples the left from the right regions of the chain. In particular, it strongly reduces the energy and mass currents fluxes. One can spot DBs simply by looking at the average mass profiles. In Figure 2, the presence of a DB is signaled by the sharp peak close to the right edge for both $T_L=9$ and $T_L=11$.

The region between the reservoir and the DB should, in principle, evolve towards an equilibrium state at temperature $T_L$. However, a close look at the $\beta$-profile in Figure 2 reveals the presence of a moderate temperature gradient that is typical of a stationary non-equilibrium state. In fact, DBs are not only born out of fluctuations, but can also collapse due to local energy fluctuations. As shown in Figure 7a, once a DB is formed, it tends to propagate towards the heat reservoir located at the opposite edge. The DB position is tagged by black dots drawn at fixed time intervals over a very long time lapse of $\mathcal{O}\left({10}^6\right)$, in natural time units of the model. This backward drift comes to a “sudden” end when a suitable energy fluctuation destroys the DB (see Figure 7a). Afterwards, mass and energy start flowing again towards the dissipator, until a new DB is spontaneously formed by a sufficiently large fluctuation close to the dissipator edge (the formation of the DB is signaled by the rightmost black dots in Figure 7a) and the conduction of mass and energy is inhibited again. The DB lifetime is rather stochastic, thus yielding a highly irregular evolution. The statistical properties of such birth/death process are discussed in the following section.

The statistical process describing the appearance/disappearance of the DB is a complex one. On the one hand, we are in the presence of a stationary regime: the mass and energy currents flowing through the dissipator are found to be constant, when averaged over time intervals much longer than the typical DB lifetime. On the other hand, strong fluctuations in the DB lifetime mean that this regime is not steady but it rather corresponds to a sequence of many different macrostates, some of which are likely to be far from equilibrium. Altogether, in this phase, the presence of long lasting DBs induces a substantial decrease of heat and mass conduction. This is clearly seen in Figure 5, where $j_a$ is plotted for different chain lengths. The two set of data corresponding to $T_L=9$ and 11 are at least one order of magnitude below those obtained in the low-temperature phase. The sharp crossover that separates the two conduction regimes is neatly highlighted in the inset of Figure 5, where the stationary mass flux is reported as a function of the reservoir temperature for the size $N=4095$.

The effect of the appearance and disappearance of the DB in the high reservoir temperature regime on the transport of heat and mass along the chain is twofold. During the fast dynamics, it produces bursts in the output fluxes of these quantities as demonstrated in Figure 7b. When the DB is present the boundary flux to the dissipator decreases, while when the DB disappears, avalanches of heat and mass reach the dissipator. In the slow dynamics obtained by averaging over many bursts, the conduction of heat and mass from the heat reservoir to the dissipator is hugely reduced with respect to the low temperature regime. We can conclude that the most important effect of the intermittent DB in the high temperature regime is to act as a thermal wall.

Finally, in Figure 6c,f we plot the heat and entropy profiles observed in the high-temperature phase, respectively. The strong fluctuations in the profiles are a consequence of the large fluctuations in the DB birth/death events and its motion. It is now necessary to average over much longer time scales to obtain sufficiently smooth profiles. It is interesting, however, to observe that the profile of $j_s$ in Figure 6f exhibits an overall shape similar to that observed in the low temperature regimes (see Figure 6d,e). This notwithstanding, there are two main differences with the low temperature behavior. First, close to the right edge of the dissipator we are now in presence of wild fluctuations of $j_s$, and second the overall scale of the entropy flux profile is heavily reduced.

5. Statistical Analysis

In order to gain information on the high-temperature regime, it is convenient to look at the fluctuations of the boundary mass flux $j_a^{\ast }$ and the boundary energy flux $j_h^{\ast }$ flowing through the dissipator edge. In Figure 8 we plot the distribution of both $j_a^{\ast }$ (panel (a)) and $j_h^{\ast }$ (panel (b)) for $T_L=11$ for different chain lengths. In both cases, power-law tails almost independent of N are clearly visible. This scenario is highly reminiscent of the avalanches occurring in sandpile models. In fact, one such analogy has been previously invoked in the context of DNLS dynamics to characterize the atom leakage from dissipative optical lattices [44].

We processed the time series of the type reported in Figure 7b to determine the duration ${\tau }_b$ of the bursts (avalanches) and ${\tau }_l$ of the “laminar” periods in between consecutive bursts (i.e., the DB life-times). In practice, we have first fixed a flux threshold ($s=4.25\times {10}^{-3}$) to distinguish between burst and laminar periods. Furthermore, a series of bursts separated by a time shorter than $d{t}_0={10}^3$ has been treated as a single burst. This algorithm has been applied to 20 independent realizations of the DNLS dynamics in the high-temperature regime. Each realization has been obtained by simulating a lattice with $N=511$ sites, $T_L=10$, and for a total integration time $t=5\times {10}^6$. In these conditions we have recorded nearly 7000 avalanches.

The probability distribution of the burst duration is plotted in Figure 9a. It follows a a Poissonian distribution, typical of random uncorrelated events. We have also calculated the amount of mass A and energy E associated with each burst, integrating mass, and energy fluxes during each burst. The results are shown as a scatter plot in the inset of Figure 9a. They display a clear (and unsurprising) correlation between these two quantities.

The time interval between two consecutive bursts is characterized by a small mass flux. Typically, during this period the chain develops a stable DB that inhibits the transfer of mass towards the dissipator. Figure 9b shows the probability distribution of the duration of these laminar periods. The distribution displays a stretched-exponential decay with a characteristic constant $\sigma =0.5$. Such a scenario is consistent with the results obtained in [26] for the Fermi-Pasta-Ulam chain and in [22] for the DNLS lattice. The values of the power $\sigma$ found in these papers are not far from the one that we have obtained here (see also [45] for similar results on rotor models).

Altogether, there is a clear indication that the statistics of the duration of avalanches and walls is controlled by substantially different mechanisms. It seems that the death of a wall/DB is ruled by rare event statistics [22,26,45], while its birth appears as a standard activation process emerging when an energy barrier is eventually overtaken. In fact, when mass starts flowing through the dissipator edge after the last death of a wall/DB, we have to wait for the spontaneous activation of a new wall/DB before the mass flux vanishes again. Conversely, the wall/DB is typically found to persist over much longer time scales and its eventual destruction is determined by a very rare fluctuation, whose amplitude is expected to be sufficiently large to compete with the energy that, in the meanwhile, has been collected by the wall/DB during its motion towards the reservoir.

6. Conclusions

We have investigated the behavior of a discrete nonlinear Schrödinger equation sandwiched between a heat reservoir and a mass/energy dissipator. Two different regimes have been identified upon changing the temperature $T_L$ of the heat reservoir, while keeping fixed the properties of the dissipator. For low $T_L$ and low chemical potential, a smooth $\beta$-profile is observed, which extends (in the central part) to negative temperatures without, however, being accompanied by the formation of discrete breathers over the time scales accessible to numerical simulations. In the light of the theoretical achievements by Rumpf [17,18,19,20], the negative-temperature regions are incompatible with the assumption of local thermodynamic equilibrium, which instead appears to be satisfied in the positive-temperature part. Therefore, despite the smoothness of the profiles and the stationarity of mass and energy fluxes, such negative-temperature configurations should be better considered as metastable states. Unfortunately, it is not easy to investigate this regime; no single heat bath can impose a fixed negative T in a meaningful way so as to be able to compare with equilibrium states. It is therefore necessary to simulate larger systems in the hope to observe the spontaneous formation of breathers, the only obvious signature of negative temperatures. We plan to undertake such a kind of numerical studies in the near future. It is nevertheless remarkable to see that negative temperatures are steadily sustained for moderately long chain lengths. As a second anomaly, we report the slow decrease of the mass-flux with the chain length: the hallmark of an unconventional type of transport. This feature is, however, not entirely new; a similar scenario has been previously observed in setups with dissipative boundary conditions and no fluctuations [42].

For larger temperatures $T_L$, we observe an intermittent regime characterized by the alternation of insulating and conducting states, triggered by the appearance/disappearance of discrete breathers. Note that this regime is rather unusual, since it is generated by increasing the amount of energy provided by the heat bath rather than by decreasing the chemical potential, as observed for example in the superfluid/Mott insulator transition in Bose-Einstein condensates in optical lattices. Although, for clarity, we referred to the low and high $T_L$ cases, there is no special difference to be expected in the intermediate regime, except for the typical timescale for DB creation that may be, still, exceedingly large. In fact, for finite N, such a timescale becomes much shorter above some typical $T_L$.

The intermittent presence of a DB/wall makes the chain to behave as a rarely leaking pipe, which releases mass droplets at random times when the DB disappears according to a stretched-exponential distribution. The resulting fluctuations of the fluxes suggest that the regime is stationary but not steady, i.e., locally the chain irregularly oscillates among different macroscopic states characterized, at best, by different values of the thermodynamic variables. A similar scenario is encountered in the XY chain, when both reservoirs are characterized by a purely dissipation accompanied by a deterministic forcing [42]. In such a setup, as discussed in Reference [42], the temperature in the middle of the chain fluctuates over macroscopic scales. Here, however, given the rapidity of changes induced by the DB dynamics, there may be no well-defined values of the thermodynamic observables. For example, during an avalanche it is unlikely that temperature and chemical potentials are well-defined quantities, as there may not even be a local equilibrium. This extremely anomalous behavior is likely to smear out in the thermodynamic limit, since the breather life-time does not probably increase with the system size, however, it is definitely clear that the associated fluctuations strongly affect moderately-long DNLS chains.