Chapter 4: Topological Phonons - Advanced Phonon Physics

This chapter explores the emerging field of topological phononics, covering Berry phase formalism for phonons, topological invariants (Chern number, Z₂), Weyl phonons, topological edge/surface states, and experimental signatures in real materials including α-quartz and Na₃Bi.

Learning Objectives

By completing this chapter, you will be able to:

✅ Calculate Berry phase and Berry curvature for phonon bands
✅ Determine topological invariants (Chern number, Z₂ index) from phonon band structures
✅ Identify Weyl phonons and phonon arc surface states
✅ Understand the role of symmetry breaking (time-reversal, inversion) in phonon topology
✅ Predict topological edge/surface states in phononic systems
✅ Analyze phonon Hall effect and its connection to Berry curvature
✅ Apply topological concepts to 2D materials (graphene, TMDs)
✅ Implement Python code for Berry phase calculations and topological classification

4.1 Introduction to Topological Phonons

The concept of topology in condensed matter physics has revolutionized our understanding of electronic systems, from the quantum Hall effect to topological insulators. In recent years, these topological concepts have been extended to bosonic systems, particularly phonons—the quanta of lattice vibrations.

What Makes Phonons Topological?

Topological phonons are phonon bands with nontrivial topological invariants that are robust against perturbations preserving certain symmetries. Unlike electronic systems, phonons are bosonic and lack fundamental time-reversal symmetry breaking (no magnetic field for phonons), making their topological properties unique.

Key Differences: Electronic vs. Phononic Topology

Property	Electronic Systems	Phononic Systems
Particle Statistics	Fermionic	Bosonic
Time-Reversal Symmetry	Can be broken by magnetic field	No fundamental breaking mechanism
Inversion Symmetry	Broken by specific crystal structures	Same mechanism
Edge States	Conduct electricity	Conduct heat
Key Observable	Hall conductivity	Thermal Hall conductivity

Historical Context

timeline title Development of Topological Phononics 2015 : First theoretical prediction of Weyl phonons 2016 : Discovery of topological phonons in α-quartz 2017 : Experimental observation of phonon surface arcs 2018 : Phonon Hall effect measurements 2019 : Topological phonons in 2D materials 2020 : Phononic topological insulators 2021-2025 : Phonon engineering and applications

4.2 Berry Phase and Berry Curvature for Phonons

4.2.1 Berry Phase Formalism

The Berry phase is a geometric phase acquired by a quantum state when adiabatically transported around a closed loop in parameter space. For phonon eigenstates $|n,\mathbf{k}\rangle$ (band $n$, wavevector $\mathbf{k}$), the Berry connection is:

$$ \mathbf{A}_n(\mathbf{k}) = i \langle n,\mathbf{k} | \nabla_\mathbf{k} | n,\mathbf{k} \rangle $$

The Berry phase accumulated along a closed path $\mathcal{C}$ in $\mathbf{k}$-space is:

$$ \gamma_n[\mathcal{C}] = \oint_\mathcal{C} \mathbf{A}_n(\mathbf{k}) \cdot d\mathbf{k} $$

4.2.2 Berry Curvature

The Berry curvature is the curl of the Berry connection, analogous to magnetic field in real space:

$$ \mathbf{\Omega}_n(\mathbf{k}) = \nabla_\mathbf{k} \times \mathbf{A}_n(\mathbf{k}) $$

For a specific component in 3D:

$$ \Omega_n^z(\mathbf{k}) = \frac{\partial A_n^y}{\partial k_x} - \frac{\partial A_n^x}{\partial k_y} $$

Using the Kubo formula, the Berry curvature can be computed from phonon eigenstates:

$$ \Omega_n^{xy}(\mathbf{k}) = -2\text{Im}\sum_{m\neq n} \frac{\langle n,\mathbf{k}|\hat{v}_x|m,\mathbf{k}\rangle \langle m,\mathbf{k}|\hat{v}_y|n,\mathbf{k}\rangle}{(\omega_{m}(\mathbf{k}) - \omega_n(\mathbf{k}))^2} $$

where $\hat{v}_\alpha = \partial H(\mathbf{k})/\partial k_\alpha$ is the velocity operator and $\omega_n(\mathbf{k})$ are phonon frequencies.

4.2.3 Practical Calculation of Berry Curvature

💻 Code Example 1: Berry Curvature Calculation for Phonon Bands

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0
# - matplotlib>=3.7.0
# - scipy>=1.11.0

"""
Example: Berry Curvature Calculation for 2D Phonon Model
Purpose: Demonstrate numerical computation of Berry curvature
Target: Advanced graduate students
Execution time: 5-10 seconds
"""

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

def phonon_hamiltonian_2d(kx, ky, t1=1.0, t2=0.5, delta=0.3):
    """
    Construct 2D phonon tight-binding Hamiltonian

    Parameters:
    -----------
    kx, ky : float
        Wavevector components
    t1, t2 : float
        Hopping parameters (nearest, next-nearest neighbor)
    delta : float
        On-site energy difference (breaks inversion symmetry)

    Returns:
    --------
    H : ndarray (2x2)
        Phonon Hamiltonian matrix
    """
    h0 = delta
    hx = t1 * np.cos(kx) + t2 * np.cos(kx + ky)
    hy = t1 * np.sin(kx) + t2 * np.sin(kx + ky)
    hz = t1 * np.cos(ky) + t2 * np.cos(kx - ky)

    H = np.array([
        [h0 + hz, hx - 1j*hy],
        [hx + 1j*hy, -h0 - hz]
    ])
    return H

def calculate_berry_curvature(kx, ky, dk=1e-4, band_index=0):
    """
    Calculate Berry curvature using finite difference method

    Ω_z(k) = ∂A_y/∂kx - ∂A_x/∂ky

    Parameters:
    -----------
    kx, ky : float
        Wavevector point
    dk : float
        Finite difference step size
    band_index : int
        Band index (0 for lower band, 1 for upper band)

    Returns:
    --------
    omega_z : float
        Berry curvature component
    """
    # Get eigenstates at four neighboring points
    H_center = phonon_hamiltonian_2d(kx, ky)
    H_right = phonon_hamiltonian_2d(kx + dk, ky)
    H_left = phonon_hamiltonian_2d(kx - dk, ky)
    H_up = phonon_hamiltonian_2d(kx, ky + dk)
    H_down = phonon_hamiltonian_2d(kx, ky - dk)

    _, v_center = eigh(H_center)
    _, v_right = eigh(H_right)
    _, v_left = eigh(H_left)
    _, v_up = eigh(H_up)
    _, v_down = eigh(H_down)

    # Extract band eigenvector (ensure gauge consistency)
    u_center = v_center[:, band_index]
    u_right = v_right[:, band_index]
    u_left = v_left[:, band_index]
    u_up = v_up[:, band_index]
    u_down = v_down[:, band_index]

    # Fix gauge: maximize overlap with center state
    for u in [u_right, u_left, u_up, u_down]:
        if np.vdot(u, u_center).real < 0:
            u *= -1

    # Berry connection via finite difference
    # A_x = i⟨u|∂u/∂kx⟩
    Ax = np.imag(np.vdot(u_center, (u_right - u_left)/(2*dk)))
    Ay = np.imag(np.vdot(u_center, (u_up - u_down)/(2*dk)))

    # Berry curvature: Ω = ∂A_y/∂kx - ∂A_x/∂ky
    Ax_right = np.imag(np.vdot(u_right, (v_right[:, band_index] - u_center)/dk))
    Ax_left = np.imag(np.vdot(u_left, (u_center - v_left[:, band_index])/dk))

    Ay_up = np.imag(np.vdot(u_up, (v_up[:, band_index] - u_center)/dk))
    Ay_down = np.imag(np.vdot(u_down, (u_center - v_down[:, band_index])/dk))

    dAy_dkx = (Ay_up - Ay_down) / (2*dk)
    dAx_dky = (Ax_right - Ax_left) / (2*dk)

    omega_z = dAy_dkx - dAx_dky

    return omega_z

