高頻半導體元件電路時域模擬之研究

全文

(1)國立交通大學電子工程學系電子研究所. 博士論文. 高頻半導體元件電路時域模擬之研究. A Time Domain Approach to High-Frequency Circuit Simulation of Semiconductor Devices. 研究生：黃堃宇指導教授：李建平中華民國九十四年九月.

(2) 高頻半導體元件電路時域模擬之研究 A Time Domain Approach to High-Frequency Circuit Simulation of Semiconductor Devices 研究生：黃堃宇. Student: Kuen-Yu Huang. 指導教授：李建平博士. Advisor: Dr. Chien-Ping Lee. 國立交通大學電子工程學系電子研究所. 博士論文 A Dissertation Submitted to Institute of electronics College of Electrical Engineering and Computer Science National Chiao Tung University in partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy in Electronics Engineering June 2005 Hsinchu, Taiwan, Republic of China. 中華民國九十四年九月 ii.

(3) 誌謝首先要感謝我的指導老師李建平教授，在過去八年的碩博士研習中與以無私的指導與幫助，老師不僅在學術研究上給予我教導指正，更是我在為人處世上的學習的目標，在此再向我敬重的老師表達誠摯的謝意。在這段求學的歷程中有許多學長同學及學弟給予許多助力，與你們的相互討論及研究，使我的論文有了新的想法與更完整的研究，特別是黃忠諤博士及廖志豪學弟，在實驗量測與元件特性上，常常向兩位請教及討論，真心感謝兩位的協助。另外，李義明博士在數學理論上對我的論文可說是有關鍵性的影響，在此也要特別說明及表達感謝之意。最後，我將我的論文獻給我的家人，有了與你們共同扶持的生活，我的生命真是豐富多彩且充滿樂趣，感謝你們，謝謝！. iii.

(4) 高頻半導體元件電路時域模擬之研究. 研究生：黃堃宇. 指導教授：李建平博士. 國立交通大學電子工程學系電子研究所. 摘. 要. 為了模擬射頻穩態電路問題 (意即：高頻週期性電路) ，在本論文中，我們發展了一個新的時域數值方法，用來解如雙調交互調頻失真 (two-tone intermodulation distortion) 之電路問題。傳統類似 SPICE 之時域電路模擬器中的暫態分析 (transient analysis) 並不適用於解高頻穩態電路問題。由於計算時域解高頻穩態電路問題的嚴苛限制條件，驅使我們發展屬於自己的時域非線性電路模擬器。此新的電路模擬器必須符合穩定、有效率的條件，並且可以處理如：半導體元件等效電路模型之強烈非線性電路問題。為此，我們成功的結合了波形分散法 (waveform relaxation method) 、單調疊代法 iv.

(5) (monotone iterative method) 與 Runge-Kutta 法 (一種常微分方程式積分法 ) ，用之於解時域電路常微分方程式 (ordinary differential equations) ，並且可以符合上述模擬高頻穩態電路之要求。此電路模擬器所用之數值方法已被證明可以收斂並於本論文中揭示其收斂特性曲線。為了要從模擬出來的時域數值資料之擷取出有用的頻域資訊，我們也做了如離散 Fourier 轉換法 (discrete Fourier transform) 等之後續分析。在本論文中，我們使用自己發展出來的數值解題法與 Gummel-Poon 大訊號電路模型來模擬異質雙接面電晶體 (heterojunction bipolar transistors, HBTs) 。于論文中，我們討論了直流電路計算、時域模擬、頻域分析與交互調頻失真特性分析 (intermodulation distortion characterization) 之結果。我們也進一步對不同的電路模擬器 (包括了我們發展的模擬器、HSPICE 與 ADS) 之模擬結果與量測數據做了相互的比較，以說明所發展的數值模擬方法是準確且有效率的。對於操作在高偏壓與高功率的 HBT 而言，熱效應主導了元件之行為特徵。為此，我們將電熱交互作用的方程式包含到所發展的數值計算方法之中，如此可以進一步增進數值模擬之真實性。在論文的內文中詳盡討論了考慮與不考慮熱效應的電路模擬之不同處。我們. v.

(6) 也另行模擬了一個多指 (multi-finger) HBT 上所發生之熱偶合效應 (thermal coupling effect) ，與此效應對 HBT 之直流、射頻功率與交互調頻失真等多項特性之影響。在附錄中，我們介紹了金氧半導體場效電晶體 (metal-oxide semiconductor field-effect transistors, MOSFETs) 之 EPFL-EKV 大訊號電路模型，並且用所發展的模擬器模擬了相關的電路。此外，我們也進一步的闡述了在本論文中所提出之數值方法的收斂特性。如本論文所言，我們所發展之數值計算方法不但可以用來解時域電路之非線性常微分方程，也可以推廣應用到包含更多數量及更多種類的半導體元件之高頻電路模擬上。. vi.

(7) A Time Domain Approach to High-Frequency Circuit Simulation of Semiconductor Devices. Student: Kuen-Yu Huang. Advisor: Dr. Chien-Ping Lee. Department of Electronics Engineering and Institute of Electronics, National Chiao Tung University. ABSTRACT In order to solve the radio-frequency (RF) steady-state, that is, high-frequency periodic circuit problem, such as intermodulation distortion, we develop a new numerical solution technique in this dissertation to simulate circuits in time-domain. Traditional transient analysis in SPICE-like time-domain solver is not suitable for RF steady-state solution. The tough criteria for solving steady-state problem in time- domain drive us to build our own time-domain circuit solver. This solver should be stable, efficient, and able to handle strongly nonlinear circuits, for instance, the large-signal models for semiconductor devices. Combining the waveform relaxation (WR) method, monotone iterative (MI) method, and Runge-Kutta (RK) method, we succeeded in solving the circuit ordinary differential equations (ODEs) in time domain. vii.

(8) that satisfies the criteria of the RF steady-state analysis. The convergence of our developed algorithms has been proved and demonstrated in this dissertation. We also perform subsequent analysis, such as the discrete Fourier transform (DFT), to extract the frequency-domain information from simulated time-domain data. In this dissertation, we use our numerical solution techniques to simulate the characteristics of heterojunction bipolar transistors (HBTs) with the Gummel-Poon (GP) model. The results of DC circuit calculation,. time-domain. simulation,. frequency-domain. analysis. and. intermodulation distortion characterization are presented. Furthermore, we compare the results of various simulators (our solver, HSPICE and ADS) with measured data to show the accuracy and efficiency of the developed method. For HBTs under high bias and high power operation, the thermal effects dominate the. behavior. of the. device.. Therefore,. we. include the. electrical-thermal interactive equations in our numerical solution technique to further improve the reality of the simulation. The difference between the simulations with and without thermal effects is well discussed in the context of this dissertation. In additional, we describe the thermal coupling effects of an multi-finger HBT and its influence on DC, RF power and intermodulation distortion characteristics. In the appendixes, the EPFL-EKV compact model of the metal-oxide semiconductor field-effect transistors (MOSFETs) is discussed and the related circuits are simulated by the proposed method and HSPICE. The convergence properties of the utilized algorithms are presented in these appendixes.. viii.

(9) According to the discussion in this dissertation, the developed approach is not only an alternative computational technique for the time-domain solution of nonlinear circuit ODEs but also can be generalized for high-frequency circuits simulation including more and variant kind of semiconductor devices.. ix.

(10) Contents Abstract (in Chinese) . . . . . . . . . . . . . . . . . . . . . . . . . . .. iv. Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1. 2. xv. Introduction. 1. 1.1. Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 2. 1.2. Background . . . . . . . . . . . . . . . . . . . . . . . . . . .. 3. 1.3. Historical Development . . . . . . . . . . . . . . . . . . . . .. 5. 1.4. Organization of the Dissertation . . . . . . . . . . . . . . . .. 6. Mathematical Models. 8. 2.1. Equivalent Circuit and Parameters of GP Model . . . . . . . .. 12. 2.2. Model Equations . . . . . . . . . . . . . . . . . . . . . . . .. 16. 2.3. Formulation of Node Equations . . . . . . . . . . . . . . . . .. 21. 2.3.1. Nodal Equations of GP Model . . . . . . . . . . . . .. 22. 2.3.2. Nodal Equations of Simulated Circuits. 23. x. . . . . . . . ..

