Physical Compact Model for Three-Terminal SONOS Synaptic Circuit Element

A well-posed physics-based compact model for a three-terminal silicon – oxide – nitride – oxide – silicon (SONOS) synaptic circuit element is presented for use by neuromorphic circuit/system engineers. Based on technology computer aided design (TCAD) simulations of a SONOS device, the model contains a nonvolatile memristor with the state variable Q M representing the memristor charge under the gate of the three-terminal element. By incorporating the exponential dependence of the memristance on Q M and the applied bias Vfor the gate, the compact model agrees quantitatively with the results from TCAD simulations as well as experimental measurements for the drain current. The compact model is implemented through VerilogA in the circuit simulation package Cadence Spectre and reproduces the experimental training behavior for the source – drain conductance of a SONOS device after applying writing pulses ranging from (cid:1) 12 V to þ 11 V, with an accuracy higher than 90%.

depicts the equivalent circuit diagram for the wellposed compact model of a three-terminal SONOS synaptic circuit element [34][35][36] for neuromorphic applications, which consists of three components. The first is found on the left vertical branch, composed of a series capacitor C M and an extended memristor [41] M, defined as M ¼ f ðQ M , VÞ and (1a) where we are not directly interested in the current I M through the memristor, which is in general too small to be measured experimentally, but rather how the memristor charge controls the longterm SONOS dynamics for synaptic weight modulation or NVM through time dependence which is discussed further. The second vertical branch on the right has a series resistor R and a capacitor C, which control the short-term dynamics represented as and in turn the drain current, i D is Equation (4) is essentially quasistatic, [42,43] without an explicit time dependence due to the small time constant τ ¼ RC ( 1ns. We performed circuit simulations based on a SONOS device with W ¼ 1.2 μm, L ¼ 7 μm, μ ¼ 350 cm 2 V À1 s À1 , and C ¼ 26 fF [18,33] as shown in Figure 1b combined with R ¼ 10 8 Ω, chosen to manage the simulation speed with a time constant τ ¼ 2.6 μs, which permits the use of a timestep as large as 1 ns in Cadence Spectre to examine both short-and long-term dynamics. For most practical circuit simulations, the short-term dynamics can be safely eliminated by making Q a parameter instead of a state variable, [18] as discussed in Section S2, Supporting Information. The last component comprises two identical variable resistors, derived from v D /i D as where χ stands for the percentage of Q M contributed by the silicon channel. The memristor conductance is a function of the Q M history as described later. We performed TCAD simulations to determine the quantitative dynamics of the state variables Q M and Q, representing charge transport via tunneling and electron drift diffusion, [44] because they are not easily experimentally accessible. A memristor circuit element was necessary to describe the dynamics of Q M , whereas a linear R-C circuit was suitable for Q. The results in Figure 2 provided the information used to derive the compact model in Figure 1. To focus on the dynamics of synaptic weight modulation, the drain and the source were grounded and denoted as "Si," which eliminated the channel current. A bias of V ¼ 8 V was applied to the gate at t ¼ 0, with Figure 2 representing two times after bias application, t ¼ 1 μs and t ¼ 200 ms. The resistor in the short-term branch (R) in Figure 1 was not included because C was fully charged within 1 μs in the TCAD simulations. The amount of charge under 8 V bias at t ¼ 1 μs (left, Figure 1 (a)) was Q ¼ 17.5 Â 10 12 e À cm À2 , corresponding to 2.8 μC cm À2 , which shows good agreement with the expected value of 2.5 μC cm À2 calculated from 8 V on a capacitor with 0.31 μF cm À2 (3 nm SiO 2 , 8 nm Si 3 N 4 , and 4 nm SiO 2 ) visualized schematically in Figure 2b. The area density of Figure 1. a) The equivalent circuit diagram for a three-terminal SONOS synaptic circuit element based on physical state variables Q M and Q. The state variable Q M is the total amount of charge that has passed through an effective memristor M, equal to the charge on the series capacitor C M (the branch on the left), and controls the slow dynamics. The memristance is a function of both the state variable Q M and the applied voltage, v G -V T -v M , with strong nonlinearity. The parallel branch on the right, composed of a series resistor (R) and capacitor (C), is responsible for the short-term dynamics with a state variable Q, which represents the typical switching behavior of a MOSFET. For simplicity, the capacitor C is described as a linear capacitor, which holds when v G -V T > 0, where the influence of charge depletion is negligible. R describes an effective resistance including the contact resistance of the gate electrode, the scattering by the ionized impurities in the depletion region, etc. The conductances and currents through the channel of the SONOS device are represented by variable resistors with identical resistances R CH and depend on the two state variables Q and Q M . b) Schematic 2D cross section perpendicular to the gate of the SONOS device used in the TCAD simulations and experiments (the study by Agarwal et al. [33]), with the doping concentration shown in different colors. Inset shows three exemplary programmed states with colors indicating the trapped electron concentration in the Si 3 N 4 layer.
www.advancedsciencenews.com www.advintellsyst.com electrons was used as the unit for Q and Q M for more intuitive analysis in this article. A slightly larger amount of charge compared with the expected capacitor value is attributed to the negative threshold voltage V T ¼ À0.01 V and nonlinear Q-V curve near the depletion region (v G ≅ V T ) of a MOSFET. The corresponding band diagram across the material stack is shown in Figure 2c. The trap layer (Si 3 N 4 ) remains electrostatically neutral until t ¼ 1 μs with Q M ¼ 0 and provides a tunneling barrier thickness of 3 nm due to the 9 eV forbidden gap of SiO 2 . After allowing enough time for tunneling, at t ¼ 200 ms (Figure 2 right column), Q M ¼ 5.8 Â 10 12 e À cm À2 , and Q G ¼ 20 Â 10 12 h þ cm À2 , hence, QÀχQ M ¼ 14.2 Â 10 12 e À cm À2 because of charge neutrality (Q G ¼ Q M þ QÀχQ M ). Consequently, the synaptic weight of the SONOS decreased based on the net charge changing from 17.5 Â 10 12 e À cm À2 at t ¼ 1 μs to 14.2Â 10 12 e À cm À2 at t ¼ 200 ms. In other words, 57% of the increased Q M (0!5.8 Â 10 12 e À cm À2 ) is contributed by the electron charge on C (¼ QÀχQ M : 17.5 Â 10 12 e À cm À2 !14.2 Â 10 12 e À cm À2 ), realizing the synaptic weight depression. Although the tunneled charge is 5.8 Â 10 12 e À cm À2 , intuitively implying Δ(QÀχQ M ) ¼ À5.8 Â 10 12 e À cm À2 , the loss of charge on C is spontaneously complemented by the gate electrode (Q G : 17.5 Â 10 12 h þ cm À2 ! 20 Â 10 12 h þ cm À2 ) due to the electrostatic force from the trapped charge located at the middle of dielectric stack (O-N-O). The Supplementary Video provides a dynamic circuit illustration of this mechanism and the relevant discussion on χ follows. The underlying physics of the memristor in our compact model is revealed in the band diagram at t ¼ 200 ms (right panel of Figure 2c). The slightly lifted conduction band edge (or electrostatic potential) due to the injected negative space charge (trapped electrons) alters the tunneling thickness so that an electron, from the conduction band edge of silicon interfaced with the SiO 2 tunneling oxide, encounters a thicker tunneling barrier of 5 nm compared with 3 nm at t ¼ 1 μs when Q M ¼ 0. Consequently, Figure 2. Synaptic weight dynamics modeling of a SONOS device with a memristor and two capacitors. The amounts of charge on the capacitors and the corresponding band diagrams are taken from TCAD simulations with a bias of V ¼ 8 V as an example. Two distinct points of time at t ¼ 1 μs (left) and t ¼ 200 ms (right) after the bias is applied are shown. Each 'À' symbol represents roughly 2.5 Â 10 12 electrons per cm 2 (e À cm À2 ). a) Equivalent circuit diagrams of a SONOS layer stack with two terminals: the gate and silicon channel (source and drain are grounded). A memristor M and a capacitor C M account for the synaptic weight modulation based on the state variable Q M . When t ¼ 200 ms, Q M increases so that Q-χQ M as well as the synaptic weight (conductance across source and drain) are reduced. b) Schematic physical structure of the SONOS stack. The thicknesses of the ONO layers are 3, 8, and 4 nm, respectively, whereas the silicon channel and gate are much larger. c) Corresponding band diagram across the SONOS layer at t ¼ 1 μs compared with t ¼ 200 ms, which exhibits a thicker effective tunneling barrier (5 nm) that is larger than the tunneling oxide (3 nm). The memristor M with the state variable Q M equivalently reflects the varying tunneling barrier thickness.
www.advancedsciencenews.com www.advintellsyst.com the tunneling current decreases, that is, Figure 2 presents results at two discrete times, but the memristance changed smoothly and continuously with the state variable Q M . Figure 3a,b shows the simulation results of Q M as a function of time under constant biases: þ9, þ8, and þ7 V in (a) and À8, À7, À6, and À5 V in (b). The dynamics are qualitatively reminiscent of an R-C circuit, where the initial current, i t¼0 , is equal to V/R and Q increases inverse exponentially until it is saturated at Q ¼ CV but is quantitatively different. For example, the two curves with þ9 and þ8 V in Figure 3a show that þ9/[dQ M /dt] t¼0 from the former is not identical to þ8/[dQ M /dt] t¼0 from the latter, immediately identified by the slope of the curves at t ¼ 0. Therefore, the resistive component in series with the capacitor has a voltage dependence. Second, the slope of the tangent decreases exponentially (not linearly) with Q M , which is evident from the replotted curves on a log scale in Figure 3c, demonstrating that the resistance is also a function of Q M , as forecast by Equation (1). The y-axis represents Q M,max À Q M for V > 0 and Supporting Information). For a usual RC circuit, ln(Q À CV ) for V < 0 (or ln(CV À Q ) for V > 0) is a linear function of time with a slope À1/ τ, where τ is the time constant . In (a), the state variable Q M increases with time and saturates at Q M ¼ Q M.max ¼ 8 Â 10 12 e À cm À2 , which is defined by the trap concentration (10 19 cm À3 ) chosen for the 8 nm-thick Si 3 N 4 layer. c) Logarithmic-scale plots of Q M for V > 0 (or Q M,max À Q M for V < 0), conveying the identical information of (a) and (b), but emphasizing the nonlinearity of the response. The decaying slopes (¼À1/τ) with respect to time show that the characteristic time is not constant. d,e) DRMs extracted from (a) and (b), respectively. The nearly linear relations between log(dQ M /dt) and Q M demonstrate that the memristance depends exponentially on Q M . The compact model linear approximations to the DRMs are denoted with dotted lines. The eccentric behaviors observed for 9 V with Q M > 5 Â 10 12 e À cm À2 and À8 V with Q M < 10 12 e À cm À2 are due to Q M approaching Q M,max ¼ 8 Â 10 12 e À cm À2 and Q M,min ¼ 0, respectively. f ) Memristance as a function of voltage. Upward and downward triangles represent the memristance values at two extremes, Q M ¼ 0 and Q M ¼ 8 Â 10 12 e À cm À2 , respectively, from (d and e). The overall trend of the memristance M versus Q M and V exhibits the dependence on the state variable (governed by a parameter γ) and a strong nonlinearity with V (governed by a parameter β) in Equation (8).
www.advancedsciencenews.com www.advintellsyst.com (¼ RC, see Equation (S5)%(S8), Supporting Information for more details). However, the observed results exhibit significant nonlinearities for all seven biases, demonstrating that the resistance is a function of Q M , which is the memristance. The dynamic route maps (DRMs) in Figure 3d,e describe the dependence of the dynamics on the state variable. The Q M dynamics results obtained from our physics simulations were fitted by a M-C circuit as shown by the dotted lines in Figure 3d,e, with the dynamical equation as follows.
This equation breaks down when Q M approaches the extrema at 0 or Q M,max , which both require extreme times to reach.  (8) where α ¼ 50, β ¼ 70, γ ¼ 50, δ ¼ 25, and n is 0 for V > 0 and 1 for V < 0. The dependences on V and Q M are described by β and γ, respectively. In addition, the sensitivity of the memristance on Q M (slopes of the dotted lines in Figure 3d,e) also exponentially increases with V, which is reflected by the fitting parameter δ. Although the curves from the analytic equation denoted by lines produce some errors compared with the physics simulation results (symbols), it can account for all possible combinations of Q M and V so that it can stand alone as a circuit element. A discontinuity in M between V ¼ 0À and V ¼ 0þ is observed when Q M ¼ Q M,max , for which the physical mechanism is attributed to a change in the charge tunneling dynamics. While a larger Q M provides a thicker tunneling barrier in the case of trapping into Si 3 N 4 (V > 0), it fosters detrapping when V < 0 because the electrostatic potential of the trap is higher compared with when Q M ¼ 0. Figure 4a,b shows the familiar pinched hysteresis loop characteristic of memristors with an applied sinusoidal voltage. In the present system, both loops are traversed clockwise, which is different from most familiar memristors that show both clockwise and counterclockwise trajectories depending on the sign of the voltage. Although the I M -V characteristic in Figure 4b appears to be that for a volatile memristor, the expanded view in Figure 4c shows that there are two different slopes at zero crossing, showing that the memristor is nonvolatile but that the states near the zero crossing have a very large resistance compared with those at high voltage amplitudes. Figure 4d shows that there is an avoided zero crossing when there is a capacitor in series with the memristor, which would be the case for an experimental measurement of this system as the capacitance is intrinsic to the structure. Blue solid curves in Figure 4a  The compact model and physics simulation present good agreement, although they have a quantitative mismatch owing to the simple analytical form of Equation (8) that can be improved at the cost of complexity. The red solid curves in Figure 4a,b represent the long-term gate current, I G,0 , which is the sole way to deduce I M through experimental measurements. The shape of I G,0 is almost identical to I M , but with different magnitudes, approximately half of I M , which is consistent with the illustration in the Supplementary Video and justifies the parameter χ ranging from 0 to 1. The detailed procedure for extracting the long-term gate current (I G,0 ) from the total gate current (I G,0 þ I G,∞ ) is available in Figure S2, Supporting Information. The memristor has both nonlinearity and asymmetry in its I-V characteristic, which leads to a "history erase effect" [45] or loss of initial condition, as shown in simulations of our compact model under sinusoidal voltage ( Figure S3, Supporting Information). We next calibrated the compact model extracted from TCAD simulations using experimental data from a SONOS device. [33] As an accurate compact model requires dynamical information from the target, we designed a protocol to extract the state variable from measurements that can actually be performed, [46] as shown in Figure 5, for a general three-terminal synaptic device. [30][31][32]47,48] For SONOS, the measured channel currents were converted to the change in threshold voltage (V T ) so that the desired state variable Q M could be obtained. Figure 5a shows the channel current of a SONOS device potentiated by À12 V and suppressed by þ11 V from the work of Agarwal et al. [33] Utilizing the rigid shape of i D -v G regardless of the shift in V T , ΔV T can be deduced from the current based on the measured i D -v G curve at an arbitrary state, as shown in Figure 5b, to obtain Figure 5c.
While a simple correlation between the amount of fixed charge (ΔQ NOT ) at the interface of the silicon channel and SiO 2 and the shift in threshold voltage (ΔV T ) holds from ΔV T ¼ ΔQ NOT /C for a MOSFET, when it comes to a SONOS device with ΔQ M (most charges exist in Si 3 N 4 ) instead of ΔQ NOT , the relation becomes ΔV T < ΔQ M /C. [19,49] Due to electrostatic equilibration, the trapped charge closer to the channel results in the larger ΔV T , mainly dictated by the amount of charge at the channel. Likewise, trapped charge closer to the gate eventually draws more charge from the gate rather than from the channel (because the gate supplies the charge to the channel, a mechanism animated in Video S1, Supporting Information), so that ΔV T is reduced. To estimate χ, ranging from 0 to 1, for a generalized correlation of ΔV T ¼ χΔQ M /C, a series of additional TCAD simulations were performed with three different scenarios for the trap concentrations in Si 3 N 4 : 1 Â 10 19 , 3 Â 10 19 , and 10 Â 10 19 cm À3 . As shown in Figure 6a, all three cases require %ΔQ M ¼ 2.5 Â 10 12 e À cm À2 to induce ΔV T ¼ 1 V; hence, the resultant χ is calculated to be 0.77, because C ¼ 0.31 μF cm À2 . For the lowest trap concentration, 10 19 cm À3 (Figure 6b), a slightly larger amount of trapped charge, ΔQ M ¼ 2.65 Â 10 12 e À cm À2 , is required to achieve ΔV T ¼ 1 V (i.e., χ ¼ 0.73) and the incremental efficiency becomes worse for the larger Q M (>4 Â 10 12 e À cm À2 ). This is because of the limited capacity of the trap layer so that for Q M > 4 Â 10 12 e À cm À2 , for instance, the trapped charges reside not only near the interface of Si 3 N 4 , but also in the bulk region of Si 3 N 4 , which is closer to the gate and induces a smaller ΔV T . With the higher trap concentrations in Si 3 N 4 (3 Â 10 19 cm À3 and 10 Â 10 19 cm À3 ), as shown in Figure 6c,d, the newly trapped electrons always occupy the interface between Si 3 N 4 and the tunneling SiO 2 layer, so that the efficiency remains the similar with χ ¼ 0.77. Based on the analysis with different trap concentration scenarios, a conversion factor of 3 Â 10 12 e À cm À2 V À1 was chosen to extract Q M from the experimental measured data of ΔV T in Figure 5c. Figure 7a,b shows the resultant DRMs from Q M of the experimentally measured channel currents versus the number of square writing pulses on the gate [33] with seven different bias conditions: þ11, þ10, and þ9 V for depression and À12, À11, 10, and À9 V for potentiation. Although the experimental points are noisy compared with those from TCAD simulations, the exponentially changing memristance with Q M and the nonlinearity with the bias voltage are present, similar to our observations from physics simulations. The triangles (blue for Q M ¼ 0 and red for Q M ¼ Q M,max ) in Figure 7c,d correspond to the bias voltages in Figure 7a,b in a similar manner to Figure 3f. The best fitting parameter set was found to be (α, β, γ, δ) ¼ (70 000, 7, 50, 1) for the experimental SONOS device compared with (α, β, γ, δ) ¼ (50,70,50,25) from the SONOS physics simulations. The significantly larger α represents a much slower tunneling by three orders of magnitude for the experimental SONOS devices at V ¼ 10 V, implying that the tunneling mass in the TCAD simulations is too low, where the default value of Synopys Sentaurus 0.36 m 0 was used. However, the memristance at V ¼ 0 (see Figure 3f ) shows the inverse trend, such that M ¼ 1.41 Â 10 20 Ω cm 2 from the TCAD model and M ¼ 1.98 Â 10 13 Ω cm 2 from the experimental SONOS data, implying that the tunneling mass may need to be smaller than Figure 6. Relationship between V T and Q M under various conditions of Si 3 N 4 trap density to obtain Q M from the experimentally obtained V T data. a) V T linearly increases by 1 V for every ΔQ M ¼ 2.5 Â 10 12 e À cm À2 . b) A detailed analysis for the case of trap density equal to 1 Â 10 19 cm À3 . The slope of 'Q-χQ M ' versus V relation is consistent with the expected capacitance of the ONO layer (%0.31 μF cm À2 ), which governs the correlation between ΔV T and χΔQ M . For a MOSFET, ΔQ M ¼ 1.93 Â 10 12 e À cm À2 is required to induce ΔV T ¼ 1 V (¼ [0.31 μF cm À2 ]/[1.6 Â 10 À19 C]) with χ ¼ 1 when the fixed charges exist at the interface between the silicon channel and SiO 2 . For the SONOS device, for which the trapped charge is dominantly in the Si 3 N 4 layer, χ ¼ 0.73 was derived from the simulation results: ΔQ M ¼ 2.65 Â 10 12 e À cm À2 for every ΔV T ¼ 1 V. c) When the trap density increases to 3Â10 19 cm À3 , the trapped charge has a narrower distribution, spatially closer to the channel rather than the gate. As a result, Δ(Q-χQ M ) contributes 77% of ΔQ M (i.e., χ ¼ 0.77) and the conversion factor slightly decreased to 2.52 Â 10 12 e À cm À2 V À1 . d) Further increase to 10 Â 10 19 cm À3 does not create a notable change compared with 3 Â 10 19 cm À3 from (c) (i.e., χ ¼ 0.77), hence the conversion factor remains as 2.52 Â 10 12 e À cm À2 V À1 .
www.advancedsciencenews.com www.advintellsyst.com 0.36 m 0 . This is likely due to the absence of trap-assisted tunneling (TAT) in the TCAD simulations that becomes prominent at smaller V. [19] The nonlinearity parameter β for the experimental SONOS device is ten times lower, so that the low-bias dynamics is faster compared with the simulation model. This may also be caused by the absence of TAT through defect states in the forbidden bands of SiO 2 and Si 3 N 4 and thus worse predictability for low-bias cases. This is a known problem for industrial SONOS TCAD models, where the pass disturb simulations with V ¼ 7-8 V typically underestimate Q M compared with fabricated devices, [19] because defects in SiO 2 and Si 3 N 4 are difficult to model. The parameter γ, responsible for the sensitivity to Q M , was found to be the same for both the physics simulations and the experimental data. Finally, δ, which handles the Q M sensitivity of M depending on V, was found to be unity, that is, there was no noticeable dependence of Q M on V from the experimental SONOS devices. Figure 7d shows the improved fits with distinct parameters for potentiation (V < 0) and depression (V > 0), such that little error remains between the extracted memristance values and the analytic equation, as also found in two-terminal compact models. [46] We deployed the compact model in a commercial circuit simulation package, Cadence Spectre, using VerilogA to assess the agreement with the measured drain current (with v D ¼ 0.1 V and v G ¼ 2.4 V during reading) [33] after various training pulses (v G ) of À12, À11, À10, À9, þ9, þ10, þ11 V, as shown in Figure 8a. Despite the remarkable simplicity of the compact model, good  [33] under v G ¼ 2.4 V and v D ¼ 0.1 V for reading the channel current. While the measurements of i D for the experimental SONOS devices were conducted after every writing pulse (À9 to À12 V for potentiation, þ9 to þ11 V for depression on v G with 10 μs width while v D ¼ v S ¼ 0), our Cadence Spectre simulations were sampled after every 100 μs. Slight errors are attributed to the deviations from the perfectly exponential increase of M with Q M assumed in our compact model. b) The corresponding evolution of the state variable Q M is available from the simulation, which demonstrates the strong correlation between ΔQ M and Δi D . . The analytical function with a parameter set (α, β, γ, δ) ¼ (70 000, 7, 50, 1) simultaneously fits both potentiation and depression, but with significant errors possibly due to different physics for charge trapping and detrapping. d) A modified function for M, where the parameters are different for potentiation (α À , β À , γ À , δ À ) ¼ (10 6 , 2.5, 3000, 1) and depression (α þ , β þ , γ þ , δ þ ) ¼ (5500, 9, 900, 1), models the experimental SONOS data with negligible error.
www.advancedsciencenews.com www.advintellsyst.com agreement between the simulation and the experiment was confirmed with the accuracy higher than 90%, proving that our compact model captures the essential physics in a simple way. Figure 8b shows the corresponding evolution of the key state variable Q M , where a linear correlation is observed between ΔQ M and Δi D as foreseen by the description of R CH in Equation (6). The linearity between ΔQ M and Δi D can be further exploited for linear and symmetric writing, [50] which is ideally desired for autonomous weight updates without the necessity for a time-consuming "write and verify" scheme in neuromorphic learning. Fuller et al. [32] reported the use of current pulses instead of voltage pulses to mitigate the nonlinearity and asymmetry in writing that is commonly observed in memristors or NVM devices including SONOS devices, as shown in Figure 3a,b, 5a,c, and 8a,b. The mathematical representation of our compact model as represented by Equation (6) is consistent with the empirical observation of Fuller et al. [32] and thus three-terminal synaptic circuit elements are expected to be promising candidates for the building unit of neuromorphic learning systems.
In conclusion, this work presented a well-posed compact model of a SONOS three-terminal synaptic circuit element that can be readily utilized by circuit designers for neuromorphic computing circuits/systems. This compact model offers the rare combination of simplicity and good predictability, stemming from physics-driven state variables. The basis of the compact model is rooted on the carrier concentration modulation in the conducting channel, so that it is readily applicable to other three-terminal synaptic circuit elements besides SONOS devices.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.