A Chain, a Bath, a Sink and a Wall

We investigate out-of-equilibrium stationary processes emerging in a Discrete Nonlinear Schroedinger chain in contact with a heat reservoir (a bath) at temperature $T_L$ and a pure dissipator (a sink) acting on opposite edges. 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 events of the breathers. The corresponding statistics exhibits 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.


Introduction
The study of nonequilibrium thermodynamics of systems composed of a relatively small number of particles is motivated by the need of a deeper theoretical understanding of 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 a layered photonic or phononic systems [8,9]. More recently, a renewed interest for this multipurpose equation emerged in the physics of gases of ultra-cold atoms trapped in optical lattices (e.g., see Ref. [10] and references therein for a recent survey). Since the seminal paper by Rasmussen et al. [11] it was realized that the presence in the DNLS of intrinsically localized solutions, named discrete breathers (DB) (see e.g. [12][13][14]), could be associated to negative absolute temperature states. In a series of important papers, Rumpf provided entropy-based arguments to describe asymptotic states above a modulational instability line in the DNLS [15][16][17][18]. It has been later found that above this line, negative temperature states can form spontaneously via the dynamics of the DNLS. They persist over extremely long time scales, that might grow exponentially with the system size as a result of an effective mechanism of ergodicity-breaking [19]. It has been recognized that most of these peculiar features of the DNLS can be traced back to the properties of its Hamiltonian that admits two first integrals of motion: the total energy H and the total number A of particles (also termed as the total norm).
A related question is how the structure of the Hamiltonian influences non-equilibrium properties when the system can exchange energy and/or mass with the environment. In a series of papers [19][20][21] it has been found that, when pure dissipators act at both edges of a DNLS chain (a case sometimes called boundary cooling [22][23][24][25]), the typical final state consists in 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 [26,27]. A similar phenomenology occurs after a quench from high to low temperatures in oscillator lattices [28]. Also, boundary driving by external forces may induce non-linear localization [29][30][31].
When, instead, the chain is put in contact with thermal baths at its edges, stationary states characterized by a gradient of temperature and chemical potential emerge [32]. The 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 [33], or persistent currents [34] are found, as well as a signature of anomalous transport in the low-temperature regime [35,36].
In this paper we consider a setup where one edge is in contact with a heat reservoir at temperature T L and chemical potential µ 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 setups [1,2], this is conceptually closer to a semi-infinite array in contact with a single reservoir. In fact, on the pure-dissipator side, energy can only flow out of the system.
A non-equilibrium steady state can be conveniently represented as a path in the (a − h)-parameter space, where a(x) is the mass density, h(x) the energy density, and x ∈ [0, 1] the rescaled position along the chain (see Fig. 1 for a few sampled paths). Making use of suitable microcanonical definitions [32], these paths can be converted into temperature (T) and chemical potential (µ) profiles. The presence of a pure dissipator in x = 1 forces the corresponding path to terminate close to the point (a = 0, h = 0), which is singular both in T and µ. Therefore, slight deviations may easily lead to crossing the β = 0 line, where β is the inverse temperature. This is indeed the typical scenario observed for small T L , when β smoothly changes sign twice before approaching the dissipator. The size of the negative-temperature region increases with T L . For high temperatures, a different stationary regime is found, characterized by strong fluctuations of mass and energy flux. In fact, upon increasing T L , 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 is controlled by the spontaneous formation (birth) and disappearance (death) of discrete breathers.
In Section 2 we introduce the model and briefly recall the main observables. Section 3 is devoted to a detailed characterization of the low-temperature phase, while the strongly non-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.

Model and observables
We consider a DNLS chain of size N and with open boundary conditions, whose bulk dynamics is ruled by the equation where (n = 1, . . . , N), z n = (p n + iq n )/ √ 2 are complex variables, with q n and p n being standard conjugate canonical variables. The quantity a n = |z n (t)| 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 iz * n , Eq. (1) can be read as the equation of motion generated by the Hamiltonian functional through the Hamilton equationsż n = −∂H/∂(iz * n ). We are dealing with a dimensionless version of the DNLS equation: the nonlinear coupling constant and the hopping parameters, that usually are 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 where we can easily check that H and the total mass A = ∑ n |z n | 2 (t) are conserved dynamical quantities.
The first site of the chain (n = 1) is in contact with a reservoir at temperature T L and chemical potential µ L . This is ensured by implementing the non-conservative Monte-Carlo dynamics described in Ref. [32]. The opposite site (n = N) interacts with a pure stochastic dissipator: the variable z N is set equal to zero 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 a dissipation process with decay rate γ ∼t −1 , wheret = (t max + t min )/2. Notice that different prescriptions, such as for example a Poissonian distribution of times with averaget, or a constant pace equal tot, 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 thermostats has been generated by implementing a symplectic, 4th-order Yoshida algorithm [37]. We have verified that a time step ∆t = 2 × 10 −2 suffices to ensure suitable accuracy.
Throughout the paper we deal with measurements of local temperature and chemical potential. Since the Hamiltonian is non separable, it is necessary to make use of the microcanonical definition provided in [38]. The general expressions are nonlocal and rather involved; we refer to [32] for details and the related bibliography.
In what follows, we consider a situation where all parameters, other than T L , are kept fixed. In particular, we have chosen µ L = 0 andt = 3 × 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 Im(z * n z n−1 ) , and the energy flux

