Modeling Electronic and Optical Properties of InAs/InP Quantum Dots

Abstract

A theoretical investigation of electronic properties of self-assembled InAs/InP quantum dots (QDs) is presented, utilizing a novel two-step modeling approach derived from a double-capping procedure following QD growth processes, a method pioneered in this study. The electronic band structure of the QD is calculated by the newly established accurate two-step method, i.e., the improved strain-dependent, eight-band k p method. The impact of various QD structural parameters (e.g., height, diameter, material composition, sublayer, and inter-layer spacer) on electronic states’ distribution and transition energies is investigated. Analysis of carrier dynamics within QDs includes intraband and interband transitions. The calculation of the carrier transitions between two atomic states, providing insights into optical gain or loss within QDs, is in terms of dipole matrix element, momentum matrix element, and oscillation strength, etc. In addition, the time-domain, traveling-wave method (i.e., rate equations coupled with traveling-wave equations) is used to investigate the optical properties of QD-based lasers. Several optical properties of the QD-based lasers are investigated, such as polarization, gain bandwidth, two-state lasing, etc. Based on the aforementioned method, our key findings include the optimization of carrier non-radiative intraband relaxation through sublayer manipulation, wavelength control through emission blue-shifting and gain bandwidth via variation of sublayer, polarization control of QDs photoluminescence via excited states’ transitions, and the enhancement of two-state lasing in InAs/InP QD lasers by thin inter-layer spacers. This review offers comprehensive insights into QDs electronic band structures and carrier dynamics, providing valuable guidance for optimizing QD-based lasers and their potential designs.

1. Introduction

Self-assembled quantum dots (QDs) have attracted much attention due to their atomic-like electronic properties [1,2,3]. Thus, QDs have been extensively studied and utilized as an active medium in various optical devices, including semiconductor lasers, and semiconductor optical amplifiers (SOAs) [4,5,6], etc. Many efforts have been made to investigate the properties of QDs, especially the InAs/InP QDs, which make the operation wavelengths of semiconductor devices fall into the C- and L-band (1.53–1.6 µm) of telecommunication [7,8,9,10,11,12]. Despite these advancements, several challenges and limitations persist, necessitating further exploration and innovation.

The Stranski–Krastanov (SK) growth mode, based on lattice-mismatched epitaxy, has been widely employed as the primary method for QD fabrication, enabling sufficient QD areal density and uniformity of QD size for QD emission [1,13,14]. However, further improvements are necessary to achieve higher-quality QD emissions, which can be widely used in optical communications. Particularly, this growth mode is still challenged to make the epitaxial growth InAs/InP QDs as a perfect material gain system, as the relatively smaller lattice mismatch leads to higher QD size dispersion and consequently lowering modal gain. Previous efforts have focused on various techniques, including adjusting barrier material compositions to control QD size [15,16,17,18], employing double-capping growth and post-growth methods to mitigate QD size dispersion in-plane [9,19,20,21], inserting GaAs or GaP inter-layers between InAs QDs and upper/lower barrier layers to regulate QD size [22,23,24,25], etc. However, these studies have often pursued specific configurations for ideal QD fabrication rather than a systematic exploration from a physical perspective. In the following sections, we will update the theoretical methods used to study QDs.

Quantum dots, being low-dimensional, quantum-confined semiconductor structures, have traditionally been studied using Bloch’s theorem, originally developed for studying the electronic properties of bulk crystals [26,27,28,29]. Three fundamental methods have emerged based on different approximations in treating the atomic potential. The first, the tight-binding (TB) method, assumes tight binding of valence electrons to atoms, considering perturbation potential from nearly adjacent atoms only [30,31]. The second and third are the pseudopotential method and the k p method, respectively, assuming loosely bound valence electrons, with perturbation potential empirically estimated from all atoms or represented by the effective mass of electrons [29]. These methods are often applied to study the energy state distribution of QDs, typically considering effects of strain and piezoelectricity for self-assembled QDs. Further, the eight-band k p method is more favored than the other two due to its balance between accuracy and computational complexity [32,33,34,35].

The conventional eight-band k p method demonstrates high accuracy in simulating quantum dot structures, particularly single QD structures, thereby providing a robust tool for the design of electro-optical devices based on QD gain media [36,37,38]. However, it has certain limitations that remain unresolved. First, the conventional eight-band k p method often relies on fundamental theoretical models that do not accommodate new and specific QD fabrication processes. Typically, this method is coupled with single-step strain analysis, treating the entire heterojunction as a unified structure. Consequently, the calculated strain profiles are symmetrical and uniformly distributed. In practical heterojunction fabrication, however, epitaxial processes are discontinuous, with techniques such as double capping procedures [8] and interrupted epitaxial growth [39] being commonly employed. These result in varied and non-uniform strain distributions that are highly dependent on the actual fabrication processes. Moreover, closely stacked QDs are increasingly preferred for enhanced optical confinement, and single mode operation, leading to a reduction in capping layer thickness from more than 30 nm to less than 10 nm [40,41]. Given this context, both strain coupling and electronic coupling become significant for multi-layer stacked QDs [42,43]. To accurately analyze strain profiles and electronic coupling effect, the eight-band k p model has been implemented in conjunction with a two-step strain analysis method, thereby improving the simulation accuracy for QD structures.

Additionally, to comprehensively understand the optical properties of QDs and optimize QD-based devices, theoretical analysis of photon and phonon dynamics is essential. QD-based lasers have been a subject of extensive research and development in recent decades [27,44,45]. The QD-based lasers have demonstrated superior performance compared with conventional quantum well lasers, showing lower threshold currents, enhanced temperature insensitivity, and other favorable characteristics, all stemming from improved carrier engineering within the QDs [44,45,46,47]. The study of carrier intraband and interband transitions within QDs has also been ongoing. Typically, the calculation of interband optical transition based on Einstein’s description of the strength of absorption and emission of atoms, it is known as Fermi’s Golden Rule, has been widely used in treatments of optical transitions between bands of states in semiconductors [6,13,48,49] Based on the above theories, two common approaches for analyzing carrier dynamics in QD-based lasers are the frequency domain approach and the time domain traveling-wave method (TDTW), both relying on rate equations to describe carrier dynamics within and outside QDs [50,51,52,53,54]. These rate equations can incorporate various factors influenced by QDs, including homogeneous [55,56] and inhomogeneous gain broadening [57], Auger recombination [58,59], and Langevin noise [53,54], etc. The rate equations method has been widely used in theoretical studies of QD-based lasers due to its flexibility. However, modifications for basic rate equations may be necessary when describing carrier dynamics of lasers based on specific QDs, such as tunneling injection QD lasers [60]. In brief, the investigation of photon and phonon dynamics in QDs demands increasingly accurate theoretical methods.

In this review article, our primary focus lies on the theoretical investigation of InAs/InP QDs and QD-based lasers. The organization of this work is as follows: in Section 2, we introduce an accurate method for modeling the structure of InAs/InP QDs, which includes a two-step strain analysis. Utilizing this method, the calculation results demonstrate higher accuracy compared to other reported methods, particularly for QDs grown using the double-capping procedure. In Section 3, the investigation of carrier dynamics within QDs is given. We approach the optical properties of QDs from a quantum mechanical perspective, where absorption or gain is viewed as the loss or gain of photons resulting from electron transitions between atomic states. By analyzing carrier transitions between these states in terms of calculation of dipole matrix elements, momentum matrix elements, and oscillation strengths, etc., we comprehensively study the effects of various fabrication parameters on the optical properties, for example, the dot height, diameter, material composition, sublayer, and inter-layer spacer etc. Based on the calculation results, we propose a series of high performances QDs laser designs including high-speed QD devices, emission blue-shifting, gain bandwidth control, polarization control, and enhancing two-state lasing. Finally, Section 4 summarizes this work.

2. Theoretical Models and Methods

The modeling method is illustrated schematically in Figure 1, which commences with the definition of a 3D QD model and then analysis through two-step strain and piezoelectricity calculations. These calculated strain and piezoelectricity values are then input into the eight-band k p Hamiltonian. Subsequently, the single particle states are determined by solving the Schrödinger equation. Additionally, carrier transitions within QDs are computed using the transition momentum matrix element derived from electronic band structures. This enables the determination of optical parameters such as carrier lifetimes, transition rates, and transition polarizations. Furthermore, the carrier dynamics of QDs is examined within the context of a typical QD-based Fabry–Perot (FP) laser. Finally, our modeling system allows for the estimation of optical properties such as photoluminescence (PL) and electroluminescence (EL).

2.1. Modeling of QDs

2.1.1. Calculation of One Particle States

