Chapter 3: Thermal Conductivity

3.1 Introduction: Thermal Conductivity in Materials

Thermal conductivity is one of the most important transport properties of materials, governing heat dissipation in electronics, energy efficiency in thermoelectrics, and thermal management in diverse applications. In insulators and semiconductors, heat is primarily carried by phonons - the quantized lattice vibrations we studied in previous chapters.

Understanding phonon thermal conductivity requires connecting microscopic physics (phonon dispersion, scattering processes) to macroscopic observables (temperature gradient, heat flux). The Boltzmann transport equation provides this connection, offering a rigorous framework for calculating thermal conductivity from first principles.

3.2 Phonon Heat Current and Thermal Conductivity Tensor

3.2.1 Heat Current Density

When phonons with different wave vectors have different populations, they carry a net energy flux. The heat current density (heat flux per unit area) is:

In thermal equilibrium, \(n_{\mathbf{k}s}^0 = 1/(e^{\hbar\omega/k_BT} - 1)\) is isotropic, and the sum over all directions gives zero heat current. However, when a temperature gradient \(\nabla T\) is applied, the distribution becomes anisotropic, creating a net heat flux.

3.2.2 Fourier's Law and the Thermal Conductivity Tensor

For small temperature gradients, the heat current is proportional to the gradient:

For cubic crystals or polycrystalline materials, symmetry makes \(\boldsymbol{\kappa}\) diagonal and isotropic:

\[ \boldsymbol{\kappa} = \kappa \mathbb{I}, \quad \mathbf{J}_q = -\kappa \nabla T \]

where \(\kappa\) is the scalar thermal conductivity. For anisotropic materials (e.g., graphite, layered materials), different components \(\kappa_{xx}\), \(\kappa_{yy}\), \(\kappa_{zz}\) can differ by orders of magnitude.

3.3 The Boltzmann Transport Equation (BTE)

3.3.1 General Boltzmann Equation

The Boltzmann transport equation describes the time evolution of the phonon distribution function \(n_{\mathbf{k}s}(\mathbf{r}, t)\) under external forces and scattering processes:

For steady-state thermal conductivity in the absence of external forces:

\[ \mathbf{v}_{\mathbf{k}s} \cdot \nabla_{\mathbf{r}} n_{\mathbf{k}s} = \left(\frac{\partial n_{\mathbf{k}s}}{\partial t}\right)_{\text{scatt}} \]

3.3.2 The Relaxation Time Approximation (RTA)

The collision term is generally complex, involving all possible scattering processes. The relaxation time approximation simplifies this:

3.3.3 Solving the BTE: Linear Response

For small temperature gradients, we assume a linearized solution:

\[ n_{\mathbf{k}s} = n_{\mathbf{k}s}^0(T) + \delta n_{\mathbf{k}s} \]

where \(\delta n_{\mathbf{k}s}\) is the deviation from equilibrium. Since \(\nabla_{\mathbf{r}} n_{\mathbf{k}s}^0 = \frac{\partial n^0}{\partial T}\nabla T\):

\[ \mathbf{v}_{\mathbf{k}s} \cdot \nabla T \frac{\partial n_{\mathbf{k}s}^0}{\partial T} = -\frac{\delta n_{\mathbf{k}s}}{\tau_{\mathbf{k}s}} \]

Solving for \(\delta n_{\mathbf{k}s}\):

\[ \delta n_{\mathbf{k}s} = -\tau_{\mathbf{k}s} \mathbf{v}_{\mathbf{k}s} \cdot \nabla T \frac{\partial n_{\mathbf{k}s}^0}{\partial T} \]

Substituting into the heat current expression and using \(\frac{\partial n^0}{\partial T} = \frac{\hbar\omega}{k_BT^2}\frac{e^{\hbar\omega/k_BT}}{(e^{\hbar\omega/k_BT}-1)^2}\):

Converting the sum to an integral over the Brillouin zone and introducing the density of states \(g(\omega)\):

\[ \kappa = \int_0^{\omega_{\text{max}}} C(\omega) \, v^2(\omega) \, \tau(\omega) \, g(\omega) \, d\omega \]

