Thermal lifetime of breathers

We are interested in a simple magnitude that approximately measures the system’s arrival to thermal equilibrium. A good candidate for this measure is the participation number defined as

where $e_{n}$ the energy of the $n^{th}$ oscillator, $E=\sum e_{n}$ is the total energy of the system, that is, e.g., the Hamiltonian in Eq. (1), and the summation is over all the oscillators, with periodic boundary conditions. The participation number $P$ takes values between $P=1$, when all the energy is concentrated in a single oscillator, and $P=N$, if it is evenly distributed among all the oscillators.

The system has a constant energy, and the statistical ensemble is, therefore, microcanonical. If the system is ergodic and large enough, we may expect to observe canonical distributions for local observables, characterized by some temperature $T$ and constant thermodynamic $\beta=1/k_{B}T$. Note that although the values of $k_{B}$ and $T$ depend on the specific units, $k_{B}T=1/\beta$ is well defined, as the statistical meaning of temperature is given by the virial theorem [29] and the average kinetic local energy at thermal equilibrium is $\langle 0.5p_{n}^{2}\rangle=0.5/\beta=0.5k_{B}T$.

Each particle can be seen as interacting with a thermal reservoir at constant temperature and the statistical ensemble of each particle is canonical. Therefore, the probability that the energy of the $n^{th}$ particle in thermal equilibrium is $e_{n}$ is given by [29]:

The partition function $Z$ can be obtained with the condition that the total probability of finding some energy within the canonical ensemble is the unity. Then:

The actual upper limit of the integral cannot be infinity but a number of the order of the total system’s energy $E\sim Nk_{B}T=N/\beta$, with $N$ the number of particles. However, the integral $\int_{E}^{\infty}\exp(-\beta e_{n})\,\textrm{d}e_{n}=\exp(-\beta E)\sim\exp(-N)$ becomes negligible as it is smaller than $10^{-10}$ for $N>21$. Thus, (4) and the derivations below should be justified good approximations for the systems with a large number of particles $N$.

The average energy of a particle becomes

as should be expected. This result is obtained with several approximations and the time averages of the local energies at thermal equilibrium differ slightly from it.

In the following, we will use the ergodic theorem [29, 30] and replace the averages within the statistical ensemble with the time averages at thermal equilibrium. This is valid only for infinite time interval averages and, of course, our simulations may be long but always finite. In addition, our approach implicitly assumes that any energy density distribution $\{e_{n}\}$ has the same probability to occur. This subtle assumption does not follow, in general, from the fundamental ergodicity property that each microstate (some point in phase space) has the same probability of occurring. Indeed, one can easily show that ergodicity results in equal probabilities of energy density distributions for small ratios of coupling to energy density $\kappa/h$ in (1) and (7). In general, however, we have to expect systematic differences and will check that these differences are smaller than the standard deviations of the temporal fluctuations.

Let us calculate the infinity time average of the inverse of the participation number $P$, that is:

The time average along any solution of the dynamical system is defined and given in the following form:

Through the ergodic theorem, we now identify $\overline{P^{-1}}=\langle P^{-1}\rangle$ and $\overline{e_{n}^{2}}=\langle e_{n}^{2}\rangle$. Then:

Numerically we can suppose and it is confirmed by the numerical simulations that $\overline{P^{-1}}\simeq(\overline{P})^{-1}$ and, therefore, $\langle P^{-1}\rangle\simeq\langle{P}\rangle^{-1}$, i.e.,

In the rest of the paper we will identify often the average over a long time with the infinitely time averages and the statistical ensemble averages.

Note that the deduction cannot be done directly for $\langle P\rangle$ from (2) as the sum of particle energies would appear in the denominator. We have used many approximations that will be ultimately confirmed by numerical simulations. It is interesting to note that the result (5) contradicts the first intuition that at thermal equilibrium $P$ would be close to $N$.

It works surprisingly well for the first three systems under study. It does not work so well for the JJN system, because as this system has no on-site potential, the approximation that we can identify the energy of an oscillator as relatively independent is not really valid. See Sec. 4.3 for details and Ref. [28] for another approach including different approximations.

We produce a random vector of momenta $p_{n}$ values of numbers between $(0,1)$ and subtract its mean $p_{n}=p_{n}-\langle\{p_{n}\}\rangle$ so that the system’s initial momentum is zero. We re-scale the momenta so that the mean kinetic energy is the desired mean thermal energy $h=E/N$, and set the initial coordinates at zero $u=[u_{1},\dots,u_{N}]=[0,\dots,0]$. After about 100$T_{0}$ time units, with $T_{0}=2\pi/\omega_{0}$, we consider the system thermalized. We let it evolve even for a longer time and compare the mean value of $P$ with the theoretical one $N/2$. For both the hard and soft potentials, $\langle P\rangle$ is generally about one or two particles above the theoretical position of $N/2$ but can also be below it for some realizations depending on the observation time window. This property does not depend significantly on the local energy $h$. Particular realizations can be seen in Fig. 1 for a lattice with $N=64$ particles and systems QH and QS.

We also obtain similar results for the FKLJ system, but for the JJN system, $\langle P\rangle$ is about 6 particles above $N/2$. However, in this case, the fluctuations of $P$ at thermal equilibrium surpass $N/2$ frequently. Being conscious that $N/2$ is, therefore, not such a good measure of thermal equilibrium for the JJN system, we still keep it as a useful measure of the proximity to equilibrium for comparison of the breather thermalization times.

We conclude that the approximate value of $P=N/2$ is an appropriate measure to indicate that the system has approached equilibrium. However, it does not guarantee it, as should not be expected of a single parameter.

There is a huge variability of the thermalization time $t_{\text{th}}$ even for the same initial value of the participation number $P$. The particular values of $u_{n}$ and $p_{n}$ are of paramount importance for $t_{\text{th}}$. However, globally, as should be expected, for a given breather energy $E_{\text{b}}$, the larger the initial value of $P$, the shorter the thermalization time $t_{\text{th}}$, where the smallest linear correlation between $P$ and $t_{\text{th}}$ is found for the JJN system. Examples of the four systems are shown in Fig. 3.

The system with the hard quartic potential is described by Eq. (1), with $\omega_{0}=1$, the frequency of isolated oscillators, and the nonlinearity parameter $s=1$. The positive coupling constant $\kappa$ may take different values. Exact breathers, their obtention method, and their properties are presented in B, particularly, their energies are of the order of a few tenths. We consider two values of the coupling parameter $\kappa=0.05$ and $\kappa=0.10$ and two values of the average initial local mean thermal energy $h=0.02$ and $h=0.04$. Breather energies are taken from $E_{\text{b}}=5h$ to $E_{\text{b}}=21h$. Smaller breather energies do not provide enough localization for breathers to form and larger values create exceptionally stable and long-lived breathers. Even during days of simulations on a supercomputer using one hundred processors, we do not achieve enough thermalization cases to provide good statistics. Generally speaking, we reproduce the method presented in the previous section $N_{r}=10^{4}$ times until thermalization is achieved. There is no time limit set for a simulation, but the supercomputer has a time limit of several days. This means that very long thermalization times (or infinite) are excluded from the statistics.

Statistical results for the quartic hard potential case (QH) can be seen in Fig. 4. The error bars are obtained by adding and subtracting the standard deviation of the mean, but they are very small and difficult to observe. In total, the results were obtained for seventeen different breather energy values. The exponential dependence of the thermalization time and averaged discrete breather lifetime can be easily seen and estimated, especially for larger breather energies, since the thermalization time $t_{\text{th}}$ axis is shown on a logarithmic scale. The obtained results can also be justified by the fact that larger energy breathers are more resilient toward background noise or thermal fluctuations. In addition, notice that the breather lifetime decreases as the coupling constant $\kappa$ increases, whereas the breather lifetime increases with the increase of the mean thermal energy value $h$.

The Hamiltonian in Eq. (1) with $s=-1$ becomes soft. The on-site potential has a potential barrier at $u_{n}=\pm 1$ that separates the potential well centered at $u_{n}=0$ from a non-physical negative infinite potential well. The quartic soft on-site potential can be considered a reasonable approximation for a given system only inside the central potential well. Then, we write the code in such a way that if during a given simulation any $u_{n}$ leaves the safe well interval, the simulation is discarded. In this way, we are selecting a more distributed localization energy and diminishing the thermalization time. Therefore, the thermalization time cannot have an approximate exponential dependence on $E_{\text{b}}/h$. For lower coupling $\kappa$, existing breathers localized at a single site have small energies, and as they are not created by an increase in $E_{\text{b}}$, the thermalization time $t_{\text{th}}$ increases linearly with $E_{\text{b}}/h$. For a larger $\kappa$, when breathers have more energy and can be created, we observe again the approximate exponential dependence of $t_{\text{th}}$ as shown in Fig. 5. Results in Figure 5 are shown only for a single mean thermal energy value while varying the coupling parameter $\kappa$ since for larger mean thermal energy, e.g., $h=0.04$, the increase in energy favors delocalization and particles are leaving the potential well towards a negative infinity well, with nonphysical results. Nevertheless, we are still observing a clear pattern that the average breather lifetime increases with the decreasing value of the coupling parameter $\kappa$ as already observed in the quartic hard potential case in Fig. 4.

In this section, we consider a Frenkel-Kontorova (FK) system [22], that is, with cosine on-site potential, and with the Lennard-Jones (LJ) interatomic interaction potential. This system provides a useful model for atoms in a crystal, taking into account the periodicity of the crystal. The on-site potential represents the interaction with other parts of the crystal, and the LJ potential provides a strong repulsion when two atoms approach each other and a potential well with a force that tends to zero when the atoms move apart, corresponding with the physical characteristics of interatomic interactions [23]. From a more technical point of view, the cosine potential provides a soft potential, without the nonphysical characteristics of the quartic soft potential of having an infinite negative well, and the need to continuously control that the coordinates do not penetrate into that well. It has also been used as a model for lattice excitations in silicates in 2D hexagonal lattices by Bajārs, Eilbeck, and Leimkhuler (BEL) [31, 32, 33], and it has been shown that it is extremely easy to generate both stationary and moving breathers in one or two dimensions. Also, polarobreathers, i.e., breathers coupled to a charge, propagate in this system extremely well [34].

To compare with the quartic soft potential, we need an appropriate scaling as commented in C. The Hamiltonian is given by:

As shown in C, we compare the FKLJ Hamiltonian with the QS, so as the linearized dynamical equations become identical. The result being that $U_{0}=\omega_{0}^{2}\sigma^{2}/(2\pi)^{2}$ and $\kappa=(\omega_{0}^{2}\sigma^{2}/72)V_{0}$. For the QS system, $\sigma=2$ and $\omega_{0}=1$, therefore, $U_{0}=1/\pi^{2}\simeq 0.103$ and $V_{0}=\kappa/18$.

As this potential is soft and at low amplitudes is fitted with the quartic soft potential, there are similar features to be expected. In particular, for low coupling $\kappa=0.05$, the thermalization time $t_{\text{th}}$ shows a linear dependence with the delivered breather energy $E_{\text{b}}$, indicating that breathers are not formed because their energy is too low. The energy is rapidly dissipated. For $\kappa=0.20$ and $\kappa=0.30$, we again obtain the exponential dependence between $t_{\text{th}}$ and $E_{\text{b}}$, indicating that breathers are formed and long-lived in the thermalized system. Differently from the QS case, we can increase the temperature or the local mean thermal energy $h$, and $E_{\text{b}}$, without the control of $u_{n}$ going outside the infinite potential well. These results are shown in Fig. 6. Interestingly, thermalization times decrease as the coupling constant $\kappa$ increases but decrease as well when the mean thermal energy $h$ is increased, which is opposite to the quartic hard potential case and results in Fig. 4.

There are other techniques to create breathers in both the QS and FKLJ systems consisting of adding and subtracting some momentum to a couple of neighboring particles, favoring the creations of $\pi$-breathers with energies above the phonon spectrum, but, in this paper, we limit our research only to a single-site excited breathers.

This system describes an array of Josephson Junctions (JJ), which is described by the Hamiltonian:

Corresponding to a dynamical equation $\dot{p}_{n}=\ddot{q}_{n}=-\partial H/\partial q_{n}$:

The variable $q_{n}$ is written differently, as it is an angle variable representing the phase of the superconducting order parameter [27, 28]. The kinetic energy $\frac{\displaystyle 1}{\displaystyle 2}p_{n}^{2}$ is the island Coulomb charging energy, and $\kappa$ is the Josephson coupling between two neighboring superconducting islands. It also represents the angle of a rotating coupled pendulum, which allows for an easier intuition of the phenomena. This system has no on-site potential, which makes it a very different system from the previous three. For example, the phonon band is not bounded from below and extends from $\omega=0$ to $\omega=2\sqrt{\kappa}$ as shown in D. When increasing a momentum $p_{n}$ for breather creation it is convenient to subtract an appropriate amount to all rotators to prevent a global rotation of the system, which is equivalent to set to zero the mean electric potential of the JJ network.