The multi-quantum well (MQW) solver is a 1D physics-based solver for calculating optical and electronic properties of multi-quantum well stacks. The solver returns the gain, or absorption, and spontaneous emission coefficients, as well as the electronic band diagram, band structure, and wave functions. The solver interface consists of script commands: buildmqwmaterial, mqwgain, mqwindex.

## Solver physics

The electronic band structure is calculated based on the \(4x4\:k \cdot p\) method, suitable for most common III-V semiconductors with zincblende crystal structure. First, the Hamiltonian is constructed, using the the \(4x4\:k \cdot p\) method and a finite difference discretization scheme, and then the Schroedinger 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 conduction band is parabolic, with the strained 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 Schroedinger 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 then the quantum well depth. 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}\).

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, where \(\sigma\) goes over the upper and lower Hamiltonian.

Now we can calculate the gain 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. $$

## Finite difference mesh

Finite difference mesh is uniform and can be defined in the smallest increment of 1 angstrom.

## Strain

Strain in each quantum well layer can be defined as biaxial with the diagonal strain tensor components given as

$$ \epsilon_{xx}=\epsilon_{yy}=\frac{a_0-a}{a}, $$

$$ \epsilon_{zz}=-\frac{2C_{12}}{C_{11}}\epsilon_{xx}, $$

where \(C_{11}\) and \(C_{12}\) are the elastic stiffness coefficients. This formulation covers one the most important strained system: that obtained by a quantum well pseudomorphically grown on a (001)-oriented substrate. This means \(R_{\epsilon}=S_{\epsilon}=0\).

## Boundary conditions and selection of bound states

There are two types of boundary conditions: hard wall and perfectly matched layer (PML). Hard wall boundary conditions force the wave function to zero at the left and right edge of the simulation domain. To reduce the error it is important to have a thick enough buffer layers at both ends of the simulation domain to ensure sufficient wave function decay. For example if the first and last layers are barriers, which is the usual case, it is important to be able to make them thick enough (e.g. 10 nm or more depending in the barrier height) regardless of the thickness of the inner barriers. PML boundary conditions are absorbing boundaries with finite thickness and they cause the wave function to decay when passing through them in order to minimize the back reflections.

Only the bound states in the quantum wells are included in the Fermi level and gain calculation. The bound state filtering can be done with user defined cut-off, and the result of the selection can be verified by inspecting the band structure and wave function outputs. In case of hard wall boundaries the cut-off applies to the wave function derivative at the boundaries, while in the case of PML boundaries it applies to the ratio of the PML to the simulation domain probability density. The cut-off should be much less than 1, typically \(10^{-2}-10^{-5}\).

## Units

Units for all input and output quantities are SI, except energy, which is in the units of eV, and unless otherwise noted in the documentation for the MQW script commands.

### References

- D. Ahn et al., J. Appl. Phys. 64, 4056 (1988)
- S. L. Chuang, Physics of Optoelectronic Devices
- Chuang, Phys. Rev. B, 43, 9649 (1991)
- Vurgaftman et al., J. Appl. Phys., 89, 5815 (2001)