Vibrations of NaCl:  A Simple Model for Lattice Dynamics

Using a simplified model of balls and springs for vibrations in a NaCl type crystal results in a non-trivial problem; however, Mathematica allows for a straightforward solution.

by E. G. Emberly, A. Kung, and P. D. Loly

The motions of atoms in a crystal play an important role in determining many physical properties of the system.  For instance the modes of vibration are important in determining such properties as sound and heat propagation within the solid.  The earliest work on the problem of vibrations of a crystal lattice date back to the early part of this century where work was done on the subject by such people as Debye and Einstein.  It was not until the 1930's where the foundation of a general theory of lattice dynamics was put forth by Born and Huang [Born and Huang 1954].  The theory is well understood and its results have been found to agree with experimental results.

In the general theory, the problem is analyzed by expanding the potential energy of the crystal in a Taylor series.  The series is truncated using the "harmonic approximation", which is valid for most crystal systems since the displacements of the ions are usually small.  By using Newton's second law this then leads to an eigenvalue problem involving Born's "dynamical matrix".  Standard treatments of the above theory can be found in most advanced solid state physics textbooks, such as Ziman [Ziman 1964].  Unfortunately, the difficulty of such an approach is evaluating the necessary sums that arise in the definition of the crystal potential energy.  For ionic crystals this involves an infinite sum over Coulomb interactions between a central cell and the rest of the crystal.  Methods have been given to evaluate such sums [Kellerman 1940].  The method taken for this article will be a simpler approach.

When the problem of vibrations in a crystal is first encountered in an undergraduate physics course, a model where the lattice is treated as a collection of balls and springs is usually used [Kittel 1986].  This model is usually only applied to a linear monatomic or diatomic chain, and at most a monatomic 2-D lattice.  Nonetheless this model is useful in solving for the vibrational modes of some complicated structures.  In a previous article, Powell and Crandall used this sort of model to solve for the dispersion curves of a DNA molecule [Powell and Crandall 1993].  The next extension is to 3 D.  This is what we undertake in this article by considering perhaps the simplest 3 D case of a NaCl type crystal which gives both optical and acoustical modes.  The algebraic complexity of even this simple problem warrants a "symbolic" setup.  The method, as suggested by DeLaunay [DeLaunay 1956], considers the interactions between ions as springs of different strength.  This method allows for an analytical solution to the Born dynamical matrix from which can be extracted the dispersion curves for the crystal.  It is still a non-trivial problem and Mathematica allows for a fairly clear and efficient implementation of the solution.

A Simple Model for Vibrations of a NaCl Type Crystal:

The problem being solved in this article is for a crystal in the NaCl structure.  The space lattice is face-centred cubic (fcc) and has a basis consisting of the two ionic species separated by a distance, a.  The reciprocal lattice (which will arise later when the vibrational modes are examined), is a body-centred cubic lattice in reciprocal or k-space. As mentioned in the introduction, the model being used treats the interactions between ions as springs.  The solution is based on taking interactions out to the second nearest neighbours for the two ions in the basis.  This allows one to write out the force acting on each ion, and then as mentioned in the general theory, Newton's second law is used set up the equations of motion for the (x,y,z) coordinates of each ion.  The displacements of the ions are then assumed to have wave-like solutions with undetermined amplitudes.  This assumption leads to an eigenvalue problem for the frequency of oscillation, ω as a function of wave-vector k.  Since in this problem the number of ions per unit cell is 2, the solution will require finding eigenvalues for a 6 x 6 matrix.  This will result in six, possibly degenerate, values for ω.  Three will be found to be proportional to k at low k, and these are refered to as the "acoustic branches".  The other three are designated as "optical branches".

Define the Atomic Positions:

The preceding paragraph described a crystal in the NaCl structure as a fcc space lattice with a diatomic basis attached to it.  A different view can also be taken.  Crystals in the NaCl structure can be viewed as a superposition of two fcc sublattices one for the Na+ and another for the Cl-, combined so that a simple cubic lattice with alternating ions at each lattice postion results.  In this geometry there are six nearest neighbours (NN) at a distance, a, and there are 12 next nearest neighbours (NNN) at a distance of (Sqrt[2]) a from any ion.  To set up the geometry of the problem we need to define nearest neighbour (nn) and next nearest neighbour (nnn) positions for a crystal in the NaCl structure.

