Two-valley Hydrodynamical Models for Electron Transport in Gallium Arsenide: Simulation of Gunn Oscillations

To accurately describe non-stationary carrier transport in GaAs devices, it is necessary to use Monte Carlo methods or hydrodynamical (or energy transport) models which incorporate population transfer between valleys. We present here simulations of Gunn oscillations in a GaAs diode based on two-valley hydrodynamical models: the classic Bløtekjær model and two recently developed moment expansion models. Scattering parameters within the models are obtained from homogeneous Monte Carlo simulations, and these are compared against expressions in the literature. Comparisons are made between our hydrodynamical results, existing work, and direct Monte Carlo simulations of the oscillator device.


INTRODUCTION
In microelectronics, near-micron n þnn þ diodes based on compound semiconductors such as GaAs, InP or GaN are currently used as generators of microwave radiation [28]. In order to optimise device performance it is very useful to have accurate numerical simulations of the transient and oscillatory behaviour of the devices [1]. When using compound semiconductors for high frequency applications, usually one deals with multi-valley band structures, and in these cases the transfer of carriers from one valley to another (which is at the root of Gunn oscillations) must be modelled. While Monte Carlo methods can be used to numerically solve the semiclassical Boltzmann equation describing carrier behaviour [17], the use of simpler transport models is generally preferred since their solution requires substantially less computational resources.
The standard drift-diffusion models [20] do not incorporate dynamical transfer of carriers between valleys and this renders them ill-suited for simulating timedependent high frequency phenomena. Generalizations of the drift-diffusion equations have therefore been sought, and these take the form of hydrodynamical models (named for their similarity with the equations of compressible fluid flow) and energy transport models.

Hydrodynamical and Energy Transport Models
Hydrodynamical semiconductor models are obtained from the infinite hierarchy of moment equations of the semiclassical Boltzmann transport equation by applying a suitable truncation procedure. This requires: (i) appropriate moments to be chosen for the expansion; (ii) closure of the hierarchy of equations by expressing the Nth-order moment in terms of lower-order ones; (iii) modelling of production terms which arise from the moments of the Boltzmann collision terms. Various closure assumptions have been suggested in the literature, leading to various hydrodynamical semiconductor models, including the Bløtekjaer [12] model which is often used in industrial simulation studies (see Ref. [2] for further references).
However, most of these closure assumptions are, at best, only phenomenological and lack a consistent physical and mathematical justification. For example, the Bløtekjaer model is inconsistent with the Onsager reciprocity principle of linear irreversible thermodynamics [2]. Also, for time-dependent phenomena such as Gunn oscillations, the energy flux relaxation time is of the order of the momentum relaxation time and therefore cannot be neglected as it is in the Bløtekjaer model.
A consistent mathematical approach for deriving hydrodynamical semiconductor models has recently been developed [2,5,7 -9,23] which utilizes a closure assumption based on the maximum entropy principle of extended thermodynamics [18,21] or equivalently the method of exponential moments [19]. In this approach, evolution equations for the heat flux and stress tensor are considered in addition to the usual balance equations for carrier density, momentum and energy. The resultant models are free from phenomenological assumptions and may enjoy important mathematical properties such as hyperbolicity. In the stationary case, linearizing the heat flux equation for small temperature gradients (Maxwellian iteration) results in an extension of Fourier's law which includes an additional convective term, and we speak of this as the Anile -Pennisi [5] model. Alternatively, extended semiconductor models have been derived [6,7] in which the stress and heat flux are treated as unknown variables, and these avoid reliance on Fourier's law altogether. A simplified model of this type, which we call the "reduced hyperbolic" model [3,4,16], proves to be particularly convenient to work with.
Multi-valley hydrodynamical models consider multiple populations of carriers, each represented by a set of transport equations. The carrier populations are coupled both through Poisson's equation for the electric field and through source terms which model the transfer of carriers between valleys. While some work has been done in deriving production terms directly from the Boltzmann collision operator [23], more often scattering processes are incorporated in hydrodynamical models through phenomenological expressions based on relaxation times for the semiconductor material. The relaxation times are typically derived as functions of carrier energy from Monte Carlo simulations in which the electric field is homogeneous and fixed [25].
An alternative to hydrodynamical models, energy transport models comprise the time-dependent balance equations for particle number and energy, supplemented by constitutive equations for the particle current and energy flux [10]. Because they make the implicit assumption that the relaxation times for velocity and energy flux are negligible compared to those for particle number and energy, these models are not appropriate for two-valley simulations such as we envisage here. In fact, as we shall see later, Monte Carlo simulations show that the momentum and energy flux relaxation times for the two populations can, under some circumstances, be of the same order as the particle number relaxation times [25]. Energy transport models will not be considered further here.

