High accuracy numerical solutions for band structures in strained quantum well semiconductor optical amplifiers

Abstract

In this paper, we have calculated the band structure of strained quantum well (QW) semiconductor optical amplifiers (SOAs) by using plane wave expansion method (PWEM) and finite difference method (FDM), respectively. The difference between these two numerical methods is presented. First, the solution of Schrödinger’s equation in a conduction band for parabolic potential well is used to check the validity and accuracy of these two numerical methods. For the PWEM, its stability and computational speed are investigated as a function of the number of plane waves and the period of QW. For FDM, effects of mesh size and QW width on its accuracy and calculation time are discussed. Finally, we find that the computational speed of FDM generally is faster than that of PWEM. However, the PWEM is more efficient than the FDM when wider SOAs are needed to be calculated. Therefore, to obtain high accuracy and efficient numerical solutions for band structures, numerical methods should be selected depending on required accuracy, device structure and further applications.

Key words： semiconductor optical amplifier; quantum well devices; plane wave expansion method; finite difference method

Introduction

Based on technological advances, quantum well [1] even quantum dot semiconductor devices [2,3] (with scales down to nanometer range) have been fabricated and exhibit very good performance. The physical phenomena such as quantum effect, nonlinear effects taking place in such nano-structures require accurate models to characterize device operation performance and optimize structures. Additionally, before fabricating the semiconductor devices, numerical models are also necessary to predict the performance of device [4,5]. Recently, the semiconductor optical amplifiers (SOAs) of strained quantum well (QW) have been extensively investigated. SOA technology is mature, and commercial devices related to this technology are available for optical communication systems. Mathematical models are required to aid in the design of SOAs and predict their operational characteristics. Band structures in the active region of the SOA of quantum well should be considered firstly, because electronic properties of QW SOA in the active region are deduced from band structures directly, such as sub-band energies and wavefunctions, gain spectrum, spontaneous emission spectrum, and refractive index-change. Additionally, the SOA of QW can be optimized for different applications. For example, the SOAs are used for all-optical signals processing, such as wavelength conversion [6,7], all-optical logic [8], where large α-factor of SOA is required. Oppositely, the α-factor of SOA should be as small as possible in the scheme of return-to-zero (RZ) differential-phase-shift keying (DPSK) signal regeneration [9]. This is because the phase information could not be preserved due to strong self-phase modulation (SPM) effect if large α-factor SOA is used.

Fortunately, the α-factor of the SOA of QW also can be optimized by band engineering [10]. However, as described above, the band structures need to be calculated firstly. Therefore, suitable numerical model is very important. In this paper, first, a comprehensive model with the spin-orbit coupling is presented for solving band structures in the strained SOA of QW. Second, the difference between two numerical methods has been given. Finally, the accuracy and computational speed of these methods are compared.

Theoretical formulations

The band structure near the zone center of a direct bandgap semiconductor is described by a k.p method [11,12]. Hamiltonian in conduction band is characterized by a parabolic-band model:

(1)H^c(k,z)={\displaystyle \frac{{\hslash }^{\mathrm{2}}}{\mathrm{2}}}\left({\displaystyle \frac{k_x^{\mathrm{2}}+k_y^{\mathrm{2}}}{m_{n,t}}}+{\displaystyle \frac{k_z^{\mathrm{2}}}{m_{n,z}}}\right)+V_e(z)+\mathrm{2}{a}_c\left(\mathrm{1}-{\displaystyle \frac{C_{\mathrm{12}}}{C_{\mathrm{11}}}}\right)ϵ,

where k is wave vector, \hslash is Planck constant, m_n,t, and m_n,z are electron effective masses perpendicular (t) and parallel (z) to growth direction, respectively, V_e(z) is the potential energy of unstrained conduction band edge, a_c is a conduction band hydrostatic deformation potential, C₁₁ and C₁₂ are stiffness constants. ϵ is given as follows:

(2)ϵ={\displaystyle \frac{a_{\mathrm{0}}-a}{a}},

where a₀ and a are the lattice constants of quantum well material and barrier material, respectively.

Electron wavefunction in conduction band can be written as

(3){\mathrm{\Psi }}_{n,k}^{c\eta }={\displaystyle \frac{e^{\mathrm{i}{k}_t\rho }}{\sqrt{A}}}{\mathrm{\Phi }}_n(z;k_t)|S,\eta \rangle ,

