    from matplotlib import pyplot
from matplotlib.patches import Circle
import numpy as np
from scipy.integrate import simpson

from common import draw_classic_axes, configure_plotting
import plotly.graph_objs as go

configure_plotting()

pi = np.pi

Interatomic interaction and molecular spectra¶

(based on chapters 6.1, 6.3, 8 of the book)

Expected prior knowledge

Before the start of this lecture, you should be able to:

Apply a Taylor expansion
Write down Newton's equations of motion of masses on springs
Write down LCAO Hamiltonian matrices describing molecular orbitals
Define and diagonalize matrices numerically (for the exercises)

Learning goals

After this lecture you will be able to:

Explain the origins of interatomic forces
Compute vibrational spectra of small molecules in 1D
Formulate Hamiltonians and equations of motion of bulk materials

In this lecture, we will study the vibrational spectra of small molecules. We start by discussing why it is reasonable to analyze the vibrations of molecules by modeling their atomic constuents as masses connected by springs. This modeling will enable us to formulate Newtonian equations of motion describing the molecular vibrations. These matrix equations form an eigenvalue problem that we will solve to find the eigenfrequencies and eigenfunctions of the molecular motion. The matrices we will find are very similar to those found in the LCAO model to describe the electronic spectrum of molecules. We will highlight the similarities and how to solve the equations.

Adding repulsion to the interatomic interaction¶

In the previous lecture, we observed that the diatomic molecule's covalent bonding energy decreases with decreasing interatomic distance. However, the trend is unphysical; the two atoms must start repelling eventually (at least when the nuclei get close, but really already much earlier). Hence from now on, we consider an interatomic repulsion between atoms that get sufficiently close. Including this repsulsion, we expect the energy of the molecule to depend on the inter-atomic distance $r$ as sketched in the figure

fig, ax = pyplot.subplots()
r = np.linspace(0.2, 3, 100)
Q = 0.01/r**3
U = np.exp(-r)
ax.plot(r, Q + U, label="$E_-$")
ax.plot(r, Q - U, label="$E_+$")
ax.set_xlabel("$r$")
ax.set_ylabel("$E$")
ax.set_xlim(-.1, 3.2)
ax.set_xticks([])
ax.set_yticks([])
ax.legend()
draw_classic_axes(ax)

svg

Towards modeling molecules as atoms connected by springs¶

Suppose we have an interatomic potential described by the blue curve in the figure below. The atoms are in equilibrium at the bottom of the interatomic potential. Let the bottom of the potential be $U = U_0$ at the interatomic distance $r = a$.

fig, ax = pyplot.subplots()
r, dr = np.linspace(0.21, 2, 100, retstep=True)
Q = 0.05 / r**2
U = np.exp(-r)
E = Q - U
x_0 = np.argmin(E)
a = r[x_0]
kappa = np.diff(E, 2)[x_0-1] / dr**2 / 2
ax.plot(r, E, label="Interatomic potential")
ax.plot(r, E[x_0] + kappa * (r - a)**2, linestyle="--", label="Parabolic approximation")
ax.set_xlabel("$r$")
ax.set_ylabel("$U$")
ax.set_xlim(-.1, 2.2)
ax.set_ylim(-.5, .3)
ax.set_xticks([a])
ax.set_xticklabels(['$a$'])
ax.legend()
ax.set_yticks([E[x_0]])
ax.set_yticklabels(['$U_0$'])
draw_classic_axes(ax, xlabeloffset=0.04)

svg

Near the minimum, the potential is approximately parabolic. We find the parabolic dependence by Taylor expanding the interatomic potential around the minimum:

\[ U = U_0 + \frac{\kappa}{2!} (r - a)^2 - \frac{\kappa_3}{3!} (r - a)^3 + ... \]

where the first derivative vanishes because we consider a minimum. Up to second order, the potential is quadratic and gives rise to a harmonic equation of motion. The higher order term $\kappa_3$ provides anharmonicity. We see that when the deviation from the equilibrium position is small, the atoms can be described by masses connected by a spring with spring constant $\kappa$.

Relation between the inter-atomic potential and rigidity¶

