Stationary states and phase diagram for a model of the Gunn effect under realistic boundary conditions

A general formulation of boundary conditions for semiconductor-metal contacts follows from a phenomenological procedure sketched here. The resulting boundary conditions, which incorporate only physically well-defined parameters, are used to study the classical unipolar drift-diffusion model for the Gunn effect. The analysis of its stationary solutions reveals the presence of bistability and hysteresis for a certain range of contact parameters. Several types of Gunn effect are predicted to occur in the model, when no stable stationary solution exists, depending on the value of the parameters of the injecting contact appearing in the boundary condition. In this way, the critical role played by contacts in the Gunn effect is clearly stablished.


I. INTRODUCTION
The Gunn effect is an ubiquitous phenomenon in many semiconductor samples presenting negative differential resistance and subject to voltage bias conditions [1][2][3][4]. In a nutshell, the negative differential resistance makes it possible the existence of a variety of pulses and wavefronts, which may be stabilized by the bias condition. Then a periodic shedding of waves by the injecting contact results in the periodic oscillation of the current through an external circuit, which constitutes the signature of the Gunn effect. There is a vast literature on this topic, despite of which different basic questions concerning the Gunn effect remain poorly understood. Paramount among these, there are the questions about which are the correct boundary conditions and, given these, how to describe all the stages of the Gunn oscillation. The lack of a precise formulation of the boundary conditions imposed by contacts on semiconductors and of a simple analytic treatment to analyze the Gunn oscillations, has not allowed to clarify in a precise way the role played by contacts in the Gunn effect. It is worth noting that clarifying this point would open, for instance, the possibility of extracting information about the contacts from the analysis of the Gunn oscillations themselves, a subject of considerable interest for applied researchers.
Recently progress has been made towards answering reasonably these two questions. On the one hand, ideas from Irreversible Thermodynamics [5] have been used to derive satisfactory boundary conditions for metalsemiconductor and other contacts in a general way [6][7][8].
Previously the usual boundary conditions used in driftdiffusion semiconductor models were: (i) periodic, [9] (ii) charge neutrality, [10] (iii) fixed field, [3,4] and (iv) control current-field characteristics of the contact [11] plus phenomenological assumptions such as the "contact length" [3]. As boundary conditions (b.c.) for a semiconductor presenting the Gunn effect, these conditions rank from clearly wrong (no current oscillation appears if the b.c. are periodic) to unsatisfactory because of their ad hoc character. Thus even when numerical simulations display the Gunn effect, the question of whether these results describe the real physical system where different contacts are present, is usually raised. In this paper, we shall present a simple derivation of appropriate b.c. for an ideal metal-semiconductor (MS) contact and use them to analyze the Gunn effect in Kroemer's model for bulk n-GaAs. Our description makes it clear which part of the derivation follows from general principles and which part includes the input from the physics of contacts.
Concerning asymptotic descriptions of the Gunn effect which delve deeper than just numerical simulations of drift-diffusion models, some progress has been made re-cently [12][13][14][15]. A detailed treatment of this topic can be found in Ref. [16].
The rest of the paper is as follows. In Section II we present our derivation of b.c. for ideal MS contacts and briefly discuss some other possibilities. Kroemer's model and its stationary solutions for these b.c. are analyzed in Section III. It is found that bistability between stationary solutions is possible for certain bias ranges depending on the values of certain dimensionless contact parameters i 0 , α 0 , which are a combination of its effective density of states, barrier height, Richardson's constant, doping and temperature. Different types of Gunn effect, namely, charge monopoles (moving charge accumulation and depletion layers), and charge dipoles, (high and low field solitary waves), are predicted to appear depending on these contact parameters, when no stable stationary solution exists. In Section V we discuss our results, whereas Appendix A is devoted to technical matters related to the main text.

