Solve phonon dispersion relation by the classical approach
Consider a crystal in which unit cells are labeled by \(\boldsymbol{l}\), and atoms in each cell are labeled by \(\boldsymbol{b}\). The total energy of the system, noted by \(\mathcal{V}\), can be expanded using Taylor's expansion (greek letters mean direction): \[ \begin{aligned} \mathcal{V} = & \mathcal{V}_0 + \underbrace{\sum_{\boldsymbol{l},\boldsymbol{b}}\sum_{\alpha} \Phi_{\alpha}(\boldsymbol{l}\boldsymbol{b}) u_{\alpha}(\boldsymbol{l}\boldsymbol{b})}_{\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} + \cdots\\ & \quad + \frac{1}{n!}\sum_{\boldsymbol{l}_1,\boldsymbol{b}_1,\boldsymbol{l}_2,\boldsymbol{b}_2,\cdots,\boldsymbol{l}_n,\boldsymbol{b}_n}\sum_{\alpha_1,\alpha_2,\cdots,\alpha_n}\Phi_{\alpha_1,\alpha_2,\cdots,\alpha_n}(\boldsymbol{l}_1,\boldsymbol{b}_1;\boldsymbol{l}_2,\boldsymbol{b}_2;\cdots;\boldsymbol{l}_n,\boldsymbol{b}_n)u_{\alpha_1}(\boldsymbol{l}_1,\boldsymbol{b}_1)\cdots u_{\alpha_n}(\boldsymbol{l}_n \boldsymbol{b}_n) \end{aligned} \] Where we note the nth order interatomic force constants (IFCs) using \(\Phi\), and we note \(u_{\alpha}(\boldsymbol{l}\boldsymbol{b})\) represent the spatial deviation of atom \(\boldsymbol{b}\) in cell \(\boldsymbol{l}\) along \(\alpha\) direction. It is clear to see that:
- The first order IFCs should be cancelled, that \(\mathcal{V}_1=0\), due to that the system is stable.
- The second order term \(\mathcal{V}_2\) represent the combination of many interatomic harmonic oscillators.
Let's ignore the higher order term, keeping only \(\mathcal{V}_2\). Using Newton's second law, \[ F_{\boldsymbol{l}\boldsymbol{b}\alpha} = m_{\boldsymbol{b}} \ddot{u}_{\alpha}(\boldsymbol{l}\boldsymbol{b}) =-\frac{\partial \mathcal{V}}{\partial \mathcal{u_{\alpha}}(\boldsymbol{l}\boldsymbol{b})}= - \sum_{\boldsymbol{l}'\boldsymbol{b}'}\sum_{\beta} \Phi_{\alpha,\beta}(\boldsymbol{l}\boldsymbol{b};\boldsymbol{l}'\boldsymbol{b}') u_{\beta}(\boldsymbol{l}'\boldsymbol{b}'),\quad \forall \boldsymbol{l},\boldsymbol{b},\alpha \] Note that the above equation exploits the symmetry of second order IFCs, that \(\Phi_{\alpha,\beta}(\boldsymbol{l}\boldsymbol{b};\boldsymbol{l}'\boldsymbol{b}') = \Phi_{\beta,\alpha}(\boldsymbol{l}'\boldsymbol{b}',\boldsymbol{l}\boldsymbol{b})\). We now seek a solution in the form of plane waves: \[ u_{\alpha} (\boldsymbol{l}\boldsymbol{b}) = \frac{1}{\sqrt{m_{\boldsymbol{b}}}} \sum_{\boldsymbol{q}} U_{\alpha}(\boldsymbol{q};\boldsymbol{b})\exp \big[ i \boldsymbol{q}\cdot \boldsymbol{l}-i\omega t\big] \]
- The purpose of seaking a solution that is multiplied by \(1/\sqrt{m_{\boldsymbol{b}}}\) is just to make the form of dynamical matrix simple and beautiful, so don't be confused.
- The value \(U_{\alpha}(\boldsymbol{q};\boldsymbol{b})\) is complex. Specifically, it contains the phase information of atom \(\boldsymbol{b}\).
A specific wavevector \(\boldsymbol{q}\) yields, \[ \omega^2 \sqrt{m_{\boldsymbol{b}}} U_{\alpha}(\boldsymbol{q};\boldsymbol{b}) \exp \big[i\boldsymbol{q}\cdot\boldsymbol{l}-i\omega t\big] = \sum_{\boldsymbol{l}'\boldsymbol{b}'}\sum_{\beta}\Phi_{\alpha,\beta}(\boldsymbol{l}\boldsymbol{b};\boldsymbol{l}'\boldsymbol{b}') \frac{1}{\sqrt{m_{\boldsymbol{b}'}}}U_{\beta}(\boldsymbol{q};\boldsymbol{b}')\exp\big[i\boldsymbol{q}\cdot\boldsymbol{l}'-i\omega t\big] \] Cancel the exponential term on both sides. Notice that due to translational invariance, \(\Phi_{\alpha,\beta}(\boldsymbol{l}\boldsymbol{b};\boldsymbol{l}'\boldsymbol{b}')=\Phi_{\alpha,\beta}(\boldsymbol{0}\boldsymbol{b};(\boldsymbol{l}-\boldsymbol{l}')\boldsymbol{b}')\), thus, \[ \begin{aligned} \omega^2 \sqrt{m_{\boldsymbol{b}}} U_{\alpha}(\boldsymbol{q};\boldsymbol{b}) &= \sum_{\boldsymbol{l}'\boldsymbol{b}'}\sum_{\beta}\Phi_{\alpha,\beta}(\boldsymbol{l}\boldsymbol{b};\boldsymbol{l}'\boldsymbol{b}') \frac{1}{\sqrt{m_{\boldsymbol{b}'}}}U_{\beta}(\boldsymbol{q};\boldsymbol{b}')\exp\big[i\boldsymbol{q}\cdot(\boldsymbol{l}'-\boldsymbol{l})\big]\\ &= \sum_{\boldsymbol{l}'\boldsymbol{b}'}\sum_{\beta}\Phi_{\alpha,\beta}(\boldsymbol{0}\boldsymbol{b};(\boldsymbol{l}'-\boldsymbol{l})\boldsymbol{b}') \frac{1}{\sqrt{m_{\boldsymbol{b}'}}}U_{\beta}(\boldsymbol{q};\boldsymbol{b}')\exp\big[i\boldsymbol{q}\cdot(\boldsymbol{l}'-\boldsymbol{l})\big]\\ &= \sum_{\boldsymbol{l}'\boldsymbol{b}'}\sum_{\beta}\Phi_{\alpha,\beta}(\boldsymbol{0}\boldsymbol{b};\boldsymbol{l}'\boldsymbol{b}') \frac{1}{\sqrt{m_{\boldsymbol{b}'}}}U_{\beta}(\boldsymbol{q};\boldsymbol{b}')\exp\big[i\boldsymbol{q}\cdot\boldsymbol{l}'\big] \end{aligned} \] Therefore, if we note a matrix called "dynamical matrix", as: \[ D_{\alpha\beta} (\boldsymbol{b}\boldsymbol{b}' | \boldsymbol{q}) = \frac{1}{\sqrt{m_{\boldsymbol{b}}m_{\boldsymbol{b}'}}} \sum_{\boldsymbol{l}'} \Phi_{\alpha,\beta}(\boldsymbol{0}\boldsymbol{b};\boldsymbol{l}'\boldsymbol{b}')\exp[i\boldsymbol{q}\cdot\boldsymbol{l}'] \] The dynamical equations are thus in the form of matrix multiplication, \[ \omega^2 U_{\alpha} (\boldsymbol{q};\boldsymbol{b}) = \sum_{\boldsymbol{b}',\beta}D_{\alpha\beta}(\boldsymbol{b}\boldsymbol{b}'|\boldsymbol{q}) U_{\beta} (\boldsymbol{q};\boldsymbol{b}') \] Solving the above equation, we find that \(\omega\) is in fact the eigenvalue of the dynamical matrix. Due to that \(D_{\cdot,\cdot}(\cdot,\cdot|\boldsymbol{q})\) is a matrix of \(3p\times 3p\), where \(p\) is the number of atoms per unit cell, and \(3\) is the directions that \(\alpha,\beta\) can take, for any given \(\boldsymbol{q}\), there are \(3p\) eigenvalues and \(3p\) eigenvectors, that correspond to different phonon branches.
Example: 1D diatomic chain
Consider the 1D diatomic chain as following:
1 | g g g g g |
The total energy is: \[ U= \frac{1}{2}g\sum_n \left[(u_{n,1}-u_{n,2})^2 + (u_{n,2}-u_{n+1,1})^2\right] \] where \(n\) labels the unit cell, \(1\) and \(2\) represent \(M\) or \(m\). The dynamical matrix is \(2\times 2\), due to that there are two atoms in each unit cell, and there is only one direction. According to the definition, the dynamical matrix is written as, \[ D_q = \begin{bmatrix} \frac{2g}{M} & -\frac{g}{\sqrt{Mm}}(1+e^{-iqa})\\ -\frac{g}{\sqrt{Mm}}(1+e^{iqa}) & \frac{2g}{m} \end{bmatrix} \] The two eigenvalues are: \[ \Rightarrow \omega^2 = g \left(\frac{1}{M} + \frac{1}{m}\right) \pm g\sqrt{\left(\frac{1}{M} + \frac{1}{m}\right)^2-\frac{4}{Mm} \sin^2\left(\frac{qa}{2}\right)} \] If we plot the two branches from \(q=-\pi/a\) to \(q=\pi/a\), the dispersion relation will be shown as:
