Introduction to the theory of bushes of nonlinear normal modes for studying large-amplitude atomic vibrations in systems with discrete symmetry

The research group from the University has been developing the theory of bushes of nonlinear normal modes (NNMs) in Hamiltonian systems with discrete symmetry since the late 90s of the last century. Group-theoretical methods for studying large-amplitude atomic vibrations in molecular and crystal structures were developed. Each bush represents a certain collection of vibrational modes, which do not change in time despite the time evolution of these modes, and the energy of the initial excitation remains trapped in the bush. Any bush is characterized by its symmetry group, which is a subgroup of the system’s symmetry group. The modes contained in the given bush are determined by symmetry-related methods and do not depend on the interatomic interactions in the considered system. The irreducible representations of the point and space groups are essentially used in the theory of the bushes of NNMs, and this theory can be considered as a generalization of the well-known Wigner classification of the small-amplitude vibrations in molecules and crystals for the case of large-amplitudes vibrations. Since using of the irreducible representations of the symmetry groups can be an obstacle to an initial familiarization with the bush theory, in the present review, we explain the basic concepts of this theory only with the aid of the ordinary normal modes, which is well known from the standard textbooks considering the theory of small atomic vibrations in mechanical systems. Our description is based on the example of plane nonlinear atomic vibrations of a simple square molecule.