As is well known, self-assembled quantum dots are formed due to strain accumulation between different materials during growth. Therefore, accurately understanding the strain distribution is crucial for studying QDs. Strain significantly impacts electrical confinement, which can be seen by the band offset caused by chemical composition variations at heterojunctions, and substantially influences electron energies and wavefunctions. Numerous studies have employed conventional continuum mechanical (CM) elasticity, often combined with the k p method, to analyze strain distribution [26,36]. In addition, the valence force field (VFF) elasticity method is mainly used in the pseudopotential and tight-binding methods [30,61]. Generally, to improve the accuracy of pseudopotential or tight-binding method of QD, the strain should be calculated carefully by the VFF method [35,62]. For example, the interaction of second or the nearest neighbor atoms and spin-orbit coupling have been taken into consideration [63]. Besides, to improve the accuracy of the k p method of QD, the VFF method is used to evaluate the strain profile instead of CM method [64], and both the strain and piezoelectric effect are taken into the total Hamiltonian [65], etc. However, for QDs, the conventional strain calculation methods may lack the necessary accuracy to describe strain distribution of QDs, especially considering various growth techniques like the double-capping procedure (i.e., alter the strain profile of QDs) used to control QD size distribution [66,67,68,69]. Hence, more precise strain calculation methods are required. To address this, we employ a two-step strain calculation method to study the strain distribution of InAs/InP QDs typically grown through the double-capping procedure [38].

In conventional strain analysis, the epitaxial self-assembled QDs arise from the lattice mismatch between the deposited InAs and the barrier InGaAsP, and the in-plane strain is defined as

$${ϵ}_0={\displaystyle \frac{a_{InAs}-a_{InGaAsP}}{a_{InAs}}}$$

(1)

Here, $a_{InAs}$ and $a_{InGaAsP}$ are the intrinsic lattice constants of InAs used in QD and the quaternary alloy $InGaAsP$ used in barrier layers, respectively. The initial strain is taken simply in all three dimensions as

$${ϵ}_{0xx}={ϵ}_{0yy}={ϵ}_{0zz}={ϵ}_0$$

(2)

$${\mathit{ϵ}}_{\mathbf{0}}=\left(\begin{array}{ccc}{ϵ}_{0xx}& 0& 0\\ 0& {ϵ}_{0yy}& 0\\ 0& 0& {ϵ}_{0zz}\end{array}\right)$$

(3)

In contrast to the above description of strain, the two-step strain, following the double-capping growth procedure, is described as follows. Initially, an epitaxial InAs film is deposited to a critical thickness under biaxial compression, leading to the random formation of quantum dots on the surface. At this stage, only in-plane strain is considered, with no initial strain accounted for in the z-direction.

$${\mathit{ϵ}}_{\mathbf{0}}=\left(\begin{array}{ccc}{ϵ}_{0xx}& 0& 0\\ 0& {ϵ}_{0yy}& 0\\ 0& 0& 0\end{array}\right)$$

(4)

In the subsequent step, a thin InGaAsP quaternary alloy is deposited as the first capping layer, partially covering the QDs then following a 30 s short growth interruption. This growth interruption process reduces the height of the QDs to a relatively uniform level due to the mixing and exchange of As/P flux [37]. The components diffusion process changes the QD geometry as well as the strain distribution. After the first-capping fabrication stabilizes, a second capping layer of the same material starts to grow and fully encapsulate the QDs. In this phase, the deformation field obtained from the initial step, along with the additional out-of-plane lattice mismatch between the fully capped layer (FCL) and the QDs, serves as the starting point for our model. Given the fabrication’s stability, no further in-plane initial strain is applied due to the minimal lattice mismatch in the x- and y-directions. The second capping layer does not affect the elastic strain analysis as its lattice constants align closely with those at the new top of the “shortened” QD. Consequently, only the extra initial strain in the z-direction is affected and considered in the calculations. The strain tensor resulting from this is described as follows,

$${\mathit{ϵ}}_{\mathbf{0}}^{\prime }=\left(\begin{array}{ccc}-{ϵ}_{0xx}^{\prime }& -{ϵ}_{0xy}^{\prime }& -{ϵ}_{0xz}^{\prime }\\ -{ϵ}_{0xy}^{\prime }& -{ϵ}_{0yy}^{\prime }& -{ϵ}_{0yz}^{\prime }\\ -{ϵ}_{0xz}^{\prime }& -{ϵ}_{0yz}^{\prime }& -{ϵ}_{0zz}^{\prime }+{ϵ}_{0zz}\end{array}\right)$$

(5)

where, ${ϵ}_0^{\prime }$ are the resulting strain tensors obtained from the first step strain calculation.

The piezoelectric effect results in a net charge distribution or electric polarization due to the loss of symmetry in severely strained nanostructures [28]. In the study of self-assembled QDs, piezoelectricity must be considered. Initially, only the linear piezoelectric component is included in QD modeling. However, later studies showed that the quadratic components are also significant, as their magnitudes are comparable to those of the linear terms [26,36,37,65]. Therefore, our QD model incorporates both linear and quadratic piezoelectric effects. The linear polarization component P₁, which arises from the calculated strain field, can be expressed as

$$P_1=\left(\begin{array}{c}{P}_x\\ P_y\\ P_z\end{array}\right)=\left(\begin{array}{ccc}2e_{14}& 0& 0\\ 0& 2e_{14}& 0\\ 0& 0& 2e_{14}\end{array}\right)\left(\begin{array}{c}{\epsilon }_{yz}\\ {\epsilon }_{xz}\\ {\epsilon }_{xy}\end{array}\right)$$

(6)

The quadratic polarization component P₂ is given by

$$P_2=2B_{114}\left(\begin{array}{c}{\epsilon }_{xx}{\epsilon }_{yz}\\ {\epsilon }_{yy}{\epsilon }_{xz}\\ {\epsilon }_{zz}{\epsilon }_{xy}\end{array}\right)+2B_{124}\left(\begin{array}{c}{\epsilon }_{yz}\left({\epsilon }_{yy}+{\epsilon }_{zz}\right)\\ {\epsilon }_{xz}\left({\epsilon }_{xx}+{\epsilon }_{zz}\right)\\ {\epsilon }_{xy}\left({\epsilon }_{xx}+{\epsilon }_{yy}\right)\end{array}\right)+4B_{156}\left(\begin{array}{c}{\epsilon }_{xz}{\epsilon }_{xy}\\ {\epsilon }_{yz}{\epsilon }_{xy}\\ {\epsilon }_{yz}{\epsilon }_{xz}\end{array}\right)$$

(7)

The piezoelectric charge density introduced by the total polarization can be derived as

$${\rho }_p\left(r\right)=-\nabla ·\left(P_1+P_2\right)$$

(8)

Then the piezoelectric potential V_p is obtained by solving Poisson’s equation:

$${\rho }_p\left(r\right)={\epsilon }_0\nabla ·\left({\epsilon }_s\left(r\right)\nabla V_p\left(r\right)\right)$$

(9)

where ε₀ and ε_s are the vacuum permittivity and relative static dielectric tensor, respectively.

The energy states and wavefunctions of bound electron and hole states can be computed by the eight-band k p method with non-periodic boundary conditions, which was developed firstly to describe electronic states in bulk materials [28,29]. Then, it was applied to heterostructures, quantum well, quantum wires and quantum dots by varying and refining the envelop function of this theory [26,70,71]. In our QD model, the total Hamiltonian includes the kinetic, the impact of two-step strain and piezoelectricity, which descripted as follows,

$$H=H_k+H_s+V$$

(10)

where the kinetic Hamilton $H_{k }$ is

$$H_k=\left(\begin{array}{cccccccc}A& 0& V^{†}& 0& \sqrt{3}V& -\sqrt{2}U& -U& \sqrt{2}{V}^{†}\\ 0& A& -\sqrt{2}U& -\sqrt{3}{V}^{†}& 0& -V& \sqrt{2}V& U\\ V& -\sqrt{2}U& -P+Q& -S^{†}& R& 0& \sqrt{3/2}S& -\sqrt{2}Q\\ 0& -\sqrt{3}V& -S& -P-Q& 0& R& -\sqrt{2}R& \sqrt{1/2}S\\ \sqrt{3}{V}^{†}& 0& R^{†}& 0& -P-Q& S^{†}& \sqrt{1/2}{S}^{†}& \sqrt{2}{R}^{†}\\ -\sqrt{2}U& -V^{†}& 0& R^{†}& S& -P+Q& \sqrt{2}Q& \sqrt{3/2}{S}^{†}\\ -U& \sqrt{2}{V}^{†}& \sqrt{3/2}{S}^{†}& -\sqrt{2}{R}^{†}& \sqrt{1/2}S& \sqrt{2}Q& -P-∆& 0\\ \sqrt{2}V& U& -\sqrt{2}Q& \sqrt{1/2}{S}^{†}& \sqrt{2}R& \sqrt{3/2}S& 0& -P-∆\end{array}\right)$$