By measuring how much a material compresses or stretches under the application of stress, we can infer the harmonic force constant $\kappa$. To that end, assume we have a 1D material with length $L$. If the equilibrium distance between atoms is $a$, then there are $N = L/a$ atoms in the material. Now lets slightly stretch the material to length $L + \delta L$ by applying a force $F$. Under such a force, the interatomic distance increases by $\delta r = \delta L/N$. We don't stretch too hard and therefore the change in the interatomic distance $\delta r $ remains small. Therefore, we can approximate the interatomic potential around the minimum by the quadratic term (Taylor expansion up to second-order). Given the quadratic interatomic potential, the material responds with a returning force:

\[ F = - \frac{\mathrm{d} U}{\mathrm{d} r} \Bigr|_{r = a+\delta r} = -\kappa a \frac{\delta L}{L}. \]

Knowing the applied force and measuring the strain $\delta L /L$, we can thus extract a macroscopic material parameter: the compressibility:

$$ \beta \equiv -\frac{1}{L} \frac{\partial L}{\partial F}= \frac{1}{\kappa a}. $$ from which we can extract the spring constant $\kappa$. Furthermore, the compressibility allows us to calculate the speed of sound by using the relation:

\[ v = \sqrt{\frac{1}{\rho \beta}}, \]

where $\rho$ is the density of the material. Because $\rho = m/a$ (in 1D), we can express the speed of sound in terms of the inter-atomic spring constant:

\[ v = \sqrt{\frac{\kappa a^2}{m}}. \]

This way we can already see how the phonon heat capacity emerges from the microscopic material structure.

Thermal expansion¶

While the quadratic approximation of the interatomic potential is certainly the most important, the anharmonic term $\kappa_3(r-a)^3/6$ also has physical significance. Let us examine its role visually by comparing the harmonic and the third order approximations in the plot below. Notice that the second-order approximation is symmetric around the minimum while the third-order term is not.

r = np.linspace(-2, 2.5, 750)

a = 0
b = 0
c = 1
d = -0.2

U_quadratic = a + b * r + c * r**2
U_cubic = U_quadratic + d * r**3

r_min = -2
r_max = 2.5
U_min = -0.2
U_max = 4

E_t_min = 0.4
E_t_max = 3.5
N_values = 20
l_width = 1.5
N_active = 0

# Create figure
fig = go.Figure()


def U_c(r):
    return a + b * r + c * r**2 + d * r**3


def line(E_t):
    right = min(np.roots([c, b, a-E_t]))
    roots = np.roots([d, c, b, a-E_t])
    roots = np.real(roots[np.isreal(roots)])
    roots.sort()
    left = roots[1]
    return [left, right]

def avg_pos_cubic(E_t):
    Z = simpson(np.exp(- U_cubic / E_t), x=r)
    r_avg = simpson(r * np.exp(- U_cubic / E_t), x=r)
    x = r_avg / Z
    return x


# Add traces, one for each slider step
for E_t in np.linspace(E_t_min, E_t_max, N_values):
    avg = avg_pos_cubic(E_t)
    fig.add_trace(
        go.Scatter(
            visible = False,
            x = r,
            y = U_quadratic,
            mode = 'lines',
            line_color = 'blue',
            name = "Quadratic potential",
        ))
    fig.add_trace(
        go.Scatter(
            visible = False,
            x = r,
            y = U_cubic,
            mode = 'lines',
            line_color = 'red',
            name = "Cubic potential",
        ))
    fig.add_trace(
        go.Scatter(
            visible = False,
            x = line(E_t),
            y = [E_t, E_t],
            mode = 'lines',
            line_color = 'black',
            name =  r'Thermal energy level'
        ))
    fig.add_trace(
        go.Scatter(
            visible = False,
            x = [0, 0],
            y = [0, E_t],
            mode = 'lines',
            line_color = 'blue',
            line_dash = 'dot',
            name =  r'$\langle r \rangle$ for the quadratic potential'
        ))
    fig.add_trace(
        go.Scatter(
            visible = False,
            x = [avg, avg],
            y = [U_c(avg), E_t],
            mode = 'lines',
            line_color = 'red',
            line_dash = 'dot',
            name =  r'$\langle r \rangle$ for the cubic potential'
        ))

# Initial starting image
N_trace = int(len(fig.data)/N_values) # Number of traces added per step
for j in range(N_trace):
    fig.data[N_active*N_trace+j].visible = True

