Modeling Emerging Semiconductor Devices for Circuit Simulation

Circuit simulation is an indispensable part of modern IC design. The significant cost of fabrication has driven researchers to verify the chip functionality through simulation before submitting the design for final fabrication. With the impending end of Moore ’ s Law, researchers all over the world are looking for new devices with enhanced functionality. A plethora of promising emerging devices has been proposed in recent years. In order to leverage the full potential of such devices, circuit designers need fast, reliable models for SPICE simulation to explore different applications. Most of these new devices have complex underlying physical mechanism rendering the model development an extremely challenging task. For the models to be of practical use, they have to enable fast and accurate simulation that rules out the possibility of numerically solving a system of partial differential equations to arrive at a solution. In this chapter, we show how different modeling approaches can be used to simulate three emerging semiconductor devices namely, silicon- on-insulator four gate transistor(G 4 FET), perimeter gated single photon avalanche diode (PG-SPAD) and insulator-metal transistor (IMT) device with volatile memristance. All the models have been verified against experimental /TCAD data and implemented in commercial circuit simulator.


Introduction
Circuit simulation is an indispensable part of modern IC design. The significant cost of fabrication has driven researchers to verify the chip functionality through simulation before submitting the design for final fabrication. With the impending end of Moore's Law, researchers all over the world are now looking for new devices with enhanced functionality and a plethora of promising emerging devices has been proposed in recent years. To leverage the full potential of such devices, circuit designers need fast, reliable models for SPICE simulation to explore different applications. Most of these new devices have complex underlying physical mechanism rendering the model development an extremely challenging task. To be of practical use, the model has to enable fast and accurate simulation, which rules out the possibility of numerically solving differential equations underlying the physics of semiconductor devices to arrive at a solution. In this chapter, we show how different approaches can be adopted to model three emerging semiconductor devices, namely silicon-on-insulator four-gate transistor (G 4 FET), single-photon avalanche diode (SPAD), and insulator-metal transistor (IMT) device with volatile memristance.

Background
The amazing technological advancements in the semiconductor industries during last 60 years have largely shaped the modern world we live in. However, bulk silicon devices are now faced with some fundamental physical limits and innovative researchers and scientists all over the world have been diligently exploring new devices and computing paradigms to continue the breathtaking pace of technology advancement. The widespread use of a technology in circuit design is heavily dependent upon good SPICE models for CAD (computer-aided design) tools that are now ubiquitous in circuit design. Sophisticated models for existing semiconductor devices integrated with CAD tools have enabled designers worldwide to design innovative and reliable circuits which are in a large part responsible for the technology boom of the last several decades. In this chapter, we will be discussing three relatively new devices highlighting the works and challenges involved in model development and SPICE simulation.
The first device is SOI (silicon-on-insulator) four-gate transistor known as G 4 FET. Here, we outline its operating mechanism, different approaches towards modeling and highlight two particular approaches, namely multivariate cubic spine interpolation model and MOS-JFET macromodel. The second device is a single photon avalanche diode (SPAD). We review and discuss the most recent developments in SPAD modeling with particular emphasis on a new class of SPAD known as PGSPAD (perimeter gated SPAD). The third device is insulator metal transition (IMT) device. Here, we focus on a very recent IMT model to demonstrate how proper simplification of device physics can lead to a suitable SPICE model facilitating circuit design.
3. Silicon-on-insulator (SOI) four gate field effect transistor (G 4 FET) G 4 FET was first reported in 2002 [1]. It has four independent gates, two of which provide vertical MOS (metal-oxide-semiconductor) field-effect action whereas the other two gates provide lateral junction field effect transistor (JFET) functionality. The unique G 4 FET structure can be leveraged to design circuits for different analog, mixed-signal and digital applications with significantly reduced transistor counts. Some of these have already been experimentally demonstrated including LC oscillators and Schmitt trigger circuit with adjustable hysteresis using negative differential resistance [2], four quadrant analog multipliers [3], adjustable threshold inverters, real time reconfigurable logic gates and DRAM cell [4], and universal and programmable logic gate capable of highly efficient full adder design [5]. G 4 FET inspired multiple state electrostatically formed nanowires have already been used for threshold logic functions [6] and high-sensitivity gas sensing [7]. Another potential application is to leverage the tunable tent map shape characteristics of the voltage controlled G 4 NDR circuit [2] to build chaotic logic gates for hardware security applications [8,9].

Device structure
One particular advantage of G 4 FET is that standard SOI process without significant modification can be used to fabricate these devices. The carrier conduction through the channel is modulated using four independent gates. The conventional source and gates become junction gates (JG1 and JG2) to provide lateral JFET functionality and the top and bottom gates provide vertical MOS functionality. Figure 1 shows the 3-D schematic structure of a p-channel G 4 FET. It is evident that no specialized fabrication procedure is necessary for this device.

