Modeling Nonlinear and Chaotic Dynamics in Semiconductor Device Structures

We review the modeling and simulation of electrical transport instabilities in semiconductors with a special emphasis on recent progress in the application to semiconductor microstructures. The following models are treated in detail: (i) The dynamics of current filaments in the regime of low-temperature impurity breakdown is studied. In particular we perform 2D simulations of the nascence of a filament upon application of a bias voltage. (ii) Vertical electrical transport in layered semiconductor structures like the heterostructure hot electron diode is considered. Periodic as well as chaotic spatio-temporal spiking of the current is obtained. In particular we find long transients of spatio-temporal chaos preceding regular spiking.


I. INTRODUCTION
Semiconductors are complex nonlinear dynamic sys- tems which give rise to a variety of current instabili- ties when they are driven by strong electric fields.These instabilities often involve switching behavior, self-generated regular or chaotic current oscillations, current filamentation and solid-state turbulence.It has been known for a long time that current filaments are formed in semiconductors with an S-shaped current 3 density-field characteristic-.The stationary struc- ture and the possible nonlinear dynamic and chaotic oscillatory behaviour of these current filaments has 4-9 been widely investigated both experimentally and theoretically 1'I0-15 in a variety of semiconducting materials, e.g.p-Ge or n-GaAs at liquid Helium tem- peratures, or layered semiconductor structures like pin, pnpn, or heterostructure hot electron diodes.
In this paper we focus on two specific model sys- tems.Firstly, the dynamics of current filaments in the regime of low-temperature impurity breakdown is studied.We present 2D simulations of the nascence of filaments.Previous theoretical attempts to model cur- rent filaments were mostly confined to one-dimen- sional simulations where only the transverse spatial coordinate perpendicular to the current flow was taken into account.Intermittent and chaotic behavior of laterally traveling filaments in crossed electric and magnetic fields has been found 11.While there has been recent progress in the microscopic analysis of low-temperature impurity breakdown in terms of sin- 16 gle-particle and many-particle Monte Carlo (MC) simulations 17 for p-Ge and for n-GaAs 8'19, the spa- tio-temporal modes of the breakdown process have so far merely been investigated in a one-dimensional longitudinal model for p-Ge, a model that neglects the transverse spatial degree of freedom and therefore cannot explain filamentation2.In order to study the nascence of current filaments 21 it is necessary to com- bine the transverse and the longitudinal degrees of freedom for a realistic 2-dimensional sample geome- try with appropriately modelled contacts, and include detailed microscopic information on the generation-recombination kinetics obtained from MC simu- lations.
Secondly, vertical electrical transport in layered semiconductor structures like the heterostructure hot electron diode (HHED) 22 is considered.The HHED has been described in terms of a specific model for a GaAs/A1GaAs heterojunction23, which exhibits an S-shaped current-voltage characterisic and stationary or spiking current filaments.In 12 a simple generic model has been proposed.Periodic as well as chaotic spatio-temporal spiking bifurcation scenarios have been obtained24.Here we focus on the transient spatio-temporal chaos preceding regular spiking25.
Other semiconductor structures where nonlinear and chaotic high-field transport phenomena have been found include parallel transport in modulation- doped heterostructures associated with real-space transfer26, and vertical transport in a superlattice3.In both systems the current-field characteristic is N-shaped and leads to the formation of electric field domains and oscillatory instabilities.