nn = {{a,0,0},{0,a,0},{-a,0,0},{0,-a,0},
	{0,0,-a},{0,0,a}};

nnn = {{a,a,0},{-a,a,0},{-a,-a,0},{a,-a,0},
	{0,a,-a},{0,a,a},{0,-a,a},{0,-a,-a},
	{a,0,-a},{-a,0,-a},{-a,0,a},{a,0,a}};

The solution to this problem is based on writing down the equations of motion for the two basis elements (namely Na+ and Cl-) in the central cell.  In our solution we are only considering interactions between one of the elements out to its next nearest neighbours.  The diagrams below depict the two types of ions and the cluster of ions with which they interact.  The diagrams are 2 dimensional cross sections of the 3 dimensional clusters being used.  In each diagram, the smaller white circles represent Na+ type ions, and the larger black disks represent the Cl- type ions.

Show[Graphics[{Table[Disk[{i, 0}, 0.25], {i, -1, 1, 2}], 
     Table[Disk[{0, i}, 0.25], {i, -1, 1, 2}], 
     Table[Circle[{i, j}, 0.2], {i, -1, 1, 2}, 
      {j, -1, 1, 2}], Circle[{0, 0}, 0.2]}], Axes -> True, 
   AxesLabel -> "Configuration 1:\n   Na Center", 
   Ticks -> None, AspectRatio -> 1, PlotRange -> All];

Show[
	Graphics[{
		Table[Circle[{i,0},0.2],
			{i,-1,1,2}],
		Table[Circle[{0,i},0.2],
			{i,-1,1,2}],
		Table[Disk[{i,j},0.25],
			{i,-1,1,2},{j,-1,1,2}],
		Disk[{0,0},0.25]
			}
		],
	Axes->True,
	AxesLabel->"Configuration 2:\n 
		Cl Center",
	Ticks->None,
	AspectRatio->1,
	PlotRange->All
	]

Set Up  the Equations of Motion:

The solution to the problem as mentioned in the introduction is based on a ball and spring model.  In essence, the problem is solved by setting up the equations of motion for the two atoms in the basis.  Since the problem is three dimensional, it will be necessary to take into account the force acting in the (x, y, z) directions for each basis atom.  Hence since the force between a basis ion and a neighbouring ion acts along the direction between the two ions, it is necessary to know the directions for all ions.  The following defines the unit vector pointing from the origin to the nth site in our cluster of ions (which will be denoted in the function as cluster[[n]] ).

e[cluster_,n_]:=
	cluster[[n]]/
		Sqrt[ Sum[ cluster[[n, j]]^2, 
			{j, 3}] ]//PowerExpand

For  example, the unit vector from the origin to the third nearest neighbour at (-a,0,0) is found using,

e[nn, 3]

The next step is to define the time dependent positions for each of the Na+ and Cl- type ions in the cluster.  The notation that will be adopted tries to alleviate some of the complexity that might arise if different symbols were used to set up the equations for the two different ions in our basis.  In the first configuration, which has the Na+ ion as the central ion, the Cl- type ions are nearest neighbours.  However, in the second configuration, where the Cl- ion is now at the center, the other Cl- ions are next nearest neighbours.  The same is true for the Na+ type ions.  To account for this we will define two sets of Na+ positions, s11 and s12.  The first index refers to the type of ion being considered ("1" refers to the Na+ ion being the first element of the basis).  The second index refers to the sets position relative to the central ion  ("1" being NN and "2" being NNN).  Using this notation, the problem is easily extendable to include third or fourth neighbour interactions by simply defining s13 and s23 and so on.   The sets themselves are composed of the position vectors for the ions they represent.  The position of the ith Na+ ion will be denoted by a vector u1[i] if it is a nearest neighbour and u2[i] if it is a next nearest neighbour.  The position of the ith Cl ion will be denoted in a similiar manner by a vector v1[i] or v2[i], again depending on whether it is a nearest or next nearest neighbour.  For example the x component of the ith nearest neighbour Na+ ion is referred to by using u1[i, 1].

s11 = Table[ u1[i, j],
		{i, 6}, {j, 3}];
		
s12 = Table[ u2[i, j],
		{i, 12}, {j, 3}];

