A Computationally Efficient Compact Model for Ferroelectric Switching With Asymmetric Nonperiodic Input Signals

In this article, we develop a Verilog-A implementable compact model for the dynamic switching of ferroelectric FinFETs (Fe-FinFETs) for asymmetric nonperiodic input signals. We use the multidomain Preisach Model to capture the saturated <inline-formula> <tex-math notation="LaTeX">$P$ </tex-math></inline-formula>–<inline-formula> <tex-math notation="LaTeX">$E $ </tex-math></inline-formula> loop of the ferroelectric capacitors. In addition to the saturation loop, we model the history-dependent minor loop paths in the <inline-formula> <tex-math notation="LaTeX">$P$ </tex-math></inline-formula>–<inline-formula> <tex-math notation="LaTeX">$E $ </tex-math></inline-formula> by tracing input signals’ turning points. To capture the input signals’ turning points, we propose an RC circuit-based approach in this work. We calibrate our proposed model with the experimental data, and it accurately captures the history effect and minor loop paths of the ferroelectric capacitor. Furthermore, the elimination of storage of each turning point makes the proposed model computationally efficient compared with the previous implementations. We also demonstrate the unique electrical characteristics of Fe-FinFETs by integrating the developed compact model of Fe-Cap with the BSIM-CMG model of the 7-nm FinFET.


I. INTRODUCTION
F ERROELECTRIC materials are widely used for memory applications due to their noncentro-symmetric structure, which leads to hysteresis in their electrical polarization (P) versus electric field (E) characteristics [1], [2], [3]. The high endurance exhibited by the ferroelectric films makes them promising candidates for storage and emerging applications such as neuromorphic computing, where a large number of history-dependent synaptic weights need to be tuned during the training process [4], [5], [6], [7], [8]. Further, the recent observations of ferroelectricity in doped hafnium oxide (HfO 2 ) materials have attracted much attention for the application of memory devices for neuromorphic applications [4], [9], [10], [11], [12]. The exhibit a limited endurance of ∼10 5 -10 7 cycles which may limit their application for neuromorphic training accelerators. However, FeFETs with such endurance are suitable for memory and inference accelerators. Thus, it is essential to have a computationally efficient compact model of the ferroelectric capacitor (Fe-Cap) for such large-scale synaptic devices by accurately capturing the history-dependent minor loops. The physics-based models for ferroelectric capacitors have been developed to capture the accurate physical behavior of the Fe-Cap [13], [14]. Recently, Wang [15] described how to build such hysteresis in the circuit simulator using internal state variables. A compact model for the Fe-Cap that captures history-dependent minor loops has already been developed in [16] and [17]. The model was based on tracing the evolution of turning points, where each turning point is stored in an array. The model further uses the memory wipe-out method [17] to reduce the storage of turning points to improve the computational efficiency of the model. However, for certain combinations of asymmetric input signals where the minor loops of P-E look like a spiral hysteretic curve, i.e., the upper and lower voltage turning points keep reducing in magnitude, the model requires large storage of each turning point, which is highly undesirable for large circuit simulation of neural network training.
In this work, we develop a computationally efficient compact model for Fe-Cap with the history-dependent minor loops, which eliminates the stringent requirement of a storage array for any combination of asymmetric input signals. To switch the polarization versus voltage (P-V) loop from forward sweep to reverse sweep or vice-versa, the model is required to trace the turning points in the input signal. In this work, we propose an R-C network to capture the turning points in the input signals. We make the approximations while developing the Fe-Cap model that is discussed in Section III, which eliminates the requirement of a storage array for any combinations of asymmetric input signals. The approximation made in the model works very well and accurately captures the experimental results for the complicated input signals. Finally, we demonstrate the Fe-FinFET characteristics by solving the self-consistent solution of the proposed model of Fe-Cap with the industry-standard BSIM-CMG model for 7-nm FinFET.
1937-4151 c 2022 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.

II. MODEL DEVELOPMENT OF FE-FINFET
A diagram of the metal-ferroelectric-metal-insulatorsemiconductor (MFMIS) type of Fe-FinFET is shown in Fig. 1. The internal metal gate present in the gate stack of the Fe-FinFET forms an equi-potential surface between the interfacial oxide and the ferroelectric layer, which helps to consider the two different circuit entities, i.e., the ferroelectric capacitor (C FE ) and the underlying conventional FinFET (M) for the modeling of the Fe-FinFET as shown in Fig. 2(a). V int is the internal metal gate voltage, and V G , V S , and V D are the gate, source, and drain terminal biases, respectively. As shown in Fig. 2(b), for the underlying conventional FinFET, we use the industry-standard BSIM-CMG Verilog-A model [18], which accurately considers the short channel effects and the quantum mechanical effects. We utilize a 7-nm open source predictive process design kit (PDK) for the model parameters of the conventional FinFET [19]. For the ferroelectric capacitor (C FE ), we propose a Verilog-A-based dynamic Preisach model including the history-dependent minor loops, which we discuss in the next section in detail. For each applied gate bias, the SPICE simulator solves self-consistently for the gate charge of the underlying FinFET (Q G ), i.e., Q G = Q FE for each V G to obtain the electrical characteristics of the Fe-FinFET.

