One-component delocalized nonlinear vibrational modes of square lattices

All possible one-component delocalized nonlinear vibrational modes (DNVMs) in a square lattice are analyzed. DNVMs are obtained taking into account exclusively the symmetry of a square lattice, and therefore, they exist regardless of the type of interactions between particles. In this work, the interactions of the nearest and next nearest neighbors are described by the β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document}-FPUT potential. For each DNVM, frequency, kinetic and potential energies are described as functions of amplitude. The mechanical stresses caused by DNVMs and the effect of DNVMs on the stiffness constants of the lattice are presented. DNVMs with higher vibration frequencies have a stronger effect on the mechanical properties of the lattice. Examples of analytical analysis of DNVMs are given. It is found that two of the sixteen one-component DNVMs can have frequencies above the phonon band in the entire range of their amplitudes. It is shown that these DNVMs can be used to construct discrete breathers by applying localizing functions. The modulational instability of these two DNVMs leads to the formation of chaotic DBs. The presented results contribute to a better understanding of the nonlinear dynamics of a square lattice by analyzing the properties of a class of delocalized exact solutions and demonstrating their connection with discrete breathers.

of the nearest and next nearest neighbors are described by the β-FPUT potential. For each DNVM, frequency, kinetic and potential energies are described as functions of amplitude. The mechanical stresses caused by DNVMs and the effect of DNVMs on the stiffness constants of the lattice are presented. DNVMs with higher vibration frequencies have a stronger effect on the mechanical properties of the lattice. Examples of analytical analysis of DNVMs are given. It is found that two of the sixteen one-component DNVMs can have frequencies above the phonon band in the entire range of their amplitudes. It is shown that these DNVMs can be used to construct discrete breathers by applying localizing functions. The modulational instability of these two DNVMs leads to the formation of chaotic DBs. The presented results contribute to a better understanding of the nonlinear dynamics of a square lattice by analyzing the properties of a class of delocalized exact solutions and demonstrating their connection with discrete breathers.
Keywords Square lattice · Nonlinear dynamics · Delocalized nonlinear vibrational mode · Exact solution

