瞬态热光栅的非扩散弛豫

电热输运

瞬态热光栅的非扩散弛豫

  1. Collins K C, Maznev A A, Tian Z, et al. Non-diffusive relaxation of a transient thermal grating analyzed with the Boltzmann transport equation[J]. Journal of Applied Physics, 2013, 114(10).
  2. 顾樵. 数学物理方法[J]. 2012.

瞬态热光栅

瞬态热光栅（Transient Thermal Grating, TTG）技术用两束激光交叉产生干涉条纹，从而在样品表面产生空间上正弦分布的温度条纹，之后用一束探测光检测温度场的演化情况。其中光栅周期\(L\)是可以调控的，当\(L\)和声子自由程相当的时候，就会产生非扩散的热输运过程，也就是弹道输运。这个实验可以用来测量声子的自由程分布。

image-20231122143237984

傅里叶解

如果这个导热过程满足傅里叶定律，那么定解问题可以写成 \[ \left\{ \begin{array}{l} \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} \quad(0 < x < L, t > 0) \\ u|_{t=0} = T_\text{max}\cos(qx), q = \frac{2\pi}{L} \quad (0 \leq x \leq L)\\ \frac{\partial u}{\partial x}\big|_{x = 0} = \frac{\partial u}{\partial x}\big|_{x = L}= 0 \quad (t>0) \end{array} \right. \] 其中\(\alpha\)是热扩散率 \[ \alpha = \frac{k}{C} \] 不论是傅里叶定律还是BTE，都可以根据解的空间周期性构造形式解，最后用本征函数法求出时间项系数。

根据空间的边界条件构造本征函数集 \[ \left\{\cos(nqx)\right\}\quad (n = 0, 1, 2, \ldots) \] 设定解问题的形式解为 \[ u(x, t) = \sum_{n=0}^{\infty} T_n(t) \cos(nqx) \] 其中，\(T_n(t)\)是时间\(t\)依赖的展开系数，把这个形式解代入控制方程中得到 \[ \sum_{n=0}^{\infty} \left[T^\prime_n(t) + \alpha n^2q^2 T_n(t) \right] \cos(nqx) = 0 \] 于是 \[ T^\prime_n(t) + \alpha n^2q^2 T_n(t) = 0 \] 它的通解为 \[ T_n(t) = C_n \exp\left(- \alpha n^2q^2 t\right) \] 将它代入到\(u(x, t)\)的形式解中 \[ u(x, t) = \sum_{n=0}^{\infty} C_n \exp\left(- \alpha n^2q^2 t\right) \cos(nqx) \] 现在利用初始条件确定展开系数\(C_n\)，设初始温度分布为\(\phi(x)\)，我们有 \[ \phi(x) = \sum_{n=0}^{\infty} C_n\cos(nqx) \] 三角函数系是线性空间中的一组正交基，于是我们有 \[ C_n = \frac{2}{L} \int_0^L \phi(x) \cos(nqx) \mathrm{d}x \] 代入初始温度 \[ C_n = \frac{2}{L} \int_0^L T_\text{max}\cos(qx)\cos(nqx)\mathrm{d}x = T_\text{max}\delta_{1n} \] 于是我们的级数解只有一级的一项 \[ u(x, t) = T_\text{max} \exp\left(- \alpha q^2 t\right) \cos(nqx) \]

假设我们关注某一个位置处的温度随时间的演化情况，可以得到 \[ \Upsilon (\zeta) = \frac{u(x, t)}{u(x, 0)} = \exp(-\alpha q^2 t) = \exp(-\beta \zeta) \] 其中\(\beta = \alpha q^2 \tau, \zeta = t / \tau\)。这是一个纯粹的指数衰减过程，光栅距离越近、材料的导热系数越高，温度场衰减得越快。为了和后面BTE的结果做对比，定义 \[ \eta = qv\tau = q\Lambda \] 根据动理学关系，有 \[ k = \frac{1}{3}Cv\tau \] 于是存在 \[ \beta = \frac{k}{C}q^2\tau = \frac{1}{3}v^2q^2\tau^2 = \frac{1}{3}\eta^2 \] Fourier

BTE

解析解

