6_bonds_and_spectra.md

from matplotlib import pyplot
from matplotlib.patches import Circle


import numpy as np

from common import draw_classic_axes, configure_plotting
from scipy.integrate import simps
import plotly.graph_objs as go
import plotly.offline as py

configure_plotting()

pi = np.pi

Interatomic interaction and molecular spectra

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

!!! success "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
- Define and diagonalize matrices numerically (for the exercises)

!!! summary "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

Adding repulsion

In the previous lecture, we observed that the diatomic molecule's 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.

First steps towards phonons

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

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

U = U_0 + \frac{\kappa}{2!} (r - a)^2 - \frac{\kappa_3}{3!} (r - a)³ + …

Up to second order, the potential is quadratic and gives rise to a harmonic equations of motion. The higher order term \kappa_3 provides anharmonicity.

Rigidity

We are interested in how materials respond to stretch and compression. To that end, assume we have a 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. Under such strain, the interatomic distance increases by \delta r = \delta L/N. We don't stretch too hard and therefore the interatomic distance $\delta r $ is small. For this reason, we approximate the interatomic potential around the minimum to be a quadratic potential (Taylor expansion 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}. This result directly leads us to a macroscopic material parameter, its compressibility:

\beta \equiv -\frac{1}{L} \frac{\partial L}{\partial F}= \frac{1}{\kappa a}.

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 finally arrive to the expression for the speed of sound, derived from the atomic properties: v = \sqrt{\frac{\kappa a^2}{m}}. This way we can already see how 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 =     a + b * r + c * r**2 + 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 = simps(np.exp(- U_cubic / E_t), r)
    r_avg = simps(r * np.exp(- U_cubic / E_t), 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)
                                                          
# Showing the figure
py.plot(fig)
fig.show();

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 focus on the mechanisms of 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. If one atom has a dipole moment \mathbf{p_1}, it creates an electric field

\mathbf{E} = \frac{\mathbf{p_1}}{4\pi \varepsilon_0 r^3}