# Calculate Berry curvature over Brillouin zone
nk = 50
kx_array = np.linspace(-np.pi, np.pi, nk)
ky_array = np.linspace(-np.pi, np.pi, nk)

berry_curvature = np.zeros((nk, nk))

print("Calculating Berry curvature over 2D Brillouin zone...")
for i, kx in enumerate(kx_array):
    for j, ky in enumerate(ky_array):
        berry_curvature[j, i] = calculate_berry_curvature(kx, ky, band_index=0)
    if i % 10 == 0:
        print(f"Progress: {i/nk*100:.0f}%")

# Visualization
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Berry curvature heatmap
im = axes[0].contourf(kx_array/np.pi, ky_array/np.pi, berry_curvature,
                      levels=20, cmap='RdBu_r')
axes[0].set_xlabel(r'$k_x/\pi$', fontsize=12)
axes[0].set_ylabel(r'$k_y/\pi$', fontsize=12)
axes[0].set_title('Berry Curvature Distribution', fontsize=13, fontweight='bold')
axes[0].set_aspect('equal')
cbar = plt.colorbar(im, ax=axes[0])
cbar.set_label(r'$\Omega_z(\mathbf{k})$', fontsize=12)

# Line cut through high-symmetry path
# Path: Γ(0,0) → X(π,0) → M(π,π) → Γ(0,0)
npts = 100
path_k = []
path_omega = []

# Γ → X
for t in np.linspace(0, 1, npts):
    kx, ky = t*np.pi, 0
    path_k.append(t)
    path_omega.append(calculate_berry_curvature(kx, ky, band_index=0))

# X → M
for t in np.linspace(0, 1, npts):
    kx, ky = np.pi, t*np.pi
    path_k.append(1 + t)
    path_omega.append(calculate_berry_curvature(kx, ky, band_index=0))

# M → Γ
for t in np.linspace(0, 1, npts):
    kx, ky = (1-t)*np.pi, (1-t)*np.pi
    path_k.append(2 + t)
    path_omega.append(calculate_berry_curvature(kx, ky, band_index=0))

axes[1].plot(path_k, path_omega, 'b-', linewidth=2)
axes[1].axhline(y=0, color='k', linestyle='--', alpha=0.3)
axes[1].axvline(x=1, color='gray', linestyle=':', alpha=0.5, label='X')
axes[1].axvline(x=2, color='gray', linestyle=':', alpha=0.5, label='M')
axes[1].set_xlabel('High-symmetry path', fontsize=12)
axes[1].set_ylabel(r'$\Omega_z(\mathbf{k})$', fontsize=12)
axes[1].set_title('Berry Curvature along Γ-X-M-Γ', fontsize=13, fontweight='bold')
axes[1].set_xticks([0, 1, 2, 3])
axes[1].set_xticklabels(['Γ', 'X', 'M', 'Γ'])
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.savefig('berry_curvature_2d_phonon.png', dpi=300, bbox_inches='tight')
plt.show()

# Calculate total Berry phase (Chern number)
chern_number = np.sum(berry_curvature) * (2*np.pi/nk)**2 / (2*np.pi)
print(f"\nCalculated Chern number: {chern_number:.4f}")
print(f"Expected (topological): integer value (0, ±1, ±2, ...)")

Key Points from the Code:

Gauge Fixing: Ensuring phase consistency of eigenvectors is critical for accurate Berry curvature calculation
Finite Difference: Numerical differentiation requires small step size (dk ~ 10⁻⁴) for accuracy
Chern Number: Integral of Berry curvature over the Brillouin zone gives a topological invariant

4.3 Topological Invariants

4.3.1 Chern Number

The Chern number $C_n$ for band $n$ is the integral of Berry curvature over the entire Brillouin zone:

$$ C_n = \frac{1}{2\pi} \int_{\text{BZ}} \mathbf{\Omega}_n(\mathbf{k}) \cdot d^2\mathbf{k} $$

For 2D systems, this is always an integer and characterizes the topological class of the band. Systems with nonzero Chern number exhibit:

Quantized Hall conductivity: $\sigma_{xy} = C \cdot e^2/h$ (for electrons)
Phonon Hall effect: Transverse heat current in response to temperature gradient
Protected edge states: Number of chiral edge modes = |C|

4.3.2 Z₂ Topological Invariant

For systems with time-reversal symmetry (TRS), Chern number is always zero. However, a different topological invariant exists: the Z₂ index (ν = 0 or 1).

The Z₂ invariant can be calculated using the parity criterion at time-reversal invariant momenta (TRIM) points:

$$ (-1)^\nu = \prod_{i=1}^{N_{\text{TRIM}}} \prod_{m=1}^{N_{\text{occ}}} \xi_{2m}(\Gamma_i) $$

where $\xi_{2m}(\Gamma_i) = \pm 1$ is the parity eigenvalue of occupied band $m$ at TRIM point $\Gamma_i$.

Time-Reversal Invariant Momenta (TRIM) in Different Dimensions

Dimension	Number of TRIM	TRIM Points (in units of reciprocal lattice vectors)
1D	2	k = 0, π
2D	4	(0,0), (π,0), (0,π), (π,π)
3D	8	(0,0,0), (π,0,0), ..., (π,π,π)

4.3.3 Mirror Chern Number

In systems with mirror symmetry $M$ (e.g., $M_z: z \to -z$), phonon bands can be classified by mirror eigenvalues $m = \pm i$. The mirror Chern number $C_M$ is defined as:

$$ C_M = \frac{C_+ - C_-}{2} $$

where $C_\pm$ are Chern numbers for bands with mirror eigenvalue $\pm i$ in the mirror-invariant plane.

4.4 Weyl Phonons and Phonon Arcs

4.4.1 Weyl Points in Phonon Spectra

Weyl points are point-like band degeneracies in 3D that act as monopoles or antimonopoles of Berry curvature. Near a Weyl point at $\mathbf{k}_W$, the phonon dispersion is linear:

$$ \omega(\mathbf{k}) \approx \omega_W + v_i (\mathbf{k} - \mathbf{k}_W)_i $$

where $v_i$ are group velocities. The Berry curvature near a Weyl point diverges as:

$$ \mathbf{\Omega}(\mathbf{k}) \approx \pm \frac{\mathbf{k} - \mathbf{k}_W}{|\mathbf{k} - \mathbf{k}_W|^3} $$

The sign (±) defines the chirality of the Weyl point.

4.4.2 Topological Charge and Protection

Each Weyl point carries a topological charge (chirality) $\chi = \pm 1$:

$$ \chi = \frac{1}{2\pi} \oint_S \mathbf{\Omega}(\mathbf{k}) \cdot d\mathbf{S} $$

where $S$ is a small sphere enclosing the Weyl point.

Nielsen-Ninomiya Theorem for Phonons

The total chirality in the Brillouin zone must sum to zero:

$$ \sum_{\text{all Weyl points}} \chi_i = 0 $$

This means Weyl points appear in pairs with opposite chirality. They can only be removed by annihilating with a partner of opposite chirality.

4.4.3 Phonon Surface Arcs

Weyl phonons exhibit characteristic surface arc states that connect the projection of Weyl points with opposite chirality on the surface Brillouin zone.

4.4.4 Example: Weyl Phonons in FeSi

💻 Code Example 2: Identifying Weyl Points in 3D Phonon Dispersion

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0
# - matplotlib>=3.7.0
# - scipy>=1.11.0

"""
Example: Weyl Point Detection in 3D Model Hamiltonian
Purpose: Demonstrate identification and characterization of Weyl points
Target: Advanced
Execution time: 10-15 seconds
"""

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

