Chaos and thermalization in a classical chain of dipoles

PHYSICAL REVIEW E 106, 014213 (2022) Chaos and thermalization in a classical chain of dipoles Manuel Iñarrea ,1Rosario González-Férez ,2J. Pablo Salas,1and Peter Schmelcher3,4 1Área de Física, Universidad de La Rioja, 26006 Logroño, La Rioja, Spain 2Instituto Carlos I de Física Teórica y Computacional, and Departamento de Física Atómica, Molecular y Nuclear, Universidad de Granada, 18071 Granada, Spain 3The Hamburg Center for Ultrafast Imaging, Luruper Chaussee 149, 22761 Hamburg, Germany 4Zentrum für Optische Quantentechnologien, Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany (Received 20 January 2022; accepted 19 May 2022; published 27 July 2022) We explore the connection between chaos, thermalization, and ergodicity in a linear chain of Ninteracting dipoles. Starting from the ground state, and considering chains of different numbers of dipoles, we introduce single site excitations with excess energy K. The time evolution of the chaoticity of the system and the energy localization along the chain is analyzed by computing, up to a very long time, the statistical average of the finite-time Lyapunov exponent λ(t) and the participation ratio (t). For small K, the evolution of λ(t)and (t) indicates that the system becomes chaotic at approximately the same time as (t) reaches a steady state. For the largest considered values of Kthe system becomes chaotic at an extremely early stage in comparison with the energy relaxation times. We find that this fact is due to the presence of chaotic breathers that keep the system far from equipartition and ergodicity. Finally, we show numerically and analytically that the asymptotic values attained by the participation ratio (t) fairly correspond to thermal equilibrium. DOI: 10.1103/PhysRevE.106.014213 I. INTRODUCTION The relationship between chaos, thermalization, and ergodicity in Hamiltonian systems with a large number of degrees of freedom is a topic of intense research with many intriguing open questions. Historically it was the pioneering study of Fermi, Pasta, Ulam, and Tsingou (FPUT) in 1953 [1,2], that initiated and opened up this field of research. Indeed, for a chain of nonlinear oscillators excited out of the equilibrium, FPUT found that the expected energy equipartition was not reached. Instead, they observed quasiperiodic energy recurrences, which are more likely to occur in integrable systems. Today, we know that these recurrences appear because the initial conditions used by FPUT were chosen near time-periodic solutions showing a strong energy localization in the normal mode space (q-breathers)[3–5]. Since then, the possibility that even in weakly nonlinear Hamiltonian systems, thermalization might not occur or be extremely slow due to the spontaneous appearance of nonergodic local fluctuations is a legitimate point of view. In complex Hamiltonian systems, the unpredictable nature of a chaotic orbit might suggest that the corresponding dynamics is ergodic and therefore such an orbit describes a thermalized system. The later implies that the chaotic trajectory is able to explore all of the available phase space. However, the combined results of the Kolmogorov-ArnoldMoser(KAM)[6–9] and the Nekhoroshev [10] theorems state that, in all weakly perturbed integrable systems it is always possible to find orbits that remain trapped close to regular phase space regions up to very long times. Furthermore, the remaining chaotic regions are connected due to Arnold diffusion, which means that, regardless of the time spent, every chaotic orbit will eventually visit every chaotic phase space region. Even though it is commonly accepted that the size of the regular islands (e.g., the KAM regime) vanishes very fast, even exponentially, for increasing number of degrees of freedom [11], ergodicity, and thermalization can only be fully developed in strongly perturbed Hamiltonian systems where there are (almost) no regular islands and phase space is then dominated by global chaos. As a consequence, although chaos always appears as the fundamental precursor of thermalization in nonlinear lattices [12–21], we also know that chaotic behavior is not always a sufficient condition to assert that a given orbit has also reached the thermalization regime [17,19,22,23]. In fact, in Refs. [17,18,20,21,24,25] we can find examples of nonlinear lattices where the time needed by the system to become chaotic is much shorter than the ergodization time. In all those systems, the large difference between the two timescales is due to the presence of chaotic breather-like excitations. In this way, the main goal of this paper is to provide more insight on a fundamental question such as the connection between chaos, thermalization and ergodicity. To this end, we use a linear chain of Nidentical rigid interacting dipoles. Starting with the system in its ground state (GS), a certain amount of energy Kis given to one of the dipoles, and we explore the transport of that excess energy with increasing time evolution. To detect chaos, we compute the maximal Lyapunov exponent λ1[26–28], being the limit t→∞of the finite-time Lyapunov exponent λ(t)=1 tlog w(t) w(0),(1) 2470-0045/2022/106(1)/014213(9) 014213-1 ©2022 American Physical Society MANUEL IÑARREA et al. PHYSICAL REVIEW E 106, 014213 (2022) i+2 θX Z a i+1 θ i θ i_1 θ i_2 θ FIG. 1. Schematic representation of the dipole chain. where w(0) and w(t) are the deviation vector of a given trajectory at t=0 and t>0. The determination of λ(t)is a common tool to estimate whether an orbit is chaotic or not. For a regular orbit it tends to zero as λ(t)∼t−1, while for chaotic orbits it reaches asymptotically a nonzero value. The inverse τ=1/λ1is the so-called Lyapunov time, which quantifies the time needed for the system to become chaotic. To measure the degree of equipartition of the initial excitation K, we use the (normalized) participation ratio (t)[29,30] (t)=1 N−1K2 N k=1Ek(t)2−1,(2) with Ek(t) being the local energy stored in each dipole, that will be defined below. The energies Kand Ekare taken with respect to the GS energy. When the excitation is completely localized, carried by a single dipole, the value of (t) is zero, while if there is complete equipartition (t)=1. When the system at hand is ergodic and fully chaotic, the equilibrium value of (t) can be computed as its average using the Boltzmann distribution. The paper is organized as follows. In Sec. II we present the Hamiltonian of the system, and a short discussion of the most relevant equilibrium configurations is provided. Starting from the GS, in Sec. III we follow the long-term time evolution of an initially localized excitation. In particular, we compute the averaged finite-time Lyapunov exponent λ(t) and the averaged participation ratio (t) for different excitation values. These quantities provide information concerning the chaoticity of the system, as well as the degree of thermalization of the dipole chain. In Sec. IV, we apply the Boltzmann statistics to determine the thermal equilibrium values of (t). We find a very good agreement between the thermal equilibrium values of (t) and the numerically calculated asymptotic values of the averaged (t), which indicates that the system has approximately reached thermal equilibrium. II. HAMILTONIAN AND EQUILIBRIUM CONFIGURATIONS The dipoles are fixed in space along the X-axis of the Laboratory Fixed Frame XYZ with a distance abetween two consecutive dipoles. They are restricted to rotate in the common XZ plane (see Fig. 1). Thence, the dipole moment of each rotor is given by the vector μi=μo(cos θi,0,sin θi), where 0 ⩽θi<2πis the angle between the dipole moment μiand the Xaxis, with i=1,...,N. Assuming periodic boundary conditions (PBC) and only interactions between nearest neighbors, the rotational dynamics of the system, as a function of the phases θi, is described by the following dimensionless Hamiltonian: H= N  i=1 p2 i 2+V(θ1,...,θ N),(3) where pi=dθi/dt and V=V(θ1,...,θ N) is the potential energy surface of the chain given by V= N  i=1 Vi= N  i=1 (sin θisin θi+1−2 cos θicos θi+1).(4) The energy E=Hin Eq. (3) is measured in units of Bχ, where B=¯h2/2Iis the molecular rotational constant of the dipoles, and χ=μ2 o/(4π0a3B) is the dimensionless dipole-dipole interaction parameter in units of B. In this formulation, the new dimensionless time is t=√χt/tBwith tB=¯h2/√2B. For more information about this reduction, we refer the reader to Ref. [30]. The GS of the system corresponds to the so-called head-tail configuration {θi=0,∀i}or {θi=π,∀i}. The minimal energy of these equilibria is Em=−2N. Besides the GS configuration, the resulting Hamiltonian equations of motion provide us with two families of equilibria that give rise to a complex choreography of equilibrium configurations. One of the families is made of alternating blocks of arbitrary number of dipoles, where all dipoles belonging to the same block are either oriented with angles 0 or π. The other family is also made of alternating blocks of arbitrary number of dipoles, but now with all dipoles belonging to the same block either oriented with angles π/2or−π/2. Determining the nature of the equilibrium configurations involves obtaining the eigenvalues of the stability matrix associated to Eq. (3) and it has been achieved in Ref. [30]. Here, we focus on the degenerate set of equilibria given by only one dipole flipped with respect to the GS configuration. Naming these equilibria as S, their energy is Es=8−2N, and they are saddle points. They are the equilibria with the closest energy to the GS. Furthermore, the energy gap between the GS and S is s=Es−Em=8, which does not depend on the chain size. These equilibria S play a very important role in the dynamics because, for energy values below Es, the phase space trajectories of the system remain trapped around the GS. Conversely, for E>Es, larger phase space regions are accessible for the trajectories, which involve also different equilibria. Then, a stronger nonlinear dynamics is expected to take place. Hereafter the energy of the system will be taken with respect to the GS energy, i.e., the total energy of the system will be shifted by Em=−2N. III. EXCITATION DYNAMICS As we mentioned before, starting form the head-tail configuration of minimal energy Em, we excite at t=0 a single dipole with an excess energy K. We use chains between N=100 and N=400 dipoles. Because PBC are assumed, without loss of generality, we excite the central dipole N/2of the lattice. Then, the initial conditions (i.c.) of the system are θi(0) =pi(0) =0,∀i= N/2, K=pN/2(0)2 2+4[1 −cos θN/2(0)].(5) 014213-2 CHAOS AND THERMALIZATION IN A CLASSICAL CHAIN … PHYSICAL REVIEW E 106, 014213 (2022) Because the energy gap s=8 between the GS and the saddle point configurations S does not depend on the chain size, we provide the Kvalues in terms of that gap. Then, because sis N-independent, for a given value of K,the larger the system’s size is, the smaller the energy per dipole (energy density) =K/Nis. In general, the influence of the system’s size on the dynamics has been studied keeping  constant and varying N(see, e.g., Refs. [13,31]). For particular values of K, we estimate λ(t) and (t) by the simultaneous numerical integration of the Hamiltonian equations of motion arising from Eq. (3) and the corresponding variational equations. More specifically, for each value of K,λ(t) and (t) are statistically determined by averaging over 20 different realizations compatible with the i.c. (5). For the numerical integration of the equations of motions we used the SABA2symplectic integrator [17,32] with fixed integration time step. On the other side, the numerical integration of the variational equations has been performed using a tangent map extension [33–35] of the SABA2algorithm. We use an integration time step h=0.1 which keeps the relative energy error less than 10−4. Because the computations of λ(t) [and so (t)] require very long integration times, the code has been parallelized. To define the local energies Ek(t) appearing in Eq. (2), we divide the potential energy in the Hamiltonian (3) between neighboring sites equally upon the involved sites (see, e.g., Refs. [30,36]), such that Ek(t)=pk(t)2 2+Uk(θk−1,θ k,θ k+1),(6) where the local potential Uk(θk−1,θ k,θ k+1) is given by Uk(θk−1,θ k,θ k+1)=1 2{sin θk(t)[sin θk+1(t)+sin θk−1(t)] −2 cos θk(t)[cos θk+1(t) +cos θk−1(t)]}+2.(7) We have performed simulations for six excitations with excess energy Kbelow the energy gap s, and two excitations with Kabove s, namely, for K=0.05s,0.1s, 0.25s,0.5s,0.75s,0.9s,1.1s, and 1.25s. Note that, when these excitations are below the energy gap s,the dipoles cannot perform complete rotations. For the computation of λ(t) we need to choose an initial deviation vector w(0). Although all choices of w(0) will converge to the same maximal Lyapunov exponent λ1, not all of them will do it at the same rate. A thorough discussion with regard to this subject can be found in Refs. [37,38]. Because a random choice of w(0) is very common, in our case we take w(0) as a 2Nnormalized vector with uniformly distributed random entries. A. Chaotic behavior of the dipole chain For fixed values of Nand K, the averaged finite-time Lyapunov exponent  λ(t) will be statistically estimated by averaging over the 20 values {λi(t),i=1, .., 20}, where λi(t) is the finite-time Lyapunov exponent of an orbit with initial conditions compatible with Eq. (5). As examples, the averaged finite-time Lyapunov exponent λ(t) for three chains with N=100, N=200,and N=400 and for the above excitations are shown in Figs. 2(a)–2(c) on a double logarithmic (a) N=100 0 123456789 -1 -2 -3 -4 -5 -6 (c) N=400 0 1 23456789 -1 -2 -3 -4 -5 -6 123456789 0 -1 -2 -3 -4 -5 (d) N=100 (b) N=200 0 123456789 -1 -2 -3 -4 -5 -6 FIG. 2. (a)–(c) Computed average  λ(t) for eight excess energies along with some error bars (standard deviation) for chains with N= 100, N=200,and N=400 dipoles, respectively. The red dashed line guides the eye for the −1 slope, which marks the expected λ(t)∼t−1behavior for regular orbits. (d) Time evolution of the finite-time Lyapunov exponent of each of the 20 initial conditions (gray lines), the corresponding averaged finite-time Lyapunov exponent  λ(t) (color lines) and some error bars for a chain of N=100 dipoles and for the excess energies K=0.05s,0.25s,0.75s, and 1.25s. scale. Some error bars computed as the standard deviation are also presented in this figure. In all cases, the small size of the error bars indicate that, in the long term,  λ(t) is insensitive to the choice of the initial conditions. This fact can be clearly observed in Fig. 2(d) where we illustrate in detail the computational procedure for  λ(t) for a chain with N=100 dipoles and for K=0.05s,0.25s,0.75s, and 1.25s. For each K,Fig.2(d) shows the finite-time Lyapunov exponent of each of the 20 initial conditions (gray lines), the corresponding averaged finite-time Lyapunov exponent λ(t) (color lines) and some error bars. From a statistical point of view, for each value of K, we observe that the values of λi(t) tend to stabilize around the same limiting value. Assuming that this limiting value is the maximal time Lyapunov exponent λ1, we conclude that, in all cases, the value of λ1do not depend on the choice of the initial conditions. In other words, for a given Kvalue, the chaotic behavior of the chain would be characterized by a single λ1. In all cases, we observe that the time evolution of  λ(t) qualitatively shows always the same behavior. Indeed, after the system is excited, there is a transient during which  λ(t) decreases in time. After that transient, there is a crossover to a plateau, and the corresponding value of the maximal Lyapunov exponent λ1is taken on. However, the timescales in the behavior of  λ(t) are very different and depend on the value of K. The larger the excess energy Kis, the shorter the transient is, and the larger the value λ1is. This hierarchy in the decay patterns of λ(t) (and so in the values of λ1) observed 014213-3 MANUEL IÑARREA et al. PHYSICAL REVIEW E 106, 014213 (2022) TABLE I. Values of the maximal Lyapunov exponent λ1obtained for five chains of N=100,150,200,300,and 400 dipoles, and for the eight considered excess energies K. For each K, the value of λ1has been obtained by averaging the corresponding  λ(t) for time beyond t=5×108. The numbers in parentheses are powers of 10. N=100 N=150 N=200 N=300 N=400 K=0.05s8.883 (−6) 6.150 (−6) 3.962 (−6) 2.336 (−6) 1.485 (−6) K=0.1s2.849 (−5) 1.798 (−5) 1.201 (−5) 6.984 (−6) 4.777 (−6) K=0.25s1.371 (−4) 8.002 (−5) 5.482 (−5) 3.131 (−5) 2.112 (−5) K=0.5s4.247 (−4) 2.457 (−4) 1.650 (−5) 9.407 (−5) 6.330 (−5) K=0.75s7.331 (−4) 4.067 (−4) 2.678 (−4) 1.495 (−4) 9.915 (−5) K=0.9s9.950 (−4) 5.303 (−4) 3.302 (−4) 1.765 (−4) 1.201 (−4) K=1.1s1.572 (−3) 8.761 (−4) 5.717 (−4) 2.7617 (−4) 1.732 (−4) K=1.25s1.928 (−3) 1.017 (−3) 6.260 (−4) 2.966 (−4) 1.855 (−4) in of Fig. 2is the manifestation of an increasingly chaotic dynamics for increasing values of K. After the system is excited with small and medium excess energies (i.e., for K=0.05s,0.1s,0.25s,0.5s), the corresponding decay pattern of  λ(t)(seeFig.2) closely follows the well-known power law  λ(t)∼t−1of regular orbits. Then, at a given time, λ(t) deviates from the regular behavior, and it tends to converge to a nonzero value which is the corresponding maximal Lyapunov exponent λ1. This behavior reveals the chaotic nature of the excitations even for very small values of K, and has been already found in different kinds of lattices such as the FPUT problem [13,31,39], disordered lattices [17,19], or in the Bose-Hubbard model [14]. In all these systems, including our dipole chain, a possible explanation of the decay pattern of λ(t) could be the existence of areas close to regular regions in phase space where, after the initial excitation, the trajectory remains trapped, possibly for a long but finite time (given by τ=1/λ1), before entering the chaotic component of the phase space. As it was pointed out in Refs. [13,39], this behavior is theoretically sustained in the KAM and Nekhoroshev theorems [6–10]. For K=0.75s, the regular decay of λ(t) reduces to a very short time after the excitation, so that for larger values of the excess energy, no trace of regular behavior can be found in the time evolution of  λ(t). From the above discussion, whenever  λ(t) has reached the plateau, we can estimate the averaged maximal time Lyapunov exponent λ1by simply averaging in time the values of  λ(t). Thus, for five chains of N=100,150,200,300,and 400 dipoles, and for the eight considered excess energies K, we average  λ(t) for time beyond t=5×108. The obtained values of λ1are given in Table I. Furthermore, the linear behavior of λ1as a function of the energy density =K/N shown on a double log scale in Fig. 3suggests a power law λ1∼a.(8) The least-squares fit, see Fig. 3, revels a weak dependence of the exponent aon N.ForlargeN,ais expected to converge to an asymptotic value [31], which is not yet obtained for the N=400 chain analyzed here. It is important to notice that, from the values of λ1in Table I, there is a fast decrease of τ=1/λ1for increasing K, which indicates that, for low excitations, the system has difficulties to find the gateway to escape from the sticky quasiregular phase space regions to the nonregular counterpart. In addition, an increase in the excess energy would also lead to a global faster dynamics. B. Energy equipartition and thermalization of the dipole chain Regarding the energy equipartition attained by the system, we ilustrate in Fig. 4the time evolution of the averaged participation ratio  (t) for the same eight excess energies Kand for N=100, N=200,and N=400. We observe that, for small excess energies (K0.25s), there is a short transient (t400) during which a fast spreading of the excitation takes place. After that transient, we always find that  (t) fluctuates around a constant value. As we can see in Fig. 4, these asymptotic values are rapidly reached for t4×102, and the amplitudes of the fluctuations around them decrease with increasing time. Comparing the time evolutions of  (t) and  λ(t), depicted in Figs. 4and 2, we can see that, for small excess energies, the participation ratio  (t) takes its asymptotic values and the system reaches thermalization when it still exhibits regular behavior because  λ(t) has not yet deviated from the power law  λ(t)∼t−1.ForN=100,150,200,300, and 400, the asymptotic values, calculated as the average N=100 a=1.658 N=150 a=1.577 N=200 a=1.566 N=300 a=1.492 0.001 0.005 0.010 0.050 0.100 5 × 10 -6 10 -5 5 × 10 -5 10 -4 5 × 10 -4 10 -3 N=400 a=1.479 FIG. 3. Maximal Lyapunov exponent λ1as a function of the energy density =K/Nfor dipole chains with N=100, 150, 200, 300, and 400. Note that a double logarithmic scale is used. The values of aare the exponents that, for each N, are obtained from the fittings of the pairs (,λ1) to the power law of Eq. (8). 014213-4 CHAOS AND THERMALIZATION IN A CLASSICAL CHAIN … PHYSICAL REVIEW E 106, 014213 (2022) 12 45678 390 0.125 0.25 0.375 0.5 0.625 0.75 12 45678 39 (a) N=100 (c) N=400 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0012 45678 39 (b) N=200 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 FIG. 4. Time averaged participation ratio  (t) for eight excess energies. The horizontal dashed lines mark the equilibrium values  assuming that it follows a Boltzmann distribution for the local energies. Chains with (a) N=100, (b) N=200,and (c) N=400 dipoles are considered. for time beyond t=106of the values of  (t)forK= 0.05s,0.1s, and 0.25s,are ≈0.5923, 0.5915, 0.5903, 0.5890, and 0.5883, respectively, which indicates that independent of Nand for small values of K, the system always reaches approximately the same degree of equipartition. For larger excess energies, the fast initial transient in the averaged participation ratio described for small Kvalues is gradually replaced by a slower increase, such that  (t) eventually reaches asymptotic values slightly smaller than those obtained for small K. Note that, for the largest excitations, it takes times t108for  (t) to reach these constant values. The very long relaxation times for the largest values of Kindicate that, despite its chaotic dynamics and before  (t) reaches the asymptotic values, the chain exhibits a long-lasting nonergodic phase. Recent studies of the FPUT, Gross-Pitaevskii, and Klein-Gordon lattices [18,20,40,41] show that this nonergodic behavior is associated with the presence of robust breather excitations that prevent the system from reaching equipartitioning. For a chain of N=200 dipoles, the color maps of Figs. 5(a) and 5(b) show the time evolution of the local energies Ek(t) [see Eq. (6)] of two excitations with K=1.1s, but with different initial conditions θ100 ≈0.6614 and θ100 ≈0.8267 and with the momenta p100 according to Eq. (5). For the excitation depicted in Fig. 5(a), we observe for an early time interval t104the presence of two chaotic breathers, such that most of the energy of the system is strongly localized in a few energy carriers (dipoles) that follow complex trajectories. For t104, the chaotic breathers collapse into a single excitation mainly localized at the central dipoles. Finally, for t3×106, this central breather disappears and the energy appears to be broadly distributed among all the dipoles. In other words, in this case the energy transfer in the lattice is to a large extent determined by the presence of breathers [12] that, for a fairly long time period, keep the system far from equipartition, and exhibiting a persistent nonergodic dynamics. However, Fig. 5(b) indicates that the energy transfer mechanism of the second excitation scheme is completely different. No breather formation is observed, and even on short timescales the energy is fairly distributed among all the dipoles. Furthermore, the different profiles of the time evolution of the local energy observed in Figs. 5(a) and 5(b) are reflected in the evolution of the corresponding participation ratio (t). Indeed, while (t) for the excitation scheme depicted in Fig. 5(b) reaches fluctuating values around a constant nonzero value when t105[see the blue line in Fig. 5(c)], we have to wait up to t3×106to find a similar behavior for the evolution of (t) for the excitation scheme depicted in Fig. 5(a) [see the red line in Fig. 5(c)]. Finally, the behavior shown in Fig. 5indicates that, besides the amount of the excess energy K, the energy transfer mechanism is highly dependent on how Kis supplied to the system, i.e., it depends on the initial conditions of the excited dipole. As a consequence, a statistical approach is necessary to obtain a general global picture of the energy transfer in the dipole chain. Thus, we find that breathers are local hot spots that destroy the global ergodic dynamics and therefore prevent the thermalization of the system. It is worth noticing that this nonergodic dynamics coexists together with the global chaotic behavior that follows from the nonzero values of the maximal Lyapunov exponent λ1. This fact ultimately implies the lack of sensitivity of λ1to detect the presence of breathers, and thence to predict thermalization [19]. In other words, although the statistical character of λ1indicates that the system exhibits a global chaotic behavior, we can not use it to assure ergodic dynamics. A similar behavior, named as weakly nonergodic dynamics, was found by Mithun etal. [20] in a Gross-Pitaevkii lattice. IV. THE DIPOLE CHAIN AS A CANONICAL ENSEMBLE: EQUILIBRIUM VALUE OF (t) For the considered chains of N=100,150,200,300,and 400 dipoles, the asymptotic values of the averaged participation ratio indicate a degree of thermalization far below the complete energy equipartition regime, for which the participation ratio takes its maximum value (t)=1 [see Eq. (2)]. At this point, we pose the question of whether the asymptotic values observed in Fig. 4indicate that the chains are in a fairly thermalized regime. Indeed, a numerical estimate of the equilibrium value of (t) can be obtained in the following way. Taking into account the participation ratio Eq. (2), the estimate of its equilibrium value can be determined using the mean values of the local energy Ekand the squared local energy E2 kat equilibrium. Assuming that the system has 014213-5 MANUEL IÑARREA et al. PHYSICAL REVIEW E 106, 014213 (2022) Site Site t t t (c) (a) (b) 0 0.05 0.1 0.15 0.2 0 0.05 0.1 0.15 0.2 0.7 0.6 0.5 0.4 0.3 0.2 0.1 FIG. 5. Color maps showing the time evolution of the local energies Ek(t)[seeEq.(6)] for a chain with N=200 and K=1.1s but with different initial conditions θ100 ≈0.6614 [map (a)] and θ100 ≈0.8267 [map (b)], and with the corresponding momenta p100 according to Eq. (5). Panel (c) shows the participation ratio (t) corresponding to the trajectories of the color maps (a) (red color) and (b) (blue color). The horizontal dashed line in panel (c) marks the equilibrium value assuming that it follows a Boltzmann distribution for the local energies. a large number of dipoles and that its dynamics is ergodic, the distribution of the local energies Ek(t)[seeEq.(6)] of the dipoles is governed by a Boltzmann distribution. Similar approaches have been carried out in different kinds of lattices as the FPUT [18,42], the sine-Gordon [43], the Klein-Gordon [18], and the Gross-Pitaevskii [44–46]. Then, the partition function Zis given by Z= exp[−H()/T]d, (9) where are the 2Nphase variables appearing in the Hamiltonian Hin Eq. (3), and Tis the temperature of the system at equilibrium. Thus, the mean values of the local energy Ek and of the squared local energy E2 kat equilibrium can be computed as Ei k=1 Z Ei k(θk−1,θ k,θ k+1,pk) ×exp[−H()/T]d, i=1,2.(10) Due to the separability between the kinetic and the potential terms in H, these mean values can be written as Ek=T 2+θUkexp[−V/T]dθ θexp[−V/T]dθ ,(11) E2 k=3T2 4+θTU k+U2 kexp[−V/T]dθ θexp[−V/T]dθ ,(12) where dθ=dθ1...dθN, and V=V(θ1, ..., θN) and Uk= Uk(θk−1,θ k,θ k+1) are defined in Eqs. (4) and (7), respectively. Assuming that K=N k=1Ek=NEkand N k=1E2 k= NE2 k, the participation ratio at the equilibrium is given by =RN −1 N−1,(13) where R=Ek2/E2 k. To obtain the expressions of Ekand E2 k,wehaveto calculate numerically the integrals in Eqs. (11) and (12)for different temperatures Tand for chains with different numbers Nof dipoles. We expect the values of Tare in a range of temperatures that corresponds to the excess energy Kadded to the system. Thus, a suitable upper limit of the temperature is estimated using the equipartition theorem in the following way. Taking into account the boundedness of the potential energy Eq. (4) of the system, we can assume that, for large enough values of K, to a certain approximation the system takes the excess energy increasing only its kinetic energy. As a result, in a perfect energy equipartition regime, the mean value of the kinetic energy of each dipole is K/N=T/2, which provides an approximate value for the temperature T. Hence, taking into account the different chains and values of Kconsidered in these study, the upper limit for Twould correspond to the case of a chain with N=100 and an excess energy K=1.25s, which yields a upper limit of T=0.2. For many-degrees of freedom systems as the dipole chain considered here, the numerical computation of the integrals in Eqs. (11) and (12) should be performed by using a Monte Carlo method (MCM) [47]. However, even for the smallest chain with N=100 dipoles considered in this work, the convergence of the MCM when applied to the integrals in Eqs. (11) and (12) is extremely slow, such that the algorithm does not provide accurate results. To circumvent the curse of dimensionality, we assume that, since the local potential Uk=Uk(θk−1,θ k,θ k+1) only involves the nearest neighbors of the dipole k, a first estimate of Eqs. (11) and (12) could be obtained by neglecting the 014213-6 CHAOS AND THERMALIZATION IN A CLASSICAL CHAIN … PHYSICAL REVIEW E 106, 014213 (2022) TABLE II. Values of the parameter R[see Eq. (13)] for temperatures Tranging between 0.1 and 0.2. Labels NN, SN, TN and FN indicate that the calculations have been extended to nearest, second, third, and fourth neighbors, respectively. T=0.1T=0.12 T=0.14 T=0.16 T=0.18 T=0.2 NN 0.643 0.643 0.643 0.643 0.644 0.644 SN 0.600 0.601 0.601 0.602 0.602 0.603 TN 0.595 0.595 0.596 0.596 0.597 0.597 FN 0.596 0.594 0.596 0.595 0.595 0.596 coupling to the dipoles with i>k+1 and i<k−1. Thus, the corresponding mean values read Ek≈T 2+kUkexp[−Vk/T]dθk−1dθkdθk+1 kexp[−V/T]dθk−1dθkdθk+1 ,(14) E2 k≈3T2 4+kTU k+U2 kexp[−Vk/T]dθk−1dθkdθk+1 kexp[−Vk/T]dθk−1dθkdθk+1 . (15) Therefore, we are actually assuming that the chain is uncoupled beyond nearest neighbors and, for a given temperature T, the three-dimensional integrals in Eqs. (14) and (15) can be easily computed. This approximation can be improved by progressively including the interactions with further neighbors, such that, at a given order converged results for Eqs. (11) and (12)are obtained. We show in Table II, for different values of T,the values of factor R,Eq.(13), obtained up to nearest, second, third, and fourth neighbors. These results in Table II have been computed using an adaptative MCM, and more specifically the VEGAS algorithm [47,48]. As we can observe in Table II, for all the values of Tand up to third neighbors, we obtain for R≈0.59. Therefore, we conclude that up to four neighbors, the equilibrium value of the participation ratio can be approximated by =RN−1 N−1≈0.595N−1 N−1.(16) For the considered chains with N=100, 150, 200, 300, and 400 dipoles, Eq. (16) gives the equilibrium values 0.591, 0.593, 0.593, 0.594, and 0.594, which are in very good agreement with the numerical asymptotic results of (t) obtained in Sec. III B. Besides these calculations, we have also computed the full expressions Eq. (11) and (12) for chains with up to N=20 dipoles with the same adaptative MCM. The results obtained are in very good agreement with those shown in Table II.As both results are rather close, one can conclude that, once (t) settles to the (still fluctuating) values observed in Fig. 4,the system has almost achieved thermal equilibrium. As we mentioned already, for larger excitations, the asymptotic values of  (t) are slightly smaller than those for small excitations, being the degree of thermalization therefore slightly smaller. However, it is clear that even in the case of very large excitations, the system is capable of reaching a degree of thermalization which is comparable to the one reached with much smaller excitations, although that requires much longer times. From the numerical results of the time evolution of  (t) depicted in Fig. 4, it is clear that the longest thermalization times are found for the largest excitations. For example from Fig. 4(c) we obtain that t≈106is a rough estimate of the thermalization time for K⩽0.75s, while for the largest values of K, the corresponding thermalization times increase up to values t107. These estimates of the thermalization time also indicate that, except for the smallest excitation value K=0.05s, the thermalization time is always much longer than the corresponding Lyapunov times τ=1/λ1that can be obtained from Table I. Similar results are found for different chain lengths. Indeed, our computations show that, except for small excitations, the system becomes chaotic before an acceptable thermalization is achieved. V. CONCLUSIONS We have explored the connection between chaos, thermalization and ergodicity in a linear chain of hundreds of interacting dipoles. Starting from the GS, the chains have been excited by supplying different excess energies Kto one of the dipoles. Our analysis tools are the finite-time Lyapunov exponent λ(t)Eq.(1) and the participation ratio (t)Eq.(2), which provide information about the chaoticity of the system and the localization of the energy. For each value of K,λ(t), and (t) are statistically determined by averaging over 20 different initial conditions compatible with initial conditions Eq. (5). It turns out that the averaged  λ(t) shows always the same behavior: Once the system is excited, there is a transient during which λ(t) decreases in time. After the transient, there is a crossover to a plateau, and the corresponding averaged maximal Lyapunov exponent λ1is reached asymptotically. However, the value of Kdictates the strongly varying times scales of the behavior of  λ(t): A larger excess energy K implies a shorter transient and a larger value of λ1.Thishierarchy indicates an increasingly chaotic dynamics for increasing values of K. When the system is excited with small and medium Kvalues the decay pattern of  λ(t) closely follows the expected power law  λ(t)∼t−1of regular orbits. Then, at a given time,  λ(t) deviates from this regular behavior, and it tends to converge to the corresponding λ1value. For the largest excitation energies considered here, there is no trace of regular behavior in the decay of λ(t) before the corresponding asymptotic value λ1is reached. For small excess energies, the averaged participation ratio  (t) shows a short transient with a fast spreading of the excitation. After that transient,  (t) fluctuates around a constant value which depends on N. For larger values of K,thefast initial transient observed for small Kvalues is replaced by a slow increase. Thence, for long times,  (t) eventually reaches an asymptotic value. For the largest values of K, we find that the extremely long relaxation times of  (t) in comparison with the values of the Lyapunov times τ=1/λ1are due to the presence of chaotic breathers that keep the system far from equipartition. Furthermore, we observe that, besides the value of K, the energy transfer mechanism is highly dependent on the initial conditions of the excited dipole. As a consequence, 014213-7 MANUEL IÑARREA et al. PHYSICAL REVIEW E 106, 014213 (2022) a statistical approach as the one carried out in this paper is necessary to obtain a proper description of the energy transfer mechanism in the dipole chain. The numerically calculated asymptotic values of  (t) indicate a degree of thermalization well below the energy equipartition. Assuming the ergodicity of the system at thermal equilibrium, we have determined the thermal equilibrium values of the participation ratio by means of the Boltzmann statistics. We find that the thermal equilibrium values are in good agreement with the asymptotic values attained by  (t). Since both values are rather close, we can assert that the asymptotic values of  (t) indicate that the system has almost achieved thermal equilibrium, which on the other side, is far from a perfect energy equipartition regime. Here we have addressed only a dynamical subspace of the system. A next step would be the study of a more general case of the same setup. A natural continuation is the extension of this work to more complex dipole systems, such as dimerized dipole chains and one-dimensional bilayers of dipoles (e.g., diamond and sawtooth arrays [49]). One exciting direction is the possibility of identifying or even designing flat bands (see, e.g., Refs. [50–52]) in such one-dimensional arrays of dipoles and to study their impact on the energy transfer mechanism of the system. ACKNOWLEDGMENTS M.I. and J.P.S. acknowledge financial support by the 560FQ Spanish Project No. MTM2017-88137-C2-2-P (MINECO). R.G.F. gratefully acknowledges financial support by the Spanish projects PID2020-113390GB-I00 (MICIN), PY20-00082 (Junta de Andalucía) and A-FQM-52-UGR20 (ERDF-University of Granada), and the Andalusian Research Group FQM-207. These work used the Beronia cluster (Universidad de La Rioja), which is supported by FEDERMINECO Grant No. UNLR-094E-2C-225. [1] E. Fermi, J. Pasta, and S. Ulam, Studies of Nonlinear Problem, Los Alamos Report LA-1940, 1955 (unpublished); later published in Collected Papers of Enrico Fermi,editedbyE. Segré (University of Chicago Press, Chicago, 1965), Vol. 2, p. 978. [2] G. P. Berman and F. M. Izrailev, Chaos 15, 015104 (2005). [3] S. Flach, M. V. Ivanchenko, and O. I. Kanakov, Phys.Rev.Lett. 95, 064102 (2005). [4] S. Flach, M. V. Ivanchenko, and O. I. Kanakov, Phys. Rev. E 73, 036618 (2006). [5] H. Christodoulidi, C. Efthymiopoulos, and T. Bountis, Phys. Rev. E 81, 016210 (2010). [6] V. I. Arnold, Uspekhi Mat. Nauk 18, 13 (1963). [7] A. N. Kolmogorov, Dokl. Akad. Nauk SSSR 98, 527 (1954). [8] J. Moser, Nachr. Akad. Wiss. Göttingen Math.-Phys. K1,1 (1962). [9] M. Tabor, Chaos and Integrability in Nonlinear Dynamics: An Introduction (Wiley, New York, NY, 1989). [10] N. N. Nekhoroshev, Funct. Anal. Appl. 5, 338 (1971). [11] C. E. Wayne, Commun. Math. Phys. 96, 311 (1984). [12] T. Cretegny, T. Dauxois, S. Ruffo, and A. Torcini, Physica D 121, 109 (1998). [13] L. Casetti, M. Cerruti-Sola, M. Pettini, and E. G. D. Cohen, Phys. Rev. E 55, 6566 (1997). [14] A. C. Cassidy, D. Mason, V. Dunjko, and M. Olshanii, Phys. Rev. Lett. 102, 025302 (2009). [15] J. D. Bodyfelt, T. V. Laptyeva, Ch. Skokos, D. O. Krimer, and S. Flach, Phys.Rev.E84, 016205 (2011). [16] M. V. Ivanchenko, T. V. Laptyeva, and S. Flach, Phys.Rev.Lett. 107, 240602 (2011). [17] Ch. Skokos, Ch. Skokos, I. Gkolias, and S. Flach, Phys. Rev. Lett. 111, 064101 (2013). [18] C. Danieli, D. K. Campbell, and S. Flach, Phys.Rev.E95, 060202(R) (2017). [19] O. Tieleman, Ch. Skokos, and A. Lazarides, Europhys. Lett. 105, 20001 (2014). [20] T. Mithun, Y. Kati, C. Danieli, and S. Flach, Phys. Rev. Lett. 120, 184101 (2018). [21] T. Mithun, C. Danieli, Y. Kati, and S. Flach, Phys. Rev. Lett. 122, 054102 (2019). [22] V. Achilleos, G. Theocharis, and Ch. Skokos, Phys.Rev.E97, 042220 (2018). [23] A. Ngapasare, G. Theocharis, O. Richoux, Ch. Skokos, and V. Achilleos, Phys.Rev.E99, 032211 (2019). [24] B. Senyange, B. M. Manda, and Ch. Skokos, Phys. Rev. E 98, 052229 (2018). [25] B. Many Manda, B. Senyange, and Ch. Skokos, Phys. Rev. E 101, 032206 (2020). [26] G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn, Meccanica 15, 9 (1980). [27] G. Benettin, L. Galgani, A. Giorgilli, and J.-M. Strelcyn, Meccanica 15, 21 (1980). [28] C. Skokos, Lect. Notes Phys. 790, 63 (2010). [29] S. Iubini, O. Boada, Y. Omar, and F. Piazza, New J. Phys. 17, 113030 (2015). [30] A. Zampetaki, J. P. Salas, and P. Schmelcher, Phys. Rev. E 98, 022202 (2018). [31] G. Benettin, S. Pasquali, and A. Ponno, J. Stat. Phys. 171, 521 (2018). [32] J. Laskar and P. Robutel, Celest. Mech. Dyn. Astron. 80,39 (2001). [33] Ch. Skokos and E. Gerlach, Phys.Rev.E82, 036704 (2010). [34] E. Gerlach, Ch. Skokos, Discr. Cont. Dyn. Sys. Supp. 2011, 475 (2011). [35] E. Gerlach, S. Eggl, and C. Skokos, Int. J. Bifurcat. Chaos 22, 1250216 (2012). [36] K. Ullmann, A. J. Lichtenberg, and G. Corso, Phys. Rev. E 61, 2471 (2000). [37] B. M. Manda, Nonlinear dynamics and chaos in multidimensional disordered Hamiltonian systems, Ph.D. Dissertation, University of Cape Town (2021), arXiv:2105.03685. [38] B. Senyange, Chaotic behaviour of disordered nonlinear lattices, Ph.D. Dissertation, University of Cape Town (2021), arXiv:2104.11696. [39] M. Pettini, L. Casetti, M. Cerruti-Sola, R. Franzosi, and E. G. D. Cohen, Chaos 15, 015106 (2005). 014213-8 CHAOS AND THERMALIZATION IN A CLASSICAL CHAIN … PHYSICAL REVIEW E 106, 014213 (2022) [40] S. Iubini, L. Chirondojan, G.-L. Oppo, A. Politi, and P. Politi, Phys.Rev.Lett.122, 084102 (2019). [41] S. Iubini and A. Politi, Chaos, Solitons Fractals 147, 110954 (2021). [42] V. V. Mirnov, A. J. Lichtenberg, and H. Guclu, Physica D 157, 251 (2001). [43] C. G. Goedde, A. J. Lichtenberg, and M. A. Lieberman, Physica D59, 200 (1992). [44] B. Rumpf, Phys. Rev. E 77, 036606 (2008). [45] B. Rumpf, Physica D 238, 2067 (2009). [46] S. Iubini, R. Franzosi, R. Livi, G.-L. Oppo, and A. Politi, New J. Phys. 15, 023032 (2013). [47] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, Cambridge, UK, 2007). [48] G. Peter Lepage, J. Comput. Phys. 27, 192 (1978). [49] A. Andreanov and M. V. Fistul, J. Phys. A: Math. Theor. 52, 105101 (2019). [50] L. Morales-Inostroza and R. A. Vicencio, Phys.Rev.A94, 043831 (2016). [51] D. Leykam, A. Andreanov, and S. Flach, Adv. Phys. X 3, 1473052 (2018). [52] D. Leykam and S. Flach, APL Photon. 3, 070901 (2018). 014213-9