(11) CONTENTS. 2.4 3. Summary of This Chapter . . . . . . . . . . . . . . . . . . . .. 28. 3.1. 30. Numerical Algorithms . . . . . . . . . . . . . . . . . . . . . Algorithm for Solving the System of Nonlinear Algebraic Equations . . . . . . . . . . . . . . . . . . . . .. 30. 3.1.2. Time Domain Solution Algorithm . . . . . . . . . . .. 36. 3.1.3. Runge-Kutta Method . . . . . . . . . . . . . . . . . .. 46. 3.1.4. Algorithm with Electrical-Thermal Feedback . . . . .. 49. 3.2. Genetic Algorithm for Parameter Extraction . . . . . . . . . .. 52. 3.3. Further Analysis of Time-Domain Results . . . . . . . . . . .. 61. 3.3.1. Transformation of the Time-Domain Data into FrequencyDomain . . . . . . . . . . . . . . . . . . . . . . . . .. 61. Calculation of Intermodulation Distortion Analysis . .. 65. 3.4. Measurement for Intermodulation Distortion . . . . . . . . . .. 68. 3.5. Summary of This Chapter . . . . . . . . . . . . . . . . . . . .. 70. 3.3.2. 5. 27. Computational Techniques and Characterization Methodology. 3.1.1. 4. xi. DC Simulation and Analysis. 72. 4.1. DC Simulation . . . . . . . . . . . . . . . . . . . . . . . . .. 74. 4.2. Self-Heating Effect . . . . . . . . . . . . . . . . . . . . . . .. 79. 4.3. Thermal Effects of Multi-Finger Device . . . . . . . . . . . .. 83. 4.4. Summary of This Chapter . . . . . . . . . . . . . . . . . . . .. 90. Time Domain Simulation and Frequency Domain Analysis. 91.

(12) CONTENTS. 6. 7. xii. 5.1. Non-Thermal-Effect Simulation . . . . . . . . . . . . . . . .. 5.2. Results with Self-Heating Effect . . . . . . . . . . . . . . . . 100. 5.3. Summary of This Chapter . . . . . . . . . . . . . . . . . . . . 106. Intermodulation Distortion and Power Characteristics. 93. 107. 6.1. Characteristics without Thermal Effects . . . . . . . . . . . . 109. 6.2. Characteristics with Self-Heating Effect . . . . . . . . . . . . 116. 6.3. Characteristics with Thermal Coupling Effect . . . . . . . . . 121. 6.4. Summary of This Chapter . . . . . . . . . . . . . . . . . . . . 129. Summary and Suggestions. 130. 7.1. Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131. 7.2. Suggestions for Further Work . . . . . . . . . . . . . . . . . . 132. References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 Appendix A The MOSFET EKV Model . . . . . . . . . . . . . . . . . . . . . . 151 A.1 Equivalent Circuit . . . . . . . . . . . . . . . . . . . . . . . . 153 A.2 Model Equations . . . . . . . . . . . . . . . . . . . . . . . . 156 Appendix B EKV Model Simulation Results . . . . . . . . . . . . . . . . . . . 162 B.1 Time Domain Simulation Results . . . . . . . . . . . . . . . . 163 B.2 Spectrum Results After FFT . . . . . . . . . . . . . . . . . . 166 B.3 Simulation of a Low Noise Amplifier . . . . . . . . . . . . . . 170.

(13) CONTENTS. xiii. Appendix C Convergence Properties . . . . . . . . . . . . . . . . . . . . . . . . 174 C.1 Convergence of Time-Domain Solution Algorithm . . . . . . 174 C.2 Convergence Properties with Thermal Effect . . . . . . . . . . 175 Appendix D VITA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181 Appendix E Publication List . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183.

(14) List of Tables 2.1. A list of BJT GP model parameters. . . . . . . . . . . . . . .. 13. 2.2. A list of parameters for thermal effects modelling. . . . . . . .. 14. 3.1. A list of parameters to be extracted. . . . . . . . . . . . . . .. 54. 4.1. A set of extracted parameters for Gummel-Poon model used in this study. . . . . . . . . . . . . . . . . . . . . . . . . . . .. 4.2. 74. A set of extracted parameters for the GP electrical and thermal models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 80. A.1 The MOSFET EKV compact model parameters. . . . . . . . . 155. xiv.

(15) List of Figures 2.1. An illustration of Gummel-Poon large-signal equivalent circuit model for the bipolar transistor. . . . . . . . . . . . . . .. 12. 2.2. A testing circuit for an HBT without thermal network. . . . . .. 23. 2.3. A simulation circuit for an HBT with thermal network. . . . .. 24. 2.4. An equivalent circuit of the investigated multi-finger HBT with high frequency input excitation. . . . . . . . . . . . . . . . .. 3.1. A flowchart of the decoupled methodology for the stationary circuit simulation. . . . . . . . . . . . . . . . . . . . . . . . .. 3.2. 25. 35. A flowchart of the proposed time-domain simulation methodology. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 45. 3.3. A illustration of the 4th-order Runge-Kutta method. . . . . . .. 48. 3.4. A flowchart of the proposed electrical-thermal simulation methodology. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 3.5. The proposed weight value W for HBT model parameter extraction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 3.6. 51. 53. The employed flowchart of parameters optimization in our study. 57 xv.

(16) LIST OF FIGURES. 3.7. A comparison of the method with and without the proposed weight value for multiple I-V curves optimization. . . . . . . .. 3.8. 58. A comparison of the method with and without the proposed dynamic mutation rate scheme. . . . . . . . . . . . . . . . . .. 3.9. xvi. 59. An illustration of the migration processes for the I-V curve optimization. . . . . . . . . . . . . . . . . . . . . . . . . . .. 60. 3.10 An illustration of discrete Fourier transformation for a given input time domain data. . . . . . . . . . . . . . . . . . . . . .. 63. 3.11 An illustration of two-tone intermodulation characteristics for a nonlinear two-port network. . . . . . . . . . . . . . . . . . .. 65. 3.12 The ideal output amplitudes at fundamental frequencies and the IM3 products versus the input amplitudes. . . . . . . . . .. 66. 3.13 A setup of on-wafer device testing with harmonic load-pull system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 68. 4.1. A cross section view of an InGaP HBT. . . . . . . . . . . . .. 75. 4.2. The used circuit in the simulation of IC −VBE and IB −VBE curves. 76. 4.3. Comparison of Gummel plot between the simulated results of HSPICE and our proposed method with the same device parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 4.4. The used circuit in the simulation and measurement of I C −VCE DC curves. . . . . . . . . . . . . . . . . . . . . . . . . . . .. 4.5. 76. 77. Comparison between our simulated results and measured results for IC −VCE curves of the InGaP HBT. . . . . . . . . . .. 78.

(17) LIST OF FIGURES. 4.6. An equivalent circuit including thermal network of an HBT for DC simulation. . . . . . . . . . . . . . . . . . . . . . . .. 4.7. xvii. 79. Comparison between our and HSPICE’s simulated results for common-emitter IC −VCC characteristics of the HBT with selfheating effect. . . . . . . . . . . . . . . . . . . . . . . . . . .. 81. 4.8. The IC −VCC characteristics with voltage source VIN . . . . . .. 82. 4.9. An equivalent circuit including thermal networks of an n-finger HBT with constant current bias. . . . . . . . . . . . . . . . .. 83. 4.10 The common-emitter I-V characteristics: (a) ICC versus VCC (b) and IC1 and IC2 versus VCC of the three-finger HBT with thermal effects being considered. . . . . . . . . . . . . . . . .. 86. 4.11 TJ1 with respect to VCC and IIN . . . . . . . . . . . . . . . . . .. 88. 4.12 TJ2 with respect to VCC and IIN . . . . . . . . . . . . . . . . . .. 89. 5.1. The time-variant VOUT , which is simulated by our numerical solver. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 93. 5.2. The time-variant VOUT , which is simulated by HSPICE solver.. 94. 5.3. The plot for the output power spectrum ,which bases on the FFT results of Fig. 5.1. . . . . . . . . . . . . . . . . . . . . .. 5.4. The plot for the output power spectrum ,which bases on the FFT results of Fig. 5.2. . . . . . . . . . . . . . . . . . . . . .. 5.5. 5.6. 95. 96. Zoom-in plots investigation for VOUT simulated by the developed solver. . . . . . . . . . . . . . . . . . . . . . . . . . . .. 98. Zoom-in plots investigation for VOUT simulated by HSPICE. .. 99.

(18) LIST OF FIGURES. 5.7. xviii. Zoom-in plots investigation for IC simulated by the developed solver. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101. 5.8. Zoom-in plots investigation of IC simulated by HSPICE. . . . 102. 5.9. Spectrum of POUT simulated by our solver including the selfheating effect. . . . . . . . . . . . . . . . . . . . . . . . . . . 103. 5.10 Spectrum of POUT simulated by HSPICE solver including the self-heating effect. . . . . . . . . . . . . . . . . . . . . . . . 104 5.11 Spectrum of POUT simulated by our solver without the selfheating effect. . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.1. Output power, simulated by our solver, at the fundamental frequencies (black-filled symbol) and the IM3 products (whitefilled symbol) versus input power. . . . . . . . . . . . . . . . 109. 6.2. Output power, simulated by HSPICE solver, at the fundamental frequencies (black-filled symbol) and the IM3 products (white-filled symbol) versus input power. . . . . . . . . . . . 110. 6.3. A deviation plot of OIP3 value versus frequency spacing ∆f . . 112. 6.4. Plots of OIP3 value versus JC , where symbols are measured data and lines are simulated results of our proposed method. . 113. 6.5. Plots of OIP3 value versus JC , where symbols are measured data and lines are simulated results of harmonic balance method.114. 6.6. A plot of OIP3 values versus IIN , which demonstrates the difference between results with and without the self-heating effect.116.