def weyl_hamiltonian_3d(kx, ky, kz, m=0.5):
    """
    3D Weyl Hamiltonian with two Weyl points

    H(k) = k_x σ_x + k_y σ_y + (k_z - m) σ_z

    Weyl points at: (0, 0, ±m)
    """
    sigma_x = np.array([[0, 1], [1, 0]])
    sigma_y = np.array([[0, -1j], [1j, 0]])
    sigma_z = np.array([[1, 0], [0, -1]])

    H = kx * sigma_x + ky * sigma_y + (kz - m) * sigma_z
    return H

def find_weyl_points(m=0.5, search_range=2.0, grid_points=20):
    """
    Search for Weyl points (band crossings) in 3D BZ

    Returns:
    --------
    weyl_points : list of tuples
        List of (kx, ky, kz, chirality) for each Weyl point
    """
    kx_array = np.linspace(-search_range, search_range, grid_points)
    ky_array = np.linspace(-search_range, search_range, grid_points)
    kz_array = np.linspace(-search_range, search_range, grid_points)

    weyl_points = []
    threshold = 0.01  # Energy gap threshold

    for kx in kx_array:
        for ky in ky_array:
            for kz in kz_array:
                H = weyl_hamiltonian_3d(kx, ky, kz, m)
                energies, _ = eigh(H)
                gap = abs(energies[1] - energies[0])

                if gap < threshold:
                    # Calculate chirality from Berry curvature flux
                    chirality = np.sign(kz - m) if abs(kz - m) > 0.1 else 0
                    weyl_points.append((kx, ky, kz, chirality))

    return weyl_points

def calculate_berry_flux_sphere(kx0, ky0, kz0, radius=0.1, npts=10):
    """
    Calculate Berry flux through sphere around point (kx0, ky0, kz0)
    This gives the topological charge (chirality)
    """
    total_flux = 0
    theta_array = np.linspace(0, np.pi, npts)
    phi_array = np.linspace(0, 2*np.pi, npts)

    for theta in theta_array:
        for phi in phi_array:
            # Point on sphere
            kx = kx0 + radius * np.sin(theta) * np.cos(phi)
            ky = ky0 + radius * np.sin(theta) * np.sin(phi)
            kz = kz0 + radius * np.cos(theta)

            H = weyl_hamiltonian_3d(kx, ky, kz)
            energies, eigvecs = eigh(H)

            # Berry curvature contribution (simplified)
            # Full calculation requires derivatives
            u = eigvecs[:, 0]  # Lower band
            # Approximation of solid angle contribution
            total_flux += np.sin(theta)

    chirality = np.sign(total_flux)
    return chirality

# Find Weyl points
print("Searching for Weyl points...")
weyl_points = find_weyl_points(m=0.5)
print(f"Found {len(weyl_points)} Weyl point candidates")

# Analytical Weyl point locations
weyl_analytic = [(0, 0, 0.5, +1), (0, 0, -0.5, -1)]

# Calculate band structure along path through Weyl points
nk = 100
kz_path = np.linspace(-1.5, 1.5, nk)
energies_path = np.zeros((2, nk))

for i, kz in enumerate(kz_path):
    H = weyl_hamiltonian_3d(0, 0, kz, m=0.5)
    energies_path[:, i], _ = eigh(H)

# Visualization
fig = plt.figure(figsize=(15, 5))

# 3D Brillouin zone with Weyl points
ax1 = fig.add_subplot(131, projection='3d')
for wp in weyl_analytic:
    kx, ky, kz, chi = wp
    color = 'red' if chi > 0 else 'blue'
    marker = '^' if chi > 0 else 'v'
    ax1.scatter(kx, ky, kz, c=color, s=200, marker=marker,
               edgecolors='black', linewidths=2,
               label=f'χ={chi:+d}' if wp == weyl_analytic[0] or wp == weyl_analytic[1] else '')

ax1.set_xlabel(r'$k_x$', fontsize=11)
ax1.set_ylabel(r'$k_y$', fontsize=11)
ax1.set_zlabel(r'$k_z$', fontsize=11)
ax1.set_title('Weyl Points in 3D BZ', fontsize=13, fontweight='bold')
ax1.legend()

# Band structure showing Weyl crossings
ax2 = fig.add_subplot(132)
ax2.plot(kz_path, energies_path[0, :], 'b-', linewidth=2, label='Lower band')
ax2.plot(kz_path, energies_path[1, :], 'r-', linewidth=2, label='Upper band')
ax2.axvline(x=0.5, color='red', linestyle='--', alpha=0.5, label='Weyl point (χ=+1)')
ax2.axvline(x=-0.5, color='blue', linestyle='--', alpha=0.5, label='Weyl point (χ=-1)')
ax2.set_xlabel(r'$k_z$', fontsize=12)
ax2.set_ylabel(r'Energy (arb. units)', fontsize=12)
ax2.set_title('Band Structure along $k_z$ axis', fontsize=13, fontweight='bold')
ax2.grid(alpha=0.3)
ax2.legend()

# Energy contour near Weyl point
ax3 = fig.add_subplot(133)
kx_grid = np.linspace(-0.5, 0.5, 50)
kz_grid = np.linspace(0, 1.0, 50)
KX, KZ = np.meshgrid(kx_grid, kz_grid)
energy_gap = np.zeros_like(KX)

for i in range(len(kx_grid)):
    for j in range(len(kz_grid)):
        H = weyl_hamiltonian_3d(KX[j,i], 0, KZ[j,i], m=0.5)
        energies, _ = eigh(H)
        energy_gap[j, i] = abs(energies[1] - energies[0])

contour = ax3.contourf(KX, KZ, energy_gap, levels=20, cmap='viridis')
ax3.plot(0, 0.5, 'r*', markersize=15, label='Weyl point')
ax3.set_xlabel(r'$k_x$', fontsize=12)
ax3.set_ylabel(r'$k_z$', fontsize=12)
ax3.set_title('Energy Gap near Weyl Point', fontsize=13, fontweight='bold')
cbar = plt.colorbar(contour, ax=ax3)
cbar.set_label('Band gap', fontsize=11)
ax3.legend()

plt.tight_layout()
plt.savefig('weyl_phonons_3d.png', dpi=300, bbox_inches='tight')
plt.show()

print("\nWeyl Point Analysis:")
print(f"Analytical Weyl points: {weyl_analytic}")
print(f"Total chirality (should be 0): {sum([wp[3] for wp in weyl_analytic])}")

4.5 Topological Edge and Surface States

4.5.1 Bulk-Boundary Correspondence

The bulk-boundary correspondence is a fundamental principle in topological physics: a nontrivial bulk topological invariant guarantees the existence of protected boundary states.

Bulk-Boundary Correspondence for Different Topological Classes

Bulk Invariant	System Dimensionality	Boundary States	Protection Mechanism
Chern number C	2D	C chiral edge modes	Topology
Z₂ index ν=1	2D	Kramers pairs at edge	Time-reversal symmetry
Weyl points	3D	Surface Fermi arcs	Chirality conservation
Mirror Chern C_M	3D	Hinge states	Mirror symmetry

4.5.2 Computing Edge States

Edge states can be computed using a ribbon geometry: finite in one direction (e.g., $y$-direction) with edges, periodic in other directions.

💻 Code Example 3: Edge States in 2D Topological Phononic System

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0
# - matplotlib>=3.7.0
# - scipy>=1.11.0

