Solve phonon dispersion relation by the classical approach

Consider a crystal in which unit cells are labeled by $\mathit{l}$, and atoms in each cell are labeled by $\mathit{b}$. The total energy of the system, noted by $\mathcal{V}$, can be expanded using Taylor's expansion (greek letters mean direction): $$\begin{array}{rl}\mathcal{V}=& {\mathcal{V}}_0+\underset{{\mathcal{V}}_1}{\underbrace{\sum _{\mathit{l},\mathit{b}}\sum _{\alpha }{\mathrm{\Phi }}_{\alpha }(\mathit{l}\mathit{b})u_{\alpha }(\mathit{l}\mathit{b})}}+\underset{{\mathcal{V}}_2}{\underbrace{\frac{1}{2}\sum _{\mathit{l},\mathit{b},{\mathit{l}}^{\prime },{\mathit{b}}^{\prime }}\sum _{\alpha ,\beta }{\mathrm{\Phi }}_{\alpha ,\beta }(\mathit{l}\mathit{b};{\mathit{l}}^{\prime }{\mathit{b}}^{\prime })u_{\alpha }(\mathit{l}\mathit{b})u_{\beta }({\mathit{l}}^{\prime }{\mathit{b}}^{\prime })}}+\cdots \\ & +\frac{1}{n!}\sum _{{\mathit{l}}_1,{\mathit{b}}_1,{\mathit{l}}_2,{\mathit{b}}_2,\cdots ,{\mathit{l}}_n,{\mathit{b}}_n}\sum _{{\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_n}{\mathrm{\Phi }}_{{\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_n}({\mathit{l}}_1,{\mathit{b}}_1;{\mathit{l}}_2,{\mathit{b}}_2;\cdots ;{\mathit{l}}_n,{\mathit{b}}_n)u_{{\alpha }_1}({\mathit{l}}_1,{\mathit{b}}_1)\cdots u_{{\alpha }_n}({\mathit{l}}_n{\mathit{b}}_n)\end{array}$$ Where we note the nth order interatomic force constants (IFCs) using $\mathrm{\Phi }$, and we note $u_{\alpha }(\mathit{l}\mathit{b})$ represent the spatial deviation of atom $\mathit{b}$ in cell $\mathit{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_{\mathit{l}\mathit{b}\alpha }=m_{\mathit{b}}{\ddot{u}}_{\alpha }(\mathit{l}\mathit{b})=-\frac{\partial \mathcal{V}}{\partial {\mathcal{u}}_{\alpha }(\mathit{l}\mathit{b})}=-\sum _{{\mathit{l}}^{\prime }{\mathit{b}}^{\prime }}\sum _{\beta }{\mathrm{\Phi }}_{\alpha ,\beta }(\mathit{l}\mathit{b};{\mathit{l}}^{\prime }{\mathit{b}}^{\prime })u_{\beta }({\mathit{l}}^{\prime }{\mathit{b}}^{\prime }),\mathrm{\forall }\mathit{l},\mathit{b},\alpha$$ Note that the above equation exploits the symmetry of second order IFCs, that ${\mathrm{\Phi }}_{\alpha ,\beta }(\mathit{l}\mathit{b};{\mathit{l}}^{\prime }{\mathit{b}}^{\prime })={\mathrm{\Phi }}_{\beta ,\alpha }({\mathit{l}}^{\prime }{\mathit{b}}^{\prime },\mathit{l}\mathit{b})$. We now seek a solution in the form of plane waves: $$u_{\alpha }(\mathit{l}\mathit{b})=\frac{1}{\sqrt{m_{\mathit{b}}}}\sum _{\mathit{q}}{U}_{\alpha }(\mathit{q};\mathit{b})\exp [i\mathit{q}\cdot \mathit{l}-i\omega t]$$

• The purpose of seaking a solution that is multiplied by $1/\sqrt{m_{\mathit{b}}}$ is just to make the form of dynamical matrix simple and beautiful, so don't be confused.
• The value $U_{\alpha }(\mathit{q};\mathit{b})$ is complex. Specifically, it contains the phase information of atom $\mathit{b}$.

