{'fraction_non_alphanumeric': 0.03646424517593644, 'fraction_numerical': 0.0093288876276958, 'mean_word_length': 4.537811824788843, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Chaotic motion of space charge wavefronts in semiconductors under time-independent voltage bias

A standard drift-diffusion model of space charge wave propagation in semiconductors has been studied numerically and analytically under dc voltage bias. For sufficiently long samples, appropriate contact resistivity and applied voltage - such that the sample is biased in a regime of negative differential resistance - we find chaos in the propagation of nonlinear fronts (charge monopoles of alternating sign) of electric field. The chaos is always low-dimensional, but has a complex spatial structure; this behavior can be interpreted using a finite dimensional asymptotic model in which the front (charge monopole) positions and the electrical current are the only dynamical variables.


I. INTRODUCTION
The dynamics of propagating solitary waves (pulses) and monotone fronts in nonlinear partial differential equations (pde's) have been the focus of extensive research. Such problems are of interest in a wide range of fields including biology (population dynamics) [1], chemical reactions and combustion [2,3], plasma physics [4], and semiconductor electronic transport [5]. A common problem is to understand and predict the form and speed of the waves, as well as the possibility of multiple excitations in a sample of finite or infinite extent, and their interactions. We focus here on a model of electrical conduction in extrinsic semiconductors (involving time and only one spatial dimension) which exhibits negative differential resistance (NDR) and moving domains of high electric field. The model is specifically relevant to experiments on cooled bulk p-type Ge under voltage bias conditions, [6,7], but much of the observed qualitative behavior applies to a broad class of semiconductor systems with space charge instabilities. Phenomena observed for the p-Ge system include time-periodic oscillation of the current in a purely resistive external circuit under dc voltage bias due to the periodic creation of a solitary wave at the injecting contact, its motion inside the semiconductor and its annihilation at the receiving contact [7]. There is some similarity to the Gunn effect in n-GaAs, except that: (i) the local current density versus field characteristics [see j(E) below] in p-Ge presents an increasing third branch after the NDR region, and (ii) the solitary waves in p-Ge move much more slowly than the carrier drift velocity (the case in the usual Gunn effect) due to the generation-recombination dynamics of ionized traps which dominate the transport properties. Other experimental observations include intermittency near the onset of the oscillatory instability [8,9], and "spatiotemporal" chaos under combined dc and ac voltage bias [10]. Another important feature is the integral constraint which corresponds to voltage bias applied across sample contacts in semiconductor problems. Integral constraints also occur in other situations, for example, expressing mass conservation in problems of phase separation in binary mixtures and in certain biological problems [1,11,12].
Many of these phenomena have been successfully explained by means of a drift-diffusion model which includes impurity trapping of mobile holes and impact ionization of neutral acceptors [6,13]. Although much work has been done on this model problem (see [14] and references therein), important basic questions concerning its asymptotic description and chaos under time independent voltage bias are still open. In this paper, we present numerical simulation results which show chaos under dc voltage bias associated with multiple shedding of wavefronts. Multiple shedding of wavefronts occurs for appropriate values of contact resistivity, and was recently predicted on the basis of asymptotic calculations [15]. Here we introduce a finite dimensional model which provides a simplified description of space charge wave dynamics in long samples. This model uses relevant information from the asymptotics in Ref. [15], although we do not rigorously derive it from such asymptotic calculations. Nevertheless, solutions of the simplified model are in good agreement with the results of direct numerical simulations.
The paper is organized as follows. In Section II we present a drift-diffusion model which accurately describes the propagation of space charge wavefronts in extrinsic semiconductors under dc voltage bias such as p-Ge. In Section III we present numerical simulation results of the drift-diffusion model which indicate that wavefronts may propagate chaotically due to multiple shedding of wavefronts for appropriate bias in sufficiently long samples. Section IV presents the asymptotic model and numerical simulations thereof which are then compared with results from the drift-diffusion model. Conclusions are finally presented in Section V.

II. THE REDUCED DRIFT-DIFFUSION MODEL
In dimensionless form, the drift-diffusion model equations for a sample of length L can be written in a form [16,17]: The first equation describes the spatio-temporal evolution of the electric field E(x, t) inside the sample, where J(t) is the total current density. The transport coefficients V , K and R are, respectively, the average (drift) velocity of microscopic charge carriers (holes in the case of p-Ge), and coefficients describing the creation of free carriers via impurity impact ionization and the destruction of free carriers by capture onto an available impurity trapping site (i.e., neutral acceptor). All are nonlinear functions of electric field, and their forms and plots have been discussed extensively in the literature [5,13,14]. In this paper we use the same forms as in Ref. [14]. Equation (2) is a global constraint which expresses the voltage bias condition, and Eq. (3) is a boundary condition which represents the Ohmic injecting contact at x = 0 with contact resistivity ρ > 0. We refer to Eqs. (1) -(3) as the reduced driftdiffusion model because they are derived from a full drift-diffusion model by systematically dropping terms that correspond to short length and time scale processes of diffusion and displacement current, respectively. For a precise derivation as well as a complete table of conversion factors to dimensional units, see [5] and [16]. Some of them are: time 2.1 10 −3 ms, length 0.01 mm, electric field 10 V/cm, density current 128.16 mA/cm 2 , cross-sectional area 0.16 cm 2 . We have adopted the same symbol for both dimensional and non-dimensional variables.
The qualitative nature of much of the dynamical behavior found in the reduced driftdiffusion model-e.g., the instability of the stationary electric-field profile and propagating high-field domains [17][18][19]-depends only on the presence of a region of negative slope of the homogeneous stationary current density [14,17,18] over an interval of positive fields, and not on the exact form of the underlying coefficients. This is particularly true when the sample is closely compensated (the ratio of the acceptor concentration to the donor concentration, α, is only slightly larger than 1). Then j(E) is N-shaped for large enough positive fields: there is an interval (E M , E m ) between the abcissas of the maximum [j(E M ) = j M ] and the minimum [j(E m ) = j m > 0] current density for which dj/dE < 0 and j(E) > 0, as shown in Fig. 1. Also shown in Fig. 1 is the injecting contact characteristic which plays a crucial role in determining when new fronts are injected into the sample. The critical current density, J c , corresponds to the intersection of the contact characteristic with the homogeneous stationary current density. The role of J c has been elucidated using a rigorous asymptotic analysis of the system Eqs.

III. NUMERICAL RESULTS
To solve the system of equations (1) -(3) for E and J, we discretize the equations using finite difference approximations to the derivatives and employ an implicit method to generate the solution. The initial condition for the electric field is spatially uniform with a value that is consistent with the global constraint, Eq. (2). In Fig. 2 we show a spacetime plot of the electric field and associated current density for φ = 6.25 V/cm, just above the threshold voltage value for which propagating domain behavior occurs. The gray scale ranges from 5.3 V/cm (black) to 14.2 V/cm (white), and this scale is used in all similar plots that follow. The dimensionless sample length is 3800 corresponding to a real p-Ge sample of length 3.87 cm, and the contact resistivity ρ is 780 Ωcm corresponding to a value of 10.0 in dimensionless units. This case corresponds well to experimental data, but published data in p-Ge were only presented for one sample and relatively low bias values and contact resistivity. We clearly see that a single domain moves across the sample at constant speed, until it reaches the receiving contact. As it disappears, a new wave is created at the injecting contact and the process repeats periodically. The current versus time plot indicates that the current is steady when the domain moves in the sample interior, while there is an increase when the domain reaches the receiving contact. It is important to note that the fronts of changing electric field (or equivalently, the regions of nonzero charge density) are sharp in space relative to other physical length scales for this problem, i.e, the extent of the flat top domains or the sample length. We have found that this separation of length scales increases with sample length and holds for most biases of interest. It is only for voltages near the onset point that one tends to observe rounded solitary waves rather than well-separated pairs of fronts; this is the dominant space charge wave structure observed in shorter samples and has been extensively reported [14,17,18].
As the bias increases, the propagating domain becomes fatter and eventually a second small domain is nucleated and propagates part way into the sample; but it dies before reaching the receiving contact or merging with the larger domain. At even larger bias values the second domain merges with the primary domain near the receiving contact and this situation is shown in Fig. 3 which corresponds to φ = 7.25 V/cm. Again the current is plotted on the far right of the figure. When the first domain reaches the receiving contact the current increases. Instead of immediately starting the nucleation of a new wave, the area lost by the dying wave is gained by the trailing wave -note that the width of the trailing wave increases after the leading wave starts to disappear. The current increases, reaching a local maximum just before the trailing domain touches the leading domain, that is, the fronts collide. The current increases abruptly after the front collision and it rises to a global maximum at which point a new domain begins to nucleate at the injecting contact. As the domain forms, the current decreases and reaches a minimum at which point a new domain detaches and begins to propagate. Then the current increases until a second smaller domain is nucleated. Finally, the current settles to a rather low constant level as the two domains move steadily and in unison across the interior region of the sample. Current behavior is apparently dominated by the major events involving the fronts: collisions with the contacts or with each other. This suggests the viability of a dynamical model that focuses on discrete front motions and the current J(t).
At larger biases, the portion of the sample occupied by the high field value E 3 is larger, reducing the separation between domains, and giving more complicated E(x, t) structure and J − t behavior. In Fig. 4 we show a space-time plot and current for what appears to be a chaotic state for an applied bias of φ = 10.0 V/cm. The spatiotemporal dynamics possess a great deal of structure and complexity. The process of multiple domain shedding is similar to that for the previous case. The large current peaks correspond to nucleation of leading domains. The leading domains cross the sample without catching up or forwardcolliding with any other high field regions, and are indicated by dark regions that extend all the way across the space-time diagram. Note also the larger spatial extent of the leading domains than in Figs. 2 and 3. The aperiodicity of the current is reflected in the irregular appearance of the maximum current peaks or, equivalently, of the dark strips that extend across the entire sample. In between them are a number of local maxima corresponding to the shedding of trailing domains. Figure 5 (a) shows a bifurcation diagram in which we plot all values of successive current maxima as a function of φ. An important feature in this diagram is the apparent presence of windows of chaotic behavior with a large number of points being visited. For periodic states we see a small number of points corresponding to perfectly repeating current maxima. Also, in the periodic regimes there are points where various branches merge or disappear and these correspond to the development or destruction of trailing domains. For the parameter values selected here we do not observe period doubling. We conjecture that the route to chaos here is of boundary crisis type in which the attractor collides with a periodic orbit on its basin boundary [20].
In Fig. 5 (b) we show the largest Lyapunov exponent λ 1 versus φ for the reduced driftdiffusion model. This unambiguously confirms the presence of chaos in the "chaotic" windows. The next two exponents have been calculated and are never positive, so that the chaos we see is of a low-dimensional variety. To compute the exponents, we used an algorithm outlined in Ref. [21] adapted for use with partial differential equations and using adaptive control of the integration time step [22]. The values of λ 1 are zero in the periodic regimes as they should be for periodic behavior. The smallness of λ 1 in the chaotic regime is easily understood by recalling the period of the system, about 1000 non-dimensional time units. This indicates that the chaos originates in processes that occur on time scales on the order of the front transit time across the sample. To our knowledge, this is the first time that chaos due to multiple shedding of wavefronts has been observed in a drift-diffusion model of this type. This type of chaotic behavior has not been reported for experiments on p-Ge with time-independent voltage bias, most likely because experimentally studied samples were too short, biases were not sufficiently large, or contact resistivity was too low. However, in early experimental studies of the Gunn effect in GaAs, J. B. Gunn [23] observed that for long samples current oscillations were almost completely random, resembling white noise. He also found that short samples produce aperiodic oscillations when circuit impedance is sufficiently large. It is plausible that Gunn may have observed a similar form of chaos to the one we have found numerically. In the next section we develop an asymptotic model in which the chaotic dynamics is understood to arise from the aperiodic nucleation of fronts at the injecting contact.

IV. ASYMPTOTIC MODEL
The asymptotic model used in this paper consists of describing the evolution of the current when all wavefronts are detached from the injecting contact by an appropriate ordinary differential equation (ode) for J, tracking the position of the wavefronts and proposing a simplified mechanism for creation and destruction of wavefronts. While our new model is compatible with the asymptotic calculations of Ref. [15], we have not rigorously derived it from these calculations. Instead, we have proposed a simplified dynamics to account for our numerical observations motivated by asymptotic results.
We begin by assuming that E M < φ < E m in Fig. 1 and that E M /j M < ρ < E m /j m . Then a Gunn effect mediated by solitary waves occurs [24], as shown by the numerical simulation in Fig. 2. For appropriate parameter values, there appears a regular oscillation of the current caused by repeated creation, motion and destruction of high-field domains in the sample. High-field domains are formed by two wavefronts separating a region where the electric field is uniform and large from regions of uniform low field. Clearly, there are positively and negatively charged wavefronts, having ∂E/∂x > 0 or ∂E/∂x < 0, respectively. Near the contacts, there are narrow boundary layer regions where the electric field changes abruptly. Creation of high-field domains occurs at the injecting contact, via an instability of the boundary layer which expels a high-field domain from the injecting contact to the bulk of the sample. Typically, the total current changes most during wavefront creation and destruction events. In the limit as L → ∞, space and time scales are x/L and t/L, respectively [15]. Then j(E) = J, except in wavefronts and boundary layers at the contacts. If the field profile consists of a single high-field domain detached from the contacts, we have E = E 3 (J) inside the domain and E = E 1 (J), outside, where E 1 < E 2 < E 3 are the three zeros of j(E) − J for j m < J < j M , [15]. High and low-field regions are joined by wavefronts which are the unique solution of Eq. (1) (with zero right hand side) in the moving coordinate χ = x − X ± (t), dX ± /dt = c ± (J) (the signs + or -refer to the charge inside the wavefront) and appropriate boundary conditions. For example, at a positively charged wavefront, E → E 1 (J) as χ → −∞ and E → E 3 (J) as χ → +∞. The numerically determined values of c ± (J) are shown in Fig. 6. Boundary layers obey (most of times) a quasistationary version of Eq. (1) with appropriate boundary conditions on a semiinfinite spatial support. The instantaneous value of the current J(t/L) determines the field profile in the low and high uniform-field regions and the velocity of the wavefronts.
Next, assume that we have an initial field profile consisting of N high-field domains (solitary waves), each formed by two wavefronts located at X where t (i) b,± denotes the time at which the i-th monopole (with positive or negative charge) was born at x = 0.
The evolution of the total current density is determined by the bias condition (2), which may be approximately evaluated as (terms of order 1/L and smaller have been ignored here; note that the X (i) ± /L are of order unity). We can get an ode for J by differentiating Eq. (5) with respect to time and then substituting dX where i goes from 1 to N. The quantities n + and n − are, respectively, the number of positive and negative monopoles detached from the contacts (i.e., excluding possible monopoles at x = 0 and x = L), while j ′ 1 and j ′ 3 denote the derivative of the static j(E) characteristic with respect to electric field, evaluated at E 1 and E 3 , respectively. Notice that the system of equations (4) to (7) completely specifies the behavior of current and field profile on the scales x/L and t/L, except that we do not have conditions for determining when new fronts are emitted from the injecting contact.
We start with the simple case of Fig. 2: the motion of a single high-field domain far from the contacts. J satisfies Eq.  Figure 6, c + (J) [resp. c − (J)] is a decreasing (resp. increasing) function of the current. Therefore, J evolves exponentially fast toward the zero of the right hand side of this equation, J = J * . When the leading wavefront at x = X − arrives at x = L, it disappears almost instantaneously in the scale s, and we obtain dJ/ds = A(J) c + (J) > 0, so that the current increases. The injecting boundary layer near x = 0 ceases to be quasistationary when J surpasses the value J c at which the line J = E/ρ intersects the second (decreasing) branch of J = j(E). The precise description of the instability, which results in expelling a narrow high-field domain from x = 0 to the interior of the sample, can be found in [15]. It is enough to say that a certain semiinfinite problem has to be solved numerically and matched to the resulting situation with a narrow high-field domain (consisting of a region of E = E 3 (J) bounded by positively and negatively charged wavefronts) near x = 0 and a high-field region from a positively charged wavefront to x = L. In the new situation, we have dJ/ds = A(J) [2 c + (J) − c − (J)], and J tries to go toward the zero J = J + of 2 c + (J) − c − (J). Depending on the resistivity ρ, J c can be larger than J + (and then J decreases toward J + ), or J c ∈ (J * , J + ) (and then J starts increasing, and a second high-field domain may be expelled from x = 0). The simplest case, J c > J + , was described asymptotically in [15]. Provided φ is large enough, J evolves exponentially fast towards J + . When the old domain leaves the sample, only two wavefronts (bounding the new high-field domain) remain, and J evolves exponentially fast toward J * , so that a period of the oscillation is completed; see Fig. 2. The second case, J c ∈ (J * , J + ) is more complicated: numerical simulations show that multiple high-field domains may coexist in the sample at the same time as in Fig. 3.
To achieve a simplified description of the current oscillation, valid for any positive value of J c , we proceed to examine further the shedding process. The simplest rule to determine when a new wavefront is shed from the injecting contact would be as follows: a positive (negative) front is emitted at the instant that J passes through J c with positive (negative) time derivative. However, this rule neglects the time needed for sufficient charge to be injected at the contact to form a propagating front. We may estimate the effective delay time by considering the time evolution of Eq. (1) evaluated at the injecting contact, where u(t) = ∂E(0, t)/∂x and the argument of V , K and R is E(0, t) = ρJ(t), i.e., the value of electric field at the injecting contact. In this equation, we can think of J(t) as driving charge injection processes which determine front launching. Based on extensive simulations of the reduced model, we have found that u(t) must attain a sufficiently large positive or negative value for the front to detach and begin to propagate. This value is generally found to lie between 50% and 90% of the steady state value that |u| would have in the case of no propagating fronts, i.e., where E near the injecting contact rapidly rises to the E 3 value or rapidly falls to the E 1 value. The asymptotic model system is fully defined once the threshold is set and consists in Eqs. (4) to (8). Thus, in the limit of an infinitely long sample, terms of order ǫ = 1/L drop and we arrive at a low-dimensional dynamical system, which consists essentially of: (i) propagating negative and positive charge points that move according to the Eq. (4), (ii) subject to the conservation law Eq. (7), (iii) which produce a measurable current according to Eq. (6), and are created according to Eq. (8). Note that in the ǫ → 0 limit the current will exhibit slope discontinuties at the formation times t (i) b,± and destruction times t (i) d , but will be otherwise continuous and governed by Eq. (6). We note that the rigorous foundations of this and similar asymptotic models have been explored recently using singular perturbation methods [15,25].
In this paper, we estimate the order-of-magnitude of the time delay associated with wavefront formation by evaluating Eq. (8) for J ≈ J c . This implies a relaxation time of i.e., the approximate time for u(t) to go from a value of u(0) = −1 to u(τ ) = 0. By dJ c /dt, we refer to the value of dJ/dt when J crosses J c . Then, we adopt the criterion that a new front is born at x = 0 at the time t + aτ where t is the time at which J = J c , and aτ is a delay time. Here a is a number of order one which is determined from simulation of the reduced model for a particular bias voltage and then assumed to apply over the complete range. Wavefront destruction is assumed occur instantaneously at times t We ignore the finite duration of (fast) monopole destruction stages which is equivalent to the well-justified approximation of neglecting the diffusive boundary layer at the receiving contact [15]. It should also be kept in mind that the index instantaneously decreases by one when wavefronts downstream collide with one another.
Let us now use the above asymptotic model to interpret the simulation results for Eqs. (1) -(3). The case of contact resistivity such that J c > J + has been explained already: we obtain the usual Gunn effect with at most one solitary wave detached from the contacts for any time [15,25]; see Fig. 2. Let us assume now that the contact resistivity is such that J c ∈ (J * , J + ). Then the current will increase after creation of a solitary wave, because 2c + − c − > 0 and multiple wave shedding is possible [15]. This situation is shown in Fig.  7, which shows simulation results of our simplified asymptotic model for similar parameters to those of Figure 3. The latter is depicted using data from direct numerical simulation of the reduced pde model Eqs. (1) -(3). To obtain Fig. 7, the values of a were set to 13.75 for positive front emission and 7.82 for negative front emission. We use these same values in the data of Fig. 8, which has same bias values as Fig. 4. The chaos appears to be closely tied to the asynchronous emission of fronts. This explains why the chaos observed for this partial differential equation system is low-dimensional. It is interesting to speculate how the maximal number of domains possible might scale with system size and contact resistivity.

V. CONCLUSION
We have utilized asymptotic analysis of a pde model (which describes the trap-dominated slow Gunn effect in a long sample) to explain the dynamics of space charge waves and current vs. time, including low-dimensional chaos which is nonetheless accompanied by spatially complex structure suggesting a loss of spatial coherence. The building blocks of this analysis are the heteroclinic orbits used to construct the typical solitary waves mediating Gunn-like oscillations. During most of the oscillation, the motion of the heteroclinic orbits and the change of the electric field inside and outside the solitary waves (enclosed by heteroclinic orbits) follow adiabatically the evolution of the total current density. When a solitary wave reaches the receiving contact, the current increases abruptly and the asymptotic model adequately approximates this as instantaneous. As an outcome, we have found a criterion that shows that single or multiple wave shedding is possible during each oscillation, depending on the resistivity of the injecting contact. While single shedding is the usual (stable) Gunn effect, multiple wave shedding may break the spatial coherence of the electric field within the sample. This new instability mechanism provides an explanation for the complicated behavior observed in experiments performed in long semiconductor samples [23,9] and in the numerical simulation of the drift-diffusion model. We have confirmed these results by direct numerical simulation of the reduced model, in particular the new predictions of multiple shedding of solitary waves in the unstable case. Although simulations and analyses have been presented here for the p-Ge model, the general approach is quite general and applies to a wide class of pde models which possess the following common properties: 1) an integral (over space) constraint; 2) standard boundary conditions which permit multiple stationary states -i.e., negative differential resistance; and 3) solitary waves (i.e., pulses) and fronts. It is interesting to speculate that different models may lead to the same class of long-sample asymptotic limiting model which completely determines and explains the long time dynam-ics (including chaotic temporal behavior and loss of spatial coherence) of the respective full models.
It is a pleasure to acknowledge beneficial conversations with P. J. Hernando and M. Kindelan. We acknowledge financial support from the Spanish DGES through grant PB98-0142-C04-01, and one of us (SWT) acknowledges support of the Fulbright Foundation.