Phonon Density of States

1. Introduction: Why Density of States Matters

In the previous chapter, we learned about phonon dispersion relations \(\omega(\mathbf{q})\), which describe how phonon frequencies vary with wavevector throughout the Brillouin zone. However, for many thermodynamic properties (like specific heat), we don't need to know the exact frequency at every wavevector—we only need to know how many phonon modes exist at each frequency.

The phonon density of states (DOS) provides exactly this information. It's a fundamental quantity that connects microscopic lattice dynamics to macroscopic thermal properties of materials.

2. Definition of Density of States

2.1 Mathematical Definition

The phonon density of states \(g(\omega)\) is defined such that \(g(\omega)d\omega\) gives the number of phonon modes with frequencies between \(\omega\) and \(\omega + d\omega\).

Formally, for a crystal with \(N\) unit cells and \(s\) atoms per unit cell (giving \(3Ns\) normal modes total):

where:

• \(j\) indexes the phonon branches (3 acoustic + \(3s-3\) optical modes)
• The integral is over the first Brillouin zone
• \(\delta(\omega - \omega_j(\mathbf{q}))\) is the Dirac delta function, which "picks out" the \(\mathbf{q}\) points where the frequency equals \(\omega\)

2.2 Normalization

The DOS is normalized such that integrating over all frequencies gives the total number of modes:

This simply reflects the fact that a crystal with \(N\) unit cells of \(s\) atoms has exactly \(3Ns\) vibrational degrees of freedom (3 for each atom).

2.3 Relationship to Dispersion Relation

The DOS and dispersion relation contain the same information, but organized differently:

graph LR
  A[Dispersion ω q] -->|"Integration over BZ"| B[DOS g ω]
  B -->|"Thermal averaging"| C[Specific Heat C T]

  style A fill:#f9f,stroke:#333,stroke-width:2px
  style B fill:#bbf,stroke:#333,stroke-width:2px
  style C fill:#bfb,stroke:#333,stroke-width:2px
      

3. Van Hove Singularities

3.1 Physical Origin

When we examine the DOS of real crystals, we often observe sharp peaks or discontinuities called Van Hove singularities. These arise from the geometry of the dispersion relation in the Brillouin zone.

Van Hove singularities occur at frequencies where the dispersion relation has critical points:

At these points, the group velocity \(\mathbf{v}_g = \nabla_{\mathbf{q}} \omega\) vanishes. When many \(\mathbf{q}\) points contribute the same frequency (flat regions in dispersion), the DOS exhibits a singularity.

3.2 Types of Van Hove Singularities

In three dimensions, there are four types of critical points:

4. The Debye Model

4.1 Motivation and Assumptions

Calculating the exact DOS requires knowing the full dispersion relation throughout the Brillouin zone, which is complex for real materials. Peter Debye (1912) introduced a simple model that captures the essential low-frequency behavior.

Key assumptions:

1. Linear dispersion: All phonon branches are approximated as acoustic modes with \(\omega = v|\mathbf{q}|\), where \(v\) is an average sound velocity
2. Isotropic medium: The crystal is treated as an elastic continuum with no directional dependence
3. Debye cutoff: A maximum frequency \(\omega_D\) (Debye frequency) is introduced to preserve the correct total number of modes

4.2 Debye Density of States

For a linear dispersion \(\omega = vq\) in 3D, the number of modes in a shell of radius \(q\) to \(q + dq\) is proportional to the volume of the shell:

where \(V\) is the crystal volume and the factor of 3 accounts for three polarizations.

The DOS grows as \(\omega^2\) (quadratic) at low frequencies, which is a universal feature of 3D systems with linear dispersion.

4.3 Debye Cutoff Frequency

The Debye frequency \(\omega_D\) is determined by requiring the total number of modes to equal \(3N\):

Substituting the Debye DOS:

Solving for \(\omega_D\):

where \(q_D = (6\pi^2 n)^{1/3}\) is the Debye wavevector and \(n = N/V\) is the number density.

4.4 Debye Temperature

