## Abstract

Asymptotic distribution of the discrete spectrum eigenvalues for a recently developed model of a double-walled carbon nanotube is presented in this paper. The corresponding initial boundary-value problem is reduced to an evolution equation, whose dynamics generator is a non-self-adjoint matrix differential operator with a purely discrete spectrum. It is shown in the paper that the entire spectrum asymptotically splits into four spectral branches. Asymptotic representation is derived for the eigenvalues along each spectral branch as the number of an eigenvalue tends to infinity. To prove the results, a two-step procedure involving the construction of the left- and right-reflection matrices has been used.

## 1. Introduction

This paper is devoted to the rigorous mathematical analysis of a specific model describing vibrations of a double-walled carbon nanotube (DWNT). The problem is represented in the form of a system of four partial differential equations modelling two coupled Timoshenko beams. This model is equipped with general linear dynamic boundary conditions. The corresponding initial boundary-value problem can be reduced to the evolution equation in the state space, which is a Hilbert space with the so-called energy norm. The dynamics generator of this evolution equation is our main object of interest. This dynamics generator is a non-self-adjoint matrix differential operator whose asymptotic and spectral properties are discussed in the paper. The results presented in this paper are the generalization of our results in Shubov & Rojas-Arenaza (2008, 2010*a*,*b*). In all of the above papers, the model is given in the form of a system of four partial differential equations equipped with the four-parameter family of physically meaningful dynamic boundary conditions. These parameters, *α*^{(i)},*β*^{(i)}, *i*=1,2, may be any complex numbers. The problem is technically so complicated that we started in Shubov & Rojas-Arenaza (2010*b*) with the case when *α*^{(1)}=*α*^{(2)} and *β*^{(1)}=*β*^{(2)}. The above restrictions on the boundary parameters yield a four-branch discrete spectrum with only two horizontal asymptotes (i.e. the eigenvalues of the problem are asymptotically of double multiplicities). In the present paper, we are able to overcome technical difficulties and eliminate the above restrictions on the boundary parameters, namely by assuming that *α*^{(1)}≠*α*^{(2)} and/or *β*^{(1)}≠*β*^{(2)}. In the case when either *α*^{(1)}=*α*^{(2)} and *β*^{(1)}≠*β*^{(2)} or *α*^{(1)}≠*α*^{(2)} and/or *β*^{(1)}=*β*^{(2)}, the discrete spectrum asymptotically splits into three sets approaching three distinct horizontal asymptotes with one of them being degenerate (i.e. containing asymptotically double eigenvalues). When *α*^{(1)}≠*α*^{(2)} and/or *β*^{(1)}≠*β*^{(2)}, there are four distinct horizontal asymptotes with asymptotically simple eigenvalues along each asymptote.

Before presenting the main results, we provide here a brief overview of contemporary literature related to the treatment of carbon nanotubes (CNTs) as continuum mechanics models (i.e. flexible beams, shells, membranes). CNTs are cylindrical structures closed by endcaps, which are hemispheres of fullerenes, with diameters in the nanometre range and lengths in the micron range. There are two basic types of nanotubes: single-walled nanotubes (SWNTs) that are exactly one atom thick and multi-walled nanotubes (MWNTs) that are essentially concentrically nested SWNTs. The individual tubes are not bonded together; they only interact through weak non-bonded van der Waals forces. The tubes remain free to slide and rotate independently with small resistive forces. The primary structural unit of a nanotube is a hexagonal ring consisting of six carbon atoms connected through the so-called σ- and π-bonds. The σ-bonds are very strong and produce tubes with excellent in-plane properties. The π-bonds are delocalized bonds and are much weaker; they are primarily responsible for the out-of-plane properties, such as wall-bending stiffness. Nanotubes have numerous atomic structures, which are determined by the roll-up angle, termed helicity or chirality and diameter of the hexagon structure (see Gibson *et al.* 2007 and references therein; Jakobson *et al.* 1996; Govindjee & Sackman 1999; Jamieson 2000).

Scientific and engineering communities have fully recognized the excellent properties of CNTs and their wide-ranging potentials. Jamieson 2000 argued that nanotechnology, mainly due to CNTs, may have a greater impact on technology than the silicon revolution. Depending on their atomic structure, CNTs have electrical properties that range from metallic to semiconducting. The mechanical properties of CNTs are also unique. They possess exceptionally high specific stiffness and strength yet are also extremely elastic, being able to bend 360^{°} without noticeable damage. The potential applications for materials with these properties are immense in fields such as nanotechnology, electronics, optics, material science and architecture (Yoon *et al.* 2002, 2003, 2004; Pantano *et al.* 2004; Wang *et al.* 2006).

Developing mathematical models of CNTs is of critical importance. These models have to be verified and quantified by performing and analysing experiments. In contemporary mathematical and engineering literature, there are two groups of models: continuum mechanics models and molecular simulation models. The continuum models are generally based on traditional engineering models such as beams, shells or membranes. They treat nanotubes as continuous materials with definite geometries and common material properties such as *Young*’*s modulus*. In contrast, molecular models consider each atom and mathematically define the interactions between the atoms. Often, the molecular models are reduced to continuum models so that the results can be compared with other continuum models and other materials directly (Arroyo & Belytschko 2003; Pantano *et al.* 2003, 2004; Yoon *et al.* 2003; Wang *et al.* 2004, 2005).

