1. Introduction
Discrete breathers (DBs), also known as intrinsic localized modes, are spatially localized nonlinear vibrational modes in defect-free discrete systems. In the past decades, huge numbers of theoretical and experimental studies had been devoted to confirming the existence and stability of DBs. Most of the theoretical studies focused on the anharmonic chains at zero temperature, from which the existence of DBs has been proved exactly (see [1, 2] and reviews [3, 4]). Since then, in other periodic systems such as cantilever arrays [5], Josephson junction arrays [6, 7], electrical lattices [8, 9], mass-spring chains [10], arrays of coupled pendula [11], chains of magnetic pendulums [12], and granular crystals [13, 14], the existence of DBs has been proved experimentally. There are also many studies in real crystals [15], for example, in ionic NaI [16], covalent Si, Ge, and diamond [17, 18], α-uranium [19], ordered alloys [20] as well as proteins [21, 22].
Standing DBs are localized modes of system energy. However, sometimes DBs can move with a small velocity. These mobile excitations can thus interact with the boundaries. It is thus interesting to study how this kind of DBs affects the energy relaxation at a finite-temperature system when the boundaries play a role. By using the cooling method, i.e. imposing the absorbing conditions on the boundaries, Aubry and Tsiron [23] first observed the stretched-exponential lattice energy relaxation induced by standing DBs in one-dimensional (1D) lattices with on-site potentials, in contrast to the standard exponential relaxation law of the corresponding linear chains. Subsequently, the counterpart two-dimensional lattices have been considered and a similar conclusion has been drawn in [24]. Further detailed studies [25, 26] of the harmonic and anharmonic chains without on-site potentials showed that the stretched-exponential energy relaxation is not universal, instead a usual exponential relaxation law is observed in a short time, followed by a power-law decay over a long time. The mechanism of the change of the power-law relaxation is induced by standing DBs. As the interactions of standing DBs with the boundaries are very weak, the energy release is actually suppressed and it is thus difficult to observe equilibrium states on typical simulation time scales. As a result, the lattice is in a metastable state, called the nonequilibrium residual state, similar to the glassy patterns in disordered systems.
In addition to the above progress, it would be worth noting that the properties of DBs have been found to be related to other physical properties. For example, by considering some complex crystal structures, such as the next-nearest-neighbor (NNN) interactions [27, 28], the lattice with altering mass or interactions [29], the on-site potentials [30–32] and the long-rang interactions [33], DBs have been studied in the perspective of energy transport. In such cases, DBs are more evident and sometimes can be moveable. This causes these kinds of DBs to be effective scatters of both phonons and boundaries [34], hence they can reduce the thermal conductivity [32, 33]. We also note that even only considering the role of nonlinearity, DBs can already affect some other macroscopic properties [35], e.g. elastic constants [36], thermal expansion [37, 38], and heat capacity [36, 38, 39].
The purpose of this paper is to study slow energy relaxation in a more complicated situation including mobile DBs. To this end, we shall focus on the Fermi–Pasta–Ulam–Tsingou-β (FPUT-β) chain with NNN interactions. By properly adjusting the ratio of the NNN coupling to the nearest-neighbor (NN) coupling, one might be able to produce the mobile DBs [27, 28]. Therefore, by employing such models we are devoted to exploring the combined effects of NNN interactions and nonlinearity (which produces moveable DBs) on the slow energy relaxation. The rest of this article is organized as follows: In section 2 we describe the reference models; section 3 presents the cooling methods to study energy relaxation and the physical quantities of interest; section 4 provides the main results for different lattice models, from which the role of NNN interactions and nonlinearity will be demonstrated; Finally, section 5 draws our conclusion.
2. Models
A general 1D lattice with both NN and NNN interparticle interactions and without on-site potential can be represented by
$\begin{eqnarray}H=\sum _{i}^{N}\left[\displaystyle \frac{{p}_{i}^{2}}{2}+V({x}_{i+1}-{x}_{i})+\gamma V({x}_{i+2}-{x}_{i})\right],\end{eqnarray}$
where pi is the ith (totally N particles and all with unit mass) particle's momentum, xi is its displacement from the equilibrium position, V(ξ) is the interparticle potential and the parameter γ specifies the comparative strength of the NNN to NN couplings. We use the FPUT-β interparticle interaction $V(\xi )=\tfrac{1}{2}{\xi }^{2}+\tfrac{1}{4}{\xi }^{4}$ to mainly consider the combined role of NNN coupling and nonlinearity in energy relaxation. For such a model, it would be worth noting that it has a special phonon dispersion relation (see figure 1 for the phonon spectrum): ${\omega }_{q}=2\sqrt{{\sin }^{2}(q/2)+\gamma {\sin }^{2}q}$, where q is the wave number for phonons and ωq is the corresponding frequency. From this dispersion relation, one obtains the phonon's group velocity ${v}_{g}={\rm{d}}{\omega }_{q}/{\rm{d}}q={\omega }^{-1}[\sin q+2\gamma \sin (2q)]$ and finds that: (i) for 0 < γ < 0.25, vg > 1, which is in contrast to vg = 1 in the case of systems with NN coupling; (ii) for γ = 0.25, vg is very close to zero in a wider k domain near the Brillouin zone boundary; (iii) for γ > 0.25, both ωq and vg are actually degenerate, i.e. one ωq/vg corresponds to two q values. These unusual properties can favor the formation of a special highly moving localized excitation (mobile DBs) in the presence of appropriate nonlinearity [27, 28] and thus one can expect that such properties would greatly influence the energy relaxation.
Figure 1. (a) Phonon frequency ωq and (b) group velocity vg versus wave number q for 1D lattice without on-site potentials. From bottom to top, the curves correspond to γ = 0, 0.25, 0.4, and 1, respectively. |
3. Method
We use the cooling method to reveal the slow energy relaxation and to obtain the information of DBs. For this purpose, we first thermalize the lattices to a fixed temperature T = 0.1 by using the Nose–Hoover [40] heat baths. Then after the systems are fully thermalized for a considerably long time, we remove the heat baths and add the absorbing boundary conditions. That is to say, we impose damping to both ends (which are free) of the lattices, and so the equation of motion of the system can be expressed by
$\begin{eqnarray}\begin{array}{rcl}\ddot{{x}_{i}} & = & V^{\prime} ({x}_{i+1}-{x}_{i})+V^{\prime} ({x}_{i}-{x}_{i-1})+\gamma V^{\prime} ({x}_{i+2}-{x}_{i})\\ & & +\gamma V^{\prime} ({x}_{i}-{x}_{i-2})-\eta {p}_{i}({\delta }_{i,1}+{\delta }_{i,N}).\end{array}\end{eqnarray}$
Here $\ddot{{x}_{i}}=\tfrac{{\rm{d}}{p}_{i}(t)}{{\rm{d}}t};$ $V^{\prime} (\xi )$ represents $\tfrac{\partial V}{\partial \xi };$ η is the dissipation exponent and we fix η = 0.1 in practice; δ is the Kronecker delta.By performing the cooling method, we are interested in the decay of the total energy. To reveal this decay, we first define the symmetrized site energies:
$\begin{eqnarray}\begin{array}{rcl}{h}_{i} & = & \displaystyle \frac{1}{2}{p}_{i}^{2}+\displaystyle \frac{1}{2}\left[V({x}_{i+1}-{x}_{i})+V({x}_{i}-{x}_{i-1})\right.\\ & & \left.+\gamma V({x}_{i+2}-{x}_{i})+\gamma V({x}_{i}-{x}_{i-2})\right].\end{array}\end{eqnarray}$
With such a definition, the total energy of the system is given by $E={\sum }_{i=1}^{N}{h}_{i}$. We can now study the decay of the normalized energy E(t)/E(0), which represents the ratio of the decayed energy to the initial energy. As it has been shown that in some cases the exponential law or the stretched-exponential law is observed. To identify what exactly the law is, it is better to plot $\begin{eqnarray}D(t)=-\mathrm{ln}\left[E(t)/E(0)\right].\end{eqnarray}$
We remind readers that since the energy localization might be strongly inhibited by fixed-end boundary conditions [25, 26], in our simulations we apply the free-end boundary conditions. The motion equations are integrated with the velocity-Verlet algorithm [41] with a time step 0.01. With both heat baths and damping, the systems are first evolved 106 times with the heat baths to reach the stationary state under the given temperature, and next evolved additional 106 times without the heat baths but with the damping presented, which enables us to reveal the relaxation process induced by the absorbing boundaries.
4. Results
4.1. Harmonic system
To show the role of NNN coupling in slow energy relaxation separately, let us first check the results of the relevant linear systems with NNN coupling. For the same harmonic system with the NN coupling only, it has been theoretically predicted that [25]
$\begin{eqnarray}\displaystyle \frac{E(t)}{E(0)}=\left\{\begin{array}{lll}{{\rm{e}}}^{-t/{\tau }_{0}} & {\rm{for}} & t\ll {\tau }_{0},\\ {\left[2\pi (t/{\tau }_{0})\right]}^{-1/2} & {\rm{for}} & t\gg {\tau }_{0}.\end{array}\right.\end{eqnarray}$
This means that one generally finds two energy relaxation processes, i.e. the exponential law and the inverse-square-root law, for short and long timescales, respectively. Such two laws have been verified in a harmonic chain with NN coupling only under a small N = 32, but some deviations have also been numerically observed [25, 26], i.e. a new exponential relaxation recovers in a long time. In view of this, we study the dependence of the energy relaxation properties on γ here in a large system size N = 4096 in figure 2(a). Surprisingly, in a time of t < 104, regardless of γ, the results of D versus t follow two general power-laws (∼tν) with ν = 0.68 and ν = 1.03, respectively. This suggests that in addition to the previously theoretically predicted exponential law, there is another stretched-exponential law ${{\rm{e}}}^{-{t}^{\nu }/{\tau }_{0}}$ with ν = 0.68 at a shorter time, independent of the NNN coupling ratio γ. Figure 2(a) also suggests that the long-time behavior of the energy relaxation should depend on γ. To see this more explicitly, figure 2(b) depicts the results of E(t)/E(0) versus t in harmonic systems with NNN coupling. A log–log plot then helps us clearly identify several power-law exponents. The best fittings give E(t)/E(0) ∼t−μ with μ = 0.54, μ = 0.52, μ = 0.16, and μ = 0.13 for γ = 0, γ = 0.025, γ = 0.25, and γ = 0.4, respectively (in a long time). It indicates that including the NNN coupling seems to further slow down the energy relaxation process. This is not trivial since from figure 1 by including the NNN coupling, on one hand the phonon group velocity vg is increased with the increase of γ. This speeds up the energy relaxation process. On the other hand, the increase of γ causes more phonons with higher frequencies to emerge, which however slows down the relaxation. Therefore, in a linear system with NNN coupling but without nonlinearity and DBs, the overall effect is to depress the energy relaxation process during the cooling.
Figure 2. (a) and (c): D versus t with several values of γ, where the two dashed lines denote the best fittings of tν and they indicate the previously predicted exponential law and the new observed stretched-exponential law decay of E(t)/E(0), respectively. (b) and (d): E(t)/E(0) versus t under several γ, where the dashed lines denote the best fittings of t−μ that indicate the power-law decay of E(t)/E(0). Figures (a) and (b) are for the harmonic potential; figures (c) and (d) are for the FPUT-β potential. |
4.2. FPUT-β system
With the above understanding, let us now consider the combined effects of NNN coupling and nonlinearity. For this purpose, as mentioned we focus on the popular 1D FPUT-β lattice. Figure 1(c) shows that the first stretched-exponential and exponential laws of E(t)/E(0) for t < 104, which are similar to that shown in figure 1(a). Therefore, these two relaxation laws seem independent of the NNN coupling, even when one includes the nonlinearity in the systems. This may also suggest that such two laws for t < 104 are general, at least for both the linear and nonlinear systems with NNN coupling. Therefore, the origin of this finding is worth studying in the future.
Again, the distinctions for different γ in the FPUT-β systems still lie in the long-time behavior. As can be seen in figure 2(d), E(t)/E(0) follows a t−μ law at long timescales, similar to those observed in the harmonic systems as shown in figure 2(b). A further comparison of figures 2(b) and (d) shows that the power-law exponent μ for harmonic and FPUT-β systems under the same γ is different, i.e. μ for the FPUT-β lattice is larger than that for harmonic lattice as long as γ ≠ 0. It presents a clue that the combined effects of NNN coupling and nonlinearity seem to speed up the relaxation during the cooling.
In order to have a close look at this trend, in figure 3 we plot μ against γ for both the harmonic and FPUT-β lattices. Therein four data points are extracted from figures 2(b) and (d), while others are calculated additionally in the same way. Indeed, in the range of γ investigated, all μ of FPUT-β systems are larger than those of the harmonic systems. Indeed, this indicates that the combined effects of the NNN coupling and nonlinearity are to make the energy relaxation relatively faster if compared with the counterpart linear systems.
Figure 3. μ versus γ, where the hollow and solid circles represent the results of harmonic and FPUT-β systems, respectively; the dashed line indicates the results of γ = 0.25. For a better illustration, we also plot D versus t for γ = 0.25, separately in the inset. |
4.3. The origin of the relatively faster energy relaxation: the role of mobile DBs
We lastly explain the observed relatively faster energy relaxation in a nonlinear system compared to the counterpart linear systems. As it has already been pointed out in the introduction, the slow energy relaxation behaviors in nonlinear lattices are usually related to the energy localization and these localization origins from the excitation of standing DBs. In the absence of nonlinearity, it was previously regarded that the energy relaxation will turn back to the traditional exponential decay [25]. However, our above results in figure 2 suggest that even for linear systems, this is not the case. The long-time behavior of the energy relaxation seems always to be a power-law decay. Combining the results of linear and nonlinear systems (see figures 2(a)–(d)), one might conclude that the power-law energy relaxation in a long time is general but the details of the power-law exponent are then affected by the localization induced by DBs. To illustrate this point, in figure 4 we plot a snapshot of the lattice displacements at the end of the simulation, i.e. by imposing the damping to both harmonic and FPUT-β systems after a time 106 (with γ = 0.25 for an example). As can be seen, only in the FPUT-β system do we identify the emergence of DBs (see the inset of figure 4(b), a profile of DBs can be found there), which is in clear contrast to that observed in the harmonic system (see the inset of figure 4(a)).
Figure 4. Snapshot of atomic displacement after imposing the damping to the free-end boundaries by a time t = 106, where the left (right) panel corresponds to a harmonic (FPUT-β) system with NNN coupling (γ = 0.25). The inset in each panel is a zoom for the reference sites. |
Since DBs can localize energies, one may wonder why the energy relaxation in the nonlinear systems with NNN coupling is relatively faster than that in the counterpart linear systems. This can be understood from figure 5 if we further examine the mobility of DBs. In figure 5 we plot the same snapshot of the lattice displacements after imposing the damping of the free-end boundaries of the nonlinear lattice (γ = 0.25) for four long times. It can be seen that due to the mobility of DBs, they can interact with both phonons and boundaries, which produces an additional passageway for speeding up the energy relaxation. This is surely different from the main passageway of phonons-boundaries scattering as shown in the relevant linear model (see figure 4(a)).
Figure 5. Snapshot of atomic displacement after imposing the damping to the free-end boundaries by a time t = 1.1 × 106 (a), t = 1.2 × 106 (b), t = 1.3 × 106 (c), and t = 1.4 × 106 (d) in the FPUT-β system with NNN coupling (γ = 0.25). |
5. Conclusion
To summarize, by applying a damping term into the system's boundaries, we have numerically studied the energy relaxation processes in both 1D linear and FPUT-β lattices when the NNN interactions are considered. As we expect that including the NNN coupling might need a longer time for systems to relax, we have considered a system size much larger than that that was taken into account of in previous studies. For both linear and nonlinear lattices, generally the short-time (t < 104) relaxation behaviors are given by a stretched-exponential law, followed by an exponential law, which is independent of the ratio of the NNN coupling. This justifies that even linear systems can also support the stretched-exponential relaxation law at some timescales. Furthermore, the long-time (t > 104) behavior is however determined by a power-law with the exponent strongly dependent on γ, i.e. as γ increases, the relaxation is slowed down. This means that including the NNN coupling makes the long-time power-law energy relaxation slower. Further introducing the nonlinearity, however, in turn, speed up the relaxation, which is induced by the emergence of mobile DBs. These mobile DBs can interact with both phonons and boundaries, which then produce an additional passageway for scattering the energies. Finally, this novel energy relaxation behavior, induced by mobile DBs in the nonlinear lattices with NNN coupling, is surely different from the stretched-exponential relaxation observed in lattices with on-site potentials that are mainly produced by the standing DBs [23, 24].