Phonons in chiral nanorods and nanotubes: A Cosserat-rod-based continuum approach

Abstract

A Cosserat-rod-based continuum approach is presented to obtain phonon dispersion curves of flexural, torsional, longitudinal, shearing, and radial breathing modes in chiral nanorods and nanotubes. Upon substituting the continuum wave form in the linearized dynamic equations of stretched and twisted Cosserat rods, we obtain an analytical expression of a coefficient matrix (in terms of the rod’s stiffnesses, induced axial force, and twisting moment) whose eigenvalues and eigenvectors give us frequencies and mode shapes, respectively, for each of the above phonon modes. We show that, unlike the case of achiral tubes, these phonon modes are intricately coupled in chiral tubes owing to extension–torsion–inflation and bending–shear couplings inherent in them. This coupling renders the conventional approach of obtaining stiffnesses from the long wavelength limit slope of dispersion curves redundant. However, upon substituting the frequencies and mode shapes (obtained independently from phonon dispersion molecular data) in the eigenvalue–eigenvector equation of the above-mentioned coefficient matrix, we are able to obtain all the stiffnesses (bending, twisting, stretching, shearing, and all coupling stiffnesses corresponding to extension–torsion, extension–inflation, torsion–inflation, and bending–shear couplings) of chiral nanotubes. Finally, we show unusual effects of the single-walled carbon nanotube’s chirality as well as stretching and twisting of the nanotube on its phonon dispersion curves obtained from the molecular approach. These unusual effects are accurately reproduced in our continuum formulation.

Keywords

Cosserat rod Phonon Chirality single-walled nanotube

1. Introduction

Phonons are linearized eigenmodes of atomic-level vibration in crystalline solids. They are used to determine thermomechanical and electrical properties of crystals as well as their stability [1 –5]. Phonon dispersion curves are typically obtained by using the translational periodicity inherent in crystals to block diagonalize the large Hessian matrix associated with the crystal’s supercell [1, 2]. Each of the smaller block matrices then yields a set of phonons corresponding to a particular wave vector. As the periodic unit cell of chiral nanotubes is usually very large, the method generates significantly large number of dispersion curves in phonon spectra. However, chiral nanotubes as well as uniformly bent and twisted achiral nanotubes possess much smaller helically repeating cell [6, 7]. Several researchers [8 –12] have used this smaller helically repeating cell and helical isometry to block diagonalize the large Hessian matrix. This drastically reduces the number of phonon dispersion curves per wave vector and also reduces the computational cost of obtaining them significantly.¹ There are two types of phonon modes for bulk crystals: acoustic and optical. The acoustic modes also have their continuum analogues and their frequencies approach zero in the long wavelength limit. The frequencies of optical modes, on the other hand, approach a finite value in the long wavelength limit. The optical phonons do not have any analogue in usual continuum theories. In the case of nanorods and nanotubes, however, both the radial breathing and shearing phonon modes have continuum analogues but their frequencies do not approach zero in the long wavelength limit [13 –15]. In fact, when helical symmetry is used in block diagonalization, the flexural phonon mode’s frequency also does not vanish at zero wave vector but at a non-zero wave vector [12].

In [16], Graff derived elastic wave solutions for longitudinal, torsional, shearing, and bending modes of vibration in isotropic rods about their stress-free configuration. In [17], Rego and Kirczenow substituted elastic constants for square GaAs nanowires in these wave solutions to obtain their phonon spectra corresponding to the above-mentioned modes. Conversely, the phonon modes having continuum analogues are also used in estimating elastic stiffnesses of nanostructures. Popov et al. [18], for example, evaluated Young’s modulus and shear modulus for achiral as well as chiral SWCNTs from the slopes of longitudinal and torsional dispersion curves in long wavelength limit. However, we show later that their calculation for chiral nanotubes is incorrect because these slopes are no longer proportional to elastic moduli owing to extension–torsion coupling present in chiral tubes. Several authors have also obtained bending stiffness of stress-free achiral nanotubes from the curvature of flexural dispersion curve [19 –22]. We again show that owing to non-zero bending–shear coupling in chiral nanotubes, the bending stiffness cannot be obtained directly from this curvature.

One can also obtain stiffnesses by directly imposing deformation. For example, one can bend isolated segment of a chiral nanotube and clamp the end atoms in a bent state. The bending stiffness can then be obtained by comparing the induced deformation or energy with classical beam theories [23 –25]. However, such methods not only require the nanotube to be long so that they behave as a continuum beam, the method also suffers from end effects as demonstrated in [24]. James [6] proposed a theory of objective structures using which one can deform achiral nanorods or nanotubes in pure bending mode without end effects and subsequently obtain their bending stiffness. However, their method cannot be employed to bend chiral nanostructures. Their method cannot be used to obtain shearing stiffness of either achiral or chiral nanostructures. Kumar et al. [7] recently proposed the helical Cauchy–Born (HCB) rule to obtain all the stiffnesses of nanorods and nanotubes modeled as special Cosserat rods. The method requires imposing uniform strain field (along the arc-length) on nanorods and nanotubes in order to obtain their stiffnesses. However, their method also cannot be applied to deform chiral nanotubes in pure bending or pure shearing mode.

To the best of the authors’ knowledge, no accurate technique exists to obtain bending, shearing, and several coupling stiffnesses of chiral nanotubes. Furthermore, the continuum wave solutions for longitudinal, torsional, flexural, shearing, and radial breathing modes in chiral rods taking into account the inherent extension–torsion–inflation and bending–shear couplings have not been obtained yet. Similarly, a continuum investigation into the effect of finitely stretching and twisting of elastic rods on these modes is also missing in the literature. There has been very limited work in this regard. For example, in [26] Aghaei et al. showed the first-order influence of stretching a nanotube on its flexural phonon dispersion curve. Chang [27] studied the effect of axial strain and twist on radial breathing frequency of SWCNTs but from purely atomistic perspective. The goal of this paper is to address the above-mentioned limitations.

The paper is organized as follows. In Section 2, we briefly discuss the theory of special Cosserat rods. In Section 3, we then discuss continuum wave solutions of longitudinal, torsional, shearing, and bending modes in stretched and twisted chiral rods and obtain a generalized eigenvalue–eigenvector equation whose solution gives us frequencies and mode shapes of these wave modes. We also obtain analytical formulas of frequencies of these modes and discuss how coupling stiffnesses as well as the induced axial force and twisting moment affect these frequencies. In Section 4, we derive continuum wave solution for longitudinal, torsional, and radial breathing modes in chiral elastic tubes using the general theory of Cosserat rods and discuss how the three modes are coupled to each other. In Section 5, we then present a methodology to obtain all the stiffnesses of chiral nanotubes by matching both frequency and mode shape information of various phonon modes obtained directly from molecular approach with the ones obtained using our continuum wave solution. In Section 6, we then obtain stiffnesses of both achiral and chiral SWCNTs using the methodology described in Section 5. In Section 7, we show interesting effects of stretching and twisting a SWCNT on its longitudinal, torsional, radial breathing, and bending phonon dispersion curves and further show that these effects are replicated accurately using our continuum approach. Finally, Section 8 concludes our paper.

Notation

The set of unit vectors ${e_{1}, e_{2}, e_{3}}$ represents a fixed, right-handed, orthonormal basis for $R^{3}$ . We further use $X = X_{i} e_{i}$ (with $X_{3} = s$ ) to denote the position vector of a typical material point of a rod in its reference configuration (Figure 1). The components $(X_{1}, X_{2})$ denote the coordinates in the rod’s undeformed cross-section ( $Ω_{0}$ ) and s denotes the coordinate along the undeformed arc-length. Repeated Latin indices always sum from 1 to 3 whereas repeated Greek indices sum from 1 to 2. The differentiation of a field variable with respect to the undeformed arclength is denoted by $\frac{d}{ds} \equiv ()'$ while that with respect to time is denoted by $\frac{d}{dt} \equiv \overset{\cdot}{(})$ .

Figure 1.

A typical rod deforming from its straight state reference configuration.

2. Brief overview of the theory of special Cosserat rods

We begin with the kinematics of special Cosserat rods.² Figure 1 shows a typical rod deforming from its straight state reference configuration. The variable $r (s)$ denotes the centerline of the deformed rod. We also have two directors ( $d_{1} (s)$ , $d_{2} (s)$ ) which are the images of two orthonormal basis vectors ( $e_{1}$ , $e_{2}$ ) in the undeformed cross-section of the rod. The third director $d_{3} (s)$ is assumed to be perpendicular to the other two. The directors rigidly rotate in this theory, i.e., $d_{i} (s) = R (s) e_{i}$ . Here, $R (s)$ denotes a three-dimensional rotation tensor and models orientation of the rod’s cross-section. The field quantities ( $r, R$ ) are thus the kinematic variables in this theory. We further introduce the following quantities:

\underline{v} = R^{T} r' = ν_{i} e_{i},