(11)

where the components in $H_{k }$ are given as

$$A=E_c-{\displaystyle \frac{\hslash }{2m^{\ast }}}\left({\partial }_x^2+{\partial }_y^2+{\partial }_z^2\right)$$

$$P=-E_v^{av}-{\displaystyle \frac{{\hslash }^2}{2m_0}}{\gamma }_1\left({\partial }_x^2+{\partial }_y^2+{\partial }_z^2\right)$$

$$Q=-{\displaystyle \frac{{\hslash }^2}{2m_0}}{\gamma }_2\left({\partial }_x^2+{\partial }_y^2+{\partial }_z^2\right)$$

$$R=\sqrt{3}{\displaystyle \frac{{\hslash }^2}{2m_0}}\left[{\gamma }_2\left({\partial }_x^2-{\partial }_y^2\right)-2i{\gamma }_3{\partial }_x{\partial }_y\right]$$

$$S=-\sqrt{3}{\gamma }_3{\partial }_z\left({\partial }_x-i{\partial }_y\right)$$

$$U={\displaystyle \frac{-i}{3}}{P}_0{\partial }_z$$

$$V={\displaystyle \frac{-i}{\sqrt{6}}}{P}_0\left({\partial }_x-i{\partial }_y\right)$$

(12)

Here, the Luttinger parameters are

$${\gamma }_1={\gamma }_1^L-{\displaystyle \frac{E_g}{3E_g+\mathrm{∆}}}$$

$${\gamma }_2={\gamma }_2^L-{\displaystyle \frac{1}{2}}{\displaystyle \frac{E_g}{3E_g+\mathrm{∆}}}$$

$${\gamma }_3={\gamma }_3^L-{\displaystyle \frac{1}{2}}{\displaystyle \frac{E_g}{3E_g+\mathrm{∆}}}$$

$$E_p={\displaystyle \frac{2m_0{P}_0^2}{{\hslash }^2}}$$

(13)

The dagger mark † here is conjugation of the matrix elements. The matrix element P₀ can be obtained by band parameters $E_p$, which can be measured by experiments.

Then, the strain-dependent Hamiltonian $H_{s }$ is

$$H_s=\left(\begin{array}{cccccccc}{a}_c{ϵ}_h& 0& -v^{†}& 0& -\sqrt{3}v& \sqrt{2}u& u& -\sqrt{2}{v}^{†}\\ 0& a_c{ϵ}_h& \sqrt{2}u& \sqrt{3}{v}^{†}& 0& v& -\sqrt{2}v& -u\\ -v& \sqrt{2}u& p-q& s^{†}& r& 0& \sqrt{3/2}s& -\sqrt{2}q\\ 0& \sqrt{3}v& -s& -p-q& 0& r& -\sqrt{2}r& \sqrt{1/2}s\\ -\sqrt{3}{v}^{†}& 0& r^{†}& 0& -p-q& s^{†}& \sqrt{1/2}{s}^{†}& \sqrt{2}{r}^{†}\\ -\sqrt{2}u& v^{†}& 0& r^{†}& s& -p+q& \sqrt{2}q& \sqrt{3/2}{s}^{†}\\ -u& \sqrt{2}{v}^{†}& \sqrt{3/2}{s}^{†}& -\sqrt{2}{r}^{†}& \sqrt{1/2}s& \sqrt{2}q& -p& 0\\ -\sqrt{2}v& u& -\sqrt{2}q& \sqrt{1/2}{s}^{†}& \sqrt{2}r& \sqrt{3/2}s& 0& -p\end{array}\right)$$

(14)

The components in $H_{s }$ are given as

$$p=a_v{\epsilon }_h$$

$$q=b{\epsilon }_b$$

$$r={\displaystyle \frac{\sqrt{3}}{2}}b\left({\epsilon }_{xx}-i{\epsilon }_{yy}\right)-id{\epsilon }_{xy}$$

$$s=-d\left({\epsilon }_{xz}-i{\epsilon }_{yz}\right)$$

$$u={\displaystyle \frac{-i}{\sqrt{3}}}{P}_0\left({\epsilon }_{xz}{\partial }_x+{\epsilon }_{yz}{\partial }_y+{\epsilon }_{zz}{\partial }_z\right)$$

$$v={\displaystyle \frac{-i}{\sqrt{6}}}{P}_0\left[\left({\epsilon }_{xx}-i{\epsilon }_{xy}\right){\partial }_x+\left({\epsilon }_{xy}-i{\epsilon }_{yy}\right){\partial }_y+\left({\epsilon }_{xz}-i{\epsilon }_{yz}\right){\partial }_z\right]$$

(15)

where ${\epsilon }_h$ and ${\epsilon }_b$ are hydrostatic strain and biaxial strain, respectively, which are defined as

$${\epsilon }_h={\epsilon }_{xx}+{\epsilon }_{yy}+{\epsilon }_{zz}$$

$${\epsilon }_b={\epsilon }_{zz}-{\displaystyle \frac{{\epsilon }_{xx}+{\epsilon }_{yy}}{2}}$$

(16)

The piezoelectric component V is

$$\mathrm{V}=\left(\begin{array}{cccccccc}-V_P& 0& 0& 0& 0& 0& 0& 0\\ 0& -V_P& 0& 0& 0& 0& 0& 0\\ 0& 0& V_P& 0& 0& 0& 0& 0\\ 0& 0& 0& V_P& 0& 0& 0& 0\\ 0& 0& 0& 0& V_P& 0& 0& 0\\ 0& 0& 0& 0& 0& V_P& 0& 0\\ 0& 0& 0& 0& 0& 0& V_P& 0\\ 0& 0& 0& 0& 0& 0& 0& V_P\end{array}\right)$$

(17)

The main material parameters used in calculation including both room temperature (RT, 300 K) and low temperature (LT, 70 K) are given in

Table 1. Material parameters used in the calculations at RT/LT [72,73].

	Unit	Temp.	InAs	GaP	In1–xGaxAsyP1–y
a	Å	RT	6.0583	5.4505	5.6533xy + 6.0583(1 − x)y + 5.4505x(1 − y) + 5.8687(1 − x)(1 − y)
		LT	6.0584	5.4512	5.6525xy + 6.0584(1 − x)y + 5.4512x(1 − y) + 5.8688(1 − x)(1 − y)
ρ	kg/m3	RT/LT	5668	4130	5317.5xy + 5668(1 − x)y + 4130x(1 − y) + 4790(1 − x)(1 − y)
εr		RT	14.55	11.1	13.18xy + 14.55(1 − x)y + 11.1x(1 − y) + 12.35(1 − x)(1 − y)
		LT	14.55	10.86	12.4xy + 14.55(1 − x)y + 10.86x(1 − y) + 11.77(1 − x)(1 − y)
C11	GPa	RT	83.29	140.5	119xy + 83.29 (1 − x)y + 140.5x(1 − y) + 101.1(1 − x)(1 − y)
		LT	83.3	143.9	122.3xy + 83.3 (1 − x)y + 143.9x(1 − y) + 102.2(1 − x)(1 − y)
C12	GPa	RT	45.26	62.03	53.8xy + 45.26(1 − x)y + 62.03x(1 − y) + 56.1(1 − x)(1 − y)
		LT	45.3	65.2	57.1xy + 45.3(1 − x)y + 65.2x(1 − y) + 57.6(1 − x)(1 − y)
C44	GPa	RT	39.59	70.33	59.4xy + 39.59(1 − x)y + 70.33x(1 − y) + 44.2(1 − x)(1 − y)
		LT	39.6	70.14	60xy + 39.6(1 − x)y + 70.14x(1 − y) + 46(1 − x)(1 − y)
e14	C/m2	RT/LT	0.045	0.1	−0.18(1 − y)2 + 0.15(1 − y) − 0.05
B114	C/m2	RT/LT	−0.5	−0.7	−0.4xy − 0.5(1 − x)y − 0.7x(1 − y) − 1.1(1 − x)(1 − y)
B124	C/m2	RT/LT	−4.1	−2.2	−3.8xy–4.1(1–x)y–2.2x(1–y)–3.8(1–x)(1–y)
B156	C/m2	RT/LT	0.2	−0.7	−0.7xy + 0.2(1 − x)y − 0.7x(1 − y) − 0.5(1 − x)(1 − y)
me	m0	RT	0.0213	Barrier	0.0632xy + 0.0213(1 − x)y + 0.158x(1 − y) + 0.077(1 − x)(1 − y)
		LT	0.023	Barrier	0.0665xy + 0.023(1 − x)y + 0.17x(1 − y) + 0.08(1 − x)(1 − y)