# Creation of the aditional images
steps = []
for i in range(int(len(fig.data)/N_trace)):
    step = dict(
        method = "restyle",
        args = [{"visible": [False] * len(fig.data)}],
        value = str(0.1*(i+1))
    )
    for j in range(N_trace):
        step["args"][0]["visible"][N_trace*i+j] = True  # Toggle i'th trace to "visible"
    steps.append(step)

# Creating the slider
sliders = [dict(
    tickcolor = 'White',
    font_color = 'White',
    currentvalue_font_color = 'Black',
    active = N_active,
    name = r'Thermal Energy',
    font_size = 16,
    currentvalue = {"prefix": r'Thermal Energy k_B T: '},
    pad = {"t": 50},
    steps = steps,
)]

# Updating the images for each step
fig.update_layout(
    sliders = sliders,
    showlegend = True,
    plot_bgcolor = 'rgb(254, 254, 254)',
    width = 700,
    height = 580,
    xaxis = dict(
        range=[r_min, r_max],
        visible = True,
        showticklabels = True,
        showline = True,
        linewidth = l_width,
        linecolor = 'black',
        gridcolor = 'white',
        tickfont = dict(size = 16)),
    yaxis = dict(
        range = [U_min, U_max],
        visible = True,
        showticklabels = True,
        showline = True,
        linewidth = l_width,
        linecolor = 'black',
        gridcolor = 'white',
        tickfont = dict(size = 16)),
    title = {'text': r'Thermal expansion of cubic potential',
        'y':0.9,
        'x':0.45,
        'xanchor': 'center',
        'yanchor': 'top'},
    xaxis_title = r'$r$',
    yaxis_title = r'$U [k_b T]$',
)

# Edit slider labels and adding text next to the horizontal bar indicating T_E
for i in range(N_values):
    fig['layout']['sliders'][0]['steps'][i]['label'] = ' %.1f ' % ((E_t_max-E_t_min)*i/(N_values-1) + E_t_min)

fig

The asymmetry due to nonzero $\kappa_3$ slows the growth of the potential when the interatomic distance increases. On the other hand, when the interatomic distance decreases, the asymmetry accelerates the growth of the potential. Therefore, stretching the material is more energetically favorable than contracting it. As a result, thermal excitations increase the interatomic distance. This gives us a simple model of thermal expansion.

Van der Waals bond¶

While we will focus on covalent bonding, let us also review another bond type.

A Van der Waals bond originates from an attraction between the dipole moments of two atoms. Suppose we have two atoms separated by an interatomic distance $r$ and that atom 1 has a dipole moment $\mathbf{p_1}$ pointing in the same direction as $\mathbf{r}$. Then, atom 1 creates an electric field

\[ \mathbf{E} = \frac{\mathbf{p_1}}{2\pi \varepsilon_0 r^3} \]

at the position of the atom 2. Atom 2 then develops a dipole moment $\mathbf{p_2} = \chi \mathbf{E}$ with $\chi$ the polarizability of the atom.

The potential energy between between the two dipoles is

\[ U(r) = -\mathbf{p_2}\cdot\mathbf{E}= -\chi |E|^2 \propto \frac{1}{r^6}. \]

The dipole attraction is much weaker than the covalent bonds but drops slower with increasing distance.

How does the strength of a covalent bond scale with distance?

The strength of the bond is determined by the interatomic hopping integral $-t = \langle 1 | H | 2 \rangle$. Since the wavefunction of a bound electron typically decays exponentially, so does the overlap integral.

Although the Van der Waals force is weak, it is the only force when there are no chemically active electrons or when the atoms are too far apart to form covalent bonds. Therefore, there are materials where Van der Waals interactions are the dominant interactions. An example of such a material is graphite. The Van der Waals bonds in graphite hold layers of covalently bonded carbon atoms together:

(image source: Wikipedia)

Looking ahead: multiple atoms¶

So far we have only considered the interatomic interactions between diatomic systems. However, our aim is to understand electrons and phonons in solids containing $N\to\infty$ atoms.

Let us see what happens when we consider more than two atoms.

Phonons¶

In order to understand phonons better, we need to understand how a vibrational motion in a solid arises. To that end, we model an array of atoms that are connected by springs with a spring constant $\kappa$. Our plan:

Consider only a harmonic potential acting between the atoms ($\kappa$ is constant)
Write down equations of motion
Compute normal modes

For simplicity we consider 1D motion, and let us start with a chain of 3 atoms:

# Defining constants
y_max = 6
max_m = 3

# Off set for the new masses
deviation_arr = [+2, -2, -3.5]

def plot_spring(x1, x2, y, annotation  = r'$\kappa$'):
    L = (x2-x1)

    width, nturns = 1, 10
    N = 1000
    pad = 200

    # Make the spring, unit scale (0 to 1)
    x_spring = np.linspace(0, L, N) # distance along the spring
    y_spring= np.zeros(N)
    y_spring[pad:-pad] = width * np.sin(2*np.pi * nturns * x_spring[pad:-pad]/L)

    x_plot, y_plot = np.vstack((x_spring, y_spring))

    # And offset it
    x_plot += x1
    y_plot += y
    ax.plot(x_plot, y_plot, c='k')
    ax.annotate(annotation, ((x2+x1)/2,y+2), fontsize = 40)


def plot_mass(x,y, annotation = r'$m$'):
    mass = Circle((x, y), .5, fc='k')
    ax.add_patch(mass)
    ax.annotate(annotation, (x-.5,y+2), fontsize = 30)


def make_plot(u, deviation_arr, max_masses = 3):
    # plotting initial masses
    plot_mass(0,0)
    plot_mass(0+deviation_arr[0],u, annotation = '')

    # Plot initial dotted lines
    pyplot.plot((0,0), (0,u-2.5), 'k--')
    pyplot.plot((deviation_arr[0], deviation_arr[0]), (u, u-2.5), 'k--')
#     pyplot.plot((0, deviation_arr[0]), (u-2.5, u-2.5), 'k--')
    pyplot.arrow(0, u-2.5, deviation_arr[0]-.5, 0, fc="k", ec="k", head_width=.25, head_length=.5)

    # Annotate the deviations
    annot_arr = [r'$u_1$', r'$u_2$', r'$u_3$']
    off_set = -1
    ax.annotate(annot_arr[0], (deviation_arr[0]/2+off_set+.4, u-3+off_set), fontsize = 30)

    # Annotate large arrow containing the interatomic distance
    pyplot.arrow( 0, u/3, 9, 0, fc="k", ec="k", head_width=.5, head_length=1)
    ax.annotate(r'$a$', (5, u/3+off_set), fontsize = 30)

    # Plotting other masses and springs
    for j in range(max_masses-1):
        # Equilibrium plot
        plot_spring(10*j,10*(j+1),0)
        plot_mass(10*(j+1), 0)

        # Plot containing deviations from it's equilibrium
        plot_spring(10*j+deviation_arr[j], 10*(j+1)+deviation_arr[j+1], u, annotation = '')
        plot_mass(10*(j+1)+deviation_arr[j+1], u, annotation = '')

        # Plot dotted lines
        pyplot.plot((10*(j+1), 10*(j+1)), (0,u-2.5), 'k--')
        pyplot.plot((10*(j+1)+deviation_arr[j+1], 10*(j+1)+deviation_arr[j+1]), (u, u-2.5), 'k--')
#         pyplot.plot((10*(j+1), 10*(j+1)+deviation_arr[j+1]), (u-2.5, u-2.5), 'k--')
        pyplot.arrow(10*(j+1), u-2.5, (deviation_arr[j+1]+.5), 0, fc="k", ec="k", head_width=.25, head_length=.5)

        # Annotations
        ax.annotate(annot_arr[j+1], (10*(j+1)+deviation_arr[j+1]/2+off_set+.4, u-3+off_set), fontsize = 30)

# Initializing figure
fig, ax = pyplot.subplots(figsize=(10,7))
pyplot.axis('off')
ax.set_xlim((-1,10*(max_m-1)+1))
ax.set_ylim((-y_max-4.5, y_max-2.5))

# Plottig system
make_plot(-y_max, deviation_arr, max_m);

svg

Let us denote the deviation of atom $i$ from its equilibrium position by $u_i$. Newton's equations of motion for this system are then given by