The Debye frequency is often expressed as a temperature via:

The Debye temperature \(\Theta_D\) is a material-specific constant that characterizes the stiffness of the lattice. Typical values:

5. The Einstein Model

5.1 Single Frequency Approximation

Albert Einstein (1907) proposed an even simpler model: all atoms vibrate independently at a single characteristic frequency \(\omega_E\) (Einstein frequency).

The Einstein DOS is simply:

This is equivalent to treating the solid as \(N\) independent quantum harmonic oscillators, all with the same frequency.

5.2 Einstein Temperature

Similar to Debye, we define:

5.3 When is Einstein Model Appropriate?

The Einstein model works best for:

• Optical phonons: These have nearly flat dispersion (weakly dependent on \(\mathbf{q}\)), so a single frequency is reasonable
• High temperatures: Where specific optical mode frequencies dominate
• Molecular crystals: Where intramolecular vibrations can be approximated as independent oscillators

6. Comparing Debye and Einstein Models

6.1 DOS Comparison

graph TB
  subgraph "Density of States"
    A[Debye Model: g ω ∝ ω²]
    B[Einstein Model: g ω = 3Nδ ω-ωE]
    C[Real Material: Van Hove singularities]
  end

  A -->|Good for| D[Low-frequency acoustic modes]
  B -->|Good for| E[Optical modes]
  C -->|Combines| F[Both acoustic and optical features]

  style A fill:#bbf,stroke:#333,stroke-width:2px
  style B fill:#fbf,stroke:#333,stroke-width:2px
  style C fill:#bfb,stroke:#333,stroke-width:2px
      

6.2 Complementary Nature

In practice, a combined approach often works best:

This uses the Debye model for the acoustic contribution (low-frequency, \(\omega^2\) dependence) and Einstein-like delta functions for each optical branch.

7. Real Materials: Beyond Simple Models

7.1 Calculating Real DOS

For real materials, the DOS is calculated from first-principles or force-constant models:

1. Generate a mesh of q-points in the Brillouin zone
2. Calculate \(\omega_j(\mathbf{q})\) at each point (from dynamical matrix)
3. Histogram the frequencies to build \(g(\omega)\)
4. Refine the mesh until converged (especially near Van Hove singularities)

7.2 Features of Real DOS

Real phonon DOS exhibits:

• Acoustic region: Near-Debye \(\omega^2\) behavior at low \(\omega\)
• Gap: Frequency gap between acoustic and optical branches in ionic/covalent materials
• Optical peaks: Sharp features from flat optical branches
• Van Hove singularities: Peaks and kinks from critical points

8. Computational Calculation of DOS

8.1 Python Implementation

Let's implement functions to calculate and visualize the Debye and Einstein models, and demonstrate how to histogram a real dispersion relation.

8.1.1 Debye Model DOS

import numpy as np
import matplotlib.pyplot as plt

def debye_dos(omega, omega_D, N=1):
  """
  Calculate Debye density of states.

  Parameters:
  -----------
  omega : array-like
    Frequency values (in arbitrary units)
  omega_D : float
    Debye cutoff frequency
  N : int, optional
    Number of atoms (for normalization)

  Returns:
  --------
  g_omega : array-like
    Debye DOS at each frequency
  """
  omega = np.asarray(omega)
  g = np.zeros_like(omega)

  # Debye DOS: g(ω) ∝ ω² for ω < ω_D
  mask = omega <= omega_D

  # Normalization factor to ensure ∫g(ω)dω = 3N
  normalization = 9 * N / omega_D**3

  g[mask] = normalization * omega[mask]**2

  return g

# Example usage
omega = np.linspace(0, 2, 1000)
omega_D = 1.5
g_debye = debye_dos(omega, omega_D, N=1)

