-- 1 -
Characteristics of Equilibrated Nonlinear Oscillator Systems
Uri Levy Weizmann Institute of Science, Rehovot 7610001, Israel E-mail: [email protected]
Keywords:
Discrete nonlinear Schrödinger equation (DNLSE), Nonlinear systems, Grand canonical formalism, Temperature, Chemical potential, Field correlations.
References [1]
Kevrekidis, Panayotis G.,
The discrete nonlinear Schrödinger equation: mathematical analysis, numerical computations and physical perspectives,
Vol. 232. Springer Science & Business Media, 2009. [2]
J.C. Eilbeck, P.S. Lomdahl, and A.C. Scott, Physica 16D (1985) 318-338. [3]
Eilbeck J. C. and Johansson M., 2003,
The discrete nonlinear Schrödinger equation - 20 years on , Proc. 3rd Conf.: Localization and Energy Transfer in Nonlinear Systems (June 17–21 2002 World Scientific) San Lorenzo de El Escorial, Madrid p 44 (arXiv: nlin/0211049). [4]
Davydov, A. S., and N. I. Kislukha,
Solitary excitons in one‐dimensional molecular chains, physica status solidi (b) 59, no. 2 (1973): 465-470. [5]
Davydov, A. S.,
The theory of contraction of proteins under their excitation,
Journal of Theoretical Biology 38, no. 3 (1973): 559-569. [6]
Malomed, Boris, and Michael I. Weinstein,
Soliton dynamics in the discrete nonlinear Schrödinger equation,
Physics Letters A 220, no. 1-3 (1996): 91-96. [7]
Rasmussen, K. Ø., T. Cretegny, Panayotis G. Kevrekidis, and Niels Grønbech-Jensen,
Statistical mechanics of a discrete nonlinear system,
Physical review letters 84, no. 17 (2000): 3740. [8]
Rumpf, Benno,
Transition behavior of the discrete nonlinear Schrödinger equation,
Physical Review E 77, no. 3 (2008): 036606. [9]
B. Rumpf,
Stable and metastable states and the formation and destruction of breathers in the discrete nonlinear Schrödinger equation , Physica D (Amsterdam) 238D, 2067 (2009). [10]
Polkovnikov, Anatoli, Subir Sachdev, and S. M. Girvin,
Nonequilibrium Gross-Pitaevskii dynamics of boson lattice models,
Physical Review A 66, no. 5 (2002): 053607. [11]
Mendl, Christian B., and Herbert Spohn,
Low temperature dynamics of the one-dimensional discrete nonlinear Schrödinger equation.
Journal of Statistical Mechanics: Theory and Experiment, no. 8 (2015): P08028. [12]
Meier, J., G. I. Stegeman, D. N. Christodoulides, Y. Silberberg, R. Morandotti, H. Yang, G. Salamo, M. Sorel, and J. S. Aitchison,
Experimental observation of discrete modulational instability , Physical review letters 92, no. 16 (2004): 163902. [13]
Pelinovsky, Dmitry E.,
Localization in periodic potentials: from Schrödinger operators to the Gross–Pitaevskii equation , Vol. 390. Cambridge University Press, 2011. [14]
Lederer, Falk, George I. Stegeman, Demetri N. Christodoulides, Gaetano Assanto, Moti Segev, and Yaron Silberberg,
Discrete solitons in optics,
Physics Reports 463, no. 1-3 (2008): 1-126. [15]
Smerzi, A., A. Trombettoni, P. G. Kevrekidis, and A. R. Bishop,
Dynamical superfluid-insulator transition in a chain of weakly coupled Bose-Einstein condensates.
Physical review letters 89, no. 17 (2002): 170402. [16]
Levy, Uri,
Ground States of Coupled Nonlinear Oscillator Systems , arXiv preprint arXiv:2009.14006 (2020). [17]
Silberberg, Yaron, Yoav Lahini, Yaron Bromberg, Eran Small, and Roberto Morandotti,
Universal correlations in a nonlinear periodic 1D system,
Physical review letters 102, no. 23 (2009): 233904. [18]
Levy, Uri, Ken Yang, Noam Matzliah, and Yaron Silberberg,
Universal correlations after thermalization in periodic nonlinear systems,
Journal of Physics B: Atomic, Molecular and Optical Physics 51, no. 3 (2018): 035401. [19]
Levy, Uri, and Yaron Silberberg,
Equilibrium temperatures of discrete nonlinear systems , Physical Review B
98, no. 6 (2018): 060303. [20]
Hermann Claudine,
Statistical Physics: including applications to condensed matter,
Springer Science & Business Media, 2006. [21]
Rebuzzini, Laura, Roberto Artuso, Shmuel Fishman, and Italo Guarneri,
Effects of atomic interactions on quantum accelerator modes , Physical Review A 76, no. 3 (2007): 031603. [22]
Roberts, B. 2020. Time Reversal.
The Routledge Companion to the Philosophy of Physics , Eleanor Knox and Alistair Wilson (eds). [23]
Wolfram MathWorld, Modified Bessel Function of the First Kind http://mathworld.wolfram.com/ModifiedBesselFunctionoftheFirstKind.html [24]
M. Popovic, arXiv preprint arXiv:1711.07326 (2017). [25]
Baierlein, Ralph,
The elusive chemical potential , American Journal of Physics , no. 4 (2001): 423-434. [26] Esteve, J., C. Gross, A. Weller, S. Giovanazzi, and M. K. Oberthaler,
Squeezing and entanglement in a Bose–Einstein condensate,
Nature 455, no. 7217 (2008): 1216-1219. [27]
Pitaevskii, Lev Petrovitch, Stringari, Sandro,
Bose-Einstein Condensation and Superfluidity , International series of monographs on physics 164, Oxford University Press, (2016). [28]
Huber, Sebastian D., Barbara Theiler, E. Altman, and Gianni Blatter,
Amplitude mode in the quantum phase model,
Physical review letters 100, no. 5 (2008): 050404. [29]
Polyanin, Andrei D., and Alexander V. Manzhirov.
Handbook of integral equations . CRC press, 2008. [30]
See http://mathworld.wolfram.com/ModifiedBesselFunctionoftheFirstKind.html
Abstract
During the evolution of coupled nonlinear oscillators on a lattice, with dynamics dictated by the discrete nonlinear Schrödinger equation (DNLSE systems), two quantities are conserved: system energy (Hamiltonian) and system density (number of particles). If the number of system oscillators is large enough, a significant portion of the array can be considered to be an “open system”, in intimate energy and density contact with a “bath” - the rest of the array. Thus, as indicated in previous works, the grand canonical formulation can be exploited in order to determine equilibrium statistical properties of thermalized DNLSE systems. In this work, given the values of the two conserved quantities, we have calculated the necessary values of the two Lagrange parameters (typically designated 𝛽, 𝜇 ) associated with the grand canonical partition function in two different ways. One is numerical and the other is analytic, based on a published approximate entropy expression. In addition we have accessed a purposely-derived approximate PDF expression of site-densities. Applying these mathematical tools we have generated maps of temperatures, chemical potentials, and field correlations for DNLSE systems over the entire thermalization zone of the DNLSE phase diagram, subjected to all system-nonlinearity levels. The end result is a rather complete picture, characterizing equilibrated large DNLSE systems. Introduction
The dynamics of a number of periodic systems - molecular, mechanical, optical, and lattice-trapped ultra-cold atoms, can be approximated by the DNLSE. This being the case, a book [1] , and a wealth of scientific papers studying the properties of DNLSE-governed systems were published over the years - [2] - [11] to cite just a few. The simplest two-term DNLSE reads [3] , [12] : 𝑖 ⋅ 𝑑𝑈 𝑚 𝑑𝜁 = 𝐶 ⋅ (𝑈 𝑚−1 + 𝑈 𝑚+1 ) + 𝛾 ⋅ |𝑈 𝑚 | ⋅ 𝑈 𝑚 (1) where 𝜁 is the evolution coordinate (distance or time), 𝑈 𝑚 (𝜁) is the complex field function of the oscillator at site 𝑚 , the parameter 𝐶 is the nearest-neighbor coupling constant and 𝛾 is the unharmonic parameter. In the next section we present a modified dynamics equation (Eq. (3) ) with the unharmonic parameter (𝛾) replaced by a nonlinearity parameter (𝛤) defined as 𝛤 ≡ 𝛾/|𝐶| . DNLSE systems are Hamiltonian systems [1] , [13] . The equivalent Hamiltonian from which the two-term DNLSE is derived (designated ℋ 𝑎 below – Eq. (5) ) is made up of two “energy” terms [2] , [7] , [14] – a tunneling energy term [8] (designated ℋ below – Eq. (6) ) and interaction energy term [8] (designated ℋ below – Eq. (6) ). The equivalent Hamiltonian ( energy ), is a conserved quantity of the (isolated) system. The other conserved quantity of the system (designated 𝒲 𝑎 below – Eq. (7) ) is density (number of particles). The Hamiltonian-derived DNLSE equation (Eq. (1) above or Eq. (3) below) has only two integrals of motion (conserved quantities) [1] , [14] , [15] and therefore, for a system of more than two sites, has no general analytic solution. Namely – no analytic expression for amplitude and phase of each and every oscillator at every distance (time) and for every possible set of initial excitation conditions. However, statistical properties of thermalized systems at equilibrium are analytically predictable, given the set of initial excitation conditions. Predictions of these equilibrium properties are the subject of the present study. A DNLSE system initially excited onto a specific point (𝓌 𝑎 = 𝒲 𝑎 𝑁 , 𝒽 𝑎 = ℋ 𝑎 𝑁 ) on a well-defined thermalization zone of the DNLSE phase diagram (cf. Figure 2 and Eq. (2) below) will thermalize [7] . Namely, at long evolution distances all systems initially excited onto the same point of the phase diagram (i.e. same (𝓌 𝑎 , 𝒽 𝑎 ) values) are doomed to reach the same equilibrium statistical properties, regardless of initial system-excitation details (cf. Figure 3 below). And, as stated above, these equilibrium statistical properties are predictable. Many published papers discuss properties of DNLSE systems initially excited onto the breather-forming zone (above the thermalization zone) of the DNLSE phase diagram. Other papers typically discuss a single property of systems initially excited onto the thermalization zone of the DNLSE phase diagram – transition behaviour [8] , entropy (with focus on breather dynamics) [9] , system instability [12] , or system ground states [16] . Characteristics of thermalized systems such as PDFs of site-densities (amplitudes squared, designated "𝐼" below), temperatures, chemical potentials, and field correlations for all initial excitation conditions and for all values of the nonlinearity parameter (𝛤) were not published in the scientific literature to-date. Previous studies predicted PDFs and temperatures for strong system nonlinearities based on the quantum phase approximation that justified system-entropy separation into the sum of density-entropy and a relative-phase entropy. Entropy separation led to the derivation of system temperatures as well as analytic expressions of equilibrium PDF (𝐼) and equilibrium PDF (𝜃) [17] - [19] . Here we have taken a rigorous approach of predicting system equilibrium statistics based on the thermodynamics formalism of grand canonical ensembles [7] (cf. Figure 2 below and related text). The thermodynamics approach allowed the extension of system-characteristics predictions to cover the entire thermalization zone – the zone of strong system nonlinearity as well as the zone of weak system nonlinearity. The grand canonical statistics is determined by two Lagrange parameters - (𝛽, 𝜇) . The value of each of these two parameters is uniquely determined by the value of the two DNLSE conserved quantities – density and energy. Once (𝛽(𝓌 𝑎 , 𝒽 𝑎 ), 𝜇(𝓌 𝑎 , 𝒽 𝑎 )) are determined, equilibrium statistical properties of the system studied such as entropy, temperature, chemical potential, nearest-neighbors field correlations, and PDF (𝐼) can be calculated. In the study presented here we have determined (𝛽, 𝜇) given (𝓌 𝑎 , 𝒽 𝑎 ) in one of two ways. First, by numerically inverting thermodynamics partial derivatives. Such 2D numerical inversion is challenging and could sometimes introduce uncertainty errors. Second, by finding 𝛽 𝑒𝑥𝑝 through a direct execution of an analytic expression. We have derived the analytic expression from an approximate analytic expression of system entropy published in [9] . Once the value of 𝛽 𝑒𝑥𝑝 is determined, the value of 𝜇 𝑒𝑥𝑝 is calculated through solving a one parameter implicit integral equation. In addition to finding equilibrium PDF (𝐼) by rigorous numerical calculations, we have gone in this study through deriving an approximate analytic PDF (𝐼) expression. Our derivation here is based on approximating the kernel associated with the partition function in the grand canonical formalism (cf. section and Appendix 1 ). The easily executed approximate analytic PDF (𝐼) expression “works” well on a large area of the thermalization zone and brings additional system-characteristics-related insights. For example - the identification of system’s temperature with the variance of the system’s PDF (𝐼) . Following is a summarizing list of the mathematical tools we have employed in the present study to quantitatively describe key equilibrium properties of thermalized DNLSE systems under all initial excitation conditions: • Simulations of system evolution: given the amplitude and phase of each and every oscillator in the array, execute the DNLSE (Eq. (3) below) to stationary-statistics distances. Next, calculate equilibrium PDFs of densities and of relative phases. Clearly, PDFs calculated this way faithfully describe the PDFs of actual thermalized systems. • The rigorous grand canonical statistics suggested in [7] . This mathematical tool requires numerical inversion of two thermodynamic equations to compute exact values of the two associated Lagrange parameters (𝛽, 𝜇) . Small numerical errors in the values determined this way often creep-in. • The grand canonical statistics with approximate values of the two associated Lagrange parameters (𝛽 𝑒𝑥𝑝 , 𝜇 𝑒𝑥𝑝 ) . The parameter 𝛽 𝑒𝑥𝑝 is determined by an analytic expression derived from an approximate analytic expression for system entropy published in [9] (hence the subscript “ exp ”, standing for expression ). The parameter 𝜇 𝑒𝑥𝑝 is determined next by numerically solving a single parameter implicit equation. • Previously published analytic expressions for PDF(densities) and PDF(relative phase-angles). These analytic expressions hold for high nonlinearity DNLSE systems (strong nonlinearity of the oscillators and/or high average of oscillator amplitudes) [18] , [19] . • Approximate analytic expression for PDF(densities). During the present study we have derived an approximate analytic expression to PDF(densities) of equilibrated DNLSE systems under a second-order approximation to the kernel of the grand canonical partition function (cf. details in
Appendix 1 ). PDFs calculated through the second-order-derived analytic expression well match the evolution-simulated PDFs for a surprisingly large class of equilibrated DNLSE systems. In our dynamics-governing DNLSE (Eq. (3) ), and unlike many DNLSE versions of other studies, we have left the oscillator fields (𝑈 𝑚 ′𝑠) dimensional and left the nonlinearity parameter (𝛤) dimensional as well. For example, if the DNLSE system is an array of proximity optical waveguides then the dimension of the complex site functions squared is typically power/volume. In this case the dimension of “energies” is also power/volume. The dimension of the nonlinearity parameter (𝛤) , and the dimension of the Lagrange parameter 𝛽 is inverse (power/volume). The dimension of Eq. (3) below is sqrt (power/volume). This amplitude-nonlinearity separation allowed us here the study of the effect of each the three variables - system density (𝓌 𝑎 ) , system energy (𝒽 𝑎 ) , and the value of the nonlinearity parameter (𝛤) independently, keeping the other two variables constant. In the next three sections we elaborate on the process of system thermalization and on the related mathematics. In the sections that follow we mathematically analyze and graphically demonstrate the characteristics of equilibrated DNLSE systems: temperature, chemical potential, field correlations and PDF(densities). Let us start then in discussing the thermalization of DNLSE systems. Thermalization of DNLSE systems
In this stage-setting section we describe the DNLSE phase diagram, discuss in more details the grand canonical formulation as applied to DNLSE systems, and present a specific thermalization example. The thermalization example consists of showing evolution-simulated and analytically predicted equilibrium PDFs for two systems. Both systems are placed onto the same location of the phase diagram (same (𝓌 𝑎 , 𝒽 𝑎 ) , see below) but each system is launched by a unique set of initial excitation conditions. The two conserved quantities, most frequently appearing in this study, are the intensive site-averaged quantities . Namely - system density (𝓌 𝑎 ) and system energy (𝒽 𝑎 ) . The two extensive conserved quantities are thus total system density - 𝒲 𝑎 = 𝑁 ⋅ 𝓌 𝑎 , and total system energy - ℋ 𝑎 = 𝒽 𝑎 ⋅ 𝑁 , were 𝑁 is the number of oscillators in the array. A DNLSE phase diagram (not to be confused with “phase space”) can now be constructed on the (𝓌 𝑎 , 𝒽 𝑎 ) plane (cf. Figure 1 ). Figure 1:
Phase diagram for DNLSE systems with a positive nonlinearity parameter (𝛤 = +1.0) . The magenta area is the thermalization zone, limited at the bottom by the zero temperature line (or the “ground state line” - 𝒽 𝑎0 ), and limited at the top by the infinite temperature line (𝒽 𝑎∞ ) . Crossing the thermalization zone is an intermediate 𝐿 𝑖 line. The colored circular disk on the 𝐿 𝑖 line represents a system with all oscillators initially excited at a uniform amplitude ( = √6 ) with random phases in the full range of [0,2 ⋅ 𝜋) . The dark magenta area above the thermalization zone is a negative temperature (or “breather-forming”) zone. The study presented here is devoted solely to the equilibrium properties of DNLSE systems initially excited onto the thermalization zone. The DNLSE phase diagram is divided into three zones – an inaccessible zone, a thermalization zone, and a negative temperature (or a “breather-forming”) zone. The thermalization zone is limited from below by a zero temperature line (𝒽 𝑎0 ) and is limited from above by an infinite temperature line (𝒽 𝑎∞ ) [7] . Crossing the thermalization zone is an intermediate 𝐿 𝑖 line [19] . The mathematical expressions for these three lines are - 𝒽 𝑎0 (𝓌 𝑎 ) = −2 ⋅ 𝓌 𝑎 + 12 ⋅ 𝛤 ⋅ 𝓌 𝑎2 𝐿 𝑖 (𝓌 𝑎 ) = 12 ⋅ 𝛤 ⋅ 𝓌 𝑎2 𝒽 𝑎∞ (𝓌 𝑎 ) = 𝛤 ⋅ 𝓌 𝑎2 (2) The 𝐿 𝑖 line is characterized by two unique properties. The first property is related to system initial excitation conditions. To be on the 𝐿 𝑖 line, systems can be initially excited with uniform amplitudes and random phases in the full range of [0,2 ⋅ 𝜋) . As the window from which random phase-angles are drawn is getting narrower and narrower, same-amplitude systems are placed closer and closer to the ground state line (𝒽 𝑎0 ) . To place a system above the 𝐿 𝑖 line, initial excitation amplitudes must be spread apart, more and more so getting closer and closer to the infinite temperature line (𝒽 𝑎∞ ) . The second property is related to two limiting cases. The 𝐿 𝑖 line becomes the lower (upper) limit of the thermalization zone as the nonlinearity parameter (𝛤) goes to infinity (zero). Following [19] , we refer to the thermalization zone below the intermediate 𝐿 𝑖 line as a cold zone and refer to the thermalization zone above the intermediate 𝐿 𝑖 line as a hot zone . In the coming sections we show quantitatively that indeed system temperatures, chemical potentials, field correlations, and PDFs of densities and relative phase-angles, change moderately for systems on the cold zone and change drastically for systems on the hot zone. The phase diagram of Figure 1 is plotted for a positive nonlinearity parameter (𝛤 > 0) . In [16] the Hamiltonian corresponding to a positive nonlinearity parameter is referred-to as a positive Hamiltonian . Here, without loss of generality, we consider only positive Hamiltonian systems, given mirror equivalence. Namely – statistical properties of a system on a point (𝓌 𝑎 , 𝒽 𝑎 ; 𝛤 > 0) are equal to the properties of a system on its mirror image position - (𝓌 𝑎 , −𝒽 𝑎 ; 𝛤 < 0) . (In [8] the nonlinearity parameter is altogether normalized away, i.e. - its value is fixed at a positive unity). Schematic of a grand canonical setting as applied to a large array of coupled oscillators is depicted by Figure 2 . A large portion of the array – the “open system” is in density contact and in energy contact with the rest of the array – the “bath”. At equilibrium, the ensemble averages (〈𝓌 a 〉, 〈𝒽 a 〉) of the open system take on the values of the two conserved quantities (𝓌 𝑎 , 𝒽 𝑎 ) of the bath, determined as soon as the entire array is initially excited. Specifically, the values of the Lagrange parameters - 𝛽(𝓌 𝑎 , 𝒽 𝑎 ), 𝜇(𝓌 𝑎 , 𝒽 𝑎 ) of the grand canonical partition function are fixed such that 〈𝓌 a (𝛽, 𝜇)〉 = 𝓌 a and 〈𝒽 a (𝛽, 𝜇)〉 = 𝒽 a [20] . Once 𝛽 and 𝜇 are determined, equilibrium statistical properties of the entire (isolated) array can be analytically predicted. Figure 2:
Grand canonical statistics applied to a large array of coupled oscillators. A significant portion of the array is (artificially) assumed to be an “open system” in contact with the rest of the oscillators in the array. At equilibrium, through diffusion-exchange of energy and density between the “bath” and the open system, the grand canonical ensemble averages (〈𝓌 𝑎 〉, 〈𝒽 𝑎 〉) settle on values equal to the known conserved density (𝓌 𝑎 ) and conserved energy (𝒽 𝑎 ) values of the entire excited array [20] . This equality dictates the values of the Lagrange parameters - 𝛽(𝓌 𝑎 , 𝒽 𝑎 ), 𝜇(𝓌 𝑎 , 𝒽 𝑎 ) . Once 𝛽 and 𝜇 are accurately determined, equilibrium statistical properties of the entire (isolated) array can be faithfully derived. Systems on the thermalization zone evolve to equilibrium [7] . Entropy rises to its maximum and stays there [16] , and other system statistical properties such as tunneling energy, interaction energy, PDFs of densities and of relative phase-angles, and field correlations, reach a stationary value or shape. These stationary statistical properties are predictable, based on the Gibbsian formalism applied to a grand canonical ensemble [7] . Figure 3:
Thermalization of DNLSE systems. A. PDF (𝐼) curves of initial excitation densities of two DNLSE systems. Considering the DNLSE phase diagram, both systems are initially excited onto the same thermalization-zone point: (𝓌 a , 𝒽 a ) =(16.0,147.2) . B. Equilibrium PDF (𝐼) curves of the two same-point-initially-excited systems. The red and blue points are respective PDFs (𝐼) calculated by evolution-to-equilibrium simulations (averaged over six realizations). The shown two practically identical red and blue simulated curves prove system thermalization, independent of initial excitation details. The continuous yellow curve and the dashed green curve are theoretically-predicted equilibrium PDF(I) curves for the given (𝓌 a , 𝒽 a ) point of the thermalization zone. The yellow curve by previously published analytic expressions assuming strong system nonlinearity [19] . The dashed green curve by the rigorous grand-canonical statistics as discussed in section below. Two PDF (𝐼) predictions of thermalized DNLSE systems are shown in panel B of
Figure 3 . One curve by previously published analytic expressions assuming strong system nonlinearity [19] . Another curve by the rigorous grand-canonical statistics as discussed in section below. Dynamics and conserved quantities
The evolution dynamics of a array of 𝑁 coupled unharmonic oscillators is given by Eq. (1) above. Throughout this work, without loss of generality, we will use the notation of optics (evolution coordinate 𝜁 as distance, or 𝑧 as a normalized distance, with coupled optical waveguides in mind). In this study we assume 𝑁 to be very large and, subject the equation to periodic boundary conditions (𝑈 𝑚+𝑁 = 𝑈 𝑚 ) . The equation consists of two terms – a linear hopping term , and a cubic on-site nonlinear term . Several options for normalizing Eq. (1) are available [1] , [3] and are often applied [10] , [17] . Here, as in [16] , we shall eliminate the coupling constant from Eq. (1) except for its sign, following a division by |𝐶| – 𝑖 ⋅ 𝑑𝑈 𝑚 𝑑𝑧 = 𝑠𝑖𝑔𝑛𝐶 ⋅ (𝑈 𝑚−1 + 𝑈 𝑚+1 ) + 𝛤 ⋅ |𝑈 𝑚 | ⋅ 𝑈 𝑚 𝑧 ≡ |𝐶| ⋅ 𝜁 ; 𝑠𝑖𝑔𝑛𝐶 ≡ 𝑠𝑖𝑔𝑛(𝐶) ; 𝛤 ≡ 𝛾|𝐶| (3) In Eq. (3) , the amplitudes 𝑈 𝑚′ 𝑠 are dimensional. The evolution coordinate (𝑧) – “distance” (or “time”) is dimensionless. The dimension of the nonlinearity parameter 𝛤 is [𝑈 𝑚 ] −2 . System nonlinearity is expressed by the dimensionless product of the nonlinearity parameter and the conserved system density (𝑖. 𝑒. 𝑎𝑠 (𝛤 ⋅ 𝓌 𝑎 ) ). If 𝑠𝑖𝑔𝑛𝐶 = 𝑠𝑖𝑔𝑛(𝛤) / −𝑠𝑖𝑔𝑛(𝛤) then Eq. (3) is a “focusing” / “defocusing” version of the DNLSE systems [1] , [21] . It is convenient at this point, and indeed done in almost every DNLSE study, to perform a Madelung transformation on the complex canonical coordinates (𝑈 𝑚 , 𝑖 ⋅ 𝑈 𝑚∗ ) into the set of density-angle real canonical polar coordinates: (𝑈 𝑚 , 𝑖 ⋅ 𝑈 𝑚∗ ) → (𝐼 𝑚 , 𝜙 𝑚 ) . The complex field functions (𝑈 𝑚 ′𝑠(𝑧)) of Eq. (3) take then the form: 𝑈 𝑚 = 𝑢 𝑚 ⋅ 𝑒 𝑖⋅𝜙 𝑚 ; 𝑢 𝑚 ≡ √𝐼 𝑚 ; 𝜃 𝑚 ≡ 𝜙 𝑚 − 𝜙 𝑚+1 (4) Equation (3) can be derived from an equivalent Hamiltonian ( ℋ 𝑎 ) which is a conserved quantity, associated with the system’s time translation invariance [13] , [22] -
10 - ℋ 𝑎 = ∑ {𝑠𝑖𝑔𝑛𝐶 ⋅ (𝑈 𝑚∗ ⋅ 𝑈 𝑚+1 + 𝑈 𝑚 ⋅ 𝑈 𝑚+1∗ ) + 𝛤2 ⋅ |𝑈 𝑚 | } 𝑁𝑚=1 (5)
The variables (𝑈 𝑚 , 𝑖 ⋅ 𝑈 𝑚∗ ) are canonical variables. Adopting the assignment 𝑞 𝑚 =𝑈 𝑚 ; 𝑝 𝑚 = 𝑖 ⋅ 𝑈 𝑚∗ , Eq. (3) is derived from the Hamiltonian (5) as 𝑑𝑈 𝑚 𝑑𝑧 = 𝜕ℋ 𝑎 𝜕(𝑖⋅𝑈 𝑚∗ ) . In the polar variables (𝐼 𝑚 , 𝜙 𝑚 ) of Eq. (4) , the DNLSE Hamiltonian (Eq. (5) ) takes on the form: ℋ (𝑧) = 𝑠𝑖𝑔𝑛𝐶 ⋅ ∑ 2 ⋅ 𝑢 𝑚 𝑢 𝑚+1 ⋅ cos 𝜃 𝑚𝑁𝑚=1 ; ℋ (𝑧) = ∑ 𝛤2 ⋅ 𝑢 𝑚4𝑁𝑚=1 ℋ 𝑎 ≡ ℋ (𝑧) + ℋ (𝑧) 𝒽 (𝑧) ≡ ℋ (𝑧)𝑁 ; 𝒽 (𝑧) ≡ ℋ (𝑧)𝑁 ; 𝒽 𝑎 ≡ ℋ 𝑎 𝑁 (6) The Hamiltonian energies - (ℋ (𝑧), ℋ (𝑧)) are the nearest-neighbor tunneling energy term and the on-site interaction energy term respectively. Obviously, both ℋ (𝑧) and ℋ (𝑧) vary with propagation distance, but their sum does not. During DNLSE evolution then, an energy diffusion process transfers energy from ℋ (𝑧) to ℋ (𝑧) or the other way around. Another conserved quantity of DNLSE systems, thanks to the system’s invariance with respect to global phase rotations [1] , [13] , is “density” (𝒲 𝑎 ) (or norm , or number of particles ) given by - 𝒲 𝑎 = ∑ 𝐼 𝑚 (𝑧) 𝑁𝑚=1 ; 𝓌 𝑎 ≡ 𝒲 𝑎 𝑁 (7) As stated above, system density and system energy (𝓌 𝑎 , 𝒽 𝑎 ) form a plane over which the DNLSE phase diagram is graphically represented (cf. Figure 1 ). Grand canonical partition function and PDF of site-densities
The grand canonical partition function (𝒵(𝛽, 𝜇)) is given as [7] : 𝒵(𝛽, 𝜇) = ∫ ∫ ∏ 𝑑𝜙 𝑚 ⋅ 𝑑𝐼 𝑚 ⋅ 𝑒 −𝛽⋅(ℋ+𝜇⋅𝒲)𝑁𝑚=12𝜋0∞0 (8)
11 -
Following integration over the phase variable (𝜙 𝑚 ) the expression for the partition function is reduced to - 𝒵(𝛽, 𝜇) = (2𝜋) 𝑁 ⋅ ∫ ∏ 𝒦(𝐼 𝑚 , 𝐼 𝑚+1 ) 𝑁𝑚=1∞0 ⋅ 𝑑𝐼 𝑚 𝒦(𝐼 𝑚 , 𝐼 𝑚+1 ) ≡ 𝑒 [−𝛽⋅(𝜇2⋅𝐼 𝑚 +𝛤4⋅𝐼 𝑚2 )−𝛽⋅(𝜇2⋅𝐼 𝑚+1 +𝛤4⋅𝐼 𝑚+12 )] ⋅ ℐ (2 ⋅ 𝛽 ⋅ √𝐼 𝑚 ⋅ 𝐼 𝑚+1 ) (9) where ℐ (∙) is the zero order of modified Bessel function of the first kind, satisfying for integer 𝑛 [23] : ℐ 𝑛 (𝑧) = ⋅ ∫ 𝑒 𝑧⋅cos 𝜃 ⋅ cos(𝑛 ⋅ 𝜃) ⋅ 𝑑𝜃 𝜋0 . The nonnegative symmetric kernel (𝒦(𝑥, 𝑦)) of Eq. (9) satisfies ∫ 𝒦(𝑥, 𝑦) ⋅ 𝑣 𝑖 (𝑥) ⋅ 𝑑𝑥 = 𝜆 𝑖 ⋅ 𝑣 𝑖 (𝑦) . If the array of oscillators is large enough then the partition function, to a very good approximation, is given by the largest eigenvalue of the kernel [7] : 𝒵(𝛽, 𝜇) ≅ [2 ⋅ 𝜋 ⋅ 𝜆 (𝛽, 𝜇)] 𝑁 (10) and the PDF (𝐼) of site-densities (𝒫 𝐼 (𝐼)) is given by the square of the eigenfunction (𝑣 (𝐼)) associated with the largest eigenvalue [7] : 𝒫 𝐼 (𝐼) = 𝑣 (𝐼) (11) For numerical calculations we treat the kernel (
𝒦(𝑥, 𝑦) ) as a 2D matrix on a Δ -spaced grid and solve for 𝜆 𝑖 and 𝑣 𝑖 : (𝒦(𝑥, 𝑦) ⋅ Δ) ⋅ 𝑣 𝑖 = 𝜆 𝑖 ⋅ 𝑣 𝑖 (12) where now 𝑣 𝑖 are the normalized eigenvectors of the matrix (𝒦(𝑥, 𝑦) ⋅ Δ) . The kernel 𝒦(𝑥, 𝑦) is written explicitly as:
𝒦(𝑥, 𝑦) = 𝑒 [−𝛽⋅(𝜇2⋅𝑥+𝛤4⋅𝑥 )−𝛽⋅(𝜇2⋅𝑦+𝛤4⋅𝑦 )] ⋅ ℐ (2 ⋅ 𝛽 ⋅ √𝑥 ⋅ 𝑦) (13) Two sets of equations relate the two conserved quantities (𝓌 𝑎 , 𝒽 𝑎 ) to the two Lagrange parameters of the grand canonical partition function (𝛽, 𝜇) . The first set is the average-calculating integrals – 𝓌 𝑎 = 〈𝓌 𝑎 (𝛽, 𝜇)〉 = ∫ 𝐼 ⋅ 𝑣 (𝛽, 𝜇; 𝐼) ⋅ 𝑑𝐼 ∞0 𝒽 𝑎 = 〈𝒽 𝑎 (𝛽, 𝜇)〉 = 〈ℋ 𝑎 (𝛽, 𝜇)〉𝑁 (14)
12 -
The second set is the set of partial derivatives of the grand canonical partition function [20] : 𝓌 𝑎 = 〈𝓌 𝑎 (𝛽, 𝜇)〉 = − 1𝛽 ⋅ 𝜕{𝑙𝑛[𝜆 (𝛽, 𝜇)]}𝜕𝜇 𝒽 𝑎 = 〈𝒽 𝑎 (𝛽, 𝜇)〉 = − 𝜕{𝑙𝑛[𝜆 (𝛽, 𝜇)]}𝜕𝛽 − 𝜇 ⋅ 𝓌 𝑎 (𝛽, 𝜇) (15) In equations (14) and (15) , 𝑣 (𝛽, 𝜇; 𝐼) , 〈ℋ 𝑎 (𝛽, 𝜇)〉 , and 𝜆 (𝛽, 𝜇) are each a 2D map calculated through Eqs. (12) and (11) . Inverting (numerically) one of these two important equation sets yields the exact value of 𝛽(𝓌 𝑎 , 𝒽 𝑎 ), 𝜇(𝓌 𝑎 , 𝒽 𝑎 ) . Once (𝛽, 𝜇) are uniquely determined, exact equilibrium statistical properties of a system excited onto a point (𝓌 𝑎 , 𝒽 𝑎 ) of the thermalization zone can be determined. In the following sections we present maps of temperatures, chemical potentials, and field correlations, for equilibrated DNLSE systems. Map calculations cover the entire thermalization zone and hold for all system nonlinearity levels. The maps reveal, for example, the universality property of field correlations, and illuminate the inverse relation between system temperatures and field correlations. Further, we present here-derived analytic expressions (cf. Appendix 1 ) that well approximate the grand canonical partition function and the PDF of site-densities on a large portion of the thermalization zone. We start with the temperature of a DNLSE system. Temperature of a DNLSE system
This temperature section consists of three parts. In the first part we rigorously derive an expression for the temperature of a DNLSE system. In the second part we present an analytic expression for the temperature following an entropy expression published in [9] . The third part is devoted to graphical illustrations.
Exact expression . Temperature of an equilibrated DNLSE system is calculated through the energy derivative of Gibbs entropy. The Gibbs entropy equation relates the entropy of a system to the probability distribution of the microstates [24] . Applied to the grand canonical ensemble considered here, continuous Gibbs entropy is given as [8] , [17] : 𝑆(𝛽, 𝜇) ≡ −𝑘 ⋅ ∫ ∫ 𝒫(𝛽, 𝜇) ⋅ 𝑙𝑛[ℎ ⋅ 𝒫(𝛽, 𝜇)] 𝑑 𝑁 𝐼 ⋅ 𝑑 𝑁 𝜃 (16) where 𝑘 is the Boltzmann constant (later on to be set at dimensionless unity). The integration sign 𝑑 𝑁 𝐼 stands for ∏ 𝑑𝐼 𝑚𝑁1 and similarly for 𝑑 𝑁 𝜃 . The value of the scale parameter ℎ , with units of 𝑒𝑛𝑒𝑟𝑔𝑦 𝑁 , is set at one. The probability density of any given state 𝒫(𝐼 , 𝐼 , … 𝐼 𝑁 , 𝜃 , 𝜃 , … 𝜃 𝑁 ; 𝛽, 𝜇) is
13 -
𝒫(𝐼 , 𝐼 , … 𝐼 𝑁 , 𝜃 , 𝜃 , … 𝜃 𝑁 ; 𝛽, 𝜇) = 1𝒵(𝛽, 𝜇) ⋅ 𝑒 −𝛽⋅(ℋ+𝜇⋅𝒲) (17) Note the use of the pair (𝐼 𝑚 , 𝜃 𝑚 ) (cf. Eq. (4) ) in the entropy expression (16) with the probability density (17) [8] , [17] and not the use of the canonical pair (𝐼 𝑚 , 𝜙 𝑚 ) . Inserting (17) into (16) and working out the integral using (10) :
1𝑘 ⋅ 𝑆(𝛽, 𝜇) = 𝑁 ⋅ 𝑙𝑛[2 ⋅ 𝜋 ⋅ 𝜆 (𝛽, 𝜇)] + 𝛽 ⋅ 𝑁 ⋅ 〈𝒽 𝑎 + 𝜇 ⋅ 𝓌 𝑎 〉 (18) Defining site-averaged system entropy 𝑠 = 𝑆/𝑁 , setting 𝑘 = 1 , and using 〈𝓌 a (𝛽, 𝜇)〉 = 𝓌 a ; 〈𝒽 a (𝛽, 𝜇)〉 = 𝒽 a , we arrive at: 𝑠(𝛽, 𝜇) = 𝑙𝑛[2 ⋅ 𝜋 ⋅ 𝜆 (𝛽, 𝜇)] + 𝛽 ⋅ (𝒽 𝑎 + 𝜇 ⋅ 𝓌 𝑎 ) (19) Since the Lagrange multipliers (𝛽, 𝜇) are each a function of (𝓌 𝑎 , 𝒽 𝑎 ) , Eq. (19) for the site-averaged system entropy (𝑠) rigorously reads: 𝑠(𝓌 𝑎 , 𝒽 𝑎 ) = 𝑙𝑛[2 ⋅ 𝜋 ⋅ 𝜆 (𝓌 𝑎 , 𝒽 𝑎 )] + 𝛽(𝓌 𝑎 , 𝒽 𝑎 ) ⋅ [𝒽 𝑎 + 𝜇(𝓌 𝑎 , 𝒽 𝑎 ) ⋅ 𝓌 𝑎 ] (20) In the present study we adhere to the definition of DNLSE system temperature as put forward in [19] : 𝑇 𝐷𝑁𝐿𝑆𝐸 (𝓌 𝑎 , 𝒽 𝑎 ) ≡ (𝛤 𝜕𝑠(𝓌 𝑎 , 𝒽 𝑎 )𝜕𝒽 𝑎 ) 𝓌 𝑎 −1 (21) Taking the derivative of 𝑠(𝓌 𝑎 , 𝒽 𝑎 ) given by Eq. (20) (multiply by 𝛤 and invert later), we arrive at the following three term expression (keeping 𝓌 𝑎 fixed): 𝜕𝑠(𝒽 𝑎 )𝜕𝒽 𝑎 = 𝜕𝑙𝑛[𝜆 (𝒽 𝑎 )]𝜕𝒽 𝑎 + 𝓌 𝑎 ⋅ 𝜕[𝛽(𝒽 𝑎 ) ⋅ 𝜇(𝒽 𝑎 )]𝜕𝒽 𝑎 + 𝛽(𝒽 𝑎 ) (22) Somewhat surprisingly, as we have shown numerically, for all (𝓌 𝑎 , 𝒽 𝑎 ) points of the thermalization zone of the DNLSE phase diagram, the first two terms of Eq. (22) cancel out: 𝜕𝑙𝑛[𝜆 (𝒽 𝑎 )]𝜕𝒽 𝑎 + 𝓌 𝑎 ⋅ 𝜕[𝛽(𝒽 𝑎 ) ⋅ 𝜇(𝒽 𝑎 )]𝜕𝒽 𝑎 = 0 (23) So that after multiplying by 𝛤 and inverting we are left with
14 - 𝑇 𝐷𝑁𝐿𝑆𝐸 (𝓌 𝑎 , 𝒽 𝑎 ) ≡ (𝛤 𝜕𝑠(𝓌 𝑎 , 𝒽 𝑎 )𝜕𝒽 𝑎 ) 𝓌 𝑎 −1 = 1𝛽(𝓌 𝑎 , 𝒽 𝑎 ) ⋅ 𝛤 (24) Down below we use just “ 𝑇 ” for the exact 𝑇 𝐷𝑁𝐿𝑆𝐸 given by Eq. (24) . The dimension of the Lagrange parameter 𝛽 is the inverse of the equivalent-Hamiltonian dimension, i.e. 𝑒𝑛𝑒𝑟𝑔𝑦 −1 . Similarly, the dimension of the nonlinearity parameter 𝛤 is also 𝑒𝑛𝑒𝑟𝑔𝑦 −1 . The dimension of the temperature of a DNLSE system is thus 𝑒𝑛𝑒𝑟𝑔𝑦 . Note that 𝑒𝑛𝑒𝑟𝑔𝑦 is also the dimension of the variance of PDF (𝐼) . Analytic expression . The author of [9] derived an expression for system entropy that is an approximation in the center of the thermalization region, and is exact at the margins. Modifying the derived expression (to explicitly include the nonlinearity parameter (𝛤) ) and taking the energy derivative (at constant density - Eq. (21) ) we arrive at the analytic expression – 𝑇 𝑒𝑥𝑝 (𝓌 𝑎 , 𝒽 𝑎 ) = 1𝛤 ⋅ 4 ⋅ (𝓌 𝑎 + 𝒽 𝑎,𝑚𝑎𝑥 (𝓌 𝑎 )4 ) − Δℎ 𝑎2 (𝓌 𝑎 , 𝒽 𝑎 )2 ⋅ Δℎ 𝑎 (𝓌 𝑎 , 𝒽 𝑎 ) 𝒽 𝑎,𝑚𝑎𝑥 (𝓌 𝑎 ) ≡ 𝛤 ⋅ 𝓌 𝑎2 ; Δℎ 𝑎 (𝓌 𝑎 , 𝒽 𝑎 ) ≡ 𝒽 𝑎,𝑚𝑎𝑥 (𝓌 𝑎 ) − 𝒽 𝑎 (25) From Eqs. (24) and (25) the Lagrange parameter 𝛽 𝑒𝑥𝑝 is determined and then the other Lagrange parameter - 𝜇 𝑒𝑥𝑝 , is calculated through solving a one parameter implicit equation given by the integral (14) : 〈𝓌 𝑎 (𝛽 𝑒𝑥𝑝 , 𝜇 𝑒𝑥𝑝 )〉−𝓌 𝑎 = ∫ 𝐼 ⋅ 𝑣 (𝛽 𝑒𝑥𝑝 , 𝜇 𝑒𝑥𝑝 ; 𝐼) ⋅ 𝑑𝐼 ∞0 −𝓌 𝑎 = 0 (26) Once (𝛽 𝑒𝑥𝑝 , 𝜇 𝑒𝑥𝑝 ) are determined, system PDF (𝐼) can be numerically computed through Eq. (12) . On the 𝐿 𝑖 line (Eq. (2) ), Eq. (25) is reduced to a particularly simple expression - 𝑇 𝑒𝑥𝑝,𝐿 𝑖 (𝓌 𝑎 ) = 4 + 2 ⋅ 𝓌 𝑎 ⋅ 𝛤𝛤 (27) Thus, if 𝑎 ⋅ 𝛤 ≫ 4 then the temperature along the 𝐿 𝑖 line is the straight line 𝑎 /𝛤 , in agreement with [19] . If, on the other hand, 𝑎 ⋅ 𝛤 ≪ 4 then the temperature along the 𝐿 𝑖 line is independent of system density (𝓌 𝑎 ) and goes to infinity with 𝛤 → 0 as 𝛤 −2 . Consulting Eq. (25) , it is clear that in general (not only on the 𝐿 𝑖 line) as 𝛤 → 0 temperatures of systems on all points of the thermalization zone (except for systems in the ground state) go to infinity. PDFs (𝐼) on the other hand are still of a finite width (see for example panel A of
Figure 17 below).
15 -
Two mathematical methods are at our disposal then to calculate the temperature of a DNLSE system initially excited onto a (𝓌 𝑎 , 𝒽 𝑎 ) point on the DNLSE phase diagram. One method – Eq. (24) -involves numerical calculation of a 𝜆 (𝛽, 𝜇) map (Eq. (12) ) and finding (𝛽, 𝜇) to satisfy the partial derivatives of Eq. (15) . Mathematically this numerical procedure should yield the exact temperature (Eq. (24) ). Practically, a certain uncertainty in the determined value creeps-in as the vicinity of the minimum of the maps involved is rather shallow. The other method is a direct execution of the analytic expression (25) . In the graphical illustrations below we use both ways as the case may be. Graphical illustrations . Temperature characteristics of equilibrated DNLSE systems such as temperature dependence on system energy (exponential-like), or dependence of temperature on the value of the nonlinearity coefficient (inverse relations) are illustrated by
Figure 4 through
Figure 10 . In
Figure 4 we show that temperatures calculated analytically (Eq. (25) ) reasonably approximate the more rigorously determined system temperatures (inversion of expressions (15) to determine 𝛽(𝓌 𝑎 , 𝒽 𝑎 ) ). Figure 4:
Temperatures of equilibrated DNLDE systems. Three of the four panels show temperatures on a log scale vs. system energy at a fixed system density as marked. Two curves are shown on each of the three panels. The discrete blue dots on the blue curve, showing exact temperatures, were calculated through inversion of the partial derivatives of Eq. (15) . The magenta curves were calculated by the analytic expression (25) . The blue-green-red vertical lines designate the crossings of the 𝒽 𝑎0 ; 𝐿 𝑖 ; 𝒽 𝑎∞ lines respectively (cf. Eq. (2) ). The lower-right panel shows three curves of 𝒫 𝐼 (𝐼) corresponding to the two marked points on the lower left panel. The red-points curve is the result of simulations, evolving the system to equilibrium (averaged over sixteen realizations). The green curve, pretty much following the red-points curve, corresponds to the green point on the left panel and was calculated through the exact procedure (Eq. (15) and Eq. (12) ). The cyan curve corresponds to the cyan dot on the left panel and was calculated by the analytic expression (Eq. (25) ). Generally, the figure indicates that whereas the temperature values predicted analytically are approximate, their general dependence on system energy and system density follows the exact-temperature dependence.
16 -
Figure 5:
Temperature vs. system energy (Eq. (25) ) at a fixed system density (𝓌 𝑎 = 2.0) , for two values of the nonlinearity parameter (𝛤 = 0.3 ; 1.0) . Energy range for both 𝛤 values is 𝒽 𝑎0 + 0.01 ⋅ Δℎ 𝑎 to 𝒽 𝑎0 + 0.95 ⋅ Δℎ 𝑎 . Vertical lines as in Figure 4 . Note the very high system temperatures at the low 𝛤 value. The curves of
Figure 5 show the rise of temperatures with system energy at a fixed system density (travelling vertically on the phase diagram). Initially temperatures rise linearly with energy, going to kind of an exponential rise to infinity for systems initially excited with higher and higher energies. Note the temperature “restriction” effect of increased system nonlinearity (𝛤 ⋅ 𝓌 𝑎 ) . Figure 6:
Temperatures and PDFs (𝐼) for a system on the 𝐿 𝑖 line having a fixed, relatively low, system density (𝓌 𝑎 = 2.0) . Temperatures and PDFs(I) are shown for three different values of the nonlinearity parameter (𝛤 = 3.0 ; 1.0 ; 0.3) . Top left panel shows analytically calculated system temperature: 𝑇 𝑒𝑥𝑝 = 1.78 ; 8.0 ; 57.8 respectively. The other three panels show simulated 𝑃 𝐼 (𝐼) curves by magenta dots corresponding to the three values of the nonlinearity parameter as indicated by the black arrows. The cyan curves were numerically calculated (Eq. (12) given
17 - (𝛽 𝑒𝑥𝑝 , 𝜇 𝑒𝑥𝑝 ) (Eqs. (24) - (26) ). The dark blue curve on the lower left panel is an exact PDF (𝐼) (Eqs. (15) , and (12) ). The figure shows the very strong dependence of system temperature on the value of the nonlinearity parameter 𝛤 , for a system on the 𝐿 𝑖 line, characterized by a relatively low system density of 𝓌 𝑎 = 2.0 . The next two figures -
Figure 6 and
Figure 7 show equilibrium temperatures for systems initially excited (at 𝑧 = 0 ) onto the intermediate (𝐿 𝑖 ) line: uniform amplitudes, 𝑢 𝑚 = √2 for all 𝑚 , and random phases in the full range of [0,2 ⋅ 𝜋) . The figures show the sharp rise of system temperature as the value of the nonlinearity parameter (𝛤) is reduced. Note that as 𝛤 → 0 the infinite temperature line (ℎ 𝑎∞ ) gets closer and closer to the intermediate 𝐿 𝑖 line (where the system in question resides), until the two lines merge at 𝛤 = 0 and system temperature is infinitely high (
Figure 7 ). It is fair then to reach the very intuitive conclusion that the energy-distance of a system from the ℎ 𝑎∞ line is a strong indication to its temperature: shorter distance – higher temperature. Figure 7:
Equilibrium distribution of densities (𝒫 𝐼 (𝐼)) for a system on the 𝐿 𝑖 line (𝑤𝑖𝑡ℎ 𝓌 𝑎 = 2.0) at zero system nonlinearity (𝛤 ⋅ 𝓌 𝑎 = 0) . The left-panel-shown densities distribution (𝒫 𝐼 (𝐼)) of the oscillator array is a decaying exponent [17] , corresponding to an infinite system temperature (cf. Eqs. (25) and (27) ). The right panel shows the position of the system on the DNLSE phase diagram. Note the merging of the 𝒽 𝑎∞ (𝓌 𝑎 ) line with the 𝐿 𝑖 (𝓌 𝑎 ) line (Eq. (2) ). Figure 8:
Equilibrium distribution of site-densities for a weak system nonlinearity (𝛤 ⋅ 𝓌 𝑎 = 0.6) and a very cold system (very close to the ground state line, see the inset). The simulated 𝒫 𝐼 (𝐼) is shown by the cyan-colored dots. Three continuous colored curves are overlaid. Blue – the exact theoretically predicted 𝒫 𝐼 (𝐼) curve for the specific point on the phase diagram [(𝓌 𝑎 , 𝒽 𝑎 ) = (0.6, −0.95)] at 𝛤 = 1 . Red – a
18 - perfect-Gaussian least-mean-squared fitted to the blue curve. The figure shows then that the equilibrium PDF (𝐼) is NOT everywhere on the phase diagram a perfect, possibly-truncated, Gaussian. Or, to make a positive statement, at weak nonlinearities and low temperatures, the shape of the distribution of DNLSE site-densities may deviate from a perfect possibly-truncated Gaussian. At strong nonlinearities, the equilibrium distribution of site-densities everywhere on the thermalization zone of the DNLSE phase diagram is a pure possibly-truncated Gaussian ( [19] and section ). Light green – a Gaussian with variance equals to the temperature of the system { 𝑇 𝐷𝑁𝐿𝑆𝐸 (𝓌 𝑎 , 𝒽 𝑎 ) = (𝛽(𝓌 𝑎 , 𝒽 𝑎 ) ⋅ 𝛤) −1 , Eq. (24) }. At weak nonlinearities then, system temperature may be significantly larger than the `variance` of the distribution of site-densities (see also Figure 9 ). At weak system nonlinearities, the shape of PDFs (𝐼) deviate from a pure Gaussian shape. In addition, the value of system temperature no longer coincides with the `variance` (the variance of an untruncated Gaussian) of the PDF (𝐼) , but rather deviates to the high side. This is shown by
Figure 8 and by
Figure 9 . Figure 9:
Equilibrium distributions of site-densities ( 𝒫 𝐼 (𝐼) ) for weak nonlinearities and very cold systems, all at system density of 𝓌 𝑎 = 1 . Value of the nonlinearity coefficient for [(A,D),(B,E),(C,F)] is sequentially increased as 𝛤 = (0.3,1.0,3.0) respectively. Each panel of the top row shows four curves: magenta – simulated-to-equilibrium PDF (𝐼) . Blue – exact theoretically predicted equilibrium PDF (𝐼) . Dashed light-blue – fit of a Gaussian curve to the theoretical PDF (𝐼) curve (very small deviation in A and close match in B and C). Light green - a Gaussian with variance equal to the temperature of the system (Eq. (24) ). Comparing the green curves to the other curves of each panel, going from A to C, we see how the large difference of temperature vs. variance of the PDF (𝐼) at low system nonlinearity gradually shrinks as the value of the nonlinearity coefficient is increased. The panels of the bottom row show the position of the analyzed systems on the respective phase diagrams.
The last figure of this temperature section -
Figure 10 – presents two maps of log-temperatures for DNLSE systems on the low-nonlinearity portion of the thermalization
19 - zone (
𝛤 ⋅ 𝓌 a of order unity). Shown temperatures were calculated analytically (Eq. (25) ). The two maps presented here are complementary to a temperature map for a stronger nonlinearity portion of the thermalization zone, presented in an already-published study [19] . Figure 10:
Maps of log temperatures on a weak nonlinearity portion (system nonlinearity of order unity) of the thermalization zone. Left:
𝛤 = 0.3 . Right:
𝛤 = 1.0 . The maps were calculated analytically (Eq. (25) ). The yellow lines crossing a green area on both maps mark the respective 𝐿 𝑖 lines. The dark lines crossing along the top of the green area on the left map and crossing the blue area on the right map are isotherm lines, having – as shown – a concave shape. These low nonlinearity maps complement a temperature map for a stronger nonlinearity portion (system nonlinearity greater then order unity) of the thermalization zone presented in an already-published study [19] . With Figure 10 our discussion of system temperatures is concluded. We proceed now to discuss the associated chemical potentials of these equilibrium-reached DNLSE systems. Chemical potentials of equilibrium-reached DNLSE systems
The multiplier 𝜇 in the grand canonical partition function (Eq. (8) ) is introduced in analogy with a chemical potential to ensure conservation of total system density (𝒲 𝑎 ) [7] . Note that here, different from the standard thermodynamic chemical potential that has a dimension of physical energy [25] , the equivalent chemical potential introduced in Eq. (8) is dimensionless. Back to conservation of total system density - indeed, we calculate the value of the density-conserving chemical potential (𝜇 𝑒𝑥𝑝 (𝓌 𝑎 , 𝒽 𝑎 )) for a system at (𝓌 𝑎 , 𝒽 𝑎 ) by solving the one parameter implicit integral equation (26) to equate the known system density (𝓌 𝑎 ) with its computed average (〈𝓌 𝑎 〉) . Equation (26) is solved with a known 𝛽 𝑒𝑥𝑝 , as determined by the analytic expression (25) (and applying Eq. (24) ). We thus designate the so-calculated chemical potential by the subscript “ 𝑒𝑥𝑝 ”.
20 -
The next three figures highlight key characteristics of the chemical potential of systems initially excited onto the thermalization zone of the DNLSE phase diagram and evolved to equilibrium.
Figure 11:
Chemical potential (𝜇 𝑒𝑥𝑝 ) as a function of system energy for a fixed system density (𝓌 𝑎 = 2.0) . Colored vertical lines as in Figure 4 . Two curves are shown, corresponding to two values of the nonlinearity parameter as marked. Somewhat like temperature, values of the chemical potential rise sharply as system energy (𝒽 𝑎 ) gets closer to the infinite temperature line (𝒽 𝑎∞ ) . Like the value of system temperature, the value of the chemical potential vary strongly with the value of the nonlinearity parameter (𝛤) . Interestingly, unlike temperature curves, the rise of the chemical potential curve with system energy is not monotonic. At small system energies, close to the ground state energy, chemical potential values are shown to decrease with system energy before the onset of a monotonic rise. Figure 12:
Chemical potential (𝜇 𝑒𝑥𝑝 ) as a function of site-averaged system density for four 𝑞 ℎ -values (the parameter 𝑞 ℎ designates a fraction of the energy span - Δℎ 𝑎 , i.e.: 𝒽 𝑎 (𝓌 𝑎 ) = 𝒽 𝑎0 (𝓌 𝑎 ) + 𝑞 ℎ ⋅ Δℎ 𝑎 (𝓌 𝑎 ) ). Bottom to top - ( 𝑞 ℎ = 0.1,0.5,0.9,0.95 ) as marked. Left: 𝛤 = 0.3 . Right:
𝛤 = 1.0 . The blue dots were numerically calculated by solving the one parameter implicit integral equation (26) ). The fitted straight green lines, each coincides with the first and last data point of the respective 𝑞 ℎ curve. It is quite illuminating to realize that for a fixed 𝑞 ℎ value, the values of the chemical potentials of DNLSE systems at equilibrium essentially fall on a straight line.
21 -
The curves of
Figure 11 show the dependence of equilibrium chemical potential on system energy for fixed system density as marked (going up vertically on the phase diagram). The two sets of curves in each of the panels were calculated for two different values of the nonlinearity parameter (𝛤 = 0.3 ; 1.0) . Moving up in energies, the curves show fast rise of chemical potential values for high system energies, but, interestingly, show slow fall of values for system energies close to the ground state energy (𝒽 𝑎0 ) . The two panels of Figure 12 reveal an interesting dependence of the DNLSE chemical potential on system density. The panels show that chemical potential values plotted against system density for systems with energies at a fixed ratio (𝑞 ℎ ) of the total energy span (𝒽 𝑎 (𝓌 𝑎 ) = 𝒽 𝑎0 (𝓌 𝑎 ) + 𝑞 ℎ ⋅ Δℎ 𝑎 (𝓌 𝑎 )) form a straight line. Not surprising, starting value (for a near-zero density) and slope of a straight line depend on the value of the nonlinearity parameter (cf. left panel vs. right panel). Such regular dependence could provide a clue for deriving an independent analytic expression to calculate equilibrium chemical potentials of DNLSE systems (not done in the present study). Figure 13:
Maps of chemical potential (𝜇 𝑒𝑥𝑝 ) shown on the thermalization zone of the DNLSE phase diagram. The green line crossing the map designates the 𝐿 𝑖 line. The concave dark line crossing each of the maps is an equipotential line. The last figure of this section -
Figure 13 , presents two maps of DNLSE chemical potentials. Note that, different from the temperature maps (
Figure 10 ), the color scale is linear. Looking at the maps, we see that values of DNLSE chemical potentials, like temperature values, sharply rise with energies if the systems are initially excited onto a high energy region of the thermalization zone. With
Figure 13 our discussion of chemical potentials is concluded. We proceed now to discuss the important dynamic phenomenon of evolving DNLSE systems – field correlation.
22 - Field correlations of DNLSE systems
A key characteristic of evolving DNLSE systems is the formation of field correlations. The formation of correlations is perhaps best viewed in terms of energy diffusion between interaction energy and tunneling energy, keeping the sum constant. Consider for example a system initially excited onto the 𝐿 𝑖 line with uniform amplitudes, say 𝑢 𝑚 = 𝑢 for all 𝑚 and random phases in the full range of [0,2 ⋅ 𝜋) . On system launch, the tunneling energy is zero - 𝒽 (0) = 0 , and interaction energy is at a minimum of 𝒽 (0) = ⋅ 𝛤 ⋅ 𝑢 . As the system evolves, amplitudes spread, amplitude-related entropy is generated, and the interaction energy goes to higher values. Necessarily, to keep the sum of energies constant, tunneling energy must go to negative values, forcing the relative phases of neighboring sites to pile-up close to 𝜋 (to get negative cos 𝜃′𝑠 as signC is positive) and thus field correlations are formed. We realize then that field correlations are intimately related to phase coherence. In fact, for nearly equal amplitudes, the normalized DNLSE field correlations ( 〈𝐶 𝑘 〉/〈𝐶 〉 , see below the defining mathematical expressions) coincide with phase correlations (〈cos 𝜃〉) , much like the phase coherence of a Bose–Einstein condensate in a lattice potential [26] . The level of the formed field correlations can be measured. In trapped ultracold atoms experiments, the degree of phase coherence is directly related to the (measurable) visibility of fringes in the interference patterns formed after expansion of Bose–Einstein condensed atomic clouds [27] . Positive Hamiltonian systems entertain two ground states (both with energy 𝒽 𝑎0 ), corresponding to the two values of signC (cf. Eq. (3) ). Correlation distance (𝑘) of the field functions of each of these two ground states extends to infinity. Absolute value of the normalized field correlations at all site-separations is unity (|𝐶 𝑘𝑁 | 𝑔𝑟𝑜𝑢𝑛𝑑 𝑠𝑡𝑎𝑡𝑒 =1.0) [16] . Most interestingly, for high enough system nonlinearity, the field correlation function and the distribution of phases assume universal forms, independent of the exact value of the nonlinearity parameter and of system density. And equilibrium field correlations can be interrogated experimentally by the study of the equivalent optical system [17] . The current “field correlations” section consists of two parts. In the first part we present field-correlations-related expressions. The second part is devoted to graphical illustrations. Expressions to calculate field correlations . Correlation of fields separated 𝑘 -sites apart, is defined as [17] - 𝐶 𝑘 (𝑧) = 1𝑁 ⋅ ∑ 𝑢 𝑚 (𝑧) ⋅ 𝑁𝑚=1 𝑢 𝑚+𝑘 (𝑧) ⋅ 𝑐𝑜𝑠[𝜃 𝑚,𝑘 (𝑧)] 𝜃 𝑚,𝑘 (𝑧) ≡ 𝜙 𝑚 (𝑧) − 𝜙 𝑚+𝑘 (𝑧) (28) Let us first consider equilibrium nearest neighbor correlation. It follows from Eqs. (6) and (28) that
23 - 𝐶 = 𝓌 𝑎 ; 𝐶 (𝑧) = 𝑠𝑖𝑔𝑛𝐶 ⋅ ℋ (𝑧)2 ⋅ 𝑁 = 𝒽 (𝑧)2 (29) and 𝐶 (𝑧) ≡ 𝐶 (𝑧)𝐶 = 𝒽 (𝑧)2 ⋅ 𝓌 𝑎 = 𝒽 𝑎 − 𝒽 (𝑧)2 ⋅ 𝓌 𝑎 (30) Equation (30) holds for any distance (𝑧) , including of course long equilibrium distances. If the Lagrange parameters (𝛽, 𝜇) are known, then the equilibrium site-averaged interaction energy (𝒽 (𝛽, 𝜇)) can be calculated as - 𝒽 (𝛽, 𝜇) = 〈ℎ (𝛽, 𝜇)〉 = 12 ⋅ 𝛤 ⋅ ∫ 𝐼 ⋅ 𝒫 𝐼 (𝛽, 𝜇; 𝐼) ⋅ 𝑑𝐼 ∞0 (31) with 𝒫 𝐼 (𝛽, 𝜇; 𝐼) given by Eq. (11) . Thus, at equilibrium, the value of nearest neighbors field correlation of DNLSE systems can be determined through finding the equilibrium interaction energy 𝒽 (𝛽, 𝜇) (cf. equations (30) , (31) , and (11) ). Alternatively, if 𝐶 is independently known, as prescribed for example in the next paragraph, then the values of both equilibrium tunneling energy (𝒽 (𝑧 𝑒𝑞𝑢𝑖𝑙 . )) and equilibrium interaction energy (𝒽 (𝑧 𝑒𝑞𝑢𝑖𝑙 . )) can be trivially determined. A second, more general approach for calculating equilibrium field correlations at any site-separation is through the value of 〈cos 𝜃〉 , 𝜃 being fields’ relative phase-angle. The Expression for field correlation at site-separation 𝑘 (〈𝐶 𝑘 〉) was derived by the authors of [18] as 〈𝐶 〉 = 〈𝐼〉 = 〈𝓌 𝑎 〉 = 𝓌 𝑎 〈𝐶 𝑘 〉 = 〈√𝐼〉 ⋅ 〈𝑐𝑜𝑠(𝜃)〉 𝑘 ; 𝑘 ≥ 1 (32) Elaborating on Eq. (32) , the authors of [18] say that the derived correlation equation is general, not limited to initial statistical-Gaussian excitation, and holds for field correlations in systems evolving under DNLSE dynamics for any long distance fields distribution where the site-densities ( 𝐼 𝑚 ’s) are not correlated and phase differences ( 𝜃 𝑚 ’s) are random (over realizations). Equation (32) shows that if the phase-differences are not flat-distributed (such that 〈𝑐𝑜𝑠(𝜃)〉 ≠ 0 ), then the fields are correlated and correlations exponentially decay with site-separation (𝑘) . If 〈𝑐𝑜𝑠(𝜃)〉 is negative, the sign of the fields’ correlations alternates with 𝑘 . To determine 〈𝑐𝑜𝑠(𝜃)〉 , the relative phase-angle density distribution (𝒫 𝜃 (𝜃)) must be known. An analytic expression for 𝒫 𝜃 (𝜃) in the hot zone near the 𝐿 𝑖 line - 𝒫 𝜃 (𝜃) ∝exp( ⋅ 𝑐𝑜𝑠 𝜃 ) – was derived in [18] , and an analytic expression for 𝒫 𝜃 (𝜃) in the
24 - entire cold zone - 𝒫 𝜃 (𝜃) ∝ exp( 𝑎 ⋅ 𝑐𝑜𝑠 𝜃 ) – was derived in [19] . The parameter 𝜂 is a Lagrange parameter with a value determined by solving an appropriate implicit equation. The parameter 𝑀 is a first-moment of an (assumed) initial statistical-Gaussian excitation. These PDFs (θ) were derived for strong system nonlinearity (𝛤 ⋅ 𝓌 𝑎 ≫ 1) justifying the quantum phase model approximation [28] , [17] under which the analytic PDFs (𝜃) expressions were derived. Calculating 〈𝑐𝑜𝑠(𝜃)〉 according to the known PDFs (𝜃) and inserting in (32) , we arrive at an analytic expression for equilibrium field correlations, at any site-separation (𝑘) for systems near and below the 𝐿 𝑖 line on the strong nonlinearity portion of the DNLSE thermalization zone: 𝐶 𝑘 = 〈𝐶 𝑘 〉 = 〈√𝐼〉 ⋅ { [( −1 ) ⋅ ℐ ( 𝓌 𝑎 ) ℐ ( 𝓌 𝑎 )] 𝑘 ; 𝒽 𝑎 (𝓌 𝑎 ) ≤ 𝐿 𝑖 (𝓌 𝑎 )[( −1 ) ⋅ ℐ ( ) ℐ ( )] 𝑘 ; 𝒽 𝑎 (𝓌 𝑎 ) > 𝐿 𝑖 (𝓌 𝑎 ) } ; 𝑘 ≥ 1 𝐶 𝑘𝑁 ≡ 〈𝐶 𝑘 〉〈𝐶 〉 = 𝐶 𝑘 𝓌 𝑎 (33) were ℐ 𝑛 (𝑧) is the modified Bessel function of the first kind. The field correlation curves and map in the illustrating figures shown below were calculated either by direct running of system-evolution simulations and applying the definition (28) , or by executing the analytic expression (33) . Graphical illustrations . The next four figures show field correlation values for systems on different parts of the DNLSE thermalization zone, along with graphical representations of the universality of field correlations formed during the evolution of DNLSE systems. Start with
Figure 14 . The red dots in all eight panels were calculated by evolution-to-equilibrium simulations and executing Eq. (28) . The red curves of all panels show high nearest-neighbors field correlations (absolute value of) for cold systems (between the vertical blue and vertical green lines) and essentially zero field correlations for hot systems close to the vertical red 𝒽 𝑎∞ line. The blue dots were calculated analytically through Eq. (33) . The blue curves show that the analytic expression (33) predicts the value of the field correlations reasonably well for high absolute correlation values or at high system nonlinearity. Consulting the top three right panels and the bottom two right panels, we see that the red curves cross the green line at essentially 𝐶 = −0.4 , indicating universality of DNLSE field correlations at high system nonlinearity [17] , [18] .
25 -
Figure 14:
Normalized nearest-neighbors equilibrium field correlations (𝐶 ) as a function of site-averaged system energy at two different values of the nonlinearity parameter (𝛤) . The red dots were obtained by simulations (averaged over several realizations) through a focusing DNLSE equation (Eq. (3) ). The blue dots were calculated analytically (Eq. (33) ). The blue dots are missing from the two left panels as the quantum phase approximation is invalid for such low system nonlinearities (𝛤 ⋅ 𝓌 𝑎 ~1) . Colored vertical lines as in Figure 4 . Clearly, normalized field correlations at all parameters are maximal (in absolute value) (𝐶 = −1) at the ground state, and are all monotonically reduced as system energy (and system temperature) goes up, reaching a zero value at the infinite temperature line. Figure 15:
Absolute value of nearest-neighbors field correlation vs. system temperature. The red dots were calculated through evolution-to-equilibrium simulations and executing Eq. (28) . System density is taken to be constant at 𝓌 𝑎 =16.0 . System energy (𝒽 𝑎 ) is a parameter, i.e. the data points of the curve are on the plane (|𝐶 (𝒽 𝑎 )|, 𝑇 𝑒𝑥𝑝 (𝒽 𝑎 )) with system energy going straight up (on the DNLSE phase diagram) from the 𝒽 𝑎0 line towards the 𝒽 line. The vertical green line designates the temperature of a system on the 𝐿 𝑖 line. The figure clearly shows the inverse relations between field correlation and system temperature. Explicitly – higher system temperature lower absolute value of the system’s field correlation.
26 -
The curve of
Figure 15 indicates the inverse relations of system temperature and field correlation (in absolute value). Explicitly – higher system temperature, weaker field correlations. Remembering that temperature is proportional to the width of PDF (𝐼) , we can make a parallel inverse-relation statement: the width of equilibrium PDF (𝐼) is inversely related to field correlation. Explicitly: wider PDF (𝐼) , weaker field correlation. Next,
Figure 16 shows the remarkable universality property of equilibrium field correlations formed during evolution of DNLSE systems at high system nonlinearities (𝛤 ⋅ 𝓌 𝑎 ≫ 1) . The straight horizontal lines of Figure 16 , each of a different energy, extend the correlation’s universality findings of [17] , and of [18] . The data points of the figure were calculated through Eq. (33) . As shown very clearly, at high system nonlinearity, normalized equilibrium values of field correlations of all DNLSE systems with energies 𝒽 𝑎 (𝑞 𝑔 ; 𝓌 𝑎 ) = 𝒽 𝑎0 (𝓌 𝑎 ) + 𝑞 𝑔 ⋅ 2 ⋅ 𝓌 𝑎 are practically equal, independent of the value of the nonlinearity parameter (𝛤) and independent of system density (𝓌 𝑎 ) . This correlation-level equality extends to each and every site-separation (𝑘) . Figure 16:
Universality of the DNLSE equilibrium field correlations. Shown are normalized values (Eq. (33) ) of the nearest-neighbors ( 𝐶 (𝑧 𝑒𝑞𝑢𝑖𝑙 ) − left) and next-nearest-neighbors ( 𝐶 (𝑧 𝑒𝑞𝑢𝑖𝑙 ) - right) field correlations as a function of the site-averaged density (𝓌 𝑎 ) and the nonlinearity parameter (𝛤 = 0.5,1.0,2.0) . The parameter (𝑞 𝑔 ) is a multiplier of the energy distance between the ground-state and the 𝐿 𝑖 line (i.e. 𝒽 𝑎 (𝑞 𝑔 ; 𝓌 𝑎 ) = 𝒽 𝑎0 (𝓌 𝑎 ) + 𝑞 𝑔 ⋅ 2 ⋅ 𝓌 𝑎 ). The remarkable universality characteristics of the field correlations in equilibrated DNLSE systems is clearly revealed. Namely, at strong enough system nonlinearity (𝛤 ⋅ 𝓌 𝑎 ≫ 1) , the values of every site-separation (𝑘) field correlation is independent of system density and of the nonlinearity parameter (𝛤) at a fixed 𝑞 𝑔 energy level across the entire DNLSE thermalization zone. Along the 𝐿 𝑖 line for example (𝑞 𝑔 = 1.0) , the value of 𝐶 is ±≈0.455 , and the value of 𝐶 is +0.215 , independent of the value of the nonlinearity parameter (𝛤) . The shown several straight horizontal correlation lines of each panel extend the correlation’s universality findings of [17] , and of [18] . The next figure -
Figure 17 – compares the characteristics of a linear system (𝛤 = 0) and the characteristics of a nonlinear DNLSE system. Both systems are initially excited onto the lower part of the cold zone, near the 𝒽 𝑎0 line (panels D,H). Key differences
27 - are seen in the PDFs (𝐼) , panels A,E and in the constant vs. decaying field correlations in panels C,G. In the extreme case of system nonlinearity approaching zero (essentially by
𝛤 → 0 ), system temperatures rise to infinity everywhere on the triangular-shaped phase diagram, with the exception of the zero-variance non-energy-diffusing ground state. Stated differently: a DNLSE system above ground state is formally approaching an infinite temperature as the nonlinearity parameter is approaching zero, yet except for systems on the 𝐿 𝑖 line (cf. Figure 7 ), the PDFs (𝐼) are of a finite variance (panel A of
Figure 17 ). Figure 17:
Characteristics of an equilibrated linear system (top row,
𝛤 = 0 ) vs. an equilibrated nonlinear system (bottom row,
𝛤 = 1.0 ). A,E: Simulated PDF (𝐼) . As shown, the shape of the 𝒫 𝐼 (𝐼) curve for the linear system does not resemble a Gaussian at all. Note the high amplitude oscillators of the stationary-statistics-reached linear system (A). B,F: PDF (𝜃) . The magenta dots are simulated. The continuous green curve is theoretical [19] . In both linear-nonlinear cases, relative phase-angles (𝜃 𝑚 = 𝜙 𝑚 − 𝜙 𝑚+1 ) are narrowly concentrated around 𝜋 , as expected for systems initially excited close to the ground state and evolving under the dynamics a focusing-type DNSLE equation. C,G: Simulated field correlations. Considering the linear system, if the system is initially excited with uniform amplitudes and with a small random-phase window, the launched fields are strongly correlated, independent of site-separation (𝑘) . These k-independent strong correlations are maintained during evolution of the single-term-Hamiltonian linear system (C). Unlike field correlations of the linear system, during evolution of the two-term-Hamiltonian nonlinear system (𝛤 > 0) , the initially-formed strong field correlations do gradually weaken (for 𝑘 ≥2) to reach a 𝑘 -dependent exponential decay at steady state distances (independent of nonlinearity value) (panel G). D,H: System position on the DNLSE phase diagram (𝓌 𝑎 = 1.75) . Looking at the panel for the linear system (D), note the coincidence of the infinite temperature line (𝒽 𝑎∞ ) with the intermediate (𝐿 𝑖 ) line, forming a triangular-shaped thermalization zone. In terms of temperature, the equilibrated linear
28 - system (and all other systems in the linear thermalization zone, except for ground -state-excited systems) is formally at infinite temperature. In contrast, the nonlinear system, initially excited near the ground state line is quite “cold”.
Figure 18:
Map of DNLSE equilibrium nearest neighbors field correlations (𝐶 ) on a strong system nonlinearity portion (𝛤 ⋅ 𝓌 𝑎 ≫ 1) of the thermalization zone. Same-color points are all at a distance of 𝑞 𝑔 (i.e. at a distance of 𝑞 𝑔 ⋅ 2 ⋅ 𝓌 𝑎 ) from the 𝒽 𝑎0 line (or “the ground-state line”, designated by the bottom orange line). The black line crossing the colored area is at a distance of 𝑞 𝑔 = 1.0 . The yellow line at the top of the colored portion is at a distance of 𝑞 𝑔 = 2.0 , showing a value of 𝐶 = −0.2 (cf. Figure 16 ). Thus, here too, universality of the DNLSE equilibrium field correlations is clearly revealed. The colored portion was calculated by the analytic expression (Eq. (33) ). The black area was left untouched. The ten isolated discrete nearly-red-colored dots are at a distance of 𝑞 𝑔 = 4.0 . Their 𝐶 values, ranging between −0.058 and - , were determined by field-propagation simulations (averaged over several realizations). The map area above these discrete point should have been covered entirely by a red color. Namely - essentially zero field correlations at this high temperature region of the DNLSE phase diagram. The last figure of this section -
Figure 18 , is a map of normalized equilibrium field correlations. The map shows a relatively strong nonlinearity region of the thermalization zone. Values for the colored part were determined analytically (Eq. (33) . Values for the discrete red dots by evolution-to-equilibrium simulations. The map area above these discrete point should have been covered red, signifying no correlations. Looking at the map, two effects become clear. First, once more – the universality of DNLSE field correlations at high system nonlinearity as discussed above. Second - the fact of essentially no field correlations at the upper region of the hot zone. Namely - no field correlations at a high temperature-high system nonlinearity region of the DNLSE phase diagram. In the next section, along with a complementary appendix (
Appendix 1 ), we derive an approximate analytic expression to the grand-canonical partition function (Eq. (9) ) and to the associated distribution function of site densities (PDF (𝐼) – Eq. (11) ). We show
29 - that the so-derived PDF (𝐼) well approximates the exact PDF (𝐼) on a large region of the DNLSE thermalization zone. Analytic PDF (𝑰) expression for systems on the hot zone
The symmetric nonnegative kernel related to the grand canonical partition function of large DNLSE systems is the product of two decaying exponential functions and a zero order modified Bessel function of the first kind (Eq. (9) ). In the “second-order” approximation, we represent the modified Bessel function by the first two terms of its series expansion (cf.
Appendix 1 below). Eigenvalues (𝜆 𝑖 ) and nonnegative eigenfunctions (𝑣 𝑖 (𝑥)) of this kernel are the solutions of the Fredholm integral equation of the second kind with symmetric kernel [29] : ∫ [𝐹(𝑥) ⋅ 𝐺(𝑦) ⋅ (1 + 𝛽 ⋅ 𝑥 ⋅ 𝑦)] ⋅ 𝑣 𝑖 (𝑥) ⋅ 𝑑𝑥 = 𝜆 𝑖 ⋅ 𝑣 𝑖 (𝑦) ∞0 𝐹(𝑥) ≡ 𝑒 [−𝛽⋅(𝜇2⋅𝑥+𝛤4⋅𝑥 )] ; 𝐺(𝑦) ≡ 𝑒 [−𝛽⋅(𝜇2⋅𝑦+𝛤4⋅𝑦 )] (34) with the largest eigenvalue (𝜆 ) directly related to the partition function (Eq. (10) ) and the associated eigenfunction (𝑣 (𝑥)) directly related to the PDF (𝐼) of the equilibrated DNLSE system in question (Eq. (11) ). It turns out that if the Bessel function of equation (9) is described by a “second-order” approximation (for arguments smaller than unity), namely by taking the first two terms of the otherwise infinite series describing the function, then the Fredholm integral equation (34) can be solved analytically (cf. Appendix 1 below). The solved expressions for the integral equation with a second-order approximated Bessel function (Eq. (34) ) are (
Appendix 1 ): 𝜆 (𝛽, 𝜇) = 𝑒 𝑞 ⋅ ν ⋅ 1 − 𝑒 𝑞 ⋅ 𝑞 ⋅ ν 𝒫 𝐼𝐸 (𝛽, 𝜇; 𝐼) = 𝒞 −1 ⋅ 𝑒 −(√𝜖⋅𝐼+𝑞 ) ⋅ (1 + 𝜅 ⋅ 𝐼) (35) where 𝜆 is the largest eigenvalue and 𝒫 𝐼𝐸 (𝛽, 𝜇; 𝐼) = 𝑣 𝐼𝐸2 (𝛽, 𝜇; 𝐼) where 𝑣 𝐼𝐸 (𝛽, 𝜇; 𝐼) is the associated eigenfunction. All constants of Eq. (35) are defined in Appendix 1 . Consulting Eq. (35) , we see that the second-order PDF (𝐼) is a product of a Gaussian and a parabola. The parabolic-parameter 𝜅 appearing in the parabolic-term is smaller than the value of the Lagrange parameter 𝛽 and is thus smaller than unity if 𝛽 is smaller than unity, i.e. if system temperature is greater than the value of .
30 -
Figure 19:
Distribution functions of site densities (𝒫 𝐼 (𝐼)) . Four curves are shown in each panel – the continuous cyan curve (1) and the dashed blue curve (3) show the 𝒫 𝐼 (𝐼) of Eq. (11) , calculated with all Bessel terms in place. The continuous thin-dark curve (2) and the dashed magenta curve (4) show the approximate analytic 𝒫 𝐼𝐸 (𝐼) of Eq. (35) , calculated with only two Bessel terms in place. The continuous cyan curve (1) and the continuous thin-dark curve (2) were calculated with exact (𝛽, 𝜇) values determined by numeric inversion of Eq. (15) . The (𝛽, 𝜇) values for the two dashed curves (3 and 4) are actually (𝛽 𝑒𝑥𝑝 , 𝜇 𝑒𝑥𝑝 ) determined through the analytic expression (25) using (24) and through the implicit integral equation (26) . The panels show that second-order PDF(I) of the “cold” system (A) does not match the exact PDF (𝐼) . However, the other two panels (B,C) show that second-order PDFs (𝐼) of Eq. (35) do match the exact PDFs (𝐼) for the two higher energy systems. The three panels of
Figure 19 show the equilibrium PDF (𝐼) of DNLSE systems at three different energies: 𝒽 𝑎 = "𝑐𝑜𝑙𝑑" ; 𝒽 𝑎 = 𝐿 𝑖 ; 𝒽 𝑎 = "ℎ𝑜𝑡" , all at the same density of 𝓌 𝑎 = 6.0 . System nonlinearity is also fixed at 𝛤 ⋅ 𝓌 𝑎 = 6.0 . The parabolic-parameter 𝜅 values for the (A;B;C) panels of Figure 19 are 𝜅 =(0.575 ; 0.038 ; 0.00057) respectively. The panels show that second-order PDF (I) of the “cold” system does not match the exact PDF (𝐼) . However, the other two panels show that second-order PDFs (𝐼) of Eq. (35) do match the exact PDFs (𝐼) for the two higher energy systems. Consulting Eq. (35) and the panels of Figure 19 , we can make the following observations: • Given high system nonlinearity (𝛤 ⋅ 𝓌 𝑎 ≫ 1) , The equilibrium PDFs (𝐼) of systems at or above the 𝐿 𝑖 line are well described by the analytic PDF (𝐼) expression of Eq. (35) . • As the parabolic curves in the analytic PDF (𝐼) expression (Eq. (35) ) for the high energy systems are very shallow, the Gaussian-parabola product is another slightly shifted Gaussian with essentially the same width. Thus, at high system nonlinearity, the shape of the equilibrium PDF (𝐼) curves for DNLSE systems at or above the 𝐿 𝑖 line is practically a pure (possibly truncated) Gaussian shape. • The parameter 𝜖 in the 𝒫 𝐼𝐸 (𝛽, 𝜇; 𝐼) expression of Eq. (35) is related to the (untruncated) Gaussian variance 𝜎 as 𝜖 −1 = 2 ⋅ 𝜎 . Thus, looking at the definition of 𝜖 in Eq. (39) of Appendix 1 (𝜖 ≡ ⋅ 𝛽 ⋅ 𝛤) we find: 𝑇 = (𝛽 ⋅ 𝛤) −1 =𝜎 . Put it all together: given high system nonlinearity, the shape of the equilibrium
31 -
PDF(I) curves for DNLSE systems at or above the 𝐿 𝑖 line is practically a pure (possibly truncated) Gaussian shape with `variance` (𝜎 ) equal to the temperature of the system (Eq. (24) ). These “second-order-related” observations end the present analytic PDF (𝐼) section and actually end our entire detailed characterization of equilibrated large DNLSE systems. Summary
The subject of the study presented here is a system consists of a long array of coupled nonlinear oscillators. The evolution dynamics of the system is taken to be the DNLSE (Eq. (3) ). A DNLSE system initially excited onto a point (𝓌 𝑎 , 𝒽 𝑎 ) of a well-defined thermalization zone of the DNLSE phase diagram (Eq. (2) ) will thermalize [7] (cf. Figure 3 ). The equilibrium statistical properties of a thermalized system are predictable. Previous studies predicted equilibrium PDFs and temperatures for strong system nonlinearities by system-entropy separation into the sum of density-entropy and a relative-phase entropy. Entropy separation was justified in the strong system nonlinearity regime by the quantum phase approximation. Entropy separation led to the derivation of system temperatures as well as the derivation of analytic expressions for equilibrium PDF (𝐼) and equilibrium PDF (𝜃) [17] - [19] . Here we have taken a rigorous approach of predicting system equilibrium statistics based on the thermodynamics formalism of grand canonical ensembles (cf. Figure 2 and related text). The thermodynamics approach allowed the extension of system-characteristics predictions to cover the entire thermalization zone. The grand canonical statistics is determined by two Lagrange parameters - (𝛽, 𝜇) .To determine the values of 𝛽 and 𝜇 we have taken here one of two routes. The first is numerically inverting the two exact thermodynamics partial derivatives of Eq. (15) . Such 2D numerical inversion is challenging and could sometimes introduce uncertainty errors. The second is finding 𝛽 𝑒𝑥𝑝 through execution of the analytic expression (25) , along with expression (24) . We have derived the analytic expression (25) from an approximate analytic expression of system entropy published in [9] . Once the value of 𝛽 𝑒𝑥𝑝 is determined, the value of 𝜇 𝑒𝑥𝑝 is calculated through solving the one parameter implicit integral equation (26) . Once (𝛽(𝓌 𝑎 , 𝒽 𝑎 ), 𝜇(𝓌 𝑎 , 𝒽 𝑎 )) are determined, equilibrium statistical properties of the system studied such as entropy, temperature, chemical potential, nearest-neighbors field correlations, and PDF (𝐼) can be calculated. To derive physics insights from numerical calculations, graphical representations are typically required. Analytic expressions on the other hand often provide deeper physics understanding by direct inspection. As part of the present study we have gone through deriving an approximate analytic PDF (𝐼) expression (section and Appendix 1 ). Probing further the derived expression enabled the insightful identification of system’s temperature with the variance of the system’s PDF (𝐼) (cf. the last part of section ). Crossing the DNLSE thermalization zone is an intermediate 𝐿 𝑖 line (Eq. (2) ). The 𝐿 𝑖 line divides the thermalization zone into a cold zone below the line and a hot zone above the line ( [19] . Characteristics of DNLSE systems can be generalized in relations to these cold-hot zones.
32 -
Some of our deduced characteristics of equilibrated DNLSE systems are the following: • Temperatures of systems on the cold zone rise approximately linearly with energy. Temperatures of systems on the hot zone rise exponentially-like with energy. • On the high system nonlinearity region of the thermalization zone (𝛤 ⋅ 𝓌 𝑎 ≫ 1) the PDF (𝐼) is of a Gaussian shape, possibly truncated. • On the high system nonlinearity region of the thermalization zone, system temperature and the `variance` of its equilibrium PDF (𝐼) are equal. • Chemical potentials of systems on the hot zone rise exponentially-like with energy. • Field correlations are inversely related to system temperature - higher system temperature, lower absolute value of the system’s field correlation. • Field correlations for all site-separations (𝐶 𝑘𝑁 ) of systems on the high system nonlinearity region of the thermalization zone at all energies are universal – their values are independent of system density or of the value of the nonlinearity parameter. At energies of 𝑎 above ground state, absolute value of nearest-neighbors normalized field correlation is already as low as (cf. Figure 18 ). • On the high system nonlinearity region of the thermalization zone, nearest-neighbors field correlations of DNLSE systems with energies (𝒽 𝑎 ) of 𝑎 or more above ground state practically vanish (cf. Figure 18 ). • Given high system nonlinearity, the equilibrium PDF (𝐼) of systems at or above the 𝐿 𝑖 line are well described by the analytic PDF ( 𝐼 ) expression of Eq. (35) . These observations, derived for an abstract DNLSE system, hold of course for the equilibrium properties of the variety of actual discrete nonlinear physical systems with evolution dynamics approximated by a DNLSE type equation as analyzed in the present study. Acknowledgement
I would like to thank Prof. David Mukamel for his many guiding comments and constructive suggestions during preparation of the manuscript.
33 -
Appendix 1:
Approximate distribution function of site densities
In this appendix we derive an approximate analytic expression for the grand canonical partition function of Eq. (9) along with an approximate analytic expression for the distribution function of site densities (𝒫 𝐼 (𝐼)) . We derive the approximate expressions through taking only two first terms of the infinite series by which the zero order modified Bessel function of the first kind is described. We refer to this approximation as a “second-order” approximation. The infinite series describing the zero order modified Bessel function of the first kind ( ℐ (𝑧) , Eq. (13) ) with 𝑧 a complex number, is given by the expression [30] : ℐ (𝑧) = ∑ (14 ⋅ 𝑧 ) 𝑘 (𝑘!) (36) If the argument of the modified Bessel function is small (2 ⋅ 𝛽 ⋅ √𝑥 ⋅ 𝑦 < 1) for (𝑥, 𝑦) in the range of interest (determined by the decay-rate of the preceding exponential functions), then the modified Bessel function can be approximated by the first two terms: ℐ (2 ⋅ 𝛽 ⋅ √𝑥 ⋅ 𝑦) ≅ 1 + 𝛽 ⋅ 𝑥 ⋅ 𝑦 ; 2 ⋅ 𝛽 ⋅ √𝑥 ⋅ 𝑦 < 1 (37) The kernel of Eq. (13) is approximated then as
𝒦(𝑥, 𝑦) ≅ 𝐹(𝑥) ⋅ 𝐺(𝑦) ⋅(1 + 𝛽 ⋅ 𝑥 ⋅ 𝑦) with 𝐹(𝑥) ≡ 𝑒 [−𝛽⋅( 𝜇2 ⋅𝑥+ 𝛤4 ⋅𝑥 )] ; 𝐺(𝑦) ≡ 𝑒 [−𝛽⋅( 𝜇2 ⋅𝑦+ 𝛤4 ⋅𝑦 )] . We are now looking for an eigenvalue (𝜆 ) and an eigenfunction (𝑣 (𝑥)) such that [29] : ∫ [𝐹(𝑥) ⋅ 𝐺(𝑦) ⋅ (1 + 𝛽 ⋅ 𝑥 ⋅ 𝑦)] ⋅ 𝑣 (𝑥) ⋅ 𝑑𝑥 = 𝜆 ⋅ 𝑣 (𝑦) ∞0 (38) The subscript “ 𝐸 ” in (38) stands for “expression”. It is convenient at this point to define a set of four parameters (𝜖; 𝛼; 𝑞 ; 𝜈 ) each depends, directly or indirectly, on the two Lagrange parameters (𝛽, 𝜇) as: 𝜖 ≡ 12 ⋅ 𝛽 ⋅ 𝛤 ; 𝛼 ≡ 𝛽 ⋅ 𝜇 ; 𝑞 ≡ 𝛼2 ⋅ √𝜖 ; 𝜈 ≡ √𝜋 ⋅ 𝑒𝑟𝑓𝑐(𝑞 ) (39) Working out the lengthy mathematics, we first determine the value of a parabolic-parameter - 𝜅 as:
34 - 𝜅 = −𝑏 ± √𝑏 + 4 ⋅ 𝛽
2 ; 𝜅 ≡ 𝜅 ; 𝑏 ≡ 𝐼 − 𝛽 ⋅ 𝐼 𝐼 𝐼 = ν = 𝑒 −q − 𝑞 ⋅ ν 𝐼 = −2 ⋅ 𝑒 −𝑞 ⋅ 𝑞 + (1 + 2 ⋅ 𝑞 ) ⋅ ν (40) And then, with 𝜅 in place, calculate the thought-for expressions for (𝒵(𝛽, 𝜇); 𝜆 (𝛽, 𝜇); 𝒫 𝐼 (𝛽, 𝜇; 𝐼)) as: 𝒵 𝐸 (𝛽, 𝜇) ≅ [2 ⋅ 𝜋 ⋅ 𝜆 (𝛽, 𝜇)] 𝑁 𝜆 (𝛽, 𝜇) = 𝑒 𝑞 ⋅ ν ⋅ 1 − 𝑒 𝑞 ⋅ 𝑞 ⋅ ν 𝒫 𝐼𝐸 (𝛽, 𝜇; 𝐼) = 𝒞 −1 ⋅ 𝑒 −(√𝜖⋅𝐼+𝑞 ) ⋅ (1 + 𝜅 ⋅ 𝐼) 𝒞 = 𝐼 + 2 ⋅ 𝜅 ⋅ 𝐼 + 𝜅 ⋅ 𝐼 (41) where 𝜆 (𝛽, 𝜇) and 𝑣 (𝛽, 𝜇; 𝐼) in 𝒫 𝐼𝐸 (𝛽, 𝜇; 𝐼) = 𝑣 (𝛽, 𝜇; 𝐼) are the solutions of Eq. (38) . As 𝒫 𝐼𝐸 (𝐼) is explicitly known, the statistical averages of system density (𝓌 𝑎𝐸 (𝛽, 𝜇)) , equilibrium interaction energy (𝒽 (𝛽, 𝜇)) , and equilibrium tunneling energy (𝒽 (𝛽, 𝜇)) can be explicitly known as well: 𝐼3 = 𝑒 −𝑞 ⋅ (1 + 𝑞 ) − 𝑞 ⋅ (1.5 + 𝑞 ) ⋅ 𝑣 𝐼4 = −𝑒 −𝑞 ⋅ 2 ⋅ 𝑞 ⋅ (5 + 2 ⋅ 𝑞 ) + (3 + 12 ⋅ 𝑞 + 4 ⋅ 𝑞 ) ⋅ 𝑣 𝓌 𝑎𝐸 (𝛽, 𝜇) = 𝒞 −1 ⋅ (𝐼 + 2 ⋅ 𝜅 ⋅ 𝐼 + 𝜅 ⋅ 𝐼 ) 𝒽 (𝛽, 𝜇) = 𝛤2 ⋅ 𝒞 −1 ⋅ (𝐼 + 2 ⋅ 𝜅 ⋅ 𝐼 + 𝜅 ⋅ 𝐼 ) 𝒽 (𝛽, 𝜇) = 𝒽 𝑎 − 𝒽 (𝛽, 𝜇) (42) As shown by
Figure 19 in the main text, the derived 𝒫 𝐼𝐸 (𝐼) expression of Eq. (41) well approximates the exact PDF (𝐼) of Eq. (11)(11)