Numerical modeling of MOS transistor with interconnections using lumped element-FDTD method

Microelectronics Journal 43 (2012) 995–1002

Contents lists available at SciVerse ScienceDirect

Microelectronics Journal journal homepage: www.elsevier.com/locate/mejo

Numerical modeling of MOS transistor with interconnections using lumped element-FDTD method Samir Labiod a,n, Saida Latreche a, Christian Gontrand b a b

´partement d’Electronique, Universite´ Mentouri, Constantine, Algeria Laboratoire Hyperfre´quences et Semiconducteurs (LHS), De Universite´ de Lyon, INSA-Lyon, INL, CNRS-UMR5270, Villeurbanne F-69621, France

a r t i c l e i n f o

abstract

Article history: Received 17 December 2011 Received in revised form 9 June 2012 Accepted 2 October 2012 Available online 24 October 2012

A detailed three-dimensional (3-D) full-wave time-domain simulation model is presented for the analysis of an active semiconductor device. Microwave analysis of the nMOS transistor using monolithic microwave integrated circuit (MMIC) is simulated. A three-dimensional finite-difference time-domain (3D-FDTD) scheme is adopted to describe the circuit passive part, whereas numerical device simulation techniques are employed for the active semiconductor devices. The proposed algorithm provides the time and space distribution of the unknown functions electrostatic potential, carriers’ concentration, and current density for the MOS transistor. Fast Fourier transform (FFT) analysis is performed on the obtained time response to calculate the scattering matrix of the S-parameters and phase frequency dependencies. Finally, the obtained results are presented and show a good agreement with those based on a hydrodynamic model. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Mos transistor Implicit backward Euler Finite difference time domain (FDTD) S-parameters

1. Introduction Electromagnetic propagation and radiative effects on semiconductor active devices become more and more important: most notably, signal reflections due to interconnecting line discontinuities, dispersion, and crosstalk phenomena. Then, active circuit design at microwave domain requires considering the electromagnetic coupling effects. This requirement can be satisfied by considering a full-wave approach, which solves Maxwell’s equations while taking into account the interaction between electromagnetic wave and active device. To this aim, FDTD method [1] is applied. FDTD method is successfully applied to the microwaves structure and optimization of linear devices (resistor, inductor and capacitor) [2]. Also, in recent years, nonlinear devices (semiconductor devices) have been analyzed at microwave frequencies, and several methods for incorporating nonlinear devices into FDTD algorithm have been proposed [3]. In this paper, we present an efficient time-domain full-wave approach for modeling an active device (2-D MOS transistor) considered for millimeter-wave applications. The proposed model couples a three-dimensional (3-D) timedomain solution of Maxwell’s equations to an active device model. In, this work, all algorithms are implemented using the MATLAB software.

n

Corresponding author. E-mail address: [email protected] (S. Labiod).

0026-2692/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.mejo.2012.10.005

The active device model is based on the drift-diffusion model (DDM) which consider carrier’s transport and Shockley–Read– Hall (SRH) model [4]. Two-dimensional finite-different (2D-FD) algorithm is proposed to solve the active device model and full implicit backward Euler’s scheme has been used to accomplish time domain. The obtained implicit matrix systems are solved by using an efficient and accurate algorithm based on the wellknown Gummel’s iterations [5]. The active devices in the microwave circuits is typically very small in size compared to a wavelength, then it can be modeled by its equivalent lumped device with a very high degree of accuracy. Thus, in the conventional lumped element-FDTD (LE-FDTD) approach, three contacts for the active device (drain, source and substrate) are considered as current sources that interact with the Maxwell’s equations, exactly with Maxwell– Ampere’s equation [6,7]. The 3D-FDTD scheme is adopted to describe the circuit passive part: interconnections and via holes. Uniaxial perfectly matched layer (UPML) conditions is used to truncate the computational domain in such way that there are no wave reflections from the boundaries. We consider a wideband Gaussian pulse as an excitation, the measured network parameters in time-domain are transformed into frequency-domain by directly using the Fast Fourier transform (FFT) technique [8]. To confirm the validity of the proposed model for the active device, we have compared it with some experimental results which were found in other works. Results of steady state regimes are presented and compared with other results based on different mathematical model. S-parameter results were carried out in the frequency range of

996

S. Labiod et al. / Microelectronics Journal 43 (2012) 995–1002

0.1–50 GHz. Finally, the calculated results for the microwave circuit are discussed and compared with those based on the hydrodynamic model.

