thermal spike hot spot phonon transport non Fourier heat

# Non-Fourier Phonon Heat Transport and Thermal Spike Kinetics in Nanoscale Devices

## 1. Introduction: Breakdown of Fourier's Law at Nanoscale

Classical Fourier's law of heat conduction:
$$q = -\kappa abla T$$

assumes that heat diffuses instantaneously through the material via local temperature gradients. This assumption fails when device dimensions become comparable to or smaller than the phonon mean free path $\lambda_{ ext{ph}}$. In modern GAA (Gate-All-Around) FinFETs and CFETs operating at sub-3nm channel lengths, the characteristic thermal diffusion length becomes shorter than typical phonon wavelengths, violating the hydrodynamic limit.

The breakdown manifests in:
- Ballistic phonon transport: Phonons traverse the device without scattering (direct transmission)
- Thermal non-equilibrium: Local temperature cannot be uniquely defined; energy is carried by directed phonon flux
- Transient overshoot: Heat fronts propagate faster than predicted by classical diffusion
- Thermal memory effects: The system "remembers" previous states over phonon relaxation times

These effects are critical during:
1. High-power RF switching (nanosecond pulse injection)
2. Ion implantation damage (thermal spike dynamics during ion stopping)
3. Joule heating in narrow current channels (hot spots exceeding equilibrium by 100s of K)

## 2. Phonon Boltzmann Transport Equation

The fundamental equation governing phonon distribution is the phonon Boltzmann transport equation (BTE):

$$\frac{\partial f_ u}{\partial t} + \vec{v}_ u \cdot abla f_ u + \frac{\vec{F}}{m^*} \cdot abla_k f_ u = \left(\frac{\partial f_ u}{\partial t} ight)_{ ext{scatt}}$$

where:
- $f_
u(\vec{r}, \vec{k}, t)$ is the phonon distribution function for mode $
u$ (phonon type: acoustic/optical)
- $\vec{v}_
u = \partial \omega_
u / \partial \vec{k}$ is the group velocity
- $\vec{F}$ is an external force (strain-induced acoustic forces)
- The scattering term accounts for phonon-phonon, phonon-defect, and boundary collisions

For heat transport without external forces, the steady-state BTE simplifies to:

$$\vec{v}_ u \cdot abla f_ u = \left(\frac{\partial f_ u}{\partial t} ight)_{ ext{scatt}}$$

The scattering term in the relaxation time approximation:

$$\left(\frac{\partial f_ u}{\partial t} ight)_{ ext{scatt}} = -\frac{f_ u - f_ u^{eq}}{ au_ u}$$

where $ au_
u$ is the mode-dependent relaxation time and $f_
u^{eq}$ is the equilibrium Bose-Einstein distribution.

## 3. Callaway Model and Effective Mean Free Path

The Callaway model provides an analytical solution to the phonon BTE by separating processes into normal (momentum-conserving) and umklapp (momentum-destroying) scattering:

$$ au_ u^{-1} = au_N^{-1} + au_U^{-1}$$

where $ au_N$ is the normal process time and $ au_U$ is the umklapp time. Umklapp processes dominate at high temperature and limit thermal conductivity.

The effective mean free path:
$$\lambda_{ ext{eff}} = \min(\lambda_{ ext{scatt}}, \lambda_{ ext{device}})$$

• When $\lambda_{ ext{scatt}} \gg L$ (device length): Ballistic limit — phonons cross without scattering, $\kappa_{ ext{eff}} \propto \sqrt{L}$
• When $\lambda_{ ext{scatt}} \ll L$: Diffusive limit — classical Fourier law applies
• When $\lambda_{ ext{scatt}} \sim L$: Transition regime — non-Fourier corrections dominate

For silicon at 300 K:
- Optical phonon mean free path: $\lambda_{ ext{opt}} \sim 10-20$ nm
- Acoustic phonon mean free path: $\lambda_{ ext{ac}} \sim 100-300$ nm

## 4. Guyer-Krumhansl Non-Fourier Heat Transport

The Guyer-Krumhansl equation extends Fourier's law to include non-local and time-dependent effects:

$$\frac{\partial Q}{\partial t} + Q = -\kappa abla T + \lambda_{ ext{ph}}^2 \frac{\partial}{\partial t}(\kappa abla^2 T) + \lambda_{ ext{ph}}^2 \kappa abla^4 T$$

where $Q$ is the heat flux, $\kappa$ is thermal conductivity, and $\lambda_{ ext{ph}}$ is the phonon mean free path.

This hyperbolic PDE replaces the parabolic heat equation:
$$ ho c_v \frac{\partial T}{\partial t} = \kappa abla^2 T$$

The Cattaneo relaxation time $ au_c \sim \lambda_{ ext{ph}} / v_s$ (sound velocity $v_s$) introduces finite heat pulse speed rather than instantaneous diffusion.