Introduction
Interest in nonlinear lattice vibrations has increased in recent decades due to the fact that crystalline materials are exposed to high-amplitude influences in many applications. One of the effects of nonlinearity in discrete periodic structures is the possibility of the existence of spatially localized oscillations of large amplitude, called discrete breathers (DBs) (or intrinsic localized modes) [1][2][3][4][5][6]. DBs contribute to the macroscopic properties of crystals [6,7], reduce thermal conductivity due to phonon scattering [8] and cause the formation of lattice defects [9][10][11].
Nonlinear lattices also support exact solutions in the form of delocalized nonlinear vibrational modes (DNVMs) [12][13][14], which are much less studied. On the other hand, there is a close relationship between DNVMs and the DBs [6,15], and this is one of the advantages of learning about the former.
In [12][13][14], a method was developed for finding exact delocalized vibrational solutions for nonlinear lattices, based on the analysis of the lattice symmetry. Such solutions were referred to in original papers as bushes of nonlinear normal modes (BNNMs), and in some later studies as DNVMs [16][17][18]. In particular, the BNNM theory was used to obtain DBs in a scalar nonlinear square lattice [15].
DNVMs can be multi-component. One-component DNVM is periodic in space and in time, having single vibration frequency, while m-component DNVM has m degrees of freedom with m incommensurate frequencies. The dynamics of the m-component DNVM can be described by m coupled dynamical equations.
BNNMs can be constructed for molecules, considering their point symmetry, and for nonlinear lattices, based on their space symmetry. Examples are one-and two-component BNNMs in the SF 6 molecule [19] and in nonlinear chains [20][21][22].
At first glance, it may seem that DBs and DNVMs are of a different nature, but it turns out that there is a relationship between them. DNVMs found for a chain and for a 2D triangular lattice were employed to construct linear and planar DBs in fcc metals [23][24][25]. DNVMs were used to find DBs of new type in monoatomic triangular Morse lattice [26,27], in triangular β-Fermi-Pasta-Ulam-Tsingou (FPUT) lattice [28] and in bcc metals [29]. DNVMs in hexagonal lattice [30] were studied in 2D materials [31], for example, in graphene [16,[32][33][34] and boron nitride [35]. Chaotic DBs can emerge in the lattices as a result of modulational instability of DNVMs [36][37][38][39][40][41][42]. The effect of chaotic DBs on macroscopic properties of triangular lattice was analysed [43]. DNVMs can be used to solve the important problem of checking the accuracy of interatomic potentials in molecular dynamics simulations [18].
This study is a continuation of the work on DNVMs in a triangular lattice [17], but here we describe onecomponent DNVMs in a square β-FPUT lattice. All DNVMs described in in this work exist in any square lattice because they are determined only by lattice symmetry. On the other hand, properties of DNVMs depend on the interparticle potential. The details of the DNVM derivation are described in [12][13][14] and outlined in Appendix A to this paper. BNNM theory is briefly presented in Appendix B.
In the rest of the paper, the model is described in Sect. 2. Then we present the phonon dispersion relation for the lattice in Sect. 3, describe one-component DNVMs in Sect. 4, give examples of analytical treatment of DNVMs in Sect. 5 and present numerical results on DNVM properties in Sect. 6. In Sect. 7 discrete breathers are obtained using localizing functions, and in Sect. 8 the appearance of chaotic discrete breathers as a result of modulation instability of delocalized modes is demonstrated. The results are summarized and future problems are outlined in Sect. 9.

Square FPUT lattice and its mechanical properties
A two-dimensional square β-FPUT lattice with nearest and next-nearest neighbor interactions is considered, see Fig. 1. The lattice with the lattice step h is a set of points in the x y plane with radius vectors where i and j are integers and the basis vectors of the lattice are e 1 = (h, 0) and e 2 = (0, h).
The lattice points are occupied by the particles of mass m, each having two degrees of freedom, the components of the displacement vector (u i, j , v i, j ). Position of the particle i, j at time t is given by the radius vector Each particle interacts with the four nearest (n) and four next-nearest (nn) neighbors via the β-FPUT potential ϕ n,nn (r ) = k n,nn 2 (r − a n,nn ) 2 + β n,nn 4 (r − a n,nn ) 4 , where r is the distance between the particles, a n = h and a nn = √ 2h for the nearest and next-nearest bonds, respectively; k n and k nn (β n and β nn ) are the coefficients for the harmonic (anharmonic) part of the potential for the nearest and next-nearest bonds, respectively. We take h = 1 and k n = 1 by choosing the units of distance and energy, respectively. Different values of k nn will be considered. For the anharmonicity coefficients we set β n = β nn = 10, then the nonlinearity effects become noticeable for particle displacements of the order of 0.1. For the particle mass m = 1 is set by choosing a unit of time.
Referring to Fig. 1, we define the following vectors connecting the nearest and next-nearest neighbors of the i, j particle: Computational cell includes I × J particles subject to the periodic boundary conditions, r i, j = r i+I, j = r i, j+J . The Hamiltonian (total energy) of the computational cell is the sum of the kinetic and potential energies of the nearest and next-nearest bonds: where overdot means differentiation with respect to time. The Hamilton's equations of motion derived from Eq. (4) are where The Störmer integrator of order six with the time step 0.002 time units is used for integrating the equations of motion Eq. (5). This symplectic integrator is efficient for the second-order ordinary differential equations not containing the first derivative of the unknown function [55].
Components of the mechanical stress in the lattice can be calculated as follows where S = h 2 I J is the area of the computational cell. Mechanical stress is related to the strain components via the Hooke's law under plane stress conditions [56] σ x x = C 11 ε x x + C 12 ε yy + C 13 ε xy , σ yy = C 21 ε x x + C 22 ε yy + C 23 ε xy , where C i j = C ji are the stiffness constants which can be calculated as follows where Excitation of a DNVM leads to the appearance of internal stresses in the lattice that are the functions of time. Elastic constants also become timedependent. We calculate the averaged over oscillation period stresses and stiffness constants, Then averaged over time mean stress (pressure) is In addition, the kinetic and potential energies per particle averaged over the oscillation period will be analyzed, where N = I J is the number of particles in the computational cell. The size of the simulation cell is taken equal to I = J = 12. With this choice, DNVMs having spatial periods with 2, 3 and 4 primitive translation cells can be described, and there are no one-component DNVMs with a different translation period, as it will be shown later.

Phonon dispersion relation
Assuming that u i j h and v i j h, we expand the right-hand side of Eq. (5) and discard terms with powers greater than linear, obtaining linearized equations of motion Looking for the solution of Eq. (14) in the form where i is imaginary unit and s and q are the wavenumbers, one finds the dispersion relation where γ = 4k n sin 2 q 2 , δ = 4k n sin 2 s 2 , Phonon frequencies at the points (q, s) = (±π, 0) and (0, ±π) are Here ω 1 is the maximal phonon frequency. Nearest and next-nearest bonds contribute equally to the maximal frequency.
Each DNVM is determined by its symmetry group G j (see Appendix A). The generators of the subgroup for each DNVM are listed in Table 1.
All nonzero initial displacement vectors in Fig. 2 have the same length A, which is the DNVM vibration amplitude. In DNVMs 1, 3, 6, 9, 15 and 16 all particles move. A half of particles in DNVMS 4,5,7,8,11,12,13 and 14 have zero initial displacements and they remain at rest at any time while other particles oscillate. The same is true for one third of particles in DNVMs 2 and 10.
In DNVMs from 1 to 8 all moving particles oscillate along a close packed direction (x or y). In DNVMs 13 and 14 moving particles oscillate in both close packed directions. In DNVMs 9-12 the particles move along y = x or along y = −x, while in DNVMs 15 and 16 they move along both these directions.

Analytical treatment of one-component DNVMs
As an example, for DNVMs 1, 6 and 16, we write down the exact Hamiltonian, the equation of motion accurate to cubic terms and obtain the dependence of the frequency on the amplitude within the framework of the cubic approximation, which is valid for not too large amplitudes.

DNVM 1
As it has been said earlier, any one-component DNVM has single degree of freedom, for example, the distance of a moving particle from its equilibrium position, a(t), with a(0) = A andȧ(0) = 0. For example, for DNVM 1 (see Fig. 2), taking into consideration the particle displacement pattern, one obtains the Hamiltonian where Here h is the equilibrium distance between lattice points, L 1 and L 4 are distances between nearest particles, L 2 , L 3 , L 5 and L 6 are distances between nextnearest particles; the interparticle potentials ϕ n (r ) and ϕ nn (r ) are defined by Eq. (2).
Expanding the potential energy Eq. (18) in Taylor series and retaining up to quartic terms one obtains the simplified Hamiltonian The equation of motion corresponding to the Hamiltonian Eq. (20) has the form Looking for the solution to Eq. (21) in the form a(t) = A sin(ωt) + A 1 sin(3ωt), where it is assumed that A 1 A, one finds the relation between frequency and amplitude in the form then DNVM 1 will demonstrate a soft type of anharmonicity even with positive values of β n and β nn . This unexpected result is explained by the competition between geometric and physical nonlinearity. Geometric anharmonicity for DNVM 1 leads to a decrease in frequency with amplitude (soft type of anharmonicity), and physical nonlinearity has the opposite effect. At small positive values of β n and β nn , the soft type of anharmonicity prevails and the oscillation frequency decreases with amplitude.
For the parameters used in our work and k nn = 1, one finds from Eq. (22)

DNVM 6
Similarly, for DNVM 6 (see Fig. 2), taking into consideration the particle displacement pattern, one obtains the Hamiltonian where h is the equilibrium distance between lattice points, L 1,2 are distances between nearest particles, and ϕ n (r ) is the interparticle potential, in our case, it is given by Eq. (2). Note that in DNVM 6 the next-nearest bonds are not deformed and the potential ϕ nn (r ) does not enter the Hamiltonian Eq. (25). Expanding the right side of Eq. (25) into a Taylor series and keeping the terms up to the fourth power, one obtains a simplified Hamiltonian in the form The corresponding cubic equation of motion is The approximate solution to Eq. (27) can be written in the form a(t) = A sin(ωt) + A 1 sin(3ωt), where it is assumed that A 1 A, and one finds the relation between frequency and amplitude in the form For the parameters used in our work and k nn = 1, one finds from Eq. (28)

DNVM 16
Similarly, for DNVM 16 (see Fig. 2), taking into consideration the particle displacement pattern, one obtains the Hamiltonian where Here h is the equilibrium distance between lattice points, L 1 ,L 3 ,L 5 and L 7 are distances between nearest particles, L 2 ,L 4 ,L 6 and L 8 are distances between nextnearest particles; the interparticle potentials ϕ n (r ) and ϕ nn (r ) are defined by Eq. (2). Expanding the potential energy Eq. (30) in Taylor series and retaining up to quartic terms one obtains the simplified Hamiltonian The equation of motion corresponding to the Hamiltonian Eq. (32) has the form Looking for the solution to Eq. (33) in the form a(t) = A sin(ωt) + A 1 sin(3ωt), where A 1 A, one finds the relation between frequency and amplitude in the form For the parameters used in our work and k nn = 1, one finds from Eq. (34) Similar analytical calculations can be easily performed for any other one-component DNVM.
The analytical results (lines) are compared to the numerically found frequencies (circles) for DNVM 1 (black) and DNVM 6 (blue) in Fig. 3. Here k nn = 1. For small amplitudes the agreement is very good and the error increases with increasing A but it does not exceed 1.3% for A < 0.3.

Numerical results
Using molecular dynamics method we analyze properties of DNVMs.

Frequency response of DNVMs
Exciting DNVM and integrating the equations of motion (5), we determine the dependence of the dis-  15) and (16). First of all, we note the general tendency of an increase in vibra- Note that only DNVM 1 and 16 have frequencies above the phonon spectrum in the entire range of vibration amplitudes for the considered parameters of the model. The rest of the DNVMs at low amplitudes have frequencies within the phonon spectrum. For the selected parameter values, all DNVMs show a hard type anharmonicity, that is, an increase in frequency with amplitude. Recall that according to the condition (23), for sufficiently small values of the parameters k n and k nn , DNVM 1 will show a soft type of anharmonicity and, therefore, its frequency will be in the phonon spectrum. Different DNVMs can have the lowest frequency, depending on the parameter k nn . So, for k nn = 0.3 this is DNVM 5, and for k nn = 3.0 this is DNVM 11.
The frequency response of DNVM 6 does not depend on k nn (see Fig. 4) in full accordance with the analytical calculations performed in Sect. 5.2, since the Hamiltonian Eq. (25) does not contain the potential φ nn .
In Fig. 4 one can see that the amplitude-frequency dependences of some modes merge as A tends to zero. With a change in the parameter k nn , different modes can have the same frequencies in the limit A → 0. However, some modes have the same frequencies in the low-amplitude limit for any values of the model parameters. Obviously, at A → 0, DNVMs are transformed into phonon modes. Let us ask ourselves what wave vectors DNVM has in the low-amplitude limit. Having the dispersion relations (15) and (16), we determine the points of the first Brillouin zone for all DNVMs in the low-amplitude limit. In addition, let us determine for each mode which of the branches of the dispersion relation it belongs to-longitudinal or transverse. Figure 5 shows the DNVM wave vectors in the limit of small amplitudes, when they transform into phonon modes. Some DNVMs are grouped in pairs. For example, Fig. 5a shows the wave vectors of DNVMs 1 and 16. DNVM 1 has a single wave vector, say, (q, s) = (π, 0), while DNVM 16 is nothing more than a superposition of two phonons with wave vectors (q, s) = (π, 0) and (q, s) = (0, π). Similarly, Fig. 5f shows the wave vectors of DNVMs 6 and 9. DNVM 9 has a single wave vector, for example, (q, s) = (−π, π), and DNVM 6 is a superposition of two phonons with wave vectors (q, s) = (−π, π) and (q, s) = (π, π ).
It is worth noting that DNVMs 6 and 9 in Fig. 5f are shown in red as longitudinal modes, but they are equally transverse modes and could be shown in blue. The point is that these modes have the same frequencies ω 1 and ω 2 at the points (q, s) = (±π, ±π) for any values of model parameters.

