## Abstract

We show, by molecular dynamics simulations, that 2:1 internal resonance may occur between a radial breathing mode (RBM) and a circumferential flexural mode (CFM) in single-walled carbon nanotubes (SWCNTs). When the RBM vibration amplitude is greater than a critical value, automatic transformations between the RBM and CFM with approximately half RBM-frequency are observed. This discovery in the discrete SWCNT atom assembly is similar to the 2:1 internal resonance mechanism observed in continuum shells. A non-local continuum shell model is employed to determine the critical conditions for the occurrence of observed 2:1 internal resonance between the RBM and CFMs based on two non-dimensional parameters and the Mathieu stability diagram.

## 1. Introduction

The vibrational characteristics of carbon nanotubes (CNTs) have been studied extensively because of their importance in nanomechanical devices and sensoring applications (Sazonova *et al.* 2004; Salvato *et al.* 2008; Wang *et al.* 2008). The vibrations of CNTs can be measured using X-ray diffraction and Raman spectroscopy, and these measurements usually focus on a radial breathing mode (RBM) and other higher frequency modes (D, G and G′ bands; Dresselhaus & Eklund 2000; Dresselhaus *et al.* 2002, 2007; Tang *et al.* 2008). Lower frequency modes than RBMs in Raman spectra have received much less attention in spectrum analysis. It has been shown that the intermittent transformation between the RBM and a lower frequency circumferential flexural mode (CFM) may happen when the amplitude of initial RBM vibration is greater than a critical value for the given geometry of single-walled carbon nanotubes (SWCNTs; Li & Shi 2008). This intermittent transformation between the RBM and CFM may be related to the 2:1 internal resonance, which has been observed in many nonlinear systems (e.g. in continuum shell structures (Nayfeh & Arafat 2006)). The proof of this automated internal resonance in CNTs is significant at least in the following two aspects. (i) It challenges the development of vibration measurements in the time domain (Gambetta *et al.* 2006), which is still not readily available in the study of CNT vibrations. (ii) It may imply other mechanisms for the development of nanomechanical devices and sensors, for example, through the interplay between the vibrational mode and electron tunnelling in CNTs (LeRoy *et al.* 2004).

In this study, we employ molecular dynamics (MD) to study the 2:1 internal resonance between the RBM and lower frequency CFMs in SWCNTs. For a given SWCNT, its vibration process with the initial, pure RBM is simulated. The displacement–time history of an SWCNT atom is obtained in the time domain. Frequency spectra are obtained from displacement–time histories using a fast Fourier transformation, which can be compared with Raman spectrum measurements. A non-local continuum model is employed to explain the observed 2:1 internal resonance in SWCNTs successfully. This study may open the door to the design of experiments to demonstrate this nonlinear internal resonance in nanostructures.

## 2. Molecular dynamics method and simulation results

In MD simulations of vibration behaviours of SWCNTs, the atomistic force-field COMPASS (condensed-phase optimized molecular potential for atomistic simulation studies) (Sun 1998) is employed to govern the motion of each carbon atom in the SWCNT. The COMPASS force field is parameterized based on *ab initio* quantum mechanics (Chen & Cao 2006). In comparing with other more accurate *ab initio* MD methods, this force field makes the MD simulations more time-efficient. Armchair (10,10) SWCNT is used to demonstrate the representative MD simulation results. For an infinitely long, standalone armchair (10,10) undergoing pure RBM vibration, due to geometry and deformation equivalence along its axis, a lattice unit of the whole structure shown in figure 1 is chosen for the simulation of a three-dimensional, constant-volume and -energy (NVE) ensemble via the Forcite Plus module in commercial MD software Materials Studio^{®;}.1