II. DYNAMICS OF CURRENT FILAMENT FORMATION
As a first model system we study current filamenta- tion in thin n-GaAs films in the regime of low temper- ature impurity breakdown.We use a model which combines a drift-diffusion approach with Monte Carlo simulations of the generation-recombination (GR) kinetics21.The carrier density in the conduction band, and hence the current density, is determined by the GR processes of carriers between the conduction band and the donor levels.The experimentally observed S-shaped current density-field relation in the regime of impurity breakdown can be explained in terms of standard GR kinetics only if impact ioniza- tion from at least two impurity levels is taken into account .Therefore we model the infinite hydrogen- like energy spectrum of the shallow donors by the ground state and an "effective" excited state close to the band edge.In this case the state of the system can be characterized by the spatial distribution of the car- rier densities in the conduction band n(x,t) as well as in the impurity ground state and excited state n l(x,t), n2(x,t), respectively, where x is the spatial coordinate and denotes time.
The temporal evolution of n, n 1, and n 2 is then gov- erned by the rate equations dn -.j+Xn2-Tnpt+Xnn +Xnn2, (1)  (2) dt -X{n2+ Tnpt-Xnn2-T*n2+X*nl, (3) where Pt ND n n2 is the density of ionized donors, N D is the total density of donors, X[ is the thermal ionization coefficient of the excited level, T[ is its capture coefficient, X, X are the impact ioni- zation coefficients from the ground and excited level, respectively, X*, T* denote the transition coefficients from the ground level to the excited level and vice versa, respectively.Within the drift-diffusion approximation the current density j can be expressed as j e(nE + DUn) with the electron charge e, the dif- fusion constant D and the mobility t. 6" is the local electric field within the sample.
The electric field is coupled to the carrier densities via Poisson's equation V" where e is the dielectric constant and N N D N A holds with the compensating acceptor concentration N A. From equations (1)-(),( 4) the charge conserva- tion equation can be derived, given by V.J=O with J=e+j, (5) where J is the total cuent density composed of dis- placement cuent and conduction cuent densities.If the initial values of n, n 1, n 2 and satisfy Poisson's equation ( 4), ( 4) and (5) are not independent and we can substitute (4) by (5) for the numerical treatment of the time-dependent problem in drift-diffusion approximation.In many cases this approach turns out to be advantageous27-3.
The essential nonlinearities of the constitutive model equations (1)-(5) in the regime of low-temper- ature impurity breakdown are contained in the dependence of the GR coefficients upon n, nl, n2, and %.In order to derive these from a microscopic theory we have performed single particle MC simulations for a spatially homogeneous steady state 18.Thermal ioni- zation of the excited donor level, acoustic phonon-assisted recombination into the excited level (Lax-Abakumov), and impact ionization from both the ground and the excited donor level were included as band-impurity processes.The relevant intraband scattering processes were elastic ionized impurity scattering (Conwell-Weisskopf approximation) and inelastic acoustic deformation potential scattering.The microscopic rates of all band-impurity processes depend upon the carrier densities in the band and impurity states, which in turn depend upon the non- equilibrium carrier distribution function.To obtain these carrier densities, the MC method has to be com- bined self-consistently with the rate equations (1)- (3)  in the homogeneous steady state, where the GR coef- ficients X 1, X, T are calculated by averaging the microscopic transition probabilities (Pii, P:iZi, Prec for impact ionization from the ground state, the excited state, and capture, respectively) over the nonequilib- rium distribution function f(k), which is extracted from the MC simulation at each step: fd3kf(k;n, nl,n2,,g)pili(k, nl), X{(n, n l , n 2 , ) X l ( n ' n l ' n 2 ' -' ) --_ T:(n, nl,n,)-- f d3kf(k;n, nl rt2,-,)Prec(k, pt) npt (6)   Note that f, and hence X 1, X and T[, in turn depend parametrically on n, n 1, n 2 and E.An iteration proce- dure, where n and n 2 are expressed by their steady-state dependence on n and E, is used to solve the .aboveproblem self-consistently.
As a result the impact ionization coefficients X and X as well as the capture coefficient T depend not only on the local electric field E, but also on the electron concentration n.This dependence on n is associated with a higher electron temperature TP on the upper branch of the S-shaped n(E) characteristic as compared to the values Tle on the lower and the middle branch.
In the following our strategy will be to insert fitted analytical representations of the MC data into the macroscopic rate equations (1)-( 5).We use this approach in order to take into account as much detailed information as possible about the micro- scopic scattering processes, while still retaining man- ageable expressions.
We have simulated a square sample with side lengths L x L z 0.02cm representing a thin GaAs film (thickness 1.4 x 10-3cm), using an implicite finite element scheme29.We model point contacts by applying Dirichlet boundary conditions to two oppo- site regions of length L c 8 x 10-4cm at the centers of the sample edges parallel to the z-axis.At the con- tacts n is fixed to a value n D 5 1015cm -3 to model Ohmic contacts.All other boundaries are treated as insulating where the components of the current den- sity j and the electric field E perpendicular to the boundaries vanish.
We study the nascence of current filaments when the applied voltage is switched rapidly to a value above breakdown threshold so that the semiconductor is forced from the nearly insulating state to a highly conducting state.Within lps the voltage is linearly increased from U= 0V to U 0.48V corresponding to an average field of E 24V/cm.Due to the Ohmic nature of the contacts we find small regions in the vicinity of the contacts in which the electron concen- tration n in the conduction band is largely enhanced compared to the bulk where it is very low (Fig. la).Practically all carriers are bound in the donor ground state.The electron temperature in the whole Sample is equal to the lattice temperature T L (Fig. 2a).Due to the assumed voltage control the electric field E reacts quasi-instantaneously forming a dipole-like electric field distribution (Fig. 3a) and inducing enlarged areas of increased electron density at both the cathode and the anode.The current density j !/I (Fig. 4a) is very low and remains so during the voltage increase.Sub- a) !, ." i .
x u.uz 0 02 FIGURE Formation of a current filament in n-GaAs at 4.2K.The temporal evolution of the simulated electron density n(x, z) is shown for a) lps, b) 0.5ns, c) 1.0ns, d) t= 2.5ns.At 0 the voltage U 0.48V is applied FIGURE 2 Density plot of the temporal evolution of the electron temperature Te(x, z) as in Fig.
sequently impact ionization multiplies the electron concentration at the cathode (at x=0.02cm) and establishes a front that moves towards the anode (Figs.b, 4b).The propagation of the front is accom- panied by a high field domain associated with a slightly increased electron temperature T e (Fig. 2b).
Although the electric field behind the front is smaller than in front of it, for reasons of current conservation the increased electron density in regions passed by the front is almost conserved because recombination is a much slower process than generation.Hence impact ionization downstream is encouraged, whereas further generation upstream is inhibited.When the front meets the region of increased carrier density around the anode, a rudimentary filament is formed, albeit with a carrier density several orders of magnitude lower than that corresponding to completely ionized donor states of density N 5 x 1015cm -3 (Fig. lc).

