Molecular Dynamical Estimates of Vibrational Frequencies of a Crystalline Sphere

Last updated: January 21, 2008
Molecular Dynamical Estimates of Vibrational Frequencies of a Crystalline Sphere

Frequencies of vibrational modes of a face centered cubic crystal sphere are obtained from direct Fourier analysis molecular dynamics simulations of clusters of up to 12000 atoms. The symmetries of the modes are classified, and the frequencies are compared to those earlier predicted for an isotropic continuum elastic sphere. The finite size dependence of the frequencies is also reported.

   The classic problem of the vibrational modes of an isotropic elastic sphere has been extensively studied. [Lamb 1882] [Love 1944] [Sato and Usami 1962] [Auld 1973] [Nishiguchi and Sakuma 1981] [Bastrukov 1994] [Graff 1994] [Bastrukov et al. 1998] [McDaniel and Holt 2000 pdf ::] [Ye 2000 pdf] [Patton and Geller 2001 pdf ::]. This problem finds diverse application in geophysics [geophysics introduction], seismology [simple introduction to earthquake waves], solid stars [Bastrukov 1996], silicon nanoparticles [Patton and Geller 2001 pdf ::] and aqueous foam spheres [McDaniel and Holt 2000 pdf ::], and an accurate list of the exact frequencies has recently been reported [Murray 2002].
   Consider a freely suspended isotropic elastic sphere not subject to body forces. The sphere has radius R, uniform density ρ, shear modulus μ and Poisson ratio ν. These parameters fully describe the problem. If the sphere is disturbed, it will begin to vibrate. A given eigenfunction has frequency ω (in radians per second). The speed of transverse waves (shear waves) is

C_s = sqrt(μ/ρ)
   It is convenient to express normal mode frequencies in terms of the dimensionless variable

η = ω R / C_s
   Elastic materials in the real world are made of conglomerations of atoms (albeit not necessarily as periodic crystals). So it is plausible to attempt to calculate the properties of an elastic object by assembling a cluster of N atoms in the shape of the object, assuming an potential energy function (such as the Lennard Jones potential) that determines the inter atomic forces, and then carrying out a molecular dynamics calculation in which a computer calculates the forces, velocities and position of each atom as a function of time.

Table I. Linear Elastic Properties of Molecular Dynamics Solid:
Property	Symbol	Average value	Unit	Min. value	Max. value	Relation:
Mass of one atom	m_atom	1.0000000	kg	(exact)
Density	ρ	1.0000000	kg/m³	(exact)
Atomic force constant	dF/dr	1.0000000	N/m	(exact)
Atom center-center distance	d_cc	1.1224621	m	(exact)		d_cc=(sqrt(2)m_atom/ρ)^1/3
Bulk Modulus	B	0.8399473	Pa	(exact)		2ρ(d_cc)²(dF/dr)/(3m_atom)
Young's Modulus	Y	1.215	Pa	1.172	1.261	(use md4ym.cpp)
Poisson's ratio	ν	0.2589	- - -	0.250	0.267	ν=(3B-Y)/(6B)
Shear Modulus	μ	0.4826	Pa	0.4623	0.5045	μ=3BY/(9B-Y)
1st Lamé constant	λ	0.5182	Pa	0.5036	0.5317	λ = B - (2/3)μ
Shear wave speed	C_s	0.6947	m/s	0.6799	0.7103	C_s=sqrt(μ/ρ)
Longitudinal wave speed	C_l	1.2179	m/s	1.2068	1.2299	C_l=sqrt((λ+2μ)/ρ)
C_s/C_l		0.5704	- - -	0.5634	0.5775

   If, while this calculation is being performed, the position or velocity of one or more of the atoms is Fourier analyzed for its frequency content, it is possible to determine the spectrum of frequencies of the vibrational modes of the object.
   This calculation is approximate in two ways. First, the number of atoms that can be stored in the memory of a computer (13000 on mine) means that a truly macroscopic object cannot be simulated. Therefore, a sequence of objects of the same shape but with increasing numbers of atoms must be studied and the limit as N becomes macroscopic must be extrapolated. The second approximation is that the simulated object is isotropic. A perfect crystal (without defects or grain boundaries) does not have isotropic elastic properties. For example, the speed of sound depends on the direction of propagation within the crystal. However, the extent to which a face centered cubic crystal deviates from isotropicity is small enough that the approximation is worthy of some hope.
   In the molecular dynamics simulation, the N atoms are initially arranged in a face centered cubic arrangement. Thinking of each atom as a sphere, it is in contact with 12 other atoms. This arrangement is conjectured to have the highest possible density of atoms per unit volume. The conventional unit cell (the basic pattern that is repeated) has a cubic shape, and there are four atoms per unit cell. (The primitive unit cell is not cubic and has one atom in it.) The symmetry of the lattice is the same as the symmetry of a cube. This configuration of atoms is common to metals of many elements.