弛豫时间近似下一维声子玻尔兹曼方程可以写成 \[ \frac{\partial g}{\partial t} + \mu v \frac{\partial g}{\partial x} = \frac{g_0 - g}{\tau} \]

\[ g_0 = \frac{1}{4\pi} \hbar\omega D(\omega) f_\text{BE}(T) \approx \frac{1}{4\pi}C_\omega T \]

\(g(x, t, \mu)\)代表单位频率单位立体角中的声子能量。假设\(g\)和\(T\)有以下形式的周期性通解 \[ g(x, t, \mu) = \tilde{g}(t, \mu)\exp(iqx), T(t, x) = \tilde{T}(t)\exp(iqx) \] 把假设的\(g\)的形式解代入到玻尔兹曼方程中 \[ \exp(iqx) \frac{\partial \tilde{g}}{\partial t} + \mu v iq\exp(iqx)\tilde{g} = \exp(iqx) \frac{\tilde{g_0} - \tilde{g}}{\tau} \] 于是 \[ \frac{\partial \tilde{g}}{\partial t} + \gamma \tilde{g} = \frac{\tilde{g_0}}{\tau} \] 其中 \[ \gamma = (1 + iq\mu v\tau) / \tau \] 这变成了一个一阶ODE，可以直接用常数变易法写出通解 \[ \tilde{g}(t, \mu) = \frac{1}{\tau} \int_0^t e^{\gamma (t^\prime - t)}\tilde{g}_0(t^\prime)\mathrm{d}t^\prime + \tilde{g}_0(0)e^{-\gamma t} \] 此外，在灰体假设下，能量守恒可以写成 \[ 2g_0(x,t) = \int_{-1}^1 g(x, t, \mu) \mathrm{d}\mu \] 代入一阶ODE得到的通解 \[ \begin{aligned} \tilde{g}_0(x, t) &= \frac{1}{2\tau}\int_{-1}^{1} \int_0^t e^{\gamma (t^\prime - t)}\tilde{g}_0(t^\prime)\mathrm{d}t^\prime \mathrm{d}\mu + \frac{1}{2}\int_{-1}^{1}\tilde{g}_0(0)e^{-\gamma t} \mathrm{d}\mu \\ &= \frac{1}{2\tau}\int_{-1}^{1} \int_0^t e^{(1/\tau + iq\mu v ) (t^\prime - t)}\tilde{g}_0(t^\prime)\mathrm{d}t^\prime \mathrm{d}\mu + \frac{1}{2}\int_{-1}^{1}\tilde{g}_0(0)e^{-(1/\tau + iq\mu v ) t} \mathrm{d}\mu \\ \end{aligned} \] 积分得到 \[ \tilde{g}_0(x, t)=\frac{1}{\tau} \int_0^te^{\left(t - t^{\prime}\right)/\tau} \tilde{g}_0\left(t^{\prime}\right) \operatorname{sinc}\left(qv \left(t-t^{\prime}\right)\right) d t^{\prime}+\tilde{g}_0(0) \operatorname{sinc}(qvt) e^{-t / \tau} \] 其中sinc函数定义为 \[ \operatorname{sinc}(x) = \frac{\sin(x)}{x} \] 代入 \[ \zeta=\frac{t}{\tau}, \quad \eta=q v \tau=\frac{2 \pi \Lambda}{L}, \quad \Upsilon=\frac{\tilde{g}_o(t)}{\tilde{g}_o(0)}=\frac{\tilde{T}(t)}{\tilde{T}(0)} \] 得到 \[ \Upsilon(\zeta)= \operatorname{sinc}(\eta \zeta) e^{-\zeta}+\int_0^\zeta \Upsilon\left(\zeta^{\prime}\right) \operatorname{sinc}\left(\eta\left(\zeta^{\prime}-\zeta\right)\right) e^{\left(\zeta^{\prime}-\zeta\right)} d \zeta^{\prime} \] 定义 \[ \operatorname{sinc}\left(\eta\left(\zeta^{\prime}-\zeta\right)\right) e^{\left(\zeta^{\prime}-\zeta\right)} = K(\zeta, \zeta^\prime) \] 可以得到 \[ \Upsilon(\zeta)= \operatorname{sinc}(\eta \zeta) e^{-\zeta}+\int_0^\zeta \Upsilon\left(\zeta^{\prime}\right) K(\zeta, \zeta^\prime) d \zeta^{\prime} \] 这是第二类沃尔泰拉积分方程（Volterra integral equation of the second kind），我们可以数值求解这个方程，把积分项挪到方程的左边，然后把\(\zeta\)在\((0, \zeta_\text{max})\)离散化，可以得到 \[ \Upsilon_i - \int_0^{\zeta_i} \Upsilon\left(\zeta^{\prime}\right) K(\zeta, \zeta^\prime)d \zeta^{\prime} = \operatorname{sinc}(\eta \zeta_i) e^{-\zeta_i} \] 左边的积分可以用梯形法用离散求和的方式表示 \[ \int_a^b f(x) d x \approx \Delta x \left(\frac{1}{2}f\left(x_0\right)+ f\left(x_1\right)+ f\left(x_2\right)+\cdots+ f\left(x_{N-1}\right)+\frac{1}{2} f\left(x_N\right)\right) \] 可以得到 \[ \Upsilon_i - \Delta x \left[\frac{1}{2}K(\zeta_i, \zeta_0)\Upsilon_0+K(\zeta_i, \zeta_1)\Upsilon_1 + \ldots + \frac{1}{2} K(\zeta_i, \zeta_{i})\Upsilon_i\right]= \operatorname{sinc}(\eta \zeta_i) e^{-\zeta_i} \] 这表示一个线性方程组 \[ (I - A)\cdot \Upsilon = B \] 其中 \[ B_i = \operatorname{sinc}(\eta \zeta_i) e^{-\zeta_i} \]