Now donor impact ionization is becoming enhanced
in the rudimentary filament because the excited level n 2 is increasingly populated.The current density (Fig. 4c) and the electron temperature are significantly growing in the rudimentary filament (Fig. 2c) The potential distribution is being deformed (Fig. 3c) by the nascent filament.Finally impact ionization leads isolines: 0.02 V d) FIGURE 3 Temporal evolution of the potential (x, z), as in Fig. Potential isolines are shown to a uniform increase of electron density until the fila- ment reaches its mature state (Fig. ld) where almost all donors are ionized, and the carrier density corre- sponds to the upper branch of the homogeneous steady state characteristic n().Within the filament the excited donor level is much more highly populated than outside, nevertheless still only about 2 per- cent of the band carriers are trapped in the excited level, while the ground level is completely depleted inside the filament.Thus the population ratio between ground and excited level is inverted in the filament.Also the current density j and the electron tempera- ture T e are much larger inside the filament than out- side (Figs.4d, 2d).The deformation of the potential distribution is completed (Fig. 3d).
Hence from our simulation we find three stages of impact ionization breakdown: a stage of front creation and propagation from the injecting contact, i.e. the cathode (stage I) followed by a stage of stagnation after the front has reached the anode (stage II) and a final stage during which the rudimentary filament grows to a mature filament (stage III).FIGURE 4 Density plot of the temporal evolution of the current density j(x, z), as in Fig.

