Asymptotic analysis of the Gunn effect with realistic boundary conditions

A general asymptotic analysis of the Gunn effect in n-GaAs under general boundary conditions for metal-semiconductor contacts is presented. Depending on the parameter values in the boundary condition of the injecting contact, different types of waves mediate the Gunn effect. The periodic current oscillation typical of the Gunn effect may be caused by moving charge-monopole accumulation or depletion layers, or by low or high-field charge-dipole solitary waves. A new instability caused by multiple shedding of (low field) dipole waves is found. In all cases the shape of the current oscillation is described in detail: we show the direct relationship between its major features (maxima, minima, plateau's, ...) and several critical currents (which depend on the values of the contact parameters). Our results open the possibility of measuring contact parameters from the analysis of the shape of the current oscillation.


I. INTRODUCTION
The Gunn effect appears in many semiconductor samples presenting negative differential resistance and subject to voltage bias conditions [1][2][3][4][5][6][7]. It consists of a periodic shedding of pulses of the electric field at the injecting contact, which then progress and are annihilated at the receiving contact. As a result there appears a periodic oscillation of the current through a passive external circuit attached to the semiconductor. Under different conditions, the current self-sustained oscillation may be caused by the motion of charge accumulation layers (charge monopoles) [2], not by the usual electric field pulses which are charge dipoles. Most of the experiments on the Gunn effect in different materials take place in samples with attached planar contacts, so that the wave motion may be safely assumed to be one-dimensional. Despite the vast literature on the Gunn effect, it is surprising that many basic questions remain poorly understood. For example, given a description of the charge transport in the bulk semiconductor (say at the level of drift-diffusion and rate equations), which are the proper boundary conditions for given contacts and how they affect the self-sustained current oscillations. The first question has been addressed in a companion paper, Ref. [8], while the second will be answered here.
Until recently, when confronted with the Gunn effect, theorists would resort to computer simulations of more or less complicated models (which were supposed to reflect the physics of a given semiconductor), and would then explain qualitatively their numerics. Resort to special solutions valid for infinite semiconductors at constant current bias conditions [3], or to extrapolations of Kroemer's NL criterion [9] were often used to interpret the simulation results. This left the processes of generation and annihilation of domains at the contacts (and in fact it also left out the dynamics of wavefronts and pulses) outside theoretical considerations. Concerning asymptotic descriptions of the Gunn effect which delve deeper than just numerical simulations of drift-diffusion models, some progress has been made recently [10][11][12][13]. These works propose asymptotic descriptions of the Gunn effect exploiting the fact that this effect is seen most clearly in semiconductors having a large value of the product of sample length times doping (basically a dimensionless length). The role of the "N-L product" in the analysis of the Gunn instability was already discovered by Kroemer [2] and exploited to study the linear stability (small signal analysis) of stationary solutions by many authors [3]. It was recognized only much later that in the limit of large dimensionless length (N-L product) it is possible to describe asymptotically both the onset [16] and the fully developed Gunn instability [10]. In this asymp-totic limit the processes of repeatedly generating a new wave (a charge monopole or dipole domain) at the injecting contact, the motion of the wave and its annihilation at the receiving contact may be well separated. Then they can be analyzed and combined to fully describe the Gunn effect. In particular, the effect of contacts on these processes and in determining the shape of the current oscillation can be clearly stated. In this paper we use our asymptotic theory to study Kroemer's model for n-GaAs under b.c. corresponding to ideal MS contacts. We find that these b.c. give rise to a multivalued control current-field characteristic at the injecting contact. The asymptotic analysis shows that the Gunn effect can be mediated by both charge monopole or dipole domains according to the values of the contact parameters. Shedding of new charge dipole waves from the injecting contact is adiabatic in clear distinction with what happens if the control characteristic of the contact is single-valued [12,13]. In the later case (analyzed for a p-Ge model in Ref. [12]) the charge dipole pulses are created very rapidly at the injecting contact and they advance and grow simultaneously (see also [10]), whereas for multivalued control characteristic the boundary layer at the injecting contact grows adiabatically to a much greater size before a new pulse can be shed. These facts may determine appreciably the shape of the current oscillations. Our analysis could be extended to more complex models also displaying the Gunn effect, [12][13][14]. Depending on the values of the parameters characterizing the injecting contact and of the dc voltage bias, we find Gunn oscillations mediated by charge accumulation and depletion monopole wavefronts, and high and low field charge dipole domains. We also find narrow regions in the parameter space where multiple shedding of dipole domains at the injecting contact occurs. We have thus found a characterization of all possible dynamic behaviors which an ideal metal-semiconductor contact would give rise to in the Kroemer model for the Gunn effect. This opens the possibility of extracting information about the contacts from the analysis of the Gunn oscillations themselves, a subject of considerable interest for applied researchers.
The rest of the paper is as follows. In Section II we present Kroemer's model and the boundary conditions for metal-semiconductor contacts discussed in the companion paper [8]. In Section III we present our asymptotic analysis of Kroemer's model and find that different types of Gunn effect are possible according to the values of the bias and of certain dimensionless parameters appearing in the b.c., i 0 and α 0 : charge monopoles (moving charge depletion and accumulation layers), high and low field solitary waves (charge dipoles), multiple (low field) charge dipoles, are predicted and confirmed by numerical simulations. Section V contains our conclusions whereas the Appendices are devoted to different technical matters related to the main text. The unipolar drift-diffusion model for the Gunn effect proposed by Kroemer [2,15], is generally accepted to provide a rather complete description of the main features of this effect. In the dimensionless units described in Ref. [8], Kroemer's model is Eq. (1) is Ampère's law which says that the sum of displacement current and drift-diffusion current is equal to the total current density J(t). It can be obtained by differentiating the Poisson equation, ∂E/∂x = n − 1, with respect to time, substituting the charge continuity equation ∂n/∂t+∂j(x, t)/∂x = 0 [the electron current density is of the drift-diffusion type: j(x, t) = nv(E) − δ ∂n/∂x], and then integrating the result with respect to x. The electron velocity is assumed to be N-shaped. For specific numerical calculations we shall use Kroemer's curve [15] v (it has a maximum v M > 0 at E M > 0 followed by a minimum 0 < v m < v M at E m > E M ) and the electron difusivity δ to be constant. The dc bias φ is the average electric field on the semiconductor sample. Eqs. (1)- (2) need to be solved with an appropriate initial field profile E(x, 0) and subject to the corresponding b.c. For ideal metal-semiconductor the following mixed boundary conditions have been derived in [8] where i i and α i , (i = 0, L), are dimensionless parameters which are a combination of the semiconductor effective density of states, contact barrier height, Richardson's constant, doping and temperature (see [8]). In what follows, i i will be assumed to be positive because the physically interesting phenomena (including the usual Gunn effect mediated by high field domains) are observed for these values of i i (see phase diagram in [8]).
For typical n-GaAs data, δ ≪ 1 and L ≫ 1 [10]. In this limit, we shall find approximate solutions to the initial boundary value problem (1)-(5) for E(x, t) and J(t). Strictly speaking, the simple asymptotic description that follows holds in the limit L → ∞, even when δ = O(1) [13]. Assuming δ ≪ 1 just simplifies the description of the traveling waves of electric field in the semiconductor through the use of characteristic equations and shock waves [10,17,18]. For example, it is shown in Appendix A (by using the method of characteristics) that the boundary condition (4) implies that the electric field at the injecting contact, E(0, t) = E 0 (t), obeys the following equation: These expressions constitute a Dirichlet boundary condition for the electric field which contains the same information as the mixed condition (4). The contact curve, j c (E), presents two extrema, a minimum, j cm = j c (E m ), and a maximum, j cM = j c (E M ), at the same field values as the electron velocity curve v(E). j c (E) tends to j sat 0 = α −1 0 + i 0 for high electric fields. As we shall see below, during most of an oscillation period, j c (E 0 ) ∼ J, and this expression yields a multivalued contact-characteristic curve relating the field at the injecting contact to the actual value of the current density.
To take advantage of the large-length limit, we will use the following rescaled time and length, Then Eqs. (1)-(2) become 1 0 E(y, s) dy = φ.
Notice that the boundary condition (6) becomes in the limit δ ≪ 1.

III. ASYMPTOTICS OF THE GUNN EFFECT
In a previous paper [8] we have analyzed the stationary solutions of Kroemer's model with metal-semiconductor contacts and discussed their stability. No stable stationary solution is expected for certain ranges of bias and the cathode contact parameters i 0 and α 0 . In these circumstances, the Gunn effect mediated by either moving charge monopoles or dipoles might appear. A rich phenomenology of propagating waves and current oscillations has been numerically observed for these parameter ranges. Among them, we have observed both high (Figs. 1, 2) and low (Fig. 3) field solitary waves (moving charge dipoles), multiple low field dipoles (Fig. 4), moving charge accumulation (Fig. 5) and charge depletion monopoles (Fig. 6). In Ref. citecompanion, we have identified the critical currents determining which type of waves mediate the current oscillation (in the limit δ ≪ 1). They are related to the boundary conditions in the following way: • If v m < j cm < j cM < v M , moving charge dipoles mediate the Gunn effect (see Ref. [8] for more details).
As shown in the figures, there are several stages in each period of the current oscillation corresponding to the processes of generation, propagation and annihilation of electric field domains. Each stage has its own time and space scales, and hence the oscillation can be suitably described by means of a matched asymptotic analysis. For instance, the annihilation of wavefronts takes place on a fast time scale compared to that governing wavefront propagation. Thus on the time scale of wavefront propagation, the annihilation of wavefronts is a quasiinstantaneous process during which the time derivative of the current density changes appreciably while the current itself, J(s), does not change. On the other hand, the generation of fronts takes place adiabatically on a much slower time scale comparable to wavefront propagation. In this stage both J(s) and its derivative change appreciably. This is quite different from the fast generation of fronts and pulses observed for other types of boundary conditions [12].
As long as these different processes take place on different time and space scales, there are different stages of the oscillation which can be analyzed separately. This happens for certain bias ranges. We shall then construct the approximate electric field and current density solutions by means of matched asymptotic expansions. A detailed description of a period of the current oscillation will be then obtained. For other bias values, several processes occurs almost simultaneously (e.g., the annihilation of a wavefront may occur during the process of detachment of another wavefront from the injecting boundary layer for low bias values). This complicates the asymptotic description without adding much to our physical understanding, so that we will omit the details.
Note that in the limit ǫ = 1/L ≪ 1, the solutions of Eqs. (9)-(10) are piecewise constant: on most of the y-interval, E is equal to one or another of the zeros of v(E) − J [notice that this equation may have three zeroes, we denote by E 1 < E 2 < E 3 ], separated by transition layers that connect them. At y = 0 and y = 1 there are boundary layers (quasi-stationary most of the time), which we call injecting and receiving layers, respectively. It can be seen that the propagation of fronts turns out to be a quasi-stationary process while the generation is not. Thus, two different asymptotic approaches will be used: one for the description of the quasistationary propagation of fronts and the other for their (non-quasistationary) generation.

A. Quasistationary propagation of wavefronts
In this section we will present the asymptotics of the quasistationary propagation of wavefronts, and of the corresponding time evolution of the density current.
Wavefronts are moving transition layers connecting regions of the sample where the electric field is spatially uniform. In order to describe their quasistationary propagation we will proceed as follows. We assume that E(y, s) is either E 1 (J) or E 3 (J) outside boundary layers and wavefronts. Let the wavefront located at y = Y (s) move with velocity c = dY /ds. For each value of J, the wavefront advances with speed c = c + (J), if it connects To find the inner structure of a wavefront and its speed, we introduce a coordinate moving with the wavefront. Then we need to solve the following problem for the equation The solution of this problem will provide both the speed and inner structure of the wave front. In terms of the phase plane (12), the previous problem is equivalent to find c = c + (J) so that there is a heteroclinic orbit connecting the saddle point The functions c ± (J) for our model are depicted in Fig. 7. Note that they intersect when J = J * , given by the so-called equal-areas rule [4].
In the limit δ ≪ 1, we can obtain explicit expressions for c ± (J) and for the corresponding wavefronts as we now recall [10,18]. The wavefront (moving depletion layer) joining (E 3 (J), 0) and (E 1 (J), 0) may be approximated by the exact solution of (12) which holds for any value of δ plus two corner layers of [18] for an explicit calculation. The width of this wavefront on the ξ length scale is ( The other wavefronts (moving accumulation layers) can be constructed by matched asymptotic expansions in the limit δ ≪ 1 and their velocities c + (J) and shapes depend on whether J is larger or smaller than J * . These wavefronts are composed of a shock wave joining two field values E − and E + (at least one of them should be equal to E i , i = 1, 3) plus a tail region which moves rigidly with the same velocity as the shock [10]. The inner structure of the shock wave (for very small but not zero δ) can be a quite complicated triple-deck set of boundary layers [18]. Let V (E + , E − ) be the velocity of the shock wave given by the equal-areas rule [10,17,18] We now have [10]: whereas E − is calculated as a function of J by imposing the condition that the tail region to the left of the shock wave moves rigidly with it Solving simultaneously (15) and (16) To the left of the shock wave (in the tail region) the field satisfies the (approximate) boundary value problem whereas E + is calculated as a function of J by imposing the condition that the tail region to the right of the shock wave moves rigidly with it Solving simultaneously (15) and (18) (with E − = E 1 ), we find both E + and c + = v(E + ) < J as functions of J. To the right of the shock wave (in the tail region) the field satisfies the (approximate) boundary value problem The functions c ± (J) are depicted in Fig. 7.
Notice that the inner structure of the shock wave has width O( √ δ) on the ξ scale while the total width of the wavefront is O(1) on the ξ scale and O(ǫ) on the y scale. In the limit δ ≪ 1, the structure of the wavefronts is thus one-sided: the wavefront is a discontinuity preceded or followed by a tail region.
In conclusion, the propagation of a single wavefront can be described by giving the position, Y (s), and velocity, c ± (J(s)), of the front and the values of the electric field on its left and right hand sides, E i (J(s)). When more than a single wavefront are co-moving inside the sample, the same description applies to each of them. All of the magnitudes involved in this description depend on time s only through the instantaneous value of the current J(s). This fact, together with the fact that the voltage must remain constant all the time, can be used to derive a simple closed equation describing the time evolution of the current density during the propagation stages. As an example of this result, let us consider the case a single propagating dipole (two fronts) (see Fig. 1b.1). In this case, neglecting transition and boundary layers, we have: E(y, s) = E 1 (J(s)) for 0 < y < Y 1 (s) and Y 2 (s) < y < 1 and E(y, s) = E 3 (J(s)) for Y 1 > y > Y 2 (s). For the voltage we have: By using the fact that the voltage is fixed, we can obtain an equation for J(s) by differentiating the bias condition (20) with respect to s. By noting that , the following simple closed equation for the current is obtained where A similar procedure can be applied to derive the corresponding closed equation for the current when different types of propagating domains are present. In general we will have n + fronts moving at velocity c + and n − moving at c − , with |n + − n − | = 0, 1. For these cases the current evolves following the equation As before, this equation holds as long as the propagating profile consists of n + fronts propagating at c + and n − at c − . Two typical cases can be considered in analyzing Eq. (24): (i) either n + or n − are equal to zero and (ii) n ± = 0, both numbers are different from zero. In the first case, we have either These results describe the quasistationary propagation of fronts and the corresponding time evolution of the current density during these stages. We will use them to interpret the results of the numerical simulations.

B. Generation of fronts
As mentioned before, for our b.c. an appreciable part of a period of the current oscillation may be spent generating new fronts non-adiabatically. Fronts of dipole domains are generated at the cathode whereas monopole wavefronts appear somewhere in the middle of the sample. In what follows, we will focus on the description of dipole domains. Monopole wavefronts have been recently described in detail elsewhere [11].
A simple rule concerning generation of dipole domains can be formulated: when the current density, J(s), crosses the maximum [resp. minimum] of the contact curve, j cM [resp. j cm ], a front moving with speed c − [resp. c + ], starts being formed when its generation is compatible with the field value at the bulk after the injecting contact.
Let us now show why the previous rule holds. Suppose that the current reaches adiabatically in the slow time scale s one of the critical values mentioned above, let us say j cM . Then the field at x = 0, given by Eq. The field in the injecting layer, in turn, increases until a moving wavefront moving with speed c − is formed. To describe this process, we adapt ideas developed for the analysis of the trap-dominated Gunn effect in p-Ge [12] to the present situation. An important difference is that the sharp increase in the contact field E 0 helps creating the new wavefront but a new pulse is not shed in the fast time scale σ = (s − s 1 )/ǫ: after the wavefront is created, the contact field varies on the third branch of j c (E) and the current has to decrease slowly [according to (24), on the s time scale] until E 0 can jump back to values on the first branch of j c (E) when J = j cm .
To leading order, the field in the injecting layer solves the following semi-infinite problem whose derivation can be found in Appendix B: x > 0, −∞ < σ < +∞ , where E 0 solves (6) with J = j cM + ǫJ (1) (σ), and J (1) (σ) is given by where h(σ) is and formulas for the positive constants α, β and γ may be found below, in Eqs. (B10). The function h(σ) is the area lost due to the motion of the old front during the time σ minus the instantaneous excess area under the injecting layer minus the constant excess area under the old wavefront at Y = Y 1 (s 1 ) [that is, under the heteroclinic orbit connecting (E 1 (j cM ), 0) and (E 3 (j cM ), 0)].
As σ → −∞, we have to impose the following matching condition on an appropriate overlap domain: Here E stat (x; J(s)) is the quasistationary injecting layer solution of (12) with c = 0 such that E stat (0; J(s)) satisfies (4) and E stat (∞; J(s)) = E 1 (J(s)) for s < s 1 , J(s 1 ) = j cM . Notice that the term ǫJ (1) (σ) is needed so as to avoid that the solution of this problem stay indefinitely in the quasistationary field E stat (x; J(s)).
The solution of the previous semi-infinite problem reveals the growth of the field at the contact and inside the injecting layer until: (i) E 0 becomes E c3 (j cM ) [E c1 (J) < E c2 (J) < E c3 (J) are the three possible solutions of J = j c (E)], (ii) E(x, σ) increases to E 3 (j cM ) as x increases from x = 0 and then (iii) it has the structure of a wavefront connecting E 3 (j cM ) to E 1 (j cM ). This wavefront advances with velocity c − (j cM ). The formation of a front when J crosses j cm can be explained similarly.

C. Putting the pieces together
In the two previous subsections, we have presented the basic features of the asymptotic description of the Gunn effect, namely, quasistationary front propagation and the generation of new fronts. Now, we are in a position to put all the pieces together, and describe a full period of current oscillation. We shall distinguish three cases: high-field dipoles, low-field dipoles and monopoles.

Dynamics of high-filed dipoles
High-field dipole domains have been observed to appear when v m < j cm and J * < j cM < v M . Then depending on the value of j cm with respect to J * , and of j cM with respect to J † ≡ J 2,1 , different situations may occur.
Let start considering the case v m < j cm < J * < j cM < J † < v M . This situation corresponds to the propagation of high-field dipole domains as shown in Fig. 1. We will assume the initial electric field configuration to correspond to a single propagating high field domain ( Fig.  1b.1). This configuration corresponds to n + = 1 and n − = 1, and hence the current satisfies Eq. (22), evolving from an initial value J(0) ∈ (j cm , J * ), towards the fixed point J * = J 1,1 . After a certain time, the wavefront located at Y 2 reaches the end of the sample and disappears, producing an abrupt change in the time derivative of the current. We have a new stage with n + = 1 and n − = 0, Fig. 1b.2, governed by Eq. (25). Its solution increases until it surpasses the value j cM . At this point the injecting layer becomes unstable and it starts shedding a new front. The formation dynamics has been explained in detail in the previous section. Then a new slowly varying stage begins. There are two leftover wavefronts, the old one located at y = Y 1 (s) [which advances toward y = 1 with speed c + (J)], and a new one located at y = Y 4 (s), moving with speed c − (J), and leaving behind a quasistationary field E 3 (J), see Fig. 1b.3. Again, the field configuration corresponds to n + = n − = 1, with J(s) following (22), and decreasing exponentially fast towards J * . Before J reaches J * , the front located at Y 1 reaches the end of the sample and disappears, thereby producing a new abrupt change in the time derivative of the current density. Then only the recently formed front [located at Y 4 (s)] is present on the sample (Fig. 1b.4), which corresponds to n + = 0 and n − = 1. J(s) decreases according to Eq. (26) until the minimum of the contact curve, J = j cm , is reached. After that, a front moving with speed c + is generated. The charge dipole wave thus created evolves adiabatically, and the current density is described again by Eq. (22), Fig.1b.1. We have come back to the initial situation, and one period of the current oscillation has been completed.
It is worth noting that by means of this analysis some of the most relevant features of the current oscillations, as the maximum and minimum currents, J max and J min , have been identified with quantities related to the contact parameters, J max ≈ j cM and J min ≈ j cm , opening then the possibility of determining the values of contact parameters from the analysis of the current oscillations.
Other situations involving dipole domains have been numerically identified. They are described by the same type of analysis. In what follows only the more relevant features of these situations will be considered. Let us as-sume now v m < J * < j cm < j cM < J † < v M . This case corresponds again to high-field solitary waves (dipole domains, see Fig. 2), but the current oscillations have a different shape. The three first stages of the oscillation, Fig.  2b.1,2,3, correspond to the propagation of a single dipole, annihilation of a front and generation of new one. They are similar to the stages described above. Now, however after the new front has been formed, Fig. 2b.3, and the current is decreasing towards J * , J reaches the critical current, j cm (because J * < j cm ). A new front moving with speed c + is then created. After the formation process has finished, we have a configuration with two fronts moving at c + and one at c − , that is, with n + = 2 and n − = 1 (Fig. 2b.4). Hence the current satisfies [see Eq.
and it starts increasing again, trying to reach J 1,2 ≡ J † . Before this value may be attained, the old front located at Y 1 , arrives at the receiving contact and disappears. We obtain again Eqs. (22)-(23) and recover the initial situation. Thus a full period of the Gunn oscillation is again completed.
There are other possible situations for propagating high field domains, but we have not found them in our numerical simulations with the curve v(E) considered in the present work. For example in the second case we have described above, v m < J * < j cm < j cM < J † < v M , after the formation of the new front, Fig. 2b.4, the current density could have reached the value J † before the old front located at Y 1 had arrived at the receiving contact. In such situation, a new front moving at c − could be formed, giving rise to an electric field configuration with n + = n − = 2. In this configuration the current would decrease again to J 2,2 = J 1,1 = J * . We would then have a more complicated situation with two pulses and the dying wavefront simultaneously present in the sample. Proceeding in a similar way, we could therefore observe simultaneous coexistence of several pulses for appropriate ranges of parameters. This situation (multiple shedding of high field domains) has been numerically observed in other models [12,13].

Dynamics of low-field dipoles
Low-field dipole domains appear when v m < j cm < j cM < J * ( fig.). Depending on the value of the applied bias, single (high bias values) or multiple (low bias values) propagating low field domains are obtained [see Figs. 3,4].
Let us start considering the case of single propagating low field dipole domains, Fig. 3. We consider an initial field configuration having a low field domain far from the contacts. Then there is a wavefront located at Y 1 , moving with speed c − , and a wavefront located at Y 2 > Y 1 , moving with speed c + , that is, n + = n − = 1, Fig. 3b.1, and J(0) ∈ (j cm , J * ). The current increases towards J * according to Eq. (22), until Y 2 = 1. Then J starts satisfying (26) (corresponding to n + = 0 and n 1 = 1, see Fig.  3b.2), and it decreases until the value j cm is reached. After this occurs, a new front (moving at speed c + ) is created, while the old wavefront Y 1 < 1. We have a stage described by Eq. (22) with n + = n − = 1, see Fig. 3b.3, during which J increases past j cM (j cM < J * ). Then a new front moving with speed c − starts being created. At the same time, the old front located at Y 1 reaches the end of the sample, giving rise to a complex stage in which the non-stationary effects can not be neglected Fig. 3b.4. After that stage, we recover the first situation Fig. 3b.1, and a complete period of the oscillation has been described.
Let us now consider an example of multiple propagating low-field dipole domains (Fig. 4). As we mentioned above, they appear for small bias values. The asymptotic analysis of this case is more complicated because there are new stages (fusion of two wavefronts inside the sample) whose detailed description is outside the scope of this paper. Let us start with the configuration shown in Fig. 4b.1, for which two low field domains coexist. This configuration corresponds to n + = n − = 2. Following Eq. (24), the current will evolve according to thereby increasing towards J 2,2 = J * . Before this value can be reached, the front located at Y 4 reaches the end of the sample. The resulting configuration (Fig. 4b.2) has now n + = 1 and n − = 2, and the current decreases according to If the curve v(E) were such that the fixed point of Eq.
(34) existed, the current would tend to such a value. For our choice of v(E), this fixed point does not exist (see Fig. 7), and therefore the current decreases all the time. During this process j cm will be crossed. According to our previous considerations, a new front moving at speed c + should then be formed. However, this does not occur because melting of the fronts Y 2 and Y 3 starts after Y 4 reaches y = 1. This melting process seems to inhibit the creation of new wavefronts until it is completed, which happens for J < j cm . Then a new front moving at speed c + is rapidly created and we have a stage, with n + = n − = 1, Fig. 4b.3, described by Eq. (22). The current increases until J = j cM , at which time a front moving at speed c − appears at y = 0, Fig. 4b.4. Then the current decreases following (34) until J = j cm . Another front moving at speed c + is then formed. We have now n + = n − = 2, Fig. 4b.5 and J increases according to (22). When J > j cM , a front moving with speed c − is created at y = 0, so that n + = 2, n − = 3, see Fig. 4b. 6. Now the current decreases towards J 2,3 according to the equation This stage lasts until Y 1 = 1, at which time we are back at the first stage, Fig. 4b.1, having completed a full period of the oscillation. More complicated examples of multiple low field domains exist and they could be described similarly.

Dynamics of charge monopoles
When j cM > v M , or when j cm < v m and j sat > v m , the Gunn effect is mediated by wavefronts which are charge monopoles. The case j cM > v M , corresponds to moving accumulation layers, Fig. 5, while moving depletion layers are observed for j cm < v m and j sat > v m , Fig.6. H. Kroemer [15] discovered numerically the Gunn effect mediated by charge accumulation monopoles for boundary conditions corresponding to the control characteristics, while its asymptotic analysis (for B ≪ 1) was performed many years later [10]. Recently a Gunn effect mediated by charge accumulation monopoles has been used to describe self-sustained oscillations of the current in GaAs/AlAs superlattices [21] in the limit of weakly coupled, low doped, long superlattices [11].
The complete asymptotic analysis found in Ref. [11] can be used to describe the present situation after a few simple changes are made. First of all the equal-areas rule in the case of the superlattice refers to 1/v(E), not v(E) [22]. The second important change is that of the injecting boundary condition in the description of the birth of a new monopole: instead of a rigid Neumann boundary condition ∂E(0, t)/∂x = c we should use (4). These two changes can be implemented without difficulty and the details will be omitted. An important difference with the Gunn effect mediated by dipole solitary waves is that the amplitude of the current oscillations is larger (its largest value is approximately v M −v m for the charge accumulation monopoles and j cM − j cm for the dipoles). The monopoles "probe" the full region of negative slope of v(E) while the dipoles "probe" a smaller region.

IV. CONCLUSIONS
We have performed an asymptotic and numerical analysis of the Gunn effect in n-GaAs under general boundary conditions for metal-semiconductor contacts. We have shown that the Gunn effect is mediated by (i) moving depletion charge monopoles, (ii) moving accumulation charge monopoles, (iii) high-field dipole solitary waves or (iv) low-field dipole solitary waves, according to whether the critical contact currents j cM = j c (E M ), . Some of these results are well-known for boundary conditions given by Kroemer's control characteristic. In addition we have shown that there are new instability mechanisms consisting of multiple generation of (low-field) charge dipole solitary waves in the region near the injecting contact if j cM is close enough to J * , and the dimensionless length is large enough. In each case we have been able to describe in detail the shape of the current oscillation, identifying some of its main features (maxima, minima, plateaus, . . . ) with critical currents appearing in our model, among them the contact currents j cm and j cM . This result opens the possibility of determining contact parameters by means of the analysis of the shape of the current oscillation. Our results might be of use in the analysis of self-sustained current oscillations in weakly coupled n-doped superlattices [21], once the role of contacts and the boundary conditions they generate are understood in these systems. Beside this, our results are not restricted to the particular model of the Gunn effect studied here, but seem to hold for a general class of models, supporting the idea of the "universality" of the Gunn effect [13].

V. ACKNOWLEDGMENTS
We thank Dr. M. Bergmann for fruitful conversations and acknowledge financial support from the the Spanish DGICYT through grant PB94-0375 and PB95-0881 and from the EC Human Capital and Mobility Programme contract ERBCHRXCT930413. One of us (G.G.) acknowledge support by Generalitat de Catalunya.

APPENDIX A: NONSTATIONARY FIELD AT THE INJECTING CONTACT
We shall now prove that the field E 0 (t) at the contact x = 0 obeys the equation: This can be seen from the method of characteristics applied to (1) with δ = 0 and the boundary condition (4). The characteristic equations are to be solved with the conditions Clearly, The last term in this equation can be obtained from the boundary condition (4) as follows. From the solution of (A3) and (A5), we obtain ∂x ∂τ (τ ; τ ) = −v(E 0 (τ )). Then the boundary condition (4) yields: Insertion of (A7) into (A6) yields (A1) (with τ = t). To obtain the equations governing the shedding stages, we need to keep terms of order ǫ in our calculations, as indicated in Ref. [12]. The outer (bulk) expansion is We can obtain easily E > follows from (B3) and these two last equations. These equations should be solved with the matching conditions that E (0) i → E i (j cM ) (i = 1, 3) and J (0) → j cM as σ → −∞. These problems have the solutions J = j cM and E (0) i = E i (j cM ). Thus we expect that the current density does not depart substantially from j cM as ǫ → 0.
The O(ǫ) corrections E (1) i (σ) obey the equations: which immediately follow from (9). To find J (1) we proceed as follows. First of all, we write the bias condition (10) including terms of order ǫ in the approximation of the electric field and the current density: Here E(ξ) is the field inside the wavefront at Y 1 and E(x, σ) is the field in the injecting layer.
We now obtain the first order differential equation (29) for J (1) (σ) by applying the operator Θ 1 Θ 3 to both sides of (B9) and then using (B7): where the coefficients α, β and γ are All functions of J in these equations are calculated at J = j cM . Solving this equation we obtain the following function The only missing function is the field at the injecting boundary layer. This field profile is the solution of the semi-infinite problem (27) to (31), with J(σ; ǫ) = j cM + ǫJ (1) (σ) given by (B11). We can write (B11) in another form that suggests a more transparent interpretation: The terms on the right side of (B12) clearly display the balance between the area lost by the motion of the old wavefront at Y 1 (s) and the excess area under the injecting layer. J (1) (σ) increases linearly with σ until the excess area under the injecting layer increases with σ at least linearly. (3) (1)