It is a great honor to write this short review in memorial of Dr. Kun Huang, one of the founding persons of Chinese semiconductor science, and a great pioneer in studying electron–phonon coupling in semiconductors. I didn't have the opportunity to meet with Huang in person, but have the fortune to learn from and collaborate with many of his students and student's students. The favorite topic of Huang is electron–phonon coupling and its consequences in semiconductors physics. In this short review, I will discuss a few recent developments in this area, mostly based on my own works, some of them in collaborations with colleagues in the Institute of Semiconductors in Chinese Academic of Science for which Huang has served as its first director. This by no mean is a complete review, but rather some personal views about this topic.

One topic is the nonradiative decay of free carriers caused by defects in semiconductor, most notably the Shockley-Read-Hall (SRH) process. This is an important process in the operation of semiconductor devices, especially for deep defect states. The electrons (or holes) hop from the band edge to the defect state through multiple phonon processes, and it can subsequently lead to carrier annihilation. The basic formalisms of such multiphonon processes have been worked out in 50's by S. K. Pekar^[1] and K. Huang, A. Rhys^[2]. These derivations are based on Franck-Condon approximation of the electron–phonon wave functions, where the electron–nuclear wave function is approximated as , here i and n are electron and phonon eigen state index, and r and R are electron and nuclear degree of freedom. While the electron wave function satisfies the electron Schrodinger's equation at a given nuclear coordinate R:

$ H(r,R) \psi_i(r,R) = \epsilon_i(R) \psi_i(r,R) . $ (1)

The nuclear wave function satisfies its own Schrodinger's equation using as its potential energy:

$ \bigg[ \sum\limits_R {-{1\over 2 M_{ R}} \nabla_{ R}^2 + \epsilon_i(R)} \bigg] \phi_{i,n}(R) = E_{i,n} \phi_{i,n}(R). $ (2)

Here is the total eigen energy of the electron–nuclear wave function , is the nuclear mass. However, the is not a true eigen state of the electron-nuclear total Hamiltonian due to a cross derivative term: as well as a second derivative term: . Usually the cross derivative term is much larger than the second derivative term, and it causes coupling between the initial electron state and the final electron state . This naturally leads to the adiabatic state coupling formalism^{[1, 2]}.

Although elegant in its derivation, very soon it was realized that the adiabatic state coupling formalism gave too small transition rates compared with experiments. To solve this problem, Kovaskiy, Sipdvskiy^{[3, 4]} and Passler^{[5, 6]} proposed to use a static formalism, where both and are replaced by their counter parts and where is chosen as the equilibrium position at one electronic state (e.g., will be either or ). Thus although the state and (the states before and after the transition) are the electron eigen states when , when deviates from this , they are no longer eigen states due to the dependence of the electron Hamiltonian on . Thus, the phonon displacement away from the equilibrium position causes the electronic state coupling.

For a long time, there were debates for the validity and the meaning of using static coupling formalism^[3–8]. Huang gave a very interesting derivation in early 1980^{[8, 9]}. He shown that, in a perturbation representation of the electronic state using the static state as the basis set, if one replaces the static eigen energies , in the perturbation theory denominator by the R dependent eigen energies: , then one can arrive at the static coupling formalism starting from the adiabatic coupling formalism. This in a sense unifies these two formalisms within the framework of perturbation theory^[8]. At the end, the static coupling formalism will be the preferred formalism.

One can view this problem from other angles. The Franck-Condon approximation is good only at R positions where there is no near degenerated i andjelectron eigen states. Unfortunately, the degeneracy (or say the energy crossing of these two states) is exactly what happens at the coordinate where the transition happens (see Fig. 1). So, the Franck-Condon approximation is invalid at the transition nuclear coordination . This is the same situation as for the conical point in quantum chemistry treatment of molecular systems. This can also been seen from Huang's perturbation treatment. If is used, then at the crossing point , the perturbation theory is diverging. Actually, if one follows exactly the adiabatic state definition, at the crossing point , the physical identity of the i and j state will change, one must be very carefully in identifying which state is which when describing the transition. For example, there will be a large and sharp spike in at due to the identity change (if i is used to indicate one continuous adiabatic state, as indicated by one solid line in Fig. 1). This makes the Franck-Condon electron and phonon wave function separation treatment a bad approximation at this , and perhaps more importantly, it makes the commonly used approximation in the adiabatic coupling treatment: replacing by , totally invalid at . All these justify the use of static formalism.