III. SPIKING AND SPATIO-TEMPORAL CHAOS IN LAYERED STRUCTURES
As a second model, we consider the following dimen- sionless reaction-diffusion system: )a(x,t) u a 2a )t (u-a) 2 q--Ta H X 2 (7) )a(x,t) 2u )t or(j0 (u a)) + D )x2 (8)   This model can be derived from the particular trans- port mechanism in the HHED by appropriate approximations 12.More generally, however, it describes charge transport in an SNDC element with an internal degree of freedom a(x,t) (In case of the HHED this is the normalized interface charge density between two layers), whose dynamics is described by eq.( 7).Eq.( 8) describes the dielectric relaxation of the normalized voltage u across the device.The device is driven under current controlled con- ditions with an external normalized current J0 and an external capacitance Cp (related to the dimensionless parameter o-1/(C + Cp),where C is the intrinsic capacitance of the element) applied parallel to the sample.In terms of nonlinear dynamics J0 and a are control parameters.T is an internal system parameter.We use Neumann boundary conditions a/x Ou/x 0 for x 0, L, where L is the size of the system.The equations are solved with finite dif- ferences and a forward Euler algorithm.For a wide range of parameters o, T, D, L the homogeneous fixed point successively exhibits a Hopf bifurcation (for homogeneous fluctuations) and a Turing instability (for inhomogeneous fluctuations) with increasing j031'32.
In part of this range periodic spatio-temporal spiking is found12, where a spatially periodic pattern forms and vanishes periodically in time as shown in Fig. 5a.
In terms of the semiconductor model the quantity u a, which we have plotted, corresponds to the current density which is the physical quantity of interest.
Throughout the following we use T=0.05, ot 0.02, D 8, J0 1.21 as parameters.For this case the Hopf and Turing instability occur at J0 1:185 and 1.204, respectively.The system size L is the only parameter which is changed.We use the initial condi- tions a(x, O)= a* + r(x), u(x, O)= u*, where r(x) is given by a random sequence from the interval [-0.5,0.5] at each position of the discretized spatial lattice.Then the periodic spiking behaviour is found after a long transient irregular behaviour.both a part of the transient (b) and the periodic behav- iour (a) for L 600.The periodic pattern is character- ized by 3 spikes which arise periodically at two different positions separated by half a period.(A spike with a maximum located at the boundary is counted as 1/2.) The transient behaviour is dominated by similar spikes occurring, however, at irregular positions.In order to visualize the dynamics more clearly we have plotted the position of the local maxima of (u(x, t)-a(x, t)) exceed- ing the value 3 in the space-time-plane (Fig. 5c).Fig. 5d shows that the same type of transient chaotic behaviour appears for L 2000, too.
We now define the transient time " where the peri- odic state is reached within a certain accuracy.Due to the approximately exponential increase of the tran- sient times with the system length 25 a large system will practically not exhibit any periodic behaviour.The spatio-temporal properties of the system in this transient chaotic state can be examined by consider- ing the temporal spreading of a localized perturbation.
To this purpose we add a perturbation to a(xc/, tel) at a distinct location at a given time d.For > td we calcu- late the trajectory for both the unperturbed and per- turbed initial condition at ta.In order to visualize the deviations resulting from the perturbation we have plotted the logarithm of the difference between these two trajectories in Fig. 6.We find that the perturbation spreads with a constant velocity which seems to be independent of the width L. There are some points, especially in the simulation shown for L= 1200, where the propagation stops and the perturbation is reduced dramatically.Afterwards the residual pertur- bation grows again.At the initial stage the propaga- FIGURE 6 Temporal evolution of a localized perturbation in the transient chaotic phase for L 600 and L 1200.The perturbation is localized at 5 000 and 15 000, respectively, in the centre of the system.The grey scale corresponds to the difference of u a for the perturbed and the unperturbed system, respectively tion is roughly parabolic rather than linear.This behaviour might be due to pure diffusion.Only when the diffusive spreading slows down, can the linear propagation be detected.Fig. 7 shows the asymptotic behavior of a slightly modified model24'33, where the local diffusive cou- pling of u is replaced by a nonlocal coupling with a(x, t)dx.In Fig. 7b the minima of u(t) cor- 0 responding to the current spikes are plotted versus the control parameter J0 exhibiting a period-doubling sce- nario to chaos.Additionally, an intermittency scenario occurs with decreasing J0 at Joe 1.267.
As further characterisation of the chaotic behaviour we have calculated the Lyapunov exponents.Fig. 7a shows the two largest Lyapunov exponents (including zero) in the interval 1.24 < J0 < 1.33.Our system dis- ECKEHARD SCH(3LL plays periodic or chaotic spiking for these parameters.
For the periodic case the maximum exponent is zero, while chaotic behaviour is characterized by at least one positive exponent.At the period-doubling bifur- cations the second Lyapunov exponent also vanishes.The Lyapunov spectrum closely matches the bifurca- tion diagram in Fig. 7b.
In conclusion, we have presented a generic model for a novel selforganized spatio-temporal behaviour of semiconductors which we call spiking.Transient and asymptotic chaotic spiking has been found.Our model reflects the nonlinear dynamic behaviour of current filaments in pin diodes 5 and other semicon- ductor devices 13.Because it is based on very simple principles of current continuity and dielectric relaxa- tion, we suggest that it might be useful in the analysis of nonlinear spatio-temporal dynamics in a wide vari- ety of semiconductor structures.
am indebted to my collaborators S. Bose, M. Gaa, R.E.Kunz, and A. Wacker who have largely contrib- uted to the results presented.I wish to thank H. Gajewski and R. Ntirnberg from the Weierstrass-Insti- tut fiir Angewandte Analysis und Stochastik, Berlin for their help in adopting our equations to their algorithms.Partial support from Deutsche Forschungsgemeinschaft is gratefully acknowledged. X*nlXlnnl,

FigFIGURE 5
FIGURE 5 Spatio-temporal spiking of current filaments.(a) asymptotic periodic spiking, (b) transient chaotic spiking, (c) positions of the spikes as a function of time t.(d) exhibits a part of the transient chaotic bahavior for a different simulation with a larger system size(L 2000)

FIGURE 7
FIGURE 7 Asymptotic chaotic spiking in a model with global cou- pling for 0.035.a) Lyapunov spectrum, b) bifurcation diagram of the minima of the voltage u(t) as a function of the applied cur- rent density J0