Transition matrix approach for Monte Carlo simulation of coupled electron/phonon/photon dynamics

A new approach for simulating the dynamics of electrons, phonons, and photons is described. The technique provides a Monte Carlo simulation of particle dynamics without the statistical noise associated with direct Monte Carlo simulation, treats physical phenomena with a wide range of time scales, and has a good computational efficiency. A transition matrix is first precomputed by direct Monte Carlo simulation. Particle populations are then updated at regular time steps by simple matrix multiplication while correcting for nonlinear effects such as carrier–carrier scattering, band filling, hot phonons, etc. The technique is well suited to studies of quantum well laser devices and pump‐probe experiments where direct Monte Carlo simulation is exceedingly difficult.

Pump-probe studies of carriers in bulk and quantum well semiconductors are widely used to study carrier-phonon and carrier-carrier scattering processes both under lasing and nonlasing conditions. 1 For photoexcited semiconductors, Monte Carlo simulations have proven to be invaluable in developing a detailed understanding of the microscopic processes in play ͑e.g., Ref. 2͒. Direct Monte Carlo techniques, however, are not well suited for treating carrier dynamics in quantum well lasers or for analyzing pump-probe experiments of laser gain dynamics ͑e.g. Ref. 3͒ where small perturbations are induced and physical phenomena with both short and long time scales are operative. In this letter, we describe an indirect Monte Carlo simulation procedure based on the use of transition matrices. In comparison to direct Monte Carlo simulation, the technique provides simulations without statistical noise, treats phenomena with both short and long time scales, and provides a much improved computational efficiency.
Consider the space-independent, field-free Boltzmann equation for a nondegenerate semiconductor where the first term on the right-hand side represents outscattering from a state at p, and the second term describes in-scattering from the states pЈ. Equation ͑1͒ could apply to bulk semiconductors, 4 and plasmas, 5 or quantum wells. Here we focus on quantum wells. Because there is symmetry in the plane, the distribution function for a given subband will depend on energy alone. The individual scattering events, however, do depend on angle. As discussed below, it is sometimes advantageous to retain the angle resolution ͑i.e., for phonon scattering events͒, but for electron-electron scattering events it is more efficient to average over angle ͑which is exact because of the symmetry͒. Consider a group of electrons located in bin j at time, t, as shown in the shaded box of Fig. 1. At time tϩ⌬t these electrons will have made transitions to other bins in the same and other subbands. We can transform Eq. ͑1͒ into an equa-tion relating the population of a bin in momentum space, N i (tϩ⌬t), to the population of the bins at time, t, by forward differencing Eq. ͑1͒ in time and by integrating over the area of each bin to obtain where P i j is the probability that an electron in bin j make a transition to bin i in a time ⌬t, and the sum over j runs over all the momentum space bins in all the subbands. It is readily shown that Eq. ͑2͒ applies only when ⌬tӶ. To evaluate P i j in Eq. ͑2͒, we perform a Monte Carlo simulation of 2D, confined electrons ͑e.g., Ref. 2͒. A large number of electrons is injected in bin j at time tϭ0 and tracked by Monte Carlo simulation for a time ⌬t as they scatter among various momentum states with probability S(p,pЈ). By noting the number of electrons in each bin at tϭ⌬t, P i j is evaluated by taking simple ratios. The transition matrix is precomputed, stored, then used repeatedly for simulation. Initial conditions are defined, and electron dynamics are then simulated according to Eq. ͑2͒ by simple matrix multiplication. Because P i j is evaluated once ͑or infrequently͒, the time-consuming Monte Carlo simulation is amortized over many simulations, and because a large number of electrons can be used to define P i j , statistical noise is virtually eliminated. As described above, the technique does not treat nonlinear effects such as band filling, electron-electron scattering, hot phonon effects, etc. without which one cannot simulate problems such as laser dynamics. We solve the problem in the following way. Consider electron-electron scattering. The transition rate for electrons in state p ͑of subband i͒ scattering to state pЈ ͑of subband m͒ from an electron initially at p 0 ͑in subband j͒ which scatters to p 0 Ј ͑in subband n͒ where បqϭ͉pϪpЈ͉, F i jmn is the overlap integral, 2 h is Plank's constant, and q s0 is the screening length. If we assume static screening, then the term in curly brackets is independent of the shape of the distribution function while the term in square brackets is shape-dependent. For quantum well lasers where the carriers have a near-Fermi-Dirac distribution, static screening is an excellent approximation, 6 although for photoexcitation experiments it is not. 2 The technique described here, however, is not restricted to static screening; the key point is to divide the scattering rate into shape-independent and shape-dependent factors. With electron-electron scattering rate term, Eq. ͑1͒ can be discretized ͓in a manner similar to Eq. ͑2͔͒ to obtain From the shape-independent factor, we precompute by direct Monte Carlo simulation a transition matrix, P i jkl , by assuming that the target electron is always present ͑ f l ϭ1͒ and that the final states are always empty ͑ f i ϭ f k ϭ0͒. During the simulation, one then uses this precomputed electron-electron transition matrix along with the dynamic estimate of f (t) to update the bin populations at tϭtϩ⌬t. Similar techniques are used to treat bandfilling. Transitions which do not occur because of the Fermi factors, renormalize P ii to P ii Ј .
Because the matrix P i jkl would be very large, we do not resolve the individual transitions in angle, but average over constant energy surfaces. As carriers relax in energy, they emit optical phonons, and the resulting nonequilibrium phonon population has a strong influence on carrier relaxation because the electronphonon scattering rate is modified. 1,2 A rate equation for 3D, unconfined phonons can be written as where t q is the optical phonon lifetime, n 0 (q) is the equilibrium phonon population, and G(q) is the phonon generation rate due to electron-phonon scattering. During the evaluation of the electron-phonon scattering process, we record the number and wavevector of the phonons emitted using a grid in phonon wavevector space similar to the one employed in Ref. 7. The population of optical phonons in each bin, N q i , is the updated every timestep ⌬t from where G q i is the generation rate for phonons with wavevector q i computed from the electron-phonon scattering events. Finally, we treat photons by a simple rate equation, where S m is the photon density in mode m ͑e.g., the longitudinal modes in a semiconductor laser͒, t p is the photon lifetime, and G op is the optical gain due to stimulated transitions as evaluated from the distribution function at time, t. Equation ͑7͒ can be discretized on a grid, which represent the various longitudinal modes, and solved by the same type of transition matrix technique; the optical gain, which depends on the carrier distribution, couples the photon transition matrix to the electron system. Having outlined the basic approach, we briefly describe the simulation process. It begins by discretizing momentum space for electrons, phonons, and photons. For the example calculation presented here, we used 200, 3 meV wide bins in energy and 30 equally spaced angle bins for electrons. The in-plane phonon wavevectors for each normal wavevector were discretized into 100 bins from 0 to q max ϭ2 ϫ10 9 m Ϫ1 . Electron-phonon ͑P i j ͒ and electron-electron ͑P i jkl ͒ transition matrices are precomputed by direct 2D Monte Carlo simulation. Initial conditions are defined; then the electron distribution is updated by splitting the effects of electron-phonon and electron-electron scattering. The electron distribution is first updated for phonon scattering using Eq. ͑2͒ adjusted for band filling, then the distribution is updated for electron-electron scattering using Eq. ͑4͒, and finally the phonon and photon populations are updated by solving Eqs. ͑6͒ and ͑7͒. The process continues at timesteps of ⌬t ͑typically 10-20 fs͒ until steady state is achieved. Note that the splitting of electron-electron and electron-phonon scattering processes is valid only if DtӶ.
To illustrate the technique, we present results for the simulation of a model optical amplifier pump-probe experiment. We consider a 150 Å wide GaAs quantum well biased to a steady-state electron density of 10 12 cm Ϫ2 . The steadystate electron population is nearly a Fermi-Dirac distribution, so the use of static screening for e -e scattering is appropriate. A LO phonon lifetime of 2 ps was assumed. The pump beam was modeled as a rectangular pulse of 150 fs duration, and the wavelength of the beam was adjusted to produce transitions in the lowest subband only. Two different wavelengths were assumed, one where the optical gain was positive, and a second in the absorption region. The results are displayed in Figs. 2 and 3.
As shown in Fig. 2 for the gain region, the pump beam depletes carriers from the lowest subband ͑subband 3͒, and there is a corresponding decrease in the gain of the probe beam. The gain recovers in about 0.7 ps, well within the experimental range. 8 Although not shown here, the population in subband 2 increases when the pump beam depletes carriers from the bottom subband. A percentage increase is expected because electrons are removed from subband 3, but an actual increase in the subband 2 population occurs because carriers cooler than average are removed by stimulated emission, so the subsequent scattering increases the electron temperature, and intervalley scattering from subband 3 to 2 increases. Figures 2͑c͒ and 2͑d͒ shows the simulated response when the wavelength is set in the absorption region. In this case, additional carriers are generated in subband 3, and the optical gain increases. It is interesting to note that the average carrier energy increases in the absorption region because the pump beam generates electrons with energies higher than the average energy of the ensemble. Note also that the optical gain recovers in about 0.4 ps, which is much shorter than the 0.7 ps observed in the gain region. This difference, which has also been observed experimentally, 8 can be attributed to electron-electron scattering. When a spectral ''heap'' is formed by absorption, e -e scattering is very effective because there are a large number of states to out-scatter to. On the other hand, when a spectral hole is burned in the gain region, e -e scattering is far less effective in filling up the hole because two empty final states are needed and the spectral hole only provides one. This also explains why it takes on the order of 1 ps to fill up a spectral hole as opposed to the 10's of femtoseconds that might be expected in a system dominated by e -e scattering.
In summary, we have described a new, Monte Carlo based simulation technique designed to examine a class of problems beyond the reach of direct Monte Carlo simulation. By using precomputed transition matrices, the timeconsuming Monte Carlo simulations can be performed once, with high precision ͑i.e., in Fig. 2, a 0.1% perturbation has been resolved͒, then used repeatedly for simulations. The technique virtually eliminates statistical noise, treats both fast and slow phenomena, and provides an excellent computational efficiency. For example, the transient simulation of Fig. 2 required less than 1 min on an IBM RS-6000/580 workstation. The matrix based technique is well suited for coupling to scattering matrix calculations 9,10 in order to address carrier transport and 3D to 2D capture issues in quantum well lasers. The technique should also prove useful for simulating pump-probe studies of quantum well dynamics and can be extended as discussed by Kuhn