Nonlinear optical vibrations of single-walled carbon nanotubes. 1. Energy exchange and localization of low-frequency oscillations.

We present the results of analytical study and molecular dynamics simulation of low energy nonlinear non-stationary dynamics of single-walled carbon nanotubes (CNTs). New phenomena of intense energy exchange between different parts of CNT and weak energy localization in the excited part of CNT are analytically predicted in the framework of the continuum shell theory. Their origin is clariﬁed by means of the concept of Limiting Phase Trajectory, and the analytical results are conﬁrmed by the molecular dynamics simulation of simply supported CNTs.


Introduction
The study of the nanoscale objects is extremely exciting not only from theoretical viewpoint, that allows to extend our knowledge of the fundamental principles, but also as a basis of development of materials and devices.New effects and properties, which may appear at the nanoscale, permit to elaborate the materials and devices with properties, which are unreachable at the macroscale.Precisely what R.Feynman kept in mind more than 50 years ago in his talk at the annual meeting of the American Physical Society at the California Institute of Technology [1].
The theoretical analysis of the various phenomena at the nanoscale is the important stage that has to stimulate the experimental and engineering progress.The study presented in this paper deals with the new phenomenon in the CNT nonlinear dynamics that is the capture of the optical vibrations energy into some domain of the single-walled carbon nanotube.This phenomenon consists in strongly non-uniform distribution of the vibration energy along the CNT axes which is preserved during the time, which exceeds significantly the period of the oscillations in the numerical calculations and forever-theoretically.The opposite process is the slow vibration energy exchange between different domains in the CNT.The latter is similar the beating in the system of two weakly coupled oscillators.Both regimes may be observed in the CNT under the change of the initial energy.The capture of vibration energy may be observed for the low-frequency optical oscillations with the amplitude that exceeds some threshold value, while at the smaller amplitudes only the slow energy migration along the CNT occurs.The derivation of nonlinear equations for the low-frequency oscillations of the CNT as well as the analysis of the bifurcations may be very useful for understanding of the energy exchange in the CNT and for the experimental revealing of transition between different oscillation regimes.
Since the discovery in 1991 [2] up to date, the CNT physical properties excite the great interest of the specialists in the micro-and nanoelectronics as well as the researches in the material science, physics, chemistry, biophysics and related fields of the science and technologies [3].The unique mechanical properties of the CNTs (unsurpassed Young module, large elastic deformation, highest tensile yield strength) and their high electrical and thermal conductivity can dramatically improve the mechanical, electrical and thermal properties of the polymer composites [4,5,6].
The great perspectives in the development of the new super-small and ultra-fast nanodevices stimulate the investigations of the electronic properties of the CNTs in the dependence on their size, structure (chirality), the presence of the impurities or structural defects and others [7] In this connection, the electron-phonon coupling is important to the properties of the CNTs [8,9,10,11].Ballistic transport, excited-state dynamics, superconductivity, Raman spectra all are also connected with electron-phonon interaction [12,13].
The phonons are the carriers in the processes of the thermal energy trans-fer in the CNTs and the composite materials containing the CNTs [14,15].In particular, the problem of the heat transfer in the CNTs [16,17] directly relates to the problem of the thermoconductivity finiteness in one-dimensional anharmonic lattices [18,19,20,21].Both the stationary and nonstationary, but non-resonant, dynamics of the CNTs can be treated in terms of linear or nonlinear normal modes (NNMs).Using the normal mode combinations, one can describe (exactly or approximately) the CNT oscillations under arbitrary initial conditions.However, the situation drastically changes if we deal with non-stationary resonance processes such as energy transfer.In the framework of the linear theory, the energy transfer requires the formation of a wave packet, the time evolution of which depends strongly on the dispersion properties of the system.The dispersion leads to the wave packet spreading that strongly affects the energy transfer efficiency.In the nonlinear systems, the dispersive spreading can be compensated by nonlinearity.
As a result, a soliton (breather) mechanism of energy transfer in the infinite quasi-one-dimensional nonlinear lattices arises.The presence of the localized excitations as well as other nonlinear phenomena may essentially influence on the thermal and electron properties of the CNTs [22].The origin of the envelop solitons (breathers) in the various physical systems has been usually associated with the phenomenon of self-modulation or modulation instability [23,24] (see also the references in [24]).This phenomenon is the interaction between a carrier wave with eigen-frequency ω and the sidebands with the nearest frequencies ω ± δω.The superposition of the sidebands and the carrier waves in the harmonic waves leads to the carrier wave modulation.The nonlinear interaction of waves results to the modulation instability, that is revealed in the amplification of modulation amplitude as well as the growth of the modulation period.The investigation of the self-modulation was initially started in the hydrodynamics and nonlinear optics and it was subsequently expanded in the numerous physical systems including water waves, plasma physics, laser beam, and the theoretical model of various nonlinear lattices [24].
However, it was recently shown [25,26,27] that the resonant interaction of NNMs in the finite lattices leads to the existence of significant non-stationary phenomena, which disappear in the infinite case: i) the intensive energy exchange between different parts of the system, which can be observed at a small enough excitation level, ii) the existence of an instability threshold for the zone-bounding mode, iii) the transition to weak energy localization in some part of the system.All these phenomena can be understood and effi-ciently described by a unified viewpoint in the framework of Limiting Phase Trajectory (LPT) concept.The LPTs correspond to strongly non-stationary processes, which are characterized by the maximum possible (under given conditions) energy exchange between different parts of the system.In the present paper we study the processes of intense energy exchange and the transition to the energy capture of the CNT low-frequency optical vibration.The appropriate description in the framework of the nonlinear thin elastic shell theory was used for the derivation of the nonlinear Schrödinger-type equation and for the investigation of the resonant NNM interaction.The analytical results were confirmed by the molecular dynamics simulation data.A brief preliminary discussion of the revealed phenomena was presented in [28].