Conduction mechanism
By putting two highly doped body contacts on the opposite sides of the traditional channel, we can convert a conventional n-channel SOI MOSFET into a pchannel G 4 FET. What used to be n + doped source and drain now act as lateral JFET-like gates. These gates can modulate the width of the conduction channel according to applied bias. The functionality of the conventional MOS gate is achieved through the top oxide gate. The substrate along the buried oxide provides an additional bottom-gate functionality. A high positive voltage on top gate creates an inversion layer of electrons near the top surface. As we reduce the top gate voltage, the top surface gradually changes from inversion to depletion, and eventually, with high negative voltage, it will enter the accumulation region. Similar effect due to bottom-gate can be observed for the surface near the buried oxide interface. In order to ensure ohmic contact with the conducting channel, body contacts on the opposite side of the channel are highly doped. Through this process, an inversionmode, n-channel MOSFET can be transformed into an accumulation/depletionmode p-channel G 4 FET. A similar mirror transformation will convert a classic SOI p-channel MOSFET into an n-channel G 4 FET. The current from drain to source follows in a direction perpendicular to a traditional MOSFET. Similarly, the geometrical features are also swapped, i.e., the length and width of the classic MOSFET now become the width and length of the new G 4 FET, respectively. As outlined here, there are three possible paths for current conduction, namely (1) top surface near gate oxide interface, mostly controlled by top gate bias, (2) bottom surface near buried oxide interface, mostly controlled by bottom gate bias, and (3) volume conduction away from vertical oxide interfaces, mostly controlled by lateral gates.

State-of-the-art G 4 FET models
Different approaches have been adopted to understand and model the operation of G 4 FET. Extraction methods for threshold voltage, mobility, and subthreshold swing in the linear region were demonstrated in [10]. The experimental results from a partially-depleted (PD)-SOI G 4 FET showed the dependence of these parameters on different gate biases. The charge coupling between front, back, and lateral junction-gates was considered and a 2-D analytical relationship for the fullydepleted body potential was derived in [11]. Here, a closed form front-interface threshold voltage expression was derived as a function of the back and the lateral gate voltages for different back interface conditions such as accumulation, depletion and inversion.
A very interesting application of G 4 FET is the formation of quantum wire. The quantum wire can be electrostatically formed when the conducting channel is surrounded by depletion regions induced by vertical MOS and lateral JFET gates [12]. In this unique conduction mechanism named depletion-all-around (DAA), majority carriers flow in the volume of the silicon film far from the silicon/oxide interfaces. The control of lateral gates on the conduction channel can be adjusted by changing biases on the vertical gates. There is a reduced sensitivity of the channel to the oxide and interface defects, low subthreshold swing, high g m /I D ratio, high mobility, low noise, and high immunity to ionizing radiation [13]. A low frequency noise model combining volume and surface noise sources was presented in [14] and effects of gate oxide and junction nonuniformity on the DC and low-frequency noise performance were investigated in [15]. A charge sheet model has been proposed to analyze the transistor characteristics of fully-depleted G 4 FETs [16]. Here, surface accumulation behavior, drain current and gate capacitance of fully-depleted G 4 FET are modeled analytically. In [17], a mathematical model is developed to determine the subthreshold swing of thin-film fully-depleted G 4 FET. Based on the exact solution of the Poisson equation, a new two-dimensional model of potential and threshold voltage for the fully-depleted G 4 FET was developed in [18]. Another mathematical model is developed in [19] to determine the 3-D potential distribution of a fully-depleted G 4 FET.
Several numerical methods have also been used for modeling G 4 FET. Multidimensional Lagrange and Bernstein polynomial approaches have been used in [20] for modeling and SPICE implementation. The authors in [21] have reported single multivariate regression polynomial approach for modeling the device operation across different operating regimes. We will focus on two new techniques in the next two sections. The first one is based on multivariate cubic spline interpolation method for DC modeling first reported in [22]. The second one is a MOS-JFET macromodel suitable for DC, AC and transient simulation [23].