A. Modeling of Switching Dynamics in Fe-Cap
To implement the switching dynamics in the Fe-Cap, we first calculate the auxiliary voltage (V aux (t)) to which the ferroelectric dipoles respond. The SPICE model to obtain (V aux (t)) across the ferroelectric layer is shown in Fig. 3. The V aux (t) can be expressed in terms of the R aux -C aux delay [14], [20] given by where V in (t) is the applied voltage across the ferroelectric layer and τ v is the relaxation time for the auxiliary voltage. For the quasi-stationary case, the auxiliary voltage equals the applied input voltage, i.e., V aux = V in . To solve the above transcendental equation in a circuit simulator, we define an extra internal node named aux in the Verilog-A code. The right-hand side of (1) is assigned to the voltage at node aux which is solved  self-consistently in the circuit simulator resulting in the actual auxiliary voltage V aux (t). Now, using the auxiliary voltage (V aux ), we calculate the auxiliary polarization (P aux ) using the Preisach model. The Preisach model can be expressed as a superposition of simple hysteresis units, each of which contains a rectangular hysteresis loop. The exact mathematical implementation of the Preisach model is as follows: is the direction operator. The solution of such equation requires numerical evaluation of double integral, which is a complex and time-consuming procedure.
In [17], a computationally efficient model for a ferroelectric capacitor was demonstrated where the behavior of the ferroelectric capacitor is mathematically described using the hyperbolic tangent function. For the forward loop, P aux↑ using hyperbolic tangent function is given by where V c and P s are the coercive voltage and the saturated polarization value of the ferroelectric material, respectively. m ↑ is the slope and P off↑ is the offset polarization value for the forward sweep. Both m ↑ and P off↑ are required to obtain the forward sweep minor loops and for the forward sweep saturated loop m ↑ = 1 and P off↑ = 0. The forward paths  (3) is given as where t fe is the ferroelectric thickness and P r is the remnant polarization of the ferroelectric material. Similarly, for the reverse sweep, the auxiliary polarization P aux↓ is written as Again m ↓ and P off↓ are required for the backward sweep minor loops and become m ↓ = 1 and P off↓ = 0 for the saturated backward sweep loop. The backward paths obtained from (3) are indicated by the red arrows in Fig. 5. Both m ↑ /m ↓ and P off↑ /P off↓ are calculated using the polarization history of the material which is discussed later in this section. Now the actual ferroelectric charge density(Q(t)) expression with the inclusion of background permittivity of the ferroelectric material [14], [20] is rearranged as where τ p is the relaxation time of polarization, k n is the coupling coefficient, and ε fe is the background permittivity of the ferroelectric material. For the quasi-stationary case, the actual ferroelectric charge in the Fe-Cap is equal to the auxiliary polarization in addition to the dielectric charge component caused due to the internal field. Still, to obtain the complete Q−V aux loop, the model needs to switch the polarization from forward sweep to backward sweep or vice-versa whenever the turning points appear in V aux .

B. Tracing Turning Points
We propose an R delay -C delay circuit to trace the turning points in V aux as shown in Fig. 4(a). We define an internal node delay in the Verilog-A code with an aim to find out voltage at node delay, where R delay and C delay are the circuit elements which are required to produce a delayed version of V aux at node delay. To achieve the delayed auxiliary voltage, V aux at each bias point is fed to the R delay -C delay circuit in the form of a current as shown in Fig. 4(a). The current source is a voltage-dependent current source with the proportionality constant α, i.e., I delay = −αV in . The unit of α is −1 . For explanation purposes, we consider the case where the input voltage is time independent. After solving Kirchhoff's current law (KCL) and differential equation in the R delay -C delay circuit, we obtain the voltage V delay for the time-independent input voltage case as We need the delayed version of V aux at node delay, thus, we fed the same V aux through the current source with α = 1 −1 .
The voltage at node delay (shown by a solid red line) is delayed by t d = R delay C delay from the actual auxiliary voltage signal (shown by a solid gray line) as shown in Fig. 4(b).
In this work, we keep R delay = 1 and C delay = 1fF, which provides a very small value of t d in the order of f s. However, to understand the topology, we have shown an exaggerated difference in the auxiliary and delayed voltages in Fig. 4(b), which is not the case during the real simulations. Note that the value t d is very small compared to the relaxation time or domain switching time of the ferroelectric material. Thus, the RC network approach does not lead to significant change in the simulated characteristics from the actual results. Note, the differential equation to solve V delay (t) for any timedependent arbitrary input signal is being solved by the SPICE simulator. Now, to obtain the turning points, we set the condition during the backward sweep as V delay (t) < V in (t). Whenever this condition is satisfied, we store the auxiliary voltage values (denoted as by A, C, E, and G) as the lower turning point of the backward sweep. Similarly, we set the condition during the forward sweep as V delay (t) > V in (t) and when the condition becomes true, we store the auxiliary voltage values (denoted as B, D, and F) as the upper turning points for the forward sweep.