\[ A_{ij} = \left\{\begin{array}{l} \Delta x \cdot \frac{1}{2} \cdot K(\zeta_i, \zeta_0), \quad j = 0\\ \Delta x \cdot K(\zeta_i, \zeta_j),\quad 0 < j < i\\ \Delta x \cdot K(\zeta_i, \zeta_i), \quad j=i \\ 0, \quad j > i \end{array}\right. \]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
import numpy as np
from scipy.linalg import solve
import matplotlib.pyplot as plt

sinc = lambda x: np.sinc(x / np.pi)

def outside_integral(zeta, eta):
    return sinc(eta * zeta) * np.exp(-zeta)

def build_matrix(zeta_values, eta):
    # Create the matrix for the system of linear equations
    n_points = len(zeta_values)
    A = np.eye(n_points)
    B = np.zeros(n_points)

    for i in range(n_points):
        zeta_i = zeta_values[i]
        B[i] = outside_integral(zeta_i, eta)
        
        for j in range(i + 1):  # Only need to go up to i because of the integral limit
            zeta_j = zeta_values[j]
            integrand = sinc(eta * (zeta_j - zeta_i)) * np.exp(zeta_j - zeta_i)
            if j == 0 or j == i:
                integrand /= 2  # Adjust endpoints for the trapezoidal rule

            A[i, j] -= integrand * (zeta_values[1] - zeta_values[0])  # Trapezoidal rule
            
    return A, B

zeta_values = np.linspace(0, 10, 200)
upsilon_values = solve(*build_matrix(zeta_values, 5))
plt.plot(zeta_values, upsilon_values)

ttg

数值解

对于Ensemble MC需要给定一个初始的声子分布按时间演化。对于Kinetic MC将每一个跟踪发射的声子时间都设置为0，再按照初始的能量分布作为权重去抽样就可以了。令初始的温度分布为 \[ T = T_\text{max} \left(1 + \cos(qx)\right), \quad x \in (0, L) \] 这是概率密度函数（PDF），积分归一化得到累积分布函数（CDF） \[ F(x) = \frac{1}{L} \left(x + \frac{1}{q} \sin(qx)\right) \] 抽样随机数\(R\in(0, 1)\) \[ R = \frac{1}{L} \left(x + \frac{1}{q} \sin(qx)\right) \] 求逆可以得到\(x\)

ttg-sample

模拟模拟可以得到和解析解一样的结果

ttgmc