"""
Example: Edge State Calculation for 2D Topological Phonon Model
Purpose: Demonstrate edge state emergence from nontrivial topology
Target: Advanced
Execution time: 15-20 seconds
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh
from scipy.sparse import diags, lil_matrix
from scipy.sparse.linalg import eigsh

def construct_ribbon_hamiltonian(kx, Ny=40, t1=1.0, t2=0.5, delta=0.3):
    """
    Construct Hamiltonian for ribbon geometry

    Parameters:
    -----------
    kx : float
        Wavevector along infinite direction (x)
    Ny : int
        Number of unit cells in finite direction (y)
    t1, t2 : float
        Hopping parameters
    delta : float
        On-site energy difference

    Returns:
    --------
    H : sparse matrix (2*Ny x 2*Ny)
        Ribbon Hamiltonian
    """
    dim = 2 * Ny  # 2 orbitals per unit cell
    H = lil_matrix((dim, dim), dtype=complex)

    for j in range(Ny):
        # On-site and intra-cell hopping
        idx_A = 2*j
        idx_B = 2*j + 1

        # On-site energies
        H[idx_A, idx_A] = delta
        H[idx_B, idx_B] = -delta

        # Intra-cell hopping
        hx = t1 * np.exp(1j * kx)
        H[idx_A, idx_B] = hx
        H[idx_B, idx_A] = np.conj(hx)

        # Inter-cell hopping (along y-direction)
        if j < Ny - 1:
            idx_A_next = 2*(j+1)
            idx_B_next = 2*(j+1) + 1

            hy = t2
            H[idx_A, idx_A_next] = hy
            H[idx_A_next, idx_A] = hy
            H[idx_B, idx_B_next] = hy
            H[idx_B_next, idx_B] = hy

    return H.tocsr()

def calculate_band_structure_ribbon(kx_array, Ny=40):
    """
    Calculate band structure for ribbon geometry

    Returns edge state dispersion
    """
    nk = len(kx_array)
    dim = 2 * Ny
    energies = np.zeros((dim, nk))
    eigenvectors = np.zeros((dim, dim, nk), dtype=complex)

    for i, kx in enumerate(kx_array):
        H = construct_ribbon_hamiltonian(kx, Ny)
        # Use sparse eigensolver for large systems
        if dim > 100:
            # Calculate only states near zero energy
            evals, evecs = eigsh(H, k=min(20, dim-2), sigma=0, which='LM')
            idx = np.argsort(evals)
            energies[:len(evals), i] = evals[idx]
            eigenvectors[:, :len(evals), i] = evecs[:, idx]
        else:
            evals, evecs = eigh(H.toarray())
            energies[:, i] = evals
            eigenvectors[:, :, i] = evecs

    return energies, eigenvectors

def calculate_edge_weight(eigenvector, Ny):
    """
    Calculate weight of state on top/bottom edges

    Returns:
    --------
    weight_top, weight_bottom : float
        Probability density on top/bottom 3 unit cells
    """
    dim = 2 * Ny
    prob = np.abs(eigenvector)**2

    # Top edge (first 3 unit cells)
    weight_top = np.sum(prob[:6])

    # Bottom edge (last 3 unit cells)
    weight_bottom = np.sum(prob[-6:])

    return weight_top, weight_bottom

# Calculate ribbon band structure
print("Calculating ribbon band structure...")
Ny = 40
nk = 100
kx_array = np.linspace(-np.pi, np.pi, nk)

energies, eigenvectors = calculate_band_structure_ribbon(kx_array, Ny)

# Identify edge states based on energy and spatial localization
edge_state_threshold = 0.3  # Energy threshold
localization_threshold = 0.3  # Edge weight threshold

edge_indices = []
for i in range(energies.shape[0]):
    for j in range(nk):
        if abs(energies[i, j]) < edge_state_threshold:
            evec = eigenvectors[:, i, j]
            w_top, w_bot = calculate_edge_weight(evec, Ny)
            if w_top > localization_threshold or w_bot > localization_threshold:
                edge_indices.append((i, j))

# Visualization
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# Band structure with edge states highlighted
for i in range(energies.shape[0]):
    axes[0].plot(kx_array/np.pi, energies[i, :], 'k-', linewidth=0.5, alpha=0.3)

# Highlight edge states
if edge_indices:
    for idx_band, idx_k in edge_indices:
        axes[0].plot(kx_array[idx_k]/np.pi, energies[idx_band, idx_k],
                    'r.', markersize=3)

axes[0].set_xlabel(r'$k_x / \pi$', fontsize=12)
axes[0].set_ylabel(r'Energy (arb. units)', fontsize=12)
axes[0].set_title('Ribbon Band Structure\n(Red = Edge States)',
                 fontsize=13, fontweight='bold')
axes[0].set_ylim([-3, 3])
axes[0].axhline(y=0, color='gray', linestyle='--', alpha=0.5)
axes[0].grid(alpha=0.3)

# Spatial profile of edge state
kx_edge = 0  # Edge state at kx=0
idx_kx = np.argmin(np.abs(kx_array - kx_edge))
energies_at_kx = energies[:, idx_kx]
idx_edge_state = np.argmin(np.abs(energies_at_kx))  # State closest to zero

evec_edge = eigenvectors[:, idx_edge_state, idx_kx]
spatial_profile = np.abs(evec_edge)**2

# Reshape to show A and B sublattices
y_positions = np.arange(Ny)
profile_A = spatial_profile[0::2]  # A sublattice
profile_B = spatial_profile[1::2]  # B sublattice

axes[1].plot(y_positions, profile_A, 'b-o', markersize=4,
            label='A sublattice', linewidth=2)
axes[1].plot(y_positions, profile_B, 'r-s', markersize=4,
            label='B sublattice', linewidth=2)
axes[1].set_xlabel('Unit cell index (y-direction)', fontsize=12)
axes[1].set_ylabel('Probability density', fontsize=12)
axes[1].set_title(f'Edge State Localization\n(E = {energies_at_kx[idx_edge_state]:.3f})',
                 fontsize=13, fontweight='bold')
axes[1].legend()
axes[1].grid(alpha=0.3)
axes[1].axhline(y=0, color='k', linestyle='-', linewidth=0.5)

# Edge state dispersion extraction
edge_energies_top = []
edge_energies_bot = []
kx_edge_points = []

for j, kx in enumerate(kx_array):
    for i in range(energies.shape[0]):
        if abs(energies[i, j]) < edge_state_threshold:
            evec = eigenvectors[:, i, j]
            w_top, w_bot = calculate_edge_weight(evec, Ny)

            if w_top > localization_threshold:
                edge_energies_top.append(energies[i, j])
                kx_edge_points.append(kx)
            elif w_bot > localization_threshold:
                edge_energies_bot.append(energies[i, j])

if edge_energies_top:
    axes[2].scatter(np.array(kx_edge_points)/np.pi, edge_energies_top,
                   c='blue', s=10, label='Top edge', alpha=0.6)
if edge_energies_bot:
    axes[2].scatter(np.array(kx_edge_points)/np.pi, edge_energies_bot,
                   c='red', s=10, label='Bottom edge', alpha=0.6)

axes[2].set_xlabel(r'$k_x / \pi$', fontsize=12)
axes[2].set_ylabel(r'Energy (arb. units)', fontsize=12)
axes[2].set_title('Edge State Dispersion', fontsize=13, fontweight='bold')
axes[2].legend()
axes[2].grid(alpha=0.3)
axes[2].axhline(y=0, color='k', linestyle='--', alpha=0.3)

plt.tight_layout()
plt.savefig('topological_edge_states.png', dpi=300, bbox_inches='tight')
plt.show()

print("\nEdge State Analysis:")
print(f"Number of edge state points found: {len(edge_indices)}")
print(f"Edge state energy range: ±{edge_state_threshold}")
print("Edge states connect bulk bands with opposite Chern numbers")

4.6 Symmetry Breaking and Topological Phase Transitions

4.6.1 Role of Time-Reversal Symmetry

Unlike electronic systems where magnetic fields break time-reversal symmetry (TRS), phonons lack a direct coupling mechanism. However, effective TRS breaking can occur through:

Coriolis force: Rotating systems (phonon angular momentum coupling)
Gyroscopic coupling: Mechanical metamaterials with rotating components
Effective magnetic field: Phonon Hall effect in temperature gradients

4.6.2 Inversion Symmetry Breaking

Inversion symmetry breaking is more natural in phononic systems and occurs in:

Non-centrosymmetric crystals: e.g., wurtzite structures, α-quartz
Heterostructures: Asymmetric layering in 2D materials
Surface terminations: Breaking of bulk inversion at surfaces

4.6.3 Topological Phase Transitions

Topological phase transitions occur when a topological invariant changes, typically accompanied by gap closing. For phonons, this can be driven by:

Pressure: Modifying lattice constants and force constants
Temperature: Thermal expansion and anharmonicity
Composition: Doping or alloying
Strain: Engineering band structure through mechanical deformation

4.7 Phonon Hall Effect

4.7.1 Thermal Hall Conductivity

The phonon Hall effect refers to the generation of a transverse heat current $\mathbf{j}_Q$ in response to a temperature gradient $\nabla T$ in the presence of an external perturbation (real or effective magnetic field).

$$ j_Q^y = \kappa_{xy} \frac{\partial T}{\partial x} $$

The thermal Hall conductivity $\kappa_{xy}$ is related to Berry curvature:

$$ \kappa_{xy} = \frac{k_B^2 T}{\hbar} \sum_n \int_{\text{BZ}} \frac{d^2k}{(2\pi)^2} c_2\left(\frac{\hbar\omega_n}{k_B T}\right) \Omega_n^z(\mathbf{k}) $$

where $c_2(x) = (1+x)[\text{Li}_2(-e^x) + \frac{x^2}{2}]$ is a thermal factor, and $\text{Li}_2$ is the polylogarithm.

4.7.2 Experimental Observations

Phonon Hall effects have been observed in:

Paramagnetic insulators: Tb₃Ga₅O₁₂ (phonon-spin coupling)
Strontium titanate (SrTiO₃): Soft phonon modes
Topological materials: Effective Berry curvature contributions

💻 Code Example 4: Phonon Hall Conductivity Calculation

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0
# - matplotlib>=3.7.0
# - scipy>=1.11.0

"""
Example: Phonon Thermal Hall Conductivity from Berry Curvature
Purpose: Demonstrate calculation of Hall conductivity
Target: Advanced
Execution time: 5-8 seconds
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import spence  # Polylogarithm Li_2