3.4 Kinetic Theory Formula

3.4.1 Debye Model Approximation

For a simple isotropic Debye model with constant velocity \(v_s\) (speed of sound) and average relaxation time \(\bar{\tau}\):

This formula reveals the three factors determining thermal conductivity:

3.5 Scattering Mechanisms

The relaxation time (or mean free path) is limited by various scattering processes. Each mechanism contributes to the total scattering rate:

3.5.1 Matthiessen's Rule

3.5.2 Phonon-Phonon Scattering: Umklapp Processes

Anharmonic interactions cause phonons to scatter with each other. The most important for thermal resistance are umklapp processes (U-processes):

At low temperatures (\(T < \Theta_D/3\)), umklapp processes are exponentially suppressed because high-energy phonons needed for momentum non-conservation are rare.

3.5.3 Boundary Scattering

In finite samples, phonons scatter from external boundaries or grain boundaries:

Boundary scattering dominates at low temperatures or in nanostructured materials where intrinsic mean free paths exceed the sample dimensions.

3.5.4 Impurity and Defect Scattering

Point defects (substitutional atoms, vacancies, interstitials) scatter phonons via mass difference and force constant difference:

The \(\omega^4\) dependence means high-frequency phonons scatter more strongly, similar to Rayleigh scattering of light (why the sky is blue).

3.5.5 Isotope Scattering

Natural isotope mixtures create mass disorder. For example, natural silicon contains ~92% ²⁸Si, ~5% ²⁹Si, and ~3% ³⁰Si:

Isotopically pure materials (e.g., ²⁸Si, isotopically enriched diamond) show significantly higher thermal conductivity by eliminating this scattering channel.

3.5.6 Summary of Scattering Mechanisms

3.6 The Callaway Model

3.6.1 Beyond Simple RTA: Normal vs. Umklapp Processes

The simple RTA treats all scattering equally. However, normal processes (N-processes) conserve total phonon momentum and don't directly resist heat flow. Only umklapp processes and other resistive scattering actually limit thermal conductivity.

Callaway (1959) developed a model that separates resistive and non-resistive scattering:

The first term \(\kappa_1\) represents thermal conductivity if all scattering were resistive. The second term \(\kappa_2\) is a correction accounting for momentum-conserving N-processes.

3.7 Mean Free Path Distribution and Spectral Analysis

3.7.1 Phonon Mean Free Path Spectrum

Different phonon modes contribute differently to thermal conductivity. The mean free path distribution reveals which phonons carry most heat:

Experiments and calculations show that in many materials:

• At room temperature: 50% of heat is carried by phonons with \(\ell < 100\) nm
• Long mean free paths: Some phonons have \(\ell > 10\) μm, contributing significantly
• Nanostructuring: Reduces \(\kappa\) by blocking long-MFP phonons via boundary scattering

3.7.2 Ballistic vs. Diffusive Transport

When the sample dimension \(L\) becomes comparable to or smaller than the phonon mean free path \(\ell\), transport transitions from diffusive to ballistic:

3.8 Size Effects in Thermal Conductivity

Nanostructuring provides a powerful knob to tune thermal conductivity for applications like thermoelectrics (low \(\kappa\) desired) or thermal management (high \(\kappa\) desired).

3.8.1 Thin Films

For a thin film of thickness \(d\):

\[ \frac{\kappa_{\text{film}}}{\kappa_{\text{bulk}}} = \frac{\int_0^{\infty} \frac{C(\omega)v^2(\omega)\tau(\omega)}{1 + \tau(\omega)v(\omega)/d} g(\omega) d\omega}{\int_0^{\infty} C(\omega)v^2(\omega)\tau(\omega) g(\omega) d\omega} \]

Phonons with \(\ell = v\tau > d\) are suppressed, reducing \(\kappa_{\text{film}}\).

3.8.2 Nanowires and Nanotubes

For nanowires with diameter \(D\), boundary scattering at the surface limits MFP:

\[ \ell_B \approx D \]

