Thin-Film Transistor Simulations with the Voltage-In-Current Latency Insertion Method

This article presents formulations for the voltage-in-current (VinC) latency insertion method (LIM) for thin-film transistors (TFTs). LIM is a fast circuit simulation algorithm that solves circuits in a leapfrog manner, without requiring intensive matrix operations present in SPICE-based simulators. This allows LIM to have a far superior scaling with respect to the size of the circuit resulting in significant time savings on large circuit networks. The VinC LIM formulation for the TFTs written in this article has the benefit of a better stability compared to the original LIM formulation which allows the use of larger time steps. The performance of the new algorithm is demonstrated through the simulation of numerical examples of large flat-panel display (FPD) circuits. It is seen that VinC LIM greatly outperforms basic LIM and commercial SPICE-based simulators, where the presented algorithm is able to simulate circuits with more than 10 million nodes or devices in a reasonable time, which is not viable in many modern day SPICE-based simulators.


I. INTRODUCTION
A S THE DRIVE for higher functionality, capacity, and efficiency continues to push the consumer industry in IC products forward, the complexity and sheer size of the analog and digital circuitry inside these devices continue to rise in tandem as well. In order for a design to succeed, a mature product with the latest process technology node requires support during every step of the development process. From front-end design steps, to back-end packaging and test procedures, all sections of the design and verification process have to be equipped with the latest technology with sufficient capabilities to avoid potential failures. This presents a big challenge in the semiconductor industry today as product designs need to be signed-off, by passing through various modeling and simulation verifications, before being sent to the foundry to reduce the time and manufacturing cycles with the ultimate goal of minimizing development costs.
In this regard, access to accurate and fast modeling and simulation tools is paramount to the success of a product.
Most of the commercial tools on circuit simulations today rely on the Simulation Program with Integrated Circuit Emphasis (SPICE) [1] algorithm, which was originally developed in 1973, based on the modified nodal analysis (MNA) framework, with a heavy reliance of sparse matrix solvers. However, with the continued scaling of process technology nodes, more and more devices are now crammed into a single product, where a full circuit design can contain tens or even hundreds of millions of device elements. This creates a simulation bottleneck, where state-of-the-art simulators are no longer able to cope with the sheer large size of the resulting MNA matrices from these circuits, thus opening the door to new methods which could potentially replace the presently aging simulation algorithms.
Recently, thin-film transistors (TFTs) have grown in prominence as the standard device in display technology particularly in flat panel displays (FPDs) based on the organic light-emitting diode (OLED) display technology [21]- [23]. In these active-matrix devices, individual active switching elements are used to control each of the pixel in the panel. As a result, the simulation of these circuits are incredibly resource consuming from a runtime perspective as they contain millions of transistors, which make up the array of pixels for a given resolution. For example, a full HD (1080p) display can contain more than 40 million TFTs, and this number goes up even higher in 4K and 8K displays. These massive circuits are far beyond the capacity of state-of-the-art commercial simulators, which often struggle in both runtime and memory requirements due to the poor scaling of their matrix based operations.
In this article, formulations for the simulations of TFTs in the LIM framework are presented. First, the background and fundamentals of the original LIM and VinC LIM are reviewed. Then, augmentations to the formulations are presented for TFTs, which include the drain-source current, and the gate-drain and gate-source intrinsic capacitance currents. Practical steps are discussed to control the accuracy and improve the runtime. Finally, numerical results are presented, including large FPD arrays of varying resolutions, to illustrate the viability of the method for the simulations of very large circuits. Results show that the proposed TFT VinC LIM simulation algorithm is able to produce accurate simulation waveforms, while having a capacity far exceeding the original LIM simulation algorithm and state-of-the-art SPICEbased commercial simulators.
The rest of the paper is organized as follows. In Section II, the fundamentals of the original LIM and VinC LIM algorithms are reviewed. Then, in Section III, new formulations of the LIM algorithms for TFTs are derived in both the original LIM and VinC LIM frameworks. In Section IV, numerical results are presented for different TFT circuits, and the performance of the developed VinC LIM TFT algorithm is compared with the original LIM algorithm and a current generation SPICE simulator in terms of accuracy and speed. Finally, Section V concludes the article and proffers some future work on the subject.

II. LIM FORMULATION
The fundamentals of basic LIM and VinC LIM are reviewed in this section.

