WHERE‘S PRQE’ERTEES 6F ELECTRGTKERR‘EAL fiEViCES .E‘gmsis {‘09 flu Deqvee cf D“. D. EiliCEiEfiéf‘é SMTE WHERE??? Bipankar Nagchoudhuri I973 0-7 639 This is to certify that the thesis entitled SWIT CHII‘I‘C} PROPERTIES O W ELECTBOTH Eiil'iAL D EV I C35 presented by Dioankar Kagohoudhuri has been accepted towards fulfillment ~ of the requirements for Ph . D 9 degree in E.E. Major professor LIB R a is. 1' MiLhig..iii 1' ’ Universrty r. I i vOVTI ‘7... . ”lectro . , 11": J""CCion In fig ' ”S f; Rik f'evlc i will-3d noglt . y , .i'i’lWI sultc.1] at 4‘. ,8 \"n rnpmt 5 ‘ F l l HT ()f “L111 ‘ In effeCt trit'fil‘at‘lr t5”? SW; ()f new“ '31 in: ABSTRACT SWITCHING PROPERTIES OF ELECTROTHERMAL DEVIC?S By Dipankar Nagchoudhuri Electrothermal semiconductor materials undergo a drastic re- duction in band gap above a critical temperature, Tcr' Two-terminal, bulk devices fabricated from these materials can exhibit current-con- trolled negative resistance I-V characteristics and high-Speed bidirec- tional switching. Switching times of the order of nanoseconds have boon reported, though associated with large storage times, often of the ordvr of hundreds of microseconds. The primary purpose of this dissertation is to study the initi- ation of the switching process in electrothermal devices (ETD'S). A computer-based model is developed incorporating the salient features of the ETD; e.g., the abrupt narrowing of the band gap and rise in car- rier mobility at the critical temperature Tcr' The model consists of a set of three coupled, nonlinear, second-order, partial differential equations. The first of these, the Temperature Equation, is arrived at from thermodynamical considerations of energy balance within the ETD. In effect, the heat energy per unit time associated with the rate of temperature rise at any point in the interior of the ETD is equated to the sum of the net electrical power dissipated and the heat input by the thermal diffusion process. The second, the Continuity Equation, is a generalized form of the continuity of charge equation, which states that «“4 \ ‘4‘ ( :be sun of the r desert and the carrier curreflt of three conpmie presence of an ‘3 the mobile charg the temperature Eaisswfs Equati field to be cons The elec refined by plac :ircdt containi: 135 a switch 8. 53.381, the two I z“.eiSSUITieCi to b {3'3 Surfaces arl 11791131 and Elec' To facu. it ‘a-antized inti :' . .1119 . ._'.._ 46.11.“ '1 ""9 1r‘terval . i 13215. The Terir c Lei‘g n -l ,. 1 r it ‘ 19.1510 T13. “10 “tees 51‘” KM “’ 312.1- Tux-rem c, <57? 'v’i 7 the sum of the rate of increase of charge in any differential volume element and the divergence of the carrier current is zero. Here, the carrier current is assumed to be due to electrons only and comprises of three components, namely, a conduction current component due to the presence of an electric field, and two diffusioi components, one due to the mobile charge carrier concentration gradients and the other due to the temperature gradients within the ETD. The third equation is the Poisson's Equation obtained from Gauss' Law by assuming the electric field to be conservative within the ETD. The electrical boundary conditions for these equations are de- termined by placing the sample between two electrodes in an electrical circuit containing an ideal voltage source Va , a sorrce resistor Rs, pp and a switch 8. The sample geometry chosen is a rectangular parallel- epiped, the two opposite surfaces being in contact with electrodes which are assumed to be ideal heat sinks and electrical conductors. The other four surfaces are in contact with air, which is assumed to be an ideal thermal and electrical insulator. To facilitate the numerical solution of the equation, the sample is quantized into 10 x 10 x 10 identical rectangular parallelepipeds, and finite dirrerence equations were developed for each volume element, thus obtaining 1331 equations for each partial differential equation for each time interval. Different algorithms are used to solve each set of equa- tions. The Temperature Equation, which is parabolicin form, is solved using Douglas' Implicit Alternating Direction (TAD) method in three Space dimensions. The Continuity Equation and the Poisson's Equation are solved using Successive Over Relaxation (SOR) methods. In addition, an overall predictor-corrector loop is employed to achieve simultaneity. Three x '.'0 as the Pr"t 2 tiassing circui ED, and by phC :Eeml charact parallelism to analysis of the 3150 carpare we stsrage times ol is shown to be t “5311? in bot fie“: as shown The model also ; iiii'rinentally i 5s in the into F345 ' ' .m-LOWitChlng ll Three methods of inducing switching in the ETD are studied using V02 as the prototype material--by Joule heating due to the external biassing circuit, by the simulation of a defect within the bulk of the ETD, and by photoswitching. The results of the simulations establish the thermal character of the preswitching region as evidenced by the close parallelism to the results obtained from simplified one-dimensional analysis of the Temperature Equation. The results of the simulations also compare well with the experimental I-V curves and the switching and storage times observed in V02 devices. The mechanism of the switching is shown to be the propagation of the narrowing of the band gap longi- tudinally in both directions in a line parallel to the applied electric field, as shown by the electric field data obtained from the simulationS. The model also provides information regarding the profiles of various eXperimentally inaccessible bulk variables like power density and heat flux in the interior of the ETD. Finally, the possibility of inducing photoswitching under Suitable bias conditions is predicted by the model. ' ""V l . SN} in Department SWITCHING PROPERTIES OF ELECTROTHERPAL DEVICES By Dipankar Nagchoudhuri A THES IS- Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Electrical Engineering and Systems Science 1973 a.” ‘I an I! A cc and encourag few names, h exceptional l littee, notal active and c. and Hellso fc ten of Comm Finally. I Hi typed Up the Press '7 sinc ACKNOWLEDGEMENTS A complete list of all those who contributed in time, effort, and encouragement would resemble, I'm sure, a telephone directory. A few names, however, cannot escape mention: Dr. P. D. Fisher, for his exceptional support and guidance and, of course, members of my com- mittee, notably Drs. Giacolletto, Ho, Nyquist, and Sinna for their very active and considerate co-operation. I also wish to mention Drs. Frame and Wellso for their invaluable help at very crucial times, and the team of computer consultants for '72-73 for their extraordinary patience. Finally, I wish to express my gratitude to the corps of friends who typed up the initial rough draft. To all of these people I wish to ex- press my sincere thanks. ..;k:~'0»L£3: LLSI 0? AP LIST OF 1111 LIST OF FIC CHAPTER: I - 1.1 1.2 1.3 1.4 1.5 CHAPTER II - TABLE OF CONTENTS ACKNONLEEJGBFIBNTS 000.000.000.00...00.00....000000000000000.0.00... LIST OF APPENDICES .00I00000000000000.0000000000.00.00000.00.00.00 LIST OF TJ\BLBS 000.00000000000000...00000....00000000000000.000... LIST OF FIGURES 00000000000.0.000000000.0...0.0.00.0..00000000000. CMPTEli I - INTRODUCTION 00000.00.0000000...0000000000000000.00.0. 1.1 1.2 1.3 1.4 1.5 CHAPTER II - 2.1 2.2. 2.3 2.4 2.5 2.6 2.7 CHAPTER III 3.1 3.2 3.3 Overview ............................................ Electrochermal Device Parameters .................... Switching in Electrothermals .....................o.. Principal ETD Models ................................ The Need for a New Model ............................ THE MODEL ooooooooooooooooooooooo00000000000000.0000. The Band Structure .................................. Donors, Electrons, and the Net Charge Density ....... Carrier Mobility and Conductivity ................... Electric Fields and Potential ....................... Electrical Current Mechanism ........................ Thermal Pr0pertiess Heat Content and Heat Transfer 00000000000000000000000000000000000000.0000. Entropy and Energy Transfer cocoon00000000000000.0000 - MATHEMATICAL FORbiUI-ATIOI‘XS OF THE P10DEL 0.0.0.0...00. Sample Geometry cooon000000000000000000000000.0000... The External Electrical Circuit Configuration ....... The Static D'C Characteristics coco-00000000000000... iii Page ii vi vii viii 10 12 13 14 17 19 23 23 25 26 3.4 3.5 3.6 3.7 3.8 3.9 3.10 CHEER IV - 4.1 4.2 4.3 4.4 4-5 4.6 WEEK v . R 5.1 5.2 5.3 5.5 3.5 3.6 3.7 3.8 3.9 3.10 3.11 3.12 3.13 Energy Storage Mechanisms and Time Constants of the ET.) 0000000000000...oooooeooo00000000000000.0000 The Thermal Boundary Conditions 00000000000000.0000. The Electric Potential at the Electrode - Semicon- dUCtOt Interface osocoo-cooooesooooecoeeoeoooooeooeo Net Charge Density at the Electrode - Semiconductor Interface 000000000000000.0000000000000000...0000000 Boundary Conditions at the Air - Semiconductor Interfaces ......................................... The External Circuit Equation ...................... The Continuity of Charge ........................... The Poisson Equation ............................... The Temperature Equation ooocoooooooo00000000000000. Summary 0000000000000000000000000000.000000000000000 CHAPTER IV - A NIMERICAL APPROACH TO THE SOLUTION OF THE PARTIAL 4.1 4.2 4.3 4.4 4.5 4.6 DIFFERENTIAL EQUATIONS ............................ The Formation of Finite Difference Equations ....... Selection of the Grid Structure .................... The Equation in G .................................. The Equation in 9 .................................. The Equation in 9 ...............................o.. GESCEIC 0f the Logic F10“ ooeoooooooo00000000000000. CHAPTER V ' RESULIS OF THE SIMULATIONS cocoon.cocoa-00000000000000 5.1 5.2 5.3 5.4 5.5 TETM1DSI I‘v curves 000000000000oooooooocooooooeoooo Steady-State Conditions in the Preswitching Region 0000000000000...ooooccoeoooooooeooooeooooocoo Spatial Distribution of Temperatures in the Pre- switching Region ocoeooeeoooeeooeeoe0000000000000... Switching Timea 00.000000000000000accesses...oeoeooo Pilament Formation in the Switching Process ........ iv 29 32 33 35 36 37 39 41 41 44 47 48 51 54 55 57 58 6O 64 65 69 73 75 5.6 5.7 cam VI - 6.1 6.2 6.3 6.4 6.5 BIBLIOG'IAPHY 5.6 Power Dissipation and Energy Flux Density ............ 5.7 Optical Switching of the ETD ......................... CHAPTER V1 - CONCLUSION ........................................... 6.1 Accomplishments of the Research ...................... 6.2 Predictive ASpects of the Model ...................... 6.3 Limitations and Deficiencies of the Model ............ 6.4 A Possible Simplification of the Model ............... 6.5 SUSBQSCiOHS for Future Research cocoooeoeoeceooeeeeooo 79 80 82 83 85 87 89 90 BIBLIOGIAPPIX' 0....0000000.00.00.0000000.00.000.0000.00.0000.0000... 141 Append ix A 5 .')EE I‘v'. THE F NORM Am; mm: LISTI? THE m THE GE Appendix A B LIST OF APPENDICES DERIVATIONS USING THE BOLTZMANN'S EQUATION ...........o. THE FINITE DIFFERENCE EQUATIONS ooooooooooooooeooooocooo NORMALIZATION coco0000000coo00000000000.000000000000000. ALGORITHM FOR SOLNING A SET OF TRIDIAGONAL EQUATION coco TREATMENT OF THE EDGE POINTS BY THE SOR METHOD ooooooooo LISTING .o........oooo..........oo.o........o.o.......o. THE DONOR IONIZATION EQUATION .............a.......o.... TIE GENERALIZED CONTINUITY EQMTION .0...‘.........0.... vi Page 93 101 106 109 111 113 136 139 Table 3.1 4.1 3.1 5.2 5.3 Bark! (kip Time Con Naturaligz ConStant. 5'2b 000 Time Con: Data Plot Cunparim mu 0 I 0 I Table 1.1 3.1 4.1 5.1 5.2 5.3 5.4 LIST OF TABLES Band Gap Data for Various Electrothermal materials ......... Time Constants of the Model ETD Using V02 as Prototype ..... Normalization CORSEBDS Fég F6, Fa, Fa 00000.0.000.0.....0... Constants of the I-V Curves Depicted in Figs. 5.2a and 5.2b 00.00.00.0000000000000.0000.0...0000.00..0....00..0000. Time Constants Computed from F130 503 000.00....0000000oooo0 Data PlOttEd in F1880 5048 and 504D 0.0000000000000000000000 Comparison of Simulated Experimental Data with Analytic mta 000.0.00.0.000000.0.0.000000.0.0.0000.000....0.000000.0 vii Page 31 58 65 66 68 72 Figure 1013 Lib L2 L3 La 2.13 2.1b 2.2 A Typic served Typical of ETD' Hystere: of ETD': Voltage Charactq Current Charactc Energy 1 Ibdfii HI ACtivati the Mode Ionized ibdel £1 Te“Pent Approxin el £1 The 383d The Plan The Sand Samp1e c Figure LIST OF FIGURES 1.1a A Typical Conductivity-vs-Temperature Characteristic 0b- 1.1b 1.2 1.3 1.4 2013 2.1b 2.2 2.3 2.4 3.1 3.2 3.3 3.4 3.5 4.1 5.1 5.2a served in ETD'S 00.0.0000...0000...00.0.0.0000000000000.000 Typical S-shaped Conductivity-vs—Temperature Characteristic 0f 81.0.8 00000000000000.0000...000.000000....0.000000000000 Hysteresis in Conductivity-vs-Temperature Characteristic 0f W's 0.0000000.000.00.000000.00.0000000.000.000.0000.00 Voltage Controlled Negative Resistance Device - I-V Characteristics 00.00.0000.000000.00....00.0000000000.0.0.0 Current Controlled Negative Resistance Device - I-V Characteristics 00.000.000.0000000oa0000.000000000000000... Energy Level Diagram at the Ambient Temperature of the ”Ode-1 ETD 00000000000000.000J0.000......00.000.000.00000... Activation Energy Plotted as a Function of Temperature of the MOCIBI ETD 0....00000000.00000000000....0000.00.00.00... Ionized Donor Density-vs-Temperature Characteristic of “Odel ETD 00.0.0.0....000000.0000000.00......0.000000000000 Tmperature Dependence Of MODLIIEY in MOdE]. ETD 0 0 0 0 0 0 0 0 0 . . Approximate Conductivity-vs-Temperature Characteristic of Model ETD ................................................. The Dead Geometry ......................................... The Planar Geometry ....................................... The Sandwich Structure .................................... Sample Geometry of the Model ETD .......................... Circuit Diagram of the External Circuit ................... Flow Chart ................................................ Typical I-V Curve of ETD .................................. I-V Curve 0f ETD from Simulation 0000.00.00000000.000000000 viii Page 4a 5a 6a 6a 10a 10a 10‘ 10b 13a 13a 23a 23a 23a 25a 25a 59a 61a 64a Sella 5.4b 5.5 5.6 5.7 5.8 5.9 5.10 5.11 50125 5.12b 5.13 Experimental I-V Characteristic of a V02 Device .......... Behavior of Bulk Variables in Low Current Region ......... tr-vs-i in the Absence of Band Gap Switching ............. Final Steady-State Temperatures-vs-Terminal Current ...... Temperature Profiles in x-direction in Low Current Region 0.0000000000000000.00000000000000000000000000.0000. Terminal Current-vs-Time in Switching Region ............. Temperature Profiles Along x Showing initiation of Switching coo-00000000000...cocococoa-0000000000000.oooooo Temperature Profiles Along y oooooeooeooooooooooeoooosocoo Temperature PrOfileB Along Z ocoooocoooooooooooeooooooocoo Spatial Distribution of Current Densities Along y ........ COded Representation 0f E'fields ggoggoogooooooooooocooeos Steady-State Spatial Distribution of Power Density ....... Power Dissipation Density-vs-x During Switching .......... Spatial Distribution 0f Heat Flux Density 0000000000000... ix 64a 65a 67a 67a 69a 74a 75a 75b 75c 78a 78b 79a 79b 79c 80a In the n. :evices was disc ~swigworkers i1: Plagued initiall base devices ex Stance, the suit resistance to lo though large del Fer “99 "ere th. ianufacturing, a CHAPTER I INTRODUCTION In the mid 1960’s, a new class of bulk semiconductor switching devices was discovered, arousing considerable curiosity and excitement among workers in the area of solid state electronics(7). Though . . l 2 plagued initially with problems of repeatability and reliab11ity< ’ ), (3'4'5'6'7). For in- these devices exhibited many desirable properties stance, the switching process was fast; switching times from the high- resistance to low-resistance state were on the order of nanoseconds even though large delay times were always present(4). Other significant pro- perties were the memory and hysteresis effects, their relative case of manufacturing, and their immunity to radiation damage(8). Since investigations indicated that electrical switching was associated with the thermal characteristics of the device material(9'10), these devices will be referred to here generically as "electrothermal devices" or simply ETD's. Initially, this switching phenomena was ob- served in some amorphous transition metal oxides and chalcogenide glasses. Thanks to researchers like Qphen, Eritzsche, and gyshinsky (founders of the C-F-O model)(11'12), Sir Neville Mott(l3'14'15), and Gubanov(l6), rapid progress was made in formulating a tranSport theory for amorphous mate- rials to be analysed using many of the well-established techniques used with conventional crystalline semiconductors. More recently, however, e1ectrotherma1 switching has also been observed in crystalline materials like CdS<17), GaAs<18'19), and Te(20). 1.1. Overview Characteristically, in ETD materials, the band gap narrows dra- matically above a critical temperature, Tcr' In this dissertation, a computer-based model is developed for simulating the tranSport and ter- minal characteristics of such semiconductors. The model is tested against experimental data and then used to predict a photo-switching phenomenon. The first few chapters are devoted to model development. The earlier sections of this first chapter enumerate some of the properties observed in electrothermal switches. Subsequently, theoretical formulations of some of the major workers in the field are discussed in order to estab- lish the need to develop a new model. In Chapter 2, the model itself is constructed. The various electrical and thermal properties are de- scribed in light of the assumptions and idealizations of the model. In Chapter 3, the model is mathematically formulated. Applying the funda- mental laws of electrodynamics and thermodynamics, a set of equations is obtained which describe the electrical and thermal transport proper- ties of the material under certain Specified conditions. In Chapter 4, each differential equation is re-expressed as a set of difference equa- tions, and then the total coupled system of equations is discussed. Chapter 5 deals with investigations of the switching phenomena. Results of various experiments are presented. These experiments were used to test the model against existing experimental data, as well as to predict some as yet unobserved phenomena. An analysis and discussion of these results are also presented in Chapter 5. In Chapter 6, some con- clusions are drawn on the basis of the results obtained, and Suggestions are made for further model development, simulations, and laboratory experiments. 1.2. Electrocher. Lntil abo charge carrier m0 reasonably $10" 3 parameters show a :eristic critical treated as step f'| bepr'marily reSpo Lnthe sections it sf Some important ”2.1. Tenperatu For most a nature, as EXP“ E T 8( ) I El 1: electrothermal certain trans it is: S‘an't' r"*UC material 1.2. Electrothermal Device Parameters Until about a decade ago, semiconductor parameters like band gap, charge carrier mobility, and carrier conductivity were presumed to be reasonably slow and continuous functions of temperature. However, ETD parameters show a very strong temperature dependence near the charac- teristic critical temperature} in fact, some of these parameters can be treated as step functions of temperature. This behavior is thought to be primarily reSponsible for the switching effects observed in ETD's(28). In the sections immediately following, observed temperature dependences of some important ETD parameters are noted. 1.2.1. Temperature Dependence of the Band Gap For most semiconductors, the band gap reduces slowly with tem- perature, as expressed by the following equation: 28(r) - 58(0) - pr (1.2.1) In electrothermal switches, however, the band gap reduces abruptly at a certain transition temperature, Tcr’ where Tcr is characteristic of the Specific material, as illustrated in Table 1.1. Table 1.1 Band GapgData for Various Electrothermal Materials E E g i Critical Temp. T < Tcr T _>_ or Material (K) (V) (V) V02 341 0.45 0.045 V0 126 0.14 0.10 T1203 450 0.04 Metallic v203 150 0.12 0.07 As might be expected, there is often evidence of a structural 22 27 change at the critical temperature( ' )3 for instance, V0 changes from an orthorhombic structure below the transition temperature, 126 K, to a rock salt structure above it; V02 changes from a monoclinic to a rutile structure at the transition temperature, 340 K. 1.2.2. Conductivity-vs-Temperature Characteristics If the proper range of temperatures is chosen, the electrical conductivity is also a strong function of temperature for most semi- conductors. This is because the conductivity is related to the mobile carrier density and carrier mobility<21), both of which are temperature dependent parameters, the relation being: 0‘ - neph + pepp (1.2.2) where 0‘: conductivity of the semiconductor, n - density of mobile electrons, p - density of mobile holes, e a charge on an electron (magnitude only), n - electron mobility which is defined as the carrier drift veloc- ity per unit applied electric field, and p a hole mobility. In electrothermal devices, the conductivity-temperature depen- dence is much stronger than in other semiconductors. For instance, in Futaki's "Critical Temperature Resister"(33), it increased by as much as three orders of magnitude at the critical temperature. Some typical conductivity-versus-temperature curves are shown in Figs. 1.1a and 1.1b. From the Eq. (1.2.2), it is apparent that a sharp increase in conductivity is possible if the mobile carrier concentration increases abruptly. In a doped, partially-ionized semiconductor, the mobile carrier 4a Conductivity Tompmt m Figure 1.1a A Typical Conductivity-vs-Temperature Characteristic Observed in ETD's Ball-ctivity hauratm Figure lalb Typical S-shaped Conductivity-vs-Temperature Characteristic of ETD's density is primarily determined by the ionized donor density. The ionized donor density is determined in part by the activation energy Ea of the ionizing impurity level of the semiconductor, the activation energy Ea for a donor impurity being the energy difference between the donor level and the conduction band edge(6). a-(r) ~Aoexp(-E3(T)/RT) (1.2.3) where A0 is a material constant, k is Boltzmanns constant, T is the temperature, and E3 is the activation energy. The conductivity can also change abruptly due to a sudden jump in carrier mobility. Many ETD's, in fact, do exhibit such jumps. These "mobility gaps“ are an integral part of the C-F-o model<11> referred to earlier, and are often associated with a structural transformation of the electrothermal material(22'27). A local structural change often introduces an increase in the free carrier concentration in the imme- diate vicinity. The screening effect on the neighbor lattice centers increases; in other words, the effective local binding potential reduces, increasing the mobility of more electrons. The effect Spreads rapidly through the material causing an abrupt increase in mobility(l3). Some ETD's also show hysteresis in their conductivity-tempera- ture plots as shown in Fig. 1.2. This feature is utilized in memory 2 2 devices( 8’ 9). The energy associated with the area of the hysteresis (27,28) loop is believed to be related to the energy involved in the structural change, that is, the latent heat of transformation. 1.3. Switching in Electrothermals The large change in conductivity is primarily reSponsible for the negative resistance region in the static I-V curves of ETD's. The 4L |l|i|ll=l=li||l| » tr «out. a u H}St Sa Czuductnity temperature Figure 1 02 Hysteresis in Conductivity-vs-Temperature Characteristic of ETD's I-V characteristics of ETD's differ markedly from junction switches in a couple of reSpects. Junction diode(21) switches are directional; current is allowed in one direction and opposed in the other. On the other hand, electrothermal switches are bidirectional; the switching is essentially a level-sensing and not a direction-sensing mechanism. As soon as a certain temperature level is reached within the device, switching will take place regardless of the direction of current in the external elec- trical circuit. Also, the electrothermal switches are current-controlled, nega- tive-resistance devices. In contrast, typical negative-resistance semi- (21) like the Gunn diode and the tunnel diode exhibit conductor devices voltage-controlled negative resistances. The difference is illustrated in Figs. 1.3 and 1.4. Because of its shape, the I-V characteristics of the ETD are also called S-curves(5). 1.4. Principal ETD Switching Models The S-curves are believed to be caused by "filament" formation, which have been actually observed in some electrothermal devices<43). The "filament" is a thin, highly localized, conductive region extending throughout the length of the material. In filamentary switching, the "off" resistance is determined by the bulk conductivity of the material. In the ”on" state, the filament has already formed. The ”on" resistance is largely determined by the conductivity of the filament. From purely (30) thermodynamic considerations, Ridley demonstrated that filament for- mation was possible in a current-controlled negative differential resis- tance device. The actual mechanism of filament formation was not dealt with. In a later chapter, an equation for energy balance will be derived in close parallelism to Ridley's development. 6a Currant Voltage Voltage Controlled Negative Resistance Device I-V Characteristics p Currant Voltage Figure 1.4 Current Controlled Negative Resistance Device I-V Characteristics A Specific equation for the system energy balance was first determined by Boer and Dohler<31)3 subsequently, Berglund(34), Fritz- 2 sche(1 ), OvShinsky(10) , and others have contributed significantly to its development. In their development, the net electrical energy input at any point is balanced by the temperature rise and the net thermal energy flowing out from the point. The equation forms a sound basis for explaining many of the features observed in electrothermal devices. A computer simulation on the basis of this equation was performed by Warren(35). The addition of the continuity equation to the model was another step forward towards the understanding of these devices. Using the Steady state continuity equation v.3 - 0 (1.4.1.) and the appropriate energy balance equation, Kaplan<36), and others were able to obtain the steady-state curves of the devices with reasonable accuracy. 1.5. The Need for a New Model DeSpite these numerous modeling efforts, many deficiencies existed. Most importantly, all the models previously developed were essentially static in nature. The only time constant involved was in the energy equation. This thermal time constant is of the order of a few hundred microseconds and cannot explain the nanosecond switching observed. Also, several significant features of the electrothermal switches have not been accounted for in models to date. For instance, the ”mo- bility gap“ and the reduction of the "band gap” don't play any part in the formulation of the models. Instead, a change in conductivity, which might include mobility in addition to other bulk parameters, is used as a lumped parameter, explicitly dependent on the temperature. Previously developed models do not take into account the diffusion of charge carriers due to the high thermal and carrier concentration gra- dients that would be inevitably set up when a thermal filament is formed. Secondarily, carrier diffusion would cause extremely high local electric fields to be set up. These effects have been completely neglected in previous models. The present model accounts for the dynamic behavior of the ETD's by incorporating many pertinent features neglected to date by other workers. Both the “mobility gap" and the reduction of the ”band gap" have been in- corporated, and conductivity is treated as a function of both the mobile carrier density and the carrier mobility. A1so, diffusion effects are included by using a more generalized form of the continuity equation. CHAPTER II THE MODEL A dynamic model suitable for observing the high-Speed switching phenomenon is developed here. Due to the complexity of the phenomenon, various simplifying assumptions are made. Each of these assumptions is either immediately or later justified. Care has been taken to retain all of the significant characteristic features of the ETD. An isotropic, nondegenerate, homogeneous, uniformly doped ”n” type electrothermal semi- conductor is considered, the sample being initially in equilibrium at the ambient temperature Tam Being ”n” type, it has a net donor con- b. centration of N atoms/m3 not all of which are ionized at room tem era- d ' P ture. 2.1. The Band Structure For simplicity, consider the semiconductor to be a direct band gap semiconductor possessing a single parabolic valence band and a single parabolic conduction band. 'Furthcr, the semiconductor is assumed to be nondegenerate in the entire temperature range of operation) in other words, though the Fermi-Dirac statistics are applicable, the Boltzmann distribution can be used as a reasonable approximation<21). In a nondegenerate semiconductor, the Fermi level Bf lies within the band gap, and, since the material is "n" type, Bf lies in the upper half of the band gap. Since it is assumed that the material is only partially ionized at the room temperature, the Fermi level if would lie approximately half way between the donor level E and the conduction d 10 13 band edge Ec( ). The energy level diagram under these conditions is shown in Fig. 2.1a, where Si is the intrinsic energy level. In an ETD, the band gap E8 reduces sharply when the critical temperature Tcr is reached. As discussed in Appendix G, the activation energy Ea will also reduce correSpondingly. Thus, the form of the energy level diagram of Fig. 2.1a remains the same after the transition tempera- ture Tcr' Since the energy gap ES and the activation energy Ea change so dramatically at the critical temperature Tcr' the continuous change with temperature is negligable in comparison. So, the activation energy Ea is simulated by a step function with respect to temperature, the step occurring at Tcr° Below and above Tcr' it is a constant. (See Fig. 2.1b.) 2.2. Donors, Electrons, and the Net Charge Density As stated earlier, the donors are assumed to be only partially ionized below the transition temperature. The donor atoms ionize ther- mally, creating one free electron per atom. It is assumed also that thermal ionization of the donors is the only significant process existing within the material which can create mobile charge carriers. If Boltzmann statistics are applicable and if Ea is the thermal energy associated with ionization, then the net density of ionized donors at any point in the material is given by: n:(;,t) - N dexpE-Ea(T)/ZRT(;,t)] (2.2.1) where Nd is the uniform donor density, ngtrat) is the ionized donor den- sity, k is the Boltzmann‘s constant, and T(§,t) is the temperature. This equation is developed in some detail in Appendix G. Since Ea(T) is a step function of temperature only (refer Fig. 2.1b), and Nb is a scalar constant for the material, n: is an explicit function of temperature, as shown in Fig. 2.20 10a -_-----—--------- E' a -.-._.-C-I_O_.-I-E' Ev Figure 2.1a Energy Level Diagram at the Ambient Temperature of the Model ETD '3" bun-v- Activation Emu h----- Tc: Tonuntm Figure 2.1b Activation Energy Plotted as a Function of Temperature of the Model ETD luim Duet Density 10b '7 I I --- .“ ‘L_.rL“”qPE15fi/‘ii' temperature Figure 2.2 Ionized Donor Density-vs-Temperature Characteristic of Model ETD 11 Like the donor atoms, the donor ions are stationary, but each ion carries a charge of +e. The ionized donors are assumed to be the only positive charges within the material; i.e., the density of holes and other positive charge carriers are negligible. Therefore, the total positive charge q+ at any point within the material is: q+(?.t) - enélfit) (2.2.2) Similarly, mobile electrons are assumed to be the only significant nega- tive mobile charge present. The total negative charge, q-, at any point is, therefore, q’('£,t) -en(i~',c) (2,2,3) where n is the free electron concentration. Thus the net charge density 9 at any point (;,t) is: «so . ma.) + q‘('r'.t) - en§(?,c) - end—3c) (2.2.2.) The net charge density if not necessarily zero everywhere within the material, because of temperature gradients and possible diffusion effects, as will be discussed in a later section. The Eqs. (2.2.1)-(2.2.4) are thus assumed to be valid for the entire operating temperature range. Again, since mobile electrons are created only by donor ioniza- tion on a one-to-one basis, the rate of generation of electrons and donor ions are equal: + 396,2) - flat) dt dt (2.2.5) So the net charge creation rate is zero everywhere as follows: Using (2e204), 931;») - gqflat) + ((5,0) dt dt Applying (2.2.2) and (2.2.3) to the above, + d_e_(i=.c) - (23236») - egnfifit) - 0 (2.2.6) dt dt dt 12 But, due to drift and diffusion of electrons, ggfl;,t) will not be zero a everywhere. Instead, the equation of charge continuity applies. (See Appendix H.) 91(33):) - 233?») + V.3'('r‘.t) (it .3: ”Sins (2.2.6) this reduces to a_e_(?,t) +v.3(i~',c) . 0 (2.2.7) at where 31;,t) in the above two equations is the net electron current den- sity. The time of ionization is assumed to be ”instantaneous": this is a reasonable assumption since this is an atomic process and involves times of the order of lO-lsse0537). 2.3. Carrier Mobility and Conductivity Since the material contains only one type of mobile carrier, namely electrons, the carrier mobility refers only to the mobility of the mobile electrons. Also, the mobility is a scalar, since the material is isotropic; and it is a function of temperature, since "mobility gaps” exist in such materials. But, the mobility change at TCr is so drastic that the mobility, like the activation energy, can be taken to be a step function of temperature as depicted in Fig. 2.3. Note that this choice is made to simplify the mathematics of the problem: any given mobility profile could have been chosen in Fig. 2.3. The mobility change is assumed to be an explicit function of tem- perature only. For instance, if the temperature of the material every- where is less than its critical temperature Tor: the mobility is scalar constant having a value Pa’ if the temperature everywhere is above Tor! the mobility is a constant of value Pb' 13 Since only n-type carriers are involved and no minority p-car- riers are present, the conductivity relation of Eq. (1.2.2) reduces to flit) - n(?.t)er(r) (2.3.1) The conductivity is thus not a function of temperature only: however, a fair approximation can be obtained from the following consideration: Though n(r,t) is not an explicit function of temperature, nglr,t) is. Also, at any point (arm; and n cannot be very radically different: a large difference in the two would imply a large local charge separation. This would result in tremendously large electric fields being set up tending to reduce the separation. Thus, for many points within the mate- rial, the conductivity-versus-temperature plot can be obtained by using a-(T)~n:(r)e’v(r) (2.3.2) The product of Figs. 2.2 and 2.3 yields Fig. 2.4 which bears a close resemblance to Futaki's experimental curve (Fig. 1.1a) where0”1, the asymptote to the first section of the curve is Nae/a and 02, the asymptote to the second section is Ndefb. For purposes of the model, however, Eq. (2.3.1) is utilized: otherwise, the effects of the large local electric fields would not be Obs erved s 2.4. Electric Fields and Potential The electric field E'within a material is defined as the force in newtons exerted on a unit charge (coulombs). In rationalized m.k.s. units, which is used consistently throughout, E'is expressed in volts/m. There are two factors that can contribute to the existence of the elec- tric field Bil Firstly, it could be due to the application of an external applied field and secondly, to the creation of local fields by internal charge separation. 13a —da—-- 3" 2" Isliflty --- Tc: human" Figure 2.3 Temperature Dependence of Mobility in Model ETD .fl I I I I I I I I I Conductivity .3 I I I I I I I I i let. Imunim Figure 2.4 Approximate Conductivity-vs-Temperature Characteristic of Model LTD 14 In the present model, two kinds of charged particles are present, namely electrons and singly ionized donor atoms. Although the donor ions are stationary, a net local charge separation may occur due to the move- ment of the mobile electrons. Gauss's Law states: v .E(¥,c) - e(¥.t)/e (2.4.1) where Q is the net charge distribution as defined earlier and E is the permittivity. The only assumption in the above equation is that is a scalar constant. Since the ETD has already been assumed isotropic, the preceding statements imply that is assumed independent of temperature. In addition to the above, an electric field can also result due to a rate of change of magnetic flux in the material. Maxwell's second equation expresses the above statement of Faraday's Law mathematically: V X E(;,t) . - 2:8:(F’t) ' (2.4.2) at (38) The rate of change of magnetic fields is associated with inductive effects; it is large if: 1. the frequency of operation is large. 2. the magnetic permeability Ffof the material is large. 3. the current through the material is large and rapidly vary- ing. 4. the path length of the current through the device is long and strongly coupled to itself, as in a coil. 5. it is coupled strongly to external magnetic fields in its neighborhood. In the system discussed in the next chapter, the dimensions are kept small and the geometry simple. The electrothermal material is non- magnetic and often disordered; so the relative permeability p} is approx- imately unity. Also, there are no strong externally applied magnetic 15 fields in proximity to the device. Even the external circuit is chosen such that the current is limited by an external resistor. Here the mag- netic field effects are neglected and Eq. (2.4.2) can be rewritten as V x E(;,t) .- 0 (2.4.3) The validity of this assumption is discussed further in Chapter 6. From the principles of vector calculus, it is known that if the curl of a vector field is zero, then the vector field is conservative and can be expressed as the gradient of a scalar potential function. The potential function correSponding to the electric field is called the electric potential V. Hence: Rat) - -VV(?.I:) (2.4.4) Combining the above equation with Gauss' Law, Eq. (2.4.1), Poisson's equation is obtained: v.vv('2':',t) = «(EU/e (2.4.5) 2.5. Electrical Current Mechanisms The electrical current mechanisms in a semiconductor can be sub- divided into two broad classes: mechanisms which involve a physical tranSport of mobile charge carriers and those which do not. In the lat- ter category is the capacitative or displacement current. The diSplacement current density (r,t) at any point is ob- diSp tained from the following relation: Jdisp(r,t) - '_B__I:_:_(r,t) (2.5.1) at It is related to the rate of change of electric field and will be large when the electric field varies rapidly, as may happen when switching takes place. The mechanisms which can result in an actual tranSport of the mobile charges, i.e., electrons, are drift and diffusion. The presence 16 of an electric field causes drift. This drift or conduction current density Econd is given by: 3.0.1.152” - «(Roman (2.5.2) where 0(E3c) is the conductivity and E'is the net electric field. Electron diffusion is produced by two methods: due to thermal diffusion caused by temperature gradients in the material and due to concentration diffusion caused by a non-uniform distribution of free electrons in the electrothermal semiconductor. The expression for the total diffusion current 3diff is: 3diff('i~’,c) =dVT(F,c) + eDVn(i~',I:) (2.5.3) and so 3631:) -o(v1*(?.t) + eDVn(?.t) + «(Bali (2.5.4) The above expression can be more rigorously obtained by a direct con- sideration of the distribution function (Appendix A1).¢Kand D are pro- portionality factors. In Appendix A2, 0(and D are evaluated, assuming Maxwell-Boltzmann statistics closely approximate the actual non-equilib- rium carrier distribution function. Hence .51ff(r,t) becomes 3diffG't) - (kle)p(T)[n(?,t)VT(?,t) + T(?,t)Vn(?,t)] which can be rewritten as: Edit-git) - (k/e)/A(T)V[n(?,t)T(F,t)] (2.5.5) where the symbols have the same meanings as already defined. The total current density 3(r,t), comprising the drift and dif- fusion components, is obtained by combining Eqs. (2.5.2) and (2.5.5): 3'5») .. Mac)? + kp(T)v[n(F,t)T(?.t)] (2.5.6) Substituting for f, the conductivity, and E, the electric field, by pre- viously obtained Eqs. (2.4.4) and (2.3.1), the above can be expressed in terms of scalar variables only: 3(2):) = -n(F,t)e/I»(T)vv + kp(T)v[n(?,c)T(F,c)] (2.5.7) 17 . /u(T){kV[n(r-,t)l‘(;,t)] -en(i7,t)VV(?,t)} (2.5.8) In summary, the model assumes that carrier current is attribut- able to electron motion only, which is made up of diffusion and conduc- tion components: 3 a 300m + 3diff (2.5.9) the conduction component being due to the electrical potential gradients in the material 3 com = -n(?,t)e,1('r)vv (2.5.10) and the diffusion resulting from both thermal and concentration gradients A third component, the diSplacement current, is also present, which is contained in the'3g_term of the continuity equation: _ 3:: a1 +V.J . 0 at or 3% +V°Jcond +V. diff .. 0 (2.5.12) 2.6. Thermal Properties: Heat Content and Heat Transfer Since the particle tranSport in the material involved the motion of the electrons only, the mechanisms were predominantly electrical. The thermal parameters play an important role, too, primarily when deal- ing with the energy transport in the material. The two important thermal parameters for this system are the heat capacity c and the thermal con- ductivity kth' The heat capacity c is defined as the amount of heat energy input (joules) required to raise unit volume of the solid (lm3) through 1 K. It is expressed in units of joules/(m3 - K). The increase in heat con- tent increases the lattice vibrational energy as well as the energy of 18 random motion of the mobile carriers--the conduction electrons. But the latter is generally negligible compared to the former unless the lattice vibrational energy is very small: e.g., when the temperature approaches 0 K. At room temperature, the electron contribution can be safely neglected. The vibrational heat capacity is usually a weak func- tion of temperature. Over the limited temperature range used in the model, it can be assumed to be independent of temperature. Also, the model assumes material isotropy and allows only temperature related in- homogeneities. Observe that the Specific heat has about the same numer- ical value for most sclids: it is very weakly dependent on the structure. Thus, if a new structure is achieved beyond Tcr' the Specific heat can be assumed to remain the same. So, the heat capacity of the solid will also remain unaltered provided no change takes place in the overall volume of the solid. This is in fact assumed to be the case in the model under discussion. This assumption is discussed further in the next sec- tion. Therefore, if the temperature of the system at any point in- creases by an amount T, the correSponding increase in internal energy u is Au :- cAT (2.6.1) 0f the three heat transfer mechanisms--radiation, convection, and thermal conduction--only thermal conduction is significant here. Radiation heat transfer is proportional to the fourth power of the dif- ference in temperatures between the two bodies exchanging heat (Stefan's Law)(39). It is insignificant if large temperature differences do not exist between the sample and a neighboring body: in other words, it is negligible if there are no heat sources in the vicinity. Convection implies fluid motion; in a solid sample, there are no convection cur- rents. 19 In thermal conduction, the thermal flux density is related to the temperature gradient by Fourier's Law(39) Eh I -ktHVT (2.6.2) kth is the thermal conductivity of the material and is expressed in units of watts/(m - K). Like a'and 6, it is, in general, a tensor. However, in the model, it is assumed to be a scalar constant. Since the ETD has already been assumed isotropic, the preceding statement implies that kth is assumed to be independent of temperature. 2.7. Entropy and Energy Transfer All material variables can be conveniently divided into two classes<39'46). The first class, called extensive variables, depend on the mass of the material, typical examples being volume V, internal energy U, number of carriers N or P, total charges Q, etc. On the other hand, intensive variables do not depend on the mass of the material. Typical examples are temperature T, chemical potential K, pressure p, and electrical potential V. Similar classifications exist in other fields of study: extensive variables are closely related to "flux" used in physics and engineering or the "through" variables of systems science, whereas intensive variables resemble the "potential“ or ”across” vari- able. A frequently used extensive variable in thermodynamics is "entropy”. It is a measure of the order of the system, and, for a totally closed system, it always increases with time. For ideal, perfectly reversible processes, the entropy for the system remains constant. Each extensive variable in the system is linked to a Specific intensive variable: for instance, the entropy S is usually linked to the temperature T, the electrical potential V with Charge Q, etc. To des- cribe a system of variables completely, a complete set of extensive 20 variables is required, together with the set of correSponding intensive variables. The total internal energy D of the electrothermal sample is a function of all other relevant extensive variables 8, V, N, P, Q. U = U(S, v, N, P, Q) (2.7.1) Taking the total time derivative of the above equation w=m §+a| fl+£ e dt S's'v,N,P,QdI-. 3v S,N,P,th 3N S,V,P,th + au 91; + 33 513 (2.7.2) "5138,V,N,th aQ S,V,N,Pdt The partials in the above expression correspond to the associated intensive variables as follows: flfl/V: N. P. Q a T (the temperature) (2.7.3a) as ED/S, N, P, Q s p (the pressure) (2.7.3b) 3V 32/3, N, P, Q = Kn (chemical potential for electrons) (2.7.30) 3N 5111/8. V. P. Q I Kp (chemical potential for + charges) (2.7.3d) 8P EE/S, V, N, Q a V (electrical potential) (2.7.38) 3Q Eq. (2.7.2) can therefore be rewritten as: (114* 199T1§+Pfl+xnfl+x d+Vgg (20704) dt dt dt dt dt dt The above equation has the general form: 19.. 9" 8P1 dqis dt 1 dt where q.1 is a general extensive or flux variable of the system and pi is the general intensive variable, or the potential reSponsible for the flux qi. In view of the various assumptions made, Eq. (2.7.4) can be sim- plified considerably. The second term is disregarded since the volume of the system i ical work done work is being d dV-O dt I150, the only electrons. Thu V _d_q I . dt Using ti t1 the volume of d_u-'rd dt 5 where u, s, n, n I electrons, and m 535916. The eneré tion is 3180 Deg Parent in a late Eq. (2.). he overall rate :hevarious Gnarg, measure of the W as a 1088 ts zero. The term K Zdu 119 to the Chen ih acreation and ' {Grim GOnOrs . 21 of the system is assumed to remain constant. The term g!_is the mechan- dt ical work done by or on the system: the model assumes that no mechanical work is being done. 9!. 0 (2.7.5) dt Also, the only kind of mobile charge carriers present in the system are electrons. Thus, the term V dQ_can be rewritten as: v g3 - -eV g_N_ dt (2.7.6) dt dt Using the above simplifying assumption and dividing the Eq. (2.7.4) by the volume of the System V, it becomes: + gurgs+ggfl+xpdndwvgfl (2.7.7) (it (it (it (it (it where u, s, n, n; refer to the internal energy, entropy, number of free electrons, and number of ionized positive charges per unit volume of the sample. The energy associated with the chemical reaction causing ioniza- tion is also neglected. The implications of this statement will be ap- parent in a later chapter. Eq. (2.7.7) is an expression of the energy balance of the system. The overall rate of change of internal energy dB is equal to the sum of dt the various energy components of the system. The first term, T g3, is dt a measure of the increase of disorder in the system; it can be looked upon as a loss term. In an idealized, reversible system, this term is zero. The term Kn g2_is the energy associated with the mobile electrons dt n due to the chemical potential Kn‘ It is thus the energy involved in the creation and tranSport of mobile electrons by thermal and concentra- + tion diffusion. Kp dnd is the energy associated with the creation of the dt ionized donors. Lastly, the term -eV dg_is the energy associated with dt 22 the tranSport of electrons due to the presence of the electrical poten- tial V. In other words, it accounts for the drifting of conduction elec- trons due to the electric field it Each of these energies contribute to the overall rate of increase of internal energy of the ETD. CHAPTER III MATHEMATICAL FORMULATIONS OF THE MODEL The previous chapter dealt with the various properties of the electrothermal material and how they were affected by the idealizations and assumptions of the model. In this chapter, a Specific geometry is chosen, and the sample is placed in an electrical circuit, comprising a battery B, a resistor R5, and a switch S. Using the material pro- perties discussed in Chapter 2, a set of equations describing the over- all system is obtained, together with the necessary initial and bound- ary conditions. 3.1. Sample Geometry, Experimenters have commonly used three kinds of geometries in studying the electrothermal switching device<41), namely: a) The Bead or Pellet Configuration; b) The Planar Structure; c) The Sandwich Con- figuration. The three types are illustrated in Figs. 3.1, 3.2, and 3.3 reSpectively. Electrothermal switching was perhaps first observed in sintered metal-oxide complexes<33). The sintering was done under pressure using two half-ellipsoidal molds, with the electrodes laid as rods across the middle and protruding at either end. The resultant beads, however, had one major practical disadvantage-~they had poor thermal contact with the heat sink, which, of course, was the electrode itself. Moreover, the mathematical analysis of such ellipsoidal structures is quite complex. 23 23a Senucswuhuior Figure 301 The Bead Geometry emicomctor ”‘7777 )77/77) //// Figure 302 The Planar Geometry Electrmic and ’ Semiconductor and Heat 630k Heat Sink figure 3.3 The Sandwich Structure 24 The planar geometry, obtained by thin-film deposition techniques and used by Berglund<34), Duchene(40), and numerous other workers<41’42’43'34), proved experimentally far more satisfactory, particularly because of the large contact area with the heat sink. But the structure was asymmetric, and so, the analytic solution of such a geometry was not very simple. The sandwich structure is the simplest to analyze, eSpecially if a Simple geometrical shape, like a cube, cylinder, or rectangular parallelepiped, is selected as the Semiconductor geometry. It also pos- sesses the advantages of the sandwich structure in that there is good thermal contact with the heat sink. The advantages of the various geom- etries have been discussed at some length by Yu<41). Nevertheless, the underlying principles governing the Operation of the ETD remain the same regardless of the geometry chosen. For purposes of this model, the material geometry is selected to be a rectangular parallelepiped with a square cross section. For mathematical analysis, a Cartesian co-ordinate system is chosen with the origin at one corner of the sample,* such that the y and z -axes lie almost along the edges forming the Square face, and the x-axis lies along the other edge. Thus, the sample is located entirely in the posi- tive octant of the co-ordinate system. (See Fig. 3.4.) The semiconductor is assumed to have the dimensions of LX'** Ly, and L2 along the x, y, and z -axes reSpectively. Since the material is assumed to have a square cross section: by - Lz. As pointed out earlier, this choice of a Square cross section has been *Actually, the origin is chosen to lie a distance within the sample where (S > 0 and $<<< I‘x' This is convenient because of the boundary conditions. ”More accurately, L + 28 where J<<< L . This is discussed in detail when discussing the boundary conditions a? the semi-electrode in- terface in Section 3.5. 25 made merely for mathematical convenience: the model would work equally well even if Ly i L2. 3.2. The External Electrical Circuit Configuration The external circuit, connected between nodes 1 and 2 of the sample (refer Fig. 3.4), comprises of the following: a) an ideal d-c voltage source of Vapp, b) a series resistor R8, assumed to be independent of tempera- ture and the current through it, c) an ideal Single pole single throw switch 8 which is initially off and which is turned on at some time t = 0, and d) two highly conducting electrodes at nodes 1 and 2 between which the sample is connected and which is turned on at some time t = 0. Also, the electrodes make ideal ohmic contacts with the sample at each node. Thus, there is no voltage drop either at the electrode or at the interface of the elec- trode and the ETD. The circuit diagram is shown in Fig. 3.5. Node 2 is grounded, and node 1 is the live terminal of the ETD. In addition, the following assumptions are made regarding the system: (a) The electrodes are assumed to be ideal heat Sinks. (b) Air, which surrounds the ETD on four sides, is assumed to be an ideal thermal insulator and dielectric. (c) Initially, before the switch S is turned on, the ETD is in equilibrium. Let the equilibrium temperature be Tamb (d) In comparison to the various circuit parameters, cable capacitances and lead inductances are assumed negligible. (e) The switch 8 is turned on at time t I 0: i.e., all times t are measured with reference to the instant the switch is turned on. (f) In this model, the various variables 25a ETD Material Electrode pl Figure 304 Sample Geometry of the Model ETD I. s , ". M. j. u 1— ' Van-e— m i Figure 3.5 Circuit Diagram of the External Electrical Circuit Used in the Simulations 26 associated with the sample are evaluated at different times, starting from a time Ato after the switch 8 is turned on. The time Ato is chosen large enough such that the initial transients due to the switching im- pulse have died down. But,‘At0 is chosen small compared to the thermal time constant‘ra. (This time constant is discussed in more detail in Section 3.4.) (g) The ETD is assumed to be electrically neutral at all points in the sample before t - O. In other words, n(;,0) - n31?;0) - NdexpE-Ea/kTamb] (3.2.1) 5(E}o) g 0 (3.2.2) 3.3. The Static D-C Characteristics When the switch 8 is turned on at t I O, the entire System is in thermal equilibrium, and the semiconductor is uniformly at the ambient temperature T If db is the conductivity of the sample at the ambient amb' temperature, then its d-c resistance R0 at this time is given by R0 I X , (30301) 2 0’L 0 y where all surface effects are assumed to be negligible. Consider the ETD at a time Ato .AtO is a time large enough such that the initial switching transients due to turning on the switch S have died out. However,1§to is much smaller than the thermal time constant I; associated with the ETD.* At this time, the direct current 10+ at node 1 is given by: + . I0 . Vapp/(RS + R0) (3.3.2) If R8 is chosen to be a fraction f(f _<_o°) of R such that: O + R8 < fRo ohms, then the current IO can be rewritten as: + . . . 10 - Vapp/RO(1 + f) (3.3.3) *This is discussed further in Sections 3.4 and 3.5. 27 Substituting the expression for R0 from Eq. (3.2.1), the current measured at node 1 at t 20+ is 2 10+ 3 Va a—OL . 1 A (3.3.4) 1.,( 1 + f The voltage appearing across the sample at t - 0+ is, therefore, 0 0 Substituting Eq. (3.2.3) in the above equation: + V v0 ”.222. (3.3.5) 1 + f Because the conductivity of the semiconducting sample is finite, power is dissipated due to the current in the semiconductor. The dis- sipated power is converted to heat: therefore, the temperature of the sample will rise. Increased temperature implies increased ionization, and so the conductivity will also increase. The conduction current will increase with time until a final steady state is reached. Hence, the voltage V0+ and current IO+ are the static maximum terminal voltage and minimum terminal current at node 1 of the sample reSpectively and so are designated hereon as vmax and 1min' To obtain the limits of maximum static current Imax and minimum static voltage vmin' it is assumed that the semiconductor is completely and uniformly ionized at some temperature T . Since the material is uni- formly ionized, no thermal or concentration gradients exist, and the only component of current in the sample is conduction current. The conduc- tivity of the sample at this time is (using Eq. 2.3.1) (rmax - Ndepb (3.3.6) Note that in the above expression, the electron concentration is taken to be the donor concentration. This is because the model has only one mechanism for creation of electrons, namely donor ionization. Conse- quently, at any instant, the total number of ionized donors exactly equal 28 the number of mobile electrons. Hence, . Ndefih ~ (3.3.7) Umaxfirb _ NdexpC¥EgalkTamb)e a If the mobility gap ratio is represented by : r T a , max Frexp(Ega/kTamb) (3.3.8) 0— 0 g 6’ Hence: dfiax OPrexP(Ega/kTamb) and so: Rmin 8 —fx—-—. R0 (30309) 0' L 2 Prexfiftga/kramb) max y Thus, the maximum external circuit current Imax is: V = B ( 3.3.10 Imax Vapp/(RS + Rmin) - ipp ( ) ROLf + l. exp(-Ega/kTamb)J 1‘ Similarly, vmin is given by: R vmin I Imamein ' app E x 08 ROD? + .L exI><._§£L)] firexp(_ae_> Pr kTamb amb . vapp (303011) [1 + fiprexp(Ega/kTamb)J Summarizing and using the symbol IO to represent current Vapp/Ro' i.e., IO = Vapp/RO, the maximum limits of the static voltage and current that are allowed by the external circuit configuration are V Vmax . . 322 , (303012) 1 + f V vmin = app , (303.13) 1 + fprexp(h8a/kT) 1min 3 Io ’ and (3.3014) 1 + f I a loprexflfiga IkT) 1 + ffrexp(88;7kT) (3.3.15) 29 The above expressions are the limits of the static terminal cur- rents and voltages. As will be discussed further in the coming section, these can be used as tests indicative of the stability of the results obtained: if the terminal variables were to lie outside the range defined, the ETD would be in an oscillatory state or even a physically unrealiz- able situation. 3.4. Energy Storage Mechanisms and Time Constants of the ETD In Spite of the apparent simplicity of the circuit, its dynamic behavior involves more than just two resistances in series. This is because of the various energy storage mechanisms present in the system, each of which can be associated with a distinct time constant. First, there is the thermal time constant, T3, due to the energy storage associated with the heat capacity of the ETD. If c is the heat capacity, kth the thermal conductivity and LR2 the cross-sectional area of a uniform cube of the material, then‘Pé is given by 2 T I x (3.4.1) kth For the ETD sample chosen, the time constant associated with unit cube of the material is computed below. The relevant material parameters are kth a 6.0 w/(m - x) 3.3 x 105 J/(m3 - K) c a Lk a 10-“ m Ta = 55038 This time constant plays a significant role in the energy balance equation which is discussed in detail in a later section. Since kth and e have been assumed constant, this time constant I; is a constant characteristic of the ETD in this model. 30 A second time constant T; is associated with the thermal and concentration diffusion of the conduction electrons, and is related to the mobility of the electrons. t a LX (3.4.2) Unlike 2;, Tb is a function of temperature and so is not a con- stant for the material. The value of Tb is evaluated at two character- istic temperatures T and Tcr for easy reference: amb a) T a 300 K amb Pa = 4.3 x 10.6 m2/V - s tbl '—' 90.8 ms b) or = 341 x fb B 4.3 X 10-) mZ/V - S tbz = 7088 1118 This time constant is an important factor in the continuity equation that is developed in Section 3.8. A third important time constant is that associated with the per- mittivity E'of the ETD giving rise to bulk capacitance effects. There again the time constant to T I E. (30403) c 5????) is not a constant for the material. Since is a function of the free electron density n and mobility/4,?c is a fairly complex function of (it). Good estimates of 2; can be obtained by using the approximations discussed in Chapter 2. y'av'n:(T)ep(T) (from 2.3.2) Thus for the two characteristic temperatures, the approximate 31 values of 2; are: a) 300 K 4.425 x 10"11 T E F/m 0b = 800 mhos/m I; ' 0055 ps 1 b) *5 I 341 K a 6.15 x 104 mhos/m a"? hi) I 0.00717 ps 2 In addition to these three, there are other time constants which have been neglected in the model. For instance, the time constants due to lead inductances and capacitances have been neglected in the formula- tion of the external circuit equation: and so has the inductive time con- stant due to the time rate of change of magnetic field in the material. The limiting values of the terminal current and voltage can be assumed to be the maximum limits, even when the dynamic characteristics of this model ETD are studied. Thus, the relationships (3.3.12) - (3.3.15) for the limiting terminal current and voltages hold for the entire time Spans of observation: they are used to test that the overall solution obtained at any t is a physically realizable one. The time constants obtained in this section are Summarized in Table 3010 Table 3.1 - 550 s '?h m Eb - 9008 m8 1 'sz - 7.88 ms It ‘ 0.055 p8 °1 'rC ‘ 0.00717 p5 32 3.5. The Thermal Boundary Conditions Looking back on to the circuit diagram of Fig. 3.5, observe that the rectangular semiconducting block, i.e., the ETD is in contact with perfectly conducting electrodes at two of its surfaces x = ~i»and x I Lx + 8, and with air, which is assumed to be a perfect insulator, on the other four Surfaces. In the above statement, the terms ”conduct- ing“ and "insulator” refer to both thermal and electrical conduction. Because the electrodes are ideal heat sinks, the temperature at each of these electrodes will remain at the ambient temperature Tan for all 1b time. Thus, for all time t, T(“8, y, Z, t) I Tam (30501) b TCLX + a, y, z, t) a Tamb (3.5.2) Since 3<< Lx' the temperature at a distance inside the material can be expected to be not significantly different from the temperature at the electrode interface at all times t. Thus T(O, y, z, t) I Tam (3.5.3) b T(Lk’ y, Z, t) . Tam (30504) b Since the surfaces y I O, y I Ly, z I O, and z = Lz are assumed to be in contact with a perfect insulator, there is no heat flowing out of these surfaces(39)3 the normal component of the heat flux density is therefore zero everywhere on these surfaces: qu I 0 at y a 0 and at y a Ly (3.5.5) qu I O at z I 0 and at z I Lz (3.5.6) Hence, applying Fourier's Law (Eq. 2.6.2) kt T I 0 at y I 0 and at y I by 9; kth:T a 0 at z I 0 and at z I L2 2 33 Thermal conductivity being a finite non-zero scalar constant, the boundary conditions become: EEK)“ 0, Z, t) I 0, (3.507) 3y 8T(x' L , Z, t) a O, (3.508) ‘—' y ay 21(x, y, 0, t) 0' and (3.5.9) 32 11$)“ y, L2, t) - 00 (305010) Dz 3.6. The Electric Potential at the Electrode-Semiconductor Interface Consider the ETD at a timeAtO after switch S has been turned on. Recall that timelito is large enough such that the initial tran- sients due to closing the switch have died down; yet, it is small com- pared to the thermal time constant a discussed in the previous section. Ato is the first observation time; for the model, this has been taken to be lO/AS. Since this time is small compared to the thermal time constant, the resistance of the ETD is assumed to be not significantly different from R0, where R0 has been defined as the ETD resistance when it is uni- formly at the ambient temperature Tamb' Hence, the voltage at node 1 at t0 is given by Eq. (3.3.5). Hence, vlmco) - v0+ . vapp/u + f) (3.6.1) Also, since node 2 is grounded, V2(At0) = V2(t) I 0 (3.6.2) Again, since the electrodes are highly conducting and make ohmic contacts with the ETD at the interfaces, the potential at each point on a face in contact with an electrode will have the node voltage corres- ponding to the electrode. Hence: 34 V(-S, y, z, Ato) = vlcoto) (3.6.3) V(l.x + 8, y, z,Atn) I 0 (3.6.4) Obviously, the above is true not only for t IAtO, but for all values of t. V(-S, y, z, t) I V1(t) (3.6.5) V(Lx +'6, y, z, t) = 0 (3.6.6) The distance 8 has been chosen such that it is smaller than the grid Spacing,* but larger than the Debye length associated with the sample. The Debye length for the semiconductor7\u is obtained from the 48 relation ( ) (3.6.7). Afinez For the present simulation, the minimum grid Spacing x and the Debye length )‘D are compared below. 0 . 3.334 x 10.6m at T = 300 K x .. 1 x 10" m A value of 8 of 5 x 10.6m is therefore appropriate. Since 5 is much smaller than the Space dimensions, i.e., 3 << Lx' l.y the potential at points meet either electrode within a S-neighborhood of the electrode interface can be assumed to be equal to the voltage at the electrode. Consequently, the Eqs. (3.6.6) and (3.6.7) for any time become V(O, y, z, t) I V1(t) (3.6.8) V(Lx, y, z, t) I 0 (3.6.9) Specifically, atlAtO, the relations become (using Eq. (3.6.l)): V(0, y, z,A to) :- Vapp/(l + f) (3.6.10) V(Lx' y, 2, Ato) ' 0 (3.6011) *The grid structure is discussed in Chapter 4. 35 V1(t) for all other times t >At0 is evaluated in a subsequent section. Because of the power dissipation in the ETD due to the electri- cal current, the temperature of the ETD will rise. However, this tem- perature rise will not be uniform due to the thermal boundary condi- tions of the ETD. 3.7. Net Charge Density at the Electrode-Semiconductor interface Once the switch 3 is turned on, a surface charge of -¢é and +0§ is built up at the semiconductor surfaces x = -6 and x I LX + 8, the magnitude of which is determined by the bulk capacitance of the ETD. If E is the permigtivity of the semiconductor, vs at Ato is given by é L 0g I [,2 Vlait) coulombs (3.7.1) x However, at a distance 5 within the semiconductor from an elec- trode, there is no surface charge density, since 6 has been assumed to be much larger than the Debye length. The volume charge density Q at such a point can be obtained from a consideration of Gauss' Law: 9 awn? Consider the semiconductor at the instantlbto where tO << 2;, the thermal time constant t0 >> 1;, the time constant associated with the turning-on of the switch. writing the continuity equation at this time: 1+v.[avE]= 0 (3-7-2) at which after a little algebraic manipulation and applying Gauss' law, yields 3 8 '9 + E.%’ 3 0 (30703) .3 a 36 Now Ema-- 30' E + a.- E + 30‘ E where E , E , and E are the electric -- x — 5.— Z X y 2 3x, 3y 2 field components along x, y, and 2 directions reSpectively. But, near the metal electrodes, the tangential components of the E-fields must be zero, because the tangential component of the E-field must be continuous across a boundary, and, in a metal, there can be no electric fields. Hence, By and El are zero and 80, Eq. (3.7.3) becomes inflatifix'o (3.7.4) 3t 6 3x Hence, at the boundaries, can be obtained by assuming Ex to be a con- stant with reSpect to t. The solution for 9 is approximately, -£’_'t 9 2‘90 + Be e Thus, for all t >>5/a’, i’ :9 at x I 0 and x = Lx' For this model, it 0 is assumed that the above statement is true for all times t >tit0. It may be noted that this is not strictly true, since, in the above argu- ment, the material was assumed to be homogeneous. Hence, €(0, Y, Z, t) ' 90 (30705) ?(LX, y, z, t) I 90 (3.7.6) 3.8. Boundary Conditions at the Air-Semiconductor Interfaces Consider now the conditions existing at the four surfaces of the semiconductor exposed to air, for instance, at the plane y = 0. At this boundary, Ex in air is non-zero since the tangential component of the electric field is continuous across the boundary. This is the "fringing" field. The normal component Ey is reSponsible for the normal conduction current, J O’E COUd y = y Since air is assumed to be a perfect insulator, there cannot be any con- duction currents if surface charge effects are neglected. Within the semiconductor,¢7, the electrical conductivity is non-zero; hence, Ey 37 at the surface y I 0 is zero. This is similarly true for all four sur- faces y I O, y = Ly, z I 0, and z I Ly. Using the relation E's -VV developed in Chapter 2, the boundary conditions can be written mathe- matically as: EEOC, 0, z, t) a 0, 3y 2109 L ' Z, t): 0, a), y §!(x, y, 0, t) I O, and 62 31(x, y, Lz’ t) - 0. 82 The volume charge density is related to the electric field by Gauss' Law f I 6V. E. The tangential components Ex and 82 are continuous across the interface y I O, and also, there is no conduction current across the interface. Constructing a pillbox of infinitesimally small thickness A, it can be shown that ‘7. E at y I 0. Hence, by Gauss' Law 9(x, O, z, t) I 0. Similarly, Q(x, Ly' z, t) I 0, 9(x, y, 0, t) I O, and “X. y. Ly. t) - 0- 3.9. The External Circuit Equations So far the terminal voltage V1(t) has only been evaluated at the first observation timeAt0 (refer to Eq. (3.6.1)). The terminal voltage V1(t) and the terminal current i(t) measured at node 1 are functions of 38 time and are related to each other by Kirchhoff’s Voltage Law: V1(t) a va p - i(t) . as (3.9.1) p This is the load line equation. Obviously, this terminal current is equal to the total current at x I 0. The current at x I O is comprised of three components: a) Conduction Current, b) Diffusion Current, and c) Diaplacement Current. The conduction current density in the x I 0 plane is: J (0, y, z, t) I -¢(0, y, z, t)7'V(0, y, z, t) (3.9.2) cond However, at x I O, the temperature T remains the same for all t and, also no significant charge separation takes place here.* Hence, ~ + _ ¢'(0, y, z, t) I n(0, y, z, t)ep(Tamb) M nd(Tamb)eu(Tamh)-<fl) Thus, the Eq. (3.9.2) reduces to Jeond(o, y, Z, t) . "f0 o'V(0, y, z, t) (30903) The diffusion current density at x I O is given by (refer Sec. 2.5): 3diff(0' y, z, t) 1/u(T)kV(n(0, y, z, t).T(O, y, z, t)) Since the temperature at x I O is Tamb for all T5/4(T) is/ga and so Fdiffu)’ y, z, t) I’uak‘WMO, y, z, t).T(0, y, z, t)) (3.9.4) The diSplacement current density Shiff is given by Jdiffu)’ y, z, t) I'eag'gm, y, z, t) (3.9.5) In Sec. 3.8, it was assumed Ex did not change with time.* Hence 3h1ff<°9 Y: z, t) ' 0 (3.9.6) The total current density 3} at x I 0 is therefore: 'J'Tw, y, z, c) . vovwo, y, z, c) 7431. (n(0, y, z,t) 0T(0, y, Z, t)) (3.907) *Surface charge effects are neglected in this model. 39 The total current can be obtained by integrating the above ex- pression (3.9.7) over the face of the semiconductor at the x I 0 plane i(t) .IF 0 d;- (3.908) surface of semi Resolving the carrier current density 3 into three components, Jx Jy, and J2, along the three axes - A A A J J +’J + J . I xx yy 22 Also, dE'I ydy x zdz I xdydz. Hence, 3 . d? . deydz. Thus, the total current i(t) is given by L L i(t) :- gylc‘) sz(0, y, z, t)dydz. (3.9.9) L Y Z _ i(t) I {L .4; (’,t(‘l‘)k[:V(nT):]x ghde. (3.9.10) The above is obtained by substituting Eq. 2.5.8) in the Eq. (3.9.1). The subscripts (x, y, z, t) have been omitted for the sake of brevity. Substituting this expression back into the load line Eq. (3.9.1) I"y Lx v1(t) - vapp - as .4 f0 .k53‘m' a) . dydz (3.9.11)* In subsequent sections, equations for the voltage V and temperature T will be obtained. 3.10. The Continuitygof Chargg The charge continuity equation was obtained for the ETD in Eq. (2.2.7). *Note that {*(T) can be treated as a constant in the above ex- pression, since the semiconductor surface, being in contact with an ideally thermally conducting electrode, is assumed earlier. 40 a_9_ + V. 3' - 0 (3.10.1) at Where 3.15 the current density, the current density 3'is given by Eq. (2.5.10): 3' - «E + kP(T)V(nT) . (3.10.2) Substituting Eq. (2.4.4) and Eq. (2.3.1) in Eq. (3.10.2) 3': kr(T)E7(nT) - neVV] (3.10.3) Computing the divergence of the above expression: V. '3' .- k/u(T)[V2(nT) - neVZV - eVn .vv] + k Vp. [mm - new] (3.10.4) where T - nit). Also, recall that'VZV I -976 (Eq. 2.4.5). Substituting this into Eq. (3.10.4) yields V. 3 I kp(T)[V2(nT) + ne 9/6 - eVn .VV] + [kvlu .V(nT) - kne Vv .vij (3.10.5) Also, recall that in Eq. (2.2.4) ' 9 I en:(T) - en. Rearranging terms in the above yields en . en:(T) - e (3.10.6) Substituting for n in Eq. (3.l0.5) and working the algebra (see Appendix C), a parabolic partial differential equation for the volume charge density is obtained as shown below: 2 1g; +V2?.(_-__kT)+V9. (vv-_2__krr~k__'_r1np)+§ kvl‘ /p&1?))3te e e L e en L: - _9_ + 1 (v kT) 1 k T)- + + T 2+ up. .V - e J + E: V 2d(n efndVV + 1h,(-endvv + kV(ndT))] . 0 (3.10.7) The differential equation obtained requires six boundary con- ditions on each face and an initial condition to be completely solvable, all of which have been arrived at in the previous sections. 41 3.11. The Poisson Equation Once the net charge density 9 is known at every point in the semi- conductor, the electrical potential V(x, y, 2) within the semiconductor can be directly evaluated by using Poisson's Eq. (2.4.5). V.VV(X, y, z) I -S’(x, y, z)/€ (3.11.1) Since €(x, y, 2) has been completely evaluated at each time step in the previous section, the above elliptic partial differential equation is solvable provided a complete set of six boundary conditions are known. Two of the boundary conditions have been obtained by consider- ing the external circuit in Sec. 3.3. The other four boundary condi- tions were arrived at when discussing the air-semiconductor interface in SECQ 3080 3.12. The Temperature Equation The equation for the temperature at any point in the sample is obtained by considering the energy balance within the sample. Rewriting the SQ. (2.7.7) Yields T93 - 93 + Kngg + K d - eVgi_n_ (3.12.1) (it (it dt dt (it (a) (b) (c) (d) (e) Each term of the above expression is considered in turn below: (a) ng -- This represents the rate of creation of entropy in the system. Sgfice the System has been idealized, such that no less terms have been included, this term is zero. (b) du -- Since u, the internal energy per unit volume, is a 2'1? function of Space and time, it can be expressed as: (see Appendix H) d -a(}'t)+V.3"'t 3.12.2 3%; S?’ . u(r, ) ( ) (f) (g) 42 Considering each of the terms separately, 22_represents the ex- 3t plicit rate of change of internal energy per unit volume at each point in the system. This must correSpond to the explicit energy increase due to the rise in temperature of the material. Hence: 136;) = c a; (3.12.3) 3 3 3; represents the energy flux crossing a unit cross-sectional area in the material per sec and is expressed as Jim2 - S). The energy flux in the semiconductor is composed of electrical and thermal, all other forms of energy flow having been assumed absent in the model. Ju 8 J; + Jelec (3.1204) where Ja-and jelec are the thermal and electrical energy per unit area per second, reSpectively. Taking the divergence of both sides, V. Ju V. “'6' + . Jelec (3.12.5) (h) (J) Recalling Fourier's Law as stated in Eq. (2.6.2), (h) reduces to gvzr (3.12.6) Again, the electrical energy flux density, Jelec' is given by V. J; ‘3 “kt Jelec : Jcondv (3.12.7) Computing the divergence of the above: ‘V ='_ '— . 5e1ec Jcond .vv + vv. Jc0nd (3.12.8) and substituting Eqs. (2.4.4) and (2.5.6) yields Vs Jelec ' ' E s E + VVo Jcond (3.12.9) Substituting Eqs. (3.12.9) and (3.12.6) back into Eq. (3.12.5) _ w \7 . Tu - -kthV2T -rb . if + v V. 30”“, (3.12.10) 43 Substituting back Eqs. (3.12.3) and (3.12.10) for terms (f) and (3), the term (b) becomes: dU'Cfi-ktVZT-Eof'i’VVQF h (3.12.11) dt at Thus Eq. (3.2.1) with terms (a) and (b) evaluated now becomes: + 2 fi' -' —‘ dn c §T_- ktfiv T - h . E +~v V. Jcond + Kp d + Kn 92. [ at J dt dt (0)7 (k) (C) (d) dt (e) (a) Consider now the last four terms in the expression (k), (c), (d), and (e). The term (d) represents the energy associated with the creation of electrons at any point in the sample due to the chemical potential Kn, which is reSponsible for thermal and concentration dif- fusion in the semiconductor. Using the now familiar technique outlined in Appendix H and remembering that the electron (particle) has a charge of -e, (d)K anKGn -1V.J. n dt n.at 8 diff ) (3.12.13) Similarly, expanding the term (e) (e) -eV in - -eV(§_r_1_ - _1_V. '3) (3.12.14) dt at e Combining Eq. (3.12.13 and (3.12.14) with (k), we obtain: - ev g2_+ V‘?. T Kn(3.’l ' lV. Jdiff) + VV' Jcond 3t e 3t . . (k) + (d) + (e) (3.12.15) Recalling that the current density 3 due to electron motion is the Sum of diffusion and conduction components, J ' Jcond + Jdiff (2'5'14) Eq. (3.12.15) reduces, after a little algebraic manipulation, to I'll w. .- cu 44 (Kn - ev>g_g_ = (k) + (d) + (e) dt But, (Kn - eV) correSponds to the net electrochemical potential acting on the electrons. Hence, combining this with term (c): + (c) + (a) + (e) + (k) - f: 93 + Kp :11 (3.12.16) dt dt Moreover, recall that the net creation rate of donor ions and electrons is everywhere equal in the model. Hence using Eq. (2.2.5): (e) + (d) + (e) + (k) a (F+K )9}; pdt i1.and Kp are opposite in Sign, and the expression gives the net energy involved in the ionization process. In the present model, this energy of ionization is neglected. The equation for temperature can therefore be written as C 3T 3 kthva +O'E o E (3.12.17) 5'1:- This equation is identical in form to the equations used by Berglund and other investigators(31’34’36). However, the approach used here closely parallels Ridley's thermodynamic argument, thus demonstra- ting the consistency between the two. 3.13. Summary In this chapter, three partial differential equations have been obtained, together with the necessary boundary conditions for T, the temperature, V, the potential, and 9, the charge density. Since all other variables EB' n3, n'f"‘E,’a—’ 3cond’ 3diff can be expressed in terms of these three, the system is completely known if these three equa- tions are solved. The numerical solution of these three partial differ- ential equations by computer simulation forms the subject of the next chapter. 45 For easy reference, the relevant equations together with their boundary conditions are recapitualated below: a) Voltage: ‘VZVCX, y, z, t) = -€(x, y, z, t)/é Boundary Conditions: 22(X, O, z, t) a 0; §!(x, Ly, z, t) a 0 3y By 3V(x, y, 0, t) a 0; 32(x, y, Lz, t) a O z 82 V(O, y, z, t) a V1(t); V(Lx’ y, z, t) a 0 Also, V(t)-V -R fly/L7“ a" 1 " app S 0 0 ak 3&222.' 0 3M)dydz 3)( 3x Att=At0, + V V (At ) = V I 822 1 0 0 1+f b) Charge: 9(Ato) a 0 1 21+vzg.(_-__}_<_l_‘_)+V9.(VV-2_1§_VT-I_<_’I_‘_Vlnf) [4(T3'3t e e e 2 en+ + 3. (iv T+ d-£+VlnP.V(V-_l_<_'_l‘_)) e E G e + (kV2(n+T) - eV +W +v1 (- +Vv + kV( +T))) d nd nr end nd Boundary Conditions: 9(0) Y: z: t) a 90! ?(LR' Y: z: t) 3 €0 9(x, 0, z, t) a O; 9(x, Ly, z, t) a O €(X. y. 0. t) =- 0: 90:, y, L2, t) . 0 c) Temperature: cBT-khvz'r+t'§.§ S? t 46 Boundary Conditions: T(O’ Y! Z, t) a Tamb; T(LX' Y. Z, t) g Tamb 22(x, 0, z, t) a O; §I(x, L , z, t) = O y y y 210‘: Y: O: t) '3 03 910‘: Y: Ly: t) ' O 32 32 CHAPTER IV A NUMERICAL APPROACH TO THE SOLUTION OF PARTIAL DIFFERENTIAL EQUATIONS In the previous chapter, we deve10ped a set of three partial differential equations, each associated with a complete set of boundary conditions. For convenience, these equations are listed again below: the subscripts (x, y, z, t) have been omitted for brevity: a) *vgv + 94. I 0 b) 1 _§_g_ +v Zeb-161‘) + V? . (we) -vz_1_<__'r - LTvlnrflD 17—“ e + +9(:_k v2 'r + end(T) - 1 +V1n’:(T) .Vw - 1(1)) e E 9 e + 1(kVZ(n:T) - e'n:(T)V¢ + V13}: . (kV(n;T) - enSVfl) 8 0 (4.0.2) SE. th The three dependent variables used to describe the system are V, the electric potential; 9 , the charge density: and T, the temperature at, any point in the semiconducting material. If V,q , and T are known for all points (x, y, z) in the semiconductor at all times t, the system is completely described. In the three equations, the only other vari- ables are n; and’I. However, both of these depend explicitly on tem- perature T, which is separately evaluated in Eq. (c). For easy refer- + ence, the expression for (T) and nd [1(1‘) 'Pa -- T < Tar [See Fig. 2.2] are recalled below: 'f'b -- T Z Ter (4.0.4) 47 Ii,“ ' 48 n+(T) I N exp(-E /2kT) - T < T [See Fig. 2.2] d d ga er a Ndexp(-Egb/2kT) - T 3. Tcr (4.0.5) The set of Eqs. (a), (b), and (c) are nonlinear and coupled, and hence are too complex to handle in a closed form analytically. An ap- proach to numerically solving these equations is presented in the re- mainder of the chapter. 4.1. The Formation of Finite Difference Equations The finite difference technique is perhaps the most convenient and frequently-used method to solve a partial differential equation (PDE) together with its associated initial and boundary conditions. In this method, a network of grid points is first established throughout the region of interest occupied by the independent variables. In this case, the independent variables are x, y, z, and t; and the region occupied by the spatial operators is the sample Space which has dimen- sions of Lx' Ly’ and L2 in the x, y, and 2 directions reSpectively. Next, the partial derivatives of the original PDE are approximated by suitable finite-difference expressions involving x, y, and z. The finite difference expressions are obtained by using Taylor's series ex- pansions of the function. The only types of Spatial partial derivatives that occur in our mathematical formulations are the gradient, diver- gence, and Laplacian operators. Thus, in Cartesian co-ordinates, no mixed partial derivatives of the form 32f are involved. The finite are involved. The finite difference exgregsions for these partial de- derivatives are developed in Appendix B. There are two major types of errors associated with this tech- nique, namely the discretization error and the round-off error. 49 The discretization error, e, is the departure of the finite dif- ference approximation from the solution of the PDE at any grid point; i.e., if f1 be the solution of the PDE and f2 the solution of the alge- braic equation, then Gsfl-fz f1 The discretization error plays an important role in the selection of the grid as discussed in the next section. Also, the computational procedure is assumed to be capable of an exact representation of the solution of the finite difference equa- tion. As only a finite number of digits can be retained by the computer, round-off error is introduced. Round-off error is minimized if a) the number of arithmetic operations at any given location is minimized, b) the numbers on which the operations are performed, as well as the result obtained, are of the same order of magnitude, and c) all variables are of the order of unity. The first criterion depends largely on the algorithm used in the solution as well as the grid size involved. The useful technique to reduce errors of the types (b) and (c) is by “normalization". Here, the dependent variable is transformed by suitable additive and scaling factors to a dimensionless quantity lying between zero and one. For instance, consider the temperature equation vz'r -P cE-T-"kth it obtained in the last chapter. In this model, we are interested prima- rily in the preswitching region, where the temperatures within the ma- terial would lie between the temperature Tamb and Tcr' Consider therefore 50 the transformation u I T - Tamb (401.1) Tcr - Tamb Here u is a dimensionless variable lying between zero and one. Using this transformation, the temperature equation becomes: ° 23 ' kthVZ'r ' P (4.1.2) at 'r - _T Ct amb Note that this transformation is useful only in the preswitching region. Once switching takes place, the temperatures in the filament region can be expected to be much larger than the critical temperature Tcr' The other transformations used are as follows: Ax-ExibAy-fiyxziz-Lz’i: 10 10 10 eV (:1 - T ) v - vappgh g- .222}; 1‘ - cr amb a, (4.1.2) L . X where , , and 9 are the normalized dimensionless dependent variables corresponding to V, , and T and x, y, and z are unit vectors in the x, y, and 2 directions reSpectively. The transformed equations thus obtained from (a), (b), and (c) are (a) Davfi?g + Dav§f + Davifi -YD55 +1 - O, (b) 31+ [F6 . y . qvue) - F31 . qVOtG) - Egg . qgvo + Fél . (nflflj :- 0, and (a) '22" PIEDaVAe‘fDaVG'i'DaV 9"? e 9].?F.9 3t 7? x 9 2 D d where 1)::fo . f(x +Ax, y, z) + f(x -Ax, y, z) Dav9f -T[f(x. Y +A)’: 2) + f(x. Y ‘AYo 3)] Davaf -Y[f(x, y, z +Az) + f(x, y, z -Az)] Delfif .- f(x +AxLJ, z)- f(x -Ax, y, z) 2 51 Delyf a J7f(x, y +Ay, z) - f(x, y -Ay, z), 2 Delzf = ‘5’“be 2 +112) - f(x, y. 2 ~02). 2 IT" _£L_: kth I . n;(T) NdexpCEgbiTcr) - Ndexp( 8é7Tém6) ‘1’. a n(T) ' NdexpCEgb7T;;) - Ndexp(E 7T ) ga amb a w .. n * [(Delx¢)2 + (mnyebz + (Delz¢)2], and WC ai§0 = ;(Tamb)' These transformations are formulated using n1 2 n2 I n3 = 10 in Appen- dix C. The choice of n1, n2, n3 is discussed in the next section. In the equations and relations listed above, the spatial subscripts have been omitted; i.e.,} , 11,, , (,1' , 6, fl, and xare all functions of space and Davx, Delx, Davy, Dely, Davz, Delz denote Special Spatial operators. (The terms Fa, Pb, PC, Fd'Y'TD' AT, and Ar are all dimen- sionless constants. 4.2. Selection of the Grid Structure In the last section, it was pointed out that the round-off error could be minimized by a suitable choice of algorithm and by normalizing the variables to dimensionless quantities of the order of unity. But this has little or no effect on the discretization error. The discretization error depends largely on the nature of the grid structure chosen. For computer applications, it is preferable to use as uniform and regular a grid structure as possible, since this re- duces the memory storage requirements and the number of operations on the computer. Thus, this entails not only a saving in terms of computer 52 time and money, it also minimizes the round-off error, which, as was pointed out earlier, depends on the number of operations performed. Hence, a cubic or a rectangular grid are the obvious choices. Recall that a rectangular parallelepiped with a Square cross section had been chosen for the semiconducting sample.* To minimize round-off error, an equal number of points are chosen along each axis. Obviously, the mesh Spacing in the three directions is not the same; in fact, it is proportional to the dimensions of the sample in the three directions. In other words, if Li, Ly' Lz are the sample dimensions and n1, n2, :13 are the number of points along the x, y, 2 directions reSpectively, then the grid Spacings Ax, Ay,.Az are given by L L Ax . E35; by = _y_; A2 a _g (4.2.1) n1 n2 n3 Since L y ' L2 and n1 - n2 8 n3 2 n x . Egg Y a El: 2 = 2 (4.202) n n n The number n is selected from the following considerations. First, it is desirable to minimize the discretization error. Since all the differential equations are of second under, consider the discretiza- tion error involved in obtaining 32f. From Appendix B, for an interior 3x2 point (x0, yo, 20) *Another logical choice, from symmetry considerations, would have been a cylindrical sample. However, this choice leads to program- ming complications because of the singularity of the Laplaciar opera- tor at r I 00 S3 82f . f("o 1’ 24’“ ’0' 20) ' 2“"0 +Ax' yo' 7‘0) 1' f("0' y0' zo) "" 2 N 3x Ax 2 3 + X a f + as. 21 3x3 The discretization error e is hence given by 33f 2 3 e m x 73x 2: a?- 8x? Thus, e is minimum when Ax ¢ 0, and hence when n 4'00. However, there is a lower limit ofiix determined by the nature of the problem. While discussing the boundary conditions for the charge distribution, it was required that the Spacing be much larger than the Debye lengthTKD to be able to neglect surface charge effects. Thus x >DVAD. Hence L LX << l..- E? D n- Since Am and D I 3.34 x 10-6m, n << 30. LXI-1x10- However, practical considerations limit the choice of n even further. Since a number of arrays (at least one for each dependent vari- able) of size n x n x n have to be simulataneously stored in the computer memory, the core memory requirements and costs increase very rapidly with n. A compromise choice of n was 10, which required an operating computer core memory of 100 k.** In Appendix C, the normalizations are performed with n I n1 I :12 I 10. **At the time, the maximum allowable core memory on the CDS 6500 computer used was 120 k. 54 4.3. The Equation in j For the equation in 15, X-is assumed to have been already eval-\ uated, and so it is a second order elliptic partial differential equa- tion, in three Space dimensions with'X as a forcing function. Typically, elliptic PDE'S are solved by some kind of relaxation method: in this case the Successive Over Relaxation (SOR)(44) method is used to Speed up the converging process. Briefly, this method consists of initializing each element of the matrix to some convenient value. In this Specific simu- lation, the value obtained at the previous time step is used. A new approximation to the value of an element is obtained using the governing Eq. (b) and the relaxation parameter w, and the current value immediately replaces the previous value of the element in the matrix. For instance, for an interior point, the equation for f is: Davx¢ + Davy¢ + Davzfi - TD” +x I 0 (4.3.1). Expanding the above by using the definitions of Dav , Davy, and Davz x ”(x +Ax , y, z) + ¢(x -Ax, y, z) + 0(x, y +Ay, z) + i(x, y -Ay, z) + “x, y, z +Az) + (Xx, y, z -Az) - {D¢(x9 y, Z) +$(x9 Y: Z) s o (4.302) Rearranging terms: ”(X, y, Z) ' L ”(X +AX, Y: z) + ”(X “AX, Y. Z) 4' ”(X9 Y +AY: 2) TD + C(x, y -Ay, z) + 0(x, y, 2 +432) 4» ¢(x, y, z -Az) +f(x, y, 2) (4.3.3) In the 80R method, a new value ¢+(x, y, z) is obtained by ¢+(x. y. z) I ¢(x, y, z) +'w * [0'(x. y. 2) - 0(x, y, 2)] (4.3.4) where I'(x, y, 2) denotes the right hand side of Eq. (4.3.3). Similar expressions for ¢'(x, y, 2) can be obtained for points lying on the edge as shown in Appendix E. 55 The relaxation parameter w determines, to a large extent, the rate of convergence. If chosen too small, the solution tends to converge slowly: if too large, the successive inerations may not converge at all. (45) Forsythe and Hasow recommends an optimal choice of w of 1.5. Convergence is determined by the relative deviation RD: RD 3 Elf}, (X: Y: z) ' “()hj: 2H (4.3.5) over 221 “x, y, 2)] all elements when RD drops below a Specified value Eps' convergence is attained, and the solution is acceptable. As before 6 and 0+ stand for the older and newer approximations to the solution reSpectively. For the complete pro- gram listing, refer to Appendix F. 4.4. The Equation is 9 The normalized equation for 9 is: as - 1 [Davifi + Davie + Davfie Cine] a PFd This is a parabolic partial differential equation in three Spa- cial dimensions. However, the equation is non-linear, the non-linearity arising from the nature of the forcing function P-FU' From a knowledge of the conductivity and the electric field E, the forcing function P is first evaluated for each time interval: the differential equation is then solved assuming the forcing function to .be invariant during the evalua- tion process. The boundary conditions on the six boundaries, as ennumerated earlier, are of the von Neumann type<45). In normalized form, these are 9 I 90 at x I O and at x I 1.x )9 I 0 at y I O and at x I L s;’ Y 22': O at z = 0 and at z I L2 32 56 Two standard methods exist for solving this kind of equation; i.e., explicit methods and implicit methods. In explicit methods, to obtain a solution at time tn + 1, all Space derivatives are evaluated at the instant tn. Though this implies that the calculation procedures are inherently simpler, it requires that very small time steps be taken relative to the Space grid size; i.e. 0 < t/ x2 5_% This is undesirable particularly because thermal time constants are usually large and hence would place severe restrictions on the size of the sample and the Space grid structure. Implicit methods overcome this difficulty at the expense of a Somewhat more complicated algorithm. It consists of representing the Spatial derivatives 5x2, Syz, etc. by a finite difference form evaluated at the advanced point of time tn + 1, instead of tn as in the explicit method. There are, however, two serious drawbacks to the standard im- plicit method in three dimensions. For the one dimensional case, the scheme is Stable and is independent of A, where A.= At/sz But, in three Space dimensions, the method is unstable for l.> 3/2(45). Furthermore, it requires an inversion of a fairly complex matrix, which takes considerable computation time. The implicit alternating direction method, developed by Peaceman and Rachford and Douglas(aa'45'49) avoids these disadvantages and manages to use a system of equations with a tridiagonal coefficient ma- trix which can be solved by fairly straightforward methods (see Appendix D). Essentially, the principle is to employ three difference equations for each time step, the solution of the first two being intermediate 57 values denoted by u* and u** as follows: (2"); )(u.'. .. ”n ) . _1_ 5x201" + on) +15y2un + Yszzun + FDP (4.4.1) 2 2 2 - I: ‘4‘ - 'k' Ax (u** on) 1 5x201 on) +1; 8y (on + u w) + 52 u zit 2 + F P (4.4.2) D Ax2(u -u)-lé 2(u*+u)+f$2(u**+u) E n + 1 n E X 1'1 ‘2" y n +r5 2(u + u ) + F 1’ (4.4.3) 2 n + l n D In the above difference equations, the time subscripts are n and n + l. The Space-subscripts have been omitted for clarity. Leaving Eq. (4.4.1) as it is, subtracting Eq. (4.4.2) from Eq. (4.4.1) and Eq. (4.4.3) from Eq. (4.4.2) and rearranging terms, we obtain 1‘cu""'--lSx 2__u*Ilun+ xzun +75}, 2u +7822 1: +FDP 2 2 n n Kuirw-TsyzuflIXu*-152u 2 Y n 2 K%+l';éyzun+IBXU**-% zun In the above set, the right side is a known column matrix, whereas the left side is a tridiagonal matrix. The algorithm for solving such a system is shown in Appendix D. The actual implementation is done in subroutines TEMPEQN and TRIDAG (see listings in Appendix F). 4.5. The Equation in 3 The equation for 9 in normalized form has been written out in Section 4.2. This also is a parabolic partial differential equation in three Space dimensions, but with Dirichlet boundary conditions at the boundaries. However, this equation is considerably more non-linear than the temperature equation, and 50, though the alternating-direction scheme deve10ped there could be applied, considerable computational time is involved in merely evaluating the tridiagonal coefficients. whim 58 Instead, an iteration scheme based on the 80R method discussed earlier is used. This method is particularly advantageous because of the numerical values of the constants Fé, F5, Fé, and F5. Typical values of these constants are illustratedin Table 4.1 for a V02 sample. Table 4.1 Normalization Constants F5, F6, F8, Pd F; 7.96 x 107 mz/s Fb 1.63 x 105 mz/s F; 6.58 x 10‘2m2/s Fa 1.29 x io'zmz/s Since F5 and F5 are so much larger than either Fé or F5, and all quantities:}$ 6, 9, etc. have been normalized, the terms containing the coefficients Pd and Fé can be neglected to a first approximation. Thus, the equation in 9 reduces to xT'xl'l -J(n - 1 F; Vs (Q‘V‘)- Ft.) Vs (Qv(§9)) (“0501) [it and hence n x g xn _ 1 + 5% [F5 v. (qivm - Ff) v. (qv(§e))] (4.5.2) The right side of the equation is known, and so the first ap- proximation to xn to start the iterative process is obtained from Eq. (4.5.2). Because of the dominance of this term, the iterative process converges very rapidly. As in the Poisson equation, a relaxation para- meter w or 1.5 is used. 4.6. Gestalt of the Logic Flow In the preceding sections, we have developed schemes to solve each of the three differential equation--in 9, ¢, and? --with the ap- propriate boundary conditions. But, the three equations are mutually it. .w 59 coupled to each other; therefore, solving each of the variables 9, fl, and Q in turn assuming the other two known does not imply simultaneous solution. To correct for this error, a predictor-corrector loop is employed, its logical flow chart being illustrated in Fig. 4.1. The method uti- lizes the feature that the thermal time constant is long; in other words, of the three dependent variables 9, fl, and? , the temperature would vary least within a time step. As the flowchart indicates, an estimate of the terminal voltageyi: t is first made by solving the n + l circuit equation using the values obtained at tn. This is then used as a boundary condition to solve for the matrix ¢ for all elements of the Space grid. Using this value ¢* and the temperature at tn, the conti- nuity equation is Solved to obtain *. The electron density n is then computed, and, applying Simpson's Rule, the new estimate of terminal current i§(t) is obtained. The circuit equation is used again to obtain a new estimate of the terminal voltage fly. The new estimate is tested against the older estimates; the looping is terminated if |¢***;*-¢:*i (EL where éL is a preset positive real number. When the above convergence test is satisfied, the current estimates of 9 and D are taken to be the true values at time t These are then used to solve the tempera- n + 1' ture equation for time tn + 1. The time is then stepped Up and the pro- cess is repeated for the next time interval. EVALUATE NORMALIZING CONSTANTS INITIALIZE ALL VARIABLES my». (we l n30 STCRE I 0 59a SET I, V TO LIMIT VALUES I Figure 4.1 Flow Chart l SOLNE CONTINUITY EQUATION FOR CHARGE [STORE I ITER SOLVE POISSONIS EQUATION FOR V IS ITER.LT. N0 SOLN NOT ::> (couvsacmc . 310? s a”. P SOLVE TEMPEQN t.“ INNER LOOP (2)45) ggran LOOP CHAPTER V RESULTS OF THE SIMULATIONS Now that the model has been translated into a computer program, it is possible to demonstrate its utility. This is best achieved by simulating the behavior of the ETD's under various biasing and other external conditions to verify its "goodness" and then making some pre- dictions regarding the device based on the outcomes of the simulations. Of course, the ”goodness" of a model is relative. It has meaning only in relation to the objectives of the model. In other words, while there are various aSpects and facets to an experiment, a model seeks to ex- plain and make understood those aSpects which are of primary interest to the investigator. Thus, the model has to fit the objectives and goals set by the investigator, and the "goodness” criteria of the model is determined by how well it is able to attain these objectives. There- fore, in order to be able to judge and criticize the efficiency of the model, it is necessary to delineate clearly the purpose of this study. From a panoramic viewpoint, the goal of this thesis is to under- stand the physical nature of the switching process in electrothermal devices. Experimental evidence to date(40'41'42) indicates that the electrical switching process must involve a complex combination of phys- ical phenomena. For instance, the thermal nature of the device has been clearly established, yet, obviously, its extremely fast switching time cannot be accounted for by thermal processes alone. So, a primary ob- jective is to study and develop a better understanding of the physical processes related to the initiation of this fast switching process. 60 61 A change in band gap and mobility gap has been presumed to be involved in the switching process. For simplicity, in this model, both of these were assumed to be abrupt step functions of temperature, the step occurring at the critical temperature, Tcr' It is proposed to Show that such a change does indeed initiate a fast switching mechanism coupled with a change of state. The ”goodness" of the model will be judged on the basis of comparisons of the switching times from actual experimental data. The terminal I-V characteristics of the device also form a focal point of interest, since it forms the interface between the device phys- icist and the applications engineer. A typical I-V curve for a current negative resistance, shown in Fig. 1.4, is reproduced here in Fig. 5.1. For convenience, this curve is Split up into four regions. The first or low current region extends from the origin to the switching or threshold voltage vthr' The associated terminal current at this point is called the threshold current Ithr' As shown in Fig. 5.1, this is a Stable region with a positive resistance coefficient. When operated here, the semiconductor will act as a non-linear resistor. The dynamic properties of a device operating in this region are Studied in the next section. Since all points in the region are stable, all bulk dynamic characterists should tend to approach a steady state value with time. The next region is termed the transition region for obvious reasons; it is a highly unstable switching region with a non-linear negative resistance coefficient. Being unstable, this region cannot be obtained as Solutions to a static model. Adler and Kaplan's solu- tion(36) of the static equations for energy flow and continuity of charge closely resembles Fig. 5.2a which was obtained as a Solution of this more dynamic model. TERMINAL CURRENT 61a 1 STAILK HBO“ CURRENT REGOON UNSTAILI HIGH CURIBNT IN!“ N I6 |ON I (JR-[INT RIGION Ii _. 733mm... venue!“ Vs»: Figure 5.1 Typical I-V Curve of ETD Showing Various Regions 62 The third important region of the curve, the post-switching region shown in Fig. 5.1, is called the unstable high-current region. It exhibits a slight negative slope; i.e., a negative resistance coef- ficient implying that the region is still unstable. In fact, as other subsequent experiments Show, the device is in a "thermal runaway" con- dition. Many of the earlier experimental failures with the ETD re- sulted from a failure to control this post-switching phenomenon. Though the temperature and hence the conductivity increases very rapidly in this region, the power dissipation, given by (ES, reduces, due to the decrease in the electric field E} This reduces the rate of temperature rise and effectively increases the terminal dynamic re- sistance of the device, until finally the resistance coefficient becomes positive. This fourth region is called the "high-current" region. Con- trolling the length of the third region such that this roll-over takes place before the temperature attained is too high and the device burns out is a current problem in device design. This dissertation, however, concentrates on the preswitching region; it is of interest here to observe the onset of the switching process rather than the postswitching stabilization process. Hence, it is the first two regions, namely the low-current and transition regions, that are of interest here. Another subject of interest are the methods of inducing switch- ing in the ETD. Four distinct types of simulated experiments are per- formed. In the first class of experiments, the external battery vol- tage V8 p is small; so that in the final steady-state condition, the p ETD remains in the low-current zone. Therefore, in this class, no switching is expected to occur. This type of simulation is useful both in the construction of the low-current zone of the I-V characteristics 63 as well as in the study of the development of steady-state conditions within the sample. ' The other three methods are attempts to induce switching by vari- ous mechanisms. The first obvious method is to apply a ”large” applied voltage such that at least some points within the sample reach or exceed the critical temperature. In the other two methods, switching is attempted by locally applying a power pulse and an optical ionization pulse. The switching characteristics observed, if any, sometimes provide a basis for correlating with existing experimental data, providing further vindica- tions for the goodness of the model. Other characteristics provide a basis for making useful pre- dictions regarding the behavior of the ETD as a switching device. After all, the model's validity and usefulness is dictated not only by its ability to explain observed phenomena, but aISo in its ability to pre- dict and foresee phenomena not yet observed. It then becomes a useful tool with which to guide future experimental research. In this set of simulations, the experiments are performed using the data of V02 as a prototype material. The selection is so partly because, in the past decade, V02 has been widely experimented with by various investigators. Thus, the documentation on them is fairly good. Correlative work is also much easier, since the switching characteristics of V02 have already been observed. Predictions regarding this material is also likely to be more than merely academically useful, since opti- cal effects in V02 is a matter of current research. In this chapter, the first few sections concentrate on establishing the goodness of the model, while the later sections deal with the predictive aSpectS of the model. 64 5.1. Terminal I—V Curves The I-V curve obtained by this simulation (Fig. 5.2a) is in agreement with published experimental results in the following sense. Although the magnitudes of currents obtained by the simulation is very much higher than experimental data would indicate, the current density normalized to a unit volume sample is of the same order of magnitude. An experimental curve obtained by Duchene et. al.(40) is reproduced in Fig. 5.2b for comparison. The discrepancy is therefore due to the large sample size chosen. Recall that the choice of a large sample size was dictated by the Debye length for the material; in order to neglect the surface effects at the metal-semiconductor interface, the grid size, and hence the sample size must be chosen much larger than the Debye length. The current densities at the threshold voltage Vt and sustaining vol- hr tage V are compared with experimentally-obtained figures in Table 5.1. sus The sustaining voltage obtained from the simulation is 0.85V which is in excellent agreement with the experimental data; the break- down or threshold voltage Vthr is, however, about an order of magnitude smaller. This is attributed to the oversimplification involved in com- puting the volume charge densities at the boundaries. In actuality, the diffusion current component due to thermal and concentration dif- fusions result in much larger charge densities near the ends. One would therefore expect large electric fields to be set up near the electrodes, and, therefore, considerable voltage drops would occur near the ends. Recall, however, that the sustaining voltage obtained compares very well with the experimental data (see Table 5.1). As borne out by the other experiments discussed later, the switching observed is indeed with fila~ ment formation. Once the filament extends from one end of the material to the other, the large electric fields at the ends can no longer exist; ' CIIIEII(A) B 5'. F 64a _;5Iv «ed—q» u 1.8 VOlIIIE(V) Figure 5023 I-V Curve of ETD from Simulation eu9e5n1(mA) VOLTA6E(V) Figure 5.2b Experimental Current-Voltage Characteristic of a V02 Device 65 thus, the error involved in the boundary conditions of the volume charge density does not have a significant contribution. Table 5.1 Constants of the I-V Curve Depicted in Figs. 5.2a and 5.2b Experimental Simulated 4 2 4 2 Isus/cross-sectional area 3'1 x 10 A/m 3'20 x 10 A/m vsus 0.8V 0.85V I 3 1 104a/ 2 1 27a/ 2 thr/cross-sectional area 0 x m - m vthr lB-V 1.39V 5.2. Steady-State Conditions in the Preswitching Region In this simulation, the voltage of the battery shown in the cir- cuit diagram of Fig. 3.5 is small. The power dissipation at any point as described by the temperature is given by r73- . E; hence, the larger the battery voltage applied, the larger is the eXpected power dissipated at any point, so the final temperature attained would be expected to be larger. By a small battery voltage it is meant that the voltage is suf- ficiently small so that no point within the material attains, in the steady-state, the critical temperature Tcr beyond which switching is expected to occur. As evidenced by the I-V characteristics discussed in the pre- ceding section, the preswitching region of the ETD is stable. So, for small voltages, the ETD should reach a steady state in a time t large compared to the thermal time constant of the system; all variables of the system--V, Q, n, T, etc.d-should tend to a final steady state value. That this is indeed the case is evidenced by the curves (Fig. 5.3) ob- tained from a simulation of the model. 65a P M II v I I do— ‘- ‘\ 0". “v \ .L I l I h l flormlml llnits Figure 5.3 Behavior of [5qu Variables with Time in Low Current iegion 66 In Fig. 5.3, the potential V, dissipated power density P, temper- ature T, and free election density n are plotted with respect to time for a representative point in the sample. As expected, all of the above vari- ables tend to approach their reSpective steady-state values. The temper- ature T starts from the initial ambient value Tamb to the final value Tcr’ the rise can be approximated by an exponewtial laws T I Tamb + (Tf - Tamb)(e-tft1) The time constant can be computed from the plot. The value of this time constant’ti is compared in Table 5.2 with the thermal time constant Tk. It is also interesting to note that the other variables V, n, etc. exhibit the same exponential nature. Whereas the potential V re- duces with time, decaying to the final steady-state voltage Vf, the power density P and electron concentration n for representative points in- crease with time. But, in each case, the approximate law is exponential; typically, for voltage V and the electron density n for a representa- tive point, reSpectively, the relations are: V=V +(V a - -c/1 f Vf)e 2 PP / ‘t t3 n a nf + (nf n0)e Note that the time constants T1, Ti, and 23 are approximately the same and close to the thermal time constant (k. Table 502 Time Constants Computed from Fig. 5.3 T1 630's t2 SBOps 7:3 720’s 1A 550,6 67 Therefore, in the preswitching region I, the thermal properties are primarily reSponsible for the tranSport phenomena in the ETD. This is why the sample Berglund model<34) , which utilizes only the equation of energy balance as discussed in the introductory chapter, works so well in the preswitching region. Consider also the temperature-current curve shown in Fig. 5.4b. This curve is obtained by varying the input current in the simulation runs. More precisely, referring to Fig. 3.5, this is the curve obtained by keeping the battery voltage constant and varying the external resis- tance Rs to alter the current through the device. The curve (Fig. 5.4a) shows that the time required to reach the switching temperature Tor de- pends on the input current. This is reasonable since the power dissi- pated in a volume element of the sample depends on the square of the current density. It was also observed that the time required tr to attain 952 of the final steady state value depended on the input current: the larger the input current, the larger the final temperature and also the shorter it took to reach the final steady state value, the latter relationship being approximately linear: i I cltr + c2 The values of tr for various input terminal currents are tabu- lated in Table 5.3. A simplistic consideration of the energy balance equation lends substance to this observation. Neglecting the diffusion term, the energy balance equation can be written as: c 33 - o-‘i . E (5.2.1) at Since the electrical conductivity U'can be written as 0"- nelu and n, the free electron density can be approximated by the ionized charge density, it can be written as an explicit function of temperature: 67a H r P H E. O I E3 a.» I- ii 1.2 1 J | I ! ' I.” ‘02 . tittout Figure 5.4a t , vs, 1 in the Absence f Band Gap Switching 't «b e 1 Inuutm (IL ”I i—-——-——-I———_I.—__._L____. 1.0 1.2 cuttlnt Figure SOAb Final Steady State Temperatures in the Absence of Band Gap Switching-vs-Terminal Current 68 e Nae/LexpC-Eg/ZkT) The above relation uses the donor ionization equation developed earlier. Table 5.3 Data P10tt€d in Ftfifle 5043 and 5.4b V t i Final app r Terminal Current Temperature Battery Voltage Rise Time (M) (K) (V) (m5) 1000 328.3 1.875 1.80 1100 333.7 1.875 1.70 1150 338.1 1.875 1.64 1300* 356.5 1.875 1.50 *No band gap switching was simulated in this case. The approximation nt~rng is particularly good in the preswitching region because of the absence of temperature and free electron concentration gradients in the y and 2 directions (see Fig; 5.6). Returning to Eq. (5.2.1), and substituting for the electrical conductivity 1'1: . NdefA exp(-E /2kT)§ . E dt c 5 Since both E3 and f.are step functions of temperature (refer Figs. 2.1b and 2.3), and assuming the electric field E'to be a constant also in this preswitching region, the Eq. (5.2.3) can be integrated. Note that, although no physical justification is given, the assumption that the electric field E.is constant in the preswitching region has been used (10,31) extensively by many experimenters with considerable success. To simplify the integration involved, consider conductivity to be a linear function of temperature given by 9’. «0(1 Mu) - W‘BKI +°tT) Extensive examples of this kind of linear approximation to o'can be found in the literature<31'34). 69 Rearranging terms: T t f dT . { noel" 32d: T _— amb l + dT 0 c Integrating: 2 E _1_1.(1+o- *HI3N7THV7.dInI —¥- tr 2.7127.3 W. ,3»- ‘9 '0 ll.“" Figure 5.12b Power Dissipation Density-vs-x During Initiation of Density - Log P(x, 34x, 4A2) - vs - x 80 switching occurs, the power dissipation increases tremendously: this follows logically, too, since the dissipated power is proportional to the square of the electric field which is also increasing rapidly as noted earlier. This, in fact, is the cause of the thermal runaway phenomenon; the rate of power dissipation increases so rapidly that it more than offsets any reduction in the rate of increase of the electric field. Outside the filament region, the power dissipation reduces slightly during this time. Along the x-direction, the profile is roughly the opposite of the temperature profile. Starting from an initially flat profile, the power density tends to Show a minima at the center and the highest power dissipations at the two ends. When a power pulse is super- posed, the power density profile tends to shift showing a maxima at both the point of impressed power (2Ax, 3Ay, 4A2) and the diSplaced ”hot” point (4Ax, 3Ay, 4A2), which is consistent with the observations made in Sec. 5.5. The heat flux resembles the power dissipation profile. (See Fig. 5.13.) The maximum heat energy flows out at either end via the conducting electrodes; the minimum occurs at the center. Once switch- ing occurs due to the occurrence of a "power pulse”, the heat flux den- sity is low both at the defect point (ZAx, 3Ay, 4A2) and the diSplaced hot point (44x, 3Ay, 4A2). 5.7. Optical Switching of the ETD In this experiment, a random interior point--chosen again to be the point (2Ax, 38y, 4Az)--is locally excited by a continuous exter- nal light source such that the point is totally ionized. In the com- puter model, this is simulated by setting the ionized donor density n+ d at the point equal to the total donor density Nd. The effect is almost [06 P 80a es ‘\ /.__._ 1...... _ 7‘ I4 .. / 6 - O- t I 14m «.6- * C 8-4". c—u-I t . 2'5“ -' '-4>-*tu.2dfrf731‘“flb ' - —x— t a ennui». #‘--.Iflagr:ISEgN.____fl.____4r____4p___.'h_———4F""e"'-" ii .- I. l :li ‘ar L5 —:L.. 02 0‘ " '. 5:0 Ill; Figure 5.13 Spatial Distribution of Heat Flux Density Showing Initiation of Switching - Log/Jq/(x, 34y, 4A2) - vs - x 81 immediate; the device switches to the unstable high-current region and proceeds, as in the previous switching methods, to experience thermal runaway. Some interesting observations can be made with regard to this experiment. In this experiment, the external bias voltage Vapp is low such that the ETD reaches a steady state just below the critical tem- perature. The light pulse is now applied, and switching is almost im- mediately observed. Note that in the case of a defect forming in the device, there is still a fairly large ( SQPS) storage time associated. In this case, however, it is much, much less than lOps, and hence opens up very interesting lines of application for the device. For one, since the thermal runaway is controlled by device design and external circuitry, there is a very fast photoswitching device available. Furthermore, the switching can be originated at will anywhere within the sample: this can be developed into a fast read/write unit which would have increasing applications in interfacing with computers and other high-Speed equip- ment. CHAPTER VI CONCLUSION The objective of this research effort was primarily to gain an understanding of the physical principles underlying the initiation of the switching phenomenon in the ETD. Experimental evidence(6’8) and literature surveys(s) indicate that this is a complex phenomenon, and a number of interdependent variables--1ike V, p, n, T--p1ay important roles in determining the switching characteristics of this device. To understand the mechanisms involved in its switching behavior, it is necessary to understand the interrelationships between these variables and be able to analyze the significance of the contribution of each. The means chosen to assist in achieving this objective was to construct a computer-based model of the ETD. The model embodies all the available physical information known or held to be true regarding the ETD, and is based on well-known physical laws. It is used to simu- late the behavior of the ETD in order to study the mechanisms involved in the switching process within the bounds of the model; in other words, the scope of the model is limited to the study of the initiation of the switching process so as to be able to develop a more comprehensive under- standing of the physical principles associated with it. Some of the results obtained from those simulations were compared with available experimental data in order to verify the goodness of the model. The model is useful in understanding the distributed bulk pro- perties of the ETD. This information is not available experimentally, R7 83 but is important in understanding the overall behavior of a given device under various electrical and thermal boundary conditions. In addition, the model possesses the capacity to simulate conditions which have not, as yet, been experimentally investigated, thus giving it a predictive aSpect. This feature enables the model to be used to design optimal devices and understand the practical limitations of existing ETD'S. 6.1. Accomplishments of this Research The major accomplishment of this research project lay in the physical principles that were uncovered and the generalizations that were made regarding the behavior of the ETD. This project has focused on the preswitching region and the processes involved in the initiation of switching. In relation to the I-V characteristics of the device, the stable low-current region and the unstable transition region described in Fig. 5.1 form the primary areas of interest. The thermal character of the preswitching region was firmly es- tablished by the model. In the stable low-current region, the thermal parameters like the thermal conductivity k and the Specific heat c, th and the thermal boundary conditions almost completely determine the tem- perature profile as was Shown by the excellent co-relation with the solution of the one-dimensional steady State equation: d2 dx2 8 =21: C where‘P, the dissipated power density of the sample is assumed constant. Also, the time constant associated with reaching a steady state was in close accord with the theoretical thermal time constant 2; as compared in Table 5.2. Therefore, both the final steady-state temperature pro- file and the time constant involved in reaching this steady state were associated with thermal parameters. Also, the other dependent variables 84 of the device show an exponential time dependence, the time constant involved being again close to the same theoretical time constant‘?;. This agrees with the expectations of numerous theoretical workers in the area like Boer(10), Dohler<31), Ovshinsky<12), and others. The thermal time constant is also associated with the large storage time observed in the device. However, the switching time in the device is very fast-of the order of nanoseconds--and cannot be attributed to thermal diffusion effects. In this transition region, the model succeeds in describing the simulated behavior of the ETD in close parallelism to actual experimental observations made on these devices. The model describes the switching mechanism as the propagation of a band gap reduction through the length of the material. Such a re- duction caused a change in the local electric field which propagated along the length as the switching progressed, as illustrated in the electric field snapshots at succeeding time Steps depicted in Fig. 5.11. The model is also able to simulate the switching behavior of the ETD induced by different methods. Consider, for example, the phenom* enon of switching induced by the introduction of a defect in the mate- rial. In the first place, the results agreed with the experimentally observed fact that the switching first occurred along an axial line containing the defect. But, that is not all. It also localizes the point of origin of the switching process at a point away from the de- fect. From hindsight, this is readily understandable and reasonable Since thermal diffusion prevents the formation of a double-humped char- acteristic with humps at the defect point and at the midpoint. Hence, the competing forces result in a maximum forming in between the two points. H 85 The combination of these two regions, the low-current and the transition region form part of the Static and dynamic I-V curves to- gether with the unstable and Stable high current regions which lie out- side the scope of this present model. The overall agreement of this curve to the experimentally observed I-V characteristics of the ETD is further testimony to the validity of the model. Thus, the model is not only capable of providing an understanding of the interactions in the interior of the device; it is also able to relate to the easily observ- able static and dynamic terminal I-V characteristics. Since the I-V characteristics form a crucial interface between the device engineer and the circuit designer, the ability of the model to depict these charac- teristics give it considerable value. 6.2. Predictive ASpectS of the Model The predictive aSpects of the model are useful guides to further material and device development. This is largely due to the capability of the model to present experimentally inaccessible data regarding the nature of the variables in the interior of the semiconductor material, and so contribute significantly to understanding the physical processes in the ETD. Thus, using the model, one can predict the nature of the electrical potential and charge distribution in the interior of the material; one can also predict the temperature profile, the power dis- sipation density and the magnitude of heat flux density at various points in the interior in both the low-current and transition regions. It is thus possible to study the interrelationships that exist between the various electrical and thermal variables. For instance, the power dissipation curve predicted by the model provides a vital link between apparently unrelated experimental observa- tions. The switching times observed in the ETD'S are of the order of 86 nanosecunds<34). This is in agreement with the simulated results of the model: in fact, the switching mechanism can actually be observed progressing through the device by studying the electric field and tem- perature profile snapshots of the interior. Also, it has been experi- mentally observed that, prior to filament formation, thermal hot Spots appear near each electrode. The ETD simulation also projects similar results as shown by the temperature profile of Fig. 5.8b. The simula- tion predicts that the power dissipation curve (see Fig. 5.12) is approxi- mately an inverted parabola. Since the thermal time constant is so much larger than the switching time, thermal diffusion is practically non- existent once the switching process has been initiated. Hence, the power dissipation accounts for the exceSsive temperature rise near the electrodes resulting in the experimentally observed hot Spots. Simulations were also made to Study various methods of initiat- ing switching within the device. Switching was initiated a) if the external circuit parameters were such that the tem- perature at a point within the sample exceeded the critical temperature, b) by simulating a defect wherein excessive power was dissipated resulting in local temperatures greater than the critical temperature, and c) by simulating a light pulse focused at a point on the ETD operating under suitable bias conditions. These first two methods agree with experimental observations; the third method has considerable predictive value. In this case, the temperatures at all points within the device were below the critical temperature Tcr' The incidence of the simulated light beam creates local donor ionization, and this is sufficient to initiate switching in the 87 device. The simulation results therefore predict that under suitable bias conditions, any energy source capable of causing local ionization is suitable to switch the ETD. Thus, radiant energy of various fre- quencies and even acoustical waves may be used to initiate switching in ETD'S under suitable bias conditions. 6.3. Limitations and Deficiencies of the Model By virtue of the Stated objectives of the model, its scope is limited to simulating the ETD behavior in only the low-current and tran- sition regions of the terminal I-V curves. The third region, namely the unstable high-current region is portrayed to some extent. Thereafter, thermal runaway occurs, and excessive temperatures are generated in the device. The fourth region, the stable high-current region, is thus never reached. Ridley<30> has pointed out that the stability of the ETD depended heavily on the external resistance Rs. This would limit the current through the device and hence the power dissipation. In the model, the iteration scheme* is such that, for large Rs' the scheme is divergent. The iteration scheme referred to is the Inner Loop* elucidated in Chapter 4. Here, the terminal current is computed by integration near the edge of the device, and the terminal voltage is computed therefrom by apply- ing Kirchhoff's voltage law. In Such a scheme, the limiting effect is not obtained since a negative voltage would still satisfy the voltage law, and the current continues to increase. Obviously, for simulation of this region, a revision in the iteration scheme is required. Although such a revision would limit the power dissipation, it still may not limit the temperature to realistic values. The crux lies *See Fig. 4.10 88 in the assumptions that the thermal variables k , the thermal conduc- th tivity, and c, the Specific heat, remain constant beyond Tcr' For sta- bility to be achieved within a reasonable temperature, both c and kth should increase: i.e., the material should require more heat input to raise its temperature and should also conduct heat away faster. In addition to this, there are other deficiencies in the model. For one, there is a discrepancy between the experimentally observed thres- hold voltage and the Vt r observed in the simulation runs. As pointed h out earlier, this is because the assumption that the net charge density remains constant is not strictly valid due to the external circuit con- figuration. One way of correcting for this would be to reconstruct the circuit so that a charge source instead of a voltage source supplies energy to the semiconductor. An alternative method would be to re- evaluate the charge density at the ends at each time step and use this variable charge density as a boundary condition. Another source of numerical disagreement is in the magnitudes of currents: however, this can be adjusted by using sample dimensions closer to the Sample sizes used in experiments. There is also a quantitative error between the switching times observed in the Simulations and experimental data. Though several reasons for this have been discussed earlier, one pos- sible explanation is that it is due to a breakdown in an assumption of the model. In the model, all magnetic field effects were neglected. However, in the transition region of the I-V curve, the current in the filament region changes very rapidly. By neglecting the magnetic field, the inductive effects of this change is neglected, too, which may affect the switching time, though the extent is difficult to determine. Incor- porating the magnetic field in the equations is difficult, as a wave equation is to be solved instead of the Poisson's equation. The wave 89 equation is a vector equation and thus is comprised of three scalar equa- tions, and both computer time and memory requirements is increased tre- mendously. Obviously, when modifying this model, a balance or compro- mise has to be readhed: there is a definite tradeoff here between the accuracy and numerical predictability on the one hand and Simplicity, ease of interpretation, computer memory requirements, and costs on the other. 6.4. A Possible Simplification of the Model Because of the existence of this tradeoff, it is worthwhile to consider schemes to simplify the model and reduce the computer costs of simulation. One Such example is outlined here. In the preceding chapter, the temperature equation, cfl-VZT'BE Est was reduced to a simple one-dimensional steady state equation dZT . is c dx2 The analytic solution of the above indicated that the tempera- ture distribution should be parabolic in Shape. Then the axial tem- perature distribution obtained by simulating the low-current Steady state condition was fitted to a parabolic curve, and the best-fit parameters obtained were very close to the analytic solution. This indicates that reasonable accuracy of simulated results (say within 52) can be obtained by considering only the temperature equation in this region of operation where the temperatures lie well below the critical temperature. Let the starting ambient temperature be Tamb' Tcr the critical temperature, and TD the temperature difference defined by T I T - T D or amb H 90 If T is the highest temperature tested at any time, then as long as 151‘ +f(Tcr-T ), amb amb where f is a preset positive real number less than unity, only this sim- plified model needs to be solved. At higher temperatures, however, the other parameters will begin to play a more significant role, and the more complete set of equations will have to be Solved. The preset factor f is a function of the material parameters of the ETD, and a good estimate can be obtained by preliminary runs with the present model. For example, for V02, a good value for f is 0.9. Using the simplified model, the temperature at the maximum is identical to that obtained from the one-dimensional analysis, as P would be assumed constant. Thus, for f I 0.9, the expected error in the maxima estimated from Table 5.5 is less than 52. 6.5. Spggestions for Future Research From the results obtained in the simulation runs thus far, some suggestions for future areas of research can be seen and are presented below: 1) An extension and broadening of the objectives of the model to include the final stable steady State of operation is a prime area for research. This would probably necessitate a revision of the computing scheme as discussed. Such a revised model would be able to simulate the behavior of an ETD in all four regions of the I-V curve. 2) Some of the assumptions of the model need to be relaxed to incorporate more of the characteristics of the real device and thus obtain closer numerical agreement. In this regard, the two suggested targets are the assumption of constant H 3) 4) S) 6) 91 charge density at the electrode interfaces and of neglect- ing the effects of magnetic field. For the former, an iterative scheme to compute the charge density at the inter- face at each time step is suggested. For the latter, the Poisson equation VZV I ~9/g which uses the relation E a -VV should be replaced by an appropriate wave equation. Since such modifications would represent a substantial increase in computer time and memory requirements, simplifying modi- fications of the type outlines in Sec. 6.4 could also be Simultaneously made. A single defect in an otherwise homogeneous material has been Simulated on this model. A possible extension would be to observe the switching phenomena when a statistical distribution of defects exist in the material. Individual material parameters of the device, like kth' c, Nd, etc., could be varied in the existing model with a view to isolating and analyzing the effect of each parameter on the device switching characteristics. This could provide considerable impetus to material development and guide ex- perimentation in that area. Environmental parameters like the ambient temperature and boundary conditions could be varied. Such studies would be of importance to both the circuit designer and the device engineer. To the circuit designer, a study of the device characteristics under various electrical circuit bias con- ditions would also be of interest. The model predicts that local ionization by the impinging of an external energy would cause the device to switch 92 under Suitable bias conditions. Optical, u-v radiation, or other radiations could be tried. Acoustic waves could also be tried. APPENDI CBS l i APPENDIX A Derivations Using the Boltzmann's TranSport Equation In Sec. A.l, we shall derive the basic charge carrier tranSport equation used in Chapter 2 from the Boltzmann's equation. Subsequently, in Sec. A.2, we arrive at expressions for the tranSport parameters, namely 03 the electrical conductivity, D, the concentration diffusion coefficient, and as the thermal diffusion coefficient. A.l. The Transport Equation Let fO be the equilibrium distribution function for the electrons in the semiconducting material for all temperatures T in the operating temperature range. Under the application of forces during operation, the perturbed distribution is f, where f is a function of ; Space and velocity Space f = f(‘f', V) Using a first-order perturbation theory, f can be written in terms of the equilibrium distribution fO as f .. £0 + Vrfo . Fe + vao . gt (11.1) where T'is the net time constant associated with the perturbation. This is known as the Boltzmann tranSport equation. Consider a unit volume of the sample. The carrier current den- sity due to the flow of electrons is given by: Ji a :ne < v > (8.2) 93 94 where < v > is the expectation value of the electron velocity. The expectation value < v > is defined, for a given statistical distribu- tion f, as: < v > I £(Vfdv (A.3) f;fdv V where v'is the total velocity Space. USing (A.1), _ r _ .9. .. < v > =[Iv(f0 + (vvrfo) + (v .vao ))dv V .5de Hence 3- : -ne(£f/f°dv +£ITWVr~f°JV7 + £ V t-vao ° t't )dv) 1:4dv Defining the root-mean-square velocity by the expression w2 8 I52 I and, in cylindrical co-ordinates, .x g fng a {:10 foz‘wzsinededadw The expression for the carrier current thus reduces to: d I .. .— - f 2A V . J +[nxe(f\2rf0dv + [0 f0 0 Cw + rfOSmededsdw “ x + 6 f0 4) 'Cv(VrfO . V)wzsin8d9d0dw] fvf 0d? For simplicity, it is assumed that, in the model, the equilibrium distribution is Naxwellian. 3/2 'fl. f0 = n( H) ) e 2kT 2KkT Since there are n free particles per unit volume, the integral of the denominator is obviously n. Also, the first integral in the numerator is zero. Thus, 95 ._ 26 J .. [foofx f -?e anrfOSinededfldw 0 0 0 (a) + ID” I; (fiscveo . V)w2.'meded¢dw] (b) Consider the term (b). The term % Simplifies as follows: % -- a7. = .1. a - 1 .- dt m dt m where F is the applied force. Since the force in this model is due to the presence of the electric field E I -eE "3| '1'? <|° e m The term (b) thus becomes: (b) - [-efowfo‘ fozxuvvfo . V)wsineded¢dw]E Obviously, this component of current is the conduction component 3cond' since it is proportional to the field E} The quantity within paren- thesis must therefore be the electrical conductivityOf: the integral is evaluated later in this section. The other term, (a), is the diffusion term 3diff; it comprises of the concentration diffusion component and the thermal diffusion com- ponent. Consider the quantity V5f0° foo I fig-V“ + inf—O.VT Br1 3 Since 2 I2 ~me fO I n( n1 )3 ZlkT Hence - I312 f 3/2 ”tr :3 0 8 ( n1 ) e I :2_ ‘arl ZKkT n 96 and 2 2 f 3/2 ’—-‘§“ 2 -m" 3 n n( m ) [-3e' RT + . e Zij er 2xx 2T 57 2kT2 T 3/2 - mw . n( m )3/2e 2n m2 _ 3 kT 2kT2 2r - f.£_mw2 - 2.1 ZRT2 2r Substituting the above relations back into termil: (a) I Fdiff a [féméx fozgwa :2 SinGdOdewJe n n (C) co 1: 2 + [f0 f0 [0 Te w4f0( mw2 - l)sin8d6d0dw]VT ZkTZ T (d) From the form of the equation, it is apparent that the quantities (c) and (d) enclosed in parenthesis are the proportionality factors D and d, the concentration and thermal diffusion coefficients reSpectively. A.2. The Proportionality Constants (a) The Electrical Conductivityp" The required integral is 3' ‘ 1— 2 _. cond ' -33 fomf; 4’2! (VVTO . V)W SinQdedeE m and so 0' a -_3_ {wag 162; (va0 . 'V)wzsin8d9d¢dw m Expanding Vyfo in Spherical co-ordinates, 3f A M .3foe+1 09+ 1 afoa 0 am- .755- .——51.9 Te Since a Spherically symmetric Maxwellian distribution has been chosen, the partials in 9 and O are zero. Hence 97 Assuming T'to be a constant, the integrations in 6 and O can be performed: K f0 SinGdQ B 2 fort“ a 2K Since the other terms are independent of 8 and 0, the integrations in 9 and 0 yields: w33f (radxezfoo-t Odw .__ 0 .___ m 8w Since 2 -mw 3 2 fo‘“( m )/e2kT 2XkT 2 -m_w__ 3:9 - n3/Ze 2” - a”. 3w 2KkT kT Therefore 2 -mw -.-. + Me fownm( m )alszae dw m kT 2KkT Substituting 1 2 .. (.11.) ’ w ZkT W4 = u4(2kT)2 m dW 3 (2kT)1/2du m Hence 2 so 3 "U 1/2 tr: + 4Ke2n f0 _‘_t_ m lze (ZkT)2 u!J . (ZkT) GU kT‘K3/2(2kr)3/2 m2 n1 2 a + A‘Kezn {000 2? uae u (111 mKT3 2 98 From Tables 2 fowure-u du : 1/2 (t ' 1) see 2 e ‘3- s l.,/r 2 2 2 2 Assuming I to be a constant, T3+4K92nXlslsloJ%-s 2? 2 2 2 n1x372 I 3nezz m Defining the average electron drift mobility by [“8 = 3 3: (Add m the expression for the electrical conductivity becomes: U“: hep. (2.3.1) This is the eXpression arrived at in Sec. 2.3. Thus : rE' 3 1'1th Jcond (b) The Concentration Diffusion Coefficient: D The current density due to concentration diffusion has been obtained earlier as 1'1 c» K 2K 4 .. : ff _1:w f Jdiffc (16 f6 0 O sin6d6d¢)eVn and so D, the concentration diffusion constant for the mobile electrons, is given by F 21' 4 D = {0°70 f0 3" fo sinededaclw n Since the Maxwellian distribution is Spherically symmetric and T, the perturbation time constant, is assumed constant, the integra*ion in 9 and 0 can be separately performed. Since 1’2 X 0 dfl I 2X 99 and r f Sin8d9 I 2 0 Hence -mw2 D I fgat4IwQC n1 )3/2e 2ET'dw ZXkT Substituting 2 2 mw-ii 2kT aniV 2kT du I 4’ m dw ZKT Thus .. - 2 1 2 D 3 f0 art<£§l)2“a(_fli_)3/2 1 e u (.231) / du m ZkT K372 m an 4 - 2 .. 4 (21.1“) foueudu "" .JK m Using the Tables as before, this integral reduces to: D'fli; legals-J? In .JK 2 2 2 =51o3 m Recalling that We, the electron mobility, had been defined as: We I 3e? m the diffusion coefficient becomes: which is the well-known Einstein relation. (c) The Thermal Diffusion Coefficient: _g The carrier current density due to thermal diffusion is -’ K 2r' 2 Jdiffd :- [fowfo ,6 eTwaf0(;mw - _l_)sinGde¢dw]VI‘ 2k'ri T 100 and so, 2x OC .. fOOOfOK/g) e‘fwhf0( mw2 - l) sinOdOdde 2kT2 T Evaluating the integrals in 0 and O as before yields 05.. 4 I Orwefwaf0( mw2 - _3_)dw ZRTZ 2T Recalling that and making the substitution w I (252)1/2u m the integral becomes 2 «=4xf0°°e‘ru“(2kr)2 n (m )3/2e“(9_2_-_3__)J'3__'kfdu m 2'K372 2kT 2 1‘ 2T m 00 _ A- fl.8kne(f u6eudu-3fwueudu) __ 0 —-0 4r m 2 Using the Tables as before I8kfl8(105.13/7‘303010‘y/FII) I8kne s_6_.\/.i'- "hfi- 16 I kn . 3e? m I‘Pkn The net diffusion current density 3diff is therefore: Jdiff I eDvn +«VT I e ekTvn +,ukn\7T e I :kV(nT) APPENDIX B The Finite Difference Equations Consider a continuously differentialbbounded function f(x, y, 2). Using a Taylor's series expansion about a point(x , 20) an adjacent 0' yo point f(xo +-Ax, yo, 20) can be expressed as: f(x0+ Ax, yo, 20) I f(xo, yo, 20) + Ax if 3x X I X0 Y'YO z I 20 3 4'sz 32f +AX3 3f + see (Bel) 7F”— 31'— 8x x I x0 'Bx x I x0 y-YO y-yo z I 20 z I zO Consider first the point (x0, yo, 20) to be in the interior of of the material. In such a case, both points f(xO +mAx, yo, 20) and f(xo -»Ax, yo, 20) exist. Expanding f(xO - Ax, yo, 20) about (x0, yo, 20) using the Taylor's series as before: 2 2 - n -A f(xO Ax, yo, zo) f(xo, yo, 20) x 2£_+-Ax a f ex 20 2 3x 3 ”AX 3f + see (B02) 3: B 3 x all partials being assumed to evaluated at (x0, yo, 20), the subscripts having been omitted for conciseness. Subtracting (8.2) from (8.1) yields the following: 3 3X 30 3 Ex 101 102 Rearranging terms : if. 2: f(x0 +Ax, yo, 20) - f(xO -Ax, yo, 20) 3x x - xo 2Ax Y ' 3’0 2 - 20 - AXZ 33f (8.3) 3! a 3 x As discussed in Sec. A.l, a discretization error of Oosxz) is tolerated. So, the truncated expression for 91 is obtained. For an interior point, a B f(x0 +Ax, yo): 20) - f(xO - Ax, yo, zo) if; axx-x0 ZKx Y'YO 2920 if; can be similarly obtained as: 3 y g f(xo, yo +Ay, 20) - f(xo, y0 - Ay, 20) 2i 3y x I x0 ZAy Y " Y0 z = 20 To obtain the second derivation 3222:, add equations (8.1) and (13.2) to obtain: Bx f(xo +AX, Yo! Z > + f(xo - AX, yo! 20) a 2f(x0) yo, 20) O x + X2 32f +2A4;af+ooo 3x2 4. 3x Rearranging terms and dividing throughout by x2, 2f . f(x0 +fx, yo, 20) + f(xo ~Ax, yo, 20) - 2f(xo, yo: 20) x2 X'X sz X4 34f + 000 (3.5) 4! X4 + N > (V Discretization error terms of order sz or higher is neglected as before: 103 3Zf a f(xo +Ax, yo, 20) + f(x0 -Ax, yo, 20) - 2f(x0, yo, 20) 3x2 x I x0 X2 y ' ’0 z - 20 (8.6) Similarly, the second order derivative in y is: . 32f . f(XO' YO +AYO' '0) + f(XO, ’0 "AYO’ Zo) " 2f(xO, yo, 20) 3y [5’2 N (3.7) Consider now a point on the edge x I 0. In this case, the point (x0 -lix, yo, 20) does not lie within the domain of the semiconductor. But, (x0 +:Ax, yo, 20) does, and f(xO + 2Ax, yo, 20) can be expanded by Taylor's Series about (x0, yo, 20) to yield: f(xo +~2Ax, yo, 20) - f(xo, yo, 20) + 44x gg_+ gZszz 32f Ex 2: 8x2 + SZAX)3 39f + 000 (8'8) 33 3 3 x Multiplying (3.1) by 4 yields: 2 2 l:f(xo +Ax, yo, 20) I 4f(x0, yo, 20) + (sax it; + 1: Ax 3 f g x + 4 AXB 33f + 000 (309) ‘1" 3. ax3 Subtracting (8.8) from (3.9) and rearranging terms: _ -3f(xo, yo, 20) -+ l:f(x0 +Ax, yo, zo) - f(xo + 2Ax, yo, zo) ii, 31: XI!) ZAx y - ’0 Z ' 20 3 + 4 x3 3 f + so. (3010) 3! 9’3 X 2 Neglecting terms containing A): or larger power of Ax, the equation be- comes I 104 9i ,3 -3f(xo, ’0' 20) + 4f(x0, yo +Ay, 20) - f(xO + 2Ax, yo, 20) 31 x I O ZAX y - yo 2 ' 2o (8.11) Similarly, the partial derivatives and at y I 0 and z I 0 respectively are: . -3f(xo’ yo, 20) + [st-(X0, ’0 + Ay, 20) " f(xové ZAY' zo) 22. 3 x I X M Y 0 y Y ' 0 z " 20 (3.12) 33' _ -3f(x0, yo, 20) + 4f(xo, yo, 20 +-Az) - f(xo, yo, zo +-2Az ) ‘3: x I x 2A2 O Y ' Yo z I. o (8.13) The second partial derivatives for this family of points is ob- tained by multiplying (8.1) by 2 and subracting it from (8.12). Per- forming the above mentioned operation and rearranging terms: .azf _ f(x0 + 2Ax, yo, zo) - 2f(x0 +uax, yo, 20) + f(xo, yo, :0) 2 2 3X x-O X y-yo z I 20 + x 3f + ... 2: 3 x Neglecting second order terms, the second partial derivations are: i I f(xo + 2.x, yo. 20) - 2f(xo +AX, yo, 20) + f(xo' YO, 20) x2 151 2f _ f(x,, ’0 + 2Ay. 20) - 2f(xo. yo Mir. 20) + f(xo. y 2 ‘Ayz f . f(xo, yo, 20 + 2A2) - 2f(xo, yo, 20 +132) + f(xo, yo, 20) :2 4:2 0’ 20) Nel The third type of point is one lying on the edge x I LX' wherein the point (x0 +-Ax, yo, 20) is not defined within the semiconductor, but (X0 - 253a yo, 20) is. Similar Taylor's series truncations can be made 105 by utilizing the expansion of f(xO - ZAX, yo, 20), the method proceeding along similar lines to the previous case. Neglecting second order terms the first and second partial derivatives are written below: 3 - - - ?_f_ - f(XOp YO, ZO) 4f(x0 AX, YO, 20) + f(XO de, yo, 20) 38 x I L 24x x Y " 3'0 2 I 20 if. a 3f(x0, yo, 20) - 4f(x0, y0 -Ay, 20) + f(xO, yO - 2A)., 20) I L y Y z =- 20 3-f_ B 3f(x0’ ’0' 20) - “f(xo' yo, 20 - Z) + f(xo, yo, 20 - 2A2) 37' x 3 X 2A2 APPEKDIX C NURMALIZATION In Chapter 1, it was pointed out that round-off error could be minimized by normalizing variables and reducing the equations to dimen- sionless form. The details of such transformations are worked out in this section. The normalizing variables are: Ax I fixzdy I 113:;st = 32;?) V I Vapgz 10 10 10 Bv 0100 .G 3 _ . V. R app 2 T (Tct Tamb)e + Tamb’ n pnb 9 L x nd pnb’ T - Tcr - Tamb where pnb ' pn0 - p0 pO ' "b/efh 4. pna - “deli ‘7 '“h {9 I lnfl Using the truncations discussed in Sec. 4.1, the relevant equa- tions are transformed as follows: 1) Poisson's Equation: vzv . J/e fi+ozv+fi+9le . 0 3x2 ’y 322 106 107 NOW 3:1” V0‘ +Ax, Y) Z) + V(X -AXI_YJ 2) ‘ 2V(XL Y: Z) 2 2 3X AX . 0 I vapp 1 O . [¢(x +~Ax, y, z) + 9(x -Ax, y, z) - 2¢(x, y, 2)] L 2 x The other derivatives can be similarly obtained. Substituting back into the Poisson Equation and using the relationships Dav“, Davy, Davz, and fidefined in the text, the final relation is obtained below: Davxfl + Davyfl + Davzfl - Yd” +X I 0 2) The Continuitwaguation: it+Vo3.o 8:: 33_ + V. (-nepVV +’AkV(nT)) I 0 at Using the charge neutrality relationship to express the equation in terms ofg , yields %’— + v. (- engvv +3:va +fek7g(n: - 9/e)T) . 0 t The above equation is now re-expressed in terms ofIX by using the trans- formations used earlier in the text. 2 a): +Paepnbvapg sz . (”Ii”) + 100,? vamp v. (mum) 5? Ioovme 100v." + ”a Paepnb sz qume) + Pa k1'cil‘ma’atpxr>"x V-QVOW = 0 10" VIP 100W 1. 2e app x Defining the following constants for conciseness 2 F' I vOP'aEPnbe looevapp 2 F' I epnbe 1006 108 0 Fe . Mavapp 0 Fa " PaVO we obtain :1 - Fav. (awn) + Fé V. mm + F1.) v. omen + F5 7.7”“) t I0 3) The Temperature Equation: ch + kt vzr arE . E 3t h ,3}. - flavzr - nefav .vv at c Applying the same transformations used earlier in this section and considering the equation term by term, we obtain 22." Tali at at kth va I ktth v26 C C Using the truncation sceme developed in Appendix B, the righthand side of the above equation reduces to 2 2 V 9 I (Davx9 + Dav 8 + Davze -‘1D9) . lOO/Lx Y and dafining Pd and 1.1. as Pd I pnbepavapploo L 2 x 7‘1? . ktthIOO cL x th 6 Temperature equation becomes 39 .. - a? (Davxs + DavyG + Davzs d9)/ 1' + Pl-‘d APPENDIX D ALGORITHM FOR SOLVING A SET OF TRIDIAGONAL EQUATIONS In Chapter 4, the partial differential equation for temperature was reduced to a set of tridiagonal equations by using finite difference techniques. To solve such a system, consider a linear set of n triI' diagonal. equations whose subdiagonal, diagonal, and superdiagonal coef- ficients are at, b1, and c1 respectively, the righthand side'is di, and 01 18 the ith variable. The ith equation, therefore, reads A recursion solution of the form ut ”‘1 - °i“1 + 1 (9.2) :‘t is found to be valid. The coefficients 0‘1 and P1 obey backward recur- sive relations of the for:- ‘1 - d1 " ‘1“1-1 (8.3) 71 p1 - bi - ‘1": - 1 (9.4) Since rewriting the first equation of the tridiagonal set yields “1 - i; " c1 “2 b1 Ff by °°lnparing coefficients, 4" and ’51 can be evaluated. ’51 " b1 “1 " d1’P1 109 110 From the recursion relations (8.3) and (0.4), all coefficientsai1 and Pi can be evaluated, and, employing the applicable boundary conditions to obtain "i and un, all ui's can be obtained. APPENDIX E TREATMENT OF THE EDGE POINTS BY THE 30R METHOD In Sec. 4.3, a finite difference scheme was outlined for solving the electrical potential distribution in the interior of the semicon- ductor. The behavior of the edge points is governed by the boundary conditions as outlined here. Consider first the edges in contact with the metal electrodes. These edges, as discussed in the boundary condi- tions in Chapter 3, are constrained to have V(Lx' y, z) I O Normalizing, ¢(Lx. y. 2) I 0 This takes care of all points on the surface x I Lx' For points on the surface x I O, the boundary conditions are obtained from a con- sideration of the circuit equation L L V(o, y, z, 1;)- Va -f yI ’(46 i! + I‘a kgn‘T) )dydz PP ZIO yIO 3x Normalizing the above equation and using the truncated grid structure discussed in Chapter 4 “0. y. z. t) ' 1 ”32,222. 101* Zigwoqm) " wanjo k LX100 " 3-1 kIl where §;'-sx:§§; 3x x I O In the computer, the double summation is performed using the Simpson's integration rule. 111 112 Consider next the surface y I 0. Here, 21-0 By and so ¢(X, 0, Z) ' ¢(X, AY; 2) where the righthand side is an interior point which is evaluated by the 80R method. The other three surfaces in contact with air can be simi- larly evaluated. APPENDIX E LISTIVG This section c-ntains a complete listing of the computer pro- gram used in the simulation runs and consists of a main program and thirteen other subroutines. The LAIN program controls the overall logic flow and the Outer Loop referred to in Chapter A. The MAIN program calls TEMPEQN which sets up the tridiagonal equations corresponding to the Temperature Equation that is solved by TRIDAG. It also calls CONTROL, which controls the Inner Loop. CONTROL, in t rn, calls the subroutines CIRCUIT, PSSN, and LMATRIX which solves the Circuit, Poisson, and Continuity of Charge Equations reSpectively. The remaining subroutines perform auxiliary functions. BYPASS sets up the initial conditions, PRINT is concerned with all output Prints, and RHP computes the power dissipation density required for the [solution of the Temperature Equation. The function subprogram SIMPS performs the integrations required in the Circuit Equation using Simpson's Rule, and DELX evaluates the truncated difference operators. Finally, BUFFIN and BUFFOUT control the flow to and from auxiliary memory locations on tape enabling data at any intermediate times to be Stored. The actual listing is given in the following pages. 113 114 PRUORAM MAINIINPUToUUTPUToTAPEOOIINPUTcTAPEbI=0UTPUToTAPEIoTAPtZs ITAPEJ) NEAL LAOLY COMMON/A/IIoJJoNloNJ COMMON [ATIIOIIOIilotbiliolloiiioLCIIIoIIoII) COMMON/LOUP/KUUNToTIHE COMMON/TsU/Pll(lloIIoliioPPTIIollvill COMMUNIRST/HSTiIIoIIoII) COHHON/UtbK/TNETAIIIQIIvii) COMMON/VOLTb/PHI(IIsIAsIII LOHMUH/LUL5/UIAO(II)ohUHIIIIOSPR(II) CUMMON/MAI/I‘Vlsvfh COMMON/FH/PAoINUAsLloPH COMMON/ATRA/FtofllolAUb COMMON/3101A CUMMUN/H/CAobAotbAofibd COMMON/VIVRQVIHOVU CUMMON/bAMS/OAMD COMMON/DUN/GAMCOALAMUA CUNMON/LIMITS/CMAXoCMINoUELVoNA COMMON/VU/VA COMMON/N/NUAIA COMMON/MU/GMUOIUUEVSDQNAX CUMMON/HEADtR/HEAUIoHtAUdoHEAUs COHMON/UUT/CO:CDENS:PF«CH COMMON/DIR/AIJIsIAIZIoJAIZIsKAIZI COMMON/VAR/VANTOVAHCOVANVOVAHEoVAHPoVARUoVARIQVAHFQINTAOINTMOINICO IINTD COMMON/ITEPS/EPSMAXIoITMAXI COMMON/CUN/GAMMAoAMDAsbAMflobAMHNoUAMMd COMMON/FACTORS/FAoffivfCoFU COMMON/Aw/ITMAAOEPSMAAQI COMMON/I/IHOOIAOTU COMMON/CONDTV/RATIU COMMON/PULSLITIMCIoTIMKF tUUIVALENCF (UELVoPK) C" SET INIYIAL CONDITIONS "' a F0NMAI¢IN0~II¢EI:.A.:K:: C 0' 1aurr=0n NU HUFFLR out 0- C9. INDEX DETERMINES INITIALIIINO INPUT " L.“ L IS TAPE NO. ON NHICH INITIAL INPUT IS HtCORUEU “I HEADIbOo717) INDEXoIHUFFoL C " COMPUTER TIME CONSTANTS '0 READIOOoTAT) dsDELIsTIMAAsEPSMAAsCX 747 FORMATISEIbo7) C 0' OHIO SETS AND UNIVERSAL CONSTANIS '0 HEADIOOOTQ?’ LloLYsUEPSoULTKsECH c '0 VARIABLE PARAMETER SET 0- KEAOID007Q7I EMUAQEMUBsébHoEGLoQI C 0' CUNSTANT PARA StT .. ' NEADC¢097~7T THKoSPHTgTAHUoTCRoSIOO C " INTEGER CONSTANTS " NEAUIbOoTITI IIoJJoNIsNJsITMAXoKMAXoNAoND READ(bso787) UELvoNxA READI6007e71IIMtloTIMtI MEADIOUO7I7)(IAIKIOJAIKICKAIKI9K=I0NUI thocso.777: EPSHAXIvllflAAIoflI 777 FORMAIIEIZsQOIIOI PRINT TTIOEPSMAAIOITMAAI AIIISIHA SA(2)IIHY $A(3)81Hl HEADI‘IOHDISTRTUUTI HEADZ'IOMUN ALONG HEADJ3IOHDIRECTION ’8’ FUHMAT‘FIOAJ-IIOD MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN SINU MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN JULY“ MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN JULY“ MAIN JULYe MAIN JULYQ MAIN MAIN MAIN MAIN MAIN MAIN MAIN csusmsuw 98M 98? C C 565 115 FORMATIMIIUI NAIIO=EMUBIEMUA PMINI 987 PRINT 9M5.NA.NAA.NINII.KFIN.IIMAA.AMAA FONMAII IHUofiISAoIbII FORMATIIM I'TME FOLLoNINu VANMLES AML NA.NAA.NINII.NIINIIIMAonuAA o.) .... CUNSTANTS 0F TEMPERATURE 00.0 TUSTCR-TAMH TAsTCM/TU IMOzIAMM/IO VUsNLIMOIU/CCM "' TIME CONSTANTS "' ALPMA=IIII-I)/LXI"Z dETA SIIJJ-II/LYI'°2 IAUAsISPHT/THKI/ALPNA IAUHaI.0/(ALDMA°CMUA0VUI TAUC80EPS/SIHU ".“ CUNSTANTS UF TfiHP Lu” "0'“ MAMMA=ILAILY)'°£. bAMZ:Z-0“GAMMA bAMMBxl.9/0AMMA 0AMM2=9.bOhAMMH 6AMD=2.0°M.U°6AMMA LbAt-EbH/IVU'ECH'ZoOI tuna-EGL/IVO':CH'£.OI NI8LA/(SIMO'LY9'2) HS”C08.RI CA=CA'LY"2 VAPP8CA'IMSOMI) VTH8VAPP/VU VR3RS'SIGU'IIl-II/IMLTA.LA) FC-LA/LY ’ FE‘FC/ZOU FU”GAMMA/Zo° FH=fE/20 POSSIOO/ILMUA'ELM) PNAapu'iAPI-LuA/TMUI'LAPIEhb/IA) PNH=PNA-PU A=POIPNM SIGM=PNH'LCM'EMUA blblzPNAAECM'EMUB UELVsPNA/PO PRINI SbSISIhIoUELv FORMATIIMOoOSIUIa OotIo.doO UELV= 0.tIh.6) uAaAOEAPI-EbA/IMUI A=UA'EAPIEGA/THUI FU=SIOB'VAPP9'2/(TD°THMI CC=-VAPP/HS ' PF=THK'TU'ALPHA 028ALPHA0DEPS'VAPP FAsTAUH/IAOTAUCI FMsFA/VIM FC=GAMU'IHO . CASMIIIFA'A-FCI AMUA=UELTITAUA FCSTAUM/DELI vA-VAPP/e. BAMCII.O°bAMMb/MMUA ALAMOAaI.oI./AMUA MMUtEUM/EbA PI'JAIMI59265358979 UtV'oOI TOUEv5082.0'DEV'DEV bMusIMMu-I.0IIISUMTIZoO'PI I'ULVI uMUsGMUOtGA MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN JULYA JULY“ JULY“ JULY“ MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN JULY» MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN 65 66 67 69 70 7| 77 7M 80 8| 83 B“ 85 86 ' 67 8M 89 9| 92 93 9h 95 96 97 9B 99 I00 IOI I02 I03 I06 'I05 I06 I07 I08 IIM IIS II6 II7 IIB II9 I2! I22 I23 I26 I25 I26 I27 I28 I29 I30 IJI l i (a. c no 556 :57 ~75 “7“ 707 272 coo us [-4 H g;\ " LIMITING VALUES UF TNL CIRCUIT CUMRENTS AND VULTAGES ... CMINaMI/IHIORbI 63P0/9NA CMAK300{/IH/NATIU°0.2I PRINT SbboLbLoEbHofiMUAIEMUHOSPHIICAIPK " PRINT INPUT DATA 9' 9. IORMATIIHIOZUAQ. INPUT UATA IN RAT MKS UNITS‘c/lo Z‘EUL 3 .Oblo. 30IOAO’LUH 3 'OGIO 39/10 3. EMUAa .obIoo399A09EMUM 8 .QGIU.3o//o 9'SPHT 8 .QGIOATO. CUMNcNT UENbITY 3‘OGIUAJO.PK= .9 bIOoJ) PHINI 5570LA0L790£PSODLIKOCCHQTHKQTANHOTCNO IIMAA FORMATIIfloo I‘LA = ‘OUIOAJOIOAI'LV I 'OGIOAJ0/lo Z‘EPSR 3 .OGIOOJ‘IOXO.HLIN I .OGIOo30IOX0.E 8 '06I0.3!I0A0'KTH 3 .9 JUIJ.39IOA‘//o M'TAMN 3 'IGIOoI0I0A9.TCK 3 .IGI0o3oI0A99TIMAA 8 .0bI0037 PRINI 559 FUNMAIIIMoo' THL FOLLUMING INIthk VARIABLES ARE ..l. d” bTAHTKoIHUFFoTAPENOOIIOJJOPSSNITMAAQKMAAQNAIMAA ITEM UF INNER LP OIo.i/O. NAAICUNVUI ITLN CRIIIARNUITMAAI.I PRINT 7I79 INUEAOIMUFFQLIIIQJJ. IIM£=A.o NOUNT39 PRINT 30VAPP0TUOSIGBQUZOPNH FORMATIINIo/lllodhAo' CONVERSIUN FAC10RS° I/lv' VOLTAGE -"9IPLI3. So. VOLTS'0 dl/o. TEMPERATUHL .OIPLIJQSQ. UEDK.O a/lo. CONDUCTIVITY -'.0IPCIJ.50. MHUS/METEN'Q d/lo. CHAMbE DENSITY '-'0IPLIJo50' COUL PER METER CUHE'Q 3 //o. CARRIER UCNbIIY".9IPCI3050. PLN MEIER CUHE'I PRINT IOCC FONHATIIHOQ'CUNHENT CDLNss SIGO0III-IIOVAPPILA L68A0CUENS CH=IMNOIUONIILA PRINT ZOCGOLH PRINT “750CULNS FORMATIIHGO. CURRENT VRINT “7999‘ FORMATIINOQ. POULR ULNSIIY PRINT 70705IUUQPOQEMUA0IAMH FUPMATIIdOo' SCALt AHUEMS 'o/llo' sluu = 'vCI2o“oIOAo' P0 = IMOIOA. Oolllo. MUA= “OEI2.498A0.TAM88 ..Elz. “I PRINT 27B F0RMATIIMI03SXI. OTHER VARIABLLS'I LIST VARIABLES‘. PRINT ZISITNOQTAITDQVKQVIMQVUotbAQEbUIPNAvoRI0M5. IUAICAOFUOAMDAIGAMCA IALAMDAQUAMMAobAMZobAMMHOUAMMZQGAMDQFUIFHOFAOFBOFCOTAUAOTAUBOTAU. FONMATIIMOO' TH0309LIZAMODAQ' TA='9LI2o905Ao/QIH00 I. Ina'ofIdofiofiAo'VH='oEI2.“95A9/9Ih09 d. VTH8.9EI2.“05AO. V0=.oEI2.“05Av/o IMO. 30 tuAsioila.“oSAoO Lufla'oEI2.~cSAI/0IH09 9. PNA8’oEIZ.“05Ao' A3'9EI2.“95A9/9IH00 ITMAAIKMAAINAINAXIITMAAI ’-'oIPLI3obo. AMPS.) DLNSITY --‘.IPLI3.509 AMPS PER MTR SQUASE'I -—9IIPLI3.bo' MATTS PER MLTLR CUHE'I .QEIZO 0' HI8.9EI20405K0. “$3.05I20905A0/0IHOO 5’ bA=.0fIZoMQSXo. CA1.9EI2099550/0IH00 0' rD=.o£IZo“ODAv' AMUA3.9EI2.“obAo/QIM09 7' 6AMC3.CEIZ.*05A9.ALAMDA8.otIZ.“obAo/0IMO. d. AMMMS'OEIZ.“¢5AO. hAMZI'oEIZ.“v5Ao/OIH09 9. UAMMH='OEI20905A0.UAMM2 3.0£I209v5A0/9IH00 I. bAMH3'otI2o“oSA99 b=°9EI2.“05A9 /0IH°0 2' FHS'AEI2.“95Ao/OIHJO 3' FA8'otIfloMOSAo' F88.OCI20995Ao/0IH09 TAUA:'9£I2.“95A0/IIH09 TAUC”.EI2.“¢SAI «0 FC8'0EI2.“05A0' 5' TAUB’.0EI2.695AN' MAIN MAIN MAIN MAIN . MAIN JULY“ JULY“ MAIN MAIN MAIN JULY“ JULY“ JULY“ JULY“ JULY“ JULY“ JULY“ JULY“ JULY“ JULY“ JULY“ MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN JULY“ MAIN MAIN SIND SINU SINU MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN I32 I33 I3“ I35. I36 II I39 I“0 I“I I2 I3 I5 I6 I7 I8 I9 20 22 I“a I“) I““ I“5 I“b I“! I“8 I“9 ISO 23 I52 I53 I5“ I55 I56 I57 158 I59 169. "II I62 I63 16‘ I» I66 I51 I68 I69 I70 I“ I72 I13 I7“ I75 116 I77 I73 I19 Iuo IBI I82 I83 (A. 95 97 99 lo PRINT MoAUUNI.|IML NUATA=I NLAU IN FROM HUFFEN UN OPTION 0' CALL BYPASS IF IINUEA) 95099097 CALL SPECINIINDEAoLI DO TO 99 CALL dUFFINIINDEAoLI CONTINUE . TIME=TIMEODELT MOUNT=KDUNT°I IFITIMCoUEoTIMEII NDATASND PRINT “oKOUNTOTIME FORMATIIHIol/lo' MOUNT 3 .0I3o' TIME 3 'OIPEI2059. SECS.) IFIMOUNT.6T.KMAAI GO TO I2 CALL. CUIITKUI CALL PNINTIPHIIVAHVI C .9 A) RIGHT HAND SIDE " POWER C I7 I6 I9 ... M) LLFT HAND SIDE --- TEMPEQN CALL NHU CALL PRINTIZAQVARLI CALL PRINTIZCIVARPI IFIIHUFFI79996 CALL SP£CUUTIIUDFFIINULAI uU TD 9 CALL UUFFUUIIIHUFF) LUNIINUE CALL IEMPEUN IFIIRUFFIIG.19.I7 CALL SPELUUIIIBUFFIINULAI 00 To I“ CALL HUFFOUIIIBUFFI CUNIINUL c0'0 SET up MAIMICES coo coo I2 bTUHE AND PRINT OUTPUTS ”'9 IFITIMfioLTATIMAAI GD TO IO CONTINUE END MAIN MAIN MAIN SIND SINO MAIN MAIN MAIN MAIN MAIN MAIN UUMMY MAIN MAIN MAIN JULY“ JULY“ MAIN MAIN MAIN JULY“ JULY“ MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN MAIN JULY“ MAIN MAIN MAIN MAIN MAIN MAIN I6“ ‘I85 I86 I88 I89 I92~ I93 I9“ I95 I97 I98 I99 26 202 203 20“ 28 205 206 207 206 209 2I0 2II 212 2I3 29 215 2I6 217 2I8 2I9 220 C 50 S“ (\h FDOY‘O SUMROUTINE TEMPLUN DIMENSION USISIIIIoIIoIII COMMON USIIIIoIIoIII.AIIIoII.III.uIII.IIoIII COMMON/vaun/ COMMON/TAU/P/IIIIOIIOII)OPIIIOIIIIII _ LOMMON/VAR/VAKTIVAHLQVAKVQVAREIVARPQVARUOVARIOVANFQINTAQINTROINTCO IINIn COMMON/dUMK/A(IIIoDIIIIoCIIIIonIIII COMMON/LOLS/UIAMIIII.SUMIIII.SPM(III LUMMUN/DUM/GAMCoALAMUA COMMON/IIMES/TAngquIAutoud UIIIOIIOIII COMMON/V/VMQVIHOJU COHHflN/A/IIOJJOHIONJ LUMHON/CUN/LIAMHA o AMIIA v (IAflgf o {IMM‘III o IIAMHt COMMON/L/LAOLG COMMON/T/THOQTAQIU tUUIVALENLEIUbToUbTbTI LUUIVALENCEIJJIKAI “0 I0 I3I0II SUHII)"U.5 SPHII73’009 CONT INUE SPNII73‘I00 5UHIIII3'IAO SOLVE IN I’DIHLLTIUN DO 50 I=I9II DO 50 J=IoJJ DO 50 NIIIJJ UIIOJONIIUIIOJQKI’TMO PII.J.KI=0II.J.AIouIIoJcNI/AMUA CONTINUE LAIZ t LMONJ OO 5“ I'IQII OIAOIIISALAMUA CONTINUL DO 5 K-IOJJ UO 5 J'IOJJ DO 6 I320NI AIIISIUIIOIoJoKIOUII'IIJIAI-2.U'UIIoJoKII‘0o5 IFIIJ.NE.II.AND.(J.NL.JJII IMIII'IUIIoJ’I'KIOUIIIJ'IOBI'2o0'UIIoJoKII'GAMMA IFIJALUoII OHIIIsIUIIQJOIOKI-UIIOJQNII. MIII‘IUIIoJ-IohI-UIIOJIKII'GAM2 IFIIKoNE-IIoANU.IfigNtoKKIICIII=IUII0JOK°II‘UIIQJOA'II-2.'UIIOJ9VII IFIJ.EO.JJI O'OAMMA IFIinfloII CIII‘IUIIOJQAOII‘UIIoJoKIIOUAM2 IFIKoEQoKNI CIIISIUIIQJOK‘II-UII0JoKII‘UAMC DIII2MMS DII73PIIOJOKI'UIIOJQ“I/LMHUA _ PIIOJOKIaPIIQJoKIOAIIIOHIII°CIII DIIIIPIIOJQKI TERM: NMAND Is COMPLETE N0 CORRLCIION NEHUIMéU AI MOUNUARILS aINCE BOUNDARY VALUL5 ARLZVRO SOLVE HY TMIUAG MAIRIA MLTMUU CALL THIDAG SOLUTION METUMNED IN DIII USTIIoJohI'0oc USTIIIOJOKI'ooc U0 5 L320NI U5TIL0J0KIIOILI CONTINUE PAN! 2 --- CONSINUCIINM USISI MAINIA LA'I 6 LM'NJol 11‘ TEMPEUN TEMPEUN TEMPEON TEMPEUN TEMPEON TEMPEON TEMPLUN TEMPEON TEMPLUN TEMPLUN TLMPLUN TEMPEON TEMPEON TEMPEON TEMPEON IEMPLON TEMPLON TEMPLON IEMPEON TEMPLON TEMPEON TEMPEON TEMPLON TEMPEUN TEMPLON TLMPEON TEMPLON TLMPLUN TEMPEON TEMPLUN TEMPLUN TEMPLUN IEMPEON TEMPLON TEMPLON TEMPLON TEMPEUN ILMPEUN TLMPEUN TLMPLON ILMPLHN IEMPEON TEMPEON IEMPEUN TEMPLUN IEMPLON TEMPEUN TEMPEUN TEMPLON IEMPEON TEMPEON IEMPLON TEMPtUN TEMfiEON TEMPEON TEMPEON TEMPLUN TEMPEON TEMPEON TEMPEON TEMPEON IEMPLUN TLMPEON TEMPLON ~C¢:N¢7UIOIHN (‘fi 119 D0 56 J=I'JJ 56 UIAb I JI =UAMC ”0 I5 I320NI U0 I5 K=I~JJ U0 Ib ngvJJ IFIIJoNfoIIoANfloIJoNLoJJII IMIJISIUIIQJOIORIOUIIOJ'I0KI'2.O'UIIOJOKII'oob IFIJ.EO.II ZMIJISUII0J0IQKI'UIIOJQKI IFIJoLNoJJI JhIJI‘UIIOJ'IQKI’UIIQJQAI AIJI'IU5TII‘lecKIAUbTII’IoJokI'200.U5TIIOJ9KII'GAHME PIIoJoKI=bANMHGPIIIJIAI‘HIJIOAIJI OIJI=PII9JIKI HMS IS CUMPLEIE CORRECT FOR BOUNDARIES CALL IRIUAG bOLUTION RETUHNEU IN UIJI U0 I5 L3IIJJ I5 U5TSTIIOL¢KI=OILI UO 25 J8IIJJ DO 25 I=ZOKI no 26 KaloJJ IF IKoLIIo I I ICIKICUIIOJOROII'UIIIJIKI IFIKACOAJJI 2CIKIIUIIoJoK-II'UIIQJQRI IFIIKQNEOIIOANOOINoNLoJJII JCIKISIUIIOJOKOIIOUIIOJQK’II‘2AO'UII9J0KII'Oo5 IFIJoifloII . IUIKI=UST5TIIQJOIOKI-UbTDTIIIJQKI IFIJofiaoJJI IuIMIsUSISIIIcJ-IoNI-USISI(IoJ.AI ITIIJoNEoIIoANDoIJoNLoJJII IDIKI8005.(U3I5TIIOJOIOKIOUbTbIIIoJ-IOKI’ZAO'USTSTIIOJOKII PIIOJOKI'PIIOJQAI‘CIKIObIKI UIBI=PII0JIKI 00” 9H5 I: CUMPLLIL “*9 .... CONNECT FUN UUUNUAklwb ...“ CALL TRIUAU 50LUTION RETURNED IN IN“) 00 25 L=IOJJ USTIIOJOL7=0IL7 (0.x dODNDARIES ... DO 220 K=I~JJ DO 220 JzIoJJ UIIoJoAI=IML - UIIIOJIKISIHU 220 CONTINUE L'. ALL INTERIOR PIS “' DO 225 I320NI DO 22) J320HJ DO 225 M320NJ 22b OIIQJIKI3U5TII0J0RI'IMU MA=2 no 7 AL=I.JJ.MJ M=MA IFIKLALUoIII M3NJ 0. ALL VACCS '° DO 253 I'ZONI U0 253 J=LONJ UIIoJoRLI=UIIoJoMI UIIoALoJI=UIIoMoJI 2‘13 CONT INHL LUNI I NUL '0 ALL LUwEb “' IEMPLON IEMPEUN ILMPLUN IEMPEON TEMPEON TEMPLON TEMPEON TEMPEON IEMPEON IEMPLON TEMPLQN TEMPLUN [EMPLON TEMPLON IEMPLON TEMPLON ILMPLON IEMPLON ILMPEUN IEMPLON IEMPtuN IEMPLUN IEMPLON TEMPLON TEMPLuN TEMPLON IEMPEUN IEMPLON TEMPLUN TEMPEON IEMPLON IEMPLON TEMPLUN IEMPLON TEMPLUN IEMPLuN ILMPEQN IEMPLUN ILMPEUN IEMPLUN IEMPLON TEMPLON IEMPLUN IEMPLUN IEMPLQN ILMPLUN TLMPEON IEMPLUN IEMPLUN IEMPLON IEMPLON ILMPEON IEMPLUN [EMPEUN ILMPLUN ILMPLUN ILMPLUN ILMPLUN IEMPLUN ILMPLUN ILMPLUN ILMPLUN ILMPLNN IIMPIUN ILMPLUN IEMPIQN IOI I02 I03 I0“ I05 100 I07 I06 I09 IIO III II? II) II“ II5 II6 II? IIb II9 I20 IZI I22 I23 I2“ I25 I26 I27 Ian I29 IJO III 235 UU 63‘.) J‘IOJJINJ “0 235 I329NI bIIoJoII‘UIIoJQZI UIIoJoJJI=UIIoJ9NJI UIIOIOJI=UIII20JI UIIIJJoJI=UIIoNJoJI CONTINUE CALL PRINTIOQVARTI RETURN LNI) ILMPLuN IEMPEON TEMPEQN TEMPEON TEMPEQN IEMPEUN TEMPEQN TEMPLON IEMPEuN TEMPLQN I32 ~I33 I34 I35 13o I37 I38 I39. 140 I“I ‘. N NUo—WO‘ 121 SUHRDUTINE PRINTIVIVARI COMMON/N/NDATA COMMON/A/IIOJJINIONJ COMMON/DIR/XIJIOIAIZIOJAIZIOKATZI COMMON/HEADER/HEAUIoHEAUZvHEAUJ DIMENSION VIIIOIIOIII DU 5 I=Io3 PRINT IoHEAUIoHEAU29AIIIIHEAU3 U0 5 K=IONDATA IX=IATKI IY=JATKI II=KAIKI 60 TU (899.10) I PRINT JOVAROATIIOJAIKIQKAIKI PRINT BIIVILOIYoIlIoL=IoIII IIO TO 6 PRINT “OVAROIAIKIOXTIIOKAIKI PRINT 29IVIIA0LOIZI9L=I9JJI UU TO 6 PRINT IOVAROIATKIOJATKIIXTII FORMATIIHOISAOAIOOZAIIIO'O'OAII.9.0III FORMATIIHOISAOAIOIZAIIII’I'IIIO'9.0AII PRINT 20IVIIX9IY0LI9L=IOJJI CONTINUE CONTINUE FORMATIIHOIROKIZAIOIABIAIOI FURHATI IHUOSAIAIOIBAIAII’I’OII9.9.9III FORMATTIHOOIIIEII0“9IAII RETURN END PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT PRINT 70 5UHROUTINE CONTROL COMMON/LIHITS/CHAAQCMINQUELVONX COMMON/NU/UHUCTUUEVSU'NKK COMMON/A/IIOJJONIONJ COMMON/RST/RSTIIITIIOIIT COMMON ZATIIOIIOII)IZUTIIOIIOIIIOZCTIIOIIOIIT ITLR=I SL=0 $N=0 . DU 70 I=IvII 'UO TO J=IIJJ 00 To R=|.JJ RbTTIvJoKI=lHTIOJOKT C“. THUS RST CONTAINS OLD CHI FOR ENTIRE INNER LOOP" 85 co. 90 can 80 IS 76 Co. TOO N=N0l H=ITER CALL LHATRIA CALL CIRCUIT CALL PSbNTITERT IFTthR-MT 80080990 " TtNOENCY TO DIVERUL on on L=L°T -- TENUb to CUchRoL " -' PRINT “.ITEROHOLON . IFIITtRoLLoNXX) GO TO 75 IFTN.oE.NAI Go To 75 00 TO 85 CONTINUE DO 76 I=I9II DO 76 J=I9JJ DO 76 K=IOJJ RSTTIoJohT=lCTIoJ9KT RbT STORES DIFFN CURRENT DENSITY 9. RETURN PRINT “'ITEROHOLON PRINT 5 CALL LXIT FORMATTIHOO' ITER = .0139' M8 .9139. L8 .OIJO. N I .OI3T FORMATTIHOQ 'THJS I5 TENOING T0 OIVEROE O) [NO . o CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL' CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL CONTROL . CONTROL CONTROL CONTROL CONTROL CONTROL OGNU‘U’TJ‘UN ‘ .2 ..HFF ‘— ¢ «I? LL-) anROUTINE RUFFOUT (THUFFI BUFFOOT COMMON lATIIoIIvII)oZUTIIOIIoIIToZCTIIoIIoIII BUFFOUT CONRON/LOOP/KOUNTOTIHE BUFFOUT COMMON/TUO/PZITIIOIIOIIIoCONDJTIIoIIoIII RUFFOUT COMMON/RbT/RSTTIITIIcIIT . BUFFOUT COMMON/OEbK/THETATIIoIIoIIT BUFFOUT COMMON/VOLTb/PHITIIOIIOIII . ‘ BUFFOUT LII . BUFFOUT «0 IFTIUUFF.LT.O) GO TO 50 BUFFOUT C '. IRUFF I: 0 ' UUFFOUT C '. BUFFER OUT ONLARELLEO ULOCK " '. .' BOFFOUT C P“ " IN ORDER. SIbMAoRHOoPOUER " BUFFOUT 73 CONTINUE A HUFFOUT IF TUNITOLT TJQOJCTR975 UUFFOUT 63 BUTFER uUTTLoITTlATIoIoITTZCTIIoIIoIITT HUFFUUT on TO 00 ‘ BUFFUUT C '. IUUFF lb - . HUFFOUT C '. " HEPFER OUT LAHLLLLO BLOCKS) bOFFOOT C " IN ORDER 9 KUUNTvTIMEoNTnJCONUoJUIfFoU9V 0' BUFFUUT 50 CONTINUE BUFFOUT IF TUNITOLT 50.55974075 BUFFOUT $5 UUVFER OUTTLoI)TKOONToPHITIIoIIoIITI BUFFOUT 60 IMUFF=-IRUFT UUFFOUT PRINT TITQIUUFFoLoKOUNT ' UUFFOUT Coo THE ABOVE ENSURES THAT LABELLED AND UNLABELLED BLOCKS ARE BuFr.o ALT BUFFUUT 90 CONTINUE BUFFUUT IFTUNITILT 9099I079975 BUFFOUT 91 CONTINUE BUPFOUT RETURN ‘ BUFFUUT 74 PRINT 30 ' UUFFOUT 36 FORMATT' EOF 0N LAST OPERATION .o/lv' lSKIPOI'I . BUFFOUT PRINT TITOIHUFFoLoKOUNT UUFFOUT CALL EXIT UUFFOUT 75 PRINT 37 HUFFOUT 37 FURMATT' PARITY ERROR 'I . BUFFOUT CALL EXIT BOFFOUT 72 FORMATT ¢ UNIT NOT READY i) BUfFOUT 717 FORMATTBIIOT HUFFOUT END BOFFOOT ccsamcuw 1:3 bUUROOTINE BUFFINIIRECNU0LI BUFFIN 2 COMMON ZATII0IIOIII0ZbTII01I0II)0ZCTII0II0III UUFFIN 3 COHMON/LOOP/KOUNT0TIHE BUFFIN 4 COMMON/TUO/Pll(IIOIIIII)0CONUJTIIOIIOIII UUFFIN 5 COMMON/RST/RSTTII0II0III I BUFFIN 6 CUMMON/DEGK/THETATII0II0III BUFFIN 7 COMMON/VOLTb/PHI(IIOIIOIII . BUFFIN 8 C P. INUEA DETER MINES THL STARTING POINTOL IS THE UUFFIN 9 C .. REQUIRED TAPE NO. RLAO IN MAIN BUFFIN IO C '. ‘° L IS THE REQUIRED TAPE NO .. BUFFIN II RERINO L UUFFIN I2 ISKIP=IRECNO'd-2 BUFFIN I3 ICOUNT =0 I BUFFIN I“ IFTIbKIP.EU.0) GO TO 527 . BUFFIN I5 C .. .. SKIP TU T"L REQUIRED RECORD 0. O. . UUTFIN IO DO 50 I=I0ISKIP BUIFIN I7 52 ILUUNT=ICUUNT0I bUFFIN I8 IF (UNTT0LI 52053074075 - BUFFIN I9 53 HUIFER INTLOIITKOKI bUFFIN 20 SO CONTINUE HUFFIN 2T 527 PRINT 7I70ICOUNT0I5KIP0IRLCNO0L MUFFIN 22 C" .. BUFFER IN ALL ARRAYS ... UUFFIN 23 73 CONTINUE BUFFIN 24 IF TUNIT0LI 73065074075 UUFFIN 25 65 CONTINUE BUFFIN 26 60 HUFFER IN (L01) IZATIOIOIIOZCIIIOIIoIIII BUFFIN 27 S CONTINUE UUTFIN 23 If TUNIT0L) 5088074075 ' . BUFFIN 29 68 BUFFER INTLOIIIKOUHT0PHITIIOII0IIII BUFFIN 30 b FORMATTIHI0' KOONT = .0140. INDEX 8 '0I40. TIME I'0EI6070. VAPP : BUFFIN 3I I'0EI407I BUFFIN 32 REAUIb007I7T LL UUFFIN 33 70 CONTINUE BUFFIN 39 IFTUNIT0LI7007707407S BUFFIN 35 77 CONTINUE BUFFIN 3b KK=KUUNT'LL UUFFIN 37 PRINT 60KOUNT0IRECNO0TIHE0PHITIDI0II HUFFIN 38 IFTKK oNEoIRECNOI CALL EXIT bUFFIN 39 RLTURN I BUFFIN 40 74 PRINT 3b ' UUFFIN “I 36 FORMATT' EOF ON LAST OPERATION I0//0' ISKIP0I'I BUFFIN 62 PRINT 7I70ISKIP0I bUFFIN 43 fl CALL EXIT UUFFIN 4“ '75 PRINT 37 ' UUTFIN 45 37 FORMATI' PARITY ERROR 0 BUFFIN ab CALL EXIT UUFFIN b7 72 FORMATT ’ UNIT NOT READY “I BUFFIN Q8 7I7 FORMATTBIIOT - UUFFIN #9 END ’ BUFFIN 50 C9. 366 I20 I30 SURROUTINE LNATRIX . .' THIS IS CALLED HY LHATRIA TO EVALUATE.RHS OF CONTINUITY EON ... COMMON/CUNIOAMRA0AHOA0GAM20OAMHB0OAMR2 COMMON/OEbK/THETATIIoIIOIII COMMON/VOLTS/PHITIIOIIOIII COMMON/LOOP/KOUNTOTINE COMMON/PULSE/TINEI0TIMEF COHHON/LINITS/CHAKQCNINOUELVONA COMMON/TIO/PllTII0II0III0TEHPIII0II0III COMMON/RaT/RSTTII0II0III COMMON PSITII0IIIIII0CHIIII0I10III0OZITII0II0IIT COMMON/A/II0JJ0NI0NJ COMMON/VAR/VART0VARC0VARV0VARE0VARP0VAR00VARI0VARF0INTA0INTROINTC0 IINTD ' COMMON/FACTORthAOFBOFCOFO COMMON/SIOIA COMMON/B/CAOOA0EGA0EOH CONNON/V/VROVTHOVO COMMON/GANS/GAHU COMMON/CONDTY/RATIO COMMON/T/THO0TA0TO COMMON/TTEPS/EPsRAA.ITMAA COMMON/VD/VA COHHON/XTRl/FCORIOTAUU Fl'2o0.FE CONTINUE DO I ISI0II DO I J8I0JJ DO I K8I0JJ IFITHETATI0J0KI0LToTAI GO TO 130 UlITIIJOKIBALOGTRATIUI EOcEOR PSITI0J0KI8‘FC/RATIO no TO_TJa ' CG‘EGA PSITI0J0KI8'FC Ull‘IQJO‘I’Uo TENPTI0J0KI‘P5IIIOJ0KT.HSTII0J0RI PZI(I0J0KI3UA.EXPTEO/THETATIOJOKII IFTTTINEobToTIHEIIoANuoITIHEoLEoTINEFTI PZI(30405I:PZIIJ060SIO IUELV CHIIchoKI=PlITIOJOKI.THETATIOJOKI CONTINUE ‘ UO 3 I3I0II UO 3 J=IoJJ DO 3 K8I0JJ ATLHP‘IOJOKI3TEMVIIOJONI9F“.IU‘VAICHIOIOJQKT.UAVTICNIOIGJOKI 'OUAVZ(CHIOIOJOKI’UANU'LHIII0J0BI I TEHPII0J0KI8TEHPII0J0KIOFR'THELAIDllvI0J0KI'DELATCHI0I0J0KI OODELYTCHI0I0JTKI‘UELYIUZI0I0J0KIOUELlICHI0I0J0KT‘UELITUZI0I0J0"II TENPIIOJOKT'TENPIIOJOKI’FA'P‘IIIOJ0KI'I UELXTUIIvIoJoKI'DELXTPHIOI O0J0KIODELYTOZI0I0J0KI'DELY(PHI0|0J0KIOUfiLZ(PRI0I0J0KI'OELZTDIIOI0J OoKII TENPTIOJOKI‘TEHPTIOJORI-FA.TDLLATPZI9IOJ9KI’UELXIPHIOIOJOKI I’DEL'TP‘IOIOJOKT.DELTIPHIOIOJOKI.UELIIPZIOI'JOKIPUELZIPNIOIOJOHII C09 PSI IS THE COEFF OF -CHI " PSITIoJoKI3UAVlTTHETA0I0J0RIOOAVYTTHLTA0I0J0KI I'DAVZITHETAol9J0KI'GAHN'ThtTATI0J0KI ‘.ZOVTH.RSTIIOJOKIOP§IIIOJOKI Z‘FAOPIITI0J0KI PSITIOJ0RI=P5ITIoJ0KI’TUELATUZI0IOJOKI'TOELXTPHI0I0J0KI'VTH 0‘DELATTHETA0I0J0KIIOULLVTuzlOIOJOKI.TUELYTPHI0I0J0KI'VTH' ODELYITHETAoloJokIIODELZTUII0I0J0KI'TUELZIPHI0I0J0KI9VTH-OELZTTHETA ‘OIOJOKIII LHATRIA LRATRIA LNATRIA LNATRIA LNATRIA LHATRIX JULVQ JULYu LHATRIX LNATRIA LHATRIA LRATRIA LRATRIA LHATRIA LNATRIX LHATRIA LMATRIA LHATRIA LHATRIA LHATRIA LNATRIA LNATRIA LHATRIA LRATRIA LNATRIA LHATRIA LRATRIA LHATRII LRATRIK LRATRIA LNATRIA, LHATRIX LHATRIA LRATRIx LRATRIx LHATRIx LHATRIA LHATRIx LHATRIA JULY“ JULY“ LNATRIA LMATRIX LRATRIx LHATRIA LRATRIx LHATRIA LRATRIA INATRIA LNATRIA LNATRIA LNATRIA LMATRIA LRATRIX LHATRIA LHATRIA LHATRIA LRATRIA LuaTRIA LHATRIx LHATRIx LMATRIA LHATRIA JULY“ Utd (TCI-GD‘IOIFCTEHV I..- m." ..m-u 3 co. 373 3 C C C 777 g .4 202 9 20) I55 0°“:TART ... FCIBOOODR'OD. .9 UR=UIFC3CA 126 CONTINUE " INITIALIZATION cf VARIABLES " UU 33 I=IOII OO 33 J=IOJJ U0 33 K3I0JJ IFT(IoEUoII.O“.TJoE00II00R.TKoEUoIII DO TO 373 IFT(I.Eu.IT).uR.TJ.Lu.JJ).uR.TK.Eu.JJ)) 00 TO 373 CHIIIOJOKI3TENPII0J0KI/P5ITI0J0KI 00 T0 33 CHIII0J0KI=Uo0 00 T0 33 CONTINUE OVERRELAKATIONo 0.0 9.. PRINT 7770EPSHA10ITMAA FORMATTIHOO. LNATNIA.9EI0070IIUI ITER=0 EP53000 ITERIITEROI SUM3°oo UO b I‘ZONI UO 6 J=20NJ DU 6 K320NJ 3T0NE3 CHIIIOJOKI CHIIIOJOKIaLHITI0J0KI’CA’T ICHITI'I9J0KI'TVA'UELATPHIOIOJQKI-UELAITNETficl.J.K)-1H£T‘([.J.K)O dtl.¢°0.5'DELxTu£10IoJoxTI) J'CHITI-IoJokI'IUELAITHETAOIOJOKI-VA h-THETATI030KT'(I.°0.5'U£LATUZI0I0J0KIII SOCHITI0JOI0KT'TFA'TVA'UELVTPHI0I0J0kI-ULLYTTHETA0I0J0KII b-THETATI0J0KT°TOAHMAOIE'UELTTUlloI0J0KIII 'OELXIPHI0I0J0KI 70CHITI0J-I0KI'TFA'IOELYTTHLTA0I0J0KI-VA' UELYTPHI0I0J0KII D'THETATI0J0KI'(GANHA‘FE.OLLYTDII0I0J0KIII 9°CHITI0J0ROII'IFA‘TVA‘UELITPHIOI0J0KI'OELITTHLTA0IOJ0KII A-THETATI0JoKI'TGAHHAOFE'OELZTOII0I0J0k)I) HOCHIII0J0K-II’TFA‘TDELITTHETA0I0J0KI-VA. OELZTPHI0I0J0KII C‘THETAIIoJoKI.IGAHNA'TE.UELZIUZIOIIJOKIII O-CHITI0J0KI.P5ITI0J0KI0 8TEHPTI0J0KII EPS'EPSOAHSTCHITI0J0KI-STORLI 5OM8§UM0ARSTCHITI0J0KII CONTINUE IFT3UH.E0.0.0I EPSsEPS/sun PRINT ZotPSCSUH IFTEP5.LE.EPSHAAT GO TO 9 IFIITER-ITNAXI “0605 PRINT ZJZOITERQEPSQSUH FORMATTIH00'NO CONVERGANCE Thu.l0'Anu >un:¢.u|l.«I 00 TU I55 , PRINT 20I01TER0EPS0SUR FORMATTIH000CUNVEROENCE IEPS ='0EH.I0 ' ANU >UR CONTINUE SUH'IOOC‘IO AFTER .0Ih0'ITERATIONS HITH £PS:O0 hAS NEACHEU AFTER 0.I«.OTT£RATIONS uTTR 3 .OUIIO“7 C ..‘OHTAIN CONDUCTIVITV0" II ISIOII J'IOJJ KI|0JJ 00 II 00 II UO II N... IFTTHETATI0J0KI00E0TAI RSRATIU PSITIOJOKIIPIIIIOJ0KI-CHITIOJ0KI/IU PlITI0J0KI=THETATI0J0KI'PSITI0J0KI PbIIIoJoKIIR'PSIIleOKI CONTINUE LRATRIx LHATRIA LHATRIA LRATRIA LNATRIA LHATRIX LHATRIA LHATRIA LHATRIA LRATRIx LHATRIA LRATRIA LHATRIA LHATRIA LHATRIA LRATRIA LNATRIK LNATRIX LHATRIA JULY» LNATRIx LNATRIA LMATRIA LRATRIA LMATRIA LRATRII LNATRIA LRATRIx LRATRIA LRATRIx LRATRIA LNATRIA LRATRIA LNATRIA LNATRIA LHATRIA LRATRIA LNATRIX LHATRIA LRATRIA LRATRIx LRATRIx LRATRIA LNATRIX LRATRIA LNATRIx LRATRIA LNATRIX LRATRII LRATRIA LHATRIX LRATRTX LRATRTA LRATRIx LHATRIL LRATRIX LNATRIA LNATRIX LHATRIA LNATRIA LHATRIA LNATRIA LRATRIx LHATRIA LHATRIK LMATRIA IIO III II2 II3 II“ IIS II6 I17 IIU II9 I20 IZI I22 I23 I26 I25 I20 I21 nnnr IZ “0P5! cumTAle SIUMA 0THLTA CONTAINS NORMALIlLD Tgup,0o ‘.'COMPUTE CURRENT IN OCURRENT OUT FOR LO0PINO" ".OZI=OIFF CURRENT DENSITY 0P2I'CDNO. CURRENT "'OZI3OIFF CURRENT DENSITY OPZI‘CONOo CURRENT DO I8 I=IOII DO I2 J=I9JJ DU IZ K=I0JJ TEHPTI0J0KI=-PSIII0J0KI'DELAIPHI0I0J0KI/A DZITI0J0KI=DELXTP11010J0KI/IA'VTHI CONTINUE FORMATIIHO0IITEII040IAII RETURN END DENSITY... DENSITY IN A DIR-0' LRATRIA LMATRIx LHATRIA LHATRIX LHATRIX LMATRIx LNATRIX LMATRIA LRATRIx LMATRIX LMATRIA LMATRIx LHATRIX I2” -IZ9 I30 131 I32 I33 I34 I35 I36 I37 I38 I39 IRO 9 co. 5 SURRUUTINE Rho COMMON/CON/OAHNNOONONOG‘NZOGANNUObflufll COMMON/A/IIOJJ0NIONJ COMMON IAIII0I|0III0lDIII0II0III0 PIIIOIIOIII COMMON/VOLTSIPHITIIOIIOIII COMMON/SID]. . CUNNON/T/IHOOT‘OTD CONHON/V/VROVTHOVO 'COHHON/FACTORS/FAOFROFCOFO COMMON/TIRES/TAUU0OI0TAUF002 C0NR0N/LUOP/ROUNT0TIRE FORRATIIR00QIIO0FI003I '0 P CONTAINS STORA-l-SDUARED DO 5 I3I0II DO 5 J3I0JJ DO 5 K'IOJJ. PIT0J0KI'FO'IDELKTPHI0I0J0KI"2 ' IODELYIPRI0I0J0KI"2 0DLL£IPHI0I0J0KI.‘2I 2.1ATIOJINI ' CONTINUE RETURN END RNO RNO RN” RNO RNO RHO RHO RRO RRO RRO RNO RNO RNO- RHO RHO RNO “HO RRO RND RHD OIICIOWIIPHHU ... t”-_fi vi can 129 SURROUTINE BYPASS COMMON/B/CA0OA0EOA0EDB COMMON/AIRS COMMON/CON/OAHHA0ARDA0UAR20GAHHB0BANR2 C0MH0N/A/II0JJ0NI0NJ COHHON/TdO/PZI(IIOIIOIII0TENPII30II0III CORHON/OEbK/THETATII0II0III COMMON/SIUIA LUNHON/VAR/VART0VARCoVARVQVAREQVARPQVAROQVARIoVARFoINTAQINTeoINTC. IINTO COMMON (ATII0II0III0IBIII0II0III0ZCIIIOIIOIII COMMON/VOLTS/PRITII0II0III CUMMUN/I/THO0TA0TD COHHUN/V/VR0VTH0VO COMMON/ATRA/FE0RI0TAUU EOUIVALLNCLI1A0PSII DIMENSION PSITII0II0III '0 CALLED FRoM BYPASS VART=IJHItHP AT VARC=IOHCHARCE AT VARV=IUMVOLTAOE AT VARE=ToHCUNUTY AT VARleonPOBER AT YARD:IOHOONORS AT VAR|=IOMCOND J AT VARF=IOHDIFF 3 AT INTA=IOHD£CHI COEF INTB=IOHDLHI COEFF INTC=IORCHT COEFFT INTD=TOHCONST TERM DO 2v T=IOII DO 20 J=IOJJ DO 20 K=I0JJ (CII0J0KI3000 THETATI0J0KI3THO C" I“ IS NET CHARhE DENSITY0ZA IS IONIZED DONOR DENSITY 0. 9° 20 lBTI0JoKI3000 PlIII0J0KI=A [A II0J0RI=A PHIII0J0KI3III-II/I200 UO 6 I=IOII DO 6 J=IOJJ O0 6 K=I9JJ TEHPTI0J0KI3-PSITI0J0KI'OELAIPHIOI9JOKIIA CALL CIRCUIT RETURN ENI) BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS ‘ BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS BYPASS oegomou~ SUBRUUTINE CIRCUIT CIRCUIT 2 C0BR0N/CUN/6ARHA0AR0A0UAHZ0GARHB0UARR2 CIRCUIT 3 COMMON/A/II0JJ0NI0NJ CIRCUIT b COMMON/V/VR0VTH0YU CIRCUIT S COMMON/SlU/A CIRCUIT b COMMON/T/THU0TA0TO CIRCUIT T CONNOR/VULTS/PRITII01I0IIT ° CIRCUIT I COMMON/BEGA/TBETATIloIloIII CIRCUIT . 9 COMMON lATII0II0III0(CIII0II0III01AIIIOII0III CIRCUIT I0 COMMON/LIN!TSICRAA.VNAA0DELY0NA CIRCUIT II COMMON/TuOIPzTTII0II0IIT0T£RPTII0II0IIT CIRCUIT I2 COMMON/CULSIFTIII0FTIIII0UIIII CIRCUIT I3 COHHON/H/CA0UA0EUA0LUB CIRCUIT IA ALso suNaNI-z CIRCUIT IS, DO 6 I3I0II0NI CIRCUIT I6 AL-KLoI CIRCUIT I1 c OBTAIN THE CURRENT AT TI030KT 0000‘ CIRCUIT In no 2T JaloJJ CIRCUIT I9 00 22 AaloJJ CIRCUIT 20 22 FTATszATI.J.AT0T£RPTI.J.KT CIRCUIT 2T c 0' 00 1A IS UTTF CURRENT DENSITY 0' 0' CIRCUIT 22 00 0' PZI IS CUNU CURRENT DENSITY '0 0' CIRCUIT 23 C" " F IS INTEGRATED TRON I0JJ FOR EACN J 9' '0 0' CIRCUIT 2b 2I rrTJTISInPSTrI CIRCUIT 25 ' FFFSSIRPSTFFI'VR CIRCUIT 26 IFTKLoEU0II TI-I.0-FFF. CIRCUIT 21 l'IKLOEUOZ’ PZ3I90‘EFP CIRCUIT 2B 6 CONTINUE CIRCUIT 29 FAso.S-TF|0F2T CIRCUIT 30 PRINT chaorz CIRCUIT JI PRINT 2 CIRCUIT 32 PRINT 30CNAA0YHAA CIRCUIT 33 ITTFA.GE.VBAAI FAIYBAA CIRCUIT 36 FFISI.0-fA CIRCUIT JS TTTTFT.UL.CNAAT FrF-CRAA CIRCUIT 36 FAaT. o-TFT CIRCUIT 31 2 FORMATIIR00O RAA LIMITS ARE OI CIRCUIT 30 PRINT 3. TTT0TA CIRCUIT 39 5 CONTINUE ' CIRCUIT A0 3 F0RHATTIB00IOA0'I- 00EI2.S0I0A0'VI '0EI2oST CIRCUIT AI 00 ST J'IOJJ CIRCUIT 02 00 5| K=I0JJ CIRCUIT 03 SI PHIII0J0KI8FA CIRCUIT A“ RETURN CIRCUIT «5 END CIRCUIT to 130 C on co. co. c no (I. c I. 131 FUNCTION SINPST'I ..P 86 CARNAHAN " HII TIUHIZ BaII A-I '0 CURRUN IN RORK ARRAY DIMENSION TTIII COMMON/A/IIOJJQNICNJ INITI‘LICE PARAMETCRS 9. .. 5UHEN03000 t SUMNIU‘Ooo " EVACUATE SUHENU AND SUMMID 9'. DO I H=I05 ' Rae'M-I bUHENUSbUHtNUOFIKI bUMHIU=SUHHID°FIROIT CONT INUL '0 R£TUNN ESTIMATED VALUE OF INTEGRAL "' 5|HPS=(2.0959MENU0R.0'SUHHID-FTIIOFIIITT/Joo NITUHN [NU SIHPS SIRPS SIHPS SINPS SIHPS SIMPS SIMPS SIHPS bIHPS >IHPS SIHP5 SIMPS SIHPS SIMPS 5IHP§ SIHPS SIHPS SIHPS OCI~iOIfll>hHU 823 lb 12 I3 I“ 'IFIJ.£O.JJI 132 SUHNOUTINE PSSNIITERI COMMON lAIIIoIIoII)oAHTIIIoIIoII)OZCIIIQIIIIII COMMON/VOLTS/AROIIIoIIoII) CONNON/CON/OAMMA.AMUAIOAHZoGANNBoOANNZ CONNON/VAN/VARToVARCoVANVoVAREvVAHPIVARDoVARIIVAN?oINTAoINTBoINTCo IINTU CUNNON/A/IIoJJcNIoNI COMMON/OAHS/OAMU COMMON/VIVRIVTHoVO CONNON/LOOP/KOUNToTINt COMMUN/AN/ITMAXQEPSMAXOUI COMMON/TITNOoTAoIU CORRUN/SIU/A EOUIVALENLETNIoNJ) FORNATIIHO.‘ Eps: IT£R=0 EPS=C.00 5UH=300 ITtR=ITEN0I UO'S J=IoJJ LI=J-I L2=JOI U0 5 K=IOJJ HI=K-I M2=KOI IFIJ.EO.II .obII.5o. EPSNAX. .OGIIQSO. TA. 'IGIIoSI GO TO I2 60 TO I3 IFIK.tOoI) 60 TO I“ IFIK.EO.JJI GO TO I5 00 T0 I7 LI32 00 TO Ib L2=NI GO TO I6 NI=Z GO TO I7 H2=NI CONTINUE UO 5 ISZONI TEMP:AR6II9J0KI ARoIIoJoKI=ANbIIkoKI'IIoO-IIIOUI' 0IANSIIOIOJIKIOAHOII-IOJOKIqGANNA. 9IANbIIoLIoKTOARbIIoLZoKTO 9AN6II9JOMIIONHOIIIJQMZIIO 9NH7TIOJ0KIIIGAND SUN=SUM0AHSIARSII0JIKII EPS=EPSOAHSIAH6IIoJoKI-TENPI EPS=EPS/SUH ' IF IEPSoLfiofiPSNAAI GO TO 6 STOP ITER IF CONVERCENT UR EAESSIVE ITER IFIITER‘ITHAXI “0“98 . NHITEIOIOZOZIITEHoEPSoTINE FORMATIIHOo'NO CONVERthCE AFTCR'OIQOP ITERATIONS UITH EP58.0 IENZIQ' AT TIME =PoGIIo4I- DO TO I55 NHITtIfiIoZOII ITENotPSoTINE FONMATIIHOo'CUNVLRUtHLE HAS HEtN "CACHED AFTER'OIQO'ITERATIONS lulIH LPs=0o IL’IOIO. AI TIME 3.0TJIII‘OI LIHITIiHIL ' HA=8 I)” I "oI—IQ’N i, . I. I . ' A , V--~o-—‘————-'— —. POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON -POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON 'POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON T"II‘.'--MI oqumauw 133 I)“ I TO I‘d. O“. 00 I70 J'ZoNJ ALL RIURT NANO PTS ARE CORRI DEFINED " ANOIIOJIKLI'ANOIIOJONT ARoIIoRLcJIcARCTIoRoJI CONTINUE R EDGES OF THE CONE REMAIN .' U0 ITI I'ZONI APNTIOIoItitARbIIoIoNI ANbIIoJJoNLISANbIIoJJoNT ANOIIOKLCII'ANSIIONOII APbIIoNLoJJItANbTIoNoNJI CONTINUE CONTINUE PRINT OZJOEPSQEPSNAAOSUN RETURN ENU POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON POISSON' 61 69 7. TI 72 13 TA 75 To 77 70 79 O. GI 82 134 SUNNOOTINE TRIOAG CONNON/AIIIOJJONIONJ CONNON/L/LAON CONNON/COLS/OHIIIToAAIIIToGOTIII CONNON/OORK/ANITIIToANJIIIIOAPSIIIIOANOIIII C SUHNOUTINE TNUG EON.SOLVING LINEAR SIRULTANEDUS EUNS. C'°"APP. NON. METHODS --PP-D . IFsLAOI ‘IRIILAI-RRIIAI - ARSILATaARJILAT/ARIILAT 00 TOT I'IEON ARIIII'OOIII-AAIII'UUII-II/ARITI-IT IfIARSIARIIIII.LE.I.0E-2JT ARIIIT-I.OE-2J IOI ARSIIIaTARJIII-AAIII'ARSII-ITIIARIIII ARJINI-ARSINI LAST-N-LA , no In: R-IoLAST c TRIOAO EON SOLVING SIRULTANEDUS E I'N-K I02 ARJIIT-ARSTII-GUTII'ARJIIOIT/ARIIII RETURN END TNIOIO TPIOAO TRIOAO TNIOAO TNTOAO TRIDAC TNIOAO TNIOAO TNIOAG TNIOAO TNIOAO TRIOAO TPIDAG TNIOAO TNIOAO TNIDAO TNIOAO TNIDAG TNIOAO TNIOAO TNIOAb TNIOAO 0.60.09” ’59” “_— OEIMR “EUR 222 333 SS 53 33 35 JT Ii QJ 65 135 IONCTION UELITUoIoJoIT UTNENSION ”IIIOIIOIIT CONNOR/A/IIIJJONIONJ CONNON/CON/OANNAOANDAOOANZOOANNOQOANNZ CONNON/SIU/A CUNNUN/V/VNOVTHOVO CONNOR/ATRA/FE.RIoIAUR ITIlotOoIT DO TO III IT IIoEOoIIT DO TO 222 OILASIUIIoIoJoKT-UII'IOJOKTT/ZoO 00 T0 JJJ DEL‘.-905.TJO°.UTSUJUKT-“OO.USO.NUJOKP.u‘NOZUJOK’. UU T0 333 OEan o.S-IJ.o-UII.J.RI-~.00UII-IoJoKIoUII-2.J.KII CONTINUE RETURN ENTRY DELV ITIJ.LO.IT 00 TO SS IEIJ.LO.JJI 00 TO 55 ULLX:EE°TUIIOJOIONT-UTIOJ'IONTT OO TO SJ UTLK‘OAO CONTINUE RE TURN ENTRY DELI ITTKoTUoIT 00 TU 55 ITIK.LO.JJT GO TO SS OtlefE'IUIIOJQNOII-UIIOJOR-ITI RETURN ENTRY DAV! IFTIAKGOIT 00 T0 N3 IEIIoEUoIIT DO TO $5 OELASUTIOIOJOKTOUII-IOJOKT RETURN ENTRY OAVT IEIJ.EO.IT GO TO 33 IFIJ.EO.JJT GO TO 35 OELAabANNA'TUIIoJ'IoKTOUIIoJ.IoKTT UO TO 37 ~ OELAIGARHA'TUIIOJOIOKTOUIIOJQNTT 00 TO 37 UELKSOARNA'TUIIOJ'IOKTOUIIOJOKTT CONTINUE RETURN ENTRY UAVl IFIK.E0.IT GO TO I2 IFIR.EO.JJI 60 TO It ULLXSOARRA'TUIIOJOKOIIOUIIOJQN'ITT RETURN UELA=OANHA'TUIIoJoKOIIOUTIoJCKTT RETURN ULLA=OANHA'IUTIoJoI-ITOUTleoKTT RETURN UELIIJoO'UTIOJOKT'doO'UIIOIOJoKT0UTI°ZQJ9KT RETURN . ‘ OELA:J.O'UTIIJ9KI'ZIO.UII-IOJOKIOUTI'ZQJQKT RETURN LNO UELA UELA UCLA UCLA UELA UCLA UELA OLLA OELA OELA UELA UELA UELA OELA UELA OELA OLLI OELA OLLA OLLA OELA OELA UCLA OLLA DELA UELA UCLA OLLA OELA OELI OELA UELA UELA DELA UELA DLLA OELA OLLA DELA UELA DLLA DELA UELA UELA OELA UELA OELA UCLA UELA DLLA UELA DLLA OELA IIELA ULLA DLLA UELA UELA OCQOGOUN PAP APPENDIX C THE DONOR IONIZATION EQUATION The n type semiconductor modelled in Chapter 2 is assumed to be only partially ionized at any temperature T in the operating tem- perature range. In this section, an expression for the ionized donor density at any specified temperature is derived using Fermi statistics and the sample band model of the semiconductor outlined in Chapter 2. The n type semiconductor chosen has been assumed to have neg- ligible acceptor attoms, i.e., Na«é<.Nd Further, the semiconductor has been assumed to be a nondegenerate semi-_ conductor. Those donors are incompletely ionized in the temperature range of operation. Thus, at T a OK, the valence band is full, and the conduction band is empty. In the operating temperature range, some of the electrons from the donor states are thermally excited into the conduc- tion band, leaving behind ionized donors. Since, as shown in Fig. 2.1a, the valence level lies much further away below the conduction band than the donor level, only a negligible number of electrons acquire sufficient energy to be excited from the valence band to the conduction band. Thus, if p(Ec) be the probability of finding an electron in the conduction band, then, applying Fermi statistics use) .- (1 + exp(Ec - 1~:f)/k1')"1 ((3.1) 137 Also, if n be defined as the number of ionized donors per unit volume and p(E ) be the probability of occupation of a donor state, then, n; a Nd(l - p(Ed)) (0.2) BLt, p(Ed), the probability of occupation of a donor state can be obtained by applying Fermi statistics: p(Ed) = (1 + exp(Ed - Ef)/kT)‘1 (0.3) From Fig. 2.1, E is greater than Ed, measured from the valence band f edge. Also, it has been assumed that the operating temperature range is low enough such that only partial ionization takes place and so kT Ed’ E The quantity (Ed - Ec)/kT is thus negligible compared to f. unity. Under equilibrium conditions, charge neutrality is maintained everywhere within the sample. Since both acceptor density and hole concentrations have been neglected in the model, then, in equilibrium, _+ n — nd Ndexp((bf - EC)/kT) - Ndexpuhd - bf)/kT) (0.4) Hence, hf = (BC + Ed)/2 Note that this is true only under equilibrium conditions. Under operating conditions, it is only approximately true. Eq. (0.2), therefore, reduces to the simple expression n; a Ndexp((Ed - Ef)/kT (6.5) Similarly, using Eq. (g.1), the number of electrons per unit volume under equilibrium conditions is given by n a Ndp(EC) n Nd (Gob) 1 + exp((Ef - HOS/kt) Also, in the selected temperature range, RT (‘1 EC - Ef 138 Hence, exp((Ec - Ef)/kT)>> 1 and so, n - Ndexp((Ef - Ec)/kT) (6.7) Defining E, the activation energy, by the relation Ea - EC " Ed (608) and substituting back in the relation Eq. (6.5) for Nd' 4. nd a Ndexp(Ea/kT) (6.9) This is the relation stated in Eq. (2.2.1). There is an approxi- mation involved since equilibrium conditions are utilized in the deri- vation. Since considerable diffusion effects are expected to be encoun- tered, the equilibrium expression for n cannot be used. However, the ionized donors are stationary, and the equilibrium expressions are here as an approximation. . _’ “- e“—Ia917 APPENDIX H THE GENERALIZED CONTINUITY EQUATION The various dependent variables, like , the charge density, are functions of space and time and often occur as total differentials. In this section, the total differential is expressed in terms of par- tial derivatives and the divergence operator, both of which can be di- rectly evaluated by finite difference techniques. Consider a real function defined in a region R. The function represents any real physical pr0perty. It is, therefore, finite, bounded, and continuously differentiable, and a function of Space and time. In cartesian co-ordinates, f 2 f(x, y, z, t) Taking the total derivative of this function with reSpect to time, df - ggh+ afdx +‘3fdz + afdz 3'1:- at axdt aydt azdt £133, dz, and 13 represent the velocities with which the property in x, dt dt dt y, and 2 directions reSpectively. Writing vx, vy, and v2 for the re- Spective velocities, £-2£+?£.V +2;v,.+?_f.v dt 3: ex x ay )2 2 v , v x , vz are the components of the velocity vector 3; Consider y the vector Emotion f3. Since V is the velocity associated with the property f, the quantity fvrrepresents the current or flux density 139 FT.“ 140 associated with the property. Performing the divergence operation on th this function, V. (f?)=fv.v'+v.Vf But V. V I 0. Also, v .‘7f I v f + v f + v f x—— y—— z— x y 2 Hence, (it at “if-+VOJf 3t where 3f is the current density associated with the property f. The above equation is utilized in arriving at the continuity equations for charge and energy flow in Chapters 2 and 3 reSpectively. BI BLIOGRAPHY BI BLIOGRAPHY Ila:- 1. 3. 4. 7. 8. 10. ll. 12. 13. A. A. (f. o E. P. BIBLIOGRAPHY D. Pearson, ”Chemical, Physical, and Electrical Properties of Some Unusual Inorganic Glasses: Part 2", Advances in Glass Technology, New York, Plenum Press, (1962) 357. D. Pearson, "Characteristics of Semiconducting Glass Switching/ hemory Diodes", IBN J Res Develop, 13 (1969) 510. R. Uvshinsky, "Amorphous Semiconductors", Sci J, 5A32 (Aug 1969) new 73. N. Boer, G. Dohler, S. K. Ovshinsky, "Time Delay for Reversible Electric Switching in Semiconducting Glasses”, J Non-Cryst r— Solids, 4 (1970) 573. A. Davis, R. F. Shaw, "Characteristic Phenomena in Amorphous Semiconductors“, J Non-Cryst Solids, 2 (1970) 107. J. Walsh et a1, ”Experimental Results in Amorphous Semicon- ductor Switching Behavior", J Non-Cryst Solids, 2 (1970) 107. Editorial Staff, ”Glassy Semiconductors Show Switching and Memory In L\ . D. k. M. H. N. Effects", Phys Today, 2231 (1969) 62. R. Shanks, "Ovonic Threshold Switching Characteristics", J Non- Cryst Solids, 2 (1970) 99. F. Neirauch, "Threshold Switching and Thermal Filaments in Amor- phous Semiconductors", Appl Phys Letters, 1635 (1970) 212. N. Boer, 3. R. Ovshinsky, "Electrothermal Initiation of an Elec- tronic Switching Mechanism in semiconducting Glasses", J Appl Phys, 4106 (1970) 2675. H. Cohen, H. Fritzsche, 3. R. Ovabinsky, “Simple Band Model for Amorphous Semiconducting Alloys", Phys Rev Letters, 22:20 (1969) 1065. Fritzsche, S. R. Ovshinsky, "Electronic Conduction in Amorphous Semiconductors and the Physics of the Switching Phenomena”, J Non-Cryst Solids, 2 (1970) 393. F. Mott, "Conduction in Non-Crystalline Systems: Part I - Localized Electronic States in Disordered Systems”, Phil Mag, 173150 (Jun 1968) 1259. 141 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. N. N. A. E. A. A. G. R. D. F. G. J. R. J. H. D. B. K. D. 142 F. Mott, "Conduction in Non-Crystalline Materials; Part III - Localized States in a Pseudogap and Near Extremities of Conduc- tion and Valence Bands", Phil Mag, 198160 (1969) 835. F. Mott, "Conduction and Switching in Non-Crystalline Materials”, Contemp Phys, 1032 (1969) 125. I. Gubanov, ”Qggntum Theory of Amorphous Semiconductors", Trans- lated by A. Tybulewicz, New York: Consultants Bureau, 1965. Harmik, A. Lapek, "Current Filaments in Acoustically Amplifying CdS", Solid State Commun, 8 (1970) 1287. M. Barnett, H. A. Jensen, "Observation of Current Filaments in Semi-Insulating GaAs”, Appl Phys Letters, 12810 (1968) 341. P. Ferro, S. K. Ghandi, ”Observations of Current Filaments in Chromium-Doped GaAs", Appl Phys Letters, 1685 (1970) 196. _OMA . - ~—.x .— Nimtz, K. Seeger, "Current Controlled Negative Differential Re- sistivity of p - Type Tellurium", Appl Phys Letters, 1431 (1969) 19. A. Smith, "Semiconductors", Cambridge (Eng), Univ Press, 1959. Adler, J. Feinleib, H. Brooks, R. Paul, “Semiconductor-to-Hetal Transitions in Transition Metal Compounds", Phys Rev, 155 (1967) 851. J. Morin, J. P. Maita, Phys Rev, 96 (Oct 1954) 28. H. Ludwig, R. L. Halters, Phys Rev, 101 (Mar 1956) 1699. Moss, "Photoconductivity in the Elements“, Butterworth, 1952. A. Smith, Physics, 20 (195a) 910. Feinleib, N. Paul, “Semiconductor-to-Metal Transition in V203", Phys Rev, 155 (1967) 841. J. Stacker, "Phenomenology of Switching and Memory Effects in Semiconducting Chalcogenide Glasses", J Non-Cryst Solids, 2 (1970) 558. _ L. Nelson, “Ovonic Device Applications”, J Non-Cryst Solids, 2 (1970) 371. K. Ridley, "Specific Negative Resistance in Solids”, Proc Phys Soc, 82 (1963) 954. H. Boer, G. Dohler, ”Temperature Distribution and Its Kinetics in a Semiconducting Sandwich", Phys Stat Sol, 36 (1969) 679. Adler, H. Brooks, "Theory of Semiconductor-to-Netal Transitions”, Phys Rev, 155 (1967) 826. 33. 34. 35. 36. 37. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48. H. A. D. C. R. P. C. B. G. A. J. B. 143 Futaki, "A New Type Semiconductor (Critical Temperature Resis- tor)", Jap J of Appl Phys, 13 (1965) 34. N. Berglund, R. H. Walden, "A Thin Film Inductance Using Thermal Filaments”, IEEE Trans on Electron Devices, ED-l7 (1970) 137. C. Warren, "Switching Mechanisms in Chalcogenide Glasses", Elec- tronics Letter, 5 (1969) 461. Adler, T. Kaplan, "Electrothermal Switching in Amorphous Semi- conductors", 8-10 (1972) 544. Kittel, "Introduction to Solid State_Physics", John Niley & Sons, 1968. W. P. King, "Fundamental Electromagnetic Theoryfi, New York: Fer Dover Publications (1963). 5. Carslow, H. C. Jaeger, "Heat Conduction in Solids", Oxford Clarendon Press, (1959). t: C. Duchene, M. M. Terruillon, N. Pailly, G. Adams, "initiation of Switching in V0 Coplanar Devices”, IEEE Trans on Electron Devices, ED-71 (1971) 1151. J. Yu, "Electrical Switchinggfiffects in Vanadium Oxide Complexgsfl, PhD Thesis, USU 1972. D. Fisher, R. J. Yu, "Inductance Properties of Oxidized Vanadium Foils", IEEE Trans on Electron Devices, ED-16 (1971) 975. N. Berglund, "Thermal Filaments in Vanadium Dioxide", IEES Trans on Electron Devices, ED-l6 (1967) 432. Car ahan, H. A. Luther, 3. 0. Wilkes, ”Applied Numerical Methods”, John Wiley & Sons, (1969). E. Forsythe, w. R. Vasow, "Finite Difference Methods for Partial Differential Equations", John Wiley & Sons, (1960). Nussbaum, "Semiconductor Device Physics“, Englewood Cliffs, Prentice Hall, (1962). Douglas, "Alternating Direction Methods for Three Space Variables”, Numerische Nathematik, 4 (1962) 41. S. Tanenbaum, "Plasma Physics", McGraw Hill, (1967). TV LIBR "'TlliliffllllififllfllfllIlllljijflflflflmuzllfllwlill“