II. BOUNDARY CONDITIONS
The aim of this section is to present a systematic procedure to derive b.c. at semiconductor contacts, established in previous works [6][7][8]. As a general rule, the method applies for non-degenerate semiconductors under moderate temperatures, that is, when thermionic emission is the dominant transport process at the contact. Hence several contacts of interest, like ideal and nonideal metal-semiconductor (MS) contacts or any type of heterojunction contact, can be modeled. Depending on the material parameters both limiting as well as ohmic contacts may then be described. It is worth noting that a precise modeling of this type of contacts may help to clarify the role played by other types of contacts used in semiconductor systems, e.g. those in which thermionic-field or field emission processes dominate [17][18][19], for which a precise description, in the sense of the present paper, is not yet available. For the sake of clarity, the method will be presented along with its application to the case of an ideal MS contact. Other contacts have been considered in previous papers [6][7][8].
Let us consider an ideal MS contact. Due to the presence of the contact, the magnitudes describing the physical properties of the system, e.g electron density, electric field, electron energy, . . . , may be discontinuous at that point. In addition, singular contributions localized at the contact itself, e.g. electron density at interface states (when they are present), may also occur. As a consequence, a given physical magnitude, d(x,t), can be decomposed as follows (1) where d n , d m , d s refer to the values in the semiconductor (n), metal (m), and surface (s) parts , respectively (when no singular contribution is present, d s vanishes). Moreover, θ(x) is Heaviside's unit step function and δ(x) Dirac's delta function. They are introduced in order to represent the discontinuity across the contact and the singular contributions, respectively. In writing Eq.(1), a one dimensional description of the system has been assumed, with the contact being located atx = 0 and the metal [semiconductor] on its left [right]. By means of this type of decomposition, b.c. can be systematically derived.
Our procedure consists of two steps: a) the identification of the relevant magnitudes describing the transport processes through the contact, and b) the derivation of precise laws describing such processes, which relate the relevant magnitudes at the contact and which constitute the desired b.c.. For the first step, use will be made of a phenomenological formulation of transport through semiconductor junctions [6] while for the second of Shockley-Read-Hall (SRH) statistics [20,21].
Let us consider a given magnitude, d(x,t), satisfying a standard balance equation of the form [5] ∂d(x,t) ∂t where J d (x,t) and σ d (x,t) refer to the current and net rate production associated to the magnitude d, respectively. It is not difficult to show that if similar balance equations were to be satisfied on each side of the junction and if surface fluxes only exist along the interface (what in a one dimensional description means J d,s (t) = 0), then the following balance equation should be satisfied at the contact, [22] ∂d s (t) ∂t Now, we can proceed to calculate the net rate of entropy production at the contact, which will allow us to identify the relevant magnitudes describing the transport processes trough the contact. To begin with, we consider the balance equation for the total energy of the system. As this is a conserved quantity, we simply have As for an ideal MS contact no interface states are present, and hence no net charge or mass is accumulated at the contact, the total energy at the contact coincides with the surface internal energy, u s = e s [22]. Hence, the balance of internal energy is described directly through Eq.(4) or alternatively through with σ u,s (t) = [J u,n (0,t) − J e,n (0,t)] − [J u,m (0,t) − J e,m (0,t)]. In the previous expression, we have introduced explicitly the flux of internal energy (equivalent to the heat flux), which is in general different from the flux of total energy. Furthermore the Gibbs equation for an ideal contact is [22] T ds s = du s (no interface states are present), where s s is the surface entropy and T the temperature. By assuming the contact to be in local equilibrium, one then has T ∂s s (t)/∂t = ∂u s (t)/∂t, which, after using Eq.(5), gives rise to the balance equation for the entropy with the entropy production given by A more explicit expression for σ s,s is obtained once the bulk expressions for the fluxes are introduced on the r.h.s. of Eq.(7). These expressions can be found elsewhere [5,6]. One has, where E F,a refers to the electron Fermi level (or chemical potential) and J a to the electron number density current. Substituting into Eq.(7), we simply have where in the second line use has been made of the continuity of the electron number density current at an ideal MS contact (this continuity follows from the corresponding balance equation for the electron number density by imposing that no carriers are accumulated (n s (t) = 0) nor created (σ n,s (t) = 0) at the contact).
The final expression to be used in what follows, is obtained by introducing the electron quasi-Fermi levels, F a (0,t) = E F,a (0,t) − eV a (0,t). Here V a (0,t) is the electric potential (which is continuous through an abrupt junction) and e > 0 is minus the charge of the electron. We then arrive at the desired expression, Eq. (10) shows directly that the relevant magnitudes describing an ideal MS contact are the electron flux (electron current density divided by e), J n (0,t), and the discontinuity in the electron quasi-Fermi levels, (F n (0,t) − F m (0,t)), which plays the role of "thermodynamic force", [5]. Both flux and force vanish at equilibrium, and we assume (in accordance with the basic tenets of Irreversible Thermodynamics [5]) that there is a relation between them. When the fundamental relation between flux and force is specified, this relation is exactly the sought-after boundary condition at the contact.
The relation between J n (0,t) and (F n (0,t) − F m (0,t)) should involve more information about the physics of the contact. First of all, let us notice that the entropy production (10) is formally equivalent to the expression corresponding to generation-recombination processes [23] (or in general, to any activated process, such as unimolecular chemical reactions [5] or surface adsorption [24]), provided J n (0,t) is identified with the net rate of the process. From this comparison we then conclude that the transport through an ideal MS contact may be described as an elementary kinetic process of the form: q m ⇀ ↽ q n , where q n and q m represent the carriers in the semiconductor and metal, respectively, with the net rate of the process being equal to J n (0,t). The kinetics of such a process can be described, for instance, by means of SRH statistics. As it is well known, this description relates the kinetic rate of the process, in our case J n (0,t), to the affinity, in our case the difference in quasi-Fermi levels, in agreement with our former general treatment. As it is shown in Appendix A, we obtain the following relation between the net rate of the process, J n (0,t), and the jump of the quasi-Fermi level at the MS contact, (F n (0,t) − F m (0,t)): Here , is the electron energy, with E 0 C the bottom of the semiconductor conduction band and V n the electric potential at the semiconductor surface. Moreover, β = (kT ) −1 , where k is the Boltzmann constant. Eq.(11) is the sought b.c. for the ideal MS contact. It can be written more explicitly by using the expressionñ = N C e β(Fn−EC ) , which holds for non-degenerate semiconductors withñ being the semiconductor electron number density and N C the effective density of states. We then have, In the non-degenerate case, λ 0 only depends on T (see Appendix A). This result ends the derivation of the b.c. for ideal MS contacts which we will use for the rest of this paper. It is worth emphasizing at this point that, as mentioned at the beginning of this section, the procedure we have sketched here for the ideal MS contact, can be applied to several other types of contacts. These include non-ideal MS contacts with the presence of interface states [7] and unipolar or bipolar heterojunction contacts with or without interface states [7]. Moreover, by adding a few assumptions one can handle non-abrupt contacts [8].
Note that for a MS contact located atx =L, that is, with the metal [semiconductor] on the right [left] hand side of the contact, the corresponding b.c. is We can compare our result, Eq.(12) for the ideal MS contact with the corresponding one reported in [25] [see also p. 261 of Ref. [26]]. Then we can identify λ i = A i T 2 /e, i = 0, L, where A i is the Richardson constant for the semiconductor in contact with the metal located at i = 0,L. Theoretically, A i , and hence λ i , would depend only on the given semiconductor but not on the metal [26]. However, in practice A i is taken as a phenomenological parameter and it can not only depend on the metal but also on the preparation procedures [27]. On the other hand, basic energetic arguments lead immediately to the following rule for the contact barrier height [26,28]: Here φ i M is the work function of the metal in the MS contact located atx = i and χ is the semiconductor electron affinity. For covalent semiconductors, the validity of this rule has been put under question for the last five decades [17,18]. However, recently it has been shown that even for this type of semiconductors, if accurately growth materials are used, good agreement is obtained with this simple rule [29,30]. Notice, that when such materials are not used, as very often happens, the contact formed turns out to be non-ideal. This is so because it is difficult to avoid that very thin insulating layers and/or interface states may be present at the contact [18,19]. Hence, a non-ideal description of the contact should be used, with for instance, a bias dependent relation for the barrier height which would include the effects of the insulating layer and of the interface states. Such contacts have been described in detail in Ref. [7,8]. In particular, it has been shown that for this non-ideal contacts there is no simple description using directly equations such as (11) or (12), for contact non-stationary effects induced by interface sates (not present for ideal MS contacts) introduce additional terms not present in Eq. (11). These terms could be of importance when describing naturally non-stationary phenomena such as the Gunn effect with non-ideal contacts, and will be considered in future works.
Eqs. (12) and (13) can be rewritten in terms of the electric field by using the Poisson equation to eliminate the electron density from them: Here the upper [lower] sign holds for i = 0 [i =L]; ε andñ 0 are the bulk semiconductor permitivity, and its doping, respectively.
We have thus shown how our procedure allows us to derive explicit and precise expressions for the b.c. imposed by a given contact. A first important consequence of this method can be drawn directly from Eq. (14), which is simply a relation between the normal derivative of the electric field and the current density at the contact. Examining our derivation shows that this result is simply a consequence of the use of kinetic models to describe the exchange of carriers through the contact. Hence, one should expect that b.c. derived in this way, will result in relationships between the normal derivative of the electric field and the current density at the contact. It is easily seen that (if diffusion effects can be neglected) our b.c. can be transformed into a Kroemer's type contact current-field control characteristics [11] (see next section). However, unlike previous models following Kroemer's approach [3,4,12,31], our control characteristics is the result of a physically precise derivation, and therefore only parameters which are physically well-defined appear in it. In particular, to use our control characteristics we do not have to invoke ad hoc assumptions involving new parameters such as the contact length [4,31].
To facilitate the analysis in the rest of the paper, it is convenient to rewrite Eq. (14) in dimensionless units. This greatly reduces the number of relevant parameters. Our dimensionless electric field E, electron densities n, current densities j(x, t), time t and position x are measured in units ofẼ 0 ,ñ 0 , eñ 0 µ 0Ẽ0 , ε/(eµ 0ñ0 ) and εẼ 0 /(eñ 0 ), respectively [32]. In these equations,Ẽ 0 and µ 0 are an electric field and the zero-field electron mobility typical of the processes occurring in the bulk of the semiconductor (see below). Then Eq. (14) becomes where we have defined As mentioned before, here the upper [lower] sign refers to i = 0 [i = L]. It is worth noting that α i is always a positive quantity (because λ i is) while i i does not have a definite sign: it depends basically on the value of the barrier height, φ i b , and on the doping value,ñ 0 . It should be noted that there are important restrictions on the possible values of the contact current density which are due to the fact that in Eqs. (12) and (13) the electron density, n, is a positive quantity: [Eqs. (16) and (17) imply that α −1 is always a positive quantity, equal to the maximum current density which the contact can provide]. These restrictions on the current are reminiscent of the rectifying properties of MS contacts. In practice, they only impose a real limitation for the case of true rectifying contacts (when one of the j sat i is small). Otherwise, i.e., for large values of j sat i , an ohmic contact is obtained which does not impose a real limitation on the current.
In order to analyze the influence of the derived b.c. on the Gunn instability, we will assume a sample formed by a certain semiconductor (able to display the Gunn effect) and by two MS contacts implemented on it. The resulting b.c. are As discussed above, for a given semiconductor the values of the contact parameters may vary somewhat depending on the metal used in the contact and on the preparation procedures. For instance α i may vary two orders of magnitude, from about 0.3 to 33.4, if we use the experimental values of Richardson's constant for GaAs reported in [27]. Similarly, the values of i i may also span two orders of magnitude, from about 0.03 to 4.01, due to the variation of α i , if we fix the barrier height, φ b ≈ 0.2 V (corresponding to Al [29]), and the donor density is 10 14 cm −3 . Thus there is a rather wide range of parameter values for the contacts, corresponding to a large variety of situations which will be described in this paper.
Lastly, the applied bias V , defined as eV (t) = F m (L, t) −F m (0, t), can be expressed as follows: In the previous expressions, V andφ i b are in units of εẼ 2 0 /eñ 0 andF m in units of Very frequently the analysis of the Gunn effect under dc voltage bias is carried out by using the opposite sign for the electric field: E = −E. Then the b.c. become With these definitions, the dc voltage is . Instead of working with the voltage V , it is more convenient to use the average electric field on the semiconductor sample, In what follows, negative voltages, V < 0, will be considered such as φ > 0. With these conventions, the carriers go from the cathode (injecting contact) at x = 0 to the anode (receiving contact) at x = L.

