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

We study nonlinear excitations propagating in a hexagonal layer which is a model for the cation layer of silicates. We consider their properties in the frequency-momentum or $\omega-k$ representation, extending the theory on pterobreathers in their moving frame for the first time to two dimensions. It can also be easily extended to three dimensions. Exact traveling waves in the $\omega-k$ representation are within {\em resonant} planes, each plane corresponding in the moving frame to a single frequency. These frequencies are integer multiples of a frequency called the fundamental frequency. A breather is within a resonant plane called the breather plane and has a single frequency in the moving frame. The intersection of the resonant planes with the phonon surfaces produce co-traveling wings with a small set of frequencies. The traveling waves obtained by perturbing the system consist of a breather and a soliton traveling together and are quasi-exact. These traveling waves can be used as seeds to obtain exact traveling waves, also formed by a breather and a soliton. The wings do exist but they are usually very small.


Introduction
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-onedimensional 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, ≃ 27 eV for crowdions and ≃ 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.

Model
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 = 1 2 (ẋ 2 +ẏ 2 ) = 1 2 (u 2 x +u 2 y ), the dot indicating the derivative with respect to time.
The particles are in their equilibrium position in a 2D hexagonal lattice corresponding 2 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 ≤ U (x, y) ≤ 1.
The interaction between particles is given by the Lennard-Jones potential: being ǫ a measure of the relative strength of the interaction potential with respect to the on-site potential. Without loss of generality we consider ǫ = 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 τ = 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 γ > 0. As noted in [20] larger values of γ produce faster moving quasi-exact breathers with larger particle displacements in the direction of propagation. See also Table 1.

The 2D hexagonal lattice
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.

The direct 2D hexagonal lattice
For a 2D hexagonal lattice of unit distance a = 1, the direct vectors that generate the lattice are (see Appendix A and Ref. [22]): with e 1 = [1, 0] and e 2 = [0, 1], the Cartesian unit vectors. Any lattice point can be obtained as a point of the Bravais direct lattice: Blue colors are lower energies, being the minima at the equilibrium particle positions. The red-yellow colors are higher energies as specified by the color bar and correspond to potential barrier between sites. The associated face-centered rectangular lattice can be visualized above abscissas 1 and 2. Note that the breather direction along a close-packed line as the horizontal axis has no perpendicular close-packed plane. See Fig. 2 for a view of U in a single primitive cell.
with n 1 , n 2 integers. Many physical properties are described by a function with the periodicity of the direct Bravais lattice, for example, the on-site potential in Eq. (1) as can be seen in Fig. 1. The minima of the potential correspond to the equilibrium position of the particles and form the hexagonal Bravais lattice. In addition, the form of the potential within a primitive cell can be seen in Fig. 2.
The hexagonal lattice can also be seen as a face centered rectangular lattice with sides 1 and √ 3, which also corresponds to two primitive simple rectangular Bravais lattices displaced half a diagonal vector equal to a 2 . Convenient labels in this description are integers l and m such that the lattice points have the coordinates: In this description l + m is even. When l and m are both odd, they describe the lattice points of a rectangular lattice, and when they are both even, they describe the displaced simple rectangular lattice. This description is used often to write the dynamical equations since the visualization is easier [23,18]. It also provides a convenient matrix form with indexes l, m, which corresponds to the lattice with a shift of alternate rows 0.5 to the left. However, it is not a Bravais lattice and therefore it is not convenient for the description in the momentum or reciprocal space. Both sets of indexes are related by x(l, m) = 1 2 l = n 1 + 1 2 n 2 and y(l, m) =

The reciprocal hexagonal lattice
The corresponding reciprocal lattice basis of the direct Bravais lattice in Eq. (5) is given by: with the properties that |b 1 | = |b 2 | = 2π 2 √ 3 , and a i · b j = 2πδ ij . The (angular) reciprocal space of Eq. (5) is the Bravais lattice expanded by {b 1 , b 2 }, that is, G = m 1 b 1 + m 2 b 3 , with m 1 and m 2 integers, with the property that G · R = 2π(m 1 n 1 + m 2 n 2 ) = 2πn, with n an integer. We can also consider the reciprocal vector Q = G/2π, which has the same meaning as G but in linear (m −1 ) instead of angular (rad·m −1 ) units, similarly to the frequency f = ω/2π and the angular frequency ω. The appropriate basis for Q are the vectorsb 1 In this basis, Q = m 1b1 + m 2b2 with the same values m 1 and m 2 of G. G and ω are sometimes called the physics definition and Q and f the crystallographic definition of the reciprocal vectors and frequency, respectively [24]. Then: e iG·R = e i2πQ·R = 1 .
Therefore, any property with the periodicity of the direct Bravais lattice can be described as a sum of exponentials exp(iG · r). For example, the on-site potential in Eq. (1) is given by: being the simplest form of a potential in the hexagonal lattice consistent with the symmetries of the lattice. The reciprocal vector b 1 + b 2 = 2π[1, 1/ √ 3] is symmetric to b 1 with respect to the horizontal lattice and it is necessary for the hexagonal symmetry of the potential.