γ1L		RT/LT	20.4	Barrier	7.1xy + 20.4(1 − x)y + 4.04x(1 − y) + 5.33(1 − x)(1 − y)
γ2L		RT/LT	8.3	Barrier	2.02xy + 8.3(1 − x)y + 0.53x(1 − y) + 1.57(1 − x)(1 − y)
γ3L		RT/LT	9.1	Barrier	2.91xy + 9.1(1 − x)y + 1.26x(1 − y) + 2.11(1 − x)(1 − y)
ac	eV	RT/LT	−10.2	−7.14	−7.6 + 1.6x + [0.85xy − 1(1 − x)y − 1.7x(1 − y) + 0.6(1 − x)(1 − y)]
av	eV	RT/LT	1	1.7	−0.85xy + 1(1 − x)y + 1.7x(1 − y) − 0.6(1 − x)(1 − y)
b	eV	RT/LT	−1.8	−1.7	−1.85xy − 1.8(1 − x)y − 1.7x(1 − y) − 1.7(1 − x)(1 − y)
d	eV	RT/LT	−3.6	−4.4	−5.1xy − 3.6(1 − x)y − 4.4x(1 − y) − 4.3(1 − x)(1 − y)
Ec	eV	RT	0.709 − 0.37y	2.43 − 0.37y	0.14(1 − y)2 + 0.46(1 − y) + 0.75
		LT	0.768 − 0.37y	2.565 − 0.37y	1.423 + 0.689x − 1.185y + 0.758x2 + 0.18y2 + 0.763xy − 1.14x2y − 0.845xy2 + 0.875x2y2
Ev	eV	RT/LT	0.35 − 0.37y	−0.33 − 0.37y	0
Δ	eV	RT	0.41	0.08	0.108(1 − y) + 0.33y − 0.06y2
		LT	0.38	0.08	0.341xy + 0.38(1 − x)y + 0.082x(1 − y) + 0.108(1 − x)(1 − y)
Ep	eV	RT/LT	21.5	Barrier	28.8xy + 21.5(1 − x)y + 31.4x(1 − y) + 20.7(1 − x)(1 − y)

.

2.1.2. Comparison of Conventional and Two-Step Strain Effects

Based on the above eight-band k p method, the 3D InAs/InGaAsP/InP QDs models are established in COMSOL 5.4 and 5.6. Further details of the models used can be found in References [60,62,63,64,65]. And the geometry of a single truncated pyramid InAs QD and the geometry of InAs QD stacks are depicted in Figure 2 [74]. The 3D models are built based on the fabricated QDs in those reported experimental works [8,9,75,76,77]. The InAs QD is deeply buried in the InGaAsP material, and the height and base diameter of the InGaAsP substrate, the InAs wetting layer (WL), the InGaAsP first capping layer (FCL) and the InGaAsP second capping layer (SCL) are set by the estimation results of fabricated QDs. All the parameters of our QD model can be found in References [38,75,78,79,80]

Using our model, we initially compare the strained confinement potential of the single QD model with conventional strain and two-step strain, as this is crucial for QD investigations. The confinement potential energy (CPE) is characterized by hydrostatic strain and biaxial strain. Hydrostatic strain primarily shifts the conduction band edge (CBE) more than the valence band edge (VBE), while biaxial strain widens the gap between the heavy hole (HH) and light hole (LH) band edges. In Figure 3, we present the distributions of hydrostatic strain and biaxial strain for a single InAs QD (buried deeply in InGaAsP (1.15 Q) with dimensions of 2.7 nm height and 30 nm diameter) along the z-axis direction through the QD center, and the strained energy band edges for electrons, heavy holes, and light holes along the (001) direction. Clearly, the two-step strain method results in increased and peaked hydrostatic strain in the lower barrier layer near the InAs and InGaAsP interface, differing from the smooth and consistent drop observed with the conventional strain (Figure 3a). Furthermore, the biaxial strain leads to a tilted top for the QD model with the two-step strain, similar to experimental strain profiles reported in the literature [81]. In contrast, a flat top is obtained from the QD model with the conventional strain. Consequently, the smaller hydrostatic strain and tilted biaxial strain of the QD model with the two-step strain result in a lower CBE and less modified CPE compared to the QD model with the conventional strain, as illustrated in Figure 3c.

Further, we investigate the strain distribution of multiple QD stacks using both the conventional strain and two-step strain analysis. Taking a five-layer QD stack as an example, we vary the spacer thickness between QDs within one column, setting it as 30 nm and 5 nm, to examine the strain profile among QD layers. In Figure 4a, we present the biaxial strain and hydrostatic strain of QD stacks with spacer thicknesses of 30 nm and 5 nm. Significant strain coupling occurs between QDs in column stacks when the spacer layer is thin (see the right panel of Figure 4a). Conversely, when the spacer layer is thick, QDs between different layers can be considered independent (see the left panel of Figure 4a). In multilayer QD stacks, the internal strain of a central QD is influenced by the strain from the QDs directly above and below it. The effect of strain distribution on neighboring QDs varies between the two-step strain analysis and conventional strain analysis. In conventional analysis, the strain profile is symmetric due to the symmetric strain distribution along the (001) direction in individual QDs. This results in a bidirectional accumulation of compressive strain with increasing spacer layer number (SLN), from the outermost (top and bottom) QDs to the innermost (central) QDs. Consequently, the innermost QD is the most strained across all layers, as illustrated in the left panel of Figure 4b. In contrast, the two-step analysis reveals a unidirectional accumulation of compressive strain from the bottom to the top QD, leading to an asymmetric strain distribution. This uneven strain growth, as evidenced by scanning tunneling microscope (STM) images of InAs/InP QD [24,74,81,82], results in a non-uniform size distribution, with the lower QD having a greater impact on the upper QD through compressive strain. This cumulative strain significantly affects the uppermost QD, as shown in the right panel of Figure 4b. The upward accumulation of strain may not fully relax in very thin barrier layers, causing surface roughness and halting QD formation after 4 or 5 stacked layers in the InAs/GaAs system [83].

Finally, we investigate the transition energy profile using both the conventional and two-step strain-based QD methods, and the results are compared. Here, the QD has two kinds of shape as described in Reference [38], i.e., the truncated pyramid shape (TP) as shown in Figure 2 and the flat lens shape (FL) as shown in Figure 1 of Reference [38]. For this simulation, QD parameters are based on measurements reported in References [8,66], with QD heights ranging from 1.5 to 2.8 nm and base diameters from 30 to 40 nm. We then compare the calculation results with photoluminescence measurement results at room temperature, as shown in Figure 5. The patterned areas indicate the tunable ranges of PL peak wavelength calculated from transitions between conduction and valence ground states (GS) at RT. Figure 5a, b illustrate comparison of results obtained with the two-step strain model for barrier materials of 1.1 Q and 1.15 Q, respectively. In these figures, the two-step strain model results are shown with green vertical stripes for TP QD and translucent blue with dotted borders for FL QD. These results align more closely with the experimental data (red dots) compared to those from the conventional strain model, which are represented by grey diagonal stripes for TP QD and translucent purple with solid borders for FL QD. The conventional model consistently predicts fundamental transition energies that are higher than the measured values, with the discrepancy widening as the thickness of the first capping layer increases. The experimental data are in close agreement with the lower and upper bounds of the tunable range shown in Figure 5a, b. This agreement supports the validity of our model, as corroborated by previous studies [8,66].

All the above and our previous investigations shown, our two-step strain analysis method is more accurate than the conventional strain analysis method. Therefore, it is reasonable to apply the two-step strain analysis method to QDs.

2.2. Modeling of Carrier Transitions of QDs

To analyze the optical properties of quantum dots, it is essential to understand the interaction between light and matter. From a quantum mechanical perspective, transitions between available energy states can be explained by the Einstein’s Golden Rule [49], which establishes the framework for calculating the probability rate of electron transitions by solving the time-dependent Schrödinger equation, considering the coherent superposition of wavefunctions of two states (e.g., ${\psi }_i$ and ${\psi }_j$ represent the initial and final electron’s envelope function). According to Bloch’s theorem, the wavefunction of a state in a QD takes a specific form [26,84]

$$\left|{\psi }_m\right.\rangle ={\displaystyle \sum }_M{\displaystyle \sum }_{f=1}^2{F}_M^f\left(r\right)\left|u_M^f\right.\rangle$$

(18)