A. Kroemer's model
The unipolar drift-diffusion model for the Gunn effect proposed by Kroemer [2,33], is generally accepted to provide a rather complete description of the main features of this effect. Yet it is simple enough to allow very detailed asymptotic analysis; other important models such as Büttiker and Thomas's [34] incorporate more detailed physics but their study is technically more demanding. In the dimensionless units described above, Kroemer's model is Eq. (23) 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 and for specific calculations we shall use [33] 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 results using other curves having the same shape are similar. If v(E) does not reach a minimum but saturates instead as E → ∞, not all the monopole and dipole waves which we have found occur. Thus we have chosen the velocity curve that yields the richest dynamical behavior. The behavior of Kroemer's model with saturating velocity will be commented upon in the discussion. The dc bias φ is the average electric field on the semiconductor sample. Eqs. (23)-(24) need to be solved with an appropriate initial field profile E(x, 0) and subject to the following mixed boundary conditions resulting from substituting j(x, t) = J(t) − ∂E/∂t [from (23)] into (21)- (22): 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 , as will be seen in the following sections.
For typical n-GaAs data, δ ≪ 1 and L ≫ 1 [12]. In this limit, we shall find approximate solutions to the initial boundary value problem (23)- (27) for E(x, t) and J(t). Strictly speaking, the simple asymptotic description that follows holds in the limit L → ∞, even when δ = O(1) [15]. 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 [12,35,36].
To take advantage of this limit, we will use the following rescaled time and length, Then Eqs. (23)-(24) become 1 0 E(y, s) dy = φ.

