Derive phonon scattering rate from anharmonic potential

Hamiltonian of a crystal

Like what we did in the previous post, we would like to do the same thing to the crystal, but in a quantum way, and also including the anharmonic potential, i.e. the third-order term in the Taylor expansion of the interatomic potential energy.

Starting from the interatomic potential, which is always expressed by the IFCs expressed by \(\Phi_{\alpha_1\alpha_2\dots\alpha_n}(\boldsymbol{l}_1,\boldsymbol{b}_1;\boldsymbol{l}_2,\boldsymbol{b}_2,\dots,\boldsymbol{l}_n,\boldsymbol{b}_n)\), with \(\alpha_i\) denotes direction, \(\boldsymbol{l}_i\) denotes the lattice site, and \(\boldsymbol{b}_i\) denotes the basis atom. The interatomic potental is: \[ \begin{aligned} \mathcal{V} = & \mathcal{V}_0 + \cancel{\mathcal{V}_1} + \underbrace{\frac{1}{2}\sum_{\boldsymbol{l},\boldsymbol{b},\boldsymbol{l}',\boldsymbol{b}'}\sum_{\alpha,\beta} \Phi_{\alpha,\beta}(\boldsymbol{l}\boldsymbol{b};\boldsymbol{l}'\boldsymbol{b}')u_{\alpha}(\boldsymbol{l}\boldsymbol{b})u_{\beta}(\boldsymbol{l}'\boldsymbol{b}')}_{\mathcal{V}_2}\\ & + \underbrace{ \frac{1}{3!}\sum_{\boldsymbol{l},\boldsymbol{b},\boldsymbol{l}',\boldsymbol{b}',\boldsymbol{l}'',\boldsymbol{b}''}\sum_{\alpha,\beta,\gamma} \Phi_{\alpha,\beta,\gamma}(\boldsymbol{l}\boldsymbol{b};\boldsymbol{l}'\boldsymbol{b}';\boldsymbol{l}''\boldsymbol{b}'')u_{\alpha}(\boldsymbol{l}\boldsymbol{b})u_{\beta}(\boldsymbol{l}'\boldsymbol{b}')u_{\gamma}(\boldsymbol{l}''\boldsymbol{b}'') }_{\mathcal{V}_3} + \cdots \end{aligned} \]

Here, we use \({u}_{\alpha}(\boldsymbol{l}\boldsymbol{b})\) to note the spatial deviation of atom \(\boldsymbol{b}\) in cell \(\boldsymbol{l}\) along \(\alpha\) direction, relative to the equilibrium position. Until here, everyone including Newton is happy, as all we wrote is classical, and all the terms are just numbers.

Classical approach is useful, as it surely can depict the oscillations in a crystal. But it is not enough to reveal the concept of phonon, which is quantized by quantum mechanics. So we have to turn to the quantum world, where everything is controlled by the Schrödinger equation, and the Hamiltonian is the key to it.

Recall the (time-independent) many-body Schrödinger equation:

\[ \hat{H} \Psi(\boldsymbol{r}_1,\boldsymbol{r}_2,\dots,\boldsymbol{r}_N) = E \Psi(\boldsymbol{r}_1,\boldsymbol{r}_2,\dots,\boldsymbol{r}_N) \]

where \(\hat{H} = \sum_{i=1}^N \frac{\hat{\boldsymbol{p}}_i^2}{2m_i} + \hat{\mathcal{V}}(\boldsymbol{r}_1,\boldsymbol{r}_2,\dots,\boldsymbol{r}_N)\) is the Hamiltonian of the system, and \(\Psi(\boldsymbol{r}_1,\boldsymbol{r}_2,\dots,\boldsymbol{r}_N)\) is the wavefunction, \(E\) is the total energy. Generally we love the following simplification to the wavefunction:

\[ \Psi(\boldsymbol{r}_1,\boldsymbol{r}_2,\dots,\boldsymbol{r}_N) = \prod_{i=1}^N \psi(\boldsymbol{r}_i) \]

that allows us to rewrite the many-body Schrödinger equation as:

\[ \hat{H}_i \psi(\boldsymbol{r}_i) = E_i \psi(\boldsymbol{r}_i), \quad\forall i=1,2,\dots,N \]

with \(\hat{H}_i = \frac{\hat{\boldsymbol{p}}_i}{2 m_i} + \hat{\mathcal{V}}\).

Using the potential energy we wrote, the Hamiltonian of a crystal is thus:

\[ \begin{aligned} \hat{H} = & \sum_{\boldsymbol{l}\boldsymbol{b}} \frac{\hat{p}^2 (\boldsymbol{l}\boldsymbol{b})}{2m_{\boldsymbol{b}}}\\ & + \underbrace{\frac{1}{2} \sum_{\boldsymbol{l}\boldsymbol{b}\boldsymbol{l}'\boldsymbol{b}'} \sum_{\alpha\beta} \Phi_{\alpha\beta}(\boldsymbol{l}\boldsymbol{b};\boldsymbol{l}'\boldsymbol{b}') \hat{u}_{\alpha}(\boldsymbol{l}\boldsymbol{b}) \hat{u}_{\beta}(\boldsymbol{l}'\boldsymbol{b}')}_{\hat{\mathcal{V}}_2}\\ & + \underbrace{\frac{1}{3!} \sum_{\boldsymbol{l}\boldsymbol{b}\boldsymbol{l}'\boldsymbol{b}'\boldsymbol{l}''\boldsymbol{b}''} \sum_{\alpha\beta\gamma} \Phi_{\alpha\beta\gamma}(\boldsymbol{l}\boldsymbol{b};\boldsymbol{l}'\boldsymbol{b}';\boldsymbol{l}''\boldsymbol{b}'') \hat{u}_{\alpha}(\boldsymbol{l}\boldsymbol{b}) \hat{u}_{\beta}(\boldsymbol{l}'\boldsymbol{b}') \hat{u}_{\gamma}(\boldsymbol{l}''\boldsymbol{b}'')}_{\hat{\mathcal{V}}_3} + \dots\\ \end{aligned} \]

Note that the deviations \(\hat{u}_{\alpha}(\boldsymbol{l}\boldsymbol{b})\), the momentum \(\hat{\boldsymbol{p}}(\boldsymbol{l}\boldsymbol{b})\) are operators rather than numbers. Their commutation relations are:

\[ [\hat{u}_{\alpha}(\boldsymbol{l}\boldsymbol{b}), \hat{p}_{\beta}(\boldsymbol{l}'\boldsymbol{b}')] = i\hbar\delta_{\alpha\beta}\delta_{\boldsymbol{l}\boldsymbol{l}'}\delta_{\boldsymbol{b}\boldsymbol{b}'} \]

where \(\hat{\boldsymbol{u}}\) is the vector form of the deviation operators.

Diagonal representation of the Hamiltonian

\((\boldsymbol{q},\boldsymbol{b})\) representation

The Hamiltonian we wrote is not diagonal, which means that the eigenstates of the Hamiltonian are not the same as the basis we used. We need to transform the Hamiltonian into a diagonal form, so that the eigenstates of the Hamiltonian are the same as the basis we used.