A. BASIC LIM FORMULATION
The original LIM or basic LIM is a circuit simulation algorithm that is devised for the solution of large and dense networks. To analyze a circuit with LIM, the circuit must be composed of nodes and branches with topologies as shown in Fig. 1, where V i and V j are node voltages at node i and node j respectively, I ij is the branch current flowing from node i to j, L ij , R ij , and E ij are the branch inductance, resistance, and voltage source respectively, C i , G i , and H i are the node capacitance, conductance, and current source respectively, and the currents labeled from I i1 to I ik are the branch currents that are connected to node i. From the circuit topology in Fig. 1, the LIM algorithm is formed by applying Kirchhoff's laws, where the inductance and capacitance in each topology results in a derivative term in their equations, which are linearized through the application of Euler's method. By solving the equations, the basic LIM formulations can be obtained as where M i indicates the total number of branches connected to node i, ∆t is the time step of the simulation, and n is the time index which marks the time point in the transient analysis. In a LIM simulation, the leapfrog concept is applied where (1) and (2)  Compared to SPICE based methods, which make use of MNA stamps and implicit simultaneous solutions of equations, the LIM simulation is performed in an explicit manner. This gets rid of heavy matrix computations and hence, allows the LIM analysis to scale better with the size of the circuit and be carried at a relatively faster speed, especially for very large networks. However, since Euler's approximation is used, the selection of the time step becomes a limiting factor for the stability of the simulation. Through the application of Lyapunov's method [24], the maximum allowable time step value for stability can be calculated by where n N is the number of nodes in the circuit, n i B is the number of branches connected to node i, C i is the value of the shunt capacitor at node i, and L im is the value of the inductor in the m th branch connected to node i. Based on (3), it can be observed that the maximum stable time step value of a LIM simulation depends on the minimum value of the latencies in the circuit. Therefore, there is always a trade-off between simulation time and accuracy in a LIM simulation if fictitious components are required to be inserted into the circuit. Bigger fictitious components lead to a drop in accuracy while smaller values lead to a significant increase in total simulation time as the maximum stable time step is reduced to very small values.

B. VOLTAGE-IN-CURRENT LIM FORMULATION
The main limitation of LIM is its conditional stability which results in the use of very small time steps. To solve this, various improved versions of LIM have been introduced. VinC LIM [8] is an improved formulation of LIM for better stability and it is formulated by rewriting equations (1) and (2) in their respective implicit versions at the same time step of n + 1. These equations are given by Then, VinC LIM is formed by substituting the voltage equation in (4) into (5) and solving the substituted equation for the I n+1 ij terms. The resulting current equation in VinC LIM is given by is the sum of currents flowing out of node i, and Mj k=1,k =i I n+p jk is the sum of currents flowing out of node j, without considering the I ij current in both summations. In order to circumvent the added complexity of the implicit formulation, a forward branch-marching scheme is applied in VinC LIM, where the sequence of the currents in the updating process is taken into consideration, and the latest available current value is always used. This is indicated by the indices p that are shown in the equation, where p is either 1 or 0, depending on whether the associated current has been solved for or not, before the present computation at this time step. Compared to the explicit derivation of the basic LIM formulation, the implicit derivation of the VinC LIM formulation leads to improvements in the stability of the algorithm. Previous results on VinC LIM have shown the method to be unconditionally stable, thus allowing the use of timesteps far greater than that which is permitted in basic LIM.

III. LIM FORMULATION FOR RPI THIN-FILM TRANSISTORS
Both the basic LIM and VinC LIM formulations discussed in the previous section are limited to the simulation of linear circuits only. Hence, more work is required to extend the coverage of the LIM formulations to other devices especially transistors in order to support the simulations of various products. In this work, the formulations for basic LIM and VinC LIM for thin-film transistors (TFTs) are presented. TFTs are the main active devices which are commonly used in display panels such as those based on the organic lightemitting diode (OLED) display technology. This section covers the formulation for the drain-source current and the gatedrain and gate-source intrinsic capacitance currents based on the Rensselaer Polytechnic Institute (RPI) TFT model using LIM's formulation.

A. BASIC LIM FORMULATION FOR RPI TFT DRAIN-SOURCE CURRENT
For circuit simulations, the operation of a TFT device is mainly classified into two parts, which are the calculation of the drain-source current, I ds , also known as the main current, and the intrinsic capacitance currents made up by the gatedrain current, I gd and the gate-source current, I gs . Fig. 2 shows a simplified equivalent circuit of a RPI TFT model where V d , V g , and V s represent the node voltages at the drain, gate, and source nodes of the TFT respectively, the three arrows represent the default directions of I ds , I gd , and I gs current flows, and C gd and C gs are the intrinsic capacitances that exist between the gate-drain and gate-source terminals VOLUME  respectively. In this part, the formulation of the I ds current in the basic LIM formulation is discussed.
In basic LIM, the voltage equations at the three different terminals of the RPI TFT model are given by where V d , V s , and V g represent the voltages at the drain, source, and gate nodes of the TFT at their respective indices of n + 1 2 and n − 1 2 . C is the shunt capacitance, G is the conductance, and H is the current source at the nodes labeled according to their respective nodes d, s, and g, and M d , M s , and M g represent the number of branches connected to nodes d, s, and g respectively.
In the basic LIM formulation, the drain-source current formula of the TFT model can be applied directly in the equations. Taking an example using version 1.0 of the RPI poly-Si TFT model, the drain-source current equation is given by where I sub is the subthreshold current, I a is the above threshold current, ∆ kink takes into account the kink effect, and I leak is the leakage current. By adding suitable indices to the current terms in the equation, the basic LIM current equation for the version 1.0 RPI poly-Si TFT drain-source current is given by Although not indicated here, the terms in (11), namely I n+1 sub , I n+1 a , ∆ n+1 kink , and I n+1 leak depend on the voltage terms, V . A detailed explanation of these terms will be shown in the next part of this section together with the explanation of the VinC LIM formulation for the drain-source current. For illustration, Fig. 3 shows the chain relationship between the terms in the version 1.0 RPI poly-Si TFT model from the terminal voltages, to the final drainsource current formula.
This method can also be extended to other RPI TFT models. For example, in the version 2.0 RPI poly-Si TFT model, the basic LIM current equation can be written as where I n+1 ds1 is the channel current, and ∆ n+1 kink and I n+1 leak are again the variables accounting for the kink effect and the leakage current respectively, and both of them are functions of V . Compared to version 1.0, this version includes more effects and is more accurate for shortchannel RPI TFT devices. The detailed formula for these terms will be shown in the next part of this section. Similar to Fig. 3, Fig. 4 shows the relationship between the voltage dependent variables in this version where more variables and more complex relationships can be observed.