\[\begin{align*} m \ddot{u}_1 &= - \kappa (u_1 - u_2), \\ m \ddot{u}_2 &= - \kappa (u_2 - u_1) - \kappa (u_2 - u_3), \\ m \ddot{u}_3 &= - \kappa (u_3 - u_2). \end{align*}\]

We write this system of equations in matrix form

\[ m \ddot{\mathbf{u}} = -\kappa \begin{pmatrix} 1 & -1 & 0\\ -1 & 2 & -1\\ 0 & -1 & 1 \end{pmatrix}\mathbf{u} \]

We are interested in phonons, patterns of motion that are periodic and have a fixed frequency $\omega$. Hence we guess that the motion of the atoms is

\[ \mathbf{u}(t) = \mathbf{u}_0 e^{i\omega t}. \]

We substitute our guess into the equations of motion to yield an eigenvalue problem:

\[ \omega^2 \mathbf{u}_0 = \frac{\kappa}{m} \begin{pmatrix} 1 & -1 & 0\\ -1 & 2 & -1\\ 0 & -1 & 1 \end{pmatrix}\mathbf{u}_0. \]

The solutions to the eigenvalue problem are phonon modes.

Electrons¶

We just looked at how a chain of atoms moves. Let us now look at how the electrons of those atoms behave. To that end, we consider a 3 atom chain without any motion.

In order to understand how electrons behave, we use the LCAO model. The LCAO model generalizes in a very simple way. Let us consider the wavefunction:

\[ \vert \psi\rangle = \varphi_1 |1\rangle + \varphi_2 |2\rangle + \varphi_3 |3\rangle. \]

\[\begin{align*} E \varphi_1 &= E_0 \varphi_1 - t \varphi_2\\ E \varphi_2 &= E_0 \varphi_2 - t \varphi_1 - t \varphi_3\\ E \varphi_3 &= E_0 \varphi_3 -t \varphi_2. \end{align*}\]

Again, we write this in a matrix form:

\[ E \begin{pmatrix} \varphi_1 \\ \varphi_2 \\ \varphi_3 \end{pmatrix} = \begin{pmatrix} E_0 & -t & 0 \\ -t & E_0 & -t \\ 0 & -t & E_0 \end{pmatrix} \begin{pmatrix} \varphi_1 \\ \varphi_2 \\ \varphi_3. \end{pmatrix} \]

Numerical test¶

Diagonalizing large matrices is unwieldy, but let's try and check it numerically to see if we notice a trend. Let us first model 3 atoms on a chain. Diagonalizing the 'spring matrix' for 3 atoms connected by springs, we find the 3 eigenfrequencies $\sqrt{\tfrac{m}{\kappa}}\omega$

def DOS_finite_phonon_chain(n):
    rhs = 2 * np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1)
    rhs[0, 0] -= 1
    rhs[-1, -1] -= 1
    pyplot.figure()
    pyplot.hist(np.sqrt(np.abs(np.linalg.eigvalsh(rhs))), bins=30)
    pyplot.xlabel("$\omega$")
    pyplot.ylabel("Number of eigenfrequencies")

DOS_finite_phonon_chain(3)

svg

Similarly, diagonalizing the LCAO matrix describing a 3-atom chain with onsite energy $E_0=2$ and nearest-neighbour hopping $-t=-1$ yields the 3 eigenenergies of the 3 molecular orbitals (calculated analytically in exercise 2 of the LCAO lecture):

def DOS_finite_electron_chain(n):
    rhs = 2 * np.eye(n, k = 0) - np.eye(n, k = 1) - np.eye(n, k = -1)
    pyplot.figure()
    pyplot.hist(np.linalg.eigvalsh(rhs), bins=30)
    pyplot.xlabel("$E$")
    pyplot.ylabel("Number of eigenenergies")

DOS_finite_electron_chain(3)

svg

However, 3 atoms are far too few to model an actual solid. Hence, we need 'many more' atoms.

From 3 atoms to 300¶

We now analyze a 300-atom chain. Numerically diagonalizing the 300x300 spring matrix yields the 300 eigenfrequencies of the phonon modes. A histogram shows the distribution of the eigenfrequencies:

DOS_finite_phonon_chain(300)

svg

We observe that when $\omega$ is small, the density of states is constant. This is in line with what we saw in the Debye model, where the density of states for a 1D system was indeed a constant! However, when the frequencies are higher, the density of states deviates from the constant behaviour.