Multivariate cubic spline interpolation model
Developing a interpolating polynomial f(x) requires fitting a polynomial to a set of chosen data points such as (x 0 , y 0 ), (x 1 , y 1 ) … (x m , y m ). It can be written as, Spline based interpolating polynomial is a piecewise low-degree polynomial which solves some of the problems associated with a single high degree interpolant.
If the dataset is regular and monotonic, this interpolation method can get rid of an important problem of single high degree interpolating polynomial known as Runge's phenomenon, i.e., oscillation between interpolation nodes, which can be catastrophic close to the boundary. With the suitable choice of spline, it is also possible to reduce error while maintain continuity of the function and its first and second order derivatives. For these reasons, multivariate cubic spline polynomials are used for the model development [24].
Given n data points, (x k , y k ), k = 1, 2,…, n with distinct x k ' s, there is a unique polynomial in x of degree less than n whose graph passes through the points. If two successive points are (x k , y k ) and (x k + 1 , y k + 1 ) then the kth interval between these two points can be interpolated using splines. The interval index k is such that, x k ≤ x < x k + 1. The local variable, s is s = xx k. The first divided difference is δ k = (y k + 1 À y k )/(x k + 1 À x k ). Let h k denote the length of kth subinterval, i.e., h k = x k + 1 À x k ; then, δ k = (y k + 1 À y k )/h k . Let d k denote the slope of the interpolant at x k , i.e., d k = P 0 (x k ). The cubic spline polynomial will be the following function on the interval x k ≤ x ≤ x k + 1 , written in terms of local variables s = xx k and h = h k.
The cubic polynomial in Eq. (2) satisfies four conditions for smooth interpolation; two on the function values and two on the derivative values so that, In order to ensure continuous second derivative, we can add an extra constraint which gives rise to the following equation, If knots are equally spaced, Eq. (4) becomes Following the above-mentioned procedure at every interior knot x k (k = 2,. ., n À 1), results in n À 2 equations with n unknown d k 's. At the interval boundary, "not-a-knot" method is used, i.e., a single cubic polynomial is used on first two subintervals, x 1 ≤ x ≤ x 3 , and last two subintervals, x nÀ2 ≤ x ≤ x n . Now, these two boundary conditions give two more equations resulting in n unknowns with n linear equations. Then d k 's are estimated by solving this system of linear algebraic equations.
For cases involving multiple independent variables, tensor product formulation can be used to extend the univariate analysis to multiple dimensions. In each respective dimension, value at a desired point is interpolated according to a cubic interpolation of the values at nearby knot points. In this section, we are using V DS (drain-source voltage), V JG (common junction-gate voltage), V TG (top-gate voltage) and V BG (bottom-gate voltage) as our independent variables for spline interpolation of drain current. Moreover, the geometric features such as W (width) and L (length) may be used as independent variables to include the effect of device geometry.
An n-channel partially-depleted SOI (PDSOI) G 4 FET was developed in a device simulator (TCAD Sentaurus from Synopsys) and a cubic spline model was formulated based on the generated data. In Figure 2, drain current versus drain voltage is plotted for different values of junction-gate voltages ranging from 0 to À4 V. The bottom gate/substrate bias, V BG is 0 V and top-gate bias, V TG is 2.25 V, respectively. Both junction gates are tied together for simplicity. From the graph and the values of corresponding mean relative error, it is clear that the model fits reasonably well for all the isolines. More details on this method can be found in [22].

MOS-JFET macromodel
G 4 FET has also been called MOSJFET [1] since it combines both metal-oxidesemiconductor field-effect transistor (MOSFET) and junction field-effect transistor (JFET) actions in a single silicon island. As mentioned earlier, the vertical gates provide MOS functionality and the lateral gates add JFET functionality. By applying different junction-gate biases, it is possible to modulate the threshold voltages of top-and bottom-gates. As a result, G 4 FET can be thought of as a combination of two MOSFETs primarily responsible for top and bottom surface conduction and a dual-gate JFET that is primarily responsible for volume conduction. The top-gate threshold voltage is V TH and the bottom-gate voltage causing the onset of accumulation and inversion at the bottom-gate are V acc BG and V inv BG , respectively. Some of the terms used in the model are introduced below: Junction-gate, top gate and bottom gate capacitances are C JG ¼ ε Si= w , C ox1 ¼ ε ox = t ox1 and C ox2 ¼ ε ox = t ox2 , respectively. Three constants based on device geometry, α, β and γ are defined as, α ¼ Si . Here, W = width of the transistor, t si = silicon film thickness, t ox1 = top oxide thickness, t ox2 = buried oxide thickness, V T = kT/q is the thermal voltage, N d = donor concentration in the body, n i = intrinsic carrier concentration, ε si = permittivity of silicon, and ε ox = permittivity of silicon dioxide.
The onset voltage of accumulation and inversion for the bottom-gate, V acc BG and V inv BG , can be expressed [27] as, Figure 2.
Comparison between isolines of test data and cubic spline model for different junction-gate voltages ranging from À4 V to 0 V in 1 V increment arranged from bottom to top [22].
The back gate may be accumulated, depleted or inverted. When the bottom-gate is in inversion, i.e., V BG < V inv BG , When the bottom-gate is depleted, i.e., V inv BG < V BG < V acc BG , When the bottom-gate is in accumulation, i.e., V BG >V acc BG , Here, V FB1 and V FB2 are the flat band voltages of the top-gate and the bottom gates, respectively.
Based on the above relationships among different gates, a macromodel subcircuit is created combining conventional SPICE Level-1 MOSFET and the JFET models which follow from the quadratic FET model of Schichman and Hodges [25]. It should be noted that the model can be improved by using higher-level MOSFET models (Level-2/3) and BSIM/BSIMSOI models for devices with deep submicron geometry. Here, the goal of this work, with simple first order MOSFET and JFET models, is to show the feasibility of the macromodel and demonstrates how it can capture the essential physics underlying the complex interaction between multiple gates. Here, the top conduction is modeled using a MOSFET and the volume conduction is modeled using a JFET. However, instead of a constant threshold MOSFET, the subcircuit allows for threshold voltage modification based on multiple gate biases using the relationship described above. The transient simulation results for the first analog multiplier circuit configuration in [3] is shown in Figure 3. The simulation result using the macromodel shows good matching with the experimental result [3].