Due to that we are dealing with a cyclic system, that we use periodic boundary condition for \(N_0 = N_1\times N_2\times N_3\) unit cells, the operators \(\hat{\boldsymbol{u}}(\boldsymbol{l}\boldsymbol{b})\) and \(\hat{\boldsymbol{p}}(\boldsymbol{l}\boldsymbol{b})\) can be represented in the Fourier's fashion:

\[ \begin{aligned} \hat{\boldsymbol{u}}(\boldsymbol{l}\boldsymbol{b}) &= \frac{1}{\sqrt{N_0}} \sum_{\boldsymbol{q}} \hat{\boldsymbol{X}} (\boldsymbol{q};\boldsymbol{b}) \exp[i\boldsymbol{q}\cdot\boldsymbol{l}]\\ \hat{\boldsymbol{p}}(\boldsymbol{l}\boldsymbol{b}) &= \frac{1}{\sqrt{N_0}} \sum_{\boldsymbol{q}} \hat{\boldsymbol{P}} (\boldsymbol{q};\boldsymbol{b}) \exp[-i\boldsymbol{q}\cdot\boldsymbol{l}] \end{aligned} \]

where \(\boldsymbol{q}\) is the wavevector in the first Brillouin zone. The inverse of the above equations are:

\[ \begin{aligned} \hat{\boldsymbol{X}} (\boldsymbol{q};\boldsymbol{b}) &= \frac{1}{\sqrt{N_0}} \sum_{\boldsymbol{l}} \hat{\boldsymbol{u}}(\boldsymbol{l}\boldsymbol{b}) \exp[-i\boldsymbol{q}\cdot\boldsymbol{l}]\\ \hat{\boldsymbol{P}} (\boldsymbol{q};\boldsymbol{b}) &= \frac{1}{\sqrt{N_0}} \sum_{\boldsymbol{l}} \hat{\boldsymbol{p}}(\boldsymbol{l}\boldsymbol{b}) \exp[i\boldsymbol{q}\cdot\boldsymbol{l}] \end{aligned} \]

Note the following properties:

\[ \begin{aligned} \hat{\boldsymbol{X}}^\dagger (\boldsymbol{q};\boldsymbol{b}) &= \hat{\boldsymbol{X}} (-\boldsymbol{q};\boldsymbol{b})\\ \hat{\boldsymbol{P}}^\dagger (\boldsymbol{q};\boldsymbol{b}) &= \hat{\boldsymbol{P}} (-\boldsymbol{q};\boldsymbol{b}) \end{aligned} \]

You may ask, that why the hell we use Fourier transform to obtain \(\hat{\boldsymbol{X}}\) and inverse Fourier transform to obtain \(\hat{\boldsymbol{P}}\). The answer is that, we want to keep the canonical commutation relation between those two. Using the relation that \([\hat{\boldsymbol{u}}(\boldsymbol{l}\boldsymbol{b}), \hat{\boldsymbol{p}}(\boldsymbol{l}'\boldsymbol{b}')] = i\hbar\delta_{\boldsymbol{l}\boldsymbol{l}'}\delta_{\boldsymbol{b}\boldsymbol{b}'}\), we know that:

\[ \begin{aligned} \left[\hat{X}_{\alpha} (\boldsymbol{q};\boldsymbol{b}), \hat{P}_{\beta} (\boldsymbol{q'};\boldsymbol{b'})\right] &= \frac{1}{N_0} \sum_{\boldsymbol{l}} \sum_{\boldsymbol{l}'} \left[\hat{u}_{\alpha}(\boldsymbol{l}\boldsymbol{b}), \hat{p}_{\beta}(\boldsymbol{l}'\boldsymbol{b}')\right] \exp[-i(\boldsymbol{q}\cdot\boldsymbol{l}-\boldsymbol{q}'\cdot\boldsymbol{l}')] \\ &= \frac{1}{N_0} \sum_{\boldsymbol{l}} \sum_{\boldsymbol{l}'}i\hbar \delta_{\alpha\beta}\delta_{\boldsymbol{l}\boldsymbol{l}'}\delta_{\boldsymbol{b}\boldsymbol{b}'}\exp[-i(\boldsymbol{q}\cdot\boldsymbol{l}-\boldsymbol{q}'\cdot\boldsymbol{l}')]\\ &= \frac{1}{N_0} \delta_{\alpha\beta}\sum_{\boldsymbol{l}}i\hbar \delta_{\boldsymbol{b}\boldsymbol{b}'}\exp[-i(\boldsymbol{q}-\boldsymbol{q}')\cdot\boldsymbol{l}]\\ &= i\hbar \delta_{\alpha\beta}\delta_{\boldsymbol{q}\boldsymbol{q}'}\delta_{\boldsymbol{b}\boldsymbol{b}'} \end{aligned} \]

Also, the reason why we multiply \(1/\sqrt{N_0}\) in the transformation is that, we want to obtain the above canonical commutation relation.

With the above equations, the Hamiltonian can be rewritten as:

\[ \begin{aligned} \text{1st term} &= \sum_{\boldsymbol{l}\boldsymbol{b}} \frac{\hat{p}^2 (\boldsymbol{l}\boldsymbol{b})}{2m_{\boldsymbol{b}}}\\ &= \frac{1}{N_0}\sum_{\boldsymbol{l}\boldsymbol{b}} \sum_{\boldsymbol{q}\boldsymbol{q}'} \frac{1}{2m_{\boldsymbol{b}}} \hat{\boldsymbol{P}} (\boldsymbol{q};\boldsymbol{b}) \hat{\boldsymbol{P}}(\boldsymbol{q}';\boldsymbol{b}) \exp[-i(\boldsymbol{q}+\boldsymbol{q}')\cdot\boldsymbol{l}]\\ &= \sum_{\boldsymbol{b}} \sum_{\boldsymbol{q}\boldsymbol{q}'} \frac{1}{2m_{\boldsymbol{b}}} \hat{\boldsymbol{P}} (\boldsymbol{q};\boldsymbol{b}) \hat{\boldsymbol{P}}(\boldsymbol{q}';\boldsymbol{b}) \delta_{\boldsymbol{0},\boldsymbol{q}+\boldsymbol{q}'}\\ &= \sum_{\boldsymbol{q}\boldsymbol{b}}\frac{1}{2m_{\boldsymbol{b}}} \hat{\boldsymbol{P}} (\boldsymbol{q};\boldsymbol{b}) \hat{\boldsymbol{P}}^\dagger (\boldsymbol{q};\boldsymbol{b})\\ \text{2nd term} &= \frac{1}{2} \sum_{\boldsymbol{l}\boldsymbol{b}\boldsymbol{l}'\boldsymbol{b}'} \sum_{\alpha\beta} \Phi_{\alpha\beta}(\boldsymbol{l}\boldsymbol{b};\boldsymbol{l}'\boldsymbol{b}') \hat{u}_{\alpha}(\boldsymbol{l}\boldsymbol{b}) \hat{u}_{\beta}(\boldsymbol{l}'\boldsymbol{b}')\\ &= \frac{1}{2}\frac{1}{N_0} \sum_{\boldsymbol{l}\boldsymbol{b}\boldsymbol{l}'\boldsymbol{b}'} \sum_{\alpha\beta} \sum_{\boldsymbol{q}\boldsymbol{q}'}\Phi_{\alpha\beta}(\boldsymbol{l}\boldsymbol{b};\boldsymbol{l}'\boldsymbol{b}') \hat{X}_{\alpha} (\boldsymbol{q};\boldsymbol{b}) \hat{X}_{\beta}(\boldsymbol{q}';\boldsymbol{b}') \exp[i\boldsymbol{q}\cdot\boldsymbol{l} + i\boldsymbol{q}'\cdot\boldsymbol{l}']\\ &= \frac{1}{2} \frac{1}{N_0}\sum_{\boldsymbol{b}\boldsymbol{b}'} \sum_{\boldsymbol{l}\boldsymbol{h}} \sum_{\alpha\beta} \sum_{\boldsymbol{q}\boldsymbol{q}'} \Phi_{\alpha\beta}(\boldsymbol{0}\boldsymbol{b};\boldsymbol{h}\boldsymbol{b}') \hat{X}_{\alpha}(\boldsymbol{q};\boldsymbol{b}) \hat{X}_{\beta}(\boldsymbol{q}';\boldsymbol{b}') \exp[i(\boldsymbol{q}+\boldsymbol{q}')\cdot\boldsymbol{l}]\exp[i\boldsymbol{q}'\cdot\boldsymbol{h}]\\ &= \frac{1}{2} \sum_{\alpha\beta} \sum_{\boldsymbol{q}} \sum_{\boldsymbol{b}\boldsymbol{b}'} \underbrace{\sum_{\boldsymbol{h}} \Phi_{\alpha\beta}(\boldsymbol{0}\boldsymbol{b};\boldsymbol{h}\boldsymbol{b}') \exp[-i\boldsymbol{q}\cdot\boldsymbol{h}]}_{\equiv \Phi_{\alpha\beta}(\boldsymbol{b}\boldsymbol{b}'|\boldsymbol{q})} \hat{X}_{\alpha}(\boldsymbol{q};\boldsymbol{b}) \hat{X}_{\beta}^\dagger (\boldsymbol{q};\boldsymbol{b}') \\ &= \frac{1}{2} \sum_{\alpha\beta} \sum_{\boldsymbol{q}\boldsymbol{b}\boldsymbol{b}'} \Phi_{\alpha\beta}(\boldsymbol{b}\boldsymbol{b}'|\boldsymbol{q}) \hat{X}_{\alpha}(\boldsymbol{q};\boldsymbol{b}) \hat{X}_{\beta}^\dagger (\boldsymbol{q};\boldsymbol{b}')\\ \text{3rd term} &= \frac{1}{3!} \sum_{\boldsymbol{l}\boldsymbol{b}\boldsymbol{l}'\boldsymbol{b}'\boldsymbol{l}''\boldsymbol{b}''} \sum_{\alpha\beta\gamma} \Phi_{\alpha\beta\gamma}(\boldsymbol{l}\boldsymbol{b};\boldsymbol{l}'\boldsymbol{b}';\boldsymbol{l}''\boldsymbol{b}'') \hat{u}_{\alpha}(\boldsymbol{l}\boldsymbol{b}) \hat{u}_{\beta}(\boldsymbol{l}'\boldsymbol{b}') \hat{u}_{\gamma}(\boldsymbol{l}''\boldsymbol{b}'')\\ &= \frac{1}{3!}\frac{1}{N_0^{\frac{3}{2}}} \sum_{\boldsymbol{l}\boldsymbol{b}\boldsymbol{l}'\boldsymbol{b}'\boldsymbol{l}''\boldsymbol{b}''} \sum_{\alpha\beta\gamma} \sum_{\boldsymbol{q}\boldsymbol{q}'\boldsymbol{q}''} \Phi_{\alpha\beta\gamma}(\boldsymbol{l}\boldsymbol{b};\boldsymbol{l}'\boldsymbol{b}';\boldsymbol{l}''\boldsymbol{b}'') \hat{X}_{\alpha}(\boldsymbol{q};\boldsymbol{b}) \hat{X}_{\beta}(\boldsymbol{q}';\boldsymbol{b}') \hat{X}_{\gamma}(\boldsymbol{q}'';\boldsymbol{b}'') \exp[i\boldsymbol{q}\cdot\boldsymbol{l}+i\boldsymbol{q}'\cdot\boldsymbol{l}'+i\boldsymbol{q}''\cdot\boldsymbol{l}'']\\ &= \frac{1}{3!} \frac{1}{N_0^{\frac{3}{2}}} \sum_{\boldsymbol{b}\boldsymbol{b}'\boldsymbol{b}''}\sum_{\boldsymbol{l}\boldsymbol{h}\boldsymbol{h}'} \sum_{\alpha\beta\gamma}\sum_{\boldsymbol{q}\boldsymbol{q}'\boldsymbol{q}''} \Phi_{\alpha\beta\gamma}(\boldsymbol{0}\boldsymbol{b};\boldsymbol{h}\boldsymbol{b}';\boldsymbol{h}'\boldsymbol{b}'') \hat{X}_{\alpha}(\boldsymbol{q};\boldsymbol{b}) \hat{X}_{\beta}(\boldsymbol{q}';\boldsymbol{b}') \hat{X}_{\gamma}(\boldsymbol{q}'';\boldsymbol{b}'') \exp[i(\boldsymbol{q}+\boldsymbol{q}'+\boldsymbol{q}'')\cdot\boldsymbol{l}]\exp[i\boldsymbol{q}'\cdot\boldsymbol{h}]\exp[i\boldsymbol{q}''\cdot\boldsymbol{h}']\\ &= \frac{1}{3!} \frac{1}{\sqrt{N_0}} \sum_{\alpha\beta\gamma}\sum_{\boldsymbol{q}\boldsymbol{q}'\boldsymbol{q}''\boldsymbol{b}\boldsymbol{b}'\boldsymbol{b}''} \underbrace{\sum_{\boldsymbol{h}\boldsymbol{h}'} \Phi_{\alpha\beta\gamma}(\boldsymbol{0}\boldsymbol{b};\boldsymbol{h}\boldsymbol{b}';\boldsymbol{h}'\boldsymbol{b}'') \exp[i\boldsymbol{q}'\cdot\boldsymbol{h}]\exp[i\boldsymbol{q}''\cdot\boldsymbol{h}']}_{\equiv \Phi_{\alpha\beta\gamma}(\boldsymbol{q}\boldsymbol{b}, \boldsymbol{q}'\boldsymbol{b}', \boldsymbol{q}''\boldsymbol{b}'') }\hat{X}_{\alpha}(\boldsymbol{q};\boldsymbol{b}) \hat{X}_{\beta}(\boldsymbol{q}';\boldsymbol{b}') \hat{X}_{\gamma}(\boldsymbol{q}'';\boldsymbol{b}'') \delta_{\boldsymbol{G},\boldsymbol{q}+\boldsymbol{q}'+\boldsymbol{q}''}\\ &= \frac{1}{3!} \frac{1}{\sqrt{N_0}}\sum_{\alpha\beta\gamma}\sum_{\boldsymbol{q}\boldsymbol{q}'\boldsymbol{q}''\boldsymbol{b}\boldsymbol{b}'\boldsymbol{b}''} \Phi_{\alpha\beta\gamma}(\boldsymbol{q}\boldsymbol{b}, \boldsymbol{q}'\boldsymbol{b}', \boldsymbol{q}''\boldsymbol{b}'') \hat{X}_{\alpha}(\boldsymbol{q};\boldsymbol{b}) \hat{X}_{\beta}(\boldsymbol{q}';\boldsymbol{b}') \hat{X}_{\gamma}(\boldsymbol{q}'';\boldsymbol{b}'') \delta_{\boldsymbol{G},\boldsymbol{q}+\boldsymbol{q}'+\boldsymbol{q}''} \end{aligned} \]

where \(\boldsymbol{G}\) is a reciprocal lattice vector. Using the above derivations, the Hamiltonian is thus expressed as,

\[ \begin{aligned} \hat{H} &= \sum_{\boldsymbol{q}\boldsymbol{b}}\frac{1}{2m_{\boldsymbol{b}}} \hat{\boldsymbol{P}} (\boldsymbol{q};\boldsymbol{b}) \hat{\boldsymbol{P}}^\dagger (\boldsymbol{q};\boldsymbol{b}) + \frac{1}{2} \sum_{\alpha\beta} \sum_{\boldsymbol{q}\boldsymbol{b}\boldsymbol{b}'} \Phi_{\alpha\beta}(\boldsymbol{b}\boldsymbol{b}'|\boldsymbol{q}) \hat{X}_{\alpha}(\boldsymbol{q};\boldsymbol{b}) \hat{X}_{\beta}^\dagger (\boldsymbol{q};\boldsymbol{b}') \\ &\quad + \frac{1}{3!} \frac{1}{\sqrt{N_0}}\sum_{\alpha\beta\gamma}\sum_{\boldsymbol{q}\boldsymbol{q}'\boldsymbol{q}''\boldsymbol{b}\boldsymbol{b}'\boldsymbol{b}''} \Phi_{\alpha\beta\gamma}(\boldsymbol{q}\boldsymbol{b}, \boldsymbol{q}'\boldsymbol{b}', \boldsymbol{q}''\boldsymbol{b}'') \hat{X}_{\alpha}(\boldsymbol{q};\boldsymbol{b}) \hat{X}_{\beta}(\boldsymbol{q}';\boldsymbol{b}') \hat{X}_{\gamma}(\boldsymbol{q}'';\boldsymbol{b}'') \delta_{\boldsymbol{G},\boldsymbol{q}+\boldsymbol{q}'+\boldsymbol{q}''} \end{aligned} \]

\((\boldsymbol{q}, s)\) representation

Now, we have just taken the first step of the Long March.

According to our experiences dealing with classical system (see the previous post that solves the classical harmonic oscillator), there are \(nd\) distinct vibration modes for a crystal with \(n\) atoms in the unit cell, and \(d\) dimensions. And at the same time, the eigenvectors that represent those modes are orthogonal to each other. For the present case, we denote those eigenvectors as \(\boldsymbol{e}(\boldsymbol{b}|\boldsymbol{q}s)\), where \(s\) is phonon branch. It is noteworthy that \(\boldsymbol{e}(\boldsymbol{b}|\boldsymbol{q}s)\) contains the magnitude, direction of vibration, as well as the phase of the vibration. And it is normalized such that,

\[ \begin{aligned} \sum_{\boldsymbol{b}} \boldsymbol{e}^*(\boldsymbol{b}|\boldsymbol{q}s) \boldsymbol{e}(\boldsymbol{b}|\boldsymbol{q}s') &= \delta_{ss'}\\ \sum_s e_{\alpha}^*(\boldsymbol{b}|\boldsymbol{q}s) e_{\beta}(\boldsymbol{b}'|\boldsymbol{q}s) &= \delta_{\alpha\beta}\delta_{\boldsymbol{b}\boldsymbol{b}'} \end{aligned} \]

That is the inner product of two eigenvectors with the same wavevector is 1 if they are the same mode, and 0 otherwise. Taken use of this set of basis, we make another transformation:

\[ \begin{aligned} \hat{X}(\boldsymbol{q}s) &= \sum_{\boldsymbol{b}} \sqrt{m_{\boldsymbol{b}}} \boldsymbol{e}^*(\boldsymbol{b}|\boldsymbol{q}s) \hat{\boldsymbol{X}}(\boldsymbol{q};\boldsymbol{b})\\ \hat{P}(\boldsymbol{q}s) &= \sum_{\boldsymbol{b}} \frac{1}{\sqrt{m_{\boldsymbol{b}}}} \boldsymbol{e}(\boldsymbol{b}|\boldsymbol{q}s) \hat{\boldsymbol{P}}(\boldsymbol{q};\boldsymbol{b}) \end{aligned} \]

Under this transformation, we still have the canonical commutation relation:

\[ \begin{aligned} \left[\hat{X}(\boldsymbol{q}s), \hat{P}(\boldsymbol{q}'s') \right] &= \sum_{\boldsymbol{b}\boldsymbol{b}'} \sum_{\alpha\beta} \sqrt{\frac{m_{\boldsymbol{b}}}{m_{\boldsymbol{b}'}}} e^*_{\alpha}(\boldsymbol{b}|\boldsymbol{q}s)e_{\beta}(\boldsymbol{b}'|\boldsymbol{q}'s')\left[\hat{X}_{\alpha}(\boldsymbol{q};\boldsymbol{b}),\hat{P}_{\beta}(\boldsymbol{q}';\boldsymbol{b}')\right]\\ &= \sum_{\boldsymbol{b}\boldsymbol{b}'} \sum_{\alpha\beta}e^*_{\alpha}(\boldsymbol{b}|\boldsymbol{q}s)e_{\beta}(\boldsymbol{b}'|\boldsymbol{q}'s') i\hbar \delta_{\alpha\beta}\delta_{\boldsymbol{q}\boldsymbol{q}'}\delta_{\boldsymbol{b}\boldsymbol{b}'}\\ &= \sum_{\boldsymbol{b}}\sum_{\alpha} e^*_{\alpha}(\boldsymbol{b}|\boldsymbol{q}s)e_{\alpha}(\boldsymbol{b}|\boldsymbol{q}'s') i\hbar \delta_{\boldsymbol{q}\boldsymbol{q}'}\\ &= i\hbar \delta_{\boldsymbol{q}\boldsymbol{q}'}\delta_{ss'} \end{aligned} \]

Annihilation and creation operators

Using \(\hat{X}(\boldsymbol{q}s)\) and \(\hat{P}(\boldsymbol{q}s)\), we make the final transformation, called by the annihilation and creation operators:

\[ \begin{aligned} \hat{a}_{\boldsymbol{q}s} &= \frac{1}{\sqrt{2\hbar\omega(\boldsymbol{q}s)}}\hat{P}(\boldsymbol{q}s) - i \sqrt{\frac{\omega (\boldsymbol{q}s)}{2\hbar}} \hat{X}^\dagger (\boldsymbol{q}s)\\ \hat{a}^\dagger_{\boldsymbol{q}s} &= \frac{1}{\sqrt{2\hbar\omega(\boldsymbol{q}s)}}\hat{P}^\dagger(\boldsymbol{q}s) + i \sqrt{\frac{\omega (\boldsymbol{q}s)}{2\hbar}} \hat{X} (\boldsymbol{q}s) \end{aligned} \]

You would ask, why the hell they are called annihilation and creation operators? The answer is that, they actually annihilate or create a phonon. To see this, I refer you to the next section.

where \(\omega(\boldsymbol{q}s)\) is the frequency of the phonon mode \((\boldsymbol{q}s)\). By such definition, using the relations that \(\hat{X}(\boldsymbol{q}s)^\dagger = \hat{X}(-\boldsymbol{q}s), \hat{P}(\boldsymbol{q}s)^\dagger = \hat{P}(-\boldsymbol{q}s)\), and that \(\omega(\boldsymbol{q}s)=\omega(-\boldsymbol{q}s)\), the inverse transformation is:

\[ \begin{aligned} \hat{X}(\boldsymbol{q}s) &= -i\sqrt{\frac{\hbar}{2\omega(\boldsymbol{q}s)}} (\hat{a}^\dagger_{\boldsymbol{q}s} - \hat{a}_{-\boldsymbol{q}s})\\ \hat{P}(\boldsymbol{q}s) &= \sqrt{\frac{\hbar\omega(\boldsymbol{q}s)}{2}} (\hat{a}_{\boldsymbol{q}s} + \hat{a}_{-\boldsymbol{q}s}^\dagger) \end{aligned} \]

In addition, the annihilation operator and the creation operator satisfy the following commutation relation:

\[ \begin{aligned} \left[\hat{a}_{\boldsymbol{q}s},\hat{a}_{\boldsymbol{q}'s'}^\dagger\right] &= \frac{1}{2\hbar\sqrt{\omega(\boldsymbol{q}s)\omega(\boldsymbol{q}'s')}}\left[\hat{P}(\boldsymbol{q}s),\hat{P}^\dagger(\boldsymbol{q}'s')\right] + \frac{\sqrt{\omega(\boldsymbol{q}s)\omega(\boldsymbol{q}'s')}}{2\hbar}\left[\hat{X}(\boldsymbol{q}s),\hat{X}^\dagger(\boldsymbol{q}'s')\right] \\ &\qquad - \frac{i}{2\hbar}\sqrt{\frac{\omega (\boldsymbol{q}s)}{\omega (\boldsymbol{q}'s')}} \left[\hat{X}^\dagger(\boldsymbol{q}s), \hat{P}^\dagger(\boldsymbol{q}'s')\right] + \frac{i}{2\hbar}\sqrt{\frac{\omega (\boldsymbol{q}'s')}{\omega (\boldsymbol{q}s)}} \left[\hat{P}(\boldsymbol{q}s),\hat{X}(\boldsymbol{q}'s') \right]\\ &= \delta_{\boldsymbol{q}\boldsymbol{q}'}\delta_{ss'} \end{aligned} \]

Looking back at the \((\boldsymbol{q},\boldsymbol{b})\) representations, we have:

\[ \begin{aligned} \hat{\boldsymbol{X}}(\boldsymbol{q};\boldsymbol{b}) &= \frac{1}{\sqrt{m_{\boldsymbol{b}}}}\sum_s \boldsymbol{e}(\boldsymbol{b}|\boldsymbol{q}s) \hat{X}(\boldsymbol{q}s) = -i \sum_s \sqrt{\frac{\hbar}{2 m_{\boldsymbol{b}}\omega (\boldsymbol{q}s)}} \boldsymbol{e}(\boldsymbol{b}|\boldsymbol{q}s) (\hat{a}_{\boldsymbol{q}s}^\dagger - \hat{a}_{-\boldsymbol{q}s})\\ \hat{\boldsymbol{P}}(\boldsymbol{q};\boldsymbol{b}) &= \sqrt{m_{\boldsymbol{b}}}\sum_s \boldsymbol{e}^*(\boldsymbol{b}|\boldsymbol{q}s) \hat{P}(\boldsymbol{q}s) = \sum_s \sqrt{\frac{\hbar m_{\boldsymbol{b}}\omega (\boldsymbol{q}s)}{2}} \boldsymbol{e}^*(\boldsymbol{b}|\boldsymbol{q}s) (\hat{a}_{\boldsymbol{q}s} + \hat{a}_{-\boldsymbol{q}s}^\dagger) \end{aligned} \]

The Hamiltonian can thus be rewritten as:

\[ \begin{aligned} \text{1st term} &= \sum_{\boldsymbol{q}\boldsymbol{b}}\frac{1}{2m_{\boldsymbol{b}}} \hat{\boldsymbol{P}} (\boldsymbol{q};\boldsymbol{b}) \hat{\boldsymbol{P}}^\dagger (\boldsymbol{q};\boldsymbol{b}) \\ &= \frac{1}{4}\sum_{\boldsymbol{q}\boldsymbol{b}} \sum_{ss'}\sum_{\alpha\beta} \hbar\omega (\boldsymbol{q}s) e^*_{\alpha}(\boldsymbol{b}|\boldsymbol{q}s) e_{\beta}(\boldsymbol{b}|\boldsymbol{q}s') \left(\hat{a}_{\boldsymbol{q}s} + \hat{a}_{-\boldsymbol{q}s}^\dagger\right) \left(\hat{a}^\dagger_{\boldsymbol{q}s'} + \hat{a}_{-\boldsymbol{q}s'}\right)\\ &= \frac{1}{4} \sum_{\boldsymbol{q}} \sum_{s} h\omega (\boldsymbol{q}s) \left(\hat{a}_{\boldsymbol{q}s} + \hat{a}_{-\boldsymbol{q}s}^\dagger\right) \left(\hat{a}^\dagger_{\boldsymbol{q}s'} + \hat{a}_{-\boldsymbol{q}s'}\right)\\ \text{2nd term} &= \frac{1}{2} \sum_{\alpha\beta} \sum_{\boldsymbol{q}\boldsymbol{b}\boldsymbol{b}'} \Phi_{\alpha\beta}(\boldsymbol{b}\boldsymbol{b}'|\boldsymbol{q}) \hat{X}_{\alpha}(\boldsymbol{q};\boldsymbol{b}) \hat{X}_{\beta}^\dagger (\boldsymbol{q};\boldsymbol{b}') \\ &= \frac{1}{2} \sum_{\alpha\beta}\sum_{\boldsymbol{q}\boldsymbol{b}\boldsymbol{b}'} \sum_{ss'} \Phi_{\alpha\beta} (\boldsymbol{b}\boldsymbol{b}'|\boldsymbol{q}) \frac{\hbar}{2\sqrt{m_{\boldsymbol{b}}m_{\boldsymbol{b}'}\omega (\boldsymbol{q}s)\omega (\boldsymbol{q}s')}} e_{\alpha}(\boldsymbol{b}|\boldsymbol{q}s) e_{\beta}^*(\boldsymbol{b}'|\boldsymbol{q}s') \left(\hat{a}_{\boldsymbol{q}s}^\dagger - \hat{a}_{-\boldsymbol{q}s}\right) \left(\hat{a}_{\boldsymbol{q}s'} - \hat{a}^\dagger_{-\boldsymbol{q}s'}\right)\\ &= \frac{1}{2} \sum_{\alpha\beta}\sum_{\boldsymbol{q}\boldsymbol{b}\boldsymbol{b}'} \sum_{ss'} D_{\alpha\beta}^\dagger (\boldsymbol{b}\boldsymbol{b}'|\boldsymbol{q}) \frac{\hbar}{2\sqrt{\omega (\boldsymbol{q}s)\omega (\boldsymbol{q}s')}} e_{\alpha}(\boldsymbol{b}|\boldsymbol{q}s) e_{\beta}^*(\boldsymbol{b}'|\boldsymbol{q}s') \left(\hat{a}_{\boldsymbol{q}s}^\dagger - \hat{a}_{-\boldsymbol{q}s}\right) \left(\hat{a}_{\boldsymbol{q}s'} - \hat{a}^\dagger_{-\boldsymbol{q}s'}\right)\\ & = \frac{1}{4} \sum_{\boldsymbol{q}s} \hbar\omega (\boldsymbol{q}s) \left(\hat{a}_{\boldsymbol{q}s}^\dagger - \hat{a}_{-\boldsymbol{q}s}\right) \left(\hat{a}_{\boldsymbol{q}s} - \hat{a}^\dagger_{-\boldsymbol{q}s}\right) \end{aligned} \]

