声子与电声耦合

一. Classical Models

单谐振子 满足运动方程：\(F=m\frac{d^2x}{dt^2}=-kt\)，解形式为\(x=Ae^{i(wt+\psi)}\)。 耦合谐振子链：

其中$x_n=na$是第n个振子的位置，位移变化可以写成\(u_n(t)=x_n(t)-na\)，可以写出运动方程： \(m\ddot{u}_n=-\lambda(2u_n-u_{n-1}-u_{n+1})\) 加入周期性边界条件$u_{n+N}=u_n$，解可以写成： $u_n=Ae^{-iwt-ikna}$,其中$k\in[-\frac{\pi}{a},\frac{\pi}{a})$ ,且由于周期性边界条件，取值为$k=\frac{2\pi}{Na}l$, with $l=-\frac{N}{2},…,\frac{N}{2}$. 代入运动方程可以解出 \(\omega=2\sqrt{\frac{\lambda}{m}}\left|\sin\left(\frac{ka}{2}\right)\right|\)

A Diatomic Chain

同样可以列出运动方程，并求解出色散关系： \(\omega^2_\pm=\frac{\lambda}{nM}\left[m+M\pm\sqrt{(m-M)^2+4mM\cos^2(ka)}\right]\)

这里如果令不同原子间的\(\lambda\)值不同，或m不同，都会出现两个不同的振动模式，利用群速度可以分辨： \(\frac{\partial \omega}{\partial k}=\frac{2k^2a\sin ka}{4k(m+M)\omega-4mM\omega^3}\) ,对 $\omega$和$k$取0极限，得到一个常数， $k$取极限理解为空间周期为无穷，即解是一个常数，晶格做整体运动，这个常数 $\frac{\partial \omega}{\partial k}\bigg|_{w,k=0}=v=\sqrt{ka^2/(m+M)}$，这和声音传播的速度的dependence是一样的。因此散射图里更低的分支被称为声学（acoustic）分支。另一个分支，当k=0但 $\omega$不等于零时，对应一个等于0的群速度和一个不为零的振动频率$\omega$，因此这是在一个原胞里的振动，被称作光学（optical）分支，这个分支通常与电磁波互相干扰。

Genearlized Coupled Oscillator Problem：

用平衡位置时得拉格朗日量描述： \(\frac{\partial L}{\partial q_{i}}-\frac{d}{dt}\frac{\partial L}{\partial \dot{q}}=0 \\ \frac{\partial L}{\partial q_{i}}\bigg |_0=0=\frac{\partial T}{\partial q_{i}}\bigg |_0-\frac{\partial U}{\partial q_{i}}\bigg |_0 \\ \frac{\partial T}{\partial q_{i}}=\frac{\partial}{\partial q_{i}}\left(\frac{1}{2}\sum_{j,k}m_{jk}\dot q_j\dot q_k\right)=0\) 因此： \(\frac{\partial L}{\partial q_i}=-\frac{\partial}{\partial q_i}U(q_1,q_2,...)\approx -\frac{\partial}{\partial q_j}\sum_{j,k}\left(\frac{\partial^2U}{\partial q_j\partial q_k}\bigg |_0\right)q_jq_k=-\frac{1}{2}\sum_{k}A_{ik}q_k\) 又有： \(\frac{\partial T}{\partial \dot q_i}=\frac{\partial}{\partial \dot q_i}\left(\frac{1}{2}\sum_{j,k}m_{jk}\dot q_j\dot q_k\right)=\frac{1}{2}\sum_km_{ik}\dot q_k\) 所以运动方程可以写作： \(\frac{\partial L}{\partial q_i}-\frac{d}{dt}\frac{\partial L}{\partial \dot q}=-\frac{1}{2}\sum_k\left(A_{ik}q_k+m_{ik}\ddot q_k\right)=0\) 带入$q_i(t)=a_ie^{i(wt-\delta)}$可以解出n个频率$w$.