Introduction
Group-theoretical methods for studying small atomic oscillations of molecules and crystals were applied for the first time by Wigner in his pioneering paper [1] published in 1930. He showed that normal modes that are introduced in the harmonic approximation for systems with discrete symmetry can be classified by irreducible representations (irreps) of the symmetry group of these systems in equilibrium. All modes related to a given irrep possess the same frequency, and eigenvectors of the matrix of the force constants corresponded to it can be considered as basis vectors of this irreducible representation.
The possibility of the above classification is based on the linearity of dynamical equations obtained in the harmonic approximation. A similar classification of quantum states is also possible due to the linearity of the Schrödinger equation (Wigner's theorem).
The theory of bushes of nonlinear normal modes (hereafter, bush theory) was developed in the 90s of the last century for studying atomic vibrations in the systems with discrete symmetry for the case of arbitrary amplitudes. As well as the Wigner theory, it essentially uses the apparatus of irreducible representations of the symmetry groups, however it allows one to take into account interactions between vibrational modes of any amplitude and type.
In the framework of this theory, the concept of bushes of nonlinear normal modes (NNMs) was introduced. They represent some exact dynamical objects in the arbitrary nonlinear systems with discrete symmetry. These dynamical systems can be described by ordinary differential equations of classical mechanics for any interatomic potentials, and also by nonlinear equations of rather different physical and mathematical nature. For example, this is the case of integrodifferential equations that arise in the Hartree-Fock method and in the density functional theory in the framework of the quantum mechanical approach.
All these problems are discussed in the present review. Let us begin with some historical remarks concerning the works of various authors, which in some sense relate to the bush theory.

Interactions of normal modes in the case of weak anharmonism
Since the harmonic approximation can be used only for studying atomic vibrations with very small amplitudes, various authors developed many different approaches for taking into account some anharmonic terms in the system's potential energy for the case of weak nonlinearity. These are different versions of the perturbation theory, as well as the continuation of the solution of linear approximation with respect to small nonlinearity parameters.
Firstly, let us refer to the works of Lyapunov [2], who proved the theorem that each normal mode ϕ can be continued with respect to a nonlinear parameter μ in such a way that a new periodic mode with the frequency ω(μ) close to that of the mode ϕ is obtained. This theorem was proven under certain conditions, which we do not specify here, while the coefficient before all terms of the Taylor series of the potential energy above harmonic approximation can be used as parameter μ. The number of these Lyapunov's nonlinear normal modes, as well as the number of conventional normal modes, is equal to the dimension of the considered dynamical system. They are constructed with the aid of a cumbersome procedure in the form of a power series with respect to the nonlinear parameter μ. However, this procedure can usually be carried out only up to sufficiently small values of the parameter μ, which is the reason that Lyapunov's modes are not widely used in solving physical problems.
Continuation of the solutions of nonlinear differential equations in the form of normal modes with respect to some nonlinear parameters has been used in many works by different authors [3]. For example, let us refer to the Lindstedt-Poincare method, where the classification of the terms of differential equations into resonant and non-resonant type is used. The former lead to a shift of the oscillation frequency to eliminate secular terms in the solution, while the latter play only the role of small perturbations. We would like to note that this method was used in the paper [4] by Flach, Ivanchenko and Kanakov to introduce the concept of q-breathers, which describe energy localization in the crystal reciprocal lattice.
Various methods of perturbation theory are also widely used in the case of weak nonlinearity of the considered system. In this regard, it is worth to mention the interactions of different types of quasiparticles in crystal physics, in particular, interactions between phonons, which are obtained as a result of quantization of normal vibrational modes [5].

Nonlinear normal modes in Hamilton systems
In the framework of advanced mathematics, nonlinear normal modes as certain time-periodic dynamical objects were considered in [6]. These works address to mathematicians dealing with non-linear systems of classical dynamics and, as a consequence, they exploit modern and quite complex mathematical techniques such as symplectic geometry, topology, singularity theory, etc.
Here, we would like to comment on the difference of these works from our approach based on the concept of bushes of nonlinear normal modes (a more detailed comparison has been presented in [7]).
First of all, it should be emphasized that our theory is applicable not only for the Hamiltonian systems only, but for a much wider class of nonlinear equations, including Onsager type equations used in the thermodynamical description of physical systems, and even nonlinear integro-differential equations arising at solving the multielectron Schrödinger equation in the framework of the density functional theory.
Secondly, bush theory can describe not only periodic regimes, but also quasiperiodic ones up to those with dynamical chaos.
Thirdly, the basic ideas of the bush theory are rather simple for understanding not only by theorists, but also by experimenters, and we hope to demonstrate this in the main text of the present review.

Rosenberg nonIinear normal modes
There are several different interpretations of the term "nonlinear normal mode" (NNM). Mentioning these dynamical objects in the present paper, we will always bear in mind the nonlinear normal modes by Rosenberg, which were introduced in [8,9].
The main property of the Rosenberg mode in the N-particle mass point mechanical system is that the evolution of all its degrees of freedom is described by the same time function. As a consequence, the ratio of the displacements of all particles to the displacement of any selected particle at any moment is a constant value. Conventional (linear) normal modes also satisfy this definition, and in this sense, Rosenberg modes are their generalization.
A detailed description of the properties of Rosenberg modes can be found in [10,11]. Rosenberg modes are widely used in studying different mechanical systems [12,13].
The value of the Rosenberg NNMs in comparison with linear normal modes is due to the fact that these modes can describe vibrations with large amplitudes. However, Rosenberg modes are rarely used in solving physical problems, because they can exist only under very severe restrictions on the interparticle interactions. Several classes of such interactions were found, the most relevant of which is the case when potential energy of the system is a homogeneous function of some order of all its variables. For physical problems, this is a very rare case.
However, it was shown [14] in the framework of bush theory that there are certain group-theoretical reasons for the existence of a certain number of Rosenberg modes in dynamical systems with discrete symmetry even in the case of arbitrary interatomic interactions. Rosenberg modes turn out to be one-dimensional bushes and they can be found with the same group-theoretical methods, as well as bushes of higher dimensions. Let us note that Rosenberg modes describe periodic regimes, while multidimensional bushes describe quasiperiodic motions.
Existence of Rosenberg modes in systems with discrete symmetry is determined by the symmetry-related causes, and we called them "symmetry determined" Rosenberg modes. Below we consider only such modes without this specification.
We would like to comment on the paper [15] by Mishra and Singh in connection with the concept of Rosenberg modes, on the one hand, and with the apparatus of irreducible representations of symmetry groups on the other hand. The authors of this paper consider oscillations in the mechanical system of weights connected by massless nonlinear springs, which is characterized by symmetry group C 4ν .
In order to find Rosenberg modes, they pass into the modal space associated with the basis vectors of all irreducible representations (irreps) of the system's symmetry group, as it is done in the Wigner theory of oscillations with small amplitudes, and then consider the potential energy in the form of homogeneous function of variables, corresponding to individual irreps separately from each other. This approach should be recognized as erroneous since the authors take into account in the decomposition of potential energy only invariants constructed from variables transforming according to individual irreps completely ignoring all mixed invariants constructed from variables belonging to the different irreps. However, mixed invariants usually play a very significant role (see below a discussion of the complete condensate of order parameters upon phase transitions in crystals).
On the other hand, taking into account all possible mixed invariants leads to the dynamical aspects of the bush theory. However, in the present review, we discuss only the geometrical aspects of this theory, which can be explained in an understandable language and in an almost elementary way.

The theory of the complete condensate of the primary and secondary order parameters upon structural phase transitions in crystals
The apparatus of irreducible representations (irreps) of the space groups is an important component of the Landau theory of the second-order phase transitions in crystals [16]. In this theory, thermodynamic potential is decomposed into series of polynomial invariants constructed from variables that are transformed according to all irreducible representations of the space group of the highly symmetric crystal phase (in the bush theory, we call it by the term "parent group"). The coefficients of this decomposition depend usually on two thermodynamic variables, temperature T and pressure p, whose change can cause a structural phase transition with decreasing of crystal symmetry from the parent group G 0 to some of its subgroup G j . At the point of the phase transition, the coefficient c(p, T) in front of the quadratic invariant of some irrep Г i of the parent group vanishes. This so-called "critical" irrep is responsible for the occurrence of the phase transition, and the set of its variables (in the general case, it is multidimensional) forms the order parameter (OP) of the considered phase transition. The appearance ("condensation") of this primary order parameter, leads to the condensation of a number of "secondary" parameters corresponding to some other irreps, due to nonlinear interactions in the crystal. These secondary order parameters are induced by mixed invariants containing variables of different irreps. Near the point of phase transition in the p-T plane, secondary order parameters are smaller in magnitude than the primary one, however as we move away from this point the role of secondary parameters increases and they can be equal in order of magnitude with the primary parameters (in the general case, all order parameters are vector quantities!).
In our works [17,18], group-theoretical algorithms were developed for constructing the "complete condensate" of all primary and secondary order parameters upon phase transitions in crystals with any space parent group. These algorithms formed the basis for development of computer software for studying structural phase transitions in crystals [19].
The bush theory represents a generalization of the static theory of the complete condensate of order parameters to the case of atomic dynamics in systems with discrete symmetries [7,14,20].
At the same time, when we worked on the theory of the complete condensate of order parameters, group-theoretical methods for studying structural phase transitions in crystals were also intensively developed at Brigham Young University (USA, Utah) by Hatch and Stokes. They developed the powerful software "Isotropy" [21], which allows one to make many group-theoretical calculations using computer.
Several group-theoretical works such as listing all possible phase transitions in structures with layer groups were performed by our and American groups almost simultaneously, which allowed us to establish a good cooperation and to carry out several joint works [22 -24]. The most ambitious was the work [23], in which all possible symmetry-determined Rosenberg nonlinear normal modes were found for all possible mechanical structures described by all 230 space symmetry groups. This work was based on the idea of "irreducible bushes of vibrational modes".

Fermi-Pasta-Ulam chains
In the scientific literature, there is a large body of works devoted to studying interactions between vibrational modes in the Fermi-Pasta-Ulam (FPU) chains.
These chains were introduced on the initiative of E. Fermi in the classical work [25] to study the possibility of thermalization of the system due to the weak nonlinear interactions between conventional normal modes. It was expected that, at sufficiently large times, the system will transfer to the state of the energy equipartition between all degrees of freedom, which are these normal modes.
To the authors' surprise, the energy from the initially excited mode with the longest wavelength has been passed to a very small number of the modes following the first one, but was not transmitted to any distant shortwave modes. Moreover, computational experiments allowed Fermi, Pasta and Ulam to reveal the so-called "return phenomenon", which consists in the fact that over time the energy from other modes that appeared earlier is again transferred almost completely to the first mode, and the whole process is repeated in a quasiperiodic manner without any trend to the energy equipartition between all vibrational modes. If we translate this phenomenon into the range of sounds heard by a human ear, then the FPU chains generated a very definite melody instead of the expected cacophony of sounds.
The nature of this phenomenon was clarified by Zabusky and Kruskal [26], who made the transition from the ODE system, describing the vibrations of the FPU chain, to the partial differential equation, which turned out to be the Korteveg -de Vries (KdV) equation, and introduced the concept of solitons. The return phenomenon in the FPUchains was explained in some a sense in this way.
FPU chains turned out to be very convenient objects for studying various types of nonlinear dynamical objects. Until now, more and more new studies of these seemingly simple systems appear in the world scientific literature [27,28].
Poggi and Ruffo in [29] revealed the possibility to excite in the FPU-β chain some sets of normal modes, which can exist without transferring the excitation to other modes. These sets are examples of bushes of vibrational modes in this one-dimensional system. However, the results of Ref. [29] were obtained by analyzing the specific type of interactions in FPU-β chain without application of any symmetryrelated methods. Such approach is rather cumbersome even for the one-dimensional FPU chain and does not allow generalization to the case of dynamical systems of a more complex structure.
Similar work on the dynamics of FPU chains were carried out later by Shinohara [30].
In Refs. [31] and [32], we presented a detailed comparison of the results of the above mentioned papers with those of the bush theory and demonstrated how much easier and more efficient can be group-theoretical methods for studying intermode interactions even in the case of one-dimensional chains, not to mention two-dimensional and threedimensional periodic structures.

Invariant manifolds and bushes of NNMs
Another approach was developed by Manevich [33], Rink and Verhulst [34 -36], who used idea of finding symmetrydetermined invariant manifolds for dynamical systems with discrete symmetry. The bush theory also begins with the idea of finding such invariant manifolds, but it goes much further.
Indeed, in the above cited papers there is no one essential element of the bush theory -decomposition of the invariant manifolds into nonlinear normal modes with the aid of the apparatus of irreducible representations of the symmetry groups. On the other hand, it is this element allows one to introduce the concept of the bush of modes as a certain fundamental dynamical object that has its own internal dynamics, since the collection of modes of a given bush preserves over time, while their amplitudes undergo temporary evolution.
The bush theory does not cease to exist when the stability of a given bush is lost because in this case it transforms into another bush of a larger dimension.
Finally, it is the using of irreducible representations of the symmetry groups that allows one to search for bushes of modes in systems with discrete symmetry of any complexity.

Bushes of vibrational modes and discrete breathers
Another turn in the history of the bush theory was associated with the theory of discrete breathers. This direction of activity is connected with the idea put forward by Dmitriev, who suggests to construct a new type of discrete breathers in crystals by applying to the vibrational bushes some "localizing functions". A series of joint works of the Ufa and Rostov research groups in the field of nonlinear dynamics was published using this idea [37,38].

Bushes of NNMs and the density functional theory
The bush theory is based on the group-theoretical methods that are applicable not only to dynamical systems with discrete symmetry described by classical differential equations, but also to a wide class of other nonlinear equations. Especially significant are the integro-differential equations of the theory of the density functional [39], since this theory enables to take into account quantum-mechanical variation of the electron shells of atoms during classical motion of their nuclei (see Sec. 7 of this review). We have published several papers on this subject [40 -43].

Invariant manifolds
We begin with studying the classical dynamics of N material points whose interactions are described by some phenomenological potentials (model 1). This is a standard but quite rough approximation for modeling dynamics of atomic systems.
In the study of dynamical systems, the concept of invariant manifolds (IMs) plays an important role. Let us explain this concept.
The dynamics of a system of N material points is described by a trajectory in the 6N-dimensional phase space of all their coordinates and momenta (velocities). By the definition, if the system's initial state, determined by the coordinates and velocities of all its particles, lies on a given IM, then the entire phase trajectory will lie on this manifold, i. e. the system does not leave this manifold at any time during its evolution. Obviously, any exact solution of the dynamical equations must lie on a certain IM, and this fact determines the importance of finding such manifolds.
Unfortunately, in mathematics, there are no general methods for searching invariant manifolds for arbitrary dynamical systems. On the other hand, if the system has some discrete symmetry, then one can propose general methods for constructing such manifolds.
The key idea is that invariant manifolds correspond to the subgroups of the system's symmetry group. The theory of bushes of nonlinear normal modes starts with this idea [7,14,20].
In Fig. 1, we show a planar square molecule consisting of four atoms located at the corners of the square. We consider this molecule in space and, therefore, its point symmetry group is C 4ν in the Schoenflis notation. This group consists of 8 elements. They are: -vertical fourth-order axis, which determines rotations by 0, 90, 180 and 270 degrees around the Z axis; -four reflection planes that pass through this vertical axis. Two of them are "coordinate planes" -they pass through the axes X and Y, respectively, while the other two reflection planes are "diagonal" -they pass through the diagonals of the square.

Harmonic approximation and normal modes
The standard approach for studying atomic vibrations with small amplitudes in molecules and crystals is the harmonic approximation [44]. In this approximation, the potential energy of the nonlinear system is expanded into the multidimensional Taylor series over all displacements of atoms from their equilibrium positions and all members of this series, starting with cubic ones, are discarded. As a result, one obtains a quadratic form. Then the classical Newton equations, determined by such quadratic potential energy, form a system of linear differential equations with constant coefficients with respect to the atomic displacements (x j ).
This system can be split into independent equations if the above quadratic form is reduced to the canonical form, in which it takes the form of a superposition of only the squares of the variables x i 2 (all terms of the type x i x k are eliminated). This can be achieved by using a linear orthogonal transformation of the force constant matrix (the matrix of the second partial derivatives of the potential energy with respect to all variables) to the diagonal form. As it is well known, this problem is reduced to finding all the eigenvalues and eigenvectors of the force constant matrix by an appropriate orthogonal transformation in the space of old variables.
The new dynamical variables y j (t), in contrast to the old variables x i (t), are normal modes of the original dynamical system. These modes are independent of each other in the harmonic approximation, and their number is equal to the number of the system's degrees of freedom.
If some weak nonlinear terms in potential energy are taken into account, the normal modes begin to interact with each other, and their interactions can be calculated in the framework of the perturbation theory.

Normal modes for the square molecule
Let us consider all normal modes for the case of atomic vibrations in the plane of the square molecule presented in Fig. 1. This molecule possesses 8 degrees of freedom because each atom has two degrees of freedom (displacements along the X and Y axes) and, therefore, 8 normal modes. A certain pattern of atomic instantaneous displacements, which are shown by arrows in Fig. 2, corresponds to each of the normal modes.
All these normal modes were obtained for the case when atoms interact via the Lennard-Jones potential in "standard form" (both phenomenological constants corresponding to the contribution from the attraction energy and repulsion energy are set equal to unity). The eigenvalues of the force constant matrix, i. e. squares of the normal modes frequencies are λ 1 = 32.575, λ 2 = λ 7 = λ 8 = 0, λ 3 = 34.541, λ 4 = −1.965, λ 5 = λ 6 = 33.986. In Fig. 2, modes φ 7 and φ 8 are not indicated (they describe the motion of the molecule as a whole along the X and Y axes). All the atomic displacement patterns were obtained from the eigenvectors of the force constant matrix.
Let us note that the form of the atomic patterns is completely independent of the interparticle interactions, while vibrational frequencies depend substantially on these interactions.
We see that the atomic patterns for all normal modes have some symmetry, and the corresponding point groups in the Schoenflies notation are also shown in Fig. 2. Let us consider the form of these patterns in more detail.
First of all, it should be emphasized that the lengths of all arrows, showing atomic displacements, are of the same length for any chosen mode! Fig. 2 a shows the instantaneous deformation of molecule corresponding to the so-called "breathing mode" φ 1 : molecule retains its square shape, which periodically increases and decreases in size. The mode in Fig. 2 b describes the rotation of our square molecule as a whole around its center according to mode φ 2 . Fig. 2 c shows the rectangular deformation of the molecule during its oscillations according to mode φ 3 : it periodically stretches along one coordinate axis, contracting along the second axis, and vice versa. Fig. 2 d shows the rhombic deformation of the molecule during its oscillations in mode φ 4 -it either stretches along one diagonal, contracting along the second diagonal, and vice versa. Finally, the trapezoidal deformation of our molecule with a different orientation along the coordinate axes is depicted in Figs. 2 e, g. It corresponds to the modes φ 5 and φ 6 .
The above atomic patterns were obtained as eigenvectors of the force constant matrix for the case of the Lennard-Jones potential. However, it is important to emphasize that the same types of patterns can be obtained with the aid of basis vectors of all irreducible representations of the molecule's symmetry group without taking into account any specific interatomic potentials [1].
It is noteworthy that eigenvalue, corresponding to the rhombic normal mode φ 4 , turns out to be a negative number and, therefore, the frequency of "oscillations" of this mode is imaginary. This means the instability of the square molecule relative to the rhombic deformation, i. e. the molecule does not oscillate with respect to the square equilibrium positions but leaves them.
Starting the study of normal modes in a square molecule, we implicitly assumed that its square configuration is stable and small atomic vibrations occur relative to this configuration. However, the results of our calculations for the case of the Lennard-Jones potential show that the square configuration of the molecule is unstable. On the other hand, it can be stable in the case of other interatomic interactions, or when one more atom is placed in the center of the square configuration. Let us note, in this regard, that centered plane square molecules do exist in nature (for example, the molecule XeF 4 ).
Our group-theoretical analysis does not depend on the centering of the square, since we exclude modes corresponding to the movement of the molecule along coordinate axes, assuming thereby the immobility of the molecule center where the additional atom can be placed.

Invariant manifolds and bushes of nonlinear normal modes for the square molecule
In the framework of classical dynamics, which is described by Newton's differential equations, the dynamical behavior of our molecule is uniquely determined by its initial state. This is a situation of so-called classical determinism -"the past completely determines the future". In such a situation, we can say that there exists a certain "law of conservation of symmetry". This means that if some element of symmetry of the system's motion exists at the initial moment (it is determined by the atomic displacement pattern), then this element exists during the whole time of movement until it is stable. Let us consider the atomic displacement pattern corresponding to the mode φ 4 , which represents a rhombus with diagonals of the same length since, as already was mentioned, all arrows showing displacements of the atoms from the corners of the square, possess the same length in the framework of the harmonic approximation. The symmetry group of this pattern is C 2ν d . However, the considered pattern is not one of the most general form, which corresponds to the point group C 2ν d ! Obviously, the rhombus with diagonals of different lengths also has the same symmetry group and it corresponds to displacements of different lengths for atoms located on different diagonals.
Thus, the invariant manifold corresponding to the point group C 2ν d is wider than that which we obtained for small oscillations in the harmonic approximation. From this point, the introduction to the bushes of vibrational modes starts.

Bushes of vibrational modes
Displacements of atoms from their equilibrium positions in the plane of the molecule can be described by a configuration vector of the form: where the first two coordinates are the x-and y-displacements of the first atom, the next two coordinates are x-and ydisplacements of the second atom, etc. (the atom's numbers are shown in Fig. 1). In these notations, the eigenvector of the force constant matrix, corresponding to the φ 4 mode and describing the rhombus of its atomic pattern, has the following form: On the other hand, the coordinate vector corresponding to the rhombic pattern of the general form with the symmetry group C 2ν d is This vector determines the rhombus with diagonals of different lengths. If we decompose it over the complete set of eigenvectors of the force constant matrix (note that by virtue of Wigner's results [1] we can look at them as basis vectors of irreducible representations of the symmetry group C 4ν ), then it is easy to verify that there are only two nonzero terms in such decomposition and they correspond to the modes φ 4 and φ 1 : If we choose initial conditions for solving nonlinear dynamical equations of our molecule in the form of (2) for some fixed values of a, b [or A, B in (3)], we find that these constants begin to change during the numerical integration of these equations, i. e. they are certain functions of time, a(t) and b(t). The configuration vector X 4 (t) will possess the form (2) for any instant t, i. e. it does not leave the invariant manifold determined by Eqs. (2) or (3): Thus, we have obtained two-dimensional bush B 4 with the symmetry group C 2ν d .
On the other hand, for vibrations corresponding to the symmetry group C 4ν , starting from the normal mode φ 1 , we will obtain the dynamical regime of the form This is the one-dimensional bush B 1 with C 4ν point group.
Using the above examples, we can introduce some notions and present some ideas of the theory of bushes of nonlinear normal modes.
-We have excited bush B 4 with the aid of excitation of normal mode φ 4 . If its amplitude is sufficiently small, we will see during integrations of nonlinear equations practically only this mode in full accordance with the theory of small atomic vibrations. However, as the initial amplitude of φ 4 (t) increases, the additional mode φ 1 (t) can be observed in the numerical solution of nonlinear equations since nonlinear terms begin to play a more and more important role. We call the primarily excited mode φ 4 (t) the "root mode" of the bush, while the mode φ 1 (t), involving into the vibration process because of its interaction with the root mode, we call the "secondary mode".
-It is important to note that our secondary mode φ 1 (t) possesses a higher symmetry group (C 4ν ) with respect to the group C 2ν of the root mode. It can be proved in the framework of the bush theory that this is a general property -symmetry of all secondary modes always is equal or higher than that of the root mode. Note that the above-discussed notions are similar to primary and secondary order parameters in the theory of phase transitions in crystals [17].
-One should not think that secondary modes always are smaller in amplitude than the root mode -they can be of the same order of magnitude for the case of large nonlinearity.
-A given bush can be certainly excited by simultaneous excitations of all its modes, not only by one root mode.
-One-dimensional bush represents a periodic motion, while the m-dimensional bush represents a quasi-periodic motion with m main frequencies (and their different integer linear combinations) in its Fourier spectrum.
Let us now consider bush B 2 generated by the root mode φ 2 (t), that describes the rotation of our square molecule as a whole (see its atomic pattern in Fig. 2 a). Being conventional normal mode, it is independent of all other modes in the harmonic approximation. But if we take into account nonlinear terms in the decomposition of the molecule Lennard-Jones potential energy, the mode φ 2 (t) involves into the vibrational process one more mode φ 1 (t) (and only this mode!). This is the same secondary mode as that in the above-considered bush B 4 . As a result, the connection between rotation and oscillation movements appears, i. e. our molecule rotates and simultaneously oscillates. The connection between these two types of motion is a well-known property of real molecules.
The complete list of all bushes, which can be excited in the square molecule, is presented in the Table 1. Each of them is determined by a certain subgroup of the molecule symmetry group C 4ν . Therefore, to find all bushes we must consider all such subgroups, taking into account their different settings into the parent group C 4ν . Each subgroup determines the certain invariant manifold, which then we must decompose into the basis vectors of the irreducible representations of the parent group to single out all modes forming the bush. It is obvious from the described procedure that bushes found in this way fully independent of any specific interatomic interactions and these group-theoretical results are exact. This is a geometrical aspect of the bush theory.
Each bush in Table 1 is presented as a linear combination of certain modes. The second column shows the symbols of irreducible representations (irreps) of the symmetry group C 4ν , which contribute to a given bush. There are five irreps in this group -four one-dimensional representations (Γ 1 ; Γ 2 ; Γ 3 ; Γ 4 ) and one two-dimensional representation (Γ 5 ). At this stage of the description of the bush theory, this column can be ignored -we will return to the discussion of connection of bushes with irreducible representations of symmetry groups later.
Above, we have briefly described some of the geometrical aspects of the bush theory, and now we would like to discuss several important issues related to their dynamics.

Rosenberg nonlinear normal modes in the framework of the bush theory
Now we should like to explain, why we use the term "bushes of nonlinear normal modes" (NNMs).
In 1962 Rosenberg introduced the concept of the "similar nonlinear normal mode" [8]. By definition, this is a dynamical regime in N-particle mechanical system for which the evolution of all degrees of freedom is described by one and the same function of time f(t): where a i (i =1..N) are constant coefficients determining the amplitudes of vibrations of all the degrees of freedom. Conventional (linear) normal modes also satisfy the definition (6) with f (t) = sin(ωt + φ 0 ). In the general case, the function f (t) for the Rosenberg mode is determined by a certain differential equation (the so-called "governing" equation).
It is very important to note that Rosenberg NNMs can exist only in the systems with very specific properties, for example in those whose potential energy represents a homogeneous function of all its variables (obviously, this is the very exotic case).
However, one can prove in the framework of the bush theory that in the systems with discrete symmetry there can exist Rosenberg NNMs fully independent of the type of interatomic interactions. The existence of these modes is dictated only by symmetry-related properties.
Comparison of Eq. (5) and Eq. (6) shows that any onedimensional bush is a Rosenberg mode. On the other hand, we can say that the two-dimensional bush B 4 , described by Eq. (3), contains two different virtual Rosenberg modes -the root mode and the secondary mode Why do we call these modes of the bush B 4 by the term virtual Rosenberg modes? The matter of the fact is that the mode (7) cannot exist without the mode (8) (the last mode automatically involves into the vibrational process by the former mode), while the mode (7) can exist as an independent Rosenberg mode only in the case when this mode will be excited at the initial time without excitation of any other vibrational modes.
Since any vibrational regime in a system with the symmetry group G 0 can be associated with a certain subgroup of this group and some bush corresponds to it, we can refer to any vibrational regime in the dynamical system with discrete symmetry as to the bush of nonlinear normal modes.

Parent symmetry group
The parent group is the group characterizing the symmetry of the considered system. Different physical ideas may be used for choosing the parent group. Wigner in his classical paper on classification of normal modes in systems with discrete symmetry [1] chose as parent group the symmetry group of the system's equilibrium state. However, most often for this role the symmetry group of the system's Hamiltonian is chosen. The symmetry group of Hamiltonian can be higher than that of the equilibrium state. Indeed, this group can dictate the existence of several symmetry-equivalent equilibrium states that transform into each other due to the action of some of its elements. For example, this is a typical situation in perovskite BaTiO 3 , where around titanium atom Number of the bush Irreps which contribute to a given bush Bush Symmetry group G⸦G 0 1 Γ 1 μ 1 φ 1 C 4ν 2 Γ 1 ; Γ 2 μ 1 φ 1 + μ 2 φ 2 C 4 3 Γ 1 ; Γ 3 μ 1 φ 1 + μ 3 φ 3 C 2ν c 4 Γ 1 ; Γ 4 μ 1 φ 1 + μ 4 φ 4 C 2ν d 5 Γ 1 ; Γ 3 ; Γ 5 μ 1 φ 1 + μ 3 φ 3 + μ 5 φ 5 C s c 6 Γ 1 ; Γ 4 ; Γ 5 μ 1 φ 1 + μ 4 φ 4 + μ 5 (φ 5 + φ 6 ) C s d 7 Γ 1 ; Γ 2 ; Γ 3 ; Γ 4 μ 1 φ 1 + μ 2 φ 2 + μ 3 φ 3 + μ 4 φ 4 C 2 8 Γ 1 ; Γ 2 ; Γ 3 ; Γ 4 ; Γ 5 μ 1 φ 1 + μ 2 φ 2 + μ 3 φ 3 + μ 4 φ 4 + μ 5 φ 5 + μ 6 φ 6 C 1 at the center of the cubic cell there are several symmetrically located local energy minima, each of which has less symmetry compared to their complete ensemble. During the phase transition of the displacement type, the titanium atom passes into one of these minimums, as a result of which the symmetry of the system decreases. In the case of the ordering type transition, the titanium atom can be found at different minima with different probabilities. If certain subgroups of the initially chosen parent group were obtained and then we find that the system has some additional symmetry elements, i. e. the parent group is higher, then the set of the initially found bushes is changed because of combining some old bushes into one larger bush with a higher own symmetry, as well as of appearing some new bushes. For example, in arbitrary monoatomic chains, there can exist only three one-dimensional vibrational bushes, but if the potential of interatomic interactions is symmetric with respect to each site of the chain, then there can exist already five such bushes, depending on the number of atoms in the chain [31,32,45,46].
To check, which choice of the parent group is the most appropriate for a given system, one can perform its modeling using methods of molecular dynamics or density functional theory methods. Based on the results of such modeling, it is necessary to reveal which secondary modes and with what amplitudes are involved into the vibrational process when the root mode of the bush is excited.

Stability of the bushes of nonlinear normal modes
As was already discussed, the phase trajectory of the considered system does not leave a given symmetrydetermined invariant manifold during its time evolution. This conclusion is based on the classical determinism provided by the uniqueness of the solution of Newton equations for the given initial conditions. If this solution does not practically change by any sufficiently small perturbations, which always present in any physical system, one can speak about the case of stability by Lyapunov [2]. However, in the case when such stability is lost, even infinitely small fluctuations can lead to fully unpredicted results.
As we have already seen, the root and secondary modes, containing in a given bush, are active modes. All other vibrational modes are "sleeping" with respect to this bush since they possess zero amplitudes. When the amplitude of the root mode increases, its frequency changes (this is the main property of nonlinear oscillations!). At some value of the root mode amplitude, the so-called "critical" amplitude, a parametric resonance can occur between the root mode frequency and the frequency of one of the sleeping modes.
As a result, this sleeping mode wakes up, i. e. becomes active. Thus, the number of active modes in the system increases, and a new bush of a larger dimension and a lower symmetry arises.

Bushes of nonlinear normal modes in the models based on the density functional theory
Our second model is the model based on the density functional theory (DFT). Quantum-mechanical models in the framework of this theory are much more adequate to reality [39]. Indeed, the typical accuracy of such models in atomic coordinates is about 1 % and about 10 % in binding energy. They take into account quantum characteristics of atoms and additional degrees of freedom corresponding to the electrons of atomic shells. The well-known software packages for treating DFT-models, such as ABINIT, Quantum ESPRESSO, VASP, allow one to consider dynamical problems in the Born-Oppenheimer approximation, and to take into account the influence of atomic-shell polarization on the classical motion of nuclei.
Since the density functional theory is based on the Schrödinger multielectron equation, which is a linear equation, one can ask what new information above that of the Wigner's classification of atomic vibrations by individual irreducible representations of the system's symmetry group can be obtained from the theory of the bushes of nonlinear normal modes?
The matter is that despite of the linearity of the exact Schrödinger equation, the well-known approximate methods for solving this equation, such as Hatree-Fock method and Kohn-Sham method of the density functional theory are based on the solving of some systems of nonlinear integrodifferential equations.
Therefore, all group-theoretical methods of the bush theory can be used for studying large-amplitude atomic vibrations in molecular and crystal models in the framework of DFT-theory. Such studies were performed in [40 -43].

Some mathematics
The key idea of the bush theory is very simple -we must find all symmetry-determined invariant manifolds from the condition of the configuration vector invariance with respect to each subgroup of the system's parent symmetry group. However, the technical realization of this idea is rather complicated, and below we very briefly and schematically outline the main points of the group-theoretical approach for studying bushes of nonlinear modes.
At the end of the 70s of the last century, we developed a set of computer programs for the group-theoretical analysis of the complete condensate of the primary and secondary order parameters corresponding to the structural phase transitions in crystals [19]. Later, this complex was generalized for studying bushes of nonlinear normal modes. Below, we consider some details of the group-theoretical methods of the bush theory.
1. Let us begin with the case when we already know a certain subgroup G j of the parent group G 0 . The invariant manifold, corresponding to this subgroup, can be found from the condition of configuration vector invariance: Here X j is the invariant vector, while Ĝ j is the group of operators, isomorphic to the symmetry group G j , which act in multidimensional vector space of all N degrees of freedom (note that elements of the group G j act in ordinary threedimensional space).
On the other hand, the vector X j can be decomposed into basis vectors φ i of all irreducible representations (irreps) Γ i of the parent space group G 0 : The invariance condition (9) corresponds to the whole N-dimensional space. The similar invariance conditions for all the individual irreps Γ i can be obtained from it: Here Γ i ↓G j is the so-called restriction of the irrep Γ i on the subgroup G j , which is the set of matrices of the representation Γ i corresponding to all elements of the group G j . For every n i -dimensional irrep Γ i , Eq. (11) represents a certain homogeneous system of n i linear algebraic equations whose general solution depends on several arbitrary constants (a, b, c, etc.). If this solution is zero, irrep Γ i does not contribute to the bush with the symmetry group G j .
The vector C ji being multiplied by the basis vectors of the irrep Γ i (see about these vectors below) determines a concrete form of the vibrational mode φ j associated with this irrep, i. e. the pattern of atomic displacements corresponding to the obtained mode.
According to Eq. (11), the vector C ji is conserved by all matrices of the irrep Γ i , which correspond to elements of the group G j . However, some other matrices of this irrep, corresponding to elements of the parent group G 0 , which are not contained in its subgroup G j , can also conserve the vector C ji . Therefore, this vector can select some additional symmetry elements of the group G 0 (except those that are already included in the subgroup G j ), i. e. it corresponds to a subgroup G j of a higher order than G j .
The subgroup G j of minimal symmetry represents the symmetry group of the whole bush B j and, at the same time, the symmetry group of its root mode (we omit here some subtleties of the discussed problem).
Eq. (11) is the source of the selection rules for excitation transfer from the root mode to the secondary ones.
At this point, it is appropriate to explain the words "the mode belonging to a given multidimensional irrep". Let us obtained a certain invariant vector C ji as the general solution of Eq. (11) for a six-dimensional irrep Γ i , which has the form [a, b, −a, c, −b, c] (six-dimensional irreps frequently occur in the cubic space groups and a few dozens of different invariant vectors with different own symmetry groups G k can be associated with each of these irreps). Then the vibrational mode, corresponding to this vector, is where φ k (k =1..6) are basis vectors of the representation Γ i .
The complete set of the bush modes with symmetry group G j is determined by the collection of all nonzero solutions of the Eq. (11).
2. In order to construct an explicit form of atomic displacement patterns, it is necessary to obtain basis vectors of those irreps of the parent group G 0 , which are contained in the decomposition of the reducible mechanical representation of this group into irreducible representations. Usually, such basis vectors are obtained by the projection operators technique. In our paper [19], a more transparent and more convenient for our purposes method was presented. This is the so-called "direct method" based on the definition of the group representation.
If the invariant vector, corresponding to the subgroup G j , and basis vectors of the irrep Γ i are already found, then the contribution of this representation to the atomic displacement pattern can be obtained by the Eq. (10).
3. A very difficult problem of the bush theory is finding of all subgroups G j of the parent space group G 0 , taking into account all their possible settings in the group G 0 . For this purpose, we developed a special group-theoretic algorithm based on the sequential finding of all possible pair-wise intersections of invariant vectors of a given representation (each of them represents some set in the space of the given irrep). We used this algorithm earlier to construct a complete condensate of order parameters in the theory of structural phase transitions in crystals [17]. 4. A separate important and cumbersome problem is obtaining the so-called "full" irreducible representations of the space groups. For this purpose, we used in our computer program the standard algorithm based on the method of projective irreducible representations, which is well described in [47]. Let us note that it is very difficult to find the full irreps of the space groups using tables from the book by Kovalev [48], due to the deep degree of embedding notations of different levels in each other. For solving specific problems associated with the construction of bushes of modes for structures with space symmetry, we use our own tables in the so-called "genesis" form (for many space groups, we published them earlier as VINITI deposited manuscripts).

Conclusion
As was noted in Introduction, Wigner published his pioneering work on the application of group-theoretical methods for studying and classification of small oscillations of systems with discrete symmetry [1] in 1930. This approach has become classical and has been included in all textbooks on the physics of molecules and crystals. Wigner showed that small vibrations within the framework of the harmonic approximation can be classified by irreducible representations of the symmetry groups of these objects. Indeed, the degeneracy degree of a given fundamental frequency is equal to the dimension of a certain irreducible representation of the symmetry group, while the corresponding normal modes are "transformed according to this representation", i. e. they can be considered as basis vectors of this irrep. Wigner's theorem [49] determines the general form of the Hamiltonian matrix corresponding to small oscillations of the system with discrete symmetry.
We emphasize once again that all these results relate only to small atomic oscillations in the framework of the harmonic approximation. In contrast, the bush theory presents a grouptheoretical classification and corresponding methods for studying dynamical regimes in the systems with discrete symmetry for any amplitudes of the atomic vibrations.
Each bush contains some modes of different irreps of the parent group of the considered systems. In Table 1, we indicated the numbers of those irreps of the group C 4ν that contribute to the bushes of the square molecule. The apparatus of the irreps of space groups, which is used in the bush theory, can be found in [19].
Let us note that the bush theory includes not only a geometrical but also a dynamical part, which the authors hope to consider later in the second part of this review.
Group-theoretical methods were used by our group not only for studying deterministic motion and delocalized modes. For example, such methods were used for studying dynamical chaos in three-dimensional systems with quadratic nonlinearities in [50] and they were applied for searching localized states in the form of discrete breathers (outside of the above-described approach based on the idea of "localizing functions") in [51]. Many similar problems were addressed in reviews [52] and [53].