Since in our research we carry out mathematical analysis of one of the continuum mechanics models, we will focus our attention on the discussion of several continuum models presented in the literature. Based on their work on atomic simulations of CNTs, Jakobson *et al.* (1996) provided justification for incorporating continuum mechanics models into the study of CNTs, stating that ‘The laws of continuum mechanics are amazingly robust and allow one to treat even intrinsically discrete objects only a few atoms in diameter’.

Numerous choices for continuum models facing the analyst are effectively reduced to those of geometry, assumptions on deformations and constitutive laws. The validity of the choice is determined by the ability of the overall model to properly predict mechanical responses. The assumptions on the nanotube deformation suggest the type of continuum model. The most commonly used models are the Euler–Bernoulli beam model, Timoshenko beam model, flexible shell model and membrane model. Typically, models for MWNTs allow for independent wall movement, and the wall interaction is a function of the local wall separation distance. In a number of research papers, the pressure is used as the most direct equivalent to the van der Waals forces. In Ru (2001), the tube walls of a DWNT were modelled as shells whose bending stiffness depends cubically on the wall thickness. A similar approach has been used in Pantano *et al.* (2003, 2004), where the walls were modelled as thin shells and the inter-wall interactions were modelled as pressures. The authors validated this model for SWNTs and MWNTs by comparing their results with molecular mechanics simulations and experimental results on bending pressure loads. They found good predictions for strain energies, local buckling strain initiations and kinking modes for SWNTs in bending. Arroyo & Belytschko (2003), rather than using thin shell theory, developed a membrane wall model based on the generalized crystal elasticity theory. The resulting model is hyperelastic and depends on the surface curvature and the crystal lattice orientation. Based on the potential for bond deformation and van der Waals forces, Li & Chou (2003) created a truss model for a DWNT. Wang *et al.* (2005) developed a multiple elastic shell model based on the Flugge elastic shell theory in order to study vibrations of MWNTs. Using a shell model, radial breathing modes were included along with the flexural modes governed by beam models, demonstrating that interlayer van der Waals forces have a significant effect on radial breathing modes of large radius MWNTs, a smaller effect on radial modes of smaller radius MWNTs, and a practically negligible effect on the torsion and longitudinal modes of MWNTs. Wang *et al.* (2004) discussed the development of a simplified elastic shell model for buckling and free vibrations of CNTs based on the shell equations of Flugge and Donnell. Classical theory of elasticity solutions for the vibrational modes of thin-walled hollow isotropic cylinders have been applied by Mahan (2002) to calculate all the low-frequency modes of single-wall CNTs, including flexural, torsional, radial breathing and longitudinal modes.

Gibson *et al.* (2007) examined several cases where modelling of nanotube vibration is carried out by using either the Euler–Bernoulli beam theory or the Timoshenko beam theory, depending on the diameter-to-length ratio. As the authors pointed out ‘it is important to have accurate theoretical models for the natural frequencies and mode shapes of carbon nanotubes…. For example, if the nanotubes are to be used as nanomechanical resonators, the oscillations frequency is a key property of the resonator. In addition, the effective elastic modulus of a nanotube may be indirectly determined from its measured natural frequencies or mode shapes if a sufficiently accurate theoretical model is used’. Yoon *et al.* (2002) have studied non-coaxial intertube resonances of MWNTs (i.e. those modes of vibration which involve relative motions or intertube displacements between the concentric tubes). They have shown that such resonances can be modelled for an *N*-wall tube by modifying the Euler–Bernoulli equation to include van der Waals interaction forces and solving *N* coupled partial differential equations. Zhang *et al.* (2004) used a similar approach to model van der Waals interactions between transversely vibrating DWNTs under a compressive axial load, showing that increasing the compressive load causes a reduction in the frequencies. Yoon *et al.* (2004) demonstrated that a double Timoshenko beam model characterizes the motion of the inner and outer tubes in a DWNT. Accordingly, the model was used to study the effects of rotary inertia and shear deformation on terahertz frequency wave propagation in DWNTs. It was concluded that, because of the relative motion between the inner and outer tubes at high frequencies, the Timoshenko beam model is more relevant than the Euler–Bernoulli model for the terahertz frequency range.

For MWNTs embedded in an elastic medium such as a polymer matrix, the effect of the matrix on the resonance frequencies of the MWNT has been estimated by using a beam-on-elastic foundation model (Yoon *et al.* 2003). In this model, the system of equations has been modified by adding a distributed elastic reaction force to the equation of motion for the outer tube in a MWNT.

