{'fraction_non_alphanumeric': 0.03309429765693907, 'fraction_numerical': 0.009510296429804043, 'mean_word_length': 4.680243955565237, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stretched--exponential relaxation in arrays of coupled rotators

We consider the non--equilibrium dynamics of a chain of classical rotators coupled at its edges to an external reservoir at zero temperature. We find that the energy is released in a strongly discontinuous fashion, with sudden jumps alternated with long stretches during which dissipation is extremely weak. The jumps mark the disappearance of strongly localized structures, akin to the rotobreather solutions of the Hamiltonian model, which act as insulating boundaries of a hot central core. As a result of this complex kinetics, the ensemble--averaged energy follows a stretched exponential law until a residual pseudo--stationary state is attained, where the hot core has reduced to a single localized object. We give a statistical description of the relaxation pathway and connect it to the properties of return periods of rare events in correlated time series. This approach sheds some light into the microscopic mechanism underlying the slow dynamics of the system. Finally, we show that the stretched exponential law remains unaltered in the presence of isotopic disorder.


Introduction
Many physical systems display unusual relaxation properties, such as slow kinetics in glasses [1] or dynamics of biomolecules [2], which are usually associated with the presence of a complex energy landscape. Long-living transient states also arise in non-linear systems, where the relaxation properties may be affected by the emergence of localized vibrations like solitons or kinks (see e.g. [3]). Of special interest is the role of a recently discovered class of spatially localized, time-periodic excitations, termed discrete breathers [4]. At variance with their counterparts in the continuum, they exist under very general conditions and their stability properties have been thoroughly investigated [4]. In particular, these objects may easily self-excite under very different physical conditions, as exemplified by several numerical and experimental studies [4,5,6,7,8,9,10].
A typical situation where breathers form spontaneously can be obtained by cooling the lattice boundaries, thus driving energy out of the system. In such conditions one is faced with a genuine nonlinear effect: the system prevents the complete dissipation of the energy by storing part of it in a residual longliving breather or multi-breather state. Such a phenomenon was first observed in a model of coupled harmonic oscillators pinned by an on-site nonlinear force [11,12] and later recovered in models of nonlinearly coupled oscillators without on-site interaction [13,14,15]. Nonetheless, the process of energy dissipation was found to exhibit quite different behaviors in these two classes of models. The main mechanism yielding a power-law decay of the energy to the residual state has been described in previous papers by the authors [13,14]. The stretched exponential behaviour, primarily observed in models with onsite interaction, was then conjectured to represent a typical feature for this class of models. Nonetheless, a satisfactory understanding of this peculiar decay process was still lacking.
The present paper deals with the nonequilibrium transient dynamics of a model of coupled classical rotators on a one-dimensional lattice. Two different cases will be considered, the one with pure nearest-neighbours coupling and the one with an additional on-site force which is a generalized version of the well-known discrete sine-Gordon chain. The quantum version of the latter model has been introduced by Fillaux and collaborators [16] to describe inelastic-neutron-scattering experiments on methyl-pyridine crystals. As this system displays evidence of extremely slow kinetics [17], it is particully interesting to investigate how cooling and nonlinear effects may affect thermalization. Although limited to the classical case, our study may provide some insight on the observed phenomenology.
The second motivation comes from the problem of stationary heat transport in this class of models. It has been recently recognized that low-dimensional lattices show anomalous properties, namely that transport coefficients (e.g. the thermal conductivity) diverge [18]. The analysis of several models clarified that this striking feature occurs generically whenever momentum is conserved. For lattice models, this amounts to requiring that at least one acoustic phonon branch is present in the harmonic limit. Remarkably, the only known exception is actually the pure nearest-neighbor model, which displays normal transport [19]. Since the transient dynamics is connected to the system's response, it is interesting to investigate the energy relaxation and to compare the peculiarities of each class of models. As we shall see, there are indeed some qualitative differences whose genuine nonlinear origin (excitation of long-lived localized structures) points at significantly differerent transport mechanism than those based on the customary phonon theory.
The plan of the paper is as follows. The details of the system of coupled classical rotators cooled at its boundaries are presented in Section 2: we consider both the cases with and without on-site interaction. The simulations reported in Section 3 indicate that the energy relaxation dynamics displays in both cases a stretched-exponential behavior. The statistical interpretation presented in Section 4 relates such a stretched-exponential behavior with rareevent statistics. To our knowledge, this correspondence has never been pointed out before in the literature. In section 5 we show that the main features and the interpretation do not change when isotopic disorder is introduced.
Two relevant results of our analysis merit to be stressed from the very beginning. The first one concerns the robustness of the phenomenon with respect to different variants of the model. The second one amounts to put such a phenomenon in a different perspective: the stretched-exponential energy decay process seems to be common to one-dimensional systems whenever the typical solutions are static breathers, rather than being a peculiarity of models with on-site interaction.

The model
We study the relaxation toward equilibrium of a chain of classical rotators coupled on a lattice with damping acting at its edges. The equations of motion are where φ i in the angle variable of the i-th rotator (with i = 1, 2, . . . N). The Hamiltonian version of (1) is sometimes referred to as the sine-lattice equations [20,21]. We impose free-end boundary conditions (φ 1 = φ 0 , φ N = φ N +1 ) since we know from previous studies that this choice favors spontaneous localization [13]. In the following we will also consider the special case G = 0, where only interactions between nearest neighbors are present and refer to it as the NN model. Furthermore, we will set I i = 1 for all i and K = 1.4911, a value that in our units corresponds to experimental data on 4-methyl-pyridine [16].
The general layout of a simulation has been detailed previously [13,14]. First, an equilibrium microstate is generated by letting the Hamiltonian (microcanonical, i.e. γ = 0) system evolve for a sufficiently long transient (typically 500 time units). We used the 5 th -order symplectic Runge-Kutta-Nyström algorithm [22]. The initial condition for the transient is assigned by setting all angles φ i to zero and by drawing velocities at random from a Gaussian distribution with standard deviation equal to the square root of twice the energy density E(0)/N. The resulting set of φ i andφ i is then used as initial condition to integrate Eq. (1) with γ > 0 with a standard 4 th -order Runge-Kutta algorithm. The subsequent dissipative dynamics is described in the following section.

The relaxation dynamics
A typical simulation is illustrated in Fig. 1. In the left panels we show the space-time contour plots of the symmetrized site energies e i , defined as In the right panels we report the time series of the total normalized energy E(t)/E(0) for the same simulation runs.
The energy density contour plot shows that the process is associated with the spontaneous creation of long-lived localized objects, characterized by a fast librating motion of a single rotator almost decoupled from its neighbors. These localized excitations (rotobreathers) have already been shown to exist in the Hamiltonian lattice [20]. The appearance of such strongly localized objects is due to the choice of the coupling which is of the same order as the on-site term (see the analysis in Ref. [20]).
The two outermost rotobreathers systematically act as barriers between a central "hot" region of the chain and its boundaries, thus blocking the energy flow toward the environment. At some stage, such dynamical barriers are spontaneously destabilized, thus allowing a portion of the trapped energy to rapidly flow away. This process goes on progressively reducing the size of the central hot patch until a single rotobreather survives in the bulk. We checked that its motion is periodic by computing the spatiotemporal spectrum S(k, ω) of the φ i s for several wave-numbers k. Indeed S(k, ω) displays sharp peaks at the fundamental frequencies plus small-amplitude harmonics.
The right panels of Fig. 1 illustrate how the energy flows toward the reservoirs.
In accordance with the above described phenomenology, the energy decay proceeds in a characteristic step-wise fashion, through a series of sudden jumps separated by long, approximately constant plateaus. The jumps correspond to the disappearance of one of the rotobreathers at the boundaries of the hot region. Between two consecutive jumps, the energy remains approximately constant with very small fluctuations.
Our simulations show that this kind of scenario persists also in the case G = 0, thus confirming previous qualitative results [23]. This finding is also illustrated in Fig. 2, where we compare the energy decay of four different initial conditions for the two kinds of potentials. In order to emphasize the slow character of the process, we rescale the time by τ 0 = N/2γ, which has been shown to be the characteristic time scale over which edge dissipation occurs [13].
It is interesting to observe that different initial conditions may give rise to substantially different sequences of jumps and plateaus, that may be regarded as different realizations of a stochastic process. For comparison, we remark that such sensitivity to the choice of the initial condition has not been observed in the relaxation of other non-linear systems, where different stories approximately resemble each other and display a much smoother decay [13,14].
In order to give a more quantitative description of the relaxation process, it is very interesting to study the average behavior of the energy decay. Figure 3 shows that the latter is characterized by three different scaling regimes. These are best visualized by plotting the quantity − ln E(t)/E(0) in loglog scale, where a law of the form exp(−t σ ) shows up as a straight line with slope σ. We find that both the systems with G = 0 and G = 1 exhibit the same properties. At short times (t ≪ τ 0 ), the energy decays exponentially as . This is what one observes for a harmonic lattice of the same size, thus signaling that the non-linear effects are negligible on such a time scale. Actually, the same holds for the Fermi-Pasta-Ulam chain [13]. At intermediate times (t ≃ 10 −2 τ 0 ), the curves undergo a change of slope and slow down to stretched exponential laws. The knee associated with such crossover may be due to power-law corrections, which are known to be present at this stage in other systems [13]. Eventually (t > 10τ 0 ), the residual quasistationary state is reached. In this last stage, a finite fraction of the initial energy is stored in a single rotobreather and remains constant over the rest of Table 1 Dependence of the exponent σ on the initial energy density for N = 100 and γ = 0.1. the simulation.
In Figure 3 we also consider the effects of the lattice size for the same value of the damping constant γ. In order to make the comparison, we rescale the time axis by the size-dependent constant τ 0 = N/2γ. We see that the data for N = 100, 200 collapse fairly well for t < τ 0 and follow approximately the same trend at later times. More marked differences are observed in the first crossover region. In particular, the range of stretched exponential decay seems to be larger for the longer chains.
To complete this analysis, we checked the dependence of the exponent σ on the initial energy. The data reported in Table 1 show that σ decreases by increasing E(0)/N. Accordingly, this suggests that the higher the initial energy density, the longer is the "lifetime" of the rotobreathers that appear in the system. There are some results in the literature that support this inference for the sine-lattice systems (with and without local coupling). In fact, it has been found that the lifetime of an excitation of energy E in a nonlinear system should be proportional to exp √ E [24].
Finally, we wish to mention that a similar stepwise decay process yielding an average stretched-exponential law has been found also in two-dimensional easy-plane ferromagnets [25]. In such case, annihilation of vortex-antivortex pairs is the basic physical mechanism.

Statistical analysis
In this section we attempt to understand better the origin of the slow kinetics. In particular, we shall focus on the statistics of plateaus and energy jumps during energy relaxation. As we already noticed above, the total energy of the lattice decreases in a stepwise fashion. Between subsequent jumps no appreciable energy exchange with the external world occurs and a plateau is observed (see again Fig. 2), whose duration may vary over different time scales. For a single realization of a given initial condition, we can thus approximate the curve E(t) as a piece-wise function through the sequences ∆E n (energy jump at the n-th step) and ∆t n (duration of the n-th plateau): where Θ(x) is the Heaviside step function. In the spirit of a statistical description we thus regard each sequence (∆E n , ∆t n ) for n = 1, . . . N p as a realization of a random process.
Since all the features of this process are unknown, we devised a numerical procedure to extract from the simulations both the ∆t n and ∆E n for each realization. The first task is to locate the beginning and end of each plateau to a reasonable level of accuracy. We chose to sample the whole time-series (with an adjustable sampling interval t s ) and look for the times where the energy differences J(t) = E(t + t s ) − E(t) exceed some pre-assigned threshold value |J(t)| > q. Physically, J(t) is the energy flux toward the external reservoirs integrated over a time t s .
In order to fix the value of the threshold, we computed the average number of plateaus N p as a function of q. We found that there is always a range of q values where N p displays a broad maximum. Therefore, we systematically chose the value of q in such a way as to minimize the variation of the number of steps. According to this prescription, our procedure identified approximately the same number of relaxation stages at each realization. The average value of the number of plateaus was found to depend primarily on the type of potential energy of the lattice, while displaying little sensitivity on the initial energy density. For example, with N = 160 the average value of the relaxation steps before reaching the quasi-stationary state in the NN lattices is N p ≈ 11 ± 1, while almost twice as much relaxation stages are needed in the presence of onsite term, N p ≈ 26±1. We infer that the latter system could be characterized by a much rougher energy landscape.
Since we do not know whether the process at hand is a stationary one, we first checked how the averages depend on time, i.e. on the index n. In Fig. 4 we see that the averages ∆t n are almost independent of the order of the relaxation n for both G = 0 and G = 0. At least to a first approximation, it is thus legitimate to assume that the process is a stationary one. It is important to stress that the above observation rules out the conjecture that the slow decay kinetics of the system may be the result of a progressive increase of the residence times in metastable states visited during relaxation.
The next step is to investigate the probability distributions of the lifetimes P(∆t) and energy jumps P(∆E). In Fig. 5 (a) we show the histograms of the energy jumps for the two cases G = 0 and G = 0. We found that the data can be fitted by the empirical normalized function Hence, the distribution of P(∆E) is exponentially small for small values of ∆E and decays following a power law of the kind 1/∆E α+1 for large values of the energy jumps. In the case shown in Fig. 5, we found α ≈ 1.6, ∆E 0 ≈ 0.0081 (G = 0), and α ≈ 1.8, ∆E 0 ≈ 0.0074 (G = 1).
We turn now to the distribution of plateaus P(∆t). In practice, it is a difficult task to construct the histograms from the measured time series, since the latter span a wide range (3 decades) yielding exceedingly large fluctuations. A way around this problem is to construct the cumulative (integrated) distributions C(∆t), that can be obtained from the numerical series by inverting the axes in a simple rank-size plot of the data. Remarkably, we found that the statistics of the time intervals follows an inherently slow stretched exponential distribution law. This finding is illustrated in Fig. 5 (b) for both the G = 0 and G = 0 cases, where we plot the quantity − log[C(∆t)]. More precisely, we see that the first part of the distribution is a pure exponential (straight lines with slope one in Fig. 5 (b)), while the tails follow a stretched exponential law. This means that the probability distribution P(∆t) (the derivative of C(∆t)) is up to power-law corrections. The stretched exponential portion of the curves sets in at times longer than ∆t * ≈ 10 3 , which is systematically below the average, approximately constant duration of each relaxation stage. This ob- servation suggests that the statistics of the plateau duration is dominated by the stretched exponential tail. Such intrinsically fat-tailed distribution of a series of relaxation events is thus directly related to the slow energy decay.
In general, the origin of a stretched exponential relaxation is traced back to a superposition of many different time scales, distributed according to given weights [12,26]. In our context, however, it is not clear how to identify such a hierarchical set of degrees of freedom, that would be responsible for the observed slow kinetics. It is thus legitimate to explore the possibility of alternative explanations. An interesting perspective has been presented by Bunde et al. in a recent paper [27], in connection with the statistics of extreme events. According to the authors, a stretched exponential law may rule the distribution of intervals between consecutive rare events provided the original time series displays power-law correlation. More precisely, they consider a stochastic time series x(t), whose autocorrelation decays by construction as t −γ with 0 < γ < 1. In this context, rare events are defined by the condition x(t) > q, where q is a pre-assigned threshold. Their numerical analysis shows that the distribution of the time intervals ∆t separating extreme events decays as exp[−∆t γ ′ ], with γ ′ ≈ γ.
It appears natural to extend this kind of description to the process of energy release in our system: each energy step ∆E, marking the collapse of a rotobreather, can be regarded as a "rare event " (sharp peak) in the series of energy flux J(t) (see Fig. 2). In order to validate this analogy, we checked that longrange correlations do exist in the series J(t) (Fig. 6). Within the statistical fluctuations, the autocorrelation function of J(t) decreases approximately as an inverse power law, whose exponent is smaller than one (≈ 0.7). This value is of the same order of the exponents reported in Fig. 5 (b). As a matter of fact, the statistical fluctuations prevent us from obtaining a better estimate. Furthermore, the distribution of plateau durations are only asymptotically stretched-exponential, thus introducing a further source of error. Therefore, to make our analogy more quantitative further numerical work should be carried out.
We note that a similar analysis of the correlations in the series of the energy flux in the Fermi-Pasta-Ulam chain does not yield a power-law decay. In fact, the autocorrelation function of J(t) vanishes exponentially. This is consistent with the absence of stretched exponential relaxation in this model [13], thus suggesting that the microscopic mechanisms operating in the chain of coupled rotators are substantially different.
In order to provide further support to our interpretation, we have set up a Monte Carlo procedure which generates stochastic series of N p energy jumps and plateaus as in Eq. (3) according to given distributions P(∆E) and P(∆t).
We have then computed the average over such a set of realizations according to Eq. (3), where N p can be fixed close to its empirical average value. We found that, by using a distribution of the form (4) and a stretched exponential law for the P(∆t), we could recover a stretched exponential decay of the average energy by just tuning at most two of the free fitting parameters ∆E 0 , α, τ 2 , β. For example, by fixing N p = 15, ∆E 0 = 0.0081 and α = 1.7 in the case G = 0, E(0)/N = 8, the best-fit value of the exponent σ from the Monte Carlo simulation (σ = 0.42) is found to be in very good agreement with the one directly fitted on the simulation data (see Table 1).

The effect of isotopic disorder
In this section we investigate the role of disorder on the relaxation dynamics of the sine-lattice system. More specifically, we consider model (1) with randomly distributed momenta of inertia, that is either I i = R (with probability ρ) or I i = 1 (with probability 1−ρ). This choice corresponds to a 4-methyl-pyridine quasi-1D chain with random insertions of deuterated methyl groups, that has been studied experimentally by Fillaux and co-workers [16,28]. Hence, we fix R = 2, that corresponds to the ratio of momenta of inertia reported in Ref. [16].
The motivation of this analysis is twofold. First, the joint effect of non-linearity and disorder on the non-equilibrium properties of discrete lattices are far from being satisfactorily understood [29]. Second, the disordered sine-lattice model has been used to interpret data from inelastic neutron scattering ex- periments [16]. More specifically, it has been shown that the effects of isotopic disorder on the dynamics rapidly disappear upon increasing the temperature. This observation has been explained quantum-mechanically as a transition to the classical regime, thus suggesting that the effects of disorder should be negligible in the latter case. By the same token, one may expect that the relaxation dynamics in the classical system should display weak sensitivity on the concentration of impurities.
Our simulations confirm the above hypothesis. In Fig. 7 we show that the energy decay is hardly affected by the value of ρ. Remarkably, the curves almost overlap and display both a wide region of stretched exponential behavior and the saturation to a residual state. Accordingly, the estimated exponents are fairly insensitive to the concentration of impurities. Moreover, we checked that the energy which is left in the residual state is approximately independent of ρ, with deviations of at most 10 % from the average.

Conclusions
The study of the energy relaxation in models of classical rotators has allowed us to identify a possible mechanism yielding the stretched exponential behavior, observed also in other one-dimensional models [11,12,29]. The static nature of localized breather states, which are spontaneously formed upon cooling the lattice from its boundaries, is one of the basic ingredients. Indeed, static localized solutions can form "dynamical barriers" that are almost decoupled from their neighbours and thus segregate an internal hot lattice region from the cooled boundaries. The energy initially fed into the system remains trapped there, until the rotobreathers are abruptly destroyed. It may be argued that this sudden death must originate from a sufficiently large resonant fluctuation emerging from the hot core. In our approach, we have employed a statistical description, which does not make any explicit reference to specific dynamical mechanisms.
Accordingly, we showed that the destruction of boundary rotobreathers may be regarded as a "rare event". The resulting statistics depends on two main factors: the energy of the rotobreathers and the probability that in a given time interval the hot region may produce an excitation capable of annihilating them.
In our numerical analysis we can just observe the effects of this complicated interaction, which results in a stretched-exponential law for the energy decay.
It may be argued that the mechanism of energy trapping must be rather peculiar to the 1D nature of the model. Indeed, we expect that in higherdimensional lattices this phenomenon may not occur at all, and that the characteristic stepwise decay of energy could be replaced by a smoother decrease. Some preliminary simulations of a 2D arrays of rotators confirm this idea.
In particular, we observe that the residual state resembles the one of the 2D Fermi-Pasta-Ulam model -a collection of static breathers randomly arranged on the lattice [14].
Altogether, we have seen that the possibility of easily exciting long-lived rotobreathers is responsible for the slow kinetics. It should be recalled that these very same objects have been invoked to explain the normal conduction in the NN model [19]. In this sense, our results complement those observations and confirm the essential role of nonlinear excitations in the out-of-equilibrium properties of many-body systems.