For the $\mathit{k}·\mathit{p}$ method, the chosen basis functions $\left|{\psi }_m\right.\rangle$ are a set of eight electron and hole states functions including spin down and up Bloch waves from the lowest conduction band and highest valence band, and the other remote states are considered via perturbation theory. M in (18) runs over electron (E), heavy-hole (HH), light-hole (LH), and split-off (SO) hole states. $F_M^f\left(r\right)$ is the corresponding envelope functions and $\left|u_M^f\right.\rangle$ is the zone-center Bloch functions expressed as [84].

$$\begin{gathered} \left|u_E^1\right.\rangle =\left|s\rangle {\chi }_{\downarrow }\right. \\ \left|u_E^2\right.\rangle =\left|s\rangle {\chi }_{\uparrow }\right. \\ \left|u_{LH}^1\right.\rangle ={\displaystyle \frac{-i}{\sqrt{6}}}\left(\left|x\right.\rangle +i\left|y\right.\rangle \right){\chi }_{\downarrow }+i\sqrt{{\displaystyle \frac{2}{3}}}\left|z\rangle {\chi }_{\uparrow }\right. \\ \left|u_{HH}^1\right.\rangle ={\displaystyle \frac{i}{\sqrt{2}}}\left(\left|x\right.\rangle +i\left|y\right.\rangle \right){\chi }_{\uparrow } \\ \left|u_{HH}^2\right.\rangle ={\displaystyle \frac{-i}{\sqrt{2}}}\left(\left|x\right.\rangle -i\left|y\right.\rangle \right){\chi }_{\downarrow } \\ \left|u_{LH}^2\right.\rangle ={\displaystyle \frac{i}{\sqrt{6}}}\left(\left|x\right.\rangle -i\left|y\right.\rangle \right){\chi }_{\uparrow }+i\sqrt{{\displaystyle \frac{2}{3}}}\left|z\rangle {\chi }_{\downarrow }\right. \\ \left|u_{so}^1\right.\rangle ={\displaystyle \frac{-i}{\sqrt{3}}}\left(\left|x\right.\rangle -i\left|y\right.\rangle \right){\chi }_{\uparrow }+i\sqrt{{\displaystyle \frac{1}{3}}}\left|z\rangle {\chi }_{\downarrow }\right. \\ \left|u_{SO}^2\right.\rangle ={\displaystyle \frac{-i}{\sqrt{3}}}\left(\left|x\right.\rangle +i\left|y\right.\rangle \right){\chi }_{\downarrow }-i\sqrt{{\displaystyle \frac{1}{3}}}\left|z\rangle {\chi }_{\uparrow }\right. \end{gathered}$$

(19)

In the above equations, $\left|s\right.\rangle$ is the s-like function, and $\left|x\right.\rangle$, $\left|y\right.\rangle$, and $\left|z\right.\rangle$ are the p-like function.${\chi }_{\uparrow }$ and ${\chi }_{\downarrow }$ are the eigenspinors of the Pauli spin matrix. The starting point for calculation of transition rates is Fermi’s Golden Rule [49]. It uses the momentum matrix element which is directly related to the dipole matrix element obtained from the time-dependent Schrodinger equation. The electronic dipole matrix elements (i.e., dipole moments) can be calculated by the following equation [49]

$$M_{ij}=|\langle {\psi }_j|-er|{\psi }_i\rangle |$$

(20)

where e denotes electron energy and r is position. And ${\mathit{P}}_{ij}$ is the momentum matrix element expressed as

$${\mathit{P}}_{ij}=\langle {\psi }_j|\mathit{P}|{\psi }_i\rangle ={\displaystyle \sum }_{M,M^{\ast }}{\displaystyle \sum }_{f, f^{\ast }=1}^2\langle F_{j{M}^{\ast }}^{f{f}^{\ast }}\left(r\right)|F_{iM}^f\left(r\right)\rangle \langle u_{M^{\ast }}^{f^{\ast }}|\mathit{P}|u_M^f\rangle$$

(21)

Here, $\mathit{P}$ is the momentum operator, for a QD structure, which can be written as [39]

$$\begin{gathered} \mathit{e}·\mathit{P}={\displaystyle \frac{m_0}{2\hslash }}\mathit{e}·\left\{2\overrightarrow{P_t}+\left[Q_t\left({\displaystyle \frac{1}{i}}{\displaystyle \frac{\overrightarrow{d}}{dz}}\right)+({\displaystyle \frac{1}{i}}{\displaystyle \frac{\overleftarrow{d}}{dz}})\ast Q_t\right]+\left[R_t\left({\displaystyle \frac{1}{i}}{\displaystyle \frac{\overrightarrow{d}}{dx}}\right)+({\displaystyle \frac{1}{i}}{\displaystyle \frac{\overleftarrow{d}}{dx}})\ast R_t\right] \\ +\left[S_t\left({\displaystyle \frac{1}{i}}{\displaystyle \frac{\overrightarrow{d}}{dy}}\right)+({\displaystyle \frac{1}{i}}{\displaystyle \frac{\overleftarrow{d}}{dy}})\ast S_t\right]\right\} \end{gathered}$$

(22)

Here, $P_t=P+p$, $Q_t=Q+q$, $R_t=R+r$, and $S_t$ = S + s are matrices given by the Hamiltonian calculation process [39], and the arrows in above equation indicate the derivatives direction. The above equation is valid for both interband and intraband transition, however the first term $2\overrightarrow{P_t}$ mainly contributes to the interband transition. So, only the $P_t$ matrix is included in interband calculations. Therefore, the oscillation strength between a state in valence band and a state in conduction band can be calculated by [84]

$$f_s^{\mathit{k}}={\displaystyle \frac{2}{m_0}}{\displaystyle \frac{{\left|e_{\mathit{k}}·{\mathit{P}}_{ji}\right|}^2}{\left(E_j-E_i\right)}}$$

(23)

where $E_j$ and $E_i$ are the final j (i.e., valence band state with subscript j = v) and initial i (i.e., conduction band state with subscript i = c) state energy, respectively; $e_{\mathit{k}}$ is the direction unit vector of electric field with $\mathit{k}=x, y, \mathrm{and} z$. $f_s^{\mathit{x}}$ and $f_s^y$ determine the x-polarized and y-polarized transverse-electric (TE) mode interband transition, respectively, while $f_s^z$ determines the transverse-magnetic (TM) polarized interband transition. To have a direct comparison to the experimental results reported, a Gaussian broadening for each calculated transition peak at energy $E_j-E_i$ is considered

$$I^k\left(E_t\right)={\displaystyle \sum }_{j,i}{\displaystyle \frac{f_s^{\mathit{k}}}{\sqrt{2\pi }\sigma }}exp[-{\displaystyle \frac{{\left\{E_t-\left(E_j-E_i\right)\right\}}^2}{2{\sigma }^2}}]$$

(24)

where $\sigma$ is the standard deviation, for comparison, which can be set according to the experimental spectrum.

In spite of the interband transition above given, the intraband transition can also be investigated. The rate of spontaneous relaxation rate ${\gamma }_{spont}$ and lifetime ${\tau }_{spont}$ is given in terms of the transition dipole moments $M_{ij}$ [85]

$${\gamma }_{spont}={\displaystyle \frac{1}{{\tau }_{spont}}}={\displaystyle \frac{{\omega }^3{n}^3}{3\pi {ϵ}_0\hslash c^3}}{\left|M_{ij}\right|}^2$$

(25)

where the finial state j and initial state i are both in conduction band. The transition frequency is denoted by $\omega$, $n$ represents the refractive index, ${ϵ}_0$ is the vacuum permittivity, c is the speed of light in a vacuum, and $\hslash$ is Planck’s constant. The envelope functions for the ground state (GS) and excited state (ES) electrons are substituted into Equation (10) and subsequently into Equation (13) to determine the relaxation rate between the ES and GS.

Lastly, the PL emission spectrum of QDs can be calculated using the function [79]

$$g\left(\omega \right)={\displaystyle \sum }_{m=1}^N{\displaystyle \frac{{\Gamma }_{xy}\omega D_{GS,ES}{G}_m{N}_d{\left|d_{GS,ES}\right|}^2}{2cn{h}_{QD}{\epsilon }_0\hslash }}Re\left\{{\displaystyle \frac{1}{\Gamma +j\left(\omega -{\omega }_m\right)}}\right\}$$

(26)

In the calculation, it is assumed that the QDs are organized into N groups, where $m$ ranges from 1 to N, and $G_m$ represents the Gaussian distribution of the QD size. The transverse confinement factor is denoted as ${\mathsf{\Gamma }}_{xy}$, $D_{GS,ES}$ refers to the degeneracy of the ground state (GS) and the first excited state (ES), $N_d$ is the dot density per unit area, $d_{GS,ES}$ is the dipole matrix element for the GS and the first ES, and $\mathsf{\Gamma }$ denotes the dephasing rate of the interband transition.