Numerical Studies of Gunn Oscillators
This work investigates the suitability of two-valley hydrodynamical models for studying non-stationary phenomena in GaAs devices. In particular, Gunn oscillations are simulated. We note that previous studies of hydrodynamical models have mainly considered stationary behaviour in silicon devices.
One of the earliest numerical simulations of oscillations in a GaAs diode is due to Curow and Hintz [14]. They adopted a single-valley model of the Bløtekjaer type to investigate Gunn devices in harmonic mode oscillator circuits.
In a series of papers, Gruzinski et al. [15] have studied oscillations in an InP diode using a single-valley hydrodynamical model. They obtain closure of the convective part of the equations by assuming the pressure and energy flux to be functions of the average electron energy only, and fit these functions to homogeneous Monte Carlo data. The source terms they use are relaxation-type expressions, where the relaxation times for momentum and energy are also taken to be functions of the average energy, again fit to homogeneous Monte Carlo data. Their results show reasonable agreement with those obtained by direct Monte Carlo simulation. We believe, however, that their approach has shortcomings: (i) it does not give a description of electron transfer between valleys, and therefore provides less detailed information than might be desired for opto-electronics applications; (ii) the closure assumptions apply only to the one-dimensional case and it is not clear how they could be generalized to higher dimensions; (iii) the mathematical nature of the model is not clear, and so the choice of a suitable numerical scheme is not obvious.
Chen, Jerome, Shu and Wang (CJSW) [13] have used a two-valley model of the Bløtekjaer type to simulate the GaAs device in a notched oscillator circuit. The relaxation terms they use in their model are derived from phenomenological considerations and cannot be considered as physically accurate. None-the-less, their simulations show sustained Gunn oscillations being produced by the device. We use the work of Chen et al. as a starting point for our studies: we consider the same basic device set-up, and we validate our numerical code against their results before performing simulations using different hydrodynamical models.

Outline of This Paper
We perform numerical simulations of GaAs devices for a range of hydrodynamical models in order to asses the validity of their underlying assumptions. In particular, results for the hydrodynamical models are compared to simulations performed using the direct Monte Carlo method, described in the second section. The same Monte Carlo code is also used to produce expressions for the relaxation times for GaAs.
For the purposes of this paper, hydrodynamical models are thought of as comprising two components: carrier transport equations with separate collision terms. We consider in the third section, three different transport models and three different collision models. Numerical solutions to the hydrodynamical models are generated using an adaptive mesh refinement code which is summarized in the fourth section. The GaAs diode simulated in this work is described in the fifth section along with the oscillator circuit to which it is coupled. In the sixth section, we present our results which include simulations of homogeneous GaAs ["Simulations of Homogeneous Samples"] and of the GaAs diode both at steady state ["Device Behaviour at Steady State"] and when producing Gunn oscillations ["Gunn Oscillations"]. We conclude in the seventh section with a discussion.

