Including Phonons and Crystal-Field Phonon interactions

In this section we discuss, how lattice dynamics may be considered in the framework of the Hamiltonian (60). We will see that this corresponds to a system of coupled Einstein oscillators. One such oscillator can be modelled by setting up a sipf file with the module phonon. Coupling has to be done in mcdisp.j. Rephrasing lattice dynamics in this way allows to couple phonons to the crystal field.

A three dimensional Einstein oscillator (for atom $n$) in a solid can be described by the following Hamiltonian

$$H_E(n)=\frac{a_0^2{\mathbf p_n}^2}{2m_n} - \frac{1}{2} {\mathbf u}^T_n \overline{K}(nn) {\mathbf u}_n$$ (118)

Here $\mathbf u$ is the dimensionless displacement vector ( $\mathbf u_n={\mathbf P}_n/a_0=\Delta {\mathbf r}_n/a_0$, with the Bohr radius $a_0=0.5219$ Å), $m_n$ the mass of the atom $n$, $\mathbf p_n=d\mathbf u_n/dt$ the conjugate momentum to $\mathbf u_n$ and $\overline{K}(nn)$ the Matrix describing the restoring force.

Coupling such oscillators leads to the Hamiltonian

$$H_{phon}=\sum_n H_E(n) -\frac{1}{2} \sum_{n\neq n'} {\mathbf u}_n^T \overline{K}(nn') {\mathbf u}_{n'}$$ (119)

Note that our coupling constants $K_{\alpha\beta}(nn')=-A_{\alpha\beta}(nn')$, where $A_{\alpha\beta}$ are the second-order derivatives of the potenatial energy as defined e.g. in [32, page 99].

In a mean field type of theory the phonon single ion module has thus to solve the Hamiltonian

$$H_E=\frac{a_0^2{\mathbf p}^2}{2m} - \frac{1}{2} {\mathbf u}^T \overline{K} {\mathbf u} - {\mathbf F}^T {\mathbf u}$$ (120)

Here the force $\mathbf F$ corresponds to the exchange field $\mathbf H_{xc}$ and $\mathbf u$ to the general operator $\mathbf I$ and $\overline{K}(nn')$ to $\mathcal J(nn')$ of equation (60), respectively. The single ion Hamiltonian (123) can be solved by transforming it to normal coordinates (main axis of the Einstein oscillator) using the transformation matrix $\overline{S}$, which diagonalises $\overline{K}=\overline{S}^T\overline{\Omega}\overline{S}$:

$$\mathbf u'=\overline{S} \mathbf u -\overline{\Omega}^{-1} \overline{S} \mathbf F$$ (121)

$$H_E=\frac{a_0^2{\mathbf p'}^2}{2m}+\frac{1}{2} {\mathbf u'}^... ...T \overline{S}^T \overline{\Omega}^{-1} \overline{S} \mathbf F$$ (122)

Due to the action of the force $\mathbf F$ the equilibrium position of the oscillator is $\mathbf u_0=\overline{S}^T\overline{\Omega}^{-1}\overline{S}\mathbf F$ (it is the task of the function Icalc to return this equilibrium position), the energies correspond to the three elements of the diagonal matrix $\overline{\Omega}$, i.e. $\Omega_{11}=m a_0^2 (\Delta_1 /\hbar)^2$, $\Omega_{22}=m a_0^2 (\Delta_2 /\hbar)^2$, $\Omega_{33}=m a_0^2 (\Delta_3 /\hbar)^2$. In order to run mcdisp we have to calculate the transition matrix elements:

The single ion susceptibility for such a transition, e.g. $\Delta_1$ - corresponds to

$$\overline{\chi}^0 = \sum_{\nu\mu}\frac{\langle \nu\vert\mathbf u\vert\mu\rangle\langle \mu \vert\mathbf u^T\vert\nu \rangle}{\Delta_1 -\hbar \omega} (p_{\nu} -p_{\mu})$$ (123)

$$= \sum_{\nu\mu}\frac{\langle \nu\vert\overline{S}^T \mathbf u'\vert... ...f u'^T \overline{S}\vert\nu \rangle}{\Delta_1 -\hbar \omega} (p_{\nu} -p_{\mu})$$ (124)

Because the different components of $\mathbf u'$ commute and the Hamiltonian (125) is separable, for the transition $\Delta_1$ only the terms with $u_1'$ in the nominator contribute:

$$\chi^0_{\alpha\beta} = S^T_{\alpha1}\sum_{\nu\mu}\frac{\langle \nu\vert u_1'\vert\mu\ran... ...ert u_1'\vert\nu \rangle}{\Delta_1 -\hbar \omega} (p_{\nu} -p_{\mu}) S_{1\beta}$$ (125)

$$= S^T_{\alpha1}S_{1\beta}\frac{\hbar^2}{2ma_0^2\Delta_1}\left(\frac{1}{\Delta_1-\hbar\omega}+\frac{1}{\Delta_1+\hbar\omega}\right )$$ (126)

In order to derive the last result we had to express $u_1'$ in terms of ladder operators $u_1'=a_0^{-1} \hbar/\sqrt{2m\Delta_1}(a+a^{\dagger})$ and apply $a^{\dagger}\vert\nu\rangle=\sqrt{\nu+1}\vert\nu+1\rangle$, $a\vert\nu\rangle=\sqrt{\nu}\vert\nu-1\rangle$ and $\sum_{\nu=0}^{\infty}(p_{\nu}-p_{\nu+1})(\nu+1)=1$, $p_{\nu}=exp(-\nu\Delta_1/kT)(1-exp(-\Delta_1/kT))$. This shows that the single ion susceptibility of our atom can be written as a sum of three effective transitions (with temperature independent susceptibility)

$$\chi^0_{\alpha\beta} = \sum_{i=1,2,3} S^T_{\alpha i}S_{i\beta}\frac{\hbar^2}{2ma_0^2\Delta_i} \left(\frac{1}{\Delta_i-\hbar\omega}+\frac{1}{\Delta_i+\hbar\omega}\right )$$ (127)

Thus the module phonon has to provide in it's function du1calc these three transitions (=number of transitions).