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. Lumerical's MQW solver provides a set of script-based functions (buildmqwmaterial, mqwgain, mqwindex) for this purpose. Additionally, the MQW solver simulation object in FE IDE provides a GUI for running the MQW simulations.
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. $$
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.
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.
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)
- C. Y.-P. Chao et al., Phys. Rev. B, 48, 8210 (1993)