MONTE CARLO SIMULATIONS
We treat the case of unipolar devices in which the current is due only to electrons. These populate two distinct energy valleys which we assume (for consistency with the hydrodynamical models of "Two-valley hydrodynamical models" Section) to be spherical and parabolic.
In the Monte Carlo method [17], charge transport in a device is simulated by tracking the position, momentum and occupied energy valley of each of a large number of electrons. Motion of the electrons comprises a series of free flights interrupted by randomly resolved scattering events. Poisson's equation for the electric field is solved consistently with the distribution of the charge carriers. While the Monte Carlo method allows very detailed physical behaviour to be incorporated in a simulation, it suffers from the drawback that a great deal of computing time is required to track the large number of particles needed to produce accurate results.
We use a modified version of Tomizawa's Monte Carlo code [25] in which carrier-phonon and carrier-impurity interactions are modelled. The Monte Carlo code is used in two ways. Firstly, simulations of a homogeneous GaAs sample are performed in which the electric field is fixed and the particle positions need not be calculated. The results are used in "New Relaxation Times from Monte Carlo Simulations" Section to determine the relaxation times (as functions of particle energy) through which collision processes may be incorporated in hydrodynamical models. Secondly, in "Numerical Results" Section, we compare Monte Carlo results to solutions of hydrodynamical models for both homogeneous GaAs samples and a one-dimensional oscillator device (described in "The GaAs Diode and Circuit" Section). For a full device simulation, close to a million electrons are used.

TWO-VALLEY HYDRODYNAMICAL MODELS
We shall consider three different two-valley hydrodynamical models for electron transport in GaAs [see Sections "The Bløtekjaer Transport Model", "The Anile-Pennisi Transport Model" and "The Reduced Hyperbolic Transport Model"], together with three different models for collision processes [see Sections "The CJSW Collision Model", "The Tomizawa Collision Model" and "New Relaxation Times From Monte Carlo Simulations"].
Variables used in the models are the number density of carriers n, their mean velocity v, the carrier energy per unit volume E (measured relative to the minimum of the valley), the pressure p ¼ ð2=3ÞE 2 ð1=3Þm * nv 2 (where m* is the effective electron mass), the temperature T ¼ p=ðk B nÞ (where k B is the Boltzmann constant), and the total energy flux q. Two carrier types are considered here, and their associated variables are identified by subscripts 1 (for G-valley electrons) and 2 (for L-valley electrons).
The electric field E and electric potential f are determined from the charge distribution through Poisson's equation, which in one dimension reads where n D (x ) is the number density of donor atoms, e is the magnitude of the electron charge, and 1 is the dielectric constant of the material. For a GaAs device the following values are assumed for the lattice temperature and other physical parameters: where 1 0 is the vacuum dielectric constant, and m e is the electron mass. Standard units are used throughout this discussion except where noted.
In the following we consider three hydrodynamical transport models: the classic Bløtekjaer model [12] [see Section "The Bløtekjaer Transport Model"], the Anile -Pennisi model [5] [see Section "The Anile -Pennisi Transport Model"], and the reduced hyperbolic model [3,4] [see Section "The Reduced Hyperbolic Transport Model"]. Three different relaxation time models for the collision processes are also examined: expressions from Chen et al. [13] [see Section "The CJSW Collision Model"], expressions from Tomizawa [25] [see Section "The Tomizawa Collision Model"], and new relaxation times derived for the Tomizawa model from homogeneous Monte Carlo simulations [see Section "New Relaxation Times from Monte Carlo Simulations"].

The Bløtekjaer Transport Model
The first model we consider is the classic two-valley Bløtekjaer hydrodynamical model [12,25]. It consists of balance equations for the number density, momentum and energy of two populations of electrons, with heat conduction effects modelled by Fourier's law for both carrier types.
In one dimension the equations are › ›t Second-order derivatives of the carrier temperatures account for the conduction of heat, and following Chen et al. [13] we take the conduction coefficients to be where t p:1 ; t p:2 are the momentum relaxation times, described later. The Bløtekjaer model is widely used for semiconductor analysis. However, the use of Fourier's law to model heat conduction prevents the transport equations from being hyperbolic, and in addition introduces a parameter, k 0 , which is not determined by any theory: its value is chosen empirically based on comparisons with Monte Carlo simulations. In response to this, various alternative hydrodynamical transport models have been proposed in which heat conduction is incorporated in a physically consistent way. Two such models are described below.