K = R^{T} R', \underline{k} = axial (K) = κ_{i} e_{i} .

(1)

The first two components of $\underline{v}$ represent shear while the third component represents axial stretch. Similarly, the first two components of $\underline{k}$ represent the components of local bending curvature while the third component represents twist. We then let $n$ and $m$ denote the internal contact force and internal moment, respectively, that act on a typical cross-section of the rod and write

n = n_{i} d_{i}, m = m_{i} d_{i} .

(2)

Here, ( $n_{1}, n_{2}$ ) are the shear forces, $n_{3}$ is the axial force, ( $m_{1}, m_{2}$ ) are the bending moments, and $m_{3}$ is the torque or twisting moment. Upon linear and angular momentum balance of an arbitrary segment of the rod, one obtains the following equations:

n' + \hat{n} = ρ_{0} A \overset{\cdot\cdot}{r},

(3a)

m' + r' \times n + \hat{m} = \overset{\cdot}{\bar{I ω}} = I α + ω \times [I ω] .

(3b)

Here, $\hat{n}$ and $\hat{m}$ are the distributed force and distributed couple, respectively, that act on the rod, $ρ_{0}$ denotes mass per unit undeformed volume while A denotes the undeformed cross-sectional area. Furthermore, $\overset{\cdot\cdot}{r}$ denotes acceleration of the mass center of the rod’s cross-section, $ω$ denotes angular velocity and $α = \overset{\cdot}{ω}$ denotes angular acceleration of the rod’s cross-section. Finally, $I = R I_{0} R^{T}$ denotes the moment of inertia tensor of the rod’s rotated cross-section. Here, $I_{0}$ denotes the same for the unrotated cross-section and is defined as

I_{0} = \sum_{j = 1}^{3} I_{j} e_{j} \otimes e_{j} where

I_{1} = \int_{Ω} ρ_{0} X_{2}^{2} d X_{1} d X_{2}, I_{2} = \int_{Ω} ρ_{0} X_{1}^{2} d X_{1} d X_{2}, I_{3} = I_{1} + I_{2} .

(4)

We have assumed the above diagonal form for $I_{0}$ because the axes $(e_{1}, e_{2})$ are assumed to be aligned along the cross-section’s principal axes of inertia.

For a hyperelastic rod, we further assume the existence of a twice-differentiable, scalar-valued function $\hat{Φ} (\underline{v}, \underline{k})$ which is strain energy per unit undeformed length in this theory. One can show that

n_{i} = \frac{\partial \hat{Φ}}{\partial ν_{i}}, m_{i} = \frac{\partial \hat{Φ}}{\partial κ_{i}} .

(5)

Using the symmetry group relevant for circular rods, Healey [29] classifies such rods into transversely hemitropic and transversely isotropic rods. The relevant symmetry group for hemitropic rods is as follows:

G \equiv {Q_{θ} : \hat{Φ} (Q_{θ} \underline{v}, Q_{θ} \underline{k}) = \hat{Φ} (\underline{v}, \underline{k}), \forall 0 \leq < 2 π, Q_{θ} \in SO (3), Q_{θ} e_{3} = e_{3}} .

(6)

Here, $θ$ denotes the angle of rotation of the rotation tensor $Q$ . Th group $G$ can make a distinction between left- and right-handed helical patterns in a rod and is useful in modeling intrinsically twisted rods, e.g., chiral nanotubes, DNA, collagen fibrils etc. Let us define the following stiffness quantities:³

\hat{A} = \frac{\partial^{2} \hat{Φ}}{\partial κ_{α} \partial κ_{α}}, \hat{B} = \frac{\partial^{2} \hat{Φ}}{\partial κ_{3}^{2}}, \hat{C} = \frac{\partial^{2} \hat{Φ}}{\partial ν_{α} \partial ν_{α}}, \hat{D} = \frac{\partial^{2} \hat{Φ}}{\partial ν_{3}^{2}}, \hat{E} = \frac{\partial^{2} \hat{Φ}}{\partial ν_{3} \partial κ_{3}}, \hat{F} = \frac{\partial^{2} \hat{Φ}}{\partial ν_{α} \partial κ_{α}} .

(7)

The coefficients $\hat{A}, \hat{B}, \hat{C}$ , and $\hat{D}$ denote the bending, twisting, shearing, and extensional stiffnesses, respectively. The coefficient $\hat{E}$ denotes coupling between axial strain and twist whereas $\hat{F}$ denotes coupling between bending and shear. Using the symmetry group (6) and certain “flip-symmetry” that models homogeneity in a rod along its arc-length, one can show that only the coefficients defined in (7) are non-zero in the rod’s stiffness tensor $\frac{\partial^{2} \hat{Φ}}{\partial [\underline{v}, \underline{k}] \partial [\underline{v}, \underline{k}]}$ when evaluated in the stress-free state (see [29] for proof). Accordingly, the most general quadratic form of $\hat{Φ}$ for hemitropic rods becomes

{\hat{Φ}}_{quad} (\cdot) = \frac{1}{2} [{\hat{A}}_{0} κ_{α} κ_{α} + {\hat{B}}_{0} κ_{3}^{2} + {\hat{C}}_{0} ν_{α} ν_{α} + {\hat{D}}_{0} (ν_{3} - 1)^{2} + 2 {\hat{E}}_{0} (ν_{3} - 1) κ_{3} + 2 {\hat{F}}_{0} ν_{α} κ_{α}] .

(8)

Here, the quantities ( ${\hat{A}}_{0}, {\hat{B}}_{0}, {\hat{C}}_{0}, {\hat{D}}_{0}, {\hat{E}}_{0}, {\hat{F}}_{0}$ ) denote the derivatives in (7) evaluated in the stress-free state. The coupling coefficients ( ${\hat{E}}_{0}, {\hat{F}}_{0}$ ) vanish in the case of transversely isotropic rods since their symmetry group also contains improper rotations in addition to (6).

3. Linearized wave solution in stretched and twisted rods

In this section, we obtain linearized wave solution in infinitely long elastic rods assuming no distributed load or distributed couple acts on the rod. We first linearize the dynamic equations of motion (3) about a uniformly stretched and twisted static configuration of the rod. Let ( $r_{ϵ}, R_{ϵ}$ ) denote an $ϵ$ -perturbed configuration of the rod, i.e.,

r_{ϵ} (s, t) = r (s) + ϵ Δ r (s, t),

R_{ϵ} (s, t) = \exp (ϵ Δ Θ (s, t)) R (s) .

(9)

Here, ( $r, R$ ) denote the configuration variables in the rod’s uniformly stretched and twisted state, i.e.,

r = s ν_{3} e_{3}, R = Q_{(θ = s κ_{3})} .

(10)

The quantity $Δ r (s, t)$ in (9) denotes a smooth, admissible perturbation in the rod’s centerline while $Δ Θ (s, t)$ denotes a smooth, admissible skew-symmetric tensor which is used to perturb $R$ . The axial vector of $Δ Θ (s, t)$ is denoted by $Δ θ (s, t)$ . Using the definition (1), the strains then get perturbed as follows:

{\underline{v}}_{ϵ} = R_{ϵ}^{T} r_{ϵ}' = \underline{v} + ϵ Δ \underline{v} + O (ϵ^{2}) where Δ \underline{v} = R^{T} (Δ r' + r' \times Δ θ),

{\underline{k}}_{ϵ} = axial (R_{ϵ}^{T} R_{ϵ}') = \underline{k} + ϵ Δ \underline{k} + O (ϵ^{2}) where Δ \underline{k} = R^{T} Δ θ' .

(11)

The linearization of equation (3a) upon imposing uniformity of strains ( $\underline{v}' = 0, \underline{k}' = 0$ ) then yields the following:

\begin{matrix} Δ Θ R (\underline{k} \times \frac{\partial \hat{Φ}}{\partial \underline{v}}) + R (Δ \underline{k} \times \frac{\partial \hat{Φ}}{\partial \underline{v}}) + R (\underline{k} \times (\frac{\partial^{2} \hat{Φ}}{\partial {\underline{v}}^{2}} Δ \underline{v} + \frac{\partial^{2} \hat{Φ}}{\partial \underline{v} \partial \underline{k}} Δ \underline{k})) \\ + R (\frac{\partial^{2} \hat{Φ}}{\partial {\underline{v}}^{2}} Δ \underline{v}' + \frac{\partial^{2} \hat{Φ}}{\partial \underline{v} \partial \underline{k}} Δ \underline{k'}) = ρ_{0} A Δ \overset{\cdot\cdot}{r} . \end{matrix}

(12)

while the linearization of equation (3b) yields

