Updated: Feb 15, 2003
Thermal Displacement of an Isotropic Elastic Sphere
Embedded in an Infinite Isotropic Elastic Matrix

........

1. Light scattering from nanoparticles
   The calculation presented here is motivated by numerous experiments involving nanospheres of various materials (Si,Ag,Au,CdS,CdSe) embedded in various kinds of glass (SiO2,GeO2). Low frequency inelastic Raman scattering (LOFIRS) reveals peaks whose frequency varies as (1/d) where d is the diameter. Pump probe femtosecond spectroscopy reveals a "ringing" at a frequency that scales as (1/d). In principle, millimeter wavelength infrared absorption could also detect these modes.

2. Raman scattering
   Raman scattering is related to time variation of the polarizability tensor, α, of the target object being studied. (Note that α is not a function of position, but is rather a property of the target object as a whole.)

p = α E
where
p = dipole moment vector induced by the electric field   (Coulomb-metres)
α = tensor of order 2    (with components αxx, αyy, αxy etc.)
E = electric field vector of incident beam      (Volts/metre)

It is the induced dipole moment p that generates the scattered radiation.
   The target object under study could be an atom, a molecule, an isolated nanoparticle or even a bulk sample. If α does not vary with time, then the scattered radiation and the incident radiation will have exactly the same frequency. This is called "Rayleigh scattering." If α varies sinusoidally in time at some frequency ω1, the scattered radiation can have two possible frequencies:
ωinc + ω1   (called "Stokes")
ωinc - ω1   (called "anti-Stokes")
   In either case, this is called "Raman scattering." An experiment to observe Raman scattering simply requires a high quality filter that rejects all scattered radiation with the same frequency as the incident beam.
   The polarizability tensor α can vary for different reasons. The simplest situation is a diatomic molecule (like N2 or O2) where α rotates as the molecule rotates. If the molecule rotates 200 times a second, a given component of α will rise and fall 400 times a second. Thus, the Raman shift is twice the frequency of rotation. Since quantum mechanics quantizes the allowed frequencies with which an object can rotate, a discrete Raman spectrum can be observed.
   If an object (atom, molecule or nanoparticle) undergoes mechanical vibrations, α will oscillate at the same frequency.
   Another important situation is where an object (atom, etc.) undergoes an electronic transition which causes α to change. The frequency of the oscillation of α is simply related to the energy difference between the two electronic levels.
   In some situations, an object can vibrate at some frequency without having a corresponding peak in the Raman spectrum. Some modes are not "Raman active". This is related to the symmetry of the situation.

3. Basic parameters
   The situation under consideration here is a sphere with isotropic elasticity (the "nanoparticle") embedded in an infinite isotropic elastic medium (the "matrix"). Their properties are specified by parameters in the table below.

Table I.     Properties of nanoparticle and matrix
propertyunitsnanoparticlematrix 
density kg/m3 ρpρm
longitudinal
sound velocity 		m/s Clp Clm Clp=sqrt((λp+2μp)/ρp)
transverse
sound velocity 		m/s Ctp Ctm Ctp=sqrt(μp/ρp)
1st Lamé constant Paλp λm
Shear modulus Pa μp μm
radius mRpRmRm → ∞
diameter mDp Dp = 2Rp

   Let this composite object take the form of a sphere of radius Rm. For simplicity, we will want to eventually take the limit Rm → ∞ and in any case assume that Rp << 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(θ).
   Mechanical vibrations 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.


Sorry, but
This page is still under construction.
   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.    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. 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&thetta; ]
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&theeta; ]
uθ = (1/r) [ (1/sin(θ)) (d/dφ) Er - (d/dr) [ r Eφ ] ]
uφ = (1/r) [ (d/dr) [ r Eθ ] - (d/dθ) Er< ]
   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.

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