C. Calculation of m ↑ /m ↓ and P off ↑ /P off ↓
In this section, we calculate m ↑ /m ↓ and P off↑ /P off↓ using the turning points obtained in the previous section, which are essential to achieve the history-dependent minor loops. The upper and lower voltage turning points are denoted by V auxu and V auxl and the corresponding polarization values are denoted by P auxu and P auxl , respectively.
In the forward bias, we write the polarization expression at the upper turning point as Similarly, at the lower turning point in forward sweep, we write the polarization as Thus, using (8) and (9), we obtain m ↑ and P off↑ as Similarly, using the backward sweep polarization expression at upper and lower turning points, we can obtain m ↓ and P off↓ for the backward sweep as Note the values m ↑ and m ↓ are equal to 1 for the saturated loop and become less than 1 for the minor loops. And the values of P off↑ and P off↓ are equal to 0 for the saturated loop and become nonzero for the minor loops.   Fig. 4(b). The Q-V aux characteristics obtained from the developed compact model is shown in Fig. 5. We initialize a backward sweep (Type = 0), and m ↑ /m ↓ = 1 and P off↑ /P off↓ = 0 for the saturated loop. Type = 0 and Type = 1 in the model flowchart imply the backward and forward sweeps, respectively. So, when the applied input signal varies from 3 to −3 V, P aux↓ is calculated from (5) with m ↓ = 1 and P off↓ = 0 until the applied bias reaches the first turning point A as shown in Fig. 4(b). Once the condition V delay (t)−V in (t) < 0 is satisfied during the backward sweep,the model now switches to the forward sweep, i.e., Type = 1. At this bias point, the model will store the polarization as P auxl↑ and voltage as V aux l ↑ as the turning points. Thus, P aux↑ is calculated from (3) using m ↑ and P off↑ from (10) and (11) using the stored P aux l ↓ and V aux l ↓ at point A, and uses P auxu↑ at V auxu↑ = 3V as another turning point. The forward sweep continues until the next the turning point B. When the model recognizes a turning point at B, it will wipe out the previously stored P auxu↑ and V auxu↑ and overwrite the new values of turning points obtained at B. As the polarization P aux↑ does not reach the saturation polarization value during the forward sweep, thus, for the next backward sweep, the model works in the minor loop with m ↓ < 1 and nonzero P off ↓ . Next, the model calculates the m ↓ and P off ↓ using new turning points obtained at B and the previously stored P auxl↓ and voltage as V auxl↓ at point A until point C. The model further experiences the turning point at point C and works in the forward minor loop, where it overwrites the previous turning points as P auxl↓ and V auxl↓ stored at A. Next, for the forward loop between point C and D, it uses the recent turning points stored at points B and C. Similarly, for the backward loop starting at D it uses the turning points stored at C and D and overwrites the previous turning points stored at B. Now, when the bias takes a path between point D and E, the model uses the turning points stored at C and D until it reaches point C. Once, the bias reaches point C, we flag the backward sweep as path = 0 and use the modified expression for m ↓ and P off ↓ using (14) and (15), respectively. For the modified m ↓ and P off ↓ , the model uses P aux and V aux as upper turning points stored at C. Further, the model takes the second turning point as a maximum negative value of applied bias, since the polarization eventually converges to the saturation polarization at the maximum negative value of applied bias. Thus, the model changes slope for the downward sweep and follows the path between the points C and E. However, in the previous models, to calculate the slope of path = 0, they use the turning point stored at A and B. Therefore, they require the storage of multiple previous turning points in the form of an array. Hence, using this approximation, we eliminate the need for storage of turning points while switching from inner loop to minor loop. The modified expressions of m ↓ and P off↓ for path = 0 are given by