(19) LIST OF FIGURES. 6.7. xix. Simulated LO-OIP3 values for an HBT under different I IN and VCC biases. . . . . . . . . . . . . . . . . . . . . . . . . . . . 117. 6.8. Simulated HI-OIP3 values for an HBT under different I IN and VCC biases. . . . . . . . . . . . . . . . . . . . . . . . . . . . 118. 6.9. A plot of OIP3 values versus IIN , where VCC is set to be 0.8 V. 119. 6.10 A plot of OIP3 values versus IIN , where VCC is set to be 6.4 V. 120 6.11 Plots of output power, PAE, and power gain versus P in , respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 6.12 Plots of PAE for Finger 1 and 2 versus Pin . . . . . . . . . . . . 122 6.13 Plots of the collector current density (JCC ) versus the input bias current (IIN ) for the cases A, B, C, and D. . . . . . . . . . . . 124 6.14 Plots of the junction temperature (TJ ) versus IIN for the cases B, C, and D. . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 6.15 A comparison of the OIP3 values under different biases I IN between the cases A and B. . . . . . . . . . . . . . . . . . . . 126 6.16 A comparison of OIP3 values of the whole device under different biases IIN between the cases C and D. . . . . . . . . . . 127 6.17 Plots of the input third-order intercept point (IIP3) values versus IIN of the case C and D. . . . . . . . . . . . . . . . . . . . 128 A.1 A EKV MOSFET equivalent circuit for DC and transit analysis.154 B.1 A MOSFET circuit for large signal time domain analysis and the EKV2.6 large signal model for MOSFET simulation. . . . 163 B.2 A plot of VOUT versus time, which is simulated by our method. 165.

(20) LIST OF FIGURES. xx. B.3 A plot of VOUT versus time, which is simulated by HSPICE simulator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165 B.4 The spectrum of output power calculated with time-domain data of our method. . . . . . . . . . . . . . . . . . . . . . . . 167 B.5 The spectrum of output power calculated with HSPICE timedomain data. . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 B.6 A plot of the output power at fundamental frequencies and IM3 products simulated by our method. . . . . . . . . . . . . 169 B.7 A plot of the output power at fundamental frequencies and IM3 products simulated by HSPICE. . . . . . . . . . . . . . . 169 B.8 The circuit of a low noise amplifier (LNA). . . . . . . . . . . 171 B.9 A plot of VOUT for an LNA simulated by the developed method. 172 B.10 A plot of VOUT for an LNA simulated by HSPICE. . . . . . . 172 B.11 The spectrum of VOUT simulated by the developed method. . . 173 B.12 The spectrum of VOUT simulated by HSPICE simulator. . . . . 173 C.1 The maximum norm error of all computed unknowns versus number of the outer iterations. . . . . . . . . . . . . . . . . . 175 C.2 The maximum norm error of the DC simulation versus number of the outer iterations. . . . . . . . . . . . . . . . . . . . . . . 176 C.3 The illustration of the chosen time steps for the convergence comparison. . . . . . . . . . . . . . . . . . . . . . . . . . . . 177 C.4 The maximum norm error versus number of the outer iterations at time steps t1 and t2 . . . . . . . . . . . . . . . . . . . . 178.

(21) LIST OF FIGURES. xxi. C.5 The maximum norm error versus number of the outer iterations at time steps t3 and t4 . . . . . . . . . . . . . . . . . . . . 179 C.6 The maximum norm error of TJ versus number of the outer iterations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180.

(22) Chapter 1 Introduction. N. umerical methods for semiconductor devices, circuits and systems provide an efficient alternative in the design and fabrication of novel. semiconductor integrated circuits (ICs). For example, one can simulate the semiconductor devices with a physical model, like drift-diffusion (DD) model or hydrodynamic (HD) model [1]. The partial differential equations (PDEs) of the devices with spatial information are solved to obtain electrical potential and flux. With the pre-process simulation results, before IC fabrication, the engineer can predict the behavior of the devices and design the desire function in the circuits. Furthermore, the designer can also utilize the simulation tools to correct the design after the characterization of the fabricated test ICs. The procedures of design, simulation, fabrication, and testing form an IC production period. Accurate and efficient computer-aided design tools can shorten these time periods and reduce production cost in the same time.. 1.

(23) 1.1 : Motivation. 2. For integrated circuits in the radio frequency (high frequency) (RFICs), such as mixers, oscillators ,low-noise amplifiers (LNAs), phase-locked loops (PLLs), and power amplifiers (PAs) [2]-[4], reliable simulation tools are necessary to shorten the design period and minimize the prime cost. Since the fabrication of RFICs is highly sensitive to semiconductor process, the designer needs more accurate simulations for circuits in order to meet the narrow design window. For RFICs applications, the physical-based simulation method is not suitable for solving the time variant problem because of its tremendous computation time. With careful physical analysis and proper approximation, the spatial variables can be pre-calculated and turned into parameters in the analytic equations. The physical based model of device structure, therefore, is simplified to model of analytic equations, called a compact model. In mathematics, the partial differential equations (PDEs) are transferred to ordinary differential equations (ODEs) for the integrated spatial terms.. 1.1. Motivation. The junction transistor and field effect transistor are the two major types of active semiconductor devices. The bipolar junction transistor (BJT) and metaloxide semiconductor field-effect transistor (MOSFET) are the examples of these two types of devices. In order to minimize the fabrication cost and improve the high frequency performance, the trend has been towards smaller device size since the invention of IC. Meanwhile, the ”compact” models for.

(24) 1.2 : Background. 3. BJTs and MOSFETs become much more complex as device size shrinks. Furthermore, the ultra high frequency operation of devices leads many additional effects for device modelling, such as distributed capacitances and lossy transmission lines. In practical cases, the circuit model for novel devices may have near or over one hundred parameters, including the noise model and parasitic terms. The complexity of model equations is a stiff problem for the circuit simulation kernel. In other words, the capability of developing a nonlinear circuit solver is necessary to establish a functional device model. In our study, we first attempt to simulate the two-tone intermodulation distortion of the heterojunction bipolar transistors (HBTs). After several failed attempts, we found that the traditional transient analysis in SPICE-like timedomain solvers [5, 6] are not suitable for the distortion analysis of strongly nonlinear circuit. The need for a more accurate and efficient circuit solver propel us to develop our own solution method. The related computation techniques and measurements for the distortion analysis are also incorporated. Finally, we extend our numerical solution method to the MOSFET compact model after the success in simulation of the HBT circuit model.. 1.2. Background. Most of the algorithms for RF nonlinear circuit simulation can be cataloged into time- and frequency-domain methods [7]-[9]. These two classes of methods are usually complementary in many aspects. Therefore, in recent years,.

(25) 1.2 : Background. 4. there have been new hybrid algorithms which combine the advantages of the two classes. Time-domain methods are much more easily perceived in mathematics. The advantages of time-domain methods are the capability to solve a strongly nonlinear problem and less approximation. The direct method, also called transient analysis, solves the ODEs in time domain until the solution arrives at a steady-state. One improved method for RF steady-state simulation is the shooting method [10]-[12]. In this method, the ODEs for the circuits are solved by forcing the constrain that the solution is periodic in the steadystate. This condition is expressed as v(t) = v(t + T ), where v is the vector of node voltage and T is the period. The algorithm determines the solution which leads to steady state as the criterion is met. The main disadvantages of the time-domain method are the stiff requirement for transient analysis and lengthy calculation time. The basis of frequency-domain methods is the harmonic balance method (HB method) [13]-[18]. The HB method is a frequency-domain technique for periodic and quasi-periodic steady-state analysis of a mildly nonlinear circuit. In this method, the key idea is the application of the Kirchhoff’s current law (KCL) [19] for each harmonic. The spectrum of all current at a node is balanced; in other words, the KCL is applied for each independent frequency. The HB method is formulated by expressing the ODEs in terms of a Fourier series [20, 21], then replacing the differentiation in time-domain with algebraic multiplication in frequency-domain. Each variable needs many Fourier.

(26) 1.3 : Historical Development. 5. coefficients to guarantee the precision of this method, hence the size of this system is much large than that of the circuit ODEs. The huge size of system requires the manipulation of relatively dense matrices and consumes a large amount of memory and CPU times. For frequency-dependent dispersive circuit elements, the HB method is capable of handling without difficulty because of the nature of the frequency-domain method. Many improvements for the above two classes of methods have been made, such as mixed frequency-time methods [22, 23] and envelope method [24][27]. These are developed to solve special types of circuit analysis. We note that the mentioned methods for nonlinear RF circuit simulation are largesignal methods. There are also small-signals method, for example, the linear time-varying analysis [28] and Volterra Series analysis [29]. In our developed algorithms, we focus on the development of time-domain direct method. Based on the waveform relaxation (WR) method [30]-[32], the monotone iterative (MI) method [33]-[36], and the Runge-Kutta (RK) method [20, 21], we succeed in solving the circuit ODEs while satisfying the criteria of the RF steady-state analysis in time domain.. 1.3. Historical Development. The circuit simulator first appeared in the late 1960’s. However, its importance arose with the dramatic growth of the IC market in the 1970’s. As the designs became larger and more complicated, the requirement for circuit simulators.