3. Results and Discussion

Based on the aforementioned methods, our investigation focuses on exploring the electronic and optical characteristics of InAs/InP QDs and understanding how they are influenced by various QD parameters, such as QD height, sublayer thickness, and QD spacer thickness, etc. By utilizing these methods, we aim to clarify the impact of these parameters on the behavior of QDs, shedding light on their fundamental properties and potential designs. Furthermore, we employ rate equations to simulate the optical properties of InAs/InP QD-based FP lasers. By applying rate equations, we can effectively model the dynamics of carriers within the QDs, allowing us to predict and analyze the performance of QD-based FP lasers under different operating conditions. This comprehensive approach enables us to gain insights into the behavior of these devices and their potential for practical applications in optoelectronics and photonics.

3.1. Intraband Relaxation Impacted by GaAs/GaP Sublayer

Semiconductor QDs hold vast potential for applications in electronic devices, information processing, and nonlinear optics, as previously mentioned [86,87,88]. Extensive research has been dedicated to studying QD electronic properties, both theoretically and experimentally, leading to a comprehensive understanding of the basic physics governing electron dynamics within QDs [27,44,45,58]. One critical aspect of this physics is the energy relaxation rate of QDs, particularly the mechanism behind electrons’ intraband relaxation. Understanding intraband relaxation is essential as it is typically characterized by the intraband relaxation rate [45,89]. Current experimental techniques, such as pump-probe spectroscopy, along with theoretical investigations, have focused on elucidating the effect of QD size, shape, material system, and initial state dependence on intraband relaxation rate [89,90,91,92,93]. It has been observed that, in both colloidal and epitaxial QDs, the intraband relaxation rate decreases with increasing QD size [93,94]. Colloidal QDs generally exhibit faster intraband relaxation compared with epitaxial QDs due to their smaller sizes.

Several intraband relaxation mechanisms have been understood, including relaxation by one longitudinal optical (LO) or one longitudinal acoustic (LA) phonon, relaxation by multiple LO phonons, and relaxation by the interaction of one LO phonon and one LA phonon, resulting in a phonon with an energy equal to the sum (LO + LA) or the difference (LO-LA) of the individual phonon energies [27,85,95,96] However, it remains unclear which of these intraband relaxation mechanisms dominates in QDs, particularly those with sublayers. Therefore, there is a pressing need for theoretical studies to investigate the various intraband relaxation mechanisms and identify the dominant relaxation process. This is especially crucial for self-assembled QDs, such as those formed from significant lattice mismatch between substrate and deposit, as seen in the InAs/InP material system [38]. Tuning and optimizing the relaxation rate in such QD systems hold great potential for specific applications.

Using the QD methods described above, we investigate the effect of adding an ultrathin sublayer, either GaP or GaAs, on the intraband relaxation of InAs/InGaAsP/InP QDs. The detailed QD models employed in this study are outlined in Reference [75]. We theoretically investigate the intraband transition between the first excited state and the ground state. We compare QDs with GaAs and GaP sublayers, where GaAs has a smaller lattice mismatch with the InAs QD material than GaP. Figure 6a shows the strain distribution along the center of the QD in the (001) direction, with and without the sublayer. The inset provides a normalized view of the strain distribution in three-dimensional space. It is observed that the GaP sublayer induces the greatest deformation in the QD. Correspondingly, Figure 6b presents the piezoelectric potential ($V_p$) along the diagonal of the QD’s base plane in the (110) direction for QDs with and without GaP or GaAs sublayers. The inset offers a normalized depiction of the piezoelectric potential distribution in three-dimensional space. The GaP sublayer produces a larger piezoelectric potential compared to the same thickness of GaAs sublayer or no sublayer. Both Figure 6a,b demonstrate how the ultrathin sublayer affects the strain and piezoelectric potential of the QD. A greater lattice mismatch between the ultrathin sublayer and the QD material results in increased QD deformation and a higher piezoelectric potential.

To examine the impact of the GaP or GaAs sublayer on the electronic band structure in greater detail, Figure 7 presents the electron state energies as a function of sublayer thickness for QDs with QD with a height of h_qd = 1.5 or h_qd = 3 nm. In Figure 7a, for a QD with a height of 1.5 nm and a GaP sublayer, the energies of states E_e3, E_e4 and E_e5 overlap or become degenerate around a sublayer thickness of 1.0 monolayer (ML). When the QD height is increased to 3 nm, as shown in Figure 7b, the degeneracy of the excited states occurs around a sublayer thickness of 1.5 ML. Conversely, when using a GaAs sublayer, Figure 7c,d illustrate that the electron state energies are relatively unaffected by the sublayer thickness, with a larger QD height resulting in more degenerate ESs. This indicates that the GaP sublayer significantly influences the electron state energies and their degeneracy, whereas the GaAs sublayer has a minimal effect. Moreover, higher degeneracy of the first ES is observed only at specific GaP sublayer thicknesses.

Building upon the aforementioned analysis, we can derive the spontaneous transition rate and lifetime using Equations (8) and (13). Figure 8a illustrates calculated electron relaxation rate $1/{\tau }_{spont}$ from the first ES to the GS for QDs with either GaP or GaAs sublayers. Observing Figure 8a, it is evident that the electron relaxation rate decreases with increasing GaP sublayer thickness, while it remains relatively constant for GaAs sublayers. This behavior can be attributed to the dependency of the relaxation rate on the energy difference between the first ES and GS, as depicted in Figure 8b. In Figure 8b, the energy difference between the first ES and GS for QDs with GaP or GaAs sublayers is presented, providing insight into the trends observed in Figure 8a. Obviously, a turning point in the spontaneous relaxation rate with the GaP sublayer occurs at a specific sublayer thickness of 0.95 MLs, indicating the existence of a critical thickness. To elucidate this critical thickness phenomenon, we consider the impact of longitudinal optical phonon relaxation induced by strained phonons in the QD, as detailed in Reference [75]. An analysis of the strained phonons and the energy gap between the first ES and GS for QDs with a height of 1.5 nm and different GaP sublayer thicknesses reveals a significant drop in the relaxation rate beyond the critical sublayer thickness. This abrupt decrease in relaxation rate can be attributed to the involvement of one LO phonon in the relaxation process. Furthermore, when the QD height is increased to 3 nm (with the QD base remaining unchanged), the relaxation rate behavior is recalculated. Detailed results for this scenario can be found in Reference [78], indicating that the abrupt decrease in relaxation rate beyond the critical thickness is determined by the involvement of two LO phonons in the relaxation process.

Briefly, the influence of an ultrathin GaP or GaAs sublayer on intraband relaxation in InAs/InGaAsP/InP QDs has been investigated. The study reveals that the presence of a GaP sublayer notably alters the electron states, whereas a GaAs sublayer has minimal impact. With the GaP sublayer, all electron energy levels experience an increase, leading to a reduction in the energy difference between the first excited state and ground state. This phenomenon potentially results in more degenerate ESs, owing to the presence of the GaP sublayer. A critical thickness is identified for the GaP sublayer, below which the spontaneous intraband relaxation is governed solely by one-LO or two-LO phonon interactions. Furthermore, our findings suggest that the critical thickness of the sublayer plays a crucial role in the design and optimization of QD devices, particularly for high-speed applications.

3.2. Wavelength Blue-Shifting and Gain Spectral Bandwidth Impacted by GaP Sublayer

Semiconductor QD lasers emitting light within the 1.55 μm wavelength range (C-band) hold significant promise for high-speed and long-range fiber optical communication systems, notably dense wavelength division multiplexing (DWDM) technology [97,98]. Various strategies have been proposed to shift the entire gain spectrum of InP-based QDs to the C-band. These include modifying compositions of barrier materials [15,16,17] adjusting the average QD height through double-capping or post-growth intermixing techniques [8,19,20,37], and introducing ultrathin GaAs/GaP sublayers [9,22,99], etc. Experimental findings have demonstrated that growing an ultrathin GaAs/GaP sublayer between InAs QDs and the lower barrier layer effectively tunes the emission wavelength of InAs/InP QD lasers to the 1.55 μm range at room temperature [8,99,100]. In order to investigate the impact of GaAs/GaP on the emission wavelength, the QD model is established based on the parameters from References [8,100]. The details are given in Reference [79] Figure 9a illustrates the effect of the GaP sublayer on the room-temperature PL characteristics of the QDs. The calculated tunable ranges for QDs with a GaP sublayer, represented by the lower shaded region, align closely with the experimental PL peak wavelengths. Simulations indicate a blue-shift of approximately 70 nm upon inserting a mere 0.28 nm GaP sublayer, demonstrating the high efficiency of wavelength blue-shifting. Interestingly, the simulation results reveal that the PL peak wavelength blue-shift still occurs even without considering As/P exchange in the calculation. This suggests that the change induced by the sublayer exclusively affects the strain field or confinement potential profile. In Figure 9a, it is observed that the QD base diameter varies between 35–40 nm without the GaP sublayer, whereas with the GaP sublayer, it ranges from 30–35 nm. This observation aligns with the findings of the statistical distribution of lateral sizes of QDs in Reference [8]. The presence of the GaP sublayer appears to decrease the dispersion of dot base sizes and mitigate lateral overgrowth.