Rearranging to standard form:
$$ au_c \frac{\partial^2 T}{\partial t^2} + \frac{\partial T}{\partial t} = D_{ ext{eff}} abla^2 T + au_c D_{ ext{eff}} abla^4 T$$

where $D_{ ext{eff}} = \kappa / (
ho c_v)$ is the effective thermal diffusivity.

## 5. Thermal Spike Kinetics: Ion Implantation Energy Deposition

During ion implantation, an energetic ion (E = 1-100 keV) dumps energy via:
1. Electronic stopping (e-e collisions): Fast, creates electron-hole pairs
2. Nuclear stopping (n-n collisions): Slower, creates atomic displacements

The thermal spike forms when:
- Electronic energy is converted to lattice vibrations (electron-phonon coupling)
- Temperature rises dramatically in $\sim 10$-100 fs timescale
- Localized region (cylinder radius ~10 nm) reaches melting or amorphous state

The peak temperature at radius $r$ during ion impact:

$$T_{ ext{peak}}(r) = T_0 + \frac{E_{ ext{dep}}}{2\pi ho c_v r L_{ ext{track}}} \exp\left(-\frac{r^2}{\lambda_{ ext{th}}^2} ight)$$

where:
- $E_{ ext{dep}}$ is deposited energy
- $L_{ ext{track}}$ is ionization track length
- $\lambda_{ ext{th}} \sim$ 10 nm is thermal diffusion length

## 6. High-Power RF Switching and Joule Heating

During RF switching at high power, current crowding in narrow channels (sub-100 nm) creates hot spots:

$$P_{ ext{Joule}} = I^2 R = \frac{V_{ ext{DS}}^2}{ ho(T) A}$$

where resistivity $
ho(T)$ increases with temperature (positive temperature coefficient in semiconductors).

The self-heating feedback loop:
$$\frac{dT}{dt} = \frac{P_{ ext{Joule}}}{ ho c_v} - \frac{\kappa}{( ho c_v) L_{ ext{th}}^2}(T - T_{ ext{ambient}})$$

For sub-3nm channels with $A \sim 10$ nm², joule power density can reach 10 MW/cm³, creating transient temperatures >500 K.

## 7. Python Implementation: Ballistic-Diffusive Heat Transport

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.sparse import diags, lil_matrix
from scipy.sparse.linalg import spsolve

def phonon_mean_free_path(T, material='Si'):
    """
    Compute phonon mean free path vs temperature.
    
    Parameters:
    -----------
    T : array (K)
        Temperature array
    material : str
        'Si' or 'SiO2' or 'TiN'
    
    Returns:
    --------
    lambda_ph : array (m)
        Phonon mean free path
    """
    
    if material == 'Si':
        # Empirical fit: λ_ph(T) ∝ 1/T at high T (umklapp dominated)
        lambda_0 = 100e-9  # Reference: 100 nm at 300 K
        T_ref = 300
        lambda_ph = lambda_0 * (T_ref / (T + 1e-10))
    elif material == 'SiO2':
        # Silica: more defects, shorter mfp
        lambda_0 = 10e-9
        T_ref = 300
        lambda_ph = lambda_0 * (T_ref / (T + 1e-10))**0.5
    else:
        lambda_ph = 50e-9 * np.ones_like(T)
    
    return np.clip(lambda_ph, 1e-10, 500e-9)  # Physical bounds

def thermal_conductivity_ballistic(T, L_device, material='Si'):
    """
    Compute effective thermal conductivity with ballistic correction.
    
    κ_eff = κ_bulk / (1 + β * λ_ph(T) / L_device)
    
    Parameters:
    -----------
    T : array (K)
        Temperature
    L_device : float (m)
        Device length scale
    material : str
        Material type
    
    Returns:
    --------
    kappa_eff : array (W/m·K)
        Effective thermal conductivity
    """
    
    if material == 'Si':
        # Bulk Si thermal conductivity (empirical)
        kappa_bulk = 150 / (1 + 0.004 * (T - 300))  # W/(m·K)
    else:
        kappa_bulk = 10 * np.ones_like(T)
    
    lambda_ph = phonon_mean_free_path(T, material)
    
    # Ballistic reduction factor
    Knudsen = lambda_ph / (L_device + 1e-12)
    reduction_factor = 1.0 / (1.0 + 2.0 * Knudsen)
    
    kappa_eff = kappa_bulk * reduction_factor
    
    return kappa_eff