\begin{matrix} Δ Θ R (\underline{k} \times \frac{\partial \hat{Φ}}{\partial \underline{k}}) + R (Δ \underline{k} \times \frac{\partial \hat{Φ}}{\partial \underline{k}}) + R (\underline{k} \times (\frac{\partial^{2} \hat{Φ}}{\partial {\underline{k}}^{2}} Δ \underline{k} + \frac{\partial^{2} \hat{Φ}}{\partial \underline{v} \partial \underline{k}} Δ \underline{v})) + Δ r' \times R \frac{\partial \hat{Φ}}{\partial \underline{v}} \\ + R [\frac{\partial^{2} \hat{Φ}}{\partial \underline{v} \partial \underline{k}} Δ \underline{v}' + \frac{\partial^{2} \hat{Φ}}{\partial {\underline{k}}^{2}} Δ \underline{k}'] + r' \times [Δ Θ R \frac{\partial \hat{Φ}}{\partial \underline{v}} + R (\frac{\partial^{2} \hat{Φ}}{\partial {\underline{v}}^{2}} Δ \underline{v} + \frac{\partial^{2} \hat{Φ}}{\partial \underline{v} \partial \underline{k}} Δ \underline{k})] = R I_{0} R^{T} Δ \overset{\cdot\cdot}{θ} . \end{matrix}

(13)

We then assume the following wave forms for ( $Δ r, Δ θ$ ):

Δ r = e^{i (ω t - ks)} Δ {\hat{r}}_{i} e_{i}, Δ θ = e^{i (ω t - ks)} Δ {\hat{θ}}_{i} e_{i} .

(14)

Here, $i = \sqrt{- 1}$ , k represents wavenumber and $ω$ denotes angular frequency. The quantities ( $Δ {\hat{r}}_{1}, Δ {\hat{r}}_{2}$ ) denote the amplitudes of displacement along transverse directions $e_{1}$ and $e_{2}$ , respectively, whereas $Δ {\hat{r}}_{3}$ denotes the axial displacement amplitude. Similarly, $Δ {\hat{θ}}_{1}$ and $Δ {\hat{θ}}_{2}$ denote the amplitudes of rotation about $e_{1}$ and $e_{2}$ axes, respectively, whereas $Δ {\hat{θ}}_{3}$ denotes the same about the tube’s axis. As we consider wave solution in stretched and twisted rods, we assume that the stiffnesses in (7) are the only non-zero terms. Upon substituting the wave forms (14) in the linearized dynamic equations (12) and (13), we finally obtain the following generalized eigenvalue–eigenvector equation for the unknown amplitudes and frequency:

Cx = - ω^{2} D x

(15)

where

C_{6 \times 6} = [\begin{matrix} A_{4 \times 4} & O_{4 \times 2} \\ O_{2 \times 4} & B_{2 \times 2} \end{matrix}], x = [\begin{matrix} Δ {\hat{r}}_{1} \\ Δ {\hat{θ}}_{2} \\ Δ {\hat{r}}_{2} \\ Δ {\hat{θ}}_{1} \\ Δ {\hat{r}}_{3} \\ Δ {\hat{θ}}_{3} \end{matrix}], D_{6 \times 6} = diag [\begin{matrix} ρ_{0} A \\ I_{1} \\ ρ_{0} A \\ I_{1} \\ ρ_{0} A \\ I_{3} \end{matrix}],

\begin{matrix} A = (\begin{matrix} - k^{2} \hat{C} & - i k (N - \hat{C} ν_{3}) & 0 & - k^{2} \hat{F} \\ i k (N - \hat{C} ν_{3}) & (- k^{2} \hat{A} + ν_{3} (N - \hat{C} ν_{3})) & - k^{2} \hat{F} & i k (M - 2 ν_{3} \hat{F}) \\ 0 & - k^{2} \hat{F} & - k^{2} \hat{C} & i k (N - ν_{3} \hat{C}) \\ - k^{2} \hat{F} & - i k (M - 2 ν_{3} \hat{F}) & - i k (N - ν_{3} \hat{C}) & - k^{2} \hat{A} + ν_{3} (N - ν_{3} \hat{C}) \end{matrix}), \\ B = (\begin{matrix} - k^{2} \hat{D} & - k^{2} \hat{E} \\ - k^{2} \hat{E} & - k^{2} \hat{B} \end{matrix}) . \end{matrix}

(16)

Here, $O$ is zero submatrix. In submatrix $A$ , $N = \frac{\partial \hat{Φ}}{\partial ν_{3}}$ and $M = \frac{\partial \hat{Φ}}{\partial κ_{3}}$ are the axial force and twisting moment, respectively, that get induced in the rod. Owing to the structure of matrix $C$ , two of the eigenmodes corresponding to submatrix $B$ get decoupled from the four eigenmodes corresponding to submatrix $A$ . The former generates longitudinal and torsional modes whereas the latter generates two bending and two shearing modes. For each of the six modes, the associated eigenvector defines the corresponding mode shape. The generalized eigenvalue problem admits simple analytical solution in special cases which we now discuss.

3.1. Longitudinal and torsional mode frequencies

The two eigenvalues corresponding to submatrix $B$ in (16) yields the following:

{(\frac{ω}{k})}^{2} = \frac{1}{2} [\frac{\hat{D}}{ρ_{0} A} + \frac{\hat{B}}{I_{3}} \pm \sqrt{{(\frac{\hat{D}}{ρ_{0} A} - \frac{\hat{B}}{I_{3}})}^{2} + \frac{4 {\hat{E}}^{2}}{ρ_{0} A I_{3}}}] .

(17)

It implies that the ratio $\frac{ω}{k}$ is independent of the wavenumber k for longitudinal and torsional modes. Accordingly, their wave speeds are independent of the wavenumber. Furthermore, owing to the presence of non-zero extension–torsion coupling stiffness $\hat{E}$ , the two modes get coupled to each other. One also cannot obtain the three stiffnesses $(\hat{D}, \hat{B}, \hat{E})$ from just two wave speed information. We show later that one also needs to use the mode shape information to obtain all the stiffnesses. In case $\hat{E}$ is zero which is typically the case for non-twisted achiral rods, the two modes decouple and the formulas for the frequencies simplify as follows:

{(\frac{ω}{k})}_{l} = \sqrt{\frac{\hat{D}}{ρ_{0} A}}, {(\frac{ω}{k})}_{t} = \sqrt{\frac{\hat{B}}{I_{3}}} .

(18)

3.2. Bending and shearing mode frequencies

The formulas for bending ( $ω_{b}$ ) and shearing mode ( $ω_{s}$ ) frequencies turn out to be very long. However, we show now their Taylor-expanded formula about long wavelength limit, i.e., $k = 0$ :

\begin{matrix} ω_{b}^{2} = \frac{N}{ρ_{0} A ν_{3}} k^{2} + \frac{M}{ρ_{0} A ν_{3}^{2}} k^{3} + \frac{1}{ρ_{0}^{2} A^{2} ν_{3}^{3} (ν_{3} \hat{C} - N)} (- ρ_{0} A M^{2} + I_{1} N^{2} \\ + ν_{3} (ρ_{0} A (2 M \hat{F} - N \hat{A}) + ρ_{0} A ν_{3} (\hat{A} \hat{C} - {\hat{F}}^{2}) - I_{1} N \hat{C})) k^{4} + \dots . \end{matrix}

(19)

We can note the influence of axial force and twisting moment on the bending frequency. Owing to the presence of $N$ , the bending dispersion curve gets non-zero slope in the long wavelength limit, i.e.,⁴

\frac{ω_{b}^{2}}{k^{2}} |_{k \to 0} = \frac{N}{ρ_{0} A ν_{3}} .

(20)

We also note that the twisting moment $M$ generates the $k^{3}$ term in formula (19) which leads to asymmetry in the bending dispersion curve. However, in the tube’s stress-free state when $N$ and $M$ both vanish, the formula (19) simplifies to

ω_{b}^{2} = \frac{(\hat{A} \hat{C} - {\hat{F}}^{2})}{ρ_{0} A \hat{C}} k^{4} + \dots .

(21)

Furthermore, in the case of achiral tubes, the bending–shear coupling stiffness $\hat{F}$ is zero and the bending stiffness $\hat{A}$ can then be obtained directly from the curvature of this dispersion curve. For the general case as in (19), deducing bending stiffness from the dispersion curve is not straightforward, which we shall address later.

Finally, we Taylor-expand the shearing mode frequency about the long wavelength limit and obtain

ω_{s}^{2} = \frac{ν_{3}}{I_{1}} (ν_{3} \hat{C} - N) + \frac{1}{I_{1}} (M - 2 ν_{3} \hat{F}) k - \frac{1}{ρ_{0} A ν_{3}^{2}} M k^{3} + \dots .

(22)

We note that in order to obtain the shearing stiffness $\hat{C}$ , one can first obtain $N$ using (20) and then use the following formula:

\hat{C} = \frac{1}{ν_{3}} (\frac{I_{1}}{ν_{3}} ω_{s}^{2} |_{k = 0} + N) .

(23)