The derivation of the 2nd term uses the fact that \(\boldsymbol{e}(\boldsymbol{b}|\boldsymbol{q}s)\) is the eigenvector of \(\boldsymbol{D}(\boldsymbol{q})\) (the dynamical matrix, see the previous post), and the corresponding eigenvalue is \(\omega (\boldsymbol{q}s)^2\).

Above, the summation over \(s'\) in the second term is cancelled due to the orthogonality of the eigenvectors for the summation over \(\alpha\) and \(\beta\). Specifically,

\[ \begin{aligned} \sum_{\boldsymbol{b}\boldsymbol{b}'}\sum_{\alpha,\beta} D_{\alpha\beta}^{\dagger} (\boldsymbol{b}\boldsymbol{b}'|\boldsymbol{q}) e_{\alpha} (\boldsymbol{b}|\boldsymbol{q}s) e_{\beta}^* (\boldsymbol{b}'|\boldsymbol{q}s') &= \omega^2 (\boldsymbol{q}s') \sum_{\boldsymbol{b}\boldsymbol{b}'}\sum_{\alpha,\beta} e_{\alpha} (\boldsymbol{b}|\boldsymbol{q}s) e_{\beta}^* (\boldsymbol{b}'|\boldsymbol{q}s')\\ &= \omega^2 (\boldsymbol{q}s') \delta_{ss'} \delta_{\boldsymbol{b}\boldsymbol{b}'} \end{aligned} \]

As for the 3rd term, we have:

\[ \begin{aligned} \text{3rd term} &= \frac{1}{3!} \frac{1}{\sqrt{N_0}}\sum_{\alpha\beta\gamma}\sum_{\boldsymbol{q}\boldsymbol{q}'\boldsymbol{q}''\boldsymbol{b}\boldsymbol{b}'\boldsymbol{b}''} \Phi_{\alpha\beta\gamma}(\boldsymbol{q}\boldsymbol{b}, \boldsymbol{q}'\boldsymbol{b}', \boldsymbol{q}''\boldsymbol{b}'') \hat{X}_{\alpha}(\boldsymbol{q};\boldsymbol{b}) \hat{X}_{\beta}(\boldsymbol{q}';\boldsymbol{b}') \hat{X}_{\gamma}(\boldsymbol{q}'';\boldsymbol{b}'') \delta_{\boldsymbol{G},\boldsymbol{q}+\boldsymbol{q}'+\boldsymbol{q}''}\\ &= \frac{1}{3!}\frac{1}{\sqrt{N_0}} \sum_{\alpha\beta\gamma}\sum_{\boldsymbol{q}\boldsymbol{q}'\boldsymbol{q}''\boldsymbol{b}\boldsymbol{b}'\boldsymbol{b}''}\sum_{ss's''} i \sqrt{\frac{\hbar^3}{8m_{\boldsymbol{b}}m_{\boldsymbol{b}'}m_{\boldsymbol{b}''}\omega (\boldsymbol{q}s)\omega (\boldsymbol{q}'s')\omega (\boldsymbol{q}''s'')}}\Phi_{\alpha\beta\gamma}(\boldsymbol{q}\boldsymbol{b}, \boldsymbol{q}'\boldsymbol{b}', \boldsymbol{q}''\boldsymbol{b}'') e_{\alpha}(\boldsymbol{b}|\boldsymbol{q}s) e_{\beta}(\boldsymbol{b}'|\boldsymbol{q}'s') e_{\gamma}(\boldsymbol{b}''|\boldsymbol{q}''s'') \delta_{G,\boldsymbol{q}+\boldsymbol{q}'+\boldsymbol{q}''} \left(\hat{a}^\dagger_{\boldsymbol{q}s}-\hat{a}_{-\boldsymbol{q}s}\right)\left(\hat{a}^\dagger_{\boldsymbol{q}'s'}-\hat{a}_{-\boldsymbol{q}'s'}\right)\left(\hat{a}^\dagger_{\boldsymbol{q}''s''}-\hat{a}_{-\boldsymbol{q}''s''}\right)\\ &= \frac{1}{3!}\sum_{\boldsymbol{q}\boldsymbol{q}'\boldsymbol{q}''}\sum_{ss's''} \underbrace{\sum_{\alpha\beta\gamma\boldsymbol{b}\boldsymbol{b}'\boldsymbol{b}''}\frac{i}{\sqrt{N_0}} \sqrt{\frac{\hbar^3}{8m_{\boldsymbol{b}}m_{\boldsymbol{b}'}m_{\boldsymbol{b}''}\omega (\boldsymbol{q}s)\omega (\boldsymbol{q}'s')\omega (\boldsymbol{q}''s'')}} \Phi_{\alpha\beta\gamma}(\boldsymbol{q}\boldsymbol{b}, \boldsymbol{q}'\boldsymbol{b}', \boldsymbol{q}''\boldsymbol{b}'') e_{\alpha}(\boldsymbol{b}|\boldsymbol{q}s) e_{\beta}(\boldsymbol{b}'|\boldsymbol{q}'s') e_{\gamma}(\boldsymbol{b}''|\boldsymbol{q}''s'') }_{\equiv \Phi (\boldsymbol{q}s,\boldsymbol{q}'s',\boldsymbol{q}''s'')} \delta_{G,\boldsymbol{q}+\boldsymbol{q}'+\boldsymbol{q}''} \left(\hat{a}^\dagger_{\boldsymbol{q}s}-\hat{a}_{-\boldsymbol{q}s}\right)\left(\hat{a}^\dagger_{\boldsymbol{q}'s'}-\hat{a}_{-\boldsymbol{q}'s'}\right)\left(\hat{a}^\dagger_{\boldsymbol{q}''s''}-\hat{a}_{-\boldsymbol{q}''s''}\right)\\ &= \frac{1}{3!}\sum_{\boldsymbol{q}\boldsymbol{q}'\boldsymbol{q}''}\sum_{ss's''} \Phi (\boldsymbol{q}s,\boldsymbol{q}'s',\boldsymbol{q}''s'') \delta_{G,\boldsymbol{q}+\boldsymbol{q}'+\boldsymbol{q}''} \left(\hat{a}^\dagger_{\boldsymbol{q}s}-\hat{a}_{-\boldsymbol{q}s}\right)\left(\hat{a}^\dagger_{\boldsymbol{q}'s'}-\hat{a}_{-\boldsymbol{q}'s'}\right)\left(\hat{a}^\dagger_{\boldsymbol{q}''s''}-\hat{a}_{-\boldsymbol{q}''s''}\right) \end{aligned} \]

To summarize the results by far, we see the Hamiltonian is expressed as:

\[ \begin{aligned} \hat{H} &= \frac{1}{4} \sum_{\boldsymbol{q}s} \hbar\omega (\boldsymbol{q}s) \left(\hat{a}_{\boldsymbol{q}s} + \hat{a}_{-\boldsymbol{q}s}^\dagger\right) \left(\hat{a}^\dagger_{\boldsymbol{q}s'} + \hat{a}_{-\boldsymbol{q}s'}\right) + \frac{1}{4} \sum_{\boldsymbol{q}s} \hbar\omega (\boldsymbol{q}s) \left(\hat{a}_{\boldsymbol{q}s}^\dagger - \hat{a}_{-\boldsymbol{q}s}\right) \left(\hat{a}_{\boldsymbol{q}s} - \hat{a}^\dagger_{-\boldsymbol{q}s}\right)\\ &\qquad + \frac{1}{3!}\sum_{\boldsymbol{q}\boldsymbol{q}'\boldsymbol{q}''}\sum_{ss's''} \Phi (\boldsymbol{q}s,\boldsymbol{q}'s',\boldsymbol{q}''s'') \delta_{G,\boldsymbol{q}+\boldsymbol{q}'+\boldsymbol{q}''} \left(\hat{a}^\dagger_{\boldsymbol{q}s}-\hat{a}_{-\boldsymbol{q}s}\right)\left(\hat{a}^\dagger_{\boldsymbol{q}'s'}-\hat{a}_{-\boldsymbol{q}'s'}\right)\left(\hat{a}^\dagger_{\boldsymbol{q}''s''}-\hat{a}_{-\boldsymbol{q}''s''}\right)\\ &= \sum_{\boldsymbol{q}s} \hbar\omega (\boldsymbol{q}s) \left(\hat{a}^\dagger_{\boldsymbol{q}s} a_{\boldsymbol{q}s} + \frac{1}{2}\right) \\ &\qquad + \frac{1}{3!}\sum_{\boldsymbol{q}\boldsymbol{q}'\boldsymbol{q}''}\sum_{ss's''} \Phi (\boldsymbol{q}s,\boldsymbol{q}'s',\boldsymbol{q}''s'') \delta_{G,\boldsymbol{q}+\boldsymbol{q}'+\boldsymbol{q}''} \left(\hat{a}^\dagger_{\boldsymbol{q}s}-\hat{a}_{-\boldsymbol{q}s}\right)\left(\hat{a}^\dagger_{\boldsymbol{q}'s'}-\hat{a}_{-\boldsymbol{q}'s'}\right)\left(\hat{a}^\dagger_{\boldsymbol{q}''s''}-\hat{a}_{-\boldsymbol{q}''s''}\right) \end{aligned} \]

Harmonic oscillator

Let's see what the present results can do. First, we ignore the already-derived third-order term in the Hamiltonian, and focus on the harmonic potential. We have:

\[ \hat{H}_{\text{harm}} = \sum_{\boldsymbol{q}s} \hbar\omega (\boldsymbol{q}s) \left(\hat{a}^\dagger_{\boldsymbol{q}s} \hat{a}_{\boldsymbol{q}s} + \frac{1}{2}\right) \]

Seems to be familiar? Yes, the energy of phonons is just derived from this equation. As well as the quantization of phonon! We recall the many-body Schrödinger equation, say there is a state \(\Psi\), whose energy is \(E\), thus \(\hat{H}_{\text{harm}}\Psi = E\Psi\). After using the annihilator or the creation operator, we have:

\[ \begin{aligned} \hat{H}_{\text{harm}}(\hat{a}_{\boldsymbol{q}s}\Psi) &= \hbar\omega (\boldsymbol{q}s) \left(\hat{a}^\dagger_{\boldsymbol{q}s} \hat{a}_{\boldsymbol{q}s} + \frac{1}{2}\right) \hat{a}_{\boldsymbol{q}s} \Psi = \hbar \omega(\boldsymbol{q}s) \left(\hat{a}_{\boldsymbol{q}s}\hat{a}_{\boldsymbol{q}s}^\dagger \hat{a}_{\boldsymbol{q}s} - \frac{1}{2}\hat{a}_{\boldsymbol{q}s}\right) \Psi\\ &= \hat{a}_{\boldsymbol{q}s}\left(\hat{H}-\hbar\omega (\boldsymbol{q}s)\right) \Psi = (E-\hbar\omega (\boldsymbol{q}s))(\hat{a}_{\boldsymbol{q}s}\Psi)\\ \hat{H}_{\text{harm}}(\hat{a}^\dagger_{\boldsymbol{q}s}\Psi) &= \hbar\omega (\boldsymbol{q}s) \left(\hat{a}^\dagger_{\boldsymbol{q}s} \hat{a}_{\boldsymbol{q}s} + \frac{1}{2}\right) \hat{a}^\dagger_{\boldsymbol{q}s} \Psi = \hbar \omega(\boldsymbol{q}s) \left(\hat{a}^\dagger_{\boldsymbol{q}s}\hat{a}_{\boldsymbol{q}s} \hat{a}^\dagger_{\boldsymbol{q}s} + \frac{1}{2}\hat{a}^\dagger_{\boldsymbol{q}s}\right) \Psi\\ &= \hat{a}^\dagger_{\boldsymbol{q}s}\left( \hat{H}+\hbar\omega (\boldsymbol{q}s)\right) \Psi = (E+\hbar\omega (\boldsymbol{q}s))(\hat{a}^\dagger_{\boldsymbol{q}s}\Psi) \end{aligned} \]

That is to say, through the action of the annihilation or the creation operator, the energy of the state is changed by \(\pm\hbar\omega (\boldsymbol{q}s)\). And the energy of the state is quantized, that is, \(E = \hbar\omega (\boldsymbol{q}s) n_{\boldsymbol{q}s} + \frac{1}{2}\hbar\omega (\boldsymbol{q}s)\), where \(n_{\boldsymbol{q}s}\) is the number of phonons in the state \((\boldsymbol{q}s)\).

If we use \(\Ket{n_{\boldsymbol{q}s}}\) (normalized eigenvector) to denote the state with \(n_{\boldsymbol{q}s}\) phonons in the state \((\boldsymbol{q}s)\), we know that,

\[ \begin{aligned} \hat{a}_{\boldsymbol{q}s}\Ket{n_{\boldsymbol{q}s}} &\propto \Ket{n_{\boldsymbol{q}s}-1}\\ \hat{a}^\dagger_{\boldsymbol{q}s}\Ket{n_{\boldsymbol{q}s}} &\propto \Ket{n_{\boldsymbol{q}s}+1} \end{aligned} \]

But the proportionality constant is not known yet. Due to Schrödinger equation,

\[ \bra{n_{\boldsymbol{q}s}} \hat{a}_{\boldsymbol{q}s}^\dagger \hat{a}_{\boldsymbol{q}s} \ket{n_{\boldsymbol{q}s}} = n_{\boldsymbol{q}s} = n_{\boldsymbol{q}s} \braket{n_{\boldsymbol{q}s}-1|n_{\boldsymbol{q}s}-1} \]

Thus \(\hat{a}_{\boldsymbol{q}s}\ket{n_{\boldsymbol{q}s}} = \sqrt{n_{\boldsymbol{q}s}}\ket{n_{\boldsymbol{q}s}-1}\). Similarly we can obtain that \(\hat{a}^\dagger_{\boldsymbol{q}s}\ket{n_{\boldsymbol{q}s}} = \sqrt{n_{\boldsymbol{q}s}+1}\ket{n_{\boldsymbol{q}s}+1}\).

Three-phonon scattering

Three-phonon scattering can be computed by the third-order term in the Hamiltonian using Fermi's golden rule, which calculates the transition rate between initial state \(i\) and final state \(j\) by:

\[ P_{i\rightarrow j} = \frac{2\pi}{\hbar} \big|\bra{i}\mathcal{V}_3\ket{j}\big|^2 \delta(E_i-E_j) \]

There are two classes of three-phonon scattering:

  • Class 1: \(\boldsymbol{q}+\boldsymbol{q}' \rightarrow \boldsymbol{q}''\)
  • Class 2: \(\boldsymbol{q}\rightarrow\boldsymbol{q}' + \boldsymbol{q}''\)

For class 1, the scattering is given by:

\[ \begin{aligned} P_{\boldsymbol{q}s,\boldsymbol{q}'s'\rightarrow \boldsymbol{q}''s''} &= \frac{2\pi}{\hbar} \big|\bra{n_{\boldsymbol{q}s}-1,n_{\boldsymbol{q}'s'}-1,n_{\boldsymbol{q}''s''}+1}\mathcal{V}_3\ket{n_{\boldsymbol{q}s},n_{\boldsymbol{q}'s'},n_{\boldsymbol{q}''s''}}\big|^2\delta (\omega(\boldsymbol{q}s)+\omega(\boldsymbol{q}'s')-\omega(\boldsymbol{q}''s''))\\ &= \frac{2\pi}{\hbar} \big|\Phi(-\boldsymbol{q}s,-\boldsymbol{q}'s',\boldsymbol{q}''s'')\big|^2 n_{\boldsymbol{q}s}n_{\boldsymbol{q}'s'} (n_{\boldsymbol{q}''s''}+1)\delta (\omega(\boldsymbol{q}s)+\omega(\boldsymbol{q}'s')-\omega(\boldsymbol{q}''s'')) \end{aligned} \]

The net scattering rate is:

\[ \begin{aligned} P_{\boldsymbol{q}s,\boldsymbol{q}'s'\rightarrow \boldsymbol{q}''s''} - P_{\boldsymbol{q}''s'',\rightarrow \boldsymbol{q}s,\boldsymbol{q}'s'} &= \frac{2\pi}{\hbar} \big|\Phi(-\boldsymbol{q}s,-\boldsymbol{q}'s',\boldsymbol{q}''s'')\big|^2 \delta (\omega(\boldsymbol{q}s)+\omega(\boldsymbol{q}'s')-\omega(\boldsymbol{q}''s''))\\ &\quad \times \Big[n_{\boldsymbol{q}s}n_{\boldsymbol{q}'s'} (n_{\boldsymbol{q}''s''}+1) - (n_{\boldsymbol{q}s}+1)(n_{\boldsymbol{q}'s'}+1) n_{\boldsymbol{q}''s''} \Big] \end{aligned} \]

Similarly, for class 2, the scattering is given by:

\[ \begin{aligned} P_{\boldsymbol{q}s\rightarrow \boldsymbol{q}'s',\boldsymbol{q}''s''} - P_{\boldsymbol{q}'s',\boldsymbol{q}''s'' \rightarrow \boldsymbol{q}s} &= \frac{2\pi}{\hbar} \big|\Phi(\boldsymbol{q}s,-\boldsymbol{q}'s',-\boldsymbol{q}''s'')\big|^2 \delta (\omega(\boldsymbol{q}s)-\omega(\boldsymbol{q}'s')-\omega(\boldsymbol{q}''s''))\\ &\quad \times \Big[n_{\boldsymbol{q}s} (n_{\boldsymbol{q}'s'}+1) (n_{\boldsymbol{q}''s''}+1) - (n_{\boldsymbol{q}s}+1)n_{\boldsymbol{q}'s'}n_{\boldsymbol{q}''s''}\Big] \end{aligned} \]

Yielding the net three-phonon scattering rate,

\[ \begin{aligned} -\left(\frac{\partial n_{\boldsymbol{q}s}}{\partial t}\right)_{\text{scatt}} &= \sum_{\boldsymbol{q}'s'\boldsymbol{q}''s''} \left[(P_{\boldsymbol{q}s,\boldsymbol{q}'s'\rightarrow \boldsymbol{q}''s''} - P_{\boldsymbol{q}''s'',\rightarrow \boldsymbol{q}s,\boldsymbol{q}'s'}) + \frac{1}{2} (P_{\boldsymbol{q}s\rightarrow \boldsymbol{q}'s',\boldsymbol{q}''s''} - P_{\boldsymbol{q}'s',\boldsymbol{q}''s'' \rightarrow \boldsymbol{q}s})\right] \end{aligned} \]

where the factor \(\frac{1}{2}\) is to avoid double counting.

References

  1. Griffiths's Introduction to quantum mechanics.
  2. Srivastava The Physics of Phonons.
  3. Garg's PhD thesis Thermal Conductivity from First-Principles in Bulk, Disordered, and Nanostructured Materials (MIT 2011).