Chapter 2: Anharmonic Phonons and Phase Transitions

1. Introduction: Beyond the Harmonic Approximation

The harmonic approximation, where atomic vibrations are treated as independent harmonic oscillators, provides an excellent starting point for understanding lattice dynamics. However, at finite temperatures and for many important physical phenomena, this approximation breaks down.

1.1 Failures of the Harmonic Approximation

The harmonic approximation fails to describe:

Thermal expansion: Harmonic crystals have zero thermal expansion
Phonon lifetimes: Harmonic phonons do not decay
Temperature-dependent frequencies: Harmonic phonons have constant frequencies
Structural phase transitions: Cannot capture soft modes and instabilities
High-temperature behavior: Underestimates entropy and heat capacity

Physical Origin of Anharmonicity

Anharmonic effects arise from the true potential energy surface deviating from a perfect quadratic form. The Taylor expansion of the potential includes cubic, quartic, and higher-order terms:

$$V = V_0 + \sum_i \frac{\partial V}{\partial u_i}u_i + \frac{1}{2}\sum_{ij}\frac{\partial^2 V}{\partial u_i \partial u_j}u_i u_j + \frac{1}{6}\sum_{ijk}\frac{\partial^3 V}{\partial u_i \partial u_j \partial u_k}u_i u_j u_k + \cdots$$

where $u_i$ represents atomic displacements. The cubic and quartic terms are responsible for anharmonic phenomena.

1.2 Temperature Effects on Phonons

At finite temperature, several phenomena emerge that require anharmonic treatment:

graph TD A[Temperature Increase] --> B[Larger Atomic Displacements] B --> C[Exploration of Anharmonic Region] C --> D[Phonon Frequency Shifts] C --> E[Phonon-Phonon Interactions] C --> F[Thermal Expansion] D --> G[Soft Mode Condensation] G --> H[Structural Phase Transition] E --> I[Phonon Linewidths] I --> J[Thermal Conductivity]

The key insight is that thermal motion causes atoms to explore anharmonic regions of the potential energy surface, leading to temperature-dependent effective force constants.

2. Anharmonic Effects at Finite Temperature

2.1 Quasi-Harmonic Approximation (QHA)

The simplest extension beyond the harmonic approximation is the quasi-harmonic approximation, which assumes:

Phonon frequencies depend on volume: $\omega_{\mathbf{q}\nu}(V)$
At each volume, the system is harmonic
Temperature dependence enters through thermal expansion

QHA Free Energy

$$F(V, T) = U_0(V) + F_{\text{vib}}(V, T)$$

$$F_{\text{vib}}(V, T) = k_B T \sum_{\mathbf{q}\nu} \ln\left[2\sinh\left(\frac{\hbar\omega_{\mathbf{q}\nu}(V)}{2k_B T}\right)\right]$$

The equilibrium volume at temperature $T$ minimizes $F(V, T)$:

$$\left.\frac{\partial F}{\partial V}\right|_{V_{\text{eq}}(T)} = 0$$

2.2 Intrinsic Anharmonicity

The QHA captures volume-dependent phonon frequencies but misses intrinsic anharmonic effects that occur even at constant volume:

Effect	Captured by QHA?	Physical Origin
Thermal expansion	✓ Yes	Volume-dependent phonons
Phonon-phonon scattering	✗ No	Cubic anharmonicity
Temperature-dependent frequencies (at constant V)	✗ No	Renormalization from thermal population
Soft mode hardening/softening	✗ No	Quartic anharmonicity
Negative thermal expansion (NTE)	Partial	Requires transverse acoustic modes

2.3 Phonon-Phonon Interactions

Cubic anharmonicity leads to three-phonon processes where one phonon can decay into two phonons or two phonons can merge into one:

Three-Phonon Scattering Rate

$$\Gamma_{\mathbf{q}\nu} = \frac{2\pi}{\hbar} \sum_{\mathbf{q}'\nu'\nu''} |V_3(\mathbf{q}\nu, \mathbf{q}'\nu', \mathbf{q}''\nu'')|^2$$