CNT low-frequency dynamics in the framework of modified Sanders-
Koiter thin shell theory The dynamics of carbon nanotubes is one of the few areas of solid state physics, in which the classical theory of thin elastic shells (TTES) can be legitimately applied [29].It is noteworthy that, in contrast to macroscopic mechanics, where the fundamental limits of TTES are restricted by the possibility of plastic deformation, this theory can also be used for description of large displacements of CNTs, even in the analysis of their collapse.The only complicating factor is the uncertainty of the parameter characterizing the thickness of the CNT [30].The applicability of a well-designed TTES allows us to obtain an effective description of the vibrational spectrum in the framework of the linear approximation [31,32,33].This can be easily performed for the simplest of the boundary conditions, when a CNT of finite size can be considered as a part of an infinite CNT.However, the modified theory presented below admits efficient study of both linear and nonlinear dynamics of CNTs under arbitrary boundary conditions.
We have found analytically the corresponding part of the oscillation spectrum and estimated the effect of the boundary conditions taking into account the existence of boundary layers [34].As for the non-stationary nonlinear problems for CNT, in particular, resonant intermodal interaction, they can be, in principle, studied by numerical methods [35,36,37].However, such approach is insufficient when dealing with prediction of new phenomena.On the contrary, we show that the analytical approach to nonlinear dynamics of the CNT turns out to be efficient from this viewpoint.
This analytical investigation is based on the reduced nonlinear Sanders-Koiter thin shell theory, and it is a far going extension of our recent study relating to the origin of the oscillation localization in the nonlinear lattices of various types [25,26,27].
Two optical-type vibration branches in the CNT spectrum (fig. 1) are of interest from the viewpoint of the energy exchange processes mentioned above: the first of them is the well-known Radial Breathing Mode (RBM), which is associated with the circumferential wave number n = 0 and corresponds to uniform radial extension-compression.The lowest optical mode in the CNT spectrum (Circumferential Flexure Mode -CFM) is specified by n = 2, and the main deformation is a deviation of the CNT cross-section from the initial circular one [38].To the best of our knowledge, the nonlinear dynamic processes on CNTs were analytically studied only on the basis of a simplest modal analysis (RBM and its parametric instability) [35,36,37].
The CFM oscillations are characterized by relative smallness of ring and shear deformations (in particular, the contour length of the lateral section is not essentially changed during deformation).So, we can assume that only the bending, torsion and longitudinal deformations contribute to the potential energy.
Considering the small-amplitude oscillations of the CNT in the limiting case of a large aspect ratios, one can write the following equation of motion in terms of the radial component of the displacement (see Appendix A for details) where W characterizes the (dimensionless) radial displacement of the shell.Independent variables ξ and τ 0 = ω 0 τ are the dimensionless coordinates along the CNT axis (0 ≤ ξ ≤ 1) and the dimensionless time reduced to the gap frequency ω 0 , respectively.The parameters ω 0 , µ, γ, κ and a 1 depend on the circumferential wave number n (= 2), the inverse aspect ratio of the CNT (i.e. the ratio of the CNT radius R to its length L) α, the ratio of the thickness of the CNT wall h to its radius -β and the Poisson ratio ν.These parameters are determined in the Appendix A (see the equations (A.11)).Introducing a small parameter ε ∼ β one can point clearly the order of the different terms in equation (1) (see the discussion in the Appendix A, equation (A.13)).
Equation ( 1) is a useful tool to analyse the effect of various boundary conditions on the spectrum of natural oscillations of CNT.It was shown in our previous paper [34] devoted to linear CNT vibrations that the solution for periodic boundary conditions can be considered as a basic one for construction of the solutions under other conditions (see brief discussion in the Appendix B).Such conditions as clamping or free edges may be taken into account in the frameworks of the dynamic boundary layer concept, which is extended to nonlinear problems also [39].The frequency spectrum in the case of simply supported edges is determined by the following expression: where k is a longitudinal wave number corresponding to the number of half-waves along the CNT axis.The conditions corresponding to simply supported CNT is a particular case of periodic boundary conditions.As for their realization, the more realistic case is the CNT bounded by hemispheres of fullerene ( in the case of macroscopic cylindrical shell this corresponds to presence of hemispherical shells at the boundaries that is frequently considered by engineers as an approximate realization of the simple supporting).
It is convenient to rewrite the equation (1) using complex variables: where the asterisk denotes the complex conjugation.Performing the multiscale expansion procedure (see Appendix C) one can get the equation for the amplitude of the main order in the "slow" time where the main order value χ 0 is coupled with the complex function Ψ = εχ 0 exp (−iτ 0 ).
First of all, equation (4) admits the plane-wave solution with the dispersion ratio where A is the amplitude.As it can be seen, this dispersion relation is in accordance with the relation (2).
Equation ( 4) is the modified Nonlinear Schrödinger Equation (NLSE), in the standard version of which the fourth derivative is absent.As it is well known, the standard NLSE admits the localized solution -the envelope soliton or the breather.The presence of fourth derivative complicates the problem, but using Pade approximation [40] one can obtain the following localized solution where Solution ( 6) describes a set of soliton-like excitations, which are parametrized by the "frequency" parameter ω.The permissible values of ω are determined by the conditions of the reality of the magnitude (λ) (the inverse soliton width) as well as of the amplitude X 0 .Therefore, these values have to be negative.It is a natural requirement because the localized solutions can exist in the gap of the vibration spectrum.
Figures (2.1) show the solution (6) and its parameters λ and X 0 at the various values of the frequency ω.Equation (4) describes the nonlinear dynamics in the asymptotic limit of the infinitely long CNT.Commonly speaking, it is the only limit when the localized soliton-like excitations occur in the  nonlinear one-dimensional systems, keeping in mind the boundary conditions in the infinity.
Nonlinear equation ( 4) can be used for the analysis of nonlinear normal modes interaction and, in particular, in order to find out the transition between two regimes -the intensive energy exchange and energy localization [28].To perform this, one should take into account that the vibration spectrum for any CNT with finite length is discrete, i.e. the longitudinal wave numbers are integers.To consider the intermodal interaction let us use the sum of the resonant NNMs with the wave numbers k 1 and k 2 .
Substituting solution (8) into equation ( 4) one should use the Galerkin procedure to obtain the equations for complex amplitudes χ 01 and χ 02 : where are the intervals between the modal frequencies.(One can estimate that the frequency shift between the lowest modes (k 1 = 1, k 2 = 2) is approximately twice smaller than that for the next pair of modes (k 2 = 2, k 3 = 3).) It is easy to see that the nonlinear terms in equations ( 9) are separated into two groups: the terms |χ 0j | 2 χ 0i (i, j = 1, 2) determine the nonlinear frequency shift, while the terms χ 2 0i χ * 0j (i = j) describe the nonlinear interaction between modes.The Hamiltonian corresponding to equations ( 9) can be written as Equations ( 9), besides the obvious energy integral (10), possess another integral which characterises the excitation level of the system, and it is an analogue of the occupation number integral in quantum-mechanical terminology.