Xu *et al.* (2006) studied the vibration of a DWNT created by nonlinear interlayer van der Waals forces. The interlayer forces are described as a nonlinear function of the interlayer spacing. The inner and outer CNTs are modelled as two individual elastic beams. Detailed results are demonstrated for DWNTs, based on the simply supported, fixed or free-end conditions, respectively. The results indicate that the nonlinear factors of the van der Waals forces have little effect on the coaxial free vibrations but greatly affect non-coaxial free vibrations. As has been indicated (Qian *et al.* 2002), although CNTs can have diameters only several times larger than the length between carbon atoms, continuum models have been found to describe their mechanical behaviour very accurately under many circumstances. Indeed, their small size and presumed small number of defects make CNTs ideal systems for the study of the links between atomic motion and continuum mechanical properties such as Young’s modulus and yield and fracture strength. Govindjee & Sackman (1999) studied the validity of modeling MWNTs with the Euler beam theory. They showed size dependence of the material properties at the nanoscale, which does not appear in classical continuum mechanics. From scaling analysis, Harik (2001) proposed three non-dimensional parameters to check the validity of the continuum assumptions.

Now, we turn to the description of the problem and our paper’s main results. We carry out a rigorous mathematical analysis of an initial boundary-value problem modelling small transversal vibrations of a DWNT. The system of equations studied in both this paper and our previous ones is similar to those studied in a number of papers (see Jakobson *et al.* 1996; Ru 2001; Qian *et al.* 2002; Pantano *et al.* 2003; Yoon *et al.* 2003; Wang *et al.* 2006; Xu *et al.* 2006; Gibson *et al.* 2007). The physical system consists of two nested nanotubes interacting through the distributed van der Waals forces; each nanotube is modelled as a Timoshenko beam with specific parameters. As pointed out by Wang *et al.* (2006), ‘Unlike the Euler–Bernoulli beam model, the Timoshenko beam model allows for the effects of transverse shear deformation and rotary inertia. These effects become significant for CNTs with small length-to-diameter ratios that are normally encountered in applications’.

We would like to emphasize that, contrary to papers using numerical analysis methods to deal with continuum mechanics models, we present purely analytical results on the asymptotic and spectral analysis of the model. The model is given in the form of two coupled Timoshenko beams (i.e. in the form of four coupled hyperbolic partial differential equations). While extensive mathematical literature exists devoted to the different aspects of the individual Timoshenko beam model, such as asymptotics, spectral and stability analysis and control, to the best of our knowledge, there is a very limited number of results on the Timoshenko beam model coupled with other types of beam models (for examples of analyses of a coupled Timoshenko and Euler–Bernoulli beam model, see Balakrishnan *et al.* (2004), Shubov & Peterson (2004) and Shubov (2006)). The system is equipped with a set of non-self-adjoint boundary conditions involving four independent complex parameters. It is convenient to rewrite this system as the first order in a time-evolution equation in the state space, which is a Hilbert space with the energy norm. Our main object of interest is the dynamics generator of the aforementioned evolution equation. The dynamics generator is an 8×8 non-self-adjoint matrix differential operator with purely discrete spectrum of the normal eigenvalues (see theorem 2.3 and corollary 2.4). By calculating asymptotic distribution of the eigenvalues of the dynamics generator, one can immediately obtain the information on the asymptotic distribution of the vibrational frequencies and the mode shapes of the actual physical system.

The present paper is organized as follows. In §2, we introduce the main Hilbert space and rewrite the system of four equations (see (2.1)) with the set of eight boundary conditions (see conditions (2.2)–(2.6)) as the first order in the time-evolution equation. This reduction allows us to introduce a matrix differential operator (see equations (2.11) and (2.16)), whose spectrum is directly connected to the vibrational frequencies of the system of coupled nanotubes. The main result of the paper is formulated in theorem 2.5.

In §3, we begin a systematic study of the spectral equation 1.1 where is the aforementioned dynamics generator. In this section, we focus on the asymptotic representation of the fundamental system of solutions of the above equation as and represent our results as theorem 3.3.

In §§4–6, we construct two solutions of equation (1.1), the first of which satisfies only the left-hand-side boundary conditions and the second one satisfies only the right-hand-side boundary conditions. To overcome technical difficulties resulting from a direct application of eight boundary conditions, we have applied the so-called two-step procedure developed in our previous works. A detailed description of this procedure is given in §4. In particular, we construct a specific matrix-valued function of *λ*, which we call *the left reflection matrix*. With the help of the left reflection matrix we derive a representation of a specific solution of equation (1.1) that satisfies the left-hand-side boundary conditions without any limitations on the behaviour at the right-hand side of the structure. We point out that such a solution is analogous to the well-known Jost solution in quantum and acoustical scattering (see Shubov 1995 and references therein).

In §5, we derive another matrix-valued function of *λ*, which we call *the right-reflection matrix*. Using this matrix, we obtain another solution of equation (1.1) that satisfies the right-hand-side boundary conditions. It is the right-reflection matrix that contains information on the boundary parameters (see conditions (2.3)–(2.6)).

By appropriately combining the reflection matrices, we obtain the key equation (see equation (5.1)), whose solutions coincide with the eigenvalues of the operator .

## 2. Statement of problem: main results from Shubov & Rojas-Arenaza (2008, 2010*a*)