D. Modeling Flow for History-Dependent Minor Loops
At point E, the model enters in forward loop and at this bias point the model notices that the backward sweep for the last bias point was flagged as a path = 0. For the path = 0 condition, the model calculates the slope using the turning points stored at E and the second turning point predicted at the maximum positive applied bias (V max ). As mentioned earlier, the polarization has to converge to the saturation polarization value for maximum positive applied bias. Thus, the modified definitions of m ↑ and P off↑ for path = 0 are given by The model uses the above definition of m ↑ and P off↑ until the bias reaches F. Again, at point F, the model faces the turning point and enters the backward sweep, it uses the turning points stored at point E and F to calculate m ↓ and P off↓ until bias reaches at G. At G, model uses the stored turning points at F and G for the forward sweep till point F. After point F, the model flags as path = 0 and hence, uses the first turning point stored at F and the second it predicts from the maximum applied bias, i.e., 3 V.

IV. MODEL VALIDATION AND DISCUSSION
Now, the compact model developed for the Fe-Cap in the last section is validated for the different complicated asymmetric input signals against the experimental results of the 10-nm HfO 2 -based Fe-Cap as shown in Fig. 7. The calibrated parameters for the Fe-Cap are listed in Table I. We validate our model for three different input signals which are shown in Fig. 7. We consider the quasi-stationary case, where the relaxation time constants are much smaller than the applied input pulse width. The model shows excellent agreement with the experimental data for all three cases. Further, the response of remnant polarization of the ferroelectric material (P r ) for different pulse widths is demonstrated in Fig. 8. The simulated P r response for different pulse widths and amplitudes shows good agreement with the experimental results [16]. For validation purposes, we consider the relaxation time of the auxiliary voltage (τ v ) as 1.5 μs whereas we ignore the relaxation time of polarization (τ p ), since the dielectric relaxation in the HfO 2 layer dominates only at very high frequencies of operation (>1 GHz) [21], [22].  Next, we demonstrate the transfer characteristics of the Fe-FinFET by jointly solving the developed compact model of the Fe-Cap with BSIM-CMG model of the FinFET. To obtain the drain current characteristics of the Fe-FinFET, we use the model parameters of a 7-nm conventional FinFET and keep t fe = 3 nm. Fig. 9 shows the history-dependent drain current characteristics of the Fe-FinFET for three different cases. The biasing scheme is shown in the inset figure of Fig. 9. Each case has the same initial and final bias point, however, it follows a different path to reach the final bias point. The drain current in Fig. 9 is obtained during the read process for all three different cases. The drain current of the Fe-FinFET heavily depends on the history-dependent polarization states and shows significant difference for each case. Fig. 10(a) shows the I DS − V GS for the program and erase voltages at V DS = 0.1 V. For programming voltages, we increase the amplitude of V PROG from 1.2 to 3 V with the increment of 0.1 V as shown by the blue color in Fig. 10(b). We keep the pulse width of each pulse as 50 μs. The drain current for each polarization is read after each programming pulse and plotted in Fig. 10(a). Similarly, we apply the erase pulse from −1.2 to −3 V in steps of −0.1 V shown by the green color in Fig. 10(b). The drain current for each polarization state after each erasing pulse is shown in Fig. 10(a).
Furthermore, to compare the efficiency of the proposed model in terms of memory requirement and computational time as compared to the prior implementation [16], we have performed device as well as circuit simulations. At the device level, we simulate the Fe-Cap model for the specific asymmetric input signal shown in Fig. 11(b). We consider the asymmetric input pulse for which the model produces the spiral hysteretic Q − V aux curve. The obtained Q − V aux from both the models are shown in Fig. 11(a). For such an input pulse, the model presented in [16] requires storage of every turning point which increases the computational time and the memory requirement. Table II shows the comparison of simulation time and memory allocation required for the applied input pulse for both the models. Our model requires significantly less computational time than the model which requires memory storage, though we do not see any significant   difference in memory allocation for the device-level simulations. Further, the performance gains while using the proposed model are more evident at the circuit level. We performed an array-level simulation for a 100 × 100 crossbar array of Fe-FinFETs; these results are shown in the bottom two rows of Table II.

V. CONCLUSION
We have formulated a computationally efficient compact model for the Fe-Cap and implemented an Fe-FinFET model with the conventional FinFET in Verilog-A code for largescale circuit simulations. The proposed model is validated for the history-dependent minor loops with the experimental results for the nonperiodic asymmetric input signals. In comparison with the previous compact models, we eliminate the need for storage of previous turning points. The model developed for the Fe-Cap further can be implemented with other industry standard models of different FET structures and may serve as the platform for design exploration of FeFETs for several unconventional applications.