LPT and localization of CNT vibrations
As it was shown in [25,28], the modal analysis becomes inadequate at the resonance conditions.Therefore we introduce new variables as the linear combinations of resonating modes with preservation the integral X: The new variables describe the dynamics of some domains of the CNT [28] (similarly to some groups of the particles in the effective discrete onedimensional chain [25,26,27]).Considering the energy distribution along the nanotube one can see that the combination φ 1 = 0 and φ 2 = 0 corresponds to a predominant energy concentration in certain domain of the CNT, while the rest of which has a lower energy density.The inverse combination (φ 1 = 0, φ 2 = 0) leads to the opposite energy distribution.In the rest of φ 1 and φ 2 values the energy distributes more uniformly along the CNT's axis.
Because of small difference between frequencies of the modes, the mentioned domains of CNT demonstrate a coherent behaviour similar to beating in the system of two weakly coupled oscillators.Therefore we can consider these regions as new large-scale elementary blocks, which can be identified as specific elements of the system -the coherence domains.The coherence domains were introduced for the nonlinear chain as the "effective particles" in [25].
The existence of integral of motion (11) allows to reduce the dimension of the phase space up to 2 variables -θ and ∆, which characterize the relationship between the amplitudes and the phase shift between the coherence domains, respectively: Substituting these expressions into equations ( 9), the equations of motion in the terms of "angular" variables (θ, ∆) one can obtain: First of all, it is easy to show that equations ( 14) have two stationary points with coordinates (θ = π/4, ∆ = 0) and (θ = π/4, ∆ = π).Taking into account the relations (13,12), one can observe that these points correspond to the stationary states χ 01 and χ 02 , respectively.All trajectories surrounding the stationary points describe the evolution of "mixed" states with different contributions of χ 01 and χ 02 .The trajectory, which separates the NNMs attraction domains, performs the function of separatrix, but it is not possessed of the main quality of the latter.Namely, the motion of imaginary point along this trajectory is continued a finite time.Such a trajectory is the most distant from the stationary points and it is named as the Limiting Phase Trajectory (LPT).As it was mentioned above the last notes the extremely nonuniform distribution of the energy (from a possible ones).
The numerical solutions of Eq. ( 14) with the initial conditions corresponding to the point (θ = 0, ∆ = π/2) for the various values of the excitation X are shown in the Fig.

