UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Aspects of modeling III-V semiconductor devices Zhou, Haosheng 1991

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-ubc_1992_spring_zhou_haosheng.pdf [ 3.3MB ]
JSON: 831-1.0064782.json
JSON-LD: 831-1.0064782-ld.json
RDF/XML (Pretty): 831-1.0064782-rdf.xml
RDF/JSON: 831-1.0064782-rdf.json
Turtle: 831-1.0064782-turtle.txt
N-Triples: 831-1.0064782-rdf-ntriples.txt
Original Record: 831-1.0064782-source.json
Full Text

Full Text

ASPECTS OF MODELING Ill-V SEMICONDUCTOR DEVICESbyHAOSHENG ZHOUB.Sc., South China Institute of Technology, 1982M.Sc., The University of British Columbia, 1987A THESIS SUBMITTED IN PARTIAL FULFILMENT OFTHE REQUIREMENT FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinTHE FACULTY OF GRADUATE STUDIESDepartment of Electrical EngineeringWe accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIADecember 1991© Haosheng Zhou, 1991In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives. It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.________________________________Department of nee.vyThe University of British ColumbiaVancouver, CanadaDate 27) /9 2.DE-6 (2/88)AbstractThis thesis is concerned with novel modeling approaches to three GaAs-based devices,namely: the heterojunction bipolar transistor (HBT), the metal-semiconductor field-effecttransistor (MESFET) and the modulation-doped field-effect transistor (MODFET).In GaAs-based HBTs, the transit time of carriers across the collector-base space-chargeregion plays an important role in determining the performance of these devices at high frequency. In this thesis, a novel phenomenological approach is developed to modeling transport phenomena in GaAs-based HBTs. The model takes the dynamics of intervalley electrontransfer into account so that effects such as velocity overshoot can be studied. The results from this computationally-effective phenomenological model are in good agreement withthose from the time-consuming Monte Carlo approach. This suggests that the phenomenological model could be a useful tool for HBT design. The model is extended to treat thepractically-important case of operation at very high collector currents. It is revealed, perhaps for the first time, that the holes in a n-p-n device play a very significant role in thedecline of the device’s speed performance at high current levels.A one-dimensional model for MESFETs is developed in which the characterizing systemof ordinary differential equations is analyzed by the phase portrait method. The model yieldstwo criteria for the formation of a stationary dipole domain in MESFETs. One of thesecriteria pertains to the product of the doping density and the gate length, and produces alength-dependent value which is different from the constant value appropriate to travellingdomains and which has been widely-used in the past for stationary domain formation. Thenew criterion should prove useful in analysing the dynamic performance of MESFETs.A new method for obtaining the electron surface charge density in the potential well ofMODFETs is presented. The surface charge density is needed in calculations of the I—Vand C—V characteristics of the device. The new method uses the WKB method to solveSchrödinger’s equation for the bound states, and introduces a modified effective densityof states function to solve Poisson’s equation for the potential profile. The method yields11excellent agreement with the much more computationally—intensive full quantum approachto the solution.111Table of contentsPageAbstract.List of FiguresList of TablesAcknowledgementDedication.1. Introduction2. Phenomenological Model for Estimating2.1 Introduction2.2 The Model2.2.1 The Basic Equations2.2.2 The Method of Solution . .2.2.3 Values for Model Parameters2.3 Programs Used in the Calculations2.4 HBT Collector StructuresTransit Times in GaAs HBTs2.5.2 Dependence of Results on the Phenomenological Parameter r22.5.3 Comparison with Monte Carlo Simulations2.6 Conclusions3. Computation of Signal Delay Times for the Base-Collector Depletion Region ofGaAs-Based HBTs3.1 Introduction3.2 Method3.3 Programs Used in the Calculation3.4 Results and Discussioniv2.5 Results and Discussion2.5.1 Comparison of Collector Structures11• viii• XiiXfll.xiv..4•.5•.71113• 16• 161819• 222323242627Current Densities4.1 Introduction4.2 Extended Model4.2.1 Equations4.2.2 Estimation of Signal Delay Time4.2.3 The Method of Solution4.2.4 Boundary Conditions4.2.5 Programs Used in the Calculations4.3 Results and Discussion4.3.1 High Collector Doping Case4.3.2 Normal Collector Doping Case4.4 Conclusions5. Analysis of the Formation of Stationary5.1 Introduction5.2 One-Dimensional Analysis . .5.2.1 The Fundamental Equations5.2.2 Phase Portraits5.3 Criteria for Stationary Domain Formation5.4 Programs Used in the Calculation . .5.5 Results and Discussion5.6 Conclusion32343435353536404042424955585859596164666872737375783.5 Conclusions4. The Estimation of Transit Times in GaAs-Based HBTs Operating at HighDomains in GaAs MESFETs6. Use of the WKB Method to Calculate the Surface Charge Density in MODFETs6.1 Introduction6.2 WKB Approximation6.3 Model for the 2DEG in MODFETsVB.1 List of DELAY.B.2 List of DRIVHBTB.3 List of EXTRAAppendix C Listing of Programs UsedOperating at High CurrentsC.1 List of SUBFINE7C.2 List of EQUI2C.3 List of DELAYC.4 List of HBTFLC.5 List of FINITE1C.6 List of ONED3CC.7 List of ASSEMC.8 List of GAUSS• 7878• 79• 808081• 8183• 84• 8692• . 92• . 96• • 99• • • 99• . . 103• • • 106• . . 107• . . 107• . . 115• . • 124• . . 132• . . 143• • . 1461561576.3.1 The Potential Energy V(x)6.3.2 The Effective Density of States6.3.3 Determination of the Factor cr(Ef)6.4 Method of Solution6.4.1 Boundary Conditions6.4.2 The Program Used in the Calculation6.5 Results6.6 Conclusions7. Summary8. BibliographyAppendix A Listing of Programs Used to Estimate Transit Times in HBTsA.1 List of PULSHBTA.2 List of TRNTCAppendix B Listing of Programs Used to Estimate Signal Delay Times in HBTsin HBTsto Estimate Signal Delay TimesvAppendix D Listing of Programs Used to Estimate Stationary Domain Parametersfor MESFETs 160D.1 List of GAAS 1 160D.2 List of PEAKS 163D.3 List of SLOPE 165Appendix E Listing of Program Used to Estimate the Energy Levels in MODFETs . 167E.1 List of MODFET 167viiList of FiguresPageFig. 2.1 Schematic of the two conduction-band valleys in GaAs and the transitions ofelectrons between them 6Fig. 2.2 The basic HBT collector-base structure. The depletion region is the shadedregion and is taken to exist wholly in the collector region 8Fig. 2.3 (a) Mean electron velocity for electrons injected into a uniform field (75 kVcm—’) in a homogeneous semiconductor. The width in A of the profiles measured atvmax/e are: 50, 100, 362 and 456 for2(ns)=1, 10, 50, 100 respectively. (b) Figurereprinted from [29] 12Fig. 2.4 Collector-base structures and field profiles: (a) conventional arrangement, (b)inverted field configuration, (c) intrinsic collector device 14Fig. 2.5 Definitions for symbols. Pb is the p-type doping in the base, P refers to thecollector doping density and N is the n-type doping in the sub-collector. Wb is thedepletion region width in the base, W is that in the sub-collector. W’ is the widthof the p-type region, d is the width of the intrinsic region. Finally, E0 is the constantfield in the intrinsic region and E is the maximum field at the p-n junction. . 15Fig. 2.6 Electron profiles for the conventional structure of Fig. 2.4a. VBC = 0, r2 = 50ns. (a) lower valley electrons, (b) upper valley electrons. The line at x = 0 marksthe base/collector metallurgical junction. The line on the right marks the edge ofthe space charge region 16Fig. 2.7 Electron profiles for the inverted structure of Fig. 2.4b. VBC = 0, r2 = 50 ns.(a) lower valley electrons, (b) upper valley electrons. The line at x = 0 marks thebase/collector metallurgical junction. The line on the right marks the edge oT thespace charge region 17viiiFig. 2.8 Electron profiles for the intrinsic structure of Fig. 2.4c. VBC 0, r2 = 50 ns.(a) lower valley electrons, (b) upper valley electrons. The line at x = 0 marks thebase/collector metallurgical junction. The line on the right marks the edge of thespace charge region 19Fig. 2.9 Comparison of collector transit velocity (XSCR/TSCR) as estimated for a conventional structure by Yamauchi and Ishibashi [32], Rockett [8] and the presentwork 20Fig. 2.10 Comparison of TSCR as estimated for conventional and inverted field structuresby Katoh et al. [10] and the present work. The labels n, 0.5 etc. refer to the dopingdensity (in units of 1017 cm3)in the collector of a conventional structure (Fig. 2.4a),whereas p, 0.5 etc. refer to the collector doping density (in units of iO’7 cm3) inan inverted field structure (Fig. 2.4b). The collector width xc is either 2000 A or1500 A, as indicated 21Fig. 3.1 Total current JT, conduction current at base-collector interface (x = 0), andconduction current at the edge of the SCR region (x = xscR) 28Fig. 3.2 The transit times and signal delays for inverted field structure. The dash linesare transit times and signal vs. the width of collector (xc). The solid line is theratio of signal delay to transit time 29Fig. 3.3 The field profile and the mean velocity profile for the inverted field structurewith xc = 0.4 1um 30Fig. 3.4 The field profile and the mean velocity profile for the inverted field structurewith xc = 0.05 im 31Fig. 3.5 The field profile and the mean velocity profile for the conventional field structurewith XC = 0.25 urn 32Fig. 3.6 The transit time and signal delay for a conventional field structure. The dashlines are transit times and signal delays vs. the width of the collector (XC). Thesolid line is the ratio of signal delay to transit time 33ixFig. 4.1 Carrier density profiles and fields for different currents. (a)—(d) are forn1(x00)n2(xcc) = 1: 10; (e)—(h) are for fli(XC) : n2(XCC) = 10 : 1 45Fig. 4.2 Base push-out. Curves (a)—(d) are the profiles of holes from Fig. 4.1(a)—(d),respectively 49Fig. 4.3 Comparison of the old model and the extended model 50Fig. 4.4 Carrier density profiles and fields for normal doping profile 52Fig. 4.5 Profiles of holes. Curves a-e corresponding to Fig. 4.4a-e, respectively. . . . 55Fig. 4.6 TB + r0 vs. J with different properties of holes. Curve a is the normal case with= 400cm2V’s’ and v8 = 1 x io cms1; curve b is with = 2000cm2V1s’and v = 1 x iO cms’; curve c is with au,,, = 2000cm2V’s’ and v8 = 5 x iO cms’;and curve d is the one-valley case with all electrons in the lower valley 56Fig. 5.1 Device cross-section, and profiles for the electron concentration and the electricfield for the case of a MESFET with a stationary dipole domain at the drain end ofthe gate. Note the assumption of a constant depth of channel under the gate. . . 60Fig. 5.2 General form of the equilibrium electron velocity versus field relationship forGaAs. The electric field is taken as negative, in accordance with Fig. 5.la. . . . 61Fig. 5.3 The zero-value isoclines for equations (5-6) and (5-7). The short line segmentson the isolines have slope equal to 0 for II’ and oo for 00’; they indicate how thesolution curves cross the isocline 62Fig. 5.4 Phase portrait for the system of equations (5-6) and (5-7) 63Fig. 5.5 Comparison of the doping density-gate length product for stationary domainformation as computed from this work with the commonly-used criterion of NDL =1012 cm2. Also shown are experimental (solid symbols) and theoretical (opensymbols) data from devices for which evidence of stationary domain formation isclaimed 71xFig. 6.1 Band diagram of the A1GaAs/GaAs interface in a MODFET. The dash line,represents the triangular potential well approximation 74Fig. 6.2 Turning points. In a classical system, a particle with energy E cannot gobeyond these two points (a x b). In a quantum system, the wave functions canbe non-zero beyond these two points 76Fig. 6.3 n3 vs. Ef, from the present method, and from [75]. The surface density is in1O’2/cm 82Fig. 6.4 The function c(Ef). The solid line is the fit (equation (6-30)) to the computedpoints 83xiList of TablesTable 2.1 TSCR (ps) For Different Collector Structures (r2 = 50 ns) 18Table 2.2 TSCR (ps) For Inverted Structures at VBC = 0 19Table 3.1 Doping profiles. All densities are in cm3. Doping density in the base is2 x 1018 cm3 and that in the sub-collector is 5 x 1017 cm3 27Table 3.2 Transit and signal delay times for uniform field of 75kV/cm 28Table 4.1 Doping densities in cm3 for the device shown in Fig. 2.1 43Table 4.2 Signal delay (TB + TC) with different boundary conditions. Estimated fromthe “box” method. Current densities are x105 A/cm2 43Table 4.3 Signal delay (TB + TC) computed using the extended model with boundarycondition ni(xC0) :n2(xCC) 10: 1, and using the old model (Chapter 3). Currentdensities are x105 A/cm2 50Table 4.4 Doping densities in cm3 for the device shown in Fig. 2.1, for normal collectordoping case 51Table 4.5 Signal delay (TB + TC) calculated using the finite element method. Currentdensities are x iO A/cm2 55Table 4.6 Signal delay (TB + Tc) with different hole properties. Current densities arex105A/cm2,mobilities are in cm2V’s’ and velocities are in cms’ 57Table 5.1 Summary of Computed Results for Various Doping Densities 69Table 5.2 Comparison of Results for Lmin from the Linear and Nonlinear Analysis . . 70xiiAcknowledgementThe author would like to express his sincere appreciation to Dr. David L. Puifrey for hissupport and guidance throughout the course of this program.The author is also grateful to Dr. Matthew J. Yedlin for many useful discussions andsuggestions on solving partial differential equations. Especially, the author would like tothank Dr. Tom Tiedje for encouraging him to pursue studies in this field.xli’To my parentsxivChapter 1. IntroductionIt is well-known that the high mobility of electrons in GaAs and its related compounds(A1GaAs, for example) makes these materials of interest for the fabrication of very highfrequency devices [1]. In field-effect devices based on GaAs, a cut-off frequency (fT) ashigh as 117 GHz has been reported for modulation-doped devices (MODFETs) [2], and a110GHz oscillator has been realized using metal-semiconductor devices (MESFETs) [3]. InGaAs-based heterojunction bipolar transistors (HBTs), devices with a maximum oscillationfrequency (fmax) of over 100GHz have been reported [1].Technologically, the intrinsic property of high mobility needs to be complemented by anability to fabricate small dimension structures, if truly high frequency or high speed devicesare to be realized. In MESFETs, electron beam lithographic techniques have enabled thefabrication of devices with gate lengths as small as 0.03 im [4]. In MODFETs, it is notonly the gate length, but also the thicknesses of the channel layers that are important. Herethe use of molecular beam epitaxy (MBE) has proved useful. Devices with a 150 A channellayer and a 20 A spacer layer have been realized [5]. MBE is also a necessary technique forbuilding high performance HBTs. Here the base width is an important property, and valuesas low as 200 A have been reported [6].For the continuing development of these devices it is necessary to be able to modelthe devices. By modeling we mean the ability to predict the effect of the device’s physicalproperties (material composition and geometry) on either the overall electrical performanceor on some specific important feature of the devices. The fact that we are considering GaAsbased materials means that any models we build should take into account the multi-valleynature of the conduction band. Exchange of electrons between the lower two bands (I’and L) is a phenomenon which greatly affects device performance and calls for a differentmodeling approach than is used in silicon devices. The fact that we are dealing with very1small devices means that parasitic resistance and capacitance values are so low that transittime phenomena may determine the high frequency/speed performance limit of the devices.The two features of intervalley electron transfer and carrier transit time are both in evidence in HBTs. The fields in the base-collector space charge region are so high that electrontransfer to the lower mobility L valley is a problem. Also, the device lateral dimensionscan be made so small with modern self-aligned technologies [7], and the base width madeso small by MBE deposition technology, that the transit time of electrons through the highfield region of the base-collector depletion region can dominate the overall emitter-collectordelay time. It follows therefore, that accurate modelling of carrier transport in this junctionspace charge region is important.Accordingly, the main part of this thesis is concerned with the development of a novelapproach to modeling transport phenomena in an environment where intervalley electrontransfers are likely. The main goal of the research is to produce a model which gives comparable accuracy to the time-consuming Monte Carlo approach [8] [9] [10], but which is sosignificantly faster that it can be used as an efficient device design tool. The model is phenomenological in nature and is described in Chapter 2. The use of the model to predict signaldelay times, as opposed to carrier transit times, is described in Chapter 3. An interestingresult of this work is that the widely-used value of 1/2 for the ratio of the delay time to thetransit time is found to be not appropriate for GaAs HBTs. To extend the phenomenologicalmethod so that the practically-important case of high-current performance can be modeled,it is necessary to forsake the depletion approximation for estimating the electrical field andto include the effect of the space charge due to the mobile electrons and holes. Our attemptsto do this, and yet still maintain the fast execution time of the phenomenological approach,are described in Chapter 4. The results presented reveal the interesting fact that it is theholes, rather than the electrons, in n-p-n devices which are responsible for the reduced highspeed performance at high collector current levels.2Turning now to modeling aspects of small and high speed GaAs FET devices, Chapter 5deals with stationary domain formation in MESFETs and Chapter 6 is concerned with thesurface electron charge in the potential well of MODFETs. The stationary domain whichcan form at the drain side of the gate in MESFETs [12][13][14][15][16][17], is particularlyimportant in small devices as the coupling capacitance it introduces tends to become thedominant parasitic capacitance. The novelty of the work presented in Chapter 5 is the use ofthe phase portrait approach to represent the solutions to the system of ordinary differentialequations which describe the device. This approach, which gives a useful graphical picture ofthe interdependence of local carrier concentration and field, yields a criterion for stationarydomain formation which is different from that pertaining to traveling domains [13]. Finally,in the MODFET modeling described in Chapter 6, we seek to solve the 1-D Schrödingerequation for the bound states in the well in a manner which is accurate, yet computationallyefficient. The approach uses the WKB method to solve the Schrödinger equation, andintroduces a modified effective density of states function to solve Poisson’s equation forthe potential profile needed in the Schrödinger equation. The new approach yields resultswhich agree very well with those from a much more computationally-intensive full quantum,self-consistent approach.3Chapter 2. Phenomenological Modelfor Estimating Transit Times in GaAs HBTs2.1 IntroductionIt is well-known [8][9][1O] that the transit time of carriers across the base-collector spacecharge region plays an important role in the microwave performance of HBTs. The relationbetween the common-emitter unity-gain cut-off frequency fT and the various transit timesand delay times is [8][20]:fT=[27r(rECBE+rB+rc+CBc(rE+REE+Rc))] (2—1)where r is the incremental emitter resistance, CBE is the total base-emitter capacitance,TB is the base transit time, rc is the base-collector depletion region signal delay time, CBCis the base-collector depletion capacitance, and REE and Rc are the emitter and collectorseries resistances, respectively. The signal delay r0 is related to the transit time TSCR ofcarriers across the base-collector space charge region (see Chapter 3). TSCR plays a moreimportant role than TB m the determination of a modern HBT’s microwave performancebecause of the narrowness of the base layer. However, there has been no direct measurementof transit time so far. All estimations of transit time are obtained by analysis of the cut-offfrequency. Theoretical studies, on the other hand, are mostly restricted to Monte Carlostudies [8][9][1Oj. It is so because in Monte Carlo studies, the velocity of individual carrierscan be traced and recorded, and thus transit time can be easily computed. Although havingshown its ability to estimate transit times, Monte Carlo methods usually take tremendousamounts of CPU time and are usually performed with a super computer. The objective ofthe work described in this chapter is to present and validate a new model of transport inGaAs-like materials from which the transit times can be estimated accurately, but with farless computational effort than is needed for Monte Carlo calculation.4The method presented here takes a phenomenological approach to carrier transport intwo-valley semiconductors by incorporating transition rates for intervalley electron transfersinto the continuity equations for upper and lower valley electrons. The transition rates arecharacterized by phenomenological transfer times. The two transfer times for upward anddownward transitions are shown to be related via the electric field in the device. Thusby specifying one transfer time and the electric field distribution (i.e., the collector-baseconfiguration) in a device, transient analysis of the system equations can yield the transittimes of carriers through specified regions of the device.2.2 The Model2.2.1 The Basic EquationsConsider a two-valley conduction band model for GaAs comprising a lower valley (F) andan upper valley (L or X, see Fig. 2.1). Let the excess electron concentrations in these twovalleys be n1 and n2, respectively. We propose the following modified continuity equations:On1 n2 n1—V. (D1Vn)+ V. (nipiV) = — — —u 2 i (2—2)On2 fi 2— V. (D2Vn)+ V (n22V) = — — —(IL T1 T2where D is the diffusivity, z the mobility, andY the electric potential. The transfer timesr1 and r2 are the phenomenological features of the model. Clearly, the physical meaning ofthe right-hand side of the expressions in (2-2) is that n1/r represents the rate of changein excess electron population due to transfer from the lower to the upper valley, and n2/rrepresents the rate of change in excess electron concentration due to transfer from the uppervalley to the lower valley (see Fig. 2.1). Note that the recombination-generation term, whichis usually found in the continuity equation for minority carriers, is omitted from (2-2). Thisis permissible for the heterojunction bipolar transistor transit-time studies of interest herebecause the base is so narrow and traversal of the base-collector space charge region is sorapid. Equation (2-2) is a balance equation of the general form found in many instances5ri2milL1 + fl2IL1= I IUpper valley (L or X)T2 EI-i-]T1(F) kFig. 2.1 Schematic of the two conduction-hand valleys in GaAs arid the transitionsof electrons between them.of irreversible thermodynamics [21]. In such an approach, r1 and T2 are phenomenologicalparameters and, for that reason, are referred to here as transfer times.To establish a link between r1 and let us consider the case of steady-state conditions,no concentration gradients in n1 and n2 and constant electric field. Under these circumstances, (2-2) reduces to(2—3)T2 T1It is under the conditions just specified that a measurement of drift velocity could be usefullymade. For two-valley semiconductors we can write(2—4)For low-field conditions, generally i > n2, ii can be taken as a constant and the value of/-2 is unimportant. At high fields we expect n2 >> n1 and the drift velocity to saturate: the6valueof is not important in this case. Summarizing these asymptotic conditions we haveP1L = constant (2—5)IL2H =where v3 is the saturation velocity and the subscripts L and H refer to the cases of low andhigh field, respectively. Combining (2-3)-(2-5) yieldsLlLIEI + ()VdV3[i+()](2-6)This equation has a similar form to the widely-used empirical relation for the steady-statedrift velocity in GaAs [15], namelyE8 \E0)Vd—V8 —1 + (0)where E0 is a material constant and E3 is related to the low-field mobility by means ofp = v3/E [15]. Taking advantage of this expression for mobility, the comparison of (2-6)and (2-7) givesr2(IEI)4 (2—8)(2—9)I’1LEquation (2-8) is a key equation in the model as it links the two phenomenological parameters r1 and r2. Detailed calculations of the actual intervalley transition times [22], [23]indicate that the transfer of electrons from the upper valley to the lower valley is practically independent of energy. Thus we make r2 a model parameter, which can be assigned aconstant value, and we allow r1 to be field-dependent and calculated via (2-8).2.2.2 The Method of SolutionThe basic device structure studied in this chapter is shown in Fig. 2.2. By specifyingr2 and the fields in the base and collector regions, r1(x) can be computed via (2-8), and theupper and lower valley electron velocities (p2V andp1V) can be computed as described7+B::E OLLETO sub Cn0 XSCR XC XCCFig. 2.2 The basic HBT collector-base structure. The depletion region is the shadedregion and is taken to exist wholly in the collector region.in Section 2.2.3. Then, subject to particular boundary and initial conditions, (2-2) can bediscretized and solved by numerical integration to yield the temporal and spatial variationof ri and n2. We are interested in the times taken for a pulse of electrons injected into thebase to cross the base (TB) and the base-collector space charge region (TSCR).The boundary conditions aren — f (no/2)(1 — cos27rt/To), 0 t To1 B, o, t>T0fl2FXB,t) =0 (2—10)ni(xc,t) = = 0where T0 is taken to be 8 ps and n0 is the peak value of the pulse of injected electrons.Equations (2-2) are solved subject to the above boundary conditions by using the program MOL1D [24] (see next subsection). Typically, a complete run for a value of r2 50 nstakes about 300s of CPU time on an IBM 3081 machine. This computation time, althoughobviously orders of magnitude shorter than that required for Monte Carlo calculations, couldbe reduced considerably if it were possible to diagonalize the spatial discretization matrix.Also, a further reduction in CPU time could be achieved by using either the “box” or “finiteelement” algorithms described in Chapter 4.SPrior to the injection of the excess electrons into the base there are no excess carriersanywhere in the device, i.e., the initial conditions aren1(x,0) = 0(2—11)n2(x,0) 0The program output is in the form of plots of n1, n2, and total excess electron concentration (n = n1 +n2) versus time and distance (see Section 2.5.1), and of printouts at 0.1-0.2ps intervals of n(t) at x = 0 (the base-collector interface) and at XSCR .(the edge of the spacecharge region in the collector). The transit times are computed fromj’tJdtfJdt(2—12)where the current density J(0,t) is used to estimate TB and the current density J(XSCR,t)is used to estimated TB + TsCR.2.2.3. Values for Model ParametersIn each of the valleys, the electron drift velocity is assumed to saturate at high fieldsand, at low and intermediate fields, to follow a relationship of the formvi = LtjE(1 — cE2) V(2—13)E 1.5—-, i = 1,2I’iLwhere -4 1RiL’2and v is the saturation velocity in the ith valley. This form of v(E) has been used for othersemiconductors [25]. Here we takev = 108 cm• s—i, v25 = iü cm• s1IL1L = 6000 cm2 V’ s (2 — 14)—*/—m1/rn2) I’iL•For the lower valley the effective mass m is taken to be 0.063 m0 [26]. The upper valleyis considered to be the X-valley, for which the conductivity effective mass m may be taken9as 0.27 m0 [26]. Relating the low-field p’s as in (2-14) should be reasonable as the reciprocalsof the relaxation times for the dominant scattering mechanisms (ionized impurity and polaroptical phonon) in GaAs can both be represented by expressions which are dependent onthe square root of the effective mass [27]. The electron diffusivities are computed from theappropriate mobilities using Einstein’s relation.Turning now to the transfer times, T1 and r2 are related through (2-8), in which we takeE0 = 3.4 kV cm’ [15]. As detailed calculations indicate [23][28j, the transfer of electronsfrom the upper valley to the lower valley is practically independent of energy. Thus it isappropriate to take T2 as a constant. It follows from (2-8) that i will decrease with increasingelectric field. This is reasonable as the upward transfer time should be strongly dependenton the energy of the electrons in the lower valley, since these electrons have to overcome anenergy barrier to be excited into the upper valley. To establish a realistic value for i2 weconsidered the case of injection of excess electrons into a homogeneous semiconductor underuniform field conditions, and observed the resulting velocity-distance profile. The field in thesemiconductor was taken to be 75 kV cm1,which was high enough towarrant the neglectof diffusion and, therefore, the use of (2-4) to estimate the average velocity. The individualvalley velocities were calculated from (2-13), and n1 and n2 were computed by the program.The results are shown in Fig. 2.3. As our model is not subject to constraints imposed bymomentum and energy considerations, the electrons attain a very high velocity immediatelyon experiencing the high field. Thus the velocity-distance profile has amuch steeper frontthan predicted by Monte Carlo calculations (see case 1, Fig. 2.3b). The distance over whichvelocity overshoot extends depends strongly on 2, and we choose as areasonable value forr2 that which gives the same width of profile (measured at (Vmax — v,)/e) as the Mojate Carlocalculations, i.e., 360A. The result is r2 = 50 ns, and we use this value, unless otherwisestated, in the calculations presented in the following sections.102.3 Programs Used in the CalculationsThe main program PULSHBT is listed in Appendix A, List 1. In this program, the subroutineMOL1D is used to solve the partial differential equations (2-2). The details of using MOL1Dare described in [24]. It is a general FORTRAN IV subroutine for numerically solving partialdifferential equations (PDEs) in one spatial dimension. The method of lines is used wherebythe space interval is discretized and spatial derivatives are approximated by second-, fourth-or sixth-order symmetric differences etc.Together with MOL1D, there are three user-supplied subroutines. BNDRY is for theboundary conditions. In our case, these are (2-10). PDE defines the partial differentialequations, as its name indicates. Finally, FUNC gives the flux functions. We chose the fluxfunctions to be the current components‘9nif1(x,t) = —D1---— +v1(---)nOX OX 2 158fl2 —f2(x,t) —D2----— +v2(---)nOX OXOther subroutines in PULSHBT are: FIELD, VDRI1, VDRI2. The subroutine FIELD givesthe field at given x. The electric field profile is defined in the input file as a piecewisefunction.The subroutines VDRI1 and VDRI2 give, respectively, the drift velocities for the lower-and upper-valley excess electrons. The variable EPS contains the criterion for error control.It was set to be 10_b. A lower value will cause the program to take extra CPU time to meetthe criterion while a larger value of EPS will increase the output error.Note that the factor given by (2-8) is very big under some circumstances, usually inhigh field regions. For instance, for El = 100 kV/cm, r2/r1 = 6.48 x iO. This causes theprogram to make extra efforts to meet the error criterion, and is the main reason for themoderate CPU time. For this reason it was decided to limit r2/r1 to values less than 106.Thus, with r2 = 50 ns, the minimum value of r1 is 0.O5ps. This is equivalent to the value atIEI = 108 kV/cm, i.e., r1 will be the same for any field higher than this value.11I I9.00j8.00G)a 7.00U6.00>a>,5.00C)C 4.00 —2100 flS€Ca) S.> •3.00 —nseca •2.00 —10 nsec1.001 nsec__________I0.00 100.00 200.00 300.00distance, nm(a) CASEL4.02 3.0a)2.0a)I :: . . .0.0 Th0 1500 225,0 3000 3750X. . .AngstrOrnS(b)Fig 2.3 (a) Mean electron velocity for electrons injected into a uniform field (75 kVcm1)in a homogeneous semiconductor. The width in A of the profiles measuredat vmar/e are: 50, 100, 362 and 456 for 2(ns)=l, 10, 50, 100 respectively.(b) Figure reprinted from [29].12Also included in Appendix A, List 1 is a sample of the input file. The output file is notincluded because of its huge size (> 100kB).Appendix A, List 2 is the program TRNTC. This program is used to process the outputfrom PULSHBT to estimate the values of r from (2-12) for specified values of x. Thevalues of ‘r at x = 0 and x = XSCR were used to estimate the transit time (TSCR = T(x =XSCR) — T(X = 0)).24 HBT Collector StructuresThe three different types of device studied are shown in Fig. 2.4. The unconventionalcollector structures shown in Fig. 2.4b and Fig. 2.4c were chosen forstudy because theyhave been predicted to give very good speed performance [9J[1O][29]. Following Hu et al.[9J,who studied the switching speed of these devices, we call the structure in Fig. 2.4b aninverted device and in Fig. 2.4c an intrinsic collector device. The doping densities indicatedin Fig. 2.4 were also taken from Hu et al.[9].The field distributions, which are needed to estimate r1 (see equation (2-8)), were computed using the depletion approximation. In the intrinsic collector device, following theresults of numerical simulations [82j, we assume that the field in the intrinsic region is a constant. The field then has the profile shown by Fig. 2.5. The parameters to be determinedare E and E0. The given parameters are P, F, and N, where P6 is thep-type doping in thebase, P refers to the collector doping density and N is the n-type doping in the sub-collector.Other unknowns are Wb and W. W6 is the depletion region width inthe base; W is that inthe sub-collector. Finally, d and W’ are also given parameters. d is thewidth of the intrinsicregion and W’ is the width of the p region in the collector.Assume the built-in potential is 4, then from Fig. 2.5, one can write,(2—16)13a (urn)0.1 0.2 0.3 0.40.D2um181 a 10p.Base Collector sub C0.lum 0.23urn 0.lSum18 172 a 10 IntrinsiC 5 a 10p•a (urn)Fig. 2.4 Collector-base structures and field profiles: (a) conventional arrangement,(b) inverted field configuration, (c) intrinsic collector device.Base Collector sub C0.lum 0.25urn 0.lSum18 16 172x10 5x10 Sa 10p•E (ky/cm)0Base Collector sub C0.lum 0.25um 0.l5um18 16 172x10 3x10 5a10p• n•B (ky/cm)-100 /(a)a (urn)-100(b)E (ky/cm)-200400(C)14Fig. 2.5 Definitions for symbols. Pb is the p-type doping in the base, P refers to thecollector doping density and N is the n-type doping in the sub-collector. Wb isthe depletion region width in the base, W is that in the sub-collector. W’ isthe width of the p-type region, d is the width of the intrinsic region. Finally, E0is the constant field in the intrinsic region and E is the maximum field at thep-n junction.The system should be neutral so one also has1’V N = P W’ + Pb H’b (2 17)Other relations areEqN—(2—18)(E-E0) qPWIAfter some algebra one obtains(i+) Wb+2 [(d+w’)+wt] Wb+ (i+) ‘W’2 = (2—19)Solve the above equation for Wb. From (2-18), it follows that E and E0 can be obtained.E015Ij)C)--Ci))t--±CDCD—(C)C))C)(<-.-H -CDCDCD-i)IIHCfr)Ci)IICDcS.;ç!jCl)s—.CD—..C))C))--_pCDCD±C —CDCDCDQ-c,°c-QCDCDCl)C-q-CC))-II-.-2-CDC))CHoCDCD-CDDICD-X2Ci)---.sCD,o)-,CDCDCDCC.C)CD-÷CDCD‘—C-±C,-c-.-c-tCDCDCDC‘•C-CDR<CDC-cCDCCD-4-—-005CD2C)CD—‘—— CCD—.o—.2—CC-)C-)—C)H5.,-CC)CC)CD2CDCD-CC-—c*.CDCl)-c-i-CD—5c---CD,CD‘-.C)-C-Cl)C))CDCi)‘Cl)—.fl0<C—.CDCDCD.-••CDCDCC—.-Cd)—5o÷CDII.CD0 - -Cl)—.(‘3CD©-CDHC)C)-C--C-.5oCDCDIICDCD-—‘CD2-C)CCDCl)CDC))aC) 0 CD C)z 0 CDCD-qoc0CD 0—.ICi)-o‘CD 0CDCDCi)-.CDC?) C)o—.CD0 c-CDCDCDCDCDC) CD‘ci)c.DCl))i-$•C)CDCi)c<CDCi)—0‘CccCD 0•—.CJC CcCDCDsCDc-—N—.—CDCD- CDCDc-cCi)c—cC)-CDCDCi)--CD--.—-—•Ci)Ci)occ—.Ci)oCi)-- CDCD)—Ci)Cl)0c—c CD—.CD—.c--cCi)p.,p.,o—CDc.-<CDCi) o c--c-c-40C)CDo,c-c-c--cCi)--CDCD><‘-cCD c-c-Qc--cc-c--CDCDCD-CDp,p.,0CDC)CDCDC)c-c---—p.,-CDo•-P)CDc—cCl)CCC)0c-cc--c-cc-cCi)CD(qc-c-—CDCDz000c-c-CD0CDc-0-CDc-c-CD-Ci)0-<0c-c-Q)CD‘—-CDCDCi)CDCD--c--c-<c-ac-iCDCD<0c-c-c—c-Ci)C)c<Ci)c-c-- CD00Ci)r;5-0c--c-cc--c-—C)pioc--c--C)P)0-c-tjCD=CD---c-C0CDCl)—--CD-.D(00-cc--i-‘—cCD—.0c-cc--c-c0.CD0-c-cc-c-CDHp.’c--cCDC)C)Ci)Ci)Ci)F‘c-rjCDIIq-—.c-c-c0)c-c- tr0—.P’CD_——c-c-Cl)•cCD_Ci)p.,—CDc-c-CD-—C)CCD0 0p.,CDc-0Ci)—c-P’0C)C) - CDCDCDCi)o —-CDc-c-•—c--c- CDp.,c-—Ci)CD0CD••HCD::‘-C)—.CD 00HCD c-I--CDp.,I-tIIr Ci)ç;\athe field in the inverted structure has an appreciable magnitude ( 110 kV. cm’) at thecollector-base boundary. This results in significant transfer of electrons to the upper valleyat the beginning of the depletion region, and a corresponding increasein TSCR• The intrinsicstructure is even more adversely affected on account of the very high field ( 300 kV cm1at VBC = —3 V) at the p/n collector boundary. For the conventional structure mostof theelectrons are in the upper valley at VBC = 0 V, so the transfer ofthe remainder, broughtTable 2.1TSCR (ps) For Different Collector Structures (r2 = 50 us).VBC=OV VBC=—3VConventional 1.35 2.19Inverted 0.78 2.37Intrinsic 0.57 2.69about by the reduction in r1 due to the increased field, is not asmarked as in the other twocases. The fact that VBC = —3 V serves to put most of the electrons into the upper valleyin all three cases, means that the width of the space-charge region essentially determinesTSCR. The respective values of XSCR for the conventional, inverted, and intrinsic structuresat VBC = —3 V are: 0.26, 0.28, and 0.32 ,um.2.5.2 Dependence of Results on the phenomenological parameterr2The phenomenological nature of the method presented here affords some latitude in theselection of the key model parameter r2. Clearly one should choose r2 so that the predictionsof the model are realistic. The results presented in the next section confirm that the valuefor ‘r2 of around 50 ns, as chosen in Section 2.2.3, is appropriate. This time is long comparedto the transit times predicted for the devices considered in this thesis. This implies that onceelectrons in the base-collector space charge region (S CR) are scattered into theupjer valley,they will not have the opportunity to relax back into the lower valley. The significanceof the actual value of r2 is, therefore, in its effect on r1 via (2-8). For r2 =50 ns and= 50 kV. cm’, r1 becomes 1 x 10_12 s. This appears to be areasonable order of18C)tCCD o)Cl)C!) CCD C •CDCD Cl)C-+-,-CC)) CD Cl) H Cl) C C) C CD C) C C C)) CD—.CCC-4-b-CD-q-CD.C••eCDCCl)—.CCDCDo—•-i,z—.Cl)aC—.-<—Cl)C--CDC<0Cl)•CC))CD<Ci)CD——cc-CCDCl) CD—Cc’- ‘-C))C))t)—.CCDCD-—C—t-CD -CDCDCD<)•)——<CD:CD‘CDCD Cl)<CDCD9)C9)-C00CD9)CD—.Cl)CCD—l•C-Cl)•C)C)CCDCi)9)C))9)C))C•CD—.CDx—•-9)czC) Cl) Iii C) CD Cl) 9)qC+—9)C±CD‘—0’)CD9)—C)CCD—C)——C-C)9)CCDCDCD-C)Cl)fr•C9’Co=C-•-CDqi_-.C)——Cl)CDC)-C’Cl)o——•_C)CDCD C)__C-C-CCl).CDC CDaq CDC-9,II9’ IICD CD, 1Cl)CDCO.çcCl)C.)100 —3BI I —I 01 .4001.301.201.10>,.2-H0 1.000‘-4S> 0.900S0.800.700.600.00 1.00 2.00 3.00 4.00V VCsFig. 2.9 Comparison of collector transit velocity (XSC.R/TSCR) as estimated for aconventional structure by Yamauchi and Ishibashi [32], Rockett [8] and thepresent work.structures, Monte Carlo simulation techniques have been used [8], [10], [29]. In this sectionwe compare the results from [8] and [10], in which values of TSCR are quoted, with theestimates of TSCR from the phenomenological model proposed here.Rockett [8] has presented Monte Carlo results for a conventional device (Fig. 2.4(a))with a built-in base field (12.5 kV. cm’) and a narrow basewidth (0.12 sum). A device withthese properties was studied experimentally by Yamauchi and Ishibashi [32]. Based on thedependence of fT on bias, the latter workers deduced that the mean electron velocity in theSCR varied substantially with reverse bias, indicating the presence of significant velocityovershoot. Contrarily, the Monte Carlo simulations suggested a mean SCR transit velocityof iO cm s which varied very little with VEC. Fig. 2.9 displays the two sets of resultsgraphically, along with our own simulations of the structure. Our results suggest that velocityovershoot is significant only at bias voltages close to 0 V. As an increasingly large reverse20ci— ref [32]ref [8] 1ci 7 =50 ns, V21.0 x 10 cm/sec>< T =50 os, v..=0.6 x 107cm/sec2cici ciB—2.0U,IU,Ci 1.00000.0Fig. 2.10 Comparison of TscR as estimated for conventional and inverted field structures by Katoh et al. [10] and the present work. The labels n, 0.5 etc. referto the doping density (in units of iO’ cm3) in the collector of a conventionalstructure (Fig. 2.4a), whereas p. 0.5 etc. refer to the collector doping density(in units of iOfl’ cm—3) in an inverted field structure (Fig. 2.4b). The collectorwidth xc is either 2000 A or 1500 A, as indicated.bias is applied, the associated large fields and depletion layer widths ensure that most ofthe electrons are scattered into the upper valley, and so the mean velocity approaches thesaturation velocity v23. Results are shown forv2’values of 1 x iOT and 6 x 106 cm s1, i.e.,the values used in [8] and [32], respectively.A detailed evaluation by Monte Carlo methods of inverted collector structures has beencarried out by Katoh et al. [10]. They investigated different doping densities and widths forthe lightly doped p-region in the inverted structure shown in Fig. 2.4(b). For comparisonpurposes they also studied the conventional structure (see Fig. 2.4a), with different dopingdensities for the n-collector. Their results, along with the results from our work, are summarized in Fig. 2.10. The agreement in absolute values of rSCR from the two methods is good.Even more gratifying is the fact that both methods of analysis predict the same trends inTSCR for both the conventional and inverted structures.ri, 0.5 n, 1.0 n, 2.0 p, 0.5 p 1.0 D, 1.02000 2000 2000 2000 2000 1500212.6 ConclusionsThe two major conclusions to be drawn from the work presented in this chapter are asfollows:1) The notion of representing the intervalley transfer times as phenomenological parametersand relating them via the electric field is reasonable. The evidence for this is the goodagreement between Monte Carlo simulations and the predictions of the proposed model.2) The proposed phenomenological model for estimating transit times appears to havepromise as an HBT design tool. This conclusion is supported by the closeness of theresults to Monte Carlo predictions and by the greatly reduced computational resourcesneeded to implement the model.In addition to the above, we would like to suggest that our two-valley phenomenologicalmodel can, perhaps, provide a theoretical base for the well-established, but empiricallyderived, relation (2-7) for the steady-state drift velocity in GaAs.22Chapter 3. Computation of Signal Delay Times for the Base-CollectorDepletion Region of GaAs-Based HBTs3.1 IntroductionIn scaled HBTs, the dominant contribution to the emitter-collector signal delay time (seeequation (2-1)), and hence to the common-emitter, short-circuit, unity gain frequency fT, isoften the signal delay time ro through the base-collector depletion region [30][32}. Usually TOis taken to be one-half of TSCR, the transit time of carriers through the space charge region.For a derivation of this relation see [33] [34] [35] [36] .Recently, Laux and Lee [38] have pointed-out that this relationship is only correct when the carrier velocity in the space charge region isuniform. In AlGaAs/GaAs HBTs, because of velocity overshoot and a rapidly varying field,the velocity profile is likely to be highly non-uniform [8], so the assumption of Tc = TSCR/2will give erroneous results. In [38] it is shown that when the velocity profile v(x) is taken tobe a simple two-step function, ro is actually less than TSCR/2 in the case when the highervelocity is nearer to the base edge of the collector-base space charge region. Such a profile-is an approximation to the situation in real A1GaAs/GaAs HBTs, suggesting that the highfrequency performance of these devices should actually be better than that predicted usingTO = TSOR/2. A similar conclusion, again using a binary model for v(x), has been reachedby Ishibashi [37].While the work cited in references [37][38] demonstrates the importance of consideringthe non-uniformity of v(x) in estimating T0, a more detailed calculation, taking into accountthe actual v(x) profile is needed for an accurate computation of Tc in real devices. Thischapter describes such an investigation for GaAs-based HBTs in which the magnitude ofthe field in the collector portion of the base-collector depletion region is taken to be eitheruniform, decreasing (conventional structure) or increasing (inverted-field structure). The23v(x) profile is computed from the E(x) profile (calculated as described in Section 2.4),taking into account inter-valley transfer.3.2 MethodThe basic base-collector structure of the devices studied is the same as shown in Fig. 2.2.As described below, for a comparison of signal delay and transit times, we are interested inthe response of the system to a sinusoidal injection of electrons into the base. We assumethat initially all the electrons are in the lower valley, i.e.,nl(—xB, t) = no(1 — cos2irt/T)fl2(XB,t)0 (3—1)ni(xcc,t) = fl2(X,t) = 0where, as before, n and n2 are the lower and upper valley excess electron concentrations,respectively, n0 is the amplitude of the waveform and T, the period, is taken to be muchlonger than the base-collector transit time TSGR. As TSCR is typically of the order of severalpicoseconds, the specification of T > TSCR, in order to obtain a stable solution to thesystem equations, does not present a serious limitation to the observation of high-frequencyphenomena. For example, here we use T = 8 ps, which corresponds to a frequency of 125Hz.The phenomenological approach, described in detail in the previous chapter, allows computation of n1 and n2 as functions of distance and time. From these concentrations theconduction current Jc (due to drift and diffusion) can be computed at any place and atany instant. Because of the sinusoidal excitation, the conduction current is sinusoidal also.Thus, observation of the phase difference between J(0, t) and JC(XSCR, t) is a convenientway of estimating directly the transit time TSCR.Similarly, the collector signal delay time can be measured from the phase differencebetween Jc(0, t) and the total current. The latter can be computed from [38]{39]1JT j Jc(x,t)dx (3—2)XSCR 024Therefore, by displaying the time dependence of Jc(x = 0), Jc(x = xsCR) and .JT, thetransit time TSCR and the signal delay time Tc can be independently evaluated.The method proposed here differs considerably from that used by Laux and Lee [38].Their method computed r from the time delay in the collector transport factor , whichcan be written as [38]1 (XSCR . cx dx’’\J dxexp(—jw , J (3—3)XSCRJO .io v(x)jThe derivation of (3-3) relies on being able to express the conduction current at any planex in the depletion region as( /__JcJ0eXp JW t_jv(x’)where J0 exp(jwt) is the current modulation at the plane of injection, i.e., at the base sideof the depletion region. Note that (3-4) is only strictly true in media for which a generalsolution to the continuity equation can be given bydx’Pv=f(t_jV(X’))(3_5)where p is the charge density. In the case of semiconductor materials with multiple conduction bands, p and v each comprise components from the various bands. The coupled natureof the particle density equations for each of the bands does not allow a simple solution of theform given in (3-5) to be written. This fact does not invalidate the important revelation ofLaux and Lee [38] concerning the deviation of (Tc/TSCR) from 0.5 when v(x) is non-uniform.However, it does emphasize that, for the purpose of ascertaining precise values of (TC/TSCR)for particular Ill-V materials and device structures, numerical methods of solution, such asthe one proposed here, are in order.253.3 Programs Used in the CalculationThe programs used which are particular to this chapter are given in Appendix B. List 1 isthe program DELAY. It is much like PULSHBT (see Appendix A) except that the boundaryconditions now are (3-1). To save CPU time and cost, the period T was chosen to be 8 Ps.The program actually simulates the transient process for 16 ps, i.e., two periods.The output of DELAY is in the same format as that of PULSHBT. It is processed bythe program DRIVHBT listed in List 2. DRIVHBT gives the currents at x = 0, x = XSCR,and the total current from (3-2).Appendix B, List 3 is the program EXTRA. This program processes the output fromDRIVHBT to estimate the time when the total excess electron density n = ni + n2 reachespeaks at x = 0 and x = XSCR. The difference of times when n reaches its maxima at x = 0and x = XSCR gives the estimate of the transit time r0. The method used is described inthe following.From the output of DRIVHBT, one can find three values of current such that J(t1) <J(t2) and J(t3) < J(t2), where < t < t3. Then the Lagrange three point interpolation-formula is used to approximate the value of J(t) at ti <t <t3:JQ2 + pat) = ‘J(t1)+ (1 —p2)J(t + P2+1)J3) (3—6)where Ij < 1 and Lt = t2 — = — t2. Standard procedures then can be applied to (3-6)to find when J(t) reaches its maximum.263.4 Results and DiscussionIn the work described in this chapter, three profiles for the base-collector depletion regionwere studied, namely: uniform, conventional and inverted. In the first case the field was setat 75 kV/cm. In the other two instances the field was calculated for a reverse-bias of 3Vusing the depletion approximation and the doping densities listed in Table 3.1. In all casesthe material was taken to be GaAs, such as would be the case for an A1GaAs/GaAs HBT,for example.Table 3.1Doping profiles. All densities are in cm3. Doping density in the base is 2 x 1018 cm3and that in the sub-collector is 5 x i0’ cm3.Collector struc. Doping type Doping densityConventional n 5 x 1OInverted 3 x 1016Considering first the uniform field case, note that, because of velocity overshoot, thevelocity of the electrons does not remain constant as they traverse the depletion region.This can be confirmed by refering to Fig. 2.3. For the case of E = 75 kV/cm, and using= 50 ns as in Chapter 2, equation (2-8) yields r1 = 0.2 ps. This is the period of timefor which entering lower-valley electrons remain in this valley before scattering to the uppervalley. As we are using a value for r2 which is far longer than TSCR, once electrons transfer tothe upper valley, they remain there. Thus, the uniform-field situation gives rise to a velocityprofile in which v(x) is high near the base side of the depletion region, and lower towards thecollector side. From the general result of [38], such a profile should result in TC/TSCR < 0.5.Our specific results confirm this (see Table 3.2), and suggest that for the field strengthconsidered here, a value of 1/3 is appropriate for TC/TSCR. The results in Table 3.2 were27400.00300.00Co200.00a,a, 3.1 Total current.JT, conduction current at base-collector interface (x = 0), andconduction current at the edge of the SCR region (x = XSCR).Table 3.2Transit and signal delay times for uniform field of 75kV/cm.XC TscR (ps) r0 (ps) TC/TSCR0.1 0.33 0.12 0.360.2 0.87 0.31 0.360.3 1.50 0.55 0.370.4 2.50 0.80 0.32extracted from the output current waveforms given by the program DRIVHBT (AppendixB, List 2). An example of the output is shown in Fig. 3.1.By way of a check on the validity of our approach, we computed TC/TSCR for the uniformfield case, but with ni(x, t) set to zero (and T1 to infinity). This simulation of an effective one-valley situation, in which a uniform velocity accompanies a uniform field, yielded TC/TSCR =0.49, i.e., essentially the expected value of 0.5.Considering now the inverted-field HBT, recall that this structure, in which the n-type28Jo X=OJo XX SCF—-. ////,,//\ /5.00 10.00- 15.00mim, ps2.00 I I I I— 0.5/ — — — — — — — ———0.4/-CSCR- 0.31.00— /—0.2—0.10.00 0.0I I I I00 01 0.2 0.3 04Xc (micron)Fig. 3.2 The transit times and signal delays for the inverted field structure. Thedash lines are transit times and signal delay times vs. the width of collector(xc). The solid line is the ratio of signal delay to transit time.collector in a conventional device is replaced by a p layer (Fig. 2.4b), is intended to improvehigh-frequency performance by prolonging the time spent by electrons in the lower valley [10],[29]. Eventually, the lower valley electrons will scatter to the upper valley, yielding a high-lowvelocity profile which should result in TC/TSCR < 0.5 [38]. Results for the specific structurestudied here are shown in Fig. 3.2. ‘JC/TSCR is, indeed, less than 0.5 and, apparently,decreases as the p layer gets wider. This tendency can be explained by referring to Fig. 3.3and Fig. 3.4 in which the electric field (as computed via the depletion approximation) andthe mean velocity are displayed. The mean velocity was estimated from= fl1V + fl2t’ (3— 7)fll + fl2where n1(x) and n2(x) are computed by the program, and v1(x) and v2(x), the electronvelocities in the lower and upper valleys, respectively, follow from the v — E profile forGaAs (see (2-13)). Equation (3-7) is approximate inasmuch as it does not include the effects2910.00 I I-0.008.00-0.406.00j-0.804.002.00 —1 .200.00— —1.600.4Fig. 3.3 The field profile and the mean velocity profile for the inverted field structurewith XC 0.4 /zm.of diffusion which are, of course, included in the full simulation. However, in view of thegenerally high fields in the depletion region, (3-7) should be adequate for the purpose ofexplaining the trends in Fig. 3.2. The collector width referred to below is XC in Fig. 2.2.For the case of a p collector width of 0.4 tim, Fig. 3.3 reveals that the inverted field structureis successful in that it keeps v high for nearly one-half of the depletion region width. For anarrower collector width of 0.05 .tm (see Fig. 3.4), the associated higher fields result in amore rapid transfer of electrons to the upper valley and, hence, an effectively more uniformvelocity profile. This tendency towards uniformity of v as XC decreases is the reason for theincrease of TC/TSCR towards the uniform-velocity value of 0.5, as depicted in Fig. 3.2.For a conventional collector structure the field is highest at the base-collector metallurgical junction (see Fig. 2.4a). Accordingly, lower-valley electrons injected into the depletionregion transfer rapidly to the upper valley and velocity overshoot is less prominent. Thevelocity profiles for all the collector widths studied in the conventional structure case resembled that for the narrow inverted structure. An example for the case of XC = 0.25 um30E00 0.2)c (mic ron)I I0.008.006.00 —1 .004.00-2.002.00-3.000.00—0.1 0.0 0.1—Fig. 3.4 The field profile and the mean velocity profile for the inverted field structurewith x0 = 0.05 m.is shown in Fig. 3.5. Thus, TC/TSCR values close to 0.4 were obtained (see Fig. 3.6). Itappears, therefore, that previous analyses based on TC/TSCR = 0.5 have overestimated thebase-collector depletion region signal delay time by about 20%Comparing the conventional and inverted device structure, at XC = 0.25 im for example,we note that the signal delay time for the inverted structure is 30% less than that for theconventional structure. The corresponding figure, as computed using the uniform velocityrelationship of TC/TSCR = 0.5 is only 17% . Thus the inverted field structure is even moresuited to high frequency operation than earlier comparisons have predicted [10], [29], [40].I I0.2,c (micron31-0.003.5 ConclusionsFrom this study of signal delay times in electric field regions of GaAs structures, the mainconclusion to be drawn is that the commonly-used practice of taking the collector signaldelay time as one-half of the collector depletion region transit time overestimates the signaldelay time. Specifically:1. for a uniform electric field of 75 kV/cm, TC is about one-third of TSCR;2. for typical inverted-field base-collector structures, TC/TSCR varies from z0.3 to 0.4, depending on the thickness of the collector;3. for typical conventional structures, TC/TSCR is relatively independent of collector thickness and has a value of about 0.4.4. the conclusion from the one-valley analysis of Laux and Lee [38] that, when the velocityof the carriers is larger near the base, then TC/TSCR < 0.5, is still valid for a two-valleysystem in which the charge and velocity are each given by two coupled partial differentialequations.E - 3.5 The field profile and thestructure with XC = 0.25 tim..4-0.40-0.80LU-1.20-1.60O3 0.4mean velocity profile for the conventional field—0_i 0_0 0_i 02,. (micron)I I32I I I -13.00—— 0.5— TC/TSCR2.00— 0.3— 0C/)I—— — 02 0‘,‘Jr,. — I.-.1.00 — —.7-I— 0.1——0.00 0.0I I I I0.0 0i 0.2 0.3 0.4Xc (micron)Fig. 3.6 The transit time and signal delay for a conventional field structure. Thedash lines are transit times and signal delays vs. the width of the collector (XC).The solid line is the ratio of signal delay to transit time.33Chapter 4. The Estimation of Transit Times in GaAsBased HBTsOperating at High Current Densities4.1 IntroductionIn silicon transistors, bipolar devices are often preferable to field-effect devices for applications involving high current densities. This has to do with the inherently different geometriesof the two types of device [41]. The high current-carrying capability of GaAs-or InP-basedHBTs can be expected to be utilized in a number of important applications, e.g. the drivertransistors in integrated optoelectronic laser transmitter circuits [42]. Such circuits mightbe expected to run at very high data rates (e.g. 1-10 GBits/sec [43]), so it is important tobe able to estimate the speed performance of the transistors in the circuit. In Si devicesfT shows a dramatic fall when a certain critical collector current density is reached [44].The mechanism, or mechanisms, which lead to this effect have been given various names:Kirk effect, base push-out, current-induced base widening and quasisaturation [45) [46]. Basically, these phenomena involve either a neutralization of the donor charge in the collectorby the electronic space charge associated with the current, and/or a forward biasing of thecollector-base junction by the induced voltage drop in the relatively high resistance collectorlayer. Whatever the actual mechanism, there will be a reduction in the electric field in theregion of the collector close to the base-collector junction. This reduction in field may leadto a difference in the high-current performance of silicon and Ill-V devices. This is because,in the case of GaAs in particular, the reduced field should reduce the scattering of electronsinto the upper valley and lead to an improvement in the space-charge region transit timeTSCR, and hence to a reduction in the space-charge region signal delay time TC. This effectwould be expected to mitigate the increase in overall emitter-collector delay time associated with the longer base transit time which results from any effective base-widening due tOcurrent-carrier induced space-charge neutralization.34To investigate this possibility by studying the effect of collector current on transit- andsignal delay-times, we need to extend our phenomenological model so that the space-chargedue to mobile carriers is considered. This involves adding the hole continuity equation andPoisson’s equation to the set used in Chapter 2 and 3. The work carried out in treating thispractically-important case of high current operation is described in this chapter.4.2 Extended Model4.2.1 EquationsAs mentioned in the last section, in high current situations, the field created by thecarrier space charge is important. Thus the extended model includes two new equationswith respect to the one described in Chapter 2: one is the transport equation for the holesand the other is Poisson’s equation. The equations for the extended model are:— V. (D2Vn)+ V. (n22V) = ——1)— V. (DVp) — V. (piVq) = 0EV2 = —q(p — — n2 + N — N)where is the potential and the other parameters are the same as described earlier. Theparameters for holes are [47]:= 10.6 cm2/sec. (4—la)= 409 cm2/V sec.where D is calculated from Einstein’s relation with T = 300 K.4.2.2 Estimation of Signal Delay TimeIn Chapter 3 the collector-base space-charge region signal delay time was computed fromthe phase difference between the conduction current entering the depletion region and thetotal current. The latter was calculated by a summation of the conduction current over thedepletion region (see equation (3-2)). Now that we are considering the free carrier space35charge, we can no longer use the depletion approximation. Therefore there is no obviouspoint at which the depletion region can be said to begin or end. It follows that a differentmethod must be used here to estimate the signal delay time. We take as the phase referencepoint that of the sinusoidal input signal applied to the base at the base-emitter boundary(x = —XB in Fig. 2.2). The total current is conveniently approximated by the conductioncurrent in the low-field regions of the device, i.e., regions where displacement current canbe expected to be very small. This was confirmed by computing the conduction currentboth in the neutral base region (—XB < x < 0) and in the neutral sub-collector region(xc < x < xcc). The results were identical. Therefore the phase difference between theinput signal at x = —XB and the conduction current at XC < x < XCC was taken as thesignal delay time. Clearly, with this method it is not possible to separately estimate TB andTC.4.2.3 The Method of SolutionIn the estimation, two different methods are used. In one of the methods, the so-called“box” method is used to discretize the system equations for the Newton-Raphson method.The other method is the finite element method which was used in the later stages of theestimation. These two methods are described separately in the following.A) The “Box” MethodThe Newton-Raphson method was chosen for the numerical calculations. This methodhas been described in detail elsewhere (see for example [48]). It is stable and suitable forthis application. A new subroutine package was written for solving the partial differentialequations. In this package, the subroutine SLE in MTS is used to solve the matrix.equation[49]a11 a12 ... X1 Yia21 a22 a2 = (4—2)a1 afl2 ... a,,, x,36Where {a3} , i,j 1,2,• n are the elements of the Jacobian matrix [47]. {x2} i 1,2,. . .are the discretized solutions to be found, and {y2} i = 1,2, n are the vectors calculatedfrom the discretized system. For more details about the physical meaning of {x2} and {y},please see [47].The following description of the “box” method is based on that given in the reference[47]. The equations in (4-1) can be cast into the following form:oF OG--+-—+H=O (43)where the functions F and G areG G(F, OF/Ox, x, t)(4 —4)H = H(F,OF/Ox,x,t)F can be n1, n2, or p; the function G usually represents the current density, e.g., —D1Vn +n1Vç; and the function H usually represents the net transition rate, e.g., Poisson’sequation can be considered as a special case of (4-3) without the first term OF/Ot. In thisspecial case, G represents the flux of the electric field and F, the charge density.In the “box” method, for the ith grid point, one has+ t) -1(t))+b(F(L + t) — F(t)) + c(F1(t+ t)—+[G+112(t+ t/2) — G1_12(t + t/2)]++h[aH+i(t + zt/2)+bH(t + Lt/2) +c2H1_(t+ L.i/2)] = 0(4—5)whereG÷112 = G+112(F,.—) , xt+12t +“ xIi+1/2 (4—6)IOF\H, = H(F,37and= (x1+i —a1 = (1 +C = (1 + rj(1+rj)12 (4-7)= 1 — a1 —xi — x1_ixi+1 —To linearize (4-6) for the Newton-Raphson method, one needs to compute the following:0G1÷12 — 0G1÷72 0F1÷12+0G1÷12 O(OF/Ox)÷112OF, — 0F112 OF3 O(OF/Ox)11j2 4 8OH, — OH, OF OH1 O(OF/Ox)1 —OF — 8F + O(OF/O) OFNow we choose a set of functions to approximate the partial differentials such that they areexplicit functions of F1. For example:() = dF1÷ + e1F +F1÷12 = (F1 + F1)/2 (4 — 9)(OF’\ 1=—(F11—F)8xJ1+i2 zx,where= xi+1 —d1= (x1_ —x1)(x_—2x1 — x — x1 (4 — 10)= (x1 —x1_)(x —fx1—xi_—(x+1 —x1_)(x+— x1)One can verify that if zx1 = Lx11 then (4-10) is reduced to d, = —0.5, ft = 0.5 and e, = 0.This will give an approximation for partial differentials with evenly-spaced grids.The combination of (4-5)—(4-7) and (4-9)—(4-l0) yields a group of equations which areonly functions of {F}. If G and H are the linear functions of F and OF/Ox, then this groupof equations will be in the form of (4-2) and can be solved for {F1} directly. But in our38case, G and H are not linear functions of F (more specifically, V), so the Newton-Raphsonmethod is applied and (4-8) is used to linearize the equation group [47]. The NewtonRaphson method usually contains an internal ioop in which the correction to the solution iscalculated and when the correction is small enough to meet the criterion, the finding of thesolution to the non-linear equation group is assumed complete. Reference [47] describes themethod in more detail.B) The Finite Element MethodThe description of the finite element method here is according to [50]. The finite elementmethod treats the system as an assembly of elements. An element is defined by a givennumber of nodes i = 1, 2,. . . , m where m is the number of nodes in the element. Thesolution of the systemd2 d 0 (4—11)for the element is approximated by a set of the so-called shape functions:y(e) = N(e)(X)y) (4 — 12)where the superscript (e) denotes the element, functions N, i = 1, 2,.. . , m are the shapefunctions. Once the shape functions are chosen, then the solution for the element is determined by the values of nodes y, i = ,m. The major objective of the finite elementmethod is to find e) = 1,2,• , m for each element.There are several methods (or points of view) to find the solutions. One of the generalschemes is the so-called Bubnov-Galerkin method [50]. In this method, the values of thenodes are determined by the following conditions:JD(e)g(,,y,x,t)N,dD = 0, i = 1,2,• (4— 13)The above equation group is usually assembled into a matrix equation and easily solved fory, i = 1,2,• . , m for each element by means of computer. For more details of the method,one can refer to [50]. The central part of the program used in the calculation is adopted39also from [50]. For the system described in this chapter, the following simple scheme is used:assume the solution at time t is known, then for the solution at time t + zt, the transportequations for electrons and holes are solved first then the solutions are used to solve Poisson’sequation. This is good enough since the time step zt used is very small (0.OOlps). The errorintroduced by such a scheme is the retard error. But from our experiences using the “box”method, this should not be serious.4.2.4 Boundary ConditionsThe current is controlled by the boundary conditions at the base-emitter interface of thedevice. In general, a steady-state solution is first found. The boundary conditions at thebase-emitter interface (x —XB in Fig. 2.2) for steady-state are:ni(x = —XB) = constant(4 — 14)n2(x = XB) = 0To estimate the signal delay, the steady-state solution is used as an initial solution, and asinusoidal component is added to the boundary conditions at the base-emitter interface. Thefollowing are the new boundary conditions for the transient analysis= XB) nio(1 + bsin(27rt/To))(4—15)fl2(XXB)0where b << 1 is a constant, and T0 is the period of the sinusoidal signal. As in Chapter 3,a value of 8ps was chosen for T0. The boundary conditions at the collector contact are notso clearly specified. We used the following at x = xcc (see Fig. 2.2)n2 = a : b(4—16)l + 2 =where a and b are given constants, and N is the doping density in the sub-collector. Thedifficulty arises because the ratio of n1 /n2 is not known. Fortunately, as is shown in Section4.3, the results for the signal delay estimation are not sensitive to the precise ratio of a : b.4.2.5 Programs Used in the Calculations40Again there are two sets of programs used in the calculations, they are described below.A) The Programs for the “Box” MethodThe subroutines developed to implement the algorithm described in the above sub-sectionare listed in Appendix C, List 1. The main subroutine is FINE. Using the solution at t, FINEgives the solution at t + it. The solution is returned by array Yl. Besides FINE are somesupporting subroutines; e.g., FINIT is a subroutine to initiate constants including a-, b,c1, d-, e, and ft; FUNB computes the right-hand side of (4-2); and FUNJ computes theJacobian.Appendix C, List 2 gives the code for EQUI2, a program to find the steady-state solutions. It specifies values for physical parameters for the system, such as D1 and D2, and1L2. Several other subroutines are also listed which provide the necessary support for FINEto solve the equations. GRIDS specifies the grid size used in the discretization, which isnot necessarily even-spaced; FLUX computes the function G in (4-3), which is the currentdensity in this case; FUNH computes the function H in (4-3), which is the right-hand side of(4-1) (with minus sign); DGDU and DHDU give ãG/ÔF and ÔH/ãF, respectively; finally,BNDRY specifies the boundary conditions.DELAY is the program which simulates the transient responses of the system (see Appendix C, List 3). It is much the same as EQUI2 except for the boundary conditions andthe format of the outputs.B) The Programs for the Finite Element MethodThe main program for the finite element method, HBTFL, is listed in Appendix C,List 4. Here are the functions of some important subroutines: NODES defines t’he nodesby specifying the x value; ELEM defines elements of the system; and FINITE is the mainsubroutine which assembles the matrices and solves the matrix equation for the values of thenodes.41In Appendix C, List 5 is the code for FINITE. FINITE includes subroutine AMATRIXwhich assembles the matrices (and vectors), BNDRY which is the subroutine to insert theboundary conditions, and SYSTEM which solves the matrix equation.Some supporting subroutines are in Appendix C, List 6, including CALEL which calculates the element matrix and SHAPE which gives the shape functions.Finally, Appendix C, List 7 is the code for a subroutine which assembles the matrix forthe system, and Appendix C, List 8 includes a package which uses Gauss elimination to solvethe matrix equation.4.3 Results and Discussion4.3.1 High Collector Doping CaseThe basic device structure is shown in Fig. 2.2. the doping densities in the base, then-type collector and the sub-collector are listed in Table 4.1. The maximum base dopingdensity was limited to this value of 1018 cm3 by convergence problems encountered whenusing higher values with the “box” algorithm described in the previous section. The relativelyhigh value of 1017 cm3 for the collector was chosen as it gave a fully depleted collector underlow injection conditions at VBC=OV, i.e., at the bias considered in this part of the work. Itwas thought that with such an arrangement, the switch from a conventional field pattern toan inverted field pattern at high current would be very clear.Two sets of results are presented here. The set A uses the following boundary conditionfor the collector contact:ni(xCC) : n2(xCC) = 1: 10 (4 — 17)whereas set B uses:ni(xCC) : n2(xCC) = 10 : 1 (4— 18)The set A represents a case in which most of the electrons are in the upper valley at thecollector end. On the other hand, the set B represents a case in which most of the electronsare in the lower valley at the collector contact. The purpose of collecting results for the42Table 4.1Doping densities in cm3 for the device shown in Fig. 2.1Doping density Width (nm)Base 1 x 1018 100Collector 1 x iO’ 250Sub-collector 4 x iO’7 100Table 4.2Signal delay (TB + TC) with different boundary conditions.Estimated from the “box” method.Current densities are x105 A/cm2.Current density n1 : = 1 : 10 n2 = 10 12.42 0.55 0.505.63 0.45 0.4510.3 0.85 0.8514.4 0.95 0.95above two sets of data is to investigate whether the ratio of ni(xCC) n2(xCc) is significant.The results for the total (base plus collector) signal delay time (TB + TC), as computed fromthe phase delay between the conduction and total currents (as described in Chapter 3), aresummarized in Table 4.2. The conclusion from these results is that the boundary conditionsdo not affect the results much. The following simple explanation of this finding is offered.In a region where the concentration of carriers is high and the electric field is low, thetype of carriers (lower-valley or upper-valley electrons) becomes unimportant be€ause thecurrent is a drift current and the drift velocity is small. Different types of carriers onlycause a small difference in the field in the region for the same current. Thus, a differencein the ratio of the types of electrons in this region causes only an insignificant difference43in the signal delay. This is fortunate because it is not known how to arrive at a precise,physical-based value for the boundary conditions.The field profiles and carrier concentrations corresponding to the current densities listedin Table 4.2 are shown in Fig. 4.la-h. As the current increases, the field at the collector-basejunction weakens and, somewhere between J = 5 and 10 x i0 Acm2 (e.g., see Fig. 4.lf.1and Fig. 4.lg.1), the field in the collector regions reaches a near-constant level. Furtherincrease in current results in the inverted-field profile shown in Fig. 4.ld.1 and Fig. 4.lh.1.This evolution of the field profile occurs in silicon devices and is well-documented [44]. Thedoping density of the collector in our case is too high for significant base push-out to beseen. This is illustrated in Fig. 4.2 where it is clear that the hole profile changes very little.What is different between the silicon and gallium arsenide cases is the effect that thechanging field profile has on the electrons. The progressively inverted nature of the fieldwith increase in current causes less and less transfer of electrons to the upper valley (see,e.g., Fig. 4.le.2 to Fig. 4.lh.2). This should lead to a reduction in electron transit time withcurrent, and would be expected to result in an improvement in signal delay time. However,as Table 4.2 indicates, this is not the case. Instead, there is only a slight decrease and then,much as in the silicon case, a rapid increase in delay time.To try and see the cause of these interesting results, the signal delay times were alsocomputed using the model described in Chapter 3. The field profiles from the extendedmodel were used in the old model, so the major difference between the two is the absenceof holes in the old model. From Table 4.3 and Fig. 4.3 one sees that the signal delay timedecreases with current density when oniy the electrons are considered. This is the expectedresult and leads us to the important realisation that the holes must play a significant role indetermining the signal delay time at high current densities.To investigate this phenomenon further, simulations were carried out with a lower collector doping density. This should lead to more pronounced base push-out effects. The resultsare described in the next section.44I I I I I•140 L 5 2/ . J=2.42x10 A/cm/ :.o.X H / -/—>40.00[-1) 0.60H I—100.00- 0.40— —-120.00 U 0.20 L —-140.00 -H 0.00 -I I II I-100.00 0.00 00.00 000.00 300.00 —100.00 0.00 100.00 200 .00 000.60distance (nm) distance (nm)(a.1) (a.2)I I‘40L J=563x10A1cm 2-80.00 -ov 0.60H4,434.4 -100.00- V0.40 -6-44-4—120,00— V (.) 0.20 II,-140.01:V 0.00— — — — -1 I I I I I—100.00 0.00 100.00 200.00 300.00 —100.00 0.00 100.00 260.00 300.60distance (nm) distance (nm)(b.1) (b.2)Fig. 4.1 Carrier density profiles and fields for different currents. (a)—(d) are forni(xcc) fl2(XCC) = 1 : 10.45E0-4SS-461)-H4-4C-)-44)-HIto061)8-401)-H0-40-4061C)Fig. 4.1 Carrier density profiles and fields for different currents. (a)—(d) are forni(xCC) fl2(XCC) = 1 : 10.r I0.00—80.00 ———12000—-140.00-L I I—100.00 0.00 100.00 200.00 300.00distance (nm)(c.1)C)-40-4>14)-H011)4-401)-H4-44-4(II01.401.201.000.800.600.400.200.00- 5 2J=103x10 A/cm-[-. H—100.00 0.00 130.30 200.00 300.00distance (nm)(c.2)//—80.00 —.-4(1)—100.03———120.10 — —-i40.12I—100.01 0.00 100.00 201.00 300.00distance (nm)(d.1)I I I I5 2 -1.40 J=1 4.4 x 10 A/cm“p1.20 —3.80 H0.60- -:::0.00 --——II I I—100.00 3.10 100.00 200.00 300.00distance (nm)(d.2)46I I I0.00/ 140[- 5 2/ — J=2.42x10 A/cm40. /I :t-100.00-- OAO fl1 _-120.0Cr- 02o[— / -1-140.0Gb I-GG.00 0.00 100.00 200.00—1.X 0O.OOdistance (nm) distance (rim)(e.1) (e.2)IL J=5.63x1Ocm 2--40.OOL \ --60.GGF o—80.00-oo—a)-100.C0r- a) fl --120.00-14C.20 H oao ----- —II I I I I J II I I-100.00 0.00 100.00 200.00 300.00 —10OO 0 O.X O.OOdistance (rim) distance (rim)(f.1) (f.2)Fig. 4.1 Carrier density profiles and fields for different currents. (e)—(h) are forn2(XCC) = 10 : 1.47I I I1A4 5 2-J=1O.3x10 AJcm120L: “p-:E1000- — — — --10000 010 1 000 30000distance (non)(g.2)-1Fig. 4.1 Carrier density profiles and fields for different currents. (e)—(h) are forni(XCC) n2(xCC) = 10 : 1.C-)C,,‘-4‘0‘-4-H4-40.00-20.00—40.00—60.00—60.00—100.02-120.00—243.00I IC-)on‘-4on>1.4J0)01)‘04-401)-Hs-I011Ci-100.00 0.00 :oo.oo 200.20 300.00distance (nm)(g.1)0.00—20.00-40.00C)-60.00on-60.00‘0‘—401)- -100.00’—12 0.00—140.00C)on-4>14)-HOn01)‘00-101)-‘-44-44-4(11C)54012011500400400400200150J=14.4x10 A/cm 2‘p-10000 0150 10000 20000 30000distance (nm)(h.2)—100.00 1.00 100.00 200.00 300.10distance (non)(h.1)48d1.40 —1.20—- 0.80base push-outF—o.___________________-100.00 0.00 100.00 200.00 300.00distance (nrn)Fig. 4.2 Base push-out. Curves (a)—(d) are the profiles of holes from Fig. 4.1(a)—(d),respectively.4.3.2 Normal Collector Doping CaseAs well as making the change in collector doping density as mentioned above, the basedoping density was also changed (to the more normal value of 1019 cm3). The use ofthe finite element method, as described in section 4.2, overcame the convergence problemsassociated with the “box” method (see also section 4.2) at high base doping densities, soallowing this case to be examined. The device parameters are listed in Table 4.4. As in theprevious section, VBC was held at OV. As the ratio ni(XCC) : n2(XQC) did not seem to affectthe results in Section 4.3.1, just the case of ni(xGc) : n2(XCC) = 10 1 was investigated here.Results for the field and doping density profiles are shown in Fig. 4.4a-e for five differentcurrent densities.The field evolves in much the same way as for the example studied in the previous section.The situation at J, = 1.68 x io Acm2 (see Fig. 4.4c.1) is particularly interesting as the field491.10 -1 .00 r—á5 0.405—0.301—....-model0.20—.0.10—0.00—I0.00 5.00 10.00 15.00Current densityJ (x 105A/cm2)Fig. 4.3 Comparison of the old model and the extended model.Table 4.3Signal delay (TB + TC) computed using the extended model with boundary conditionfli(XCC) : l2(’cC) = 10: 1, and using the old model (Chapter 3).Current densities are x105 A/cm2.Current density Extended model (ps) Old model (ps)0.0263 0.60 0.580.196 0.60 0.612.42 0.50 0.525.63 0.45 0.2910.3 0.85 0.2214.4 0.95 0.19is almost constant, indicating the collapse of the space-charge region at the base-collectorjunction. Thereafter, a nearly-neutral region appears (and grows) on the collector side of the50Table 4.4Doping densities in cm3 for the device shown in Fig. 2.1, for normal collector doping case.Doping density Width (nm)Base 1 x 100Collector 2 x 1016 400Sub-collector 5 x 1018 100junction. This signifies a high injection condition in which holes from the base are injectedacross the junction which has now become forward biassed by the voltage drop across thecollector n region. The hole portion of this electron-hole plasma can be seen developing inFig. 4.4d.2. It is very pronounced at the highest current studied (7.2 x iO Acm2 (seeFig. 4.4e.2) Also clearly visible in this figure is, on the collector side of the thickened quasi-neutral region, a thin electron accumulation layer. According to [45] this layer is due to thereduction in hole density caused by the effectively increasing potential as one moves towardsthe sub-collector.It is quite clear in Fig. 4.4e.2, and at slightly lower current (Fig. 4.4d.2) that in advanceof the electron-hole plasma is a region which has gone above the critical condition at whichthe equilibrium electron density determined by the dopant concentration in the collector canno longer support the current. Thus electrons must accumulate in this layer as more andmore current is forced through it. This space charge layer attracts holes from the base. Sothe quasi-neutral hole-electron plasma grows, effectively thickening the base region. Thephenomenon is clearly seen in Fig. 4.5.The above description of high current effects follows closely that given by Liou, Lindholmand Wu [45] in their qualitative preamble to a paper on the development of J-V models forHBTs operating at high current densities. The novel feature of our work is the estimationof signal delay times for these conditions. The results are shown in Table 4.5. As was thecase for the device studied in Section 4.3.1, the diminishing field at the collector side of51Fig. 4.4 Carrier density profiles and fields for normal doping profile.52° L- LC-)> -•1aih— -— -I-jü) -14O.-- --H-l6O. - --l8O.- -HOf3distance (nm)(a.1)J=0.0478 x 105A/cm23p1e18 / H> I’clel7L I’ci) Ji•03 /ci) /1ielci[distance (nm)(a.2)cnr--ie- -l1- -Hci) -140C0 - --r-1-l6O - -distance (nm)(b.1)lel 5 2J0.480 x 10 A/cm3-p 1/1e18 - / — —:E_\I1e1531e14distance (nm)(b.2)Fig. 4.4 Carrier density profiles and fields for normal doping profile.53>.sorH-i‘-Ici) .I----44-116o4 I4distance (nm)(c.1)1e19’5 5 2J=1.68 x 10 A/cm3-pI11::n21e17 -NS-1ei6rA’Ie153 j1e4‘Idistance (nm)(c .2)C)C)>1-H9:5ci)1e4-4C-)9:,—Ici, -14O.3—distance (nm)(d. 1):ei9-----‘j 5 2I J=4.83 x 10 A/cmI_fl’/e18-1e17\\ i n2j1e1 / j1 I Idistance (nm)(d.2)I IFig. 4.4 Carrier density profiles and fields for normal doping profile.the base-collector junction first results in an improvement in signal delay time, due to thereduced transfer of electrons from the lower to the upper valley. However, once the holeinjection into the collector becomes significant, the signal delay time starts to rise. Thus,once again, we see the importance of the holes.To emphasize the importance of the holes we present results in Table 4.6 for differenthole parameters. These results are also shown graphically in Fig. 4.6. Increasing the holemobility from 400 cm2V’s to 2,000 cm2V’s’ does lower rB + TC in the hole-transportdominated regime. However, increasing the hole saturation velocity from 1 x i07 cms1 to5 x 1W cms1 has virtually no effect. This confirms that it is the transport of holes acrossthe low field portion of the device (from the base into the quasi-neutral collector) that islimiting the signal delay time. To our knowledge it is the first time that this important roleof holes has been so identified.QLU>-G) •I4O--d14•16 --J=7.20 x 105A/cm2 -n -1e193C-) 1e18>1•1—) 31e171e16C-)31e153/ —,--distance(e.1)(nm) distance (nm)(e .2)54::::>le-4-17 —...-3 I1c161-.-15 H—3 F—\\1e14 a\C—100.00 0.00 100.00 200.00distance x (nrn)Fig. 4.5 Profiles of holes. Curves a-e corresponding to Fig. 4.4a-e, respectively.Table 4.5Signal delay (TB + r0) calculated using the finite element method.Current densities are x iO A/cm2.Current density Signal delay (ps)0.0478 0.550.480 0.421.68 0.344.83 1.36‘ 7.20 3.754.4 ConclusionsFrom this study of the signal delay times in HBTs operating at high current densities, it canbe concluded that:1) the transport properties of both electrons and holes are important in determining the554OO I I3.50[—300—2.5O a —-H-0.000.00 2.00 4.00 6.005 2Current Density Jc ( x 10 A/cmFig. 4.6 TB + TC vs. J with different properties of holes. Curve a is the normalcase with ,u = 400cm2V’s’ and u.. = 1 x 10 cms’; curve b is with =2000 cm2Vs’ and v8 = 1 x 10 cms1; curve c is with jz = 2000cm2V’sand u5 = 5 x 1OT cms1; and curved is the one-valley case with all electrons inthe lower valley.signal delay time;2) once significant base push-out has occured it is the holes, because of their larger effectivemass, which make the greater contribution to the signal delay time;3) the inversion of the field in the collector region as the current density increases doeslead to an improvement in the signal delay of the electronic contribution to the current.However, this effect, which might have been expected to lead to a different behaviour thanis seen in single-valley semiconductors, is hardly significant because of the dominance ofthe heavy holes;4) to improve the high-speed performance of n-p-n HBTs it appears that semiconductorsin which both types of carriers have low effective mass are desirable.56Table 4.6Signal delay (TB + ‘7-c) with different hole properties.Current densities are x i0 A/cm2,mobilities are in cm2Vs1andvelocities are in cms1.Current density Delay (ps) Delay (ps) Delay (ps)= 400 = 2000 = 2000v3=1x107 v=1x107 v=5x1070.048 0.55 0.55 0.550.480 0.42 0.42 0.421.68 0.34 0.32 0.324.83 1.36 1.04 1.037.20 3.75 1.89 1.8257Chapter 5. Analysis of the Formation of Stationary Domainsin GaAs MESFETs5.1 IntroductionThe nature of the electron drift velocity versus field relationship in GaAs dictates that, undercertain operating conditions, doping densities, and channel dimensions, a stationary dipoledomain will form at the drain side of the gate in MESFETs [12][13][14][15J[16][17]. Thedomain is stable in those circumstances where the conducting channel beyond the gate onthe drain side widens, causing the field to drop below the value necessary to sustain travelingGunn domains. The presence of a stationary domain in a MESFET influences the saturationcurrent and the capacitive coupling between the device’s terminals and the current gain cutoff frequency fT [51]. While the formation of a stationary domain is beneficial in that itreduces the drain-gate feedback capacitance, it is undesirable as regards fT as it causes thisparameter to decrease strongly with drain bias. These important practical phenomena makeit of interest to develop a criterion for the formation of stationary domains. There are somewidely accepted criteria for the formation of traveling domains, namely: [26][52]NDL> 1 x 1012 cm(5—1)NDa > 2 x 1011 cm2where ND is the doping density (assumed uniform), a is the active layer thickness, and Lis the gate length. But these would seem to be unsatisfactory as criteria for stationarydomain formation because two-dimensional modeling [13] reveals that stable operation ofMESFETs, in which the domain is stationary, is possible even in devices with parametervalues exceeding these criteria. Therefore, in the work described in this chapter, we seekto develop a criterion for stationary domain formation. Specifically, we propose a onedimensional model in which the analytical approach is based on the phase portrait methodof representing the solutions to ordinary differential equations. The analysis yields two58necessary conditions for the formation of stationary domains. One of these is material basedand relates to the velocity-field relationship. The other is device based and provides a newcriterion for the minimum gate length at which stationary domain formation is possible. Thelatter may prove useful in the design and analysis of MESFET’s.5.2 One-Dimensional Analysis5.2.1 The Fundamental EquationsAs Fjeldly has pointed out [15], the cross-section of the conducting channel varies littleover the source side and the accumulation region of the domain. Therefore, for an approximate analysis, it is justifiable to use one-dimensional forms of Poisson’s equation and thetransport equation, i.e.,(5—2)J— —qnv+qD (5—3)where n is the local electron density, q is the absolute value of the electronic charge andp and D are the field-dependent carrier velocity and diffusivity, respectively. Note that, inaccordance with the directions specified in Fig: 5.1, v> 0 and E < 0. The assumption of aconstant channel cross-section, in conjunction with the need for current continuity, impliesthat J is constant along the channel. If electrons accumulate in the channel such thatthe field at some point becomes larger than E (see Fig. 5.2), v may decrease to such anextent that the flow of electrons by diffusion must augment the drift flow in order to keepJ constant. This situation will lead to a decrease in n. If the gate is long enough, n maybecome less than ND, and a depletion region will appear. The juxtaposition of accumulationand depletion regions constitutes the stationary domain. It is the goal of the work in thischapter to compute the minimum gate length Lmjn at which a stationary domain can form.59fl-NDND(b)DEPLETION REGIONGATE / DOMAINDRA__1___ ___IL-—zFig. 5.1 Device cross-section, and profiles for the electron concentration and theelectric field for the case of a MESFET with a stationary dipole domain at thedrain end of the gate. Note the assumption of a constant depth of channel underthe gate. [81]Introducing the dimensionless variablesdY1 — qr---_--(1-Y2)dY2 — v(Ei) fv(Ei}”i)dxD(E1)l v(Ei)SOURCE I : : : I‘I(a) fl—GaAs— :_____E ::(C)0allows (5-2) and (5-3) to be rewritten as(5—4)= n/Ne (5— 5)where E1 is the lesser of the two fields which satisfy the equation= —qNDV(E)(5—6)(5—7)(5—8)60VEp E2Fig. 5.2 General form of the equilibrium electron velocity versus field relationshipfor GaAs. The electric field is taken as negative, in accordance with Fig. 5.la.where J is the appropriate constant value of J. E1 can be thought of as the field inthe channel which would produce a current density J if the channel were entirely neutral.Fig. 5.1(c), which is based on the results of two-dimensional modeling [13][16j, indicates thatthe actual field under the gate is larger than E1. This is because of the diffusion currentthat opposes the drift current in this region. The field at the source end of the gate, i.e.,at x = 0, is labeled IEA’I to be consistent with the notation used in the phase portrait ofFig. 5.4. Note also from Fig. 5.1(c) that 1E > EC where EC is the field in the neutralregion between the source contact and the source end of the gate. This is simply due to thedifferences in cross-sectional area due to the presence of a depletion region under the gate.5.2.2 Phase PortraitsThe differential equations (5-6) and (5-7) are autonomous. They do not explicitly containthe parameter x. Their solutions, therefore, can be represented in a phase plane (Y12) asa family of curves, directed with increasing x, called trajectories. The resulting figure iscalled a phase portrait (see, for example, [53}). A similar representation, but using nonnormalized variables, has been employed in studies of the effect of boundary conditions onthe electrical behavior of homogeneous GaAs diodes [54][55j. For a qualitative examination of61V2I,0 1 V1Fig. 5.3 The zero-value isoclines for equations (.5-6) and (5-7). The short line segments on the isolines have slope equal to 0 for II’ and cc for 00’; they indicatehow the solution curves cross the isoclines.the MESFET situation, we look first for the fixed points of the system, i.e., the solutions forthe case when dY1/dx = dY2/dx = 0. Inspection of (5-6) and (5-7) reveals that (Y1,Y2)=(1,1)is a fixed point. From Fig. 5.2 we havev(E1) = v(E2). (5—9)Using (5-9) in (5-6) and (5-7) uncovers the other fixed point, (Y1,Y2) = (E2/1,1). Thecurves corresponding to either dYi/dx = 0 or dY2/dx = 0 are called isoclines in the phaseplane (Y1}). The fixed points, together with the isoclines 00’ (for dYi/dx = 0) and II’(for dY2/dx = 0) are shown in Fig. 5.3. 00’ is simply Yj = 1 from (5-6). II’ is a solutionof v(E1Y)Y2= v(E1) (from (5-7) at dY2/dx = 0) for the v versus E relationship of Fig. 5.2.The short line segments cross the isoclines at dY2/d1 = (d}/dx)/(dYi/dx) 0 for the caseof II’, and at d)’/dr = (dY/dx)/(dY/dx) = 0 for the case of 00’. The arrows indicatethe directions of the solution curves as they cross the isoclines. For example, for dYi/dx = 062I dYr20dV1—00 dxI I“20 viFig. 5.4 Phase portrait for the system of equations (5-6) and (5-7).(i.e. }‘ = 1) and Y < 1, (5-7) indicates that dY2/dx < 0, whereas for dYi/dx = 0 andY1 > 1 we have dY2/dx > 0.Knowing the directions of the solution curves and the positions of the zero-value isoclineswe can sketch the phase portrait. Fig. 5.4 gives the qualitative picture. The actual portrait,drawn to scale, would be much compressed in Y2.As Fig. 5.1(b) indicates, when a stationary domain is present in the channel, n NDat the source end of the gate and at the maximum field position in the domain. Thisphysical situation is represented on the phase portrait by a trajectory which crosses the 00’isocline (i.e., n = ND) twice. In accordance with the field profile of Fig. 5.1(c), the firstcrossing must be in the direction of increasing Y1 and the second crossing must be in thedirection of decreasing Y1. Inspection of Fig. 5.4 shows that only the trajectory A’A”A”meets these conditions. The neighboring trajectory B’B” does not cross the isocline 00’twice, indicating that no domain is present for this solution. Clearly, somewhere between A’and B’ there is a critical point (Y1c, 1) which marks the boundary between starting pointsfor the two classes of trajectory. However, as is demonstrated below, it is not necessary to6310Aknow the actual position of the critical point in order to arrive at the criterion for stationarydomain formation we are seeking. This is fortunate as there does not appear to be a simpleanalytical method of determining the critical point.5.3 Criteria for Stationary Domain FormationA necessary condition for the formation of a stationary domain in a real device is that thegate must be long enough for the electron concentration to decrease to ND. Beyond the gaten must then fall below this value, as required in the creation of the depletion zone of thedomain. On the phase portrait this is manifest by a solution that follows trajectory A’A”out to, and beyond, A”. (Recall that movement along a trajectory is in the direction ofincreasing x.)Often in the solution of a set of nonlinear differential equations it is convenient to considerthe first-order approximations to the equations about a fixed point. In the present caseit turns out that consideration of the linearized form of (5-6) and (5-7) about the fixedpoint (E2/1,1) leads to a straightforward solution for the minimum gate length for domainformation. The first order equations aredZ1 qN 5—10dZ2 — v(Ei) dv(E2)z v(E2) Z 5 11dxD(E) dE 1+D(E2 2 —whereZ1 = Y1 — (5 — 12a)andZ2=Y—1. (5—12b)Equations (5-10) and (5-11) comprise a linear homogeneous system of differential equationswith constant coefficients. The eigenvalues ) can be computed from the determinant—AEi dv(E2) v(E2) —D(E2) dE D(E2)64which yieldsv(E2) 11 ± 1 — 4qND(E2)dv(E2)1/25 13—2D(E) ev2(E) dE —Recalling that v> 0 and E <0 (see Fig. 5.1), then for El > lEl (see Fig. 5.2), we havedv(E2)> (5 — 14)which allows for the possibility of ) being a complex number. This is the situation we areseeking as the magnitude of the imaginary part of .A can be defined as the reciprocal of thecharacteristic gate length L0, whereupon a solution for Z2 can be writtenZ2 =A1e°+A2e_fl0 (5 — 15)where c is the real part of the eigenvalue, and A1 and A2 are constants. It follows from theboundary condition n ND at x = 0 (i.e., at the source end of the gate) that A1 = —A2.We are interested in the condition when n = ND at high fields as this would indicate theformation of a domain (see Fig. 5.1(b)).When n = ND, Z2 = 0; therefore, (5-15) gives= 1 (5 — 16)for which two solutions are x = 0 and x = irL0. The latter can be viewed as the minimumlength of gate required for formation of a stationary domain. The true minimum valueoccurs when L0 is a minimum, i.e., 1m2(k) is a maximum. It follows from (5-13) that theseconditions are met at the value of E2 which maximizes 1m2(.X(E)) in the expressionJm2p(E))= eD(E2)dv(E2)—> 0. (5 — 17)For a given ND, and known relationships for v versus E and D versus E, (5-17) can be solvednumerically to find the maximum. Labeling E as the value of field at which this occurs, theminimum gate length is given byLmin lIm((E))V(5-18)65Some results are given in Section 5.5The necessary condition for stationary domain formation described by (5-17) can beexpressed asdv(E2) ev2(E) 5 19dE 4qND(E —Note that the right-hand side of this inequality is positive and that we have defined E asbeing negative. It follows that (5-19) can only be satisfied at values of E( in excess of(E( (see Fig. 5.2). This implies, on the basis of one-dimensional considerations alone, thatstationary domain formation is not possible in MESFET’s made from material such as siliconin which v increases monotonically with E. Another physical meaning is that the v versus Erelation for the semiconductor material must decrease at above a certain rate in order for adomain to be able to form. The criterion (5-19) has previously been discussed in the contextof stationary high-field domain formation at the anode of Gunn diodes [56].5.4 Programs Used in the CalculationTo determine E2, velocity v(E2), D(E2), and Lmjn for given ND, a simple algorithm wasused. This is a general method to find the optimal parameter for a target function, knownas the Fibonacci Search [57].Suppose that the target function is f(x) where x is a parameter to be optimized. Tofind the maximum of f(x), one starts from two points a and b. The target function f(x) hasone and only one maximum in the region (a, b).Now compute the values of f at x1 = a + 0.382(b — a) and x2 = a + O.618(b — a), wherethe number 0.618 is the value of the ratio of successive Fibonacci numbers [57] at the limitk —f oc, i.e.,Fk1 1urn = = 0.618k—oo 2where Fk, k 0 are the Fibonacci numbers defined asF0=O, F1=1, Fk=Fk_1+Fk_266Now if f(x1) > f(x2) then b is replaced by x2, otherwise a is replaced by x1. Repeat theabove process until the value of b — a is less than a given number, say, 1 x i0. Then thevalue (a + b)/2 is a good estimate for the optimal parameter. VThe first program, GAAS1, is used to find the value of E2 which optimizes the targetfunction (5-17). It is listed in Appendix D, List 1. Also listed in Appendix D, List 1 are theinput and output files of the program. The function f in the program GAAS 1 is the targetfunction. The subroutine MINLEN is used to calculate v, D, and Lmin.The second program listed, PEAKS, is used to find the value of E at which the velocityreaches the maximum. The target function in this case is (5-20). The program is listed inAppendix D, List 2 together with sample input and a output files.The third program, SLOPE, is used to calculate E1 from E2. It is listed in Appendix D,List 3. The target function isf(E) = abs(v(E) — v(E2))Instead of finding the maximum of f, SLOPE finds the zeros of f, i.e., the minima. Finally,the subroutine VELO in SLOPE is used to calculate v(E2).675.5 Results and DiscussionTo evaluate Lmin requires knowledge of E, v versus E and D versus E. We take the lattertwo relationships from Fjeldly [15], i.e.,(E)— Em/Es+(Em/Eo)4E 5 20V— —vs 1 + (ErniE0)4 Emwhere Em is the magnitude of the field Ev8 = 1 x lO7cm/s= 3.4kV/cmE3 = E30[1 + (ND/No)”2]= 1.2kV/cmN0 = 5.3 x 10’7cm3andD E — f D0[1 -t- (E/E1)2], E ET 5 21( ) — D[1 + (Eh/E)25], E ET —whereD0 = 200cm/s= 1.66kV/cmEh = 6.5kV/cmET = 3.5kV/cmUsing (5-20) and (5-2 1) in (5-17) allows numerical computation of the value of field E whichrenders 1m2(E)a maximum. This field is tabulated for various channel doping densities inTable 5.1. Also shown are the values of v(E) = v(Ej) from (5-20), D(E) from (5-21), andLmin from (518). Lmin is somewhat sensitive to the choice of D(E), and it is noteworthythat values of D lower than used here have been used by other workers [14], [58]. However,it transpires that Lmin is not strongly dependent on D. For example, by reducing D0 by afactor of 2 to 100 cm2/s causes Lmin to decrease by only 6 percent at ND = 10’5cm3andby 29 percent at ND = 5 x 10’7cm3.68Table 5.1Summary of Computed Results for Various Doping DensitiesND v(Ej) = v(E) D2 Lmin(cm) (kV/cm) (cm2/s) (pm)1 x iO 1.58 vs 4.74 641 1.725 x 10’s 1.49 v5 5.01 584 0.681 x 1016 1.46 v 5.07 573 0.495 x io’ 1.36 vs 5.21 548 0.241 x 1017 1.30 v5 5.30 533 0.185 x iO’ 1.16 v5 5.77 469 0.11To assess the validity of the results from the linearized analysis, we first note that it isimplicit in the first-order analysis that Z1,Z2 << 1. For Z2 this inequality is likely to beeasily satisfied as carrier concentrations of interest in GaAs MESFETs have been estimatedto be such that (n — ND)/ND is much less than unity [15][59]. For Z1 the situation is not soclear cut, but some evaluation can be made by noting that physically meaningful solutionsfor domain formation are represented on the phase portrait (see Fig. 5.4) only by thosetrajectories that start between the critical point (Y1c, 1) and (Y1p, 1), where the latter is thepoint corresponding to E = E at the sources end of the gate. For a trajectory startingat (Y1p,1), Z1 = (E — E2)/E1. Taking the case of ND = 10’7cm3 as an example, withE evaluated for the minimum gate length case (see Table 5.1), gives Z1 0.8. Becausethis value is not much less than unity, some additional error analysis is needed to allowassessment of the usefulness of the linearized solution. To obtain this information, thenon-linear equations ((5-6) and (5-8)) were solved numerically using a fifth-order RungeKutta-Fehlberg approximation [60]. The field E1 was assigned the same values as used inthe linearized method, and the initial conditions were E = E,, and n = ND at x = 0. Theshooting method used to determine Lmtn employed variable step size error control. Theresults are shown in Table 5.2, from which it can be seen that the differences in Lmin for69the two methods are less than 20 percent. This is quite acceptable for the purpose of aone-dimensional analysis. The Runge-Kutta-Fehlberg method failed to yield a solution forthe ND = 10’5cm3case, presumably because the field at the critical point is, in this case,greater than EI.Table 5.2Comparison of Results for Lmin from the Linear and Nonlinear AnalysisND IEI Eu Error(cm) (kV/cm) (kV/cm) (zm) (itm) %1 x iO’ 3.13 2.09 - 1.72 -5 x iO’ 3.16 2.05 0.79 0.68 141 x 1016 3.19 2.07 0.56 0.49 135 x i0 3.31 2.23 0.28 0.24 141 x 1017 3.41 2.36 0.22 0.18 18s x 1017 3.90 3.90 0.13 0.11 15* from the non-linear analysis (Runge-Kutta-Fehlberg)f from the linearized analysisThe product NDLmin is displayed in Fig. 5.5. It is not a constant, and, at long and veryshort gate lengths, it deviates markedly from the value of 10’2cm,which is often used asa criterion for traveling domain formation [131.Qualitatively, the fact that Lmin decreases with ND can be explained by examining thephase portrait and recalling the comments following (5-8) concerning the fact that a diffusioncurrent opposes the drift current in the channel under the gate. This diffusive flow is drivenby a concentration gradient dn/dx, so n must increase on moving away from the source endof the gate (point A’ in Fig. 5.4). This causes E to increase in the channel. As long asEl < lEl, the drift velocity also increases on moving further from the source end of thegate. The associated increase in drift current demands an increase in the opposing diffusioncurrent in order to hold the total current constant. Thus an accumulation region forms.70ECEzFig. 5.5 Comparison of the doping density-gate length product for stationary domain formation as computed from this work with the commonly-used criterionof ADL = 1012 cm2. Also shown are experimental (solid symbols) and theoretical (open symbols) data from devices for which evidence of stationary domainformation is claimed.The increases in E and n are clearly visible on the phase portrait. The turning point for n(i.e., for Y2), is reached when the trajectory A’A” crosses the isocline II’. Physically, thedecreases in electron concentration after this point arises because v now decreases with E,causing the drift current to fall to such an extent that it must now be augmented by thediffusion current to maintain current continuity. If the gate is long enough, n can decline toequal ND. This occurrence appears on the phase portrait as point A” and corresponds tothe drain side edge of the gate. Moving beyond the gate (toward A”) signifies the formationof the depletion zone of the domain and a decrease in electric field. If the doping density isincreased, so too is the drift current. Thus, the constraint of current continuity will demandthat the diffusion current near to the source end of the gate also increase. This means thatthe sequence just described will evolve over a shorter distance.To compare the predictions of this analysis with experimental results requires, ideally,data from a set of devices in which the NDL product varies through the NDLmin value.71N0 Lmin10’0Aoo[62JIlsi. (1 2jI 1171A f611.0.5 1.0 1.5 2.0L (pm)Confirmation of the validity of the proposed criterion would be indicated by the loss ofevidence for stationary domain formation as NDL was reduced below NDLmin. Unfortunatelysuch a data set does not appear to exist. The experimental results shown by the solid symbolsin Fig. 5.5 are from devices claimed to exhibit stationary domain formation [12][16][17][61].They are in agreement with the present analyses inasmuch as they all lie above the boundarydefined by NDLmin.Comparison of the proposed criterion with the results of 2-D analyses [13j[14][62] can bemade by examining the open symbol data points in Fig. 5.5. The result of Yamaguchi etal. [13] is for the application of their domain formation criteria to a 1-nm device. The othertheoretical data do not represent boundary values for NDL, but instead refer to devices thathave been shown by 2-D simulation to exhibit the effects of stationary domain formation.Again, the comparison is satisfactory as the data points lie above the NDLmin boundary.5.8 ConclusionThe main conclusion from this model is that a one-dimensional representation of conductionin a GaAs MESFET can lead to a new criterion for stable, stationary domain formation. Thecriterion sets a minimum length to the gate. The associated doping density-minimum gatelength product is not a constant and differs significantly from the commonly used criterionof NDLmin > 1012 cm. Moreover, the phase portrait is shown to be effective in qualitativeanalysis due to its graphical properties. The linearization of the equations (5-2) and (5-3)and the utilization of the fixed points are two other effective tools in qualitative analysis,although the linearization technique should be used with caution.72Chapter 6. Use of the WKB Method to Calculate the Surface Charge Densityin MODFETs6.1 IntroductionGaAs-based MODFET devices are suitable for high-speed and high-frequency applicationsbecause a very high charge carrier mobility is possible in the two-dimensional electron gas(2DEG) at the interface between doped AIGaAs and undoped GaAs (see Fig. 6.1). Dueto the spatial separation of the ionized donors and the electrons, Coulomb scattering isstrongly reduced, and therefore high electron mobilities are achieved [63]. Due to the lightmass of the electrons in this material, the energy band of electrons in the potential wellnear the A1GaAs/ GaAs interface splits into subbands [64]. Generally, in MODFET devicemodeling, it is important to know the position of these subbands and the correspondingelectron concentrations (surface charge densities). From this information a charge controlmodel can be developed to determine the I—V and C—V characteristics of the device [65].In its most complete form, the problem of estimating the subband characteristics (poition and electron concentration) involves the solution of the equations of Poisson andSchrodinger in a self-consistent manner. Such solutions have been obtained (see, for example, [66] [67] [68] [69]). They involve the iterative solution of non-linear differential equationsand are computationally intensive.Perhaps the most popular approach to the problem is to approximate the potentialprofile at the junction as a triangular well [63][70j[71][72][73], (see Fig. 6J). With such apotential profile, an analytical solution of Schrödinger’s equation for the subbani energylevels is possible [74]. The value of the constant field F on the GaAs side of the well ischosen so that the total surface charge density qn8 computed fromcF=qn3 (6—1)73Tringu1r Po1enti1 W11Fig. 6.1 Band diagram of the A1GaAs/GaAs interface in a MODFET. The dash linerepresents the triangular potential well approximation.is consistent with the total surface charge density computed from the subband energy levelsyielded by the Schrödinger equation [63]. in (6-1) is the permittivity of GaAs.A comparison of the results for the variation of the surface concentration with Fermi level,as computed by the full self-consistent approach and the triangular potential well (TPW)approach has been published [75]. Only by introducing an empirical factor to modify thefield F can good agreement be achieved between the two approaches.Here we propose a new approach to the problem which avoids the computational complexities of the full, self-consistent, quantum approach, yet yields more accurate results thanthe TPW approach. The two novel features of the proposed model are:(1) the use of the WKB approximation to solve the Schrödinger equation for the boundstates in the potential well;E(x=O)-EN174(2) the use of a classical approach involving a modified effective density of states functionto solve Poisson’s equation for the potential energy profile in the well.The application of this approach to computing the surface charge density in a MODFETis described in this chapter. The method may also be useful in the analysis of other quantumwell devices.6.2 WKB ApproximationIn this section we briefly review the WKB approximation [76]. The WKB method is atruncation approximation of the exponential power expansion. The general 1-D Schrödingerequation has the formh2 a2Eb(x)=+ V(x)] (x). (6-2)where V(x) is the potential energy term, rn is the mass of an electron (or effective mass),iJ’ is the wavefunction, and finally, E is the eigenvalue of energy. It can be shown that thesolution of (6-2) can be expressed as an exponential power series of h, the reduced Planck’sconstant [77]:(x) = exp { ( hSn(x)) } (6-3)Where S, n = 0, 1,2,•. are functions to be determined by substituting (6-3) into (6-2)and comparing both sides of (6-2) with respect to the power of 1. The simplest WKBapproximation results from ignoring any term with power of h higher than 0 in the series,i.e.,h (6—4)This approximation is valid under the following two conditions:a) the potential energy V(x) should be a smooth function of x except for some specialdiscontinuity such as an infinite barrier [77];b) the region of interest should be large compared with the size of an atom (—‘1 A).In all available quantum devices, the above conditions are well satisfied.75V( x)Fig. 6.2 Turning points. In a classical system, a particle with energy E cannot gobeyond these two points (a x < b). In a quantum system, the wave functionscan be non-zero beyond these two points.For bound states, like the two-dimensional electron gas (2DEG) in MODFETs (in onespatial direction) [64], the approximation (6-4) can be rewritten as [77]=or)[hJXPH+] (6—5)where c and c’ are constants, a and a’ are phase factors determined by the boundaryconditions (for more details on this, see [77]or [76]), a and b are the so-called “turning points”(Fig. 6.2) where V(a) = V(b) = E, and E is the energy level. Function p(x) is given byp(x) = [2rn(E — V(x))]”2 (6 — 6)where m is the mass of the particles (carriers). The quantization condition depends on theboundary conditions. In the case shown in Fig. 6.2, V(x) has finite value and is continuous76Ea b xat both turning points. The quantization condition in this case is [77]b 1j p(x)dx= (i+)h, i=0,1,2,•• (6—7a)In our case, we consider the potential at the interface to be that of an infinite barrier, thiscauses the wavefunction to vanish at this boundary and to introduce an additional phasechange of ir/4. It can be shown that the energy levels are determined by the followingequation [76][78]:b 3j p(x)dx= (i+)irh, i=0,1,2,•• (6—7b)As a check on this WKB approach, let us consider a simple case: the triangular potentialwell.Assume that the triangular potential has the formV(x) = qF3x (6 — 8)where q is the value of electric charge for an electron and F8 is a constant which representsthe electrc field. Then from (6-6),p(x) = [2m(E — qFx)]”2 (6 — 9)Therefore the integraljP(x)dx = j[2m(E — qF8x)]2dx(2rnEj3I (6 — 10)-3rnqFwhere E, i = 0, 1,2,... are the energy levels for electrons. x = a is the normal turning pointwhere V(a) = E1. Substituting the above result into the quantization condition (6-7b), onehas(2mE)312—(. 3’ . ——1z + — z —3mqF5 \ 4,orh21/3 3 2/3= (_) [qFs7r ( + (6-11)This is precisely the same as can be obtained directly from an analytical solution to theSchrödinger equation for a triangular potential well [74].776.3 Model for the 2DEG in MODFETs6.3.1 The Potential Energy V(x)As an application of the WKB method, let us consider a simple model for the 2DEGin MODFETs. In this model, we consider only the relation between the Fermi level E1and the electron surface densityn3(Ef). Fig. 6.1 shows the profile of V(x) in GaAs. TheA1GaAs/GaAs interface is on the left-hand side (x < 0). For a given Fermi level Ef, thepotential energy V(x) is given by the solution of Poisson’s equation&V(x)q22 — —[p—n—ND] (6—12)dxwhere p is the density of holes, n is the density of electrons, and ND is the total dopingdensity Na — Nd.6.3.2 The Effective Density of StatesThe effect of the appearance of the subbands can be described within the classical modelby the concept of the effective density of states. Here we propose the following expressionfor n:n = NF112(7) (6 — 13)where N is the modified effective density of states for the conduction band, F112(q) is theFermi-Dirac integral of order one-half and= E1 —V(x) (6— 14)where kB is Boltzmamn’s constant. We relate N to the more usual parameter N, which isthe effective density of states in the classical limit (no subbands), byNc*•=aNc (6—15)where c is a factor which should be at least a function of Ef:= t(E) (6 — 16)78Moreover, because the subbands make the effective density of states less than in the trueclassical case, a satisfies the relation0<a<1 (6—17)The density of states for holes, on the other hand, stays intact at its dassical value, i.e.,p = NF112(—ii—(6— 18)where E9 is the band gap of GaAs and the Fermi-Dirac integral is:2 00= J 1 + (6— 19)In the calculation, the function F112 is approximated byFi12(i) = [e + (ii)]—’ (6 — 20)where= 3\/7[(’1 + 2.13) + (I — 2.1312.4 + 9.6)5/12]3h12 (6 — 21)This approximation has a maximum error of ±0.5 % for all values of [79].6.3.3 Determination of the factor cr(E1)To determine c(Ej), we assume that the surface density calculated from the classicaltreatment equals that calculated from the WKB method, i.e.,00J ndx=Enj (6—22)0 1=0where n, i = 0, 1,... 00 is the surface density at the ith energy level and can be calculatedfrom: [64]= mkT in {i + exp (Ef_EJ] (6— 23)where E1 i = 0, 1,” are the energy levels in the potential well. The left-hand side of (6-22)refers to the surface density as computed from the classical approach describe above. Theupper limit for the integral was taken as 0.3 tm, rather than infinity, as the conduction bandis essentially fiat at this point. n, was computed from edV/dxIo, making due allowance forthe doping density ND and the surface density of holes. The employment of (6-22) gives aself-consistency to our method of estimating the surface density of electrons.796.4 Method of Solution6.4.1 Boundary ConditionsWith reference to Fig. 6.1, the boundary conditions for solving (6-12) with a given valueof E1 are:V(x = 0) = given value of E(x = 0) — Ef (6 — 24)V(x—+oo)=E9—çbwhere is the value of Ef — E(x =(6—25)NDThe boundary conditions (6-24) are sufficient for solving (6-12) for V(x). However, in thecalculation, it is desirable that the derivative of V(x) with respect to x at x 0 is knownbeforehand. This is because the Runge-Kutta-Fehlberg method can then be used mostefficiently for obtaining V(x). Without knowing the exact form of V(x), one can calculatedV/dxl=o by changing the variable. Let(6—26)then one obtains from (6-12)dw qWJ-J=—{P—Ti—ND] (6—27)where p and n are given by (6-13) and (6-18), respectively, and are only functions of V.Integrating (6-27), one obtains2qkT (Eg—)/k8Tw2(x 0) =V(=O)-Ef)/kBT[ND + aNFi12(—i) — NF112(q — Eg/kBT)] d(6—28)It follows that the boundary condition dV/dx_0 can be calculated.6.4.2 The Program Used in the CalculationThe program used in the calculation, MODFET, is listed in the Appendix E. The mainprogram searches for the value of the factor a which satisfies the self-consistency condition(6-22). The method used here is the Fibonacci Search described in Chapter 5. The function80target in the main program is only a function of c and returns the square of the differencebetween the surface density calculated by both the classical and WKB methods. The functiontarget includes the subroutines field and level. The subroutine field uses the Runge-KuttaFehlberg method to solve (6-12), with the proper boundary conditions at z = 0; and thesubroutine level uses the WKB method to get the energy levels from (6-6) and (6-7b) usingthe V(x) calculated as just described. Other functions denhole and denelec calculate thecontributions of holes and electrons to the surface charge density.6.5 ResultsTo compare the results from our method with those of the full, self-consistent, quantum andTPW approaches reported in [75], we use the same parameters as Yoshida, namely,E9 = 1.425eVm*= 0.067mN = 4.350 x i0’’ cm3N,, = 9.313 x 1018 cm3 (6 — 29)e=12.9e0= 0.33eVNa — Nd = 1 x 1014 cm3where Na — Nd is the residual acceptor density in the nominally undoped GaAs and isthe conduction band discontinuity at the interface.The results are shown in Fig. 6.3, where the Fermi energy is taken with respect to thepotential energy at the bottom of the well (see Fig. 6.1). As can be clearly seen, the proposedWKB method yields results which are in excellent agreement with those from the full, self-consistent, quantum approach. As the computational effort involved in implementing theWKB approach is considerably less than that required by the full quantum approach, it812.50present result> /2.00 /Cfull treatment [75]/1.50- /approx.[75)0.00— — I I I 11111111 I I—0.10—0.00 0.10 0.20Fermi Lev& (eV)Fig. 6.3 n3 vs. E, from the present method, and from [75]. The surface density isin 10’2/cmappears that the new method may make a useful addition to the design tools used forquantum-well devices.To this end, we include a plot of a vs. E in Fig. 6.4. The best fit to this data is givenbya(Ej) = 0.6.54887— 0.3319.SEf— 13.214E + 266.776E(6—30)_-1663.83E + 3558.21Ewhere Ef is in eV and is measured with respect to the bottom of the well. By utilizingFig. 6.4, and so estimating N via (6-15), the problem of finding the surface charge densityreduces to that of the straightforward, common classical situation.820.800.600-c-0.400.200.00—0.10Fig. 6.4 The function u(Ej). The solid line is the fit (equation (6-30)) to thecomputed points.ConclusionsThe success of the proposed method of calculating the surface charge density in the potentialwell of a MODFET leads to three conclusions:1) the WKB method is highly appropriate for solving the Schrödinger equation for thebound states in the potential well;2) the use of a modified effective density of states function to allow a classical solution ofPoisson’s equation for the potential profile in the well is very reasonable;3) the combination of the WKB method for solving Schrödinger’s equation and a modifiedeffective density of states function to solve Poissons equation yields a procedure which isvery accurate but much less computationally-intensive than the full quantum approach.Fermi Leve’83Chapter 7. SummaryIn this thesis three novel approaches to modeling certain important phenomena in GaAselectronic devices have been developed. The specific conclusions drawn from the study ofGaAs HBTs, MESFETs, and MODFETs have been stated at the end of the respectivechapters, i.e., Chapter 2, 3, 4, 5 and 6.The main contribution has been the proposal of a phenomenological approach to carriertransport in GaAs, and the implementation of it in estimating transit times and signal delaytimes in the base and base-collector space charge regions of HBTs. The good agreementof the results with those from Monte Carlo methods, as published in the literature, suggestthat the phenomenological approach is accurate. The fact that it achieves this accuracy withconsiderably less computational effort than is needed for Monte Carlo calculations, suggeststhat the phenomenological method may be a useful design tool.The work on HBTs has:1) confirmed that HBTs with novel collector structures (np+p devices or n+p+in devices)have shorter base-collector space charge region transit times than conventional iip+structures [40]. This is important because transport through this space charge regionoften dominates the overall emitter-collector delay.2) predicted that the widely-used value of 0.5 for the ratio of signal delay time to transittime in the collector-base space charge region does not hold for GaAs HBTs [80]. This isbecause the assumption of uniform velocity, which leads to the value of 0.5, does not holdin GaAs devices due to the presence of velocity overshoot. The ratio for np+n structuresis around 0.4, whereas that for typical np+p structures can be as low as 0.3. Therefore,the latter structures (inverted field devices) should be even faster than previous analyseshave indicated.3) revealed that the principal cause of the fall-off in high-speed performance at high collectorcurrent densities is the low mobility of the holes. This new finding is significant becausea consideration of electron transport alone would suggest that the collector-base space84charge signal delay time should actually decrease with increasing J. This is becausethe neutralization of the collector donor charge by the injected electrons leads to aprogressively more inverted nature of the field. Our results show that this effect isswamped by the the behaviour of the holes, suggesting that one look for materials inwhich both electrons and holes have high mobility in order to maintain good speedperformance (e.g., high fT) at high collector current densities.Concerning MESFET devices, a contribution to the literature on the formation of stationary domains has been made [81]. As these domains can significantly affect the device’scapacitance, it is useful to know under what conditions they may form. Previously, thesame condition as for the formation of travelling domains had been widely-used. The workdescribed in Chapter 5 uses a phase portrait approach to indicate that the actual criterionfor the formation of stationary domains can be either higher or lower than the travellingdomain criterion, depending on the length of the gate of the device.Finally, this time in the context of MODFET devices, a new method for computing thesurface charge density in the potential well at the A1GaAs/GaAs interface has been developed. Knowledge of the surface charge density is useful as it provides a starting point forcalculations of the device’s current-voltage and capacitance-voltage characteristics. The proposed method invokes the WKB method to solve Schrödinger’s equation for the bound states.The potential energy profile on the GaAs side of the well is a required input to the waveequation. In Chapter 6 a quasi-classical approach is proposed to obtain this profile. Thisprocedure leads to a very rapid calculation of the surface charge density, and yields excellentagreement with published results from the full, quantum self-consistent approach. Thus, thenew method may prove to be a useful component in the design tool set for MODFET andother quantum well devices.858. Bibliography[1] P. M. Asbeck, M. F. Chang, K. C. Wang, D. L. Miller, G. J. Sullivan, N. H. Sheng,E. A. Sovero and J. A. Higgins, “Heterojunction bipolar transistors for microwave andmillimeter-wave integrated circuits,” IEEE 1987 Microwave and Millimeter- Wave Monolithic Circuits Symposium, pp. 1-5.[2] G.-W. Wang, Y.-K. Chen, W. J. Schaff, and L. F. Eastman, “A 0.1-nm gate Alo.sIno.sAs/Ga0•5In5As MODFET fabricated on GaAs substrates,” IEEE Trans. Electron Devices,vol. ED-35, pp. 818-823, 1988.[3] H. Q. Tserng and B. Kim, “110 GHz GaAs FET oscillator,” Electron. Lett., vol. 21, pp.258-259, 1985.[4] J. A. Adams, I. G. Thayne, N. I. Cameron, M. R. S. Taylor, S. P. Beaumont, C. D.Wilkinson, N. P. Johnson, A. H. Kean, C. R. Stanley, “Fabrication and high frequencycharacterisation of GaAs MESFET with gate length down to 30 nm,” Microelectron.Eng., vol. 11, pp. 65-68, 1990.[5] N. Nishiyama, H. Yano, S. Nakajima, and H. Hayashi, “n-AlInAs/(InAs)(Ga s)1superlattice modulation-doped field effect transistor grown by molecular beam epitaxy,”Electron. Lett., vol. 26, pp. 885-886, 1990.[6] C. Takano, K. Taira, and H. Kawai, “Improving collector-current uniformity in emitter-graded AlGaAs/GaAs heterojunction bipolar transistors,” IEEE Trans. Electron DeviceLett., vol. 9, pp. 125-127, 1988.[7] N. Hayama, M. Madihian, A. Okamoto, H. Toyoshima, and K. Honjo, “Fully self-alignedA1GaAs /GaAs heteroj unction bipolar transistor for high-speed integrated-circuits application,” IEEE Trans. Electron Devices, vol. 35, pp. 1771-1777, 1988.[8] P. I. Rockett, “Monte Carlo study of the influence of collector region velocity overshooton the high frequency performance of AlGaAs/GaAs heterojunction bipolar transistors,”IEEE Trans. Electron Devices, vol. ED-35, pp. 1573-1579, 1988.[9] J. Hu, K. Tomizawa and D. Pavlidis, “Monte Carlo approach to transient analysis ofHBT’s with different collector designs,” IEEE Trans. Electron Device Lett., vol. 10. pp.55-57, 1989.[10] R. Katoh, M. Kurata, and J. Yoshida, “Self-consistent particle simulation for (A1Ga)As/GaAs HBT’s with improved base-collector structures,” IEEE Trans. Electron Devices,vol. ED-36, pp. 846-853, 1989.[12] R. W .H. Engelmann and C. A. Liechti, “Gunn domain formation in the saturated currentregion of GaAs MESFETs,” in IEDM Tech. Dig., Dec. 1976, pp. 351-354.[13] K. Yamaguchi, S. Asai, and H. Kodera, “Two-dimensional numerical analysis of stabilityof GaAs FET’s,” IEEE Trans. Electron Devices, vol ED-23, pp. 1283-1290, 1976.86[14] T. Wada and J. Frey, “Physical basis of short-channel MESFET operation,” IEEE Trans.Electron Devices, vol ED-26, pp. 476-490, 1979. /[15] T. A. Fjeldly, “Analytical modeling of the stationary domain in GaAs MESFET’s” IEEETrans. Electron Devices, vol. ED-33, pp. 874-880, 1986.[161 M. S. Shur and L. F. Eastman, “Current-voltage characteristics, small-signal parameters,and switching times of GaAs MESFET’s,” IEEE Trans. Electron Devices, vol. ED-25,pp. 606-611, 1978.[17] T. A. Fjeldly and J. S. Johannessen, “Negative differential resistance in GaAs MESFETs,” Electron. Lett., vol. 19, pp. 649-650, 1983.[18] H. Kroemer, “Negative conductance in semiconductors,” IEEE Spectrum, vol. 5, p. 47-56,1968.[20] S. M. Sze, Physics of Semiconductor Devices, New York: Wiley, 1969, p. 282.[21] H. J. Kreuzer, Nonequilibrium Thermodynamics and its Statistical Foundations, Oxford,UK: Clarendon, 1981, p. 19.[22] J. A. Copeland, “LSA oscillator diode theory,” J. Appl. Phys., vol. 38, pp. 3096-3103,1967.[23] E. M. Conwell and M. 0. Vassell, “High-field distribution function in GaAs,” IEEETrans. Electron Devices, vol. ED-13, pp. 22-27, 1966.[24] J. M. Hyman, “MOL1D: A general purpose subroutine package for the numerical solutionof partial differential equations,” Manual LA-7595-M, Los Alamos Sci. Lab. (P.O. Box1663, Los Alamos NM 87545).[25] R. A. Smith, Semiconductors, 2nd. ed. London, UK: Cambridge Univ. Press, 1978, p.389.[26] D. E. Aspnes, “Gallium arsenide and related compounds,” in Inst. Phys. Conf. Ser., vol.33b, p. 100, 1977.[27] M. Shur, GaAs Devices and Circuits, New York, NY: Plenum, 1986, p. 39.[28] P. Das and R. Bharat, “Hot electron relaxation times in two-valley semiconductors andtheir effect on bulk-microwave oscillators,” Appl. Phys. Lett., vol. 11, pp. 386-388, 1967.[29] C. M. Maziar, M. E. Klausmeier-Brown, and M. S. Lundstrom, “A proposed tructurefor collector transit-time reduction in A1GaAs/GaAs bipolar transistors,” IEEE Trans.Electron Device Lett., vol. EDL-7, pp. 483-485, 1986.[30] T. Ishibashi and Y. Yamauchi, “A possible near-ballistic collection in an AlGaAs/GaAsHBT with a modified collector structure,” IEEE Trans. Electron Devices, vol. 35, pp.401-404, 1988.87[31] J. E. Carroll, Hot Electron Microwave Generators, London, UK: Edward Arnold, 1970,p. 61.[32] Y. Yamauchi and T. Ishibashi, “Electron velocity overshoot in the collector depletionlayer of A1GaAs/GaAs HBT’s,” IEEE Trans. Electron Device Lett., vol. EDL-7, pp 655-657, 1986.[33] J. M. Early, “P-N-I-P and N-P-I-N junction transistor triodes,” Bell Syst. Tech. J., vol.33, pp. 517-533, 1954.[34] F. N. Trofimenkoff, “Collector depletion region transit time,” Proc. IEEE, vol. 52, pp.86-87, 1964.[35] R. G. Meyer and R. S. Muller, “Charge-control analysis of the collector-base space-charge-region contribution to bipolar-transistor time constant TT,” IEEE Trans. ElectronDevices, vol. ED-34, pp. 450-452, 1987.[36] A. van der Ziel and T. G. M. Kleinpenning, “High-frequency response of microwavetransistors,” Solid-State Electron., vol. 30, pp. 771-772, 1987.[37] T. Ishibashi, “Influence of electron velocity overshoot on collector transit times ofHBT’s,” IEEE Trans. Electron Device, vol. 37, pp. 2103-2105, 1990.[38] S. E. Laux and W. Lee, “Collector signal delay in the presence of velocity overshoot,”IEEE Trans. Electron Device Lett., vol. EDL-11, pp. 174-176, 1990.[39] D. R. Hamilton, J. K. Knipp and J. B. H. Kuper, Klystrons and Microwave Triodes,(MIT Radiation Laboratory Series), New York: McGrawHill, 1948, ch. 3.[40] H. Zhou, D. L. Puifrey and M. J. Yedlin, “A phenomenological approach to estimatingtransit times in GaAs HBTs,” IEEE Trans. Electron Devices, ED-37, pp. 2113-2120,1990.[41] D. L. Puifrey and N. G. Tarr, introduction to Microelectronic Devices, New Jersay:Prentice-Hall, 1989, pp. 382-383.[42] J. Shibata, I. Nakao, Y. Sasai, S., Kimura, N. Hase, and H. Serizawa, “Monolithic integration of an InGaAsP/InP laser diode with heterojunction bipolar transistors,” Appl.Phys. Lett., vol. 45, pp. 191-193, 1984.[43] K.-Y. Liou, S. Chandrasekhar, A. G. Dentai, E. C. Burrows, G. J. Qua, C. FL. Joyner,and C. A. Burrus, “A 5 Gb/s monolithically integrated lightwave transmitter with 1.5zm multiple quantum well laser and HBT driver circuit,” IEEE Photonics Technol. Lett.,vol. 3, pp. 928-930, 1991.[44] H. C. Poon, H. K. Gummel, and D. L. Scharfetter, “High injection in epitaxial transistors,” IEEE Trans. Electron Devices, vol. ED-16, pp. 455-457, 1969.88[45] J. J. Liou, F. A. Lindholm and B. S. Wu, “Modeling the cut-off frequency of singleheterojunction bipolar transistors subjected to high collector-layer current,” J. Appl.Phys., vol. 67, pp. 7125-7131, 1990.[46] R. M. Warner Jr. and B. L. Grung, “Transistors: Fundamentals for the Integrated-CircuitEngineer,” Florida: Robert Kreiger Inc., 2nd Edn., 1990, pp. 626-667.[47] Z.-P. Yu and R. W. Dutton, SEDAN III - A Generalized Electronic Material DeviceAnalysis Program, Integrated Circuits Laboratory, Stanford University, Stanford, CA94305, 1985.[48] M. R. Pinto, C. S. Rafferty, and R. W. Dutton, PISCES-Il User’s Manual, StanfordElectronics Laboratories, Department of Electrical Engineering, Stanford University, Stanford CA 94305, 1984.[49] UBC MATRIX, edited by Tom Nicol, Computing Centre, The University of BritishColumbia, 6356 Agricultural Road, Vancouver, B.C., Canada, V6T 1W5.[50] K. H. Huebner and E. A. Thornton, The Finite Element Method For Engineers, NewYork: John Wiley & Sons, 1982, Ch. 2,4,5,10, and 11.[51] R. W. H. Engelmann and C. A. Liechti, “Bias dependence of GaAs and InP MESFETparameters,” IEEE Trans. Electron Devices, vol. ED-24, pp. 1288-1296, 1977.[52] G. S. Kino and P. N. Robson, “The effect of small transverse dimensions on the operationof Gunn devices,” Proc. IEEE, vol. 56, pp. 2056-2057, 1968.[53] D. K. Arrowsmith and C. M. Place, Ordinary Differential Equations, London: Chapmanand Hall, 1982, chs. 1 and 2.54] K. W. Boer and G. Dohler, “Influence of boundary conditions on high-field domains inGunn-diodes,” Phys. Rev., vol. 186, pp. 793-800, 1969.[55] E. M. Conwell, “Boundary conditions and high-field domains in GaAs,” IEEE Trans.Electron Devices, vol. ED-17, pp. 262-270, 1970.[56] P. Gueret and M. Reiser, “Switching behavior of over-critically doped Gunn diodes,”Appl. Phys. Lett., vol. 20, pp. 60-62, 1972.[57] Richard L. Branham, Jr., Scientific Data Analysis, New York: Springer-Verlag, 1990, pp.175- 177.[58] J. A. Copeland, “Stable space charge layers in two-valley semiconductors,” J. Appi.Phys., vol 37, pp.3602-3609, 1966.[59] B. L. Gelmont and M. S. Shur, “Analytical theory of stable domains in high-doped Gunndiodes,” Electron. Lett., vol. 6, pp. 385-387, 1970.[60] H. Flanders, Scientific PASCAL, Reston: Reston Publishing, 1984, pp. 273-283.89[61] J. H. Abeles, R. F. Leheny, G. K. Chang, and S. J. Allen, “Experimental measurementof a high-field dipole in GaAs field-effect transistors,” Appi. Phys. Lett., vol. 49, pp.1387-1398, 1987.[62] R. W. McColl, R. L. Carter, J. M. Owens, and T.-J. Shieh, “GaAs MESFET simulationusing PISCES with field-dependent mobility-diffusion relation,” IEEE Trans. ElectronDevices, vol. ED-34, pp. 2034-2039, 1987.[63] S. Kugler, “Generation-Recombination Noise in the Saturation Regime of MODFETstructures,” IEEE Trans. Electron Devices, vol. ED-34, pp. 623-628, 1988.[64] D. Delagebeaudeuf and N. T. Linh, “Metal-(n) A1GaAs-GaAs two-dimensional electrongas FET,” IEEE Trans. Electron Devices, vol. ED-29, pp. 955-960, 1982.[65] M. Shur, GaAs Devices and Circuits, New York, NY: Plenum, 1986, p. 545.[66] Y. Ando and T. Itoh, “Analysis of charge control in pseudomarphic two-dimensionalelectron gas field-effect transistors,” IEEE Trans. Electron Devices, ED-35, pp. 2295-2301, 1988.[67] B. Vinter, “Subbands and charge control in a two-dimensional electron gas field-effecttransistor,” Appi. Phys. Lett., vol. 44, pp. 307-309, 1984.[68] S. Voinigescu and A. Muller, “Charge dynamics in heterostructure Schottky-gate capacitors and their influence on the transconductance and low-frequency capacitance ofMODFET’s,” IEEE Trans. Electron Devices, vol. ED-36, pp. 2320-2327, 1989.[69] T. Ando, “Self-consistent results for a GaAs/A1Gai_As heterojunction. I. Subbandstructure and light-scattering spectra,” J. Phys. Soc. Japan., vol. 5, pp. 3893-3899, 1982.[70] S. P. Svensson, “Theoretical analysis of the layer design of inverted single-channel heterostructure transistors,” IEEE Trans. Electron Devices, vol. ED-35, pp. 992-998, 1987.[71] W. A. Hughes and C. M. Snowden, “Nonlinear charge control in A1GaAs/GaAsmodulation-doped FET’s,” IEEE Trans. Electron Devices, vol. ED-34, pp. 1617-1625,1987.[72] M. L. Majewski, “An analytical dc model for the modulation-doped field-effect transistor,” IEEE Trans. Electron Devices, vol. ED-34, pp. 1902-1910, 1987.[73] A. Eskandarian, “Determination of the small-signal parameters of an A1GaAs/GaAsMODFET,” IEEE Trans. Electron Devices, vol. ED-35, pp. 1793-1801, 1988.[74] F. Stern, “Self-consistent results for n-type Si inversion layers,” Phys. Rev., B5, pp.4891-4899, 1972.[75] J. Yoshida, “Classical versus quantum mechanical calculation of the electron distributionat the n-A1GaAs/GaAs heterointerface,” IEEE Trans. Electron Devices, vol. ED-33, pp.154-156, 1986.90[76] Leonard I. Schiff, Quantum Mechanics, New York: McGraw-Hill Book Company, 1968,pp.268-277.[77] David Bolim, Quantum Theory, New York: Prentice-Hall, Inc., 1951, ch.12.[78] A. B. Migdal and V. Krainov, Approximation Methods in Quantum Mechanics, NewYork: W. A. Benjamin, Inc., 1969, p. 121.[79] Robert F. Pierret, Advanced Semiconductor Fundamentals, Reading, Massachusetts:Addison-Wesley Publishing Company, 1987, p. 118.[80] H. Zhou and D. L. Puifrey, “Computation of transit and signal delay times for thecollector depletion region of GaAs-based HBTs,” Solid-State Electron., accepted in July,1991.[81] H. Zhou and D. L. Puifrey, “A criterion for stationary domain formation in GaAs MESFET’s,” IEEE Trans. Electron Devices, vol. ED-36, pp. 872-878, 1989.[82] J. Hu, D. Pehike, K. Tomizawa and D. Pavildis, “Heterojunction bipolar transistor studies,” Bell Northern Research Report, p.53, 1988.91Appendix A Listing of Programs Used to Estimate Transit Times in HBTsA.1 List of PULSHBTC PROGRAM PULSHBTCIMPLICIT REAL*8(A—H,O-Z)DIMENSION UZ(2,1O1),XM(1O1),MORD(2,3),TOUT(30)DIMENSION CFE(100) ,CFN(100) ,CFX(100)COMMON/MOLPLT/XAL,NPLT,MG,UAL,XAXIS(2),UAXIS(2,2)COMMQN/CONST/AMU1 ,AMU2,EO,D1 ,D2COMMON/BLOCK/TAU,VS ,T0COMMON/CONF/CFE, CFN , CFX , NCONREAD(5,2)TINT,TLAST,TAUWRITE(6 , 1)TINT,TLAST,TAU1 FORMAT(’ TINT= ‘,D10.3,’ TLAST= ‘,D1O.3,/ ‘ TAU2= ‘,D10.3,/)2 FORMAT(3D10.3)READ (5 , 3) NCON3 FORMAT(14)DO 4 I=1,NCONREAD(5,5)CFX(I) ,CFE(I),CFN(I)WRITE(6,5)CFX(I) ,CFE(I) ,CFN(I)5 FORMAT(3D10.3)4 CONTINUEXLEFT=CFX(1)XRIGHT=CFX (NCON)E0=3 4D-2D1=1 554D00D2=O 751D00AMU1=60 ODOOAMU2=29.ODOOVS=0 6D00T0=0. 1ODOONPDE=2NPTS= 101KEQN=2KBC=2DO 10 I=1,NPDEMORD(I,1)4MORD(I,2)=410 MORD(I,3)=4METH=22EPS=1 OD—lOMOUT=1KMOL=ODT=(TLAST-TINT)/30 ODOODO 20 1=1,30TOUT( I) 1*DT+TINT20 CONTINUEDX (xRI GHT-XLEFT) /DFLOAT (NPTS- 1)DO 50 11,NPTSXM (I) =DFLOAT (I—i) *DX+XLEFTUz(1,I)=0.oDoo92so Uz(2,I)=0.ODO0NPLT=214G=—1XAL=10 ODOOUAL=10 . 01300XAXIS(1)=XM(1)XAXIs(2)xM(NPTs)DO 70 I=1,NPDEUAXIS(I,1)=O.ODOOUAXIS(I,2)=1 .000070 CONTINUECALL MOL1D(NPDE,NPTS ,KEQN,KBC ,METH,EPS,MORD,/ TINT,TLAST,MOUT,TOUT,UZ,XM,KMOL)STOPENDCSUBROUTINE PDE(UT,U,UX,UXX,FX,T,X,IX,NPDE)IMPLICIT REAL*8(A—H , O-z)DIMENSION U(NPDE) ,UT(NPDE) ,UX(NPDE) ,UXX(NPDE) ,FX(NPDE)COMMON/CONST/AMU1 , AMU2 , EQ .01 , D2COMMON/BLOCK/TAU,VS ,T0DATA TNAX/1.0D6/DATA QE/0. 828792522000/TEMP 1FIELD(X)TEMP2=TEMP1/E0TEMP2=TEMP2*TEMP2*TEMP2*TEMP2IF(TEMP2.LT.TMAX) GO TO 10TEMP2TMAX10 TEMP3(U(2)—U(1)*TEMP2)/TAUUT(1)=—FX(1)+0. 1DOO*UXX(1)/ +TEMP3UT (2 ) —FX (2 ) +0 . 5D—0 1*UXX (2)/ —TEMP3RETURNENDCSUBROUTINE FUNC(F,U,UX,UXX,T,XM,IX,NPDE)IMPLICIT REAL*8 (A—H, O—z)DIMENSION F(NPDE) ,U(NPDE) ,UX(NPDE) ,UXX(NPDE)COMMON/CONST/AMU1 , AMU2 , E0 ,D1 ,D2DATA CKTQ/0. 0259000/V1=VDRI1(XM)V2=VDRI2(XM)AD 1=D 1AD2=D2E=FIELD (XM)IF(E.EQ.0.ODOO) GO TO 10AD1=—V1*CKTQ/EAD2=—V2*CKTQ/E10 F(1)=—AD1*UX(1)+U(1)*VDRI1(XM)+0.1DOO*UX(1)F(2)=—AD2*UX(2)+U(2)*VDRI2(XM)+0. 5D—01*UX(2)RETURNEND93CSUBROUTINE BNDRY(T,UL,AL,BL,CL,UR,AR,BR,CR,NPDE)IMPLICIT REAL*8 (A-B, O—Z)DIMENSION UL(NPDE),AL(NPDE),BL(NPDE),CL(NPDE),UR(NPDE),/ AR(NPDE) ,BR(NPDE) ,CR(NPDE)COMMON/CONST/AMU1 , AMU2,EQ ,D1 , D2COMMON/BLOCI{/TAU,VS,TODATA P1/3.141592664/DO 10 1=1,2AL(I)0.ODOOBL(I)=O.ODOOCL(I)=O.ODOOAR(I)=0.ODOOBR(I)0.0000CR(I)=O.ODOO10 CONTINUEAL(1)=1 .0000AL(2)1 .0000IF(T.GT.T0) GO TO 20CL(1)=0.SDOO*(1.ODOO—DCOS(2.0000*PI*T/TO))20 AR(1)=1.ODOOAR(2)1.0000RETURNENDCFUNCTION FIELD(X)IMPLICIT REAL*8 (A-H,O--Z)DIMENSION CFE(100) ,CFN(100) ,CFX(100)COMNON/CONST/AMU1 , AMU2 ,E0 ,D1 ,D2COMMON/BLOCK/TAU ,VS ,TOCOMMON/CONF/CFE, CFN • CFX , NCONDO 10 I=1,NCONIF(X.LT.CFX(I)) GO TO 2010 CONTINUE20 R=(CFE(I)—CFE(I—1))/(CFX(I)—CFX(I—1))FIELD=CFE(I—1)+R*(X—CFX(I—1))RETURNENDCFUNCTION VDRI2(X)IMPLICIT REAL*8 (A—H,O—Z)COMMON/CONST/AMU1 ,AMU2,EO,D1 ,D2COMMON/BLOCK/TAU,VS,TOTEMP1=FIELD(X)BX=ANU2*TEMP 1IF(DABS(BX).GT.(1.5DOO*VS)) GO TO 10TEMP2=BX/VSVDRI2=—BX* (1 . 0000—TEMP2*TEMP2*4 . ODOO/27. ODOO)RETURN10 VDRI2=—VS*TEMP1/DABS(TEMP1)RETURNENDC94FUNCTION VDRI1(X)IMPLICIT REAL*8 (A—H,O—z)COMI4ON/CONST/AMU1 , AMU2, E0 ,D1 ,D2COMMON/BLOCK/TAUVS ,TOV1=10. 0D00*VSTEMP 1FIELD(X)BX=ANU 1 *TEMP 1IF(DABS(BX).GT.(1.5D00*V1)) GO TO 10TEMP2=BX/V1VDRI1=—BX*(1 .0000—TEMP2*TEMP2*4 .ODOO/27 . 0000)RETURN10 VDRI1=—V1*TEMP1/DABS(TEMP1)RETURNEN])CC A sample of input file. It defines T00.OPS, T16.OPS, TAU2SONS,C AND A UNIFORM FIELD PROFILE.CC+O. 0000+00+1.6000+01+5.0000+04C6,C—i. 00000+00+0.00000+00+0. 00000+00C—O. 10000+00+0. 00000+00+0 . 00000+00C+0 .00000+00—0.75000+00+0. 0000D+O0C+5. 0000D+00—0 . 75000+00+0 . 00000+00C-i-S. 10000+00+0. 0000D+00+0 . 00000+00C+6.S0000+00+0.00000+00+0.00000+0095A.2 List of TRNTCIMPLICIT REAL*8(A—H ,O—z)DIMENSION CFE(100),CFN(100),CFX(100)COMMON/MOLPLT/XAL,NPLT,MG,UAL,XAXIS(2),UAXIS(2,2)COMMON/CONST/AMU1 , AMU2 , EO ,D1, D2COMMON/BLOCK/flu ,VS ,TO, TMIDCOMZ1ON/CONF/CFE,CFN,CFX,IJCONREAD(3 , 2)TINT,TLAST,TAUWRITE(6, 1)TINT,TLAST,TAU1 FORMAT(’ TINT= ‘,D1O.3,’ TLAST= ‘,D1O.3,/ ‘ TAU2= ‘,D1O.3,/)2 FORMAT(3D1O.3)READ (3 , 3) NCON3 FORMAT(I4)DO 4 11,NCONREAD(3,5)CFX(I),CFE(I),CFN(I)WRITE(6,5)CFX(I) ,CFE(I) ,CFN(I)5 FORMAT(3D10.3)4 CONTINUEXLEFT=CFX (1)XRIGHT=CFX(NCON)EO=3 40—2D1=1.554DO0D2=O 751000ANU1=60. ODOOAMU2=29 0000VS=1 .0000CALL VELOSTOPENDCFUNCTION FIELD(X)IMPLICIT REAL*8 (A—H,O—Z)DIMENSION CFE(100),CFN(100),CFX(100)COMMON/CONST/AMU1,AMU2,EO,D1 ,D2COMMON/BLOCK/TAU,VS,T0,TMIDCOMMON/CONF/CFE, CFN, CFX , NCONDO 10 I=1,NCONIF(X.LT.CFX(I)) GO TO 2010 CONTINUE20 R=(CFE(I)—CFE(I—1))/(CFX(I)--CFX(I—1))FIELD=CFE(I—1)+R*(X—CFX(I—1))RETURNENDCFUNCTION VDRI2(X)IMPLICIT REAL*8 (A—H,O—z)COMMON/CONST/AMU1 , AMU2 , EO ,Dl , 02COMMON/BLOCK/TAU ,VS ,TO, TMIDTEMP1=FIELD(X)BX=AMU2*TEMP 1IF(DABS(BX).GT.(1.S000*VS)) GO TO 10TEMP2=BX/VS96VDRI2=-BX*(1 ODOO—TEMP2*TEMP2*4 ODOO/27 ODOO)RETURN10 VDRI2=—VS*TEMP1/DABS(TEMP1)RETURNENDCFUNCTION VDRI1(X)IMPLICIT REAL*8 (A—H,o—z)COMMON/CONST/AMU1 , AMU2 , EO ,D1 ,02COMMON/BLOCK/TAU ,VS ,TO, TMIDV110 ODOO*VSTEMP1=FIELD(X)BX=AMU 1 *TEMP 1IF(DABS(BX).GT.(1.5DOO*V1)) GO TO 10TEMP2=BX/V1VDRI1=—BX* (1 ODOO—TEMP2*TEMP2*4 ODOO/27 ODOO)RETURN10 VDRI1=—V1*TEMP1/DABS(TEMP1)RETURNENDCSUBROUTINE VELOIMPLICIT REAL*8 (A—H, O—z)COMMON/CONST/AMU1 , AMU2 , EO ,Dl , D2DATA CKTQ/0 0269D00/DIMENSION X(151),AN(151),ANT(151),AT(151),U1(1S1),U2(1S1)NSET=30WRITE(6, 100)100 FORMAT(’INPUT Ml: ( FORMAT 13 ) ‘)READ(5,25)M1DO 110 I=1,M1READ(4, 120)110 CONTINUE120 FORNAT()LNUM=100DO 200 I=1,LNUMAN(I)=0.OEOOANT(I)=O.OEOO200 CONTINUE26 FORMAT(13)DO 20 I=1,NSETIF((I.EQ.1).OR.(I.EQ.NSET)) GOTO 210F=1.OEOOGO TO 220210 F=0.SEOO220 READ(4,60)T60 FORMAT(20X,G11.4)DO 230 K=1,6READ (4, 120)230 CONTINUEDO 300 J=1,LNUMREAD(4,70)NUM,X(J) ,U1(J) ,U2(3)300 CONTINUE97D=X(2)—X(1)DO 400 J=1,LNUMCC CALCULATE DERIV.CIF(J.EQ.LNUM) GO TO 430IFCJ.EQ.1) GO TO 410UX1=O.5*(U1((J+1))—U1((J--1)))/DUX2=O.5*(U2((J+1))—U2((J—1)))/DGO TO 420410 Uxl=(ul(2)—U1(1))/DUX2=(U2(2)—U2(1))/DGO TO 420430 Ux1=(U1(LNUM)—U1((LNUM—1)))/DUX2=(U2(LNUM)—U2((LNUM—1)))/DCC DRIFT VELOCITYC420 V1VDRI1(X(3))V2=VDRI2(X(J))AD 1=D 1AD2D2E=FIELD(X(i))IF(E.E.0.0DOO) GO TO 350AD1-V1*CKTQ/EAD2-V2*CKTQ/E350 AN(J)=AN(J)+F*(—AD1*UX1+V1*U1(J))+F*(—AD2*UX2-f-V2*U2(J))ANT(J)=ANT(J)+F*(—AD1*UX1+V1*U1(J))*T+F*(—AD2*UX2+V2*U2(J))*T400 CONTINUEMTEMP=LNUM+ 1DO 450 J=MTEMP,101READ (4, 460)450 CONTINUE460 FORMAT()70 FORMAT(14,D10.3,2(1X,D11.4))20 CONTINUEAT(1)=0.OEOODO 500 I=2,LNUNAT(I)=ANT(I)/AN(I)500 CONTINUEDO 600 I=1,LNUMWRITE(7,30)X(I),AT(I)600 CONTINUE30 FORMAT(2(1X,D11.4))RETURNEND98Appendix B Listing of Programs Used to Estimate Signal Delay Times in HBTsB.1 List of DELAYC PROGRAM DELAYCIMPLICIT REAL*8 (A-il, O—Z)DIMENSION Uz(2,1O1),XM(1o1),MORD(2,3),TOUT(3o)DIMENSION CFE(100) ,CFN(100) ,CFX(100)COMMON/MOLPLT/XAL,NPLT,MG,UAL,XAXIS(2),UAXIS(2,2)COMMON/CONST/AMU1 , AMU2 , EO , Dl, D2COMMON/BLOCK/TAU ,VS ,TOCOMNON/CONF/CFE, CFN, CFX • NCONREAD(5 ,2)TINT,TLAST,TAUWRITE(6, 1)TINT,TLAST,TAU1 FORMAT(’ TINT= ‘,D1O.3,’ TLAST= ‘,D1O.3,/ ‘ TAU2= ‘,D10.3,/)2 FORMAT(3D1O.3)READ (5 , 3) NCON3 FORMAT(I4)DO 4 I=1,NCONREAD(5,5)CFX(I) ,CFE(I) ,CFN(I)WRITE(6,5)CFX(I) ,CFE(I) ,CFN(I)5 FORMAT(3D1O.3)4 CONTINUEXLEFTCFX (1)XRIGHT=CFX(NCON)EO=3.4D-2D1=1 554D00D2=O 751D00AMU16O ODOOAMU229 ODOOVS=1 ODOOTO=4. 00000NPDE=2NPTS= 101KEQN=2KBC=2DO 10 I=1,NPDEMORD (I, 1) =4MORD(I,2)=410 MORD(I,3)=4METH=22EPS=1 OD—lOMOUT=1KMOL=ODT=(TLAST—TINT)/30 .0000DO 20 1=1,30TOUT(I)=I*DT+TINT20 CONTINUEDX (XRIGHT—XLEFT ) /DFLOAT (NPTS- 1)DO 50 I=1,NPTSXM(I)=DFLOAT(I-1)*DX+XLEFTUz(1,I)=o.oDoo99so Uz(2,I)=O.oDOONPLT=2MG=—1XAL=1O 0DO0UAL1O. ODOOXAXIS(1)=XM(1)XAXIS(2)=XM(NPTS)DO 70 I=1,NPDEUAXIS(I,1)=O.ODOOUAXIS(I,2)=1 .ODOO70 CONTINUECALL MOL1D(NPDE,NPTS ,KEQN ,KBC ,METH ,EPS,MORD,/ TINT,TLAST,MOUT,TOUT,UZ,XM,KMOL)STOPENDCSUBROUTINE PDE(UT,U,UX,UXX,FX,T,X,IX,NPDE)IMPLICIT REAL*8 (A-H, O-Z)DIMENSION U(NPDE),UT(NPDE),UX(NPDE),UXX(NPDE),FX(NPDE)COMMON/CONST/AMU1 , AMU2 , E0 ,D1 ,02COMMON/BLOCK/TAU ,VS ,TODATA TMAX/1.0D6/DATA QE/0.828792522D00/TEMP 1=FIELD(X)TEMP2=TEMP1/E0TEMP2=TEMP2*TEMP2*TEMP2*TEMP2IF(TEMP2.LT.TNAX) GO TO 10TEMP2=TMAX10 TEMP3=(U(2)—U(1)*TEMP2)/TAUUT(1)=—FX(1)+O. 1DOO*UXX(1)/ +TEMP3—F1(X)*QE*(U(1)+U(2))UT(2)=—FX(2)+O. SD—O1*UXx(2)/ -TEMP3RETURNENDCSUBROUTINE FUNC(F,U,UX,UXX,T,XM,IX,NPDE)IMPLICIT REAL*8(A—H ,O-Z)DIMENSION F(NPDE) ,U(NPDE) ,UX(NPDE) ,UXX(NPDE)COMMON/CONST/AMU1 , AMU2 , EO ,Dl,D2DATA CKTq/O 0259D00/V1=VDRI1(XM)V2=VDRI2(XM)AD1=D1AD2=D2E=FIELD (xM)IF(E.EQ.O.ODOO) GO TO 10AD1=-V1*CKTQ/EAD2—V2*CKTQ/E10 F(1)—AD1*UX(1)+U(1)*VDRI1(XM)+O. 1DOO*UX(1)F(2)=—AD2*UX(2)+U(2)*VDRI2(XM)+O SD—O1*UX(2)RETURNEND100CSUBROUTINE BNDRY(T,UL,AL,BL,CL,UR,AR,BR,CR,NPDE)IMPLICIT REAL*8(A—H , O—Z)DIMENSION UL(NPDE) ,AL(NPDE) ,BL(NPDE) ,CL(NPDE) ,UR(NPDE),/ AR(NPDE) ,BR(NPDE) ,CR(NPDE)CONMON/CONST/AMU1 , AMU2 ,E0 ,D1 ,D2COMMON/BLOCK/TAU ,VS • TODATA P1/3.141592654/DO 10 1=1,2AL(I)0.ODOOBL(I)=O.ODOOCL(I)0.ODOOAR(I)0.ODOOBR(I)0.ODOOCR(I)0.ODOO10 CONTINUEAL(1)1 .ODOOAL(2)1 .ODOOC IF(T.GT.TO) GO TO 20C CL(1)0.SDOO*(1.ODOO—DCOS(2.ODOO*PI*T/T0))CL(1)0.2D00*(1.ODOO—DCOS(2.ODOO*PI*T/T0))20 AR(1)=1.ODOOAR(2)1 .ODOORETURNENDCFUNCTION FIELD(X)IMPLICIT REAL*8 (A-H,O-Z)DIMENSION CFE(100),CFN(100),CFX(100)COMMON/CONST/AMU1 , AMU2 ,EO ,D1 ,D2COMMON/BLOCK/TAU ,VS ,TOCOMMON/CONF/CFE, CFN , CFX , NCONDO 10 I=1,NCONIF(X.LT.CFX(I)) GO TO 2010 CONTINUE20 R=(CFE(I)—CFE(I—1))/(CFX(I)--CFX(I--1))FIELD=CFE(I—1)+R*(X—CFX(I—1))RETURNENDCFUNCTION VDRI2(X)IMPLICIT REAL*8 (A-H,O—Z)COMMON/CONST/AMU1 , ANU2, E0 ,D1 ,D2COMMON/BL.OCK/TAU ,VS ,T0TEMP1=FIELD(X)BXAMU2*TEMP 1IF(DABS(BX).GT.(1.5D00*VS)) GO TO 10TEMP2BX/VSVDRI2=—BX* (1 . ODOO—TEMP2*TEMP2*4. ODOO/27 . ODOO)RETURN10 VDRI2=—VS*TEMP1/DABS(TEMP1)RETURNEND101CFUNCTION VDRII(X)IMPLICIT REAL*8 (A-H,O-Z)COMMQN/CONST/AMU1 ,AMU2,EO ,D1 ,D2COMMON/BLOCK/TAU,VS ,TOV1=1O ODOO*VSTEMP1=FIELD(X)BXAMU1*TEMP1IF(DABS(BX).GT.(1.S000*V1)) GO TO 10TEMP2=BX/V1VDRI 1=—BX*( 1 0000—TEMP2*TEMP2*4 0000/27.0000)RETURN10 VDRI1=—V1*TEMP1/DABS(TEMP1)RETURNENDCFUNCTION F1(X)IMPLICIT REAL*8 (A-H,o-z)DIMENSION CFE(100),CFN(100),CFX(100)COMMON/CONST/AMU1 ,AMU2,E0,D1 ,D2COMMON/BLOCK/TAU ,VS ,TOCOMMON/CONF/CFE, CFN , CFX , NCONDO 10 I=1,NCONIF(X.LT.CFX(I)) GO TO 2010 CONTINUE20 F1=CFN(I—1)RETURNEND102B.2 List of DRIVHBTC PROGRAM DRIVHBTCIMPLICIT REAL*8 (A—H, O-Z)DIMENSION CFE(100) ,CFN(100) ,CFX(100)COMNON/MOLPLT/XAL,NPLT,MG,UAL,XAXIS(2),UAXIS(2,2)COMMON/CONST/AMU1 , AMU2 • EO ,D1 , D2COMMON/BLOCK/TAU ,VS ,TO, TMIDCOMMON/CONF/CFE, CFN , CFX , NCONREAD(3,2)TINT,TLAST,TAUWRITE(6, 1)TINT,TLAST,TAU1 FORMAT(’ TINT= ‘,D1O.3,’ TLAST= ‘,D1O.3,/ ‘ TAU2= ‘,D1O.3,/)2 FORMAT(3D10.3)READ(3,3)NCON3 FORMAT(14)DO 4 I=1,NCONREAD(3,S)CFX(I) ,CFE(I),CFN(I)WRITE(6,5)CFX(I) ,CFE(I) ,CFN(I)S FORMAT(3D10.3)4 CONTINUEXLEFT=CFX (1)XRIGHT=CFX (NCON)E0=3.4D—2D1=1 554D00D2=O 751D00AMU1=60 0000AMU2=29 ODOOVS=1 ODOOCALL VELOSTOPENDCFUNCTION FIELD(X)IMPLICIT REAL*8 (A—H,O—Z)DIMENSION CFE(100) ,CFN(100) ,CFX(100)COMMON/CONST/AMU1 , AMU2 ,EO ,D1 ,D2COMMON/BLOCK/TAB ,VS ,TO, TMIDCOMNON/CONF/CFE , CFN , CFX , NCQNDO 10 I=1,NCONIF(X.LT.CFX(I)) GO TO 2010 CONTINUE20 R=(CFE(I)—CFE(I—1))/(CFX(I)—CFX(I—1))FIELD=CFE(I—1)+R*(X—CFX(I—1))RETURNENDCFUNCTION VDRI2(X)IMPLICIT REAL*8 (A-H,o—Z)COMMON/CONST/AMU1 , AMU2 , EO ,Dl , 02COMMON/BLOCK/TAU,VS ,T0,TMIDTEMP 1=FIELD(X)BX=AMU2*TEMP 1103IF(DABS(BX).GT.(1.5D00*VS)) GO TO 10TEMP2=BX/VSVDRI2=—BX* (1 . ODOO—TEMP2*TEMP2*4 . ODOO/27 ODOO)RETURN10 VDRI2=—VS*TEMP1/DABS(TEMP1)RETURNENDCFUNCTION VDRI1(X)IMPLICIT REAL*8 (A—li,O—Z)COMMON/CONST/AMU1 ,AMU2,E0,D1 ,D2COMMON/BLOCK/Thu ,VS ,TO , TMIDV110 0000*VSTEMP1=FIELD(X)BX=AMU1 *TEMP 1IF(DABS(BX) GT (1 SDOO*V1)) GO TO 10TEMP2=BX/V1UDRI 1=—BX* (1 ODOO—TEMP2*TEMP2*4 ODOO/27 ODOO)RETURN10 VDRI1=—V1*TEMP1/DABS(TEMP1)RETURNENDCSUBROUTINE VELOIMPLICIT REAL*8(A—H,O-Z)COMNON/CONST/AMU1 , AMU2, EO ,D1 , D2DATA CKTQ/O. 0259D00/DIMENSION X(151),U1(1S1),U2(151)DIMENSION TP1(30),TP2(30),TP3(30),TT(30)NSET=30WRITE(6, 100)100 FOR!4AT(’INPUT 141, 142, AND L: ( FORMAT 313 ) ‘)READ(S,25)M1 ,M2LNUMDO 110 I1,M1READ (4 , 120)110 CONTINUE120 FORMAT()25 FORMAT(313)DO 20 11,NSETREAD(4,60)TT(I)60 FORMAT(20X,G11.4)DO 230 K=1,6READ(4, 120)230 CONTINUEDO 240 K=1,M2READ(4, 120)240 CONTINUEDO 300 J=1.LNUMREAD(4,TO)NUM,X(J) ,U1(J) ,u2(J)300 CONTINUED=X(2)—X(1)AN=O ODOODO 400 J=1,LNUM104IF((I.EQ.1).OR.(I.EQ.L.NUM)) GO TO 210F=1.OEOOGO TO 220210 F0.5E00CC CALCULATE DERIV.C220 IF(J.EQ.LNUM) GO TO 430IF(J.EQ.1) GO TO 410Ux1=0.S*(U1((3+1))—U1((J—1)))/DUX2=0.S*(U2((3+1))—U2((J—1)))/DGO TO 420410 UX1=(U1(2)—U1(1))/DUX2=(U2(2)—U2(1))/DGO TO 420430 UX1=(U1(LNUM)—U1((LNUN-1)))/DUX2=(U2(LNUM)-U2((LNUM-1)))/DCC DRIFT VELOCITYC420 V1=VDRI1(X(J))V2=VDRI2(X(J))AD 1=0 1AD2=D2E=FIELD(X(J))IF(E.EQ.0.ODOO) GO TO 350AD1-V1*CKTQ/EAD2=-V2*CKTQ/E350 TEMP=(—AD1*UX1+V1*U1(J))+(—A02*UX2+V2*U2(J))AN=AN+TEMP*FIF(J.NE.1) GO TO 360TP1(I)=TEMP360 IF(3.NE.LNUN) GO TO 400TP2(I)=TEMP400 CONTINUEMTENP=LNUM+M2+1DO 450 J=MTENP,101READ (4 460)450 CONTINUE460 FORMAT()70 FORNAT(14,D10.3,2(1X,D11.4))TP3(I)=AN20 CONTINUEDO 500 I=1,NSETWRITE(7,30)TT(I) ,TP1(I) ,TP2(I) ,TP3(I)500 CONTINUE30 FORMAT(4(1X,D11.4))RETURNEND105B.3 List of EXTRAC PROGRAM EXTRACWRITE(6, 10)10 FORMAT(’Fl ‘)READ (S , 20) FlWRITE(6,30)30 FORMAT(/,’F2= ‘)READ (S , 20 ) F2WRITE(6,40)40 FORMAT(/,’F3 ‘)READ(5,20)F3WRITE(6 ,S0)50 FORMAT(/,’T ‘)RE.4D(S,20)TWRITE(6,60)60 FORMAT(/’LINE #= ‘)READ (S , 70 )N20 FORMAT(Fl1.4)70 FORMAT(13)T=T/30.0P=o.S*(Fl—F3)/(Fl—2*F2+F3)IF(ABS(P).GT.l.0) GOTO 100p=(p+N*l.0/2.0)*TWRITE(6,80)P80 FORMAT(’TMAX ‘,Ell.4)STOP100 WRITE(6,l10)110 FORMAT(’ABS(P)>l.O’)STOPEND106Appendix C Listing of Programs Used to Estimate Signal Delay Times in HBTsOperating at High CurrentsC.1 List of SUBFINE7C FINE INTEGRATE SYSTEM ONE STEP FORWARD.C THE SYSTEM IS SPECIFIED BYCC DF/DT + DG/DX + II 0CC EVEN TIME STEP IS ASSUMED.CC MAXIMUM NUMBER OF POINTS, NPT*N, IS 512C MAXIMUM NUMBER OF EQUATIONS, N, IS 4C MAXIMUM NUMBER OF GRIDS IS 150, INCLUDING START AND END POINTSCC WHEN ENTER, N NUMBER OF EQUATIONSC NPT NUMBER OF GRIDSC TO TIMEC T TO+DELTC YO SULOTION AT TO, MUST BE YO(512)C Yl GUESS SULOTION AT T, MUST BE Y1(512)C IFLAG INDICATE IF THIS IS THE FIRST TIMEC THE SUB. IS CALLED. IF 0, THE FIRST CALL.C MUST BE A VARIABLE.C ITYPE TYPE OF SYSTEM. IF 0, LINEARC WHEN EXIT Yl SULOTION AT TO+DELTC ITYPE SET TO 1CSUBROUTINE FINE(N,NPT,TO,T,YO,Y1,IFLAG,ITYPE)IMPLICIT REAL*8 (A—H,O-Z)DIMENSION YO(512),Y1(512)CC ARRAYS FOR SLE, A SUBROUTINE TO SOLVE A*X=BCDIMENSION DT(512,512) ,IPERM(1024) -CC COMMON AREA MUST BE SPECIFIED IN MAIN PROGRAM ASC AS FT(9*15O+512*512+4)C THE HIGHEST DIMENSION FOR THE SYSTEM, THEN, IS 148CCOMMON/COMM/GRID(150),ACOE(150),BCOE(150) ,CCOE( 150)/ ,TDX(150),RC(3,lso),TH(15o),DDT(4),DJACO(512,512)CC INTEGERS RESERVED FOR FINE, MUST BE SPECIFIED IN MAIN PROGRAM ASC IFT(4)CCOMMON/ICOM/NLAST, NLAST3 , NLAST4 , NLASTSCC COMMON AREA FOR ALL SUBROUTINES WITHIN FINECCOMMON/FINE1/Y2(512) ,Y3(512) ,DELY(512) ,F(512)COMMON/FINE2/DELT,T1 ,T2C START OF FINE107DELT=1 ODOO/(T—T0)C CHECK IF THIS IS THE FIRST CALLIF(IFLAG.NE.0) GO TO 20C SET THE FLAGIFLAG=1CALL GRIDS((NPT+2),GRID)CALL FINIT(N,NPT)CALL DEQDDT(N,DDT)DO 40 I=1,NLAST3Y2(I)=Y0(I)40 Y3(I)=Y1(I)T2=T0T 1=TC CALCULATE JACOBIANCALL FUNJ(N,NPT)GO TO 3020 CONTINUEDO 50 I=1,NLAST3Y2(I)=Y0(I)50 Y3(I)=Y1(I)T2=T0T1=TIF(ITYPE.EQ.0) GO TO 30C CALCULATE JACOBIANCALL FUNJ(N,NPT)30 CONTINUEC CALCULATE BCALL FUNB(N,NPT,F)C SOLVE A*X=BCALL SLE(NLAST3,512,DJACO, 1,512,F,DELY,IPERM,512,DT,DDET,JEXP)DO 10 I=1,NLAST310 Y1(I)=Y1(I)+DELY(I)RETURNENDCC FINIT CALCULATES VARIOUS CONSTANTSCC WHEN ENTER, N THE NUMBER OF EQUATIONSC NPT THE NUMBER OF GRIDSC WHEN EXIT, CONSTANTS SUCH AS ACOE, BCOE WILL BE SETC AND RETURNED THROUGH COMMON AREACSUBROUTINE FINIT (N, NPT)IMPLICIT REAL*8 (A-H,O—Z)COMMON/COMM/GRID(150) ,ACOE(150) ,BCOE(150) ,CCOE(150)/ ,TDX(150) ,RC(3, 150) ,TH(150) ,DDT(4) ,DJACO(512,512)COMMON/ICOM/NLAST, NLAST3 , NLAST4 , NLASTSNLAST=NPTNLAST3=NPT*NNLAST4=NPT+ 1NLAST5=NPT+200 20 I=1,NLAST420 TDX(I)=1 .0000/(GRID((I+1))—GRID(I))108DC 10 11,NPTRAT=(GRID((I+1))—GRID(I))/(GRID((I+2))—GRID((I+1)))ACOE(I)=(1 ODOO+RAT*(1 ODOO—2 . ODOO*RAT)/ /(1.0000+RAT))/12.oDOoCCOE(I)=(1.0D00+(RAT—2.ODOO)/(RAT*(1.ODOO+RAT)))/12.D00BCOE(I)=1 ODOO—ACOE(I)—CCOE(I)TH(I)=(GRID((I+2))—GRID(I))*0.SDOORC(1 , I)=(GRID((I+1))—GRID(I))*TDX((I+1))*0.5D00/TH(I)RC(2,I)=(GRID(I)+GRID((I+2))—2.ODOO*GRID((I+1)))*TDX(I)/ *TDX((I+1))RC(3,I)=(GRID((I+1))—GRID((I+2)))*0.SDOO*TDX(I)/TH(I)10 CONTINUEDO 510 L=1NPTDO 520 K=1,NPTDJACO(L,K)=0.ODOO520 CONTINUE510 CONTINUERETURNENDCC FUNB RETURNS MATRIX BCC WHEN ENTER, N THE NUMBER OF EQUATIONSC NPT THE NUMBER OF GRIDSC WHEN EXIT, B MATRIX BCSUBROUTINE FUNB(N ,NPT,B)IMPLICIT REAL*8 (A—H,O-z)DIMENSION B(512)COMMON/COMM/GRID(150) ,ACOE(150) ,BCOE(150) ,CCOE(150)/ ,TDX(150),RC(3,150),TB(150),DDT(4),D3ACO(512,512)COMMON/ICOM/NLAST , NLAST3 , NLAST4 , NLAST5COMMON /FINE1/Y0(512),Y1(512),DELY(512),F(512)COMMON /FINE2/DELT,T,TODIMENSION G1(4),G2(4)DIMENSION H1(4) ,112(4) ,113(4)DIMENSION UT1(4) ,UT2(4) ,UT3(4)DIMENSION HU1(4) ,HU2(4)DIMENSION HUX1(4) ,HUX2(4)DIMENSION UX1(4) ,UX2(4) ,UX3(4)DIMENSION U10(4),U20(4),U30(4),U40(4),U50(4)DIMENSION U1(4),U2(4),U3(4),U4(4),U5(4)DO 10 I=1.,NPT11=1—212=1—113=1+114=1+2X1GRID(I)X2=GRID((I+1))X3=GRID((I+2))X4=0 5D00*(X1+X2)X50 . 5D00*(X2+X3)T2=0.5D00*(T+T0)109CALL UVALUE(I1,N,U10,U1)CALL UVALUE(12,N,U20,U2)CALL UVALUE(I,N,U30,U3)CALL UVALUE(13,N,U40,U4)CALL UVALUE(14,N,TJSO,U5)CALL DUDX(12,N,U1,U2,U3,UX1)CALL DUDX(I,N,U2,U3,U4,UX2)CALL DUDX(13,N,U3,U4,U5,UX3)CALL HUVL(N,U2,U3,HU1)CALL HUVL(N,U3,U4,HU2)CALL RDUDX(I,N,U2,U3,JIUX1)CALL HDUDX(13,N,U3,U4,HUX2)CALL FLUX(N,X4,T2,HU1,HUX1,G1)CALL FLUX(N,X5,T2,HU2,HUX2,G2)CALL FUNH(N,X1,T2,U2,UX1,H1)CALL FUNR(N,X2,T2,U3,UX2,H2)CALL FUNH(N,X3,T2,U4,UX3,113)DO 20 J=1,NIT=I+(J—1)*NPTB(IT)=—TH(I)*DELT*DDT(J)*(ACOE(I)*U40(J)/ +BCOE(I)*U30(J)+CCOE(I)*U2O(J))/ .-(G2(J)—G1(J))—TH(I)*(ACOE(I)*H3(J)+BCOE(I)*H2(J)/ +CCOE(I)*H1(J))20 CONTINUE10 CONTINUERETURNENDCC FULl RETURNS JACOBIANCC WREN ENTER, N NUMBER OF EQUATIONSC NPT NUMBER OF GRIDSCSUBROUTINE FUNJ(N,NPT)IMPLICIT REAL*8 (A—H,O—z)COMMON/COMM/GRID(150) ,ACOE(15O) ,BCQE(1S0) ,CCOE(1bO)/ ,TDX(150),RC(3,1S0),TH(1S0),DDT(4),DJACO(512,512)COMMON/ICOM/NLAST , NLAST3 , NLAST4 , NLASTSCOMMON /FINE1/Y0(512) ,Y1(512) ,DELY(512) ,F(512)COMMON /FINE2/DELT, T,TODIMENSION G1(2,4,4),G2(2,4,4)DIMENSION H1(2,4,4),R2(2,4,4),113(2,4,4)DIMENSION RU1(4) ,HU2(4)DIMENSION UX1(4),UX2(4),UX3(4),HUX1(4),HUX2(4)DIMENSION U10(4) ,U20(4) ,U30(4) ,U40(4) ,U50(4)DIMENSION U1(4) ,U2(4) ,U3(4) ,U4(4) ,U5(4)DO 10 J=1,NPTDO 20 I=1,NPT11=1—2121—113=1+114=1+2DELTA1=DEL(I1 ,J)110DELTA2=DEL(12, 3)DELTA3=DEL(I ,3)DELTA4=DEL(13, 3)DELTAS=DEL(14, 3)IF((DELTA1+DELTA2+DELTA3+DELTA4+DELTA5) .EQ. 0. ODOO)/ GO TO 500X1=GRID(I)X2=GRID((I+1))X3=GRID((I+2))X4=O . 5000*(X2+X1)X5=0 . 5D00*(X3+X2)T20 . 5D00*(T+TO)CALL UVALUE(I1,N,U1O,U1)CALL UVALUE(12,N,U20,U2)CALL UVALUE(I,N,U30,U3)CALL UVALUE(I3,N,U40,U4)CALL UVALUE(14,N,U50,U5)CALL DUDX(12,N,U1,U2,U3,UX1)CALL DUDX(I,N,U2,U3,U4,UX2)CALL DUDX(13,N,U3,U4,U5,UX3)CALL HUVL(N,U2,U3,HUI)CALL HUVL(N,U3,U4,HU2)CALL HDUDX(I,N,U2,U3,HUX1)CALL HDUDX(13,N,U3,U4,HUX2)CALL DGDU(N,X4,T2,HU1,HUX1,G1)CALL DGDU(N,X5,T2,HU2,HUX2,G2)CALL DHDU(N,X1,T2,U2,UX1,H1)CALL DHDU(N,X2,T2,U3,UX2,H2)CALL DHDU(N,X3,T2,U4,UX3,113)DO 30 L=1,NDO 40 K1,NDELKL=DEL (K, L)IT=I+(K—1)*NPT.JT=J+(L—1)*NPTDJACO(IT, JT)=TH(I)*DELT*DDT(K)*DELKL*(ACOE(I)*DELTA4/ +BCOE(I)*DELTA3/ +CCOE(I)*DELTA2)+0.5D00*((G2(1,K,L)*(DELTA4+DELTA3)*0.5D00/ +G2(2,K,L)*TDX(13)*(DELTA4—DELTA3))—(G1(1,K,L)*0.5D00/ *(DELTA3+DELTA2)+G1(2,K,L)*TDx(I)*(DELTA3—DELTA2)))/ +0.5D00*TH(I)*(ACOE(I)*(H3(1 ,K,L)*DELTA4+H3(2,K,L)*DUXDU(13,3))/ +BCOE(I)*(R2(1,K,L)*DELTA3+H2(2,K,L)*DUXDU(I,J))/ +CCOE(I)*(H1(1 ,K,L)*DELTA2+H1(2,K,L)*DUXDU(I2,J)))40 CONTINUE30 CONTINUE500 CONTINUE20 CONTINUE10 CONTINUERETURNENDCC UVALUE RETURNS U AT (T+T0)/2CC WHEN ENTER, I I-TH POINT111C N NUMBER OF EQTJATIONSC UO DIFFERENCE OF U(T)-U(TO)C WHEN EXIT, U U AT T+DELT/2CSUBROUTINE UVALUE(I ,N,UO,U)IMPLICIT REAL*8 (A—H,O—Z)DIMENSION UO(4) ,U(4)COMMON/COMM/GRID(150) ,ACOE(150) ,BCOE(150) ,CCOE(150)/ ,TDX(150),RC(3,150),TH(150),DDT(4),DJAC0(512,512)COMMON/ICOM/NLAST , NLAST3 , NLAST4 , NLASTSCOMMON /FINE1/YO(512) ,Y1(512) ,DELY(512) ,F(512)COMMON /FINE2/DELT,T,T0DIMENSION BV1(4) ,BV2(4)IF(I.LT.1) GO TO 20IF(I.GT.NLAST) GO TO 30DO 10 .31,NIO=I+(J—1)*NLASTU(J)=O .5D00*(Y1(IO)+YO(IO))UO(J)=Y1(IO)—YO(IO)10 CONTINUERETURN20 CONTINUECALL BNDRY(N,0,T,BV1)CALL BNDRY(N,O,T0,BV2)DO 40 .J=i,NU(J)=O. 5D00*(BV2(J)+BV1(J))UO(J)=BV1(J)—BV2(J)40 CONTINUERETURN30 CONTINUECALL BNDRY(N,1,TO,BV1)CALL BNDRY(N,1,T,BV2)DO 50 3=1,NU(J)=0.5Doo*(Bvl(J)+Bv2(J))U0(J)=BV2(J)—BV1(J)50 CONTINUERETURNENDCC DUDX RETURNS DU/DX AT (T+T0)/2CC WHEN ENTER, I I-TB POINTC N NUMBER OF EQUATIONSC U1—3 INPUT OF UC WHEN EXIT, UX VALUE OF DU/DXCSUBROUTINE DUDX(I,N,U1,U2,U3,UX)IMPLICIT REAL*8 (A—H,O—z)DIMENSION U1(4) ,U2(4) ,U3(4) ,UX(4)COMMON/COMM/GRID(150) ,ACOE(150) ,BCOE(150) ,CCOE(150)/ ,TDX(150),RC(3,150),TH(150),DDT(4),DJACO(512,512)COMMON/ICOM/NLAST , NLAST3 , NLAST4 , NLAST5COMMON /FINE1/YO(512) ,Y1(512) ,DELY(512) ,F(512)112COMMON /FINE2/DELT , T ,TOIF(I.EQ.0) GO TO 20IF(I.EQ.(NLAST+1)) GO TO 30DO 10 J=1,NUX(3)=RC(1,I)*U3(J)+RC(2,I)*U2(3)+RC(3,I)*U1(J)10 CONTINUERETURN20 DO 40 i=1,N40 UX(J)=(U3(J)—U2(3))*TDX(1)RETURN30 DO 50 J=1,Nso UX(J)=(U2(J)—U1(J))*TDX(NLAST4)RETURNENDCC HUVL RETURNS U VALUE AT (T+T0)/2, X(I+1/2)CC WHEN ENTER, N NUMBER OF EQUATIONSC U1—2 INPUT OF UC WHEN EXIT, HU VALUE OF U AT X(I+1/2)CSUBROUTINE HUVL(N,U1,U2,HU)IMPLICIT REAL*8 (A—H,O-z)DIMENSION U1(4) ,U2(4) ,BU(4)COMMON/COMM/GRID(150) ,ACOE(150) ,BCOE(150) ,cCoE(1S0)/ ,TDX(150),RC(3,150),TH(150),DDT(4),D.JACO(512,512)COMMON/ICOM/NLAST , NLAST3 , NLAST4 , NLASTSCOMMON /FINE1/Y0(512) ,Y1(512) ,DELY(512) ,F(512)COMMON /FINE2/DELT , T ,TODO 10 J=1,N10 HU(J)=O.SDOO*(U1(J)+U2(J))RETURNENDC HDUDX RETURNS DU/DX AT (T+T0)/2, X(I+1/2)CC WHEN ENTER, I I— AND (I—1)—TH POINTSC N NUMBER OF EQUATIONSC U1—2 INPUT OF UC WHEN EXIT, HUX VALUE OF DU/DX AT X(I+1/2)CSUBROUTINE HDUDX(I,N,U1,U2,HUX)IMPLICIT REAL*8 (A-H,O-z)DIMENSION U1(4) ,U2(4) ,RUX(4)COMMON/COMM/GRID(15O),ACOE(150),BCOE(15O),CCOE(150)/ ,TDX(1SO),RC(3,150),TH(1SO),DDT(4),DJACO(512,512)COMMON/ICOM/NLAST, NLAST3 , NLAST4 , NLAST5COMMON /FINE1/YO(512) ,Y1(512) ,DELY(512) ,F(512)COMMON /FINE2/DELT , T ,TODO 10 J=1,N10 HUX(J)=TDX(I)*(U2(J)—U1(J))RETURNEND113CC DUXDU RETURNS D(UX)/DUCC ON ENTRY I,J INDEXC EXIT VALUE OF D(UX)/DUFUNCTION DUXDU(I,J)IMPLICIT REAL*8 (A—H,O—Z)DIMENSION U1(4) ,U2(4) ,HU(4)COMMQN/COMM/GRID(160),ACOE(150),BCOE(150),CCGE(1S0)/ ,TDX(1S0),RC(3,150),TH(1S0),DDT(4),DJACO(512,612)COMMON/ICOM/NLAST , NLAST3 , NLAST4 , NLASTSCOMMON /FINE1/Y0(512) ,Y1(512) ,DELY(512) ,F(512)COMMON /FINE2/DELT,T,T0IF(I.EQ.O) GO TO 20-IF(I.EQ.(NLAST+1)) GO TO 30C DO 10 J=1,NDUXDU=RC(1,I)*DEL((I+1),J)+RC(2,I)*DEL(I,J)+RC(3,I)*DEL((I—1),J)C CONTINUERETURNC DO 40 J=1,N20 DUXDU=(DEL((I+1),J)—DEL(I,J))*TDX(1)RETURNC DO 50 J=1,N30 DUXDU=(DEL(I,J)—DEL((I—1),J))*TDX(NLAST4)RETURNENDCC DEL RETURNS 1 IF I=J, OTHERWISE 0.CFUNCTION DEL(I,J)IMPLICIT REAL*8 (A—H,o—z)DEL=0 . ODOOIF(I.EQ.J) DEL=1.ODOORETURNEND114C.2 List of EQUI2IMPLICIT P.EAL*8 (A—E,O-z)COMMON/COMM/FT (263498 )COMMON/ICOM/IFT (4)CDIMENSION Yo(512) ,Y1(512) ,GPID(1S0) E(512)COMMDN/IEL/NPTCOMMDN/PARA/VELO(3) ,AMU(3) ,DIFF(3),/ E0,TAU,C1,C2,C3,QEPS,DOP1,DOP2,DOP3POT,BVL,RATIOCOMI1ON/COOR/XCON 1, XCON2 , XCON3XCON4COMMON/BCOND/RIGBT1, RIGRT2DATA IIJMST/20/NPT= 110N=4NTOTN*NPTIFLAG=0EPS1 . OD—2READ(5,2)NSTEP,TAU,DOP1,DOP2,DOP3,PEIWRITE(6,1)NSTEP ,TAU,DOP1 ,DOP2,DDP3,PHI1 FORI4AT(’ NSTEP= ‘,I5,/ ‘ TAU2= ‘,D11.4,/,/ ‘DOPINGl= ‘, D11.4, ‘ DOPING2= ‘, D11.4, ‘ DOPING3=’,D11.4,/ ‘ BIAS (V)= ‘, D11.4)2 FOR14AT(15,5D11.4)READ(5.3)XCON1 ,XCON2,XCDN3.XCON4WRITE(6,3)XCON1 ,XCON2,XCON3,XCON43 FORMAT(4D11.4)READ(5,4) BVL,RATIOWRITE(6,4) BVL,P.ATIO4 FORMAT(2D11.4)CC INITIALIZE PROGRAMCT0.ODOOCC RECIPRECAL OF E0CEO=1.ODOO/3.4D—2VELO(1)=10.ODOOAMU(i)=60 . ODOODIFF(1)=1 .554D00VELO(2)=1 .0000AMU(2)=29 . ODOODIFF(2)=0 . 751D00AMU(3) 4. 0000DIFF(3)=1 . 0340—01C1=4.ODOO/27.0000C2=2*C 1C3=4.0000/9.0000POT=0. 0259D00*DLOG((DABS(DOPI*DOP3)/5.0625D—18))+PEIQEPS1 . 38132087D—02CC REDEFINE ThU TO BE 1. ODOO/TAU115CTAU=1 . ODOO/TAUCALL GRIDS((NPT+2) ,GRID)CC INITIAL CONDITIONCTEMP2=D0P3/(1 . ODOO+RATIO)TEMP1=DOP3—TEMP2TEMP3=BVL/ (GRID (29) —GRID (1))C RIGHT1=TEMP1C RIGHT2=TEMP2DO 10 1=1,29Y0(I)=BVL—TEMP3*(GRID(I)—GRID(1))Y0((I+NPT))=0. ODOOYo((I+2*NPT))=—DOPIY0((I+3*NPT))=O.ODOOYO((NPT-I+1))=TEMP1YO((2*NPT—I+1))=TEMP2Yo((3*NPT—I+1))=o.ODO0Y0 ( (4*NPT—I+1) ) =POT10 CONTINUEDP=POT/53 . ODOODO 30 1=1,52Y0((I+29))=0.ODOOY0((I+29+NPT))0 . ODOOY0 C (I+29+2*NPT) ) =0. ODOOY0((I+29+3*NPT))=DP*DFLOAT(I)30 CONTINUEDO 20 11,NTOT20 Yi(I)=Y0(I)DELTO.SDOOTOTT=T+DELTDDOP=DOP2/NUMSTDO 40 11,NUMSTDOP2=DDOP*DFLOAT (I)250 K1100 CONTINUEC RIGET1=Y1(NPT)C RIGET2=Y1 C (2*NPT))C DO 270 L=1,2CALL SPOTEN(N,NPT,T0,T,Y0,Y1,IFLAG,1)270 CALL SPART(N,NPT,T0,T,Y0,Y1,IFLAG,1)CALL FINE(N,NPT,T0,T,Y0,Y1,IFLAG,1)CC UPDATE YO AND YlCDO 90 J=1,NTOT90 Yo(J)=0.sDoo*(n(J)+Y0(3))DO 95 3=1,NTOP95 Yi(J)=Yo(J)CIF(K.NE.1) GO TO 200116K=2DO 210 L=1,NTOT210 E(L)=Y1(L)GO TO 100200 ERROF=DERR(Y1 ,E,NTOT)WRITE(6,201)ERROR201 FORMAT(D11.4)DO 230 3=1,NTOT230 E(J)=Y1(J)K=K+1IF(K.GT.60) GD TO 300IF(ERROR.GT.EPS) GO TO 100IF(I.EQ.(NUMST—1)) EPS=1.OD—0840 CONTINUECALL PRINT(N,NPT,GRID,Y1)STOP300 WRITE(6,310)310 FORMAT(’TOO MANY ITERATIONS’)STOPENDSUBROUTINE GRIDS (NPT2, GD)IMPLICIT REAL*8 (A—H,O—z)DIMENSION GD(150)COMMON/CO OR/XCON 1, XCON2 , XCON3 , XCON4TEMP 1=(XCON2—XCON1—0. 1SDOO)/26. ODOOTEMP2= (XCON3—xCON2—0 . 3D00) /41 ODOOTENP3=(XCON4—XCON3—0. 1SDOO)/26 ODOODO 10 1=1,2510 GD(I)=XCON1+(I—1)*TEMP1DO 35 1=26,3535 GD(I)=GD(25)+(I—25)*0.03D00DO 20 1=36,7620 GD(I)=GD(35)+(I—35)*TEMP2DO 45 1=77,8645 GD(I)GD(76)+TEMP2+CI—76)*0.03D00DO 30 1=87,11230 GD(I)GD(86)+(I—86)*TEMP3RETURNENDSUBROUTINE DEQDDT (N, DOT)IMPLICIT REAL*8 (A-E,O-Z)DIMENSION DOT(4)DDT(1)1.ODOODDT(2)1 .0000DDT(3)=1 . ODDODDT(4)0. 0000RETURNENDSUBROUTINE FLUX(N,X,T,U,UX,G)IMPLICIT REAL*8 (A—B,O—Z)DIMENSION 13(4) ,Ux(4) ,G(4)COMMON/PARA/VELO(3) ,AMIJ(3) ,DIFF(3),/ E0,TAU,C1,C2,C3,QEPS,DQPI,DOP2,D0P3,POT,BVL,R.ATIO117DIMENSION v(3) ,B(3)CALL VMU(X,U,UX,V,B)G(1)=—B(1)*UX(1)+V(1)*U(1)G (2) =—B (2) *Ux(2) (2) *U(2)G(3)=—B(3)*UX(3)+V(3)*U(3)G (4) =UX (4)RETURNENDSUBROUTINE FUJR(NX,T,U,UX,R)IMPLICIT REAL*8 (A—H,O—Z)DIMENSION U(4) ,T3X(4) ,E(4)COMMON/PARA/VELO(3) ,AMU(3) bDIFF(3),/ E0,TAU,C1,C2,C3,QEPS,DOP1,00P2,DOP3,POT,BVL,RATIODIMENSION v(3) ,B(3)CALL P.ATE(X,U,UX,R)N(1)=TAU*(R*U(1)—U(2))(2) =—E (1)N(3)0. ODOOU(4)=—QEPS*(U(1)+U(2)—U(3)—DOP(x))RETURNENDSUBROUTINE DGDU(N,X,T,U,UX,DG)IMPLICIT REAL*8 (A—B,O--z)DIMENSION U(4),UX(4),DG(2,4,4)COMMON/PARA/VELO(S) ,AKU(3) ,DIFF(3),/ E0,TAU,C1,C2,C3,QEPS,DOP1,DOP2,DOP3,POTBVL,RATIODIMENSION V(3) ,B(3) ,DV(3) ,DB(3)CALL VMU(X,U,UX,V,B)CALL DVMU(X,U,UXDV,DB)DO 10 K=1,2DO 10 1=1,4DO 10 J1,410 DG(K,I,J)0.ODOODG (1, 1, 1) =V( 1)DG (1, 2, 2) =V (2)DG(1,3,3)=V(3)DG (2 1, 1) =—B (1)DG(2,2,2)=-B(2)DG(2,3,3)=—B(3)DG(2,4,4)1 .ODOODG(2,1,4)DB(1)*UX(1)—DV(1)*U(1)DG(2,2,4)=DB(2)*UX(2)—DV(2)*U(2)DG(2,3,4)DE(3)*UX(3)—DV(3)*U(3)RETURNENDSUBROUTINE DRDU(N,X,T,U,UX,DR)IMPLICIT REAL*8 (A—R,0—Z)DIMENSION U(4),UX(4),DE(2,4,4)COMMON/PARA/VELO(3) ,ANU(3) ,DIFF(3),/ E0,TAU,C1 ,C2,C3,QEPSJDOP1,DOP2,DOP3,POT,BVL,RATIODIMENSION V(3),B(3),DV(3),DB(s)CALL RATE(X,U,UX,R)CALL DRATE(X,U,UX,DR)118DO 10 K=1,2DO 10 1=1,4DO 10 J=1,410 DE(K,I,J)=0.ODOODR(14 1)=—QEPSDE(1,4,2)=—QEPSDH(1,4,3)QEPSDR (1,1 , 1) =TAU*RDR(1, 1 ,2)=-TAUDE(1,2,1)=—DE(1,1,1)DR(1,2, 2) =—DH(1, 1, 2)DE(2, 1,4)—TAU*DR*U(1)DE(2,2,4)=—DE(2,1,4)RETURNENDSUBROUTINE BNDRY(N, IFLAG,T,BV)IMPLICIT REAL*8 (A—H,O-z)COMMON/IEL/NPTCOMMON/PARA/VELO(3) AMU(3) ,DIFF(3),/ E0,TAU,C1,.C2,C3,QEPS,DOP1,DOP2,DOP3,POT,BVL,PATIOCOMMON/BCDND/RIGHTI • RIGRT2DIMENSION BV(4)IF(IFLAG.EQ.1) GO TO 10BV(1)=BVLBV(2)=0.ODOOBV(3)—DOP1BV(4)0 .ODOORETURN10 BV(2)D0P3/(1 .0000+RATIO)BV(1)=DOP3-BV(2)C 10 BV(1)=RIGHT1C BV(2)=RIGET2BV(3)0.ODOOBV(4)=POTRETURNENDSUBROUTINE VMU(X,U,UX,V,B)IMPLICIT REAL*8 (A—E,O-Z)DIMENSION U(4),UX(4),v(3),B(3)COMMON/PARA/VELO(3) ,AMU(3) ,DIFF(3),/ E0,TAU,C1,C2,C3.QEPS,DOP1 ,DOP2,DOP3,POT,BVL,RATIODATA CKTQ/0. 0259D00/FIELD—UX (4)DO 10 I—1b2BX=AMU(I)*FIELDIF(DABS(BX).GT.(1.SDOO*VELO(I))) GO TO 20TEMP=BX/VELO (I)V(I)=—BX*(1 ODOO—TEMP*TEMP*C1)GO TO 3020 V(I)VELO(I)IF(FIELD.GT.0.ODOO) v(I)=—vELO(I)30 B(I)=DIFF(I)IF(FIELD.EQ.0.ODOO) GO TO 10119B ( I) =—v (I) *CKTQ/FIELD10 CONTINUEv(3)=FIELD*AMu(3)B (3) =DIFF (3)RETURNENDSUBROUTINE RATE(X,U,UX,R)IMPLICIT REAL*8 (A—B,o-z)DIMENSION U(4) ,Ux(4)COMMON/PARA/VELO(3) ,AMU(3) ,DIFF(3),/ E0,TAU,C1 ,C2,C3,QEPS,DOP1 ,DOP2,DOP3,POT,BVL,RATIODATA CKTQ,RCKTQ/0 .02&9D00, 38. 61003861D00/DATA TMAX/31 . 6227766D00/FIELD=-UX(4)TEMP=FIELD*E0IF(DABS(TEMP).GE.TMAX) GO TO 10R=TENP*TEMP*TEMP*TEMPRETURN10 R=1.0D6RETURNENDSUBROUTINE DVMU(X,U,UX,DV,DB)IMPLICIT REAL*8 (A—B,O—z)DIMENSION U(4) ,UX(4) ,Dv(3) ,DB(3)COMMON/PARA/VELO(3) ,AMU(3) ,DIFF(3),1 EO,TAU,C1,C2,C3,QEPS,DOP1,DOP2,DOP3,POT,BVL,RATIODATA CKTQ/0 . 0259D00/FIELD=-UX (4)DO 10 1=1,2BX=FIELD*AMU(I)IF(DABS(BX) . GT. (1. SDOO*VEL0(I))) GO TO 20TEMP=BX/VELO (I)DV(I)=—AMU(I)*(1 .ODOO—TEMP*TEMP*C3)DB(I)=—DIFF(I)*AMU(I)*AMU(I)*C2*FIELD/(VELO(I)*VELO(I))GO TO 1020 DV(I)=0.ODOODB(I)=—VELO(I)*FIELD*CKTQ/(FIELD*FIELD*DABS(FIELD))10 CONTINUEDV(3)=AMU(3)DB(3)=0.ODOORETURNENDSUBROUTINE DRATE(X,U,UX,DR)IMPLICIT REAL*8 (A—N,0—z)DIMENSION U(4),UX(4)COMMON/PARA/VELO(3) ,AMU(3) ,DIFF(3),/ E0,TAU,C1,C2,C3,QEPS,DOP1 ,DOP2,DOP3,POT,BVL,RATIODATA CKTQ ,RCKTQ/0 .0259D00, 38. 61003861D00/DATA TMAX/31 . 6227766D00/FIELD=-UX(4)TEKP=FIELD*E0IF(DABS(TEMP).GE.TMAX) GO TO 10DR=4 . OODOO*TEMP*TEMP*TEMP*E0120RETURN10 DR=O.ODOORETURNENDFUNCTION DOP(X)IMPLICIT REAL*8 (A—R,O—Z)COMMON/PARA/VELO(3) ,AMU(3) ,DIFF(3),/ E0,TAU,C1,C2,C3,QEPSDOP1,DOP2,DOP3,POT,BVL,R.ATIOCOMMON/COOR/XCON 1, XCON2,ICON3 , XCON4IF(X.LE.XCDN2) GO TO 10IF(X.LE.XCON3) GO TO 20DOP=DOP3IF(X.LT. (XCON3+0. 1SDOO)) DOP=0.SDOO*(DOP3+DOP2)+/ 0. 5D00* (DQP3—DOP2) * (x—XCON3) *6. 66666666TDOORETURN10 DOP=DOP1IF(X. GT. (XCON2—0. 1SDOO)) DOPO .5D00*(DOP1+D0P2)+/ 0. 5D00*(DOP1—DOP2)*(XCON2—X)*6 .666666667D00RETURN20 DOP=DOP2IF(X.LT. (xC0N2+0. 1SDOO)) DOP=0. SDOO*(DOP2+DOP1)+/ 0 .SDOO*(DOP2—DOP1)*(X—XCON2)*6 .666666667D00IF(X. GT. (xC0N3—0. 1SDOO)) DOPO .5D00*(DOP2+DOP3)+/ 0 .5D00*(DOP2—DOP3)*(XCON3—X)*6.666666667D00RETURNENDFUNCTION DERR(Y,EN)IMPLICIT REAL*8 (A—H,O—z)DIMENSION Y(512),E(512)DEL=0.ODOODMOD=0.ODOODO 10 I=1,NDMOD=DMOD+Y (I) *y(I)TEMP=Y(I)—E(I)10 DEL=DEL+TEMP*TEMPDERR=DSQRT (DEL/DMOD)RETURNENDSUBROUTINE SPOTEN(N,NPT,T0,T,Y0,Y1,IFLAG,ITYPE)IMPLICIT REAL*8 (A—R,O—Z)COMMON/COMM/FT (263498)COMMON/ICOM/IFT(4)CDIMENSION Y0(512) ,Y1(512)C DO 20 K1,4CALL FINE(NNPT,TO,T,Y0,Y1,IFLAG,ITYPE)CC UPDATE POTENTIAL ONLYCDO 10 I=1,NPTY1(I)=Y0(I)• Y1((I+NPT))=Y0((I+NPT))Y1((I+2*NPT))=Y0((I+2*NPT))121Y0((I+3*NPT))=0.SDOO*(yl((I+3*NpT))+yo((x+3*NpT)))10 Y1((I+3*NPT))=Y0((I+3*NPT))RETURNENDSUBROUTINE SPART(NNPT,T0TJ0,Y1,IFLAG,ITYPE)IMPLICIT REAL*8 (A-R,o-z)COI4MON/COMI1/T1’ (263498)CONMON/ICOM/IFT (4)CDIMENSION Y0(512) ,Y1(512)C DC 20 K=1,4CALL FINE(N,NPT,T0,T,Y0,Y1,IFLAGITYPE)CC UPDATE DENSITIES ONLYCDO 10 I=1,NPTY0(I)=0.SDOO*(Y0(I)+Y1(I))Y1(I)=YO(I)Yo((I+NPT))=0.SD0O*(Yo((I+NPT))+Y1((I+NPT)))Y1((I+NPT))=Y0((I+NPT))Y0((I+2*NPT))=0.5Doo*(yo((I+2*NpT))+yl((I+2*NpT)))Y1((I+2*NPT))=Y0((I+2*NPT))10 Y1((I+3*NPT) )=Y0((I+3*NPT))RETURNENDSUBROUTINE PRINT(N,NPT,GRID ,Y1)IMPLICIT REAL*8 (A—fi,O-z)DIMENSION GRID(1S0) ,Y1(512)DIMENSION U(4) ,Ux(4) ,G(4)DO 10 J1NPTWRITE(6,20)J,GRID((i+1)),Y1(J)Y1((J+NPT)),Y1((J+HPT*2)),/ Y1((3+NPT*3))10 CONTINUE20 FORMAT(IS,5(1X,D11.4))T=0.SDOON1NPT—1DO 30 I=1N1X=0.SDOO*(GRID((I+2))+GRID((I+i)))U(1)=0.SDO0*(Y1((I+1) )+fl(I))U(2)=0.SDO0*(Y1((I+NPT+1))+Y1((I+NPT)))U(3)=0 .5D00*(Y1((I+NPT*2+1))+Y1((I+2*NPT)))U(4)=O.5D00*(Y1((I+NPT*3+1))+Y1((I+3*NPT)))DELX=1 .ODOO/(GRID((I+2))—GRID((I+1)))UX(1)=(Y1((I+1) )—Y1(I))*DELXUX(2)=(Y1((I+NPT+1))—Y1((I+NPT)))*DELXUX(3)=(Y1((I+2*NPT+1))—Y1((I+2*NPT)))*DELXUx(4)=(Y1((I+3*NPT+1))—Y1((x+3*N?T)))*DELxCALL FLUX(3,I,T,U,UX,G)CRNT=G(3)—G(1)—G(2)FIELD=-UX(4)WRITE(6,40)I,X,CILNT,FIELD30 CONTINUE40 FORMAT(IS,3(1X,D11.4))122CAZ9C.3 List of DELAYC BOUNDARY CONDITION CHANGED TO —DOP1+BVL. OTHERWISE THE SAMEC AS DELAY6. 9/4/91CIMPLICIT REAL*8 (A—H,O—Z)COMMON/COMM/FT (263498)COMMON/ICOM/IFT(4)CDIMENSION YO(512) ,Y1(512) ,GRID(1SO) ,E(512)COMMON/IEL/NPTCOMMON/PARA/VELO(S) ,AMU(3) ,DIFF(3),/ EO,TAU,C1 ,C2,C3,QEPS,DOP1,DOP2,DOP3,PDT,BVL,R.ATIOCOMMON/COOR/XCON 1, XCON2 , XCON3.XCON4NPT 110N4rrOT=N*NPTIFLAG=OEPS1 OD-8NSTARTONEND=0READ(5,2)NSTEP,TAU,DOP1 ,DOP2,DOP3,PHIWRITE(6, 1)NSTEP ,TAU,DOP1 ,DOP2,DOP3,PHI1 FORMAT(’ NSTEP= ‘,IS,/ ‘ TAU2= ‘,D11.4,/,/ ‘DOPINGl= ‘, D11.4, ‘ DOPING2= ‘, D11.4, ‘ DOPING3=’,D11.4,/ ‘ BIAS (V)= ‘, D11.4)2 FORMAT(I5,SD11.4)READ(5 ,3)XCON1 ,XCON2 ,XCON3 ,XCON4WRITE(8, 3) XCON1,XCON2 , XCON3 , XCON43 FORMAT(4D11.4)READ(5,4) BVLRATIOWRITE(6,4) BVL,RATIO4 FORMAT(2D11.4)CC INITIALIZE PROGRAMCT=0.ODOODELT=0.05D00CC RECIPRECAL OF EOCEO=1 ODOO/3 . 4D—2VELO(1)10.ODOOAMU(1)=6O .ODOODIFF(1)=1.554D00VELO(2)=1 ODOOAMU(2)29 ODOODIFF(2)0.751D00AMU(3)=4 ODOODIFP(3)=1.034D—O1C1=4 . ODOO/27 ODOOC22*C1C3=4.ODOO/9.ODOO124POTO . 0259D00*DLOG ((DABS (DDPI*DOP3) /5. 0625D—18) ) +PEIQEPS=1 . 38132087D—02CC REDEFINE TAU TO BE 1. ODOO/TAUCTAU=1 . ODOO/TAUCALL GRIDS((NPT+2) ,GRID)CC INITIAL CONDITIONCDO 10 I=iNPTREAD(4,30) J,TEI4P, YO(I),Y0((I+NPT)),YO((I+2*NPT)),/ YO((I+3*HPT))10 CONTINUE30 FORi4AT(I5,5(1XDi1.4))DO 40 I=1,NTOT40 Yi(i)=YO(I)CALL PCUR(N,NPT,NSTART,NEND,T,GRID,Y0)WRITE (6,4S )NSTART , NEND45 FOPJ(AT(2(IX,I5))TOTDO 500 I=1NSTEPT0 (T0+T) *0. SDOOT=T0+DELT*2. ODOOK=1510 CONTINUECALL FINE(N,NPT,TO,T,YO,Y1,IFLAG,1)IF(K.NE.1) GO TO 520K=2DO 210 L=1,NTOT210 E(L)Y1(L)DO 215 L=1,NPTL1L+3*NPTY0 (LI) =0 5D00* (YO (LI) +Y1 (Li) )215 Yi(Li)=Y0(L1)GO TO 510520 ERROIt=DERR(Y1,E,NTOT)C WRITE(6,201)ERRORC 201 FORMAT(D11.4)DO 60 L=1,NTOT60 E(L)=Y1(L)K=K+iDO 65 L=1,NPTL1L+3*N?TY0(L1)=0.5D00*(YO(LI)+Y1(L1))65 Y1(L1)=Y0(L1)IF(K.GT.8) GO TO 300IF(ERROLGT.EPS) GO TO 510CALL PCUR(N,NPT,NSTART,NEND,T,GRIDJO)DO 70 L=1,NTOTY1(L)=O .5D00*(Y1(L)+YO(L))70 YO(L)=Y1(L)500 CONTINUE125STOP300 WRITE(6,310)310 FORMAT(’TOO MANY ITERATIONS’)STOPENDSUBROUTINE GRIDS (NPT2 ,GD)IMPLICIT REAL*8 (A—H,0-Z)DIMENSION GD(150)COMMON/COOR/ICON1 , XCON2 , XCON3 , XCON4TEMP1= (XCON2—XCON1—0. 15D00)/26 . ODOOTEMP2= (XCQN3—XCON2—0. 3D00) /41. ODOOTEMP3=(XCON4—XCON3—0. 15D00)/26. ODOODO 10 1=1,2510 GD(I)=IC0N1+(I—1)*TEMP1DO 35 1=26,3535 GD(I)=GD(25)+(I—25)*0.03D00DO 20 1=36,7620 GD(I)=GD(35)+(I—35)*TEMP2DO 45 1=77,8645 GD(I)=GD(76)+TEMP2+(I—76)*0.03D00DO 30 1=87,11230 GD(I)=GD(86)+(I—86)*TEMP3RETURNENDSUBROUTINE DEQDDT(N ,DDT)IMPLICIT REAL*8 (A-E,O-z)DIMENSION DDT(4)DDT(1)=1.ODOODDT(2)=1 . ODOODDT(3)=1. ODOODDT(4)=0.ODOORETURNENDSUBROUTINE FLUX(N,X,T,U,UX,G)IMPLICIT REAL*8 (A—R,o—z)DIMENSION U(4) ,uX(4) ,G(4)COMMON/PARA/VELO(S) ,AMu(3) ,DIFF(3),/ E0,TAU,C1 C2,C3,QEPS,DOP1,DOP2,DOP3,POT,BVL,RATIODIMENSION V(3),B(3)CALL VMU(X,U,UX,VB)G(1)—B(1)*UI(1)+V(1)*U(1)G(2)—B(2)*UX(2)+V(2)*U(2)G(3)=—B(3)*UX(3)+V(3)*U(3)G(4)=UX(4)RETURNENDSUBROUTINE FUNR(N,X,T,U,UX,H)IMPLICIT REAL*8 (A—H,O—z)DIMENSION U(4) ,UX(4) ,E(4)COMMON/PAB.A/VELO(3) ,AMU(3) ,DIFF(3),/ E0,TAU,C1 ,C2,C3,QEPS,DOP1,DOP2,DOP3,POT,BVL,RATIODIMENSION V(3),B(3)CALL RATE(X,U,UX,R)126B(1)=TAU*(a*U(1)—U(2))H(3)=0.0DOOE(4)=—QEPS*(U(1)+U(2)—U(3)—DOP(X))RETURNENDSUBROUTINE DGDU(N,X,T,U,UX,DG)IMPLICIT REAL*8 (A—E,O—Z)DIMENSION U(4),UX(4).,DG(2,4,4)COMMON/PARA/VELO(3) ,AMU(3) ,DIFF(3),/ E0,TAU,C1,C2,C3,QEPS,DOP1,DOP2,DOP3,POT,BVL,RATIODIMENSION v(s) ,B(3) ,DV(3) ,DB(3)CALL vMU(x,U,Ux,v,B)CALL DVMU(X,U,U1,DV,DB)DO 10 K=1,2DO 10 1=14DO 10 3=1,410 DG(K,I,J)=O.ODOODG(1,1,1)V(1)DG(1,2,2)V(2)DG(1,3,3)=V(3)DG (2 , 1, 1) =—B (1)DG(2,2,2)—B(2)DG(2,3,3)=—B(3)DG(2,4,4)=1 .ODOODG(2, 1,4)=DB(1)*UI(1)—DV(1)*U(1)DG(2,2,4)=DB(2)*UX(2)—DV(2)*U(2)DG(2,3,4)=DB(3)*UX(3)—DV(3)*U(3)RETURNENDSUBROUTINE DHDU(N,X,T,U,UX,DB)IMPLICIT REAL*8 (A—R,O—Z)DIMENSION U(4),Ux(4),DM(2,4,4)COMMON/PARA/VELO(3) ,AMIJ(3) ,DIFF(3),/ E0,TAU,C1 ,C2,C3,QEPS,DOP1 ,DOP2,DOP3,POT,BVL,RATIODIMENSION v(3),B(3) ,DV(3) ,DB(3)CALL RATE(X,U,UX,R)CALL DRATE(X,U,UI,DR)DO 10 K1,2DO 10 1=1,4DO 10 .7=1,410 DH(K,I,3)=0.ODOODE(1,4, 1)=—QEPSDR(1 ,4,2)=—QEPSDH(1,4,3)=QEPSDE(1, 1, 1)=TAU*RDE (1, 1, 2) =—TAUDE(1,2,1)—DE(1, 1,1)DE(1,2, 2) =—DH(1, 1, 2)DR(2,1 ,4)=—TAU*DR*U(1)DR(2,2,4)=—DR(2, 1,4)RETURNEND127SUBROUTINE BNDRY(N,IFLAG,T ,BV)IMPLICIT P.EAL*8 (A—H10—Z)CDMMON/IEL/NPTCOMMON/PARA/VELO(3) ,ANU(3) ,DIFF(3),/ E0,TAU,C1,C2,C3,QEPS,DQP1,DOP2,DOP3,POT,BVL,RATIODATA PI,T0/3. 14159265359D00,8.ODOO/DIMENSION BV(4)IF(IFLAG.EQ.1) GO TO 10BV(1)=BVL*(1 .0+0 .2D+0O*DSIN(2.ODOO*PI*T/T0))C BV(1)=5.ODOO*(1 .ODOO—DCOS(2.ODOO*PI*T/T0))BV(2)=O.ODOOBV(3)=-DOP1+BVLBV(4)0.ODOORETURN10 BV(2)DOP3/(1 . ODOO+RATIO)BV(1)=DOP3—BV(2)BV(3)=O.ODOOBV(4)=POTRETURNENDSUBROUTINE VMU(X,U,UX1,B)IMPLICIT REAL*8 (A—N,O—Z)DIMENSION U(4) ,Ux(4) ,V(3) ,B(3)COMNON/PARA/VELO(3) ,AMU(3) ,DIFF(3),/ E0,TAU,C1 ,C2,C3,QEPS,DOP1 ,DOP2,DOP3,POT,BVL,RATIODATA CKTQ/0 0259D00/FIELD-UX (4)DO 10 1=1,2BX=AMU(I)*FIELDIF(DABS(BX).GT.(1.5D00*VELO(I))) GO TO 20TEMP=BX/VELO (I)v(I)=—Bx*(1 .ODOO—TEMP*TEMP*C1)GO TO 3020 v(I)=VELO(I)IF(FIELD.GT.0.ODOO) V(I)=—VELO(I)30 B(I)=DIFF(I)IF(FIELD.EQ.0.ODOO) GO TO 10B(I)=—V(I)*CKTQ/FIELD10 CONTINUEv(s)=FIELD*Axu(3)B(3)=DIPF(3)RETURNENDSUBROUTINE RATE(X U,UX,R)IMPLICIT REAL*8 (A—E,O—Z)DIMENSION U(4) ,Ux(4)COMMON/PARA/VELO(3) ,AMU(3) ,DIFP(3),/ E0,TAU,C1,C2,C3,QEPS,DOP1,DOP2,DOP3,POT,BVL,RATIODATA CKTQ,RCKTQ/0 .0259D00, 38. 61003861D00/• DATA TMAX/31.6227766D00/FIELD=-UX (4)TEMP=FIELD*EOIF(DABS(TEMP).GE.TT4AX) GO TO 10128R=TEMP*TEMP*TEMP*TEMPRETURN10 R=1.0D6RETURNENDSUBROUTINE DVMU(X,U,UI,DV,DB)IMPLICIT REAL*8 (A—R,O—z)DIMENSION U(4) ,UX(4) ,DV(3) ,DB(3)COMMON/PARA/VELO(3) ,AMU(3) ,DIFF(3),/ E0,TAU,C1 ,C2,C3,QEPS,DOP1 ,DOP2,DOP3,POT,BVL,R.ATIODATA CKTQ/0 0259D00/FIELD=-TJX (4)DO 10 1=1,2BXFIELD*AMU(I) VIF(DABS(BX).GT.(1.SDOO*VELO(I))) GO TO 20TEMP =BX/VELO (I)DV(I)=—AMU(I)*(1.ODOO—TEMP*TEMP*C3)DB(I)=—DIFF(I)*AMU(I)*AMU(I)*C2*FIELD/(VELO(I)*VELO(I))GO TO 1020 DV(I)0.ODOODB(I)=—VELO(I)*FIELD*CKTQ/(FIELD*FIELD*DABS(FIELD))10 CONTINUEDV(3)=AMU(3)DB(3)=0.ODOORETURNENDSUBROUTINE DRATE(X,U,UX,DR)IMPLICIT REAL*8 (A—B,O-z)DIMENSION U(4) ,Ux(4)COMMON/PARA/VELO(3) ,AMU(3) ,DIFF(3),/ E0,TAU,C1 ,C2,C3,QEPS ,DOP1,DOP2,DOPS,POT,BVL,R.ATIODATA CKTQ,RCKTq/0 .0259D00, 38.61003861D00/DATA TMAX/31 6227766D00/FIELD =—UX (4)TEMP=FIELD*E0IF(DABS(TEMP).GE.TMAX) GO TO 10DR=4 00000*TEMP*TEMP*TEMP*E0RETURN10 DRO.ODOORETURNENDFUNCTION DOP(X)IMPLICIT REAL*8 (A-E,O-Z)COMMON/PARA/VELO(3) ,AMU(3) ,DIFF(3),/ E0,TAU,C1 ,C2,C3,QEPS,DOPI,DOP2,DOP3,POT,BVL,RATIOCOMMON/COOR/XCON1 XCDN2 , ICON3 XCON4IF(X.LE.XCON2) GO TO 10IF(X.LE.XCON3) GO TO 20 VDOP=DOP3IF(X.LT. (XCONS+O. 1SDOO)) DOP=0.SDOO*(DOP3+DOP2)+/ 0.SDOO*(DOP3—DOP2)*(X—XCON3)*6 .666666667D00RETURN10 DOP=DOP1129IF(X.GT. (xC0N2—0 1SDOO)) DOP=DOP1*(XCON2—X)*6.666666667D00RETURN20 DOP=DOP2IF(X.LT. (XCON2+0. 15D00)) DOP=D0P2*(X—XCON2)*6. 666666667D00IF(X GT. (XCON3—0. ISDOO)) DOP=0. SDOO*(DOP3+DOP2)+/ 0 .5D00*(D0P2—DOP3)*(XCON3—X)*6 .666666667D00RETURNENDFUNCTION DERR(Y,E,N)IMPLICIT REAL*8 (A—E,O—Z)DIMENSION Y(512) ,E(512)DEL=0.ODOODMOD=0.ODOODO 10 11,NDMOD=DMOD+Y (I) *y (I)TEMP=E(I)—Y(I)10 DEL=DEL+TEMP*TEMPDERR=DSQRT (DEL/DMOD)RETURNENDSUBROUTINE PCUR(N,NPT,NSTART,NEND,T,GRID)Y1)IMPLICIT REAL*8 (A—H,O—z)COMMON/PAP.A/VELO(3) ,AMIJ(3) ,DIFF(3),/ E0,TAU,C1,C2,C3,QEPSDOP1,DOP2,DOP3,POT,BVL,RATIQDIMENSION GRID(150) ,Y1(512)DIMENSION U(4),Ux(4) ,G(4)DIMENSION x(150) ,CRNT(1S0) ,FIELD(150)N1=NPT—1DO 30 I=1N1.x(I)=0.SDOO*(GRID((I+2))+GRID((I+1)))U(1)=0.5D00*(Y1((I+1) )+Y1(I))U(2)=0 .SDOO*(Y1((I+NPT+1))+Y1((I+NPT)))U(3)=0 LDOO*(Y1((I+NPT*2+1))+Y1((I+2*NPT)))U(4)=o 5D00*(Y1((I+NPT*3+1))+Yi((I+3*NPT)))DELX=1.ODOO/(GRID((I+2))—GRID((I+1)))UX(1)=(Y1((I+1))—Y1(I))*DELXUX(2)=(Y1((I+NPT+i))—Y1((I+NPT)))*DELXUx(s)=(Y1((I+2*NPT+1))—Y1((I+2*NPT)))*DELIUX(4)=(Y1((I+3*NPT+1))—Y1((I+3*NPT)))*DELXCALL FLUX(3,X,T,U,UX,G)CRNT(I)=G(3)—G(1)—G(2)FIELD(I)=—UX(4)30 CONTINUECC SELECT POINTS NSTART AND lENDCIF(HEND.NE.0) GO TO 60NSTART=1SNEND=19GO TO 60I=0FRAC=—DOP1*0 O1DOO40 CONTINUE1301=1+1ITEMP=NPT*2+IIF(Y1(ITEZ4P).GT.FRAC) GO TO 40NSTART=II=NPTFR.AC=DOP3*O. 0 1DOO/ (1 ODOO+RATIO)50 CONTINUE1=1—1ITENP=NPT+IIF(Y1(ITEMP).GT.FRAC) GO TO 50NEND=I60 CONTINUECC CALCULATE TALTOAL CURRENTSCANO.0DO 70 I=NSTART,NENDAN=AN+0 .5D00*(X((I+1))—X(I))*(CRNT((I+1))+CRJT(I))70 CONTINUEAN=AN/(X((NEND+1))—X(NSTART))WRITE(6,80)T,CR.JT(NSTART) ,CP.JT((NEND+1)),AN80 FORI1AT(4(1X,D14.7))RETURNEND131C.4 List of HBTFLC USED FOR HIGH CURRENT SIMULATION 10/9/91CC 10/12/91 FIX OUTPUT SUBROUTINECC 10/23/91 MODIFIED BOUNDARY CONDITIONS FOR HOLESCIMPLICIT REAL*8 (A—B, O-z)DIMENSION IEL(240,3) ,ID(480)DIMENSION A(480,20) ,R(480),/IE(3),EM(6,6),F(6)COMMON/NEL/NPT,NFREECOMMON/EL/PHIO(240) ,PHII(240) ,ELCO(480) ,ELC1(480) ,X(240),BV(480)/,HOLO(240) ,HOL1 (240)CC COMMON BLOCK FOR PARAMETERSCCOMMON/PARAM/VELO , AMU , VELOP , AMUP,VELO2 , AMU2,/ E0,TAU,C1,C2,C3,QEPS,DOP1,DOP2,DOP3,POT,BIAS,BVL,RATIOCOMMON/COOR/XCON1 , XCON2,XCON3, XCON4COMMON/PART/TO, T ,DELTCC INPUT DATACCALL DATAINCC INITIALIZE PARAMETERSCCALL INITI(NEQ,NPT,NFREE,ND,NE,NN)CC DEFINE NODESCCALL NODES(NEQ,NPT,X,ID)DO 40 I=1,NPTDPP=DOP(X(I))IF(DPP.LT.0.ODOO) DPPO.ODOOI11*2—11211+1ELCO(I1)DPP/(1 .ODOO+RATIO)ELCO(12)=DPP—ELCO(I1)ELC1(I1)=DPP/(1 .0DO0+PATIO)ELC1(12)=DPP—ELC1(I1)40 CONTINUEDO 50 I=1,NPTDPP=DOP(X(I))IF(DPP.GT.0. ODOO) DPP=0.ODOOBOLO(I)=—DPPHOL1(I)=—DPP50 CONTINUECC DEFINE ELEMENTCCALL ELEM(ND,NN,IEL)132CC CALCULATE THE SEMIBAND,WIDTH AND NCOLSCCALL SMBAND(ND,NN,NFREE,MBAND0,NCOLSO,IEL)NFREE=2CALL SMBAND(ND,NN,NFREE,MBAND1,NCOLS1,IEL)CC SOLVE THE SYSTEMCEPS=1 OD—4DO 10 1=1,8000ITYPE=ONFREE=1NEQ=NPT*NFREENENN*NFREEMBAND=MBANDONCOLS=NCOLSOID(1)=1ID(NEQ)=1CALL FINITE(T0,T,NPT,NEQ,NCOLS,MBAND.ND,NE,NN,NFREE,/IEL,ID,A,R,PHIO,PHI1,X,BV,IE,EM,F,ITYPE)ER1=ERROR (NEQ , PHIO , PHIl)DO 30 J=1,NEQIF(I.EQ.1) GO TO 35PRI0(J)=(PHII(J)+PfiI0(J))*0. SDOOPHI1(J)=PHIO(J)GO TO 3035 CONTINUEPHIO(J)PHI1(J)30 CONTINUEID(1)=0ID(NEQ)=0ITYPE=1NFREE=2NEQ=NPT*NFREENE=NN*NFREEMBAND=MBAND 1NCDLS=NCOLS 1ID(1)=1ID(2)=1ID((NEQ—1) )=1ID(NEQ)=I.CALL FINITE(TO,T,NPT,NEQ,NCOLS,MBAND,ND,NE,NN,/NFREE,IEL,ID,A,R,ELCO,ELC1,X,BV,IE,EM,F,ITYPE)ER2=ERROR(NEQ , ELCO , ELC1)DO 20 J=1,NEQELCO(J)=0.5D00*(ELC1(J)+ELCO(3))ELC1(J)=ELCO(J)C ELCO(J)=ELC1(J)20 CONTINUEID (1) =0ID (2 ) =0ID((NEQ—1))=0133ID(NEQ)=0ITYPE=2NFREE 1NEQ=NPT*NFREENE=NN*NFREEMBAND=MBANDONCOLS=NCOLSOID (1) = 1ID(NEQ)=1CALL FINITE(T0,T,NPT,NEQ,NCOLS,MBAND,ND,NE,NN,/NFREE,IEL,ID,A,R,ROLO,MOL1,X,BV,IE,EM,F,ITYPE)ER3=ERROR C NEQ , BOLO , ROLl)DC 60 J=1,NEQROLO(J)0. SDOO*(HOL1(J)+BOLO(J))HOL1(J)110L0(J)C ROLO(J)=ROL1(3)60 CONTINUEID(1)0ID(NEQ)0T00 6D00*(T0+T)C TOTT=T0+DELTIF((ER1.LT.EPS).AND.(ER2.LT.EPS).AND.(ER3.LT.EPS)) GO TO 7010 CONTINUE70 CONTINUEWRITE(6,90)I,ER1 ,ER2,ER390 FORMAT(15,1X,3(D11.4,1X))CC PRINT SOLUTIONCCALL OUTPUT(NEQ,NPT,NFBEE,ND,NN,IEL)CCTERMINATE PROGRAMCSTOPENDCCCSUBROUTINE INITI(NEQ,NPT,NFREE,ND,NE,NN)IMPLICIT REAL*8 (A—H, O-z)COMMON/PARAM/VELO,AMU,VELOP,AMUP,VELO2SAMU2,/ E0,TAU,C1,C2,C3,QEPS,00P1,DOP2,DOP3,POT,BIAS,BVL,RATIOCOMMON/COOR/XCON 1 • XCON2 , XCON3 , XCON4COMMON/PART/TO ,T,DELTCOMMON/EL/PHIO(240) ,PHI1(240) ,ELCO(480) ,ELC1(480) ,x(240) ,Bv(480)/,ROLO(240) ,ROL1(240)NPT239NFREE=1NEQ=NPT*NFREEND(NPT—1)/2NN=3134NE=NN*NFREEVELO=1 . ODO1AMU=60.ODOOVELO2=1 . ODOOAMU2=29.ODOOVELOP=1 . ODOOAMUP=4.09D00EO=1 . ODOO/3. 4D—2TO=O.ODOOT=O . 0020D00DELT=T-TOTAU=1 . ODOO/TAUPOT=O . 0259000*DLOG((DABS(DOP1*DOP3)/5. 0625D—18))+BIASC 1=4. ODOO/27 . ODOOC2=2*C1C3=4.ODOO/9.ODOOQEPS=1 . 38132087D—02CALL ZARRAY(NEQ,PHIO)CALL ZARRAY(NEQ,PHI1)CALL ZARRAY(NEQ ,BV)RETURNENDCCCSUBROUTINE DATAINIMPLICIT REAL*8 (A—H, O—z)COMMON/PARAM/VELO , AMU ,VELOP,AMUP,VELO2 , AMU2,/ EO,TAU,C1,C2,C3,QEPS,DOP1,DOP2,DDP3,POT,BIAS,BVL,RATIOCOMMON/COOR/XCON1 , XCON2 , XCON3 , XCON4READ(5,2)NSTEP,TAU,DOP1 ,DOP2,DOP3,BIASWRITE(6, 1)NSTEP,TAU,DOP1 ,DOP2,DOP3,BIASI FORMAT(’ NSTEP= ‘,IS,/ ‘ TAU2= ‘,D11.4,/,/ ‘DOPINGl= ‘, D11.4, ‘ DOPING2= ‘, D11.4. ‘ DOPING3=’,D11.4,/ ‘ BIAS (V)= ‘, D11.4)2 FORMAT(I5b5D11.4)READ(5, 3)XCON1 ,XCON2 ,XCON3 ,XCON4WRITE(6,3)XCON1 ,XCON2 ,XCON3,XCON43 FORMAT(4D11.4)READ(5,4) BVL,RATIOWRITE(6,4) BVL,RATIO4 FORMAT(2D11.4)RETURNENDCC DEFINE NODES AND SET FLAGS FOR BOUNDARY CONDITIONSCSUBROUTINE NODES (NEQ , NPT , X ,ID)IMPLICIT REAL*8 (A—H, O-z)DIMENSION X(NPT) ,ID(480)COMMON/COOR/XCON1 ,XCON2 ,XCON3,XCON4DO 60 1=1,480135ID(I)=060 CONTINUEID(1)1ID(NEQ)=1x(1) =XCON 1x (NPT) =XCON4N1NPT-1DX=(X(NPT)—X(1))/(NPT—1)DO 50 I=2,N1x(I)=X(1)+(I—1)*DX50 CONTINUEC XX1=ICON1C XX2=XCON2—0.2D00C 1X3=XCON2+0. 2D00C 1X4=ICON3-0 2D00C XXS=XCON3+0 . 2D00C XX6=XCON4C DX=(XX2—XX1)/20 ODOOC DO 10 1=1,21C x(I)=xCONI+(I—1)*DxClO CONTINUEC DX= (XX3—Xx2 ) /29. ODOOC DO 20 1=22,50C x(I)=x(21)+(I—21)*DXC20 CONTINUEC DX=(XX4—XX3)/21 .ODOOC DO 30 1=51,71C x(I)=X(50)+(I—50)*DXC30 CONTINUEC DX=(XX5—XX4)/29 . 0000C DO 40 1=72,100C X(I)=X(71)+(I—71)*DXC40 CONTINUEDX= (xx6—Xx5 ) /20. 0000C DO 50 1101,240C x(I)=X(100)+(I—100)*DXC50 CONTINUERETURNENDCCCSUBROUTINE FUNCF(MNF,NFREE,T,X,IE,U,UX,F,ITYPE)IMPLICIT REAL*8 (A—H, O—Z)DIMENSION IE(3) ,U(MNF) ,Ux(MNF) ,F(MNF)COMMON/PARAM/VELO, AMU , VELOP, AMUP , VELO2 , AMU2,/ E0,TAU,C1 ,C2,C3,QEPS,00P1,00P2,DQP3,POT,BIAS,BVL,R.ATIOCOMMON/COOR/XCON1 ,XCON2,XCON3,XCON4COMMON/PART/TO ,T1 ,DELTIF(ITYPE.EQ.0) GO TO 10IF(ITYPE.EQ.1) GO TO 20F(1)0.ODOORETURN13610 CONTINUEF(i)=—QEPS*(DENS(MNF,X,IE)—DOP(X))RETURN20 CONTINUEF(1)=0.ODOOF(2)=0.ODOOENDCCCSUBROUTINE FUNCG(MNF,NFREE,T,X, IE,U,UX,G,ITYPE)IMPLICIT REAL*8 (A—H ,O—Z)DIMENSION IE(3) ,U(MNFY,UX(MNF),G(MNF)COMMON/PARAM/VELO , AMU, VELOP, AMUP , VELO2 , AMU2,/ E0,TAU,C1,C2,C3,QEPS,DOP1,DOP2,DOP3,POT,BIAS,BVL,RATIOCOMMON/COOR/XCON1 , XCON2 , XCON3 , XCON4COMMON/PART/TO, Ti, DELTG(1)=0.ODOOIF(ITYPE.EQ. 1) G(2)=0. ODOORETURNENDCCCSUBROUTINE DFDU(MNF,NFREE,T,X,IE,U,UX,F, ITYPE)IMPLICIT REAL*8 (A—H, 0-z)DIMENSION IE(3) ,U(MNF) ,UX(MNF) ,F(2,MNF,NNF)COMMON/PARAM/VELO , AMU ,VELOP , AMUP , VELO2 , AMU2,/ EO,TAU,C1,C2,C3,QEPS,DOP1,DOP2,DOP3,POT,BIAS,BVL,RATIOCOMMON/COOR/XCON1 ,XCON2,XCON3 ,XCON4COMMON/PART/TO ,T1 ,DELTIF(ITYPE.EQ.1) GO TO 10F(i,i,1)O.ODOOF(2,i,1)0.ODOORETURN10 CONTINUER=RATE(MNF,X,IE)F(1,i,1)=R*TAUF(1,1,2)=—TAUF(i,2,1)=—F(i,i,1)F(1,2,2)—F(i, 1,2)F(2,i,i)=O.ODOOF(2,1,2)0.0000F(2,2, 1)=O.ODOOF(2,2,2)=O.ODOORETURNENDCCCSUBROUTINE DGDU(MNF,NFREE,T,X,IE,U,UX,G, ITYPE)IMPLICIT REAL*8 (A—H, O—Z)DIMENSION IE(3) ,U(MNF) ,UX(MNF) ,G(2,MNF,MNF)137COMMON/PARAM/VELO , AMU ,VELOP, AMUP , VELO2 , AMU2,/ E0,TAU,C1,C2,C3,QEPS,00P1,DOP2,DOP3,POT,BIAS,BVL,P.ATIOCOMMON/COOR/XCON1 , XCON2 • XCON3 , XCON4COMMON/PART/TO, Ti , DELTIF(ITYPE.EQ.O) GO TO 10CALL VMU(ITYPE,MNF,X,IE,D,V,D2,V2)IF(ITYPE.EQ.i) GO TO 20G(1,1,i)=VG(2,i, i)=—DRETURN10 CONTINUEG(1,i,i)=O.ODOOG(2, 1,1)=1 .ODOORETURN20 CONTINUEG (1, 1 , 1) =VG(2,i,i)=—DG(i,2,2)=V2G(2,2.2)=—D2G(1,1,2)=0.ODOOG(1,2,i)0.ODOOG(2,1,2)=O.ODOOG(2,2,i)0.ODOORETURNENDCCCSUBROUTINE BCOND(T,NEQ ,NPT,NFREE,ID,BV,ITYPE)IMPLICIT REAL*8 (A—H, O—z)COMMON/PARAM/VELO , AMU , VELOP,AMUP , VELO2 , AMU2,/ E0,TAU,C1,C2,C3,QEPS,DOP1,DOP2,DOP3,POT,BIAS,BVL,RATIOCOMMON/COOR/XCON 1, XCON2, XCON3 , XCON4COMMON/PART/T2 ,Ti , DELTCOMMON/EL/PHIO(240),PHI1(240),ELCO(480)fELC1(480),X(240),BV1(480)/,HOLO(240) ,HOL1 (240)DIMENSION ID(NEQ) ,BV(NEQ)DATA P1/3. 141592654000/IF(ITYPE.EQ.O) GO TO 10IF(ITYPE.EQ.2) GO TO 20C BV(1)=1.ODOO—DCOS(T*PI)C BV(i)=TBV(i)=BVL*T/2. ODOOIF(T.GT.2.ODOO) Ev(i)=BVLBV(2)=0.ODOOBV((NEQ—i))D0P3/(1 .ODOO+RATIO)BV(NEQ)=DOP3—BV((NEQ—i))C BV((NEQ—i))=ELCO((NEQ—3))C BV(NEQ)=ELCO((NEQ—2))RETURN10 CONTINUEBV(1)=0.ODOOBV(NEQ)=POT*T/2 ODOO138CCC1211+lWRITE(6,20)I,X(I),ELCO(I1),ELCO(12),HOLO(I),PHIO(I)10 CONTINUEDO 40 11,NDDO 50 J=1,NNIE(J)=IEL(I,J)50 CONTINUEITYPE= 1NFREE=2NEQ=NPT*NFREENENN*NFREECALL CURRENT(TO,T,I,NEQ,NPTNE,NN,NFREE,IE,EM,F,X,ELC0,ELC1,ITYPE)CNT1=—F(1)CNT2=—F(2)CNT3=F(5)CNT4=F(6)ITYPE=2NFREE1NEQ=NPT*NFREENE=NN*NFREECALL CURRENT(T0,T,I NEQ,NPT,NE,NN,NFREE,IE,EM,F,X,HOL0,HOL1 ,ITYPE)CNT5=-F(1)CNT6=F(3)TEMP 1=CNT5— (CNT1+CNT2)TEMP2=CNT6- (CNT3+CNT4)FIELD1=FVAL(X(IE(1)) ,IE)FIELD2=FVAL(X(IE(2)) , IE)WRITE(6,20)I,X(IE(1)) ,TEMP1,FIELD1WRITE(6,20)I ,x(IE(2)),TEMP2,FIELD240 CONTINUERETURN20 FORMAT(15,1X,10(D14.7,1X))IF(T.GT.2.ODOO) BV(NE)=POTC BV(1)=PHIO(2)C BV(NEQ)=PHI0((NEQ—1))RETURN20 CONTINUEBV(1)=—DOP1+BVL*T/2 0000IF(T.GT.2.ODOO) BV(1)=—DOP1+BVLBV(NEQ)0. ODOOENDOUTPUT HANDLES UP TO 9 PDE’SSUBROUTINE OUTPUT(NEQ,NPT,NFREE,ND ,NN,IEL)IMPLICIT REAL*8 (A—H • o-z)DIMENSION IEL(ND,NN)COMMON/EL/PHIO(240) ,PHI1(240) ,ELCO(480) ,ELC1(480) ,X(240) ,BV(480)/,UOLO(240) ,EOL1(240)DIMENSION IE(3) ,EM(12,12),F(12)NM,T,DELTDO 10 I=1,NPT11=1*2—113930 FORMAT(/)ENDCCCFUNCTION DENS (MNF , XX, IE)IMPLICIT REAL*8 (A—H, O—z)DIMENSION IE(3)CONJ4ON/EL/PHIO(240),PHI1(240) ,ELCO(480),EL.C1(480),X(240),BV(480)/,HOLO(240) ,HOL1 (240)I1=IE(1)*2—1I21E(2)*2—113=IE(3)*2—114=11+115=12+116=13+1DENS=(ELCO(I1)+ELCO(14))*SHAPE(1,X(IE(1)),X(IE(2)),X(IE(3)),XX)+/(ELCO(12)+ELCO(I5))*SHAPE(2,X(IE(1)),JC(IE(2)),X(IE(3)),XX)-f/(ELC0(I3)+ELCO(16))*SEAPE(3,X(IE(1)) ,X(IE(2)),x(IE(3)) ,xX)DENS=DENS—HOLO(IE(1))*SHAPE(1,X(IE(1)),X(IE(2)),X(IE(3)),XX)—/ROLO(IE(2))*SHAPE(2,X(IE(1)),X(IE(2)),X(IE(3)),XX)—/HOLO(IE(3))*SHAPE(3,X(IE(1)),X(IE(2)),X(IE(3)),XX)RETURNENDCCCSUBROUTINE VMU(ITYPE,MNF,XX,IE,D,V,D2,V2)IMPLICIT REAL*8(A-R, O-z)DIMENSION IE(3)COMMON/PARAM/VELO , AMU , VELOP , AMUP , VELO2 , AMU2,/ E0,TAU,C1 ,C2,C3,QEPS,DOPI ,DOP2,00P3,POT,BIAS,BVL1RAT OCOMMON/COOR/XCON 1, XCON2 , XCON3 , XCON4COMMON/EL/PHI0(240),PHI1(240),ELCO(480),ELC1(480),X(240),BV(480)/,HOLO(240) ,HOL1(240)DATA CKTQ/0 . 0259D00/FIELD=FVAL(XX,IE)IF(ITYPE.EQ.2) GO TO 40BXAMU*FIELDIF(DABS(BX) . GT. (1. 5D00*VELO)) GO TO 20TEMP=BX/VELQV=-BX* (1. ODOO—TEMP*TEMP*C 1)GO TO 3020 V=VELOIF(FIELD. GT .0. ODOO) V=-VELO30 D=CKTQ*AMUIF(FIELD.EQ.0.ODOO) GO TO 10D=-V*CKTQ/FIELD10 CONTINUEBXAMU2*FIELDIF(DABS(BX).GT.(1.SDOO*VELO2)) GO TO 80TEMP=BX/VELO2V2—BX*(1 . ODOO—TEMP*TEMP*C1)140GO TO 9080 V2=VELO2IF(FIEL.D. GT.0. 0D00) V2=-VELO290 D2=CKTQ*AMU2IF(FIELD.EQ.0.ODOO) GO TO 100D2=-V2*CKTQ/FIELD100 CONTINUERETURN40 CONTINUEBXAMUP*FIELDIF(DABS(BX) . GT. (1. 5D00*VELOP)) GO TO 60TEMP=BX/VELOPV=BX* (1. ODOO—TEMP*TEMP*C1)GO TO 5060 V=VELOPIF(FIELD.LT. 0. ODOO) V=-VELOP50 DCKTQ*AMUPIF(FIELD.EQ.0.ODOO) GO TO 70D=V*CKTQ/FIELD70 CONTINUERETURNENDCCCFUNCTION RATE (MNF ,XX, IE)IMPLICIT REAL*8(A—H,O—Z)DIMENSION IE(3)COMMON/PARAM/VELO, AMU , VELOP , AMUP , VELQ2 , AMU2,/ E0,TAU,C1,C2,C3,QEPS,DOP1 ,00P2,DOP3,POT,BIAS,BVL,RATIOCOMMON/COOR/XCON1 , XCON2 , XCON3 , XCON4COMMON/EL/PHIO(240),PHI1(240),ELCO(480),ELC1(480),X(240),BV(480)/,HOLO(240) ,MOL1 (240)DATA TMAX/31.6227766D00/FIELD=FVAL(XX,IE)TEMP=FIELD*E0IF(DABS(TEMP).GE.TMAX) GO TO 10RATE=TEMP*TEMP*TEMP* TEMPRETURN10 RATE=1.0D6RETURNENDCCCFUNCTION DOP(XX)IMPLICIT REAL*8 (A—H, O—z)COMMON/PARAM/VELO , AMU , VELOP , AMUP, VELO2 , AMU2,/ E0,TAU,C1 ,C2,C3,QEPS,DOP1 ,DOP2,DOP3,POT,BIAS,BVL,RATIOCOMMON/COOR/XCON1 ,XCON2,XCON3,XCON4XX2XCON2—0. 2000XX3=XCON2+0 2000XX4XCON3—0. 2000141XXS=XCON3+0.2D00IF(XX.GT.XX5) GO TO 10IF(XX.GT.XX4) GO TO 20IF(XX.GT.XX3) GO TO 30IF(XX.GT.XX2) GO TO 40DOP=DOP 1RETURN10 CONTINUEDOP=DOP3RETURN20 CONTINUEDOPDOP2+ (DOP3—DOP2) * (XX—xX4) *2. 5D00RETURN30 CONTINUEDOP=DOP2RETURN40 CONTINUEDOP=DOP 1+ (DOP2—DOP 1) * (XX—XX2) *2. 5D00RETURNENDCCCFUNCTION ERROR(NEQ,Y0,Y1)IMPLICIT R.EAL*8(A-R,O—Z)DIMENSION Y0(NEQ) ,Y1(NEQ)E=0.ODOODO 10 I=1,NEQTEMPYO(I)—Y1(I)E=E+TEMP*TEMP10 CONTINUEERROR=DSQRT (E/NEQ)RETURNENDCCCFUNCTION FVAL(XX • IE)IMPLICIT REAL*8 (A—H, o—z)DIMENSION IE(3)COMMON/EL/PHIO(240),PHI1(240),ELCO(480),ELC1(480),X(240),BV(480)/,HOLO(240) ,EOL1 (240)FVAL=—PHIO(IE(1))*DSHAPE(1 ,X(IE(1)) ,X(IE(2)) ,X(IE(3)) ,xx)—/PHIO(IE(2))*DSHAPE(2,X(IE(1)) ,X(IE(2)) ,x(IE(3)),xx)—/PHIO(IE(3))*DSHAPE(3,X(IE(1)) ,x(IE(2)) ,X(IE(3)) ,XX)RETURNENDCCC$INCLUDE ASSEMC$INCLUDE GAUSSC$INCLUDE FINITE1C$INCLUDE ONED3C142C.5 List of FINITE1CC * INTERFACE TO PARTIAL DIFFERENTIAL EQUATIONS *CC SUBROUTINE CALEL RETURNS ELEMENT MATRIX AND VECTOR.C SUBROUTINE BCOND RETURNS BOUNDARY CONDITIONS FOR NODESC WITH FREEDOM NFREE.CC THE REST OF SUBROUTINES IN THIS FILE SHOULD BE INDEPENDENT OFC PDE’S.CCC ENTRY TO TIME TOc T TIME T (>To)C NPT NUMBER OF NODESC NEQ NUMBER OF EQUATIONSC NCOLS NUMBER OF COLUMNSC MBAND SEMIBANDWIDTHC ND NUMBER OF ELEMENTSC NE DIMENSION OF ELEMENT MATRIX AND VECTORC NN NUMBER OF NODES IN THE ELEMENT (<=4)C NFREE FREEDOMS OF NODESC IEL DEFINITION OF ELEMENTS (NODES IN THEM)C ID FLAGS FOR BOUNDARY CONDITIONSC A SYSTEM MATRIXC R SYSTEM VECTORC YO SOLUTION AT TOC Yl SOLUTION AT T (GUESS SOLUTION)C X COORDINATES OF NODESC BV BOUNDARY CONDITIONSC IE ELEMENT NODESC EM ELEMENT MATRIXC F ELEMENT VECTORC ITYPE TYPE OF EQUATIONS (GROUPS)CC EXIT Yl SOLUTION AT TCCSUBROUTINE FINITE(T0 ,T,NPT,NEQ ,NCOLS,MBAND,ND,NE,NN,/NFREE,IEL,ID,A,R,YO,Y1 ,X,BV,IE,EM,F,ITYPE)IMPLICIT REAL*8 ,O—Z)DIMENSION IEL(ND,NN) ,ID(NEQ)DIMENSION A(NEQ,NCOLS),R(NEQ),YO(NEQ),Y1(NEQ),X(NPT),BV(NEQ)DIMENSION IE(NN) ,EM(NE,NE) ,F(NE)CCDO 20 K=11CC ZEROES A,R, AND BVCCALL ZERO(NEQ,NCOLS,A,R)CC ASSEMBLE MATRIX k VECTOR143CCALL AMATRIX(T0,T,NPT,NEQ,NCOLS,MBAND,ND,NE,NN,/NFREE,IEL,ID,A,R,YO,Y1,X,IE,EM,F,ITYPE)CC INSERT BOUNDARY CONDITIONSCCALL BNDRY(T0,T,NPT,NMAX,NEQ,NCOLS,MBAND,NE,NN,NFREE,/ID,BV,A,R,ITYPE)CC SOLVE SYSTEM AND PRINT RESULTSCCALL SYSTEM(NEQ,NPT,NFREE,NCOLS,MBAND,A,R,X)CC UPDATE Yl-CDO 10 I=1,NEQY1(I)=R(I)10 CONTINUE20 CONTINUERETURNENDCCCSUBROUTINE AMATRIX(T0,T,NPT,NEQ,NCOLS ,MBAND,ND,NE,NN,/NFREE,IEL,ID,A,R,Y0,Y1,X,IE,EM,F,ITYPE)IMPLICIT REAL*8 (A—B,O—Z)DIMENSION IEL(ND,NN) ,ID(NEQ)DIMENSION A(NEQ,NCOLS),R(NEQ),YO(NEQ),Y1(NE),X(NPT)DIMENSION IE(NN) ,EM(NE,NE)DIMENSION F(NE)DO 10 I=1,NDDO 20 J=1,NNIE(J)=IEL(I,J)20 CONTINUECALL CALEL(T0,T,I,NEQ,NPT,NE,NN,NFREE,IE,EM,F,x,Yo,Y1,ITYPE)CALL ASSEM(NEQ,NCOLS ,MBAND,NEQ,NCOLS ,NE,NN,NFREE,IE,EM,F,A,R)10 CONTINUERETURNENDCCCSUBROUTINE SYSTEM(NEQ , NPT ,NFREE, NCOLS ,MBAND, A ,R,X)IMPLICIT REAL*8 (A—H, O—z)DIMENSION A(NEQ,NCOLS) ,R(NEQ) ,x(NPT)CALL FACTOR(NEQ ,NCOLS ,MBAND ,NEQ , NCOLS ,A, IPUNT)CALL SOLVE(NE,NCOLS,MBAND,NEQ ,NCOLS,A,R, IPUNT)RETURNENDCCC144SUBROUTINE BNDRY(T0,T, NPT,NMAX ,NEQ,NCOLS,/I4BAND,NE,NN,NFREE,ID,BV,A,R,ITYPE)IMPLICIT REAL*8(A—H,O—Z)DIMENSION ID(NEQ),BV(NEQ),A(NEQ,NCOLS),R(NEQ)C TT=0.5D00*(T0+T)CALL BCOND(TNEQ,NPT,NFREE,ID,BV,ITYPE)DO 30 I=1,NEQIF(ID(I).EQ.0) GO TO 30DO 20 JJ=1,NCOLSII=I+MBAND—JJIF(II.LT.1) GO TO 10IF(II.GT.NEQ) GO TO 10IF(A(II,JJ).EQ.0.ODOO) GO TO 10R(II)=R(II)—BV(I)*A(II,JJ)A(II,JJ)=o.ODOO10 CONTINUEIF(A(I,33).EQ.0.ODOO) GO TO 20A(I,J3)=0.ODOO20 CONTINUER(I)=BV(I)A(I ,MBAND)=1 ODOO30 CONTINUERETURNENDCCCSUBROUTINE ZERO(NEQ ,NCOLS ,A,R)IMPLICIT REAL*8 (A—B, O—Z)DIMENSION A(NEQ,NCOLS) ,R(NEQ)DO 10 I=1,NEQR(I)=0.0000DO 10 J=1,NCOLSA(I,J)=0.ODOO10 CONTINUERETURNEND145C.6 List of ONED3CSUBROUTINE EL.EM(ND,NN, IEL)IMPLICIT REAL*8 (A—H, o—z)DIMENSION IEL(ND,NN)DO 10 I=1,NDDO 10 3=1,NNJ 1 = (NN— 1) * (I—i) +JIEL(I,J)J110 CONTINUERETURNENDCCcSUBROUTINE SMBAND(ND ,NN,NFREE,MBAND ,NCOLS, IEL)IMPLICIT REAL*8 (A—H, o—z)DIMENSION IEL(ND,IIN)MBAND=0DO 10 I=1,NDMIN=100000MAX=0DO 20 J1,NNIF(IEL(I,J).GT.MAX) MAX=IEL(I,3)IF(IEL(I,J).LT.MIN)MIN=IEL(I,J)20 CONTINUENDIF=MAX-NIN+ 1IF (NDIF . GT MBAND )MBAND=NDIF10 CONTINUEMBAND=MBAND*NFREENCOLS=2*MBAND— 1RETURNENDCCCSUBROUTINE ZARRAY(NEQ ,A)IMPLICIT REAL*8 (A—H, O-Z)DIMENSION A(NEQ)DO 10 I=1,NEQA(I)=0.ODOO10 CONTINUERETURNENDCC CALEL HANDLES NFREE<=4CSUBROUTINE CALEL(T0,T,M,NEQ ,NPT,NE,NN ,NFREE,/IE,EM,F,X,Y0,Y1,ITYPE)IMPLICIT REAL*8(A—H, O—Z)DIMENSION IE(NN) ,EM(NE,NE)DIMENSION F(NE)DIMENSION X(NPT) ,Yo(NEQ) ,Y1(NEQ)DIMENSION G1(2,4,4),G2(2,4,4),146/G3(2,4,4) ,G4(2,4,4) ,G5(2,4,4),F1(2,4,4) ,F2(2,4,4),/F3(2,4,4),F4(2,4,4),F5(2,4,4),H1(4),H2(4),H3(4),/H4(4) ,H5(4) ,H6(4) ,H7(4) ,H8(4) ,H9(4) ,H1O(4)DIMENSION U1(4) ,U2(4) ,U3(4) ,U4(4) ,U5(4)DIMENSION UX1(4) ,UX2(4) ,UX3(4) ,UX4(4) ,UX5(4)DIMENSION EM1(12,12),EM2(12,12)CC MATRIX DEFINED IN A REVERSED WAYC SO THE MATRIX EQUATION IS (—A)*U=(—B).C WHEN ADDING NEW FEATURES TO THIS SUBROUTINE,C BE AWARE OF THISCI4NF=4M1=IE(1)M2=IE(2)M3=IE(3)TT=O 5D00*(TO+T)X1=X(M1)X5=X(M3)DX=(X5—X1)*O 25D00X2=X1+DXX3=X2+DXX4=X5-DXDELX1=X(M2)-X(M1)DELX2=X(M3)—X(M1)DELX3=X (M3) -X(M2)RDELT=1 .ODOO/(T—TO)TEMP2=(DELX1+DELX3)/90. ODOOCALL UVAL(MNF,NEQ,NFREE,IE,X1,X2,X3,X4,X5,/YO,Y1 ,U1 ,U2,U3,U4,U5)CALL UXVAL(MNF,NEQ,NPT,NFREE,IE,X1,X2,X3,X4,X5,/YO,Y1 ,UX1 ,UX2,UX3,UX4,UX5)CALL DGDU(MNF,NFREE,TT,X1,IE,U1,UX1,G1,ITYPE)CALL DGDU(MNF,NFREE,TT,X3,IE,U3,UX3,G3,ITYPE)CALL DGDU(MNF,NFREE,TT,X2.IE,U2,UX2,G2,ITYPE)CALL DGDU(NNF,NFREE,TT,X4,IE,U4,UX4,G4,ITYPE)CALL DGDU(MNF,NFREE,TTX5,IE,U5,UX5,GS,ITYPE)CALL DFDU(MNF,NFREE,TT,X1,IE,U1,UX1,F1,ITYPE)CALL DFDU(MNF,NFREE,TT,X3,IE,U3,UX3,F3,ITYPE)CALL DFDU(MNF,NFREE,TT,X2,IE,U2,UX2,F2,ITYPE)CALL DFDU(MNF,NFREE,TT,X4,IE,U4,UX4,F4,ITYPE)CALL DFDU(MNF,NFREE,TT,X5,IE,U6,UXS,F5,ITYPE)CALL FUNCF(MNF,NFREE,TT,X1,IE,U1,UX1,H1,ITYPE)CALL FUNCF(MNF,NFREE,TT,X2,IE,U2,UX2,112,ITYPE)CALL FUNCF(MNF,NFREE,TT,X3,IE,U3,UX3,H3ITYPE)CALL FUNCF(MNF,NFREE,TT,X4,IE,U4,UX4,H4,ITYPE)CALL FUNCF(MNF,NFREE,TT,XS,IE,U5,UXS,HS,ITYPE)CALL FUNCG(MNF,NFREE,TT,X1,IE,U1,UX1,H6,ITYPE)CALL FUJCG(MNF,NFREE,TT,X2,IE,U2,UX2,H7,ITYPE)CALL FUNCG(MNF,NFREE,TT,X3,IEIU3,UX3,H8,ITYPE)CALL FUNCG(MNF,NFREE,TT,X4,IE,U4,UX4,H9,ITYPE)CALL FUNCG(MNF,NFREE,TT,X5,IEU5,UX5,H1O)ITYPE)DO 10 I=1,NN147SI1=SHAPE(I,X(M1),X(M2),X(M3) ,x1)S12=SEAPE(I ,x(M1) ,x(M2) ,x(M3) ,X2)S13=SEAPE(I ,x(M1),X(M2) ,X(M3) ,X3)S14=SEAPE(I,X(M1).X(M2),X(M3),X4)SI5=SHAPE(I,X(M1),X(M2) ,X(M3) ,x5)DSI1=DSHAPE(I,X(M1) ,x(M2) ,x(M3) ,xi)DSI2=DSHAPE(I,X(M1) ,x(M2) ,X(M3) ,X2)DSI3=DSHAPE(I,X(M1) ,X(M2) ,x913) ,X3)DSI4=DSHAPE(I,X(M1) ,X(M2) ,X(M3) ,X4)DSI5=DSHAPE(I,X(M1) )X(M2) x(M3) ,X5)DO 10 3=1,NNS31=SHAPE(3,X041),X(142) ,X(M3) ,X1)SJ2=SEAPE(J,X(M1).,X(M2) ,X(M3) ,X2)S33=SRAPE (3 , X (Ml) X 042) , X (M3 ) , X3)S34=SBAPE(J,X(M1)X(M2) ,X(M3) ,X4)SJ5=SEAPE(J,X(M1),X(M2) ,x(M3) ,x5)DSJ1=DSEAPE(J,X(M1) ,x(M2) ,X(M3) ,xl)DS32=DSHAPE(3,X(M1) ,x(M2),x(M3) ,X2)DSJ3=DSHAPE(J,X(Ml),X(M2) ,x(M3) ,X3)DSJ4=DSHAPE(J,X(M1) ,x(M2) ,x(M3) ,x4)DSJ5=DSHAPE(iX041) ,X(M2) ,X(M3) ,XS)DO 10 K=1,NFREEDO 10 L=1,NFREEI1NFREE*(I-1)+KJ1NFREE*(J—1)+LEMTMP1=(7.ODOO*(G1(1 ,KL)*SJ1+Gl(2,K,L)*DSJ1)*DSI1/+32.0D00*(G2(1,KL)*S32+G2(2,K,L)*DS32)*DSI2/+12.ODOO*(G3(1,K)L)*S33+G3(2,K,L)*DSJ3)*DSI3/+32.ODOO*(G4(1,K,L)*SJ4+G4(2,K)L)*DSJ4)*DS14/+7.ODOO*(G5(1,K,L)*SJS+G5(2,K,L)*DSJ5)*DSI5)—/(7.ODOO*(F1(l ,K,L)*S31+F1(2,K,L)*DSJ1)*SI1/+32.ODOO*(F2(1 ,K,L)*SJ2+F2(2,K,L)*DSJ2)*S12/+12.ODOO*(F3(1 ,K,L)*S33+F3(2,Kj.i*DSJ3)*S13/+32.ODOO*(F4(l,K,L)*S.J4+F4(2,K,L)*DSJ4)*S14/+7.ODOO*(F5(1 jC,L)*SJ5+FS(2,K,L)*DSJ5)*SI5)EM(I1 • J1)=TEMP2*EMTMP110 CONTINUEIF(ITYPE.EQ.0) GO TO 45DO 30 11,NFREE12=NFREE+I13=NFREE+12EM1(I,I)=DTCON(1 4 ,DELX1 ,DELX2,DELX3)EMI(12, 12)DTCON(2,2,DELX1 ,DELX2,DELX3)EM1(13,13)DTCON(3,3,DELX1 ,DELX2,DELX3)EM1(I ,12)=DTCON(1 ,2 ,DELX1 ,DELX2,DELX3)EM1(I , 13)=DTCON(1 3 ,DELX1 ,DELX2,DELX3)EM1(12,13)=DTCON(2,3,DELX1 ,DELX2,DELX3)EM1(12,I)=EMI(I,12)EM1(13,I)=EM1(1,13)EM1(13,12)=EM1(12,13)DO 30 31,NFREE32=NFREE+ 333=NFREE+32148IF(J.EQ.I) GO TO 30EM1(I,J)=0.ODOOEM1(12,J2)=O.ODOOEM1(13,J3)=O.ODOOEM1(I,J2)=0.ODOOEM1(I,J3)=0.ODOOEM1(12,J3)0.ODOOEM1(12,J)=O. ODOOEMI(13,J)=0.ODOOEM1(13, J2)0.ODOO30 CONTINUEDO 40 I=1NEDO 40 .J=1,NECC NOTE: THE EQUATIONS ARE DEFINED IN REVERSE WAY SOC THESE ARE THE RIGHT RELATIONS.CEM2(I,J)=—0.5D00*EM(I,J)—EM1(I,J)*RDELTEM(I,J)=0.SDOO*EM(I,J).-EM1(I ,J)*RDELTCC40 CONTINUE45 CONTINUEDO 20 I=1,NNSIi=SHAPE(I ,x(M1),x(M2) ,x(M3) ,xi)S12=SHAPE(I,X(M1),X(M2) ,x(M3) ,X2)S13=SHAPE(I,X(M1),X(M2).,X(M3),X3)S14=SHAPE(I,X(Mi),X(M2) ,X(M3) ,X4)SI5=SHAPE(I,X041),X(M2) ,x(M3) ,x5)DSI1=DSHAPE(I,X(M1) ,x(M2) ,x(M3) ,x1)DSI2=DSHAPE(I,X(M1) ,X(M2) ,X(M3) ,X2)DSI3=DSHAPE(I,X(M1) ,X(M2) ,x(M3) ,X3)DSI4=DSHAPE(I,X(M1) ,X(M2) ,X(M3) ,X4)DSIS=DSHAPE(I,X(M1) x(M2) ,x(M3) ,X5)DO 20 J=1,NFREEIi =NFREE* (I—i) +JFTMP1=7 .ODOO*(H1(J)*SIi—H6(J)*DSI1)+/32. ODOO*(H2(J)*SI2—E7(J)*DSI2)+/12. ODOO*(H3(J)*S13—H8(J)*DSI3)+/32. ODOO*(H4(J)*SI4—H9(J)*DSI4)+/7 .ODOO*(H5(J)*S15—Hi0(J)*DSI5)F(I1)FTMP1*TEMP220 CONTINUEIF(ITYPE.EQ.0) GO TO 55I1=NFREE* (Mi—i)DO 50 I=1,NEDO 50 J=i,NE12=I1+JF(I)=F(I)+EM2(I ,J)*Y0(I2)50 CONTINUE55 CONTINUERETURNEND149CC CURRENT HANDLES NFREE<=4CSUBROUTINE CURRENT(TO,T,M,NEQ,NPT,NE,NN ,NFREE,/IE,EM,F,X,YO,Y1,ITYPE)IMPLICIT REAL*8 (A—H, O—z)DIMENSION IE(NN) ,EM(NE,NE)DIMENSION F(NE)DIMENSION X(NPT) ,YO(NEQ) ,Y1(NEQ)DIMENSION G1(2,4,4),G2(2,4,4),/G3(2,4,4) ,G4(2,4,4),G5(2,4,4),F1(2,4,4),F2(2,4,4),/F3(2,4,4) ,F4(2,4,4) ,F5(2,4,4),H1(4) ,H2(4),H3(4),/114(4) ,H5(4) ,116(4) ,117(4) ,118(4) ,119(4) ,H1O(4)DIMENSION U1(4),U2(4),U3(4),U4(4),U5(4)DIMENSION UX1(4) ,UX2(4) ,UX3(4) ,Ux4(4) ,UXS(4)DIMENSION EM1(12,12),EM2(12,12),C(12)CC MATRIX DEFINED IN A REVERSED WAYC SO THE MATRIX EQUATION IS (—A)*U=(-B).C WHEN ADDING NEW FEATURES TO THIS SUBROUTINE,C BE AWARE OF THISCMNF=4M1=IE(1)M21E(2)M3=IE(3)TT=O SDOO*(TO+T)X1=X(M1)X5=X(M3)DX=(X5—X1)*O .25D00X2=X1+DXX3=X2+DXX4=X5-DIDELX1=X(M2)-X(M1)DELX2=X(N3)—X(M1)DELX3=X(M3)—X(M2)RDELT=1 ODOO/(T-TO)TEMP2=(DELX1+DELX3)/90. 0000CALL UVAL(MNF,NEQ,NFREE,IE,X1,12,X3,X4,XS,/YO,Y1 ,U1,U2,U3,U4,US)CALL UXVAL(MNF,NEQ,NPT,NFREE,IE,X1,X2,X3,X4,XS,/YO,Y1 ,UX1,UX2,UX3,UX4,UX5)CALL DGDU(MNF,NFREE,TT,X1,IE,U1,UX1,G1,ITYPE)CALL DGDU(MNF,NFREE,TT,X3,IE,U3,UX3,G3,ITYPE)CALL DGDU(MNF,NFREE,TT,X2,IE,U2,UX2,G2,ITYPE)CALL DGDU(MNF,NFREE,TT,X4,IE,U4,UX4,G4,ITYPE)CALL DGDU(MNF,NFREE,TT,X5,IE,US,UX5,G5,ITYPE)CALL DFDU(MNF,NFREE,TT,X1,IE,U1,UX1,F1,ITYPE)CALL DFDU(MNF,NFREE,TT,X3, IE,U3,UX3,F3,ITYPE)CALL DFDU(MNF,NFREE,TT,X2,IE,U2,UX2,F2,ITYPE)CALL DFDU(MNF,NFREE,TT,X4,IE,U4,UX4,F4,ITYPE)CALL DFDU(MNF,NFREE,TT,X5,IE,U5,UX5,F5,ITYPE)CALL FUNCF(MNF,NFREE,TT;X1,IE,U1,UX1,R1,ITYPE)150FUNCF(MNF,NFREE,TT,X2, IE,U2,UX2,R2,ITYPE)FUNCF(MNF,NFREE,TT,X3,IE,U3,UX3,H3,ITYPE)FUNCF(MNF,NFREE,TT,X4,IE,U4,UX4,E4,ITYPE)FUNCF(MNF,NFREE,TT,X5 , IE,U5,UX5,R5, ITYPE)FUNCG(MNF,NFREE,TT,X1 ,IE,U1 ,UX1,E6,ITYPE)FUNCG(MNF,NFREE,TT,12, IE,U2,UX2,R7,ITYPE)CALL FIJNCG(MNF,NFREE,TT,X3,IE,U3,UX3,M8,ITYPE)CALL FUNCG(MNF,NFREE,TT,X4,IE,U4,UX4,U9,ITYPE)CALL FUNCG(MNF,NFREE,TT,XS,IE,U5,UXS,R10,ITYPE)DO 10 I=1,NNSI1=SHAPE(I,X(M1),X012),X(M3),X1)S12=SEAPE(I,X(M1),X(M2),X(M3),X2)S13=SIIAPE(I ,x(M1),x012) ,X(M3),X3)SI4SHAPE(I ,X(M1),x(M2) ,x(M3) ,X4)SI5=SHAPE(I ,x(M1),x(M2) ,X(H3) ,X5)DSI1=DSHAPE(I,X(M1) X(M2) ,X(M3) ,xi)DSI2=DSIIAPE(I.X(M1) ,X(M2) ,x(M3) ,X2)DSI3=DSHAPE(I.X(M1) ,X(M2) x(M3) ,X3)DSI4=DSBAPE(I,X(M1) ,X(M2) ,X(M3) ,X4)DSIS=DSHAPE(I,X011) ,x(M2) ,x(M3) ,X5)DO 10 J=1,NNSJ1=SEAPE(J,X(M1),X(M2) ,X(M3) ,X1)SJ2=SEAPE(J,X(M1) ,x(M2) ,X(M3) ,x2)SJ3=SUAPE(J ,X(M1),x(M2) ,X(M3) ,X3)SJ4=SEAPE(3,X(M1).,X(M2),X(M3) ,X4)SJ5=SHAPE(J,X(M1) ,X(M2) ,x(M3) ,X5)DSJ1=DSHAPE(J,X(M1),X(M2) ,X043) ,x1)DSJ2=DSHAPE(3.X(M1),X(M2),X(M3) ,X2)DSJ3=DSEAPE(3,X(N1) ,X(M2) ,X(M3) ,X3)DSJ4=DSHAPE(J,X(M1) ,x(M2).,x(N3),x4)DSJ5=DSIIAPE(J,X(M1) ,X(M2) ,x(M3) ,Xs)DO 10 K=1,NFREEDO 10 L=1,NFREEI1=NFREE*(I—1)+KJ1=NFREE*(J—1)+LEMTMP1=(7.ODOO*(G1(1 ,K,L)*SJ1+G1(2,K,L)*DSJ1)*DSI1/+32.ODOO*(G2(1 ,K,L)*SJ2+G2(2,K,L)*DSJ2)*DSI2/+12.ODOO*(G3(1,K,L)*SJ3+G3(2,K,L)*DSJ3)*DSI31+32. ODOO*(G4(1 ,K,L)*SJ4+G4(2,K,L)*DSJ4)*DSI4/+7.ODOO*(G5(1 ,K,L)*SJS+G5(2,K,L)*DSJ5)*DSI5)—/(7.ODOO*(F1(1,K,L)*SJ1+F1(2,K,L)*DSJ1)*SI1/+32.ODOO*(F2(1 ,K,L)*S32+F2(2,K,L)*DSJ2)*S12/+12.QDQO*(F3(j ,K,L)*SJ3+F3(2,K,L)*DS33)*S13/+32.ODOO*(F4(1 ,K,L)*SJ4+F4(2,KL)*DS34)*SI4/+7.ODOO*(F5(1,K,L)*SJ5+F6(2,K,L)*DSJ5)*S15)EM(I1 , J1)=TEMP2*EMTMP110 CONTINUEIF(ITYPE.EQ.0) GO TO 46DO 30 I=1,NFREE12=NFREE+II3NFREE+12EM1CI,I)=DTCON(1 ,1 ,DELX1 ,DELX2,DELX3)E141(12, 12)=DTCON(2,2,DELX1 ,DELX2,DELX3)CALLCALLCALLCALLCALLCALL151EM1(13,13)=DTCON(3,3,DELX1 ,DELX2,DELX3)EM1(I ,12)DTCON(1,2,DELX1 ,DELX2,DELX3)EM1(I, 13)=DTCON(1 ,3,DELX1 ,DELX2,DELX3)EM1(12, 13)=DTCON(2,3,DELX1 ,DELX2,DELX3)EM1(12,I)=EM1(I ,12)EM1(I3,I)EM1(I,I3)EM1(13,12)=EM1(12,13)DO 30 J=1,NFREEJ2=NFREE+J33=NFREE+32IF(J.EQ.I) GO TO 30EM1(I ,J)=O.ODOOEM1(12, 32)=0.ODOOEM1(13,J3)=0.ODOOEM1(I,J2)=0.ODOOEM1(I,J3)=O.ODOOEM1(12,33)=0 .ODOOE!41(12,3)0.ODOOEM1(13.J)=0.ODOOEM1(13,32)=0.ODOO30 CONTINUEDO 40 11,NEDO 40 .J=1,NECC NOTE: THE EQUATIONS ARE DEFINED IN REVERSE WAY SOC THESE ARE THE RIGHT RELATIONS.CEM2(I ,J)=—0. 5D00*EM(I,J)—EM1(I ,3)*RDELTEM(I ,J)=0.5D00*EM(I,J)—EM1(I )J)*RDELTCC40 CONTINUE45 CONTINUEDO 20 I=1,NNSI1=SHAPE(I,X(M1),X(M2),X(M3),X1)S12=SHAPE(I ,X(M1),X(M2) ,x(M3) ,X2)S13=SHAPE(I,X(M1) ,x(M2) ,x(M3) ,X3)S14=SEAPE(I,X(M1) ,x(M2) ,x(M3) ,X4)S15=SHAPE(IX(M1),X(M2),X(M3),X5)DSI1=DSHAPE(I,X(M1) ,x(M2) )x043) ,X1)DSI2=DSHAPE(I,X(M1) ,X(M2) ,x(M3) ,X2)DSI3=DSHAPE(I,X(M1) ,x(M2) ,x(M3) ,X3)DSI4=DSHAPE(I,X(M1) x(M2) ,x(M3) ,X4)DSI5=DSHAPE(I,X(M1) ,X(M2) ,X043) ,X5)DO 20 J=1,HFREEI1NFREE*(I—1)+JFTMP1=7 .ODOO*(H1(J)*SI1—H6(J)*DSI1)+/32. ODOO*(E2(J)*S12—H7(J)*DSI2)+/12. ODOO*(H3(J)*SI3—H8(J)*DSI3)+/32. ODOO*(H4(J)*S14-H9(3)*DSI4)+/7 .ODOO*(B5(J)*SI5—H10(J)*DSI5)F(I1)=FTMP1*TEMP220 CONTINUE152IF(ITYPE.EQ.0) GO TO 55I1=NFREE*O.11—1)DO 50 1=1 NEDO 50 J=i,NE12=11+3F(I)F(I)+EM2(I,J)*YO(12)50 CONTINUE55 CONTINUEDO 60 11,NEC (I) =—F (I)DO 60 3=1,NE12=11+3c(I)=c(I)+EM(I,J)*Y1(I2)60 CONTINUEDO 70 I=1NEF(I)C(I)70 CONTINUERETURNENDCCC12=1 1+NFREE13=12+NFREEXX1=0.5D00*(YO(I1)+Y1(I1))XX20 .5D00*(Y0(12)+Y1(12))XX3=0.5D00*(Y0(13)+Y1(13))U1(I)=XX1*SRAPE(1,X1 ,X3,X5,X1)+/XX2*SHAPE(2,X1 ,X3,X5,X1)+/XX3*SHAPE(3,X1 ,X3,X5,X1)U2(I)=XX1*SHAPE(i ,X1 ,X3,X5X2)+/XX2*SRAPE(2,X1 ,X3,X5,X2)+/XX3*SBAPE(3 ,X1 ,X3,X5,X2)U3(I)=XX1*SRAPE(1,X1 ,X3,X5,X3)+/XX2*SBAPE(2,X1 ,X3,X5,X3)+/XX3*SHAPE(3,X1 ,X3 ,X5 ,X3)U4(I)=XX1*SRAPE(i ,X1 ,X3,X54)+/XX2*SHAPE(2,X1 ,X3,X5 ,X4)+/XX3*SRAPE(3,X1,X3,X5,X4)U5(I)=XX1*DSRAPE(i ,X1 ,X3,X5,X5)+/XX2*SEAPE(2,X1,X3,X5,X5)+/XX3*SHAPE(3 lxi ,X3 ,X5 ,X5)10 CONTINUERETURNENDCSUBROUTINE UVAL(MNF,NEQ,NFREE, IE,X1 ,X2,X3,X4,X5,/Y0,Y1 ,U1 ,U2,U3,U4,US)IMPLICIT REAL*8 (A-H, O-Z)DIMENSION IE(3),YO(NEQ),Y1(NEQ),U1(MNF),U2(MNF),/U3(MNF) ,U4(MNF) ,U5(MNF)DO 10 I1,NFREEI1=NFREE*(IE(1)—1)+I153CCSUBROUTINE UXVAL(MNF,NEQ,NPT,NFREE,IE,X1 ,X2,X3,X4,X5,/Y0,Y1 ,UX1 ,UX2 ,UX3 ,UX4,UX5)IMPLICIT REAL*8(A-H , O-z)DIMENSION IE(3),Y0(NEQ),Y1(NEQ),UX1(MNF),UX2(MNF),UX3(MNF)/,UX4(MNF) ,uxs(MNF)DO 30 I=1,NFREEI1=NFREE*(IE(1)—1)+I12=1 1+NFREE13=12+NFREEXX1=O.5D00*(Y0(I1)+Y1(I1))XX2=0. &DO0*(Y0(12)+Y1(I2))XX30 5D00*(Y0(13)+Y1(13))UX1(I)=XX1*DSHAPE(1,X1,X3,XS,X1)+/XX2*DSEAPE(2,X1,X3,X5,X1)+/XX3*DSHAPE(3.X1 ,X3,X5,X1)UX2(I)=XX1*DSBAPE(1 ,X1,X3,X5,X2)+/XX2*DSHAPE(2,X1 ,X3,X5,X2)+/XX3*DSHAPE(3,X1 X3,X5,X2)UX3(I)=XX1*DSEAPE(1 ,X1,x3,X5,X3)+/XX2*DSEAPE(2,X1 ,X3,X5,X3)+/XX3*DSEAPE(3,X1 ,13,X5,X3)UX4(I)=XX1*DSEAPE(1 ,X1,X3,X5,X4)+/XX2*DSHAPE(2,X1 X3,X5,X4)+/XX3*DSHAPE(3,X1 ,13,X5,X4)UX5(I)=XX1*DSHAPE(1 ,X1,X3,X5,XS)+/XX2*DSRAPE(2,X1 ,X3,X5,X5)+/XX3*DSHAPE(3,X1 ,X3,X5,X5)30 CONTINUERETURNENDCCCFUNCTION SRAPE(I,X1,X2,X3)X)IMPLICIT REAL*8 (A—H, O-z)IF(I.EQ.2) GO TO 10IF(I.EQ.3) GO TO 20SRAPE=(X—X2)*(X—X3)/((X1—X2)*(X1—X3))RETURN10 CONTINUESHAPE=(X—X1)*(X--X3)/((X2—X1)*(X2—X3))RETURN20 SHAPE=(X—X1)*(X—X2)/((X3—X1)*(X3—X2))RETURNENDCCCFUNCTION DSHAPE(I,X1,X2,X3,X)IMPLICIT REAL*8 (A-H, O—z)IF(I.EQ.2) GO TO 10154IF(I.EQ.3) GO TO 20DSHAPE=(2. ODOO*X—X2—X3)/((X1—X2)*(X1—X3))RETURN10 DSHAPE=(2 . ODOO*X—X1--X3)/((X2—X1)*(X2—13))RETURN20 DSMAPE=(2. ODOO*X—X1—X2)/((X3—X1)*(X3—X2))ENDCCCFUNCTION DTCON(I,J,DELX1 ,DELX2,DELX3)IMPLICIT REAL*8 (A-H,O-Z)IF(I.EQ.2) GO TO 10IF(I.EQ.3) GO TO 20IF(J.Eq.2) GO TO 30IF(J.EQ.3) GO TO 40DTCON= (DELX2*DELX1*DELI1+2 . ODOO*DELX2*DELX2* (0. 25D00*/DELX3—0. 2D00*DELX2))/(3. ODOO*DELX1*DELX1)RETURN30 CONTINUEDTCON=(0 . 4D00*DELX2-0 . 25D00*(DELX2+DELX3))*DELX2*DELX2*DELX2/1(3. ODOO*DELX1*DELX1*DELX3)RETURN40 CONTINUETEMP =DELX 1*DELX 1 *DELX 1*DELX1TEMP 1 =DELX3*DELX3*DELX3*DELX3DTCON=(0. 25D00*(DELX1—DELX3)*(TEMP—TEMP1)/—0 .4D00* (TEMP1*DELX3+TEMP*DELX1))/(3. ODOO*DELX1*/DELX2*DELX2*DELX3)RETURN10 CONTINUEIF(J.EQ.1) GO TO 30IF(J.EQ.3) GO TO 50DTCON=DELX2*DELX2*DELX2*DELX2*DELX2/(30. ODOO*/DELX 1 *DELX 1 *D ELX3 *DELX3)RETURN50 DTCON=DELX2*DELX2*DELX2* (0. 4D00*OELX2-0 . 25D00*/(DELX2+DELX1) )/(3. ODOO*DELX3*DELX3*DELX1)RETURN20 CONTINUEIF(J.EQ.1) GO TO 40IF(J.EQ.2) GO TO 50DTCON=DELX2* (DELX3*DELX3—2 . ODOO* (0. 2D00*DELX2*DELX2—0. 25D00*/DELX1*DELX2) )/(3. ODOO*DELX3*DELX3)RETURNEND155C.7 List of ASSEMCC ** FOR TRIAGLE ELEMENT**CC NE=NFREE*NNCC ENTRY NEQ NUMBER OF EQUATIONSC NCOLS NUMBER OF COLUMNSC MBAND SEMIBAND WIDTHC NMAX 1ST DIMENSION OF MATRIX A AND VECTOR R,C NMAX>=NEQC MCOL 2ND DIMENSION OF MATRIX A, MCOL>NCOLSC NE DIMENSION OF ELEMENT MATRIX AND VECTORC NN NUMBER OF NODES IN THE ELEMENTC NFREE FREEDOMS OF NODESC IE ELEMENT NODE NUMBERS (I,J,K)C EM ELEMENT MATRIXC F ELEMENT VECTORC A SYSTEM MATRIXC R SYSTEM VECTORC EXIT A ASSEMBLED SYSTEM MATRIXC R ASSEMBLED SYSTEM VECTORCSUBROUTINE ASSEM(NEQ ,NCOLS ,MBAND,NMAX,MCOL,NE,NN,NFREE,/IE,EM,F,AR)IMPLICIT REAL*8 (A—H, O—Z)DIMENSION IE(NN) ,EM(NE,NE)DIMENSION F(NE) ,A(NMAX,MCOL) ,R(NMAX)DO 10 K1,NFREEDO 10 I=1,NNI1=(IE(I)—1)*NFREE+KR(I1)=R(I1)+F(((I—1)*NFREE+K))DO 10 L=1,NFREEDO 10 J=1,NN31((IE(J)—1)*NFREE+L—I1)+MBAND -A(I1,J1)=A(I1,J1)+EM(((I—1)*NFREE+K),((J—1)*NFREE+L))10 CONTINUERETURNEND156C.8 List of GAUSSCCCC ENTRY NEQ NUMBER OF EQUATIONS (RAws)C NCOLS NUMBER OF COLUMNSC MBAND SEMIBANDWIDTRC NMAX 1ST DIMENSION OF MATRIX A, NMAX>=NEQC MCOL 2ND DIMENSION OF MATRIX A, MCOL>=NCOLSC A MATRIX TO BE REDUCEDC EXIT A REDUCED MATRIXC IPUNT ERROR FLAG. IF SET, ERRORCSUBROUTINE FACTOR(NEQ ,NCOLS ,MBAND, NMAX,MCOL,A, IPUNT)IMPLICIT REAL*8 (A—H,O—z)CC REDUCE MATRIX BY GAUSS ELIMINATIONCC A(I,J) HAS SEMI-BANDWIDTH MBAND AND MAY BE ASYMMETRICCDIMENSION A(NMAX,MCOL)IPUNT=0KMIN=MBAND+1DO 50 N=1,NEQIF(A(N,MBAND).EQ.O.ODOO) GO TO 60IF(A(N,MBAND).EQ.1.ODOO) GO TO 20C=1 ODOO/A(N ,MBAND)DO 10 K=KMIN,NCOLSIF(A(N,K).EQ.O.ODOO) GO TO 10A (N ,K) C*A (N ,K)10 CONTINUE20 CONTINUEDO 40 L=2,MBANDJJ=MBAND-L+1I=N+L—1IF(I.GT.NEQ) GO TO 40IF(A(I,JJ).EQ.0.ODOO) GO TO 40KI=MBAND+2-LKF=NCOLS+1-LJ=MBANDDO 30 K=KI,KFJ=J+1IF(A(N,J).EQ.0.ODOO) GO TO 30A(I,K)=A(I,K)—A(I,J3)*A(N,i)30 CONTINUE40 CONTINUE50 CONTINUERETURN60 CONTINUEIPUNT=1WRITE(6,70) N,A(N,MBAND)RETURNC15770 FORMAT(1H1,SX,34H SET OF EQUATIONS MAY BE SINGULAR,//,5X,12SRDIAGONAL TERM OF EQUATION,I5,1311 IS EQUAL TO ,D15.8)ENDCCCCC ENTRY NEQ NUMBER OF EQUATIONS (RAws)C NCOLS NUMBER OF COLUMNSC MBAND SEMIBANDWIDTUC NMAX 1ST DIMENSION OF MATRIX A, NMAX>=NEQC MCDL 2ND DIMENSION OF MATRIX A, MCOL>NCOLSC A REDUCED MATRIXC B LOAD VECTORC EXIT B SOLUTION OF A*X=BC IPUNT ERROR FLAG. IF SET, ERRORCSUBROUTINE SOLVE(NEQ,NCOLS,MBAND ,NMAX,MCOL,A,B,IPUNT)IMPLICIT REAL*8 (A-H,O—z)CC REDUCTION OF A LOAD VECTOR B(I)CDIMENSION A(NMAX,MCOL) ,B(NMAX)IPUNT=ODO 30 H=1,NEQIF(A(N,MBAND).EQ.O.0000) GO TO 60IF(A(N,MBAND).EQ.1.ODOO) GO TO 10B(N)=B(N)/A(N,MBAND)10 CONTINUEDO 20 L=2,MBIJDJJ=MBAND—L+ 1I=N+L—1IF(I.GT.NEQ) GO TO 20V IF(A(I,JJ).EQ.O.ODOO) GO TO 20B(I)=B(I)—A(I,JJ)*B(N)20 CONTINUE30 CONTINUECC BACKSUBSTITUTIONCLL=MBAND+1DO 50 M=1,NEQN=NEQ+1—MDO 40 L=LL,NCOLSIF(A(N,L).EQ.O) GO TO 40K=N+L—MBANDB(N)=B(N)—A(N,L)*B(K)40 CONTINUE50 CONTINUERETURN60 CONTINUEIPUNT= 1WRITE(6 ,70)N ,A(N,MBAND)158RETURNC70 FORMAT(1H1,5X,3111 SET OF EQUATIONS ARE SINGULAR 71.51,12SHDIAGONAL TERM OF EQUATION,I5,1311 IS EQUAL TO ,D15.8)END159Appendix D Listing of Programs Used to Estimate Stationary Domain Parameters for MESFETs#def me#def me#def me#define#de:f me#define#define#defineIP1vseO 3.4esO 1.2nO 5.3e17el 1.66eh 6.5et 3.5} while(l>e);field2=(a0+bO)/2 .0;minlen(field2,nd,&v,&d,&length);prints (“%12. lle’hlO . 2l’/.13 . 3leY.13. OlfI.13 . 41f Y.13 . 41e\n”,D.1 List of GAAS 1/* Program GAAS1 */#include <stdio . h>#include <math.h>3. 1415926541.0e7double dO;double f(double t, double nd);void iuinlen(double field2,double nd,double *v,double *d,double *length);void main(void){jut i;double ahold,bhold, e;double nd,aO,bO,l,tl,t2,f1,f2;double field2,length,d,v;scax(”/.le’/.le’/.le’/.le\n” ,&ahold,&bho].d,&e ,&dO);priutl(”DO’/.le cm2/sec\n”,d0);printf(”Y.12s/.lOs’/.l3sY.l3sY.13s’/.13s\n” ,“Nd cm (—3)” ,“E2 kV/cm”,“V cm/sec” , ‘D cm2/sec” , “Lmin le—6m” , “Nd*L cin(—2)”);for(i0;i<6;i++) {scanf(”Y.le\n” ,&nd);aOa.hold;b0bhold;do {l=bO—aO;t 1=. 382*1+aO;t2= . 618*l+aO;f1=(t1,nd);f2f(t2,nd);if(!1>2) {b0t2;} else {aOtl;160nd,field2,v,d,length,nd*length*1.Oe—4);}}void minlen(double I ield2,double nd,double *v,double *d,double *length){double e;double temp,templ ,temp2,dv,es;efield2;teiupe/eO;temp2temp*t emp*temp;temptemp2*telup;temp2temp2/eO;tentpll.O+temp;esesO*(1 .O+sqrt (nd/no));*v=vs*(e/es+temp)/templ;dv=(vs/es+4. O*teinp2*(vs—(*v)))/templ;*v*0.O1;dv*1.Oe—7;if(e<et) {tempe/el;*ddO*(1 O+temp*temp);}else{tempeh/e;*d=dO*(1.O+pow(temp,2.5));}(*d)*1 .Oe—4;*length(—4.O*1.40273412e—3*dv*(*d)*nd/((*v)*(*v))—1.O)*(*v)*(*v)/(4.O*(*d)*(*d));*lengthPl*1 .0e6*sqrt(1 .O/(*length));(*d)*1 .0e4;}double I(double t,double nd){double e;double temp,teiup1,temp2,dvv,es,d;e=t;tenip=e/eO;temp2=temp*temp*temp;temptemp2*temp;temp2temp2/eO;templl .O+temp;esesO*(1 .O+sqrt (nd/no));v=vs*(e/es+temp)/templ;dv=(vs/es+4 O*temp2*(vs—v))/templ;v*0.O1;dv*1 .Oe—7;iI(e<et) {teinpe/el;ddO*(1 O+teinp*temp);161}J else {tempeh/e;d=dO*(1 .O+pow(temp,2. 5));}d*1 .Oe—4;return((—4.O*1.40273412e—3*dv*d*nd/(v*v)—j.o)*v*v/(4.o*d*d));1* Sample of input data1.1 20.1 1.Oe—5 200.01.0e155. Oe 151. Oe 165.0e161.0e175. Oe 17/* Sample of outputD02 . 000000e+02 cm2/secNd cm(—3) E2 kit/cm v cm/sec D cm2/sec Lmin le—6m Nd*L cm(—2)1.Oe+15 4.74 1.584e+05 641 1.7210 1.7210e+115.Oe+15 5.01 1.492e+05 584 0.6813 3.4066e+111.Oe+16 5.07 1.457e+05 573 0.4869 4.8690e+115.Oe+16 5.21 1.357e+05 548 0.2411 1.2057e+121.Oe+17 5.30 1.301e+05 533 0.1850 1.8497e+125.Oe+17 5.77 1.lSSe+0S 469 0.1147 5.7347e+12162D.2 List of PEAKS/* Program PEAKS *//*#include <stdio . h>#include <inath.h>#define P1 3.141592664#define vs 1.0e7#define eQ 3.4#define esO 1.2*define nO 5.3e17#define el 1.66#define eh 6.5#define et 3.5double dO;double !(double t, double nd);void main(void){mt i;double ahold,bhold,e;double nd,aO,bO,l,tl,t2,fl,!2;double ield2,length,d,v;scanf(”/.le’/.le’/.leY.le\n” ,&ahold,&bhold&e,&dO);print(”DO’/.le cnt2/sec\n” ,dO);printf(”/.12sY.13s\n” ,“Nd cnr(—3)” , “Ep kV/cm”);for(i0;i<6;i++) {sca.nf(”/.le\n” ,&nd);aOahold;bObhold;do{lbO—aO;tl=. 382*l+aO;t2= . 618*l+aO;fl=f(tl,nd);f2=f(t2,nd);if (11>1 2) {b0t2;} else {aOtl;}} while(l>e);field2(aO+bO)/2 .0;priutf(”/.12. lle’/.13 .41f\n” ,nd,field2);}}double f(double t,double nd){163double e;double temp,templ,temp2,dv,v,es,d;et;tempe/eO;temp2=temp*temp*temp;temp=temp2*temp;temp2=tenip2/eO;templl .0+temp;esesO*(1 . 0+sqrt (nd/no));v=vs*(e/es+temp)/templ;return v;}1* Sample of input file1.1 3.1 1.Oe—5 200.01.Oe+15 4.745.Oe+15 5.011.Oe+16 5.075.Oe+16 5.211.Oe+17 5.305.Oe+17 5.77/* Sample of output fileD02 . 000000e+02 cm2/secNd cm(-3) Ep kV/cm1.Oe+15 3.12675.Oe+15 3.16301.Oe+16 3.19095.Oe+16 3.31461.Oe+17 3.41425.Oe+17 3.9032164D.3 List of SLOPE1* Progr SLOPE#IZClUde <stdjo .#include <math.h>efIne#definc#defme#defi5 esO 1.2#defj nO S.3e17#def ins ci 1.66#defj5 eh 6.5et 3.5double dO;double f(doubi t, double fld.doUbl5 ‘12);double velo(double t,double nd);void main(vojd){mt i;double ahold bhold e;double fld,aObO2tl2double field2,lenhdVSCa.nf(H1/.ley.ley.ly.\do)Printf(11y.1231/35\ ,, “Nd (—3)”, “Elfor(i=O;i<6;j÷) {sca.nfQel/.le•/le\,,&‘12=veio(e2nd)aO=ahold•bObhOid;do {lbO—ao;tl.382*i+aOt2.6l*l+aO;f2=f(t2fld)if(fl<f2) {J else {) Vhi15(>).field2(aO+bO), 0;printf(’ø/12 lle133if\fl,,d)Idouble f(doUbl5 t,double Ud,dOUbl5 v2)P1 3.141592654vS 10e7sO 3.4b0t2;aO=tl;165{}double e;double temp,templ.temp2,dv,v,es,d;et;tempe/eO;temp2=temp*temp*teinp;temp=tezup2*temp;temp2=temp2/eO;tenipl=1 .O+temp;eseso*(1 . 0+sqrt (nd/no));vvs*(e/es+temp)/templ;return(Iabs(0.01*v—v2));double velo(double t,double nd){double e;double temp,templ,temp2.dv,v,es,d;et;tempe/eO;temp2=temp*tenip*temp;temp=temp2*temp;temp2=temp2/eO;templl .O+tezup;esesO*(1 .0+sqrt (nd/no));vvs*(e/es+temp)/temp1;return(0.O1*v);}/* Sample input file1.1 3.1l.Oe+155.Oe+151. Oe+ 165. Oe+ 161.Oe+175. Oe+ 171.Oe—6 100.04.484.945.* Sample output fileD02 . 000000e+02Nd ciu(-3)1.Oe+155. Oe+ 151. Oe+ 165.Oe+161. Oe+ 175. Oe+ 17cm2/secEl kV/cm2.0862.0482.0762.2312.3602.935166Appendix E Listing of Program Used to Estimate the Energy Levels in MODFETsE.1 List of MODFETProgram modfet.c#include <stdio.h>#include <string.#include <math.h>/* Constant for MODFET *7#define nc2 4.350e17 7* Effective density of states for electrons *7#define nv2 9.313e18 7* Effective density of states for holes *7#define kt 0.02685 7* k..B*T at 300K (eV) */#define epsilon2 1.142166e—10 /* Dielectric constant for GaAs 12.9*epsilono *7#define eg2 1.425 7* Energy gap in GaAs *7#define dec 0.3 /* barrier high (eV) *7#define eg2kt 55.12572634 7* eg2/kt */#define qepsilon 1.40259e—9 7* q/epsilon2 *7/* Constants for RKF *7#define MAXODE 10#define SCAMIN 0.1#define SCAMAX 10.0#define RMAX 1.Oe-4#define h 1.Oe—2/* Variables for MODFET */double phi2; /* E_c2 — E_f at x+oo *7double nd; /* Dopping density in GaAs*/double eg2phi2kt; /* (eg2—phi2)/kt *7double eta; 7* E_c2(x=0)—E...f *7double etakt; 7* eta/kt *7double quan; 7* Quantum number *7double qhbar, energy;double quano,quanl;double sdent; /* Surface density */double levelsEl000]; /* energy levels *7double poten[500],dpotenCsoo]; 7* potential profile *77* Variables for RKF *7mt n;double t,tf,hl;double y[MAXODE] ,kEMAXODE] [6];double c,refe0;void mnitial(void);void output(double alpha);double target(double alpha);void level(void);double denhole(void);167double denelec(void);double fhole(int n);double fe(double x);double fl(int a);double iutfe(double xl,double x2);/* Functions ±or RKF */void field(double nsdouble alpha);void rkf4(double rmax,double scainin,double scaiuax,double alpha);void suntk(int j,double p,double ql,double q2,double q3,double q4,double q5,double alpha);void evaif (double t ,double y [MAXODE] ,double f [MAXODE] ,double alpha);void main(void){mt i;double ahold,bhold, e;double aO,bO,l,tl,t2,fl,f2;initialO;e1.Oe—3;for(i0;i<31;i+4-) {eta(double)i*O.O1—O.2;etakt=eta/kt;printfQ’nd=%le, Ec—Ef=’/.le\u”,nd,eta);a00.O;b01.O;do {lbO—aO;t 1. 382*l+aO;t2 618*l+aO;fl=target(tl);f2target(t2);if(I1<f2) {b0t2;} else {aOtl;}} while(l>e);printf(1Y.le Y.le\n”,(aO+bO)*O.S,fl);output((aO+bO)*O.5);printfQ’ \n”);}}void initia].(void){mt i;double xl,x2;scanf(”Y.le Y.le\n”,&nd,&eta);phi2kt*log(nv2/nd);168eg2phi2kt = (eg2—phi2) /kt;etakteta/kt;qhbar=163.O*sqrt(O.067);I1* function target gives the difference of the surface densitycalculated from two method */double target(double alpha){inti;double nsl,ns2;1* surface density is in cm{—2) *1/s calculate the surface density from classical model */nsl=2.0e6*epsilon2*kt*(nd*(eg2phi2kt—etakt)+alpha*nc2*intfe(—eg2phi2kt ,—etakt)—nv2*intfe(etakt--eg2kt,eg2phi2kt—eg2kt));usl/1.602e—11; /* scale to cm{—4} */usl=sqrt(nsl);/* solve Poisson’s equation */field(nsl ,alpha);/s get energy levels */levelO;/* calculate surface density */ns2denholeO; 7* surface density of holes 5/sdentdenelecO; /* surface density of electrons */ns2+sdent;7* return value of target *7nslnsl—ns2;return nsl*nsl;}void field(double ns,double alpha){mt i;double delt;x is in le—7 m *//* Ec is in eV *7n=2;yCO]=eta;yEl]ns*qepsilon*1 .Oe—3;169t0.0;pot en [0] =y [0];dpoten[0]=y[1];delt=3.0/499;tf=delt;for(i=0;i<499;i++) {rkf4 (RMAX,SCAMIN,SCAMAX ,alpha);poten[i+1]=y[O];dpoten[i+1]y[1];tf+=delt;}}void level(void)mt j,n,alt,q;double x,s,delta,a,al;double dx,de ,levell ,lastx,laste,lastp,lastq;delta3.0/499;quan=0.0; 1* The first quantum number *1quanoo.0;quanllOO.0; 7* signal that this is the first time call *1for(ni:6 ;n<=498 ;n+2) {energy=poten[n];a1t4;sfl(0)+f 1(n);Ior(j1;j<n;j++) {s+=alt*fl(j);a1t6—alt;}s*delta/3 .0;s—=0.75; 1* s is the quantum number *1if(quanl!100.0) { 1* if not the first time *1quano=quanl;quanl=s;} else { /* the first time *1quanl=s;}if(s>quan) { /* if a level is found *1dx=delta*(quan-quano)/(qua.ni-quano);de(poten[n]—poten[n—1] )*(quan—quano)/(quanl—qua.uO);levell9oten tn—i] +de;x= (n—i) *delta+dx;if(levell<dec) {q=(int)floor(quan);levels [q] leveli;lastedpoten [n—i] +(dpotentn] -dpoten[n—1] )*(quan-quano)/(quanl-quano);lastxx;lastplevell;lastqquan;} else break;170quan+=1 .0;}}quanfloor(lastq); 1* The maximum quantum number *1}void rkf4(double rmax,double scamin,double scamax ,double alpha){mt i,iter;double hinin;double ratio,errest,scale,temp;hlpov(rmax,0.25);if(tf—t<O.0) h1-hl;hmino.5e—4*hl;iter0;do {if(hl*(t+hl—tf)>0.0) h1tf—t;suink (1,0.0,0.0,0.0,0.0,0.0,0.0, alpha);sumk(2,0.25,0.25,0.0,0.0,0.0,0.0,alpha);suink(3,0.375,3.0/32.0,9.0/32.0,0.0,0.0,0.0,alpha);suink(4,12.0/13.0,1932.0/2197.0,—7200.0/2197.0,7296.0/2197.0,0.0,0.0,alpha);suink(5,1.0,439.0/218.0,—8.0,3680.0/513.0,—845.0/4104.0,0.0,alpha);sumk(6,0.5,—8.0/27.0,2.0,—3544.0/2565.0,1859.0/4104.0,—11.0/40.0,alpha);ratio0.0;for(i0;i<n;i++) {errestk[i] [0]/360.0—128.0*k[iJ t2]/4275.0—2197.0*kEi] [3]/75240.0+k [i] [4]/50 .0+2. 0*k[i] [5]/55.0;tempfabs (errest/hi);if (temp>ratio) ratiotemp;}if(ratio<rmax) {t+h1;suink(0,0.0,25.0/216.0,0.0,1408.0/2565.0,2197.0/4104.0,—0.2,alpha);iter++;}scaleo.84*pow(rmax/ratio,0. 25);if(scale<scamin) scalescamin;if (scale>scaznax) scale=scaiuax;hlscale*h;it(ttf) return;} while(fabs(hl)>hmin);printf (“Apparent singularity near t = ‘/.e\n”, (float)t);}void sumk(int j,double p,double ql,double q2,double q3,double q4,double qS,double alpha){mt i;double fsuxn[MAXODE];double sum [MAXODE];171for(i0;i<n;i++) {sujn[i] y Ci] +ql*k Ci] [0] +q2*k Li] [1] +q3*k Ci] [2]+q4*k[i] [3]+qS*k[i] [4];if(j=0) y[i]=sum[i];}if(j0) return;evaJ.f(t+p*hl )sum,fsum,alpha);for(i0;i<n;i++) {k[i] [j—1]=h1*sum[i];}}void eval! (double t ,double y [MAXODE] ,double f [MAXODE] ,double alpha){double ktl;t[O]=yEl];ktll.0/kt;f [1]1 . Oe—8*qepsilon*(nv2*!e((y[0J—eg2)*ktl)—aJ.pha*nc2*fe(—y [0]*ktl)—nd);}double fl(int n){double y;yqhbar*sqrt(energy-poten[n]);return y;}double denhole(void){mt j.,n,aJ.t,q;double x,s,delta,a,al,levell;double dx,de,level,lastx,laste,lastp,lastq;alt4;delta3 .0/499;n498;sfhole(0)+Thole(n);or(j1;j<(n—1);j++) {s+alt*fhole(j);alt6—alt;}s*delta*1 .Oe—5/3.0;returns;}double denelec(void){jut j,n,alt;double x,s,delta,a,al,levell;double dx,de,level,lastx,laste,lastp,lastq;172s0.O;a4. 174e14*O .067*kt;for(n=O;nc=(int)qua.n;n-i-i-) {s+log(1 . O+exp(—levels[nJ/kt));}return s;}double fe(double x){double s;s=3.75994212*pow((x+2.13)+pow(pow(fabs(x—2.13),2.4)+9.6,O.416666666),—1.5);s1.O/(ezp(—x)+s);return s;}double fhole(int n){return (nd—nv2*fe((potenEn]—eg2)/kt));}void output(double alpha){mt i;double x,dx;printf(”Ec—Ef’/.le eV\n”,eta);printf(”Alpha/.le\n” ,a].pha);printf (“Effective density of states: ‘/.le\n”,aJ.pha*nc2);printf (“Surface density’/.le\n” ,sdent);printf (“Potential profile:\n\n”);printf(”x V(x) dV/dx\n”);dx3.O/499;x0.O;for(i0;i<499;i++) {printf(”Y.le ‘/,le ‘/,le\n” ,x,potenEi] ,dpoten[i]);x+dx;}printf(”Number of levels: ‘/.d\n”, (int)quan+1);prmntf(”List of levels:\n\n”);printf(”i Ei\n”);for(iO;i<(mnt)qua.n;i++) {prmntf(”/,d Y.le\n”,i,levels[i]);}}double intfe(double xl,double x2){mt j,n,alt;double x,s,delta;173n=100;xxl;a1t4;delta (x2—xl ) / (double )n;se(x1)+fe(x2);for(j0;j<n—1;j++) {x+delta;s+alt*fe(x);a1t6—alt;}s*delta/3 .0;return s;}174


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items