B. Stationary states and its stability
Before describing the Gunn effect in the present model, it is convenient to discuss how to construct the stationary solutions of the model in the limit ǫ ≪ 1 and δ ≪ 1.
(In the case δ = O(1) the procedure is slightly more complicated and we shall omit the corresponding details; see [32]. In this section we shall analyze the stationary states of Kroemer's model in n-GaAs [2,33] under dc voltage bias with the new boundary conditions (21)- (22). Our work is based upon previous asymptotic and numerical studies of this and related models [12][13][14][15].
In this asymptotic limit, any stationary solution can be described as composed of outer and inner solutions: the outer bulk solution is a piecewise constant field profile valid everywhere except for two narrow boundary layers located at the contacts and, for particular values of the current density, a narrow transition layer somewhere in the middle of the sample (see Fig. 1 and explanation  below). First of all, if we ignore inner solutions, the stationary state solves the equations except for particular values of J which will be specified below. These equations result from retaining only order-one terms in Eqs. (29) and (30), and assuming E(y) =const. Then for those values of φ such that the outer solution (31) is compatible with the boundary conditions, we have J = v(φ) + O(ǫ). Let us denote by Then the outer (bulk) field profile will be E(y) = E i , i = 1, 3, depending on the value of the bias φ.
At y = 0 and y = 1 there are boundary layers, which we will call injecting and receiving layers, respectively. E(y; J), the field at the injecting boundary layer of width O(ǫ) at y = 0, obeys [we ignore the O(δ) diffusive term]: The analysis of Eqs. (32) and (33) is more easily carried out if we express the derivative b.c. Eq.(33) in terms of a b.c. for the electric field at the contact, E(0; J) . We can obtain E(0; J) from (33) by using (32) to eliminate ∂E(0; J)/∂y. The result is that E = E(0; J) solves Notice that the contact curve j c (E) has the same extrema as the velocity curve v(E) and saturates for high electric fields to the value j sat 0 , defined in Eq. (18). Kroemer's contact characteristic for shallow-barrier metalsemiconductor contacts presented in [33], corresponds to a particular case of our model in which the electrons in the metal are assumed to be in equilibrium with those of the semiconductor near the contact. For this case, one would take α 0 → 0 with |α 0 i 0 | < ∞, so that j c (E) would be then proportional to v(E). In contradistinction with Kroemer's contact characteristic, the general curve j c (E) may intersect the bulk velocity curve v(E). The main difference between these two cases is that a Gunn effect mediated by charge dipole waves is seen only if j c (E) intersects the bulk velocity curve v(E). If (34) has a solution, Eq. (32) indicates that E(y; J) approaches one of the solutions of (31) as we leave the boundary layer. The boundary layer at the receiving contact y = 1 is a much narrower diffusive boundary layer of width O(ǫδ). The field there is [32] where The idea now is to fix J and to discuss for which values of J the above construction yields a stationary solution. Additionally, its stability will be considered. Clearly we may distinguish different cases according to the values of the contact parameters (i i , α i ), i = 0, L. In what follows we shall assume for the sake of simplicity that the boundary layer equations at the receiving contact, Eqs. (37) and (38), always have a solution, and hence only the parameters of the cathode, (i 0 , α 0 ), need to be considered.
The general situation encountered when constructing the stationary solutions is the following. For each value of J, there are one or three values of the contact electric field at x = 0, which are solutions of Eq. (34). We shall denote these field values by E ci (J), with E c1 (J) < E c2 (J) < E c3 (J). Then the field profile in the injecting boundary layer is a monotonic solution of Eq. (32), which joins E ci (J) (i = 1, 2, 3) to one of the solutions of J = v(E) (outer solution). Furthermore, the outer solution may be a constant field profile given by E(y) = E l (J) (l = 1, 3) which extends to the end of the sample, where a narrow receiving boundary layer exists [see Fig. 1]. For such an electric field profile, the bias is φ ≈ E i (J). The corresponding J-φ characteristics satisfies J ≈ v(φ). By this construction, we identify the portions of the J-φ characteristics which follow the first or the third branch of v(E) (see the details below). Other portions of the J-φ characteristics are flat, with J = J f for certain constant values of the current. The corresponding outer field profile is step- Fig. 1]. ∆Y is chosen so as to satisfy the bias condition φ ≈ E 2 (J f )∆Y + (1 − ∆Y )E i (J f ). The flat part of the J-φ characteristics corresponds to a bias range E 1 (J f ) < φ < E 3 (J f ) (see the details below). Finally, when J is near the value j sat 0 , the field at the injecting contact is very large, the contact region is almost depleted of electrons and its extension, = 1, 3), may be comparable to the length of the sample. Assuming the extension of the depletion layer at the injecting contact is less than the sample length, the corresponding bias is 1, 3). The characteristics tends to saturate at j sat 0 .
Following this general scheme, different possibilities may be distinguished according to the relative values of j cm < j cM < j sat 0 with respect to v m < v M . Here, and v k , with k = m, M , refer to the minimum (k = m) and maximum (k = M ) of the contact, j c (E), and velocity, v(E), curves, respectively . We now discuss the different cases which appear for our velocity curve.
In this case, see Fig. 2 , we have for 0 < J < J cM a class of solutions joining E c1 (J) and E 1 (J), with voltage φ ≈ E 1 (J). For 0 < φ < E 1 (j cM ) the curve J-φ then follows the first branch of v(E). Furthermore, a second class of solutions joining E c3 (J) and E 1 (J), will exist for j cm < J < j sat 0 . In this case for J not near j sat 0 , the voltage is given by φ ≈ E 1 (J), and for J near j sat Then, in the characteristics the third branch starting at φ ≈ E 1 (j cm ) follows the first branch of v(E) at the beginning, till it tends to saturate to j sat 0 for larger voltages (see Fig. 2 ). Joining these two classes of solutions, there exists a third class for j cm < J < j cM which joins E c2 (J) to E 1 (J), with φ ≈ E 1 (J). These solutions are unstable, and they give rise to the second branch in the characteristics, Fig.  2, which also tend to follow the first branch of v(E). Note that for voltages E 1 (j c m) < φ < E 1 (j c M ), the two classes of stable stationary solutions coexist [see inset at the bottom of Fig. 2]. Hysteresis between them is then possible.
2. jcm < vm < jcM < j sat 0 < vM In this case, see Fig. 3 , the description is very similar to the previous case, except on what concerns the third branch of the J-φ characteristics. Now, this branch is composed of two types of solutions: i) for j cm < J < v m there is a class of solutions joining E c3 (J) and E 1 (J). Most of the time one has φ ≈ E 1 (J), except for J near v m that the solution is step-like with φ ≈ E m ∆Y + (1 − ∆Y )E 1 (v m ). We expect that these solutions become unstable on a bias range which is a subinterval of E M < φ < E m [37,38]. Then a Gunn effect mediated by moving charge monopoles (which are charge depletion layers) might appear (see companion paper [16]). For v m < J < j sat 0 , there is a class of solution joining E c3 (J) and E 3 (J), with φ ≈ E 3 (J) for J not near j sat Then, the third branch of the J-φ curve starts following the first branch of the v(E) curve for E 1 (j cm ) < φ < E 1 (v m ), then it presents a flat region for E 1 (v m ) < φ < E m with J = v m , corresponding to the step-like solutions, and finally a region for E m < φ which starts following the third branch of the v(E) curve and tends to saturate to j sat 0 for larger voltages. Note the presence again of bistability for voltages E 1 (j cm ) < φ < E 1 (j cM ) [inset at the bottom of Fig. 3] and hence the possibility of hysteresis.
A similar situation to the one depicted above would appear for j cm < j cM < v m < j sat The main difference for this case with respect to the previous ones, relies on the second and third branches (see Fig. 4). Now, the third branch of the J-φ curve involves one single class of solutions joining E c3 (J) and E 3 (J), with voltages φ ≈ E 3 (J) for J not near j sat 0 , and The second (unstable) branch is formed of two classes of solutions: one class for i 0 < J < j cM starting at E c2 (J) and ending at E 1 (J), and another class for j cm < J < i 0 , starting at E c2 (J) and ending at E 3 (J). For J ≈ i 0 , these solutions are step-like with the voltage given through φ ≈ E c2 (i 0 )∆Y + (1 − ∆Y )E i (J), with i = 1, 3 depending on the class of solutions considered. The (unstable) branch in the J-φ curve then starts following the first branch of the v(E) curve for E 1 (j cM ) < φ < E 1 (i 0 ), then it presents a flat portion for E 1 (i 0 ) < φ < E 3 (i 0 ) with J ≈ i 0 , and it ends following the third branch of the v(E) curve for E 3 (j cm ) < φ < E 3 (i 0 ). Note that for voltages E 1 (j cM ) < φ < E 3 (j cm ) there is no stable stationary solution [see inset at the bottom of Fig. 4]. Thus we expect that the usual Gunn effect (mediated by moving charge dipoles) will be present for these values of the bias (see the companion paper [16]).
A similar situation appears for v m < j cm < j cM < v M < j sat 0 .
4. vm < jcm < vM < jcM < j sat 0 In this case (see Fig. 5), the third branch of the J-φ characteristics is described as in the previous case. The first branch is composed of two types of solutions: one class, for 0 < J < v M , joining E c1 (J) and E 1 (J), and the other, for v M < J < j cM , joining E c1 (J) and E 3 (J). For the first type of solutions, one has φ ≈ E 1 (J) except . These step-like solutions are expected to become unstable in a subinterval of E M < φ < E m [37,38]. Then a Gunn effect mediated by moving charge monopoles (which are charge accumulation layers) might appear (see the companion paper [16]). Thus, the first branch starts following the first branch of the v(E) curve for 0 < φ < E M , then it presents a flat portion for E M < φ < E 3 (v M ), with J = v M , and ends following the third branch of the v(E) for E 3 (v M ) < φ < E 3 (j cM ). The second (unstable) branch of the J-φ curve is formed by a class of solutions that starts at E c2 (J) and ends at E 3 (J), with φ ≈ E 3 (J). Then, this branch follows the third branch of the v(E) curve, for E 3 (j cm ) < φ < E 3 (j cM ). Note that for this range of bias, two stationary stable solution coexists [see inset at the bottom of Fig. 5] and hysteresis may appear.