(27) 1.4 : Organization of the Dissertation. 6. increased. The development of the modern circuit simulator is credited in large part of the SPICE group at the University of California, Berkeley. The major advantage of computer-aided circuit simulation tools is formation of the matrix of nodal equations automatically by computer programs with the modified nodal analysis [37]. The modest beginning resulted in the simulation programs being developed and culminated in the release of SPICE in 1972 and SPICE2 in 1975 [5]. Because of the propagation of the standard simulator, the SPICE became very important. The source code for SPICE is available at a nominal cost. In the late 1980’s, Berkeley upgraded SPICE to SPICE3 [6]. SPICE3 was written in C langue and is easier to add new component models. Meanwhile, Berkeley also released a new circuit simulator named Spectre [38, 39]. Spectre utilized the HB method to compute steady-state analysis of nonlinear RF circuits in the frequency domain. Spectre was further developed by Hewlett-Packard, becoming the Agilent Advanced Design System (ADS) [40]. Also, Cadence Design Systems replaced the HB algorithms with transient analysis algorithms and developed SpectreRF [41]. A more comprehensive historical development of the circuit simulator can be found in [39].. 1.4. Organization of the Dissertation. The dissertation is organized as follows: In Chapter 2, we state the studied mathematical models, which include equations of equivalent circuit model,.

(28) 1.4 : Organization of the Dissertation. 7. node equations, and electrical-thermal feedback equations. Then, in Chapter 3, we describe the computational techniques for the simulation of DC, time-domain, and electrical-thermal interaction. The characterization methodology for frequency-domain analysis used in this work is also introduced in this chapter. Besides, the genetic algorithm (GA) for parameter optimization and measurement for intermodulation distortion are briefly explained in Chapter 3. To prove the advantages of our algorithms, we demonstrate the simulated and measured results for HBTs in subsequent chapters: DC simulation and analysis (Chapter 4), time domain simulation and frequency domain analysis (Chapter 5), and intermodulation distortion and power characteristics (Chapter 6). Among these chapters, we discuss the results with and without thermal effects. Our simulation results are compared with the outcomes of other simulators, such as HSPICE and Agilent ADS, and measured data. The DC, power and distortion characteristics of an multi-finger HBT are also shown in Chapter 4 and Chapter 6. After the discussion of results, we give the conclusion and some suggestions for the future work in Chapter 7. In addition, we represent the MOSFET EKV model and related simulation results in Appendix A and Appendix B to show the capability of our kernel in solving different kinds of devices. Finally, to complete the description of our developed numerical solution algorithms, we further remark on the convergence properties of these algorithms in Appendix C..

(29) Chapter 2 Mathematical Models. A. mong modern semiconductor devices, the most widely studied types are bipolar junction transistor (BJT) and field effect transistor (FET).. Though both of these two types of devices are commonly used, the physical principle of them are quite different [42]. Therefore, the formations of equivalent circuit model equations have almost nothing in common for these two types of devices. To demonstrate the adaptability of our proposed simulation kernel, we simulate HBT and MOSFET devices with the Gummel-Poon BJT model and the EKV MOSFET model, respectively. Equivalent circuit models for BJT have been greatly improved since the middle of the twentieth century. The models change with the shift in structure and material of devices [43]-[48]. The Ebers-Moll (EM) model is a well understood large-signal model [49]. This model include main mechanisms for middle level current injection operation of bipolar transistors. By including. 8.

(30) 9. the concept of internal charge control, the Gummel-Poon (GP) model was introduced in 1970. The GP model improves the simulation for low and high level current injection regions. This model also considers the effects such as leakage current at collector-emitter and base-emitter junctions, β-falloff mechanism, and Kirk effect [50]-[52]. Based on EM and GP models, many complex new models have been developed in the past few decades. The researchers focus on improving performance for particular characteristics: dependency on temperature, 2-D or 3-D geometry factors, substrate and leakage current terms, physical structure and material related parameters, and parasitic terms for high frequency operation. For the GP model, there are three internal nodes in the large-signal model. The models developed in recent years have more internal nodes for better simulation of new devices. For examples, the High-Current Model (HICUM) [53]-[55] and Most Exquisite Transistor Model (MEXTRAM) [56, 57] have five internal nodes, and the Vertical Bipolar Inter-Company model (VBIC) [58]-[60] has six internal nodes. The increment of internal nodes can meet the requirement for better fit, but also raises difficulties in nonlinear circuit computation and model parameter extraction. In order to compare with other commercial simulators, such as HSPICE [61] and ADS [40, 62, 63], we use the most popular GP large-signal model in this work. We can still fit the measurement data well by GP model with proper extracted parameters [64]. Thermal effects influence the behavior of semiconductor devices during.

(31) 10. DC or RF operation. Both the ambient temperature and heat generated by device power dissipation can induce thermal effects [65]-[72]. The linearity of the device, which strongly depends on bias condition, can also be affected by the thermal-effects [73]-[77]. We introduce the temperature-dependent terms into the conventional GP model. In the meanwhile, the equations, which express the relation between power dissipation and junction temperature, are formed to construct the thermal-electrical iteration loops. With above additional thermal network, our circuit model can simulate the thermal effects of the multi-finger HBT. Furthermore, the field effect transistor, such as metal-oxide semiconductor field-effect transistor (MOSFET), is another important kind of semiconductor device. Though the MOSFETs are not the earliest of semiconductor transistors, they are the mostly fabricated devices today. There are several famous and popular MOSFET compact models. Two of them are the BSIM and EKV models. The BSIM model was developed by the BSIM Research Group in the Department of EECS at the University of California, Berkeley. In the BSIM model series [78, 79], BSIM3 is the most widely applied MOSFET model in industry today. The latest version of this model is BSIM3v3.2.4, released in December, 2001. However, as the device line width shrinks to 100 nm or narrower, the BSIM3 model becomes less precise. The BSIM4, the extension of the BSIM3 model, was developed to address the MOSFET physical behavior in the sub-100 nm regime. Though accuracy is improved, the complexity of.

(32) 11. the model greatly increased. This situation challenges the corresponding circuit simulator and extraction procedure of model parameters. A similar trend happens to the EKV model. The EKV MOSFET model was developed by the Electronics Laboratories, Swiss Federal Institute of Technology (EPFL), in Lausanne, Switzerland [80]. The newest version is EKV3.0, announced in 2004, and is based on the surface potential model combined with inversion charge linearization. The number of model parameters in EKV3.0 is less than that in BSIM4, but it is still difficult to solve this model. For the simulation demonstrated in Appendix B, we list model equations of the EKV model version 2.6 in Appendix A. The EKV2.6 is the previous version of EKV model. We use it in order to compare our solver with HSPICE solver, which has the EKV2.6 model. In this chapter, the equivalent circuit of GP model and parameters of both electrical and thermal models will be described first. Then, we write the equations of the electrical and thermal models. The internal and external node equations of circuits used in this dissertation are formed in the third section. Finally, we give a brief summary of this chapter..

(33) 2.1 : Equivalent Circuit and Parameters of GP Model. 2.1. 12. Equivalent Circuit and Parameters of GP Model CX RC. C CJCX. CDR. CJCI. I2. IBL2. RB. CDF. ICT /qb. B. CJE. IBL1. I1. BX. E. RE EX. Figure 2.1: An illustration of Gummel-Poon large-signal equivalent circuit model for the bipolar transistor.. Figure 2.1 shows the equivalent circuit of a GP large-signal bipolar junction transistor. The model includes internal nodes: C, B, and E, and external nodes: CX, BX, and EX. The capacitance terms of this model are: 1. Base-emitter junction capacitance, CJE , 2. Intrinsic portion of collector-emitter junction capacitance, C JCI , 3. Extrinsic portion of collector-emitter junction capacitance, C JCX ,.

(34) 2.1 : Equivalent Circuit and Parameters of GP Model. Notation IS BF NF BR NR ISE NE ISC NC IKF IKR RB RE RC CJEO VJE MJE CJCO VJC MJC XCJC TF XTF VTF ITF TR FC M. 13. Description Unit Transport saturation current A Ideal maximum current gain in forward-active mode Ideality factor of the forward current Ideal maximum current gain in reverse-active mode Ideality factor of the reverse current Base-emitter leakage current A Ideality factor of the base-emitter leakage current Base-collector leakage current A Ideality factor of the base-collector leakage current Corner for the forward beta high-current roll-off A Corner for the reverse beta high-current roll-off A Zero bias base resistance Ω Emitter resistance Ω Collector resistance Ω Base-emitter zero bias junction capacitance F Base-emitter junction built-in potential V Base-emitter junction exponential factor Base-collector zero bias junction capacitance F Base-collector junction built-in potential V Base-collector junction exponential factor Factor for intrinsic part of the base-collector capacitance Ideal forward transit time Sec Pre-coefficient for bias dependence of TF Coefficient of VBC dependence of TF V Coefficient of IC dependence of TF A Reverse transit time Sec Coefficient for forward-bias capacitance formula Multiplier factor for the transistor connection -. Table 2.1: A list of BJT GP model parameters.. 4. Diffusion capacitance of the charge due to forward active current, C DF , 5. Diffusion capacitance of the charge due to reverse active current, C DR ..