Another observation to note from formula (22) is that the shearing dispersion curve has non-zero slope at $k = 0$ . However, this slope vanishes in the stress-free state of achiral tubes.

4. Continuum wave solution for radial breathing mode

In order to understand the radial breathing mode from continuum perspective, we need to explicitly account for the tube’s inflation in the continuum model. As the cross-sectional deformation is not an independent measure in the special Cosserat rod model, we use the general Cosserat rod model [30] for this purpose and restrict ourselves to its coupled extension–torsion–inflation deformation. Kumar and Mukherjee [30] derived the quadratic form of strain energy for this model whose restriction to extension–torsion–inflation deformation is as follows:

{\tilde{Φ}}_{quad} (ν_{3}, κ_{3}, β, β') = \frac{1}{2} ({\tilde{B}}_{0} κ_{3}^{2} + {\tilde{D}}_{0} {(ν_{3} - 1)}^{2} + {\tilde{I}}_{0} β^{2} + {\tilde{K}}_{0} {β'}^{2} + 2 {\tilde{E}}_{0} κ_{3} (ν_{3} - 1) + 2 {\tilde{G}}_{0} β (ν_{3} - 1) + 2 {\tilde{H}}_{0} κ_{3} β) .

(24)

Here, $β$ denotes inflation strain which we define as the ratio of radial displacement to the tube’s undeformed radius. The coefficients appearing in the above expression are the stiffnesses in the stress-free state which are defined as follows at a general state of strain $(ν_{3}, κ_{3}, β, β')$ :

\tilde{B} = \frac{\partial^{2} \tilde{Φ}}{\partial κ_{3}^{2}}, \tilde{D} = \frac{\partial^{2} \tilde{Φ}}{\partial ν_{3}^{2}}, \tilde{I} = \frac{\partial^{2} \tilde{Φ}}{\partial β^{2}}, \tilde{K} = \frac{\partial^{2} \tilde{Φ}}{\partial {β'}^{2}}, \tilde{E} = \frac{\partial^{2} \tilde{Φ}}{\partial ν_{3} \partial κ_{3}}, \tilde{G} = \frac{\partial^{2} \tilde{Φ}}{\partial ν_{3} \partial β}, \tilde{H} = \frac{\partial^{2} \tilde{Φ}}{\partial κ_{3} \partial β} .

(25)

Here, $\tilde{I}$ denotes inflation stiffness, $\tilde{G}$ denotes the stiffness corresponding to extension–inflation coupling while $\tilde{H}$ denotes the stiffness corresponding to torsion–inflation coupling. The quantity $\tilde{K}$ denotes the stiffness corresponding to $β'$ while the other stiffnesses have the same meaning as in equation (7).⁵ In addition to the dynamic equations (3) of the special Cosserat rod theory, we have the following additional equation for inflation (see Appendix A for its derivation):

(\frac{\partial \tilde{Φ}}{\partial β'})' - \frac{\partial \tilde{Φ}}{\partial β} = (\overset{\cdot\cdot}{β} - (1 + β) {\overset{\cdot}{κ}}_{3}^{2} s^{2}) I_{3} .

(26)

To obtain the linearized wave solution, we linearize equations (3) and (26) with the following perturbation forms:

r_{3, ϵ} = s ν_{3} + ϵ Δ r_{3} (s, t), θ_{3, ϵ} = s κ_{3} + ϵ Δ θ_{3} (s, t), β_{ϵ} = β + ϵ Δ β (s, t) .

(27)

Assuming the energy form as in (24), the linearized dynamic equations become

\tilde{D} {Δ r ″}_{3} + \tilde{E} {Δ θ ″}_{3} + \tilde{G} Δ β' = ρ_{0} A Δ {\overset{\cdot\cdot}{r}}_{3},

\tilde{E} {Δ r ″}_{3} + \tilde{B} {Δ θ ″}_{3} + \tilde{H} Δ β' = I_{3} Δ {\overset{\cdot\cdot}{θ}}_{3},

\tilde{K} Δ β ″ - \tilde{G} Δ r'_{3} - \tilde{H} Δ θ'_{3} - \tilde{I} Δ β = I_{3} Δ \overset{\cdot\cdot}{β} .

(28)

Upon further substituting the following wave forms for perturbations:

Δ r_{3} = Δ {\hat{r}}_{3} e^{i (ω t - ks)}, Δ θ_{3} = Δ {\hat{θ}}_{3} e^{i (ω t - ks)}, Δ β = Δ \hat{β} e^{i (ω t - ks)},

(29)

we obtain the following generalized eigenvalue problem for the frequencies and mode shapes of longitudinal, torsional, and radial breathing modes:

(\begin{matrix} - k^{2} \tilde{D} & - k^{2} \tilde{E} & - i k \tilde{G} \\ - k^{2} \tilde{E} & - k^{2} \tilde{B} & - i k \tilde{H} \\ i k \tilde{G} & i k \tilde{H} & - k^{2} \tilde{K} - \tilde{I} \end{matrix}) [\begin{matrix} Δ {\hat{r}}_{3} \\ Δ {\hat{θ}}_{3} \\ Δ \hat{β} \end{matrix}] = - ω^{2} (\begin{matrix} ρ_{0} A & 0 & 0 \\ 0 & I_{3} & 0 \\ 0 & 0 & I_{3} \end{matrix}) [\begin{matrix} Δ {\hat{r}}_{3} \\ Δ {\hat{θ}}_{3} \\ Δ \hat{β} \end{matrix}]

(30)

Owing to the presence of coupling stiffnesses, the three modes get coupled to each other. This complicates the procedure to obtain all the stiffnesses from phonon dispersion curves directly. At $k = 0$ however, we note that equation (30) gives only one non-zero eigenvalue which corresponds to the radial breathing mode frequency ( $ω_{r}$ ) and the same can be used to obtain inflation stiffness as follows:

\tilde{I} = I_{3} ω_{r}^{2} |_{k = 0} .

(31)

The characteristic equation of the eigenvalue problem (30) can be solved numerically to obtain dispersion curves of all the three phonon modes. However, in case of stress-free achiral tubes, we obtain the following analytical formulas for the three frequencies:

\begin{array}{l} ω_{(l, r)}^{2} = \frac{k^{2} I_{3} \tilde{D} + ρ_{0} A \tilde{ℐ} + k^{2} ρ_{0} A \tilde{K} \mp \sqrt{{(- k^{2} I_{3} \tilde{D} - ρ_{0} A \tilde{ℐ} - k^{2} ρ_{0} A \tilde{K})}^{2} - 4 I_{3} ρ_{0} A (k^{2} \tilde{D} \tilde{ℐ} + k^{4} \tilde{D} \tilde{K} - k^{2} {\tilde{G}}^{2})}}{2 I_{3} ρ_{0} A}, \\ ω_{t}^{2} = k^{2} \frac{\tilde{ℬ}}{I_{3}} . \end{array}

(32)

Note that the torsional mode is decoupled from the remaining two since the coupling constants $(\tilde{E}, \tilde{H})$ are both zero. However, the longitudinal and radial breathing modes are coupled due to non-zero $\tilde{G}$ . Another point to observe is that the formula for torsional frequency is linear in wavenumber just as in the special Cosserat rod model (see equation (18)). However, the longitudinal frequency is no longer linear in wavenumber. In fact, we shall see later that equation (32) for longitudinal frequency represents the longitudinal phonon dispersion molecular data more accurately. In the special case when $\tilde{G} = 0$ which is typically the case when the Poisson’s ratio is zero, equation (32a) reduces to

ω_{l}^{2} = \frac{k^{2} \tilde{D}}{ρ_{0} A}, ω_{r}^{2} = \frac{\tilde{I} + k^{2} \tilde{K}}{I_{3}} .

(33)

5. Deducing stiffnesses of chiral nanotubes using phonon dispersion data from molecular simulation

We saw in the preceding section that due to the structure of matrix $C$ in equation (15), the longitudinal and torsional modes get decoupled from the bending and shearing modes. However, the bending and shearing modes themselves are coupled in general and their frequencies can be obtained by solving the eigenvalue problem associated with the submatrix $A$ in (16). In fact, we get the following characteristic equation (in terms of five parameters: three stiffnesses ( $\hat{A}, \hat{C}, \hat{F}$ ) and ( $N, M$ )) whose roots yield bending and shearing frequencies for every wavenumber k:

\begin{matrix} (k^{2} N^{2} + ρ_{0} A ω^{2} (k (M + \hat{A} k) - I_{1} ω^{2}) + 2 k^{3} N \hat{F} + k^{4} {\hat{F}}^{2} - k^{2} (k (M + \hat{A} k) - I_{1} ω^{2}) \hat{C} \\ + ν_{3} (- ρ_{0} A ω^{2} (N + 2 k \hat{F}) - k^{2} N \hat{C} + ρ_{0} A ω^{2} \hat{C} ν_{3})) = 0 . \end{matrix}

(34)

We shall use this equation in an inverse way to obtain the five parameters. We first generate phonon dispersion curves from molecular simulation as described in Appendix B. We then take three frequency data corresponding to three different wavenumbers of bending mode and two data from the shearing mode and further substitute them in the above characteristic equation to get five nonlinear equations for the five unknown parameters.⁶ The nonlinear equations are then solved using the Newton–Raphson method.

The characteristic equation corresponding to longitudinal and torsional modes cannot be used to obtain the three stiffnesses ( $\hat{B}, \hat{D}, \hat{E}$ ) since the two dispersion curves associated with them are straight lines. Hence, they will at most generate only two independent equations for the three unknown stiffnesses. However, we can generate sufficient equations by matching the mode shape of the phonon mode with the eigenvector from the continuum approach. This procedure is discussed now in detail.

5.1. Extracting ( $Δ {\hat{r}}_{3}, Δ {\hat{θ}}_{3}, Δ \hat{β}$ ) from atomic displacements

In the case of longitudinal, torsional, and radial breathing modes, only ( $Δ {\hat{r}}_{3}, Δ {\hat{θ}}_{3}, Δ \hat{β}$ ) are non-zero. We first present a way to obtain ( $Δ {\hat{r}}_{3}, Δ {\hat{θ}}_{3}, Δ \hat{β}$ ) from displacement data of atoms in these phonon modes. As described in Appendix B, the displacement of an atom ( $q, l$ ) in the phonon mode ( $p, v$ ) is written as (see equation (67) for the formula of ${\hat{u}}_{(q, l)}^{[p, v]}$ ):

u_{(q, l)}^{[p, v]} = {\hat{u}}_{(q, l)}^{[p, v]} \exp (- i ω t) .

(35)

The quantity ${\hat{u}}_{(q, l)}^{[p, v]}$ is usually complex and contains information of both amplitude and phase of displacement. We drop the superscript $[p, v]$ from now on since we always look at a particular phonon mode while extracting ( $Δ {\hat{r}}_{3}, Δ {\hat{θ}}_{3}, Δ \hat{β}$ ). Let us denote the three components of ${\hat{u}}_{(q, l)}$ by

{\hat{u}}_{(q, l)} = [\begin{matrix} {({\hat{u}}_{1})}_{(q, l)} \\ {({\hat{u}}_{2})}_{(q, l)} \\ {({\hat{u}}_{3})}_{(q, l)} \end{matrix}] .

(36)

Here, ${\hat{u}}_{1}$ and ${\hat{u}}_{2}$ denote displacement amplitudes in transverse directions $e_{1}$ and $e_{2}$ , respectively, whereas ${\hat{u}}_{3}$ denotes the same along the axis of the nanorod/nanotube. In order to obtain $Δ {\hat{r}}_{3}$ , we simply take the projection of displacement of each atom along the tube’s axis and sum over all the atoms within the nanorod’s “cross-section,”⁷ i.e.,

Δ {\hat{r}}_{3} = \frac{1}{m} \sum_{l = 0}^{m - 1} {(\hat{u} \cdot e_{3})}_{(0, l)} = \frac{1}{m} \sum_{l = 0}^{m - 1} {({\hat{u}}_{3})}_{(0, l)} .

(37)

Here, “ $m$ ” denotes the total number of atoms within the nanorod’s cross-section.

To extract $Δ {\hat{θ}}_{3}$ , which denotes rotation of a cross-section about the nanorod’s axis, we take the projection of atomic displacement in the direction perpendicular to the radial line and then sum over all “ $m$ ” atoms as follows:

Δ {\hat{θ}}_{3} = \frac{1}{m} \sum_{l = 0}^{m - 1} [{(\hat{u} \times X)}_{(0, l)} \cdot e_{3}] / R_{(0, l)}^{2} = \frac{1}{m} \sum_{l = 0}^{m - 1} [{({\hat{u}}_{2} X_{1})}_{(0, l)} - {({\hat{u}}_{1} X_{2})}_{(0, l)}] / R_{(0, l)}^{2} .

(38)

Here, $X_{(0, l)}$ denotes the relaxed position of the atom $(0, l)$ in the stretched and twisted nanorod and $R_{(0, l)} = \sqrt{X_{1}^{2} + X_{2}^{2}}$ denotes its radial coordinate.

Finally, to obtain $Δ \hat{β}$ , we first obtain the projection of atomic displacement in the radial direction and then sum as follows:

Δ \hat{β} = \frac{1}{m} \sum_{l = 0}^{m - 1} [{({\hat{u}}_{1} X_{1})}_{(0, l)} + {({\hat{u}}_{2} X_{2})}_{(0, l)}] / R_{(0, l)}^{2} .

(39)

5.2. Deducing stiffnesses ( $\hat{B}, \hat{D}, \hat{E}$ ) of chiral nanotubes

As only the submatrix $B$ influences longitudinal and torsional modes, the eigenvalue–eigenvector equation associated with this submatrix in equation (15) can be rearranged as a linear problem for the unknown stiffnesses as follows:

(\begin{matrix} k^{2} Δ {\hat{r}}_{3} & k^{2} Δ {\hat{θ}}_{3} & 0 \\ 0 & k^{2} Δ {\hat{r}}_{3} & k^{2} Δ {\hat{θ}}_{3} \end{matrix}) [\begin{matrix} \hat{D} \\ \hat{E} \\ \hat{B} \end{matrix}] = ω^{2} [\begin{matrix} ρ_{0} A Δ {\hat{r}}_{3} \\ I_{3} Δ {\hat{θ}}_{3} \end{matrix}] .

(40)

Finally, we take one frequency data each from longitudinal and torsional phonon dispersion data, extract the corresponding mode shape and substitute them into the above equation to get a total of four equations for the three unknowns. We then solve this system using a least-squares method to obtain all the stiffnesses.

5.3. Deducing stiffnesses ( $\tilde{B}, \tilde{D}, \tilde{E}, \tilde{G}, \tilde{H}, \tilde{I}, \tilde{K}$ )

The inflation stiffness $\tilde{I}$ can be obtained independently using equation (31). We then rearrange the generalized eigenvalue problem (30) in terms of the remaining unknown stiffnesses as follows:

(\begin{matrix} k^{2} Δ {\hat{r}}_{3} & k^{2} Δ {\hat{θ}}_{3} & 0 & i k Δ \hat{β} & 0 & 0 \\ 0 & k^{2} Δ {\hat{r}}_{3} & k^{2} Δ {\hat{θ}}_{3} & 0 & i k Δ \hat{β} & 0 \\ 0 & 0 & 0 & - i k Δ {\hat{r}}_{3} & - i k Δ {\hat{θ}}_{3} & k^{2} Δ \hat{β} \end{matrix}) [\begin{matrix} \tilde{D} \\ \tilde{E} \\ \tilde{B} \\ \tilde{G} \\ \tilde{H} \\ \tilde{K} \end{matrix}] = [\begin{matrix} ω^{2} ρ_{0} A Δ {\hat{r}}_{3} \\ ω^{2} I_{3} Δ {\hat{θ}}_{3} \\ (ω^{2} I_{3} - \tilde{I}) Δ \hat{β} \end{matrix}] .

(41)

From the phonon solution, we then take one frequency data each from longitudinal, torsional, and radial breathing modes and also extract the corresponding mode shapes. These frequencies and mode shapes are then substituted in the above equation to get a total of nine equations for the six unknowns. Finally, we solve this system using a least-squares method to obtain all the stiffnesses. Note that one cannot obtain the six unknown stiffnesses by forming a system of six equations using just the longitudinal and torsional mode data. This is due to the fact that these two modes are insensitive to $\tilde{K}$ in the long wavelength limit.

6. Deducing stiffnesses of SWCNTs from phonon dispersion curve

We now deduce stiffnesses of both achiral and chiral SWCNTs using the methodology presented in the preceding section. The phonon dispersion curves are obtained using the SWCNT’s two-atom helically repeating cell. The readers are requested to go through Appendix B in order to get familiar with the notation used in this section. Owing to the presence of a group generator g (see Appendix A) in longitudinal, torsional, and radial breathing phonon modes, the summation in Section 5.1 could be done over just the two-atom repeating cell of the SWCNT instead of all atoms in the cross-section. The Tersoff–Brenner potential [31] is used to model interactions between carbon atoms. For the stiffnesses of the special Cosserat rod model, we compare our results with those obtained using the HCB rule of [7]. Similarly, for the general Cosserat rod stiffnesses, we compare our results with those obtained in [32].

6.1. (10,10) SWCNT

To generate a (10,10) stress-free SWCNT using two-atom basis, the following group generators are used: $h = (R_{θ_{0}} | τ_{0} e_{3})$ and $g = (R_{ψ} | 0)$ . Here, $τ_{0} = 0.1257$ , $θ_{0} = π / 10$ , and $ψ = π / 5$ .

6.1.1. Stiffnesses in the stress-free configuration

Let $(k_{1}, k_{2})$ denote the components of the wavevector $k$ and P denote perimeter of the nanotube’s cross-section. The three lowest curves in Figure 2(a) corresponding to $k_{2} P / π$ =0 are the longitudinal, torsional, and radial breathing modes, respectively, while the two lowest curves in Figure 2(b) corresponding to $k_{2} P / π = 0.2$ are the bending and shearing modes.⁸ The mode shapes of bending and shearing phonons are also presented in Figure 3(a) and (b), respectively.⁹

Figure 2.

Phonon dispersion curves of a (10,10) stress-free SWCNT depicting (a) longitudinal, torsional, and radial breathing modes and (b) shearing and bending modes.

Figure 3.

(a) Bending mode observed in the lowest branch of Figure 2(b) at $k_{1} τ_{0} = π / 10$ . (b) Shearing mode observed in the second lowest branch of Figure 2(b) at $k_{1} τ_{0} = π / 10$ .

In the case of longitudinal, torsional, and radial breathing modes, the continuum wavenumber k is the same as $k_{1}$ in phonon calculations. In the case of bending and shearing modes however, $k = k_{1} - θ_{0} / π$ . Finally, we obtain all the stiffnesses in the stress-free state using the methodology presented in Section 5 and present the same in Table 1. We can note that the stiffnesses from the phonon data match very well with the stiffnesses calculated from other techniques.

Table 1.

(I) Stiffnesses of a stress-free (10,10) SWCNT obtained from phonon data and HCB rule [7]. (II) Stiffnesses from phonon data and the scheme of [32]. “×” indicates the corresponding quantity cannot be obtained.

		${\hat{A}}_{0}$	${\hat{B}}_{0}$	${\hat{C}}_{0}$	${\hat{D}}_{0}$	${\hat{E}}_{0}$	${\hat{F}}_{0}$	$N$	$M$
I	Phonon data	1547.068	1067.827	1120.654	6362.219	0	0	0	0
	HCB rule	1544.278	1067.817	1085.856	6362.567	0	0	0	0
		${\tilde{B}}_{0}$	${\tilde{D}}_{0}$	${\tilde{I}}_{0}$	${\tilde{E}}_{0}$	${\tilde{G}}_{0}$	${\tilde{H}}_{0}$	${\tilde{K}}_{0}$
II	Phonon data	1067.827	7681.928	7677.448	0	3192.261	0	−21.377
	Singh et al. [32]	1067.817	7681.493	7690.887	0	3184.919	0	×

There is a slight difference in only shearing stiffness of about 3% (marked in bold in Table 1). To investigate this discrepancy, we obtained shearing stiffness using our approach for arm-chair SWCNTs of different radii, normalized these values with those obtained from HCB rule and present in Figure 4(a). It can be observed that with increasing nanotube radius, the difference between the two approaches becomes insignificant. This difference may be existing because the phonon approach is based on linearized dynamics while the HCB rule is based on static energy minimization.¹⁰

Figure 4.

Normalized shearing stiffness versus radius in arm-chair SWCNTs of different radii.

6.1.2. Stiffnesses in uniformly twisted configuration

We now evaluate various stiffnesses for the twisted (10,10) nanotube at $κ_{3} = 0.2$ . Twisting of a nanotube makes it chiral. Accordingly, the longitudinal, torsional, and radial breathing modes all get coupled. Similarly, the bending and shearing modes also get coupled and their mode shapes look completely different too (see Figure 5). The long wavelength limits of bending and shearing modes now occur at $k_{1} = \frac{θ_{0} + κ_{3} τ_{0}}{π}$ . We again obtain the stiffnesses as well as the induced axial force and twisting moment using our approach and present them in Table 2. Owing to the inherent difference in shearing stiffness, we also observe a large difference in the bending–shearing coupling stiffness $\hat{F}$ (marked in bold in Table 2). However, all other values match pretty well. We also note that all coupling stiffnesses are now non-zero since the tube has now become chiral.

Figure 5.

(a) Out-of-plane bending mode and (b) shearing mode in a twisted (10,10) SWCNT.

Table 2.

(I) Stiffnesses, axial force, and twisting moment of a twisted (10,10) SWCNT obtained from phonon data and HCB rule at $κ_{3} = 0.2$ . (II) Stiffnesses from phonon data and the scheme of [32]. “×” indicates the corresponding quantity cannot be obtained.

		$\hat{A}$	$\hat{B}$	$\hat{C}$	$\hat{D}$	$\hat{E}$	$\hat{F}$	$N$	$M$
I	Phonon data	1495.088	785.280	868.074	6303.956	−258.776	65.537	−28.556	194.463
	HCB rule	1499.3777	785.273	843.849	6304.307	−258.752	26.096	−28.556	194.439
		$\tilde{B}$	$\tilde{D}$	$\tilde{I}$	$\tilde{E}$	$\tilde{G}$	$\tilde{H}$	$\tilde{K}$
II	Phonon data	834.123	7649.071	7479.076	−2.363	3181.490	605.958	−35.951
	Singh et al. [32]	834.210	7648.650	7493.529	−2.259	3173.937	605.569	×

6.2. Stress-free chiral (10,7) SWCNT

We finally take up the example of a (10,7) chiral nanotube and evaluate its stiffnesses in the stress-free configuration. The group parameters in this case are $τ_{0} = 0.0147$ , $θ_{0} = 1.8506$ , and $ψ = 2 π$ . We then evaluate the stiffnesses using the methodology presented in Section 5 and further present in Table 3. Again, owing to inherent chirality, all the coupling stiffnesses are non-zero. The application of HCB rule to stress-free chiral nanotubes is not straightforward. We accordingly do not have data corresponding to HCB rule in Table 3. However, we have compared our data with the data obtained using the theory of objective structures (see [6]). However, the latter approach cannot deduce bending stiffness, shearing stiffness, and bending–shearing coupling stiffnesses of chiral tubes.

Table 3.

(I) Stiffnesses of a (10,7) SWCNT obtained from phonon data. (II) Stiffnesses from phonon data and the scheme of [32]. “×” indicates the corresponding quantity cannot be obtained.

		${\hat{A}}_{0}$	${\hat{B}}_{0}$	${\hat{C}}_{0}$	${\hat{D}}_{0}$	${\hat{E}}_{0}$	${\hat{F}}_{0}$	$N$	$M$
I	Phonon data	962.849	663.522	954.404	5418.947	9.642	−2.504	0	0
	James [6]	×	663.929	×	5417.246	9.690	×	0	0
		${\tilde{B}}_{0}$	${\tilde{D}}_{0}$	${\tilde{I}}_{0}$	${\tilde{E}}_{0}$	${\tilde{G}}_{0}$	${\tilde{H}}_{0}$	${\tilde{K}}_{0}$
II	Phonon data	663.525	6550.894	6550.122	8.509	2734.638	−2.796	−17.568
	Singh et al. [32]	663.524	6548.928	6565.888	8.524	2723.65	−2.696	×

7. Unusual effects of extension and torsion on phonon dispersion curves of (10,10) SWCNTs

We now present some interesting effects of extension and torsion on phonon dispersion curves of (10,10) SWCNTs and further show that they are accurately replicated using our continuum formulas.

7.1. Coupling of longitudinal, torsional, and radial breathing modes

We note from equation (30) that the longitudinal, torsional, and radial breathing modes get coupled when the coupling stiffnesses $\tilde{E}, \tilde{H}, \tilde{G}$ are all non-zero which is the case for stress-free chiral nanotubes or twisted chiral/achiral nanotubes. However, only $\tilde{G}$ is non-zero in stress-free achiral nanotubes. Accordingly, their torsional mode gets decoupled from longitudinal and radial breathing modes (see Figure 6(a)). The latter two are coupled to each other though and it becomes difficult to distinguish between the two modes. One can then look at the ratio of the amplitudes of radial and longitudinal displacements ( $| \frac{Δ \hat{β} R}{Δ {\hat{r}}_{3}} |$ ) from phonon data of these modes to distinguish. A ratio less than unity would imply the mode to be longitudinal. We can also obtain the formula for this ratio from the eigenvectors of our continuum formulation using which we find that this ratio becomes unity at

k = \sqrt{\frac{\tilde{I}}{- \tilde{K} + R^{2} \tilde{D}}},

(42)

which characterizes the wavenumber at which the two modes transition into one another. This transition phenomena has also been reported earlier [35]. However, the analytical formula at which this transition occurs has been missing. We note from Figure 6(a) that the dispersion curves as well as the transition point are very well replicated using our general Cosserat rod continuum model. In Figure 6(b), we again look at these dispersion curves but for the twisted state of the SWCNT. Owing to twist, all three modes get coupled to each other now. Again, our general Cosserat rod continuum model accurately replicates the phonon dispersion molecular data. We are not reporting the continuum formulas for the two transition points (longitudinal↔radial and radia↔torsional) though since they turn out to be very long. However, their numerical values match with the ratio obtained from molecular data. We find that the lowest branch transitions from torsional to radial breathing mode; the second lowest branch transitions from longitudinal to radial breathing mode and then finally to torsional mode while the third branch transitions from radial breathing to longitudinal mode. Another point to note is that the special Cosserat rod model (as per equation (17)) fails to capture the trend from molecular data at high wavenumber.

Figure 6.

Longitudinal, torsional, and radial breathing modes of a (10,10) SWCNT: (a) stress-free configuration; (b) twisted configuration ( $κ_{3} = 0.2$ ).

To get an idea of the three frequencies and the transition point in terms of material constants, we think of stress free achiral SWCNTs as a linearly isotropic thin tube for which the formulas of various stiffnesses are already available [32]. We report the same here:

{\tilde{B}}_{0} = 2 π R^{3} TG, {\tilde{D}}_{0} = {\tilde{I}}_{0} = \frac{2 π RTE}{(1 - ν^{2})}, {\tilde{G}}_{0} = ν {\tilde{I}}_{0}, {\tilde{E}}_{0} = 0, {\tilde{H}}_{0} = 0 .

(43)

Here, $ν$ is the Poisson’s ratio, E is the Young’s modulus, G is the shear modulus, and T denotes the tube’s thickness. The formula for ${\tilde{K}}_{0}$ is not available in literature but its numerical value is very low compared with other stiffnesses (see Table 1). Furthermore, it is found to have insignificant effect on dispersion curves in Figure 6. Accordingly, we set it to zero. Upon substituting these stiffnesses in equation (32), we finally obtain the following formulas for the frequencies of the three modes in stress-free state:

ω_{l, r}^{2} = \frac{E (1 + k^{2} R^{2} \mp \sqrt{k^{2} R^{2} (k^{2} R^{2} + 4 ν^{2} - 2) + 1})}{2 ρ_{0} (1 - ν^{2}) R^{2}}, ω_{t} = k \sqrt{\frac{G}{ρ_{0}}} .

(44)

The radial breathing mode frequency further simplifies to the following in long wavelength limit:

ω_{r} |_{k = 0} = \frac{1}{R} \sqrt{\frac{E}{ρ_{0} (1 - ν^{2})}} .

(45)

The inverse dependence of radial breathing mode frequency on tube’s radius is well known [20]. We also find the transition wavenumber in equation (42) to be inverse of the tube’s radius, i.e.,

k = \frac{1}{R} .

(46)

7.2. Effect on bending phonon frequency

We now discuss the effect of extension, compression, and twist on bending phonon dispersion curve. We first consider the effect of extension. In Figure 7, we present this dispersion curve obtained using atomistic calculation as well as by solving equation (34). Both continuum and atomistic data match pretty well. We also observe that the dispersion curve has non-zero slope at $k = 0$ when the nanotube is axially strained. Furthermore, the magnitude of the slope increases with strain since the induced axial force also increases. This is in complete agreement with formula (20) derived earlier.

Figure 7.

Bending frequency in the stress-free and axially stretched state.

We then consider the bending dispersion curve in compressed and twisted state. We again note from Figure 8 that continuum and atomistic results match pretty well. We are presenting the square of bending frequency here which exhibits unusual pattern. In compressed state, the dispersion curve is symmetric about $k = 0$ . Furthermore, it first decreases to become negative and then increases monotonically. The initial decreasing trend is due to the coefficient of $k^{2}$ term being negative in formula (19) since the axial force is compressive now. Another interesting point to note is the asymmetry in dispersion curve in the nanotube’s twisted state. This asymmetry is due to non-zero coefficient of $k^{3}$ in formula (19) since the twisting moment $M$ is non-zero now.

Figure 8.

Bending frequency in the stress-free, axially compressed ( $ν_{3} = 0.99$ ) and twisted state ( $κ_{3} = 0.2$ )

Finally, using the continuum formula for bending frequency, we also obtained the wavenumber at which $ω^{2}$ crosses zero to become positive. It turns out to be

\begin{matrix} k = \frac{- 2 N \hat{F} + M \hat{C} \pm \sqrt{\hat{C} (4 N^{2} \hat{A} - 4 M N \hat{F} + M^{2} \hat{C} + 4 N ν_{3} ({\hat{F}}^{2} - \hat{A} \hat{C}))}}{2 ({\hat{F}}^{2} - \hat{A} \hat{C})} \\ = \pm \sqrt{\frac{N (N - \hat{C} ν_{3})}{\hat{A} \hat{C}}} (at zero twist) . \end{matrix}

(47)

All bending phonons having wavenumber less than the above value are unstable. In fact, this instability is related to Euler buckling of clamped–clamped rods. As the nanotube is compressed and twisted, all phonons having wavelength larger than the critical buckling length corresponding to the induced compressive force and twisting moment will automatically become unstable owing to Euler buckling. Conversely, the formula (47) can also be used to find the buckling load of compressed and twisted continuum rods or nanotubes of length $L$ : we just need to set $k = 2 π / L$ on the left-hand side. Note that this formula is more general than the classical buckling load formula of $4 π^{2} \hat{A} / L^{2}$ because the classical formula holds only for unshearable and inextensible Kirchhoff rods at zero twist whereas the formula (47) holds for special Cosserat rods which can shear as well as stretch.

7.3. Instability in torsional mode at large twist

We now show in Figure 9(a) the effect of twisting a nanotube on its torsional dispersion curve. The continuum solution here corresponds to torsional wave solution of the general Cosserat rod model. We note that the continuum and atomistic solutions match pretty well. They start to deviate at higher wavenumber though. An interesting point to note is the negative values for $ω^{2}$ at large twist which signals instability. Indeed, the expression $\hat{B} \hat{D} - {\hat{E}}^{2}$ becomes negative in this case which renders the elasticity tensor of special Cosserat rod theory negative definite [29]. This instability due to large twist was also reported in [36], but they found it by tracking eigenvalues of the Hessian matrix of the full length SWCNT. For readers’ reference, we have also presented the phonon mode shape of this case in Figure 9(b).

Figure 9.

(a) Torsional dispersion curve in twisted state. (b) Phonon mode shape for $κ_{3} = 0.5$ .

8. Conclusions

We have presented a novel technique to obtain stiffnesses of chiral nanorods and nanotubes modeled as Cosserat rods from their phonon dispersion curves. The approach is straightforward having no end effects as prevalent in several existing techniques. We have also derived continuum wave solutions in chiral Cosserat rods for longitudinal, torsional, bending, shearing, and radial breathing modes and obtained formulas for their frequencies in terms of the rod’s stiffnesses as well as the induced axial force and twisting moment. Along the way, we have derived the dynamic equation corresponding to inflation of tubes in Appendix A. Our study allows a clearer understanding of how the rod’s coupling stiffnesses (extension–torsion–inflation and bending–shearing couplings) as well as twisting of the rod lead to coupling of various wave modes in continuum tubes and that of phonon dispersion modes in nanorods and nanotubes. We also saw unusual effects of stretching and twisting a nanotube on its phonon dispersion curves which we were able to reproduce accurately using our continuum formulas. Finally, the continuum formulas derived could be useful in tuning the dispersion curves by changing the nanostructures’s stiffness properties as well as by stretching and twisting a nanostructure. Such tunings could be potentially useful in blocking or facilitating waves of certain frequency range to pass through these structures.

Footnotes

Appendix A. Derivation of the dynamic equation corresponding to inflation

Appendix B. Generating phonon dispersion curves for SWCNTs

In this appendix, we present how to generate phonon dispersion curves of stretched and twisted SWCNTs by using helical symmetry inherent in them. The whole SWCNT can be generated using its two-atom basis and applying helical symmetry operations to it. If we denote by $x_{i, l}$ the position of the ith image of atom l in the two-atom basis unit cell, then

(61)

x_{i, l} \equiv x_{(i_{1}, i_{2}), l} = h^{i_{1}} g^{i_{2}} (x_{(0, 0), l}) : i_{1} \in Z, i_{2} = 0, 1, 2, \dots, c - 1, l = 0, 1 .

Here $c = \gcd (n, m)$ for a $(n, m)$ SWCNT. The mathematical form of the generating elements h and g are as follows:

h = {R_{θ} | τ e_{3}},

(62)

g = {R_{ψ} | 0}, ψ = 2 π / c .

Here $R_{θ}$ and $R_{ψ}$ denote rotations by $θ$ and $ψ$ , respectively, about the nanotube’s axis. The generator h rotates and translates an arbitrary point $x$ as follows: $h (x) = R_{θ} x + τ e_{3}$ . The generator g on the other hand does pure rotation. Basically, the generator g when applied to the two-atom basis unit cell generates the nanotube’s “cross-section” (ring of atoms). We then apply h to this cross-section to generate the whole nanotube.¹¹ The parameters ( $θ, τ$ ) in h change as follows when the tube is stretched and twisted:

(63)

τ = ν_{3} τ_{0}, θ = θ_{0} + κ_{3} τ_{0} .

Here ( $τ_{0}, θ_{0}$ ) are the reference state group parameters (see Section 7 of [32]). To generate the whole SWCNT, one just needs to know the position of the two atoms ${x_{(0, 0), l}}_{l = 0}^{1}$ which is obtained by minimizing the unit cell energy (see [32]). We then follow the methodology of [12] to obtain phonon dispersion curves which we summarize below.

For illustration, we present the phonon dispersion curves obtained using the two-atom basis as well as using the periodic unit cell for a (10,10) SWCNT in Figure 10. The periodic unit cell methodology has 40 atoms per unit cell and accordingly it has 120 branches for every wave vector. It is just that in case of periodic unit cell, the wave vector is a scalar whereas in case of two-atom basis, the wave vector has two components ( $k_{1}, k_{2}$ ): $k_{1}$ is the component in the continuous direction and $k_{2}$ is the component in the discrete direction.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: P Gupta acknowledges the financial support received from a DST-INSPIRE fellowship and A Kumar acknowledges the support from SERB, India (grant number YSS/2014/000023).

ORCID iD

Ajeet Kumar

Notes

References

Dove

. Introduction to Lattice Dynamics, Vol. 4. Cambridge: Cambridge University Press, 1993.

Born

Huang

. Dynamical Theory of Crystal Lattices. Oxford: Oxford University Press, 1998.

Wallace

. Thermodynamics of Crystals. Courier Corporation, 1998.

Dresselhaus

Eklund

. Phonons in carbon nanotubes. Adv Phys 2000; 49(6): 705–814.

Elliott

Triantafyllidis

Shaw

. Stability of crystalline solids—I: Continuum and atomic lattice considerations. J Mech Phys Solids 2006; 54(1): 161–192.

James

. Objective structures. J Mech Phys Solids 2006; 54(11): 2354–2390.

Kumar

Gupta

. A helical Cauchy–Born rule for special Cosserat rod modeling of nano and continuum rods. J Elast 2016; 124: 81–106.

Popov

Van Doren

Balkanski

. Lattice dynamics of single-walled carbon nanotubes. Phys Rev B 1999; 59(13): 8355.

Milošević

Dobardžić

Damnjanović

. Phonons in narrow carbon nanotubes. Phys Rev B 2005; 72(8): 085426.

10.

Gunlycke

Lawler

White

. Lattice vibrations in single-wall carbon nanotubes. Phys Rev B 2008; 77(1): 014303.

11.

Jiang

Tang

Wang

. A lattice dynamical treatment for the total potential energy of single-walled carbon nanotubes and its applications: relaxed equilibrium structure, elastic properties, and vibrational modes of ultra-narrow tubes. J Phys Condens Matter 2008; 20(4): 045228.

12.

Aghaei

Dayal

Elliott

. Symmetry-adapted phonon analysis of nanotubes. J Mech Phys Solids 2013; 61(2): 557–578.

13.

Jishi

Venkataraman

Dresselhaus

. Phonon modes in carbon nanotubules. Chem Phys Lett 1993; 209(1–2): 77–82.

14.

Kürti

Kresse

Kuzmany

. First-principles calculations of the radial breathing mode of single-wall carbon nanotubes. Phys Rev B 1998; 58(14): R8869.

15.

Kuzmany

Burger

Hulman

Kürti

Rinzler

Smalley

. Spectroscopic analysis of different types of single-wall carbon nanotubes. Europhys Lett 1998; 44(4): 518.

16.

Graff

. Wave Motion in Elastic Solids. New York: Dover Publications, 1975.

17.

Rego

. Quantized thermal conductance of dielectric quantum wires. Phys Rev Lett 1998; 81(1): 232.

18.

Popov

Van Doren

Balkanski

. Elastic properties of single-walled carbon nanotubes. Phys Rev B 2000; 61(4): 3078.

19.

Saito

Takeya

Kimura

Dresselhaus

. Raman intensity of single-wall carbon nanotubes. Phys Rev B 1998; 57(7): 4145.

20.

Sánchez-Portal

Artacho

Soler

Rubio

Ordejón

. Ab initio structural, elastic, and vibrational properties of carbon nanotubes. Phys Rev B 1999; 59(19): 12678.

21.

Mahan

Jeon

. Flexure modes in carbon nanotubes. Phys Rev B 2004; 70(7): 075405.

22.

Zimmermann

Pavone

Cuniberti

. Vibrational modes and low-temperature thermal properties of graphene and carbon nanotubes: Minimal force-constant model. Phys Rev B 2008; 78(4): 045410.

23.

Iijima

Brabec

Maiti

Bernholc

. Structural flexibility of carbon nanotubes. J Chem Phys 1996; 104(5): 2089–2092.

24.

Cao

Chen

. The effects of chirality and boundary conditions on the mechanical properties of single-walled carbon nanotubes. Int J Solids Struct 2007; 44(17): 5447–5465.

25.

Chandraseker

Mukherjee

Paci

Schatz

. An atomistic-continuum Cosserat rod model of carbon nanotubes. J Mech Phys Solids 2009; 57: 932–958.

26.

Aghaei

Dayal

Elliott

. Anomalous phonon behavior of carbon nanotubes: First-order influence of external load. J Appl Phys 2013; 113(2): 023503.

27.

Chang

. Radial breathing mode frequency of single-walled carbon nanotubes under strain. Appl Phys Lett 2008; 93(6): 061901.

28.

Antman

. Nonlinear Problems of Elasticity. New York: Springer-Verlag, 1995.

29.

Healey

. Material symmetry and chirality in nonlinearly elastic rods. Math Mech Solids 2002; 7: 405–420.

30.

Kumar

Mukherjee

. A geometrically exact rod model including in-plane cross-sectional deformation. J Appl Mech 2011; 78: 011010.

31.

Brenner

. Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films. Phys Rev B 1990; 42: 9458.

32.

Singh

Kumar

. Effect of intrinsic twist and orthotropy on extension–twist–inflation coupling in compressible circular tubes. J Elast 2017; 128(2): 175–201.

33.

Cowper

. The shear coefficient in Timoshenko’s beam theory. J Appl Mech 1966; 33(2): 335–340.

34.

Hutchinson

. Shear coefficients for Timoshenko beam theory. J Appl Mech 2001; 68(1): 87–92.

35.

Ghavanloo

Fazelzadeh

. Vibration characteristics of single-walled carbon nanotubes based on an anisotropic elastic shell model including chirality effect. Appl Math Model 2012; 36(10): 4988–5000.

36.

Gupta

Kumar

. Effect of material nonlinearity on spatial buckling of nanorods and nanotubes. J Elast 2017; 126(2): 155–171.

Phonons in chiral nanorods and nanotubes: A Cosserat-rod-based continuum approach

Abstract

Keywords

1. Introduction

Notation

2. Brief overview of the theory of special Cosserat rods

3. Linearized wave solution in stretched and twisted rods

3.1. Longitudinal and torsional mode frequencies

3.2. Bending and shearing mode frequencies

4. Continuum wave solution for radial breathing mode

5. Deducing stiffnesses of chiral nanotubes using phonon dispersion data from molecular simulation

5.1. Extracting ( Δ r ^ 3 , Δ θ ^ 3 , Δ β ^ ) from atomic displacements

5.2. Deducing stiffnesses ( B ^ , D ^ , E ^ ) of chiral nanotubes

5.3. Deducing stiffnesses ( B ~ , D ~ , E ~ , G ~ , H ~ , I ~ , K ~ )

6. Deducing stiffnesses of SWCNTs from phonon dispersion curve

6.1. (10,10) SWCNT

6.1.1. Stiffnesses in the stress-free configuration

6.1.2. Stiffnesses in uniformly twisted configuration

6.2. Stress-free chiral (10,7) SWCNT

7. Unusual effects of extension and torsion on phonon dispersion curves of (10,10) SWCNTs

7.1. Coupling of longitudinal, torsional, and radial breathing modes

7.2. Effect on bending phonon frequency

7.3. Instability in torsional mode at large twist

8. Conclusions

Footnotes

Appendix A. Derivation of the dynamic equation corresponding to inflation

Appendix B. Generating phonon dispersion curves for SWCNTs

Funding

ORCID iD

Notes

References

5.1. Extracting ( $Δ {\hat{r}}_{3}, Δ {\hat{θ}}_{3}, Δ \hat{β}$ ) from atomic displacements

5.2. Deducing stiffnesses ( $\hat{B}, \hat{D}, \hat{E}$ ) of chiral nanotubes

5.3. Deducing stiffnesses ( $\tilde{B}, \tilde{D}, \tilde{E}, \tilde{G}, \tilde{H}, \tilde{I}, \tilde{K}$ )