Figs. 2.2(a, b)
show the evolution of θ(τ 2 ) and ∆(τ 2 ) for small value of X, when the system is close to the linear one.In this case one can see the nonsmooth behavior of the relative amplitudes as well as of the phase shift of the coherence domains φ 1 and φ 2 .Such a behaviour correlates with that the any states belonging to the lines θ = 0 or θ = π/2 with ∆ = (π/2 ± mπ), in fact, are some "virtual" ones, and they must be passed in the infinitesimal time.The respective solutions should be described in the terms of "non-smooth" functions [41].
Figs. 2.2(c,d) demonstrate the behaviour of the functions θ and ∆ if the excitation level X is large enough and is extremely close, but smaller than the threshold value X loc , the origin of which we will consider below.These figures show that this behaviour qualitatively does not different from that in fig.2.2(a,b).
However, Figs.2.2(e,f) exhibit the drastic changes in the evolution of the functions θ and ∆ in spite of excitation of the system was increased of 0.5% only.First of all the variation range of the function θ becomes twice less.It is the most important fact, which shows, that the state with θ = π/2 is inaccessible, if the initial conditions correspond to θ = 0 and vice versa.The second feature in the fig.2.2(f) is the unlimited growth of the function ∆.Such a behaviour corresponds to the transit-time trajectories.
To clarify the bifurcations of the solution of equations ( 14) let us rewrite the Hamilton function (10) in terms of angular variables θ and ∆ and consider the topology of the phase space.
Fig. 4 shows the phase portraits for various values of the parameter X.The initial structure of the phase space for the small X, when the system is close to the linear, is clearly shown in the left panel of Fig. 4. The representative domains of the phase space are bounded by the intervals 0 ≤ θ ≤ π/2 and −π/2 ≤ ∆ ≤ 3π/2.Two stable stationary points correspond to the normal modes χ 01 and χ 02 .The LPT contains two lines θ = 0 and θ = π/2, and two fragments, which connect the pairs of points: The analogous trajectory rounds the stationary state χ 02 .The motion along the LPTs leads to the non-smooth behaviour as it is shown in the Fig. 2.2 However, the stationary state χ 01 becomes unstable if the parameter X exceeds some threshold.Its value X ins can be calculated from the instability condition: Two new stationary points arise after loosing the mode χ 01 its stability.They are new NNMs.The distance between them grows while the parameter X increases.These new stationary points correspond to some non-uniform distribution of the energy along the CNT, however, this non-uniformity is weak.The main feature of these states consists in that any trajectory surrounding them cannot attain the separatrix, which passes through the unstable stationary state χ 01 .Therefore, the non-uniformity of energy distribution remains for the infinite time.Nevertheless, any trajectories, which are situated in the gap between the separatrix and the LPT, preserve the possibility to pass from the vicinity of φ 1 state (θ = 0) into the vicinity of φ 2 state (θ = π/2) (see Fig. 2.2(c, d).This process is accompanied with the slow energy transfer from one part of the CNT to another one and vice versa.
However, the behavior of the solution of equations ( 14) is changed drastically if the value of X overcomes next threshold X loc .The existence of this threshold results from that the new stationary states move away from the unstable state and the separatrix grows while the LPT moves to the unstable state in the vicinity of θ = π/4.The principal changes happen when the LPT reaches the point (θ = π/4, ∆ = 0).At this instant the gap between the LPT and the separatrix disappears and the only trajectory passing from θ = 0 to θ = π/2 is the LPT.The further increasing of the parameter X leads to that new separatrix passing through the unstable stationary points (θ = π/4, ∆ = 0) and (θ = π/4, ∆ = 2π), arises (see Fig. 4(c)).It separates the phase space of the system into uncoupled parts and any trajectories, which start near the θ = 0, cannot attain the value θ = π/2 and vice versa.It means that the energy originally given to in a part of CNT is kept in this part.The new LPTs enclose the stationary points, which correspond to the stable NNMs.  the domain between LPTs and separatrix.Therefore the solution, which is shown in the Fig. 2.2(f), demonstrates the infinite rise of the variable ∆.
The condition of the bifurcation discussed is the degeneration of the energy of the states χ 01 , φ 1 and φ 2 , i.e.
So, the value of the localization threshold turns out to be Fig. 2.2 demonstrates the dependence of the localization threshold in the terms of the radial displacement w from the inverse aspect ratio of the CNT.4) and analogous procedure, based on the direct modal analysis [28], respectively.The dashed red curve shows the threshold value estimated by the numerical method (see Appendix D for detail).Periodic boundary conditions are considered.

