Frequency-momentum representation of moving breathers in a two dimensional hexagonal lattice

The study of nonlinear excitations as breathers in a hexagonal lattice is motivated by the experimental evidence of the existence of nonlinear excitations, also known as intrinsic localized modes (ILM), that transport energy and momentum in a quasi-one-dimensional form along the chains of the cation layer of muscovite mica and other layered silicates [1, 2, 3, 4]. These lattice excitations were called quodons making reference to the quasi-one-dimensional propagation and the particle-like behaviour identified in the fossil dark tracks in muscovite [5, 6]. It is worth noting that muscovite has been demonstrated able to record the passage of swift particles as positrons, protons and antimuons [7, 8]. The reader is referred to three complete reviews [9, 10, 11].

There have been several attempts to model these nonlinear excitations in one dimension as kinks or crowdions [12, 13, 14, 15] and breathers [16]. These models use a substrate potential and interaction potential that have been deduced from physical properties and empirical potentials, and therefore are able to provide physical values of the magnitudes involved as, for example, energies, $\simeq 27$ eV for crowdions and $\simeq$ 0.3 eV for breathers. The latter article  [16] also develops the theory of exact breather solutions, which are generally coupled to an extended plane wave called wing and are therefore called pterobreathers, “ptero” meaning “wing” in Greek. However, for some specific velocities the wings disappear and the exact pterobreathers become exact breathers. Their stability can be determined by a variant of the Floquet method. The same article also characterized exact pterobreathers and breathers by their unique frequency in the moving frame. These models have the obvious limitation of being one-dimensional but this limitation also allows for easier mathematical and numerical insight into the physical and mathematical properties.

Other model, inspired by the same phenomena, was proposed [17, 18] using a more generic model in a two dimensional (2D) hexagonal lattice also with a substrate potential. They obtained that kinks or crowdions can propagate in a quasi-dimensional form [19] and also breathers [17, 18], depending on the specifics of the potential. Breathers are also scattered by other breathers and can migrate to other close-packed chains [20].

In this paper we extend the theory developed in [16] to two dimensions and apply it to the model in Ref. [18] in order to describe moving breathers in the moving frame in a hexagonal lattice with a substrate potential, to obtain their frequency-momentum representation, their frequencies in the moving frame and the existence or absence of wings. This is fundamental to be able to interpret possible signatures of localized nonlinear waves in physical spectra of real crystals [21].

The paper is written illustrating the theory together with particular analytical and numerical results. The mathematical model is briefly described in Sect. 2. The on-site potential used to describe the direct hexagonal lattice and the phonon bands to describe the reciprocal and momentum lattices are presented in Sect. 3. The phonon frequencies and polarization are obtained in Sect. 4. The theory of exact traveling waves in 1D and 2D is developed in Sect. 5. In Sect. 6 quasi-exact solutions are used to obtain important characteristic parameter values for the exact propagating modes. The procedure to obtain exact traveling wave solutions is explained in Sec. 7 and the corresponding exact solutions are presented in Sect. 8. The paper is concluded with the Conclusions and an appendix containing some details of the reciprocal basis.

We consider a system of $N$ particles with mass $m=1$ with their position described by coordinates $(x,y)$ in a Cartesian system of reference and $(u_{x},u_{y})$ with respect to their equilibrium positions. Their kinetic energy is given by $K=\tfrac{1}{2}(\dot{x}^{2}+\dot{y}^{2})=\tfrac{1}{2}(\dot{u}_{x}^{2}+\dot{u}_{y}^{2})$, the dot indicating the derivative with respect to time.

The particles are in their equilibrium position in a 2D hexagonal lattice corresponding to the minima of the on-site potential given in the dimensionless form by:

where the unit distance has been chosen $a=1$ and $0\leq U(x,y)\leq 1$.

The interaction between particles is given by the Lennard-Jones potential:

being $\epsilon$ a measure of the relative strength of the interaction potential with respect to the on-site potential.