In Figure 9b, the calculated ultra-broadband tunable range of gain at RT is presented. Without the sublayer, it is evident that the maximum shifting range towards shorter peak wavelengths is approximately 130 nm. The peak wavelengths vary from 1512 to 1642 nm as the dot height is modified from 1.5 to 2.7 nm. Although the GaP sublayer thickness is varied as depicted in Figure 9a, it is apparent that a greater shifting range towards shorter peak wavelengths is achieved. With a GaP sublayer thickness of 1.09 monolayers, there is an additional shifting range of approximately 65 nm towards shorter peak wavelengths. This results in a total shifting range of 195 nm for the peak wavelengths, with the peak wavelength reaching as low as 1447 nm at a dot height of 1.5 nm, corresponding to a 50% enhancement.

A broadband optical gain can be achieved by employing multiple quantum dot active layers with varying thicknesses, a structure known as a chirped structure [101,102]. Conventional chirped structures have been based solely on variations in dot height [103]. However, our findings suggest that incorporating additional GaP sublayers into the structure can further enhance wavelength tuneability and blue-shifting, particularly advantageous in the case of InAs/InP self-assembled QDs. Based on the simulation results of RT PL peak wavelength for both cases (without and with a GaP sublayer) as shown in Figure 9, an improved chirped structure consisting of five non-identical dot layers can be constructed. Details of the QD structures are provided in Reference [79].

Subsequently, the gain spectra of conventional chirped QD structures and the improved QD structures are calculated and compared. Figure 10a displays the gain spectra for each layer in the conventional chirped structure, with five distinct curves showing peak wavelengths at 1512, 1547, 1580, 1611, and 1642 nm. The full width at half maximum (FWHM) of the total gain, indicated by the bold black solid line, is 184.6 nm, centered at a wavelength of 1580 nm. In contrast, Figure 10b shows the gain spectra for the enhanced structure, featuring peak wavelengths at 1447, 1483, 1535, 1580, and 1642 nm. The FWHM of the gain spectral bandwidth extends to 245.7 nm, reflecting a 30% improvement. Additionally, the peak wavelength is blue-shifted to around 1510 nm, aligning more favorably with the application requirements of the C-band.

In summary, the analysis indicates that the sublayer enhances quantum confinement rather than strengthening Ga-P bonds, resulting in smaller QD sizes. This introduces an additional degree of freedom for engineering the emission peak wavelength and the total gain spectral bandwidth by precisely controlling the sublayer thickness.

3.3. Polarization of QDs

In self-assembled quantum dots, optical transitions between conduction and valence band states exhibit an anisotropic polarization, typically with the TE mode being stronger than the TM mode in material gain and total emitted light [15]. This polarization anisotropy poses challenges for achieving polarization independence in semiconductor devices like QD-SOAs [13]. To address this, extensive researches have focused on managing polarization anisotropy through various methods, including modifying QD shapes [104,105], aspect ratios [47], injection current densities [106], stacking layer numbers (SLN) [4,25], stacking spacer thicknesses [5], indium composition profiles [17], and growth conditions [16]. Previously, theoretical studies suggested that achieving high TM polarization in QD stacks with a small SLN would be challenging due to low mixing between heavy-hole (HH) and light-hole (LH) bands. However, experimental observations contradicted this, showing high TM polarization in QD stacks with small SLNs [77]. To understand this phenomenon, the physic mechanism behind TM polarization in QD stacks with small SLNs is studied using the above QD methods. Specifically, a three-layer QD stack is analyzed to clear the details. The analysis includes examining the energy state distribution, wavefunctions of GS and the first ES energy states, and the mixing between HH and LH bands. This comprehensive analysis aims to clarify the polarization differences between GS- and ES-dominated photoluminescence (PL). The details of the considered QD models for our simulations are given in Reference [80].

Based on previous studies, the high TM polarization is usually attributed to the high HH and LH mixing, resulting from closely stacked QDs with larger SLNs. In order to understand the high TM polarization in QD stacks with small SLNs, the HH and LH mixing in QD stacks with small SLN (from 2 to 5) is investigated, the component C_M of the electron (E), heavy-hole (HH), light-hole (LH), and split-off (SO) hole states in the GS is calculated [105] ($C_M = ʃʃʃ\sum _{j=1}^2{\left|F_M^j\left(r\right)\right|}^2{d}^{3}r,$ where M notes E, HH, LH, and SO and $F_M^j\left(r\right)$ is the corresponding envelope functions in (6)). Figure 11a illustrates the variations of $C_{HH}$, $C_{LH}$, and $C_{SO}$ as a function of spacer layer number (SLN), ranging from 2 to 5. It shows that the contributions from light hole (LH) states increase, while those from heavy hole (HH) states decrease as the SLN rises from 2 to 5. This trend indicates that the LH states are not sufficiently high to drive strong transverse magnetic (TM) polarization transitions [104]. So, the high TM polarization of the closely stacked QDs (i.e., 4 nm spacer layer thickness) with small SLNs cannot be attributed to the HH and LH mixing.

Indeed, to understand the TM polarization transition, it is crucial to consider the wavefunction types of each energy state in both the GS and the first ES. Figure 11b provides insight into this by depicting the probability density isosurfaces (wavefunctions) of the lowest electron energy states included in the GS and the first ES for closely stacked three-layer QDs with a spacer thickness (T_s) of 4 nm. The s-type and p-type denote the types of orbitals associated with the bound electrons at each energy state. Figure 11b displays the probability shapes, illustrating the dominance of specific Bloch functions from Equation (7) among the eight Bloch functions. In the case of closely stacked three-layer QDs, the ground state and excited state energy levels include two s-type wavefunctions each, while the remaining six ES energy states are characterized by p-type wavefunctions. This results in an ES to GS degeneracy ratio of 4 for these QD stacks, indicating that the transitions between ES and GS have a significant impact in the closely stacked three-layer QD configurations.

To further understand the PL emission dominated by either the first ES or the GS transitions, the PL emission spectrum of closely stacked three-layer QDs is calculated and depicted in Figure 12a. This spectrum shows the PL of the closely stacked three-layer QDs, with the peak wavelengths of the ES and GS PL around 1557 nm and 1679 nm, respectively. Notably, the PL is largely dominated by the ES, given the higher degeneracy of the ESs compared to the GSs, as revealed in Figure 11a. Moving forward, the transition intensities of TE and TM polarizations for the closely stacked three-layer QDs are computed using Equations (8)–(12). The oscillation strength $f_s^k$ is calculated for all transitions between the 8~10 electron-states and 15 hole-states, including all GS and the first ES energy states. A standard deviation of 30 meV is set in Equation (12) to align with the PL emission calculation settings. The results of these calculations are presented in Figure 12b.

It is evident that the contributions of TE and TM polarizations are almost similar, exhibiting the same overall trend. Initially, the TM polarization slightly surpasses the TE polarization, but beyond a critical point, the TE polarization begins to slightly exceed the TM polarization, particularly in the vicinity of the peak region. Notably, the peak wavelength of the ES-dominated TE and TM polarization transitions in Figure 12b is approximately 1553 nm, which closely aligns with the peak wavelength of the PL spectrum in Figure 12a at around 1557 nm. This observation underscores that the TM polarization is predominantly contributed by the ES transitions in the closely stacked three-layer QDs. Thus, it can be concluded that the ES-dominated PL is specific to the closely stacked three-layer QD stacks, where the p-type wavefunction exists solely in the ES, indicating that the TM polarization is predominantly influenced by the ES transitions.

In summary, PL emission and polarization characteristics (TM/TE) in closely stacked quantum dot (QD) structures with minimal SLNs show that the pronounced TM polarization is mainly due to strong transitions from ES rather than the mixing of HH and LH states. This work offers a framework for understanding and predicting the high TM polarization observed in experimental studies of closely stacked QDs with small SLNs. This finding offers a promising avenue for achieving high TM polarization in QD-based devices, potentially paving the way for novel methods in device design and fabrication.

