Chaotic dynamics of electric-field domains in periodically driven superlattices

Self-sustained time-dependent current oscillations under dc voltage bias have been observed in recent experiments on n-doped semiconductor superlattices with sequential resonant tunneling. The current oscillations are caused by the motion and recycling of the domain wall separating low- and high-electric- field regions of the superlattice, as the analysis of a discrete drift model shows and experimental evidence supports. Numerical simulation shows that different nonlinear dynamical regimes of the domain wall appear when an external microwave signal is superimposed on the dc bias and its driving frequency and driving amplitude vary. On the frequency - amplitude parameter plane, there are regions of entrainment and quasiperiodicity forming Arnol'd tongues. Chaos is demonstrated to appear at the boundaries of the tongues and in the regions where they overlap. Coexistence of up to four electric-field domains randomly nucleated in space is detected under ac+dc driving.

Self-sustained time-dependent current oscillations under dc voltage bias have been observed in recent experiments on n-doped semiconductor superlattices with sequential resonant tunneling. The current oscillations are caused by the motion and recycling of the domain wall separating low-and high-electric-field regions of the superlattice, as the analysis of a discrete drift model shows and experimental evidence supports. Numerical simulation shows that different nonlinear dynamical regimes of the domain wall appear when an external microwave signal is superimposed on the dc bias and its driving frequency and driving amplitude vary. On the frequency -amplitude parameter plane, there are regions of entrainment and quasiperiodicity forming Arnol'd tongues.
Chaos is demonstrated to appear at the boundaries of the tongues and in the regions where they overlap. Coexistence of up to four electric-field domains randomly nucleated in space is detected under ac+dc driving.

