Quantum Theory of Polarized Superlattice Optical Response: Faithful Reproduction of Nakamura’s Blue Laser Spectra

1. Introduction

The discovery in the 1990s of blue laser diodes based on ${\mathrm{In}}_x$${\mathrm{Ga}}_{1-x}$N superlattices in the active region [1,2] triggered rapid progress in solid-state lighting, optoelectronics, and related technologies. Since then, intensive efforts have been devoted to improve epitaxial growth, structural quality, and device performance of GaN-based heterostructures and superlattices (SLs)[3,4,5,6,7,8,9,10,11]. Subsequent research concentrated on Ga-rich InGaN and AlGaN alloys, whose bandgaps enable photon emission across the blue–green spectral range [12,13].

Despite extensive work on carrier transport, effective masses, device efficiency, and brightness, the detailed optical response of InGaN/GaN superlattice diodes has received comparatively less attention. High-resolution spectra reveal groups of narrow emission peaks, but their theoretical explanation has long been regarded as overly demanding due to the quantum-mechanical complexity of the problem. In particular, interpreting the spectra first reported in Nakamura’s pioneering experiments has remained a challenge. Several theoretical and experimental studies have addressed partial aspects of the emission process. Examples include the joint density-of-states model for multi-quantum-well heterostructures [14] and works comparing polar and non-polar quantum wells, which showed that structural non-uniformity broadens optical transitions [15]. Conventional approaches based on effective-mass and $k·p$ models, however, rely implicitly on Bloch’s theorem, which is rigorously valid only for infinite crystals. For finite periodic systems such as superlattices, Bloch-based descriptions become inadequate: they replace a quantum-discrete spectrum with a classical continuum approximation, and in doing so eliminate precisely the discrete features observed experimentally.

In the 1980s, several independent works [35] applied the transfer-matrix method to multilayered structures and, without apparent knowledge of earlier studies, effectively rediscovered the compact formulas of Jones [36] and Abelès [37]. These results were used mainly to analyze transport within the one-propagating mode approximation. In parallel, and starting in the 1990s, the Theory of Finite Periodic Systems (TFPS) was further developed to determine not only resonant energies and transmission coefficients but also energy eigenvalues and eigenfunctions. This approach combines exact single-cell quantum solutions with transfer-matrix properties and provides a rigorous framework that overcomes the limitations of Bloch-based models [16,17,18,19,20]. The TFPS restores the quantum nature of finite systems, yielding discrete eigenvalues and eigenfunctions that naturally reproduce the spectral features of superlattices. Applications of TFPS to polarized structures [21] and to resonant features of Nakamura’s blue laser spectra [22,23] demonstrated its explanatory power. More recently, TFPS has been applied to ${\mathrm{In}}_x$${\mathrm{Ga}}_{1-x}$N blue-emitting SLs, successfully reproducing both the narrow emission peaks and their characteristic clustering, as seen in Nakamura’s spectra [24]. In this framework, optical transitions occur between quantized conduction- and valence-band subbands, while surface states—slightly shifted from the bulk subbands—account for the observed groups of emission lines.

Nevertheless, discrepancies between theory and experiment remain. For example, the calculated spectrum in Figure 1 spans about 3 nm in wavelength, whereas the measured spectrum extends over nearly 4 nm. Moreover, theoretical resonance peaks do not capture the gradual decrease in intensity at higher photon energies. These differences suggest that additional physical effects—cladding-layer asymmetry, local Stark fields, variations in effective masses, and recombination dynamics—must be included. Addressing these effects is the main purpose of the present work.

In our earlier calculations we assumed symmetric cladding layers and neglected polarization. Under these conditions, conduction- and valence-band surface states are nearly degenerate, producing essentially a single isolated transition peak. Introducing asymmetry, however, splits these peaks, as demonstrated in Ref. [17]. For the case considered in Figure 1, the cladding layers have Al contents of 0% and 20%, leading to an observable separation of the surface-state resonances.

A distinctive feature of nitride-based heterostructures and SLs is the presence of spontaneous and piezoelectric polarization [25]. Experiments by Kozodoy et al. [26] showed two-dimensional electron and hole accumulation at opposite interfaces of $Mg$-doped ${(A{l}_{x}G{a}_{1-x}N/GaN)}^n$ and ${(I{n}_{x}G{a}_{1-x}N/GaN)}^n$ superlattices. The associated charge distributions bend the conduction- and valence-band edges linearly, giving rise to the quantum-confined Stark effect (QCSE). QCSE modifies the subband structure, eigenfunctions, and transition probabilities, as demonstrated in numerous experimental and theoretical works [27,28,29,30,31,32]. Within TFPS, the Stark effect can be incorporated by transfer matrices W, expressed in terms of Airy functions and their derivatives, as discussed in Ref. [20].

Another factor shaping the emission spectra is recombination dynamics. Injection currents can suppress some transitions, while energy-dependent lifetimes cause peak intensities to decrease with increasing transition energy. This effect, reported for instance by Anikeeva et al. [33], resembles the lifetime scaling of harmonic oscillator states. To capture these features we employ the time distribution of transition probabilities [34], which naturally accounts for the energy dependence of recombination rates.

To calculate the optical response we use the optical susceptibility:

$$\begin{array}{c}\hfill {\chi }_{op}=\sum _{\nu ,{\nu }^{\prime },\mu ,{\mu }^{\prime }}{f}_{eh}{\displaystyle \frac{{\displaystyle {\left|\int dz{\left[{\psi }_{{\mu }^{\prime },{\nu }^{\prime }}^v\left(z\right)\right]}^{*}\frac{\partial }{\partial z}{\psi }_{\mu ,\nu }^c\left(z\right)\right|}^2}}{{(\hslash \omega -E_{\mu ,\nu }^c+E_{{\mu }^{\prime },{\nu }^{\prime }}^v+E_B)}^2+{\Gamma }_{\nu }^2}},\end{array}$$

(1)

where $E_{\mu ,\nu }^c$ and ${\psi }_{\mu ,\nu }^c\left(z\right)$ are the eigenvalues and eigenfunctions in the $\mu$th conduction subband, with $\nu =1,2,\dots ,n+1$ (n being the number of unit cells). $E_{{\mu }^{\prime },{\nu }^{\prime }}^v$ and ${\psi }_{{\mu }^{\prime },{\nu }^{\prime }}^v\left(z\right)$ are the corresponding valence-subband quantities, $E_B$ is the exciton binding energy ($\approx 10$ meV), and ${\Gamma }_{\nu }$ is the level broadening, inversely proportional to the lifetime of the $\nu$th state.

The remainder of the paper is organized as follows. In Section 2 we outline the TFPS framework including the Stark effect. In Section 3 we present subband structures, eigenfunctions, and optical spectra for different electric-field strengths and asymmetry conditions. Finally, Section 4 discusses the influence of recombination lifetimes on the spectral line shapes.

2. Materials and Methods

2.1. Outline of the Transfer Matrix Approach

Solving the Schrödinger equations of biased superlattices is a prerequisite for calculating the optical response in Equation (1). The spontaneous polarization of charges in ${(I{n}_{x}G{a}_{1-x}N/GaN)}^n$ superlattices produces a potential profile at the conduction- and valence-bands edges, as shown in Figure 2. This type of Schrödinger equation can be solved exactly using the transfer matrix method and the Theory of Finite Periodic Systems, briefly outlined in this section. More detailed analysis are provided in Refs. [16,18,20].

For clarity in labelling wells, barriers, and unit cells, we introduce two indices: the index $j=0,1,2,...,n$ labels the coordinates $z_0=0$, $z_1$, $z_2$, ..., $z_n$, that define the unit-cell boundaries, and the index $l=0,1,...,2n+2$ labels the layers. The layers $l=0$ and $l=2n+2$ are the left and right cladding layers, odd l correspond to wells and even l to barriers.

We define a unit cell as one barrier flanked by two half-wells. Thus, the j-th unit cell contains half of layer $2j-1$, the full layer $2j$, and half of layer $2j+1$. The unit-cell length is $l_c=a+b$. The electric fields in the barriers and wells have opposite orientations and their magnitudes are such that $|F_b|b=|F_w|a$. In the following, we will assume that $F_w$ is positive. When the polarization is in the opposite direction, one must change, accordingly, the signs.

The Schrödinger equations in the wells (l odd) and barriers (l even) are

$$\left(-{\displaystyle \frac{{\hslash }^2}{2m_w^{*}}}{\displaystyle \frac{d^2}{d{z}^2}}+V_w\left(z\right)\right){\psi }_l\left(z\right)=E{\psi }_l\left(z\right),l=odd,$$

(2)

with $V_w\left(z\right)=c_w+F_{w}z$, and

$$\left(-{\displaystyle \frac{{\hslash }^2}{2m_b^{*}}}{\displaystyle \frac{d^2}{d{z}^2}}+V_b\left(z\right)+V_o\right){\psi }_l\left(z\right)=E{\psi }_l\left(z\right),l=even.$$

(3)

with $V_b\left(z\right)=c_b-F_{b}z$. In these equations, $m_w^{*}$, $c_w$, $m_b^{*}$ and $c_b$ are the effective masses and ordinates at the origin, in wells and barriers, respectively, $V_o$ is the barrier height.

With the changes of variables

$$\begin{array}{c}\hfill \xi ={\eta }_w\left(z-{\displaystyle \frac{E-c_w}{F_w}}\right)\mathrm{and}\zeta ={\eta }_b\left(z-{\displaystyle \frac{V_o-E+c_b}{F_b}}\right)\end{array}$$

(4)

where

$$\begin{array}{c}\hfill {\eta }_w={\left({\displaystyle \frac{2m_w^{*}{F}_w}{{\hslash }^2}}\right)}^{1/3}\mathrm{for}{\eta }_b=-{\left({\displaystyle \frac{2m_b^{*}{F}_b}{{\hslash }^2}}\right)}^{1/3},\end{array}$$

(5)

Equations (2) and (3) reduce to Airy equations, whose solutions are

$${\psi }_l(z,E)=a_{l}Ai\left(\xi \right)+b_{l}Bi\left(\xi \right),\mathrm{for}l=\mathrm{odd},$$

(6)

and

$${\psi }_l(z,E)=a_{l}Ai\left(\zeta \right)+b_{l}Bi\left(\zeta \right),\mathrm{for}l=\mathrm{even}.$$

(7)

These Airy-function solutions form the basis for defining state vectors, transfer matrices and ultimately the energy eigenvalues, the eigenfunctions and optical response. Each layer solution contains two coefficients, $a_l$ and $b_l$, which in principle must be determined by enforcing continuity, finiteness and physical boundary conditions. For the 22-layer system considered here, this would involve solving for 44 coefficients directly. Here lies the first major advantage of the transfer matrix method: by constructing a transfer matrix that relates the wavefunction (and its derivative) between two points, e.g. $z_i$ and $z_f$, one automatically ensures continuity throughout the interval without explicitly solving for all coefficients.

Two types of transfer matrices are commonly used: the transfer matrix M, which connects wavefunction components (two for one-dimensional systems), and the transfer matrix W, which connects wavefunctions and their derivatives. The transfer matrices M and W are related by a similarity transformation, as explained in Ref. [20]. For our problem, the W- matrices are the most convenient. Given the solutions ${\psi }_l(z,E)$ in wells and barriers, we define the vectors

$$\begin{array}{c}\hfill f_w\left(z\right)=\left(\begin{array}{c}{\psi }_l(z,E)\\ {{\psi }^{\prime }}_l(z,E)\end{array}\right)=\left(\begin{array}{cc}Ai\left(\xi \right)& Bi\left(\xi \right)\\ {\eta }_{w}A{i}^{\prime }\left(\xi \right)& {\eta }_{w}B{i}^{\prime }\left(\xi \right)\end{array}\right)\left(\begin{array}{c}{a}_l\\ b_l\end{array}\right)\end{array}$$

(8)

for wells, and for barriers,

$$f_b\left(z\right)=\left(\begin{array}{c}{\psi }_l(z,E)\\ {{\psi }^{\prime }}_l(z,E)\end{array}\right)=\left(\begin{array}{cc}Ai\left(\zeta \right)& Bi\left(\zeta \right)\\ {\eta }_{b}A{i}^{\prime }\left(\zeta \right)& {\eta }_{b}B{i}^{\prime }\left(\zeta \right)\end{array}\right)\left(\begin{array}{c}{a}_l\\ b_l\end{array}\right).$$

(9)

If $z_i$ and $z_f$ lie in the same well or barrier, there is no discontinuity point in between. In that case, the coefficients remain the same, and the vectors $f\left(z_i\right)$ and $f\left(z_f\right)$ can be written as ( $s=w,b$)

$$f_s\left(z_i\right)=\left(\begin{array}{cc}Ai\left({\varsigma }_i\right)& Bi\left({\varsigma }_i\right)\\ {\eta }_{s}A{i}^{\prime }\left({\varsigma }_i\right)& {\eta }_{s}B{i}^{\prime }\left({\varsigma }_i\right)\end{array}\right)\left(\begin{array}{c}{a}_l\\ b_l\end{array}\right)$$

(10)

where

$$\begin{array}{c}\hfill {\eta }_s={\mathrm{sign}}_s{\left({\displaystyle \frac{2m_s^{*}{F}_s}{{\hslash }^2}}\right)}^{1/3}\mathrm{for}s=w,b,\mathrm{and}{\mathrm{sign}}_{w/b}=+/-,\end{array}$$

(11)

with ${\varsigma }_i=\xi$ when $z_i$ is a well, and ${\varsigma }_i=\zeta$ when it is in a barrier. Similarly

$$f_s\left(z_f\right)=\left(\begin{array}{cc}Ai\left({\varsigma }_f\right)& Bi\left({\varsigma }_f\right)\\ {\eta }_{s}A{i}^{\prime }\left({\varsigma }_f\right)& {\eta }_{s}B{i}^{\prime }\left({\varsigma }_f\right)\end{array}\right)\left(\begin{array}{c}{a}_l\\ b_l\end{array}\right)$$

(12)

The transfer matrix $W_s(z_f,z_i)$ satisfies the relation

$$\begin{array}{ccc}\hfill f_s\left(z_f\right)& =& W_s(z_f,z_i)f_s\left(z_i\right)\hfill \\ \hfill \left(\begin{array}{cc}Ai\left({\varsigma }_f\right)& Bi\left({\varsigma }_f\right)\\ {\eta }_{s}A{i}^{\prime }\left({\varsigma }_f\right)& {\eta }_{s}B{i}^{\prime }\left({\varsigma }_f\right)\end{array}\right)\left(\begin{array}{c}{a}_l\\ b_l\end{array}\right)& =& W_s(z_f,z_i)\left(\begin{array}{cc}Ai\left({\varsigma }_i\right)& Bi\left({\varsigma }_i\right)\\ {\eta }_{s}A{i}^{\prime }\left({\varsigma }_i\right)& {\eta }_{s}B{i}^{\prime }\left({\varsigma }_i\right)\end{array}\right)\left(\begin{array}{c}{a}_l\\ b_l\end{array}\right)\hfill \end{array}$$

(13)

since the coefficients are the same in both sides of this equation, it is easy to show that the transfer matrix $W_s(z_f,z_i)$ is given by

$$W_s(z_f,z_i)=\left(\begin{array}{cc}Ai\left({\varsigma }_f\right)& Bi\left({\varsigma }_f\right)\\ {\eta }_{s}A{i}^{\prime }\left({\varsigma }_f\right)& {\eta }_{s}B{i}^{\prime }\left({\varsigma }_f\right)\end{array}\right){\left(\begin{array}{cc}Ai\left({\varsigma }_i\right)& Bi{\varsigma }_i)\\ {\eta }_{s}A{i}^{\prime }\left({\varsigma }_i\right)& {\eta }_{s}B{i}^{\prime }\left({\varsigma }_i\right)\end{array}\right)}^{-1}$$

(14)

This matrix propagates the solution from point $z_i$ to point $z_f$. Now consider the case where $z_i$ and $z_f$ belong to adjacent layers l and $l+1$, for example, a well and in its neighboring barrier. Let $z_d$ denote the coordinate of the discontinuity point. Strictly speaking, the points inmediately to the left and right of $z_d$ should be written as $z_d^{\pm }=z_d\pm ϵ$. However, this notation quickly becomes cumbersome. For simplicity, when we write $f_l\left(z_d\right)=W_l(z_d,z_i)f_l\left(z_i\right)$, we mean that $z_d$ lies in the same layer l as $z_i$. Similarly, when we write $f_{l+1}\left(z_f\right)=W_{l+1}(z_f,z_d)f_{l+1}\left(z_d\right)$, we understand that $z_d$ lies in layer $l+1$, the same as $z_f$. The continuity of the wave function and its derivative at $z_d$ requires

$$\begin{array}{ccc}\hfill f_{l+1}\left(z_d\right)& =& f_l\left(z_d\right)\hfill \\ & =& W_l(z_d,z_i)f_l\left(z_i\right)\hfill \end{array}$$

(15)

On the other hand, the vector $f_{l+1}\left(z_f\right)$ can be written as

$$\begin{array}{ccc}\hfill f_{l+1}\left(z_f\right)& =& W_{l+1}(z_f,z_d)f_{l+1}\left(z_d\right)\hfill \end{array}$$

(16)

Substituting $f_{l+1}\left(z_d\right)$ from Equation (15), we obtain

$$\begin{array}{ccc}\hfill f_{l+1}\left(z_f\right)& =& W_{l+1,l}(z_f,z_i)f_l\left(z_i\right)\hfill \end{array}$$

(17)

where

$$\begin{array}{c}\hfill W_{l+1,l}(z_f,z_i)=W_{l+1}(z_f,z_d)W_l(z_d,z_i).\end{array}$$

(18)

Here $W_l(z_d,z_i)$ and $W_{l+1}(z_f,z_d)$ are as defined earlier in Equation (14), with the appropriate modifications. Thus, the transfer matrix $W_{l+1,l}(z_f,z_i)$ propagates the solution from one layer to the next, ensuring continuity without explicitly solving for the coefficients.

The second advantage of the transfer matrix method is its multiplicative property: $W(z_3,z_1)=W(z_3,z_2)W(z_2,z_1)$. This property allows the continuity conditions to be imposed in extended systems, such as a multiple quantum wells (MQWs) or superlattices. In practice, one first identifies all discontinuity points of the potential $V\left(z\right)$ between the initial and final positions, and then applies the continuity conditions at the interfaces between wells and barriers.

For a periodic system like that in Figure 2, the variables, $\xi$ and $\zeta$, at the discontinuity points of the j-th unit cell, bounded by $z_{j-1}$ and $z_j$, are denoted as

$${\xi }_{jl}\left(z_{jl}\right)=\xi {|}_{z_{j-1}+a/2},{\zeta }_{jl}\left(z_{jl}\right)=\zeta {|}_{z_{j-1}+a/2},{\zeta }_{jr}\left(z_{jr}\right)=\zeta {|}_{z_j-a/2},{\xi }_{jr}\left(z_{jr}\right)=\xi {|}_{z_j-a/2}.$$

(19)

The indices $jl$ and $jr$ refer to the coordinates at the left and right boundaries of the barrier in the j-th unit cell.

As explained before, the transfer matrices that connect the state vectors at the ends of the wells, satisfy

$$\begin{array}{ccc}\hfill f_w\left(z_{lj}\right)& =& W_w(z_{lj},z_{j-1})f_w\left(z_{j-1}\right),\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill f_w\left(z_j\right)& =& W_w(z_j,z_{rj})f_w\left(z_{rj}\right),\hfill \end{array}$$

(21)

while for the barriers we have

$$\begin{array}{ccc}\hfill f_b\left(z_{rj}\right)& =& W_b(z_{rj},z_{lj})f_b\left(z_{lj}\right).\hfill \end{array}$$

(22)

Thus, the unit-cell transfer matrix is given by

$$\begin{array}{c}\hfill W(z_j,z_{j-1})=W_w(z_j,z_{jr})W_b(z_{jr},z_{jl})W_w(z_{jl},z_{j-1}),\end{array}$$

(23)

with

$$\begin{array}{c}\hfill W_w(z_q,z_p)=\left(\begin{array}{cc}Ai\left({\xi }_p\right)& Bi\left({\xi }_p\right)\\ {\eta }_{w}A{i}^{\prime }\left({\xi }_p\right)& {\eta }_{w}B{i}^{\prime }\left({\xi }_p\right)\end{array}\right){\left(\begin{array}{cc}Ai\left({\xi }_q\right)& Bi{\xi }_q)\\ {\eta }_{w}A{i}^{\prime }\left({\xi }_q\right)& {\eta }_{w}B{i}^{\prime }\left({\xi }_q\right)\end{array}\right)}^{-1},\end{array}$$

(24)

and

$$\begin{array}{c}\hfill W_b(z_q,z_p)=\left(\begin{array}{cc}Ai\left({\zeta }_p\right)& Bi\left({\zeta }_p\right)\\ {\eta }_{b}A{i}^{\prime }\left({\zeta }_p\right)& {\eta }_{b}B{i}^{\prime }\left({\zeta }_p\right)\end{array}\right){\left(\begin{array}{cc}Ai\left({\zeta }_q\right)& Bi{\zeta }_q)\\ {\eta }_{b}A{i}^{\prime }\left({\zeta }_q\right)& {\eta }_{b}B{i}^{\prime }\left({\zeta }_q\right)\end{array}\right)}^{-1}\end{array}$$

(25)

In a more compact notation, we can write

$$\begin{array}{c}\hfill W_w(z_q,z_p)=\left(\begin{array}{cc}{\theta }_w({\xi }_p,{\xi }_q)& {\mu }_w({\xi }_p,{\xi }_q)\\ {\nu }_w({\xi }_p,{\xi }_q)& {\chi }_w({\xi }_p,{\xi }_q)\end{array}\right),\end{array}$$

(26)

and

$$\begin{array}{c}\hfill W_b(z_q,z_p)=\left(\begin{array}{cc}{\theta }_b({\zeta }_p,{\zeta }_q)& {\mu }_b({\zeta }_p,{\zeta }_q)\\ {\nu }_b({\zeta }_p,{\zeta }_q)& {\chi }_b({\zeta }_p,{\zeta }_q)\end{array}\right).\end{array}$$

(27)

Accordingly, the unit-cell transfer matrix takes the form

$$W=W_{w_{r}j}{W}_{bj}{W}_{w_{l}j}=\left(\begin{array}{cc}\theta & \mu \\ \nu & \chi \end{array}\right).$$

(28)

A key property of periodic systems is that all unit cells share the same transfer matrix. Thus, the transfer matrix connecting two points separated by exactly j unit cells, say the points $z_0$ and $z_j$, denoted as $W_j$, is simply

$$W_j=W^j$$

(29)

More explicitly,

$$W_j=W_j(z_j,z_0)=\left(\begin{array}{cc}{\theta }_j& {\mu }_j\\ {\nu }_j& {\chi }_j\end{array}\right)={\left(\begin{array}{cc}\theta & \mu \\ \nu & \chi \end{array}\right)}^j.$$

(30)

At first sight, one might attempt to evaluate $W_j$ analytically by repeated matrix multiplication. However, for large j, this procedure becomes impractical, as the algebra quickly grows in complexity. Numerical multiplication can be used instead, but this approach tends to obscure valuable analytical information. The theory of finite periodic systems offers a more elegant and accurate alternative.

The compact closed-form expressions for the matrix elements of $W_j$ in the single-mode case (one propagating mode to the left and one to the right) were first derived independently by Jones [36], Abelès [37] and later rediscovered in several works [35]. These results, which greatly simplify the analysis of multilayer systems can be written as1,

$$\begin{array}{ccc}\hfill {\theta }_j& =& U_j-\chi U_{j-1}{\mu }_j=\mu U_{j-1}\hfill \\ \hfill {\nu }_j& =& \nu U_{j-1}{\chi }_j=U_j-\theta U_{j-1}.\hfill \end{array}$$

(31)

where $U_j$ is the Chebyshev polynomial of the second class and order j, evaluated at $(\theta +\chi )/2$.

A few remarks are in order. An important property of transfer matrices, established in Ref. [38], is that they form group whose structure depends on the Hamiltonian symmetry class (or universality class). It was shown also that the transfer matrices factorize into the product of two matrices belonging to compact and noncompact subgroups, respectively. Over the last two or three decades, analyses of periodic systems have often conflated the Abéles transfer matrix with Bloch-Floquet theory, which as mentioned before applies strictly to infinite periodic systems. As clarified in Ref. [20], the only transfer matrices consistent with Bloch-Floquet theory are those belonging to the compact subgroup. These matrices faithfully capture phase evolution but lose another essential quantum property: tunneling probability.

2.2. Energy Eigenvalues and Eigenfunctions

Although the transfer matrices W describe how wave functions evolve inside the superlattice, our main objective here is to determine the eigenvalues and eigenfunctions of the confined ${(I{n}_{y}G{a}_{1-y}N/I{n}_{x}G{a}_{1-x}N)}^n$ superlattice. As in standard quantum mechanics, the eigenvalue equations arises when boundary conditions for the confined system are enforced. This requires including the cladding layers on both sides of the superlattice. Given the transfer matrix $W(z_R,z_L)$, that connects the state vectors at $z_L=0$ and $z_R=n{l}_c+a$, still inside the superlattice, one has first to transform the transfer matrix $W(z_R,z_L)$ into the transfer matrix $M_S(z_R,z_L)$, which connects wave vectors defined in the cladding regions [20]. Specifically,

$$\begin{array}{ccc}\hfill {\left(\begin{array}{c}{a}_R{e}^{q_R(z-z_R)}\\ b_R{e}^{-q_R(z-z_R)}\end{array}\right)}_{z=z_R}& =& M_S(z_R,z_L){\left(\begin{array}{c}{a}_L{e}^{q_{L}z}\\ b_L{e}^{-q_{L}z}\end{array}\right)}_{z=z_L}\hfill \\ & =& M_R{W}_r(z_R,z_n)W_n(z_n,z_0)W_l(z_0,z_L)M_L{\left(\begin{array}{c}{a}_L{e}^{q_{L}z}\\ b_L{e}^{-q_{L}z}\end{array}\right)}_{z=z_L}\hfill \end{array}$$

(32)

where $M_R$ and $M_L$ are the transition transfer matrices[20]

$$\begin{array}{c}\hfill M_R={\left(\begin{array}{cc}1& 1\\ q_R& -q_R\end{array}\right)}^{-1}\mathrm{and}{M}_L=\left(\begin{array}{cc}1& 1\\ q_L& -q_L\end{array}\right)\end{array}$$

(33)

with

$$q_L={\left({\displaystyle \frac{2m_L^{*}}{{\hslash }^2}}(V_{wL}-E)\right)}^{1/2}{q}_R={\left({\displaystyle \frac{2m_R^{*}}{{\hslash }^2}}(V_{wR}+F_{w}a-E)\right)}^{1/2}.$$

(34)

The finiteness of the wave functions $\psi (z<z_L)=a_L{e}^{q_{L}z}+b_L{e}^{-q_{L}z}$ and $\psi (z>z_R)=a_R{e}^{q_R(z-n{l}_c+a)}+b_R{e}^{-q_R(z-n{l}_c+a)}$, at the left and right sides of the superlattice, imply that $b_L=a_R=0$. Therefore the relation between the wave vectors at the right and the left are expressed as

$$\left(\begin{array}{c}0\\ b_R{e}^{-q_R(z-n{l}_c+a)}\end{array}\right)=M_R{W}_r(z_R,z_n)W_n(z_n,z_0)W_l(z_0,z_L)M_L\left(\begin{array}{c}{a}_L{e}^{q_L\left(z\right)}\\ 0\end{array}\right)$$

(35)

From this condition we have the energy eigenvalues equation. If the transfer matrices $W_r(z_R,z_n)$, $W_n(z_n,z_0)$ and $W_l(z_0,z_L)$ are represented as

$$W_r=\left(\begin{array}{cc}{\theta }_r& {\mu }_r\\ {\nu }_r& {\chi }_r\end{array}\right),W_n=\left(\begin{array}{cc}{\theta }_n& {\mu }_n\\ {\nu }_n& {\chi }_n\end{array}\right),\mathrm{and}{W}_l=\left(\begin{array}{cc}{\theta }_l& {\mu }_l\\ {\nu }_l& {\chi }_l\end{array}\right)$$

(36)

The energy eigenvalues equation is

$$h_w{U}_n+g_w{U}_{n-1}=0$$

(37)

where

$$\begin{array}{c}\hfill h_w={\theta }_l{\nu }_r+q_r({\theta }_l{\theta }_r+{\mu }_r{\mu }_l+q_l({\theta }_r{\mu }_l+{\mu }_r{\chi }_l))+{\nu }_l{\chi }_r+q_l({\mu }_l{\nu }_r+{\chi }_l{\chi }_r)\end{array}$$

(38)

and

$$\begin{array}{ccc}\hfill g_w& =& q_R\left({\theta }_l{\mu }_r\nu +{\theta }_r\mu {\nu }_l-\theta {\mu }_r{\nu }_l-{\theta }_l{\theta }_r\chi +q_L\left({\mu }_l{\mu }_r\nu -{\theta }_r{\mu }_l\chi +{\theta }_r\mu {\chi }_l-\theta {\mu }_r{\chi }_l\right)\right)\hfill \\ & & +\mu {\nu }_l{\nu }_r-{\theta }_l{\nu }_r\chi +{\theta }_l\nu {\chi }_r-\theta {\nu }_l{\chi }_r+q_L\left(-{\mu }_l{\nu }_r\chi +\mu {\nu }_r{\chi }_l+{\mu }_l\nu {\chi }_r-\theta {\chi }_l{\chi }_r\right)\hfill \end{array}$$

(39)

The solutions of (37) are the discrete energy eigenvalues $E_{\mu ,\nu }$, with two indices; the index $\mu =$1,2,... labels the subbands, in the conduction or the valence bands, and the index $\nu =$ 1, 2, ..., $n+1$, labels the quantized levels within each subband.

The next step is to compute the eigenfunctions ${\psi }_{\mu ,\nu }\left(z\right)$. Before doing so, it is convenient to obtain the general wave function at any point z inside the superlatice, as a function of both z and energy E. This is achieved from the transfer matrix relation

$$\left(\begin{array}{c}\psi (z,E)\\ {\psi }^{\prime }(z,E)\end{array}\right)=W(z,z_L)M_L\left(\begin{array}{c}{a}_L\\ 0\end{array}\right)$$

(40)

Thus, $\psi (z,E)$ at z is fully determined for any E. The eigenfunctions ${\psi }_{\mu ,\nu }\left(z\right)$ follow from evaluating $\psi (z,E)$ at the eigenvalues $E=E_{\mu ,\nu }$.

The relations presented here will be used in the next section to determine the energy eigenvalues and the eigenfunctions in the conduction and valence bands for a specific system. As emphasized in the Introduction, our aim is to apply this formalism to refine the theoretical description of the optical response of blue emitting ${(I{n}_{x}G{a}_{1-x}N/GaN)}^n$ superlattices reported by Nakamura[2]. For this purpose we will introduce the cladding layers asymmetry and the local Stark effect, which directly affect the wave numbers $q_L$ and $q_R$, and the arguments of the Airy functions.

3. Results

The theoretical formalism outlined above provides exact formulas that solve the Schrödinger equation for the periodic potential describing the nitride-based superlattice $GaN{(I{n}_{y}G{a}_{1-y}N/I{n}_{x}G{a}_{1-x}N)}^n$$I{n}_{y}G{a}_{1-y}N/A{l}_{y}G{a}_{1-y}N$ within the effective mass approximation. It has recently been shown [47] that the main assumptions of this approximation remain valid for thin layers.

Before analyzing the effects of the confining potential and local fields on the optical response of the blue-emitting superlattice, we first examine how these parameters affect the fundamental quantities required for the calculation of the optical susceptibility: the energy eigenvalues and the corresponding eigenfunctions.

3.1. Effect of Symmetry and Polarization on Energy Eigenvalues and Eigenfunctions

For concreteness, we focus on the ${(I{n}_{y}G{a}_{1-y}N/I{n}_{x}G{a}_{1-x}N)}^{n}I{n}_{y}G{a}_{1-y}N$ superlattice, with $x=0.05$ and $y=0.2$, and lateral layers $GaN$ and $A{l}_{y^{\prime }}G{a}_{1-y^{\prime }}N$, with $y^{\prime }=0.2$, whose optical response is shown in Figure 1. The heights of the confining potentials in the lateral layers,2$V_R$ and $V_L$, define the asymmetry $a_s=V_R/V_L$.

We selected this system to apply the theory developed in the previous section because its optical spectrum encompasses the full set of interband optical transitions between quantized conduction- and valence-subband levels, including those involving surface states. Moreover, the optical response was measured with the highest available resolution (∼0.016 nm). From earlier quantum-theoretical calculations, where asymmetry and local polarization were neglected, we found that three parameters play a decisive role in shaping the main spectral features: the carrier effective masses, the width and depth of the local potential, and the repulsion of the surface energy levels.

A key spectral feature, closely tied to the eigenvalue spectrum, is the width of the resonance groups: approximately 2 nm for the nine resonances in the long-wavelength (low-energy) region, and about 0.5 nm for the second group of resonances between $\lambda =$419 nm and $\lambda =$420 nm. The resonance spacing is also an important characteristic. At higher energies, isolated peaks correspond precisely to the separation energies of repelled levels, identified as surface states. Their exact positions provide additional constraints on the quantum system, which our formalism reproduces accurately. The local barrier heights (or well depths) and the lateral confining potentials also play an essential role in reproducing the observed spectrum.

We now examine how the local Stark field and the asymmetry affect the eigenvalues and eigenfunctions. Figure 3 and Figure 4 show the influence of polarization and confining-layer symmetry on the eigenvalues and eigenfunctions of the superlattice ${(I{n}_{y}G{a}_{1-y}N/I{n}_{x}G{a}_{1-x}N)}^{n}I{n}_{y}G{a}_{1-y}N$.

Figure 3 illustrates the effect of polarization. While it has only a minor influence on the symmetry of the eigenfunctions, its impact on the energy eigenvalues is significant: the levels shift upward under a positive electric field F, and downward under a negative field. A careful analysis of the spacing between eigenvalues reveals the well-known tendency of the Stark effect to redistribute energy levels toward an almost equidistant spectrum. Overall, the field effect is modest, acting mainly as a tuning parameter when reproducing the optical response.

Figure 4 compares the eigenfunctions and eigenvalues for symmetric and asymmetric confining layers. Here, the main effect is on the surface state: the surface energy level separates from the subband energy levels, and the symmetry of the eigenfunctions is visibly altered near the surface. The loss of symmetry increases with increasing asymmetry $a_s$. The surface states (not shown here) become localized at the surface. In the calculation of transition matrix elements, those involving surface states—either as initial or final states—are the most important, and are responsible for the appearance of groups of resonances, as shown in Ref. [24]. While the surface state is strongly affected, shifting the energy level and modifying the overall spectrum width, the bulk subband energy levels remain essentially unchanged.

There are both differences and similarities between the effects of polarization and symmetry. The local field leaves the overall spectrum width essentially unchanged, while it slightly reduces the level spacing. The main effect of confining-layer asymmetry is to widen the spectrum by shifting the surface energy level. In both cases, however, the bulk subband width defined by the difference $\Delta E_s=E_{\mu ,n-1}-E_{\mu ,1}$ remains nearly constant. In fact, for the first conduction subbands shown in Figure 3 and Figure 4, we find $\Delta E_s=0.0147841$ eV, $\Delta E_s=0.014705$ eV, and $\Delta E_s=0.01486$ eV for the symmetric, polarized, and asymmetric cases, respectively.

The excited electron energies must also include the exciton binding energies. These are small (∼10 meV) compared with the uncertainty in determining the energy gaps, but they are strongly affected by dimensionality and will be accounted for together with the gap energy $E_g$. Our calculations indicate that in superlattices—quasi-two-dimensional systems—the carrier effective masses are likewise affected by dimensionality, yielding values smaller than those commonly reported. For example, if the widely quoted electron effective mass $m_e^{*}=0.2$ is used, the calculated width of the low-energy resonance set is only about half of the observed value, unless unrealistically small barrier heights and widths are assumed. A similar discrepancy arises for the second resonance group when the literature value $m_{hh}^{*}=2$ is used for the heavy-hole effective mass. Better agreement with experiment is obtained using electron effective masses of about $m_e^{*}=0.15$ in the superlattice wells and hole effective masses of about $m_h^{*}=1.2$–1.4.

Figure 5 shows the subband width as a function of the effective mass. The graph illustrates the strong influence of the effective mass on the subband width $\Delta E_s$. The observed spectrum of the ${(I{n}_{y}G{a}_{1-y}N/I{n}_{x}G{a}_{1-x}N)}^{n}I{n}_{y}G{a}_{1-y}N$ superlattice (with $y=0.2$) corresponds to an effective mass of approximately $m_e^{*}\sim 0.135+0.09\times 0.2\sim 0.153$.

3.2. Faithful Reproduction of Nakamura’s Blue Laser Optical Response

Our purpose in this section is to reproduce the optical response of Nakamura’s blue laser diode. Reproducing the optical spectrum requires more than an accurate formalism: the quantum solutions, with appropriate parameter choices, must yield an optical susceptibility that faithfully reproduces the experimental response shown in the upper part of Figure 1. In the next section we will analyze the effect of each parameter separately.

Figure 6 presents several best reproductions of the optical response, obtained by adjusting the confining potential asymmetry $a_s=V_R/V_L$, the polarization field $F_w$, and the electron and heavy-hole effective masses.

We present multiple reproductions because of two reasons: first, the polarization symmetry is slightly broken by confining asymmetry; second, uncertainties and discrepancies in reported parameter values allow multiple parameter sets that describe the experimental results equally well. Figure 6 shows two such sets of parameters (barrier heights, effective masses, and local fields). In the left column, the spectra correspond to conduction-band barrier heights of $0.1995$ eV, while in the right column they correspond to $0.2086$ eV. In (b) and (d) the fields are assumed negative, and in (c) and (e) positive. As can be seen, the predictions are practically indistinguishable, reproducing most features of the experimental spectrum, with slightly different—but consistently smaller than reported—effective masses.

Among the parameter values indicated in Figure 6, those of the effective masses for both electrons and holes are especially noteworthy. As mentioned above, the effective masses required to reproduce the spectrum are significantly smaller than those reported for bulk semiconductors. Literature values exhibit large discrepancies in effective masses, bowing parameters, and related quantities [39,40,41,42,43,44,45,46]. In this context, the quantum-theoretical predictions presented here provide guidance in determining the parameter values that best reproduce the experimental optical response.

In earlier work [19], transition selection rules were derived for symmetric superlattices. With the introduction of asymmetry, one might expect those rules to break down. Our present results, however, show that the same rules remain applicable to an excellent approximation, reducing the number of allowed transitions by about a factor of two. While this reduction is not critical when the number of levels is on the order of dozens, it demonstrates the robustness of the theoretical description.

It is important to emphasize that these results stem from exact quantum-mechanical calculations, based on analytical solutions of the Schrödinger equation, with the effective mass approximation as the only assumption. Therefore, the degree of agreement with experiment directly reflects the adequacy of the chosen parameters for the system under study.

Given the quantum solutions and the optimum parameter sets, we will choose the set corresponding to panel (b) in Figure 6. In the following sections, we analyze how each parameter affects the optical spectrum by varying one parameter while keeping the others fixed.

3.3. Optical Response as a Function of Asymmetry, Polarization, and Effective Mass

Figure 7 and Figure 8 show the effects of asymmetry, polarization, and effective mass. In each case, one parameter is varied while the others are kept fixed at their optimum values. As noted earlier, and visible in the left column of Figure 7, increasing asymmetry shifts one of the surface-state levels upwards, thereby widening the entire spectrum. In the right column, the effect of different polarization fields is shown: less negative fields shift the spectrum upward. The field effect is subtle—larger $|F_w|$ shifts the lower edge of the spectrum to smaller energies, decreases the intra-subband spacing, and slightly increases the asymmetry-induced repulsion of surface states.

As already mentioned, smaller effective masses are needed to reproduce the optical response. The lower panels of Figure 8 show the effect of increasing the effective masses. On the left, the electron mass controls the width of the low-energy resonance group ($\lambda =420$–422 nm). On the right, the hole mass controls the higher-energy resonance group ($\lambda =419.5$–420 nm). In both cases, larger effective masses shrink the group widths. Clearly, adopting bulk effective masses would yield optical responses much narrower than observed. The quantum theory presented here demonstrates that effective masses in quasi-two-dimensional systems are smaller.

3.4. Influence of the Mean Lifetime of Energy Levels

In the previous sections we analyzed the effects of effective mass, symmetry, and polarization on the spectrum width and the relative positions of the resonances. Another important factor, often independent of the injected current, is the mean lifetime of the energy levels within each subband. Not all levels are equally populated at a given time; once populated, the recombination dynamics, together with the data-collection window, determine the observed spectral appearance.

In a recent publication, one of us (PP) derived from time-dependent perturbation theory the time distribution of the transition probability [34] from an energy level $E_{\mu ,\nu }$ with mean lifetime ${\tau }_{\nu }=1/2{\gamma }_{\nu }$:

$$\varrho (t,{\gamma }_{\nu })=2e^{-2t{\gamma }_{\nu }}\left(e^{t{\gamma }_{\nu }}-1\right),$$

(41)

where ${\gamma }_{\nu }={\Gamma }_{\nu }/\hslash$ and ${\Gamma }_{\nu }$ is the level broadening. For simplicity, all level broadenings ${\Gamma }_{\nu }$ in the independent resonances of the optical susceptibility are taken equal; however, the time distribution $\varrho (t,{\gamma }_{\nu })$ is retained because it has a stronger influence on peak heights. Since no direct data exist for subband lifetimes in semiconductor heterostructures, we use as a guide the scaling of harmonic-oscillator lifetimes, given by

$${\tau }_n={\displaystyle \frac{3{ϵ}_{0}mc{\lambda }^2}{2e^2\pi }}{\displaystyle \frac{1}{n}}.$$

(42)

In Figure 9 (left) we plot the time distribution $\varrho (t,{\gamma }_{\nu })$ for the energy levels $E_{1,1}$, $E_{1,2}$ and $E_{1,11}$, with decay times ${\tau }_1=7.933$ns, ${\tau }_2={\tau }_1/2$ and ${\tau }_{11}={\tau }_1/11$. Differences in decay times are crucial for the observed spectral shape: at an observation time $t_{ex}=10$ ns, the transition probability from higher-energy levels is already negligible. In the right-hand column of Figure 9, we plot the time-dependent optical susceptibility

$$\begin{array}{c}\hfill {\chi }_{op,\varrho }\left(t\right)=\sum _{\nu ,{\nu }^{\prime },\mu ,{\mu }^{\prime }}{f}_{eh}{\displaystyle \frac{{\displaystyle |\int dz{\left[{\psi }_{{\mu }^{\prime },{\nu }^{\prime }}^v\left(z\right)\right]}^{*}\frac{\partial }{\partial z}{\psi }_{\mu ,\nu }^c\left(z\right){|}^2}}{{(\hslash \omega -E_{\mu ,\nu }^c+E_{{\mu }^{\prime },{\nu }^{\prime }}^v+E_B)}^2+{\Gamma }_{\nu }^2}}\varrho (t,{\gamma }_{\nu }),\end{array}$$

(43)

for two values of $t=t_{ex}$. In the upper panel at $t_{ex}=$0.5 nm $\sim {\tau }_{11}$, the heights of the high-energy resonances are larger. These resonances practically disappear when $t_{ex}=$8 nm $\sim {\tau }_1$, as can be seen in the optical spectrum of the lower panel.

Reported spectra, however, are not snapshots but the result of data integrated over a time interval $[t_i,t_f]$. Given the instantaneous probability $\varrho (t,{\gamma }_{\nu })$ the probability that a transition from level $E_{\mu ,\nu }$ occurs during the elapsed time $[0,t_{ex}]$ is

$$\pi (t_{ex},{\gamma }_{\nu })=\gamma {\int }_0^{t_{ex}}\varrho (t,{\gamma }_{\nu })dt={\left(1-e^{-t_{ex}{\gamma }_{\nu }}\right)}^2.$$

(44)

This probability approaches unity faster for shorter lifetimes ${\tau }_{\nu }$ (i.e. larger ${\gamma }_{\nu }$. The probability that the transition occurs specifically within the window $[t_i,t_f]$ is then

$$\Delta \pi (t_i,t_f,{\gamma }_{\nu })=\pi (t_f,{\gamma }_{\nu })-\pi (t_i,{\gamma }_{\nu }),$$

(45)

with $\pi (t_{ex},{\gamma }_{\nu })={(1-e^{-t_{ex}{\gamma }_{\nu }})}^2$. In the left-hand panel of Figure 10 we plot $\Delta \pi (t_i,t_f,{\gamma }_{\nu })$ as function of the energy level index $\nu$, for fixed $t_i=$3 ns and different $t_f$. The magnitude decreases steadily with $\nu$, with small differences at low $\nu$. If $t_i=0$ all probabilities tend to one as $t_f$ grows, mirroring Equation (46). In the right-hand panel we show the optical spectra for two time-interval windows. The initial time $t_i=0.3{\tau }_1$ is the same in both panels, but $t_f$ is ${\tau }_1$ in the lower panel an $5{\tau }_1$ in the upper panel. It can be shown that starting the collection at smaller $t_i$ allows larger contributions from high-energy levels, while increasing $t_f$ enhances contributions from long-lived, low-energy states.

To conclude the analysis based on the time distribution of the transition probability, we summarize in Figure 11 the experimental, theoretical without time factor and the optical response evaluated with the optical susceptibility

$$\begin{array}{c}\hfill {\chi }_{op,\pi }(t_i,t_f)=\sum _{\nu ,{\nu }^{\prime },\mu ,{\mu }^{\prime }}{f}_{eh}{\displaystyle \frac{{\displaystyle |\int dz{\left[{\psi }_{{\mu }^{\prime },{\nu }^{\prime }}^v\left(z\right)\right]}^{*}\frac{\partial }{\partial z}{\psi }_{\mu ,\nu }^c\left(z\right){|}^2}}{{(\hslash \omega -E_{\mu ,\nu }^c+E_{{\mu }^{\prime },{\nu }^{\prime }}^v+E_B)}^2+{\Gamma }_{\nu }^2}}\Delta \pi (t_i,t_f,{\gamma }_{\nu })\end{array}$$

(46)

Figure 11 shows the experimental and the theoretical spectra with and without the time factor. In the upper panel the experimental response. In panel (b) the same spectrum as in panel (b) of Figure 6, and in panel (c) the optical spectrum of Equation (45), which takes into account the polarization, symmetry and time effects. The combined effects of asymmetry, polarization, effective mass, and recombination lifetimes provide a faithful description of the blue laser diode optical response.

4. Discussion

The results presented demonstrate that including polarization fields, asymmetry of the confining layers, and recombination lifetimes leads to a faithful reproduction of Nakamura’s blue laser spectra. Importantly, the analysis shows that effective masses in quasi-two-dimensional superlattices are significantly smaller than commonly reported bulk values. This reduction in effective mass, while not the central motivation of this work, emerges naturally from the requirement to reproduce experimental data.

The role of polarization is primarily as a fine-tuning mechanism: it shifts energy levels slightly and redistributes level spacing, consistent with Stark-effect tendencies, but does not strongly alter the overall spectrum width. By contrast, asymmetry in the confining layers has a pronounced effect on surface states, broadening the spectral features and modifying peak positions. Finally, incorporating recombination lifetimes provides a more realistic description of peak intensities and spectral evolution. Together, these results underscore the importance of including both structural asymmetry and dynamical effects when modeling superlattice-based optoelectronic devices.

5. Conclusions

We have developed and applied an exact quantum theory for polarized and asymmetric superlattices, incorporating piezoelectric polarization, confining-layer asymmetry, and recombination lifetimes. This approach, based on analytical solutions of the Schrödinger equation within the Theory of Finite Periodic Systems, allows precise calculation of eigenvalues, eigenfunctions, and optical spectra.

The theory reproduces Nakamura’s blue laser spectra with high fidelity. The analysis demonstrates that realistic reproduction requires adopting smaller effective masses than typically assumed for bulk materials, unless unphysical barrier parameters are used. While polarization primarily plays the role of a fine adjustment, asymmetry significantly affects spectral width and surface-state contributions. Moreover, recombination lifetimes, when treated through their time-dependent distribution, account for the energy dependence of peak intensities and spectral evolution.

These findings not only clarify the role of polarization and asymmetry in superlattice devices but also highlight the need to revisit effective-mass assumptions in quasi-two-dimensional systems. The methodology and insights provided here may guide the design and optimization of next-generation optoelectronic devices.