The same is done for the Cl- ions, except now the first index will be ``2'' referring to the fact that Cl- is defined as the 2nd element of the basis.

s21 = Table[ v1[i, j],
		{i, 6}, {j, 3}];
		
s22 = Table[ v2[i, j],
		{i, 12}, {j, 3}];

As mentioned above, the forces between atoms will be considered to obey Hooke's Law.  The total force on one of the basis atoms will be the sum of the individual forces due the surrounding 18 ions.  The coupling constants are different for different pairs of ionic species.  The coupling between Na+ - Cl- ions is governed by g1, the coupling between Na+ - Na+ ions is governed by g2 and the coupling between Cl-  - Cl- is governed by g3.  The force due to the nth particle acting on the origin takes on the form F = -k x and is given by:
	F(n) = g [e(n).( (Subscript[s, 0]) - (Subscript[s, n]))] e(n)
where e(n) is the previously defined unit vector from the origin to the nth particle, and (Subscript[s, 0]) is the time dependent position of the particle at the origin.

Now define the time dependent locations of the Na+ ion at the origin in the first configuration, s10, and the Cl- ion considered at the origin in the second configuration, s20.

s10 = Table[ u0[j], {j, 3}];
s20 = Table[ v0[j], {j, 3}];

It is now possible to define the interaction force between the nth particle acting on the ion located at the origin.  The "1" denotes the case when the Na+ ion is at the origin and the "2" denotes the case when the Cl- ion is at the origin.  The arguments to the function are to be used in the following way.  The argument, s, can be any of s11, s12, s21 or s22.  The argument, n, refers to the nth ion in the specified set, s.  Based on the choice of s, the cluster will be either nn or nnn, depending on whether the chosen s refers to a nearest or next nearest neighbour cluster.  Of course, g is chosen depending on which interaction is involved.

fn10[cluster_, s_, n_, g_] := 
	- g ( e[cluster,n].( s10 - s[[n]] ) ) e[cluster,n]
	
fn20[cluster_, s_, n_, g_] := 
	- g ( e[cluster,n].( s20 - s[[n]] ) ) e[cluster,n]

For example the force between the the Cl- ion at (-a, 0, 0) and the central Na+ ion is given by

fn10[ nn, s21, 3, g1]

The total force acting on the Na+ ion located at the origin is the sum of all of its NN interactions (being Cl- ions)  plus all of its NNN interactions (being Na+ ions).  By using the command Collect,terms in the equation involving the coupling constants g1 or g2 are being collected so that the result for the force in one particular direction will be of the form (terms)*g1 + (terms)*g2.

f10 = Collect[
	Sum[ fn10[ nn, s21, i, g1], {i, 6}] + 
		Sum[ fn10[ nnn, s12, i, g2], {i, 12} ]//Simplify,
			  {g1, g2}];

The total force acting on a Cl- ion located at the origin is the sum of its NN interactions (being Na+ ions) and its NNN interactions (being Cl- ions, which is governed by g3).

f20 = Collect[
	Sum[ fn20[ nn, s11, i, g1], {i, 6} ] +
		Sum[ fn20[ nnn, s22, i, g3], {i, 12} ]//Simplify,
			  {g1, g3}];

Now the assumption is made that the solutions for the displacements are wave like (i.e. Exp[ i k r ]).  The amplitudes A are for the Na+ ions and the Amplitudes B are for the Cl- ions.

First define the wave vector.

k = {kx, ky, kz};

Now assign wave-like solutions to the already defined diplacements.  Here are the Na+ nearest neighbour displacements defined for the nth ion:

u1[n_, j_]:= A[j] Exp[I k.nn[[n]] ];

Cl- nearest neighbour displacements:

v1[n_, j_]:= B[j] Exp[I k.nn[[n]] ];

Na+ next nearest neighbour displacements:

u2[n_, j_]:= A[j] Exp[I k.nnn[[n]] ];

Cl- next nearest neighbour displacements:

v2[n_, j_]:= B[j] Exp[I k.nnn[[n]] ];

Define the position of the Na+ ion at origin in the first configuration. The Cl- ion is at the orgin of the second configuration.

r10 = {0,0,0};
r20 = {0,0,0};

Assume displacements of these ions are also wave-like.

u0[j_]:= A[j] Exp[I k.r10 ];

v0[j_]:= B[j] Exp[I k.r20 ];

Now use Newton's 2nd law F = ma where m = M1 is the mass of the Na+ type ion, and m = M2 is the mass of the Cl- type ion.  The acceleration, a, for a solution which is harmonic in time is equal to -(w^2) (Subscript[s, 0]) where w is the frequency of oscillation.  The following equations are being collected in terms of the unknown amplitudes (A[1], A[2], A[3], B[1], B[2], B[3]).

eqnsNa = Table[
		-w^2 M1 s10[[i]] == 
			Collect[ f10[[i]],
				{A[1],A[2],A[3],B[1],B[2],B[3]} ],
		{i, 3}];

eqnsCl = Table[
		-w^2 M2 s20[[i]] ==
			Collect[ f20[[i]],
				{A[1],A[2],A[3],B[1],B[2],B[3]} ],
		{i, 3}];

Form the Eigenvalue Problem:

The resulting equations involve sums of complex exponentials so it will be necessary to convert the complex exponentials into trig functions.  The following package needs to be loaded to accomplish this task.

<<Algebra`Trigonometry`

For purposes of being able to extract the required matrix, the right hand side of the above equations is brought over to the left hand side (-w^2 M - rhs ==0), and the coefficients of the aformentioned unknowns are collected again.

eqnsNa2 = Table[
	Collect[
		eqnsNa[[i, 1]] - 
			ComplexToTrig[ eqnsNa[[i, 2]] ]//Simplify,
		{A[1],A[2],A[3],B[1],B[2],B[3]}
			] 
			 	== 0,
		{i, 3}];

eqnsCl2 = Table[
	Collect[
		eqnsCl[[i,1]] - 
			ComplexToTrig[ eqnsCl[[i, 2]] ]//Simplify,
		{A[1],A[2],A[3],B[1],B[2],B[3]} 
			] 
				== 0,
		{i, 3}];

Here is the equation of motion for the x component of the central Na+ ion in terms of the 6 unknowns.  (Some substitutions are made to make the equation more legible).

eqnsNa2[[1]]/.Cos[a kx]->Cx/.Sin[kx a]->Sx/.
	Cos[ky a]->Cy /.Sin[ky a]->Sy /.
	Cos[kz a]->Cz /.Sin[kz a]->Sz//Simplify

These two sets of three equations form a complete set of six equations in the six unknowns, A[1], A[2], A[3], B[1], B[2], and B[3].  For there to be a solution the determinant of the coefficients must vanish.  This provides an eigenvalue problem in ω, and allows for the determination of the frequencies for a given k.  The following extracts the matrix from the above six equations.  First of all collect both sets into one big set of 6 equations.

eqns = Join[eqnsNa2,eqnsCl2];

The next set of code does the actual extraction of the coefficients for the six unknowns.

Alist = Table[
	Table[
		Coefficient[ eqns[[i,1]], A[j]],
			{i, Length[eqns]} ],
		{j, 3}];

Blist = Table[
	Table[
		Coefficient[ eqns[[i,1]], B[j]],
			{i, Length[eqns]} ],
		{j, 3}];

Now form the 6 x 6 matrix.  (The dynamical matrix is the following matrix without the -M (ω^2) along the diagonal).

M = Transpose[
		Join[Alist,Blist]
			];

Dispersion Curves:

With the matrix set up above, it is now possible to solve for the dispersion relations along different directions of the crystal lattice.  Two of these directions will be explored below, namely the [100] and [110] direction in k-space.  (Note: to simplify the computation, the determinants will be solve for (ω^2) = W, hence (ω^2) is replaced by W, and then afterwards one can simply take the square root to get the angular frequency).  It should be mentioned that the following procedure was also attempted for the [111] direction, but the computation ran for a long time on a Sun Sparc 20 and no results were obtained.

Along the [100] direction ky = kz = 0.

M100 = M /. ky->0 /. kz-> 0/.w^2->W//Simplify ;

Along the [110] direction, kx = ky; kz=0.

M110 = M /. kx->K /. ky->K /.kz->0/.w^2->W//Simplify;

The following values for the masses and the spring constants will be used.  The masses do correspond to the masses of Na+ and Cl-, but the spring constants have been chosen arbitrarily - although the g1 constant (i.e. Na+ - Cl-) interaction will be the strongest.  There is a vast amount of data achieved from scattering experiments on the lattice dynamics of various crystals, so finding reasonable values for the spring constants is possible.  To see specific experimental results for dispersion curves of real NaCl one can look at the paper by W. J. L. Buyers and T. Smith [1966].

m1 = 3.8 10^-26;
m2 = 5.9 10^-26;

G1 = 6;
G2 = 3;
G3 = 2;

Now solve for the 6 possible values of (ω^2) for the [100] direction.

sol100 = Solve[
	Det[M100]==0,
			W];

Substitute the solutions into the 6 squared frequencies:

W100 = Table[
			sol100[[i,1,2]],
				{i,1,6}];

The square roots of the above 6 equations will be plotted to produce the dispersion curves. We change variables to kx a -> π/2 X, where X is the reduced wave vector.  The extra factor of 2 in the above change occurs because of our use of a to specify the spacing between Na+ and Cl-, whereas the conventional definition of a as the lattice parameter is twice this.   In the [100] direction the zone boundary for a bcc (the bcc lattice is reciprocal to a crystal with a fcc lattice) lattice ends at k = 2π/a, so X ranges from 0 to 2.  A similiar approach was used by Keeler [Johnston, Keeler, Rollins and Spicklemire 1996] in his Turbo Pascal program to generate dispersion curves for NaCl.

W100 =
	W100/.kx a->(Pi/2)X;

Here's a plot of the frequencies versus the reduced wave vector in the [100] direction.  (Note the branching into acoustic and optical modes and also the existance of degeneracy in the travnsverse vibrational modes).

Plot[
	Sqrt[
W100/.M1->m1/.M2->m2/.g1->G1/.g2->G2/.g3->G3]//Evaluate,
		{X,0,2}];

Try [110] direction:

sol110 = Solve[
	Det[M110]==0,
			W];

W110 = Table[
			sol110[[i,1,2]],
				{i,1,6}];

Again make a change of variables.

W110 =
	W110/.K a->(Pi/2)X;

For the bcc reciprocal lattice, the zone boundary in the [110] direction ends at ( (3/2)π/a, 
(3/2)π/a ), so X goes from 0 to 3/2.

Plot[
	Sqrt[
W110/.M1->m1/.M2->m2/.g1->G1/.g2->G2/.g3->G3]//Evaluate,
	{X,0,3/2}];

Thus in the [110] direction all of the degeneracies are broken so that the transverse modes are separate along with the longitudinal vibrational mode, and again optical and acoustic splitting is found.

Conclusions:

Using a simplified model of balls and springs for vibrations in a NaCl type crystal has resulted in a non-trivial problem.  However using Mathematica allowed for a straightforward solution to the problem without having to do the necessary algebra by hand.  Although it was not displayed in the paper, exact analytical results for the dynamical matrix and the dispersion curves can be extracted using the above code.  The dispersion curves were found and it was interesting to find that such a simple model would yield rather intricate results.  This method can be easily extended to include third neighbour interactions, however this adds to the complexity of the problem and makes the analytical results less tractible.  It is possible to study different crystal geometries also with only minor changes to the code.

References:

Born, M. and K'un Huang, Dynamical Theory of Crystal Lattices, Oxford Clarendon Press, 1954.

Buyers, W. J. L, and Smith, T., Phys. Rev., 150, 758 (1966).

DeLaunay, J., Solid State Physics:  Advances in Research and Applications, (ed. Seitz, F., Turnbull, D., and Ehrenreich, H.), 2, 219,    New York Academic Press, 1956.

Johnston, Keeler, Rollins, and Spiklemire, Solid State Physics Simulations:  Consortium for Upper-level Physics Software, John Wiley & Sons, 1996.

Kellerman, E. W., Phil. Trans. Roy. Soc., A238, 513 (1940).

Kittel, C., Introduction to Solid State Physics, 6th edition, Wiley & Sons, 1986.

Powell, J. W. and Crandall, R. E., Mathematica in Education, Spring 1993, pp. 15-20.

Ziman, J. M., Principles of the Theory of Solids, Cambridge University Press, 1964.

About the Authors

E. G. Emberly, A. Kung, and P. D. Loly
		
Department of Physics,
University of Manitoba,
Winnipeg, MB
R3T 2N2