Experiments on silicon nanowires show:

• \(\kappa\) can be reduced by 100× compared to bulk at \(D < 100\) nm
• Surface roughness enhances boundary scattering, further reducing \(\kappa\)
• ZT (thermoelectric figure of merit) dramatically improved

3.9 Minimum Thermal Conductivity

Is there a lower limit to thermal conductivity? Cahill and Pohl (1992) proposed a minimum \(\kappa\) based on the idea that phonon MFP cannot be shorter than half a wavelength (\(\ell_{\min} \sim \lambda/2\)):

Amorphous materials (glasses, amorphous silicon) and highly disordered crystals approach this minimum. It represents the thermal conductivity when phonons scatter after each lattice spacing - the most disordered possible state.

3.10 Python Implementation: Thermal Conductivity Calculation

Let's implement a thermal conductivity calculator using the Callaway model with realistic scattering mechanisms:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.constants import hbar, k as kB, pi

class PhononThermalConductivity:
    """
    Calculate thermal conductivity using Callaway model with multiple
    scattering mechanisms.
    """

    def __init__(self, material_params):
        """
        material_params: dict with keys
            - theta_D: Debye temperature (K)
            - v: average sound velocity (m/s)
            - M: average atomic mass (kg)
            - a: lattice constant (m)
            - gamma: Grüneisen parameter
            - V_atom: atomic volume (m³)
        """
        self.theta_D = material_params['theta_D']
        self.v = material_params['v']
        self.M = material_params['M']
        self.a = material_params['a']
        self.gamma = material_params['gamma']
        self.V_atom = material_params['V_atom']

    def tau_umklapp(self, omega, T):
        """Umklapp scattering relaxation time"""
        if T < 1e-6:
            return np.inf

        # Simplified umklapp rate: B*omega^2*T*exp(-theta_D/bT)
        b = 2.0  # typical value
        B = self.gamma**2 / (self.M * self.v**2 * self.theta_D)
        rate = B * omega**2 * T * np.exp(-self.theta_D / (b * T))
        return 1.0 / (rate + 1e-20)  # avoid division by zero

    def tau_boundary(self, omega, L):
        """Boundary scattering relaxation time"""
        # tau_B = L / v
        return L / self.v

    def tau_impurity(self, omega, Gamma_impurity):
        """Point defect scattering (Rayleigh scattering)"""
        # rate ~ Gamma * omega^4
        if Gamma_impurity < 1e-10:
            return np.inf
        rate = Gamma_impurity * omega**4 / (4 * pi * self.v**3)
        return 1.0 / rate

    def tau_isotope(self, omega, Gamma_isotope):
        """Isotope scattering"""
        if Gamma_isotope < 1e-10:
            return np.inf
        rate = Gamma_isotope * omega**4 / (4 * pi * self.v**3)
        return 1.0 / rate

    def tau_normal(self, omega, T):
        """Normal process rate (approximate)"""
        # Normal processes ~ omega^2 * T but don't conserve momentum
        if T < 1e-6:
            return np.inf
        B_N = self.gamma**2 / (self.M * self.v**2 * self.theta_D) * 0.5
        rate = B_N * omega**2 * T
        return 1.0 / (rate + 1e-20)

    def tau_combined(self, omega, T, L=np.inf, Gamma_imp=0, Gamma_iso=0):
        """Combined relaxation time (Matthiessen's rule)"""
        tau_U = self.tau_umklapp(omega, T)
        tau_B = self.tau_boundary(omega, L)
        tau_I = self.tau_impurity(omega, Gamma_imp)
        tau_iso = self.tau_isotope(omega, Gamma_iso)

        # Resistive scattering rate
        rate_R = 1/tau_U + 1/tau_B + 1/tau_I + 1/tau_iso

        return 1.0 / rate_R

    def thermal_conductivity_RTA(self, T, L=np.inf, Gamma_imp=0, Gamma_iso=0):
        """
        Calculate thermal conductivity using simple RTA
        """
        omega_D = self.theta_D * kB / hbar

        def integrand(x):
            # x = hbar*omega / (kB*T)
            omega = x * kB * T / hbar

            if omega > omega_D:
                return 0.0

            # Heat capacity factor
            C_factor = x**4 * np.exp(x) / (np.exp(x) - 1)**2

            # Relaxation time
            tau = self.tau_combined(omega, T, L, Gamma_imp, Gamma_iso)

            return tau * C_factor

        # Debye integral
        x_max = self.theta_D / T
        integral, _ = quad(integrand, 0, x_max, limit=100)

        # Prefactor: kB/(2*pi^2*v) * (kB*T/hbar)^3
        prefactor = kB / (2 * pi**2 * self.v) * (kB * T / hbar)**3

        kappa = prefactor * self.v**2 * integral

        return kappa

    def thermal_conductivity_Callaway(self, T, L=np.inf, Gamma_imp=0, Gamma_iso=0):
        """
        Calculate thermal conductivity using Callaway model
        (includes normal process correction)
        """
        omega_D = self.theta_D * kB / hbar

        def integrand_1(x):
            omega = x * kB * T / hbar
            if omega > omega_D:
                return 0.0

            tau_C = self.tau_combined(omega, T, L, Gamma_imp, Gamma_iso)
            C_factor = x**4 * np.exp(x) / (np.exp(x) - 1)**2

            return tau_C * C_factor

        def integrand_2_num(x):
            omega = x * kB * T / hbar
            if omega > omega_D:
                return 0.0

            tau_N = self.tau_normal(omega, T)
            C_factor = x**4 * np.exp(x) / (np.exp(x) - 1)**2

            return tau_N * C_factor

        def integrand_2_denom(x):
            omega = x * kB * T / hbar
            if omega > omega_D:
                return 0.0

            tau_N = self.tau_normal(omega, T)
            tau_R = self.tau_combined(omega, T, L, Gamma_imp, Gamma_iso)

            C_factor = x**4 * np.exp(x) / (np.exp(x) - 1)**2

            return (tau_N / tau_R) * C_factor

        x_max = self.theta_D / T

        # First term (RTA-like)
        I1, _ = quad(integrand_1, 0, x_max, limit=100)

        # Second term (normal process correction)
        I2_num, _ = quad(integrand_2_num, 0, x_max, limit=100)
        I2_denom, _ = quad(integrand_2_denom, 0, x_max, limit=100)

        # Avoid division by zero
        if I2_denom < 1e-20:
            I2_denom = 1e-20

        prefactor = kB / (2 * pi**2 * self.v) * (kB * T / hbar)**3

        kappa_1 = prefactor * self.v**2 * I1
        kappa_2 = prefactor * self.v**2 * I2_num**2 / I2_denom

        return kappa_1 + kappa_2