plt.figure(figsize=(10, 6))
plt.plot(omega, g_debye, 'b-', linewidth=2, label='Debye model')
plt.axvline(omega_D, color='r', linestyle='--', label=f'$\omega_D$ = {omega_D}')
plt.xlabel('Frequency $\omega$ (arbitrary units)', fontsize=12)
plt.ylabel('DOS $g(\omega)$', fontsize=12)
plt.title('Debye Density of States', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

8.1.2 Einstein Model DOS

def einstein_dos(omega, omega_E, N=1, width=0.01):
  """
  Calculate Einstein density of states.

  Parameters:
  -----------
  omega : array-like
    Frequency values
  omega_E : float
    Einstein frequency
  N : int, optional
    Number of atoms
  width : float, optional
    Gaussian broadening width (to visualize the delta function)

  Returns:
  --------
  g_omega : array-like
    Einstein DOS (delta function approximated as Gaussian)
  """
  omega = np.asarray(omega)

  # Approximate delta function with narrow Gaussian
  # δ(ω - ω_E) ≈ (1/√(2πσ²)) exp(-(ω-ω_E)²/(2σ²))
  sigma = width
  g = (3 * N / (sigma * np.sqrt(2 * np.pi))) * np.exp(-0.5 * ((omega - omega_E) / sigma)**2)

  return g

# Example usage
omega_E = 1.0
g_einstein = einstein_dos(omega, omega_E, N=1, width=0.02)

plt.figure(figsize=(10, 6))
plt.plot(omega, g_einstein, 'r-', linewidth=2, label=f'Einstein model ($\omega_E$ = {omega_E})')
plt.xlabel('Frequency $\omega$ (arbitrary units)', fontsize=12)
plt.ylabel('DOS $g(\omega)$', fontsize=12)
plt.title('Einstein Density of States', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

8.1.3 Combined Model

def combined_dos(omega, omega_D, omega_E, N=1, width=0.02):
  """
  Combined Debye (acoustic) + Einstein (optical) model.

  Parameters:
  -----------
  omega : array-like
    Frequency values
  omega_D : float
    Debye cutoff frequency
  omega_E : float or list of float
    Einstein frequency(ies) for optical modes
  N : int
    Number of atoms
  width : float
    Gaussian width for Einstein peaks

  Returns:
  --------
  g_omega : array-like
    Combined DOS
  """
  # Debye contribution (acoustic modes)
  # Assume 3 acoustic modes out of 3s total (for s atoms per unit cell)
  g = debye_dos(omega, omega_D, N=N * 0.5)  # Reduced weight for acoustic

  # Einstein contribution (optical modes)
  if not isinstance(omega_E, (list, tuple)):
    omega_E = [omega_E]

  for w_E in omega_E:
    g += einstein_dos(omega, w_E, N=N * 0.5 / len(omega_E), width=width)

  return g

# Example: Material with acoustic + 2 optical branches
omega = np.linspace(0, 2.5, 1000)
omega_D = 1.0
omega_E_list = [1.5, 2.0]

g_combined = combined_dos(omega, omega_D, omega_E_list, N=1)

plt.figure(figsize=(12, 7))
plt.plot(omega, g_combined, 'g-', linewidth=2.5, label='Combined model')
plt.plot(omega, debye_dos(omega, omega_D, N=0.5), 'b--', alpha=0.6, label='Debye (acoustic)')
for i, w_E in enumerate(omega_E_list):
  plt.plot(omega, einstein_dos(omega, w_E, N=0.25, width=0.02),
       'r--', alpha=0.6, label=f'Einstein (optical {i+1})' if i == 0 else '')
plt.xlabel('Frequency $\omega$ (THz)', fontsize=12)
plt.ylabel('DOS $g(\omega)$', fontsize=12)
plt.title('Combined Debye-Einstein DOS Model', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

8.1.4 DOS from Dispersion Relation (Histogram Method)

def dos_from_dispersion(q_points, omega_q, num_bins=100):
  """
  Calculate DOS from a dispersion relation by histogramming.

  Parameters:
  -----------
  q_points : array-like
    Wavevector sampling points
  omega_q : array-like
    Frequencies at each q-point (can be 2D for multiple branches)
  num_bins : int
    Number of frequency bins

  Returns:
  --------
  omega_bins : array
    Frequency bin centers
  dos : array
    Density of states at each frequency
  """
  omega_flat = np.asarray(omega_q).flatten()

  # Create histogram
  hist, bin_edges = np.histogram(omega_flat, bins=num_bins, density=False)
  omega_bins = 0.5 * (bin_edges[1:] + bin_edges[:-1])

  # Normalize to number of modes
  bin_width = bin_edges[1] - bin_edges[0]
  dos = hist / bin_width / len(q_points)

  return omega_bins, dos

# Example: 1D monatomic chain dispersion
def monatomic_dispersion(q, C=1.0, M=1.0, a=1.0):
  """1D monatomic chain: ω(q) = 2√(C/M)|sin(qa/2)|"""
  return 2 * np.sqrt(C / M) * np.abs(np.sin(q * a / 2))

# Calculate dispersion
q = np.linspace(-np.pi, np.pi, 5000)
omega_1d = monatomic_dispersion(q)

# Calculate DOS
omega_bins, dos_1d = dos_from_dispersion(q, omega_1d, num_bins=100)

# Plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Dispersion
ax1.plot(q, omega_1d, 'b-', linewidth=2)
ax1.set_xlabel('Wavevector $q$ ($\pi$/a)', fontsize=12)
ax1.set_ylabel('Frequency $\omega$', fontsize=12)
ax1.set_title('1D Monatomic Chain Dispersion', fontsize=14)
ax1.grid(True, alpha=0.3)

# DOS
ax2.plot(dos_1d, omega_bins, 'r-', linewidth=2)
ax2.set_xlabel('DOS $g(\omega)$', fontsize=12)
ax2.set_ylabel('Frequency $\omega$', fontsize=12)
ax2.set_title('1D Monatomic Chain DOS', fontsize=14)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print("Note the Van Hove singularity at ω_max (divergent DOS)!")

9. Experimental Measurement of DOS

9.1 Inelastic Neutron Scattering (INS)

Principle: Neutrons scatter from atomic nuclei, exchanging energy and momentum with phonons. By measuring the scattered neutron energy and angle, both \(\mathbf{q}\) and \(\omega\) can be determined.

Advantages:

• Direct access to full dispersion \(\omega(\mathbf{q})\)
• Can measure all phonon branches
• DOS obtained by integrating over all \(\mathbf{q}\) at fixed \(\omega\)

Challenges:

• Requires neutron source (reactor or spallation source)
• Large sample volumes needed
• Expensive and time-consuming

9.2 Inelastic X-ray Scattering (IXS)

Similar to INS but uses high-energy X-rays. Synchrotron facilities enable high-resolution measurements, especially for heavy elements where neutron scattering is weak.

9.3 Raman Spectroscopy

Principle: Inelastic scattering of photons from phonons, typically measuring modes near the Brillouin zone center (\(\mathbf{q} \approx 0\)).

Advantages:

• Non-destructive, requires small samples
• Good for optical phonon frequencies
• Room-temperature measurements

Limitations:

• Only certain modes are Raman-active (symmetry selection rules)
• Typically limited to zone-center modes
• Does not give full DOS

9.4 Specific Heat Measurements

While not a direct DOS measurement, temperature-dependent specific heat \(C_V(T)\) can be inverted to extract \(g(\omega)\) using:

This integral equation can be solved (approximately) to extract \(g(\omega)\) from measured \(C_V(T)\) data across a wide temperature range.

10. Summary

graph TD
  A[Dispersion Relation ω q] --> B[Calculate DOS g ω]
  B --> C[Van Hove Singularities]
  B --> D[Debye Approximation Low ω]
  B --> E[Einstein Approximation Optical]
  D --> F[Specific Heat C T³]
  E --> G[High-T Specific Heat]
  C --> H[Experimental: INS/IXS]

  style A fill:#f9f,stroke:#333,stroke-width:2px
  style B fill:#bbf,stroke:#333,stroke-width:3px
  style F fill:#bfb,stroke:#333,stroke-width:2px
  style G fill:#bfb,stroke:#333,stroke-width:2px