Without loss of generality we consider $\epsilon=0.05$. For short and long-range interactions we consider smooth cut-off of the Lennard-Jones potential in Eq. (2) with a cut-off radius $r_{c}$ described in detail in [18]. In this paper we consider $r_{c}=3a$. Derived Hamiltonian equations of dynamics are integrated in time with the second order time reversible symplectic Verlet method. In the following, all numerical examples are performed with time step $\tau=0.04$ and periodic boundary conditions. The dimensions of the lattice are $N_{1}=64$ and $N_{2}=32$, i.e., $N=N_{1}N_{2}$, unless stated otherwise.

To produce quasi-exact discrete breather solutions, see Sec. 6, we excite three neighboring particle velocities with the pattern:

where $\gamma>0$. As noted in [20] larger values of $\gamma$ produce faster moving quasi-exact breathers with larger particle displacements in the direction of propagation. See also Table 1.

For understanding the frequency-momentum representation it is convenient to review the properties of the hexagonal lattice. In this section we perform that revision while at the same time we present important properties of the system, such as its potential energy and linear spectrum.

In Ref. [18] the equations that determine the eigenvectors and eigenvalues were deduced. Here we reproduce the deduction but also obtain the polarization of the phonons. The phonon dispersion relation is plotted along special paths in Fig. 4 and the phonon surfaces are displayed in a primitive reciprocal cell in Fig. 5. The polarization is described in Fig. 4 and displayed in Fig. 6.

Linear modes or phonons are given by solutions (to the linearized equations) of the form:

The variables are as follows: $\mbox{\boldmath${w}$\unboldmath}=[u_{x},u_{y}]^{T}$ is the vector of displacements in the $x$ and $y$ directions of a particle of the lattice and $\mbox{\boldmath${A}$\unboldmath}=[a,b]^{T}$ is its polarization vector. Without loss of generality we can choose $a$ and $b$ such that $\sqrt{a^{2}+b^{2}}=1$ and $a\geq 0$. In principle $a$ and $b$ could be complex, but in fact, in our system, they can always be chosen real. $\mathbf{k}$ is a wavevector and $\mathbf{R}=n_{1}\mathbf{a}_{1}+n_{2}\mathbf{a}_{2}$ is the equilibrium position corresponding to a given particle. The angular frequency $\omega$ can be chosen positive without loss of generality because the propagation direction of the wave is given by $\mathbf{k}$.

In terms of these indexes the equation above becomes:

leading to the following equation [18]:

where $\kappa=16\pi^{2}/9=17.5460$ is the minimal squared frequency, derived from the linearization of the on-site potential in Eq. (1), and $\beta=V_{LJ}^{{}^{\prime\prime}}(1)=72\epsilon=3.6$.

We obtain the homogeneous system of linear equations:

where

Note that $A\geq\kappa$ and $B\geq\kappa$ and therefore positive, while the sign of $D$ depends on $k_{x}$ and $k_{y}$.

The necessary condition for the existence of nonzero solutions of the system in Eq. (13) is the characteristic equation: $(\omega^{2}-A)\times(\omega^{2}-B)-D^{2}=0$ or $(\omega^{2})^{2}-(A+B)\omega^{2}+(AB-D^{2})=0$, that is:

As $\omega>0$, there are two phonon bands corresponding to the two signs of $\pm$. Let us suppose initially that $D\neq 0$, which also implies that $a\neq 0$. For the frequencies $\omega_{1,2}$, the two equations in Eq. (13) are equivalent and only one is needed, for example, the first one. Let us denote $\alpha=b/a$, then

The normalized polarization vector becomes $a_{1,2}=1/\sqrt{1+\alpha_{1,2}^{2}}$ and $b_{1,2}=\alpha_{1,2}/\sqrt{1+\alpha_{1,2}^{2}}$.

It turns out that $D=0$ is an interesting case. In this case, there are two possibilities, either $\sin(k_{x}/2)=0$ or $\sin(\sqrt{3}/2k_{y})=0$. In the first one $\sin(k_{x}/2)=0$, which implies that $k_{x}=0$ as $k_{x}=2\pi m$ is equivalent, then:

These functions reproduce the curves in Fig. 4-right and also identify that the upper curve corresponds to $a=0$, that is, there is no vibration in the $u_{x}$ coordinate, while the lower curve corresponds to $b=0$, that is, there is no vibration in the $u_{y}$ coordinate.

For the second possibility, $\sin\left(\tfrac{\sqrt{3}}{2}k_{y}\right)=0$, then $\tfrac{\sqrt{3}}{2}k_{y}=m\pi$ with $m\in\{0,1\}$ and:

where the $+$ or $-$ in $\pm$ corresponds to $m=0$ and $m=\pi$, respectively.

The curves corresponding to $(\pm)=+1$ correspond to the curves in Fig. 4-left. The upper one have polarization $a=1$ for $q_{1}<2/3$ at point $K$ and $a=0$ for $q_{1}>2/3$. The lower curves have the opposite polarization. The curves with $(\pm)=-1$ can also be observed in the phonons that are in the phonon band depicted in Fig. 9 right. The phonon surfaces can be seen in Fig. 5.

The breathers described in this work derive from the first maximum at point K in the 1st Brillouin zone and, therefore, the central mode of the components of a breather propagating in the $k_{x}$ direction is a mode with polarization $b=u_{y}=0$, that is, it has no perturbation in the $y$ direction. This is a very favorable circumstance because, in this way, it tends not to perturb the adjacent chains. Note also that for $D\neq 0$, $D$ and therefore $\alpha$ changes sign with the sign of $k_{y}$, therefore, the sum of plane waves with equal amplitude, the same $k_{x}$ and opposite $k_{y}$ and $b$ will cancel out leading to polarization $u_{y}=0$ for the breather, that is, with vibration only in the $x$ directions. This is what happens with the breathers observed in our system.

In this section we review the theory of moving nonlinear excitations developed in Ref. [16] and extend it to two dimensions. It is presented in a heuristic way together with results for traveling waves in our system. The ideas developed here use the fact that a periodic function in the supercell described in Sect. 3 can be expressed as a sum of plane waves with wavevectors in the corresponding momentum space. The amplitudes of these plane waves are obtained by the Fourier Transform (FT). We will specify often the variables that are considered, as, for example, XYTFT indicates the Fourier transform on the two spatial variables and time.

We explore different simulations for good traveling solutions that can be the seed to obtain exact traveling solution as explained in the next section. In general, the procedure is to find a long-lived one, add dissipation at the borders of the simulation cell parallel to the close-packed line where the initial perturbation has been provided and thereafter let it propagate some time without dissipation. In our study we initiate traveling solutions with the velocity pattern in Eq. (3) and different values of $\gamma$. Such excitation produces only small amount of phonon background.