MD simulation
According to equation ( 12) the energy distribution corresponding to the first or second effective particles is the predominantly concentrated in the one or other part of the CNT.The motion along the LPT, as it is shown in Fig. 4(a), leads to slow energy transfer from the one effective particle to another one.These processes may be observed in the numerical simulation of the CNT vibration.Also the transition from the regime of energy exchange to the energy capture should be revealed while the amplitude of initial excitation grows.
To verify the results of analytical model the simulation of the low-frequency vibrations of CNTs was performed by molecular dynamics (MD) techniques using the realistic inter-atomic potential functions (force fields).The force fields contain the energies of the covalent C-C bonding, valence and torsion angels deformation, and Van-der-Waals interaction.The potential energy of CNT may be written in the form: The energy of covalent bonding U 1 is written as follows: where R j is the coordinate of j-th atom and ρ 0 = 1.42,Å is the equilibrium length of C − C bond.Only three nearest neighbours are taken into account for each j-th atom in equation (20).The force constant K 1 = 20.564,eV / Å2 .
The potential energy of valence angle deformation is where the equilibrium valence angle φ 0 = 120 0 and the force constant K 2 = 3.852 eV /degree 2 .The second sum in equation ( 19) consists of the triads of the atoms, which form the valence angels.
The energy of torsion angles contains two terms: where the force constants are K 3,1 = 0.1285 eV and K 3,2 = 0.01585 eV and the equilibrium dihedral angle Θ 0 = 0.The torsion terms corresponding to the equilibrium dihedral angle Θ 0 = π are negligible and they were not taken into account in the current simulation [42].The last term in equation ( 19) corresponds to the Van der Waals interaction of the non-covalent bonding carbon atoms: where σ = 1.750Å and the force constant = 1.414 10 −3 eV / Å.One should note, that in spite of that the potential functions (equations ( 20) -( 21)) contain the quadratic forms only, the nonlinearities of geometrical origin are taken into account in accordance with the analytical model of the thin elastic shell ( Appendix A).
Because the analytical model discussed above does not taken into account the chirality of the CNT, we tested the processes of the energy exchange and its localization during the MD simulation of the zigzag CNTs (m, 0) with indexes m = 10 and m = 20, and with the different aspect ratios.
The CNT vibrations in the current system were modelled using the program complex "PUMA".This program uses the Verlet algorithm for the integration of the equations of motion with the time step δt = 0.5 f s.The temperature was controlled by the collision-type thermostat [43].
The typical MD experiment consisted of several stages.At the first stage the CNT was kept at high temperature ( 400K) for structural relaxation.The duration of this stage was approximately 20-25 ps.Then the thermostat temperature was decreased with a constant rate (∼ 1 K/ps) down to approximately 1 K with a subsequent low-temperature relaxation.The third stage dealt with CNT deformation according to analytical solution with subsequent relaxation.The relaxation times for second and third stages were approximately 15-20 ps.(Another version of initial conditions was given by initial velocities of atoms at zero initial displacements.)The short relaxation times during the MD simulation concern to the relaxation of the special prepared state of the CNT, when the initial deformations correspond to a linear combination of the modes considered.In such a case, the differences between the CNT profiles which are described by the analytical expressions and those defined by the MD potentials, were minimal.Therefore, the relaxation times were smaller than those which may be observed in the realistic experiments.Nevertheless, the relaxation processes were continued until the essential parameters of the system became stationary (corresponding to thermodynamic equilibrium).After that the external field was turned off, and the free natural oscillations of CNT with the fixed boundary conditions were realized.To assert the boundary conditions according to the analytical model, the atoms at the edges of CNT were fixed by the force field against any radial displacements (W (0, t) = W (1, t) = 0) that corresponds to the boundary conditions similar to simply supported shell.At that time any longitudinal and circumferential displacements of the edge atoms were not suppressed.The consequent analysis of MD simulation data included the control of natural frequencies and energy distribution along the CNT axis via variation of the oscillation amplitude.The "maps" of energy distribution along the CNT axis measured during the MD simulations are shown in the figures 6 (a-c).Figure 6(a) corresponds to the small value of the excitation level X, when the nonlinearity does not appreciably influence on the modes behaviour and the periodic process that is similar to the beating in the system of two weaklycoupled oscillators, results in the slow energy exchange between different parts of the CNT.Near the localization threshold the process of the energy transfer is changed appreciably (Fig. 6(b)).The period of it essentially grows and the maximum of the energy distribution slowly drifts along the CNT in contrast to the case of small excitation, where the change of the energy location likes the series of jumps (Fig. 6(a)).After overcoming the localization threshold the energy of CNT vibration turns out to be captured in the initially excited part of the CNT (Fig. 6(c)).The reason of that is obvious from considering the phase portrait on figure 4(c).One should notice of that the energy localization is weak enough, because only two modes form the energy distribution.However, for the best of our knowledge, the presence of additional modes leads to the greater energy localization [25].So we can conclude that the bifurcation of the LPT at high excitation level brings to the energy capture in some domain of the CNT.Fig. 2.3 shows the variation of the CNT vibration spectrum with changing of the initial excitation level.The dot, black and red curves correspond to the excitation X < X loc .
One should note, that the large narrow peak near the frequency ω 0.05 corresponds to the lowest eigenvalue of the system under considera- tion.Therefore no vibration with the smaller frequency exists in the gap for any initial excitation, if it does not exceed the localization threshold X loc .However, the overcoming of the localization threshold changes the spectrum drastically.From one side, one can see in Fig. 6(c), that the shape of the initial excitation is deformed essentially during the MD simulation and the local temperature reaches a great value (∼ 1000K).Therefore, the onset of high-frequency modes is naturally sufficient.On the other side, the intensive oscillations in the gap of the spectrum can not be explained by the increasing of the local temperature.To fill this gap, the excitation of another low-frequency modes is required (they may be acoustic type modes like the bending or the longitudinal stretching modes -see Fig. 1).Another possibility is forming of the localized excitation, because the Fourier spectrum is wide enough.Unfortunately, the two-mode approximation used in our analysis can not answer in this question.The accurate study of this problem needs in the consideration of the inter-branch mode interaction.This problem will be formulated in our future studies.