Holes Electrons

2. Active device model The semiconductor models used are usually obtained by applying approximations and simplification of the Boltzmann transport equation (BTE). These assumptions can result in a number of different transport models such as the drift diffusion model (DDM). This later formulates the problem using three dependent variables V, n, and p. Poisson’s equation relates the electrostatic potential to the space charge density: q

r2 V ¼  ðpðx,yÞnðx,yÞ þDopðx,yÞÞ e

ð1Þ

The charge conservation equations formulated for electrons and holes, respectively: @n 1 !!  r J n ¼ r SRH @t q

ð2Þ

@p 1 !! þ r J p ¼ r SRH @t q

ð3Þ

where V, n and p denote the electron density, hole density and electrostatic potential respectively. From Eq. (1), Dop accounts for the net ionized impurity concentration, and rSRH represents the Shockley–Read–Hall recombination, which is a general recombination process using traps in the forbidden band gap of the semiconductor. npn2i r SRH ¼ tp ðn þni Þ þ tn ðp þ ni Þ

Table 1 Constants values for mobility model.

ð4Þ

where q is the electronic charge, ni is the intrinsic carrier concentration and tn, tp are the carrier lifetimes of electron and hole respectively. ! ! Electron and hole currents densities J n and J p are expressed here through the ‘‘drift-diffusion’’ considerations respectively: ! ! ! J n ¼ qnmn r V þ qDn r n

ð5Þ

! ! ! J p ¼ qpmp r VqDp r p

ð6Þ