Performing a similar analysis for a 300-atom LCAO matrix, we compute a histogram with the distribution of the 300 electron eigenenergies:

DOS_finite_electron_chain(300)

svg

The numerical results once again agree with the model we developed earlier: In the Sommerfeld free-electron model, the density of states in 1D was proportional to $\frac{1}{\sqrt{E}}$. The above histogram reflects this proportionality at small energies $E$. However, when $E$ is higher, we observe a significant deviation from the $\frac{1}{\sqrt{E}}$ behavior.

For both the vibronic modes and the electronic orbitals, we find that our models agree with the numerical results whenever the frequencies/energies are small. However, when the frequencies/energies are high, we find that there is a significant deviation from the Debye/Sommerfeld models. The nature of this deviation is the subject of the next lecture!

Summary¶

The interatomic potential describes the rigidity and the thermal expansion of a material.
The approximately harmonic (quadratic) interatomic potential motivates describing atomic chains by masses connected by springs.
In the small $\omega$ limit, the density of states of the phonon modes in the Debye model is justified by modeling the atoms as masses connected by springs.
In the small $E$ limit, the density of states of electrons in the Sommerfeld model is justified by modeling the electrons as particles that can hop between atoms.

Exercises¶

Warm-up questions*¶

Discuss the interatomic potential sketched in this lecture. How do you determine which part of the interatomic potential is attractive and which is repulsive?
Write down the equations of motion describing two masses connected by a spring in matrix form. Argue (no calculations!) what should be the form of the eigenvectors. Compute the eigenfrequencies.
Write down the LCAO matrix for three atoms. Assume that each atom has a single orbital, onsite energy $E_0$ and the hopping between neighboring atoms $-t$. In addition, include a next-nearest neighbour hopping $-\tilde{t}$.
Write down the LCAO matrix for six atoms. You may assume that each atom has a single orbital, onsite energy $E_0$ and the hopping between neighboring atoms $-t$.

Exercise 1*: The vibrational modes of a linear triatomic molecule¶

In this exercise, we analyze the vibrational modes of a 3-atom molecule. We do so by modeling the molecule as 3 atoms connected by springs. We will compute the eigenfrequencies and eigenvectors analytically, making use of the symmetry of the system to simplify the analysis. This analysis is very similar to that used in exercise 2 of the LCAO lecture, where we analyzed the eigenenergies of electrons on a 3-atom molecule.

Consider carbon dioxide (C0$_2$), which is a linear, triatomic molecule shown below

carbon dioxide

source

By Jasek FH. - Own work, CC BY-SA 4.0, Link

Excluding the modes corresponding to a translation of the entire molecule, how many normal modes does this molecule have assuming motion in only 1D? And if the atoms can move in all three dimensions?
For simplicity, we only consider 1D motion of the atoms. Write down Newton's equations of motion for the atoms, you may assume that the spring constant is the same for both bonds.
Consider the symmetric vibrational mode, for which the displacements of the oxygen atoms are equal in magnitude and have an opposite direction. Write down the eigenvector (no calculations) of this mode and find its eigenfrequency. Discuss the eigenfrequency in relation to that for a single mass on a spring connected to a wall.
Now consider the antisymmetric mode, in which both oxygen atoms move in phase and have the same displacement. Considering that the center of mass of the molecule should remain at rest, write down the eigenvector describing this motion (no calculations).
Compute the eigenfrequency of the antisymmetric mode. Discuss if the eigenfrequency behaves as you expect when you would let the mass of the center atom $M$ go to infinity: $M \rightarrow \infty$.
The vibrational modes of molecules such as CO$_2$ can absorb light when the frequency of the light matches the eigenfrequency of the mode AND the atomic motion mode leads to a net, oscillating electric dipole. Do you expect the non-vibrating molecule to have an electric dipole moment? Do you expect any inhomogeneous electron distribution on the molecule? Argue if you would expect the symmetric or the antisymmetric mode to absorb light.

Hint

Compare your answers with Wikipedia, which also offers insight into the important role of absorption of solar radiation by the vibrational modes of molecules in our atmosphere, mainly occuring in the infrared part of the spectrum.

