MQW Solver Physics

This page gives a detailed description of the physics of the MQW solver model.

MQW Solver Calculations

Internally, the gain of the MQW structure is calculated using the following algorithm:

1. Electronic band structure calculation: The electron band structure is calculated using the k.p method. The result of this step is the electronic wave functions and sub-bands in the conduction and valence bands of the quantum wells. These represent states confined in the quantum wells.
2. Fermi level calculation: This energy level controls the occupation probability for each confined state through the Fermi-Dirac distribution function. Its value depends on the band structure (calculated in step 1) and on the charge density (input parameter).
3. Stimulated and spontaneous emission calculation: These emissions depend on the occupation of quantum well levels (step 2) and the overlap of wave functions between conduction and valence bands. The stimulated emission corresponds to gain, while the spontaneous emission mainly contributes to optical loss (it also seeds lasing when gain becomes larger than all loss mechanisms). Negative gain is equivalent to absorption, which can be used in the simulations of EAMs based on the quantum confined Stark effect. Exciton effects can also be included in the absorption calculation for EAM simulations.

The rest of this page gives the details for each of these calculations.

Electronic Band Structure

The electronic band structure is calculated based on the \(\:k \cdot p\) method, suitable for most common III-V semiconductors with zincblende and wurtzite crystal structures. First, the Hamiltonian is constructed, using the the \(\:k \cdot p\) method and a finite difference discretization scheme, and then the Schrodinger equation in the form of the eigenvalue problem is solved for each transverse electronic wave vector to obtain the band structure and wave functions in the quantum wells.

The \(\:k \cdot p\) method is currently supported on three levels (orders), i.e. 4x4, 6x6, and 8x8. Under the lowest 4x4 order, the conduction band is assumed parabolic, with its Hamiltonian defined as:

$$ H_c = \frac{\hbar^2}{2m_e^*}\left(k_x^2+k_y^2+k_z^2\right) + V_e(z) + a_c\left(\epsilon_{xx}+\epsilon_{yy}+\epsilon_{zz}\right),$$

where \(m_e^*\) is the electron effective mass, \(k_x\) and \(k_y\) are the transverse wave vectors, \(k_z\) is the wave vector in the growth direction, \(a_c\) is the conduction band deformation potential, \(\epsilon_{ij}\) are the components of the strain tensor, and \(V_e(z)\) is the electronic quantum well potential in the growth direction including the external electric field. The parabolic conduction band approximation is justified in most commonly used III-V semiconductors, due to the relatively large band gaps that minimize conduction and valence band mixing. The wave function in the z direction can be separated from the wave function in the x and y directions and we can obtain the quantum well subbands n from the following Schrödinger equation

$$ \left(-\frac{\hbar^2}{2}\frac{\partial}{\partial z}\frac{1}{m_e^*(z)}\frac{\partial}{\partial z}+V_e(z)\right)\psi^n_c(z) = \left(E(k_{\rho})-\frac{\hbar^2k_{\rho}^2}{2m_e^*(z)}\right)\psi^n_c(z), $$

where we used \(k_z = -i\partial/\partial z\). We assume the x-y plane to be isotropic, so that it can be described in cylindrical coordinates with \(k_{\rho}\) only.

The valence band is non-parabolic with 4 bands: heavy hole and light hole bands, including spin. This method is sufficient for most common III-V zincblende semiconductors where the spin-orbit split-off band lies several hundred meV below the heavy and light hole bands, which is usually more than the quantum well depth. Under the 4x4 \(\:k \cdot p\)-order, the strained Hamiltonian in the \(\{|3/2,3/2\rangle,|3/2,1/2\rangle,|3/2,-1/2\rangle,|3/2,-3/2\rangle\}\) basis can be defined as

$$ H_v = - \begin{pmatrix} P_t+Q_t-V_h(z) & -S_t & R_t & 0 \\ -S_t^{\dagger} & P_t-Q_t-V_h(z) & 0 & R_t \\ R_t^{\dagger} & 0 & P_t-Q_t-V_h(z) & S_t \\ 0 & R_t^{\dagger} & S_t^{\dagger} & P_t+Q_t-V_h(z) \end{pmatrix}, $$