A system of two Timoshenko beams coupled through the van der Waals force of strength can be given in the following form (Yoon *et al.* 2003, 2004; Wang *et al.* 2006; Xu *et al.* 2006):
2.1
Here *x*∈[0,*L*], , *t*≥0. In system (2.1), the following notation has been used: *W*^{(i)} is the transverse displacement of the *i*th tube, *i*=1,2; *Φ*^{(i)} is the slope of the *i*th tube due to bending deformation alone; *I*_{i} is the second moment of area of cross section; *A*_{i} is the cross-sectional area; *ρ* is the mass density per unit volume; *k*_{i} is shear correction factor; *E* is Young’s modulus; *G* is the shear modulus.

### Remark 2.1.

Remark on van der Waals intertube interaction.

Since the inner and outer tubes of a DWCNT are originally concentric and the van der Waals interaction is determined by interlayer spacing between two tubes, the net van der Waals force remains zero if the tubes vibrate coaxially (i.e. they share the same deflection curve). Hence, the van der Waals interaction plays no role in the simple beam model for isolated double-walled nanotubes. For the double-beam model (Yoon *et al.* 2003, 2004; Xu *et al.* 2006), however, the inner and outer tubes are described by two individual deflection curves and they are not necessarily coincident. Therefore, for small-amplitude non-coaxial linear vibrations, the van der Waals force at any position *x* between the two tubes depends linearly on the difference of their deflection curves at that position.

Together with the system of equations, we introduce the following set of boundary conditions. The boundary conditions at *x*=0 are
2.2
At the right-hand side of the structure, we impose the most general set of the boundary conditions involving four numerical parameters , , *i*=1,2, i.e.
2.3
2.4
2.5and
2.6

### (a) Operator setting of problem: dynamics generator

Let us reformulate problem (2.1)–(2.6) as an evolution problem in a Hilbert space, which is a state space of the system.