def solve_guyer_krumhansl_1d(x, t, T_init, T_boundary, power_density_array, 
                              L_device, dt_out, material='Si'):
    """
    Solve Guyer-Krumhansl non-Fourier heat equation in 1D.
    
    ∂T/∂t = D_eff ∂²T/∂x² + τ_c D_eff ∂⁴T/∂x⁴ + P(x,t)/ρc_v
    
    Parameters:
    -----------
    x : array (m)
        Spatial grid
    t : array (s)
        Time output points
    T_init : array (K)
        Initial temperature
    T_boundary : float (K)
        Boundary temperature
    power_density_array : array (W/m³)
        Time-dependent power input
    L_device : float (m)
        Device length
    dt_out : float (s)
        Output time step
    material : str
        Material type
    
    Returns:
    --------
    T_history : array (N_space, N_time)
        Temperature evolution
    t_out : array
        Output time array
    """
    
    # Material properties
    if material == 'Si':
        rho = 2330  # kg/m³
        c_v = 700   # J/(kg·K)
    else:
        rho = 2200
        c_v = 600
    
    dx = x[1] - x[0]
    dt = 0.1 * dx**2 / (100e-9 / (rho * c_v))  # Stable timestep
    
    N_space = len(x)
    N_time_steps = int(t[-1] / dt)
    t_out = np.linspace(0, t[-1], len(t))
    
    T = T_init.copy()
    T_history = np.zeros((N_space, len(t_out)))
    T_history[:, 0] = T
    
    # Finite difference stencils
    L_laplacian = diags([-2, 1, 1], [0, -1, 1], shape=(N_space, N_space)) / dx**2
    L_laplacian[0, :] = 0
    L_laplacian[-1, :] = 0
    L_laplacian = lil_matrix(L_laplacian)
    
    L_bilaplacian = L_laplacian @ L_laplacian
    L_bilaplacian = L_bilaplacian.tocsr()
    
    output_idx = 0
    
    for n in range(N_time_steps):
        t_current = n * dt
        
        # Temperature-dependent properties
        kappa_eff = thermal_conductivity_ballistic(T, L_device, material)
        lambda_ph = phonon_mean_free_path(T, material)
        D_eff = np.mean(kappa_eff) / (rho * c_v)
        tau_c = np.mean(lambda_ph) / 3000  # Sound velocity ~3000 m/s
        
        # Power input (interpolate from array)
        P_idx = min(int(t_current / (t[-1] / len(power_density_array))), 
                    len(power_density_array) - 1)
        P = power_density_array[P_idx]
        
        # RHS: D_eff ∇²T + τ_c D_eff ∇⁴T + P/(ρc_v)
        grad2T = L_laplacian @ T
        grad4T = L_bilaplacian @ T
        
        dT = D_eff * grad2T + tau_c * D_eff * grad4T + P / (rho * c_v)
        
        # Forward Euler
        T_new = T + dt * dT
        T_new[0] = T_boundary
        T_new[-1] = T_boundary
        T = T_new
        
        # Store output
        if abs(t_current - t_out[output_idx]) < dt:
            T_history[:, output_idx] = T
            output_idx = min(output_idx + 1, len(t_out) - 1)
    
    return T_history, t_out

# Simulation setup
x = np.linspace(0, 100e-9, 500)  # 0-100 nm channel
t = np.logspace(-12, -9, 100)    # 1 ps to 1 ns

# Initial condition: room temperature
T_init = 300 * np.ones_like(x)

# Power density: thermal spike from ion impact
# Gaussian pulse centered at t=100 fs, width ~50 fs
t_pulse = 100e-15
sigma_t = 50e-15
power_scale = 1e18  # 1 EW/m³ peak power

power_density = power_scale * np.exp(-((t - t_pulse)**2) / (2 * sigma_t**2))

# Solve GK equation
T_history, t_out = solve_guyer_krumhansl_1d(x, t, T_init, 300, power_density, 
                                            100e-9, 1e-12, material='Si')

# Plotting
fig, axes = plt.subplots(2, 2, figsize=(13, 10))

# Panel A: Temperature profiles at different times
time_indices = [0, 20, 50, 80, 99]
colors = plt.cm.hot(np.linspace(0.2, 1, len(time_indices)))

for idx, i in enumerate(time_indices):
    axes[0, 0].plot(x*1e9, T_history[:, i], color=colors[idx], 
                   label=f't={t_out[i]*1e15:.0f} fs', linewidth=2)

axes[0, 0].set_xlabel('Position (nm)', fontsize=11)
axes[0, 0].set_ylabel('Temperature (K)', fontsize=11)
axes[0, 0].set_title('Thermal Spike Evolution', fontsize=12, fontweight='bold')
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)

# Panel B: Phonon mean free path vs T
T_range = np.linspace(300, 1000, 100)
lambda_ph = phonon_mean_free_path(T_range, 'Si')
axes[0, 1].semilogy(T_range, lambda_ph*1e9, 'b-', linewidth=2)
axes[0, 1].set_xlabel('Temperature (K)', fontsize=11)
axes[0, 1].set_ylabel('λ_ph (nm)', fontsize=11)
axes[0, 1].set_title('Phonon Mean Free Path', fontsize=12, fontweight='bold')
axes[0, 1].grid(True, alpha=0.3, which='both')

