Problem 1) Molecular Potentials and Crystal Binding

Problem Set 9 Solutions

Spring 2003

PHYS 7440

Problem 1) Molecular potentials and Crystal binding

a) Use Mathematica to plot the potential above at ten different R positions and graphically find the minimizing values of Z. Plot the approximate values of Z(R) and produce a table of the minimum potential vs. proton-proton separation. This table is a numerical approximation to the actual effective proton-proton potential.

Often, the use of such numerical potentials is concidered not worth the effort for the precision desired. An alternative approach is to use a simple function that reproduces the numerical data reasonably well. The Lennard-Jones 6-12 potential is one example. Another is the Morse potential, which looks like this:

One of the cool things about this potential is that you can actually solve the Schroedinger Eq. exactly!

b) Use Mathematica to plot the Morse potential along with the points you found above for the numerical approximation to the actual effective proton-proton potential. Adjust and list the parameters of the potential to get close to the data. If you go nuts, you can use the NonlinearFit function of Mathematica to optimize the parameters. Now you have an approximate potential in closed form

Here’s the Mathematica stuff to do parts a) and b)

Fitting to a Morse potential

In the section above, we looked for the equilbrium or minimum value of the energy. Here, let's try to see how the total ground state energy is predicted to behave for various proton-proton separations. Eventually, this type of information allows us to estimate molecular and solid vibrational energies.

The problem is that we actually have an energy function of two parameters, Z and R, whereas we want a potential that is just a function of position. In the variational calculation, Z is considered a parameter that must be optimized at each position. Therefore, it is explicitely a function of position. Here, I simply sample the ground state energy function at a discrete set of points around the known equilibrium position of R=2 and Z=1.23 and solve for Z(R) graphically. I'll choose say ten points from R=0.3 to R=3, evenly spaced by 0.2 Bohr radii and maybe a few at larger separation:

In[1]:=

s[z0_,R_]:= Exp[-z0*R]*(1+(z0*R)+(((z0*R)^2)/3))

In[2]:=

int2[z0_,R_]:=s[z0,R]-(2*Exp[-z0*R]*(1+(z0*R)))

In[3]:=

int4[z0_,R_]:=Exp[-2*z0*R]*(1+(1/(z0*R))) -

(1/(z0*R))

In[4]:=

int5[z0_,R_]:=-Exp[-z0*R]*(1+(z0*R))

In[5]:=

egs[z0_,R_]:=(1/(1+s[z0,R]))*

( (z0^2)*(1-int2[z0,R]) -

2*z0*(1-int4[z0,R]-2*int5[z0,R]))

In[6]:=

Plot[(egs[z,r]+(2/r))/.r->0.3,

{z,1.5,2},

Frame->True]

Out[6]=

-Graphics-

Graphically, we find that the potential can be listed as a set of (energy, Z, R) with this particular point given by (2.938,1.88,0.3). Similarly, we can generate these graphs and find e.g.:

In[7]:=

Plot[(egs[z,r]+(2/r))/.r->4.9,

{z,0.8,1.3},

Frame->True]

Out[7]=

-Graphics-

Together, the results are:

(0.535,1.78,0.5), (-0.352,1.68,0.7), (-0.756,1.56,0.9), (-0.968,1.50,1.1), (-1.076,1.43,1.3), (-1.13,1.36,1.5), (-1.16,1.305,1.7), (-1.172,1.26,1.9), (-1.172,1.22,2.1), (-1.164,1.18,2.3), (-1.156,1.15,2.5), (-1.145,1.13,2.7), (-1.132,1.12,2.9), (-1.08,1.03,3.9), (-1.04,1.1.01,4.9).

I group these results as a list of energy vs. position. The only additional point is that the Morse potential that we aim to fit is referenced to an energy zero when the particles are at infinite separation. In the hydrogen molecule problem, we have an energy of -1Rydberg at infinite separation due to our inclusion of the hydrogen atom ground state energy. Therefore, in my list, I have offset the energy by 1Rydberg (remember, we always have this freedom):

In[9]:=

energyList={ {0.3,2.938},

{0.5,0.535+1},

{0.7,-0.352+1},

{0.9,-0.756+1},

{1.1,-0.968+1},

{1.3,-1.076+1},

{1.5,-1.13+1},

{1.7,-1.16+1},

{1.9,-1.172+1},

{2.1,-1.172+1},

{2.3,-1.164+1},

{2.5,-1.156+1},

{2.7,-1.145+1},

{2.9,-1.132+1},

{3.9,-1.08+1},

{4.9,-1.04+1} };

In[10]:=

energyPlot=ListPlot[energyList,

PlotRange->{-0.2,0.5}]

Out[10]=

-Graphics-

Here, I've cut off the highest energy point to accent the potential well section.

Next, we can try to construct a Morse potential of the same basic shape. From the calculation, we know that the well minimum is at 2 Bohr radii and that the well depth is -0.173. Then, the potential only has one free parameter:

In[11]:=

morsePot[kappa_,r_]:=0.173*

(Exp[-2*kappa*(r-2)]-

2*Exp[-kappa*(r-2)])

Now let's plot it up for a few kappa values and see how it compares:

In[12]:=

morsePlot1=Plot[morsePot[kappa,r]/.kappa->1,

{r,0,5},

PlotRange->{-0.2,0.5}]

Out[12]=

-Graphics-

In[13]:=

Show[morsePlot1,energyPlot]

Out[13]=

-Graphics-

Seems to blow up too fast. Let's try a smaller kappa:

In[14]:=

morsePlot3=Plot[morsePot[kappa,r]/.kappa->0.8,

{r,0,5},

PlotRange->{-0.2,0.5}]

Out[14]=

-Graphics-

In[15]:=

Show[morsePlot3,energyPlot]

Out[15]=

-Graphics-

Pretty good in the area of the well. Good enough for me. Thus, the Morse potential for the ionized hydrogen molecule would probably give a fair estimate of the vibrational properties for the parameters V0=0.173, Lambda=2, and Kappa=0.8

The Morse potential sees a fair amount of use in this role. In the next few parts, let’s imagine that we use it to model the binding potentials in a simple cubic solid.

Assume that each atom only interacts with its closest neighbor atoms.

c) Calculate the equilibrium lattice spacing, a, for the system. What is the total binding energy of the crystal (assuming N atoms and zero kinetic energy)?

In a cubic crystal, there are 6 nearest neighbors, all at the same distance. Thus, the total potential is just the sum of six equal Morse terms:

You find the equilbrium spacing by minimizing this result with respect to R:

Thus, the Morse potential is nice in that one of the parameters directly determines the equilibrium spacing.

d) Do a second order expansion of the total potential about the equilibrium lattice position. Determine how the spring constant, C, which describes harmonic motion of the atoms depends upon the parameters of Morse potential. (Recall that C is the coefficient in Hooke's Law for a spring such that the force on a spring when stretched to position, x, is F=Cx i.e., the potential is Cx2/2.)

This part is just a matter of expanding for small deviations around the equilibrium point i.e., we let:

with the deviation assumed small. Then,

Therefore, the approximate spring constant (which arises from the 6 different spings is just:

e) Typically, the sound velocity is related to the spring constant, lattice spacing and the mass of the atoms by:

Assume that the interacting atoms are Sodium atoms: Then, the mass is MNa=3.7x10-23gm. Futher, Sodium has a sound velocity of roughly 100,000cm/sec. and a lattice spacing of 3.5x10-8cm.

Use these numbers to produce an estimate of the value of C.

Plug and chug to get

Problem 2) Singularities in density of states.

This problem is designed to investigate the singularities that can develop in the density of phonon modes due to the nature of the phonon dispersion relationship. In fact, for essentially all problems in solid state, whether for phonons or electrons, densities of states are guaranteed to have singularities somewhere in the first Brillouin Zone. These singularities are generally referred to as Van Hove singularities.

To see why these occur, start with the general result for the phonon density of modes. For a 3-dimensional system, the mode density per unit volume is given by a result that is very similar to the electronic case:

I apologize to anyone who got confused, but in Kittel’s book (from which I selected these problems) he defines a density of modes which is not per unit volume:

In any case, to see why the density of states (either g or D) could have singularities, recall that the function is only uniquely given in the first Brillouin Zone. For k-vector values outside the zone, we just periodically reproduce the first zone values. Thus, you can think of as a periodic function of k-vector with a periodicity of the first Brillouin Zone length. Further, the values of  are bounded i.e., there is some minimum and some maximum value of frequency. Together, for reasonable functions, these two behaviors will lead to the function having points (usually at the maximum and minimum values) where is zero. Therefore, at these points, the integrand will diverge and the density of states will develop a range of singularity types depending on the dimensionality of the system. For propogation of plane-wave, the gradient of the dispersion relation is also referred to as the group velocity.

a) 1-Dimensional chain.

Here, we consider a 1D problem, namely the monatomic chain. The dispersion relation is:

First note that this dispersion relation is a maximum at the zone boundary where . Then,

so that:

Given this dispersion relation, we can generate the group velocity:

The next step in determining the density of modes is to evaluate the integral. In class, I considered the 3-dimensional case where the integral is over a 2-dimensional surface. By following the argument in class by which we derived the density of states, you can see that for a 2D system, the integral is over a 1D surface i.e., a line or several line segments. Finally, for a 1D system, the integral is over a 0D surface i.e., is just the evaluation of the integrand at one or more points. The integrand of interest is generally:

where d is the dimensionality of the system. The 'surface' of constant frequency for the 1D problem is actually the two points, +/-k, which yield the same frequency from the dispersion relation. Therefore, the integration is just the integrand at these two points, added together:

Finally, we must rewrite this result in terms of the frequency. We do this via the original dispersion relation:

or

Then, we find the desired result:

This result shows that the density of modes diverges at the zone boundary due to the fact that the group velocity goes to zero. Thus, there are asymptotically an infinite number of modes per unit frequency out at the zone edge.

Plotting of this function is straight-forward: The natural units for the frequency scale are clearly in terms of the maximum frequency. Rewriting the result above in terms of /m yields:

Then, it is reasonable to measure D in units of . Note that the units are then (number/frequency) as expected. Here are the Mathematica lines:

density[omega_]:=1/Sqrt[1-omega^2]

Now plot this up:

Plot[density[omega],{omega,0,0.99},

PlotRange->{0,10},

AxesOrigin->{0,0},

Frame -> True,

FrameLabel -> {"w/wm","D(w)/(2N/Pi)"}]

b) 3D Optical branch density of modes.

The first part of this problem is to see what the approximate optical phonon dispersion relation is like for the dimerized 1-D chain in the case where one spring is very much stronger than another i.e., K>G . Let’s start by rewriting the dispersion relation:

To first order in the small quantity G/K, we can expand the square root to find:

or

Next, consider the small ka expansion of this relation:

Therefore, the dispersion relation can be approximated as:

where:

Density of States

For this problem, we return to the full 3D case:

where for the 3D optical branch, we approximate the dispersion relation near k=0 by an inverted parabola:

In this approximation, it should be clear that there are no modes available for frequencies above 0. Therefore, in this region, the density of modes is zero.

However, below 0 there are modes and a certain density of them. In this region, the group velocity is given by:

Further, from the dispersion relation we can see that the surfaces of constant  correspond to the surfaces of constant k2 i.e., to spherical surfaces in k-space. Under these circumstances, the integral for the density of states becomes:

Using the dispersion relation to rewrite this result as a function of  gives us:

Thus, we have:

Here is the Mathematica plot using dimensionless units such that frequency is measured in units of 0 and density is in units of :

optic[omega_]:=If[ omega<1,Sqrt[1-omega],0]

Plot[optic[omega],{omega,0.5,1.5},

PlotRange->{0,1},

AxesOrigin->{0.5,0},

Frame -> True,

FrameLabel -> {"w/w0","D(w)/D(0)"}]

-Graphics-

Despite the statement in Kittel, a plot of this density of modes shows that it is not discontinuous. Rather, its first derivative is discontinuous. In general, it has been found that in 1D problems it is possible to have densities of modes which diverge and are discontinuous at certain points in the Brillouin Zone. However, in 2 or 3D, the integration changes divergences in the integrand into discontiunities in the slope of the density of modes.

Problem 3) Heat capacity of layered lattices.

a)Rigid layer coupling.

We consider a solid composed of layers that are rigidly coupled so that the only possible motion is in the plane of the layers. Then, each layer is effectively a 2D system described by lattice vibrations only in the plane and k-vector solutions that live in a 2D first Brillouin Zone. Therefore, this problem is an example that requires the calculation of the density of modes to be done as an integration of the inverse group velocity over a 'surface' of constant  that is actually just a line.