I. INTRODUCTION
Negative differential conductivity (NDC) in weaklycoupled narrow-miniband semiconductor superlattices (SL's) results in the formation of electric-field domains, which have been studied both experimentally 1 and theoretically. [2][3][4][5][6] The domains are stable if the doping or photoexcitation are large enough to form a stationary charge accumulation layer (the domain wall). The domain wall moves from well to well as the bias increases and gives rise to the jumps (discontinuities) of the current in the stationary I-V characteristics.
When the carrier density is not sufficiently large to form stable domains, but it is large enough for the uniform field distribution to be unstable, periodic timedependent oscillations of the current under fixed dc bias are possible. These oscillations have been observed in recent experiments on GaAs/AlAs SL's; they are damped for undoped photoexcited samples 7 and undamped for doped SL's without photoexcitation. 8,9 According to a discrete drift model, 5 the current oscillations are caused by a periodic motion of the domain wall over a few periods of the SL, 10 , and this is confirmed by photoluminescence measurements. 8,9 Note that damped oscillations (with frequencies up to 20 GHz) have been observed experimentally by Le Person et al. 11 in a photoexcited wide-miniband SL with strong interwell coupling. These oscillations have been interpreted in terms of dipole charge waves. We do not consider the miniband transport regime in the present paper.
This situation is reminiscent of that found in bulk GaAs with NDC caused by the intervalley transfer of electrons, where the well-known Gunn oscillations mediated by high-field domain dynamics may appear under dc voltage bias. 12,13 There are several important differences between the Gunn oscillations and those observed in the SL's. A significant difference is that the space charge wave is a dipole in the case of the Gunn oscillations (a domain of high electric field with two charge layers: accumulation and depletion) and a charge monopole in the SL current oscillations (a wave of one charge accumulation layer, or the domain wall, connecting two electric-field domains). Another difference is that the Gunn waves are generated close to the injecting contact whereas the domain walls appear clearly inside the SL 10 . Notice finally that the different transport mechanisms determining the velocity of waves and the characteristic frequency of the oscillations give rise to different limitations in the performance of possible devices. While for the Gunn effect the oscillation frequency is limited by a parameter of the material (the intervalley scattering time), for the SL's it is determined by the tunneling time. Hence, the oscillation frequency can be varied by tuning the growing parameters (barrier widths, etc.) and/or the dc bias. The bias regions between different resonance peaks give rise to quite different frequencies, ranging between hundreds of KHz and several GHz over wide temperature ranges including room temperature. 9 In our previous paper 14 we predicted the appearance of chaos by applying an additional small ac signal to the dc voltage bias in the regime where the SL exhibits self-sustained oscillations. In particular, when the ratio between the natural and the driving frequencies is fixed at the golden mean, chaotic current oscillations and strange attractors with a multifractal measure have been obtained. 14 Here we present a detailed numerical study of nonlinear current oscillations in SL's over the whole driving frequency -driving amplitude parameter plane. Based upon a relatively simple (but otherwise self-consistent) discrete drift model, 5 the shape of frequency-locked regions (Arnol'd tongues), chaotic and quasiperiodic regions is determined, thereby providing a global bifurcation picture. The largest Lyapunov exponents are calculated to identify chaotic solutions and to characterize their fractal dimension. The critical line where the Arnol'd tongues begin to overlap and chaos appears is found at very low values of the amplitude of the driving force, unlike other periodically driven semiconductor systems which also display complex nonlinear interaction between internally generated oscillations and an external ac signal (e.g., GaAs Gunn diodes, 15 p-Ge, 16,17 ).
An additional point worth noticing is that the physical mechanism of the NDC in the SL system (responsible for chaos) is based on a pure quantum-mechanical effect (resonant tunneling) which is absent in the classical limit. Thus we may have a case of quantum chaos without classical counterpart 18 which might be easier to detect experimentally than that proposed by Jona-Lasinio et al. 18 . In both cases the self-consistent mean-field potential created by the charges is instrumental in inducing chaos, which is also the case when a quantum system is coupled to purely classical subsystems having a widely different time scale. 19 All these phenomena are different from the so-called quantum chaos, 20 i.e., the behavior of quantum systems whose classical counterpart is chaotic.
Recently chaotic dynamics in semiconductor superlattices with miniband transport has also been considered within a classical balance-equation approach. 21 In this case the chaotic dynamics is due to the negative effective mass of the electrons in the minibands for the regime of the Bloch oscillations under an external ac electric field (field self-consistency is also taken into account). In this completely different physical situation, the field was assumed to be spatially homogeneous 21 which precludes consideration of the spatial structure of chaos, in contrast with our work. 14 The paper is organized as follows. In Sec. II we describe the physical model and its governing dynamical equations for a dark doped SL. In Sec. III we recall the main qualitative features of the self-sustained oscillations in a dc-voltage biased SL. The general behavior of the self-oscillating SL subjected to external ac and dc biases is discussed in Sec. IV. Frequency-locked solutions (Arnol'd tongues) are analyzed and compared with those for the Gunn diode. The chaotic dynamics at the golden mean ratio between the natural and the driving frequencies is described in Sec. V. Different techniques to analyze chaotic behavior (Poincaré mapping, bifurcation diagram, Lyapunov exponents, phase plots, Fourier spectra, first return map, etc.) are discussed. Section VI draws the main conclusions of this work. In the Appendix we present the algorithm to calculate the Lyapunov exponents for our system and discuss the Lyapunov dimension calculated by the Kaplan-Yorke formula.

II. THE PHYSICAL MODEL
We consider a set of weakly interacting quantum wells (QW's) under high voltage biases, so that the electron states are localized in the wells and the tunneling process is sequential. The relaxation from excited levels to the ground state within each QW (∼ 1 ps) is considered to be much faster than the tunneling process between adjacent QW's (∼ 1 ns). Therefore, a single QW reaches a local equilibrium between two consecutive tunneling processes, and the state of the system can be characterized by values of the electric field E i (t), and the electron density n i (t), with i = 1, . . . , N denoting the QW index. 5 The one-dimensional equations governing the dynamics of the doped SL are the Poisson equation averaged over one SL period l, Ampère's equation for the balance of current density, and the voltage bias condition 14 Here ǫ, e, and N D are the average permittivity, the electron charge, and the average doping density, respectively. The total current density J(t) is the sum of the displacement current and the electron flux due to sequential resonant tunneling en i v(E i ). The effective electron velocity v(E) (proportional to the tunneling probability) exhibits maxima at the resonant fields for which the adjacent levels of neighboring QW's are aligned 22 (Fig. 1). Notice that the mechanism of sequential resonant tunneling is responsible for the NDC region in the velocity curve which in turn gives rise to the current instabilities in the SL. The voltage V (t) in (3) is the sum of a dc voltage V b and an ac microwave signal of relative amplitude A and driving frequency f d : The boundary condition at the first contact ǫ(E 1 − E 0 )/(el) = n 1 − N D = δ allows for a small (δ ≪ N D ) negative charge accumulation in the first well. The physical origin of δ is that the n-doped SL is typically sandwiched between two n-doped layers with an excess of electrons, thereby forming a n + -n-n + diode. 8 Then some charge will be transferred from the contact to the first QW creating a small dipole field that will cancel the electron flow caused by the different concentration of electrons at each side of the first barrier.
Note that the rate equation for the electron density for our model can be derived by differentiating (1) and using (2): This is the equation of charge conservation under sequential resonant tunneling between neighboring QW's. Finally, the current density J(t) in the external circuit under voltage bias condition can be obtained from the time-derivative of (3) and Ampère's law (2) in the form where c V = ǫ/(lN ) is the intrinsic SL capacitance per unit cross-sectional area.
To study our equations, it is convenient to render them dimensionless, using characteristic physical quantities. As the unit of the electric field E = E/E 1−2 , we adopt E 1−2 = ∆/(el), with ∆ being the energy separation between the first and the second electron subbands [the maximum of v(E)]. This yields the characteristic charge density n 0 = ǫ E 1−2 /(el) used to normalize the doping ν = N D /n 0 and the electron density. The dimensionless velocity v(E) is obtained normalizing v(E) by its value at E = E 1−2 , where it has a local maximum due to resonance in tunneling (see Fig. 1). The others dimensionless quantities are defined as follows: the time τ = t/t tun where t tun = l/v(E 1−2 ) is the characteristic tunneling time; the dc bias V = V b /E 1−2 lN ; the ac bias amplitude a = AV; the driving frequency ω = 2πf d t tun . Now substituting (1) and (5) into (2) we obtain a system of N equations for the electric-field profiles with the boundary condition E 0 = E 1 − νδ and initial conditions E i (0) = V, ∀i.
The system (6) is solved numerically by the fourthorder Runge-Kutta method using 4000 time steps per one oscillation period. We start with a uniform initial field profile and solve the equations for dc bias. After a short transient, the self-sustained oscillations set in and we switch on the ac part of the bias. Every time step we calculate the electric-field distribution over SL E i (t) and the total current density J(t) [from Eq. (5)]. The main features of our numerical results are as follows.

III. SELF-SUSTAINED OSCILLATIONS UNDER DC BIAS
For pure dc case (a=0) the SL exhibits undamped time-periodic current oscillations when the doping density is between some critical values ν * and ν * * . In Fig. 2 the temporal behavior of the current starting from t=0 is shown for different doping densities. Below ν * ≈ 0.066 the electric-field distribution over the SL remains almost uniform, so that the current exhibits no oscillations. Above ν * * ≈ 0.175, after some transient period of the order of t tun N , stable stationary electric-field domains are formed. This results in an increase of the current and its saturation with time.
For in-between values ν * < ν < ν * * and appropriate dc voltage bias, we find undamped oscillations with frequency slightly dependent on the doping. In Fig. 3 the frequency of the natural oscillations versus doping concentration is shown for different boundary parameters δ. One can estimate the natural frequency to be approximately equal to the inverse total tunneling time 1/(t tun N ), which gives, after substitution of our value for t tun ≈ 2.7 ns and N = 40, a frequency about 10 MHz in close agreement with the value observed in experiment. 8 This estimation improves as N increases. All the results presented below were obtained for the value of δ = 10 −3 .
The spatial structure of the oscillations can be seen from the electron density and the electric-field distributions of Fig. 4. Starting from uniform distribution at time τ =0, when the dc bias is turned on, after a transient period τ ≈ 3, a charge accumulation wave (monopole) is created and then it moves toward the correspondent contact. Depending on the applied voltage, it may or may not reach the end of the SL before it disappears and a new monopole is formed starting a new period of the oscillation. Simulations clearly show monopole recycling with two monopoles coexisting during some part of one current oscillation period [see Fig. 4(b)]. The period of the oscillation is mainly determined by the total travelling time of the monopole across the SL t tun N , as was shown in Fig. 3.
The total current density J(t) determines the average state of the system (see Eq. 5), nevertheless the information on the spatial dynamics is presented in the time dependences of the current. In the curves of Fig. 2 small current spikes can be seen during the initial stage of the domain formation, which are correspondent to the wellto-well jumping of the charge accumulation layer. The similar, but less clear, fine structure can be resolved in the periodic part of the current oscillations (become more pronounced after taking the derivative dJ/dt). From the number of spikes per period (∼6 for our curves) it is possible to estimate the number of QW's the monopole moves across.
The existence of the threshold value for the doping ν * , above which current instabilities take place, can be understood from the charge conservation law (4) rewritten in the continuum limit as Taking spatial derivative and using the Poisson law for ∂E ∂x one gets where the sign of the velocity derivative over the field v ′ E is crucial. In the NDC region where v ′ E < 0 the last term in (8) gives rise to an exponential growth of the charge with the characteristic time τ ǫ ∼ ǫ/(eN D |v ′ E |). For v ′ E > 0 the same term gives charge relaxation to quasineutrality n ≈ N D . Taking into account the characteristic transit time τ t ∼ lN/v from the convection term of (8) we can get significant growth of charge if τ t > τ ǫ , that is N D N > ǫv/(el|v ′ E |). This condition is similar to the critical nL-product for the Gunn effect. 12 The treatment above is of course too rough, but it clarifies the influence of different parameters on the instability. Rewritten in our dimensionless units as this inequality implies: (i) for fixed number of the SL periods N there should exist the critical value for the doping ν * above which instabilities of the current may occur, and (ii) for any fixed doping ν one can bring the system to instability by increasing the length of the SL N . Precise bounds for ν * and ν * * may be found elsewhere. 23 Notice the importance of the NDC region in the velocity curve v(E) and the positive sign of ∂E ∂x at the boundary, which are both necessary for monopole recycling. Besides the oscillations of the domain wall due to the travelling charge wave, the field in the low-and highfield domains also oscillates [see Fig. 4(a)]. The field behind the monopole is uniform in space, up to a small correction of the order of δ near the boundary. When the current reaches its maximum then the field behind the monopole takes values on the NDC region, and those corrections increase exponentially in time, nucleating a new monopole. 24 The condition (9) in the strong limit νN ≫ v/|v ′ E | implies a large separation between the characteristic times τ t and τ ǫ . This can be exploited to perform an asymptotic analysis of the current oscillations, which gives a reasonable approximation to the numerical simulations for long SL's with N > 100. 24 A key ingredient of the analysis is that the domain walls become shock waves (i.e., discontinuities moving towards the right and separating quasineutral regions of low and high electric field), whose velocity obeys an equal-area rule. 10

IV. FREQUENCY-LOCKING UNDER DC AND AC BIAS
Self-oscillating systems that are forced by an external oscillating signal represent an important class of coupled oscillators. An inherent feature of periodically forced nonlinear system is that the actual oscillation frequency depends on the amplitude of the forcing. Therefore both the frequency f d and the amplitude of the driving a can be used as control parameters to study nonlinear dynamics of our system.
For relatively low amplitudes we could expect two possible effects of the forcing according to the ratio between the driving frequency f d and the natural frequency f 0 . If this ratio is a rational number the regime of entrainment or mode-locking will be observed, if it is irrational the regime of quasiperiodicity will occur because the frequencies are incommensurate. 25 The statement above is confirmed by our numerical calculations. The phase diagram of Fig. 5 shows the distribution of frequency-locked solutions (Arnol'd tongues) over the frequency -amplitude parameter plane. Within each tongue, the current is a periodic function of time whose actual frequency f s does not coincide with f 0 nor with f d in general. Instead it is usually related to one of their harmonics. Frequency-locking can be considered as the trend of the driven system to keep its frequency close to that of the unforced oscillation, f 0 , by taking the value of the harmonic of f d which is closest to f 0 . The periods of the actual oscillations are shown by numbers in Fig. 5. Notice that the frequency-locked solutions are interspersed with regions of quasiperiodicity where both frequencies f 0 and f d coexist. At a=0 the Arnol'd tongues arise in intervals around rational ratios of f d /f 0 , and here f s = f 0 . As a > 0 increases, the tongues open up and eventually intersect. Then more complicated dynamic behaviors (e.g., chaos) usually occur as the solution jumps between the various overlapping modes in an erratic way.
Frequency locking is a typical nonlinear phenomenon that can be found in different periodically driven systems. 25 In our model the Arnol'd tongues begin to overlap already at very low driving amplitudes (a ∼ 0.01). This is quite different of the case of the Gunn diode studied by Mosekilde et al. 15 In addition, we observe 'the red shift' in the dependence of the frequency f s with the amplitude a, which manifests itself in the bending of the locked regions. The last phenomenon is also absent in the bulk Gunn effect and it might be a consequence of the delay of the natural oscillations when the forcing is increased. We have a qualitative argument to understand why the overlapping of Arnol'd tongues (and therefore the critical line for chaos) occurs at relatively low ac-voltages. Notice that our system has important voltage and time scales due to its discreteness: recall that there are spikes superimposed to the oscillation of the current tracking the jumps of the domain wall from one QW to the next one. It is plausible that overlapping of Arnol'd tongues occurs when the ac voltage amplitude is equal to the time-averaged voltage drop per period, V b /N , needed to balance the typical discrete scale of voltage. The latter is equal to the energy separation between subbands divided by the charge of the electron, E 1−2 l = ∆/e. 26 Unlike the ac-driven bulk semiconductors, one could thus expect overlapping of resonances for a SL driven with an ac amplitude V b A ∼ ∆/e ∼ V b /N , which is much smaller than the dc bias V b . This value seems to give a reasonable estimation of the lowest critical line for chaos as shown in Fig. 5.
Let us now further describe the Arnol'd tongue diagram of Fig. 5. Between the main tongues correspondent to n : 1 locking, smaller n : m intervals are found, representing more complex entrainments. For instance, between the 3 : 1 and 4 : 1 tongues one can see 10 : 3, 7 : 2, and 11 : 3 ratios, between the 4 : 1 and 5 : 1 tongues there are 13 : 3, 9 : 2, and 14 : 3 ratios, and so on. Scanning over both parameters a and f d to analyze the system of 40 equations was very time-consuming, because we calculated the largest Lyapunov exponent (see below) to distinguish between quasiperiodic and chaotic solutions. Thus we restricted ourself to scanning with a frequency step ∆f d = 0.05f 0 , which allowed us to detect the tongues n : m with at most m=3. Higher order locking and much smaller new chaotic regions at the boundaries of the tongues are expected to appear if a smaller step ∆f d is used. An additional important result in the mode-locking diagram of Fig. 5 is that the richest structure of different mode mixing and chaos is concentrated between the 2 : 1 and 3 : 1 tongues. In the next section we consider in more detail the results obtained when the driving frequency is fixed within that region, namely, at the inverse golden mean ratio.

A. Bifurcation diagram
Here we shall consider the driving amplitude a as the control parameter, fix the driving frequency as f d = Gf 0 (where G = ( √ 5 + 1)/2 ≈ 1.61803... is the inverse golden mean ratio 29 ), and calculate the current J as a function of time.
To detect and visualize the chaotic regions in parameter space, we need to define a Poincaré mapping. The current is a good measure of the amplitude (norm) of the solutions, which is illustrated by the use of current versus voltage characteristics as bifurcation diagrams 13 . Let T d = 1/f d be the driving period. We adopt as our Poincaré mapping (for each value of a) the current at times τ m = mT d , m = 0, 1, . . ., (after waiting enough time for the transients to have decayed). 14 The result is the bifurcation diagram in Fig. 6(a) which is constructed as follows. For each a we compute J m = J(mT d ) until the solution becomes periodic within a 10 −5 accuracy. At that time, we stop the simulation and depict all the J m corresponding to one period of the solution. If the solution is not periodic, we eliminate the first 500 transient points J m and depict the next 200 points. Thus, aperiodic solutions can be very easily distinguished from periodic ones in Fig. 6(a) by the large number of points (periods) corresponding to each value of a. Notice that period-2 orbits span the widest parameter region from a ≈ 0.01 up to a ≈ 0.085. Period-doubling cascades can be seen, which points out the existence of chaos near their accumulation points.
The next important consideration is to discriminate chaotic from quasiperiodic regions, both distinguished by having a large number of points in the bifurcation diagram. This can be achieved by computing the largest Lyapunov exponent λ 1 . For chaotic regions λ 1 > 0, which indicates exponential divergence of nearby trajectories. For periodic solutions of our nonautonomous system λ 1 < 0, while for quasiperiodic solutions λ 1 = 0.
In Fig. 6(b) we present the first and the second Lyapunov exponents calculated by the algorithm explained in the Appendix. We see that the system starts being quasiperiodic (λ 1 = 0) at a = 0, as it should be for an irrational frequency ratio. Then at a ≈ 0.005 it locks to a period-5 orbit (λ 1 < 0) terminated by some chaotic windows at a ≈ 0.01 (λ 1 > 0). Notice the first appearance of chaos at relatively small driving amplitude (∼ 1% of V b ). After a > 0.085 several chaotic windows can be seen, and then the solution becomes again quasiperiodic (as it was at a = 0) before locking to the driving frequency f d at a ≈ 0.145.
More insight into the transition between chaotic and nonchaotic regions of the bifurcation diagram can be obtained by sweeping-up (or sweeping-down) through it, using the electric-field profile stored at the end of each simulation (for a = a 0 ) as an initial condition for the next simulation (for a = a 0 + ∆a). This approach is close to the process of experimental sweeping-up measurements. By this technique we investigated some critical points in the bifurcation diagram and found different pictures when sweep-up and sweep-down runs were made, thus demonstrating hysteresis. For example, such hysteretic behavior we obtained for the transition point between period-2 orbit and chaos at a ≈ 0.085, which points out to the existence of a subcritical bifurcation.   Fig. 7(b),(d)] the FS be-come very irregular with a large number of peaks, which could be considered as an additional method to detect chaos. The FS for chaotic solutions should not be necessarily continuous. In our system the sharp frequency components are also present in the FS as in the case of the familiar Rössler attractor. 30 The Poincaré mapping used to obtain the bifurcation diagram of Fig. 6(a) can be understood from Fig. 7 as the successive crossing of the orbit through the line V=1.2, where ac part of the voltage crosses zero (at upper values of the current corresponding to the increasing voltages). For aperiodic solutions those crossing points are distributed over some interval (or intervals) of the current.

C. First return map
One of the main quantities observable in experiment is the current density in the external circuit J(t). By sampling the current J(t) each driving period T d , one obtains the data set J m . This set can be analyzed by means of the first return map plotting J m+1 as a function of J m . After some transient time the resultant attractor for a n-period solution will be just the n separate points (0-dimensional object), while for an aperiodic (chaotic or quasiperiodic) solution it will be represented by a higher dimensional object.
A closed smooth loop with regular distribution of the points indicates quasiperiodicity (Fig. 8a). The chaotic attractor is a layered (sometimes folding) structure with varying density of the points on different regions (Fig.  8b), and it can be characterized by the multifractal dimension D q . 31 Sometimes the chaotic attractor contains several separate branches almost continuously filled as in Fig. 8c. In this particular case the attractor has five branches, indicating that the solution is close to period-5 locking for this value of the control parameter a. The chaotic solution has the same symmetry as the frequencylocked solution which lies nearby.

D. Spatiotemporal aspects of chaos
So far we have only characterized the temporal aspects of our chaotic current oscillations leaving aside the spatial dependence of electric field and charge density. This dependence is a very characteristic feature of our system because the oscillations of the current are due to the dynamics of nonlinear travelling charge waves. Thus we analyze here what are the field and charge density profiles in the dynamical regimes of interest described previously.
Under pure dc voltage bias, the SL exhibits timeperiodic current oscillations accompanied by a periodical recycling of the monopole charge wave (domain wall) in space as described in Sec. III. When an ac signal is superimposed on the dc voltage bias, the current can become chaotic for particular values of the control parameters. For that case the motion of the charge waves becomes chaotic too, as shown in Fig. 9(a). The ac part of the voltage causes the electric field in different parts of the SL to take on values in the NDC region at different instants of time. This results in irregular amplification or damping of the charge disturbance at different QW's.
There are two significant new features for a SL under ac and dc voltage bias as compared to the pure dc case. The first distinctive feature is that the monopoles may be generated closer to the beginning of the SL, thereby leaving more room to their motion towards the end of the SL; see Fig. 9(a). Then the peak-to-valley ratio in the current oscillations increases as it can be appreciated in Fig. 7. Amplification of the ac signal is not linked to chaos since it can also be observed in cases where there is frequency locking and an uniform electric field profile inside the SL. 32 In fact, the amplification is most pronounced at the instability threshold for the doping ν * , where the current oscillations appear. Then the electric field is almost uniform all the time.
The second new feature is the qualitative change of the traveling wave picture for longer SL. For short SL (N < 80), at most two domain walls (separating three domains) can coexist at a given time (see Fig. 9(a)), whereas up to three coexisting domain walls (separating four domains) can be observed for long SL's during certain time intervals (for the case of Fig. 9(b) at τ ≈ 41.1; 46.5). In Fig. 10(a)-(c) we show the presence of one, two or three electric charge monopoles at different times for a 200-well SL under ac and dc voltage bias. We can understand this as follows. The width of a domain wall (≈6 wells for our particular set of parameters) is determined by the velocity profile and it is approximately independent of the SL length. 5 Then the monopole has more room to move on a longer SL (under pure dc bias, monopole recycling takes place on about 12 periods of a 40-well SL as compared to about 90 periods of a 200-well SL). The disturbance caused by the ac bias on monopole motion is thus much larger for a longer SL, which results in quite different spatiotemporal structures under ac and dc bias [compare the charge densities of Figs. 9(a) and (b)]. For a short SL the chaotic behavior can be associated with a chaotic domain wall dynamics that resembles the dc voltage bias case: most of the time there is only one domain wall which moves to the end of the SL and disappears, and about that time another monopole is generated, sometimes quite close to the beginning of the SL [see Fig.  9(a)]. On the other hand, for a long SL the graph of the electron density is disjoint: In addition to long-living waves traveling over almost the whole SL, there are shortliving waves existing only at the beginning of the SL. The two types of waves are distributed chaotically in space.
We have not found more than three monopoles by increasing further the SL length N . This may be due to the fact that there is not enough charge in the SL (determined by doping) to provide more than three jumps in the electric field. It might be possible to observe more complicated field structures in parameter regions where the current instabilities exist at much higher electron densities.
The spatially chaotic nature of the solutions can be illustrated by picking two far-away QW's and depicting the simultaneous values of the electric field at them after each period of the driving force T d . The resultant attractors will be very similar to those obtained for the first return map J m+1 − J m described in Sec. V C (cf. Fig.2 in Ref. 14 with Fig. 8 of the present paper).

VI. CONCLUSIONS
Dynamic properties of the high-field transport in weakly-coupled superlattices under ac and dc biases has been studied numerically within the simple selfconsistent discrete drift model. A nonlinear interaction between an internally generated periodical motion of the accumulated charge wave and an external microwave signal gives rise to a frequency-locking, quasiperiodicity or chaos depending on the external driving parameters. The calculated phase diagram of the frequency-locked solutions (Arnol'd tongues) shows the following new features: (i) the first overlapping of tongues giving transition to chaos occurs under very weak driving (a ∼ 0.01); (ii) there is a bending of tongues under strong driving ('red shift' in frequency) which we associate with the delay in natural oscillation frequency, and that is a consequence of our spatially extended system.
Besides the chaotic dynamics, the microwave forcing was found to give an amplification of the self-sustained current oscillations, which is most pronounced at the instability threshold for the doping ν * , where the oscillations appear.
It would be of interest to verify our predictions by making measurements in currently available n-doped GaAs/AlAs samples forming n + -n-n + diodes. These samples exhibit self-sustained oscillations under pure dc voltage bias, 8 and are thus suitable candidates for observing chaos when an appropriate microwave signal is superimposed on the dc voltage bias. More highly doped samples present multistable stationary solution branches (corresponding to coexisting static electric field domains 1 ) in their current -voltage characteristics under dc voltage bias. Studying the response of these samples to ac and dc voltage bias would also be of interest. Preliminary calculations show rich spatiotemporal behavior including chaos.

APPENDIX: LYAPUNOV EXPONENTS
Consider the N -dimensional phase space for the electric-field vector E i . The evolution of its trajectory E i (τ ) is described by the set of nonlinear equations (6). Since we want to know how trajectories infinitesimally close to the real trajectory of our system diverge from it, we can use a linearized version of the equations (6): Hereê i , with i = 1, . . . , N , is the vector of the disturbances of the electric field. The velocity v(E τ i ), its derivative v ′ (E τ i ), and the electron density n τ i = E τ i − E τ i−1 + ν are calculated on the fiducial trajectory E i (τ ). The boundary condition for the system (A1) isê 0 =ê 1 .
To calculate the largest Lyapunov exponent λ 1 is very easy. We take an arbitrary initial vectorê i (0). Then both the nonlinear system (6) and the linearized system (A1) are advanced forward in time simultaneously.
Denoting by ê i (τ ) the euclidean norm of the vector e i at time τ , the largest Lyapunov exponent is given by The linear system without nonlinearities will give an exponential growth (relaxation) of the solution for a positive (negative) λ 1 . Therefore, the perturbations will need to be renormalized from time to time to prevent overflow (underflow), because they are only represented with a finite floating point numbers in the computer. To avoid this we use the algorithm of Benettin et al. 33 and renormalize our perturbation vector after each half period T d /2 of the ac voltage bias. Then λ 1 will be given by where d j is the vector perturbation growth during the jth renormalization period. Fig. 11 shows the results for λ 1 (t k ) obtained after computing of k renormalization intervals. In the limit k → ∞ λ 1 converges to its limiting value (positive, zero, or negative) according to the control parameter a. Notice a little difference between a=0.01 and a=0.0102 gives a drastic difference in values of the Lyapunov exponents from negative to positive.
To look at our chaotic solutions from the point of view of their dimensionality we have to compute the next Lyapunov numbers. Of course, we don't need all the Lyapunov spectrum of 40 numbers to characterize our system. Kaplan and Yorke defined a quantity called the Lyapunov dimension D L , given by a formula 34 where K is the largest integer such that K i=1 λ i ≥ 0 (all λ i are arranged in decreasing sequence). Since different Lyapunov exponents characterize the stretching and contracting of phase space in different directions, the Lyapunov dimension is then the number of vectors in phase space needed to describe an infinitesimal volume that remains constant on the average. D L is related to the information dimension of the system D 1 34 and can be used to characterize the fractal dimension of the associated chaotic attractor.
To calculate the first n exponents λ 1 , . . . , λ n we define n initial perturbation vectorsê (n) i (0), which are orthonormal, so that they form n-dimensional sphere. As the system is integrated forward in time, the sphere of trajectories will become an ellipsoid since different directions in phase space will expand or contract at different rates. The Lyapunov exponents λ i are determined by the rate of expansion or contraction of the principal axes of the ellipsoid averaged over the entire attractor. One additional problem arises here. Since all the vectorŝ e (n) i (τ ) tend to line up in the direction of the greatest growth we have to use the Gram-Schmidt orthonormalization procedure to separate the vectors into orthogonal components. The Lyapunov exponents λ i are then calculated by the same formula (A3) as for λ 1 , with d (n) j being the orthonormalized perturbation growth of ê (n) i during jth orthonormalization period. Fig. 12 shows the temporal convergence of the first three Lyapunov exponents for the particular chaotic solution (a=0.09). The estimated from those values Lyapunov dimension D L ≈ 2.208 indicates low-dimensional chaos for chaotic dynamics in our SL. By further increasing the length of the SL N the characteristic Lyapunov dimension does not grow. Since we did not see a transition from low-to high-dimensional chaos, it could be concluded that our system is not extensively chaotic. We have not observed hyperchaos either, where the second Lyapunov exponent is positive. Temporal convergence of the first three Lyapunov exponents for chaotic solution at a=0.09. Lyapunov exponents are scaled on the driving period T d .