with \(V_h(z)\) the hole quantum well potential with an external electric field, \(P_t=P+P_{\epsilon}\), where \(P_{\epsilon}\) is the strain contribution, and similarly for \(Q_t\), \(R_t\), \(S_t\), and the basis vectors are defined as

$$ \begin{array}{l}{|\frac{3}{2},\frac{3}{2}\rangle=-\frac{1}{\sqrt{2}}|\left(X+iY\right)\uparrow\rangle} \\ {|\frac{3}{2},\frac{1}{2}\rangle=-\frac{1}{\sqrt{6}}|\left(X+iY\right)\downarrow\rangle+\sqrt{\frac{2}{3}}|Z\uparrow\rangle} \\ {|\frac{3}{2},-\frac{1}{2}\rangle=\frac{1}{\sqrt{6}}|\left(X-iY\right)\uparrow\rangle+\frac{\sqrt{2}}{3}|Z\downarrow\rangle} \\ {|\frac{3}{2},-\frac{3}{2}\rangle=\frac{1}{\sqrt{2}}|\left(X-iY\right)\downarrow\rangle} \end{array}$$

Explicit expressions for Hamiltonian matrix terms are

$$ \begin{array}{l}{P=\frac{\hbar^2\gamma_1}{2m_0}\left(k_x^2+k_y^2+k_z^2\right)} \\ {Q=\frac{\hbar^2\gamma_2}{2m_0}\left(k_x^2+k_y^2-2k_z^2\right)} \\ {R=-\frac{\hbar^2\gamma_2}{2m_0}\sqrt{3}\left(k_x^2-k_y^2\right)+i\frac{\hbar^2\gamma_3}{2m_0}2\sqrt{3}k_xk_y} \\ {S=\frac{\hbar^2\gamma_3}{2m_0}2\sqrt{3}\left(k_x-ik_y\right)k_z} \\ {P_{\epsilon}=-a_v\left(\epsilon_{xx}+\epsilon_{yy}+\epsilon_{zz}\right)} \\ {Q_{\epsilon}=-\frac{b}{2}\left(\epsilon_{xx}+\epsilon_{yy}-2\epsilon_{zz}\right)} \\ {R_{\epsilon}=\frac{b}{2}\sqrt{3}\left(\epsilon_{xx}-\epsilon_{yy}\right)-id\epsilon_{xy}} \\ {S_{\epsilon}=-d\left(\epsilon_{xz}-i\epsilon_{yz}\right)} \end{array}.$$

where \(a_v\), \(b\), \(d\) are the valence band deformation potentials and \(\gamma_i\) are Luttinger parameters. Using a basis transformation, the 4x4 valence band Hamiltonian can be block diagonalized and cast as two 2x2 Hamiltonians, termed upper and lower:

$$ H_v = \begin{pmatrix} H^U & 0 \\ 0 & H^L \end{pmatrix} = -\begin{pmatrix} P_t+Q_t & \tilde{R}_t & 0 & 0 \\ \tilde{R}_t^{\dagger} & P_t-Q_t & 0 & 0 \\ 0 & 0 & P_t-Q_t & \tilde{R}_t \\ 0 & 0 & \tilde{R}_t^{\dagger} & P_t+Q_t \end{pmatrix},$$

where \(\tilde{R}=|R_t|-i|S_t|\). Using the effective mass theory, we substitute \(k_z = -i\partial/\partial z\), discretize the z-direction using a suitable finite-difference scheme, and then solve the one dimensional Schroedinger equation at each \(k_{\rho}\), again assuming isotropy in the x-y plane (axial approximation). In the case of parabolic conduction band it was sufficient to solve only one eigenvalue problem, but for the non-parabolic valence band we have to solve an eigenvalue problem for each \(k_{\rho}\).

The 8x8 \(\:k \cdot p\) Hamiltonian, when formulated with the following eight basis functions

$$ |S\uparrow\rangle, |S\downarrow\rangle, |X\uparrow\rangle, |Y\uparrow\rangle, |Z\uparrow\rangle, |X\downarrow\rangle, |Y\downarrow\rangle, |Z\downarrow\rangle ,$$

takes the block-matrix form

$$H_{k \cdot p} = \begin{pmatrix}
H_{cc} & 0 & H_{cv} & 0 \\
0 & H_{cc} & 0 & H_{cv} \\
H_{vc} & 0 & H_{vv} & 0 \\
0 & H_{vc} & 0 & H_{vv} \\
\end{pmatrix} $$

The diagonal block of conduction band reads

$$H_{cc}= E_{c} + {k_{x}A}_{c2}k_{x} + {k_{y}A}_{c2}k_{y} + {k_{z}A}_{c1}k_{z}$$

\(E_{c}\) is the conduction band minimum of bulk material, and

$$A_{c1} = \frac{\hslash^{2}}{2m_{0}}\left( \frac{m_{0}}{m_{e,z}} - E_{P}\frac{E_{g} + 2\Delta_{2}}{\left( E_{g} + \Delta_{1} + \Delta_{2} \right)\left( E_{g} + 2\Delta_{2} \right) - 2\Delta_{3}^{2}} \right)$$

$$A_{c2} = \frac{\hslash^{2}}{2m_{0}}\left(\frac{m_{0}}{m_{e,xy}} - E_{P}\frac{\left( E_{g} + \Delta_{1} + \Delta_{2} \right)\left( E_{g} + \Delta_{2} \right) - \Delta_{3}^{2}}{E_{g}\left\lbrack \left( E_{g} + \Delta_{1} + \Delta_{2} \right)\left( E_{g} + 2\Delta_{2} \right) - 2\Delta_{3}^{2} \right\rbrack}\right)$$

where \(E_{P}\) is the energy parameter for optical matrix elements, and \(m_{e,xy}\) and \(m_{e,z}\) are electron effective masses in and out of the quantum well plane. \( \Delta_1, \Delta_2, \Delta_3 \) are related to the band splitting parameters through \(\Delta_{1} = \Delta_{cr}\) and \(\Delta_{2} = \Delta_{3} = \Delta_{so}/3\), where \(\Delta_{cr}\) is termed the crystal field splitting energy and \(\Delta_{so}\) is the spin-orbit coupling energy.

The off-diagonal blocks of \(H_{k \cdot p}\) reads \(H_{cv} = \begin{pmatrix}
iPk_{x} & iPk_{y} & iPk_{z} \\
\end{pmatrix}\ ,\ \)where \(P = \sqrt{\frac{ħ^{2}}{2m_{0}}E_{P}}\). When \(H_{cv} \) and \(H_{vc}\) (Hermitian conjugate to each other) are neglected, the conduction and valence bands become decoupled, and hence the k.p order reduces to 6x6.

The diagonal block of valence bands reads

$$H_{vv} = \begin{pmatrix}
{\overline{E}}_{v} + k_{x}L_{1}k_{x} + k_{y}M_{1}k_{y} + k_{z}M_{2}k_{z} & k_{x}N_{1}^{+}k_{y} + k_{y}N_{1}^{-}k_{x} & k_{x}N_{2}^{+}k_{z} + k_{z}N_{2}^{-}k_{x} \\
k_{y}N_{1}^{+}k_{x} + k_{x}N_{1}^{-}k_{y} & {\overline{E}}_{v} + k_{x}M_{1}k_{x} + k_{y}L_{1}k_{y} + k_{z}M_{2}k_{z} & k_{y}N_{2}^{+}k_{z} + k_{z}N_{2}^{-}k_{y} \\
k_{z}N_{2}^{+}k_{x} + k_{x}N_{2}^{-}k_{z} & k_{z}N_{2}^{+}k_{y} + k_{y}N_{2}^{-}k_{z} & {\overline{E}}_{v} + k_{x}M_{3}k_{x} + k_{y}M_{3}k_{y} + k_{z}L_{2}k_{z} \\
\end{pmatrix} $$

Here \({\overline{E}}_{v}\) stands for the average of valence band edges of the underlying bulk material. The valence band maximum \(E_{v}\) is related to \({\overline{E}}_{v}\) through the following relation:

$${\overline{E}}_{v} = E_{v}^{} - \max\left( \Delta_{1} + \Delta_{2},\ \ \frac{\Delta_{1} - \Delta_{2}}{2} + \sqrt{\left( \frac{\Delta_{1} - \Delta_{2}}{2} \right)^{2} + 2\Delta_{3}^{2}} \right)$$

The \(L,M,N\) parameters are computed differently between zincblende and wurtzite materials. In the zincblende case, they are computed from Luttinger parameters\(\ \gamma_{1}\) to \(\gamma_{3}\), while in the wurtzite case, they are computed from \(A_{1}\) to \(A_{6}\), which are effective mass parameters specific to wurtzite materials.

The strain Hamiltonian under the 8x8 basis takes the following block-diagonal form:

$$H_{\epsilon} = \begin{pmatrix}
H_{c}^{(\epsilon)} & 0 & 0 & 0 \\
0 & H_{c}^{(\epsilon)} & 0 & 0 \\
0 & 0 & H_{v}^{(\epsilon)} & 0 \\
0 & 0 & 0 & H_{v}^{(\epsilon)} \\
\end{pmatrix}$$

where \(H_{v}^{(\epsilon)}\) is a three-by-three matrix block

$$H_{v}^{(\epsilon)} = \begin{pmatrix}
l_{1}\varepsilon_{xx} + m_{1}\varepsilon_{yy} + m_{2}\varepsilon_{zz} & n_{1}\varepsilon_{xy} & n_{2}\varepsilon_{xz} \\
n_{1}\varepsilon_{xy} & m_{1}\varepsilon_{xx} + l_{1}\varepsilon_{yy} + m_{2}\varepsilon_{zz} & n_{2}\varepsilon_{yz} \\
n_{2}\varepsilon_{xz} & n_{2}\varepsilon_{yz} & m_{3}\left( \varepsilon_{xx} + \varepsilon_{yy} \right) + l_{2}\varepsilon_{zz} \\
\end{pmatrix}$$

and \(H_{c}^{(\epsilon)} = a_{z}\varepsilon_{zz} + a_{xy}(\varepsilon_{xx} + \varepsilon_{yy})\). For zincblende, the \(l,m,n\) parameters are computed from the valence band deformation potentials \(a_{v}\) and \(b\), whereas for wurtzite, they are computed from \(D_{1}\) to \(D_{6}\), which are deformation potentials specific to wurtzite materials. As for the deformation potentials of conduction band, \(a_{z} = a_{xy} = a_{c}\) for zincblende, while \(a_{z} = \left( a_{1} + D_{1} \right)\) and \(a_{xy} = \left( a_{2} + D_{2} \right)\) for wurtzite.

On top of the k-dependent \(H_{k \cdot p} \) and the strain contribution \( H_{\epsilon}  \), high order \(\:k \cdot p\) methods (including 6x6 and 8x8) also take into account the band splitting effects due to built-in crystal field and relativistic spin-orbit coupling. These effects are characterized by two respective parameters named "delta crystal" and "delta spin-orbit", which can be set through the material database. They result in a modification to the \(\:k \cdot p\) Hamiltonian in a k-independent manner.

Fermi Levels

Given the input charge density n=p, averaged over the thickness t as

$$ n=p=\frac{1}{t}\int_0^t n(z)dz, $$

and the band structure and wave functions, we can calculate the electron and hole quasi-Fermi levels in the quantum wells using

$$ n=\sum_l \int \rho_c^{2D}(E^l)f(E^l)dE^l $$

in the conduction band, where the summation is over subbands, \(\rho_c^{2D}\) is the subband density of states, and f is the Fermi-Dirac distribution function, and using

$$ p=\sum_l \sum_{\sigma} \sum_{k_x} \sum_{k_y} f(E(l,\sigma,k_x,k_y)) $$

in the valence band. Here the \(\sigma\) index goes over the upper and lower Hamiltonian when \(\:k \cdot p\)-order is 4x4; otherwise \(\sigma\) is neglected.

Stimulated and Spontaneous Emission

Now we can calculate the stimulated and spontaneous emission coefficients using the following equations

$$ g_{st}\left(\omega\right)=\frac{\pi e^2}{n_rc\varepsilon_0m_0^2\omega}\frac{2}{t}\sum_{\sigma}\sum_{n,m}\int_{0}^{\infty}{\frac{k_{\rho}dk_{\rho}}{2\pi}M_{nm}^\sigma(k_{\rho})\times\frac{\Gamma/2\pi}{{(\Gamma/2)}^2+\left(E_{hm}^{en}(k_{\rho})-\hbar\omega\right)^2}\left[f_c^n\left(k_{\rho}\right)-f_v^{\sigma m}(k_{\rho})\right]}, $$

$$ g_{sp}\left(\omega\right)=\frac{\pi e^2}{n_rc\varepsilon_0m_0^2\omega}\frac{2}{t}\sum_{\sigma}\sum_{n,m}\int_{0}^{\infty}{\frac{k_{\rho}dk_{\rho}}{2\pi}M_{nm}^\sigma(k_{\rho})\times\frac{\Gamma/2\pi}{{(\Gamma/2)}^2+\left(E_{hm}^{en}(k_{\rho})-\hbar\omega\right)^2}f_c^n\left(k_{\rho}\right)\left[1-f_v^{\sigma m}(k_{\rho})\right]}, $$

where \(n_r\) is the effective refractive index, n is the conduction band subband, m is the valence band subband, \(M_{nm}^{\sigma}\) is the momentum matrix element, \(\Gamma\) defines Lorentzian broadening of energy levels due to intraband scattering, and \(f_c\) and \(f_v\) are the conduction and valence band Fermi-Dirac distribution functions.

The complex index \(\bar{n}=\Delta n+i\kappa\), where \(\Delta n\) is the change in the input refractive index due to MQW gain (or absorption) and \(\kappa\) is the corresponding attenuation coefficient, can be obtained from the linear electric susceptibility

$$ \chi\left(\omega\right)=\frac{2}{\varepsilon_0t}\sum_{\sigma}\sum_{n,m}\int_{0}^{\infty}{\frac{k_{\rho}dk_{\rho}}{2\pi}M_{nm}^\sigma(k_{\rho})\times\frac{f_v^{\sigma m}\left(k_{\rho}\right)-f_c^n\left(k_{\rho}\right)}{E_{hm}^{en}\left(k_{\rho}\right)-\hbar\omega-i\frac{\Gamma}{2}}}, $$

where now \(M_{nm}^{\sigma}\) is the electric dipole matrix element instead of the momentum matrix element. The relation between the complex index and susceptibility is given by the following approximate relations:

$$ \bar{n} = \frac{1}{2n_r}\chi. $$

Please note that due to a different formulation in the complex susceptibility formula, where the dipole matrix element is used, compared to the gain formula, where the momentum matrix element is used, there may be a slight difference in the values of gain (or absorption) calculated from the complex susceptibility compared to the values calculated from the gain formula directly.

Exciton Effects

There are two methods for modeling exciton effects, the original method is the direct method and a newly added method in release 2023 R2.1 is the variational method. Variational method is an alternative for evaluating the binding energies of ground state excitons.

Direct Method

For modeling exciton effects, the Hamiltonian has an additional term for Coulomb interaction potential between an electron in the conduction band and a hole in the valence band [5]:

$$H_{\nu, \nu^{\prime}}=H_{c} \delta_{\nu, \nu^{\prime}}+H_{v}^{\nu, \nu^{\prime}}+V_{\text {Coul }} \delta_{\nu, \nu^{\prime}}$$

where \(H_c\) and \(H_{v}^{\nu , \nu^{\prime}}\) are the electron and hole Hamiltonians defined previously, \(\nu\) represents the label of the basis functions in the valence band {|3/2,3/2⟩,|3/2,1/2⟩,|3/2,−1/2⟩,|3/2,−3/2⟩}, and the Coulomb potential is given by

$$V_{\text {Coul }}=-\frac{e^{2}}{4 \pi \epsilon_{0} \epsilon} \frac{1}{\mid \boldsymbol{r}_{\boldsymbol{e}}-{\boldsymbol{r}_{\boldsymbol{h}} \mid}}$$

Now, the Schrödinger equation for the exciton becomes

$$\sum_{\nu^{\prime}} H_{\nu, \nu^{\prime}} \Psi_{\nu^{\prime}}^{X}\left(\boldsymbol{r}_{\boldsymbol{e}}, \boldsymbol{r}_{\boldsymbol{h}}\right)=E_{X} \Psi_{\nu^{\prime}}^{X}\left(\boldsymbol{r}_{\boldsymbol{e}}, \boldsymbol{r}_{\boldsymbol{h}}\right)$$

where the exciton wave function can be expanded in terms of the electron and hole wave functions \(f_{n \sigma}(\boldsymbol{k}, z_{e})\) and \(g_{m \nu}(\boldsymbol{k}, z_{h})\) as

$$\Psi_{\nu}^{X}\left(\boldsymbol{\rho}, z_{e}, z_{h}\right)=\sum_{n, m} \int \frac{d^{2} \boldsymbol{k}}{(2 \pi)^{2}} \phi_{n m}^{X}(\boldsymbol{k}) f_{n \sigma}\left(\boldsymbol{k}, z_{e}\right) g_{m \nu}\left(\boldsymbol{k}, z_{h}\right) \exp (i \boldsymbol{k} \cdot \boldsymbol{\rho})$$

where \(n\) and \(m\) are conduction and valence subband indices, respectively, \(k\), \(\rho\) are the 2D wave vector and coordinate in the plane of quantum wells, and z is the coordinate in the quantum wells thickness direction. The Schrödinger equation can now be formulated in the momentum space as

$$T_{n m}(\boldsymbol{k}) \phi_{n m}^{X}(\boldsymbol{k})+\sum_{n^{\prime}, m^{\prime}} \frac{d^{2} \boldsymbol{k}^{\prime}}{(2 \pi)^{2}} V_{n m n^{\prime} m^{\prime}}\left(\boldsymbol{k}, \boldsymbol{k}^{\prime}\right) \phi_{n^{\prime} m^{\prime}}^{X}\left(\boldsymbol{k}^{\prime}\right)=E_{X} \phi_{n m}^{X}(\boldsymbol{k})$$

where \(V_{n m n^{\prime} m^{\prime}}(\boldsymbol{k}, \boldsymbol{k}^{\prime})\) is the matrix element of the Coulomb potential and \(T_{n m}(\boldsymbol{k})\) is the kinetic energy term. The exciton wave functions at the output of the solver are expressed in this momentum basis. Applying the axial approximation in the plane of quantum wells due to usually low warping of the band structure in that plane:

$$
\begin{array}{l}
f_{n \sigma}\left(\boldsymbol{k}, z_{e}\right)=f_{n \sigma}\left(k, z_{e}\right) e^{-i \sigma \theta} \\
g_{m \nu}\left(\boldsymbol{k}, z_{h}\right)=g_{m \nu}\left(k, z_{h}\right) e^{-i \nu \theta} \\
T_{n m}(\boldsymbol{k})=T_{n m}(k) \\
V_{n m n^{\prime} m^{\prime}}\left(\boldsymbol{k}, \boldsymbol{k}^{\prime}\right)=\sum_{l} V_{n m n^{\prime} m^{\prime}}^{l}\left(k, k^{\prime}\right) e^{i l \theta}
\end{array}
$$

where \(\theta\) is the in-plane angle, the Schrödinger equation with excitons can be written as

$$T_{n m}(k) \phi_{n m}^{X l}(k)+\sum_{n^{\prime}, m^{\prime}} \int \frac{k^{\prime} d k}{2 \pi} V_{n m n^{\prime} m^{\prime}}^{l}\left(k, k^{\prime}\right) \phi_{n^{\prime} m^{\prime}}^{X l}\left(k^{\prime}\right)=E_{X l} \phi_{n m}^{X l}(k)$$

where l becomes a well-defined angular momentum quantum number.

The absorption coefficient (negative gain) including exciton effect is

$$\alpha(\omega)=\frac{2 \pi e^{2} \hbar}{n_{r} \epsilon_{0} m_{0}^{2} c t} \sum_{X, l, \nu, n, m} \frac{1}{E_{X l}} \int \frac{k d k}{2 \pi} M_{n m}^{X l \nu}(k) \frac{\frac{\Gamma}{2 \pi}}{\left(\frac{\Gamma}{2}\right)^{2}+\left(E_{X l}(k)-\hbar \omega\right)^{2}}$$

where \(M_{n m}^{X l \nu}(k)\) is the magnitude squared of the momentum matrix element in the dipole approximation (long optical wavelength compared to the material lattice constant) and the assumption is that the conduction band state is empty, while the valence band state is full, which is a reasonable assumption in reverse biased pin junctions used for electro-absorption modulator applications, where the quantum wells are almost fully depleted. Similarly, the full complex susceptibility can be written as

$$\chi(\omega)=\frac{2}{\epsilon_{0} t} \sum_{X, l, \nu, n, m} \int \frac{k d k}{2 \pi} M_{n m}^{X l \nu}(k) \frac{1}{E_{X l}(k)-\hbar \omega-\frac{i \Gamma}{2}}$$

where \(M_{n m}^{X l \nu}(k)\)  is now the magnitude square of the electric dipole matrix element in the dipole approximation.

Variational Method:

In summary, the 'best' estimate of the binding energy of the exciton is obtained by employing trial wavefunctions for the relative motion of the electron and hole, then invoking the variational principle. Therefore, we start with the Hamiltonian representing the interaction of the electron-hole as a sum of three terms:

$$H_{c\nu}=H_{e} + H_{h} + H_{X}$$

where \(H_{e}\) and \(H_{h}\) are the one-particle Hamiltonians appropriate to the conduction and valence bands, respectively. The third term, \(H_{X}\), represents the electron–hole interaction, and is composed of two terms. The first term corresponds to the kinetic energy of the relative motion of the electron and hole in the \(x\)–\(y\) plane, while the other represents the Coulombic potential energy.

$$H_{X} = \frac{\hbar^2\nabla_{\rho}^2}{2m_r} - \frac{e^2}{4\pi\varepsilon\sqrt{\boldsymbol{\rho}^2 + \left|z_e - z_h\right|^2}}$$

with in plane effective mass

$$\frac{1}{m_r} = \frac{1}{m_e^*} + \frac{1}{m_h^*}$$

We construct the Schrodinger equation for quantum well system and find the eigenfunctions \(\psi\) and eigenvalues \(E\), by employing an exciton wavefunction in the form:

$$\psi = \psi_e(z_e)\psi_h(z_h)\phi_{e-h}(\rho,\left|z_e - z_h\right|)$$

where \(\psi_e(z_e)\) and \(\psi_h(z_h)\) are the eigenfunctions of the one-particle Hamiltonians, and \(\phi_{e-h}\) represents the electron and hole relative motion. It will be the variational wave function employed to minimize the total energy \(E\) of the system. The wavefunction is a hydrogenic type, given by after using \(a = \left|z_e - z_h\right|\) [6]:

$$\phi_{e-h}(r\prime) = e^{-\displaystyle\frac{\left|r\prime\right|}{\lambda}} = e^{-\displaystyle\frac{\sqrt{\rho^{2} + \gamma^{2}a^{2}}}{\lambda}}$$