First, we return to problem (2.1) and rewrite it in the form
2.7
where
2.8
and, without misunderstanding, we will use the same notation *E* for *E*/*ρ*.

We rewrite system (2.7) as the first order in the time-evolution equation for the following eight-component vector-valued function:
2.9
where the superscript ‘T’ stands for transposition. System (2.7) can be written in the form
2.10
with being given as the following matrix differential expression:
2.11
We assume that differential expression (2.11) is defined on the functions subject to the boundary conditions obtained from equations (2.2)–(2.6),
2.12
2.13
Here we have coupled beams subject to a four-parameter family of dynamical boundary conditions at the right end. Now we are in a position to give a Hilbert space setting. For each *t*, we assume that the eight-component vector of equation (2.10) belongs to a state space (see equation (2.15)). We consider the above vector-valued function of *x*
2.14
as an element of an energy space , which is a Hilbert space obtained as a closure of smooth functions in the following energy norm:
2.15
It is clear that the components of belong to the following Sobolev spaces:
In , we consider a matrix differential operator , where *α* and *β* are the two-component vectors, *α*=(*α*^{(1)},*α*^{(2)}), *β*=(*β*^{(1)},*β*^{(2)}), the operator defined by matrix differential expression (2.11) on the following domain:
2.16
In what follows, it is convenient to represent the operator as the sum, , where *L*_{αβ} coincides with if one sets . So both operators and *L*_{αβ} have the same domains. The operator is defined by the formula
2.17
where is the notation for a 4×4 matrix with the zero entries.

Our next result, shown in Shubov & Rojas-Arenaza (2008, 2010*a*), is the following statement.

### Lemma 2.2.

*given by formula (2.17) defines a bounded operator in* *. Therefore, the operator* *is a relatively compact perturbation of the operator L*_{αβ}*.*

### (b) Spectral properties of operator *L*_{αβ}

It is shown in our paper (Shubov & Rojas-Arenaza 2008) that *L*_{αβ} is an unbounded non-self-adjoint operator with compact resolvent . This means that the spectrum of *L*_{αβ} is purely discrete and each eigenvalue is normal. (We recall (Gohberg & Goldberg 1981; Gohberg & Krein 1996) that an eigenvalue *λ*_{0} of a bounded operator *A* in a Hilbert space *H* is normal if (i) the algebraic multiplicity of *λ*_{0} is finite and (ii) the space *H* breaks up into the direct sum of subspaces , where is the root subspace of the operator A corresponding to *λ*_{0}, and is an invariant subspace of A in which the operator is invertible.)

By lemma 2.2, the operator from formula (2.17) is relatively compact with respect to the operator *L*_{αβ}. (We recall this notion: let *A* be a linear operator having at least one regular point *λ*. An operator *B* is said to be compact relative to *A* if its domain contains and the operator is compact. Obviously, this definition does not depend on the choice of *λ* (Gohberg & Goldberg 1981; Marcus 1998).) Therefore, the sum has the same structure of the spectrum as *L*_{αβ} (i.e. the spectrum of consists of *normal* eigenvalues).

Our main result of Shubov & Rojas-Arenaza (2008) is the following statement.

### Theorem 2.3.

*Let L_{00} be a particular case of the operator L_{αβ} corresponding to α*

^{(1)}=

*α*

^{(2)}=

*β*

^{(1)}=

*β*

^{(2)}=0.

*Then both the operator*2.18

*L*_{αβ}and the operator*L*_{00}have compact inverses and the following relation holds:*where is a finite-dimensional operator in , whose rank is at most 4.*

*For , is given by the formula*
2.19

### Corollary 2.4.

—

*The operator is a relatively compact perturbation of the operator**L*_{αβ}, which means that has a purely discrete spectrum of normal eigenvalues.—

*The operator**L*_{αβ}is dissipative, i.e. for any and ℜ*α*^{(i)}≥0, ℜ*β*^{(i)}≥0,*i*=1,2. (We recall (Marcus 1988; Gohberg & Krein 1996 ) that a linear operator*A*with the domain is called dissipative if ℑ(*Af*,*f*)≥0 for—

*The operator**L*_{αβ}generates a strongly continuous semi-group of contractions, and generates a strongly continuous semi-group of bounded operators (Pazy 1983 ).

The rest of the paper is devoted to the derivation of the asymptotic formulae for the eigenvalues of the operator . As we know, the eigenvalues are directly related to the natural frequencies of the corresponding flexible structure (i.e. a double-walled nanotube). Having a model of the two coupled Timoshenko beams, we are dealing with a system involved in four different types of motion simultaneously. (A detailed asymptotic analysis of the spectrum for a single Timoshenko beam model with variable parameters can be found in Shubov (2002).) In what follows, we prove that the discrete spectrum splits asymptotically into four spectral branches and derived asymptotic representation for each spectral branch.

Our main result is the following statement.

### Theorem 2.5.

*Assume that the structural parameters are such that and Assume that the boundary parameters α^{(i)} and β^{(i)}, i=1,2, are such that the following conditions hold: λ_{1}≠α^{(2)}, β^{(1)}≠β^{(2)}*,
2.20

*The set of the eigenvalues of the operator asymptotically splits into four branches denoted by ,*.

*j*=1,…,4, . All four branches are located in a strip of the complex plane parallel to the real axis, i.e. and as*The following asymptotic approximation holds for the first spectral branch as*
2.21
*where, under the symbol ‘ln’, the principal value of the logarithm is understood.*

*The following asymptotic approximation holds for the second spectral branch as*
2.22
*The following asymptotic approximation holds for the third spectral branch as*
2.23

The following asymptotic approximation holds for the fourth spectral branch as 2.24

## 3. Left reflection matrix

### (a) Fundamental system of solutions of the spectral equation

In this section, we focus on the spectral equation 3.1 Rewriting equation (3.1) component-wise, we obtain the following system of four equations: 3.2

For the rest of §3 and in §4, we are interested in the asymptotic representation of the fundamental system of solutions of system (3.2) as . It is convenient to reserve the notation ‘*Ψ*’ for a solution of system (3.2) satisfying the boundary conditions as well, i.e. . Since we are interested in a solution for system (3.2), which does not necessarily satisfy the boundary conditions, we will use another notation, ‘*F*’. If *F*≡(*f*_{0},*f*_{1},*f*_{2},*f*_{3})^{T}, then system (3.2) can be written as the following matrix differential equation:
3.3
Introducing the following notation:
3.4
we obtain from equation (3.3) a system of two matrix equations with respect to the unknown vectors *F*^{(1)} and *F*^{(2)} defined in notation (3.4)
3.5
Denoting , and eliminating *F*^{(1)} from system (3.5), we transform this system to the desired form,
3.6
Let us look for a solution of equation (3.6) in the form *F*^{(2)}(*λ*,*x*)=*e*^{μ(λ)x}*F*_{0}, where *F*_{0} is a two-component vector with constant components. Substituting *F*^{(2)}(*λ*,*x*) into equation (3.6), we get
3.7
Equation (3.7) has a non-trivial solution if and only if the determinant of the matrix in equation (3.7) acting on *F*_{0} is zero. It is convenient to introduce new notation,
3.8

### Theorem 3.1.

*The matrix differential equation (3.6 ) has eight linearly independent solutions given by the following formulae:*
3.9
*with functions μ_{j}(λ) being given explicitly as*
3.10
3.11

*and*3.12

*where*

*δ*=

*δ*

_{1}+

*δ*

_{2}.

*The constant vectors are defined by the formulae*
3.13

In the sequel, we need asymptotic approximations for all eight roots as .

### Lemma 3.2.

*The following asymptotic approximations are valid for* *as*
3.14
3.15
3.16*and*
3.17

We do not provide the proof of this lemma since it is based on straightforward usage of standard methods of asymptotic analysis, and can be easily reconstructed.

### (b) Reflection matrices

The general solution of the differential equation (3.6) can be written in the form 3.18 where and are unknown coefficients.

If we substitute *F*_{g}(*λ*,⋅) into boundary conditions, we will have a linear system of eight equations in eight unknowns , , *n*=1,…,4. Let *K*(*λ*) be the 8×8 matrix of coefficients from the aforementioned system. We can write the set of boundary conditions in a matrix form, , where
3.19
Finding approximations for the solution of the equation is an extremely difficult problem. So we take an alternative approach, which we call the two-step procedure. This approach was instrumental in our study of the spectral asymptotics of a non-homogeneous Timoshenko beam model (Shubov 2002), coupled bending-torsion vibration model (Shubov 2001, 2005; Shubov & Peterson 2004), and a system of coupled Timoshenko beam models in Shubov & Rojas-Arenaza (2010*b*).

Let us introduce two four-component vectors and by the following formulae: 3.20

— Using the left-hand-side boundary conditions, we will find a relation between and in the form 3.21 We call

*the left-reflection matrix*. We notice that if the vectors and are related through*the left-reflection matrix*then the solution*F*_{g}(*λ*,⋅) of equation (3.6) satisfies only the left-end boundary conditions without any restrictions of its behaviour at the right end.— Using the right-hand-side boundary conditions, we will find a relation between and in the form 3.22 We call

*the right-reflection matrix*. We notice that if the vectors and are related through*the right-reflection matrix*then the solution*F*_{g}(*λ*,⋅) of equation (3.6) satisfies only the right-end boundary conditions without any restriction on its behaviour at the left end.— Having relations (3.21) and (3.22), we conclude that they are consistent if and only if 3.23

Provided that exists, we obtain from equation (3.23) that
3.24
It is clear that homogeneous equation (3.24) has non-trivial solutions only at those points *λ* that satisfy the following equation:
3.25
Equation (3.25) will be called the main spectral equation. It will be clear from the proof below that this equation has a countable set of solutions. Our goal is to find an asymptotic approximation for this set. Using our two-step procedure and working with the left- and right-reflection matrices, we are able to reduce asymptotic analysis of an 8×8 matrix to the analysis of 4×4 matrices. This makes our problem tractable.

To carry out our programme, we need asymptotic approximations for both reflection matrices. In this section, we will derive asymptotic approximation for the *left-reflection matrix* . However, we begin with one preliminary step.

### (c) New representation of boundary conditions

We recall that our initial point was to study the eigenvalue problem (3.1). As the first step, we have focused on the equation , disregarding the boundary conditions. We have shown that the above equation involving the 8×8 matrix differential expression (2.11) can be written as a system of two coupled differential equations (3.5) and each equation contains only 2×2 matrices. By eliminating *F*^{(1)} from system (3.5), we have arrived at the fourth-order differential equation for unknown vector-function *F*^{(2)}. Since our goal in §3 is to describe a solution of equation (3.6) satisfying the left-hand-side boundary conditions, we have to rewrite these conditions in terms of *F*^{(2)} only.

Let *G*_{α} and *G*_{β} be two matrices defined by
3.26
The boundary conditions at *x*=0 can be given in the form
3.27
and
3.28
The boundary conditions at *x*=*L* can be given in the form
3.29
and
3.30

### (d) Left-reflection matrix

As we know, the general solution *F*_{g}(*λ*,⋅) of differential equation (3.6) is given by formula (3.18). Let us substitute *F*_{g}(*λ*,⋅) into the boundary conditions (3.27) and (3.28) to have
3.31
and
3.32
It is convenient to introduce new notation (see equation (3.6) for )
3.33
With equation (3.33), system (3.31) and (3.32) can be written as the following matrix equation:
3.34
where the matrices and are given by the formulae
3.35
and
3.36
The 2×2 matrices and are defined as follows:
3.37
In terms of and , the *left-reflection matrix* introduced in relation (3.21) can be given by
3.38
Comparing formulae (3.35) and (3.36), we notice that
3.39
where
3.40
Using relation (3.39) and assuming that exists, we derive that
3.41
We point out here that the question of existence of the matrix has been addressed in remark 4.3 in Shubov & Rojas-Arenaza (2010*b*). Since the left-hand-side boundary conditions are similar in the present paper and in Shubov & Rojas-Arenaza (2010*b*), all the arguments of remark 4.3 of Shubov & Rojas-Arenaza (2010*b*) are valid for the present work as well.

### Theorem 3.3.

*The left-reflection matrix that is defined as exists for any and has the following asymptotic representation as*
3.42

In the present paper, we do not give any proof of theorem 3.3 since a detailed proof of a similar result is given in Shubov & Rojas-Arenaza (2010*b*). We recall that the present work differs from Shubov & Rojas-Arenaza (2010*b*) only with respect to the right-hand-side boundary conditions (i.e. with respect to the representations and properties of the right-reflection matrix). All the results on the right-reflection matrix are rigorously stated and proved in §4.

## 4. Right-reflection matrix

To evaluate the right-reflection matrix, we take the general solution of equation (3.6), *F*_{g}(*λ*,⋅), given by equation (3.18) and substitute into the right-hand-side boundary conditions (3.29) and (3.30). First, we modify equation (3.29). Since and *G*_{β} are diagonal matrices, this equation can be reduced to
4.1
Leaving all terms of the highest order in *λ* at the left-hand side of equation (4.1) and moving all terms of the lower order in *λ* to the right-hand side, we obtain the following equation:
4.2
We observe that each term of equation (4.2) is an exponential–polynomial function of *λ*; the order of each polynomial standing at the left-hand-side terms is 4, while the order of each polynomial standing at the right-hand-side terms is at most 2.

Similarly, rearranging terms in equation (3.30), we obtain the second condition in the form 4.3

### Remark 4.1.

As will be clear from calculations below, to evaluate the asymptotic representation for the right-reflection matrix, we can use the reduced equations obtained from equations (4.2) and (4.3) by replacing the right-hand sides with zeros. Such a reduction will change neither the leading asymptotic terms in the entries of the right-reflection matrix nor the orders of all remainder terms. Thus, in what follows, we deal with the reduced system
4.4
and
4.5
Substituting *F*_{g}(*λ*,⋅) of equation (3.18) into both equations (4.4) and (4.5) yields
4.6
and
4.7
Taking into account explicit formulae for all 2×2 matrices *A*, , *G*_{α} and *G*_{β} (see formulae (3.4) and (3.26)) and expressions for , *n*=1,…,4 (equation (3.13)), we can see that the system of two equations (4.6) and (4.7) is equivalent to the one matrix equation
4.8
where
4.9
and and are defined in formula (3.20). From equation (4.8), it follows that the right-reflection matrix can be represented in terms of as
4.10
In what follows, it is convenient to use new notation for the matrix ,
4.11
The main result of this paragraph, theorem 4.2, deals with the asymptotic approximation for the matrix , as . Its proof can be split up into the proofs of five technical lemmas. These lemmas with the detailed proofs can be found in the electronic supplementary material for the paper. To present the main results, we introduce new notation. Let , and be functions defined by the formulae
4.12
Before proving theorem 4.2, we make the following remark. It turns out that the proof of the main result depends on whether the boundary parameters are critical or non-critical. We say that one of the parameters, *α*^{(i)} or *β*^{(i)}, *i*=1,2, has a critical value if one of the following conditions holds (*γ* is from (3.8)):

The proof of the results for non-critical values of the boundary parameters is technically more complicated. The results of theorem 4.2 and of §5 are proved for non-critical cases. The adjustment of the proofs to the case of the critical values of the boundary parameters is addressed in remark 5.2 at the end of §5.

### Theorem 4.2.

*Let the boundary parameters satisfy the following set of conditions*
4.13
*In addition, there exists an absolute constant C_{0}>0 such that the structural and boundary parameters satisfy the estimates*
4.14

*where*4.15

*γ*is defined in (3.8) and*C*_{0}is an absolute constant whose precise value is immaterial for us. Then the entries of the matrix from (4.11) can be approximated for as follows:*where the diagonal entries are given by the formulae*4.16

*The non-diagonal entries are given by the formulae*4.17

*The notation*.

*ω*_{ij}(*λ*) in formulas (4.16) and (4.17) stands for the factor (1+*O*(*λ*^{−2})) occurring on the intersection of the ith row and jth column

### Proof.

Taking into account that *α*^{(1)}≠*α*^{(2)}, *β*^{(1)}≠*β*^{(2)} and the results of lemmas 4.3, 4.4 and 4.7 from the electronic supplementary material, we obtain the following approximations as for all entries of the matrix :
4.18
Rewriting formulae (4.18) in terms of the notation introduced in formula (4.12) and using conditions (4.13) and (4.14), we arrive at the desired result. ■

The theorem is shown.

### Remark 4.3.

If we consider a particular case of the parameters when *α*^{(1)}= *α*^{(2)} and *β*^{(1)}=*β*^{(2)}, then the approximation of the matrix obtained in theorem 4.2 coincides with the approximation for the matrix obtained in Shubov & Rojas-Arenaza (2010*b*).

## 5. Spectral asymptotics

We recall that *the main spectral equation*, equation (3.25), yields the asymptotic distribution of the eigenvalues of the problem. Since , we rewrite this equation in the equivalent form
5.1
where the asymptotic representations for the matrices and are given by formulae (3.42) and ((0.19) in the electronic supplementary material) and theorem 4.2, respectively, (for , see formulae (4.16)–(4.18)). Let us denote the exponential factors of the matrix from equation (4.9) by
5.2
It is convenient to rewrite the approximation for the matrix using the following set of notation:
5.3
and
5.4
Then the matrix can be written in the form (see formulae (4.16)–(4.18))
5.5

Combining equations (3.42), (5.2) and (5.5), we rewrite equation (5.1) as follows:
5.6
The factor has a meaning similar to the meaning of the factor *ω*_{ij}(*λ*) introduced in equations (4.16) and (4.17). We already know that equation (5.6) has a countable set of roots, each of finite multiplicity. Let , , be the aforementioned set of roots counting multiplicities.

Our first result is the following statement.

### Lemma 5.1.

*The set of roots of equation (*5.6*) is located in a strip on the complex plane* *parallel to the real axis*
5.7
*where C*_{0} *is a positive constant.*

We do not give the proof since it is similar to the proof of lemma 6.5 from Shubov & Rojas-Arenaza (2010*b*).

Now we are in a position to prove the main result of the paper, theorem 2.5.

### Proof of theorem 2.5.

Let us rewrite matrix (5.6) using new notation. Let be the leading terms of the diagonal entries of the matrix from equation (5.6), i.e.
5.8
Using lemma 5.1 and notation (5.8), we can represent equation (5.6) in the form
5.9
where . Consider the ‘model equation’ obtained from equation (5.9) by omitting the term *O*(*λ*^{−2}) and have
5.10
Evaluating both determinants in equation (5.10), we obtain the following equation:
5.11
Based on notation (5.2), equation (5.11) yields the following system of two equations:
5.12

Let us consider the first equation from system (5.12) and rewrite it using approximations (3.14) and (3.16) for *μ*_{1,3}(*λ*) as . We have
5.13
If we denote
5.14
then omitting the asymptotically small term, we obtain the following quadratic equation:
5.15
Taking into account the formula for *η*_{α} from notation (5.4), we simplify equation (5.15) and obtain a new equation
5.16
whose roots, and , are (see formulae (5.3))
5.17
Each equation in (5.17) has a countable set of roots that can be given by explicit formulae. Let us take, for example, the equation for , then the set of its roots is
5.18
If instead of the ‘model equation’ (5.16), we examine the complete equation (5.13) containing the asymptotically small term *O*(*n*^{−1}), then, by using Rouche’s theorem (Jeffrey 2006), we derive the approximation for the first spectral branch eigenvalues as follows:
5.19
Formula (5.19) coincides with formula (2.21) of theorem 2.5. We do not present all estimates necessary for the application of Rouche’s theorem since the detailed derivation of similar estimates can be found in Shubov & Rojas-Arenaza (2010*b*). Formula (2.22) of theorem 2.5 can be obtained from the second equation from (5.17).

Finally, the second equation of system (5.12) yields formulae (2.23) and (2.24) for the approximations of the third and fourth spectral branches, respectively.

The theorem is completely shown.

### Remark 5.2.

Now we present modifications of the proofs for the case of a critical value boundary parameter. Assume that one of the conditions from estimates (4.14) fails, i.e. assume that and the remaining three conditions are valid. In this case, called a critical case of the boundary parameters, the diagonal entries of matrix (5.5) can be given by the formulae 5.20 with the entries and being unchanged. As a result, the second equation of system (5.12) is the same as before and the first equation has a new form 5.21 Taking into account approximation (3.16), we reduce equation (5.21) to an equation similar to equation (5.13), i.e. 5.22 Using notation (5.14), this equation can be written as 5.23 This equation has two distinct roots: and . Obviously, the first root generates the asymptotic representation for the first spectral branch given precisely by formula (2.21) of theorem 2.5. The second root yields the formula for the second spectral branch corresponding to the case . It can be easily seen that the second spectral branches corresponding to the critical and non-critical cases of the boundary parameters are in agreement, i.e. the asymptotic representation for the eigenvalues of the third branch for the critical value case can be obtained from formula (2.23) by replacing with .

Completing similar analysis for the remaining three cases of the critical values of the boundary parameters, we conclude that the results of theorem 2.5 are valid regardless of whether the boundary parameters are critical or non-critical.

The proof is complete. ■

## 6. Concluding remarks and future research

In this paper, we presented a set of formulae describing the behaviour of vibrational frequencies for a DWNT model. The model has recently appeared in multiple research papers (see Jakobson *et al.* 1996; Ru 2001; Qian *et al.* 2002; Pantano *et al.* 2003; Yoon *et al.* 2003; Wang *et al.* 2006; Xu *et al.* 2006; Gibson *et al.* 2007 and references therein). However, to the best of our knowledge, rigorous asymptotic and spectral analysis resulting in formulae similar to formulae (2.21)–(2.24) had not been accomplished.

— Regarding this direction of research, we plan to perform a numerical analysis of the system to check whether the numerical approach yields a similar four-branch spectrum. In addition, we will continue our investigation to address questions concerning: (i) properties of the mode shapes such as minimality, completeness and the Riesz basis property in the energy space , (ii) controllability of the model using the boundary parameters

*α*^{(i)},*β*^{(i)},*i*=1,2, as the control parameters, and (iii) stability and the pole assignment problems.— Examining our analytical formulae, we can address an issue that has been discussed in multiple sources on nanotubes. Namely, there is a statement indicating that nanotube vibration goes on ‘frictionless’, without energy dissipation. The following significant question arises: Do mathematical formulae (2.21)–(2.24) reflect this ‘frictionless’ motion? The answer is most likely affirmative if one pays attention to the fact that there is no van der Waals force strength in the leading asymptotic terms. Certainly, this constant ‘’ is present in the remainder terms

*O*(|*n*|^{−1}), which means that even though the problem is not conservative due to the existence of the van der Waals interaction, this interaction is not strong enough to affect the leading asymptotic terms describing the vibrational spectrum.— As is well known for any mechanical system of macro-size, it is important to investigate a couple of dozen of the lower vibrational frequencies and mode shapes since higher energy modes are heavily damped. Contrary to this case, for a nanoscale system, owing to the aforemention ‘frictionless’ motion, one needs to know hundreds or even thousands of the lower vibrational frequencies. This is exactly the point where all the advantages of a continuous model become obvious.

— An extremely important research direction will be to extend the model to nonlinear cases both for van der Waals force nonlinearity and for structural nonlinearities.

— Finally, developing a rigorous model of DWCT combining a Timoshenko beam model for the inner tube and an elastic shell model for the outer tube and performing its mathematical analysis will be a challenging step in our programme.

## Acknowledgements

Partial support by the National Science Foundation grant DMS no. 0604842 is highly appreciated by the first author.

- Received March 29, 2010.
- Accepted May 26, 2010.

- This journal is © 2010 The Royal Society