B. VINC LIM FORMULATION FOR RPI TFT DRAIN-SOURCE CURRENT
To derive the VinC LIM formulation, the voltage equations at the drain, source and gate terminals of the TFT can first be written as where the same time index of n + 1 is used in the voltages and the currents, and the I ds terms have been extracted out from the total current summations for emphasis. The signs of the I ds currents are different in the two equations to signify a current direction flowing in or out of a node. The I ds current is not shown in the gate voltage equation as in most cases, the gate node is independent of the I ds current, unless the gate and drain, or gate and source nodes of a TFT is shorted,  cases, equation (13) or (14) will be used to represent the gate voltage equation depending on the respective situation. Then, a general equation of I ds can be written as for a RPI TFT model where its drain-source current depends on the voltages at all the terminals of the TFT. Since the terminal voltages of the TFT can also be functions of I ds , this equation can be further written as where the time indices have been omitted to avoid convoluting the expressions.
To apply the concept of VinC LIM into the I ds current formulation, (17) needs to be solved for the I ds term. This can normally be done by moving all the I ds terms in the right-hand side of the equation to the left-hand side of the equation, and solving for I ds . However, since the TFT is a non-linear device, this is not a straightforward process, and linearizing it would cause accuracy drops especially when larger time step values are applied. To overcome this, a Newton-Raphson iteration is introduced to the process of solving the drain-source current. Such an idea has also been applied in the formulation of VinC LIM for diodes [25]. The general formula of the Newton-Raphson method is given by and by substituting the x terms into I ds , the equation can be written as where I ds,new and I ds,old are the present and previous I ds values in the Newton-Raphson loop, F (I ds,old ) is the value of the differentiable function F (I ds ) when substituted with the I ds,old value, and F (I ds,old ) is the derivative of F (I ds,old ) with respect to I ds,old . The differentiable equation F (I ds ) and its derivative F (I ds ) can be obtained from (17), where they are described by (20) and are the resulting derivatives with respect to I ds of the drain voltage, V d , and the source voltage, V s , from the LIM node equations in (13) and (14) respectively. Note that V d (I ds ) and V s (I ds ) can be zero if the node is independent of the I ds current. The exact equations use in F (I ds ) and F (I ds ) are dependent on the types and properties of the TFT models. In this formulation, the version 1.0 and version 2.0 RPI poly-Si TFT models are considered. For the version 1.0 RPI poly-Si TFT with the I ds equation shown in (10), the differentiable equation for this TFT model is given by (22) and the derivative of the equation is given by where I a , I sub , ∆ kink , and I leak are the derivatives of I a , I sub , ∆ kink , and I leak respectively with respect to I ds . Then from the RPI TFT modeling equations [26], the subthreshold current, I sub is given by where u s is the subthreshold mobility, f cox is the ratio of the gate insulator permittivity to the thin-oxide thickness, w ef f is the effective width, l ef f is the effective length, and v sth = η · v th , where η is the subthreshold ideality factor and v th is the thermal voltage at the device temperature. These four parameters are model parameters that are independent of I ds .
Conversely, the terms V gt and V ds are variables that depend on I ds . V gt can be expressed as with where a t is the first drain induced barrier lowering (DIBL) parameter, b t is the second DIBL parameter, v si is the first parameter for V gs dependence, v st is the second parameter for V gs dependence, v to is the threshold voltage, and The formula for V ds is simply Then, the derivatives of (24) to (28) with respect to I ds are given by are temporary variables used to simplify the equation, by assuming that the gate terminal is independent of I ds , and (33) For the above threshold current, I a , the equation is given by where a sat is the proportionality constant of V dsat , and V dsat , µ f et , and V gte are I ds dependent variables with equations where u 0 and u 1 are the high and low field mobility parameters respectively, m u is the low field mobility exponent applied in (36), and where δ is the transition width parameter. The differentiation of (34) to (37) then yields The kink effect is observed during large drain biasing in the TFT model and the equation is written as , v kink is the kink effect voltage, l kink is the kink effect constant, m k is the kink effect exponent, and V dsep is the inner variable which depends on I ds . It is given by where V dsat is given in (35). The differentiation of (42) and (43) with respect to I ds through the chain rule, yields and where are temporary variables used to simplify the equation. Finally, for the leakage current, I leak , its effect to the overall drain-source current is minimal and its equation is described by where i 0 is the leakage scaling constant, b lk is the leakage barrier lowering constant, x te is a temperature-dependent variable, and X tf e and I diode are functions of I ds , where X tf e also depends on X tf e,lo , X tf e,hi , P f , and F f . Their detailed formulations can be referred from [26] and the relationship between the variables can be referred from Fig.  3. The derivative of (46) is given by VOLUME 4, 2016 This completes the formulation for the drain current of the TFT in the version 1.0 RPI poly-Si TFT model. If intrinsic capacitances are not considered in the TFT model, then at this point, a complete VinC LIM simulation for the TFT circuit can be performed by using (19), (22) and (23) to solve for the I ds currents and (13), (14) and (15) for the node voltages. To ease the convergence of the Newton-Raphson iterations, the I ds value at the previous time step can be used as the initial guess in I ds,old , and the converged value will be taken as the solution.
Next, the formulation for the drain-source current of version 2.0 of the RPI poly-Si TFT model will be discussed using the VinC LIM approach. Its drain-source current equation is given in (12) where the subthreshold current and above threshold current terms have been replaced with the channel current term. The differentiable equation and its derivative for this TFT model are described by and F (I ds ) = 1 − (I ds1 · (1 + ∆ kink ) + I ds1 · ∆ kink ) − I leak (49) where I ds1 , ∆ kink and I leak are derivatives with respect to I ds of I ds1 , ∆ kink and I leak . The equation for the channel current, I ds1 , is given by where λ is the channel length modulation parameter, m e is the long channel saturation transition parameter, and G ch , V dsat , and V ds are functions of I ds . The derivative of (50) with respect to I ds can be written as are temporarily variables used to simplify the equation. V ds and its derivative are given in (28) and (33) respectively, V dsat and its derivative are given by and where G ch and G ch are described as where r s and r d are the effective access resistances. Proceeding further, I sat is given by G chi · r s + 1 + α + 1 + 2 · G chi · r s + (1 + α) 2 (56) with G chi , V gte , and α are variables changing according to the terminal voltages and hence, depend on I ds . Differentiating (56) with respect to I ds through the chain rule yields where L = G chi · V gte , L = G chi · V gte + G chi · V gte , P = G chi · r s + 1 + α + 1 + 2 · G chi · r s + (1 + α) 2 , and P = G chi · r s + α + G chi ·rs+α ·(1+α) √ 1+2·G chi ·rs+(1+α) 2 . G chi is given by where N s and µ ef f are voltages dependent variables and q is the electron charge. The derivative of G chi with respect to I ds gives For N s and N s , the equations are given by where S = Vgt η f ·v th and S = where V l = v max · l ef f µ ef f and v max is the saturation velocity, and the derivative of α is given by . For µ ef f , its expression is described as where u s is the subthreshold mobility, θ is the mobility degradation parameter, t ox is the thin-oxide thickness, and µ f et is a function of the voltages and hence the drain-source current. The differentiation of (64) with respect to I ds results in Then, the equations for µ f et and µ f et are given by and where all the model parameters and variables have been introduced above. Next, the expression for η f is given by where m η is the η floating-body parameter, and ∆ kink1 is similar to the equation described in (42) but by changing V dsep to V dse and multiplying by the ratio of w ef f /l ef f . The derivative of (68) is given by where ∆ kink1 is also similar to (44), but by changing the V dsep and V dsep terms to V dse and V dse respectively and multiplying by the ratio of w ef f /l ef f . V dse is expressed as and its derivative is given by where m ss is the V dse transition parameter, V sat = V gte , and hence, V sat = V gte . In conjunction with this, and where β = a sat · Vgt 2·v sth and β = a sat · V gt 2·v sth . The expressions for V gt and V gt can be referred from (25) to (28) and (30) to (33) as they are the same in both version 1.0 and 2.0 RPI poly-Si TFT.
The formula for the kink effect in the version 2.0 model is written as where V dsenew is a function of I ds and is similar to (70) by substituting V ds with V ds0 and V sat with V satnew . Differentiating ∆ kink with respect to I ds yields where V dsenew is similar to (71) by changing V ds to V ds0 and V sat to V satnew . The equations for V ds0 and V ds0 are given by and while the equations for V satnew and V satnew are given by where the variable Q s and its derivative are defined as VOLUME 4, 2016 and Finally, the formula for the leakage current and its derivative are similar to (46) and (47). The exact relationship between the variables that make up the leakage current can be referred from [26] and Fig. 4. This completes the formulation for the drain current of the TFT in the version 2.0 RPI poly-Si TFT model. A VinC LIM simulation for the TFT circuit in the model can be done similar to that in version 1.0, but by using the appropriate equations in (48) to (81) for the version 2.0 model. Do note that all the TFT equations presented thus far are for the n-type TFT model. For a p-type TFT model, the polarity of the voltage values should be inverted and the current value calculated is flipped in its polarity as well at the end of its calculation. As an illustration, the current equation for a p-type TFT model in basic LIM is given as while in the VinC LIM formulation it is given as , −V s (I ds ) = 0 (84) which are then used in (19).
This completes the formulation for the drain-source current for TFTs in both the basic LIM and VinC LIM formulations. The effects of the intrinsic capacitance and its current formulation will be discussed next.