Moblities mn, mp and diffusion coefficients Dn, Dp of electrons and holes are related by: ( Dn ¼ mn U T ð7Þ Dp ¼ mp U T where UT ¼0.026 V and present the thermal voltage at temperature T ¼300 K. The variation of the mobility (m) with the electric field (E) and the ionized impurity concentration (Dop) is approximated by the following expression [5]:  2 !  2    2 E=A m0 Dop E þ ¼ 1þ ð8Þ þ Dop=S þ N B m E=A þ F

2.1. Boundary conditions for the active device The boundary potential at an ohmic contact is the sum of the potential VA, which is externally applied to the contact, and the so

N

S

A

F

B

480 1400

4  1016 3  1016

81 350

6.1  103 3.5  103

1.6 8.8

2.5  104 7.4  103

called built-in-potential, which is produced by the doping:   DopðsÞ DopðsÞ  log V ðsÞ ¼ V A þU T  ni DopðsÞ Dirichlet boundary conditions for electrons and holes densities can be obtained by assuming charge neutrality and thermal equilibrium. 8 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 > > > nðsÞ ¼ ðDopðsÞ þ DopðsÞ2 þ4ni 2 Þ < 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 > 2 2 > > : pðsÞ ¼ 2 DopðsÞ þ DopðsÞ þ 4ni where S indicates the value specified along the boundary surface. Neumann boundary conditions is used at free surfaces, the normal derivatives of electrostatic potential, electron and hole densities are set to zero. ( @V @N ¼ 0 @n @N

¼

@p @N

¼0

where N indicates the normal with respect to the boundary surface.

3. Electromagnetic model Formulation of the FDTD method begins by considering the differential form of Maxwell’s equations that govern the propagation of fields in the computational domain. These equations are first-order linear coupled differential equations relating the elec! ! tric field E , magnetic field H and current densities at any point in space and at any time. Maxwell’s curl equations may be written as: ! ! ! dH : ð9Þ r  E ¼ m dt !

!

r  H ¼e

! ! dE ! þ J media þ J Device : dt

ð10Þ

! Here e is the dielectric constant and J media accounts for the contribution of the current flowing along the distributed media.

4. Coupling the two models The coupling between the active device and the electromagnetic model can be done using LE-FDTD method. The discretization uses a first and second order of (1) in 2-Dfinite difference mesh, leads to have V t ði þ 1,jÞ þ V t ði1,jÞ þ V t ði,j þ1Þ þV t ði,j1Þ2V t ði,jÞ ¼

The values of m0, N, S, A, F and B are presented in Table 1.

m0

2 qh 

e

nt ði,jÞpt ði,jÞDopði,jÞ



ð11Þ

Euler implicit method [9,10] seeks to approximate the derivatives in (2) and (3) with regard to the discrete solutions points defined by spatial and temporal cells. Then the carrier continuity equations may be discretized in implicit form as follows: a1 ði,jÞnt þ 1 ði1,jÞ þ a2 ði,jÞnt þ 1 ði þ 1,jÞ þ a3 ði,jÞnt þ 1 ði,j1Þ

S. Labiod et al. / Microelectronics Journal 43 (2012) 995–1002

þ a4 ði,jÞnt þ 1 ði,j þ 1Þ½a5 ði,jÞ þ 1nt þ 1 ði,jÞ ¼

tuT h2 m0  Rði,jÞnt ði,jÞ 2 DtLD mn

Etyþ 1 ði,j,kÞ ¼ Ety ði,j,kÞ þ ð12Þ

a6 ði,jÞpt þ 1 ði1,jÞ þ a7 ði,jÞpt þ 1 ðiþ 1,jÞ þa8 ði,jÞpt þ 1 ði,j1Þ þ a9 ði,jÞpt þ 1 ði,j þ 1Þ½a10 ði,jÞ þ 1pt þ 1 ði,jÞ ¼

tuT h2 m0  Rði,jÞpt ði,jÞ DtL2D mp

ð13Þ

a1  a10 are the matrix coefficients which depend with the Vt vector [11]. The space step (h) and time step (Dt), are limited by Debye length (LD) and dielectric relaxation time (t) respectively. At each time step, Eqs. (11), (12) and (13) are solved using Gummel’s iterations [12]. Poisson equation is solved at all grid points, followed by electron continuity equation, and then by hole continuity equation, each matrix system is solved using the successive over-relaxation (SOR) method [13]. To preserve the computational efficiency of 3D-FDTD solution scheme, we assume that device equations are solved at even time-steps; this makes the averaging of current density contribution between time steps t and tþ1. The explicit FDTD algorithm is obtained using the centered difference approximation on both the time and space first-order partial differentiations of all the components of electric and magnetic fields. The electric field components are computed at ‘‘integer’’ time-steps and the magnetic field at ‘‘half-integer’’ time-step [14]. We consider a locally uniform medium in the vicinity of an x–y directed lumped element Fig. 1. A lumped component is assumed to coincide with the E field components. The current along the element and the rate of change of the E field across the element determines the values of the H field components surrounding the E field. The discrete form of such an additional equation is found by discretizing and projecting (10) along the x and y directions Etxþ 1 ði,j,kÞ ¼ Etx ði,j,kÞ þ 



997

Dt h

eDz

   i Hxt þ 1=2 i,j,kþ 1=2 Htxþ 1=2 i,j,k1=2

   i Htz þ 1=2 i þ 1=2,j,k Hzt þ 1=2 i1=2,j,k    Dt t þ 12 1 tþ1 1 JDrain,y i,j þ ,k J media,y ði,j,kÞ þ  2 2 e   1 ð15Þ þJ tDrain,y i,j þ ,k 2 ! where Dx, Dy and Dz are spatial discretization and J Drain,Source is an impressed current density through which the MOS transistor will be incorporated. In our model, such a contribution comes from the time domain solution of the related device described by the model’s equations [15]. The active device contacts are oriented in y and x directions, where the currents densities are given by: 8 t > < JtDrain,y ¼ DIDrain zDx ð16Þ It > : J tSource,x ¼ DSource z Dy 

Dt h

eDx

The voltage drop DV across the MOS transistor (Gate contact) is simply obtained at (t þ1) time step by integrating the electric field Ey component across the y direction.   1 ð17Þ DV Gate ¼ DyEyt þ 1 1,j þ ,k : 2 The numerical algorithm for Maxwell’s curl equations and semiconductor equations requires that the time increment Dt must have a specific bound relative to the spatial discretization Dx, Dy, and Dz. The time discretization interval, Dt is determined by the Courant condition [16], given by 1

Dt r qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C

1 ðDxÞ2

þ ðD1yÞ2 þ ðD1zÞ2

ð18Þ

where C is the velocity of light in considered medium. The whole numerical procedure to describe the coupling between the two models can eventually be summarized as follows:

  i Dt h t þ 1=2  Hz i þ1=2,j,k Htz þ 1=2 i1=2,j,k eDx

  i Dt h t þ 1=2  Hy i,j þ1=2,k Htyþ 1=2 i,j1=2,k e:Dy

     1 tþ1 1 1 t þ 12 J Source,x i þ ,j,k þJ tSource,z i þ ,j,k J media,x ði,j,kÞ þ 2 2 2 e ð14Þ

Dt

5. Results and discussions The device considered in this analyze is shown in Fig. 2 The active device is a MOS transistor with a silicon (Si) substrate doping with 2  1016 At.cm  3. Source and drain junctions are doped as Gaussian function with a maximum and minimum densities of 2  1020 and 2  1018 At.cm  3, respectively. The work function difference between the gate and the silicon interface is jms ¼0.2 V corresponding an aluminum contact metal. In this simulation all calculations were done at temperature T¼300 K. The dimensions of the MOS transistor are given as:

Fig. 1. Yee’s FDTD cell by incorporating lumped element.

– Oxide thickness TOX ¼5 nm. – Source-drain separation Leff ¼0.5 mm. – Drain and source junctions depth, Xj ¼0.05 mm.

998

S. Labiod et al. / Microelectronics Journal 43 (2012) 995–1002

Fig. 2. Schematic representation of a bidimensional MOS transistor.

Table 2 Numerical methods comparisons for the active device model resolutions. Method

Iterations number

Residual potential error

CPU time (s)

NR SOR MINRES (minimal residual) LU decomposition

4600 3500 2822

2.23 10  5 1.06 10  5 1.14 10  5

125.36 93.04 109.78

126

1.35 10  5

372.12 s

Fig. 4. Calculated potential distribution.

Fig. 5. Electron density distribution in logarithmic scale.

Fig. 3. Potential error in logarithmic scale using the SOR method.

The 2D numerical simulation uses a non uniform mesh with 170  70 nodes. 5.1. Active device simulation The device is biased with a direct potential of VG ¼0.8 V and VD ¼1 V. The DC parameter distributions are obtained by solving the drift-diffusion model where time derivatives for the three unknown parameters are directly set to zero. Different numerical methods are presented and compared in Table 2. We notice that some numerical methods converge for a little iteration however takes the maximum CPU time, whereas other numerical methods takes the minimum CPU time for several iterations. One of the original reasons is that in the iterative methods likes NR (Newton Raphson) and SOR the matrix requires less storage in memory than direct methods likes LU decomposition. Fig. 3 presents the norm of the residual error for the variations in potential vector using the SOR method. A good convergence of

Fig. 6. Contour plot of current density in logarithmic scale.

the proposed method appears for a little time. In this simulation, each iteration takes approximately 0.0265 s. Fig. 4 shows the potential distribution inside the device. Fig. 5 shows the computed electrons and holes density profile in the entire device in logarithmic scale. Important quantity of electrons which attracted at oxide–semiconductor interface can be observed. The figure also shows the onset of the pinch-off effect which indicates the lack of channel region near the drain

S. Labiod et al. / Microelectronics Journal 43 (2012) 995–1002

999

Fig. 7. Calculated drain current characteristics: (a) output characteristic and (b) transfer characteristic.

Fig. 10. Calculated transient currents. Fig. 8. Calculated I–V characteristics for the proposed device.

Fig. 9. Applied gate bias versus time. Fig. 11. Electron charge quantity variation in the channel for different time values.

junction. Contour plot of current density lines with 0.65 steps in logarithmic scale is presented in Fig. 6. Notice the correct flow of electrons which injected from source into the drain contacts. The current density matrix has been obtained by using Eqs. (5) and (6). The static current-potential characteristics for the n-channel MOS transistor are finally presented in Fig. 7. For a given bias, drain current is calculated by integrating current density over drain contact.

The obtained result for the I–V characteristic of the MOS transistor is compared with other simulator which is based on the hydrodynamic model. On the other hand, the proposed model is compared with other measured characteristics. Therefore we use the same parameters which are found in [17]. The device is0.35 mm. Drain and source contacts have maximum, minimum and substrate doping in the order of 1020, 1018, and 1016 At./cm3,

1000

S. Labiod et al. / Microelectronics Journal 43 (2012) 995–1002

respectively. The oxide thickness was 0.0096 mm, and source/ drain junction depth of approximately 0.17 mm. The deduced values presented in Fig. 8 agree well with the independent measurement presented in this work. Transient simulations were simulated at Vd ¼1 V with Vg switched from 0 V to 0.8 V with a non linear variation of duration 0.01 ps (see Fig. 9). Drain, source and substrate current variations versus time are shown in Fig. 10. When the channel is clearly formed, the stability of

the drain and source current appears above 10 ps whereas substrate current above 1 ns. We can clearly show that the transient solution converge towards the results obtained in the steady state solution and by this we can confirm the validity of the proposed method. Fig. 11 shows the variation of electron concentration in the channel–oxide interface versus time, we note that electrons are injected from substrate into the channel region. These charges do not recombine immediately which results in the instability of currents for the three contacts—drain, source and substrate.

5.2. Microwave circuit simulation A simulation example is introduced in Fig. 12. At high frequency, the analysis of the MOS transistor is considered by introducing the interconnections effect. So the proposed circuit consist of a MOS transistor, connections (stripe lines), via holes and two load resistors The 3D microwave circuit would be built using the monolithic microwave integrated circuit (MMIC) technology with 100 mm thick substrate and 3.9 dielectric constant. The proposed microwave circuit was designed as follows:

Fig. 12. The simulated structure.

Fig. 13. Electric field distribution for 700 time steps (a, b).

– The active device (MOS transistor) is supposed like a two current sources oriented following x and y directions. The source contact is related with ground plane (commune source).

Fig. 15. S-parameters of nMOS transistor from 50 MHz to 6 GHz.

Fig. 14. Frequency analysis of nMOS transistor from 0.1 GHz to 50 GHz: (a) S-parameters and (b) phase analysis of the output wave at drain contact.

S. Labiod et al. / Microelectronics Journal 43 (2012) 995–1002

1001

Fig. 16. Set of calculated S-parameters for the microwave circuit.

– Metallic strips are used to interconnect the circuit components, these interconnections are modeled as transmission lines with 35 mm thickness. – Source contact and load resistor are connected to ground plane by using via-holes in the proposed circuit. – Input and output ports for RF signals were added and connected in order to have 50 O transmission lines. – The excitation pulse that is used at the input port is a Gaussian pulse with a maximum frequency of 50 GHz.

LE-FDTD simulation is performed with a uniform grids of space steps Dx¼0.189 mm, Dy¼ 0.2 mm and Dz¼0.02 mm. To take, propagation along an infinite medium, the proposed circuit has been surrounded by adopting uniaxial perfectly matched layer (UPML) boundary conditions [18]. Perfect electric condition (PEC) is applied to define the metallic surfaces and viaholes. The active device is biased with a direct potential of Vg ¼0.8 V and Vd ¼1 V. A view of the computed electric field vertical component is shown in Fig. 13 (the normal component of the electric field with respect to the circuit surface), from which the nonuniform potential distribution along the conductive lines can be clearly observed. S parameters and phase analysis for the nMOS transistor are obtained and displayed in Fig. 14 for the interval frequency (0.1–50 GHz). The scattering matrix has been extracted from the time responses for gate and drain contacts. To estimate the S parameters for the microwave circuit, the output port 2 has been selected (open) [19], when a Gaussian pulse has been input at port 1. Signal frequency spectra have been extracted, through FFT analysis, from the transient response predicted by time-domain analysis. S parameters for the active device displayed in Fig. 15 are compared with experimental measure which are found in [20] for the interval frequency (50 MHz to 6 GHz). Here the channel length is 0.25 mm and the device is biased at Vg ¼ Vd ¼ 1.2 V. Fig. 16 shows the reflection coefficient and transmission coefficient computed at the two active ports (S11 and S21). To validate such results, the same structure has been simulated using the two commercial simulators (HFSSþISE-TCAD). In this case, we have extracted the S-parameters matrix device from ISETCAD simulator which is based on the hydrodynamic model, whereas the passive circuit has been analyzed by using the HFSS simulator software. The slight differences between LE-FDTD and finite element simulators outcome can hence be ascribed to the different approaches such as method resolution, and absorbing boundary conditions type.

6. Conclusion In this work, we have presented an efficient time-domain fullwave simulator for modeling an nMOS transistor used for millimeterwave applications. In the first part, we have developed a two-dimensional timedependent simulator that solves the three semiconductor equations using backward Euler scheme for a MOS transistor. In order to incorporate the active device into the electromagnetic model, drain, source and substrate currents versus time are calculated and presented. Drift-diffusion model and Maxwell’s equations are combined in the 3D-LE-FDTD algorithm. The scattering matrix of the S-parameters and phase frequency dependencies are extracted from the FDTD algorithm using the FFT method. The simulation results have brought a number of interesting findings which are summarized as follows: – From the transient simulation, the current stability of the active device appears above 10 ps. – At high frequency, connections between components (metallic strips) and via holes should be considered for the electromagnetic simulation. – The critical point frequency 12.3 GHz of the microwave circuit, present a brutal variation in phase and a minimum gain. – The simulation confirms that the proposed active device can be used in the high frequency domain. The obtained results were presented and show a good agreement with those obtained using commercial software (ISE-TCAD and HFSS software) which is based on the hydrodynamic model. The model presented in this paper represents an interesting contribution to the full-wave characterization of active device. Also, the proposed analysis can be extended to more complex microwave circuits.

Acknowledgments This work is supported by UPM (Union Pour la Me´diterrane´e). References [1] K.S. Yee, Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media, IEEE Trans. Antennas Propag. (1966) 302–307. [2] X. Wenbin, M. Laiqing, Y. Wenhua, M. Raj, A novel expression of lumped element in FDTD simulation, Int. J. Infrared Millimeters Waves (2006) 985–993.

1002

S. Labiod et al. / Microelectronics Journal 43 (2012) 995–1002

[3] M. Picket-May, A. Taflove, J. Baron, FDTD-modeling of digital signal propagation in 3-D circuits with active and passive loads, IEEE Trans. Microwave Theory Tech. (1994) 1514–1523. [4] S.M. Sze, Kwok K.Ng, Physics of Semiconductor Devices, John Wiley & Sons, Inc., Hoboken, New Jersey, 2007. [5] D.L. Schafetter, H.K. Gummel, Large-signal analysis of a silicon read diode oscillator, IEEE Trans. Electron Devices (1969) 64–77. [6] C.C. Wang, C.W. Kuo, An efficient scheme for processing arbitrary lumped multiport devices in the finite-difference time–domain method, IEEE Trans. Microwave Theory Tech. (2007) 958–965. [7] O. Gonzalez, J.A. Pereda, A. Herrera, A. Vegas, An extension of the lumpednetwork FDTD method to linear two-port lumped circuits, IEEE Trans. Microwave Theory Tech. (2006) 3045–3051. [8] X. Zhang, K.K. Mei, Time–domain finite difference approach to the calculation of the frequency-dependent characteristics of microstrip discontinuities, IEEE Trans. Microwave Theory Tech. (1988) 1775–1787. [9] J.C. Butcher, Numerical Methods for Ordinary Differential Equations, 2rd Edition, John Wiley & Sons Ltd., England, 2008. [10] D. Vasileska, S.M. Goodnick, Computational Electronics, 1st Edition, Morgan & Claypool, USA, 2006. [11] R. Mirzavand, A. Abdipour, G. Moradi, Full-wave semiconductor devices simulation using ADI–FDTD method, Prog. Electromagn. Res. M (2010) 191–202.

[12] A. Quarteroni, R. Sacco, F. Saleri, Numerical Mathematics, Springer-Verlag press, New York, 2000. [13] A. Quarteroni, R. Sacco, F. Saleri, Numerical Methods, Springer-Verlag Press, Milano, 2007. [14] D.M. Sullivan, Electromagnetic Simulation Using the FDTD Method, IEEE Press Series on RF and Microwave Technology, New York, 2000. [15] P. Ciamolini, L. Roselli, G. Stopponi, Mixed-mode circuit simulation with fullwave analysis of interconnections, IEEE Trans. Electron Devices (1977) 2098–2105. [16] A. Taflove, S. Hagness, Computational Electrodynamics: The Finite-Difference Time–Domain Method, 3rd Edition, Artech House Inc., 685 Canton Street Norwood, 2005. [17] W. Liang, N. Goldsman, I. Mayergoyz, J. Oldiges, 2-D MOSFET modeling including surface effects and impact ionization by self-consistent solution of the Boltzmann, Poisson, and Hole-continuity equations, IEEE Trans. Electron Devices (1997) 257–267. [18] J. Li, C. Miao, A new implementation of the uniaxial perfectly matched layer, Microwave Millimeter Wave Technol. (2008) 770–773. [19] X. Ye, J.L. Drewnaik, Incorporating two-port networks with S-parameters into FDTD, IEEE Trans. Microwave Wireless Components (2001) 77–79. [20] S.D. Wu, G.W. Huang, L.P. Cheng, C.Y. Chang, Characterization of 2-port configuration mosfet’s amplifiers by 4-port measurement, in: Proceedings of the Asia-Pacific Microwave Conference, (2003).