Updated: January 23, 2008
Thermal Amplitudes of Spherical Phonon Modes

   Thermal average squared amplitudes of phonon modes in spherical coordinates are found (ignoring quantum mechanical corrections), with the following results:
Longitudinal spheroidal modes: |A|2 = dω kBT (2 l + 1) / (π2 ρm clm ω2)
Torsional (transverse) modes:    |A|2 = dω kBT (2 l + 1) / (l(l+1) π2 μm ctm)
Transverse spheroidal modes: |A|2 = dω kBT (2 l + 1) / (l(l+1) π2 ρm ctm ω2)
   At 300 K, phonon modes can be treated classically for wavenumbers below 200 cm-1, in other words, phonons of frequency f ≤ 6×1012 Hz.

   Apart from quantum mechanical corrections, a linear mechanical system in thermal equilibrium at temperature T (absolute temperature, Kelvin degrees) will have average energy kBT for each mechanical degree of freedom. kB is Boltzmann's constant, 1.381×10-23 J/K. This is commonly called the "equipartition theorem". Apply this to a linear elastic continuum solid (in particular, the glass "matrix" surrounding a nanoparticle) with density ρm, shear modulus μm, longitudinal speed of sound clm and transverse speed of sound ctm. Realistic properties of the material, such as anisotropy of its elasticity or dispersion of phonon modes at higher frequencies, are ignored in everything that follows. Let this object take the form of a sphere of radius Rm. For a coordinate system, choose the origin to be at the center of the sphere, and use spherical coordinates (r,θ,φ), related to Cartesian coordinates by x = r sin(θ) cos(φ), y = r sin(θ) sin(φ), z = r cos(θ).
   Disturbances in this solid are described in terms of a displacement field u(r,θ,φ,t) (more conveniently denoted as "u(r,t)"). The units of u are metres. u is the displacement r' - r of a point of material whose equilibrium position is r. For greatest simplicity, the boundary conditions are assumed to be that
u(Rm,θ,φ,t) = 0 for all θ and φ
that is, the surface of the sphere is rigidly fixed.
   Since the system is linear, such a general disturbance is the linear combination of normal modes of vibration of the form
u(r,t) = cos(ωt) u(r)
where ω is the angular frequency of the vibration (in radians per second).
   The next job is to come up with a orthogonal basis of functions corresponding to vibrational modes with angular frequency ω. Vibrational normal modes in an isotropic continuous solid can be longitudinal or transverse. Transverse vibrations have two possible directions of polarization.

Case I: Longitudinal ("Spheroidal") Modes
    The longitudinal modes (often called "spheroidal" modes) have the form (for integer l = 0, 1, 2, 3...)
u(r,t) = cos(ωt) ∇( A jl (klm r) Pl (cos(θ)) )
ω = angular frequency (radians/second)
∇ = gradient differential operator
A = amplitude (real number, in metres squared)
klm = ω / clm    ("l " = "longitudinal", "m" = "matrix")
clm = longitudinal speed of sound in matrix material
jl (.) = spherical Bessel function of first kind of order l       (See: mathworld.wolfram.com)
Pl (.) = Legendre polynomial of order l

Evaluating the gradient,
ur = cos(ωt) A (d/dr) jl (klm r) Pl (cos(θ))
uθ = cos(ωt) A (1/r) (d/dθ) jl (klm r) Pl (cos(θ))
In the large-r limit, only ur makes a significant contribution to the energy, so uθ will be ignored.
The velocity is
vr = - ω sin(ωt) A (d/dr) jl (klm r) Pl (cos(θ))
   A useful identity for spherical Bessel functions of the first kind is (from: mathworld.wolfram.com)
jn(x) = (-1)n (x)n ( (1/x) d/dx )n (sin(x)/x)
while for spherical Bessel functions of the second kind
nn(x) = (-1)n (x)n ( (1/x) d/dx )n (-cos(x)/x)
(In particular, j0(x) = sin(x)/x, j1(x) = sin(x)/x2 - cos(x)/x, n0(x) = -cos(x)/x, n1(x) = - cos(x)/x2 - sin(x)/x. )
   If x is large, so that the leading term goes as 1/x, then approximately
