Carbon Nanostructures as an Electromechanical Bicontinuum
Carbon Nanostructures as an Electromechanical Bicontinuum
Cristiano Nisoli∗, Paul E. Lammert∗, Eric Mockensturm† and Vincent H. Crespi∗
Department of Physics and Materials Research Institute
Department of Mechanical and Nuclear Engineering
The Pennsylvania State University, University Park, PA 16802-6300
(Dated: November 2, 2018)
A two-field model provides an unifying framework for elasticity, lattice dynamics and electrome-
chanical coupling in graphene and carbon nanotubes, describes optical phonons, nontrivial acoustic
branches, strain-induced gap opening, gap-induced phonon softening, doping-induced deformations,
and even the hexagonal graphenic Brillouin zone, and thus explains and extends a previously dis-
parate accumulation of analytical and computational results.
PACS numbers: 62.25.+g, 81.05.Tp, 63.22.+m, 77.65.-j, 46.05.+b
Vibrations in carbon nanostructures such as tubes,
fullerenes, or graphene sheets [1, 2, 3] have a ubiquitous
influence on electronic, optical and thermal response:
scattering from optical phonons limits charge transport
in otherwise ballistic nanotube conductors [4, 5]; twist
deformations gap metallic tubes [6, 7]; ballistic phonons
transport heat in nanotubes with great efficiency [8, 9,
10]; resonant Raman spectroscopy can unambiguously
identify a tube’s wrapping indices (n,m) [11, 12, 13, 14];
electron-phonon interactions may ultimately limit the
electrical performance of graphene [15, 16]. Computa-
tionally intensive atomistic models of lattice dynamics
often lack simplified model descriptions that can facili-
tate insight, yet traditional analytical continuum mod-
els [1, 2, 17, 18], while very useful and important, cannot
describe atomistic phenomena without phenomenologi-
cal extensions [19, 20, 21]. Although continuum models
are restricted to long-wavelength physics, they have been
used to describe atomic-scale phenomena in bulk binary
compounds by incorporating a separate continuum field
for each sublattice [23]: in graphene, two fields are neces-
sary. Here we present an analytical “bicontinuum” model
that represents the full atomistic detail of the graphenic
lattice, including optical modes, nonlinear dispersion of
in-plane phonons, electromechanical effects and even the
hexagonal graphenic Brillouin zone, a construct generally
held to be exclusively atomistic.
Graphene decomposes into the two triangular sublat-
tices of Fig. 1. We describe in-plane deformations of the
sublattices via two fields, ui(x), vi(x), i = 1, 2, and their
strain tensors uij = ∂(iuj) and vij = ∂(jvi). The density
of elastic energy contains direct and cross terms:
V [u, v] = d[u] + d[v] + c[u, v]. (1)
Six-fold symmetry of the sublattices implies isotropy of
the direct terms [24]:
d[u] = µ′ uijuij +
j. (2)
Symmetry dictates the form of the cross term
FIG. 1: The two sublattices (circles and squares) of graphene
and the three unit vectors ê(l) used in the text. φ, z are
cylindrical coordinates of a tube, while Ψ = π/6− θc with θc
the chiral angle. Also, anisotropic (uxx = uxy = 0, uyy = 2γ,
qx = ℓ γ), shear (uxx = uyy = 0, uxy = η, qy = −ℓ η) strains.
c[u, v] = 2 µuijvij + λu
+ α (u− v)2 (3)
− β eijk
uij + vij
uk − vk
The tensor eijk, which is invariant under C3v, can be
represented by the three unit vectors {ê(l)} of Fig. 1:
eijk =
k . (4)
Only the last term in Eq. 4 is not invariant under gen-
eral rotation. (In nanotubes, it depends on the helical
angle θc: eφφφ = −eφzz = − sin(3θc), ezzz = −eφφz =
− cos(3θc), where φ, z are defined in Fig. 1). This elastic
energy density, the lowest-order approximation in both
derivatives and fields, contains six parameters: µ′ and λ′,
being confined to one sublattice, describe next-neighbor
interactions; the cross terms µ and λ describe nearest-
neighbor interaction; α describes the stiffness against rel-
ative shifts of the sublattices; β determines the strength
of rotational symmetry breaking and so carries the point
group symmetry of graphene. These parameters are nor-
malized to the sublattice surface density σs, so that the
elastic energy is W =
σs V d
Taking 1
u̇2 + v̇2
as the surface density of kinetic
http://arxiv.org/abs/0704.1305v2
energy, the equations of motion read
üi = ∂jσ
ui − vi
+ β e ilm
vlm + ulm
v̈i = ∂jσ
ui − vi
− β e ilm
vlm + ulm
with the sublattice 2-D stress tensors
= 2µ′ uij + λ′ δijukk + 2µ v
ij + λ δijvkk
−β eijk
uk − vk
= 2µ′ vij + λ′ δijvkk + 2µu
ij + λ δijukk
−β eijk
uk − vk
As expected, α determines the frequency of two degener-
ate k = 0 optical modes: ωΓ
2 = 4α.
First, we briefly show that the usual macroscopic elas-
tic energy of graphene and its Lamé coefficients can be
obtained from V by considering a static, uniform solu-
tion of Eqs. 5 with identical deformations on both lattices
with an internal displacement 2qi ≡ ui − vi :
2qi = ℓ e ilm u
lm = ℓ e ilm v
lm, (7)
where ℓ = β/α is a characteristic length. Anisotropic
(2γ = uxx−uyy) and shear (η = uxy) strains produce in-
ternal displacements qx = ℓγ and qy = −ℓη (Fig. 1). The
elastic energy for uniform deformations Wu =
Vuσg d
then simplifies to
Vu[u, q] =
uijuij +
+ 2αq2 − 2 β eijkuijqk, (8)
where σg = 2σs = 2.26 g cm
−2 is the surface density
of graphene, µR ≡ µ + µ′ − β
, λR ≡ λ + λ′ + β
the measurable Lamé coefficients [24]. Macroscopic prob-
lems do not distinguish between the two sublatices; elim-
inating qi in Eq. 8 through Eqs. 4 and 7 we obtain
the familiar, isotropic, macroscopic energy for graphene,
Vu = µRu
ijuij + λRu
j/2. In the long wavelength
limit Eqs. 5 returns the familiar longitudinal and trans-
verse speeds of sound in terms of the Lamé coefficients:
v2L = 2µR + λR, v
T = µR.
The out-of-plane displacements u⊥(x) and v⊥(x) do
not couple with the in-plane ui, vi in the harmonic limit:
invariance under simultaneous sign change of u⊥ and v⊥
prevents it, for flat sheets. Introducing 2p⊥(x) = u⊥(x)+
v⊥(x) and 2q⊥(x) = u⊥(x)−v⊥(x), V⊥ must be invariant
under p⊥ → p⊥ + L(x), L(x) a linear function in the
plane, and thus, can contain only second (and higher)
derivatives in p⊥. Symmetry dictates (cf. Appendix)
V⊥ = 4α⊥q
⊥ − 4α′⊥∂iq⊥∂iq⊥ + 4β⊥eijk ∂kq⊥∂ijp⊥
+ 2µ+
∂ijp⊥∂
ijp⊥ + λ
∂iip⊥∂
− 2µ−
∂ijq⊥∂
ijq⊥ − λ−⊥∂
iq⊥. (9)
The frequency of the k = 0 out-of-plane optical mode is
α⊥, and the out-of-plane acoustic branch is quadratic
at small wave-vector, as expected.
FIG. 2: Bicontinuum phonons compared to EELS data (dia-
monds [26] and squares [27]), fitting either to the entire Bril-
louin zone (top) or just around Γ along Γ → M .
The bicontinuum phonons are much more richly struc-
tured than in a traditional continuum model: they in-
clude all the optical branches, show nonlinear dispersion
at large wavevector, and even display the main features
of the Brillouin zone, all without sacrificing the advan-
tages of a continuum framework. Plane-wave solutions of
Eqs. 5 returns an analytically solvable fourth-order secu-
lar equation in ω(k), yielding two acoustic and two opti-
cal branches. The longitudinal branches cross at the ver-
tices of a hexagon. Since the two-field elastic energy den-
sity respects the point group symmetry of the graphene
lattice, this hexagon is oriented just as the graphene Bril-
louin zone; although the model, unlike in the envelope
function approach [25], has no built-in length scale, the
elastic parameters can be constrained so that the crossing
point coincides with the K point of graphene. A simi-
lar argument holds for the out-of-plane modes: strikingly
one can construct the correct Brillouin zone within a con-
tinuum model. Fig. 2 shows the bicontinuum phonons fit
to electron-energy-loss spectroscopy (EELS) data [26, 27]
for parameters fitted either to the full Brillouin zone or
just around Γ [28].
The bicontinuum provides a unified framework for nan-
otube mechanics which can describe all current compu-
tational results on the coupling of nanotube phonons
to static structural distortions, to each other (e.g.
breathing-to-Raman or longitudinal-to-transverse modes
in helical tubes) and to the tube electronic structure.
In a cylindrical geometry with coordinates {r, φ, z}, a
coupling between the tangential displacements ui, and
the radial ur = u⊥ appears in V of Eq. 1 via u
φ + ur
/r (and similarly for v); this accounts for
the emergence of the Radial Breathing Mode (RBM) [29].
We consider uniform solutions: u = uoe
−iωt, v = voe
−iωt.
The tube’s helicity can be subsumed into new axes {ξ, ζ}
(ξ = φ cos 3θc + z sin 3θc ζ = −φ sin 3θc + z cos 3θc) ro-
tated by an angle 3θc with respect to the base of the
tube. In terms of p, q we obtain pξ, pζ = 0 and
ω2 − 4α
pr = 0
ω2 − v
+β2/α
qζ = 0
ω2 − 4α
ω2 − 4α⊥ + 2µ−2µ
′+λ−λ′
. (10)
Unlike standard elasticity [17], which cannot describe op-
tical modes, or standard atomistic descriptions, which
cannot be solved analytically, the two-field continuum
model enables an exact analytical solution for the cou-
pling between the RBM and the graphite-like optical
mode through the first two of Eqs. in (10); the RBM
induces a shear in the sublattices, uφφ = vφφ = ur/r,
which couples with the internal displacement through
β, and vice versa. Thus, the RBM is not purely ra-
dial, but has a longitudinal component qzB ∼ ℓ2r cos 3θc,
as previously seen in a numerical calculation[30]. Ex-
pansion of the RBM frequency in powers of l/r re-
veals a correction to the the standard continuum re-
sult vL/r [17]: ωB =
. The
graphite-like optical modes of chiral tubes are ωξ =
ωζ/ωξ = 1 +
, also of mixed longitudi-
nal/transverse character except for armchair and zig-zag
nanotubes, while the out-of-plane optical mode ω⊥ =
4α⊥ − 2µ−2µ
′+λ−λ′
is purely radial. A density
functional theory calculation of the breathing mode [31]
reports different frequencies with (ωB) and without (ω̃B)
coupling to optical modes. We predict r2
ω̃2B − ω2B
β2/α as r → ∞: using ref [31] data for ω̃B, ωB we ob-
tain ℓ ≡ β/α = 0.25 Å (0.27 Å) for non metallic zig-zag
(armchair) tubes, in good agreement with the parameters
from our fit to the graphene phonons [28].
The bicontinuum can also describe electron-lattice cou-
pling to both acoustic and optical modes, by incorporat-
ing a tight-binding model whose nearest neighbor hop-
ping integrals t(1), t(2), t(3) are modulated by the in-plane
elastic deformations:
dt(l) = −τ ê(l)i ê
ij + τ ê
i/e (11)
where e is the inter-atomic distance and τ a parameter
to be determined [32]. For example, lattice deformations
open gaps in metallic tubes, and these gaps in turn affect
vibrational frequencies. If ǫc, ǫv are the conduction and
valence bands, we have to nearest neighbors
ǫc(k)
2−ǫv(k)2 =
t(l)+2
t(l)t(m) cos(k·a(n)), (12)
where a(n) ≡ e(l) − e(m), n(l,m) is cyclic in {1, 2, 3} (e.g.
a(3) ≡ e(1) − e(2)) and {e(i)} connects nearest neighbors.
From Eqs. 11,12 we find the band gap opened by strain
in a metallic nanotube to be
uijuij −
eijku
eijkφ̂
kuij φ̂hq
h − 1
eijku
ij φ̂k
In the second line of equation (13) the symmetry of the
honeycomb lattice is broken by the unit vectors φ̂i, ẑi of
the cylindrical coordinates. In terms of 2γ′ ≡ uφφ − uzz,
η′ ≡ uφz, qz , equation (13) reads
∆ = 3τ |qz/e+ γ′ cos(3θc) + η′ sin(3θc)| , (14)
which corrects and extends a well known previous result
within a one-field continuum model [7] that neglected the
inner displacement (i.e. qi = 0).
Opening bandgaps in metallic nanotubes causes several
shifts in observed quantities. The term proportional to
q2z in Eq. (13) show that longitudinal optical modes open
a bandgap in metallic tubes of any helicity; the elastic en-
ergy lowers by a term proportional to the square of the
bandgap, leading to a the softening of longitudinal optical
frequency in metallic nanotubes, as revealed by a recent
DFT study [33]. Eq. (13) predicts also a softening of the
RBM in metallic nanotubes δωB
= −A cos2 (3θc), high-
est for zig-zag tubes as seen in DFT [31], and relates it to
the optical softening, with A = (1−ℓ/e)ωoptδωopte2/4v2L,
ωopt the graphite-like optical mode, and δωopt its soften-
ing in metallic tubes (A ≃ 2%). Other shifts can be
predicted: the speed of sound for the twist mode softens
by ∆ct
= − v
A sin2 (3θc), or ≃ 2.2% in armchair tubes.
Doping-induced structural deformations can also be
studied by minimizing the total energy (elastic plus
doped electrons). Subtle phenomena absent in other
models [22] can be accessed within the bicontinuum
framework. Going to next-nearest-neighbor in the hop-
ping integrals (dt
1 = −τ1 â
ij [32]), we find that
at first order in both a/r and the number of dopant elec-
trons per atom ρe, semiconducting (n, 0) nanotubes show
doping-induced changes in tube length (dL/L = uzz) and
axial bond-length (dbax = eu
zz − qz):
dL/L = ρeτ
+ 3τ1
2µR+λR
µR+λR
dbax = ± ρeτ2mCω2opte
. (15)
where mC is the mass of the carbon atom. The sign is
positive (negative) for r = n mod 3 = 2 (n mod 3 = 1).
Recent DFT results [34] indeed show shrinking or stretch-
ing of bax for n = 16, 13 or n = 14, 11 tubes respec-
tively, as predicted by Eq. 15. In DFT, the overall tube
lengthens in the second case (n = 14, 11), again in ac-
cord with the bicontinuum; the lengthening found for
r = 2, is less than for r = 1, perhaps a consequence of
the change in sign in Eqs. 15. Finally the shrinking of
the axial bond determines an up-shift in the longitudinal
graphite-like optical mode and might explain recent Ra-
man results that point toward anomalous bond contrac-
tion under doping in semiconducting nanotubes [35, 36].
In summary, a symmetrized two-field continuummodel
of graphene and carbon nanotubes provides the first uni-
fied analytical treatment for a wide range of vibrational
and electromechanical phenomena including nonlinear
dispersion of in-plane phonons, zone-edge degeneracies
and optical modes. A full range of vibrational-electronic-
mechanical couplings, which were absent from previous
continuummodels or happened upon in an ad hoc fashion
in computational work, can now be understood within a
single unified analytical framework. Extending the for-
malism to include higher-order effects arising from curva-
ture or metallic character (i.e. symmetry breaking terms
containing φ̂i, ẑi, as in Eq. 13), anharmonicity (terms
higher order in uij , vij), or long-distance interactions
(higher partial derivatives) is straightforward. An ex-
tension to boron nitride nanotubes, with different coef-
ficients for each sublattice in the direct terms of Eq. 2,
might prove useful to study their piezoelectricity.
Appendix: Derivation of Eq. 3
The term c[u, v] must be invariant under the combina-
tion of 2π/6 rotations and the exchange of fields u ↔ v.
Adding reflection through the x axis (Fig. 1) then im-
plies C3v invariance. There is also a field translation
invariance: u(x) → u(x) + p, v(x) → v(x) + p. The
objects ui, vj , uij , and vij can be combined pairwise
only into tensors of rank two, three and four; thus c[u, v]
decomposes into three parts. The first part has terms
like uivj ; symmetry then implies the form α(u − v)2
with α > 0 to ensure an energy minimum. The sec-
ond part has terms like uijvkl; the only admissible form
is 2µuijvij + λu
j . The third part contains only rank
three terms such as uijvk contracted with a C3v invari-
ant tensor eijk, giving eijku
ijvk. By requiring invariance
under 2π/6 rotations conjugated with sublattice switch-
ing, and also the field translation invariance, we obtain
the form eijku
ij(uk − vk) + e∗ijkvij(vk − uk), where the
star means a 2π/6 rotation. Since C3v invariance implies
e∗ijk = −eijk we finally obtain the third row of Eq. 3.
[1] S. Iijima, Nature (London) 354, 56 (1991).
[2] R. Saito, G. Dresselhaus and M. S. Dresselhaus, Physical
properties of Carbon Nanotubes (Imperial College Press,
London 1998).
[3] M. S. Dresselhaus, G. Dresselhaus and P. C. Eklund, Sci-
ence of Fullerences and Carbon Nanotubes (Academic,
New York, 1996).
[4] S. J. Tans et al., Nature 386, 474 (1997).
[5] C. L. Kane et al., Europhys. Lett. 41, 683 (1998).
[6] A. Rochefort, P. Avouris, F. Lesage, D. R. Salahub, Phys.
Rev. B 60, 13824 (1999).
[7] L. Yang and J. Han, Phys. Rev. Lett. 85, 154 (2000)
[8] S. Berber, Y-K. Kwon, and D. Tománek, Phys. Rev. Lett.
84, 4613 (2000).
[9] H.-Y. Chiu et al., Phys. Rev. Lett. 95, 226101 (2005).
[10] P. Kim, L. Shi, A. Majumdar and P. L. McEuen Phys.
Rev. Lett. 87, 215502 (2001).
[11] E. Richter and K. R. Subbaswamy Phys. Rev. Lett. 79,
2738 (1997).
[12] R. Saito et al., Phys. Rev. B 64, 085312 (2001).
[13] A. Jorio et al., Phys. Rev. Lett. 86, 1118 (2001).
[14] A. Jorio, R. Saito, G. Dresselhaus and M. S. Dresselhaus,
Phil. Trans. R. Soc. Lond. A 362, 2311 (2004).
[15] K. S. Novoselov et al., Nature 438, 197 (2005).
[16] S. Y. Zhou et al., Nature Phys. 2, 595 (2006).
[17] G. D. Mahan, Phys. Rev. B 65, 235402 (2002).
[18] H. Suzuura and T . Ando, Phys. Rev. B 65, 235412
(2002); A. Raichura et al., J. Appl. Phys. 94, 4060 (2003).
S.V. Goupalov, Phys. Rev. B 71, 085420 (2005)
[19] F. Comas et al., Phys. Rev. B 47, 7602 (1993).
[20] L. Chico and R. Pérez-Álvarez Phys. Rev. B 69, 35419
(2004). L. Chico and R. Pérez-Álvarez, Phys. Rev. B 73,
075425 (2006)
[21] Y. N . Gartstein et al., Phys. Rev. B 68, 115415 (2003)
[22] M. Verissimo-Alves, B. Koiller, H. Chacham, and
R. B. Capaz, Phys. Rev. B 67, 161401 (R) (2003).
[23] H. Deresiewicz et al. The collected papers of Raymond D.
Mindlin, Springer-Verlag New York (1989).
[24] L. D. Landau, E. M. Lifshitz “Theory of Elasticity” Perg-
amon Press Oxford (1986). The density of elastic energy
for an isotropic system has the form f = µuijuij+
where λ, µ are the Lamé coefficients. Here we renormalize
the coefficients to σg = 1.
[25] B. A. Foreman, Phys. Rev. B 52, 12260 (1995).
[26] C. Oshima et al., Solid State Commun. 65, 1601 (1988).
[27] S. Siebentritt, R. Pues, K. H. Rieder, and A. M. Shikin,
Phys. Rev. B 55, 7927 (1997).
[28] The fit around Γ returns (in Km s−1) vL = 21, vT = 14,
−2µ+ 2µ′ − λ+ λ′ + β2/α
= 4.4, (µ− µ′)
15, 6 and ℓ ≡ β/α = 0.3 Å [28]. The fit to the
full zone uses vL = 16.5 Km s
−1, vT = 10.8 Km
2µ− 2µ′ + λ− λ′ − β2/α
= 8.7 Km s−1,
(µ− µ′)
= 6, 6 Km s−1, ℓ ≡ β/α = 0.24 Å; an exten-
sion to higher derivatives would improve the agreement.
[29] The chiral vector of the tube breaks the hexagonal sym-
metry and allows for new terms to be introduced in V
as curvature corrections, which for simplicity we won’t
consider here. Different problems will suggest different
leading corrections.
[30] E. Dobardžić et al., Phys. Rev. B 68, 045408 (2003).
[31] J. Kürti et al., New J. Phys.5, 125 (2003).
[32] t(l), t
1 are the absolute values of the hopping integrals
for nearest and next nearest neighbors. We assume they
depend only on distance, and thus dt = −τde/e, dt1 =
−τ1da/a.
[33] O. Dubay et al., Phys. Rev. Lett. 88, 235506 (2002).
[34] R. E. Margine et al. submitted to Phys. Rev. Lett.
[35] G. Chen, C. A. Furtado, U. J. Kim, and P. C. Eklund
Phys. Rev. B 72 155406 (2005).
[36] G. Chen, C. A. Furtado, S. Bandow, S. Iijima and
P. C. Eklund Phys. Rev. B 71 045408 (2005).