where A is the area of the quantum well, k_t is wave number, \rho =\sqrt{x^{\mathrm{2}}+y^{\mathrm{2}}}, Φ_n is the envelope function of the nth conduction sub-band, |S,\eta \rangle is the basis operator of wavefunction.

Conduction band structure E_n^c(k_t) can be calculate by solving

(4)H^c{\mathrm{\Phi }}_n(k,z)=E_n^c(k){\mathrm{\Phi }}_n(k,z).

In a valence band, under an axial approximation, the 6×6 Hamiltonian can be simplified into a block-diagonal form [5,13]

(5)H(k,z)=\left[\begin{array}{cc}{H}_{\mathrm{3}\times \mathrm{3}}^U(k,z)& \mathrm{0}\\ \mathrm{0}& H_{\mathrm{3}\times \mathrm{3}}^L(k,z)\end{array}\right],

(6)H_{\mathrm{3}\times \mathrm{3}}^{\sigma }=-\left[\begin{array}{ccc}P+Q-V_h(z)& R_k\mp \mathrm{i}{S}_k& \sqrt{\mathrm{2}}{R}_k\pm {\displaystyle \frac{i}{\sqrt{\mathrm{2}}}}{S}_k\\ R_k\pm \mathrm{i}{S}_k& P-Q-V_h(z)& \sqrt{\mathrm{2}}Q\pm \mathrm{i}\sqrt{{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}}{S}_k\\ \sqrt{\mathrm{2}}{R}_k\mp {\displaystyle \frac{\mathrm{i}}{\sqrt{\mathrm{2}}}}{S}_k& \sqrt{\mathrm{2}}Q\mp \mathrm{i}\sqrt{{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}}{S}_k& P+\mathrm{\Delta }(z)-V_h(z)\end{array}\right],

where

(7)\begin{array}{l}P=P_k+P_{ϵ},\hfill \\ Q=Q_k+Q_{ϵ},\hfill \\ P_k={\displaystyle \frac{{\hslash }^{\mathrm{2}}}{\mathrm{2}{m}_{\mathrm{0}}}}{\gamma }_{\mathrm{1}}(K_t^{\mathrm{2}}+\mathrm{2}{K}_z^{\mathrm{2}}),\hfill \\ Q_k={\displaystyle \frac{{\hslash }^{\mathrm{2}}}{\mathrm{2}{m}_{\mathrm{0}}}}{\gamma }_{\mathrm{2}}(K_t^{\mathrm{2}}-\mathrm{2}{K}_z^{\mathrm{2}}),\hfill \\ R_k={\displaystyle \frac{{\hslash }^{\mathrm{2}}}{\mathrm{2}{m}_{\mathrm{0}}}}\sqrt{\mathrm{3}}{\displaystyle \frac{{\gamma }_{\mathrm{2}}+{\gamma }_{\mathrm{3}}}{\mathrm{2}}}{K}_t^{\mathrm{2}},\hfill \\ S_k={\displaystyle \frac{{\hslash }^{\mathrm{2}}}{\mathrm{2}{m}_{\mathrm{0}}}}\mathrm{2}\sqrt{\mathrm{3}}{\gamma }_{\mathrm{3}}{K}_t{K}_z,\hfill \\ P_{ϵ}=-a_x({ϵ}_{xx}+{ϵ}_{yy}+{ϵ}_{zz}),\hfill \\ Q_{ϵ}=-{\displaystyle \frac{b}{\mathrm{2}}}({ϵ}_{xx}+{ϵ}_{yy}-\mathrm{2}{ϵ}_{zz}),\hfill \end{array}

where σ in Eq. (6) indicates U and L corresponding with the upper and lower signs in Eq. (5), respectively. m₀ is the electron effective mass,γ₁, γ₂, γ₃ are Luttinger parameters, k_t^{\mathrm{2}}=k_x^{\mathrm{2}}+k_y^{\mathrm{2}}, k_z is interpreted as -\mathrm{i}\partial /\partial z and {ϵ}_{xx}, {ϵ}_{yy}, {ϵ}_{zz} are