Single photon avalanche diode (SPAD)
Single photon avalanche diodes (SPADs) based optical detectors have gained interest for use in a wide range of applications such as biochemical analysis, imaging and light ranging applications [26][27][28][29][30][31][32]. The main causes of the increased popularity are the exceptional level of miniaturization and portability, overall high performance and low costs due to the integration of SPADs with mixed-signal readout circuits in standard complementary metal-oxide-semiconductor (CMOS) technology [26][27][28][29][30][31][32]. The capability to detect single photon and provide subnanosecond time resolution along with the low-power and high-speed CMOS readout circuits have made SPADs superior to other optical detectors in high performance weak optical signal detection applications. In order to correctly simulate the behavior of the SPAD device and the performance of its readout circuit during the design phase, a suitable and comprehensive model of the SPAD is required. In this work, we review and discuss the most recent developments in SPAD modeling to predict the real behavior of the CMOS SPAD systems before fabrication reducing the design cycle.

SPAD theory and operation
A single photon avalanche diode (SPAD) is a p-n junction, which is biased above the breakdown voltage. The device operates in the so-called Geiger mode. With the increase of reverse bias voltage beyond a critical potential, the current increases rapidly. The rapid increase in current is due to the free carriers, which have gained enough energy to ionize the fixed lattice atoms and free other carriers through high-energy collisions. This rapid multiplication process results in a sudden large avalanche current. The reverse voltage above which this multiplication process occurs is called breakdown voltage [26][27][28][29][30][31][32]. In order to respond to an incident photon, impact ionization is necessary for a SPAD. A high electric field is required to break the covalent bonds and generate the free carriers, which take part in the avalanche process. The probability that carriers gain the threshold energy is a function of the local electric field and previous states of the carrier. This probability determines the ionization rates. Figure 4 illustrates impact ionization causing an avalanche current in the presence of a high electric field. A number of free carriers is generated at each step, quickly multiplying and generating a large current [10].
Free carriers are also generated in absence of photon due to inherent noise processes such as thermal generation, minority carrier diffusion, and band-to-band tunneling. Dark count rate is the effective noise for SPADs caused by these non-photon driven carrier generations. The large avalanche current must be quenched to protect the long-term use of the device. An external circuit, typically a passive ballast resistor or a transistor, is employed to quench the large avalanche current. As the current increases the voltage drop across the quenching resistor increases. As a result, the bias voltage across the SPAD reduces below the breakdown voltage ending the avalanche process. The SPAD is then recharged and ready to detect photon event induced avalanche. The avalanche gain is large enough to be almost infinite in this Geiger mode of operation. The generation of an electron-hole pairs causes the device to be in an "ON" or "OFF" state. Since the generation of an avalanche gives a current or voltage spike, the device exhibits an inherently digital operation [26][27][28][29][30][31][32]. Figure 5 shows the Geiger mode operation of SPAD and quenching of the avalanche current. Here, V br represents the breakdown voltage and V s is the applied voltage. During the current flow, some carriers may be trapped in the deep-levels and are released after a characteristics lifetime later. The released carrier may trigger a new avalanche process when the SPAD is ready to detect another photon. The new avalanche process associated with trapped carriers during the first avalanche pulse causes after-pulses and it reduces detector performance [26].
SPADs fabricated in CMOS provides the benefits of commercial CMOS process such as low cost, miniaturization, and improved performance [27][28][29][30][31][32]. However, due to the planar nature of the junction in standard CMOS process, electric field distributions exhibit maxima at the diodes' edge [27][28][29][30][31][32]. Since SPADs operate in a high reverse bias region, they are more vulnerable to breakdown at the edges. Therefore, full volumetric breakdown cannot be achieved as the diode periphery undergoes breakdown earlier compared to the other parts of the diode. As a result, the active area is decreased and the photon detection efficiency is significantly compromised. This premature breakdown is one of the major problems of SPADs implemented in standard CMOS process [27][28][29][30][31][32]. It has been previously proved that the perimeter gated SPAD (PGSPAD), a p-n junction incorporating a poly-silicon gate surrounding the periphery, suppresses the premature breakdown effectively [29,31,[33][34][35][36][37][38][39]. The applied voltage at the gate terminal modulates the electric field making it uniform throughout the junction preventing the premature breakdown.