C. Phase diagram
By collecting the information obtained in the previous subsections, the phase diagram describing the different behaviours of the system can be sketched, in terms of the injecting contact parameters, i 0 and α 0 , Fig.  6. Stable, non-oscillatory stationary solutions are expected for values of i 0 and α 0 such as j sat Otherwise, some kind of oscillatory solution should be found. Charge accumulation monopoles appear for j cM > v M (or equivalently for i 0 > v M ), charge dipoles for v m < j cm < j cM < v M , that is, for v m < i 0 < v M , and charge depletion monopoles for j sat 0 > v m with j cm < v m (i 0 < v m ). For completeness, also the separation between low field and high field dipoles, discussed in the companion paper [16], has been depicted. It is worth noting that this rich phenomenology of oscillatory states appears just by changing the value of the contact parameters. This fact should be taken into account in analyzing the Gunn effect in real systems, where, as mentioned before, a wide range of values for the contact parameters, depending on the metal used and preparation procedures, may appear.

IV. DISCUSSION
We have presented a general formulation for the derivation of the boundary conditions imposed by metalsemiconductor contacts on semiconductor systems. According to this general formulation, the appropriate boundary conditions for ideal metal-semiconductor contacts are linear relations between the normal derivative of the electric field at the contacts and the electron current there. For the classical unipolar drift-diffusion Kroemer's model of the Gunn effect, these boundary conditions are of mixed type. In this paper, we have investigated how the boundary conditions for ideal metal-semiconductor contacts affect the stationary solutions of the Kroemer model, and their stability. Depending on the values of the contact parameters, bistability and hysteresis may appear. Moreover, for some range of parameters no stable stationary solution is expected to occur. In those parameter ranges we expect to find the Gunn effect. Numerical simulations show that different types of Gunn effect appear, mediated by a variety of waves: (i) charge monopole accumulation wavefronts, (ii) monopole depletion wavefronts, or (iii) charge dipole waves (high and low electric field domains). Why these types of Gunn effect appear in the simulations will be explained by the asymptotic theory of the companion paper [16]. It suffices to say that without this theory we would have missed significant possible instabilities. For example (ii) seems to have been missed by earlier workers, in spite of past extensive simulations of Kroemer's model [3]. With our boundary conditions, the previously described types of Gunn effect are found in the following ranges of dimensionless critical contact currents: (i) j cM > v M (ii) j cm < v m and j sat 0 > v m , and (iii) v m < j cm < j cM < v M . Here j cM = j c (E M ), j cm = j c (E m ), j sat 0 are the critical currents, and v m and v M are the minimum and maximum values of the electron drift velocity v(E), E > 0. When we want to characterize experimental samples displaying the Gunn effect, it is important to bear in mind the great influence of the contact parameters on the type of wave mediating the Gunn effect. A wide range of values for these contact parameters may be obtained depending on the type of metal used or on the contact preparation procedure followed.
N-shaped velocity curves occur naturally in recently observed self-sustained oscillations in weakly-coupled ndoped GaAs/AlAs superlattices (see Ref. [39] for the most complete data so far). In these superlattices there is strong indirect evidence for a Gunn effect mediated by charge accumulation monopoles through photocurrent and photoluminescence measurements [40]. It is hard to say at this point which other possibilities of those found in our analysis might be realizable in these systems. An important issue to be decided is the form of the boundary conditions. Our analysis needs to be modified in order to be extended to these systems, as quantum tunneling plays an essential role in the injection of carriers through contact regions.
The most used velocity curves v(E) for the classical Gunn effect in bulk n-GaAs lack the third branch after v m . The reason is that avalanche breakdown appears at electric fields smaller than E m . The avalanche field is smaller for the longer samples needed to observe the Gunn effect and this precludes reaching the high fields on the third branch of v(E). Then low field dipole domains and charge depletion monopoles are not observed in the usual bulk samples or in strongly coupled superlattices with wide minibands, which are analogous to them [34].