(8)\begin{array}{l}{ϵ}_{xx}={ϵ}_{yy}={\displaystyle \frac{a_{\mathrm{0}}-a}{a}},\hfill \\ {ϵ}_{zz}=-{\displaystyle \frac{\mathrm{2}{C}_{\mathrm{12}}}{C_{\mathrm{11}}}}{ϵ}_{xx}.\hfill \end{array}

Hole wavefunction can be written as

(9){\mathrm{\Psi }}_{m,k}^{v,\sigma }={\displaystyle \frac{e^{\mathrm{i}k\rho }}{\sqrt{A}}}{\displaystyle \sum _{i=hh,lh,so}{g}_{m,i}^{\sigma }}(z,k)|u_i^{\sigma }\rangle .

The valence band structure can be obtained by solving

(10){\displaystyle \sum _{i=hh,lh,so}{H}_{\mathrm{3}\times \mathrm{3}}^{\sigma }}{g}_{m,i}^{\sigma }(z,k)=E_{\sigma ,m}^v{g}_{m,i}^{\sigma }(z,k).

Numerical simulation methods

Finite difference method

The finite difference method (FDM) is used to solve Eqs. (4) and (10). All differential operators in these equations can be written as A(z){\partial }^{\mathrm{2}}/\partial z^{\mathrm{2}}, B(z)\partial /\partial z and C(z). Hermitian prosperities of the Hamiltonian are obtained by the following finite difference formulas [4]:

(11){A(z){\displaystyle \frac{\partial g}{\partial z^{\mathrm{2}}}}|}_{z=z_i}\to {{\displaystyle \frac{\partial }{\partial z}}\left(A(z){\displaystyle \frac{\partial g}{\partial z}}\right)|}_{z=z_i}={\displaystyle \frac{A(z_{i-\mathrm{1}})+A(z_{i+\mathrm{1}})}{\mathrm{2}{(\mathrm{\Delta }z)}^{\mathrm{2}}}}g(z_{i+\mathrm{1}})-{\displaystyle \frac{A(z_{i-\mathrm{1}})+\mathrm{2}A(z_i)+A(z_{i+\mathrm{1}})}{\mathrm{2}{(\mathrm{\Delta }z)}^{\mathrm{2}}}}g(z_i)+{\displaystyle \frac{A(z_i)+A(z_{i-\mathrm{1}})}{\mathrm{2}{(\mathrm{\Delta }z)}^{\mathrm{2}}}}g(z_{i-\mathrm{1}}),