In the first step of simulation, the lattice unit is geometrically optimized by minimizing enthalpy of the whole system to reach a stable state. This is achieved by using the function ‘Geometry Optimization’ of the Forcite Plus module, which reduces the radius of the SWCNT from 6.78 nm to approximately 6.50 nm. In the next step, all atoms are imposed an initial speed of 1 nm ps^{−1} in their respective radial directions to initiate RBM vibration for the unit structure. The following vibration process is simulated through the ‘Dynamic Analysis’ function from the same module of Materials STUDIO. The simulation time step is 1 fs and the initial temperature is set to room temperature, i.e. 298 K. It is observed that RBM vibration continues for about 4 ps and then one CFM is gradually excited. This CFM dominates vibration for nearly 2 ps before the vibration is switched back to the initial RBM. The lattice unit experiences this intermittent vibration mode transformation in the subsequent vibration period (the total simulation time is 20 ps). Seven snapshots of the lattice unit perpendicular to its cross-section are shown in figure 2 in the duration of the first two vibration mode transformations. It can be seen that when the CFM is excited, the cross-section of the lattice unit changes between the RBM (circles) and CFM (triangles). Variation of the radial displacement of one atom versus time is plotted in figure 3, where displacement is normalized by *a**v*_{0}/*c*, in which *a* is the radius of the SWCNTs, *v*_{0}=1 nm ps^{−1} is the amplitude of the initial RBM vibration velocity and is the phonon velocity, where *E*, *ν* and *ρ* are the nominal elastic modulus, Poisson’s ratio and mass density of the SWCNT, respectively. Here, *h* is a representative effective thickness of the SWCNT. Vibrations dominated by the RBM and CFM are indicated by solid and dashed arrows, respectively, in figure 3. It is evident that the RBM vibration becomes unstable after about 4 ps, and the CFM and RBM have intermittently dominated the vibration since then. More simulations on the same lattice unit with different initial RBM velocities show that only when the amplitude of the initial RBM velocity is sufficiently large, the vibration mode transformation may happen after a finite period of RBM-dominated vibration. It is found that when *v*_{0}≤0.22 nm ps^{−1}, stable RBM vibration continues in the whole simulation period of 50 ps. When *v*_{0}=0.23 nm ps^{−1}, the CFM is excited at about 46 ps. Therefore, the required initial velocity for the occurrence of vibrational mode transformation between the RBM and CFM of armchair (10,10) can be estimated to be 0.23 nm ps^{−1}. Frequency analyses of vibrations with different initial velocities of 0.22, 0.5, 0.75 and 1.0 nm ps^{−1} in a 20 ps duration are shown in figure 4. It can be seen that when the initial velocity is small (0.22 nm ps^{−1}), only one resonance peak appears on the frequency spectrum that corresponds to a stable RBM vibration. However, when the initial velocity exceeds a critical velocity value, RBM vibration becomes unstable after a finite time period, which is shown in the spectrum by two or more groups of frequency peaks, one group is around the original RBM frequency and the other group is around half the RBM frequency corresponding to CFMs excited from the initial RBM.

For armchair (10,10), its critical initial RBM vibration velocity (*v*_{c}≈0.23 nm ps^{−1}) is initially estimated based on several simulations with different initial RBM vibration velocities in a given time period. Then a more accurate prediction of *v*_{c} is obtained using the following procedure. A limited number of simulations (five in the present study) with different initial velocities close to the estimated critical velocity *v*_{c}≈0.23 nm ps^{−1} are imposed on the armchair (10,10). The initiation time (*t*_{0}) of the first mode transformation for each of the five simulations can then be obtained, which are plotted inversely against the imposed *v*_{0} values in figure 5. As it is naturally expected that *t*_{0} approaches infinity or equivalently 1/*t*_{0} approaches zero as *v*_{0} approaches *v*_{c}, the dependence of 1/*t*_{0} on *v*_{0} can be approximated by 1/*t*_{0}=*a*(*v*_{0}−*v*_{c})^{b}, which can be determined by fitting the simulation data as 1/*t*_{0}=0.252(*v*_{0}−0.213)^{0.595}. Therefore, the simulation data are extrapolated to give *v*_{c}=0.213 nm ps^{−1}, which is close to the estimated value of *v*_{c}=0.23 nm ps^{−1} with a relative error of 8.0 per cent.

The initial RBM vibration velocities cannot be perfectly uniform. The imperfections in initial RBM velocities are introduced to the pure RBM vibration velocity field through a perturbation velocity field with stochastic directions and amplitudes in two different ways. First, for an initial RBM vibration with velocity amplitude *v*_{0}=1.0 nm ps^{−1}, a perturbation velocity field of fixed stochastic directions and amplitudes, but with varying maximum perturbation velocity amplitude *v*_{p} (*v*_{p}/*v*_{0}=0.001, 0.01 and 0.1, respectively), is introduced. Its effects on the radial displacement of an atom is plotted in figure 6, where the case *v*_{p}/*v*_{0}=0 corresponds to pure RBM vibration. With the increase of perturbation strength, the transformation from RBM to CFM occurs at an earlier time. Second, for an initial RBM vibration with velocity amplitude *v*_{0}=0.25 nm ps^{−1}, perturbation velocity fields with varying stochastic directions and amplitudes, but fixed maximum perturbation velocity amplitude *v*_{p}/*v*_{0}=0.1, are introduced and their effects on radial displacement of an atom are shown in figure 7. Once again, we can observe the change of the transformation time from RBM to CFM due to the variation of the perturbation velocity field. It is concluded that the first transformation time from RBM to CFM depends on the distributions of stochastic direction and amplitude of the perturbation velocity field. Therefore, the required initial RBM velocity amplitude of 0.213 nm ps^{−1} for the occurrence of mode transformation, which is based on pure initial RBM vibration, is an idealized prediction. The actual value is directly related to the degree of imperfections of the initial RBM velocity field. In practice, there are many factors that may affect the uniformity of the initial RBM velocity field, e.g. excitation methods, environmental media effects, etc.

A large amount of simulations on a selection of armchair lattice units from (5,5) to (11,11) with different initial RBM velocity fields indicate that, when the initial RBM velocity is sufficiently large, the RBM can become unstable after a finite period of RBM vibration and the CFM will be excited. Frequencies of the RBM and three CFMs (the excited one and the nearest two) can be obtained by linear modal analysis of individual lattice unit structure, and the results are presented in table 1. Though the third flexural mode (i.e. CFM-3) is excited for all armchair lattice units (5,5) to (11,11), it does not imply that CFM-3 will still be excited for larger armchairs. For example, CFM-4 is excited for armchair (15,15). It should be pointed out that the cross section of larger armchair lattice units will not remain circular after geometry optimization and, therefore, simulations on them are directly implemented without the geometry-optimization step. This may influence the mode transformation time and the critical initial RBM velocity, but not the excited CFM.

It is shown in table 1 that the CFM that possesses a frequency value closest to half of the RBM frequency will be excited when the initial RBM vibration is strong enough. This indicates that the well-known 2:1 internal resonance mechanism observed in many nonlinear dynamic systems still happens in nanostructures, at least for the cylindrical-shell-like structures such as SWCNTs.

## 3. Non-local elastic shell model

Continuum thin-shell model, which neglects non-local atomic interactions, can successfully predict the RBM frequency of SWCNTs (Mahan 2002; Wang & Hu 2005; Chico *et al.* 2006). When flexural waves are involved, non-local continuum model may take the long-range forces of nanotube atoms into account (Zhang *et al.* 2004, 2005; Wang & Hu 2005; Hu *et al.* 2008). The non-local elastic constitutive equation includes strain gradient effects by introducing a non-local length *l*, i.e.
3.1
where *E*_{1}=*E*/(1−*ν*^{2}); *E* and *ν* are the elastic modulus and Poisson’s ratio, respectively; ∇^{2} is the Laplace operator; and *σ*_{θ} and *ϵ*_{θ} are the circumferential stress and strain in the SWCNT, respectively. Equation (3.1) is the non-local constitutive form for plane–strain elasticity problems (Eringen 1983).