The Anile -Pennisi Transport Model
Described in Ref. [5], this model uses the principles of extended thermodynamics [21] to replace Fourier's law for the heat flux with an alternative diffusion-like expression which is free of arbitrary parameters. Numerical simulations using this model have previously been performed by Romano and Russo [24].
Generalized to two carrier types, the evolution equations in one dimension are the same as for the Bløtekjaer model (2) where t q:1 ; t q:2 are the heat flux relaxation times, described later. Note that the expression for the heat flux in this model is determined uniquely by the theory and no arbitrary parameters (such as k 0 above) remain to be specified.

The Reduced Hyperbolic Transport Model
The third model we use is a generalization to two carrier types of another in the family of hydrodynamical transport models introduced by Anile et al. [2 -5,7,9]. Further, application of extended thermodynamics principles to the Anile -Pennisi model results in evolution equations in which thermal diffusion is not described by diffusion-like terms but instead by additional balance equations for the energy flux q (arising from higher moments of the Boltzmann equation).
The particular model adopted here, which we call the reduced hyperbolic model, is obtained as a limit of a more complete extended model when quadratic terms in the deviations from local thermal equilibrium are neglected. This model is described and used in numerical simulations by Anile, Nikiforakis and Pidatella [3,4].
In one dimension the equations are › ›t

The CJSW Collision Model
We consider here (and in the following two sections) three models for collision processes in GaAs to complete the transport models described above.
Firstly, in the article by Chen, Jerome, Shu and Wang (CJSW) [13] the two-valley Bløtekjaer model, Eq. (2), is completed by expressing the collision terms as giving values in seconds. The transfer rate t 21 n:12 from the first to the second valley is made "smooth" by inserting a seventh-order polynomial section into the function such that it is continuous up to the third derivative. For the threshold energy in this function, CJSW use various values, a ¼ 3:66 £ 10 220 ; 2:40 £ 10 220 ; 1:44 £ 10 220 ; 0:48 £ 10 220 J; while the width of the transition region is fixed by The crude relaxation time expressions used by CJSW enable them to reproduce in simulations the basic oscillatory behaviour of a GaAs device. In the following two sections we consider more realistic models based on Monte Carlo calculations.

The Tomizawa Collision Model
The second collision model we use comes from the book by Tomizawa [25] where it is incorporated in the two-valley Bløtekjaer transport model (2). The collision terms are ð›n i v i =›tÞ C ¼ 2v i ðn i =t p:i þ 2n i =t n:ij 2 n j =t n:ji Þ; ð8bÞ A simple form for the energy flux relaxation times is assumed: In "New Relaxation Times from Monte Carlo Simulations" Section, it is seen that this relationship between the relaxation times for energy flux and momentum is reasonable.

New Relaxation Times from Monte Carlo Simulations
Later in this work ("Numerical Results" Section), we compare numerical results for hydrodynamical models to Monte Carlo simulations for a GaAs device. In order to make these comparisons as meaningful as possible, we have calculated new relaxation times for the Tomizawa collision expressions ["The Tomizawa Collision Model" Section] using the same Monte Carlo code ["Monte Carlo Simulations" Section] that we use for the device simulations. As far as possible, the assumptions made in these Monte Carlo calculations are the same as those used for the hydrodynamical models.
In particular, parabolic energy bands are assumed in both cases. Our collision model thus comprises the Tomizawa expressions (8a) -(8d) and (10) An impurity concentration of 10 22 m 23 in the GaAs sample has been assumed. These relaxation time expressions are based on the original Tomizawa ones, but the coefficients have been altered so that a good fit is achieved with the Monte Carlo results. Note that we have also determined relaxation times t q:i for the energy fluxes. Figure 1 shows how these expressions for the relaxation times compare with the Monte Carlo data.