The supercell lattice and the wavevector or momentum space
However, as phonons, breathers or kink solutions do not have the symmetry of the lattice (although the whole set of them has it) and we have Born-von Karman periodic conditions, we have to consider the Bravais direct lattice expanded by N 1 a 1 and N 2 a 2 , where N 1 and N 2 are the dimensions of the lattice. It is called the supercell lattice as its primitive cell is usually the computational lattice. The corresponding reciprocal space is the momentum or wavevector space and it is expanded by the basis vectors {b 1 /N 1 , b 2 /N 2 }, or in other form they correspond to: Any wavevector k is equivalent to k + G, with G a vector in the reciprocal lattice. The primitive cell of the wavevector space is visualized by the phonon modes obtained in Sect. 4 as can be seen in Fig. 3. By translations of vectors G, it is often constructed the first Brillouin zone, that is, a zone in momentum space which is closer to a given point than to any other and it is also visualized in the same figure. Important points and paths are often used to visualize the phonon bands as also depicted in the same figure. We can also define the crystallographic wavevectors q = k/(2π) to get rid of the factor 2π, then q = q 1b1 + q 2b2 = q x e 1 + q y e 2 . Note that q x = k x /2π and q y = k y /2π, but q 1 and q 2 are the same in both representations. As they have the reciprocal basis vectors as unit of length they are the reduced wave vectors [25]. Any wavevector k represents a set of planes 1 perpendicular to k separated by a distance 2π/|k| = 1/|q| or a plane wave with plane wavefronts separated by the same distance. Therefore, it is better seen as defining a plane than as a direction. As k + G represents the same set of planes, it is enough to consider a primitive reciprocal lattice with non equivalent vectors k, either within the primitive reciprocal lattice, the Brillouin zone or any other construction through translations by vectors G.

Phonons and the 1st Brillouin zone
In this subsection we will work in the q wavevector space for convenience to avoid a factor of 2π everywhere. Primitive vectors of the reciprocal lattice expanded byb 1 and b 2 represent planes of atoms separated by a distance 1/|b i | = √ 3/2, and perpendicular to them, i.e., at −60 • with the x-axis and horizontal, respectively. Corresponding plane waves propagate in the directions ofb i . Plane waves that propagate parallel to the xaxis do so in the direction of the smallest reciprocal wavevector 2e 1 = 2[1, 0] = 2b 1 +b 2 defining vertical planes separated by a distance 1/2. Therefore, they have wavevectors with q 2 = 2q 1 . The 1st Brillouin zone is a regular hexagon centered at Γ = (0, 0) with sides 2/3, aphotem 1/ √ 3, maximal radius 2/3, with two horizontal sides parallel to the horizontal axis cutting the vertical axis at 1/ √ 3 and two vertices at the x-axis at ±2/3. Important critical points are vertices M = (±2/3, 0) and middle points of the two horizontal sides K = (0, ±1/ √ 3) and their equivalents through several ±60 • rotations. It is convenient to plot the wavevectors in the direction [1,0] with components along that direction. However, (1, 0) is outside the 1st Brillouin zone because planes that bisectb 1 andb 2 cut at vertex K = (2/3, 0). The wavevector [1,0] , but the meaning of the wavevectors are not so clear. By a 2 × 60 • rotations M is also equivalent tob 1 /2.
Plots can be done in the Cartesian coordinates both in the direct and reciprocal space because e 1 = [1, 0] and e 2 = [0, 1] form its own reciprocal basis. We write them as r = xe 1 + ye 2 or in the direct lattice coordinates r = r 1 a 1 + r 2 a 2 . Similarly, wavevectors can be written as k = k x e 1 + k y e 2 = q 1b1 + q 2b2 . Given the values ofb 1 andb 2 above, it is easy to obtain k x = 2πq 1 and k y = 2π(q 2 − 1 2 q 1 ) 2 √ 3 .