In the calculation, distance units are chosen so that the density of atoms is 1 atom per cubic metre. The center to center distance between two atoms is 1.1224621 (i.e. 2^1/6). The mass of each atom is assumed to be 1 kilogram. The strength of the neighbour to neighbour atom force is such that dF/dr is 1 N/m. The forces between atoms that are initially not nearest neighbours are ignored to allow the calculation to proceed quickly. In any case, those forces would be small.
   If the velocities of the atoms became too large, the forces would become nonlinear. Since we are interested in linear elastic behavior, the velocities of the atoms are kept small enough that the forces are linear. Normally, in these simulation, the peak speed of an atom is on the order of 1x10^-5 m/s. If the speed is made too small, there are numerical precision problems.
   The sharpness of the Fourier spectrum peaks depends on the amount of time over which the integrations are made.

Figure 1. Dimensionless frequency parameter η for vibrational modes of a face centered cubic crystal sphere are shown. The vertical axis is the reciprocal of the sphere's radius. The upper edge of the chart (above the black band) corresponds to infinite R. The bottom edge of the chart corresponds to R=5. The maximum R and maximum number of atoms used are plotted at the lower right. (2ran.gif)

From Figure 1, it is apparent that there are no vibrational modes with η less than about 2.3 (after extrapolating to infinite R at the top). The white line at η=0 corresponds both to translational and rotational zero frequency modes. As radius varies, a given mode appears as a sloped but approximately straight line on the graph. In every case, η increases as R gets bigger. This is likely a result of decreased phase velocity of phonons as the boundary of the Brillouin zone is approached, in other words, when the wavelengths are comparable to atomic distances. Thus, the finite size dependence of the mode frequencies is of the form

η(R) = η(infinity)(1 - D / R )
ω(R) = (η(infinity) C_s / R)( 1 - D / R )
where constant D depends on the mode, and C_s is the speed of transverse waves. Some modes have stronger finite size dependence than others. Apparently, D is always positive. For most modes, D is on the order of 1 or less. But remember that I am using microscopic distance units. If R was in metres, then D would be on the order of 0.3 nm or so.
For R=5 (a cluster of 520 atoms) the lowest vibrational mode is at η=2.10. Researchers studying nanoparticles normally adapt results from continuum elasticity theory without considering finite size effects [Cerullo, De Silvestri and Banin 1999 pdf ::] [Geller et al. 2001 pdf ::] [Patton and Geller 2001 abstract, preprint pdf ::] [Portales et al. 2002 pdf ::]. For a nanoparticle of 520 atoms, the resulting error in frequency would be on the order of 10%.
Based on the approximate spherical symmetry (approximate because a cubic crystal is not truly isotropic), the modes can be classified [Ye 2000]. Modes are of two types: spheroidal and torsional (seismology literature calls the latter "toroidal"). In spheroidal modes, the velocity field is of the form [Ye 2000]

u(r,θ,φ,t) exp(-i η C_s t / R ) =

A_l ∇

[

j(l, k_l r ) P(l,cos(θ))

]

+ B_l ∇× ∇×

[

r j(l, k_s r ) P(l,cos(θ))

]

where r is the position vector.
In order to identify the spheroidal modes of order n, the initial velocities of the atoms should be assigned according to the velocity field expected. It is sufficient and more convenient simply to use either the zero curl or the zero divergence term alone, since B_l/A_l for the lattice mode is unknown at the outset. The value of k_l is chosen to match the frequency predicted by continuum theory. The results are not affected as long as k_l is approximately in the correct range. The orientation of the crystal lattice is randomized each time. The subsequent motion is then Fourier analyzed to give the frequency. In some cases, three peaks are seen, instead of one. The three peaks are due to the broken spherical symmetry of the fcc lattice.