C. BASIC LIM AND VINC LIM FORMULATIONS FOR RPI TFT INTRINSIC CAPACITANCE CURRENTS
While the main current in a TFT is the drain-source current discussed in the previous section, there are also minor current flows associated with the gate-drain terminals and gatesource terminals which can affect the accuracy of the simulation. These currents are modeled by intrinsic capacitances which exist between the two terminals, as shown in Fig. 2, and are contributed by the overlap capacitances which are formed according to the physical properties such as the thinoxide thickness, permittivity of the gate insulator, fringing factor, width of the model, etc., and the parasitic capacitances which depend on the model used for the TFT. These non-zero capacitances, labelled as C gd and C gs , connect the terminals of the TFT and hence, currents, labelled as I gd and I gs , are able to flow across the terminals when there are changes in the terminal potentials. In this section, the basic LIM and VinC LIM approach for these intrinsic capacitance currents are discussed. The intrinsic capacitances in a TFT can be represented as a branch capacitor connection in the LIM topology. In the basic LIM formulation, a fictitious inductor is inserted into the branch with the capacitor. This is illustrated in Fig. 5, where L ij is the inserted latency to the branch and V o is the voltage at the node between the branch capacitor and the fictitious inductor. Based on this, it can be written that (87) Substituting (87) into (86) yields which is the branch capacitance current updating equation for basic LIM. Thus, in a basic LIM simulation, (85) is applied first followed by (88) to solve for the branch capacitor current. In a TFT model, i and j can simply be replaced by the respective terminal nodes (e.g. g, d or s), and the value of C ij can be obtained from the model. The drawback of this approach is the additional latency component that is inserted to the respective branch, which may affect the simulation accuracy and stability.
For the VinC LIM formulation, the current equation for the branch capacitor can be derived by substituting the respective voltage terms with their LIM voltage equations and solving for the I n+1 ij current. Note that the insertion of fictitious inductances are not required for this approach. The branch current VinC LIM equation is given by and the other terms are as described in Section II. Equation (89) is an implicit VinC LIM branch capacitance current formulation. The original paper showing the complete derivation steps for a semi-implicit VinC LIM branch capacitance formulation can be referred to in [27].
Equation (89) can be used to solve for the intrinsic capacitance currents in a TFT by substituting i and j with the respective terminal nodes (e.g. g, d or s) and the value of C ij can be obtained from the model. It can be noted that the equation can be simplified into in a situation where terminal g is connected to a voltage source, or in a situation where terminal d or s is connected to a voltage source.
Considering a RPI poly-Si TFT model (for both versions 1.0 and 2.0), the overlap capacitances are given by where c gsos is the gate-source overlap capacitance, c gdos is the gate-drain overlap capacitance, w ef f overlap is the effective overlap width which depends on the area calculation method used in the TFT model, and c gso and c gdo are the source overlap capacitance factor and the drain overlap capacitance factor respectively, which can be calculated from the formula where l f , l d , t ox , and ε gate are the fringing factor, lateral diffusion into channel from source and drain, thin-oxide thickness, and permittivity of the gate insulator respectively. Then, for the parasitic capacitance, it is evaluated depending on the selected capacitance model. For example, a constant capacitance model which only depends on the device parameters can be calculated from the model equation as The total intrinsic capacitance of the RPI poly-Si TFT model is the summation of the overlap capacitance and the parasitic capacitance and can be written as and where C gs is the gate-source intrinsic capacitance and C gd is the gate-drain intrinsic capacitance.