Kinetic and potential energies of DNVMs
The energy of low-amplitude oscillations is equally shared between the kinetic and potential components; however, an increase in the amplitude of oscillations in nonlinear systems leads to a deviation from the equality of these energies. Thus, the difference between kinetic and potential energies serves as a measure of the anharmonicity of oscillations.
In Fig. 6, as the functions of the oscillation amplitude, the kinetic (red lines) and potential (blue lines) For all sixteen DNVMs, as expected, the difference between kinetic and potential energies increases with increasing amplitude A. For all modes, the kinetic energy is higher than the potential energy, which reflects the hard type of mode nonlinearity. Indeed, an increase in the vibration frequency is associated with an increase in the average particle velocities and, consequently, the average kinetic energy.
The rate of growth of energy with amplitude is higher for modes with higher frequencies, as can be seen from the comparison of Figs. 4a and 6. The highest-frequency DNVMs 1 and 16 have the highest energy. The lowest energy, as well as the lowest frequency, has DNVM 5. Since 6 represents energy per particle, its growth rate also depends on the fraction of particles that oscillate. For example, in DNVM 6 all particles vibrate, while in DNVM 7 only half, therefore, despite the fact that the frequency of DNVM 7 is slightly higher than the frequency of DNVM 6, the latter has a higher energy per particle.

Mechanical stresses produced by DNVMs
The stresses calculated from Eq. (7) and averaged according to Eq. (11) are shown in Fig. 7 for sixteen DNVMs for the case k nn = 0.3. The angle brackets are omitted for brevity. Stresses σ x x are shown with red solid lines, σ yy with blue dashed lines, and σ xy with green solid lines.
Shear stress (σ xy ) is not equal to zero only for DNVMs from 9 to 12. This is due to the fact that in these modes all particles vibrate in the same direction, along the bonds connecting the second neighbors. Interestingly, the shear stress in these modes differs only slightly from the normal ones (σ x x and σ yy ).
Equality of normal stress components is observed only for DNVMs from 13 to 16. This is a consequence of the fact that the symmetry groups of these DNVMs (see Fig. 2 and Table 1) contain a fourfold rotation axis (φ), and therefore, corresponding symmetry-invariant second-rank tensors must be proportional to the unit tensor.
DNVMs 5, 11 and 13 create very low lattice stresses due to their low vibration frequencies, see Fig. 4a. The stiffness constants of the equilibrium square lattice, calculated by Eq. (9), are given Table 2 for three values of the stiffness parameter of the bonds between the next-nearest neighbors, k nn .
As can be seen from Table 2, C 11 = C 22 due to the symmetry of the square lattice. The shear stiffness of the lattice is lower than the tensile stiffness: C 33 < C 11 = C 22 . With an increase in the stiffness of the bonds connecting the second neighbors, the difference between the shear and tensile stiffness decreases.
Let us consider how excitation of DNVM affects the lattice stiffness constants averaged over the oscillation period according to Eq. (11).
In Fig. 8, the dependencies of the stiffness constants averaged over the oscillation period on the DNVM amplitude are plotted. The red, blue dashed and green curves show C 11 , C 22 and C 33 , respectively. The angle brackets denoting time averaging are omitted for brevity. In this example, k nn = 0.3.
As follows from Fig. 8, all DNVMs increase the stiffness of the lattice, and their effect on stiffness increases Table 2 The stiffness constants of the equilibrium square lattice given by Eq. with the DNVM amplitude. The only exception is the value of C 22 for DNVM 6, which slightly decreases with an increase in A, however, this is accompanied by a fairly rapid increase in stiffness C 11 . A decrease in C 22 with increasing amplitude is observed for DNVM 6 not only for k nn = 0.3, but also for k nn = 1.0 and 3.0. Note that the shear stiffness C 33 in the entire studied range of amplitudes for all DNVMs remains below the stiffness C 11 and C 22 .
For DNVMs from 1 to 8 for A > 0, C 11 > C 22 is fulfilled, which is due to the fact that atoms in these modes vibrate along the x axis. However, they can fluctuate along the y axis, and then the relationship between the constants C 11 and C 22 will be reversed.
For DNVMs from 9 to 16, C 11 = C 22 for any A. This is because atoms in these modes vibrate equivalently with respect to the x and y axes.
DNVMs 5, 11 and 13 have very little effect on the stiffness of the lattice, which is associated with their low vibration frequencies, see Fig. 4a. High-frequency DNVMs 1 and 16 have the greatest influence on the lattice rigidity.