(12){B(z){\displaystyle \frac{\partial g}{\partial z}}|}_{z=z_i}\to {{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(B(z){\displaystyle \frac{\partial g}{\partial z}}+{\displaystyle \frac{\partial B(z)g}{\partial z}}\right)|}_{z=z_i}={\displaystyle \frac{B(z_{i+\mathrm{1}})+B(z_i)}{\mathrm{4}\mathrm{\Delta }z}}g(z_{i+\mathrm{1}})-{\displaystyle \frac{B(z_i)+B(z_{i-\mathrm{1}})}{\mathrm{4}\mathrm{\Delta }z}}g(z_{i-\mathrm{1}}),

where A(z) and B(z) are position-dependent parameters, Δz is step size. From these equations, we can find that the accuracy of the solution depends on the step size (Δz).

Plane wave expansion method

The other numerical method used in our numerical model is the plane wave expansion method (PWEM). We put the quantum well into an assumed large box with the well located at the center. For convenience, a periodic boundary condition is imposed. Orthogonal and normalized plane waves are introduced. The envelope function can be expanded as [5]

(13)\begin{array}{l}{g}_v^{m,\sigma }(k,z)={\displaystyle \sum _{p=-\infty }^{\infty }{a}_v^{m,\sigma }(k)|p\rangle },\\ |p\rangle ={\displaystyle \frac{\exp (\mathrm{i}\mathrm{2}p\mathrm{\pi }z/l)}{\sqrt{l}}},\end{array}

where |p\rangle is the plane wave, p is a positive integer and l is the period. Substituting Eq. (13) into Eq. (10), we can obtain

(14){\displaystyle \sum _{v=\mathrm{1}}^{\mathrm{3}}{\displaystyle \sum _{p=-\infty }^{\infty }\langle q\left|H_{\mathrm{3}\times \mathrm{3}}^{uv}\right|p\rangle }}{a}_m^{v,p}(k_t)=E_m(k_t)a_m^{v,p}(k_t).

The accuracy of the solution depends on the number of plane waves. If the number of plane wave is infinite, the obtained wavefunctions and corresponding eigenenergies are identical to the true ones. However, a finite plane waves still leads to very good accuracy.

Results and discussion

Parabolic potential well solution

The solution of Schrödinger’s equation for parabolic potential well is analytically known [11]. Therefore, the analytical solution of the parabolic potential well can be used to check the validity and the accuracy of the two numerical methods.

The eigenenergies of the parabolic well is shown in Fig. 1. In this figure, the black lines are wavefunctions of the lowest four sub-bands in the conduction band. The blue line indicates the parabolic potential. Note that, the ground state has an energy level \hslash \omega /2 at the absolute zero of temperature [14]. Above this, it has equal-spaced energy step \hslash \omega.

Fig.1 The lowest four confined eigenenergies of a finite parabolic quantum well, black lines from bottom to top: the lowest four confined eigenenergies, blue line: parabolic potential

Full size|PPT slide

Table 1 displays the energies of the lowest ten sub-bands within a parabolic quantum well. The total width of parabolic quantum is 24 nm. Mesh size is 0.1 nm (FDM). For PWEM, 60 plane waves are used. Clearly, both numerical methods have very high accuracy. Although as the levels of sub-bands in the conduction band approach the top of the barrier, the eigenenergies experience the finiteness of the potential, and hence some discrepancy with the infinite parabola solution arises. However, as seen the value of E_n/(2E₁) in Table 1, the magnitude of the error is not significant.

Tab.1 Convergence test of two numerical methods at г point

	FDM	PWEM	accurate solution
n	En/eV	En/(2E1)	En/eV	En/(2E1)	En/eV	En/(2E1)
1	0.118488	0.500	0.118169	0.500	0.118209	0.500
2	0.355527	1.500	0.354589	1.500	0.354629	1.500
3	0.592525	2.500	0.591009	2.500	0.591049	2.500
4	0.829482	3.500	0.827429	3.501	0.827469	3.500
5	1.066400	4.499	1.063849	4.501	1.063889	4.500
6	1.303277	5.499	1.300269	5.501	1.300309	5.500
7	1.540114	6.498	1.536689	6.502	1.536729	6.500
8	1.776910	7.498	1.773109	7.502	1.773149	7.500
9	2.013666	8.497	2.009529	8.502	2.009569	8.500
10	2.250382	9.496	2.245949	9.503	2.245989	9.500

The percentage error with respect to the analytical solutions for the eigenenergies of the lowest ten sub-bands is shown in Fig. 2. Results solved by PWEM are more accurate than those solved by finite-difference method in the case of the number of the plane wave with 60 and the mesh size with 0.1 nm. The errors of both two numerical results are quite small. The maximum of error with FDM is as small as 0.25%, which can be tolerated in real numerical model.

Fig.2 Percentage error in the first ten eigenenergies solved by two numerical methods

Full size|PPT slide

Figure 3 shows that calculation speed with two numerical methods as a function of parabolic quantum well period. We can find that the CPU time increases with the parabolic quantum well period increasing. It should be noted that, when the parabolic quantum well is larger than 13 nm, the calculated speed of FDM is slower than that of PWEM. This is because the impact of QW width on FDM is stronger than that on PWEM. Figure 3(a) shows the CPU time of the FDM as a function of mesh size and Fig. 3(b) depicts the CPU time of the PWEM as a function of the number of plane wave.

Fig.3 CPU time with two different numerical methods as a function of parabolic quantum well period

Full size|PPT slide

Rectangle potential well solution

In this section, we considered the SOA of QW with rectangle potential well, in which the quantum well material is In_0.53Ga_0.47As and the barrier material is In_0.73Ga_0.27As_0.58P_0.42. The substrate material is InP. The Hamilton in conduction is simpler than that in the valence band. The extremely accurate solution of the conduction band structure can be obtained by using Newton-Raphson method [11]. These solutions will also be used to check the validity and accuracy of our numerical methods. In this case, the width of quantum well is set to 12 nm. The barrier width is set to 12 nm.

Based on the PWEM, in Eq. (14), the valence Hamiltonian can be expressed as the sum of several generalized items such as k_{z}A(z)k_z, (B(z)k_z+k_{z}B(z))/\mathrm{2} and C(z), where A(z), B(z) and C(z) are the position-dependent parameters. Taking the center of the well at z=0, we can calculate \langle q\left|H_{\mu \nu }^{\sigma }\right|p\rangle based on the following formulations:

(15)\begin{array}{l}\langle q\left|k_{z}A(z)k_z\right|p\rangle ={\left({\displaystyle \frac{\mathrm{2}\mathrm{\pi }}{l}}\right)}^{\mathrm{2}}pq\left[{\delta }_{pq}{A}_b+(A_w-A_b){\displaystyle \frac{\sin \left[(p-q+\xi )\mathrm{\pi }{\displaystyle \frac{l_w}{l}}\right]}{(p-q+\xi )\mathrm{\pi }}}\right],\hfill \\ \langle q\left|(k_{z}B(z)+B(z)k_z)/\mathrm{2}\right|p\rangle ={\displaystyle \frac{\mathrm{2}\pi }{l}}{\displaystyle \frac{p+q}{\mathrm{2}}}\left[{\delta }_{pq}{B}_w+(B_w-B_b){\displaystyle \frac{\sin \left[(p-q+\xi )\mathrm{\pi }{\displaystyle \frac{l_w}{l}}\right]}{(p-q+\xi )\mathrm{\pi }}}\right],\hfill \\ \langle q\left|C(z)\right|p\rangle ={\delta }_{pq}{C}_b+(C_w-C_b){\displaystyle \frac{\sin \left[(p-q+\xi )\mathrm{\pi }{\displaystyle \frac{l_w}{l}}\right]}{(p-q+\xi )\mathrm{\pi }}}.\hfill \end{array}

Figure 4 shows the eigenenergies of the rectangular potential well in the conduction band. The first two wavefunctions are calculated by FDM (Fig. 4(a)) and the PWEM (Fig. 4(b)), respectively. The difference between two numerical methods is quite little.

Fig.4 Rectangular potential well and the first two wavefunctions in conduction band. (a) FDM method; (b) PWEM method (solid line: the first sub-band in the conduction, dotted line: the second sub-band)

Full size|PPT slide

The envelope functions of the top valence sub-band for the three quantum well structures are shown in Fig. 5. The results shown in Fig. 5(a) are solved by the FDM, and the results presented in Fig. 5(b) are solved by the PWEM. It can be found that the difference between two numerical method results is not significant.

Fig.5 Envelope wavefunction for HH1 sub-band of rectangular potential well with different numerical methods. HH is heavy hole sub-band, LH is light hole sub-band, SO is spin-orbit sub-band. (a) FDM method; (b) PWEM method

Full size|PPT slide

In Fig. 6, we show the valence sub-band structures for unstrained quantum well semiconductor as a function of the transverse wave vector k_t solved by two different numerical methods. For the PWEM, the results are obtained under the conduction that 60 plane-waves used for the expansion. The period is set to 24 nm. For FDM, the mesh size is set to 0.1 nm. As it is shown in Fig. 6(a), when we solve energy band structures in the conduction band by using two numerical methods, it is found that the difference between two numerical method results is very little. However, when the valence band structures are calculated based on two numerical methods, the results from two methods are significantly different.

Fig.6 Energy dispersions in conduction band (a) and valence band (b) solved by two different numerical methods, inset (a1): the detail of conduction band

Full size|PPT slide

As shown in Fig. 7, the percentage error in the solution of conduction band eigenenergies is analyzed. For the PWEM, Fig. 7(a) shows the percentage error as a function of quantum well period. The quantum well period should be chosen large enough to ensure the envelope functions vanish at the boundary. Figure 7(b) shows the percentage error as a function of the number of plane waves. It can be seen that the percentage error of energy eigenenergies keep unchanged within 1×10^-4% for N larger than 30. Figure 7(c) shows the impact of mesh size on the accuracy based on the FDM. The percentage error will be augmented with the mesh size increasing.

Fig.7 Percentage error in the first eigenenergy (compared with analytic solution) as a function of quantum well period (a), the number of plane waves (b) and mesh size (c), respectively

Full size|PPT slide

Conclusion

In this paper, different methods used to solve band structures are compared. Considering the parabolic potential well, when the quantum well period is smaller than 13 nm, the FDM is faster than the PWEM. But the PWEM has higher calculation efficiency for solving wider quantum well period. The accuracy of the finite-difference method is affected by the mesh size while the number of functions and the quantum well period are the mainly considered parameters using the PWEM in the case of the rectangle potential well solution. Therefore, an efficient combination of the solution methods depends on the required accuracy, on the device structure and on further applications of the solution.

Acknowledgements

This paper was supported in part by the National Natural Science Foundation of China (Grant Nos. 61007042 and 60877056), the Doctoral Program Foundation of Institutions of Higher Education of China (No. 20090142110052), and the National Basic Research Program of China (No. 2011CB301704).

References

1 Huang X, Qin C, Huang D X, Zhang X L. Local carrier recovery acceleration in quantum well semiconductor optical amplifiers. IEEE Journal of Quantum Electronics, 2010, 46(10): 1047-1013

2 Meuer C, Schmidt-Langhorst C, Schmeckebier H, Fiol G, Arsenijević D, Schubert C, Bimberg D. 40 Gb/s wavelength conversion via four-wave mixing in a quantum-dot semiconductor optical amplifier. Optics Express, 2011, 19(4): 3788-3798

PMID

3 Meuer C, Schmidt-Langhorst C, Bonk R, Schmeckebier H, Arsenijević D, Fiol G, Galperin A, Leuthold J, Schubert C, Bimberg D. 80 Gb/s wavelength conversion using a quantum-dot semiconductor optical amplifier and optical filtering. Optics Express, 2011, 19(6): 5134-5142

PMID

4 Tan G L, Xu J M. Modeling of gain, differential gain, index change, and linewidth enhancement factor for strain-compensated QW’s. IEEE Photonics Technology Letters, 1998, 10(10): 1386-1388

5 Chang C S, Chuang S L, Lien C S. Modeling of strained quantum-well lasers with spin-orbit coupling. IEEE Journal on Selected Topics in Quantum Electronics, 1995, 1(2): 218-229

6 Liu Y, Tangdiongga E, Li Z, de Waardt H, Koonen A M J, Khoe G D, Shu X W, Bennion I, Dorren H J S. Error-free 320-Gb/s all-optical wavelength conversion using a single semiconductor optical amplifier. Journal of Lightwave Technology, 2007, 25(1): 103-108

7 Liu Y, Tangdiongga E, Li Z, Zhang S X, de Waardt H, Khoe G D, Dorren H J S. Error-free all-optical wavelength conversion at 160 Gb/s using a semiconductor optical amplifier and an optical bandpass filter. Journal of Lightwave Technology, 2006, 24(1): 230-236

8 Dong J J, Zhang X L, Xu J, Huangg D X. 40 Gb/s all-optical logic NOR and OR gates using a semiconductor optical amplifier: experimental demonstration and theoretical analysis. Optics Communications, 2008, 281(6): 1710-1715

9 Mørk J, Mecozzi A. Response function for gain and refractive index dynamics in active semiconductor waveguides. Applied Physics Letters, 1994, 65(14): 1736-1738

10 Nielsen M L, Mørk J, Suzuki R, Sakaguchi J, Ueno Y. Experimental and theoretical investigation of the impact of ultra-fast carrier dynamics on high-speed SOA-based all-optical switches. Optics Express, 2006, 14(1): 331-347

PMID

11 Yamanaka T, Yoshikuni Y, Yokoyama K, Lui W, Seki S. Theoretical study on enhanced differential gain and extremely reduced linewidth enhancement factor in quantum-well lasers. IEEE Journal of Quantum Electronics, 1993, 29(6): 1609-1616

12 Harrison P. Quantum wells, wires, and dots. Theoretical and Computational Physics. West Sussex, UK: John Wiley & Sons, 2001

13 Chuang S L. Physics of Optoelectronic Devices. New York: Wiley-Interscience, 1995

14 Chao C Y P, Chuang S L. Spin-orbit-coupling effects on the valence-band structure of strained semiconductor quantum wells. Physical Review B: Condensed Matter and Materials Physics, 1992, 46(7): 4110-4122

PMID