Exercise 2: Analyzing the Lennard-Jones potential¶

A simple model approximating the interaction between a pair of noble gas atoms such as Argon is the Lennard-Jones potential, in which the potential energy as a function of interatomic distance is

\[ U(r) = 4\epsilon\big[\big(\frac{\sigma}{r}\big)^{12}-\big(\frac{\sigma}{r}\big)^6\big] \]

where $r$ is the distance between two atoms, $\epsilon$ is the depth of the potential well, and $\sigma$ is the distance at which the inter-particle potential is zero.

Sketch $U(r)$ as a function of interatomic distance and mark the regions of repulsive and attractive forces acting between the atoms.
Find the distance, $r_0$ (bond length), at which the potential energy is minimal and find the value of the potential energy at this distance (binding energy of the molecule).
Expand $U(r)$ in a Taylor series around $r_0$ up to second order. By considering a second-order (=harmonic) potential approximation around the minimum ($r_0$), find an expression for the spring constant, $\kappa$, in terms of $\epsilon$ and $\sigma$.
Using the spring constant $\kappa$ you found earlier, find the ground state energy of the molecule by comparing the molecule to a quantum harmonic oscillator. What is the energy required to break the molecule apart?

Hint

Because the diatomic molecule is modeled as a one-body problem (in the center of mass rest frame of the molecule), the mass should be replaced by the reduced mass.
What is the approximate number of phonons that can occupy this mode before the potential becomes anharmonic?

Exercise 3: Numerical simulation of the vibrational modes of an $N$-atom chain¶

For this exercise, use Python (or any other programming environment you are comfortable with).

Define a matrix that relates forces to displacements in a linear 1D chain containing $N=5$ atoms. Repeat for $N=200$. You may assume that the masses and spring constants are equivalent throughout the chain.

Hint

In Python use the function numpy.diag.
Using numerical diagonalization (numpy.linalg.eigvalsh), compute the eigenfrequencies of this atomic chain. Plot a histogram of these eigenfrequencies.
Make the mass of every even atom different from the mass of every odd atom. Compute the eigenfrequencies of this atomic chain and plot a histogram.
Now make the masses of even and odd atoms equivalent again. Furthermore, make the spring constants of every even spring different from the odd spring. Compute the eigenfrequencies of this atomic chain and plot a histogram.

Exercise 4: Applying the LCAO model to a ring of 4 atoms¶

In this question we analyze the electronic spectrum of 4 atoms forming a ring-shaped molecule using the LCAO model. The ring is formed by the existence of a nearest-neighbour hopping $-t$ between atoms 1 and 2, between atoms 2 and 3, between atoms 3 and 4, and between atoms 4 and 1. The goal is to study the eigenstates of the system and the role of the symmetry of the system. The onsite energy is zero for all atoms.

Consider the wavefunction written as a linear combination of atomic orbitals $|\psi\rangle = \sum_i \phi_i |i\rangle$, where $i=1,2,3,4$. Argue why we should expect that $|\phi_1|^2 = |\phi_2|^2 = |\phi_3|^2 = |\phi_4|^2$ for all eigenvectors.
From the previous subquestion we conclude that the $\phi_i$ can only differ by a phase. Argue if you expect the phase difference between $\phi_1$ and $\phi_2$ to be the same or different from that between $\phi_2$ and $\phi_3$ for an eigenstate of the system.
Writing $\phi_{i+1} = e^{i\alpha}\phi_i$, identify the possible values of $\alpha$. Use these values to write down the 4 eigenvectors of the system.

Hint

Make use of the fact that the total change in phase should be an integer multiple of $2\pi$ going around the ring. In other words, we have periodic boundary conditions. (Large-$N$ atomic chains with periodic boundary conditions are the topic of the next lecture.)
Write down the LCAO Hamiltonian of the 4-atom ring, check if your eigenvectors are indeed eigenvectors, and find the corresponding eigenvalues. You should find 3 distinct eigenvalues, of which one has a double degeneracy. Argue why this degeneracy is to be expected given the form of the associated eigenvectors.
Now assume that the hopping between atoms 1 and 2 becomes $- t' $, and that the hopping between atoms 3 and 4 also becomes $-t'$. Which eigenenergies do you expect if $-t'=0$? What is the degeneracy of these eigenvalues?