Fig. 1 also connects the nonradiative state transition to the Landau-Zener transition^{[10, 11]} and Marcus theory^[12]. The Landau-Zener transition is mostly a 1D model. It has been used by Lang et al.^[13] to treat the nonradiative transition. Similar approach has also been used by Alkauskas et al.^[14] more recently in an ab initio treatment. However, the multiphonon process is not necessarily a 1D process, in the sense that the transition cross point (the energy conservation modes) and the coupling (the promoting modes) are not caused by the same phonon degree of freedoms^[15]. More specifically, the energy conservation (accepting) modes (or say the phonon modes which lead to the transition point ) come from the displacement . This displacement comes from the diagonal elements of the electron–phonon coupling: (together with the phonon frequencies). On the other hand, the promoting modes comes from the cross terms of the electron–phonon coupling: . The and the cross coupling can both be viewed as multidimension vectors, and the directions of these two vectors can be rather different. In the high dimensional space of R, their dot product could be close to zero^[16]. For example, in the Zn_Ga+ V_N defect complex in GaN, the energy accepting modes are consisted with acoustic phonon modes, while the promoting modes come from optical phonons^[15].

Marcus theory is also often used to describe electron transitions from one state to another, especially for electronic states located at different positions^[17]. While the energy barrier of the Marcus theory also comes from atomic displacements (hence the diagonal elements of electron–phonon coupling), the electronic coupling can come from different sources, e.g., an applied electric field, or more intrinsic coupling. In the case the coupling is also caused by phonon modes (e.g., the pure multiphonon process as described above), Huang has derived a very interesting high temperature classical approximation formula^[9], which essentially gives the coupling constant of Marcus theory based on the electron–phonon coupling constants and phonon mode frequencies. Even when the main coupling is not caused by electron–phonon coupling, it might be possible that the coupling constant provided by Huang's formula is related to the environmental effects during the Landau-Zener transition^{[18, 19]}. This will be an interesting research topic for the future.

Modern theoretical defect studies are mostly depending on ab initio calculations, especially using density functional theory (DFT)^[20–22]. In recent years, the use the hybrid functional^[23] has made the defect level calculations more accurate, and can yield energy level results consistent with experiments^[24]. One current frontier is to study the dynamic properties of such defects, including the nonradiative decay rates. As discussed above, the basic formalism (e.g., the static coupling formalism) already existed. Thus a contemporary challenge is to use ab initio method to calculate the electron–phonon coupling, and compared the results with experiments. The straight forward calculation of requires the calculation of the change of (the total selfconsistent potential in DFT method) due to the displacement of one atomic coordinates. Thus, if there are N atoms in a supercell (as in a defect calculation), there need to be 3N self-consistent field (SCF) calculations. As N can be in the order of 1–2 hundreds, this makes the calculation rather expensive. We have developed an approach which can yield all the electron–phonon coupling constants for a given , pair in one SCF calculation^[15]. This makes the multiphonon process calculation more practical.

More specifically, to calculate , during the SCF calculation, one can add one extra term in the total electronic charge density: , here is the occupation number of the Kohn-Sham eigen state . The first term is the usual formula to calculate the charge density, while the second term is the additional charge density with being a small parameter. Besides the charge density, all other formulas are kept the same as in a normal Kohn-Sham calculation. Furthermore, the atomic force is calculated using Hellmann-Feynman formula as: ( is the nuclear ionic potential). Then we can shows that^[15]:

$ C(i,j,R) = {{\rm d} F_{R,\alpha} \over {\rm d} \alpha } . $ (3)

The derivative regarding to can be done numerically, by using (the original SCF calculation results) and a small like 0.1. In doing so, we can get all the coupling constant in one extra SCF calculation. The same variational approach can be applied to hybrid functional calculations. In that case, during the SCF calculation, besides the above additional term in the charge density, one also needs to add one extra term in the Fock exchange integral as:

$ {\hat P} \psi_k = \alpha \int [\psi_i(r) \psi_j^*(r')+\psi_j(r) \psi_i^*(r')] f(|r-r'|) \psi_k(r'){\rm d}^3r', $ (4)