Definition of Stress and Strain The generalization of Hooke’s law \(F=-Kx\), for a 3D continuous medium, is: \(\sigma_{ij}=C_{ijlm}\epsilon_{lm}\) $\tilde\sigma$is a 3x3 matrix called stress tensor, corresponding to $\vec F=\tilde\sigma\cdot \hat nA$. The displacement in responding to stress is also a 3x3 matrix\(\tilde \epsilon\)called strain tensor. $\widetilde C$ is a 3x3x3x3 double tensor consisting of crystal elastic constant Any uniaxial strain can be decomposed into a hydrostatic strain (has equal term along the diagonal) and shear strain (has zero trace).

2. Computing Phonon Spectrum

Expanding Hamiltonian $H=T+V=\sum_i\frac{1}{2}m_I\dot{u_I}^2+U(u_0,u_1,u_2,..)$ at equilibrium position: \(H\approx E_0+\frac{1}{2}\sum_{I,J}\frac{\partial^2U}{\partial u_I\partial u_J}u_I,u_J\) since we know: $\dot{q}=-\frac{\partial H}{\partial u}$, so: \(m\ddot{u}_I=\sum_{J}\frac{\partial^2U}{\partial u_I\partial u_J}u_J\) We again using one dimensional system, but the system here is infinite large. To include the 3d cases, one should add index $\alpha, \beta \in{x,y,z}$.

Again we can set \(u_I=\frac{C_I}{\sqrt{m_I}}e^{-i\omega t}\) and applied the periodic boundary condition $u_{i+R}=u_i$ with period $R$ to make it looks more like a solid system, this will gives: \(u_{iv}=\frac{C_i(k)}{\sqrt{m_i}}e^{-i(\omega t+kR_v)}\) where $i$ is in the primitive cell, and $kR_v=2\pi$, bring it to the equation of motion we have:

\(\omega^2\sqrt{m_i}C_i(k)e^{-i(\omega t+kR_v)}=\sum_j\frac{1}{\sqrt{m_j}}\frac{\partial^2U}{\partial u_I\partial u_J}\sum_\mu C_i(k)e^{-i(wt+kR_\mu)}\) here sum over $J$ becomes sum over all repetitive unit by lattice vector $R_v$ and atom in unit cell, i.e. by summing over $j,v$. Rewrite it we have: \(\begin{align} \omega^2 C_i(k)&=\sum_{j}\left[\sum_\mu\frac{1}{\sqrt{m_im_j}}\frac{\partial^2U}{\partial u_{iv}\partial u_{j\mu}}e^{-ik(R_\mu-R_v)}\right]C_i(k)\\ &=\sum_j \mathcal{D}_{ij}C_i(k) \end{align}\)

Dynamic matrix of solid is \(\mathcal{D}_{ij}=\sum_\mu\frac{1}{\sqrt{m_im_j}}\frac{\partial^2U}{\partial u_{iv}\partial u_{j\mu}}e^{-ik(R_\mu-R_v)}\) (note: think about why $\mathcal{D}$ doesnot dependent on $v$).

Computing $\mathcal{D_{ij}}$(direct method)

This is the same procedure as the k-point sampling in electronic structures:

We can notice that $\mathcal{D_{ij}}$ has its period, depends on $k$. So we should write $\mathcal{D_{ij}}(k+G)=\mathcal{D_{ij}}(k)$ for $G(R_\mu-R_v)=2\pi$. So this can be write as a fourier integral, where $\frac{1}{\sqrt{m_im_j}}\frac{\partial^2U}{\partial u_{iv}\partial u_{j\mu}}$ is the fourier spectrum: \(\frac{1}{\sqrt{m_im_j}}\frac{\partial^2U}{\partial u_{iv}\partial u_{j\mu}}=\frac{L}{2\pi}\int dk \mathcal{D_{ij}}(k)e^{ik(R_\mu-R_v)}\)