(35) 2.1 : Equivalent Circuit and Parameters of GP Model. 14. The current terms, I1 , I2 and ICT are the intermediate current variables. With these terms, we can calculate the base, collector and emitter currents. The remaining current terms, IBL1 and IBL2 , are the leakage currents of the base-emitter and base-collector junctions, respectively. Furthermore, q b is the ratio of the base charge QB to the equilibrium base charge QB0 .. Notation Ea Eb XTI XTB BB. Description First bandgap correction factor Second bandgap correction factor Temperature exponent for IS Temperature exponent for BF and BR Fitting parameter for the thermal conductivity. Unit eV/K K -. Table 2.2: A list of parameters for thermal effects modelling.. The GP large-signal model used in this work involves 28 parameters. The purpose of these parameters can be separated into several groups. The IS, BF, NF, BR and NR are the main parameters for the forward and reverse current terms. The leakage current related parameters include ISE, NE, ISC and NC. IKF and IKR are introduced to model the high current region of BJT operation. The internal resistance terms include RB,RC and RE. There are seven parameters: CJEO, VJE, MJE, CJCO, VJC, MJC and XCJC, used for junction capacitance of the model. Finally, the fitting parameter FC and the multiplier factor M are also included in this model. For M connected.

(36) 2.1 : Equivalent Circuit and Parameters of GP Model. 15. identical devices, the parameters with current or capacitance units are multiplied by M. On the other hand, the parameters with resistance units should be divided by M. A detailed description of parameters and their units are shown in Table 2.1 We also include five important parameters to simulate thermal effects. They are Ea , Eb , XTI, XTB and BB. Descriptions for these are listed in Table 2.2..

(37) 2.2 : Model Equations. 2.2. 16. Model Equations. In order to implement the GP model to perform the DC and time-domain characteristic simulations, we must solve the nonlinear large-signal equations. I 1 , I2 and ICT , the intermediate current variables which dominate current behavior, can be written as follows:. ICT. VB − V E IS −1 , · exp I1 = NF · VT BF VB − V C IS − 1 , and · exp I2 = NR · VT BR VB − V C VB − V E . − exp = IS · exp NR · VT NF · VT. VT represents the thermal voltage (VT =. kTJ ). q. (2.1). (2.2). (2.3). For a device operated at 300K,. VT is equal to 0.0259V. In the above equations, VB , VC and VE represent the voltage, which are variables to be solved, at internal node B, C and E, respectively. Leakage currents IBL1 and IBL1 are written as:. VB − V E −1 and IBL1 = ISE · exp NE · VT VB − V C −1 . IBL2 = ISC · exp NC · VT . . (2.4). (2.5). Then, the base charge related term qb is q qb = 1 + 2. r. q1 2 + q2 . 2. (2.6). We note that the q1 and q1 in equation (2.6) are q1 = 1 +. VB − V E VB − V C + VAR VAF. and. (2.7).

(38) 2.2 : Model Equations. 17. VB − V C IS VB − V E IS −1 . · exp −1 + · exp q2 = NR · VT IKR NF · VT IKF (2.8) Equation (2.7) models the Early and reverse Early effect of a BJT device. Unlike the Si BJTs, the base-width modulation phenomenon is rarely observed in InGaP HBTs, due to heavy base doping. Therefore, the forward and reverse Early voltage, VAF and VAR are large enough in magnitude to approximate q1 as 1.0 for the HBTs. The differential voltages among node C, E, B, and BX are defined as the follows. VBC ≡ VB − VC ,. (2.9). VBE ≡ VB − VE , and. (2.10). VBXC ≡ VBX − VC .. (2.11). With the simplified notation, we can write the equations of all diffusion capacitances, CDR =. i h VBC ∂ TR · IS · exp NR·V T. CDF =. ∂VBC ∂ τF ·. Ibf qb. ∂VBE. . ;. and. (2.12). (2.13).

(39) 2.2 : Model Equations. 18. and junction capacitances,   VBE −MJE  CJEO · 1 −   VJE       if VBE ≤ FC · VJE    CJE = CJEO · (1 − FC)−MJE ·      MJE  · V 1 − FC · (1 + MJE) +  BE VJE       if VBE > FC · VJE. CJCX. (2.14). ,.   VBXC −MJC  (1 − XCJC) · CJCO · 1 −   VJC       if VBXC ≤ FC · VJC    = (1 − XCJC) · CJCO · (1 − FC)−MJC ·       · VBXC 1 − FC · (1 + MJC) + MJC  VJC       if VBXC > FC · VJC , and   VBC −MJC  XCJC · CJCO · 1 − VJC         if VBC ≤ FC · VJC    CJCI = XCJC · CJCO · (1 − FC)−MJC ·      MJC  · V 1 − FC · (1 + MJC) +  BC VJC       if VBC > FC · VJC .. (2.15). (2.16). In equation (2.13), two terms can be described more clearly as . Ibf = IS · exp. and. . . VBE NF · VT. 2. τF = TF · 1 + XTF · x · exp. . . −1. . VBC 1.44 · VTF. (2.17). . ,. (2.18).