$$\times \left[(n_{\nu'}+n_{\nu''}+1)\delta(\omega-\omega'-\omega'') + (n_{\nu'}-n_{\nu''})\delta(\omega-\omega'+\omega'')\right]$$

where $V_3$ is the cubic anharmonic coupling and $n_{\nu} = 1/(e^{\hbar\omega/k_B T}-1)$ is the Bose-Einstein distribution.

The phonon linewidth is $\Gamma_{\mathbf{q}\nu}$, and the phonon lifetime is $\tau_{\mathbf{q}\nu} = 1/\Gamma_{\mathbf{q}\nu}$. These quantities are crucial for thermal conductivity.

3. Self-Consistent Phonon Theory

3.1 Basic Concept

Self-consistent phonon (SCP) theory addresses the limitation of the harmonic approximation by recognizing that thermal atomic vibrations renormalize the effective force constants. The key idea is:

Instead of expanding around the equilibrium lattice positions, expand around the thermally averaged positions, accounting for the thermal spread of atomic positions.

3.2 Theoretical Derivation

The SCP method constructs an effective harmonic Hamiltonian that best approximates the true anharmonic system at temperature $T$:

SCP Effective Hamiltonian

$$H_{\text{eff}} = \sum_i \frac{p_i^2}{2m_i} + \frac{1}{2}\sum_{ij}\Phi_{ij}^{\text{eff}}(T) u_i u_j$$

The effective force constants include thermal averaging over anharmonic terms:

$$\Phi_{ij}^{\text{eff}}(T) = \left\langle \frac{\partial^2 V}{\partial u_i \partial u_j} \right\rangle_T$$

where $\langle \cdots \rangle_T$ denotes thermal averaging.

For a potential including up to quartic terms:

$$V = \frac{1}{2}\sum_{ij}\Phi_{ij}^{(2)} u_i u_j + \frac{1}{6}\sum_{ijk}\Phi_{ijk}^{(3)} u_i u_j u_k + \frac{1}{24}\sum_{ijkl}\Phi_{ijkl}^{(4)} u_i u_j u_k u_l$$

The effective force constants become:

$$\Phi_{ij}^{\text{eff}}(T) = \Phi_{ij}^{(2)} + \frac{1}{2}\sum_{kl}\Phi_{ijkl}^{(4)}\langle u_k u_l \rangle_T$$

3.3 Self-Consistency Condition

The thermal average $\langle u_k u_l \rangle_T$ depends on the phonon frequencies, which in turn depend on $\Phi_{ij}^{\text{eff}}(T)$. This leads to a self-consistency equation:

SCP Self-Consistency Loop

Start with initial guess for $\Phi_{ij}^{\text{eff}}(T)$ (e.g., harmonic values)
Solve eigenvalue problem to get $\omega_{\mathbf{q}\nu}(T)$
Calculate thermal averages: $$\langle u_k u_l \rangle_T = \sum_{\mathbf{q}\nu} \frac{\hbar}{2m\omega_{\mathbf{q}\nu}}\left(2n_{\mathbf{q}\nu}+1\right) e_{k,\nu}(\mathbf{q}) e_{l,\nu}^*(\mathbf{q})$$
Update $\Phi_{ij}^{\text{eff}}(T)$ using the formula above
Repeat until convergence: $|\omega_{\mathbf{q}\nu}^{\text{new}} - \omega_{\mathbf{q}\nu}^{\text{old}}| < \epsilon$

3.4 Physical Interpretation

The SCP effective force constants have a clear physical interpretation:

graph LR A[Anharmonic Potential] --> B[Thermal Fluctuations] B --> C[Renormalized Force Constants] C --> D[Temperature-Dependent Frequencies] D --> E[Updated Thermal Fluctuations] E --> C style C fill:#e1f5ff

The quartic anharmonicity $\Phi_{ijkl}^{(4)}$ acts as a "correction" to the harmonic force constants, weighted by the thermal population of phonons. At high temperatures, larger thermal fluctuations $\langle u_k u_l \rangle_T$ lead to stronger renormalization.

4. Modern Anharmonic Methods

4.1 Self-Consistent Ab Initio Lattice Dynamics (SCAILD)

SCAILD combines DFT with SCP theory to compute temperature-dependent phonons from first principles without empirical parameters.

SCAILD Methodology

Generate displaced configurations: Create supercells with atomic displacements representing thermal motion
DFT force calculations: Compute forces for each configuration
Extract force constants: Fit anharmonic force constants $\Phi^{(3)}, \Phi^{(4)}$ from forces
SCP iteration: Solve self-consistency equations for $\Phi^{\text{eff}}(T)$
Phonon spectrum: Diagonalize dynamical matrix to get $\omega_{\mathbf{q}\nu}(T)$

4.2 Stochastic Self-Consistent Harmonic Approximation (SSCHA)

SSCHA is a variational method that goes beyond perturbative SCP theory. It determines the best harmonic approximation to the anharmonic free energy.

SSCHA Variational Principle

The trial harmonic free energy is:

$$F_{\mathcal{R}}(T) = \langle H \rangle_{\mathcal{R}} - T S_{\mathcal{R}}$$

where $\mathcal{R}$ represents the trial harmonic potential, $\langle H \rangle_{\mathcal{R}}$ is the expectation value of the true anharmonic Hamiltonian with respect to the harmonic density matrix, and $S_{\mathcal{R}}$ is the entropy of the harmonic system.

SSCHA minimizes $F_{\mathcal{R}}(T)$ with respect to the trial force constants and equilibrium positions.

Advantages of SSCHA over SCP:

Fully variational (satisfies thermodynamic inequalities)
Handles strong anharmonicity (e.g., near phase transitions)
Includes quantum nuclear effects naturally
Can optimize both force constants and lattice structure
Provides access to free energy derivatives (stress tensor, heat capacity)

4.3 Temperature-Dependent Effective Potential (TDEP)

TDEP extracts temperature-dependent force constants from molecular dynamics or Monte Carlo simulations:

TDEP Workflow

Ab initio MD/MC: Run first-principles simulations at target temperature
Collect snapshots: Extract atomic positions and forces from trajectory
Covariance analysis: Compute position covariance matrix $\langle u_i u_j \rangle$
Force constant fitting: Determine $\Phi_{ij}^{\text{eff}}(T)$ that best reproduces forces and covariances
Phonon calculation: Diagonalize to obtain temperature-dependent spectrum

TDEP is particularly powerful for systems where perturbative methods fail, such as superionic conductors or highly anharmonic materials.

4.4 Method Comparison

Method	Computational Cost	Accuracy	Best For	Software
QHA	Low	Good for weak anharmonicity	Thermal expansion, simple materials	Phonopy, QE
SCP/SCAILD	Medium	Good for moderate anharmonicity	Temperature-dependent phonons	ALAMODE, Phono3py
SSCHA	Medium-High	Excellent for strong anharmonicity	Phase transitions, soft modes	python-sscha
TDEP	High	Very high (includes all effects)	Complex anharmonicity, liquid-like	TDEP package

5. Temperature-Dependent Phonon Frequencies

5.1 General Temperature Dependence

Phonon frequencies typically exhibit temperature dependence of the form:

$$\omega_{\mathbf{q}\nu}(T) = \omega_{\mathbf{q}\nu}(0) + \Delta\omega_{\mathbf{q}\nu}(T)$$

where the frequency shift $\Delta\omega_{\mathbf{q}\nu}(T)$ can be decomposed into:

$$\Delta\omega_{\mathbf{q}\nu}(T) = \Delta\omega_{\mathbf{q}\nu}^{\text{V}}(T) + \Delta\omega_{\mathbf{q}\nu}^{\text{anh}}(T)$$

$\Delta\omega^{\text{V}}$: Volume contribution (captured by QHA)
$\Delta\omega^{\text{anh}}$: Intrinsic anharmonic contribution

5.2 Grüneisen Parameters

The mode Grüneisen parameter quantifies how phonon frequency changes with volume:

$$\gamma_{\mathbf{q}\nu} = -\frac{V}{\omega_{\mathbf{q}\nu}}\frac{\partial \omega_{\mathbf{q}\nu}}{\partial V}$$

Positive $\gamma_{\mathbf{q}\nu}$ means frequency decreases with volume expansion (most common). Negative $\gamma_{\mathbf{q}\nu}$ can lead to negative thermal expansion.

5.3 Experimental Observation

Temperature-dependent phonon frequencies can be measured by:

Inelastic neutron scattering (INS): Direct measurement of $\omega_{\mathbf{q}\nu}(T)$
Raman spectroscopy: Zone-center optical phonons ($\mathbf{q} = 0$)
Infrared spectroscopy: IR-active optical modes
Inelastic X-ray scattering (IXS): High-energy resolution for acoustic phonons

Example: Silicon

In silicon, the optical phonon at the zone center ($\Gamma$ point) softens with temperature:

$$\omega_{\Gamma}(T) \approx 520 \text{ cm}^{-1} - 0.02 \text{ cm}^{-1}\text{/K} \times T$$

This softening is primarily due to volume expansion (positive Grüneisen parameter $\gamma \approx 0.98$) and cubic anharmonicity (phonon-phonon scattering).

6. Soft Modes and Structural Instabilities

6.1 Definition and Physical Meaning

A soft mode is a phonon whose frequency decreases with temperature and potentially vanishes at a critical temperature $T_c$:

$$\omega_{\mathbf{q}_0\nu}(T) \to 0 \quad \text{as} \quad T \to T_c$$

Physical Significance

A soft mode indicates a structural instability. When $\omega^2 < 0$ (imaginary frequency), the crystal structure is unstable and will undergo a phase transition to a lower-symmetry structure.

6.2 Origin of Soft Modes

Soft modes arise when anharmonic interactions cause effective force constants to decrease with temperature:

$$\Phi_{ij}^{\text{eff}}(T) = \Phi_{ij}^{(2)} + \frac{1}{2}\sum_{kl}\Phi_{ijkl}^{(4)}\langle u_k u_l \rangle_T$$

If $\Phi_{ijkl}^{(4)} < 0$ (negative quartic coupling) and $\langle u_k u_l \rangle_T$ increases with $T$, then $\Phi_{ij}^{\text{eff}}(T)$ can become negative, leading to imaginary phonon frequencies.

6.3 Landau-Peierls Instability

For a mode at wavevector $\mathbf{q}_0$, the free energy can be expanded as:

$$F = F_0 + \frac{1}{2}A(T) Q^2 + \frac{1}{4}B Q^4 + \frac{1}{6}C Q^6 + \cdots$$

where $Q$ is the amplitude of the soft mode, and:

$A(T) = a(T - T_c)$ with $a > 0$: quadratic coefficient
$B$: quartic coefficient (can be positive or negative)
$C > 0$: sextic coefficient (ensures stability)

The equilibrium displacement $Q_{\text{eq}}$ minimizes the free energy:

$$\frac{\partial F}{\partial Q}\bigg|_{Q_{\text{eq}}} = 0 \quad \Rightarrow \quad A Q_{\text{eq}} + B Q_{\text{eq}}^3 + C Q_{\text{eq}}^5 = 0$$

6.4 Identifying Soft Modes Computationally

Calculate phonon spectrum at 0 K: Check for imaginary frequencies
Perform SCP/SSCHA at finite T: Track temperature-dependent softening
Identify critical temperature: Find $T_c$ where $\omega_{\mathbf{q}_0\nu}(T_c) = 0$
Analyze mode eigenvector: Determine atomic displacement pattern
Predict new structure: Apply displacement pattern to find low-symmetry phase

Example: BaTiO₃ Ferroelectric Transition

In cubic BaTiO₃, a soft transverse optical (TO) mode at the $\Gamma$ point drives the ferroelectric phase transition at $T_c \approx 403$ K:

High-T phase (T > 403 K): Cubic, $Pm\overline{3}m$, $\omega_{\text{TO}} > 0$
Transition: TO mode softens, $\omega_{\text{TO}}(T) \propto \sqrt{T - T_c}$
Low-T phase (T < 403 K): Tetragonal, $P4mm$, Ti displacement along [001]

The soft mode corresponds to Ti displacement relative to the oxygen cage, creating spontaneous polarization.

7. Landau Theory of Phase Transitions

7.1 Order Parameter

Landau theory describes phase transitions using an order parameter $\eta$ that:

Is zero in the high-symmetry phase
Becomes non-zero in the low-symmetry phase
Transforms according to the irreducible representation of the soft mode

For structural phase transitions driven by phonon softening, the order parameter is the amplitude of the soft mode displacement pattern.

7.2 Landau Free Energy

General Landau Expansion

$$F(\eta, T) = F_0(T) + \frac{1}{2}a(T-T_c)\eta^2 + \frac{1}{4}b\eta^4 + \frac{1}{6}c\eta^6$$

Equilibrium condition: $\frac{\partial F}{\partial \eta} = 0$

$$a(T-T_c)\eta + b\eta^3 + c\eta^5 = 0$$

7.3 Second-Order Phase Transitions

When $b > 0$ (positive quartic coefficient), the transition is second-order (continuous):

For $T > T_c$: $\eta = 0$ (high-symmetry phase)

For $T < T_c$: $\eta = \pm\sqrt{\frac{a(T_c - T)}{b}}$ (low-symmetry phase)

The order parameter grows continuously from zero as $\eta \propto \sqrt{T_c - T}$.

7.4 First-Order Phase Transitions

When $b < 0$ (negative quartic coefficient), the sextic term $c\eta^6$ with $c > 0$ is needed for stability. The transition becomes first-order (discontinuous):

The transition occurs at $T_t$ where the free energies of the two phases are equal:

$$F(\eta=0, T_t) = F(\eta=\eta_{\text{eq}}, T_t)$$

The order parameter jumps discontinuously: $\eta = 0 \to \eta_{\text{eq}} \neq 0$

7.5 Critical Exponents

Near second-order transitions, physical quantities exhibit power-law behavior characterized by critical exponents:

Quantity	Critical Behavior	Landau Theory Exponent
Order parameter	$\eta \propto (T_c - T)^{\beta}$	$\beta = 1/2$
Susceptibility	$\chi \propto \|T - T_c\|^{-\gamma}$	$\gamma = 1$
Heat capacity	$C \propto \|T - T_c\|^{-\alpha}$	$\alpha = 0$ (jump)
Correlation length	$\xi \propto \|T - T_c\|^{-\nu}$	$\nu = 1/2$

Note: Landau theory is mean-field and does not capture fluctuation effects. Actual critical exponents can differ (e.g., 3D Ising model: $\beta = 0.326$, $\gamma = 1.237$).

8. Displacive vs Order-Disorder Transitions

8.1 Classification of Structural Phase Transitions

Structural phase transitions can be classified based on the microscopic mechanism:

Feature	Displacive	Order-Disorder
Mechanism	Continuous softening of phonon mode	Ordering of pre-existing local distortions
High-T phase	Dynamically stable (positive $\omega^2$)	Locally distorted, time-averaged symmetric
Low-T phase	Static displacement along soft mode	Frozen-in ordering of local distortions
Order parameter	Collective displacement amplitude	Correlation of local displacements
Temperature dependence	Soft mode frequency $\omega(T) \to 0$	Relaxation time $\tau(T) \to \infty$
Typical examples	BaTiO₃, SrTiO₃, quartz	KH₂PO₄, Rochelle salt, some perovskites

8.2 Displacive Transitions

In a displacive transition, the high-temperature phase has a single potential energy minimum (undistorted structure), but the curvature (force constant) decreases with temperature due to anharmonicity.

graph TD A[High Temperature] --> B[Soft phonon mode] B --> C[ω decreases with T] C --> D[ω → 0 at Tc] D --> E[Potential energy minimum splits] E --> F[Atoms displace to new equilibrium] F --> G[Low-symmetry phase]

Characteristics:

Continuous softening of a specific phonon mode
Well-defined soft mode above $T_c$
Often second-order transitions (but can be first-order if $b < 0$)
Central peak absent or weak in scattering experiments

8.3 Order-Disorder Transitions

In an order-disorder transition, the potential energy surface has multiple local minima even above $T_c$. Atoms hop between these minima, and the high-temperature phase is symmetric only on average.

Order-Disorder Mechanism

Above $T_c$: Atoms rapidly hop between equivalent off-center positions. Time average gives symmetric structure.

Below $T_c$: Hopping slows down ($\tau \to \infty$). Atoms freeze into one of the equivalent positions, breaking symmetry.

Characteristics:

Local distortions present above $T_c$
No well-defined soft mode (overdamped central peak in scattering)
Transition driven by entropy (ordering reduces entropy)
Strong central peak in neutron/Raman scattering

8.4 Intermediate Cases

Many real materials exhibit mixed character:

BaTiO₃: Primarily displacive (well-defined soft mode) with some order-disorder character
PbTiO₃: Strongly displacive
KNbO₃: Mixed character depending on temperature range
Relaxor ferroelectrics: Strongly order-disorder (nanoscale polar regions)

The distinction can be quantified by the ratio of the soft mode frequency to the damping rate:

$$\text{Displacive: } \omega_0 \gg \Gamma \quad \text{(underdamped)}$$

$$\text{Order-Disorder: } \omega_0 \ll \Gamma \quad \text{(overdamped)}$$

9. Case Studies: Perovskites and NTE Materials

9.1 Perovskite Phase Transitions: BaTiO₃

BaTiO₃ is the prototypical ferroelectric perovskite with multiple phase transitions:

Temperature	Phase	Space Group	Ti Displacement
$T > 403$ K	Cubic (paraelectric)	$Pm\overline{3}m$	None (centered)
278 K < $T$ < 403 K	Tetragonal (ferroelectric)	$P4mm$	Along [001]
183 K < $T$ < 278 K	Orthorhombic (ferroelectric)	$Amm2$	Along [011]
$T < 183$ K	Rhombohedral (ferroelectric)	$R3m$	Along [111]

Computational Analysis of BaTiO₃

First-principles calculations reveal:

Cubic phase instability: Phonon calculation at 0 K shows imaginary frequencies (soft modes) at $\Gamma$, $X$, $M$, and $R$ points
Dominant instability: The $\Gamma$-point TO mode (Ti displacement) has the largest imaginary frequency
Anharmonic stabilization: SSCHA calculations show that the cubic phase is stabilized at high temperature by anharmonic effects
Temperature-dependent soft mode: $\omega_{\text{TO}}(T) \approx A\sqrt{T - T_c}$ above $T_c$

9.2 SrTiO₃: Quantum Paraelectric

SrTiO₃ is a fascinating case where quantum fluctuations prevent ferroelectric ordering:

Quantum Paraelectric Behavior

At low temperature, the soft mode frequency saturates at a finite value ($\sim 40$ cm$^{-1}$) instead of going to zero, due to quantum zero-point motion:

$$\omega_{\text{TO}}(T \to 0) = \omega_{\text{quantum}} > 0$$

This prevents the phase transition. SrTiO₃ remains cubic down to 0 K (hence "quantum paraelectric").

SrTiO₃ exhibits other interesting phenomena:

Antiferrodistortive transition at 105 K: Oxygen octahedra rotate, cubic → tetragonal
Strain-induced ferroelectricity: Epitaxial strain can induce ferroelectric ordering
Isotope effect: Replacing $^{16}$O with $^{18}$O increases the soft mode frequency (stronger quantum effects)

9.3 Negative Thermal Expansion: ZrW₂O₈

ZrW₂O₈ exhibits negative thermal expansion (NTE) over a wide temperature range (0.3–1050 K):

$$\alpha_V = \frac{1}{V}\frac{dV}{dT} \approx -9 \times 10^{-6} \text{ K}^{-1}$$

Phonon Mechanism of NTE

NTE arises from low-frequency transverse acoustic phonons with negative Grüneisen parameters:

The volumetric thermal expansion coefficient is:

$$\alpha_V = \frac{1}{BV}\sum_{\mathbf{q}\nu} \gamma_{\mathbf{q}\nu} C_{\mathbf{q}\nu}$$

where $B$ is bulk modulus, $\gamma_{\mathbf{q}\nu}$ is the mode Grüneisen parameter, and $C_{\mathbf{q}\nu}$ is the mode heat capacity.

For NTE: dominant modes have $\gamma_{\mathbf{q}\nu} < 0$

Physical picture: Transverse vibrations of Zr-O-W linkages cause the framework to contract on average. As temperature increases, stronger transverse vibrations lead to greater contraction.

Computational Prediction of NTE

Calculate full phonon spectrum and eigenvectors
Compute mode Grüneisen parameters: $\gamma_{\mathbf{q}\nu} = -\frac{V}{\omega_{\mathbf{q}\nu}}\frac{\partial \omega_{\mathbf{q}\nu}}{\partial V}$
Identify modes with large negative $\gamma_{\mathbf{q}\nu}$
Calculate thermal expansion using QHA or anharmonic methods
Verify against experimental dilatometry data

Other NTE materials: ScF₃, Ag₃[Co(CN)₆], MOF-5 (metal-organic framework)

10. Computational Implementation

10.1 Workflow for Anharmonic Phonon Calculations

graph TD A[DFT Ground State] --> B[Harmonic Phonons] B --> C{Stable?} C -->|Yes| D[QHA Calculation] C -->|No - Imaginary ω| E[Anharmonic Calculation] D --> F[Temperature-dependent properties] E --> G{Method Choice} G --> H[SCAILD - Moderate anharmonicity] G --> I[SSCHA - Strong anharmonicity] G --> J[TDEP - Complex systems] H --> K[SCP iterations] I --> L[Free energy minimization] J --> M[MD/MC simulations] K --> F L --> F M --> F

10.2 Python Implementation: SCP Calculation

Here is a simplified implementation of the self-consistent phonon method:

Python scp_phonons.py - Self-Consistent Phonon Calculation

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh

class SCPCalculator:
    """
    Self-consistent phonon (SCP) calculator for anharmonic systems.

    This implementation demonstrates the core SCP algorithm for a
    1D chain with quartic anharmonicity.
    """

    def __init__(self, natoms, mass, phi2, phi4):
        """
        Initialize SCP calculator.

        Parameters
        ----------
        natoms : int
            Number of atoms in unit cell
        mass : float
            Atomic mass (in atomic units)
        phi2 : ndarray
            Harmonic force constant matrix (natoms x natoms)
        phi4 : ndarray
            Quartic anharmonic coupling tensor (natoms x natoms x natoms x natoms)
        """
        self.natoms = natoms
        self.mass = mass
        self.phi2 = phi2
        self.phi4 = phi4
        self.kB = 8.617e-5  # Boltzmann constant in eV/K

    def thermal_displacement_covariance(self, omega, eigvec, temperature):
        """
        Calculate thermal average <u_i u_j> for each mode.

        Parameters
        ----------
        omega : ndarray
            Phonon frequencies (nmodes,)
        eigvec : ndarray
            Phonon eigenvectors (natoms, nmodes)
        temperature : float
            Temperature in Kelvin

        Returns
        -------
        cov_matrix : ndarray
            Displacement covariance matrix <u_i u_j> (natoms, natoms)
        """
        nmodes = len(omega)
        cov_matrix = np.zeros((self.natoms, self.natoms))

        hbar = 0.6582119569  # eV·ps

        for mode in range(nmodes):
            if omega[mode] < 1e-6:  # Skip zero or imaginary modes
                continue

            # Bose-Einstein distribution
            if temperature > 0:
                n_BE = 1.0 / (np.exp(hbar * omega[mode] / (self.kB * temperature)) - 1)
            else:
                n_BE = 0.0

            # Quantum + thermal contribution
            prefactor = (hbar / (2 * self.mass * omega[mode])) * (2 * n_BE + 1)

            # Add contribution from this mode
            for i in range(self.natoms):
                for j in range(self.natoms):
                    cov_matrix[i, j] += prefactor * eigvec[i, mode] * eigvec[j, mode]

        return cov_matrix

    def effective_force_constants(self, cov_matrix):
        """
        Calculate temperature-dependent effective force constants.

        Φ_eff = Φ^(2) + (1/2) Σ_kl Φ^(4)_ijkl <u_k u_l>

        Parameters
        ----------
        cov_matrix : ndarray
            Thermal displacement covariance (natoms, natoms)

        Returns
        -------
        phi_eff : ndarray
            Effective force constant matrix (natoms, natoms)
        """
        phi_eff = self.phi2.copy()

        # Add quartic contribution
        for i in range(self.natoms):
            for j in range(self.natoms):
                for k in range(self.natoms):
                    for l in range(self.natoms):
                        phi_eff[i, j] += 0.5 * self.phi4[i, j, k, l] * cov_matrix[k, l]

        return phi_eff

    def solve_phonons(self, phi_eff):
        """
        Solve phonon eigenvalue problem.

        Parameters
        ----------
        phi_eff : ndarray
            Effective force constant matrix

        Returns
        -------
        omega : ndarray
            Phonon frequencies (in THz or energy units)
        eigvec : ndarray
            Phonon eigenvectors (normalized)
        """
        # Construct dynamical matrix
        dyn_matrix = phi_eff / self.mass

        # Solve eigenvalue problem
        eigenvalues, eigenvectors = eigh(dyn_matrix)

        # Convert eigenvalues to frequencies
        # ω² = eigenvalue, so ω = sqrt(eigenvalue)
        omega = np.sqrt(np.abs(eigenvalues))
        omega[eigenvalues < 0] *= -1  # Keep track of imaginary frequencies

        return omega, eigenvectors

    def run_scp(self, temperature, max_iter=100, tol=1e-6):
        """
        Self-consistent phonon iteration loop.

        Parameters
        ----------
        temperature : float
            Temperature in Kelvin
        max_iter : int
            Maximum number of iterations
        tol : float
            Convergence tolerance for frequencies

        Returns
        -------
        omega : ndarray
            Converged phonon frequencies
        eigvec : ndarray
            Converged phonon eigenvectors
        iterations : int
            Number of iterations required
        """
        print(f"\nSCP calculation at T = {temperature} K")
        print("=" * 50)

        # Initial guess: harmonic phonons
        omega_old, eigvec = self.solve_phonons(self.phi2)

        for iteration in range(max_iter):
            # Step 1: Calculate thermal displacement covariance
            cov_matrix = self.thermal_displacement_covariance(
                omega_old, eigvec, temperature
            )

            # Step 2: Update effective force constants
            phi_eff = self.effective_force_constants(cov_matrix)

            # Step 3: Solve for new phonon frequencies
            omega_new, eigvec = self.solve_phonons(phi_eff)

            # Check convergence
            max_diff = np.max(np.abs(omega_new - omega_old))

            if iteration % 10 == 0:
                print(f"Iteration {iteration:3d}: max Δω = {max_diff:.6e}")

            if max_diff < tol:
                print(f"\nConverged in {iteration + 1} iterations")
                print(f"Final frequencies: {omega_new}")
                return omega_new, eigvec, iteration + 1

            omega_old = omega_new

        print(f"\nWarning: Did not converge in {max_iter} iterations")
        return omega_new, eigvec, max_iter


def example_1d_chain():
    """
    Example: 1D harmonic chain with quartic anharmonicity.

    Model potential: V = (1/2)k₂(u_{i+1} - u_i)² - (1/4)k₄(u_{i+1} - u_i)⁴
    Negative k₄ causes softening at high temperature.
    """
    print("SCP Example: 1D Chain with Anharmonicity")
    print("=" * 60)

    # System parameters
    natoms = 1  # Single atom per unit cell for simplicity
    mass = 12.0  # Carbon-like mass
    k2 = 5.0     # Harmonic force constant (eV/Å²)
    k4 = -0.5    # Quartic anharmonicity (eV/Å⁴), negative causes softening

    # Construct force constant matrices
    phi2 = np.array([[k2]])
    phi4 = np.array([[[[k4]]]])

    # Initialize SCP calculator
    scp = SCPCalculator(natoms, mass, phi2, phi4)

    # Temperature range
    temperatures = np.linspace(0, 1000, 11)
    frequencies = []

    # Run SCP at each temperature
    for T in temperatures:
        omega, _, _ = scp.run_scp(T, max_iter=50, tol=1e-6)
        frequencies.append(omega[0])

    # Plot results
    plt.figure(figsize=(10, 6))
    plt.plot(temperatures, frequencies, 'o-', linewidth=2, markersize=8)
    plt.axhline(y=0, color='r', linestyle='--', label='Instability threshold')
    plt.xlabel('Temperature (K)', fontsize=14)
    plt.ylabel('Phonon Frequency (THz)', fontsize=14)
    plt.title('Temperature-Dependent Phonon Softening (SCP)', fontsize=16)
    plt.grid(True, alpha=0.3)
    plt.legend(fontsize=12)
    plt.tight_layout()
    plt.savefig('scp_temperature_dependence.png', dpi=300)
    plt.show()

    print("\n" + "=" * 60)
    print("Results:")
    print(f"Harmonic frequency (T=0): {frequencies[0]:.4f} THz")
    print(f"Frequency at 1000 K: {frequencies[-1]:.4f} THz")
    print(f"Softening: {frequencies[0] - frequencies[-1]:.4f} THz")


if __name__ == "__main__":
    example_1d_chain()

10.3 Using ALAMODE for SCP Calculations

ALAMODE (Anharmonic Lattice MODElling) is a powerful software for anharmonic phonon calculations.

Bash ALAMODE workflow for SCP

# Step 1: Generate displaced configurations
# Create ALAMODE input file (alm_input.in)

# Step 2: Run DFT calculations for each displaced configuration
# (Use VASP, Quantum ESPRESSO, etc.)

# Step 3: Extract force constants using ALM
alm alm_input.in

# Step 4: Run self-consistent phonon calculation with ANPHON
anphon anphon_input.in

# anphon_input.in should specify:
# - Temperature range
# - SCP mode (SELF_CONSISTENT = 1)
# - Convergence criteria
# - Output options (phonon DOS, thermal properties)

# Step 5: Analyze results
# - Temperature-dependent phonon dispersion
# - Free energy
# - Heat capacity
# - Thermal expansion

10.4 Using SSCHA

Python SSCHA basic workflow

import cellconstructor as CC
import cellconstructor.Phonons
import sscha, sscha.Ensemble, sscha.SchaMinimizer

# Load harmonic dynamical matrices
harmonic_dyn = CC.Phonons.Phonons("harmonic_dyn", nqirr=8)

# Define temperature
temperature = 300  # Kelvin

# Create SSCHA ensemble
ensemble = sscha.Ensemble.Ensemble(harmonic_dyn, temperature)

# Generate configurations
ensemble.generate(N_configs=1000)

# Compute energies and forces with DFT
# (This step requires interfacing with DFT code)
ensemble.get_energy_forces(calculator=my_dft_calculator)

# Initialize SSCHA minimizer
minimizer = sscha.SchaMinimizer.SSCHA_Minimizer(ensemble)

# Set minimization parameters
minimizer.min_step_dyn = 0.01
minimizer.kong_liu_ratio = 0.5
minimizer.gradi_op = "all"

# Run SSCHA minimization
minimizer.run()

# Extract results
sscha_dyn = minimizer.dyn
sscha_free_energy = minimizer.get_free_energy()

# Compute phonon dispersion with SSCHA correction
sscha_dyn.Symmetrize()
sscha_dyn.save_qe("sscha_dyn")

10.5 Practical Recommendations

Scenario	Recommended Method	Software
Simple material, weak anharmonicity	QHA	Phonopy + VASP/QE
Temperature-dependent phonons, moderate anharmonicity	SCP/SCAILD	ALAMODE, Phono3py
Soft modes, near phase transition	SSCHA	python-sscha
High-temperature, liquid-like behavior	TDEP	TDEP package
Thermal conductivity	Phono3py	Phono3py + VASP/QE

11. Summary

Key Concepts

Anharmonic effects become important at finite temperature, leading to thermal expansion, phonon-phonon scattering, and temperature-dependent frequencies
Self-consistent phonon (SCP) theory renormalizes force constants by thermally averaging anharmonic terms, capturing intrinsic temperature dependence
Modern methods (SCAILD, SSCHA, TDEP) extend beyond perturbative SCP to handle strong anharmonicity and quantum effects
Soft modes signal structural instabilities; when $\omega \to 0$, a phase transition occurs
Landau theory describes phase transitions using order parameters and free energy expansions, predicting first- or second-order behavior
Displacive vs order-disorder transitions differ in their microscopic mechanisms: continuous softening vs ordering of local distortions
Computational tools like ALAMODE, SSCHA, and TDEP enable quantitative prediction of anharmonic properties from first principles

Practical Implications

Understanding anharmonic phonons is essential for:

Predicting structural phase transitions in functional materials (ferroelectrics, thermoelectrics)
Calculating thermal conductivity and designing thermal management materials
Understanding negative thermal expansion and engineering low-expansion materials
Interpreting temperature-dependent spectroscopy experiments
Designing materials with tailored temperature-dependent properties

12. Exercises

Exercise 1: Conceptual Understanding

Question: Explain why a purely harmonic crystal cannot exhibit thermal expansion. What anharmonic term (cubic or quartic) is responsible for thermal expansion?

Hint

Consider the symmetry of the potential energy and how the average position changes with temperature.

Exercise 2: SCP Self-Consistency

Problem: For a 1D harmonic oscillator with potential $V = \frac{1}{2}k_2 x^2 + \frac{1}{4}k_4 x^4$, derive the SCP equation for the effective force constant $k_{\text{eff}}(T)$.

Show that:

$$k_{\text{eff}}(T) = k_2 + \frac{3k_4}{2}\langle x^2 \rangle_T$$

Solution Outline

1. The effective potential is $V_{\text{eff}} = \frac{1}{2}k_{\text{eff}} x^2$

2. For a harmonic oscillator: $\langle x^2 \rangle = \frac{\hbar}{2m\omega}(2n+1)$ where $n$ is the Bose-Einstein distribution

3. Take the second derivative of $V$ and average: $\langle \partial^2 V/\partial x^2 \rangle = k_2 + 3k_4\langle x^2\rangle$

4. Self-consistency: solve $k_{\text{eff}} = k_2 + \frac{3k_4}{2m\omega_{\text{eff}}}(2n+1)$ iteratively

Exercise 3: Landau Theory Analysis

Problem: Consider a Landau free energy:

$$F = \frac{1}{2}a(T-T_c)\eta^2 + \frac{1}{4}b\eta^4$$

with $a = 0.5$ eV/Å² K, $T_c = 400$ K, $b = 2.0$ eV/Å⁴.

(a) Calculate the equilibrium order parameter $\eta_{\text{eq}}(T)$ for $T = 300$ K.

(b) Compute the free energy barrier for the phase transition.

(c) Determine the order parameter susceptibility $\chi = \partial^2 F/\partial \eta^2$ at $T = 410$ K.

Solution

(a) Minimize $F$: $\partial F/\partial \eta = a(T-T_c)\eta + b\eta^3 = 0$

For $T < T_c$: $\eta_{\text{eq}} = \pm\sqrt{-a(T-T_c)/b} = \pm\sqrt{0.5 \times 100 / 2.0} = \pm 5.0$ Å

(b) Barrier is zero for second-order transition (b > 0)

(c) $\chi^{-1} = a(T-T_c) = 0.5 \times 10 = 5.0$ eV/Å², so $\chi = 0.2$ Å²/eV

Exercise 4: Negative Thermal Expansion

Problem: A material has three phonon modes with Grüneisen parameters:

Mode 1 (acoustic): $\gamma_1 = +2.0$, $C_1 = 3k_B$
Mode 2 (acoustic): $\gamma_2 = +1.5$, $C_2 = 3k_B$
Mode 3 (optical): $\gamma_3 = -5.0$, $C_3 = 3k_B$

The bulk modulus is $B = 100$ GPa and volume is $V = 50$ Ų.

Calculate the volumetric thermal expansion coefficient $\alpha_V$.