State-of-the-art SPAD models
Various SPAD models have been presented in literature to predict the actual behavior and choose suitable parameters during the design phase to ensure the optimal performance for CMOS SPAD systems in existing and new applications [29,[35][36][37][38][39][40][41][42][43][44][45][46][47]. State of the art SPAD modeling approaches can be grouped in three sub groups named as model based on device physics, circuit simulation model, and model based on information theory. The main features of the models of each group along with the limitations are discussed below.

Models based on device physics
The SPAD models in this group have been developed to simulate the SPADs at the semiconductor level based on purely quantum mechanical or semi-classical representation of the transport properties of semiconductor structure underlying the SPAD device [29,[35][36][37]. The key operational features such as breakdown voltage, avalanche gain, and noise have been generated through the model.
In [35], a new compact numerical model to simulate the forward and reverse dc diode characteristics has been presented. The model is based on the solution of 1D steady-state hole continuity equation in the depletion layer of a p-n junction. The developed dc model includes the physical mechanisms such as band-to-band tunneling, trap-assisted tunneling, Shockley-Read-Hall recombination, and avalanche breakdown.
The model presented in [36], focuses on breakdown state stability studying the SPAD transient behavior based on physical and numerical approach. At first, the spatial steady-state free charge carrier distribution during breakdown is computed. Then the temporal change of the distribution is considered. It explores the dependence of the probability for the avalanche process on the quenching resistance. It aids designers to accurately choose the optimal passive quenching resistance for minimizing the overall dead time value.
The reported model in [37] investigates the lateral spreading of the avalanche in Geiger mode reach-through avalanche photo diode, which is the main cause of avalanche dynamics and the time jitter. The model estimates the photon assisted effect based on the Poisson equation and the time-dependent semiconductor continuity equation to ensure the avalanche dynamics. A Monte Carlo simulation of the avalanche current rise has been developed in the model.
The effect of structural effect on the breakdown voltage and the localization of the avalanche region within p+/n-well junction was determined in model [29]. The reported 2-D numerical model is based on Poisson, hole and electron current continuity equation coupled with carrier generation rate equation due to impact ionization. The model highlights the regions most susceptible to breakdown through the determination of the spatial maxima of carrier generation rates and establishes that the generation rate reaches its peak at the perimeter and the surface of the junction causing perimeter breakdown. Finally, it demonstrates that full volumetric breakdown can be achieved using perimeter breakdown suppression techniques.

Circuit simulation model
Circuit simulation models aim to represent the equivalent circuit representation of the SPAD for circuit level designs and simulation. These models provide a first order estimation of performance at the circuit level and allow simulation of readout circuits based on the electrical characteristics of SPAD [31,33,[38][39][40][41]. The readout delays, time constant for both quenching and reset can be calculated using the circuit simulation model. Traditional SPAD basic model reported in [38] includes voltage sources, resistors and capacitors to represent the diode characteristics ( Figure 6). In the model, the diode resistance, R D , includes both space charge resistance and resistance of neutral regions crossed by the avalanche current. C AC represents the junction capacitance and C AS and C CS represents the stray capacitances from anode and cathode respectively. The closure of NMOS switch simulates the triggering of an avalanche due to photon absorption or other phenomenon including thermal generation or bandband-band tunneling. Figure 7 shows the improved version, which includes the triggering, the selfsustaining process, and the self-quenching of the avalanche by incorporating of current-voltage controlled switches [39]. In order to represent the nonlinear I-V characteristics above breakdown, the model includes nonlinear voltage generator, V SPAD , to generate piecewise linear curve with different slopes.
An accurate model for SPADs [40] exploits the behavioral description of the sensor with biasing circuits and aims to provide optimal biasing circuit design through SPAD simulation. In order to represent the SPAD above breakdown and with the avalanche triggering and self-quenching mechanisms, a complete SPICE circuit model for SPAD has been presented in [41]. The SPAD is implemented using a piecewise non-linear approximation and modeled through a non-linear series resistance above breakdown. The model also includes the forward region and the secondary breakdown due to edge-junction or punch-through effects (Figure 8).
SPAD stochastic phenomena affect the switching behavior of the device and define the transition from no-avalanche to avalanche (turn-on) and vice versa. In order to provide a more realistic simulation platform, more comprehensive models [34,[42][43][44][45] include stochastic nature of Geiger mode operation based on theoretical equations of dark count rate, breakdown probability, signal-to-noise ratio and afterpulsing probability.
The behavioral model presented in [42], simulates the noise model due to dark count rate and after pulsing based on fabricated component for the first time. The model improves the previous circuit model [41] by adding the noise sources of avalanche photodiode. The model has been developed using specter active and passive components to represent the unwanted after-pulse and dark count events. At first, the parasitic capacitors of the sensors (anode-cathode capacitor, anodebulk and the cathode-bulk capacitors) are included to represent the dynamic behavior. Three different branches representing forward, reverse and effects of noise behavior are incorporated to simulate the I-V behavior. Each branch is controlled by an ideal voltage controlled switch. Finally, two noise sources, represented with two corresponding branches are added to simulate the noise effect ( Figure 9). An accurate behavioral model reported in [43] models the major statistical behaviors of SPADs including the turn-off probability, dark count rate and afterpulse phenomena. One the major focus of the model was to explore the dependence of the device capacitance on the reverse bias. Complete SPAD model [41] with switches S TRIG and S SELF representing avalanche triggering and selfsustaining; two additional branches are included to represent forward and reverse biasing. Improved SPAD model [39]; switch S TRIG is used to trigger the detector through "photon" input, S SELF is used to include self-sustaining and self-quenching of the avalanche, and V SPAD represents the nonlinear I-V curve.
SPAD is unresponsive to any subsequent photons after detecting the first photon (during avalanche, quenching and resetting cycle) until it is reset back to undergo another avalanche. This time duration is known as dead time. In order to accurately interpret the SPAD measurements, dead time needs to be taken into account especially for large arrays of CMOS SPAD based detectors. In [44], a dead time model was reported for externally reset SPADs. Here, time distribution of incidents counts is assumed to follow Poisson distribution.

Model based on information theory
In order to simulate the behavior of SPAD's array and predict the behavior at system level design, model based on information theory include a combination of the previously mentioned models [47,48]. These system levels models are helpful to simulate the operation of multiple SPADs including both high-density and lowdensity arrays. Various operational aspects of SPAD arrays such as array breakdown voltage non-uniformity, pixel cross-talk can be investigated using these models. Some models are able to study the performance at single-pixel-level [48]. They consider the SPAD front-end as an information channel. Then the information capacity and mutual information can be calculated from the information channel as a function of excess bias voltage, dark count rate, and after pulsing probability. The information theoretic model reported in [47] includes a channel-capacity metric to evaluate the performance of SPADs as a detector in laser communication systems for optimizing its performance with respect to the device structure and operating voltage.

State-of-the-art PGSPAD models
A perimeter gated SPAD (PGSPAD) is a SPAD with an additional polysilicon gate terminal and its effects need to be properly modeled. PGSPAD models were reported in the literatures to predict the Geiger mode response to a photon; simulate the static and dynamic I-V characteristics considering the effect of applied voltage at additional gate terminal [31,33,34,48]. The equivalent circuit models for PGSPAD [31,33] simulate the static and dynamic behavior of PGSPADs with the additional gate terminal. These models are based on the improved SPAD equivalent model [41]. The PGSPAD model simulates the avalanche in similar way to SPAD model [41]. In addition, they include a dependent voltage source to represent the dependence of breakdown voltage (V br ) on gate voltage (V G ) and the reverse biased resistance as a variable resistance depending on gate bias voltage Figure 9. SPAD behavioral model including I-V behavior and noise sources for dark counts and after-pulsing [42]. and reverse bias voltage ( Figure 10). Two parallel branches are incorporated to simulate the forward and reverse bias operation. Model parameters and threshold voltage of the switches are extracted from the fabricated PGSPAD device. From the empirical data, it was found that the breakdown voltage depends on zero-bias breakdown voltage and the gate voltage. The empirical relation [31] of V br, reverse branch resistance (R R ), forward branch resistance (R F ) with gate voltage (V G ) and anode-cathode biasing (V CA ) are, where γ is a fitting constant obtained from experimental perimeter gate voltage versus breakdown voltage, P, Q , R, and S are constants determined experimentally, V br0 is the zero bias breakdown voltage.
An improved PGSPAD SPICE model was reported in [34] to simulate the stochastic behavior of PGSPAD. The model correctly simulates the dark count rate with variation of gate voltage and excess bias voltage, the applied reverse biased voltage in excess of breakdown voltage. The model includes the dark current noise due to non-photon driven carrier generation mechanisms and the spectral response with the additional gate terminal. It uses Shockley-Hall-Read equation and tunneling current equation to calculate the thermal carrier generation rate and carrier generation rate due to band-to-band tunneling respectively. Finally, the model established that the band-to-band tunneling mechanism depends on the gate bias voltage through the electric field. The PGSPAD model presented in [48] evaluates the performance of PGSPAD using the information theory principles. The model estimates the optimization of the information rate of the PGSPAD by applying the asymmetric communication channel model to the device.

IMT
Insulator Metal Transition (IMT) devices are volatile memristors that have recently spurred significant interest in the research community [49]. In crossbar memory arrays, IMT devices are often leveraged as selector devices owing to their Back-End-Of-Line (BEOL) compatibility and high ON/OFF ratio [50,51]. On the neuromorhpic front, the inherent switching dynamics of IMT devices have been shown to mimic the Integrate-And-Fire neurons which alleviates the need for complex CMOS circuitry, thus, enabling higher integration density [52,53]. In addition, volatile memristors with short-term memory [54][55][56] can also be used as synapses for neuromorphic application.

Device operation
Several studies have shown that temperature is the main driver of phase transition [57,58] while other works argue that the electric field is the main driver of resistive switching while the Joule heating plays a secondary role [59]. A more thorough study about IMT's switching mechanism is reported in [60,61] which show that Joule heating is not sufficient to instigate switching and that an Electric field assisted phase transition is more likely. The work in [60] proposes that a certain threshold voltage is required to effect a phase transition which decreases as temperature increases.

Motivation behind compact model
In [52], an IMT device model that captures the positive feedback between the temperature and electric field was presented. Here, the relationship between the device resistance and the device temperature was implemented using a look up table. In [62], another device model was proposed based on band theory. The IMT device is modeled as a low bandgap semiconductor. Increasing the device temperature results in decreasing the bandgap, thus, increasing carrier concentration. This increase in carrier concentration ultimately results in decreasing the device's resistance. A model is also presented which captures the change in the thermal conductivity with temperature. The thermal conductivity model along with the band gap model are then solved in a self-consistent fashion to effect a temperature dependent phase transition. The model in [62] was implemented in Sentaurus TCAD simulator wherein the built-in electrothermal models and finite element drift-diffusion model where leveraged. This model, similar to the previous one, is best used in a TCAD simulation flow and not SPICE level simulators. The lack of a physics-based compact model, however, has hindered circuit designers from exploring the full potential of IMTs in circuit applications which is the motivation behind the development of the following simplified compact model.

Simplified compact model
This work focuses on an IMT compact model originally proposed in [63]. The model describes the IMT device as a memristor with the device's local temperature being the state variable. The proposed model is simulated using Specter from Cadence and exhibits a close match to experimental data and electro-thermal simulations based on the models in [52,62].
As the current flows through the device, the device's temperature increases until it reaches a critical temperature at which point it transitions from a high resistance state to a low resistance state. As the temperature decreases, the resistance relaxes back to its default high resistance state. The memristive dynamics of the IMT can be expressed as follows [64].
where (13) and (14) represent the port and state equations, respectively, and x is the state variable. The proposed model consists of two governing equations: (a) the resistance change equation that corresponds to the port equation and (b) the temperature evolution equation which corresponds to the state equation, with the temperature being the state variable such that x = T(t).
The behavior of the resistance change versus temperature is described by two thermistor states for the high resistance state and the low resistance state and a sigmoid function capturing the transition from a high resistance state to a low resistance state. Since the thermistance behavior depicts a linear relationship between the resistance and the temperature in the Log-Linear domain, one can readily model the two thermistors as exponential functions of the temperature such that R LRS ¼ R LRS F e ÀB LRS T t ð ÞÀT F ð Þ and R HRS ¼ R HRS 0 e ÀB HRS T t ð ÞÀT 0 ð Þ . R LRS F is the low resistance state defined at temperature T F (a reference temperaure) and R HRS 0 is the high resistance state defined at the ambient temperature T 0 . B LRS and B HRS are the negative temperature coefficients of the thermistors which can be readily extracted from the slope of the thermistance versus temperature plot. This modeling framework, however, requires clipping of the R LRS and R HRS at some minimum and maximum values to avoid any unphysical behavior in circuit simulation. Clipping, however, often uses conditional statements, which hampers the "smoothness" of the model yielding potential convergence difficulties in circuit simulation. Hence, the model equations are reformulated such that R LRS and R HRS smoothly approach to R LRS F and R HRS 0 at high and low temperatures, respectively.
This relationship between the resistance and the temperature can be described as follows: where The fitting parameter, T x , captures the sharpness of the resistance transition. T c is the critical temperature which is around 340 K in the case of VO 2 devices [53]. R LRS and R HRS are the low resistance state and high resistance state, respectively. A is a control parameter that governs how the two thermistors approach the asymptotes [65]. While the model might seem complicated at first glance, the principal equations are simple exponential functions. This formulation is employed to abide by compact modeling practices as suggested in [65,66].
The temperature evolution dynamics are described by the compact thermal model in [67] as shown in (4): where V IMT I IMT is the Joule heating, C th and R th are the effective thermal capacitance and the effective thermal resistance, respectively, and T 0 is the ambient temperature. This model assumes that the device stays at an effective temperature T(t) and exchanges heat with the ambient environment at an ambient temperature T 0 .
The model predictions closely follow the experimental data as well as electrothermal simulation [52] as shown in Figures 11-15. Figures 11 and 12 capture the resistance transition around the critical temperature which is about 340 K in VO 2 devices fitted against experimental data from [52]. Figure 13 depicts the hysteresis in the I-V domain (a characteristic of memristors) exhibited by the IMT device as shown in [52,53] and fitted against the experimental data from [52]. Figures 14 and  15 capture the time dependence of temperature and resistance evolution, respectively, fitted against electrothermal simulations from [52]. Three voltage levels, based on the values used in [52] were applied across the IMT device: 1.4 V (blue), 1.6 V (red) and 1.8 V (green). One can readily observe in Figure 14 that the local temperature of the device saturates at a higher temperature value for higher     voltages due to increased Joule heating. In Figure 15, higher voltages result in a faster transition due to higher rate of joule heating. Figure 16 depicts a simple neuron based on the proposed circuit in [52]. First, the input voltage spikes are fed through the synaptic network. The dot product multiplication is then executed between the input spikes and the synaptic weights via Ohm's law such that the accumulated sum follows I sum ¼ ∑ i V i G i . I sum is then fed to the neuron which is represented by the parallel combination of the IMT device and the capacitor. The core of the neuron is the IMT device which switches from R HRS to R LRS when the device's temperature exceeds the critical temperature resulting in a current spike. The current spike is converted to a voltage pulse via a CMOS inverter. A spike generator is then added to generate an output voltage spike with a pulse width controlled by the RC time constant of the RC network preceding the output buffer. Figure 17 depicts the simulation results of the neuron circuit.  The essence of neuron oscillation rests in the IMT device alternating between R HRS and R LRS . At R HRS , the steady state temperature exceeds the critical temperature and, accordingly, the neuron fires. At R LRS , the steady state temperature drops below the critical temperature, the neuron resets and the process is repeated for the next inputs.

Application in neuron design
At steady state ( dT dt ¼ 0Þ, the solution to the differential equation in (18) can be described as follows: where T ss is the steady state temperature of the IMT device. Hence, according to the aforementioned explanation, the oscillation condition can be expressed as follows: R th I 2 IMT R LRS < T c À T 0 < R th I 2 IMT R HRS Inequality (7) establishes the oscillation condition for the IMT-based neuron as a function of device parameters such as R th , R LRS , R HRS and T c and circuit variables such as I IMT which is a function of the amplitude of the voltage spike and the series resistance with the IMT device including the synapse resistance, diode ON resistance and the IMT series resistance. The neuron typically operates in three phases: (I) accumulation wherein the inputs from the synapse networks are summed, (II) firing when the accumulated value reaches the neuron's threshold and (III) refractory period wherein the neuron is idle.
In CMOS neurons, two operational amplifiers are required for accumulation (integration) and firing (comparison). In addition, a feedback circuit is often employed to implement the refractory period. These operational amplifiers, besides entailing all the complexities of analog design, are also area consuming and power hungry. On the other hand, the IMT device can provide the accumulation function through the heating of the device, fire through device transition and a refractory period during device cooling should the device be placed in such configuration. A typical CMOS neuron, such as the work in [68], requires more than 20 transistors while the proposed IMT-based neuron only requires 7 transistors.

Conclusion
Computer aided design and simulation plays a major role in the advancement of semiconductor industry. A multitude of exciting new devices has emerged in recent years. There is a growing need in the industry and academia for fast and reliable compact models of these emerging devices to enable useful circuit design leveraging their unique capabilities. In this chapter, we show how different approaches can be adopted to model three emerging semiconductor devices namely, silicon-on-insulator four gate transistor (G 4 FET), single photon avalanche diode (SPAD) and insulator-metal transistor (IMT) device with volatile memristance. All the models have been verified against experimental/TCAD data and implemented in commercial circuit simulator. The ideas developed in this chapter can also be transferred to model other emerging devices.