Figure 2. Dimensionless frequency parameter η for spheroidal modes of a face centered cubic crystal sphere are shown. The vertical axis is the reciprocal of the sphere's radius. The upper edge of each chart (above the black band) corresponds to infinite R. The maximum R and maximum number of atoms used are plotted at the lower right.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)

The above frequencies may be compared to the corresponding exact values obtained for an isotropic elastic continuum sphere [Murray 2002 link].

Table II: Spheroidal Modes
l	η	η (2nd peak)	η (3rd peak)	Exact η	Error	figure
0	4.55	-	-	4.528	+0.5%	see Fig 2(a)
0	10.77	-	-	10.633	+1.3%	see Fig 2(b)
1	3.43	4.10	4.19	3.444	-0.4%	- - -
1	6.69 (?)	-	-	6.836	-2.2%	- - -
2	2.48	2.63	3.00	2.642	-0.5%	- - -
2	4.73	4.91	5.26	4.892	+0.4%	- - -
2	-	-	-	8.369	-	- - -
3	3.75	4.06	4.24	3.921	+3.5%	- - -
4	5.08	5.22	-	5.017	+1.3%	- - -
5	5.98	6.07 (?)	-	6.043	-1.1%	- - -

Earlier [Murray 2002a] the torsional vibrational mode of a sphere of atoms was excited by giving atoms their initial velocities (v_x,v_y,v_z) based on their initial position (x,y,z) as follows:

v_x = 0.0001 z y
v_y = - 0.0001 z x
v_z = 0

This mode could be called the "wet dog" mode because of its similarity to the motion that a dog uses to shake water off of its fur. The head and body of the dog rotate in opposite directions. In the "wet dog" mode, the z>0 end of the sphere rotates clockwise while the z<0 end rotates counterclockwise. The z=0 plane remains approximately stationary. (An animated GIF image illustrates this [Murray 2002a] ).
The exact expression for the displacement field u for torsional modes is

u(x,y,z,t) exp(-i η C_s t / R ) =

C_l ∇×

[

(x,y,z) j(l, k_s r ) P(l,cos(θ))

]

In Figure 3 the frequencies that result from this initial excitation are shown.

Figure 3. Dimensionless frequency parameter beta for torsional modes of a face centered cubic crystal sphere are shown. The vertical axis is the reciprocal of the sphere's radius. The upper edge of the chart (above the black band) corresponds to infinite R. The bottom edge of the chart corresponds to R=5 in each case. The maximum R and maximum number of atoms used are plotted at the lower right. (md4mod3.cpp)

(a)
(2wetdog.gif)
dbeta=0.05

Table III: Toroidal Modes
l	η	η (2nd peak)	η (3rd peak)	Exact η	Error	figure
2	2.36	2.45	2.59	2.502	-2.1%	--

Appendix: Practical Application of Results:

   To read off frequencies from the above charts, you must first know (1) the radius of the sphere, R, in meters; (2) the density of the material, ρ, in kg/m³; (3) the shear modulus of the material, μ, in Pascals; (4) the Poisson ratio of the material. From the red scale, note the values of η at which the white peaks cross the top scale. Then, convert to frequency using,

f = 0.15915 η sqrt(μ/ρ) / R (f in Hertz)
ω = η sqrt(μ/ρ) / R (ω in radians/second)
   For practical engineering situations involving ordinary materials for which the Young's modulus (Y or sometimes E, also known as tensile modulus) is known, and only an approximate frequency is needed, then the following formula can be used, together with the charts:

f[Hz] = 100 η sqrt(Y[GPa] / density[g/cc] ) / radius[m]
(This formula makes the reasonable assumption that the Poisson ratio is 0.30.)
   For example, for an aluminum ball bearing of diameter 2 cm (R=0.01 m), elastic modulus 73 GPa, Poisson ratio 0.33 and density 2.7 g/cc, we refer to Figure 2(a) and find a frequency mode at η = 4.55. (Interpolate to the top edge of the chart since R in this case is approximately one hundred million, and approximately infinite.) The frequency in Hertz is