def thermal_factor_c2(x):
    """
    Thermal factor for thermal Hall conductivity

    c_2(x) = (1 + x)[Li_2(-e^x) + x^2/2]

    where Li_2 is the dilogarithm (Spence function)
    """
    # Avoid overflow for large x
    if x > 10:
        return 0

    # Li_2(-e^x) = -Li_2(-exp(x))
    # Use spence: spence(x) = Li_2(1-x)
    # So Li_2(-e^x) = spence(1 + e^x)
    li2_term = -spence(1 + np.exp(x))

    c2 = (1 + x) * (li2_term + x**2/2)
    return c2

def phonon_thermal_hall_conductivity(berry_curvature, frequencies,
                                    temperature, kx_array, ky_array):
    """
    Calculate thermal Hall conductivity κ_xy

    Parameters:
    -----------
    berry_curvature : ndarray (nk_x, nk_y)
        Berry curvature Ω^z(k) in 2D BZ
    frequencies : ndarray (nk_x, nk_y)
        Phonon frequencies ω(k) [in THz]
    temperature : float
        Temperature [K]
    kx_array, ky_array : ndarray
        k-point grids

    Returns:
    --------
    kappa_xy : float
        Thermal Hall conductivity [W/K·m]
    """
    # Constants
    hbar = 1.054571817e-34  # J·s
    kB = 1.380649e-23  # J/K

    # Convert frequency to angular frequency [rad/s]
    omega = frequencies * 2 * np.pi * 1e12  # THz → rad/s

    # Dimensionless frequency
    x = hbar * omega / (kB * temperature)

    # Calculate thermal factor for each k-point
    c2_array = np.zeros_like(x)
    for i in range(x.shape[0]):
        for j in range(x.shape[1]):
            if x[i, j] < 50:  # Avoid numerical overflow
                c2_array[i, j] = thermal_factor_c2(x[i, j])

    # Integrand: c_2(x) * Ω(k)
    integrand = c2_array * berry_curvature

    # BZ integration (trapezoidal rule)
    dkx = (kx_array[-1] - kx_array[0]) / len(kx_array)
    dky = (ky_array[-1] - ky_array[0]) / len(ky_array)

    integral = np.sum(integrand) * dkx * dky / (2*np.pi)**2

    # Thermal Hall conductivity
    kappa_xy = (kB**2 * temperature / hbar) * integral

    return kappa_xy

# Example: Calculate for model 2D phonon system
# Use Berry curvature from previous calculation (Example 1)
nk = 50
kx_array = np.linspace(-np.pi, np.pi, nk)
ky_array = np.linspace(-np.pi, np.pi, nk)

# Create mock Berry curvature and frequency data
# (In real calculation, use actual phonon dispersion)
KX, KY = np.meshgrid(kx_array, ky_array)

# Model: Gaussian Berry curvature peak
berry_curvature = 2.0 * np.exp(-((KX-0.5)**2 + (KY-0.5)**2)/0.5)

# Model: Simple phonon dispersion [THz]
frequencies = 3.0 + 2.0 * (np.cos(KX) + np.cos(KY))
frequencies = np.abs(frequencies)  # Ensure positive

# Calculate thermal Hall conductivity vs temperature
T_array = np.linspace(10, 300, 50)  # Temperature range [K]
kappa_xy_array = []

print("Calculating thermal Hall conductivity vs temperature...")
for T in T_array:
    kappa = phonon_thermal_hall_conductivity(
        berry_curvature, frequencies, T, kx_array, ky_array
    )
    kappa_xy_array.append(kappa)

kappa_xy_array = np.array(kappa_xy_array)

# Visualization
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# Berry curvature distribution
im1 = axes[0].contourf(KX/np.pi, KY/np.pi, berry_curvature,
                       levels=20, cmap='RdBu_r')
axes[0].set_xlabel(r'$k_x/\pi$', fontsize=12)
axes[0].set_ylabel(r'$k_y/\pi$', fontsize=12)
axes[0].set_title('Berry Curvature $\Omega^z(\mathbf{k})$',
                 fontsize=13, fontweight='bold')
axes[0].set_aspect('equal')
cbar1 = plt.colorbar(im1, ax=axes[0])
cbar1.set_label('Berry curvature', fontsize=11)

# Phonon frequency
im2 = axes[1].contourf(KX/np.pi, KY/np.pi, frequencies,
                       levels=20, cmap='viridis')
axes[1].set_xlabel(r'$k_x/\pi$', fontsize=12)
axes[1].set_ylabel(r'$k_y/\pi$', fontsize=12)
axes[1].set_title('Phonon Frequency $\omega(\mathbf{k})$',
                 fontsize=13, fontweight='bold')
axes[1].set_aspect('equal')
cbar2 = plt.colorbar(im2, ax=axes[1])
cbar2.set_label('Frequency [THz]', fontsize=11)

# Thermal Hall conductivity vs temperature
axes[2].plot(T_array, kappa_xy_array * 1e3, 'b-', linewidth=2)
axes[2].set_xlabel('Temperature [K]', fontsize=12)
axes[2].set_ylabel(r'$\kappa_{xy}$ [mW/K·m]', fontsize=12)
axes[2].set_title('Phonon Thermal Hall Conductivity',
                 fontsize=13, fontweight='bold')
axes[2].grid(alpha=0.3)
axes[2].axhline(y=0, color='k', linestyle='--', alpha=0.3)

plt.tight_layout()
plt.savefig('phonon_hall_effect.png', dpi=300, bbox_inches='tight')
plt.show()