3.4. Two-State Lasing of QD Laser

Two-state lasing, characterized by the simultaneous emission of coherent light at two distinct wavelengths corresponding to the QD GS and ES, represents a unique capability exclusive to QD-based semiconductor lasers [107]. Consequently, two-state lasing has emerged as a distinctive characteristic of QD semiconductor lasers, which has garnered significant research attention over the past decade. On one hand, two-state QD lasers have been investigated for their applicability in all-optical memory components and self-mixing velocimetry [108]. On the other hand, the effects of two-state lasing have been extensively explored to enhance QD lasers for fiber-optical communications [109,110,111]. In particular, QD lasers operating ES emission demonstrate advantages over GS emission, including reduced relative intensity noise (RIN) and phase noise [112]. These lasers also exhibit a nearly zero linewidth enhancement factor (LEF) [113] and offer a wider modulation bandwidth [111]. However, two-state lasing is typically observed in InAs/GaAs QD lasers with either a single-section or two-section structure, not specifically optimized for long-haul optical transmission. Conversely, achieving stable two-state lasing experimentally in the InAs/InP system remains a significant challenge. Previous research demonstrated that two-state lasing in InP-based QDs was only attainable under optical injection while sole GS lasing observed under continuous wave (CW) electrical bias at RT [114]. This limitation was attributed to non-uniform carrier injection and poor hole mobility within the barrier layers. However, the ES emission was not observed even with p-doping in the barrier layers [115,116]. There have been theoretical studies on two-state lasing primarily focused on optical injection, so that the reasons leading to the absence of ES lasing in InAs/InP QDs are still unclear.

The two-state lasing behavior is commonly associated with the phonon-bottleneck effect [27]. This bottleneck effect can arise from both the random population effect at low temperatures and the presence of a large density of available states (i.e., high degeneracy) at the upper energy level, as clarified in our previous work [78]. The theoretical possibility of achieving high degeneracy of ES due to the reduction of spacer thickness between QD layers at RT has also been discussed. Details of these calculations can be found in the referenced works [38,78]. Specifically, we have calculated the six lowest electron energy levels (E_1–6) as a function of the spacer thickness (T_S) between QD layers, as illustrated in Figure 13. A notable observation is the significant energy splitting between binding states (solid lines) and anti-binding states (dotted lines) when dot layers are closely stacked, attributed to resonant electronic coupling. This splitting diminishes with increasing Ts, accompanied by a rise in the energy levels of the binding s-like (square marks) and p-like (circular marks) states, while the energy of the anti-binding counterparts decreases. Eventually, the splitting vanishes at T_S = 30 nm. Additionally, the anti-binding s-orbital (represented by the red dotted line with empty square marks) intersects with two binding p-orbitals (depicted by blue solid lines with circular marks) at approximately 10 nm, leading to a mixture of states or a “quasi-continuum band” at approximately 1.049 eV (assuming a double spin degeneracy for each state). This approximately six-fold degenerate excited state enhances the bottleneck effect, thereby playing a crucial role in the manifestation of two-state lasing at RT.

To clarify the impact of high degeneracy ES on the two-state lasing, the lasing behaviors in single-section FP QD lasers for the two specific cases (T_S = 10 and 30 nm) are calculated and compared. Here, the time domain traveling wave model is used, it is a simplified four-level, rate-equation model coupled with a pair of forward/backward propagating optical fields, and the entire details of this model can be found in References [78,118].

The simulated optical spectra, obtained at various levels of CW bias currents ranging from 1 to 6 times the threshold currents (I_th), are depicted in Figure 14a–d for both the T_S = 10 and 30 nm cases, considering inhomogeneous gain bandwidth (IGB) levels of full width at half maximum (FWHM) ΔE = 35 and 40 meV. The calculated threshold currents are 23, 15, 24, and 16 mA, respectively. It is obvious that the presence of two-state lasing (approximately ES: λ = 1560 nm and GS: λ = 1640 nm) is evident only in the 10 nm case (Figure 14a,c). It can be seen as evidence of the bottleneck effect caused by the increased ES to GS degeneracy ratio (i.e., 3:1). Conversely, sole GS lasing around 1600 nm is observed for the 30 nm case (Figure 2b), reproducing the electroluminescence (EL) results from prior studies [114]. This can be attributed to a relatively lower ES to GS degeneracy ratio (i.e., 2:1). Furthermore, two-state lasing at two different inhomogeneous gain bandwidths (Figure 14a,c) have distinct behaviors. For inhomogeneous gain bandwidth ΔE = 35 meV, the optical behavior is similar to the observations of several previous experiments [114,119,120] and the PL results in Reference [114]. GS lasing is achieved under the injection of I_th and then ES + GS lasing (patterned in red) is observed with two I_th injection currents (Figure 14a) and, finally, the peak power of ES emission starts to exceed GS when the injection current is larger than three I_th (patterned in purple). In contrast, for an inhomogeneous gain bandwidth ΔE = 40 meV, stable lasing starting on ES is presented (Figure 14c), and the peak power of the ES emission consistently surpasses that of GS. It is similar to the previous findings [121]. This reversal lasing might be related to the inhomogeneous gain bandwidth level, with a higher inhomogeneous gain bandwidth representing a larger dot-size dispersion, and it seems that GS is more adversely affected by the larger dot-size dispersion in providing a sufficient gain to overcome intra-cavity loss, reducing the probability of GS lasing.

To further clarify the contribution of emission from each state, it is essential to isolate ES and GS emissions and calculate the output power for each state separately. The calculated power-current characteristics for both the T_S = 10 and 30 nm cases with inhomogeneous gain bandwidth ΔE = 35 and 40 meV are depicted in Figure 15. This can be achieved by filtering the output optical field of stable ES + GS emissions at the uncoated facet for each injection level, employing Hanning windows of suitable width (e.g., eight times the free spectral range (FSR)) centered at ES and GS lasing frequencies, followed by time-averaging of the state-resolved temporal power $P_m\left(t\right)=\left(1-r_L^2\right){\left|E_m^{+}\left(L,t\right)\right|}^2$ for m = ES, and GS. Compared to the 10 nm case with smaller IGB (blue line with marks in the lower inset), the larger dot-size dispersion scenario (red line with marks in the upper right inset) exhibits much smaller I_th for ES lasing (25 mA versus 42 mA) and suppressed P_GS (almost becoming half under approximately 3 I_th injection), suggesting that QDs with larger ΔE are more favorable for ES lasing.

In summary, the exclusive occurrence of two-state lasing in the 10 nm case underscores the importance of the bottleneck effect, which amplifies the ratio of ES to GS degeneracy. Furthermore, a modest increase in the inhomogeneous gain bandwidth facilitates the initiation of preferred ES lasing. Most importantly, the utilization of small spacer thickness presents a novel avenue to attain two-state lasing in InAs/InP QD lasers under CW injection. This finding holds significant implications for QD device design, offering potential advancements in performance and functionality.

4. Summary

In this article, we have offered a comprehensive overview of investigation methods for QDs. Firstly, a novel two-step, strain-based, eight-band k·p method is introduced for more accurate electronic structure calculations of QDs. Subsequently, the optical properties of QDs from a quantum mechanical perspective are investigated, describing absorption or gain as the loss or gain of photons through upward or downward transitions of electrons between two allowed states of an atom. The carrier (electron and hole) transitions between these atomic states are explored, providing insight into gain and recombination in QDs through factors such as dipole matrix element, momentum matrix element, and oscillation strength. Furthermore, four compelling research subjects utilizing the aforementioned theoretical methods are explored. Firstly, the electronic and optical properties of InAs/InP QDs with ultrathin GaP sublayers are analyzed, which has revealed potential benefits for optimizing high-speed QD devices, as the ultrathin GaP sublayer offers control over non-radiative transitions between the GS and the first ES. Additionally, it may facilitate emission blue-shifting and gain bandwidth control in QD lasers by altering the carrier confinement potential, leading to potential applications in chirped QD stacks and ultra-broadband PL emission. Next, analysis of the electronic and optical properties of InAs/InP QDs with a small spacer thickness have shown that such structures can aid in polarization control of QD devices, as the high TM polarization of closely stacked QDs with small SLNs is attributed to strong transitions between first excited states rather than significant HH and LH mixing. Additionally, this structure may facilitate two-state lasing in InAs/InP QD lasers under CW injection, as the phonon-bottleneck effect could be mitigated by the higher degeneracy ratio between ES and GS, resulting from thin-spacer-induced electronic coupling.