Low-temperature regime: coupled transport and negative temperatures
In the left panel of Fig. 2 we report the average profile of the inverse temperature β(x) 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 [32,39,40]. In the present setup, this second thermodynamic observable is the chemical potential µ(x), set equal to zero at the left edge. Rather than plotting µ(x), in the right panel of Fig. 2, we have preferred to plot the more intuitive mass density a(x). There we see that the profile for T L = 1 deviates substantially from a straight line, 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 Fig. 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 Fig. 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 Fig. 2. The crossings of the infinite temperature points (β = 0) at the boundaries of the negative temperature region correspond to infinite (negative) values of the chemical potential µ in such a way that the product βµ 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 Fig. 4 we report the results obtained for T L = 3 and N ranging from 511 to 4095. In panel (a) 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 [41], it is tempting to test to what extent an anomalous scaling (say n/ √ N) can account for the observations. In the inset of Fig. 4(a), we have rescaled the position along the lattice by √ N. For relatively small but increasing values of n/ √ N we do see a convergence towards an asymptotic curve, which smoothly crosses the β = 0 line. This suggests that close to the left edge, positive temperatures extend over a range of order O( √ N), 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/ √ N (see the panels (b) and(c) in Fig. 4, respectively).  Further information can be extracted from the scaling behavior of the stationary mass flux j a . In Fig. 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 [42]).
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 (x) = j h − j a µ(x) [42] Since j a and j h are constant, the profile of the heat flux j q is essentially the same of the µ profile (up to a linear transformation). In Fig. 6(a) we report the heat flux for T L = 1. It is similar to the temperature profile displayed in Fig. 3. It is not a surprise to discover that j q is larger where the temperature is higher. Panel (b) in the same Fig. 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 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 panel (b) representing the β 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 Fig. 6(d)) More important is to note that the singularities displayed by j q for T L = 3 are almost removed since they occur where T → ∞ (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

High temperature regime: DB dominated transport
Let us now turn our attention to the high-temperature case. As shown in Fig. 2(a), 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 destruction of large DBs close to the dissipator. On average, such process gives rise to locally steep amplitude profiles that are reminiscent of barriers raised close to the right edge of the chain, see Fig. 2(b). 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 [19][20][21]. 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 [19]. 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 Fig. 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 β-profile in Fig. 2  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 Fig. 7 (a), 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 O(10 6 ), in natural time units of the model. This backward drift comes to a "sudden" end when a suitable energy fluctuation destroys the DB (see Fig. 7(a)). 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 Fig. 7(a)) 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, the 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 Fig. 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 Fig. 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 burts in the output fluxes of these quantities as demonstrated in Fig. 7(b). 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 Fig. 6 (panels (c) and (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 Fig. 6 (f) exhibits an overall shape similar to that observed in the low temperature regimes (see Fig. 6 panels (d) and (e)). This notwithstanding, there are two main differences with the low temperature behaviour. 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.

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 and the boundary energy flux j * h flowing through the dissipator edge. In Fig. 8 we plot the distribution of both j * a (panel (a)) and j * h (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 [43].
We processed the time series of the type reported in Fig. 7 (b) to determine the duration τ b of the bursts (avalanches) and τ 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 × 10 −3 ) to distinguish between burst and laminar periods. Furthermore, a series of bursts separated by a time shorter than dt 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 × 10 6 . In these conditions we have recorded nearly 7000 avalanches. The probability distribution of the burst duration is plotted in Fig. 9(a). It follows a a Poissonian distribution, typical of random uncorrelated events. We have also calculated the amount of mass A and energy E associated to each burst, integrating mass and energy fluxes during each burst. The results are shown as a scatter plot in the inset of Fig. 9(a). 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. Fig. 9 (b) shows the probability distribution of the duration of these laminar periods. The distribution displays a stretched-exponential decay with a characteristic constant σ = 0.5. Such a scenario is consistent with the results obtained in [24] for the FPU chain and in [20] for the DNLS lattice. The values of the power σ found in these papers are not far from the one that we have obtained here (see also [44] 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 [20,24,44], 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.

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 β-profile is observed, which extends (in the central part) to negative temperatures, without, however, being accompanied by the formation of discrete breathers. In the light of the theoretical achievements by Rumpf [15][16][17][18], we can say that such regions at negative temperature are certainly incompatible with the assumption of local thermodynamic equilibrium, usually invoked in linear response theory of transport processes. In this sense, 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. Lacking any numerical evidence of their final evolution even for moderate chain lengths, we can just speculate that they could last over astronomically long times. Further studies are necessary to clarify whether and how this regime can persist in the thermodynamic limit. 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 [41].
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. 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 [41]. In such a setup, as discussed in Ref. [41], 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.