D. PRACTICAL STEPS FOR ACCURACY IMPROVEMENT IN VINC LIM
One of the main differences between the basic LIM and VinC LIM algorithms is the utilization of a forward branchmarching scheme when solving for the branch currents in VinC LIM. Unlike the basic LIM branch equation that depends only on the voltages at the nodes and the current through that branch at the previous time step, the VinC branch equation depends also on the currents through the other branches which are solved for at the same time step. The forward branch-marching scheme alleviates this by using the most recently available current based on the order in which the currents are solved. However, this creates a situation where the solution is dependent on the order in which the branches are stepped through in the simulation. In a physical simulation, this order can be selected to start from the voltage and current sources, and propagating outwards to the rest of the circuit. However, determining this order is not a trivial task. To reduce the reliance on any particular order, in this work, a random order is used, but the order is alternated each time the branches are evaluated. In other words, if initially the order is selected as I a to I z , then on the next time step, the order will be flipped to be I z to I a . This idea is similar to other alternating direction algorithms, most notably to the one presented in [7] for transmission lines. Besides that, both basic LIM and VinC LIM are still inherently explicit methods, where the solution of the current time point is based on the solution at the previous time point. While this avoids the solution of large systems of simultaneous equations, such as those present in SPICE, the accuracy of the solution depends on the size of the time step used. In basic LIM, this is less of an issue, since the maximum time step to ensure stability is normally much less than that which is needed to obtain an accurate result, but since VinC LIM relaxes the stability criteria, accuracy can degrade when significantly larger time steps are used. In our experiments, we found that the accuracy can be improved by repeating the branch evaluations to obtain a better convergence during the forward branch-marching scheme. This can provide a better trade-off between speed and accuracy, compared to using a smaller time step, since only part of the solution process needs to be repeated. Numerical examples are presented in Section IV which illustrate this trade-off.

IV. RESULTS AND DISCUSSION
In this section, numerical examples are presented which show the application of the developed formulations in the simulations of TFT circuits. Silvaco's SmartSpice™ will be used as a benchmark. All simulations are performed on a Linux server with an Intel(R) Xeon(R) CPU E5-2699 v3 @ 2.30 GHz with 528 GB of RAM.