NUMERICAL SOLUTION OF HYDRODYNAMICAL MODELS
Approximate numerical solutions to the hydrodynamical models of "Two-Valley Hydrodynamical Models" Section are obtained using a computer code described in full in [16]. The code is designed to evolve solutions to a wide range of hydrodynamical models, and it has been validated using single-valley versions of the Bløtekjaer and reduced hyperbolic models [see Sections "The Bløtekjaer Transport Model" and "The Reduced Hyperbolic Transport Model"] in the case of a onedimensional silicon diode. The inclusion of two-valley models for GaAs within the code is straightforward. The main features of the code are summarized below.
The code solves the full time-dependent equations for the hydrodynamical models. Solutions are advanced in time using a first-order splitting approach which allows specialized numerical methods to be applied to the different components of the evolution equations. The source terms of the equations (the collision and forcing terms) are treated using the explicit second-order Runge -Kutta method for solving systems of ordinary differential equations [22]. While more sophisticated methods are available, this method is found to be adequate for the problem at hand. The left-hand sides of Eqs. (2), (4) and (5) (omitting the heat diffusion terms in the Bløtekjaer and Anile -Pennisi models) form hyperbolic systems of conservation laws, and are solved using the Slope Limiter Centred (SLIC) scheme [27]. This is a high resolution scheme which is second-order accurate when the solution is smooth and which resolves discontinuities without introducing unphysical oscillations. For the evolution of hydrodynamical models, the SLIC scheme has the advantage that it does not rely on the general solution being known to the Riemann problem for the system, and indeed for the Anile -Pennisi model this information is not available. Finally, for the Bløtekjaer and Anile -Pennisi models, the parabolic heat conduction terms must be treated. Uncoupled from the rest of the system, these take the form of linear scalar diffusion equations, and they are solved using the forward Euler method [22].
The code uses the adaptive mesh refinement (AMR) technique of Berger and Oliger [11] to enable the resolution of a simulation to vary in space and time in response to the behaviour of the evolving solution. This can greatly improve the speed with which numerical results are obtained. For most of the GaAs device simulations we describe here, a base resolution of Dx ¼ L/200 was used (where L ¼ 2 mm is the length of the device). In regions surrounding jumps in the doping profile n D , the resolution was increased to Dx ¼ L=800 via buffer regions of intermediate resolution Dx ¼ L/400. A minimum resolution of Dx ¼ L=400 was also used in regions where the gradient or the curvature of either of the velocity variables v i was judged by the code to be high. These criteria for dynamically varying the resolution have previously been found to be effective in producing accurate results for time-dependent simulations of a ballistic silicon diode [16].