We can set a grid of $k$ and evaluate the integral, and for each $k$, only a specific set of $R_\mu$.

The only remaining problem is how to evaluate $\mathcal{D_{ij}}$ at large distance atom, i.e. $\frac{\partial^2U}{\partial u_{iv}\partial u_{j\mu}}$ when $iv,j\mu$ is far away from each other: Here, two assumption is used:

• first, this them decay for long term atom. this brings throuble when long-term interaction, such as in polar insulator, or medal.

• second, the small displacement in other supercell doesnot interupt the one you are considering. Since if you doing pbc, once a atom is perturbed in unit cell, the corresponding atom is also displaced from their eq position.

Then \(\frac{\partial^2U}{\partial u_{iv}\partial u_{j\mu}}=\frac{\partial F_{iv}}{\partial u_{jv}}\approx\frac{F_{iv}(\delta u_{jv})}{\delta u_{jv}}\)

It feels like a little self contradictory here, but I haven’t figure out where. (This is vital).

Density Functional Perturbation Theory

二. Quantum Waves

量子单谐振子： 对于固体中的两个原子，建模成谐振子后，哈密顿量可以写成： \(H=\frac{p^2}{2M}+\frac{1}{2}Kx^2=\frac{p^2}{2M}+\frac{1}{2}M\omega_0^2x^2\) 定义： \(\hat{x}=\sqrt{\frac{M\omega_0}{h}}x, \quad \hat{p}=\sqrt{\frac{1}{Mh\omega_0}p}\) 带入H, 得到： \(H=\frac{1}{2}(\hat{x}^2+\hat{p}^2)h\omega_0\) 定义 \(a=\frac{1}{\sqrt{2}}(\hat{x}+i\hat{p}), \quad a^\dagger=\frac{1}{\sqrt{2}}(\hat{x}-i\hat{p})\) \(\left[a,a^{\dagger}\right]=1, \quad H=\left(a^{\dagger}a+\frac{1}{2}\right)h\omega_0\) Its eigenstates have energies:\(E=(N+\frac{1}{2})h\omega_0\)， and