jn(large x) =   sin(x) / x n = 0, 4, 8, ...
jn(large x) = -cos(x) / x n = 1, 5, 9, ...
jn(large x) = -sin(x) / x n = 2, 6, 10, ...
jn(large x) =   cos(x) / x n = 3, 7, 11, ...
   In what follows, the sign of the velocity is not important. That is because we are only interested in the energy. Also, whether it is cos(x) or sin(x) is not important. That is because of the large r limit being taken. In this spirit, the approximate form for the velocity field r component is
vr = - ω sin(ωt) A (d/dr) ( cos(klm r)/(klm r) ) Pl (cos(θ))
Taking the r derivative, and keeping the dominant term only
vr = ω sin(ωt) A klm sin(klm r)/(klm r) Pl (cos(θ))
vr = ω sin(ωt) A sin(klm r)/(r) Pl (cos(θ))
   Kinetic energy density (joules per cubic metre) is uKE = ½ ρm v2. The time-averaged, local-space-averaged kinetic energy density is
(0.5) ρm ω2 0.5 A2 0.5 r -2 (Pl (cos(θ)))2
The time-averaged, local-space-averaged total energy density (including kinetic and potential energy) is
(0.5) ρm ω2 0.5 A2 r -2 (Pl (cos(θ)))2
The total energy of the vibrational mode is obtained by integrating this over the entire volume (0≤r≤Rm), which involves introducing the differential volume factor, dV = r 2 sin(θ) dθ dφ dr.
Let Fl represent the integral over θ, from 0 to π, of
Fl = ∫ sin(θ) ( Pl (cos(θ)) )2 dθ
Based on standard properties of Legendre polynomials (Introduction to Electrodynamics 3rd Edition, D. Griffiths, page 140), Fl = 2/(2l+1). I also verified this numerically with a short C++ program finfl.cpp.
The φ part of the integral simply introduces a factor of 2π.
The energy of the vibrational mode then is the integral from 0 to Rm of
U = ∫ (0.5) ρm (0.5) ω2 2π A2 Fl r 2 r -2 dr
which equals the integral from 0 to Rm of
U = ∫ (0.5) ρm π ω2 A2 Fl dr
which is:
U = (0.5) ρm π ω2 A2 Fl Rm
In thermal equilibrium, (ignoring quantum effects) this equals kBT, so therefore,
kBT = (0.5) ρm π ω2 A2 Fl Rm
and consequently,
A2 = kBT / ( (0.5) ρm π ω2 Fl Rm )
In particular, A varies as one over the frequency of the mode.
   The frequency of individual modes is determined by the boundary condition, u = 0 at r = Rm. This means that
sin(klm Rm) = 0      (if l is even)
cos(klm Rm) = 0      (if l is odd)
   In any case, the allowed frequencies are approximately of the form
ω Rm / clm = n π
where n is a positive integer.
ω = n π clm / Rm
   Consider a range of frequencies of width dω. The number of modes with frequencies in this range is
n = dω Rm / (π clm )
Suppose I want to represent all normal modes in the frequency range dω by a single function with amplitude A. I want this function to have a total energy n kBT. (Note that this function represents a combination of normal modes). The amplitude of this function is obtained from:
A2 = (dω Rm / (π clm )) kBT / ( (0.5) ρm π ω2 Fl Rm )
   It is at this point that Rm cancels out, and we are easily able to take the macroscopic (i.e. thermodynamic) limit that Rm is large. In that case,
A2 = dω kBT / ( (0.5) π2 ρm clm Fl ω2)
A2 = dω kBT (2 l + 1) / (π2 ρm clm ω2)
Note that there is frequency dependence of this amplitude. Also note that the spherical symmetry of the problem means that there are actually 2l + 1 vibrational modes of this type, corresponding to different values of the azimuthal quantum number m = -l, .., l.

Case II: Torsional Modes
    There are two polarizations of transverse modes. The displacement can be polarized along the θ direction or the φ direction. One class of transverse modes ("torsional modes", polarized along the φ direction) has the form (for integer l = 1, 2, 3...)
u(r,t) = cos(ωt) ∇×( r A jl (ktm r) Pl (cos(θ)) )
A = amplitude (real number, in metres)
ktm = ω / ctm
ctm = transverse speed of sound in matrix material

   u has components ur, uθ, and uφ. Making use of the formula for curl in spherical coordinates, ur = 0, uθ = 0 and