V. ACKNOWLEDGMENTS
We thank Dr. M. Bergmann for fruitful conversations, and acknowledge financial support from the the Spanish DGICYT through grants PB94-0375 and PB95-0881, and from the EC Human Capital and Mobility Programme contract ERBCHRXCT930413. One of us (G.G.) acknowledge support by the Generalitat de Catalunya. Let us consider the elementary kinetic process: q m ⇀ ↽ q n describing the charge transport through the junction. By assuming the validity of the SRH statistics to describe this process, the following expression for its kinetic rate, J n , can be obtained [23]: where D a (E a ), a = n, m, is the density of states of system a, f a (E a ) its occupation function, given through the FD distribution f a (E a ) = (1 + e β(Ea−EF a ) ) −1 , with E F,a the corresponding Fermi level. γ mn (E m , E n ) [resp. γ nm (E m , E n )] is the probability per unit time for the transition between states of energy E m and E n [resp. E n and E m ]. At equilibrium, we must have J n = 0 and F m = F n , with F a = E F,a − eV a being the corresponding quasi-Fermi levels. This implies E F,n −E F,m = e (V n −V m ) = 0 (using the assumed continuity of the electric potential at the contact). These equations follow from Eq. (A1) if the latter is supplemented with the following detailed balance relation: A term β (V n − V m ) has to be added to the argument of the exponential in (A2) if V n = V m ; see [8] for a more general case. We now substitute Eq. (A2) into Eq. (A1) and use the equations C is bottom of the semiconductor conduction band and eφ 0 b is the height of the contact barrier). After straightforward manipulations, we derive (11), in which the transition coefficient λ 0 is: When the semiconductor is non-degenerate, we may approximate 1 − f n (E n ) ≈ 1, whereas for a metal we may approximate f m (E m ) by its equilibrium value. Then, for this case, λ 0 is a function of T only.