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:

- Electronic band structure calculation: The electron band structure is calculated using the 4x4 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.
- 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).
- 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 \(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 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 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}\).

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

## 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

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.

## 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)