The generalization of the monatomic 1-D chain to a 2-dimensional square sheet yields the dispersion relation:

First, let's see what this dispersion relation looks like. I am choosing to define the function such that frequency is measured in units of , the maximum frequency (found at the zone corner) and each of the kvector components is in units of the inverse of the lattice spacing. Then, the dispersion relation is just:

omega[kx_,ky_]:=Sqrt[(2-Cos[kx]-Cos[ky])]/2

We can then plot it up:

Plot3D[omega[kx,ky],

{kx, -Pi, Pi},

{ky, -Pi, Pi}]

-SurfaceGraphics-

Thus, you can see that for small k-vectors near the zone center that the dispersion relation is approximately linear.

In the Debye model, we approximate the dispersion relation as a simple linear acoustic relation throughout the zone. Then,

so that the group velocity is just the sound velocity:

Then, the density of modes is given by a line integral of the inverse group velocity over the line of constant frequency. From the dispersion relation, it is clear that constant frequency implies constant k, so that the integration takes place around a circular path (similar to the spherical integral in 3D).

Then, the density of modes in this Debye approximation is given by:

This result is linear in the frequency as compared to the 3D case that goes as the second power of frequency.

To proceed with the Debye description, we calculate the approximate internal energy and specific heat by approximating the integrals over the first Brillouin Zone with an integral over a circular disk of k-space area sufficient to contain the same number of modes for each phonon branch as the original Brillouin Zone. This defines a Debye k-vector such that each branch of the phonon spectrum (each has N modes) is approximated by integrating the acoustic approximation out to:

or

With this definition for the Debye k-vector, we proceed to define the Debye frequency (the maximum limit of the frequency integration) as:

Then, the internal energy is given by (remember, there are two acoustic modes in 2D):

and the heat capacity is:

Plugging in the result above for the density of modes gives the following for the specific heat (heat capacity per unit area):

Change variables to so that we have:

Again, substitute into the denominator the value of kD and make the usual definition of the Debye temperature and we can rewrite this result as:

At very low temperatures, the integral can be approximated by letting the upper limit go to infinity so that the integral is just some definite number. In this case, the integral can be evaluated in terms of the Riemann Zeta Function as 6(3) or approximately 7.2. In any case, the specific heat goes as T2. The low temperature result is then:

Note also that at high T we can again approximate the integral by expanding the integrand for small x so that we get the correct 2D Dulong-Petite value:

b) Weak layer coupling.

Now think about what happens in a solid where the inter-layer coupling is very weak i.e., where the springs which connect the atoms between different layers are very soft so that the spring constant, C, between layers has a very small value. Recall that in our treatment of phonon dispersion relations that the spring constant and the mass of the atoms combine as something like to determine the maximum value of the frequency in the disperson relation. For example, in the simple 1D monatomic case, the dispersion relation is:

In the Debye model, we approximate the dispersion relation by its acoustic behavior. Therefore, in this limit,

and the very weak spring constant corresponds to very low sound velocity. Because the springs between the layers are so weak, it takes a very long time for energy on one layer to propagate to the next and so-forth, so the sound waves travel very slowly from layer to layer. In the limit of no coupling, the sound waves never leave the layer and we have a set of independent 2D layers.

Now, in terms of the specific heat behavior of such a system, imagine what the phonon branches look like in the zone: Two of them (the acoustic branches within the layers) have some comparable and rather large (compared with the weak-spring direction) sound velocity. The third branch has low sound velocity and therefore very shallow slope compared to the other two. The branches then look like this in the Debye approximation:

Given this picture, we can immediately provide an approximate description of how the specific heat is likely to behave at different temperatures:

High T:

At very high temperatures, the thermal energy is sufficient to excite all the vibrational modes of the lattice and we expect to see the high temperature Dulong-Petite value of:

Intermidiate T:

Each of the branches can be characterized by a Debye temperature below which the thermal energy is no longer sufficient to excite the higher frequency modes. This Debye temperature is directly related to the sound velocity since:

Therefore, the two in-layer branches will begin to enter the quantum regime at a much higher temperature than will the layer-to-layer branch. Because only two of the three branches are in the quantum limit, the specific heat begins dropping with a T2 dependence, essentially identical to that found above in the rigid layer case. The third mode simply continues to contribute its temperature independent Dulong-Petite value.