A. NAND CIRCUIT
In the first example, a NAND circuit is constructed using TFTs modeled with the version 1.0 RPI poly-Si model. This example is simulated as a verification of the basic LIM and VinC LIM formulations presented in the previous section. Fig. 6 shows the circuit schematic for this example, where Q 1 to Q 4 are TFTs, V ss = 5V is the power supply voltage, V g1 and V g2 are the input signals, and V 1 and V 2 are the output nodes of the circuit. In this example, the intrinsic capacitance currents are ignored and only the drain-source current of each TFT in the circuit is considered. Fig. 7 shows the simulated waveforms obtained from SmartSpice, basic LIM, and VinC LIM at the input and output nodes for the circuit. Fictitious capacitances of 10 −18 F are inserted at the two LIM nodes of the circuit for the LIM simulations. The simulation using basic LIM is performed at a time step of 1 ps, which is its maximum stable time step, while the simulation using VinC LIM can be performed at a timestep of 50 ps which is 50× the maximum stable time step in basic LIM, while still being both stable and accurate. This clearly shows the advantage of the VinC LIM formulation over the basic LIM formulation. All simulations are comparable in terms to accuracy to SmartSpice.

B. 7T1C PIXEL CELL CIRCUIT
In the second example, a single 7T1C cell, which makes up one of the color cells in an RGB pixel of an OLED display   [25] and [27]. Two simulations are performed using this circuit. In the first simulation, a version 2.0 RPI poly-Si model is used for the drain-source current and the intrinsic capacitance currents are ignored. Fig. IV-B shows the results obtained from SmartSpice, basic LIM, and VinC LIM at three distinct nodes in the circuit, V 1 , V 2 , and V 3 . In addition, Table 1 tabulates the number of time points, runtime, and RMS error in the simulations for different time step values for basic LIM and VinC LIM, based on the maximum stable time step in basic LIM, which is 0.224 ps. Note that some additional breakpoints are added by the simulator in order to capture the changing edges of the input waveforms accurately in all simulations.
From Fig. IV-B and Table 1, the result from basic LIM is only stable at its maximum stable time step value which is 0.224 ps for this circuit. On the other hand, the VinC LIM simulations remain stable even though the time step is increased up to 200 times the initial value. Comparing the runtime, for the same time step, basic LIM is faster than VinC LIM, due to the simplicity of its formulation. However, as VinC LIM is not limited by the stability criteria of basic Next, the same simulation using this circuit is repeated, but with the intrinsic capacitance currents taken into account. Fig. IV-B shows the simulation results from SmartSpice, basic LIM and VinC LIM. Additional branch loops are performed in VinC LIM as described in Section III-D to improve the accuracy. Table 2 shows a comparison of the runtime and accuracy of basic LIM and VinC LIM under different conditions. It is observed that the VinC LIM simulations are able to outperform basic LIM in terms of runtime at larger time steps, but with some degradation in accuracy. By adjusting the number of branch loops, the accuracy of VinC LIM can be improved, while still retaining its runtime advantage over basic LIM.

C. FULL TFT FLAT-PANEL DISPLAY CIRCUIT
In this third example, full TFT flat-panel display (FPD) circuits are simulated in VinC LIM and SmartSpice, to show the advantage of VinC LIM on large circuits compared to stateof-the-art commercial simulators. First, a full RGB pixel cell is simulated, where it consists of three 7T1C cells, each representing the red, green, and blue color in a pixel. A version 2.0 RPI poly-Si model is used for the drain-source current and the intrinsic capacitance currents are all taken into account.   and VinC LIM where a similar level of accuracy can be observed. It is to be noted that the simulation in VinC LIM was carried out using a time step 5000 times the maximum stable time step of basic LIM without any extra loop on the branch simulations. If necessary, the simulation speed can be further boosted by using a larger time step and a proper control on the number of branch loops. Then, large FPD circuits are constructed by using the RGB pixel cell to form the appropriate resolutions as in Table 3. Each circuit is simulated in VinC LIM and SmartSpice and the time spent per iteration is compared. It can be seen that VinC LIM greatly outperforms SPICE based simulators such as SmartSpice, especially on larger circuits. The runtime of VinC LIM scales almost linearly with the number of nodes in the circuit, while SPICE and its reliance on matrix based operations start to show very poor scaling on very large circuits. For circuits with more than 10 million nodes or devices, SmartSpice simulations were either unable to be completed even after seven days of runtime, or exceeded the total memory available on the systems. This is indicated by an entry of "dnc." in the table.