The RBM and CFM frequencies based on non-local continuum shell model are
3.2
and
3.3
where *α*^{2}=*h*^{2}/12*a*^{2},*c*_{0}=2.998×10^{10} cm s^{−1} is the speed of light and *n* is the wavenumber of the inextensional flexural vibration wave.

The nonlinear differential equations in Li & Shi (2008) are further extended to include non-local effects described by equation (3.1), as shown in appendix A, which can be reduced to Mathieu equations to describe the coupling between the RBM and inextensional flexural modes (Goodier & McIvor 1964; Li & Shi 2008). The two non-dimensional parameters in the Mathieu stability diagram with the consideration of non-local effects are
3.4a
and
3.4b
where . If (*Ω*_{n},*μ*_{n}) falls in the unstable region in the Mathieu stability diagram, vibration-mode transformations between the RBM and those involved flexural modes will occur. Based on the Mathieu stability diagram and the (*Ω*_{n},*μ*_{n}) values given in equations (3.4a) and (3.4b), we examine the possible 2:1 internal resonance observed in MD simulations.

We determined the thickness of SWCNTs using in-plane stiffness *E**h*≈360 J m^{−2} and in-plane bending stiffness *D*=*E**h*^{3}/12(1−*ν*^{2})=2.0 eV with Poisson’s ratio *ν*=0.19 (Yakobson *et al.* 1996; Wang & Hu 2005), which gives an approximate effective thickness of 0.1 nm for SWCNTs and phonon velocity of *c*≈22 nm ps^{−1}. For SWCNTs, the non-local length *l* is related to the carbon bond length and the wavelength of the flexural deformation mode (Duan *et al.* 2007). Frequently, its value is fixed for a particular study, e.g. it was suggested to be 0.116 nm in the study of buckling strain in multi-walled CNTs (Zhang *et al.* 2005). The non-local length values in this study were determined by comparing the calculated frequencies of different vibration modes with the corresponding frequencies obtained from MD simulations, which are shown in table 2.

