Phonon Transport at Interfaces
1. Introduction
In solids, atoms are not perfectly still; they vibrate. These vibrations carry heat, much like data packets travel across a network. At ann interface(boundary between two materials), the way vibrations pass through dertermines how well heat flows.
Our goal: simulate this “vibration transport” using only linear algebra and Python- no advanced physics background needed.
2. Equations of Motion -> Matrix Form
Imagine each atom is a point mass attached to its neighbours by springs. The physics law is just newtons second law: \(M\ddot{u}(t) = -K u(t)\). * M : diagonal matrix of atomic masses * K : stiffness (spring constant) matrix * u(t) : displacement vector (positions)
this is just like solving a second-order system of ODEs. Instead of time-stepping, we assume vibrations oscillate like \(u(t) = U e^{-i \omega t}\) and reduce the ODE to an eigenvalue problem \(\big( K - \omega^2 M \big) U = 0\)
This tells us the frequencies \(\omega\) (the signal frequencies) and shapes \(U\) (like signal modes).
3. The interface problem
- Think of a chain of atoms: left half are “Material A”, right half are “Material B”.
- They are joined in the middle -> the interface.
- Each material supports certain vibration patterns(“phonons”).
- At the interface, do vibrations transmit or reflect?
This is what we want to compute.
4. Lattice Dynamics + Green’s Functions
We use a standard trick: * Transform the equations into frequency domain. * Replace solving Ordinary differential equations with matrix inversion.
The central object is the Green’s function: \[ G(\omega) = \Big[ \, (\omega^2 + i\eta) I - D - \Sigma_L(\omega) - \Sigma_R(\omega) \, \Big]^{-1} \]
Where:
- \(D = M^{-1/2} K M^{-1/2}\) is the mass-normalized dynamical matrix
- \(\Sigma_L, \Sigma_R\) are the self-energies of the left/right leads
- \(i\eta\) is a small numerical damping term (keeps the matrix invertible)
Transmission Function
Once \(G(\omega)\) is known, the phonon transmission is:
\[ T(\omega) = \mathrm{Tr} \Big[ \, \Gamma_L \, G \, \Gamma_R \, G^\dagger \, \Big] \]
where
\[ \Gamma = i \, \big( \Sigma - \Sigma^\dagger \big) \]
5. Implementing in Python
First, lets set up a toy 1D model with two materials. Each atom is a mass, each bond a spring.
Build the device region: 8 atoms in a chain. Left half: light atoms with strong springs. Right half: heavy atoms with weak springs.
5a) Mass matrix and mass-normalized dynamical matrix (D)
We convert the force-constant matrix (K) into the mass-normalized dynamical matrix
\[ D \;=\; M^{-1/2}\, K\, M^{-1/2}, \]
where \(M=\mathrm{diag}(m_1,\dots,m_N)\) comes from the per-atom masses.
This converts \((K-\omega^2 M)U=0\) into a standard eigenproblem \(D U=\omega^2 U\).
device_force_constants
device_force_constants (N:int, k_left:float, k_right:float)
Build the NxN force-constant (Hessian) matrix K for a 1D chain with an interface. - First half bonds use k_left - Second half bonds use k_right Nearest-neighbor springs only (tridiagonal K).
6) Leads and surface Green’s functions (Sancho–Rubio)
We treat the left/right semi-infinite materials (“leads”) via their **surface Green’s functions\[ g_L, g_R\] For a 1D monatomic chain (1 DOF per cell), the mass-normalized on-site and coupling blocks are
\[ H_0 = \frac{2k}{m}, \qquad H_1 = -\frac{k}{m}. \]
We compute the surface Green’s function using the Sancho–Rubio decimation (a fixed-point block-matrix iteration).
surface_g_mono
surface_g_mono (w:float, m:float, k:float, eta:float=1e-06)
Retarded surface Green’s function g^r(ω) for a semi-infinite 1D monatomic chain in the mass-normalized* basis.
Chain blocks: H0 = 2k/m, H1 = -k/m (scalar). g satisfies: (h1^2) g^2 - (z - h0) g + 1 = 0, z=ω^2+iη. Pick the root with Im(g) < 0 (retarded condition).*
sancho_rubio
sancho_rubio (H0:numpy.ndarray, H1:numpy.ndarray, w2:float, eta:float=1e-12, max_iter:int=200)
Surface Green’s function g^r for a semi-infinite 1D chain (1x1 blocks).
lead_blocks
lead_blocks (mass:float, kspring:float)
7) Self-energies, device Green’s function, and transmission
Couple the device’s end atoms to the leads. In the mass-normalized basis, the coupling vectors are
\[ V_L[0]=-\sqrt{k_L/m_L}, \quad V_R[-1]=-\sqrt{k_R/m_R}. \]
Then compute: - Lead self-energies \(\Sigma_{L,R} = V_{L,R}\, g_{L,R}\, V_{L,R}^T\) - Broadening matrices \(\Gamma = i(\Sigma - \Sigma^\dagger)\) - Device Green’s function \(G = \big[(\omega^2+i\eta)I - D - \Sigma_L - \Sigma_R\big]^{-1}\) - Transmission \(T(\omega) = \mathrm{Tr}\big[\Gamma_L\, G\, \Gamma_R\, G^\dagger\big]\)
transmission_at
transmission_at (w:float, eta:float=1e-06)
gamma_from_sigma
gamma_from_sigma (Sigma:numpy.ndarray)
8) Frequency sweep and plot \(T(\omega)\)
We scan $ $ and visualize the transmission spectrum.
9) System eigenmodes (the Gordiz–Henry “correct modes”)
Diagonalize the device’s dynamical matrix (D) to obtain eigenpairs \((\omega_n^2, U_n)\).
These are the system modes spanning the interface region (not bulk modes).
They provide a physically correct basis to interpret which patterns couple to the leads.
A quick heuristic: modes with larger amplitude on the end atoms couple better to the leads.
10) (Optional) Spectral channels via the spectral function
For deeper mode analysis, compute the spectral function $ A()=i(G-G^) $ and diagonalize it to get frequency-resolved transport channels.
spectral_channels
spectral_channels (G:numpy.ndarray)