uφ = cos(ωt) (-1/r) (d/dθ) ( r A jl (ktm r) Pl (cos(θ)) )
= cos(ωt) (-1/r) r A jl (ktm r) (d/dθ) Pl (cos(θ))
= - cos(ωt) A jl (ktm r) (d/dθ) Pl (cos(θ))
   In thermal equilibrium, the energy of any single mode of this type (on average) will be kBT. This condition will ultimately be shown here to determine the rms value of the magnitude of A.
   The first step is to calculate the total kinetic energy associated with this mode. (the potential energy is identical) The velocity field for this mode is v = (d/dt) u. The φ component of velocity is
vφ = ω sin(ωt) A jl (ktm r) (d/dθ) Pl (cos(θ))
   A useful identity for spherical Bessel functions of the first kind is (from: mathworld.wolfram.com)
jn(x) = (-1)n (x)n ( (1/x) d/dx )n (sin(x)/x)
while for spherical Bessel functions of the second kind
nn(x) = (-1)n (x)n ( (1/x) d/dx )n (-cos(x)/x)
(In particular, j0(x) = sin(x)/x, j1(x) = sin(x)/x2 - cos(x)/x, n0(x) = -cos(x)/x, n1(x) = - cos(x)/x2 - sin(x)/x. )
   If x is large, so that the leading term goes as 1/x, then approximately
jn(large x) =   sin(x) / x n = 0, 4, 8, ...
jn(large x) = -cos(x) / x n = 1, 5, 9, ...
jn(large x) = -sin(x) / x n = 2, 6, 10, ...
jn(large x) =   cos(x) / x n = 3, 7, 11, ...
   In what follows, the sign is not important. Also, whether it is cos(x) or sin(x) is not important. In this spirit, the approximate form for the velocity field φ component is
vφ = ω sin(ωt) A (1/(ktm r)) cos(ktm r) (d/dθ) Pl (cos(θ))
vφ = sin(ωt) A (ctm/ r) cos(ktm r) (d/dθ) Pl (cos(θ))
vφ = sin(ω t) A (ctm/ r) cos(ktm r) (d/dθ) Pl (cos(θ))
   Kinetic energy density (joules per cubic metre) is ½ ρm v2. The time and local space averaged kinetic energy density is
(0.5) ρm (0.5) (0.5) A2 ctm2 r -2 ( (d/dθ) Pl (cos(θ)) )2
The total time and local space averaged energy density which includes potential energy, is simply twice this:
(0.5) ρm (0.5) A2 ctm2 r -2 ( (d/dθ) Pl (cos(θ)) )2
The total energy of the vibrational mode is obtained by integrating this over the entire volume (0≤r≤Rm), which involves introducing the differential volume factor, dV = r 2 sin(θ) dθ dφ dr.
Let El represent the integral over θ, from 0 to π, of
El = ∫ sin(θ) ( (d/dθ) Pl (cos(θ)) )2 dθ
Numerically, I found using a C++ program findel.cpp that El = 2 l (l+1)/(2l+1).
The φ part of the integral simply introduces a factor of 2π.
The energy of the vibrational mode then is the integral from 0 to Rm of
U = ∫ (0.5) ρm (0.5) 2π A2 ctm2 El r 2 r -2 dr
which equals the integral from 0 to Rm of
U = ∫ (0.5) ρm π A2 ctm2 El dr
which is:
U = (0.5) ρm π A2 ctm2 El Rm
Finally, since ctm = sqrt(μm/ρm), where μm is the shear modulus of the matrix material (in Pascals), the energy of the mode is
U = (0.5) μm π A2 El Rm
In thermal equilibrium (ignoring quantum effects), U equals kBT, so therefore,
kBT = (0.5) μm π A2 El Rm
and consequently,
A2 = kBT / ( (0.5) μm π El Rm )
In particular, A does not vary depending on the frequency of the mode.
   The frequency of individual modes is determined by the boundary condition, u = 0. This means that
sin(ktm Rm) = 0      (if l is even)
cos(ktm Rm) = 0      (if l is odd)
   In any case, the allowed frequencies are approximately of the form
ω Rm / ctm = n π
where n is a positive integer.
ω = n π ctm / Rm
   Consider a range of frequencies of width dω. The number of modes with frequencies in this range is
n = dω Rm / (π ctm )
Suppose I want to represent all normal modes in the frequency range dω by a single function with amplitude A. I want this function to have a total energy n kBT. (Note that this function represents a combination of normal modes). The amplitude of this function is obtained from:
A2 = (dω Rm / (π ctm )) kBT / ( (0.5) μm π El Rm )
   It is at this point that Rm cancels out, and we are easily able to take the macroscopic (i.e. thermodynamic) limit that Rm is large. In that case,