print(f"\nThermal Hall conductivity at T=100K: {kappa_xy_array[20]*1e3:.4f} mW/K·m")
print(f"Peak value: {np.max(kappa_xy_array)*1e3:.4f} mW/K·m at T={T_array[np.argmax(kappa_xy_array)]:.1f} K")

4.8 Topological Phonons in 2D Materials

4.8.1 Graphene Phonons

Graphene, despite having inversion symmetry, exhibits rich phononic topology when symmetry is broken. The acoustic phonon branches near the Dirac point show:

Dirac cones: Linear dispersion at K and K' points
Valley-contrasting physics: Opposite Berry curvature at K and K'
Edge states: Zigzag edge phonon modes with valley polarization

4.8.2 Transition Metal Dichalcogenides (TMDs)

TMDs (MoS₂, WSe₂, etc.) naturally break inversion symmetry in monolayer form, leading to:

Valley-selective phonons: Optical phonon modes at K/K' with Berry curvature
Phonon Hall effect: Measurable at room temperature
Coupled electron-phonon topology: Both electronic and phononic edge states

graph TB subgraph "TMD Monolayer (e.g., MoS₂)" A[Broken Inversion
Symmetry] B[Valley-Dependent
Berry Curvature] C[Chiral Phonon
Modes] D[Phonon Hall
Effect] end A --> B B --> C B --> D C --> D style A fill:#ff6b6b style B fill:#4ecdc4 style C fill:#ffe66d style D fill:#95e1d3

4.8.3 Phononic Crystals and Metamaterials

Phononic crystals are periodic structures engineered to control phonon propagation, analogous to photonic crystals for light. They can realize:

Topological band gaps: Frequency ranges with no bulk propagation
Protected waveguides: Topological edge states for robust phonon transport
Acoustic Chern insulators: Classical analogs of quantum Hall systems
Higher-order topology: Corner and hinge states in 2D/3D structures

💻 Code Example 5: Valley-Polarized Phonons in Honeycomb Lattice

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0
# - matplotlib>=3.7.0

"""
Example: Valley-Polarized Phonon Modes in Graphene-like Honeycomb Lattice
Purpose: Demonstrate valley-contrasting Berry curvature
Target: Advanced
Execution time: 8-12 seconds
"""

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

def honeycomb_phonon_hamiltonian(kx, ky, t=1.0, delta=0):
    """
    Tight-binding model for honeycomb lattice phonons

    Parameters:
    -----------
    kx, ky : float
        Wavevector components
    t : float
        Nearest-neighbor hopping
    delta : float
        Sublattice energy difference (breaks inversion symmetry)

    Returns:
    --------
    H : ndarray (2x2)
        Hamiltonian at k-point
    """
    # Nearest-neighbor vectors in honeycomb lattice
    a = 1.0  # Lattice constant
    delta1 = np.array([a, 0])
    delta2 = np.array([-a/2, a*np.sqrt(3)/2])
    delta3 = np.array([-a/2, -a*np.sqrt(3)/2])

    # Structure factor
    f_k = (np.exp(1j * (kx * delta1[0] + ky * delta1[1])) +
           np.exp(1j * (kx * delta2[0] + ky * delta2[1])) +
           np.exp(1j * (kx * delta3[0] + ky * delta3[1])))

    # Hamiltonian
    H = np.array([
        [delta, t * f_k],
        [t * np.conj(f_k), -delta]
    ])

    return H

def calculate_valley_chern_numbers(delta=0.1, nk=30):
    """
    Calculate Chern numbers around K and K' valleys

    Returns:
    --------
    C_K, C_Kprime : float
        Chern numbers for K and K' valleys
    """
    # K and K' points in honeycomb Brillouin zone
    K = np.array([4*np.pi/(3*np.sqrt(3)), 0])
    Kprime = np.array([-4*np.pi/(3*np.sqrt(3)), 0])

    # Small region around K point
    dk = 0.5
    kx_K = np.linspace(K[0]-dk, K[0]+dk, nk)
    ky_K = np.linspace(K[1]-dk, K[1]+dk, nk)

    berry_curv_K = np.zeros((nk, nk))

    for i, kx in enumerate(kx_K):
        for j, ky in enumerate(ky_K):
            # Simplified Berry curvature calculation
            H = honeycomb_phonon_hamiltonian(kx, ky, delta=delta)
            energies, eigvecs = eigh(H)

            # Berry curvature approximation from band structure
            # (Full calculation requires derivatives as in Example 1)
            # This is a simplified pedagogical version
            berry_curv_K[j, i] = np.imag(
                eigvecs[1, 0] * np.conj(eigvecs[0, 0])
            ) * 2

    C_K = np.sum(berry_curv_K) * (2*dk/nk)**2 / (2*np.pi)

    # Similarly for K' (opposite chirality due to time-reversal)
    C_Kprime = -C_K

    return C_K, C_Kprime

# Calculate band structure along high-symmetry path
# Path: Γ → M → K → Γ
npts = 100

# High-symmetry points
Gamma = np.array([0, 0])
M = np.array([np.pi/np.sqrt(3), 0])
K = np.array([4*np.pi/(3*np.sqrt(3)), 0])

path_k = []
path_labels = []
k_ticks = [0]

# Γ → M
for i, t in enumerate(np.linspace(0, 1, npts)):
    kx = Gamma[0] + t * (M[0] - Gamma[0])
    ky = Gamma[1] + t * (M[1] - Gamma[1])
    path_k.append([kx, ky, i/npts])
k_ticks.append(1)

# M → K
for i, t in enumerate(np.linspace(0, 1, npts)):
    kx = M[0] + t * (K[0] - M[0])
    ky = M[1] + t * (K[1] - M[1])
    path_k.append([kx, ky, 1 + i/npts])
k_ticks.append(2)

# K → Γ
for i, t in enumerate(np.linspace(0, 1, npts)):
    kx = K[0] + t * (Gamma[0] - K[0])
    ky = K[1] + t * (Gamma[1] - K[1])
    path_k.append([kx, ky, 2 + i/npts])
k_ticks.append(3)

path_k = np.array(path_k)

# Calculate energies along path
energies_pristine = []
energies_broken = []

for kpt in path_k:
    H_pristine = honeycomb_phonon_hamiltonian(kpt[0], kpt[1], delta=0)
    H_broken = honeycomb_phonon_hamiltonian(kpt[0], kpt[1], delta=0.3)

    E_pristine, _ = eigh(H_pristine)
    E_broken, _ = eigh(H_broken)

    energies_pristine.append(E_pristine)
    energies_broken.append(E_broken)

energies_pristine = np.array(energies_pristine)
energies_broken = np.array(energies_broken)

# Calculate valley Chern numbers
C_K, C_Kprime = calculate_valley_chern_numbers(delta=0.3)

# Visualization
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Band structure comparison
axes[0].plot(path_k[:, 2], energies_pristine[:, 0], 'b-',
            linewidth=2, label='Lower band (pristine)', alpha=0.7)
axes[0].plot(path_k[:, 2], energies_pristine[:, 1], 'r-',
            linewidth=2, label='Upper band (pristine)', alpha=0.7)
axes[0].plot(path_k[:, 2], energies_broken[:, 0], 'b--',
            linewidth=2, label='Lower band (δ=0.3)', alpha=0.7)
axes[0].plot(path_k[:, 2], energies_broken[:, 1], 'r--',
            linewidth=2, label='Upper band (δ=0.3)', alpha=0.7)

for tick in k_ticks:
    axes[0].axvline(x=tick, color='gray', linestyle=':', alpha=0.5)

axes[0].set_xlabel('Wave vector', fontsize=12)
axes[0].set_ylabel('Energy (arb. units)', fontsize=12)
axes[0].set_title('Honeycomb Phonon Band Structure\n(Inversion Symmetry Breaking)',
                 fontsize=13, fontweight='bold')
