Thermalization in the one-dimensional Salerno model lattice

The Salerno model constitutes an intriguing interpolation between the integrable Ablowitz-Ladik (AL) model and the more standard (non-integrable) discrete nonlinear Schr{\"o}dinger (DNLS) one. The competition of local on-site nonlinearity and nonlinear dispersion governs the thermalization of this model. Here, we investigate the statistical mechanics of the Salerno one-dimensional lattice model in the nonintegrable case and illustrate the thermalization in the Gibbs regime. As the parameter interpolating between the two limits (from DNLS towards AL) is varied, the region in the space of initial energy and norm-densities leading to thermalization expands. The thermalization in the non-Gibbs regime heavily depends on the finite system size; we explore this feature via direct numerical computations for different parametric regimes.


I. INTRODUCTION
Enlightening of the thermalization properties of lattice dynamical models and complex networks is a crucial issue in understanding and exploiting the transport and localization phenomena of relevance to a wide range of physical problems. In spite of substantial research, the coexistence of diverse physical processes and correlations among them developing on different spacetime scales, lacks a conclusive interpretation [1][2][3][4][5][6][7][8][9][10][11]. It is therefore natural to expect that the application of statistical and thermodynamical approaches and related mixing, ergodicity, and energy equipartition concepts to such problems is a topic of substantial ongoing interest. Among the numerous prototypical nonlinear physical examples are: the statistics of the Fermi-Pasta-Ulam-Tsingou chains [11], chains of Josephson-junctions [12], Gross-Pitaevskii/discrete nonlinear Schrödinger lattices with different types of nonlinearities [13][14][15][16][17][18][19], Toda and Morse lattices [20], etc. The localized wave patterns that naturally emerge in such systems as a result of the interplay between lattice dispersion and nonlinearity play a crucial role in the thermalization and lattice dynamics [21].
Our study is partially motivated by recent findings regarding the statistics of different discrete nonlinear physical systems [12,19] in which, as a core mechanism, the relaxation of nonlinear localized excitations and related ergodization is considered. In the many-body systems, the ergodization demands infinite-time averages of an observable during a microcanonical evolution to match with their proper phase space averages. In a class of dynamical systems where ergodization timescales sensitively depend on the control parameters, dynamical glass behavior is postulated to be a generic system property on the route towards the integrable limits [12]. This glassy behavior is further attributed to the short range network in the action space. Hence it is interesting to explore the thermalization in a system where both an integrable a non-integrable (yet physically relevant) limits can exist as a suitable parameter is varied. One such model is the Salerno model (SM) [22]. On the one hand, and despite the two decades that have ensued since the attempt to thermodynamically describe the discrete nonlinear Schrödinger (DNLS) model [23,24], the problem continues to attract significant attention; see, Refs. [25][26][27] for recent studies. On the other hand, the SM has provided an excellent platform for exploring the interplay between nonlinear localized structures and near-linear extended ones, between nonlinearity and dispersion, on the path between integrability and non-integrability; see, e.g., [28] for a recent review.
Bearing the above features in mind, our aim here is to explore the statistical mechanics and thermalization properties of the Salerno model. In addition to interpolating between the fully integrable Ablowitz-Ladik (AL) and the DNLS models [6,22,[29][30][31], the SM incorporates coexistence and competition of nonlinear dispersion and nonlinear local interactions. In the present work, we adopt a grand-canonical description of SM, decomposing the parameter space of energy and norm densities (corresponding to the conserved quantities of the energy and total norm, respectively) into Gibbs and non-Gibbs regimes. In the former, we expect "regular" thermalization. In the latter non-Gibbs regime, we expect to encounter energy localization in the form of long-living nonlinear excitations in line with the corresponding DNLS prediction of [23]. One part of the intriguing story is that as the SM parameters are varied homotopically interpolating between the DNLS and the AL models, we should be able to observe that regions that used to belong in the non-Gibbs regime will now "regularly" thermalize as we approach the AL limit. That is, we should be able to observe a key change in behavior parametrically along the paths of parametric variation considered. However, there are also additional features that add to the complexity of the story.
While the general expectation previously was that the non-Gibbs regime is non-ergodic, a recent study of the DNLS lattice using the statistics of excursion times of equilibrium Poincaré manifolds and finite time average distributions of an observable has shown that a part of the non-Gibbs regime is weakly ergodic [11,12,19] in the setting of finite lattices. The non-ergodicity may be a feature of the thermodynamic limit of infinite lattices. Motivated by this array of recent developments and challenging observations, we numerically investigate the thermalization in the SM with respect to both the Gibbs and non-Gibbs regimes. We will utilize a combination of thermodynamic tools, including the transfer integral approach, and direct numerical simulations, and measures of both statistical properties (e.g., the probability distribution function (PDF) of different amplitudes) and dynamical properties (such as Lyapunov exponents) in order to characterize the system, including for different lattice sizes, so as to address the fundamental question of the behavior of the SM under parametric, initial condition and lattice size variations.
Our presentation is structured as follows. In section II, we present the fundamentals of the model. In section III, we lay the theoretical foundations for the statistical mechanical analysis of the SM. In section IV, we present the corresponding numerical analysis (both through statistical and dynamical diagnostics) and finally in section V we summarize our findings, and present our conclusions, as well as some directions of future research. The Appendix presents details of our theoretical analysis.

II. MODEL DESCRIPTION
The SM can be considered as a dicretization of the continuous fully integrable nonlinear cubic Schrödinger equation. Its equations of motion [22,30,31], upon suitable rescaling of the inter-site coupling, read: where ψ n is the complex wave function at site n, γ = 2(1 − µ) and µ ≥ 0. The parameter γ represents the strength of local nonlinear interaction and µ represents the strength of nonlinear dispersion (inter-site nonlinearity). In the limit µ → 0 the model reduces to the standard DNLS equation with on-site (local) cubic nonlinearity. On the contrary, the limit γ = 0 corresponds to the completely integrable AL model [24,32,33]. A simple transformation ψ n → ψ n e i2t leads to The full set of invariants of motion in the completely integrable, AL limit is considered in [34], while the nonintegrable DNLS limit is characterized by two integrals of motion. Therefore, regardless of the limits, Eq. (2) can be characterized by two conserved quantities: norm A and Hamiltonian H [35]: where N is the total number of lattice nodes and periodic boundary conditions are used. The SM equations (Eqs. (2)) can be derived from the with respect to the canonically conjugated pairs of variables ψ n and iψ * n defining the deformed Poisson brackets [36]

III. STATISTICAL MECHANICS OF THE SALERNO NETWORK
Here we attempt to clarify the thermalization properties of the SM starting from the DNLS limit µ = 0, in which the thermalization and statistical properties are extensively investigated [2,19,23,36,37]. After a brief remark on findings in the DNLS limit we probe the extension of the Gibbs approach to the SM with competing local and nonlocal nonlinearities.
Applying the canonical transformation ψ n = √ A n exp (iφ n ), where A n and φ n denote the amplitude and phase, we obtain from Eq. (3) the following expressions for the conserved quantities: The corresponding grand-canonical partition function of the SM can then be presented in a form where parameters α and β are introduced in analogy with the chemical potential and the inverse temperature [23] (i.e., they are the corresponding Lagrange multipliers). This expression can be reduced to the integral form after the integration over the phase variables φ n : where I 0 stands for the modified Bessel function of the first kind (with index 0). In the thermodynamic limit of large systems N → ∞, the integral can be evaluated exactly using the eigenfunctions and eigenvalues of the transfer integral operator (TIO) (Appendix A). From the latter calculation, in the infinite temperature limit β → 0, the following relation between the energy (h = H/N ) and norm (a = A/N ) density can be derived which for µ = 0 reduces to the corresponding relation of the DNLS lattice [6,19,23]. In terms of the largest eigenvalue λ 0 of the kernel Eq. (A10), the norm density, a can be expressed as a = ∞ 0 y 2 0 (A)AdA, where y 2 0 (A) represents the probability distribution of amplitudes P (A) corresponding to the largest eigenvalue.
Following the statistical mechanical analysis of the Appendix, we distinguish the Gibbs regime in the (a, h) parameter space by determining the characteristic phase curves β → 0 and β → ∞. While formally the first one separates the microcanonically inaccessible regime from the Gibbs region of phase space, the second one separates regions characterized by positive temperature ( 1 β > 0) from those with negative temperature ( 1 β < 0) whose accessibility is experimentally and numerically proven and which leads to the prolonged emergence of coherent structures.
The transition curve β → ∞ can be determined from by minimizing the Hamiltonian Eq. (3) with the planewave solution in a form ψ n = √ de inθ and taking θ = π. To clarify the thermalization properties we calculate the amplitude probability density function (P (A)) and the excursion time probability (P + ). The P (A) obtained from a direct numerical simulation of Eq. (1) is compared with the one calculated via the TIO approach [23]. On the other hand, the time intervals which the local norms spend between two consecutive intersections of the plane A n = a (the Poincaré section) form the excursion time distribution P + (τ ), where + denotes A n > a and τ n (i) = t i+1 n − t i n . The distribution has the average µ τ and the standard deviation σ τ . The value of σ τ can be associated with the divergence of the average excursion time and weak nonergodicity in the lattice [11,19]. A third measure of thermalization is the finite time average (FTA) of the observable, The distribution of the FTA for a set of trajectories is characterized by the first moment m 1 (T ) and the second moment m 2 (T ). For an ergodic regime, at large time m 2 (T → ∞) → 0 [12,38].
As an additional diagnostic, we estimate the maximal Lyapunov Characteristic Exponent (mLCE) Λ [39][40][41] which is in general a measure of the degree of chaos in the system. More concretely, we derive the evolution equations for small perturbations χ n (0) of the initially injected plane-wave profile ψ n (0) adopting the standard procedure based on the linearization in the presence of small perturbations and numerically solve the obtained variational equations The mLCE is then obtained as with λ(t) denoting the so-called finite time mLCE (ftmLCE), χ χ χ(t) = (χ 1 , χ 2 , ..., χ N ) being the deviation vector and || · || the usual Euclidean norm. In order to investigate thermalization in the SM we perform numerical experiments and base our corresponding analysis on the system's phase diagram illustrated in Fig. 1, in the parameter space (a, h). The red curve corresponds to the zero temperature limit β → ∞, while black curves correspond to the infinite temperature limit, β → 0 for a few values of parameter µ in the interval 0 to 0.5. These are based on the analytical predictions of these limits given above (and derived in the Appendix). The region between β = ∞ and β = 0 lines denotes the Gibbs regime of the SM where the approach based on the grand-canonical-Gibbs statistics (Eq. (7)) is applicable [19,23]) and thus the model is expected to thermalize. Outside this region, a non-Gibbs regime featuring the presence of coherent structures is identified where, however, for finite domains only weak non-ergodicity may be present as discussed in [19]. We now proceed to analyze our numerical computations at different selected points within this parameter space bearing in mind that a key feature of the SM case is that regimes identified as non-Gibbsian for the DNLS model of µ = 0 can be Gibbsian for larger values of µ.

IV. NUMERICAL ANALYSIS
In this section, we describe the numerically obtained results for various µ values corresponding to both the Gibbs and non-Gibbs regimes. This section aims to explore the thermalization properties in these regimes induced by the competing nonlinearities, local nonlinear interaction and nonlinear hopping, which cause selftrapping and nonlinear dispersion, respectively. We solve Eq.

A. Gibbs regime
We consider the parameter set (a = 1.5, h = 3) which is in the Gibbs regime irrespective of the µ value considered (see Fig. 1). The numerically calculated PDFs, P (A) (solid curves in Fig. 2) show that the amplitude A increases with the increase of µ. In order to ensure that the obtained P (A) represents a thermalized state, we compare them with the probability distributions obtained by the transfer integral operator (TIO) approach based on the corresponding dominant (squared) eigenvector of the TIO approach (see, e.g., Appendix I and also [23]). Our results indicate that TIO solutions (dashed curves in Fig. 2) and numerical results match very well, which corroborates the anticipated thermalization in this regime.
The amplitude profiles corresponding to the case (a = 1.5, h = 3) for three different µ values are shown in Fig. 3. Though high amplitude nonlinear localized excitations emerge in the system as a result of the modulational instability of the initial condition, they are all rather short lived and the long time evolution of the system appears to be thermalized into a phononic bath and without evidence of any kind of persistent localization. The thermalization is further verified from the second moment, m 2 (T ) of the FTA of local integral norms (Fig. 4), which decays as 1/T at large times for all values of the considered µ. Additionally, the calculation of the mLCE shows that Λ > 0 (red points in Fig. 5) which is a signature of chaoticity in the system. Interestingly, we find that the latter grows exponentially as a function of µ.

B. Non-Gibbs regime
To investigate the thermalization in the non-Gibbs regime, we consider the parameter set (a = 1.5, h = 5), the blue colored point in Fig. 1. For this parameter set, the critical value µ c ≈ 0.17 sets the transition point from Gibbs to non-Gibbs regime. That is, for µ > µ c , (a = 1.5, h = 5) is in the Gibbs regime, where TIO solutions match exactly with the numerically calculated P (A) as shown in Fig. 6.
For µ < µ c , the tail of P (A) develops a bump (at suitably large amplitudes) at initial times, which can be associated with the accumulation of large amplitude nonlinear excitations. In spite of it, we observed an exponential cutoff that might be related to the finite size of the system, which has been shown to affect all statistical measures we calculated in the non-Gibbs regime. It is worthwhile to observe that this bump is no longer present in the cases of µ > µ c (while the cutoff is still featured)  and at large times for µ < µ c . To corroborate this observation in the system, but also to explore the role of the finite size effects, we first plot the amplitude profiles as a function of time for three different values of µ, Fig. 7.
Interestingly, for all values of µ, both those that belong to the Gibbs regime (µ = 0.5 and µ = 0.2 of panels (c) and (b); the latter is near and above the critical value) and even in the non-Gibbs regime (µ = 0.05), the decay of initially generated large amplitude structures is eventually observed. Nevertheless, the coherent structures are significantly more prominent in the non-Gibbs case. The latter case of panel (a) represents the possibility of a non-Gibbsian regime which, however, within a finite lattice manifests the features of quasi-ergodic behavior as has been discussed in [19]. The second moment of the FTA of the integral local norms, shown in Fig. 8 Gibbs regimes (µ < µ c ). A similar result (in the sense of not distinguishing between the Gibbs and non-Gibbs regimes) is also found for the mLCEs in Fig. 5.
In line with recent developments in the field [12,19], it is quite relevant to explore the effects of the finite size of our computation and their implications in connection with the infinite domain thermodynamic analysis, e.g., implicit in the TIO calculations. To explore the role of finite system size, we calculate the variance, σ 2 τ of the excursion time distribution for different system sizes as depicted in Fig. 9. Upon increasing the system size N , the variance σ 2 τ increases indicating the significance of the finite size effect in the thermalization of the non-Gibbs regime for h → h β=0 . The somewhat non-smooth nature of the growth might be related to the small number of initial conditions used for the averaging. Nevertheless, it can be conjectured that in the limit N → ∞ the excitations in the non-Gibbs regime will be persistent. It is also worthwhile to note that in the calculation of the variance, σ 2 τ , the excitations whose life time is higher than the total integration time, T = 10 7 , are not considered. On the other hand, the decay of the second moment, m 2 of the FTA distribution shown in Fig. 8 indicates that there are no such long-living excitations. Since we observed that the finite size plays a crucial role in the non-Gibbs regime for h → h β=0 , we next consider a point in the parameter space (a, h), that is far from h β=0 . This higher ratio of h h β=0 can be more straightforwardly obtained for a small norm and we fix (a = 0.2, h = 0.4) (orange symbol in Fig. 1). As the ratio h h β=0 becomes larger, the contribution of the nonlinear interaction term in Eq. (3) is higher due to the boundedness of the kinetic energy [23]. The corresponding amplitude evolutions are shown in Fig. 10 for three different values of µ. For all these µ values, the parameter set (a = 0.2, h = 0.4) is in the non-Gibbs regime and accordingly, the amplitude evolution shows the presence of at least one excitation with life time greater than the total computation time. This hints at a non-ergodic behavior for the considered long but finite time.
In summary, our analysis has illustrated that the Gibbs regime for the DNLS model remains Gibbsian for the SM of the present considerations in line with the corresponding analytical results. The situation becomes considerably more interesting beyond the thermalization limit of β = 0 for the DNLS; progressively the thermalization region of the (a, h) space expands upon increase of µ rendering non-thermalized parameter ranges for smaller µ thermalized as µ grows past a certain threshold. However, there exist additional features to consider. On the one hand, the finiteness of the lattice plays a considerable role as to whether non-ergodicity will be preserved and it is indeed found that finite lattices may lead to an apparent ergodization. Nevertheless, as the domain size grows to infinity, so does the life time of the highamplitude, nonlinear excitations in line with the thermodynamic model prediction. At the same time, the thermalization in the non-Gibbs regime heavily depends on the h h β=0 ratio. When the latter becomes larger, then the system may find itself in a non-ergodic regime even in the case of the finite lattice (long-time) computations. Admittedly, these features are the ones that are qualitatively reported here and merit additional quantification through extensive and highly-demanding (in their computation time and accuracy) computations. Nevertheless, we trust that the above results offer useful insights in this direction.

V. CONCLUSIONS & FUTURE CHALLENGES
The thermalization in the nonintegrable regime of the Salerno model (SM) has been investigated by using a combination of analytical and numerical techniques. More specifically, we have complemented the transfer integral operator (TIO) analysis by performing relevant long-time numerical simulations. In the latter, we have used a set of diagnostics such as the probability distribution of amplitudes and finite time averages of local probability density (and its moments) as well as, e.g., Lyapunov characteristic exponents. A key feature of the model is the coexistence of local nonlinearity and  nonlinear dispersion. The competition between these two effects sets two limits, a nonintegrable and an integrable one in the system, namely the discrete nonlinear Schrödinger (DNLS) and the Ablowitz-Ladik (AL) models. One can then mono-parametrically interpolate between these limits. Our study has mainly focused on the non-integrable limit, but providing a sense of the qualitative variation of the thermodynamic properties as the integrable limit is gradually approached. In the former limit the statistical mechanics of the system yields a phase diagram in the parameter space of energy and norm densities. An infinite temperature line decomposes the phase space into Gibbs and non-Gibbs regimes. The freedom afforded by the Salerno model is that we can achieve a transition from a non-Gibbs to a Gibbs regime for the same norm and energy parameters, upon variation of the strength of the nonlinear dispersive term (µ) interpolating between the two limits. Our analysis has shown that in the Gibbs regime, the SM is ergodic, and that as µ is increased towards the integrable limit, the region of the two-parameter space that thermalizes progressively expands. Nevertheless, there are some additional important features. More specifically, in the non-Gibbs regime, the ergodic properties heavily depend on the initial conditions and the ratio h/h β=0 . Additionally, the finite system size plays a crucial role and the TIO predictions are (expected to be) genuinely valid in the thermodynamic limit.
We believe that the present work provides insights into the thermalization of lattice systems and especially as we start approaching the integrable limit, including the role of parameters such as the lattice size and how "deep" in the non-Gibbs regime the initial conditions are. However, our results also raise a number of significant questions for future studies. Specifically, it is important to understand to what degree we can extend the present picture further towards the integrable limit and what happens in its im- mediate vicinity. Here, an important issue that arises is that additional conservation laws come into play [34] and how these can be incorporated in thermodynamic considerations. The role they play in modifying (or dynamically constraining) the picture is something that is especially relevant to understand, in the immediate vicinity of the integrable limit and then further away from it. In that vein, revisiting also related studies exploring the creation and disappearance, as well as mobility of discrete breathers as they interact with the phonon bath [45] (and how these mechanisms change while approaching integrability) would be an especially interesting direction. Additionally, while we have illustrated the relevance of finite size and of ratios such as h/h β=0 , obtaining a quantitative characterization of their role and of, e.g., the scaling dependence on the coherent structure (average) lifetime on them emerges as an especially relevant problem. This would greatly help appreciate the influence of quasi-ergodicity ideas such as those put forth in [12,19]. Lastly, all of the above features have been explored in one-dimensional contexts yet it would be rather natural to extend considerations to higher-dimensional one. These topics are presently under consideration and findings will be presented in future publications. To obtain the β → ∞ line we substitute the planewave like solution ψ n = √ de inθ with θ = π into Eq. (3) and obtain After integration over the phase variable φ m (see Eq. 6) the following expression is obtained: From this expression we find in the limit µ = 0 (no nonlocal nonlinearity) the whole set of equations derived for the DNLS model with only local nonlinearity [23,45].
The line β = 0 which separates the Gibbsian from the non-Gibbsian regime for the Salerno lattice can be obtained in two ways: a. Method I: Analytical In the limit β → 0, I 0 (2β A m A m+1 ) ≈ 1. Now we take βγ = x and βα = y, A = z. Then Eq. (A3) can be expressed as where E(x, y) = (ln |1+µ Am|+ln |1+µ Am+1|)+ γ 2µ (Am+Am+1)) . (A9) In order to evaluate the integral, we consider the thermodynamic limit N → ∞ of the system and evaluate using the integral equation is the kernel of the integral operator of Eq. (A10).
Here K(x, z) is symmetric and in the limit z → ∞, K(x, z)dx dz should converge. In the TIO, the partition function can be expressed in a simple form: Z ≈ (2πλ 0 ) N , where the largest eigenvalue (λ 0 ) of the kernel function is only taken into account. Therefore, approximately the norm and energy densities are