A2 = dω kBT / ((0.5) π2 μm ctm El)
A2 = dω kBT (2l+1) / (l(l+1) π2 μm ctm)
Unlike Case I and Case III, there is no frequency dependence of this amplitude. Note that there are 2l + 1 modes of this type because of spherical symmetry.

Case III: Transverse Spheroidal Modes
    The other class of transverse modes ("transverse spheroidal modes") has the form (for integer l = 1, 2, 3...)
u(r,t) = cos(ωt) ∇×∇×( r A jl (ktm r) Pl (cos(θ)) )
"∇×∇×" signifies taking the curl twice in succession
A = amplitude (real number, in metres squared)
ktm = ω / ctm

First, define D (with components Dr, Dθ, Dφ) as
D(r) = r A jl (ktm r) Pl (cos(θ))
and also E(r) = ∇ × D
so that u(r) = ∇ × E
Dr = r A jl (ktm r) Pl (cos(θ))      Dθ = 0      Dφ = 0
The formula for curl in spherical coordinates is, respectively:
Er = (1/(r sin(θ)) [ (d/dθ) [ ( sin(θ) Dφ ] - (d/dφ) D&theeta; ]
Eθ = (1/r) [ (1/sin(θ)) (d/dφ) Dr - (d/dr) [ r Dφ ] ]
Eφ = (1/r) [ (d/dr) [ r Dθ ] - (d/dθ) Dr< ]
ur = (1/(r sin(θ)) [ (d/dθ) [ ( sin(θ) Eφ ] - (d/dφ) E&thheta; ]
uθ = (1/r) [ (1/sin(θ)) (d/dφ) Er - (d/dr) [ r Eφ ] ]
uφ = (1/r) [ (d/dr) [ r Eθ ] - (d/dθ) Err ]
   In this case, since only Dr is nonzero, and depends on r and θ only,
Eφ = (-1/r) (d/dθ) Dr
Eφ= (-1) A jl (ktm r) (d/dθ) ( Pl (cos(θ)) )
   Again making use of the formula for curl in spherical coordinates, uφ=0, and
ur = cos(ωt) (1 / (r sin(θ)) ) (d/dθ) [ sin(θ) (-1) A jl (ktm r) (d/dθ) ( Pl (cos(θ)) ) ]
ur = - cos(ωt) ( r -1 A jl (ktm r) (1 / sin(θ) ) (d/dθ) [ sin(θ) (d/dθ) Pl (cos(θ)) ] )
uθ = cos(ωt) (-1/r) (d/dr) ( r (-1) A jl (ktm r) (d/dθ) ( Pl (cos(θ)) ) )
uθ = cos(ωt) (1/r) (d/dr) [ r A jl (ktm r) ] (d/dθ) [ Pl (cos(θ)) ]
   Note that ur depends on r, θ and t. The r-dependent factor is r -1 jl (ktm r). When r is large, this approaches
(±1)(1/ktm) sin(ktm r) r -2 (if l is even)
or
(±1) (1/ktm) cos(ktm r) r -2 (if l is odd).
   uθ also depends on r, θ and t. The r-dependent factor is (1/r) (d/dr) [r jl (ktm r)]. When r is large, this approaches
(±1) cos(ktm r) r -1 (if l is even)
or
(±1) sin(ktm r) r -1 (if l is odd).
   Clearly the ur component dies out more quickly and only the uθ component makes a significant contribution to the energy of the mode after the "thermodynamic limit" of large Rm is taken. This also illustrates the θ-polarized transverse character of these modes. The velocity is,
vθ = - ω sin(ωt) (1/r) (d/dr) [ r A jl (ktm r)] (d/dθ) [ Pl (cos(θ)) ]
Which, in the case where r is large, is "approximately"
vθ = (±1) A ω sin(ωt) (1/r) (d/dr) [ r sin(ktm r) (1/(ktm r)) ] (d/dθ) [ Pl (cos(θ)) ]    (even l )
vθ = (±1) A ω sin(ωt) (1/r) (d/dr) [ r cos(ktm r) (1/(ktm r)) ] (d/dθ) [ Pl (cos(θ)) ]    (odd l )
For brevity, I will only consider the even l case in the steps that follow.
vθ = (±1) A (1/(ktm)) ω sin(ωt) (1/r) (d/dr) [ sin(ktm r) ] (d/dθ) [ Pl (cos(θ) ]
vθ = (±1) A (1/(ktm)) ω sin(ωt) (1/r) ktm cos(ktm r) (d/dθ) Pl (cos(θ))
vθ = (±1) A ω sin(ωt) (1/r) cos(ktm r) (d/dθ) Pl (cos(θ))
   Kinetic energy density (in J/m3) is uKE = ½ ρm v2. The time-averaged, local-space-averaged kinetic energy density (denoted by < uKE > ) is
< uKE > = (0.5) ρm A2 ω2 (0.5) r -2 (0.5) [(d/dθ) Pl (cos(θ))]2
The time-averaged, local-space-averaged total energy density is twice < uKE >.
< uE > = (0.25) ρm A2 ω2 r -2 [(d/dθ) Pl (cos(θ))]2
   The total energy of the vibrational mode is obtained by integrating < uE > over the entire volume (0≤r≤Rm), which involves introducing the differential volume factor, dV = r 2 sin(θ) dθ dφ dr.
Let El represent the integral over θ,
∫0 → π   [ (d/dθ) Pl (cos(θ)) ]2 sin(θ) dθ
which comes out to 2 l (l+1) / (2l+1).
   The φ part of the integral simply introduces a factor of 2π.
   The energy of the vibrational mode, U, then is the integral from 0 to Rm of
U = ∫0→Rm (0.25) ρm A2 ω2 r -2 El r 2 (2π) dr
which equals the integral from 0 to Rm of
U = ∫0→Rm(π/2) ρm A2 ω2 El dr
which is easy to integrate, since the function being integrated is constant:
U = (π/2) ρm A2 ω2 El Rm
   In thermal equilibrium, (and assuming T is high enough to ignore quantum effects) U equals kBT, so therefore,
kBT = (π/2) ρm A2 ω2 El Rm
and consequently,
A2 = kBT / ( (π/2) ρm ω2 El Rm )
   The frequency of individual modes is determined by the boundary condition, u = 0. This means that,
cos(ktm Rm) = 0      (if l is even)
sin(ktm Rm) = 0      (if l is odd)
   In any case, the allowed frequencies are determined by
ω Rm / ctm = ( n + ½ ) π      (if l is even)
ω Rm / ctm = n π      (if l is odd)
where n is a positive integer.
   Consider a range of frequencies of width dω. The number of modes n with frequencies in this range is
n = dω Rm / (π ctm )      (l even or odd)
   Suppose (for example, for the purposes of doing a Fourier transform) I want to represent all normal modes in the frequency range dω by a single function with amplitude A. I want this function to have a total energy n kBT. (Note that this function represents a combination of normal modes). The amplitude of this function is obtained from:
A2 = (dω Rm / (π ctm )) (kBT / ((π/2) ρm ω2 El Rm))
   It is at this point that Rm cancels out, and we are easily able to take the macroscopic (i.e. thermodynamic) limit that Rm is large. In that case,
A2 = dω kBT / (π2 ctm ω2 (0.5) ρm El)
A2 = dω kBT (2 l + 1) / (l(l+1) π2 ρm ctm ω2 )
There are 2l + 1 modes of this type due to spherical symmetry.

Lattice and Quantum Mechanical Effects:
    The applicability of these results is certainly limited to frequencies sufficiently low that dispersion is not significant. Dispersion (variation of the sound velocity) happens when the wavelength of the vibrations becomes so small that it is comparable to the distance between atoms in the crystal. Furthermore, at high frequencies quantum mechanical effects will become important. Classical statistical mechanics is valid as long as kBT is much greater than hω/(2π), where h is Planck's constant (h=6.626×10-34 J-s). For a sample at room temperature (T=300K), kBT = 4×10-21 J, so quantum mechanical effects become important for ω > 4×1013 rad/s, or 200 cm-1. Since Raman scattering experiments with nanoparticles typically deal with vibrational frequencies much less than this, a quantum mechanical treatment of the phonon amplitudes is not necessary.

Daniel Murray
Associate Professor
Math, Stats & Physics Unit
University of British Columbia - Okanagan
Kelowna, BC, Canada
daniel "dot" murray "at" ubc "dot" ca

For a list of related articles click here.

Hosted by www.Geocities.ws

1