It has been shown in MD simulations (table 1) that most vibration mode transformations between the RBM and CFM belong to 2:1 internal resonance. We can use armchair (10,10) to explain the 2:1 internal resonance observed in SWCNTs in MD simulations. According to table 2, the non-local length values vary from 0.112 to 0.449 nm for the vibration modes from the RBM to CFM-4 in armchair (10,10). We used non-local length values from 0.10 to 0.25 nm with increment of 0.05 nm in the non-local model to find out which flexural modes are likely to be excited from the RBM, as shown in figure 8. Since *μ*_{n} increases with *Ω*_{n} monotonically and the unstable regions in the Mathieu stability diagram become narrower with the increase of *Ω*_{n}, it is clearly shown that 2:1 internal resonance is the most likely mechanism for the transformation between the RBM and other CFMs. The normalized initial radial velocity in figure 8 is *v*_{0}/*c*=0.042 nm ps^{−1}, which is related to the RBM excitation power intensity. If this is reduced, those points in unstable regions will move down vertically and may enter the stable region. The radial displacement history in time domain and the corresponding frequency spectrum for *v*_{0}/*c*=0.04 nm ps^{−1} and *l*=0.1 nm are shown in figure 9, which clearly shows the appearance of the RBM and CFM-3. On the other hand, figure 8 also demonstrates the necessity to include the non-local effects in the continuum model for the correct prediction of 2:1 internal resonance. In Li & Shi (2008), where a local continuum model is employed using the old thickness parameter 0.066 nm (Yakobson *et al.* 1996), mode transformation between the RBM and CBM-4 was predicted. If the excitation power intensity is further reduced, CBM-5 will be involved under 1:1 internal resonance mechanism, as shown in figure 10. These local model predictions do not agree with the present MD simulation result where only CMF-3 was involved in mode transformation with the RBM, and it seems that the inclusion of non-local length in the continuum model strongly supports the MD simulation results.

## 4. Conclusions

Based on MD and non-local continuum model, it is shown that the 2:1 internal resonance in SWCNTs controls non-linear coupling between the RBM and other circumferential flexural vibration modes. The initiation of this 2:1 internal resonance is controlled by two non-dimensional numbers in the Mathieu stability diagram. Due to the significances of this 2:1 internal resonance phenomenon on the applications of the SWCNT’s dynamic behaviour, further numerical and experimental verifications are necessary.

## Acknowledgements

M. X. S. acknowledges the financial support from Dorothy Hodgkin Postgraduate Award (EPSRC-Shell DHPA). The authors also thank Dr Gábor Csányi for suggesting the calculation method of critical velocity of mode transformation.

## Appendix A

Take the SWCNT as a thin cylindrical shell where *a* and *h* are the radius and thickness, respectively. Assume that the shell undergoes small plane–strain deformation, while surface normals are reserved and thickness remains unchanged during deformation. A point on initial middle surface (*a*,*θ*) will move to a point represented by (*r*,*ϕ*), where *r* and *ϕ* are functions of initial angle *θ* and time *t*, respectively. Two non-dimensional radial and tangential deflection variables *ξ*=1−*r*/*a* and *ψ*=*ϕ*−*θ* are introduced and expressed as
A 1
and
A 2
In the above expansions, the term *a*_{0} represents RBM vibration. The other term *c*_{n} refers to inextentional flexural deformation, which corresponds to the deformation caused by CFM vibration. Using the same governing equations in Li & Shi (2008), the unknown coefficient functions *a*_{0} and *c*_{n} in equations (A 1) and (A 2) can be obtained by means of the Lagrange equation as
A 3
and
A 4
where and correspond to the additional contributions due to the non-local term in non-local constitutive equation
A 5
and are given as
A 6
and
A 7
In equations (A 3)–(A 7), dots above symbols indicate partial derivatives with non-dimensional time variable *τ*=(*c*/*a*)*t*, [···] in equation (A 7) represents higher-order coupled terms between flexural modes *c*_{m} and *l* is the non-local length, *E*_{1}=*E*/(1−*ν*^{2}),*E* and *ν* are the elastic modulus and Poisson’s ratio of material, respectively.

Consider the modal expansion of *ξ* and *ψ* given in equations (A 1) and (A 2), the initial conditions of RBM vibration with imperfections for an SWCNT can be given as
A 8
where the term *v*_{n} represents the imperfections to the uniformity of the initial RBM vibration field *v*_{0} and are normally at least two orders smaller than *v*_{0}. Correspondingly, the initial conditions of *a*_{0} and *c*_{n} are
A 9
Initially, the amplitudes of *c*_{n} (*n*≥2) are small because the non-uniformity is small in initial condition equations (A 9). Then, when terms *c*_{n} are neglected in equation (A 3), it will be reduced to
A 10
where
A 11
A harmonic vibration solution of equation (A 10) with initial conditions in equation (A 9) is
A 12
With solution *a*_{0} given in equation (A 12) and discarding those small terms in equation (A 4), it can then be reduced to a Mathieu-type equation
A 13
where dots above *c*_{n} refer to the partial derivative with respect to a new non-dimensional time variable *η*=*k**τ*, and *Ω*_{n} and *μ*_{n} are given in equations (3.4a) and (3.4b).

## Footnotes

↵1 http://accelrys.com/products/materials-studio/(accessed on April 30, 2009).

- Received March 22, 2009.
- Accepted June 26, 2009.

- © 2009 The Royal Society