Discrete breathers obtained by using localizing functions
Let us show that DBs can be obtained by applying the localizing functions on DNVMs 1 and 16 with frequencies above the phonon spectrum, in the same way as it was done for a triangular lattice [28]. The localizing function is taken in the form where a i j is the magnitude of the initial displacement of particle with coordinates (x i j , y i j ), A is the DNVM amplitude, β x and β y are the localization parameters, and the coordinates (x 0 , y 0 ) determine the position of

Chaotic discrete breathers
DNVMs with sufficiently large amplitude are unstable [34,[36][37][38][39][40][41][42]. If the DNVM frequency is outside the phonon band of the lattice, its energy cannot be directly transferred to delocalized phonons, and the development of instability occurs through energy localization on chaotic DBs.
The degree of energy localization in the lattice can be measured by the localization parameter where e n is the energy of the nth particle and N is the number of particles in the computational cell. If all particles have the same energy, L = 1, and if all the energy of the system is concentrated on one particle, then L = N . In Fig. 11 L(t) is shown for a DNVM 1 and b DNVM 16. Three values of DNVM amplitude are considered: A = 0.02 (red), 0.03 (blue) and 0.04 (green). All curves show a sharp increase in L after a period of time that increases with decreasing A. This increase in L is associated with the formation of chaotic DBs in the lattice. After L reaches a maximum, which also increases with decreasing A, the localization parameter begins to decrease. The decrease in L is associated with the emission of DB energy in the form of low-amplitude waves. Eventually, the system reaches a state of thermal equilibrium, which is characterized by a low value of L.
The energy distribution in the lattice is shown in Fig. 12 for a,b DNVM 1 and c,d DNVM 16 with amplitude A = 0.03. The particles are colored according to their energy. White (black) color corresponds to the minimum (maximum) lattice energy. Panels a and c show the initial stage of energy localization. Panels b and d show the energy distribution when the localization parameter reaches its maximum. One can see strongly localized chaotic DBs.

Conclusions and future challenges
In this study, the properties of all 16 possible onecomponent DNVMs in a square lattice with nearest and next nearest interactions via the β-FPU potential were analyzed using molecular dynamics and analytics. The methods of derivation of DNVMs are outlined in Appendices A and B. DNVMs exist as exact solutions in square lattices with any interparticle interactions and for any amplitudes, since only the symmetry of the lattice was taken into account in their construction. However, the properties of DNVMs depend on the interaction potential.
Examples of analytical analysis of DNVMs are presented in Sect. 5, where the equations of motion and amplitude-frequency relations were derived for DNVMs 1, 6 and 16.
DNVMs with frequencies above the phonon spectrum are important in the study of discrete breathers [26,27,41]. Numerically (see Fig. 4) and analytically (see Sect. 5) it was found that only DNVMs 1 and 16 have frequencies above the phonon spectrum in the entire range of amplitudes, and for DNVM 1 this is true only if the condition Eq. (23) is not satisfied.
Wave vectors of DNVMs in the limit of small amplitudes were identified, see Fig. 5.
Kinetic and potential energies per particle, averaged over the DNVM oscillation period were calculated and presented in Fig. 6. The difference between kinetic and potential energies characterizes the degree of anharmonicity of the DNVMs. As expected, the difference between kinetic and potential energies increases with the amplitude. For all DNVMs, the kinetic energy is higher than the potential energy due to the hard type nonlinearity of DNVMs.
The stress tensor components averaged over the oscillation period are presented in Fig. 7. The relationship between normal stresses σ x x , σ yy and shear stress σ xy reflects DNVM symmetry. DNVMs with higher frequencies cause higher stresses in the lattice.
Stiffness constants of the lattice averaged over the oscillation period are reported in Fig. 8. DNVMs with higher frequencies stronger affect the lattice stiffness.
Examples of constructing discrete breathers by imposing a localizing function on the DNVMs with frequencies above the phonon spectrum are presented in Sect. 7.
The appearance of chaotic DBs as a result of the modulation instability of DNVMs with frequencies above the phonon spectrum was modeled in Sect. 8.
In the present study, only one-component DNVMs were analyzed. It is important to study the properties of multicomponent DNVMs in square lattices with different interaction potentials. DNVMs with transverse vibrations of particles (normal to the plane of the lattice) will also be analyzed.

Data availability
The data that support the findings of this study are available on request from the corresponding author, S.V.D.

Conflict of interest
The authors declare that they have no conflict of interest.

Appendix A: Some comments on methods for studying nonlinear lattice dynamics
We would like to comment on several papers related to the study of the nonlinear dynamics of a square lattice.
An important concept in nonlinear dynamics is the concept of an invariant manifold (or invariant set). This is such a set of points in the phase space of the system, starting from which the phase trajectory does not leave this set. Any exact solution of dynamical equations must belong to some invariant manifold, so finding invariant manifolds is an important step in finding exact solutions to such equations.
Currently, there are no methods for finding invariant manifolds for dynamical systems of general form. However, in the theory of BNNMs, such a method was developed for arbitrary dynamical systems with discrete symmetry [14]. This is a specific group-theoretical method based only on the knowledge of the symmetry of the system, regardless of the specific type of its interparticle interactions.
In Ref. [15], this method was applied to find scalar DBs in the square lattice. However, the same method can be used to search for DBs of other types, in particular, those with particle oscillations in the lattice plane.
For a Hamiltonian dynamical system, the space symmetry is determined by the set of all spatial transformations that do not change its Hamiltonian. The complete set of such transformations forms a certain symmetry group G 0 . In the case of molecules, this is one of 32 point groups, and in the case of three-dimensional crystals, it is one of 230 space groups.
G 0 is the symmetry group of the system's equilibrium state, while its subgroups G j ⊂ G 0 determine the symmetry of particle displacements at general instant, which take place during atomic vibrations around their equilibrium positions. The main proposition of the theory of BNNMs is the assertion that each symmetry group G j determines some invariant manifold of the studied dynamical system [13,14]. Thus, the problem of finding all symmetry-determined invariant manifolds reduces to finding all atomic configurations that correspond to each of the subgroups of the original symmetry group G 0 .
The symmetry group G 0 of the square lattice is the plane space group C 1 4v = P4mm. Invariant manifolds, on which stationary DBs lie, must have point symmetry group, i.e., such group that remains a fixed point relative to breather oscillations, which determines the lattice site of the breather localization. In crystallographic handbooks (see for example [57]), one can find lists of so-called Wyckoff positions (WPs) for all space groups. Particles of the crystal are located at these positions. Wyckoff positions correspond to the lattice in equilibrium. However, during vibrations, particles are displaced from these equilibrium positions, and instantaneous configuration of the moving atoms is already described by some subgroup G j of the original space group G 0 .
Hereafter, we will refer to these Wyckoff positions as WP1, WP2 and WP3, respectively. For each WP, the above formulas give the following information: multiplicity, label, local point group and localization of the fixed point of this group in a square primitive cell.
The primitive cell of the lattice is a translation cell of the minimum volume. An extended primitive cell, consisting of M × N primitive translational cells, is considered to find the particle displacement pattern Q M N corresponding to the exact solution of the dynamical equations of motion. We must consider all subgroups of local symmetry groups of WP1, WP2 and WP3 demanding that the pattern Q M N of dynamical variables q i j (t) is invariant with respect to these subgroups. As a result, we will obtain all possible symmetry-determined invariant manifolds. The restriction on Q M N to be invariant under the action of a symmetry group G j ⊆ G 0 leads to certain relations between the dynamical variables q i j (t): some of them can be equal or differ from each other only in sign.
In [15], we studied only strongly localized DBs. This allowed us to deal with relatively small fragments Q M N of each invariant manifold. Specific values of M and N were chosen from the condition that amplitudes q i j (t) for the peripheral sites must be much smaller than those for the breather core.
The explicit forms of all possible symmetrydetermined patterns of scalar DBs in the square lattice, as well as the explicit forms of breathers for several dynamical models, can be found in Ref. [15].
There is a deep connection between localized and delocalized symmetry-related invariant manifolds, which is determined by the so-called induced representations of space symmetry groups [57]. Discussion of this issue goes far beyond the scope of the Appendix A.
In the theory of BNNMs, symmetry-determined dynamical objects, both localized and delocalized, are studied by the same group-theoretical method.
Discrete symmetry of the system makes it possible to find a wide class of invariant manifolds, which contain some exact dynamical objects.
Of course, the method under discussion cannot give us all possible exact solutions of the considered nonlinear dynamical system, but by applying this method we obtain a wide class of exact solutions. The structure of these solutions is universal in the sense that it does not depend on the specific type of interparticle interactions in the system under consideration. However, the dynamical properties of the considered dynamical objects, such as frequencies and stability, depend significantly on interparticle interactions.
Despite the fact that only oscillatory bushes were considered in this paper, we emphasize that one can also speak about bushes of a different physical nature. For example, there exist bushes of spin modes and even mixed (multiferroic) bushes, some of whose modes are vibrational, while others are spin modes [14]. Bushes can exist not only in molecules and crystals, but also in other structures, such as different materials described in Ref. [58].
Thus, finding all possible patterns of atomic displacements, similar to those shown in Fig. 2, reduces to choosing invariant manifolds corresponding to different subgroups of the symmetry group of the crystal under consideration.
When selecting subgroups of the space symmetry group, it is necessary to consider extended cells of a sufficiently large size so as not to lose some bushes of the specified dimension. In the present work, we study only one-dimensional bushes (Rosenberg modes), and to achieve this goal, it is necessary to consider extended cells, the edges of which are 2, 3 and 4 times larger than those of the primitive cell of the square lattice.
Group-theoretical methods for finding subgroups of space groups have been described in [59], and an elementary introduction to the problems under discussion can be found in [60].

Appendix B:Theory of BNNMs for large-amplitude atomic vibrations
Two aspects of the theory of BNNMs should be discussed, namely the geometric and dynamical aspects.
Geometric aspect Let us start with conventional normal modes (NMs), which are introduced in the harmonic approximation and are independent of each other in this approximation. When small anharmonic terms are taken into account in the expansion of the potential energy of the system in terms of atomic displacements from equilibrium positions, this independence is lost, and NMs begin to interact with each other.
It turns out that in systems with discrete symmetry, there are some selection rules for the transfer of energy between NMs of different symmetry. These selection rules lead to the possibility of the existence of dynamical objects consisting of a certain number of vibrational modes, the collection of which is preserved in time, and which constantly exchange energy with each other. For a number of reasons, discussed below, these objects have been called BNNMs.
Each bush singles out some invariant manifold, which is a subspace of the phase space spanned by some eigenvectors of the force constants matrix. These vectors define an atomic displacements pattern corresponding to some normal mode. Each invariant manifold corresponds to a symmetry subgroup G j of the symmetry group G 0 of the considered dynamical system (G j ⊆ G 0 ). As a result, we are able to find the geometric structure of the bush by solving systems of linear algebraic equations of the form which represent conditions of invariance of the bush, as a pattern of atomic displacements X, with respect to a subgroup G j of the group G 0 . The number of modes of a given bush defines its dimension. Figure 2 of the main text of this paper shows patterns of atomic displacements for all one-dimensional bushes in a plane square lattice, where the arrows show the displacements of individual atoms from their positions corresponding to the equilibrium state.
Below we would like to answer some questions that readers may have.

Why do we claim that exactly 16 one-dimensional bushes can exist in the system under consideration?
The point is that the system Eq. (B1) is a homogeneous system of linear algebraic equations whose general solution depends on a certain number (m) of arbitrary constants. In our case, this number determines the dimension of the bush.
If we consistently handle all subgroups G j of G 0 in Eq. (B1), we will obtain all bushes of different dimensions (m j ). The value of this dimension can be an integer ranging from one to the full dimension of the considered system. In this way we find all bushes of all possible dimensions.
In the case of a plane square lattice, it turned out that there are exactly 16 one-dimensional bushes, because all other bushes possess higher dimensions. This is the answer to the above question.
It should be emphasized once again that in the case of arbitrary dynamical systems, there are no general methods for searching invariant manifolds, while for systems with discrete symmetry such methods were developed in the framework of the theory of BNNMs.
It should be noted that the apparent simplicity of Eq. (B1) is illusory. It is rather easy to find all subgroups of any point group of molecule, while for crystals the situation becomes much more complicated. Indeed, in this case, we have to deal with space symmetry groups whose subgroups G j can be embedded in the space group G 0 in different ways. This fact makes the problem rather non-trivial, especially for the case of complex space groups G 0 . For this reason, we have developed some special group-theoretical methods, in particular, the method of invariant (stationary) vectors of irreducible representations (irreps) of space groups [59].

How does the Wigner classification of normal modes relate to the theory of BNNMs?
For the first time, the classification of small atomic vibrations of molecules and crystals using irreps of their symmetry groups was carried out by Wigner in his classical work [61]. However, this classification is valid only within the framework of the harmonic approximation, i.e., only for the case of conventional NMs. It turns out that each NM corresponds to some irrep of the symmetry group of the system in the equilibrium state. The dimension of such irrep determines the fre-quency degeneracy of the given mode, and its basis vectors determine atomic displacements pattern corresponding to this mode.
Since these results apply only to NMs, the following natural question arises: What will happen with increasing amplitudes of atomic vibrations, when it is necessary to take into account the anharmonic terms in the expansion of the potential energy in a multidimensional power series?
The answer to this question is proposed by the theory of BNNMs, which states that at large amplitudes of atomic vibrations, each BNNM corresponds not to one irrep, but to some set of irreps of the group G 0 . In this sense, the theory of BNNMs is a generalization of Wigner's theory of the classification of conventional NMs. This is the novelty of the theory of BNNMs in comparison with the traditional approach to the study of the dynamics of molecules and crystals.

Dynamical aspect
The geometric properties of bushes, such as the number of bushes of a given dimension and their structure, do not depend on amplitudes of their modes. However, the dynamical properties of bushes depend on these amplitudes in essential way. Let us follow the possible scenarios of dynamical behavior of the system with a discrete symmetry group G 0 with increasing amplitudes of atomic vibrations.
The initial configuration of our system is described by the coordinates of all its atoms. One can transit from such a description to a description of the system in the modal space whose dimension coincides with the system's dimension (the number of NMs is equal to this dimension).
In this case, the new dynamical variables are NM amplitudes. The anharmonic terms that appear in the expansion of the potential energy into power series and describe the interaction of modes with amplitudes x and y can be, for example, x 2 y, x 2 y 2 , x 3 y 2 , etc. Hereafter, we call modes belonging to this bush by term "active modes" (they possess nonzero amplitudes), and all other modes will be called "sleeping modes" (they possess zero amplitudes).
Let us consider the anharmonic term of the form x 2 y. Such a term in potential energy generates the equations of Newton's second law of the form: (the masses of particles is assumed to be unit). It can be seen from the first equation that the first mode (x) can be at rest (x = 0) for any values of the second mode (y). In contrast, the second mode (y) must be nonzero if the first mode has been excited (i.e., if x = 0). In other words, the mode y is involved into the oscillatory process due to the direct (force) interaction with the mode x, but not vice versa, the mode x is independent of the mode y. This is characteristic example of a two-dimensional bush with two modes x and y, which demonstrates the difference between these two modes.
In terms of the theory of BNNMs, the mode x is called the root mode, if only it is excited at the initial instant, and then mode y is called the secondary mode. But in the expansion of the potential energy, there may be other anharmonic terms, for example, the term x 2 y 2 . Then it follows from Newton's equations that the excitation of one of the modes does not necessarily imply the excitation of the second mode, they are independent of each other.
However, for appropriate values of the amplitudes of the considered modes, it may turn out that, for example, the initially excited mode x still leads to the excitation of the previously sleeping mode y due to internal parametric resonance. (We fall into the region of instability of the zero solution of the Mathieu equation for the sleeping mode y). Such processes was considered in detail in our work [62], in which the concept of parametric interaction between NMs was introduced. Now it is important for us that, as a result of internal parametric resonance, some dormant nonlinear mode can "wake up," as a result of which the number of active modes in the system increases and, thus, the dimension of the initial bush increases. At the same time, its symmetry group automatically reduces: It turns into a new bush of higher dimension. But even in this case, the dynamics of the system is described by some bush, and we remain in the framework of the theory of BNNMs.
Above we discussed the first scenario of the behavior of the dynamical system with increasing amplitudes of atomic vibrations, which leads to expansion of the bush, i.e., to increasing its dimension. In particular, the bush can expand up to a general bush, which includes all N NMs of the considered system. But at the same time, we still remain within the framework of the theory of BNNMs.

Expansion of bushes due to awakening previously sleeping modes
Below, as an illustration of this process, we consider the stability diagrams of one-dimensional and twodimensional bushes in the chains of α-FPU and β-FPU types given in our work [62].
Each one-dimensional bush in the FPU chain generates some pattern of atomic displacements (similarly to that shown in Fig. 2 of this paper). In [62], the following patterns of atomic displacements were found for one-dimensional bushes in the α-FPU chain: • |x, −x| This is the π -mode or zone-boundary mode associated with doubling the primitive cell of the one-dimensional crystal. The variable x denotes an arbitrary displacement of an atom from its position in the equilibrium state of the chain. • |x, 0, −x|, |0, x, −x|, |x, −x, 0|. These three bushes are associated with a tripling the size of the primitive cell of the one-dimensional crystal. • |0, x, 0, −x|, |x, 0, −x, 0|. These two bushes are associated with quadrupling the size of the primitive cell.
Each bush corresponds to a certain subgroup G j of the group G 0 of the chain in equilibrium. In our case, G j are subgroups of the dihedral symmetry group D n . We define subgroups G j by the list of their generators. Let us emphasize that there are no other onedimensional bushes in one-dimensional chains with any potential of the interatomic interactions (not only with FPU potentials). As we already explained above, the reason lies in the fact that for all other subgroups G j of the dihedral group D n , the solution of Eq. (B2) depends on 2, 3, 4, etc., arbitrary parameters.
In [62], for all above-listed one-dimensional bushes in the FPU chains, we have presented stability diagrams with respect to their transformation into bushes of higher dimension. Note that Lyapunov's stability of motion is not lost in this case.