A specific wavevector $\mathit{q}$ yields, $${\omega }^2\sqrt{m_{\mathit{b}}}{U}_{\alpha }(\mathit{q};\mathit{b})\exp [i\mathit{q}\cdot \mathit{l}-i\omega t]=\sum _{{\mathit{l}}^{\prime }{\mathit{b}}^{\prime }}\sum _{\beta }{\mathrm{\Phi }}_{\alpha ,\beta }(\mathit{l}\mathit{b};{\mathit{l}}^{\prime }{\mathit{b}}^{\prime })\frac{1}{\sqrt{m_{{\mathit{b}}^{\prime }}}}{U}_{\beta }(\mathit{q};{\mathit{b}}^{\prime })\exp [i\mathit{q}\cdot {\mathit{l}}^{\prime }-i\omega t]$$ Cancel the exponential term on both sides. Notice that due to translational invariance, ${\mathrm{\Phi }}_{\alpha ,\beta }(\mathit{l}\mathit{b};{\mathit{l}}^{\prime }{\mathit{b}}^{\prime })={\mathrm{\Phi }}_{\alpha ,\beta }(\mathbf{0}\mathit{b};(\mathit{l}-{\mathit{l}}^{\prime }){\mathit{b}}^{\prime })$, thus, $$\begin{array}{rl}{\omega }^2\sqrt{m_{\mathit{b}}}{U}_{\alpha }(\mathit{q};\mathit{b})& =\sum _{{\mathit{l}}^{\prime }{\mathit{b}}^{\prime }}\sum _{\beta }{\mathrm{\Phi }}_{\alpha ,\beta }(\mathit{l}\mathit{b};{\mathit{l}}^{\prime }{\mathit{b}}^{\prime })\frac{1}{\sqrt{m_{{\mathit{b}}^{\prime }}}}{U}_{\beta }(\mathit{q};{\mathit{b}}^{\prime })\exp [i\mathit{q}\cdot ({\mathit{l}}^{\prime }-\mathit{l})]\\ & =\sum _{{\mathit{l}}^{\prime }{\mathit{b}}^{\prime }}\sum _{\beta }{\mathrm{\Phi }}_{\alpha ,\beta }(\mathbf{0}\mathit{b};({\mathit{l}}^{\prime }-\mathit{l}){\mathit{b}}^{\prime })\frac{1}{\sqrt{m_{{\mathit{b}}^{\prime }}}}{U}_{\beta }(\mathit{q};{\mathit{b}}^{\prime })\exp [i\mathit{q}\cdot ({\mathit{l}}^{\prime }-\mathit{l})]\\ & =\sum _{{\mathit{l}}^{\prime }{\mathit{b}}^{\prime }}\sum _{\beta }{\mathrm{\Phi }}_{\alpha ,\beta }(\mathbf{0}\mathit{b};{\mathit{l}}^{\prime }{\mathit{b}}^{\prime })\frac{1}{\sqrt{m_{{\mathit{b}}^{\prime }}}}{U}_{\beta }(\mathit{q};{\mathit{b}}^{\prime })\exp [i\mathit{q}\cdot {\mathit{l}}^{\prime }]\end{array}$$ Therefore, if we note a matrix called "dynamical matrix", as: $$D_{\alpha \beta }(\mathit{b}{\mathit{b}}^{\prime }|\mathit{q})=\frac{1}{\sqrt{m_{\mathit{b}}{m}_{{\mathit{b}}^{\prime }}}}\sum _{{\mathit{l}}^{\prime }}{\mathrm{\Phi }}_{\alpha ,\beta }(\mathbf{0}\mathit{b};{\mathit{l}}^{\prime }{\mathit{b}}^{\prime })\exp [i\mathit{q}\cdot {\mathit{l}}^{\prime }]$$ The dynamical equations are thus in the form of matrix multiplication, $${\omega }^2{U}_{\alpha }(\mathit{q};\mathit{b})=\sum _{{\mathit{b}}^{\prime },\beta }{D}_{\alpha \beta }(\mathit{b}{\mathit{b}}^{\prime }|\mathit{q})U_{\beta }(\mathit{q};{\mathit{b}}^{\prime })$$ 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 |\mathit{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 $\mathit{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:

 g     g     g     g     g
---[m]---[M]---[m]---[M]---
     |<--- a --->|

The total energy is: $$U=\frac{1}{2}g\sum _n[(u_{n,1}-u_{n,2}{)}^2+(u_{n,2}-u_{n+1,1}{)}^2]$$ 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=\left[\begin{array}{cc}\frac{2g}{M}& -\frac{g}{\sqrt{Mm}}(1+e^{-iqa})\\ -\frac{g}{\sqrt{Mm}}(1+e^{iqa})& \frac{2g}{m}\end{array}\right]$$ The two eigenvalues are: $$\Rightarrow {\omega }^2=g(\frac{1}{M}+\frac{1}{m})\pm g\sqrt{{(\frac{1}{M}+\frac{1}{m})}^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: