Analytical Expressions for Numerical Characterization of Semiconductors per Comparison with Luminescence.

Luminescence is one of the most important characterisation tools of semiconductor materials and devices. Recently, a very efficient analytical set of equations has been applied to explain optical properties of dilute semiconductor materials, with an emphasis on the evolution of peak luminescence gain with temperature and its relation to sample quality. This paper summarizes important steps of the derivation of these expressions that have not been presented before and delivers a theoretical framework that can used to apply exactly solvable Hamiltonians for realistic studies of luminescence in various systems.

1. Introduction

Materials development requires characterization techniques. Among them, photoluminescence, or in more general terms, radiation emission due to different excitation mechanisms [1], is a very powerful tool to study semiconductor materials and map specific characteristics of new devices [2,3] for applications from the THz-Mid Infrared (TERA-MIR) to ultraviolet ranges [4,5,6]. A one-to-one correspondence between measured spectra and fundamental materials properties requires a clear theoretical model, ideally easy to understand and to programme, but at the same time with microscopic information for conclusive interpretation and as free as possible from phenomenological parameters.

A recent theoretical effort led to the development of analytical solutions for the interband polarization, which plays the selfenergy role in the Dyson equation for the Photon Green’s functions [7], which have been applied them to study photoluminescence of Coulomb-correlated semiconductor materials. The accuracy of the resulting easily programmable solutions has been demonstrated by consistently explaining the low temperature s-shape of the luminescence peak of dilute semiconductors, such as ternary GaAsSb, InAsN, and quaternary InAs(N,Sb) [7,8,9]. The interplay of homogeneous versus inhomogeneous broadening at low and high temperatures are described, together with the relevance of many body effects, which are in very good agreement with experiments [10,11,12,13]. A similar set of equations was also used to study nonlinearities in GaAs–AlGasAs and GaAsN–AlGaAs superlattices [14,15,16]. The superlattice case is particularly noteworthy for room temperature GHz nonlinear multiplication into the THz range [17,18]. This paper has two objectives: to show hitherto unpublished details of the mathematical steps that lead to the equations used in Refs. [7,8,9,14,15,16], and to draw a bridge between the luminescence and nonlinear absorption calculations in superlattices. It is organized as follows: the main steps involving manipulations of hypergeometric functions that characterize the solution of the Hulthén potential problem are delivered. Next, a direct connection between luminescence and absorption equations is given together with a connection with generalized semiconductor Bloch equations, with potential for the study of polaritons in superlattices within a dielectric approach such as that used to investigate valence band THz polaritons and antipolaritons [6]. A brief conclusion follows.

2. Mathematical Model

In order to make this paper self-contained, some of the steps shown in Ref. [7] are followed to guide the presentation towards the more complete derivation presented here.

2.1. Integro-Differential Equation for the Power Spectrum

Luminescence, or equivalently the optical power density spectrum , is described quantum mechanically by the Poynting vector, which can be directly related to the transverse polarization function , which is the selfenergy in the Dyson equations for the transverse photon Green’s function components in the Keldysh formalism [7,19,20]. All the quantities presented in this paper are considered in frequency space, i.e., evaluated at steady state.

The free photon Green’s function represents the photons propagating without any interaction with the medium. When carriers are injected the transverse polarization function , which is the selfenergy in photon Green’s function Dyson equation, determines how the excited medium modifies the photon propagation. The lesser Keldysh component is proportional to the carriers recombination rate and yields the number of emitted photons per unit area. It thus governs the power emission spectrum, as seen in Equation (1). The imaginary and real parts of are, respectively, proportional to absorption and gain and refractive index changes, since the dielectric function of the medium reads as shown in Ref. [19]. The starting points for the results derived here are the equations in Refs. [7,19,20].

Here, , and denote, respectively, the electron charge, the speed of light, the sample volume and the velocity matrix element, which is the expectation value of the velocity operator, i.e., the momentum operator divided by the electron mass. It stems directly from the fact that current is charge times velocity. The formal definition of the transverse polarization function selfenergy in terms of functional derivatives is , where and are, respectively, expectation values of the induced current and vector potential operators. The full expression involves labels along the Keldysh contour and is tensorial. A complete discussion is beyond the scope of this paper. For details see Refs. [19,20].

The crystal momentum is a consequence of Fourier transforming from real space. Likewise, ω and thus the photon energy stem from a corresponding Fourier transformation from time to the frequency domain. The matrix element satisfies the integro-differential equation [20], where is the screened Hulthén potential [21,22,23,24]. Furthermore,

Electrons or holes are labelled, respectively, by , the renormalized energies , and dephasing are calculated from the real and imaginary parts of the selfenergy in the Dyson equation for the retarded carriers Green’s functions. This paper focuses on quasi-equilibrium luminescence and on three dimensional (bulk semiconductors) with one conduction and one valence band.

Under these conditions, denotes a Fermi function characterized by a chemical potential and the spectral function in Equation (4) for each particle, derived from components of the carriers’ Green’s function in the Keldysh formalism, reads

The next step is to re-write the last term in Equation (4) by means of the identity and to approximate this factor by , where is the total chemical potential, where T is the temperature in Kelvins and is the Boltzmann constant. Different versions of this approximation has been used before in phenomenological approaches for absorption Refs. [21,22,23] and delivered good agreement with experiments (see details and further references in Ref. [23]). Within the Keldysh Green’s functions, context, a detailed derivation of its application is given in Ref. [20]. The fully numerical solutions of the equations that use this version of the approximation have given very good agreement with both single beam and pump-probe luminescence [20,25]. Its usefulness has been further confirmed recently by the good agreement between the analytical solutions shown here and the experimental luminescence of dilute semiconductors [7,8,9,14].

Figure 1 and Figure 2 depict typical ranges of parameters, showing that is an excellent approximation. Note that this theory is applied for photon energies around the semiconductor bandgap and thus in Figure 1 and Figure 2, .

Figure 1

, the term in curly braces in Equation (6) which has been approximated by , evaluated at for bulk GaAs at low temperatures, where the dephasing is typically small . In this case, only a small range of detunings contribute to the integral in Equation (4). Thus the range chosen in the x-axis, from zero to approximately twice the exciton binding energy (), is even larger than necessary. (a) T = 10 K; (b) T = 20 K.

Figure 2

, the term in curly braces in Equation (6) which has been approximated by , evaluated at for bulk GaAs at higher temperatures and high densities, where the dephasing is larger than in the low temperature case. In this case a wider range of detunings contribute to the integral in Equation (4). Thus, the range chosen in the x-axis is even larger than necessary. (a) T = 150 K, (b) T = 300 K.

Low temperature luminescence is typically performed with a small density of injected carriers. Very good agreement of this theory with results from different experimental teams for a variety of materials has been obtained with carrier densities around 1015 carriers/cm3 [7,8,9], further justifying the range of densities in the y-axis. The theory has also been used for high temperatures and high densities to investigate optical nonlinearites [14,15,16], and this range is illustrated in Figure 2.

Relevant dephasing mechanisms such as electron-electron, electron-phonon and electron-impurity scattering can be added to the selfenergy [17,18], and the resulting is frequency and momentum dependent. However, in what follows, it is replaced by averaged values, leading to a simple approximation for consistent with the Ansatz solution, where and . In 3D, the material resonance energy is: , where The bandgap is given by the sum of the fundamental band gap , and a many body renormalisation term where denote, respectively the electron and hole effective masses. The equation for simplifies to:

The total dephasing Γ will determine the luminescence linewidth. Thus, it can be treated as a phenomenological parameter used to interpret data, and at the same time estimate the strength of the scattering and dephasing processes [7,8,9] by comparison of adjusted data with microscopic calculations derived from the relevant selfenergies [17,18]. At this point, the Kubo–Martin–Schwinger (KMS) relation under the form derived in Ref. [20] can be applied to Equation (8), together with the auxiliary variable: , leading to the relation:

Expressing from Equation (8) in terms of the corresponding integro-differential equation becomes where .

Before proceeding, the Hulthén potential [21,22,23,24] should be revised. The usual approximation for a static 3D screened potential is the Yukawa potential, . However, the corresponding Schrödinger equation does not have known analytical solutions. In contrast, the Hulthén potential: , has known analytical solutions that have proven to be very useful for the description of bulk absorption [22]. Recent applications have confirmed its relevance to explain experimental luminescence studies [7,8,9]. Figure 3 shows that, in the range of carrier densities and temperatures of interest, the Yukawa potential can be replaced by the Hulthén potential with negligible differences in numerical values.

Figure 3

Comparison of the Hulthén (dashed lines) and Yukawa (thin solid lines) potentials as a function of distance. The black curves are for T = 10 K and N = 1010 carriers/cm3. The blue curves are for T = 10 K and N = 1018 carriers/cm3 and the red curves are for T = 300 K and N = 1018 carriers/cm3. Both cases depend on the temperature T and carrier density N through the inverse screening length and the inset explains the results, because increases with increasing carrier density and with decreasing temperature. The dot-dashed (cyan) curve is for T = 10 K, while the double-dot-dashed (orange) curve is for T = 300 K.

The Fourier transform of the Hulthén potential has an analytical expression, where is the sample volume, is the Trigamma function [26], is the screening wavenumber and by including in at Equation (11), . Analytical approximations for and are given in Ref. [7]. Note that the bandgap renormalization including Coulomb hole and screened exchange corrections reads

The Fermi functions are evaluated at the peak of the spectral function for each particle, i.e., in Equation (5), . More details are given in Ref. [7]. Equation (13) goes beyond phenomenological term for the bandgap shift [21,22,23], and also, in contrast to those, here we can in principle take into account a reduction in the Coulomb interaction due to phase space filling through the factor . Note however that in the range of carrier densities and temperatures of interest i.e., , as shown in Figure 4.

Figure 4

Inversion factor, for GaAs as function of and ranges of densities consistent with low temperature luminescence, (a) T = 10 K [7,8,9] and (b) high temperature nonlinear absorption T = 300 K [14,15,16,22].