# Example: Silicon thermal conductivity
print("=" * 60)
print("Phonon Thermal Conductivity Calculator")
print("=" * 60)

# Silicon parameters
Si_params = {
    'theta_D': 645,      # Debye temperature (K)
    'v': 6400,           # Average sound velocity (m/s)
    'M': 28.085 * 1.66e-27,  # Atomic mass (kg)
    'a': 5.43e-10,       # Lattice constant (m)
    'gamma': 1.0,        # Grüneisen parameter
    'V_atom': (5.43e-10)**3 / 8  # Atomic volume (m³)
}

calculator = PhononThermalConductivity(Si_params)

# Calculate thermal conductivity vs temperature
T_range = np.logspace(1, 3, 50)  # 10 K to 1000 K
kappa_pure = []
kappa_isotope = []
kappa_nano = []

# Natural silicon isotope parameter
Gamma_Si_nat = 2e-4  # Natural isotope mixture

for T in T_range:
    # Pure silicon (no impurities, bulk)
    k1 = calculator.thermal_conductivity_Callaway(T)
    kappa_pure.append(k1)

    # Silicon with natural isotope mixture
    k2 = calculator.thermal_conductivity_Callaway(T, Gamma_iso=Gamma_Si_nat)
    kappa_isotope.append(k2)

    # Silicon nanowire (100 nm diameter)
    k3 = calculator.thermal_conductivity_Callaway(T, L=100e-9, Gamma_iso=Gamma_Si_nat)
    kappa_nano.append(k3)

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

# Linear scale
ax1.plot(T_range, kappa_pure, 'b-', linewidth=2, label='Pure Si (isotopically pure)')
ax1.plot(T_range, kappa_isotope, 'r--', linewidth=2, label='Natural Si (isotope scattering)')
ax1.plot(T_range, kappa_nano, 'g:', linewidth=2, label='Si nanowire (100 nm)')
ax1.set_xlabel('Temperature (K)', fontsize=12)
ax1.set_ylabel('Thermal Conductivity (W/m·K)', fontsize=12)
ax1.set_title('Silicon Thermal Conductivity (Callaway Model)', fontsize=14)
ax1.grid(True, alpha=0.3)
ax1.legend()
ax1.set_xlim(0, 1000)

# Log-log scale
ax2.loglog(T_range, kappa_pure, 'b-', linewidth=2, label='Pure Si')
ax2.loglog(T_range, kappa_isotope, 'r--', linewidth=2, label='Natural Si')
ax2.loglog(T_range, kappa_nano, 'g:', linewidth=2, label='Si nanowire (100 nm)')
ax2.set_xlabel('Temperature (K)', fontsize=12)
ax2.set_ylabel('Thermal Conductivity (W/m·K)', fontsize=12)
ax2.set_title('Log-Log Plot (Power Law Analysis)', fontsize=14)
ax2.grid(True, alpha=0.3, which='both')
ax2.legend()

plt.tight_layout()
plt.show()

# Print specific values
print("\nThermal Conductivity Values:")
print(f"{'Temperature (K)':<20} {'Pure Si':<15} {'Natural Si':<15} {'Si Nanowire':<15}")
print("-" * 65)
for i, T in enumerate([50, 100, 200, 300, 500]):
    idx = np.argmin(np.abs(T_range - T))
    print(f"{T:<20} {kappa_pure[idx]:<15.1f} {kappa_isotope[idx]:<15.1f} {kappa_nano[idx]:<15.1f}")

print("\nKey Observations:")
print("1. Low T (< 100 K): κ ~ T³ (boundary/defect scattering dominates)")
print("2. High T (> 300 K): κ ~ 1/T (umklapp scattering dominates)")
print("3. Isotope scattering reduces κ by ~10% at room temperature")
print("4. Nanostructuring (100 nm) reduces κ by ~10× (boundary scattering)")

3.11 Comparison with Experimental Data

The Callaway model with proper scattering parameters reproduces experimental thermal conductivity remarkably well:

Key experimental observations:

• Peak at low T: Most crystals show a peak where phonon population (C_V ~ T³) increases faster than scattering rates destroy MFP
• High-T 1/T decrease: Umklapp scattering grows linearly with T
• Isotope effects: Isotopically enriched samples show 50% higher κ at low T
• Size effects: Nanostructures show dramatic κ reduction and peak suppression

3.12 Advanced Topics and Current Research

3.12.1 Beyond the RTA: Full Solution of BTE

Modern computational methods solve the BTE exactly without RTA approximations:

• Iterative methods: Solve BTE by iteration until self-consistency
• Variational methods: Minimize functional to find optimal distribution
• Ab initio calculations: Calculate scattering rates from first principles (DFT + phonon code)

3.12.2 Hydrodynamic Phonon Transport

In ultra-pure crystals at intermediate temperatures, normal processes can dominate over umklapp. In this regime, phonons behave like a viscous fluid with collective flow:

Recently observed in graphite, graphene, and boron nitride at ~100 K.

3.12.3 Four-Phonon Scattering

Recent calculations show that four-phonon scattering (two phonons → two phonons) can be significant in materials with weak three-phonon scattering:

• Diamond: Four-phonon reduces κ by ~30% at 300 K
• Cubic BN: Similar effects observed
• Higher-order anharmonicity: Important for accurate predictions

3.13 Summary

3.14 Exercises