Conclusion
We demonstrate that instability of edge-spectrum optical modes of CNT vibrations is only the preliminary condition of non-stationary (in slow time scale) energy localization in some domain of CNT.The energy capture in one of the CNT parts can be achieved, if the excitation level exceeds the specific threshold X = X loc , which corresponds to merging two trajectories, which are the LPT and the separatrix appeared at X = X ins .When this threshold is exceeded the phase portrait of the system under consideration changes drastically: the separatrix passing through the unstable stationary point (θ = π/4, ∆ = 0) (see fig. 4(b)) encircles the stable stationary point (θ = π/4, ∆ = π) and prevents full energy exchange between effective particles ψ 1 and ψ 2 .Simultaneously a set of transit-time trajectories, which involve any values of phase difference ∆, is created.It means that initial conditions corresponding to identical velocities or displacements of both modes lead to the energy capture by the coherence domains.Then only a partial energy exchange becomes possible along the trajectories, surrounding the stable stationary point and situated inside the separatrix.One should keep in mind that the process of energy capture does not suggest the creation of strongly localized solutions whose formation requires a participation of more components of the spectrum.This can be achieved for CNT with larger aspect ratio.
It should be noted once more that the development and the use of the analytical framework based on the LPT concept is motivated by the fact that resonant non-stationary processes occurring in a broad variety of finite dimensional physical models are beyond the well-known paradigm of nonlinear normal modes (NNMs), fully justified only for quasi-stationary processes and non-stationary processes in non-resonant case.While the NNMs approach has been proved to be an effective tool for the analysis of instability and bifurcations of stationary processes (see, e.g., [44]), the use of the LPTs concept provides the adequate procedures for studying strongly modulated regimes as well as the transitions to energy localization and chaotic behaviour [25].Such an approach clarifies also the physical nature of the breathers formation in infinite discrete or continuum systems.
As a conclusion we would like to note that the phenomenon of energy localization considered above has universal character and it is the common peculiarity of the systems possessing the optical-type branches of vibrational spectrum.However, as it was studied early [25,27,26], the occurrence of the localization depends on the types of nonlinearity as well as on the relations between coefficients σ i,j in the Hamiltonian (10).If some ratios between coefficients σ i,j are satisfied, an additional integral of motion arises that leads to the effective linearization of the equations of motion and, as consequence, to the absence of the localization processes [41].So, in spite of that the interaction of resonating nonlinear modes described by the Hamiltonian (10) is valid for a wide class of the nonlinear systems, the results of this interaction may vary considerably.In any case, the analysis of the Hamiltonian (10) in combination with the LPT concept gives us a useful tool for the study of nonlinear systems.

Acknowledgments
The work was supported by Russia Basic Research Foundation (grant 14-01-00284a ) and Russia Science Foundation (grant 14-17-00255).The authors are grateful to N.K. Balabaev for the assignment of the code of the program complex "PUMA".

Appendix A. The reduced Sanders-Koiter thin shell theory
It is convenient to use the dimensionless variables into expression of the elastic deformation of the circular thin shell.In such a case all the components of the displacement field (u -longitudinal along the CNT axis, v -tangential and w -radial displacement, respectively) are measured in the units of CNT radius R. The displacements and respective deformations refer to the middle surface of the shell.The coordinate along the CNT axis ξ = x/L is measured via the length of nanotube and varies from 0 up to 1, and θ is the circumferential angle.
One can define the dimensionless energy and time variables, which are measured in the units E 0 = Y RLh/(1 − ν 2 ) and t 0 = 1/ Y /ρR 2 (1 − ν 2 ), respectively.Here Y is the Young modulus of graphene sheet, ρ -its mass density, ν -the Poisson ratio of CNT, and h is the effective thickness of CNT wall.There are two dimensionless geometric parameters which characterize CNT: the first one is inverse aspect ratio α = R/L and the second -effective thickness shell β = h/R.
The energy of elastic deformation of CNT in the dimensionless units is written as follows: where ε ξ , ε ϕ and ε ξϕ are the longitudinal, circumferential and shear deformations, and κ ξ , κ ϕ and κ ξϕ are the longitudinal and circumferential cur-vatures, and torsion, respectively.The respective forces and momenta may be written in the physically linear approximation: One should note that both curvatures and torsion are the dimensionless variables in accordance with our definition of displacement field (u, v, w).
Taking into account the relationships (A.2) one can rewrite the expression for the elastic energy deformation as The Sanders-Koiter approximation of defectless thin shell allows to write the nonlinear deformations (ε) and curvatures (κ) in the following form One should make some physically grounded relationships between the displacement components to simplify the description of the CNT nonlinear dynamics.We consider the low-frequency optical-type vibrations which are specified by circumferential wave number n = 2.The deformation of the CNT cross-section for the amplitude of the radial displacement W 0 = 0.2 is shown in figure Appendix A (a).
This branch is characterized by relatively small circumferential and shear deformations, while the displacements themselves may not be small.The comparison of elastic deformation energies associated with different type of deformations is shown in figure Appendix A (b).One can see that the bending deformation (defined by curvatures (A.5)) contains the majority of elastic deformation energy, while the energies of circumferential and shear deformations may be negligible.The energy of longitudinal deformation is considerably smaller than the bending one only for the longitudinal wave numbers k = 1, 2.
In such a case one can formulate the hypotheses: However, these hypotheses don't mean that the displacements included to circumferential and shear deformations are small.
The components of displacement field should be written as follows Relations (A.6) allow us to express the longitudinal and tangential components, and axially symmetric part of the radial displacement, via the radial one.Corresponding relationships can be written as follows: Because the kinetic energy contains the inertial terms corresponding to all components of the deformation field we need in taking into account the relations (A.8) also.
Omitting the calculation details one can write the final equation of motion in terms of radial displacement W (ξ, t): where Eq. (A.10) allows us to calculate the eigenfrequencies in the linear approximation as well as to estimate the effect of nonlinearity on these frequencies at different boundary conditions.The analysis of the different terms of equation (A.10) shows that the essential contribution is given by the first nonlinear term only.In further we skip the rest of nonlinear terms.
The applicability of the thin elastic shell theory requires the smallness of the effective wall thickness β.Therefore, the characteristic frequency of the gap ω 0 is small.One should note, that the parameter α is small enough in most cases.
It is convenient to introduce the 'new' time, which is scaled by the gap frequency ω 0 : τ 0 = ω 0 τ .Then Let us introduce a small parameter ε ∼ β and assume that the inverse aspect ratio α has the same order of smallness.Then, taking into account expressions (A.11) one can evaluate the orders of smallness for the parameters of equation (A.13).
We suppose that the parameters α 1 and β 1 have the value of order of unity.In further we can clearly show the smallness of the parameters of equation (A.13) by assigning the respective order of ε.