If we observe the particles in the central line, we can label them with just an index $n\in\{0:N_{1}-1\}$, i.e., $u_{n}(t_{l})$, with $t_{l}=l/N_{t}$, and $l\in\{0:N_{t}-1\}$. From the breather line in the XTFT we obtain the breather velocity $V_{b}=\Delta\omega/\Delta k$. An exact solution has the property that $V_{b}=s/T_{F}$, with $s$, the step, and $T_{F}$, the fundamental time. Then $T_{F}=s/V_{b}$ is a multiple of the velocity period $1/V_{b}=1/f_{b}$. We consider times $T_{s}=s/V_{b}$ with $s=1,2,..$ which are candidates to be $T_{F}$. We find a time $t_{0}$ for which the breather amplitude is at maximum and compare the functions $u_{n}(t_{0})$ and $u_{n+s}(t_{0}+t_{s})$, being $t_{s}=\mathrm{round}(T_{s}/\tau)\tau$, where $\mathrm{round}(x)$ gives the closest integer to $x$ as $t_{s}$ is the closest integration time to $T_{s}$. Usual steps are $s=1$ or $s=2$ [16] and we check the overlapping of the functions. Of course, it is possible to automatize the process, but in practice it is not worth it. Figure 12 shows a good example. Velocity was observed to be $V_{b}=0.3147$, therefore, candidates for fundamental time are $T_{1}=1/V_{b}=3.1776$ and $T_{2}=2/V_{b}=6.1369$, with integration step $\tau=0.04$, $t_{1}=\mathrm{round}(T_{1}/\tau)=76\tau=2.9600=3.1600$ and $t_{2}=\mathrm{round}(T_{2}/\tau)=152\tau=2.9600=6.3600$. Mathematically, $s=2$ results in a good agreement, but $s=1$ is good enough and if $s=2$ and $T_{F}=T_{2}$, then at midtime $T_{F}$ the coordinates would be inverted, which does not happen. Therefore, $s=1$ and $T_{F}$ is close to $T_{1}$. The other adjustment for $T_{F}$ is that the breather line for $k=0$ is the breather frequency in the moving frame is multiple of the fundamental frequency $\omega_{MF}=m_{b}\omega_{F}=m_{b}\,2\pi/T_{F}$. In this way a small adjustment of $T_{F}$ is possible. Note that the agreement is not perfect, first, due to the fact that the breather is not exact and, second, due the finite time sampling interval $\tau$.

Breathers obtained by providing some recoil velocity are interesting as this is a likely mechanism from $\beta$ decay, as in ⁴⁰K, or after the impact of incoming radiation. As they have long lives, they are quite good solutions, however, they are only an approximation to an exact solution. Exact solutions are interesting as they are amenable to mathematical and numerical methods and their lives are theoretically infinite. The generic solutions are breathers accompanied by an extended wing, that is, pterobreathers [16]. In the following, for completeness, we recall some concepts and the method used to obtain exact pterobreathers.

An exact traveling solution is a solution that repeats itself after some time $T_{F}$, the fundamental time, displaced a lattice step ${\bf s}=[s_{1},s_{2}]$, where $s_{1}$ and $s_{2}$ ar integers. Let ${\bf n}=[n_{1},n_{2}]$, where $n_{1}$ and $n_{2}$ are integers that represent a lattice site $n_{1}\mathbf{a}_{1}+n_{2}\mathbf{a}_{2}$, and ${\bf U_{n}}=\{{\bf u_{n}}(t),{\bf v_{n}}(t)\}$ represents the positions and velocities of the particles of the system at site $\bf n$.

Let us define the maps $\hat{L}[{\bf m}]$, with ${\bf m}=[m_{1},m_{2}]$, and $\hat{T}[T]$ as the operators of translation in the lattice and time, respectively:

The action of the operator $\hat{T}$ is obtained by integrating the equations of dynamics of the system.

The composite translation map in the lattice and time $\hat{S}[{\bf m},T]$ becomes:

Therefore an exact traveling solution with fundamental time $T_{F}$ and step $\bf s$ is given by the equation:

The election of $t=0$ is irrelevant as the Hamiltonian and dynamical equations are time invariant. Let as define the function $\bf f$:

The exact solution is obtained by the Newton method. Let as suppose that ${\bf U}^{0}$ is an approximate solution, a seed, that is, ${\bf f}({\bf U}^{0})\neq 0$ but small. The approximate ${\bf U}^{0}$ can be identified by its initial positions and velocities or by other means as the Fourier coefficients, but this does not affect what follows.

If ${\bf U}={\bf U}^{0}+\delta{\bf U}$ is an exact solution close to ${\bf U}^{0}$, i.e., ${\bf f}({\bf U})={\bf f}({\bf U}+\delta{\bf U})=0$:

where $\nabla{\bf f}$ is the Jacobian of ${\bf f}$ calculated at the approximate solution ${\bf U}^{0}$.

The new initial variables become ${\bf U}={\bf U}^{0}+\delta{\bf U}$ and they will be closer to an exact traveling solution and becomes a new seed. Repetition of the operation above until the desired accuracy is achieved leads to an exact solution.