THE GAAS DIODE AND CIRCUIT
Our studies are based on a GaAs notched diode, as used by Chen et al. [13]. This is coupled to an RLC tank circuit which stimulates Gunn oscillatory effects. The same configuration is used for both hydrodynamical and Monte Carlo simulations.
The (one-dimensional) diode has length L ¼ 2 mm and is doped with donor atoms in a profile n D ðxÞ ¼ In our studies, the transitions in the doping profile at the device junctions are discontinuous (in contrast to Ref. [13]. For both hydrodynamical and Monte Carlo simulations, the following initial data at time t ¼ 0 is used:  For the Monte Carlo simulations, this initial state is implemented by starting all particles off in the first valley with velocities randomly picked to give a Maxwellian distribution.
The following boundary conditions at x ¼ 0 and x ¼ L are used for the hydrodynamical simulations: Implementation of equivalent boundary conditions in the Monte Carlo simulations is not entirely straightforward. The approach we use is to extend the Monte Carlo domain to include short "buffer" regions beyond the specified edges of the device. During the course of a simulation the number density of electrons in these buffer regions is made to be constant, equal to the doping density n D at the contacts, by removing at random excess electrons, or creating replacement electrons in the first valley with velocities fitting a Maxwellian distribution consistent with Eq. (15). Simulation results are not found to be overly sensitive to the conditions at the boundaries.
Boundary conditions for Poisson's equation (1) are imposed by specifying the electric potential at the device contacts: Here the potential difference V is either a prescribed function of the simulation time (typically a constant value of 2 V), or it is determined by coupling the device to a system of ODEs modelling a simple circuit.
For simulating Gunn oscillations we use the following circuit model given by Chen et al. [13]. The voltage V across the device and the current I through the battery satisfy differential equations where the particle current in the device is calculated as where 1 is the dielectric constant of GaAs, L is the diode length, and its cross-sectional area is The bias voltage of the circuit is and the oscillator equations (17) are given the initial state when the circuit is engaged at time t 0 .

NUMERICAL RESULTS
We present numerical results for the two-valley hydrodynamical models of "Two-Valley Hydrodynamical Models" Section. Comparisons are made with the results of Chen et al. [13] [see "Comparison with Existing Results" Section], and between hydrodynamical and Monte Carlo simulations for a homogeneous sample [see "Simulations of Homogeneous Samples" Section], and for a diode both at steady state [see "Device Behaviour at Steady State" Section] and when coupled to an oscillator circuit [see "Gunn Oscillations" Section].

Comparison With Existing Results
To validate our numerical approach for two-valley hydrodynamical models ["Numerical Solution of Hydrodynamical Models" Section], we have reproduced the GaAs diode simulations of CJSW [3]. The semiconductor model used combines the Bløtekjaer transport equations of "The Bløtekjaer Transport Model" Section with the CJSW collision expressions of "The CJSW Collision Model" Section. For the electron effective masses, CJSW use values m * 1 ¼ 0:065 m e and m * 2 ¼ 0:222 m e which are different from the values we use for the other simulations described here. The GaAs diode and oscillator circuit are as described in "The GaAs Diode and Circuit" Section with the difference that for the CJSW simulations the doping profile is smooth across the device junctions.
Numerical results are presented by CJSW for the steady state of the device under a fixed voltage and for its oscillatory behaviour when coupled to a circuit. We obtain results which are in good agreement in all cases. What differences there are may be plausibly attributed to the different resolutions and numerical methods used. (The CJSW simulations use a third-order ENO scheme with Runge -Kutta time stepping.) We note that the differences in the effective masses and the doping profile mentioned above are observed to have significant effects on simulation results.

Simulations of Homogeneous Samples
The first of our comparisons between hydrodynamical models and Monte Carlo results considers the case of homogeneous GaAs samples when fixed electric fields are applied. For the hydrodynamical models, the evolution equations reduce to systems of ODEs which are easily integrated using the Runge -Kutta method to obtain steady solutions. For the Monte Carlo simulations, a simplified version of the code described in "Monte Carlo Simulations" Section is used in which a single particle is tracked and its behaviour timeaveraged to determine macroscopic properties at steady state. Figure 2 shows steady state results over a range of electric field strengths, with Monte Carlo results being contrasted against hydrodynamical solutions. The reduced hyperbolic transport model (5)  It may be observed from Fig. 2 that the new relaxation times (12a) -(12h) give a better agreement overall with the Monte Carlo results than the Tomizawa expressions do. This is to be expected, of course, since these relaxation times are derived using the same Monte Carlo code and the same physical assumptions. However, in neither case is the agreement between hydrodynamical and Monte Carlo results perfect, and this highlights the limitations of the relaxation time approach for modelling collision processes: more detailed models may be needed to accurately capture semiconductor behaviour, even in the homogeneous case.

Device Behaviour at Steady State
The GaAs diode described in "The GaAs Diode and Circuit" Section is evolved to steady state under a fixed applied voltage to provide another comparison between semiconductor models.
The initial state of the device is given by Eqs. (13) and (14), and the boundary conditions by Eq. (15). A voltage ]. It has previously been observed [16] for a silicon diode that the steady state solutions for the Bløtekjaer and the reduced hyperbolic transport models (which differ significantly only in their treatment of heat conduction) are very similar if both use the same collision terms. This has been found to be the case also for the GaAs diode (and also when the Anile -Pennisi transport model is considered), and consequently results for other combinations of transport and collision models are not presented in Fig. 3. For the CJSW collision model ["The CJSW Collision Model" Section] the numerical results differ substantially from the Monte Carlo calculations, as would be expected.
It may be seen from Fig. 3 that the hydrodynamical models give solutions which are broadly in agreement with the Monte Carlo results. However, the agreement is not perfect, and, perhaps surprisingly, the results obtained using the Tomizawa relaxation times may be judged to be closer to the Monte Carlo results than the ones obtained using relaxation times derived from the same Monte Carlo code. This suggests that the "effective" relaxation times for a device may be different from the relaxation times for a homogeneous sample.