At quasi-equilibrium, used in Refs. [7,8,9,14,15,16,22], the total chemical potential is calculated self-consistently with the many body renormalization of the bandgap and can be written exactly as , where is the total free carrier chemical potential calculated from the bottom of each band. In other words, the inversion factor can be equivalently written as , and it is now clear why Figure 4 has the detuning in the x-axis. Furthermore, the 3D exciton binding energy for GaAs is 4.2 meV and there is no luminescence of absorption for a detuning below 4.2 meV, unless there are deep levels due to impurities and defects, which are not considered here. Thus the approximation, i.e., for the dielectric constant used in the Hulthén potential is clearly excellent in the low power luminescence case. Nonlinear absorption studies are only meaningful away from population inversion leading to optical gain, i.e., the studies are in the range . Thus a decreasing occupation reflects phase space feeling and even for the approach is valid. In order to study the gain regime, the approximation used in the literature is to make at the Hulthén potential and consider the inversion factor only on the right hand side of Equation (8). In other words, in the traditional “plasma theories” for bulk semiconductor absorption and gain, phase space filling ( is not taken into account. Since in the high density case where gain develops, the Hulthén potential decreases due to screening, which described by large in Equation (12), there is still good agreement with experiments. See e.g., Refs. [22,23].

Equation (11) can now be Fourier-transformed

Here is the sample volume, denotes the Dirac delta function. Expanding in the basis of eigenstates of the Hamiltonian: ,

Thus, Equation (15) can be rewritten as

Projection onto state yields

Substitution into Equation (16)

Fourier-transforming back to k-space

From and , a closed expression can be obtained. where a factor 2 for spin has been explicitly written out of the summation over all quantum numbers. Introducing and combining Equations (1), (2), (10) and (21) leads to where reduces to a Dirac delta function for . The velocity matrix element is expressed in terms of the dipole moment matrix element and the fundamental bandgap as . Next, the Schrödinger Equation for the Hulthén potential must be solved, so that can be inserted in Equation (23). The first step is to separate the wavefunction in radial and angular parts. The label thus spans the set ,

The corresponding Schrödinger Equation, which is a generalized Wannier equation [23] can be cast in the form:

The energy eigenvalues depend only on the quantum numbers, and thus we can replace by . Introducing the 3D Rydberg and Bohr Radius , as well , and , leads to

Note that the angular momentum operator has been applied to the wavefunction directly from Equation (25) to Equation (26), i.e., Only solutions that do not vanish at contribute to the emitted power, so is selected. The labels will be dropped at the moment to simplify the notation,

Introducing , and Equation (27) is transformed into

The auxiliary variables, and , lead to the equation which reduces to the Hypergeometric Equation [26,27],

The generalized Wannier Equation, Equation (25), has two types of solutions: bound states for and unbound solutions for . The wavefunctions and eigenvalues are thus different and it makes sense to study each case separately and then add all contributions when a sum over all possible as required from Equation (23).

2.2. Bound States

For bound states, , and thus is real. The function must be normalized, so . This is achieved for either or a negative integer. Set e.g., , in Equation (31), leading to

Thus, this selected solution of the Hypergeometric equation that does not vanish at is where the normalization factor is defined by the normalization condition:

Since the wavefunction is normalized to 1 over angular integration. Setting or equivalently , the other relations become: or . Equation (8.962.1) of Ref. [27] can be used leading to where the Legendre Polynomial and the Gamma function have been used. Introducing the normalization constant thus reads

However, , see Equation (8.960) of Ref. [27] leading to

The solutions for the integrals above are given by Equations (7.391.1) and (7.391.5) of Ref. [27],

Which combined with Equation (6.1.15) of Ref. [26], yields

However, Consequently,

The intensity spectrum requires , where . Note that ,

Note that and . Furthermore, the bound state energies are Thus, in summary

Note that . In other words, the number of bound states that contribute to the power spectrum are determined by the size of the inverse screening length. For zero screening, .

2.3. Continuum States

The unbound solutions that make a continuum have positive eigenvalues, , and thus imaginary . Dropping labels to simplify the development in the next equations yields which can be written for simplicity as

Note that the transformation is being used. The solution that will be later inserted in Equation (23), will be normalized in a sphere of radius and asymptotic solutions, obtained a large radius will be investigated. Next, Equation (15.3.6) from Ref. [26] is used, i.e.,

Furthermore, note that and for all values of . Thus,

, which, combined with (Equation (6.1.2.3) of Ref. [26]), gives in the asymptotic limit

Leading to asymptotic forms of where . The normalization constant is thus given by ,

However, see e.g., Equations (6.1.29) and (6.1.31) of Ref. [26], plus a little algebra deliver the continuum normalization constant

The required value of the wave function at the origin can thus be expressed as

Next, note that , but , leading to

The sum of continuum states becomes an integral, Introducing 2 for spin and changing variables, the continuum contribution becomes which combines with the bound states to deliver the power spectrum

Here, , , and the square of the velocity matrix element is , where the spin orbit shift, the free-carrier bandgap, and renormalized bandgap are given by , and .

Note that this approach does not include cavity effects, which can be introduced in the Photon Green’s functions solution following Ref. [19]. Quasi-periodic structures can also be addressed by a Green’s functions formalism as shown in Ref. [28].

3. Numerical Application

The goal of this section is to illustrate the approach and the many quantities and parameters used making reference to published material, where the equations are used delivering very good agreement with experimental data. Photoluminescence is a very powerful tool to characterize semiconductor materials and map specific characteristics of new devices. Equation (58) is the reference, since it delivers the emission spectrum, which can be directly compared with experimental data. The carriers generated by the photo excitation process modify the spectrum and these modifications are described in Equation (58) approximately by the corrections induced by the (screened) potential, bandgap renormalization, and changes in linewidth governed by the dephasing or scattering . The temperature T can be measured and used as an input parameter. The carrier density N can be estimated by measuring the input power and its spot size when focused on the sample, but in our recent investigations, where this theory has been very successfully compared with experiments [7,8,9], it has been treated as a free parameter, which has been globally adjusted. The other parameters that characterize the material, i.e., the fundamental band gap , the electron and hole effective masses , the static dielectric constant , and the spin-orbit shift can either be found in the literature or robust numerical methods such as simulated annealing can be used to determine these parameters by direct comparison with experiments.

As a reference for the material parameters recently used for dilute nitrides and bismides and the corresponding bandstructure calculations that lead to the material parameters, see: Ref. [7] for GaAs1−Bi; Ref. [9] for InAs1−N and Ref. [8] for more complex quaternary materials, such as InAs1−NSb. The “s-shape” in the luminescence profiles as a function of temperature for these materials have been well explained in Refs. [7,8,9]. However, in the case of completely new materials, expecting to have general characteristics as the ternary or quaternary above, the parameters leading to the bandstructure may be unknown. This theory can be used as a numerical characterization tool as follows.

The corresponding bandstructure can depend on a number of unknown parameters for new compounds, but the approach used in Refs. [7,8,9] can be extended in the following way to extract these parameters by a systematic comparison between theory and experiments. For fixed excitation power, the luminescence can be measured for a number of different temperature points. The dephasing corresponding to different excitation processes can be calculated or taken also as a parameters. Thus, at each Temperature T, there is an ensemble of parameters, such as

Experiments provide a series of data points measured at T = (T1, …, T). The calculated luminescence spectrum will be a function of T and will depend on the ensemble of parameters, denoted The least squares method leads to estimates of the parameter ensemble by minimizing the residual between the theoretical function and the experiments. Therefore, the problem becomes

Trust Region-Reflective (TRR) methods deliver an efficient solution for this numerical problem and Ref. [8] gives further details of their application.

Figure 5 depicts a numerical example to further illustrate choices for the main input parameters. Short period superlattices with strong delocalization of the electron and hole wavefunctions can be described in many cases by anisotropic 3D media, characterized by in-plane and transverse (along the growth direction) effective masses and dielectric constants. The anisotropy parameter is given by the ratio between the in-plane and perpendicular reduced effective masses, , with and which are calculated from the non-interacting superlattice Hamiltonian , , , for i = e, h.

Figure 5

Luminescence for a GaAs–Al0.3Ga0.7As superlattice with repeated barrier and well widths equal to 2 nm. All curves have been calculated with the same broadening meV and a carrier density N = 1014 cm−3 From left to right, the orange, black, blue, violet, green and turquoise curves have been calculated for temperatures given respectively by T = 150 K, 170 K, 200 K, 250 K, 270 K, and 300 K. The corresponding free carrier bandgaps needed for input are = 2.021 eV, 2.015 eV, 2.003 eV, 1.982 eV, 1.974 eV, and 1.960 eV.

These can be calculated from the corresponding free carrier Hamiltonian, and full details of the method, which has led to good agreement with experimental data, can be found in Refs. [15,16,21,29]. Figure 5 shows calculated luminescence using the modified parameters determined by anisotropic medium theory for a short period GaAs–Al0.3Ga0.7As superlattice with repeated barrier and well widths equal to 2 nm. The resulting effective masses are 0.08; 0.12, 0.08, 0.53. These lead the anisotropy parameter 0.67. The resulting exciton binding energy and Bohr radius are given respectively by meV and 11.04 nm.

Except of course for the actual value of the bandgap, which is larger due to quantum confinement, the strong delocalization of electrons and holes in this short period superlattice make the evolution of the luminescence with temperature look qualitatively similar to a three dimensional (bulk) semiconductor, notably the evolution of the line-shape. This is quite similar to the calculations presented in Ref. [9], which are in very good agreement with the experiments discussed in Ref. [13].

4. Conclusions

Photon Green’s function techniques have been used to study different types of luminescence over the years, but the solutions are typically numerically intensive and not accessible for experimentalists or non-specialists. To bridge this gap, a simple analytical solution of the relevant Green’s functions was necessary and the approach described here meets those needs. Notably, the evolution of luminescence with temperature has been successfully compared for GaAsBi [7] and InAs(N,Sb) [8] and InAsN [9] dilute semiconductors. Furthermore, the approach has been used to predict nonlinearities in short period superlattices treaded as anisotropic three-dimensional media [14,15,16]. However, details of the mathematical steps needed to achieve the final formulas used in these publications have not been previously presented and they are given in this review, complemented by numerical results demonstrating the range of validity of the main approximations used in the development of the approach. To complete the picture, an application section shows how to use the method as a numerical characterization machine, and the main parameters needed in typical simulations are illustrated with results for luminescence of short period superlattices. Screening of the Coulomb interaction between electrons and holes is discussed by means of the Hulthén potential and the steps provided are suitable as a guideline to the study of other interacting potentials of interest. They can be followed for further development of a suite of algorithms for efficient and easily programmable numerical characterization tools for a host of new bulk materials or superlattices that can described be as effective 3D media using anisotropic medium approximations.