here is the truncated Coulomb interaction kernel used in the hybrid functional. Besides the electron–phonon coupling constant, one also needs to calculate the phonon spectrum of the defect system in order to use the analytical formulas derived by Huang et al. One way to calculate the phonon spectrum of the supercell is to calculate the dynamic matrix . Once again, this requires the numerical displacement of all N atoms within the supercell. Thus the benefit of the above variational calculation of the electron–phonon coupling constant will be lost if such dynamic matrix needs to be calculated directly. Fortunately, we found that^[15] if both and are away from the point defect, beyond a cutoff radius , then can be approximated by the counter part from the perfect crystal. For crystal, due to the translational symmetry, one only needs to displace the atoms within one primary cell, instead of all the atoms within a supercell. As a result, we only need to calculate for within . Typically 10–50 numerical displacements are needed to carry out the SCF calculations to yield all needed . It has been shown this procedure yield very accurate phonon spectrum for the supercell system containing one point defect^[15].

After we obtain both the electron–phonon coupling constant and phonon spectrum, the multi-phonon nonradiative process can be calculated directly. There are different ways to derive the analytical formula, all based on the use of Fermi Golden rule between the initial and final states with Franck-Condon separation of the electron wave function and phonon wave function. The coupling Hamiltonian in the static coupling approximation is a perturbation term proportional to (note, the R in the summation and in is used as an index, while in it is a vector). is the starting point for the static coupling calculation. There is an ambiguity for what to choose for . Usually one either choose it as or . However, for most cases, since the promoting mode direction of is almost perpendicular to the accepting mode direction of , the choice of or does not really matter (as well be shown later in Fig. 2). What left is to evaluate under a thermodynamic assembly. The harmonic approximation is taken to describe phonon modes and . Very often it is assume they are the same phonon modes (and frequences) but with a zero point displacement of , although analytical equations can also be obtained if different harmonic phonon modes at and are used^[25]. To evaluate the under a thermodynamic assembly and a delta function for energy conservation, one first converts the delta function for energy conservation into an integral using Dirac distribution function:

$ \delta(\omega) = {1\over 2{\text{π}}} \int_{-\infty}^{\infty} {\rm e}^{i \omega t} {\rm d}t . $ (5)

This always leads to an integration for the final formula. After subsequent derivations based on matrix manipulations^[26–28], one can obtain a formula for the i to j nonradiative transition rate as:

$ W_{ij} = {2{\text{π}}}\sum\limits_{k_1,k_2} C_{i,j}^{k_1} C_{i,j}^{k_2} A_{ij}^{k_1,k_2} , $ (6)

here is just the electron–phonon coupling constant , however, converted into the phonon mode coordinate . The matrix can be evaluated as:

$ A_{ij}^{k_1,k_2} = {1\over 2{\text{π}} Z} \int_{-\infty}^{\infty} \chi_{ij}^{k_1,k_2}(t,T) {\rm e}^{-i (E_i-E_j) t} {\rm d}t, $ (7)

here is the phonon partition function and , the and are the defect state energies of the electronic state i and j at their equilibrium atomic positions and respectively (Fig. 1). The expression for the matrix is a bit complicated^[16]. It is expressed by several other matrices with phonon frequencies and the atomic displacement as their variables.

In our previous work^[16], we have adopted the above formalism by Borrelli et al.^[28]. However, we later found that Huang gave a different derivation earlier in 1981^[9]. He explicitly integrated out the harmonic phonon wave functions using their Gaussian representations. He arrived at a more concise formalism:

$ \begin{split} W_{ij} =\; & {2{\text{π}}}\!\!\int_{-\infty}^{\infty}\!\!\bigg\{ \Big[\!\sum\limits_k \! C_{i,j}^k \Delta Q_{ij}^k \big( {\rm{cos}}(\omega_k t) + i{\rm{coth}}(\beta \omega_k/2) {\rm{sin}}(\omega_k t)\big)\!\Big]^2\\ & {+\, {1\over2} \sum\limits_k |C_{i,j}^k|^2 {1\over \omega_k} \big( {\rm{coth}}(\beta \omega_k/2) {\rm{cos}}(\omega_k t)+i {\rm{sin}}(\omega_k t)\big) \bigg\} }\\ & \times {1\over 2{\text{π}}} {\rm{exp}}\big[-i t (E_j-E_i) - \sum\limits_s {\omega_s\over2} |\text{Δ} Q_{ij}^s|^2\big( {\rm{coth}}(\beta \omega_s/2)\\ & \times(1-{\rm{cos}}(\omega_s t)-i {\rm{sin}}(\omega_s t)\big) \big] {\rm d}t. \end{split}$ (8)

Here and are the atomic displacements between state and converted into phonon mode coordinates k and s respectively. We have numerically tested this concise equation versus the more complicated equation derived by Borrelliet al.^[28], they give the exactly same results.

The concise formalism allows Huang to apply the steepest decent approximation to get an closed analytical formula (without the dt integration) for high temperature approximation. To do that, one finds the maximum of the exponent as a function of t, and expands it with a second order approximation. The resulting Gaussian exponential can be integrated over t, yielding in a closed analytical formula. Such closed analytical formulas exist for adiabatic approximation of multiphonon transition^[7], as well as for quantum mechanical treatment of the Marcus theory^[29] where the transition coupling between the two electron states and is assumed to be an constant, independent of the phonon degree of freedom. For the static coupling approximation, we only found such closed formula in Huang's work^[9]. We can thus call it the Huang's formula, which is:

$ W_{ij} = \left({{\text{π}} kT\over S \bar \omega}\right)^{1/2} \left( \sum\limits_k {1\over \omega_k^2} |C_{i,j}^k|^2\right) {\rm{exp}}\left(-{(E_i-E_j -S \bar \omega)^2\over 4 kT S \bar \omega}\right) . $ (9)

Here S is the Huang-Rhy's factor^[2], and . The is nothing but the reorganization energy used in classical formula like the Marcus theory. This reorganization energy equals as shown in Fig. 1. The above formula can be directly compared with the Marcus theory which is:

$ W_{ij} = \left({{\text{π}} kT\over \lambda}\right)^{1/2} {1\over kT} |V_{\rm c}|^2 {\rm{exp}}\left(-{(E_i-E_j-\lambda)^2\over 4 \lambda kT}\right) . $ (10)

Thus, Huang's formula provides an expression for the coupling constant in the Marcus theory as: this is valid for the phonon induced coupling between the electron states and . As discussed before, one can think about other causes for the coupling in a more general case. Nevertheless, the Huang's formula can be used to estimate the phonon contribution to such coupling.

We have tested Huang's formula of Eq. (9) versus direct dt integration of Eq. (8). As shown in Fig. 3, the high temperature formalism is valid when the temperature T is higher than 300 K for the case of Zn_Ga–V_N defect in GaN and a hole transition from the valence band edge to the defect state.

We can now numerically calculate the multiphonon transition rate using ab initio methods, either using the direct integration as in Eq. (8), or using the high temperature Huang's formula of Eq. (9). Such calculation procedure has been implemented within the PWmat code package^[30], in an automatic way to calculate the electron-phonon coupling, defect phonon spectrum, and the static coupling formula for the nonradiative transition rate. The calculated results are usually in good agreement with the experiments. Fig. 2 shows a comparison of different calculation methods, both with the explicit multiphonon mode calculation, and 1D models along the direction. As we can see that, in this case, the 1D model all give very different results compared with the explicit multiphonon static coupling formula. The static coupling results are in good agreement with the experiment. Furthermore, the 1D model sensitively depend on which equilibrium point is chosen for : or . On the other hand, the static coupling result is insensitive to such choices.

There are still challenges for accurate prediction of the nonradiative decay rate and the related studies. The first is an experimental one, as there is a scarcity of the experimentally measured nonradiative decay rate for different defects. The commonly used method: deep level transient spectroscopy (DLTS) can miss some deep defect with low concentrations, and results can also be influenced by factors like Coulomb repulsion between the charged defect and the band edge free carriers. Although there are other alternative techniques, both optical measurements and non optical measurements^[32], due to possible multichannel competition, the interpretation of the results can still be challenging. In terms of calculation, the anharmonicity of the phonon oscillation is one major uncertainty. In ab initio calculation, it is often found that the directly calculated the atomic relaxation energy after the electron transfers from i to j state (e.g. from band edge state to the defect state), can be larger than the relaxation energy as calculated by . If the relaxation and the transition coupling happens in the same direction, then an 1D approximation can be used, and direct numerical calculation can be used to study the transition process as a Landau-Zener transition^[33]. However, if that is not the case, such 1D approximation cannot be used. In such cases, one approximation is to rescale all the phonon frequency so that will yield the same results as the numerically calculated reorganization energy . This also allows us to evaluate the Eq. (8) in the low temperature situation. Another challenge is that, the phonon modes at and can be different. For Huang's formula, while one can use the average formula for the effective coupling constant, what more critical is the exponential term, which is , here is the barrier between i and j, as shown in Fig. 1. One approximation is to re-evaluate this , the lowest valley crossing point for the two multi-dimension parabolas betweeni and j, and using that to replace the exponential term in Huang's high temperature formula. As a matter of fact such correction can be even applied to the low temperature integration formula. We found that, after such correction, the calculated transition rate can increase by almost an order of amplitude in some cases, bringing the result further closer to the experimental values.

We also like to mention that sometime the SRH electron-hole recombination can happen through a multi-step process, with several intermediate transitions corresponding to different occupations and charge states for the defect level^[34]. Although analytical multi-phonon formula can still be used in such cases, one alternative approach could be to do direct real-time time dependent density functional theory (rt-TDDFT) simulations^[35]. In the rt-TDDFT calculation, the nuclear movement follow the Ehrenfest dynamics. It is a classical description for the nuclear movement, thus it is likely adequate in high temperature limit. Since the electron movement follows the time dependent Schrodinger' equation, it can be used to describe the Landou-Zener transition. Compared to the analytical formalism, one advantage of direct rt-TDDFT is its ability to describe anharmonic nuclear movement and strong electron-phonon coupling. On the other hand, due to the classical description of the nuclear movement, it lacks the detail balance between the i to j transition and j to i transition. Recently, we have added such detailed balance within rt-TDDFT, as a result, it can be used to describe multiphonon nonradiative decay. For example, we have used such direct simulation to study the molecule dissociation caused by electron ionization. This is particularly useful to describe very strong electron-phonon coupling, e.g., inside a small molecule, or for a carrier on a localized defect state and going through multiple stages in such defect state. Future investigation of such problems will be interesting. For larger systems, one can ignore the back reaction from the electron movement to the nuclear movement, thus can use the ground state ab initio molecular dynamics trajectory to describe a time dependent Hamiltonian, then use such Hamiltonian to study the carrier dynamics. This approach is called nonadiabatic molecular dynamics, (NAMD), which allows the simulations of much larger systems (e.g., a few hundred atoms) for much longer times^[36] (e.g., 10 ps). All these provide alternative approaches to study the electron-phonon coupling and its related carrier dynamics behavior. Compared with the analytical formula, one limitation is the simulation time. If the decay lifetime is much longer than 10 ps, direct simulation will become rather difficult. So they are only good for strongly electron-phonon coupled systems (for rt-TDDFT simulations) or for problems involve many electronic state transitions (for NAMD simulations).

In summary, in this short review, some current developments in calculating nonradiative decay rates of defects in semiconductors are presented, as well as a brief review for the derivation of Huang's formula. In particular, one procedure was developed to calculate all the electron-phonon coupling constant within one self-consistent field calculation. Another approximated way is introduced to calculate all the phonon modes within a supercell containing a point defect. Such developments, together with Huang's formalism, allow us to calculate the nonradiative decay rate at ab initio level, and the results agree well with the experiment. We have also discussed some of the remaining challenges and possible approaches to overcome them. These include the anharmonic phonon effects and different phonon modes at electronic stateiandj. Finally, we discussed modern direct simulation methods, either rt-TDDFT, or NAMD, which can be used to study problems with strong electron-phonon coupling and strong anharmonicities of the phonon modes, or to study carrier dynamics involving many electronic states.

The author likes to thank many collaborators from the Institute of Semiconductors, including Prof. Jianbai Xia, Prof. ShuShen Li, Prof. Jingbo Li, Prof. Jun-Wei Luo, Prof. Xiangwei Jiang, for many helpful discussions. Many of the works discussed in this paper were done in collaboration with Dr. Lin Shi. The author also likes to thank Dr. Danylo Zherebetskyy for helping find some of the literatures. This work is supported by the Director, Office of Science (SC), Basic Energy Science (BES)/Materials Science and Engineering Division (MSED) of the U.S. Department of Energy (DOE) under the Contract No. DE-AC02-05CH11231 through the Theory of Material project.