Note that it is fundamental to have a good seed and an accurate guess of the step and fundamental time. We are interested, in principle, in traveling wave solutions along close-packed lines, which are most likely to occur in a mathematical and physical system. As all close-packed lines are equivalent, the simplest one is a line parallel to $\mathbf{a}_{1}$ and the step becomes $\mathbf{s}=[s_{1},0]$. For simplicity, we will refer often to the step just as the scalar $s\equiv s_{1}$. Typical values of $s$ are very low integer numbers 1,2,…

With the velocity pattern in Eq. (3) we find different quasi-exact traveling waves composed of a breather and a soliton. The soliton is always present and can be reconstructed by filtering out the frequencies above the soliton line and performing the inverse Fourier transform. Figure 7 shows both profiles of the soliton-breather and the soliton. The soliton has a “S”-shape, a compression followed by a decompression. This is coherent with the physical mechanism for producing propagating lattice excitations, i.e., the impact of a swift particle or the recoil from a decay event as the $\beta$-decay of ⁴⁰K. However similar methods did not produce a soliton-breather but a pure breather in the related model by Archilla et al. [16]. The soliton-breather has a length of 5-6 particles.

Many solutions have a long life, specially, if the initial phonons are absorbed by the borders by computational means. These traveling waves are quasi-exact, meaning that the properties of exact traveling breathers hold very well, although not exactly. Using them as a seed and through the Newton method described above, it is possible to obtain exact traveling waves also composed of a breather and soliton and with small amplitude wings. In Fig. 13 the plot of the particle coordinates $u_{x}$ is presented for three times separated by $20T_{F}$. Both the localization and the exactness can be appreciated. Figure 14 demonstrates the $u_{y}$ coordinates. They have similar spectral properties as $u_{x}$ but their amplitudes are much smaller.

In Table 1 some examples of parameters of exact soliton-breathers are shown. Particular exact traveling wave solutions were obtained on the lattice with dimensions $N_{1}=32$ and $N_{2}=16$. Unfortunately, we have been able to find only steps $s=1$ and we presently do not know if this is a characteristic of the system or it shows the need for very different generation methods to produce different seeds.

We have expanded the theory of exact traveling waves developed for one dimension in a previous publication to two dimensions in a hexagonal lattice, applying it to a model for silicate layers in which some authors have previously found breathers with very long life. The theory is based in the representation in the frequency-momentum space and has allowed the determination of the structure of the propagating waves in the model. They are formed by a breather and a soliton traveling together, that is, a soliton-breather or solbreather. We have described the system in terms of the direct and reciprocal hexagonal lattice, finding the structure of the 1st Brillouin zone including the polarization of the phonon surfaces. The theory describes traveling waves in the $\omega-k$ representation and shows that exact traveling waves lie within parallel planes in that space, each one corresponding to a specific frequency in the moving frame. These frequencies are integer multiples of a minimal one, called the fundamental frequency. One of these planes is the breather plane where the breather lies. The others produce the wings at their intersection with the phonon surfaces. The soliton corresponds to a resonant plane with zero frequency. We have developed a method for obtaining the fundamental time of the traveling waves and used it to obtain exact traveling waves. A variety of them have been found, all of them with small wings and step $s=1$. We are exploring methods to obtain other steps or to find out if they are not possible in our system. The main conclusion is that the theory of exact traveling waves and their spectral representation are powerful means to observe the structure of traveling waves and obtain exact ones. From the physical point of view exact solutions are more likely to appear in physical processes than approximate solutions and their properties can be studied more easily.

J Bajārs acknowledges support from PostDocLatvia grant No.1.1.1.2/VIAA/4/20/617. JFR Archilla thanks projects PAIDI 2021/FQM-280 and MICINN PID2019-109175GB-C22. He also acknowledges a travel grant from VIPPITUS-2020 and the University of Latvia for hospitality. Both authors aucknowledge Prof. M.E. Manley for useful information and discussions about neutron scattering spectra of some materials.