# Panel C: Contour plot of thermal evolution
X_plot, T_plot = np.meshgrid(x*1e9, t_out*1e15)
contour = axes[1, 0].contourf(X_plot, T_plot, T_history.T, levels=20, cmap='hot')
axes[1, 0].set_xlabel('Position (nm)', fontsize=11)
axes[1, 0].set_ylabel('Time (fs)', fontsize=11)
axes[1, 0].set_title('Thermal Spike Contour (Space-Time)', fontsize=12, fontweight='bold')
plt.colorbar(contour, ax=axes[1, 0], label='T (K)')

# Panel D: Peak temperature vs time
T_max_history = np.max(T_history, axis=0)
axes[1, 1].plot(t_out*1e15, T_max_history, 'r-o', linewidth=2, markersize=4)
axes[1, 1].set_xlabel('Time (fs)', fontsize=11)
axes[1, 1].set_ylabel('Peak Temperature (K)', fontsize=11)
axes[1, 1].set_title('Thermal Spike Decay Kinetics', fontsize=12, fontweight='bold')
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('thermal_spike_phonon.png', dpi=150, bbox_inches='tight')
plt.show()

print("=== Guyer-Krumhansl Non-Fourier Heat Transport Simulation Complete ===")
print(f"Peak temperature: {np.max(T_history):.1f} K")
print(f"Time to cool to <400 K: {t_out[np.where(T_max_history < 400)[0][0]]*1e15:.1f} fs")

## 8. Phonon Drag and Dember Effect

At high thermal gradients, the phonon drag effect causes an enhancement of the Seebeck coefficient:

$$S_{ ext{total}} = S_0 + S_{ ext{drag}}$$

where $S_{ ext{drag}} \propto (\partial au / \partial E)
abla T$ represents momentum transfer from phonons to electrons.

The Dember effect (thermal voltage without applied heat gradient in semiconductors) arises from differential diffusion of electrons and holes under photoexcitation:

$$V_{ ext{Dember}} = \frac{k_B T}{q} \ln\left(\frac{\mu_e}{\mu_h} ight)$$

Both effects are relevant in multi-layered heterostructures where phonon populations in different regions become thermally decoupled.

## 9. Transient Hot Carrier Relaxation

Hot carriers (electrons/holes with $E > E_F + k_B T$) relax via:

1. Optical phonon scattering (picosecond timescale)
2. Acoustic phonon scattering (shorter mean free path)
3. Impact ionization (creates additional carriers)

The hot carrier cooling time:
$$ au_{ ext{cool}} = \hbar \omega_{ ext{opt}} / (E - E_c)$$

For GaAs: $ au_{ ext{cool}} \sim 100$ fs for carriers 100 meV above bandedge, rising to picoseconds for cooler carriers.

## 10. Two-Temperature Model

In metals and heavily doped semiconductors, electrons and lattice temperatures may differ:

$$ ho c_{e} \frac{\partial T_e}{\partial t} = \kappa_e abla^2 T_e - G(T_e - T_l) + P_{ ext{ext}}$$
$$ ho c_l \frac{\partial T_l}{\partial t} = \kappa_l abla^2 T_l + G(T_e - T_l)$$

where $G$ is the electron-phonon coupling constant (~10¹⁷ W/m³·K in metals).

This model predicts transient overshoot of electron temperature far above lattice temperature during ultrafast pulse excitation.

## 11. Thermal Boundary Resistance (Kapitza Resistance)

At interfaces, phonons scatter inelastically, creating a thermal boundary resistance $R_B$:

$$R_B = \frac{1}{C_l v_l} \sqrt{\frac{\pi m_1 m_2}{k_B T(m_1 + m_2)}}$$

where $C_l$ is phonon heat capacity, $v_l$ is sound velocity, and $m_{1,2}$ are reduced masses.

For Si/SiO₂ interface: $R_B \sim 10^{-8}$ m²·K/W, limiting thermal extraction in SOI devices.

## 12. Applications to GAAFET/CFET Thermal Management

Sub-3nm GAA transistors suffer from:
- Reduced contact area for heat extraction
- High power density ($>10$ W/μm²)
- Thermal runaway risk at high bias

Mitigation strategies:
1. Thermal interface layers (TiN, W) with high thermal conductivity
2. Back-gate cooling (heterostructure contacts)
3. Phononic crystals (periodic defects to scatter phonons below critical frequency)
4. Operation at lower $V_{DS}$ (reduces Joule heating, limits frequency)

The balance between performance and thermal stability defines the power-temperature envelope.