(40) 2.2 : Model Equations. 19. where x=. Ibf . Ibf + ITF. (2.19). These represent the forward-active current and forward transit time of this device. We also include thermal effects of HBT in this work [77]. To simulate the thermal-electrical feedback mechanism, the temperature dependant equations of some physical parameters are introduced to the GP model: IS(TJ ) = IS · (. ISE(TJ ) = ISE·(. Eg (TJ ) Eg (TN ) TJ XTI )], )−( ) · exp[( k · TJ k · TN TN. (2.20). Eg (TJ ) Eg (TN ) TJ XTI −XTB )], (2.21) )−( ·exp[( ) NE NE · k · TJ NE · k · TN TN. ISC(TJ ) = ISC · (. Eg (TJ ) Eg (TN ) TJ XTI −XTB )], )−( · exp[( ) NC NC · k · TJ NC · k · TN TN (2.22) TJ (2.23) BF (TJ ) = BF · ( )XTB , and TN TJ (2.24) BR(TJ ) = BR · ( )XTB . TN. In the equations (2.20)-(2.22), the energy band gap with temperature dependency is shown as. Ea · T2J Ea · T2N , + Eg (TJ ) = Eg (TN ) + TN + E b TJ + E b. (2.25). where TJ and TN are junction and nominal temperatures, respectively. In our case, TN is set to be the ambient temperature, which is 300K. We need to note that TN is the temperature on the back of the substrate for the highly powered devices. Above equations include the temperature-dependent parameters, not.

(41) 2.2 : Model Equations. 20. involved in the original GP model. The parameters are transport saturation current IS(TJ ), base-emitter leakage current ISE(TJ ), base-collector leakage current ISC(TJ ), ideal forward current gain BF (TJ ), and ideal reverse current gain BR(TJ ). The thermal model also expresses the relation between power dissipation and junction temperature. In consideration of temperature-dependent thermal conductivity, the junction temperature for an n-finger HBT is    TJ1       TJ2  −1 (BB − 1)  [RTH · PD ]} BB−1 =  TJ = TN {1 −  . , TN  ..      TJn. (2.26). where TJn is the junction temperature of n-th finger [42, 52, 71, 72]. R TH ·PD is given by .     RTH · PD =     .  . RT 11 RT 12 . . . RT 1n       RT 21 RT 22 RT 2n   ·   .. . . . ..   . .     RT n1 RT n2 RT nn. . PD1    PD2  .  ..  .   PDn. (2.27). Here RT nn and RT nm denote the self-heating thermal resistance of n-th finger and the coupling thermal resistance which includes the coupled heat from mth finger to n-th finger. Furthermore, the power dissipation of n-th finger is denoted by PDn . We note here that in general the values of RTH must be computed in advance with the help of a three-dimensional thermal analysis method (for example, a finite element simulation of the structures)..

(42) 2.3 : Formulation of Node Equations. 2.3. 21. Formulation of Node Equations. The overall node equations for circuit simulation can be formulated by the Kirchhoff’s current law (KCL) [19]. For each node KCL states:. For all lumped circuits, the algebraic sum of the currents leaving any node is equal to zero at all time t.. For circuits which consist of both linear passive and nonlinear active elements, one should formulate the nodal equations at both the external nodes and the internal nodes of active devices. It is necessary to note that KCL is employed to formulate equations for lumped circuits in classic circuit theory. For a lumped circuit, voltages and currents need to be well defined everywhere in the circuit. Clearly, microwave circuits cannot treated as lumped. At microwave frequencies, distributed-circuit and radiation related to propagation of electromagnetic energy become important. In order to derive KCL from Maxwell’s equations, one needs to assume slow time variations in operation [9, 17, 18]. In this regime, classic circuit theory is valid as a quasi-static approximation of general electromagnetic theory. However, Kirchhoff’s laws are still applicable at high frequencies provided two requirements are met: First, distributed-circuit and radiation effects must be modeled within circuit elements since these effects are ignored by Kirchhoff’s laws. Second, voltages and currents be adequately defined at terminals of every circuit element.

(43) 2.3 : Formulation of Node Equations. 22. including those introduced to model distributed-circuit and radiation effects. For instance, a section of transmission line has to be treated as a two-port. Electromagnetic coupling between two parts of the circuit connected by a transmission line has to be modeled through circuits elements like sections of coupled transmission lines or combinations of lumped elements.. 2.3.1 Nodal Equations of GP Model The GP model has three internal nodes and three external nodes for the HBTs simulated here. The substrate portion of the model is neglected because of the vertical structure of the HBTs. For the circuit shown in Fig. 2.1, we formulate the node equation at node C: CJCX ( dVdtBX −. I2 + IBL2 −. dVC ) dt. ICT qb. +. B − + CDR ( dV dt. VCX −VC RC. dVC ) dt. B − + CJCI ( dV dt. dVC )+ dt. (2.28). =0,. node E: CDF (. ICT VEX − VE dVB dVE dVB dVE = 0 , and + ) + I1 + IBL1 + − ) + CJE ( − RE qb dt dt dt dt (2.29). node B: B − CDR ( dV dt. dVC ) dt. B − + CJCI ( dV dt. I1 + IBL1 + I2 + IBL2 +. VB −VBX RB. dVC ) dt. B − + CDF ( dV dt. dVE ) dt. B − + CJE ( dV dt. dVE )+ dt. =0. (2.30).

(44) 2.3 : Formulation of Node Equations. 23. VCC. RCCS CX. B2. VIN. BX RB2. Vin. HBT. VOUT. EX. REE. Figure 2.2: A testing circuit for an HBT without thermal network.. 2.3.2 Nodal Equations of Simulated Circuits The external nodes of the active device are connected by lumped elements or other active devices. At these connecting nodes, the external nodal equations are formulated with KCL. For the circuit as shown in Fig. 2.2, four node equations are given by the following, node BX: CJCX (. VB − VBX VB2 − VBX dVC dVBX =0, + )+ − RB2 RB dt dt. (2.31). VB2 = VIN + Vin ,. (2.32). VC − VCX VCC − VCX = 0 , and + RCCS RC. (2.33). node B2:. node CX:.

(45) 2.3 : Formulation of Node Equations. 24. VCC. IC CX. IB. Pin. BX. IIN. HBT. PD. TJ. EX. Figure 2.3: A simulation circuit for an HBT with thermal network.. node EX: VE − VEX VEX =0. − REE RE. (2.34). In this circuit, VCC and VIN are the DC bias voltages. Time-variant voltage input excitation is denoted by Vin . This circuit also includes resistors which serious connect node BX to B2, node EX to ground, and node CX to voltage source VCC . RB2 , REE , and RCCS represent these series resistors, respectively. To simulate the single-finger HBT with self-heating effect, the electrical circuit and thermal network is plotted in Fig. 2.3. P D and TJ represent the power dissipation and junction temperature, respectively. V CC and IIN , are collector bias voltage and base injection current. The input excitation is a sinusoidal power signal Pin . External nodal equations at node BX, CX, and.

(46) 2.3 : Formulation of Node Equations. 25. VCC. ICC. IB1. M1. CX. IC1 IB2. M2. IC2. IBn. Mn. ICn. BX. EX1. Pin. PD1 TJ1. EX2. PD2 TJ2. EXn. PDn TJn. IIN Figure 2.4: An equivalent circuit of the investigated multi-finger HBT with high frequency input excitation.. EX are written as CJCX (. Pin VB − VBX dVC dVBX =0, + IIN + )+ − VBX RB dt dt. (2.35). VCX = VCC , and. (2.36). VEX = 0 .. (2.37). For an n-finger HBT, we draw the equivalent circuit with thermal networks in Fig. 2.4. The whole device is represented by a circuit of n shunted transistors. The k-th finger is denoted by M k , k = 1, 2, ..., n. The voltage, current, power dissipation, and junction temperature of each distinct finger can be written as VXk , IXk , PDk and TJk . Each finger is excited by the same time-variant.

(47) 2.3 : Formulation of Node Equations. 26. power source Pin and biased with the same DC voltage source VCC at collector and DC current source IIN at base. Now we list the node equations of external nodes: at node BX: n X Pin VBk − VBX dVCk dVBX =0, + IIN + )+ − CJCXk ( V RB dt dt BX k=1 k=1. n X. (2.38). at node CX: VCX = VCC ,. (2.39). and at Node EX of k-th finger: VEXk = 0 ,. k = 1, 2, ...., n .. (2.40). Finally, the input excitations for above circuits need to be defined more rigorously. For the circuit in Fig. 2.2, the two-tone input excitation is a timevariant voltage source Vin = Vm [cos(ω1 t) + cos(ω2 t)],. (2.41). where ω1 = 2πf1 and ω2 = 2πf2 . f1 and f2 are fundamental frequencies of the tow-tone signal. The amplitude of this signal is denoted by V m . On the other hand, for the circuit in Fig. 2.3 and Fig. 2.4, the input power signal of two-tone modulation is written as Pin = Pm [cos(ω1 t) + cos(ω2 t)],. (2.42).

(48) 2.4 : Summary of This Chapter. 27. where Pm represents the amplitude of the time-variant input signal. Here we define the input signal as power sources because of the convenience in identifying the input power level in the frequency-domain.. 2.4. Summary of This Chapter. In this chapter, we discussed the following issues in detail: • the equations of GP large-signal model, • a brief description for the parameters of GP model, • the testing circuits utilized in this study, • the internal and external node equations for the testing circuits, and • the equations for the thermal-electrical interaction. All the equations and parameters mentioned in this chapter will be used in the numerical methods to be presented in Chapter 3. In addition, we will also demonstrate the simulation and measurement results of the tested circuits in this dissertation..

(49) Chapter 3 Computational Techniques and Characterization Methodology. C. ircuit simulation solves a system of coupled ordinary differential equations (ODEs) in time- or frequency-domain. These ODEs are the. node or loop equations formulated by electrical circuit analysis. The nodal equations of linear lumped elements can be formulated in matrix form automatically by computer programs [9]. The major difficulty arises from the nonlinear portion of the circuit. Here, our efforts focus on solving the nonlinear elements in the time domain without neglecting details. We introduce the waveform relaxation (WR) method [30]-[32] and monotone iterative (MI) method [33]-[36] to solve the ODEs in the time domain directly. With the proposed algorithms, we overcome error propagation, the major problem of two-tone analysis in the time domain. The solution maintains accuracy after. 28.

(50) 29. many simulated time steps. This chapter is organized as follows. In section 3.1, we state the numerical algorithms for both stationary and time-variant circuit simulation. In the second section, we introduce the genetic algorithm (GA) for parameter extraction and optimization. For further analysis, additional calculation techniques for the simulated time-domain data are presented in Sec. 3.3. These methods provide useful information in frequency-domain, such as the intermodulation distortion characteristic. In the fourth section, the measurement for intermodulation characterization is described briefly. Finally, a brief summary of this chapter is given..

(51) 3.1 : Numerical Algorithms. 3.1. 30. Numerical Algorithms. After the formulation of node equations are , the circuit simulation still need to solve them correctly. In this section, we will present a time-domain circuit simulation method. In order to solve the system of ODEs in the large-scale time period efficiently, we propose a decoupled and globally convergent simulation technique for solving the system ODEs. First, under the steady state condition, we find the DC solutions, which provide a starting point to compute the time dependent solutions.. 3.1.1 Algorithm for Solving the System of Nonlinear Algebraic Equations DC equilibrium is an important parts of circuit analysis. In a traditional RLC circuit, the circuit equations form a system of algebraic-differential equations represetned by: ∂q(x) + F (x) − i(t) = 0 . ∂t. (3.1). Under certain conditions, the circuit will reach an equilibrium state, where all circuit variables such as currents, voltages, and charges are time-invariant. In this situation, the time varying function. ∂q(x) ∂t. will become zero, and i(t) is. replaced by the constant i0 . Hence the equation (3.1) becomes F (x0 ) − i0 = 0 .. (3.2).

(52) 3.1 : Numerical Algorithms. 31. The unknown to be solved is x0 . Due to the strongly nonlinear function F , the DC analysis problem is hard to solve and can even challenge the most robust algorithms, such as Newton’s method. Under the DC condition, the capacitor can be seen as an open circuit. By applying KCL to the circuit shown in Fig. 2.2, we can formulate a system of nonlinear algebraic equations for DC simulation. In other words, the equations (2.28)-(2.30) and (2.31)-(2.34) become I2 + IBL2 −. ICT VCX − VC =0, + RC qb. (3.3). I1 + IBL1 +. ICT VEX − VE =0, + RE qb. (3.4). I1 + IBL1 + I2 + IBL2 +. VB − VBX =0, RB. (3.5). VB2 − VBX =0, RB2. (3.6). VB2 = VIN ,. (3.7). VC − VCX VCC − VCX = 0 , and + RCCS RC VE − VEX VEX =0. − REE RE. (3.8). (3.9). The unknowns are VC , VE , VB , VBX , VCX , and VEX . If we let f be a general form of these nonlinear equations, f : D → Rn .. (3.10).

(53) 3.1 : Numerical Algorithms. 32. The node equations (3.3)-(3.9) can be rewritten as ICT VCX − VC = 0, + RC qb ICT VEX − VE = 0, + fE (VC , VE , VB , VBX , VEX , VCX ) = I1 + IBL1 + RE qb VB − VBX = 0, fB (VC , VE , VB , VBX , VEX , VCX ) = I1 + IBL1 + I2 + IBL2 + RB VIN − VBX = 0, fBX (VC , VE , VB , VBX , VEX , VCX ) = RB2 VC − VCX VCC − VCX = 0, and + fCX (VC , VE , VB , VBX , VEX , VCX ) = RCCS RC VE − VEX VEX = 0. − fEX (VC , VE , VB , VBX , VEX , VCX ) = REE RE fC (VC , VE , VB , VBX , VEX , VCX ) = I2 + IBL2 −. (3.11) Involving a notational change, .         Z=        . . VC    VE     VB  ,  VBX     VCX    VEX. (3.12). and decoupled Scheme shown in Fig. 3.1, the iteration can be written as.

(54) 3.1 : Numerical Algorithms. 33. (k+1). ) = 0,. (k+1). ) = 0,. (k+1). ) = 0,. fC (Z (k) → VC (k+1). → VE. (k+1). → VB. (k+1). fE (Z (k) , VC (k+1). , VE. , VE. (k+1). , VB. → VBX ) = 0,. (k+1). (k+1). → VCX ) = 0, and. fB (Z (k) , VC (k+1). fBX (Z (k) , VC (k+1). fCX (Z (k) , VC. (k+1). fEX (Z (k) , VC. (k+1). , VE. (k+1). , VE. , VB. , VBX. (k+1). , VBX , VCX. , VB. (k+1). (k+1). (k+1). (k+1). (k+1). → VEX ) = 0. (3.13). (0). (0). (0). (0). (0). (0). Giving initial guesses VC , VE , VB , VBX , VCX , and VEX , we first solve the nonlinear equation fC and substitute the new VC into fE to approximate a new VE . Then the algorithm proceeds by repeating the above procedure for VB , VBX , VCX , and VEX , respectively. Each decoupled algebraic equation is solved with the MI iterative method. When the overall outer loop error is within the error tolerance, the solution algorithm stops. The steps of the solution algorithm shown in Fig. 3.1 are described as follows: • Step 1: With the given bias conditions VIN and VCC , we properly input a set of values V(0) X to all unknowns as initial guesses. • Step 2: The system of nonlinear algebraic equations are decoupled through WR method. • Step 3: With VE , VB , and other variables kept constant, the decoupled algebraic equation of unknown VC are solved with MI.

(55) 3.1 : Numerical Algorithms. . formula to get new V(n+1) C − V(n) • Step 4: If the norm error |V(n+1) C | is less than the preset tolC is used to substitute the VC in the next erance (TOL), V(n+1) C algebraic equations of unknown VE , VB and etc. If not, Step 3 is repeated until the norm error meets the prescribed criterion. • Step 5: The new VC solved in previous iterative loop is substituted. Then, the decoupled equation of unknown VE are solved with other variables kept constant. • Step 6: If the norm error |V(n+1) − V(n) E E | is less than the preset tolis used to substitute the VE in the next erance (TOL), V(n+1) E algebraic equations of unknown VB , VBX and etc. If not, Step 5 is repeated until the norm error meets the prescribed criterion. • Step 7: Similar iterative loops as that described in Step 5 and Step 6 are performed for the remaining unknowns. • Step 8: If all norm errors (|VC − VCO |, |VE − VEO |, ...) are less than TOL, we arrive at a set of solution for the DC simulation. If not, the data of VXO is updated with newly solved VX , and the procedure is repeated from Step 3.. 34.

(56) 3.1 : Numerical Algorithms. 35. Given bias conditions VIN and VCC. Giving the initial condition for VX , where X = C, E, B, BX, CX, EX, etc.. VCO <− VC(0) Solve the corresponding equation with the MI formula VC(n+1) = VC(n) + F1(VC(n), VEO , VBO , ....) / λ. |VC(n+1) − VC(n) | < TOL. No. Update the result VC(n) <− VC(n+1). Yes. VEO <− VE(0) Solve the corresponding equation with the MI formula VE(n+1) = VE(n) + F2(VC , VE(n), VBO , .... ) / λ. Update all data VEO = VE , VBO = VB , VCO = VC , ..... |VE(n+1) − VE(n) | < TOL. No. Update the result VE(n) < − VE(n+1). Yes . . .. Yes. No. |VC − VCO | < TOL, |VE − VEO | < TOL, etc. Valid ?. Yes Obtain all computed data for stationary simulation. Figure 3.1: A flowchart of the decoupled methodology for the stationary circuit simulation..

(57) 3.1 : Numerical Algorithms. 36. 3.1.2 Time Domain Solution Algorithm In order to simulate the circuit shown in Fig. 2.2, we formulate a system of ODEs and algebraic equations by using the KCL for this circuit (see Chapter 2, equations (2.28)-(2.30) and (2.31)-(2.34)). For a specified time period T , these nonlinear ODEs needs to be solved in the time domain. The following steps are taken to accomplish this [35, 82]: • Step 1: The initial time step t, total time period T , and time step 4t. are determined first. • Step 2: The decoupling method is used to decouple all equations (2.28)-(2.30) and (2.31)-(2.34). • Step 3: Each decoupled ODE is solved sequentially with the MI and Runge-Kutta methods. • Step 4: Each MI loop is convergence tested. • Step 5: Overall outer loop is convergence tested. • Step 6: If the specified stopping criterion is not met for the outer loop, all data is updated with newer results and recalculate from Step 3. • Step 7: Step 3 - 6 are repeated until the time step meets the specified time period T . First, all coupled ODEs are decoupled by the waveform relaxation (WR) method [30]-[32]. The WR method starts with a general nonlinear system of.

(58) 3.1 : Numerical Algorithms. 37. D ODEs with associated initial conditions, dV = f (V, t) , dt. V (0) = V0 ,. t ∈ [0, T ] ,. (3.14). where T > 0, f : RD × [0, T ] −→ RD , V0 = [V1,0 , V2,0 , · · · , VD,0 ] ∈ RD is the initial vector of V , and V (0) = [V1 (t), V2 (t), · · · , VD (t)] ∈ RD is the solution vector. The system can be written as follows:    d  V = f1 (V1 , V2 , · · · , VD , t), V1 (0) = V1,0  dt 1         dtd V2 = f2 (V1 , V2 , · · · , VD , t), V2 (0) = V2,0  ..    .         d VD = fD (V1 , V2 , · · · , VD , t), dt. VD (0) = VD,0. (3.15). .. The WR method for solving Eq. (3.14) is a continuous-time iterative method. Therefore, given a function which approximates the solution, it calculates a new approximation along the whole time-interval of interest. Clearly, it differs from most standard iterative techniques in that its iterations are functions in time instead of a scalar value. The iteration formula is chosen in such a way that one avoids having to solve a large system ODEs. A particularly simple, but often very effective iteration scheme which maps the old iterate V (n−1) is written here:.

(59) 3.1 : Numerical Algorithms. 38.   (n) (n) (n−1) (n−1) (n−1)  d  , t), V = f1 (V1 , V2 , V3 , · · · , VD  dt 1         dtd V2(n) = f2 (V1(n) , V2(n) , V3(n−1) , · · · , VD(n−1) , t),  ..    .        d V (n) = fD (V1(n) , V2(n) , V3(n) , · · · , V (n) , t), D dt D. (n). V1 (0) = V1,0 (n). V2 (0) = V2,0. (n). VD (0) = VD,0. . (3.16). This is similar to the Gauss-Seidel method for iteratively solving linear and nonlinear systems of algebraic equations; the so-called Gauss-Seidel waveform relaxation scheme. This scheme converts the task of solving a differential equation in D variables into a task of solving a sequence of differential equations in a single variable. A closely related iteration is the Jacobi waveform relation scheme, the iteration formula is given by:   (n) (n) (n) (n−1) (n−1) (n−1)  d  , t), V1 (0) = V1,0 V = f1 (V1 , V2 , V3 , · · · , VD  dt 1       (n)   dtd V2(n) = f2 (V1(n−1) , V2(n) , V3(n−1) , · · · , VD(n−1) , t), V2 (0) = V2,0  ..    .         d V (n) = fD (V1(n−1) , V2(n−1) , V3(n−1) , · · · , V (n) , t), D dt D. (n). VD (0) = VD,0. .. (3.17). Note that the Jacobi algorithm is fully parallel. The equations can be solved simultaneously. The Gauss-Seidel waveform relaxation algorithm is formulated in the following. Algorithms:.

(60) 3.1 : Numerical Algorithms. 39. n = 0; choose Vx0 (t) for t ∈ [0, T ], x = 1, 2, · · · , D do n=n+1 for x = 1, 2, · · · , D solve. (n) d V dt x. (n). (n). (n). (n). (n−1). (n−1). = fx (V1 , V2 , · · · , Vx−1 , Vx , Vx+1 , · · · , VD. , t). with Vxn (0) = Vx,0 end for while |Vxn − Vxn−1 | < TOL, for x = 1, 2, · · · , D. Each decoupled ODE is solved with the MI algorithm. To clarify the MI algorithm for the numerical solution of the decoupled nonlinear ODEs, we write the ODEs (2.28)-(2.30) and (2.31) in the following form dVx = f (Vx , t) , dt. (3.18). where Vx is the unknown to be solved and f is the collection of nonlinear functions from each node. We define the MI parameter λ =. ∂f ∂Vx. and insert λ. into Eq. (3.18), arriving at the MI equation dVx = f (η, t) − λ(Vx − η) , dt. (3.19). where v0 ≤ η ≤ ω0 is a value in [0, T ], and v0 and ω0 are the lower and upper solutions of Eq. (3.18), respectively. Based on nonlinear behaviors of each decoupled circuit ODEs,we will show mathematically that the solution algorithm exhibits monotone convergence..

(61) 3.1 : Numerical Algorithms. 40. To clarify the MI algorithm [35, 36, 81] for the numerical solution of the decoupled nonlinear ODEs, we rewrite the decoupled ODEs as the following form: (g). dVX dt. (g). = f (VX , t). (g) VX (0) (g). =. (g) V X0. ,. (3.20). is the unknown to be solved, g is the decoupling index (g =. where VX. 0, 1, 2, · · · ). We note that the f is a collection of nonlinear functions with f ∈ C[I × R, R] and I = [0, T ]. For a fixed index g and X, since the upper (g). (g). (g). (g). and lower solutions, V X and V X , exist in the circuit and V X ≥ V X , we (g). (g). (g). can prove the solution existence in the set Ω = {(t, VX ) | V X ≥ VX. ≥. (g). V X , ∀t ∈ I} for each decoupled circuit ODE. (g). (g). Theorem 3.1.1 Let V X and V X be the upper and lower solutions of Eq. (3.20) (g). (g). in C 1 [I × R, R] such that V X ≥ V X in the time interval I and f ∈ C[I × (g). (g). (g). R, R]. Then there exists a solution VX of Eq. (3.20) such that V X ≥ V X in. the time interval I. Proof 3.1.1 It is a direct result with the continuous property of f , here the comparison theorem is applied [35]. Remark 3.1.1 We note that for each decoupled ODE, the nonlinear function (g). f is nonincreasing function of the unknown VX and the upper and lower (g). (g). solutions V X (0) and V X (0) of Eq. (3.20) in I can be found. We can further (g). (g). prove there exists a unique solution VX of Eq. (3.20) in I and V X (0) ≥ (g). (g). VX ≥ V X (0)..

(62) 3.1 : Numerical Algorithms. 41. We see that the Theorem 3.1.1 provides an existence condition of the problem, and we can now describe a monotone constructive method for the computer simulation of the circuit ODEs. The constructed sequences will converge to the solution of Eq. (3.20) for all decoupled ODEs in the circuit simulation. In this condition, instead of original nonlinear ODE to be solved, a transformed ODE is solved with such as the RK method. Now we state the main result for the solution of each decoupled circuit ODEs. (g). (g). Theorem 3.1.2 Let the f ∈ C[I × R, R], V X (0) and V X (0) be the upper (g). (g). and lower solutions of Eq. (3.20) in I. Because f (t, VX ) − f (t, VeX ) ≥. (g) (g) (g) (g) (g) (g) −λ(VX − VeX ), V X (0) ≥ VX ≥ VeX ≥ V X (0) and λ ≥ 0, sequences. (g) unif.. (g). (g) unif.. (g). V Xn −→ V X and V Xn −→ V X as n −→ ∞ monotonically in I [35]. (g). (g). Proof 3.1.2 For V ∈ C[I × R, R] such that V X ≥ V ≥ V X , we consider. the following transformed ODE equation for the fixed g and X dVX dt. (g). = f (V, t) − λ(VX − V) (g) VX (0). =. (3.21). ,. (g) V X0 (g). (g). then ∀V, ∃!VX of Eq. (3.20) in I. Define ΘV = VX and we can prove: (i) (g). (g). (g). (g). ΘV X (0) ≥ V X (0) and ΘV X (0) ≤ V X (0); (ii) Θ is a monotone operator (g). (g). (g). (g). (g). (g). in [V X (0), V X (0)] ≡ [VX ∈ C[I, R] | V X (0) ≥ VX ≥ V X (0)]. (g). Now we construct two sequences by using the mapping Θ : ΘV Xn = (g). (g). (g). V Xn+1 and ΘV Xn = V Xn+1 and by above observations, the following rela-. tion holds (g). (g). (g). (g). V X 0 ≥ · · · ≥ V Xn ≥ V X n ≥ · · · ≥ V X 0.

(63) 3.1 : Numerical Algorithms. 42. (g). (g) unif.. (g) unif.. (g). in I. Hence V Xn −→ V X and V Xn −→ V X as n −→ ∞ monotonically in (g). (g). I. Furthermore the V Xn and V Xn satisfy  (g)  dV  (g) (g)   Xdtn+1 = f (V (g) Xn , t) − λ(V Xn+1 − V Xn ).   (g)  V (g) Xn (0) = VX0. and.  (g)   dV (g) (g) (g)   Xdtn+1 = f (V Xn , t) − λ(V Xn+1 − V Xn )  (g)   V Xn (0) = VX(g) 0 (g). ,. (3.22). ,. (3.23). (g). respectively. Thus V X and V X are the solution of Eq. (3.20).. Theorem 3.1.3 For decoupled ODEs. the nonlinear function f is nonincreas(g). (g). (g). (g). (g). (g). (g). ing in VX and f (t, VX1 ) − f (t, VX2 ) ≥ −λ(VX1 − VX2 ), ∀VX1 ≥ VX2 . o∞ n n (g) o∞ (g) converge uniformly and monotonically to and V Xn Thus V Xn n=1. n=1. the unique solution. (g) VX. of Eq. (3.20) [35].. Proof 3.1.3 By using Theorem 2 and note the nonincreasing property of f , the result follows directly..

(64) 3.1 : Numerical Algorithms. 43. Fig. 3.3 shows the flowchart of the proposed time domain simulation methodology. In order to solve the equations (2.28)-(2.30) and (2.31)-(2.34) in the time domain, we solve VC , VE , VB , VBX , VCX , and VEX sequentially by the decoupled scheme of the system ODEs. Here we state the computational steps of the methodology: • Step 1: For a given time interval T , we obtain the DC solution as an initial guess and start the time evolution processes. • Step 2: Giving a time step (t), we can apply the MI and RungeKutta (RK) method [20, 21] to solve the decoupled ODE in (t−1). Eq. (2.28) with the solved solution (VC (t−1). (t−1). (t−1). VBX , VCX ). The notation VX. (t−1). , VE. (t−1). , VB. ,. represents the solution. of voltage VX at time step (t − 1). In the monotone loop, (t−1). VC. has been solved with the RK formula (fourth-order. RK method is used in practical case). The estimation time step of RK method is adaptive to fit the error tolerance criterion in this work [20]. (t). (t). • Step 3: If the norm error |VC − VCO | is less than the given error tolerance value (TOL), the solution algorithm will exit the monotone loop and enter the next decoupled ODE Eq. (2.29). (t). (t). Otherwise we update VCO with VC and continue the monotone loop described in Step 2 until the tolerance condition has been achieved..

(65) 3.1 : Numerical Algorithms. 44. (t). • Step 4: With the newly solved VC , we can solve the ODE Eq. (2.29) with the same method as described in Step 2. (t). (t). • Step 5: When the norm error |VE −VEO | is less than TOL, we solve (t). (t). the next ODE Eq. (2.30) with the new results (VC , VE ). Otherwise we return to the Step 4. • Step 6: Similar iterative loops, such as that formed by Step 2 and (t). Step 3, are executed for all unknown variables. As (VC , (t). (t). (t). (t). VE , VB , VBX , VCX ) converge in their respective inner loops, we apply the convergence test for the overall outer loop. If the outer loop error is less than the TOL, this calculated set (t). (t). (t). (t). (t). of (VC , VE , VB , VBX , VCX ) is the solution for this given time step (t). The process is then repeated for the next time (t). step (t + 1). If the TOL is not satisfied, we update all (VCO , (t). (t). (t). (t). VEO , VBO , VBXO , VCXO ) and continue the inner loop for each ODE until the outer loop tolerance condition is met. • Step 7: When t exceeds the given time interval T , the simulation will terminate and the results be ready for further processes, such as fast Fourier transform (FFT). The convergence properties of this time-domain solution algorithm are presented in Appendix C..