Appendix B. Influence of boundary conditions
The presence of boundary conditions different from the simple supporting affects on the NNMs and their frequency.In this Appendix we consider the effective method for the solution of boundary problem and demonstrate the procedure of the normal mode construction on the example of CNT with free edges.
It is intuitively clear that the strong boundary conditions like the clamping lead to frequency growth while the more "soft" conditions can decrease the frequencies.To estimate the variation of normal modes we used the linear approximation of RSKTST (reduced Sanders Koiter thin shell theory).
Let us assume that the solution of linearized equation of the the CNT vibrations is represented as the periodic process Taking into account expression (B.2) explicitly, one can rewrite equation (B.1) with the help of the product of two differential operators: where the parameters µ, γ, and ω are linked by the relationships Because the operators (d 2 /dξ 2 + k 2 ) and (d 2 /dξ 2 − λ 2 ) are commutative ones, any function f (ξ), which satisfies one of the equations where k, λ, C j , j = 1, 2, 3 and ξ j , j = 0, 1 are the constants determined by the boundary conditions and the symmetry of solution.Equation (B.7) shows that the proposed approach allows clearly single out the exponential boundary layer as a part of the solution.One should note that expressions (B.4) provide the coupling between the parameters of the solution.The estimation of the parameters of solution (B.7) is needed for formulation of the boundary conditions in terms of the radial displacements.
Let us consider the vibrations of CNT under condition of free edges.One can show that the free edges boundary conditions correspond to two equations for the radial displacement W (ξ, τ ):