Frequency and polarization of linear modes
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: w = [u x , u y ] T is the vector of displacements in the x and y directions of a particle of the lattice and A = [a, b] T is its polarization vector. Without loss of generality we can choose a and b such that √ a 2 + b 2 = 1 and a ≥ 0. In principle a and b could be complex, but in fact, in our system, they can always be chosen real. k is a wavevector and R = n 1 a 1 + n 2 a 2 is the equilibrium position corresponding to a given particle. The angular frequency ω can be chosen positive without loss of generality because the propagation direction of the wave is given by k.
In terms of these indexes the equation above becomes: leading to the following equation [18]: where κ = 16π 2 /9 = 17.5460 is the minimal squared frequency, derived from the linearization of the on-site potential in Eq. (1), and β = V ′′ LJ (1) = 72ǫ = 3.6. We obtain the homogeneous system of linear equations: where Note that A ≥ κ and B ≥ κ 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: As ω > 0, there are two phonon bands corresponding to the two signs of ±. Let us suppose initially that D = 0, which also implies that a = 0. For the frequencies ω 1,2 , the two equations in Eq. (13) are equivalent and only one is needed, for example, the first one. Let us denote α = b/a, then The normalized polarization vector becomes a 1,2 = 1/ 1 + α 2 1,2 and b 1,2 = α 1,2 / 1 + α 2 1,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( √ 3/2k y ) = 0. In the first one sin(k x /2) = 0, which implies that k x = 0 as k x = 2πm is equivalent, then: If k x = 0, then q 1 = 0 and 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 where the + or − in ± corresponds to m = 0 and m = π, respectively. The curves corresponding to (±) = +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 (±) = −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 10 this way, it tends not to perturb the adjacent chains. Note also that for D = 0, D and therefore α 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.

Basic theory
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.

Exact traveling waves in one dimension
The concepts explained below are illustrated in Figs. 7 and 8, where the traveling wave along the central close-packed chain is represented in the real and frequency-momentum spaces, respectively.

Traveling wave
It is represented by u(n, t) = f (n − V b t, ω MF t), being 2π periodic in the second argument. If it is partially localized in the first argument, it becomes a traveling localized wave. The frequency ω MF is the frequency in the moving frame, where it becomes the only frequency of a stationary profile.

Solitons, kinks and breathers
If ω MF = 0 and if f (±∞, 0) = 0, then u represents a soliton. If ω MF = 0 and f (·, 0) is only zero at +∞ or −∞ and a constant value at the other infinity, u represents a kink. If ω MF = 0 and f is localized in the first argument, u represents a breather. The profiles of a breather and a soliton are represented in Fig.7.

Exact traveling wave
If there is a minimal time T F and integer s such that u(n+ s, t+ T F ) = u(n, t), then u is an exact traveling wave. A quasi-exact traveling wave can be seen in Fig. 7 Fundamental time T F and step s The parameters T F and s are denominated fundamental time and step, respectively. They are related by the velocity of propagation V b = s/T F .

Fundamental frequency
It is defined as ω F = 2π/T F . The condition of an exact traveling wave applied to u(n, t), that is, u(n + s, t + T F ) = u(n, t), implies that ω MF = mω F , where m is an integer. More generally, an exact traveling wave is given by a sum of functions u with moving-frame frequencies that are integer multiple of ω F .

Resonant plane waves and resonant lines
For a given step s and T F and therefore velocity V b = s/T F , resonant plane waves are plane waves u = exp(i[kn − ω L t]) that are exact with the same step s and T F . Therefore, they can be cast in the form u = exp(ik(n − V b t)) exp(−imω F ). The moving frame frequencies mω F are related to the laboratory frequencies ω L by: In the ω − k representation the above equation represents straight lines called resonant lines that cut the vertical axis k = 0 at the moving frame frequencies ω MF at integer multiples of ω F . They can be seen in Fig. 8.

Breather line
An exact traveling solution is a sum of resonant plane waves. Therefore, in its XT Fourier transform (XTFT) most of the intensities lie on a resonant line, called the breather line or on a few of them. If there is only one, then the breather line cuts the axis k = 0 at the moving frame frequency ω MF = m b ω F and we recover a single frequency of a breather as in the common stationary case. See Fig. 8.

Wings
Often, there are intensities of the XTFT of an exact traveling wave at the crossing points of a resonant line and the phonon dispersion relation. This means that the traveling waves travel together with one or several resonant phonons. They are called wings, and should not be confused with tails, that are the diminishing amplitudes from the core of a traveling wave, because the wing amplitudes tend to be constant far from the core of the traveling wave. See Fig. 8.
Rest or π frequency Breathers in hard potentials typically derive from a maximum of ω at k = π. The breather mode at k = π does not translate and, therefore, its frequency is called the rest frequency or π-frequency. Often, moving breathers are obtained by perturbing a stationary breather with wavenumbers centered at π and they are a nonlinear perturbation of the π phonon. Its value is ω π = (m b + s/2)ω F , and therefore, if s is odd, it is a semi-integer multiple of ω F . For breathers derived from a ω(k) minimum at k = 0, the rest frequency coincides with the frequency in the moving frame.

