www.elsevier.com/locate/cpc
A novel numerical approach to heterojunction
bipolar transistors circuit simulation
Yiming Li
a,b,∗, Kuen-Yu Huang
c aNational Nano Device Laboratories, Hsinchu 300, TaiwanbMicroelectronics and Information Systems Research Center, National Chiao Tung University, Hsinchu 300, Taiwan cInstitute of Electronics, National Chiao Tung University, Hsinchu 300, Taiwan
Received 9 August 2002; accepted 25 November 2002
Abstract
In this paper we present a novel computational method for calculating the heterojunction bipolar transistor (HBT) physical characteristics in the time domain. To calculate the HBT high frequency properties, the Gummel–Poon equivalent circuit model is applied to replace the HBT in the circuit and a set of governing ordinary differential equations (ODEs) is formulated. We directly decouple the system ODEs and solve each decoupled ODE with the monotone iterative method in the time domain. This solution methodology proposed here has been applied to semiconductor device simulation by us earlier, and we find this method for the HBT simulation has good accuracy and converges globally. Compared with the HSPICE circuit simulator results, our results present the accuracy, efficiency, and robustness of the method.
2002 Elsevier Science B.V. All rights reserved.
PACS: 73.40.Kp; 72.30.+q; 07.05.Tp; 02.60.Cb; 02.30.Hq
Keywords: HBT; ODE; Computer simulation; Monotone iterative method
1. Introduction
Due to the extremely high linearity of HBTs operating at high frequency, these semiconductor devices and structures have been of great interest for physical applications, such as wireless communications in recent years [1–4]. The harmonic balanced method [5–8] is a way for the solution of the HBT ODEs in the frequency domain. This frequency domain approach has some merits, but it has also limitations in studying the physical properties of HBTs with time evolution. Contrary to those solution methods in the frequency domain, another computational alternative to analyze the physical behavior for an HBT device is to solve a set of equivalent circuit ODEs in the time domain directly. Together with the fast Fourier transformation (FFT), the time domain approach can be extended to study the HBT high frequency properties in the time and frequency domain at the same time. The discretized ODEs in circuit simulation are solved traditionally with the Newton’s iterative (NI) method [9]. It is
* Corresponding author.
E-mail address: ymli@mail.nctu.edu.tw (Y. Li).
0010-4655/02/$ – see front matter 2002 Elsevier Science B.V. All rights reserved. doi:10.1016/S0010-4655(02)00849-4
method to solve the ODEs. The MI method has been successfully developed and applied to semiconductor device simulation by us earlier [12–17]. Based on the decoupling procedure (so-called the waveform relaxation (WR) method [20]) and MI method, the HBT circuit ODEs are directly solved in the time domain. The computed time domain results are then analyzed with the FFT to obtain the necessary information in frequency domain. This approach is strongly dependent on the robustness of the nonlinear ODE solver; compared with the conventional NI method, this method converges globally and is inherently parallel. First of all, a set of ODEs are decoupled with the decoupling algorithm. Each decoupled nonlinear ODE is then solved directly with the MI technique and the Runge–Kutta (RK) method. The proposed computational approach has been successfully implemented on a PC-based cluster with message passing interface (MPI) library [21]. The primary parallel results show that a well-designed parallel algorithm can reduce the execution time up to an order of magnitude.
This paper is organized as follows. Section 2 states the HBT ODE model, and introduces computational procedure. Section 3 presents the simulation results. Section 4 draws the conclusions.
2. An HBT circuit model and simulation method
As shown in Fig. 1(a), based on the node current flow conservation (the well-known Kirchhoff’s current law) and utilize the Gummel–Poon large signal equivalent circuit model (Fig. 1(b)) for the HBT device [9,22,23], the complete simulation model can be formulated with nodal equations. The system of node equations for time dependent HBT circuit is a set of nonlinear coupled ODEs. At nodes C, E, and B we have the following differential equations (1)–(3), respectively. CJCX dVBX dt − dVC dt + CDR dVB dt − dVC dt + CJCI dVB dt − dVC dt + I2+ IBL2− ICT q1 2 + (q1 2)2+ q2 +VCX− VC RC = 0, (1) CDF dVB dt − dVE dt + CJE dVB dt − dVE dt + I1+ IBL1 + ICT q1 2 + (q1 2)2+ q2 +VEX− VE RC = 0, (2) CDR dVB dt − dVC dt + CJCI dVB dt − dVC dt + CDF dVB dt − dVE dt + CJE dVB dt − dVE dt + I1+ IBL1+ I2+ IBL2+ VB− VBX RB = 0. (3)
Similarly, at nodes BX, CX, and EX we formulate, respectively, the equations as follows: CJCX dVC dt − dVBX dt +VB− VBX RB + (VIN+ Vin)− VBX RB2 = 0, (4)
(a) (b)
Fig. 1. (a) The HBT circuit for the high frequency simulation, (b) a plot of the Gummel–Poon equivalent circuit model for the HBT device.
VC− VCX RC + VCC− VCX RCCS = 0, (5) VE− VEX RE − VEX REE = 0. (6) Eqs. (1)–(4) are the ODEs, and Eqs. (5) and (6) are the algebraic equations. These equations are subject to proper initial values at time t= 0 for all unknowns to be solved. All currents I and capacitances C above are nonlinear functions of unknown variables. These nonlinear terms are
I1= IS BF exp VB− VE NFVT − 1 , (7) I2= IS BR exp VB− VC NRVT − 1 , (8) IBL1= ISE exp VB− VE NEVT − 1 , (9) IBL2= ISC exp VB− VC NCVT − 1 , (10) ICT= −IS exp VB− VC NRVT − exp VB− VE NFVT , (11) CDR= ∂ ∂(VB− VC) TRISexp VB− VC NRVT , (12)
CJE= − VJE if (VB− VE) FCVJE, CJEO(1− FC)−MJE 1− FC(1+ MJE)+ MJE VJE (VB− VE) , if (VB− VE) > FCVJE, (14) CJCX= (1− XCJC)CJCO 1−VBX− VC VJC −MJC , if (VBX− VC) FCVJC, (1− XCJC)CJCO(1− FC)−MJC 1− FC(1+ MJC)+ MJC VJC (VBX− VC) , if (VBX− VC) > FCVJC, (15) CJCI= XCJCCJCO 1−VB− VC VC −MJC , if (VB− VC) FCVJC, XCJCCJCO(1− FC)−MJC 1− FC(1+ MJC)+ MJC VJC (VB− VC) , if (VB− VC) > FCVJC, (16)
where VT =kTq is thermal voltage; q1, q2and τF are as follows:
q1= 1 + VB− VE VAR + VB− VC VAF , (17) q2= IS IKF · exp VB− VE NFVT − 1 + IS IKR exp VB− VC NRVT − 1 , (18) τF = TF 1+ XTF Ibf Ibf+ ITF 2 exp VB− VC 1.44VTF . (19)
The model parameters above for the In0.49Ga0.51P/GaAs HBT in our simulation are listed in Table 1. There
are 4 coupled ODEs with the nonlinear current and capacitance models have to be solved and the unknowns to be calculated in the system of ODEs are VC, VE, VB, VBX, VCX, and VEX, respectively. We note that the system consists of strongly coupled nonlinear ODEs, due to the exponential dependence of current and capacitance models [1,20,22].
We propose here a decoupled and globally convergent simulation technique to solve the system ODEs in the large-scale time domain directly. Firstly, under the steady state condition, we find the DC solution as the starting point to compute other time dependent solutions. For a specified time period T , to solve these nonlinear ODEs in the time domain, our approach consists of following steps:
(i) Let an initial time step t be given.
Table 1
A list of the parameters used in the modeling and simulation Parameters Numeric value Unit
BF 86.95 – BR 1.47 – CJCO 24.27E–15 F CJEO 130.0E–15 F FC 0.50 – IKF 0.1815 A IKR 1.032E–3 A IS 2.85E–24 A ISC 2.142E–14 A ISE 2.34E–18 A ITF 419.80E–3 A MJC 0.266 – MJE 0.1188 – NC 1.954 – NE 1.910 – NF 1.068 – NR 1.060 – RB 48.130 Ohm RB2 8.750 Ohm RC 6.750 Ohm RE 1.256 Ohm TF 2.680E–12 second TR 350.0E–12 second VJC 0.7161 V VJE 1.3670 V VTF 66.0000 – XCJC 0.3428 – XTF 275.6000 –
(iii) Each decoupled ODE is solved sequentially with the MI and RK methods. (iv) Convergence test for each MI loop.
(v) Convergence test for overall outer loop.
(vi) If the specified stopping criterion is reached for the outer loop, then go to step (vii), else update the newer results and back to step (iii).
(vii) If t < T , t= t + t and repeat the steps (iii)–(vi) until the time step meets the specified time period T . For a given specified time step t and the previous calculated results (i.e. results at time step t− 1), Fig. 2 is the flowchart of the proposed simulation procedure for the solutions (VC, VE, VB, VBX, VCX, VEX) in the large-scale time domain simulation. In step (ii), the decoupling algorithm solves the circuit equations sequentially, for instance, the VC in Eq. (1) is solved for given the previous results (VE, VB, VBX, VCX, VEX). The VE in Eq. (2)
is solved for newer given VC and (VB, VBX, VCX, VEX). The VB in Eq. (3) is solved for newer given (VC, VE)
and (VBX, VCX, VEX). We have similar procedure for other unknowns [13,17,19,20]. We use the DC solutions and previous calculated results as the starting point for all unknowns to be solved at present time step t . This method converges globally and is not sensitive to the initial guesses in the solution process [19].
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 above ODEs as the following form
dVX
Fig. 3. The maximum norm error versus number of outer iterations. All computed unknowns have the same strictly convergent behavior.
Table 2
A comparison of simulation CPU time between our developed circuit simulator and HSPICE
Number of points per cycle 200 400 800 1600 CPU time of our method (sec.) 2670 5450 11,030 22,600 CPU time of HSPICE (sec.) 2490 6020 14,570 27,020
where VX is the unknowns to be solved and f is the collection of the nonlinear functions. We define the MI
parameter λ=∂V∂f
X [13,17–19] and inserted the λ into Eq. (20), then we have the MI equation dVX
dt = f (η, t) − λ(VX− η), (21)
where v0 η w0 is a value in [0, T ], and v0 and w0 are the lower and upper solutions of Eq. (20),
respectively [18,19]. Based on the nonlinear behaviors of each decoupled circuit ODEs, mathematically we have shown that the solution algorithm has monotone convergence property [19]. As shown in Fig. 2, we solve Eq. (21) with respect to each unknown using the RK method, where the notation subscript o is the previous iteration results, and the superscripts (t− 1) and (t) represent the previous and present time steps. The nonlinear functions Fi(.), i= 1, . . ., 6 are the corresponding nonlinear terms of the ODEs above, where the TOL is the error tolerance,
h is the increment in the RK formula of fourth order [24].
Fig. 3 shows the achieved convergence properties (maximum norm error versus number of outer iterations) on the proposed method for the HBT simulation, where the VIN= 1.42 V and VCC= 5 V. In our calculation expe-rience, a convergence criterion for all quantities (maximum norm error < 10−10 and 10−7 for the inner and the overall outer loops) can be reached by only 8−12 MI loops and 25−30 outer loops, respectively. The simulation time is an important consideration in circuit simulation. Under the same convergence criterion above, a comparison of the simulation CPU time between our developed simulator and HSPICE simulator is performed on a PC-based Linux workstation (CPU is with Pentium-IV 1 GHz and RAM= 512 MBytes). In the comparison, 40 cycles are computed, the input amplitude of two-tone signal Vm= 0.002 V, and the frequencies f1 and f2 are 1.71 and
1.89 GHz, respectively. As shown in Table 2, we find the simulation time of two simulations is in the same order magnitude. For a fixed convergence criterion, simulation time of our method increases linearly when the number of simulation points per cycle are increased. Furthermore, our method can be easily implemented on a Linux-cluster
3. Results and discussion
The DC characteristics of the HBT circuit shown in Fig. 1 is performed firstly, and the computed results are inputted as starting points in high frequency simulation. In our simulations, computed family DC curves with our method demonstrate very good consistence with the results from the well-known HSPICE simulation.
For the nonlinearity analysis of the equivalent HBT circuit, we calculate the two-tone input signal in the time domain with intermodulation. The simulation was done directly in the time domain relies on the robustness of the developed circuit simulator with the MI method. From the computed periodic output results, we found the output signal has good repetition. Fig. 4 shows the calculated output voltage (VOUT) in the time domain with a
Fig. 4. The VOUT(t) versus the time simulated with our developed simulator.
Fig. 6. The output power versus frequency. The results are calculated with our time domain results and FFT.
Fig. 7. The frequency domain results from the HSPICE simulation. The IM3 products cannot be identified from this calculation.
two-tone input signal. As shown in Fig. 1, the input signal VIN is the DC bias and the expression of Vin is as follows: Vin= Vmsin(2πf1t)+Vmsin(2πf2t), where the Vm= 0.002 V, and the frequencies f1and f2are 1.71 and
1.89 GHz, respectively. Figs. 4 and 5 are our results and the HSPICE results, respectively. Contrary to the HSPICE results initially have some unstable calculations, our simulator presents its robustness in the large signal time domain analysis. We convert the time domain results into the frequency with the FFT directly. Figs. 6 and 7 are the corresponding results with Figs. 4 and 5, respectively. We find the products of the third-order intermodulation (IM3) at 2f1−f2and 2f2−f1are clearly observed in Fig. 6. However, as shown in Fig. 7 it is difficult to identified the two
IM3 products. Our methodology for large scale time domain analysis and two-tone intermodulation demonstrates its superiority over some approaches. Compared to conventional NI-based circuit simulation methodology (e.g.,
A novel circuit simulation based on the waveform relaxation and monotone iterative methods has been proposed. With this approach, the DC and two-tone AC characteristics for an HBT circuit were directly computed from time domain results. Compared with the HSPICE simulator, our solution method for an HBT high frequency simulation has been presented to show the accuracy and efficiency of the method. This method provides an alternative for the time domain solution of circuit ODEs. This solution methodology can be generalized for larger circuit simulation including more transistors.
Acknowledgements
This work was supported in part by the National Science Council of Taiwan, under contract number NSC 91-2112-M-317-001 and the 2002 Research Fellowship Award of the Pan Wen-Yuan Foundation in Taiwan.
References
[1] S.M. Sze, Physics of Semiconductor Devices, 2nd edn., Wiley-Interscience, New York, 1981. [2] M. Yanagihara, et al., Solid-State Electron 41 (1997) 1615.
[3] S. Sheu, et al., IEEE Trans. Electron Devices 45 (1998) 326. [4] T. Oka, et al., IEEE Trans. Electron Devices 48 (2001) 2625.
[5] K.S. Kundert, J.K. White, A. Sangiovanni-Vincentelli, Steady-state Method for Simulating Analog and Microwave Circuit, Kluwer Academic Publishers, 1990.
[6] B. Li, S. Prasad, IEEE Trans. Microwave Theory Tech. 45 (1997) 1135.
[7] B. Troyanovsky, Z. Yu, R.W. Dutton, Comput. Methods Appl. Mech. Engrg. 181 (2000) 467. [8] M. Gourary, et al., Comput. Methods Appl. Mech. Engrg. 181 (2000) 451.
[9] G.F. Carey, et al., Circuit, Device and Process Simulation—Mathematical and Numerical Aspects, Wiley, New York, 1996. [10] Ron M. Kielkowski, SPICE: Practical Device Modeling, McGraw-Hill, New York, 1995.
[11] R.C. Jaeger, Computerized Circuit Analysis Using Spice Programs, McGraw-Hill College Div, 1997. [12] Y. Li, et al., Comput. Phys. Commun. 142 (2001) 285.
[13] Y. Li, WSEAS Trans. on Systems 1 (2002) 68. [14] Y. Li, et al., Comput. Phys. Commun. 147 (2002) 697.
[15] Y. Li, et al., in: Tech. Proc. 5th Internat. Conf. Modeling and Simulation of Microsystems, Puerto Rico, 2002, p. 376.
[16] K.-Y. Huang, et al., in: D. Tsoukalas, C. Tsamis (Eds.), Simulation of Semiconductor Processes and Devices, Springer, New York, 2001, p. 226.
[17] Y. Li, et al., in: Proc. IEEE Internat. Symp. of VLSI Tech., Systems, and Applications, Taipei, 1999, p. 27.
[18] V. Lakshmikantham, A.S. Vatsala, Generalized Quasilinearization for Nonlinear Problems, Kluwer Academic Publishers, Boston, 1998. [19] Y. Li, WSEAS Trans. Math. 1 (2002) 159.
[20] E. Lelarasmee, IEEE Trans. CAD. IC Systems 1 (1982) 131.
[21] P.S. Pacheco, Parallel Programming with MPI, Morgan Kaufmann, San Francisco, 1997. [22] W. Liu, Handbook of III–V Heterojunction Bipolar Transistors, John Wiley & Sons, 1998. [23] L.E. Getreu, Modeling the Bipolar Transistor, Elsevier, New York, 1984.