axes[0].set_xticks(k_ticks)
axes[0].set_xticklabels(['Γ', 'M', 'K', 'Γ'])
axes[0].legend(fontsize=9)
axes[0].grid(alpha=0.3)
axes[0].axhline(y=0, color='k', linestyle='-', linewidth=0.5)

# Valley Chern number illustration
valley_data = {
    'Valley': ['K', "K'"],
    'Chern Number': [C_K, C_Kprime],
    'Berry Curvature': ['+', '-']
}

axes[1].bar(valley_data['Valley'], valley_data['Chern Number'],
           color=['red', 'blue'], alpha=0.7, edgecolor='black', linewidth=2)
axes[1].axhline(y=0, color='black', linestyle='-', linewidth=1)
axes[1].set_ylabel('Valley Chern Number', fontsize=12)
axes[1].set_title('Valley-Contrasting Topology\n(δ = 0.3)',
                 fontsize=13, fontweight='bold')
axes[1].grid(axis='y', alpha=0.3)
axes[1].set_ylim([-1, 1])

# Add annotations
for i, (valley, C) in enumerate(zip(valley_data['Valley'], valley_data['Chern Number'])):
    axes[1].text(i, C + 0.05, f'{C:.3f}',
                ha='center', va='bottom', fontsize=11, fontweight='bold')

plt.tight_layout()
plt.savefig('valley_polarized_phonons.png', dpi=300, bbox_inches='tight')
plt.show()

print("\nValley Chern Number Analysis:")
print(f"C_K = {C_K:.4f}")
print(f"C_K' = {C_Kprime:.4f}")
print(f"Total Chern number = {C_K + C_Kprime:.4f} (should be 0 by symmetry)")
print("\nValley-polarized phonons enable:")
print("• Valley-selective excitation with circularly polarized light")
print("• Topological valley transport (valley Hall effect)")
print("• Protected phonon edge states")

4.9 Real Material Examples

4.9.1 α-Quartz (SiO₂)

α-Quartz was one of the first materials where topological phonons were predicted and observed. Key features:

Crystal structure: Trigonal (P3₁21 or P3₂21), non-centrosymmetric
Topological features: Weyl phonon points at ~25 THz
Experimental signature: Phonon surface arcs observed via inelastic X-ray scattering
Chirality: Left- and right-handed crystals have opposite Weyl point chiralities

Experimental Detection of Topological Phonons in α-Quartz

Technique	Observable	Key Finding
Inelastic X-ray Scattering	Phonon dispersion	Weyl points at ~25 THz
Surface phonon spectroscopy	Surface states	Phonon arcs connecting Weyl projections
Thermal transport	Hall conductivity κ_xy	Weak but measurable signal
First-principles (DFT)	Berry curvature	Confirms topological classification

4.9.2 Na₃Bi (Sodium Bismuthide)

Na₃Bi is a Dirac semimetal in its electronic structure and also hosts topological phonon modes:

Crystal structure: Hexagonal, space group P6₃/mmc
Electronic Dirac points: Protected by crystal symmetry
Phononic topology: Nodal lines and Weyl points in phonon dispersion
Coupled topologies: Electron-phonon interaction enhanced by both being topological

4.9.3 Other Material Candidates

FeSi: Cubic, topological phonon bands with large Berry curvature
CoSi: Chiral crystal structure, both electronic and phononic Weyl points
Half-Heusler compounds: RPtBi (R = rare earth), topological phonons relevant for thermoelectrics
2D materials: Monolayer MoS₂, WSe₂ (valley-polarized phonons)

💻 Code Example 6: Analyzing Phonon Topology from DFT Data

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0
# - matplotlib>=3.7.0
# - ase>=3.22.0 (Atomic Simulation Environment)

"""
Example: Loading and Analyzing Phonon Data from First-Principles Calculations
Purpose: Demonstrate realistic workflow for topological analysis
Target: Graduate researchers
Execution time: Variable (depends on data)
"""

import numpy as np
import matplotlib.pyplot as plt

def load_phonopy_band_structure(filename='band.yaml'):
    """
    Load phonon band structure from Phonopy output

    Returns:
    --------
    qpoints : ndarray
        q-point coordinates
    distances : ndarray
        Distance along path
    frequencies : ndarray (nqpts, nbands)
        Phonon frequencies [THz]
    eigenvectors : ndarray (complex)
        Phonon eigenvectors
    """
    # This is a simplified version
    # In practice, use: from phonopy import load

    # Simulated data for demonstration
    nqpts = 100
    nbands = 12  # Example: 4 atoms × 3 directions

    qpoints = np.random.rand(nqpts, 3)
    distances = np.linspace(0, 5, nqpts)
    frequencies = 5 + 10 * np.random.rand(nqpts, nbands)
    frequencies = np.sort(frequencies, axis=1)  # Sort bands

    eigenvectors = (np.random.rand(nqpts, nbands, nbands) +
                   1j * np.random.rand(nqpts, nbands, nbands))

    # Normalize eigenvectors
    for i in range(nqpts):
        for j in range(nbands):
            eigenvectors[i, :, j] /= np.linalg.norm(eigenvectors[i, :, j])

    return qpoints, distances, frequencies, eigenvectors

def identify_band_crossings(frequencies, threshold=0.1):
    """
    Identify potential Weyl points or nodal lines

    Parameters:
    -----------
    frequencies : ndarray (nqpts, nbands)
        Phonon frequencies
    threshold : float
        Energy difference threshold for crossing [THz]

    Returns:
    --------
    crossings : list of tuples
        (q_index, band1, band2, gap)
    """
    nqpts, nbands = frequencies.shape
    crossings = []

    for iq in range(nqpts):
        for ib1 in range(nbands - 1):
            for ib2 in range(ib1 + 1, nbands):
                gap = abs(frequencies[iq, ib2] - frequencies[iq, ib1])
                if gap < threshold:
                    crossings.append((iq, ib1, ib2, gap))

    return crossings

def calculate_wilson_loop(eigenvectors, k_plane='kz', nk1=20, nk2=20):
    """
    Calculate Wilson loop for Z2 invariant determination

    The Wilson loop is:
    W(k2) = ∏_{k1} ⟨u_n(k1, k2) | u_n(k1 + dk1, k2)⟩

    Returns:
    --------
    wilson_eigenvalues : ndarray (nk2,)
        Eigenvalues of Wilson loop operator
    """
    # Simplified calculation for demonstration
    # Real implementation requires careful gauge fixing

    wilson_eigenvalues = np.exp(2j * np.pi * np.random.rand(nk2))

    return wilson_eigenvalues

# Load phonon data (simulated for this example)
print("Loading phonon band structure data...")
qpoints, distances, frequencies, eigenvectors = load_phonopy_band_structure()

print(f"Loaded {len(qpoints)} q-points with {frequencies.shape[1]} phonon bands")

# Identify band crossings
crossings = identify_band_crossings(frequencies, threshold=0.5)
print(f"\nFound {len(crossings)} potential band crossing points")

if crossings:
    print("\nFirst 5 crossings:")
    for i, (iq, ib1, ib2, gap) in enumerate(crossings[:5]):
        print(f"  {i+1}. q-index {iq}: bands {ib1}-{ib2}, gap = {gap:.4f} THz")

# Visualization
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# Phonon band structure
for ib in range(frequencies.shape[1]):
    axes[0].plot(distances, frequencies[:, ib], 'b-', linewidth=1, alpha=0.7)

# Mark band crossings
if crossings:
    crossing_q = [distances[c[0]] for c in crossings]
    crossing_freq = [frequencies[c[0], c[1]] for c in crossings]
    axes[0].scatter(crossing_q, crossing_freq, c='red', s=50,
                   zorder=10, label='Potential Weyl/Nodal points')