Exact traveling waves in two dimensions
Most of the theory in two dimensions is similar to one dimension, with some obvious changes. The variable u may have two components (u x , u y ) and depends on two indexes in the plane, i.e., u = u(n 1 , n 2 , t), with n 1 and n 2 indexes in a sublattice of the Bravais lattice. These concept are illustrated in Fig. 9 as a projection on the ω − k x plane and in three dimensions in Fig. 10.

Exact traveling waves in 2D
If we suppose that the direction of a propagation is along n 1 , then a traveling wave is represented by u = f (n 1 − V b t, n 2 , ω MF t), being f localized in the first two variables and 2π periodic in the third variable. It is exact if u(n 1 +s, n 2 , ω(t+T F )) = u(n 1 , n 2 , t) for some step s and fundamental time T F .

Resonant planes and breather plane
Resonant lines become resonant planes in the ω − k space. The breather line becomes the breather plane which is parallel to the n 2 direction. The intensities 14 of the XY T F T are within the breather plane(s). If the traveling wave is very much localized along a given direction, the plane may be formed from parallel lines perpendicular to it in k-space. They can be seen in Fig. 10. If the Fourier transform is done only on the n 1 variables, then the breather lines reappear as can be seen in Fig. 11.

2D wings
Wings may appear at the intersections of the resonant planes and the phonon surfaces. They are therefore more complex than in 1D, since they tend to involve more wavevectors, but the frequencies in the moving frame are still a discrete small set corresponding to the few resonant planes that cut the phonon surfaces. See Fig. 10.

Travelling localized waves in three dimensions and in physical crystals
There is no special difficulty to extend the theory to three dimensions, although the visualization and understanding becomes problematic as the frequency-momentum space has four dimensions. Resonant planes becomes hyperplanes and so on. Visualization in 3D will be projections of the 4D ω − k space.
Then, we can think at how quasi-exact travelling localized waves might appear in actual spectra of physical crystals. Most likely, they will be not so much localized as in simulations, therefore they will be some platelet within a plane which is close to tangent to some slope around a maximum or minimum of the phonon bands. There will be not a single travelling wave but the whole spectrum of them according to their probability of formation. With different velocities and directions they will form a thick truncated cone half filling a valley or a mountain-top in the phonon hypersurfaces. Similar structures have been observed in spectra of materials as SePb at high temperature [21] (See Fig. 2 in that reference), which seems promising.
6. Quasi-exact soliton-breather in the central row. Finding T F .
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 γ. 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 ∈ {0 : N 1 − 1}, i.e., u n (t l ), with t l = l/N t , and l ∈ {0 : N t − 1}. From the breather line in the XTFT we obtain the breather velocity V b = ∆ω/∆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 = round(T s /τ )τ , where 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 τ = 0.04, t 1 = round(T 1 /τ ) = 76τ = 2.9600 = 3.1600 and t 2 = round(T 2 /τ ) = 152τ = 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 ω MF = m b ω F = m b 2π/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 τ .

Exact traveling solutions and the Newton method
Breathers obtained by providing some recoil velocity are interesting as this is a likely mechanism from β decay, as in 40 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 s = [s 1 , s 2 ], where s 1 and s 2 ar integers. Let 16 The action of the operatorT is obtained by integrating the equations of dynamics of the system. The composite translation map in the lattice and timeŜ[m, T ] becomes: Therefore an exact traveling solution with fundamental time T F and step s is given by The election of t = 0 is irrelevant as the Hamiltonian and dynamical equations are time invariant. Let as define the function f : The exact solution is obtained by the Newton method. Let as suppose that U 0 is an approximate solution, a seed, that is, f (U 0 ) = 0 but small. The approximate 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 U = U 0 + δU is an exact solution close to U 0 , i.e., f (U) = f (U + δU) = 0: where ∇f is the Jacobian of f calculated at the approximate solution U 0 . The new initial variables become U = U 0 + δ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 a 1 and the step becomes s = [s 1 , 0]. For simplicity, we will refer often to the step just as the scalar s ≡ s 1 . Typical values of s are very low integer numbers 1,2,...

Exact soliton-breathers
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 β-decay of 40 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.

Conclusions
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 ω − 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.