Stress and Strain Definition: 声子（Phonon）： 依然用一个一维原子链的模型，哈密顿量写作： \(H=\sum_n{\frac{p_n^2}{2m}}+\frac{\lambda}{2}\sum_n(x_n-x_{n-1})^2\) (1) 用傅里叶变换解耦 \(\sum_n(x_n-x_{n-1})^2=\sum_n(x_n^2+x_{n-1}^2-2x_nx_{n-1})\\=\sum_n(2x_n^2-2x_nx_{n-1})\) parseval‘s identity: \(\sum_nx_n^2=\sum_k\hat{x}_k^2\) and \(x_{n\pm1}=(\mathcal{F}^{-1}\circ\mathcal{F})[ x_{n\pm1}]=\mathcal{F}^{-1}[x_ke^{\pm i ka}]\) \(\sum_nx_nx_{x\pm1}=\sum_n\left[\left(\sum_kx_ke^{ikan}\right)^\dagger\left(\sum_{k'}x_{k'}e^{ik'a(n\pm1)}\right)\right]\\=\sum_{k,k',n}x_k^\dagger x_{k'}e^{i(k'-k)an}e^{\pm ika}=\sum_k{|x_k|^2e^{\pm ika}}\) 因此哈密顿量可以写成： \(H=\sum_k\left(\frac{|p_k|^2}{2M}+\frac{1}{2}M\omega_k^2|x_k|^2\right)\) 量子化后变为： \(H=\sum_k\left(a^\dagger_k a_k+\frac{1}{2}\right)h\omega_k\) 可以通过傅里叶逆变换将表示成产生和湮灭算符的k空间displacement算符 \(\hat{x}_k=\frac{1}{\sqrt{x}}(a_k+a_{-k}^\dagger), \quad \hat{p}_k=\frac{-i}{\sqrt{2}}(a_k-a_{-k}^\dagger)\) 变换为实空间对于displacement算符（Local amplitude operators）： \(x(r)=\sum_k\sqrt{\frac{h}{2\rho V\omega_k}}\left(a_ke^{ikr}+a_k^\dagger e^{-ikr}\right)\) \(\dot x(r)=-i\sum_k\sqrt{\frac{h\omega_{k}}{2\rho V}}\left(a_ke^{ikr}-a_k^\dagger e^{-ikr}\right)\) In 3D cases, a polarization vector $\eta_{\vec k\lambda}$ is inserted, with: \(\vec x(\vec r)=\sum_{\vec k,\lambda}\hat \eta_{\vec k\lambda}\sqrt{\frac{h}{2\rho V\omega_{\vec k\lambda}}}\left(a_{\vec k\lambda}e^{i\vec k\cdot \vec r}+a_{\vec k\lambda}^\dagger e^{-i\vec k \cdot \vec r}\right)\) \(\dot{\vec x}(\vec r)=-i\sum_{\vec k,\lambda}\hat\eta_{\vec k\lambda}\sqrt{\frac{h\omega_{\vec k\lambda}}{2\rho V}}\left(a_{\vec k\lambda}e^{i\vec k\cdot \vec r}-a_{\vec k\lambda}^\dagger e^{-i\vec k \cdot \vec r}\right)\)

Coherence State

Spatial Field Operator For a 3D system, a field operator is the Fourier inversion transform of create and annihilation opeartor: \(\psi^\dagger(\vec r)=\sum_{\vec k}\frac{e^{-i\vec k\cdot \vec r}}{\sqrt{L}}a^{\dagger}_{\vec k}, \quad \psi(\vec r)=\sum_{\vec k}\frac{e^{i\vec k\cdot \vec r}}{\sqrt{L}}a_{\vec k}\) Which creates and destroys a particle at exactly the point r. Single and many particle with wave function now is represented as: \(|\phi\rangle=\int{d^3r\phi(\vec r)\psi^\dagger(\vec r)|0\rangle}, \quad |\psi\rangle=\prod_i\frac{1}{\sqrt{N_i!}}\left(\int{d^3r_i\phi_i(\vec r)}\psi^\dagger(\vec r_i)\right)^{N_i}|0\rangle\) they satisfied normalization condition. Electron Fermi Field Operators

三. 声子电子相互作用 Deformation Potential Scattering Deformation Potential is a parameter that can be measured by the shift of the electron bands under hydrostatic pressure. Strain: \(\epsilon_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\) \(\vec u(\vec r)=\sum_{\vec k,\lambda}\hat \eta_{\vec k\lambda}\sqrt{\frac{h}{2\rho V\omega_{\vec k\lambda}}}\left(a_{\vec k\lambda}e^{i\vec k\cdot \vec r}+a_{\vec k\lambda}^\dagger e^{-i\vec k \cdot \vec r}\right)\) Then we have hamiltonian: \(H_{int}=DTr(\tilde\epsilon)=D\left(\epsilon_{xx}+\epsilon_{yy}+\epsilon_{zz}\right)\\=D\left(\frac{\partial u_x}{\partial x}+\frac{\partial u_y}{\partial y}+\frac{\partial u_z}{\partial z}\right)\) Integrated with specific spatial field operator, that: \(H_D=\int d^3x\Psi_n^\dagger(\vec x)H_{int}(\vec x)\Psi_n(\vec x)\\=\sum_{\vec k.\vec k_1}Dk\sqrt{\frac{h}{2\rho V\omega_k}}i\left(a_{\vec k}b^\dagger_{n,\vec k_1}b_{n,\vec k_1}-a^\dagger_{\vec k}b^\dagger_{n,\vec k_1-\vec k}b_{n,\vec k_1}\right)\) This expression can be generalized to account for shear deformation, replacing D with operator acting on band states, we have a new Hamiltonian: \(H_{n,n'}=\sum_{\vec k,\lambda,\vec k_1}\langle n'|D|n\rangle k\sqrt{\frac{h}{2\rho V\omega_k}}i\left(a_{\vec k\lambda}b^\dagger_{n,\vec k_1}b_{n,\vec k_1}-a^\dagger_{\vec k \lambda}b^\dagger_{n,\vec k_1-\vec k}b_{n,\vec k_1}\right)\) When\(\langle n'|D|n\rangle=D\delta_{n,n'}\), it represents a hydrostatic strain, but the off-diagonal elements are for shear strains.

用电子声子相互作用矩阵$M$，电声相互作用哈密顿量可以写成：

\(\hat{H}_{el-ph}=\sum_{k,v,k',v',\lambda}M_{k,v,k',v',\lambda}a_{k',v'}^\dagger a_{k',v'}(b_{q,\mu,\lambda}+b_{-q,-\mu,\lambda}^\dagger)\) 2104.08998.pdf (arxiv.org) 里，M可以由DFPT计算 \(M^v_{mn}(k,q)=\langle\psi_{mk+q}(r)|\delta_{vq}V(r)|\psi_{nq}(r)\rangle\)

四. 声子自能

Pourfath, Mahdi. The non-equilibrium Green’s function method for nanoscale device simulation . Vol. 3. Vienna: Springer, 2014.

分为三种情况考虑，Acoustic，Optical，和Polar Optical。前两种只有相同原子时，自能才有元素，所以不会影响哈密顿量的稀疏性质（分块三对角）。但是Polar Optical和声子波矢有关，推导时会产生非对角项，让哈密顿矩阵变得稠密，只能直接求逆。这是造成含声子NEGF计算开销大很重要的原因。

计算时，先计算出lesser/greater self-energy, 接着retarted self-enrgy可以近似为\(\Sigma_{i_1i_2}^r(k_t,E)\approx -\frac{i}{2}\Gamma_{i_1,i_2}(k_t,E), \quad \Gamma_{i_1i_2}(k_t,E)=2Im[\Sigma_{i_1i_2}^>(k_t,E)]\) Acoustic Phonon

用electron-phonon interation matrix 写：

Optical Phonon

用electron-phonon interation matrix 写：

Polar Optical Phonon 这里就包含非对角项的元素了：

用effective mass model，上面的积分可以等价于：

其中 $\rho_{2D}$是二维的态密度，

\(\epsilon_{q_t}=h^2q_t^2/(2m^*)\) \(\rho_{2D}=m^*/(2\pi h^2)\)

Anantram, M. P., Mark S. Lundstrom, and Dmitri E. Nikonov. “Modeling of nanoscale devices.” Proceedings of the IEEE 96.9 (2008): 1511-1550. 应该只考虑了对角元的声子自能，输入是不同mode\(\eta\)的声子的phonon deformation potential\(D^\eta\)，计算方法为： \(\Sigma^{in}_{q,q}(E)=\sum_{\eta}D_q^\eta \left[n_B(h\omega_{phonon})G^n_{q,q}(E-h\omega_{phonon})+(n_B(h\omega_{phonon})+1)G_{q,q}^n(E+h\omega_{phonon})\right].\)

Direct Method for Calculating Temperature-Dependent Transport Properties 先用“Force Constant Approach”计算声子谱，后将所有phonon modes求和，得到原子核做扰动的位移函数$u_l$： \(u_l(T,t)=\frac{1}{\sqrt{N_qM_l}}\sum_{sq}\epsilon_{sq}A_{sq}(T)e^{i(q\cdot R_l-\omega_{sq}t+\phi_{sq})}\) $N_q$是散射区域的横向超胞的波矢数，$M_l$是第l个原子的质量，$\phi_{sq}$是随机相位，用来模拟thermal disorder。$A_{sq}$是扰动振幅，$\omega_{sq}^2A_{sq}^2/2$应该等于sq phonon mode的总能量。