where \(\lambda\) is a variational parameter and will be used to minimize the total energy \(E\), and \(\gamma\) is the second variational parameter and represents the dimensionality of the system. The total exciton energy of the system is given by [6]:

$$E = \frac{\langle{\psi|H|\psi}\rangle}{\langle{\psi|\psi}\rangle} = \frac{\langle{\psi|H_e|\psi}\rangle + \langle{\psi|H_h|\psi}\rangle + \langle{\psi|H_X|\psi}\rangle}{\langle{\psi|\psi}\rangle} $$

Details of sweeping the variational parameters \(\lambda\) and \(\gamma\) to minimax the exciton binding energy can be found in [6].

The excitonic contribution to the complex susceptibility due to each pair of electron and hole sub-bands is then given by [7]:

$$\chi_{ex}(\omega) = \frac{1}{\epsilon_{0}t}\sum_{\sigma}\sum_{n,m}M_{nm}^{\sigma}(k_{\rho}=0)\times \frac{\left|\phi_{e-h}(\rho=0,a=0)\right|^2}{E_x - \hbar\omega - i\frac{\Gamma}{2}},$$

where \(M_{nm}^{\sigma}\) is in the 4x4 basis, which already includes the spin so prefactor 2 is not necessary. The continuum contribution is given by:

$$\chi_{\mathrm{cont}}(\omega) = \frac{1}{\epsilon_{0}t}\sum_{\sigma}\sum_{n,m}\int_{0}^{\infty}\frac{k_{\rho}dk_{\rho}}{2\pi}S_{2D}\left(E_{hm}^{en}(k_{\rho})\right)M_{nm}^{\sigma}(k_{\rho})\times\frac{1}{E_{hm}^{en}(k_{\rho}) - \hbar\omega - i\frac{\Gamma}{2}}$$

where \(M_{nm}^{\sigma}\) is in the 4x4 basis and \(S_{2D}\) is the Sommerfeld enhancement factor defined for the 2D case as [2]:

$$S_{2D}\left(E_{hm}^{en}(k_{\rho})\right) = \frac{2}{1 + e^{-2\pi/\sqrt{\displaystyle\frac{E_{hm}^{en}(k_{\rho}) - E_{hm}^{en}(0)}{R_y}}}}$$

and \(R_y\) is the Rydberg energy, given by:

$$R_{y} = -\frac{m_{r}e^{4}}{2\hbar^{2}(4\pi\varepsilon)^{2}}$$

Finally, the absorption spectrum is calculated as:

$$\alpha(\omega) = \frac{1}{2n_r}(\mathrm{Imag}(\chi_{\mathrm{cont}}) + \mathrm{Imag}(\chi_{ex}))$$

References

1. D. Ahn et al., J. Appl. Phys. 64, 4056 (1988).
2. S. L. Chuang, Physics of Optoelectronic Devices (1996).
3. Chuang, Shun-Lien, et al. “Exciton Green’s-Function Approach to Optical Absorption in a Quantum Well with an Applied Electric Field.” Physical Review B, vol. 43, no. 2, Jan. 1991, pp. 1500.
4. Vurgaftman et al., J. Appl. Phys., 89, 5815 (2001).
5. C. Y.-P. Chao et al., "Momentum-space solution of exciton excited states and heavy-hole-light-hole mixing in quantum wells," Phys. Rev. B, 48, 8210 (1993).
6. P. Harrison et al., Quantum wells, wires and dots: theoretical and computational physics of semiconductor nanostructures. John Wiley & Sons (2016).
7. Mares, P. J., and S. L. Chuang. “Modeling of Self‐electro‐optic‐effect Devices.” Journal of Applied Physics 74, no. 2 (July 15, 1993): 1388–1397.
8. V. A. Fonoberov and A. A. Balandin, "Excitonic properties of strained wurtzite and zinc-blende GaN/AlxGa1−xN quantum dots," Journal of Applied Physics, vol. 94, p. 7178–7186, 2003.
9. S. Birner, "Modeling of semiconductor nanostructures and semiconductor–electrolyte interfaces," doctoral thesis, 2011.