axes[0].set_xlabel('q-point path', fontsize=12)
axes[0].set_ylabel('Frequency [THz]', fontsize=12)
axes[0].set_title('Phonon Band Structure\n(with crossing points)',
                 fontsize=13, fontweight='bold')
axes[0].legend()
axes[0].grid(alpha=0.3)
axes[0].set_ylim([0, frequencies.max() * 1.1])

# Phonon density of states (mock)
freq_bins = np.linspace(0, frequencies.max(), 100)
dos, _ = np.histogram(frequencies.flatten(), bins=freq_bins)

axes[1].fill_between(dos, freq_bins[:-1], alpha=0.5, color='steelblue')
axes[1].plot(dos, freq_bins[:-1], 'b-', linewidth=2)
axes[1].set_xlabel('Phonon DOS (states/THz)', fontsize=12)
axes[1].set_ylabel('Frequency [THz]', fontsize=12)
axes[1].set_title('Phonon Density of States', fontsize=13, fontweight='bold')
axes[1].grid(alpha=0.3)

# Wilson loop eigenvalue evolution (for Z2 calculation)
nk2 = 50
wilson_phases = []
for i in range(nk2):
    # Simplified: just show phase evolution
    wilson_phases.append(np.angle(np.exp(2j * np.pi * i / nk2)))

axes[2].plot(range(nk2), wilson_phases, 'g-', linewidth=2, marker='o',
            markersize=4)
axes[2].axhline(y=np.pi, color='r', linestyle='--', label='π (topological)')
axes[2].axhline(y=0, color='b', linestyle='--', label='0 (trivial)')
axes[2].set_xlabel('Wilson loop step', fontsize=12)
axes[2].set_ylabel('Phase [rad]', fontsize=12)
axes[2].set_title('Wilson Loop Evolution\n(Z₂ Invariant)',
                 fontsize=13, fontweight='bold')
axes[2].legend()
axes[2].grid(alpha=0.3)
axes[2].set_ylim([-np.pi - 0.5, np.pi + 0.5])

plt.tight_layout()
plt.savefig('dft_phonon_topology_analysis.png', dpi=300, bbox_inches='tight')
plt.show()

print("\n" + "="*60)
print("TOPOLOGICAL ANALYSIS SUMMARY")
print("="*60)
print(f"Material: [Load from DFT calculation]")
print(f"Number of atoms in primitive cell: {frequencies.shape[1]//3}")
print(f"Phonon frequency range: {frequencies.min():.2f} - {frequencies.max():.2f} THz")
print(f"Potential topological features: {len(crossings)} band crossings")
print("\nNext steps for full topological characterization:")
print("1. Calculate Berry curvature over full 3D BZ")
print("2. Compute Chern/Z2 invariants")
print("3. Analyze surface/edge states with slab geometry")
print("4. Check symmetry protection of topological features")
print("5. Predict experimental signatures (ARPES, IXS, thermal transport)")

4.10 Summary

This chapter introduced the emerging field of topological phononics, extending concepts from electronic topological materials to bosonic lattice vibrations. Key takeaways:

Essential Concepts

Berry Phase & Curvature: Geometric phase accumulated by phonon states, sourced by band structure geometry
Topological Invariants: Integer quantities (Chern number, Z₂ index) that classify phonon band topology
Weyl Phonons: Point degeneracies acting as Berry curvature monopoles, appearing in pairs with opposite chirality
Surface States: Topologically protected boundary modes (arcs, edge states) guaranteed by bulk-boundary correspondence
Symmetry Breaking: Inversion symmetry breaking enables nontrivial topology; TRS breaking less accessible for phonons
Phonon Hall Effect: Transverse heat current from Berry curvature, measurable in experiments
Material Realizations: α-quartz (Weyl phonons), TMDs (valley topology), metamaterials (engineered topology)

Connections to Broader Physics

graph LR A[Topological
Phonons] --> B[Thermal
Transport] A --> C[Phonon
Engineering] A --> D[Quantum
Materials] A --> E[Metamaterials] B --> F[Thermoelectric
Efficiency] C --> G[Thermal
Management] D --> H[Electron-Phonon
Coupling] E --> I[Acoustic
Devices] style A fill:#ff6b6b style B fill:#4ecdc4 style C fill:#ffe66d style D fill:#95e1d3 style E fill:#ffd93d

Future Directions

Higher-order topology: Corner and hinge phonon modes in 2D/3D
Non-Hermitian phonons: Dissipation-induced topology
Floquet phonons: Time-periodic driving to engineer topology
Phonon spintronics: Coupling phonon topology to spin degrees of freedom
Applications: Topological thermal diodes, robust waveguides, enhanced thermoelectrics

Exercises

Exercise 4.1: Berry Phase Calculation (Intermediate)

Consider a 1D phonon band with dispersion $\omega(k) = \omega_0 + 2t\cos(ka)$ and eigenvector $|u(k)\rangle = [\cos(\theta(k)/2), e^{i\phi(k)}\sin(\theta(k)/2)]^T$ where $\theta(k) = \pi/4$ and $\phi(k) = k$.

Tasks:

Calculate the Berry connection $A(k)$ as a function of $k$
Compute the Berry phase $\gamma$ for a loop from $k=0$ to $k=2\pi/a$
Discuss whether this Berry phase is gauge-invariant

Exercise 4.2: Chern Number from Band Structure (Advanced)

Using the 2D phonon model from Code Example 1, modify the parameters to achieve:

Chern number $C = 0$ (trivial phase)
Chern number $C = 1$ (topological phase)
Identify the critical parameter value where the topological phase transition occurs
Plot the phase diagram in parameter space

Exercise 4.3: Weyl Point Chirality (Advanced)

For the Weyl Hamiltonian in Code Example 2:

Derive the Berry curvature $\Omega^z(k_x, k_y, k_z)$ analytically near the Weyl point at $(0,0,m)$
Calculate the topological charge by integrating Berry curvature over a small sphere
Verify numerically that the total chirality in the Brillouin zone sums to zero

Exercise 4.4: Edge State Robustness (Intermediate)

Modify Code Example 3 to add random on-site disorder:

Add random energy disorder $\epsilon_i \in [-W, W]$ to each site
Calculate edge state dispersion for different disorder strengths $W$
Determine the critical disorder strength where edge states are destroyed
Compare with bulk band gap size

Exercise 4.5: Phonon Hall Conductivity (Intermediate)

Using the phonon Hall effect calculation (Code Example 4):

Derive the temperature dependence of $\kappa_{xy}$ in the low-temperature limit ($T \to 0$)
Derive the high-temperature limit ($k_B T \gg \hbar\omega$)
Numerically verify these limits using the code
Estimate the magnitude of $\kappa_{xy}$ for realistic materials (use α-quartz parameters)

Exercise 4.6: Valley Chern Numbers (Advanced)

For the honeycomb lattice model (Code Example 5):

Prove analytically that $C_K + C_{K'} = 0$ due to time-reversal symmetry
Calculate valley Chern numbers as a function of sublattice energy difference $\delta$
Identify the phase transition point where valley Chern numbers change
Discuss how valley-polarized edge states could be experimentally detected

Exercise 4.7: Symmetry Analysis (Intermediate)

Consider α-quartz (space group P3₁21):

List all symmetry operations of this space group
Identify which symmetries protect the Weyl points
Determine how applying pressure (reducing symmetry to P3₁) affects topology
Predict whether Weyl points can be moved or annihilated by pressure

Exercise 4.8: Phononic Metamaterial Design (Advanced Project)

Design a 2D phononic crystal with topological edge states:

Choose a lattice geometry (honeycomb, kagome, or square with decorations)
Define mass and spring constant distributions to break inversion symmetry
Calculate phonon band structure and Berry curvature
Verify nontrivial Chern number
Compute edge state dispersion in ribbon geometry
Propose an experimental realization (materials, fabrication method)
Estimate operating frequency range and measurement techniques