Appendix C. The multiscale expansion
Because we consider the small-amplitude oscillations, one can represent the complex amplitude Ψ as a series of small parameter ε: Next we should introduce the 'time' series: τ 1 = ετ 0 , τ 2 = ε 2 τ 0 , . . .and the respective time derivatives: Substituting the expansion (C.1) into (1) and taking into account the hierarchy of the times, we get the equations in the different orders by small parameter ε.
Then we get: Taking into account the relations (C.3), one can integrate the equation (C.4) with respect to τ 0 and τ 1 .Then the condition of the secular terms absence give us the equation for the main approximation amplitude χ 0 (4).the threshold values of oscillation amplitude.However, the influence of the other part of the spectrum is very important for the estimation of the reliability of the obtained results.Therefore, they should be verified by numerical methods.One of the approaches consists in the direct numerical integration of the modal nonlinear equations of the Sanders-Koiter thin shell theory.
In order to carry out the numerical analysis of the CNT dynamics, a twostep procedure was used: i) the displacement field was expanded by using a double mixed series, then the Rayleigh-Ritz method was applied to the linearized formulation of the problem, in order to obtain an approximation of the eigenfunctions; ii) the displacement field was re-expanded by using the linear approximated eigenfunctions, the Lagrange equations were then considered in conjunction with the nonlinear elastic strain energy to obtain a set of nonlinear ordinary differential equations of motion.
So, to satisfy the boundary conditions the displacement field was expanded into series r(ξ, ϕ, t) = The maximum number of variables needed for describing a general vibration mode with n = 2 nodal diameters (Circumferential Flexural Mode) is obtained by the relation (N max = M u + M v + M w + 3 − p ), where ( M u = M v = M w ) denote the maximum degree of the Chebyshev polynomials and p describes the number of equations for the boundary conditions to be respected.
A specific convergence analysis was carried out to select the degree of the Chebyshev polynomials: degree 11 was found suitably accurate, ( M u = M v = M w = 11).
In the case of a SWCNT with simply supported edges ( p = 8), the maximum number of degrees of freedom of the system is equal to ( N max = 33 + 3 − 8 = 28).
The equations (D.1) were inserted into the expressions of the potential energy E el and kinetic energy T to compute the Rayleigh quotient R(q) = E el,max /T * , where E el,max = max(E el ) is the maximum of the potential energy during a modal vibration, T * = T max /ω 2 , T max = max(T ) is the maximum of the kinetic energy during a modal vibration, ω is the circular frequency of the synchronous harmonic motion and q = [..., U m,2 , V m,2 , W m,2 , ...] T represents a vector containing all the unknown variables.
After imposing the stationarity to the Rayleigh quotient, it was obtained the eigenvalue problem (−ω 2 M + K)q = 0 (D.2) which gives approximate natural frequencies (eigenvalues) and modes of vibrations (eigenvectors).The results of performed calculation show that the eigenspectrum values are in the good accordance with the estimations made in the framework of reduced Sanders-Koiter theory discussed above.The specific difference between eigenvalues amounts to the values 2 − 4% for the long-wave modes and reachs up to 20% while the longitudinal wavenumber grows [34].
In the nonlinear analysis, the full expression of the dimensionless potential energy E el containing terms up to the fourth order (cubic nonlinearity), was considered.
In the case of simply supported boundary conditions, the two low-frequency optical-type circumferential flexure modes (m = 1, n = 2) and (m = 2, n = 2) were considered.
Using the Lagrange equations d dτ ∂T ∂q i + ∂E el ∂q i = 0, (D.3) a set of nonlinear ordinary differential equations was obtained; these equations must be completed with suitable initial conditions on displacements and velocities.This system of nonlinear equations of motion was finally solved by using the implicit Runge-Kutta numerical methods with suitable accuracy, precision and number of steps.The solution of nonlinear equations with initial conditions in the vicinity of the bifurcation threshold shows the coincidence of the threshold values in the analytical model and the numerical one for the wide interval of aspect ratios (see Fig. 2.2).The procedure and results will be discussed in the nearest future.

Figure 1 :
Figure 1: (Color online) The example of the CNT vibration spectrum according to the exact Sanders-Koiter thin shell theory: solid curves correspond to circumferential wave number n = 0, dashed ones -to n = 1 and dot-dashed ones -to n = 2.The insert shows the small longitudinal wave number part of the CFM branch.All the frequencies ω are measured in dimensionless units and k -denotes the number of longitudinal half-waves along the CNT.The aspect ratio of the CNT λ = 20 and the effective wall thickness β = 0.08
Fig. 4(c)  shows the transit-time trajectories, which are in

Figure 5 :
Figure 5: (Color online) Radial displacement W at the localization threshold for the (20,0) zigzag CNT vs inverse aspect ratio α: solid black and doted blue curves are the analytical predictions based on the Eqs.(4) and analogous procedure, based on the direct modal analysis[28], respectively.The dashed red curve shows the threshold value estimated by the numerical method (see Appendix D for detail).Periodic boundary conditions are considered.

Figure 6 :
Figure 6: (Color online) (a) The "map" of the energy distribution along the (20,0) zigzag CNT during the MD-simulation with the small excitation level X = 0.15X loc , that corresponds to phase portrait in the Fig. 4(a).(b, c) The same as in the panel (a) with initial excitation level X = 0.995X loc and X = 1.25X loc , respectively.Insets show the energy measured in kelvins.

Figure A. 8 :
Figure A.8: (Color online) (a) CNT cross-section for circumferential flexure mode with circumferential wave number n = 2. Solid black and red dashed curves show the undeformed and deformed cross-section of the CNT, respectively.R is the radius of the CNT and ϕ is the circumferential angle.(b) The comparison of elastic energy components for the circumferential flexure deformation of the CNT with aspect ratio λ = 20.Dashed, dotted, dot-dashed, and solid curves show the energies of longitudinal, circumferential, shear and bending deformations, respectively.The effective wall thickness β = 0.1 and the Poisson ratio ν = 0.2.The energies are measured in the dimensionless units.

6 )
is a solution of equation (B.3).So a general solution of equation (B.3) is a linear combination of the solutions of equations (B.5), (B.6):

8 )
Solving the equations (B.4,B.8), one can estimate all the parameters of the solution (B.7).Fig. B.9 shows the spectra for the system under different boundary conditions -the simply supported, free and clamped edges.

Figure B. 9 :
Figure B.9: The comparison of CNT vibration spectra for the different boundary conditions.CNT parameters -L = 10 nm and R = 0.79 nm.
n T * m (ξ) cos nϕ f (t) (D.1)where function r(ξ, ϕ, t) substitutes the displacements u, v or w.In equations (D.1) T * m (ξ) = T m (2ξ − 1) are the Chebyshev orthogonal polynomials of the m − th order, n is the number of nodal diameters, and f (t) describes the time evolution of the CNT vibrations.