Gunn Oscillations
For our final comparison between models, we consider the time-dependent behaviour of a GaAs diode subject to Gunn oscillations.
"The GaAs Diode and Circuit" Section describes the circuit used to determine the voltage across the diode. In our simulations, the circuit is engaged at the time when the GaAs diode is judged to have reached the steady state discussed in "Numerical Results" Section. (The applied voltage V ¼ 2 V for the steady state is the initial voltage used in the circuit model.) Figure 4 shows the voltage V(t ) across the device for Gunn oscillations in seven hydrodynamical models and also for a Monte Carlo simulation. It is apparent from the figure that significantly different results are produced by the different models: the oscillations can initially grow or decay, and may reach a steady amplitude or eventually vanish. The hydrodynamical results closest to the Monte Carlo ones are for the reduced hyperbolic transport model when the Tomizawa collision terms are used [see Sections "The Reduced Hyperbolic Transport Model" and "The Tomizawa Collision Model"], although the decay rates of the oscillations in the two cases are clearly different. (The small oscillations seen in the Monte Carlo results at late times are believed to be spurious: as the number of particles used in simulations increases, the amplitude of these oscillations decreases, and it is likely that they would not be present in a "converged" solution produced using sufficiently many particles.) It may be noted from Fig. 4 that there are significant differences between results for the three different transport models we use (the Bløtekjaer, Anile-Pennisi and reduced hyperbolic models). This is in contrast to the steady state results obtained for the device ["Device Behaviour at Steady State" Section] and reinforces the tentative conclusion of Ref. [16] (for a silicon diode) that the form of heat conduction used in a hydrodynamical model has more of an effect on time-dependent behaviour than on steady state behaviour.

CONCLUSIONS
We have compared numerical results for a variety of hydrodynamical models against Monte Carlo simulations for a GaAs diode both at steady state and when producing Gunn oscillations. The usefulness of hydrodynamical models in semiconductor device design is apparent: using modern numerical methods, accurate solutions for the models have been obtained in minutes, compared to Monte Carlo running times of possibly several days. However, while the hydrodynamical results are reasonably close to the Monte Carlo ones for the diode at steady state, only one of the hydrodynamical models tested produced Gunn oscillations similar to those of the Monte Carlo code.
The results of this study suggest (very tentatively) that the "reduced hyperbolic" transport model coupled to the Tomizawa collision model gives the most realistic results for GaAs devices. It is perhaps surprising that the Tomizawa relaxation expressions produce better results than alternative expressions derived specifically to be consistent with the Monte Carlo code in the homogeneous case. However, as has already been remarked, it could be that the "effective" relaxation times for the device are different from the homogeneous ones due to the influence of higher moments caused by, for example, jumps in the density. The Tomizawa expressions might, by chance, provide a better fit to these inhomogeneous relaxation rates. This calls for further work to be done in deriving semiconductor models within the moments framework, with higher moments being systematically included to capture complicated time-dependent behaviour such as Gunn oscillations.