The Dielectric Continuum model , describes optical phonons in bulk materials or heterostructures. As presented
below, by solving for the associated phonon potential, it will find only modes which produce a potential, ie modes which
are purely transverse will not be identified. In wurtzite materials, “ordinary” modes are purely transverse, but the
“extraordinary” modes, even those known as transverse optical (TO), do have some non-vanishing longitudinal component and
will appear from the following math.
Equations of Motion
If the system is described by a driven oscillator equation
\[\newcommand{\abs}[1]{\left|#1\right|}
\begin{equation}
\mu\partial_t^2u_i =-\mu\omega_{0i}^2u_i+e^* E_i
\end{equation}\]
where \(\mu\) and \(i\) runs over \(\perp, \parallel\). Then \(P_i=\frac{e^*}{V_u}u_i\), giving
\[\begin{align}
P_i =\varepsilon_0\chi_iE_i, \quad \chi_i=\frac{e^{*2}}{\varepsilon_0\mu V_u}\frac{1}{\omega_{0i}^2-\omega^2}
\end{align}\]
Defining an \(\varepsilon_i\) so that \(D_i=\varepsilon_i E_i=\varepsilon^\infty E_i +P=(\varepsilon^\infty +\varepsilon_0\chi_i) E_i\) leads to
\[\begin{equation}
\varepsilon_\perp=\varepsilon^\infty\frac{\omega_{LO\perp}^2-\omega^2}{\omega_{TO\perp}^2-\omega^2},
\quad\quad
\varepsilon_\parallel=\varepsilon^\infty\frac{\omega_{LO\parallel}^2-\omega^2}{\omega_{TO\parallel}^2-\omega^2}
\end{equation}\]
where \(\omega_{TO,i}^2=\omega_{0i}^2\) and \(\omega_{LO, i}^2=\omega_{0i}^2+\frac{e^{*2}}{\mu V_u\varepsilon^\infty}\)
The phonons are found by solving
\[\begin{equation}
-\nabla\cdot\overset\leftrightarrow{\epsilon}\nabla\phi=0
\end{equation}\]
which can be framed as an eigenvalue problem in \(q\) at fixed \(\omega\)
\[\begin{equation}
\partial_z\epsilon_\parallel\partial_z\phi=q^2\epsilon_\perp\phi
\end{equation}\]
Normalization
The \(\phi\) must then be normalized by the prescription \(\int dz\frac{\mu}{V_u}\abs{\vec u}^2=\frac{\hbar}{2\omega}\). To do this, we’ll need to extract \(\vec u\) from the \(\phi\):
\[\begin{align}
\vec u = \frac{V_u}{e^*}\vec P=\frac{V_u}{e^*}\varepsilon_0\chi\vec E=-\frac{V_u}{e^*}\chi\nabla\phi
\end{align}\]
So
\[\begin{align}
\vec u = -\frac{\varepsilon_0V_u}{e^*}\left( \chi_\parallel \partial_z\phi \hat{z} +iq\chi_\perp\phi\hat{q} \right)
\end{align}\]
\[\begin{align}
\frac{\mu}{V_u}\abs{\vec u}^2 = \frac{\varepsilon_0\mu V_u}{e^{*2}}\left( (\chi_\parallel \partial_z\phi)^2 +q^2(\chi_\perp\phi)^2 \right)
\end{align}\]
\[\begin{align}
\frac{\mu}{V_u}\abs{\vec u}^2 = \frac{e^{*2}}{\mu V_u}\left( \left(\frac{\partial_z\phi}{\omega_{TO\parallel}^2-\omega^2}\right)^2 +\left(\frac{q\phi}{\omega^2_{TO\perp}-\omega^2}\right)^2 \right)
\end{align}\]
\[\begin{align}
\frac{\mu}{V_u}\abs{\vec u}^2 = \varepsilon_\infty(\omega^2_{LO}-\omega^2_{TO})\left( \left(\frac{\partial_z\phi}{\omega_{TO\parallel}^2-\omega^2}\right)^2 +\left(\frac{q\phi}{\omega^2_{TO\perp}-\omega^2}\right)^2 \right)
\end{align}\]
where, since \(\omega_{LO\perp}\) et al are taken from experiment, we’ll average the \((\omega_{LO}^2-\omega_{TO}^2)\) factor over \(\perp, \parallel\).
Thus the normalization condition is written in terms of the potential as
\[\begin{align}
\frac{\hbar}{2\omega} = \int dz \varepsilon_\infty(\omega^2_{LO}-\omega^2_{TO})\left( \left(\frac{\partial_z\phi}{\omega_{TO\parallel}^2-\omega^2}\right)^2 +\left(\frac{q\phi}{\omega^2_{TO\perp}-\omega^2}\right)^2 \right)
\end{align}\]
Binary Wurtzite Heterojunction
Following the quantum well example of Komirenko1999, DielectricContinuum_SWH solves a
single heterojunction structure with uniform Wurtzite regions 1/2, assuming the bottom region is semi-infinite.
The top surface at \(z=0\) is assumed Dirichelet.
The thickness of the top layer is \(t_1\) and the normalization thickness for the bottom layer which will be used to
break the continuum of states by an artificial Dirichelet condition is \(t_2\).
It is required that the frequencies be in this order:
\[\omega_{TO\perp, s}<
\omega_{TO\parallel, s}<
\omega_{TO\perp, f}<
\omega_{TO\parallel, f}<
\omega_{LO\perp, s}<
\omega_{LO\parallel, s}<
\omega_{LO\perp, f}<
\omega_{LO\parallel, f}\]
where s/f refer to the slower/faster material. This is satisfied for GaN/AlN, and for GaN/AlGaN if the Aluminium
composition is sufficiently high (see Komirenko1999). The class checks this assertion so it is safe to experiment.
In a given region, solutions are oscillating if \(\varepsilon_\perp\varepsilon_\parallel<0\)
and exponential if \(\varepsilon_\perp\varepsilon_\parallel>0\).
At an interface, the derivative switches signs iff \(\varepsilon_{1\parallel}\varepsilon_{2\parallel}<0\).
We will use the following convenient definitions, similar to the notation of Komirenko1999
but for a factor of two in \(\alpha\)
\[\xi_{i}=\sqrt{\abs{\varepsilon_{i\perp}\varepsilon_{i\parallel}}},\quad
\alpha_i=\sqrt{\abs{\varepsilon_{i\perp}/\varepsilon_{i\parallel}}}\]
Then the vertical wavevector of a mode in a given region is \(k_i=q\alpha_i\).
If the solution is written in Region 1 with a normalization constant \(A\) and Region 2 with a normalization constant \(B\), then the first matching condition gives us \(A/B\), and the normalization condition will be written
\[\begin{equation}
A^2=\frac{\hbar}{2\omega}\big/\left[ \beta^2_{\parallel 1}\gamma_{\parallel 1}^2+\beta^2_{\perp 1}\gamma_{\perp 1}^2 + \left(\frac{B}{A}\right)^2\left(\beta^2_{\parallel 2}\gamma_{\parallel 2}^2+\beta^2_{\perp 2}\gamma_{\perp 2}^2 \right) \right]
\end{equation}\]
with
\[\begin{equation}
\beta_{\parallel i}^2=\varepsilon^\infty_i\left( \omega_{LOi}^2-\omega_{TOi}^2 \right)\left( \frac{k_i}{\omega_{TO\parallel i}^2-\omega^2} \right)^2
\end{equation}\]
\[\begin{equation}
\beta_{\perp i}^2=\varepsilon^\infty_i\left( \omega_{LOi}^2-\omega_{TOi}^2 \right)\left( \frac{ q}{\omega_{TO\perp i}^2-\omega^2} \right)^2
\end{equation}\]
and
\[\begin{equation}
\gamma_{\parallel i}^2=\int_{i} dz \left(\frac{\partial_z \phi}{A k_i}\right)^2, \quad
\gamma_{\perp i}^2=\int_{i} dz \left(\frac{\phi}{B}\right)^2
\end{equation}\]
Confined to Region 1
If the solution is oscillating in Region 1 and decaying in Region 2, we can write
\[\begin{align}
\phi_1=A\sin(k_1 z), \quad \phi_2=Be^{-k_2z}
\end{align}\]
Matching interface conditions gives
\[\begin{equation}
q=\frac{1}{\alpha_1 t_1}\left[ \tan^{-1}\left( \xi_1/\xi_2 \right) +\pi n \right]
\end{equation}\]
with \(B/A=\sin(k_1 t_1)e^{k_2t_1}\) and
\[\begin{equation}
\gamma_{\parallel 1}=\frac{1}{2}(t_1+\frac{1}{2k_1}\sin(2k_1t_1)), \quad
\gamma_{\perp 1}=\frac{1}{2}(t_1-\frac{1}{2k_1}\sin(2k_1t_1))
\end{equation}\]
\[\begin{equation}
\gamma_{\parallel 2}=\frac{1}{2k_2}e^{-2k_2t_1}, \quad
\gamma_{\perp 2}=\frac{1}{2k_2}e^{-2k_2t_1}
\end{equation}\]
Note that if the energy region adjoins \(\omega_{LO\perp 1}\), the above does admit one peculiar \(n=0\) mode which has a very small \(q\) for all \(\omega\). This mode thus does not contribute to in-plane scattering, but is terrible for the math, so I’d filter it out.
Confined to Interface
If the solution is decaying in both regions, we can write
\[\begin{align}
\phi_1=A\sinh(k_1 z), \quad \phi_2=Be^{-k_2z}
\end{align}\]
Matching interface conditions gives
\[\begin{equation}
q=\frac{1}{2\alpha t_1}\log\left[ \frac{\xi_2+\xi_1}{\xi_2-\xi_1} \right]
\end{equation}\]
with \(B/A=\sinh(k_1t_1)e^{k_2t_1}\) and
\[\begin{equation}
\gamma_{\parallel 1}=\frac{1}{2}\left(\frac{1}{2k_1}\sinh(2k_1t_1)+t_1\right), \quad
\gamma_{\perp 1}=\frac{1}{2}\left(\frac{1}{2k_1}\sinh(2k_1t_1)-t_1\right)
\end{equation}\]
\[\begin{equation}
\gamma_{\parallel 2}=\frac{1}{2k_2}e^{-2k_2t_1}, \quad
\gamma_{\perp 2}=\frac{1}{2k_2}e^{-2k_2t_1}
\end{equation}\]
For large \(q\) mode approaches an resonant energy of \(\omega_{IF}\) defined by \(\xi_1=\xi_2\).
Confined to Region 2
If the solution is decaying in Region 1 and oscillating in Region 2, we can write
\[\begin{align}
\phi_1=A\sinh(k_1 z), \quad \phi_2=B\sin(k_2 z+\theta)
\end{align}\]
Matching interface conditions gives
\[\begin{align}
\theta=\tan^{-1}\left( \frac{\xi_2}{\xi_1}\tanh(k_1t_1) \right)-k_2t_1
\end{align}\]
with \(B/A=\sinh(k_1t_1)/\sin(k_2t_1+\theta)\). The \(t_2\) thickness normalization gives \(k_2=\pi (n+1)/t_2\), so
\[\begin{equation}
q=\frac{\pi(n+1)}{\alpha_2t_2}
\end{equation}\]
Normalization is accounted for via
\[\begin{equation}
\gamma_{\parallel 1}=\frac{1}{2}\left(\frac{1}{2k_1}\sinh(2k_1t_1)+t_1\right), \quad
\gamma_{\perp 1}=\frac{1}{2}\left(\frac{1}{2k_1}\sinh(2k_1t_1)-t_1\right)
\end{equation}\]
\[\begin{equation}
\gamma_{\parallel 2}=\frac{t_2}{2}, \quad
\gamma_{\perp 2}=\frac{t_2}{2}
\end{equation}\]