f = (100)(4.55) sqrt(73/2.7) / 0.01 = 236587 Hz
   The crystalline structure of many elements is face centered cubic such as Aluminum, solid Argon, Calcium, Copper, Gold, solid Krypton, Lead, solid Neon, Nickel, Platinum, Silver, Strontium, and solid Xenon. Many other metals are body centered cubic, or hexagonal close packed, which would have somewhat different frequencies, for example Iron is body centered cubic with lattice parameter 0.287 nm [C. Kittel, "Introduction to Solid State Physics" Wiley 4th Ed. 1971, see table 4 on page 38].

References:

H. Lamb, "On the vibrations of an elastic sphere" Proc. London Math. Soc. 13, 189 (1881-1882).
A. E. H. Love, A Treatise on the Mathematical Theory of Elasticity, (Dover, New York, 1944).
Y. Sato and T. Usami, Geophys. Mag. 31, 15 (1962)
C. Kittel, "Introduction to Solid State Physics" Wiley 4th Ed. 1971, see table 4 on page 38
B. A. Auld, Acoustic Fields and Waves in Solids, (John Wiley & Sons, New York, 1973).
N. Nishiguchi and T. Sakuma, Sol. Stat. Comm. 38, 1073 (1981).
Sergey I Bastrukov, Phys. Rev. E 49, 3166 (1994).
F. K. Graff, Wave Motion in Elastic Solid, (Ohio University Press, Ohio, 1994).
S. Bastrukov, Phys. Rev. E (1996)
Bastrukov et al., Physica A 250, 435 (1998).
Zhen Ye, "On the Low Frequency Elastic Response of a Spherical Particle" Chinese Journal of Physics Vol. 38, pages 103-110, (2000). pdf
J. Gregory McDaniel and R. Glynn Holt "Measurement of Aqueous foam Rheology by acoustic levitation" Phys. Rev. E Rapid Communications volume 61 R2204 (2000). pdf ::
Kelly Patton & Michael Geller "Phonon Spectrum in a Nanoparticle mechanically coupled to a substrate" Journal of Luminescence 94-95 (2001) 747-570 pdf ::
Michael R. Geller et al. "Theory of electron-phonon dynamics in insulating nanoparticles" Preprint 11 July 2001 t / d where v_t is the bulk transverse sound velocity." > pdf ::
Kelly R. Patton, Michael R. Geller, "Phonon Spectrum in a Nanoparticle Mechanically Coupled to a Substrate" (Preprint Tues, June 26, 2001) abstract, preprint pdf ::
G. Cerullo, S. De Silvestri and U. Banin "Size-dependent dynamics of coherent acoustic phonons in nanocrystal quantum dots" Phys. Rev. B volume 60 (July 15, 1999) pdf ::
H. Portales, L. Saviot, E. Duval, M. Gaudry, E. Cottancin, M. Pellarin, J. Lermé and M. Broyer "Resonant Raman Scattering by Quadrupolar Vibrations of Ni-Ag Core-shell nanoparticles" (Preprint Mar 22 2002) pdf ::
Daniel B. Murray, "Molecular Dynamics Simulation of an Elastic Material" 2002 (link to article)
Daniel B. Murray, "Vibrational Frequencies of an Elastic Sphere" 2002 (link to article)

Properties of InAs: Source http://www.ioffe.rssi.ru
Bulk modulus 5.8�1011 dyn cm-2
Density 5.68 g cm-3
Elastic constants at 293 K
C11=8.34�1011 dyn/cm2
C12=4.54�1011 dyn/cm2
C44=3.95�1011 dyn/cm2
Temperature dependences of elastic constants.
(Burenkov et al. [1975]).
For T= 200K Bulk modulus (compressibility-1) Bs= 5.81�1011 dyn/cm2
Shear modulus C'= 1.90�1011 dyn/cm2
[100] Young's modulus Yo= 5.14�1011 dyn/cm2
[100] Poisson ratio so = 0.35
Acoustic Wave Speeds
Wave propagation Direction Wave character Expression for wave speed Wave speed (in units of 105 cm/s)
[100] VL (C11/? )1/2 3.83
VT (C44/? )1/2 2.64
[100] Vl [(C11+C12+2C44)/2?]1/2 4.28
Vt|| Vt||=VT=(C44/?)1/2 2.64 Vt? [(C11-C12)/2?]1/2 1.83
[111] Vl' [(C11+2C12+4C44)/3?]1/2 4.41
Vt' [(C11-C12+C44)/3?]1/2 2.13

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 on the link.

Hosted by www.Geocities.ws