0 Introduction

When the digital circuits are operated near the threshold/moderate inversion region, they are highly sensitive to process variations. Intrinsic and extrinsic factors contribute to process variation. Intrinsic differences are often caused by the limits of the device physics, while extrinsic variations are produced by wafer-to-wafer and die-to-die variations. Although extrinsic differences may be minimized via improved process control, intrinsic variations are uncontrollable and serve as a significant barrier to device scalability. Significant causes of intrinsic variation include the following^[1-4]:

1) Random dopant fluctuation: It has been shown in the literature that the quantity and location of dopant atoms implanted are random in nature at a scaled technology node. Due to the extremely small number of dopant atoms in scaled devices, there is a significant fluctuation in the threshold voltage (which is a function of dopant concentration). When the dopant atoms are big, the Poisson distribution may be used to estimate the dopant fluctuations, as opposed to the Gaussian distribution^[5-6].

2) Roughness of the line edge: This refers to the variance in the form of the gate along the channel width direction. It happens more often near the mask borders, where bigger material grains dissolve more readily than smaller material grains. The roughness of the line edge does not scale with technology and cannot be decreased via lithographic process control. However, this has an effect on the threshold voltage and leakage current of the device^[7-8].

3) Oxide thickness fluctuations: When the oxide layer thickness is a few Si atom layers, the atomic scale roughness at the interface of Si - SiO₂ and gate SiO₂ produces oxide thickness variations in scaled devices. This has an effect on the carrier mobility and thus on the device's threshold voltage^[9-10].

All of the above-mentioned fluctuations become more significant as the technology node scales down, but random dopant fluctuations (which result in threshold variation) become increasingly crucial as the supply voltage drops. It is because the variation in the drain current is asymmetric at the ultradeep submicron channel length and at lower supply voltage. Statistical metrics like mean and standard deviation cannot be used to identify the shape of the distribution^[9-10].

Thus, the probability distribution of drain current has been studied using Skewness, a measure of asymmetry. Fig. 1 illustrates the skewness of a CMOS inverter's current distribution with regard to the supply voltage. The findings were acquired using a Monte Carlo simulation with 20000 HSPICE nodes at the 32 nm technology node. As the supply voltage drops, the drain current varies asymmetrically owing to the change in threshold voltage, which is considerably greater than at higher supply voltage.

In general, the variability of any random variable may be expressed as its standard deviation normalized by its mean (σ/μ) value, which indicates its global deviation^[11]. This leads to a statistical distribution of performance in comparison with the circuits operating in the strong inversion region. Low voltage operation (in the moderate inversion region) exhibits non-Gaussian delay distribution, enabling utilization of advanced statistical models^[12].

Statistical based solutions are classified as either fitting-based or analytical-based. Examples of fitting based non-Gaussian solutions are based on Monte Carlo simulations such as log skew normal distribution, inverse Gaussian distribution, and Weinbull distribution. These distributions have a long simulation time and do not provide any physical insight into the design parameters^[13-14]. Many analytical models have been proposed in literature to study the impact of process variations on circuit performance^[13-14]. In Ref. [15], log normal distribution for gate delay variation is applied to circuits operating in the subthreshold region. It has been assumed that the threshold voltage follows a Gaussian distribution under process variation. The analysis has been extended to stacked transistors and applied to the NAND gate as an example. In Ref.[16], the gate delay variability (σ/μ) model for subthreshold circuits has been presented analytically with the condition that the correlation between the stacked transistors is derived from generic logic gates under specific corner conditions.

The transregional MOSFET model with closed form equations was developed by Ref. [17]. The model is utilized to generate a log normal delay distribution under variation in threshold voltage by Ref. [18]. Further improvement of the model was given by Ref. [19] that described the micro-architectural statistical timing performance. All the above models followed log normal distribution for delay under threshold voltage variation for low supply voltages. But, none of the models was specific to the moderate inversion region. The following are the major contributions:

1) The unified transregional MOSFET model applicable in the entire range of moderate inversion has been presented to find the delay model for ultradeep submicron technology up to 22 nm.

2) It has been shown that the TRM model^[17] verification with BSIM4^[20] shows lower errors at the edge of weak inversion, which increases as the transistor moves towards the strong inversion region. But the proposed model is applicable for the entire range of moderate inversions.

3) In order to study variability at the 22 nm technology node, an analytical delay model under threshold voltage (dominant parameter) variation was developed using a log skew normal distribution.

4) The moment matching function has been utilized to calculate an analytical log skew normal based delay model. The distribution parameters are represented as a function of threshold voltage, supply voltage and the number of transistors that are stacked.

5) The proposed model is validated in the entire moderate inversion range using a CMOS inverter as a sample gate.

6) The proposed methodology has been compared with log normal and Monte Carlo simulations done using predictive technology model 22 nm (PTM)^[21]. The obtained results were better than that with both the techniques.

7) The proposed model is capable of finding the delay variability, maximum and minimum delay for CMOS inverters at 22 nm.

1 CMOS Inverter Delay Derivation

The velocity versus electric field for electrons in the inversion layer is described as in Ref. [22] where U_sat is saturation velocity (m/s), E represents a lateral electric field, while E_c represents the critical field at which the carriers become velocity saturated:

$ v(E)=\frac{\mu_{\mathrm{eff}} E}{\left[1+\left(\frac{E}{E_{\mathrm{c}}}\right)^M\right]^{\frac{1}{M}}} $

(1)

where

$ E_{\mathrm{c}}=\frac{v_{\text {sat }}}{\mu_{\text {eff }}} $

Here, v_sat is the saturation velocity, μ_eff is the effective mobility and M is constant. It is commonly acknowledged in literature that M=2 with electrons and while M=1 for holes. Because the mathematics required to calculate the drain current is difficult and time-consuming, the modified first degree U-E model inserting M=1 with fitting constants a, b, and c is presented below, where $ E_{\mathrm{c}}=\frac{b+c}{a} \cdot \frac{v_{\mathrm{sat}}}{\mu_z}, \mu_z$ is the mobility dependent on the vertical applied field. a, b, c are the temperature and channel length dependent parameters as given in Table 1.

$ v\left(E_x\right)=\frac{a \mu_z E_x}{\left[b+c\left(\frac{E_x}{E_{\mathrm{c}}}\right)\right]} $

(2)

and

$\mu_z=\frac{\mu_0}{1+\theta\left(V_{\mathrm{gs}}-V_{\mathrm{th}}\right)} $

(3)

where θ represents coefficient of mobility degradation, E_x is the lateral electric field, μ_z represents the degradation effect due to vertical field, and μ₀ denotes mobility due to low electric field; V_gs is the bias voltage between gate and source; V_th represents the threshold voltage. The relationship between μ and E can be expressed as^[23-24]

$ \mu_{\mathrm{eff}}\left(E_x\right)=\frac{a \mu_z}{\left[b+c\left(\frac{E_x}{E_{\mathrm{c}}}\right)\right]} $

(4)

In the presence of parasitic resistance, the mobility degradation coefficient can be expressed as that in Ref. [25]:

$ \theta_m=\theta+\beta_0 R_{\mathrm{s}} $

(5)

where β₀ denotes the maximum value of device transconductance, θ represents coefficient of mobility degradation, and R_s represents the value of device series resistance.

PTM^{[21, 26]} describes the effect of velocity overshoot as

$ V_{\mathrm{sat}}=V_{\mathrm{sato}}+0.13 \mu_{\mathrm{eff}} \sqrt{\tau \mu_{\mathrm{eff}} \frac{k T}{q}}\left(\frac{V_{\mathrm{ds}}}{L^2}\right) $

(6)

where V_sato represents the carrier velocity saturation without overshoot; τ denotes the transition time of the charge carriers to travel from source to drain; V_ds represents the applied bias potential from drain to source; L is the effective channel length, k is the Boltzmann constant, T is the temperature and q is the unit charge.

Normalizing every quantity to make it unitless^[27]:

$ q_{\mathrm{i}}=\frac{-Q_{\mathrm{i}}}{2 n \phi_{\mathrm{t}} C_{\mathrm{ox}}} $

(7)

where q_i indicates the normalized value of the inversion charge Q_i, n is the slope factor, C_ox is the oxide capacitance and ϕ_t is the thermal voltage.

$ I_0=2 n \mu \frac{W}{L} \phi_{\mathrm{t}}^2 C_{\text {ox }} $

(8)

where I₀ indicates technological current, W is the channel width and μ is the carrier mobility.

$ v_{\mathrm{c}}=\frac{V_{\mathrm{c}}}{\phi_{\mathrm{t}}} $

(9)

where v_c indicates the normalized value of channel potential, V_c.

$ i(\text { drain current })=\frac{I}{I_0} $

(10)

where i indicates the normalized value of drain current I.

$ u=\frac{\mu_{\text {eff }}\left(E_{x}\right)}{\mu_{z}} $

(11)

where u indicates mobility.

$ e=\frac{\left|E_{x}\right|}{E_{\mathrm{c}}}=-\lambda_{\mathrm{c}} \frac{\mathrm{d} q_{\mathrm{i}}}{\mathrm{d} \xi} $

(12)

where e represents horizontal electric field, λ_c is the coefficient velocity saturation and is a function of effective channel length and ξ is the normalized position along x-axis. $ \lambda_{\mathrm{c}}=\frac{2 \phi_{\mathrm{t}}}{E_{\mathrm{c}} L}$. The expression of drain current as a function of vertical and horizontal electric field is expressed as^[27]

$ i_{\mathrm{ds}}=-\left(2 q_{\mathrm{i}}+1\right) u \frac{\mathrm{d} q_{\mathrm{i}}}{\mathrm{d} \xi} $

(13)

Putting u from Eq. (11) in Eq. (13) and rearranging,

$i_{\mathrm{ds}} b \mathrm{~d} \xi=-\left(a\left(2 q_{\mathrm{i}}+1\right)-c \lambda_{\mathrm{c}} i_{\mathrm{ds}}\right) \mathrm{d} q_{\mathrm{i}} $

(14)

where i_ds is the device drain to source current.

Integrating Eq. (14) from ξ=0 and q_i=q_s to ξ=1 and q_i=q_d,

$ i_{\mathrm{ds}}=\frac{a\left(q_{\mathrm{s}}^{2}-q_{\mathrm{d}}^{2}+q_{\mathrm{s}}-q_{\mathrm{d}}\right)}{b+c \lambda_{\mathrm{c}}\left(q_{\mathrm{s}}-q_{\mathrm{d}}\right)} $

(15)

where q_s indicates normalized inversion charge at the source end, and q_d indicates inversion charge at the drain end.

In order to further simplify Eq. (15), assume a=b=c=1

$ i_{\mathrm{ds}}=\frac{q_{\mathrm{s}}^{2}-q_{\mathrm{d}}^{2}+q_{\mathrm{s}}-q_{\mathrm{d}}}{1+\lambda_{\mathrm{c}}\left(q_{\mathrm{s}}-q_{\mathrm{d}}\right)} $

(16)

The inversion charge density in the presence of velocity saturation, i.e., when q_i=q_dsat^[27],

$ i_{\mathrm{dsat}}=\frac{2 q_{\mathrm{dsat}}}{\lambda_{\mathrm{c}}} $

(17)

$v_{\text {dsat }}=v_{\mathrm{p}}-2 q_{\text {dsat }} $

(18)

$ \begin{gathered} q_{\mathrm{d}}= \begin{cases} \frac{v_{\mathrm{p}}-v_{\mathrm{d}}}{2}, & \text { if } v_{\mathrm{d}} <v_{\mathrm{dsat}} \\ q_{\mathrm{dsat}}, & \text { if } v_{\mathrm{d}} \geqslant v_{\mathrm{dsat}}\end{cases} \end{gathered} $

(19)

where q_dsat, v_dsat, v_p, v_d, are the normalized values of saturated drain charge, saturated drain voltage, pinch-off voltage and drain voltage.

For q_dsat, utilizing Eq.(16) and Eq.(17) at i_ds=i_dsat,

$ \frac{2 q_{\mathrm{dsat}}}{\lambda_{\mathrm{c}}} =\frac{q_{\mathrm{s}}^{2}-q_{\mathrm{dsat}}^{2}+q_{\mathrm{s}}-q_{\mathrm{dsat}}}{1+\lambda_{\mathrm{c}}\left(q_{\mathrm{s}}-q_{\mathrm{dsat}}\right)} $

(20)

$ \begin{aligned} q_{\mathrm{dsat}} & =\frac{1}{2 \lambda_{\mathrm{c}}}\left[\lambda_{\mathrm{c}}\left(1+2 q_{\mathrm{s}}\right)-\right. \\ & \left.\sqrt{\lambda_{\mathrm{c}}^{2}+4 \lambda_{\mathrm{c}}\left(1+2 q_{\mathrm{s}}\right)+4}+2\right] \end{aligned} $

(21)

$ \begin{aligned} i_{\mathrm{dsat}}= & \frac{2}{\lambda_{\mathrm{c}}^{2}}\left[1+\frac{1}{2} \lambda_{\mathrm{c}}\left(1+2 q_{\mathrm{s}}\right)-\right. \\ & \left.\frac{1}{2} \sqrt{\lambda_{\mathrm{c}}^{2}+4 \lambda_{\mathrm{c}}\left(1+2 q_{\mathrm{s}}\right)+4}\right] \end{aligned} $

(22)

The relationship of inversion charge density and channel potential can be expressed following the Ref. [17] and expressed in Eq. (23) where k_a, k₂ and k₃ denotes fitting coefficients; v_y=v_p-v_c, v_p represents normalized value of pinch off potential, v_c denotes the voltage across the channel and v_y is the voltage across channel y-axis, $V_{\mathrm{p}}=\frac{V_{\mathrm{g}}-V_{\mathrm{th}}}{n \phi_{\mathrm{t}}} $ :

$ q_i \approx k_a \mathrm{e}^{k_2 v_y+k_3 v_y^2} $

(23)

where V_p is the unnormalized pinch-off voltage, V_g is the unnormalized gate voltage, V_th is the unnormalized threshold voltage, and k_a is the fitting coefficient.

Unnormalize Eq. (22), and substitute q_i from Eq. (23). Consider source reference and deep saturation, considering

$A=k_a \mathrm{e}^{k_2 \frac{\mathrm{VG}}{n \phi_{\mathrm{t}}}+k_3\left(\frac{\mathrm{VG}}{n \phi_{\mathrm{t}}}\right)^2}, \mathrm{VG}=V_{\mathrm{gs}}-V_{\mathrm{th}} $

(24)

$\begin{aligned} & I_{\mathrm{dsat}}=2 n \mu_z \frac{W}{L} \phi_t^2 C_{\mathrm{ox}}\left[\frac{1}{4}\left(\frac{E_{\mathrm{c}} L}{\phi_{\mathrm{t}}}\right)^2 .\right. \\ & \left.\left(\sqrt{\left(\frac{\phi_{\mathrm{t}}}{E_{\mathrm{c}} L}\right)^2+\frac{2 \phi_{\mathrm{t}}}{E_{\mathrm{c}} L}(1+2 A)+1}-1\right)^2-\frac{1}{4}\right] \end{aligned} $

(25)

where VG is the overdrive voltage, V_gs is the gate to source voltage, A is the multiplying factor and I_dsat is the drain saturation current.

Eq. (25) corresponds to drain current for short channel devices considering both the valid drift and diffusion component in all three operating regions. For the moderate inversion region, when the field is low, expand the expression in Eq. (25) through binomial expansion and take initial terms:

$ \begin{gathered} I_{\mathrm{dsat}}=n \mu \frac{W}{L} \phi_{\mathrm{t}}^{2} C_{\mathrm{ox}} \frac{1}{2}\left[\left(\frac{\phi_{\mathrm{t}}}{4 E_{\mathrm{c}} L}\right)^{2}+(1+2 A)^{2}+\right. \\ \left.\left(\frac{\phi_{\mathrm{t}}}{2 E_{\mathrm{c}} L}\right)(1+2 A)-1\right] \end{gathered} $

(26)

$ I_{\mathrm{dsat}}=2 n \mu \frac{W}{L} \phi_{\mathrm{t}}^{2} C_{\mathrm{ox}}\left[A^{2}+A\right] $

(27)

Substituting A from Eq. (24),

$ \begin{gathered} I_{\mathrm{dsat}}=2 n \mu \frac{W}{L} \phi_{\mathrm{t}}^{2} C_{\mathrm{ox}}\left[\left(k_{a} e^{k_{2} \frac{\mathrm{VG}}{n \phi_{\mathrm{t}}}+k_{3}\left(\frac{\mathrm{VG}}{n \phi_{\mathrm{t}}}\right)^{2}}\right)^{2}+\right. \\ \left.k_{a} \mathrm{e}^{k_{2} \frac{\mathrm{VG}}{n \phi_{\mathrm{t}}}+k_{3}\left(\frac{\mathrm{VG}}{n \phi_{\mathrm{t}}}\right)^{2}}\right] \end{gathered} $

(28)

At low supply voltage, the first exponent term is neglected:

I_ON is the current that is described as drain current when V_gs=V_ds=supply voltage V_dd. Therefore, in above equation, defining technology current I_S=2nμ ϕ_t²C_ox.

$ I_{\mathrm{ON}}=I_{S} \frac{W}{L} K_{a} e^{Z} $

(29)

$ Z=\mathrm{e}^{k_{2} \frac{\mathrm{VG}}{\phi_{\mathrm{t}}}}+k_{3}\left(\frac{\mathrm{VG}}{n \phi_{\mathrm{t}}}\right)^{2} $

(30)

where Z is the constant factor.

The RC delay of CMOS inverter can be expressed as t_d as described in Eq. (31) where K_d is the fitting parameter, and C_L denotes the load capacitance. The average value of K_d shall be utilized to reduce the fitting error:

$ t_{\mathrm{d}}=K_{\mathrm{d}} \frac{V_{\mathrm{dd}} C_{L}}{I_{\mathrm{ON}}} $

(31)

The logarithmic expression of Eq. (31) can be written as

$ Y=\ln \left(t_{\mathrm{d}}\right)=-Z+C_{t_{\mathrm{d}}} $

(32)

The delay of the gate is assumed to follow a log skew normal distribution. Therefore, its logarithmic form expressed in Eq. (32) follows a skew normal distribution whose probability density function can be written as^[28]

$f_{\mathrm{SN}}(Y)=\frac{2}{\omega} \phi\left(\frac{Y-\epsilon}{\omega}\right) \varPhi\left(\lambda \frac{Y-\epsilon}{\omega}\right) $

(33)

Here,

$ \epsilon=\mu-\omega \beta \sqrt{\frac{2}{\pi}} $

(34)

$ \omega=\sqrt{\frac{\sigma^{2}}{1-\frac{2}{\pi} \beta^{2}}} $

(35)

$ \beta=\frac{\pi}{2} \frac{\left(\frac{2 \gamma_{1}}{4-\pi}\right)^{\frac{2}{3}}}{1+\left(\frac{2 \gamma_{1}}{4-\pi}\right)^{\frac{2}{3}}} $

(36)

where $\epsilon $ is location parameter, ω is scale parameter, β is shape parameter, ϕ(.) denotes the PDF and Φ(.) denotes the CDF of the distribution of type standard normal (SN), and γ₁ denotes the skewness. The distribution parameters $\epsilon $, ω and λ are represented in terms of mean denoted by μ, variance denoted by σ² and skewness denoted by β. The relationship of β and λ can be represented as

$ \lambda=\frac{\beta}{\sqrt{1-\beta^{2}}} $

(37)

From Eq. (30) and Eq. (32), it can be seen that Y is dependent on threshold voltage. Y can be calculated from the mean (first moment), variance (second moment) and standard deviation (third moment) assuming a Gaussian distribution of threshold voltage. The following are the expressions:

$ \mu=-\frac{k_{2}}{n \phi_{\mathrm{t}}} E\left[V_{\mathrm{d} t}\right]-\frac{k_{3}}{\left(n \phi_{\mathrm{t}}\right)^{2}} E\left[V_{\mathrm{d} t}^{2}\right]+C_{t_{\mathrm{d}}} $

(38)

$\sigma^{2}=\frac{k_{3}^{2}}{\left(n \phi_{\mathrm{t}}\right)^{4}}\left(E\left[V_{\mathrm{d} t^{\prime}}^{4}\right]-E\left[V_{\mathrm{d} t^{\prime}}^{2}\right]^{2}\right) $

(39)

$ \begin{gathered} \gamma_{1}=\frac{-\frac{k_{3}^{3}}{\left(n \phi_{\mathrm{t}}\right)^{6}}\left(E\left[V_{\mathrm{d} t^{\prime}}^{6}\right]+\frac{3}{4} \frac{k_{2}^{2} k_{3}}{\left(n \phi_{\mathrm{t}}\right)^{4}} E\left[V_{\mathrm{d} t^{\prime}}^{4}\right]\right)}{\sigma^{3}}+ \\ \frac{\frac{3}{16} \frac{k_{2}^{4}}{\left(n \phi_{\mathrm{t}}\right)^{2} k_{3}} E\left[V_{\mathrm{d} t^{\prime}}^{2}\right]}{\sigma^{3}}+\frac{\frac{k_{2}^{6}}{64 k_{3}^{3}}-3 \mu \sigma^{2}-\mu^{3}}{\sigma^{3}} \end{gathered} $

(40)

V_dt and V_dt' are the expressions that follow Gaussian distribution and expressed in Eq. (41) and Eq. (42), $C_{t_{\mathrm{d}}}=\ln \left(K_{\mathrm{d}} \frac{V_{\mathrm{dd}} C_{L}}{I_{S} \frac{W}{L} k_{a}}\right) $ and E[.] calculate the mean value in Eq. (40). Also, μ₀ and σ₀ represents mean and standard deviation respectively. Thus,

$ V_{\mathrm{d} t}=V_{\mathrm{dd}}-V_{\mathrm{th}} \sim N\left(V_{\mathrm{dd}}-\mu_{0}, \sigma_{0}\right) $

(41)

$ \begin{aligned} V_{\mathrm{d} t^{\prime}}=V_{\mathrm{dd}}-V_{\mathrm{th}}+\frac{n \phi_{\mathrm{t}} k_{2}}{2 k_{3}} \sim N\left(V_{\mathrm{dd}}-\mu_{0}+\right. \\ \left.\quad \frac{n \phi_{\mathrm{t}} k_{2}}{2 k_{2}}, \sigma_{0}\right) \end{aligned} $

(42)

where V_dt and V_dt^′ are the expressions that follow Gaussian distribution. Further, in Eq. (38), Eq. (39) and Eq. (40), several terms can be described as

$ E\left[V_{\mathrm{d} t^{\prime}}^{2}\right]=\sigma_{0}^{2}+\left(E\left[V_{\mathrm{d} t^{\prime}}\right]\right)^{2} $

(43)

$ E\left[V_{\mathrm{d} t^{\prime}}^{4}\right]=3 \sigma_{0}^{4}+6 \sigma_{0}^{2}\left(E\left[V_{\mathrm{d} t^{\prime \prime}}\right]\right)^{2}+ \\ \left(E\left[V_{\mathrm{d} t^{\prime \prime}}\right]\right)^{4} $

(44)

$ E\left[V_{\mathrm{d} t^{\prime}}^{6}\right]=15 \sigma_{0}^{6}+45 \sigma_{0}^{4}\left(E\left[V_{\mathrm{d} t^{\prime}}\right]\right)^{2}+ \\ 15 \sigma_{0}^{2}\left(E\left[V_{\mathrm{d} t^{\prime \prime}}\right]\right)^{4}+\left(E\left[V_{\mathrm{d} t^{\prime \prime}}\right]\right)^{6} $

(45)

where

$ E\left[V_{\mathrm{d} t}\right]=V_{\mathrm{dd}}-\mu_{0} $

(46)

$ E\left[V_{\mathrm{d} t^{\prime}}\right]=V_{\mathrm{dd}}-\mu_{0}+\frac{n \phi_{\mathrm{t}} k_{2}}{2 k_{3}} $

(47)

$ E\left[V_{\mathrm{d} t^{\prime \prime}}\right]=V_{\mathrm{dd}}-\mu_{0}+\frac{\frac{k_{2}}{n \phi_{\mathrm{t}}}}{\frac{k_{3}}{\left(n \phi_{\mathrm{t}}\right)^{2}}} $

(48)

Eqs.(38)-(40) show that most of the parameters can be obtained from fitting, whether they are process dependent or process independent, except the mean and standard deviation of the device threshold voltage. The mean value of the Gaussian-distributed threshold voltage can be utilized as the nominal value for various device sizes in SPICE simulations. Monte Carlo simulations are performed to find the standard deviation for the smallest transistor at each process node. For long channel transistors, Pelgrom's Law which states that the standard deviation is inversely to find the standard deviation^[29]. The amount of effort required for Monte Carlo simulations is reduced to one at different values of the transistor size and the load capacitance because of the model explained in Eqs.(38)-(40). By combining Eqs.(34)-(36) with Eqs.(38)-(40), the distribution parameters $\epsilon $, ω and β can be found.

Now, T_d follows log skew normal distribution with distribution parameters the same as skew normal distributed Y^[30]. Thus,

$ T_{\mathrm{d}} \sim \operatorname{LSN}\left(\epsilon, \omega^{2}, \lambda\right) $

(49)

Its PDF and CDF can be expressed as

$ f_{\mathrm{LSN}}\left(t_{\mathrm{d}}\right) =\frac{2}{\omega y} \phi\left(\frac{\ln \left(t_{\mathrm{d}}\right)-\epsilon}{\omega}\right) \varPhi\left(\lambda \frac{\ln \left(t_{\mathrm{d}}\right)-\epsilon}{\omega}\right) $

(50)

$\begin{aligned} F_{\mathrm{LSN}}\left(t_{\mathrm{d}}\right) =\varPhi\left(\frac{\ln \left(t_{\mathrm{d}}\right)-\epsilon}{\omega}\right)-2 T\left(\lambda \frac{\ln \left(t_{\mathrm{d}}\right)-\epsilon}{\omega}, \lambda\right) \end{aligned} $

(51)

Here, Owen's function T can be expressed as

$ T\left(\frac{\ln \left(t_{\mathrm{d}}\right)-\epsilon}{\omega}, \lambda\right)=\frac{1}{2 \pi} \int_{0}^{\lambda} \cdot \frac{\mathrm{e}^{\frac{-\ln \left(t_{\mathrm{d}}\right)-\epsilon^{2}}{\omega}}\left(1+t^{2}\right)}{1+t^{2}} \mathrm{~d} t $

(52)

where t indicates real variable (the frequency parameter).

The mean and variance of the gate delay can be expressed as (after finding PDF and CDF)

$ E_{t_{\mathrm{d}}}=2 \mathrm{e}^{\epsilon} \mathrm{e}^{\frac{\omega^{2}}{2}} \phi(\beta \omega) $

(53)

$ D_{t_{\mathrm{d}}}=2 \mathrm{e}^{2 \epsilon} \mathrm{e}^{\omega^{2}}\left(\mathrm{e}^{\omega^{2}} \phi(2 \beta \omega)-2 \phi^{2}(\beta \omega)\right) $

(54)

2 Analysis of Worst and Best Case Delay

For the log skew normal distributed model for delay, the points of ±3σ percentile can be obtained by setting the CDF given in Eq. (51) to be equal to Φ(±3) represented as

$ \begin{array}{r} F_{\mathrm{LSN}}(Y)=\varPhi\left(\frac{\ln (Y)-Y}{\omega_{Y}}\right)-2 T\left(\frac{\ln (Y)-\epsilon_{Y}}{\omega_{Y}}, \lambda_{Y}\right)= \\ \varPhi( \pm 3) \end{array} $

(55)

where λ_Y denotes the shape parameter.

Now, in order to solve the above equation, the shape distribution parameter of delay t_d in log skew normal distribution is approximately taken to a value close to 1. Using the Owen's function,

$ T(h, 1)=\frac{1}{2} \varPhi(h)(1-\varPhi(h)) $

(56)

where h is a numeric vector.

The above equation can be derived from

$ \varPhi^{2}\left(\frac{\ln (Y)-\epsilon_{Y}}{\omega_{Y}}\right)=\varPhi( \pm 3) $

(57)

The maximum and minimum 3σ or -3σ percentile points are located at

$ t_{\mathrm{dmax}}=\mathrm{e}^{\epsilon_{Y^{+}} \varPhi^{-1}(\sqrt{\varPhi(3)}) \omega_{Y}} \approx \mathrm{e}^{\epsilon_{Y}+3.21 \omega_{Y}} $

(58)

$ t_{\mathrm{d} \min }=\mathrm{e}^{\epsilon_{Y^{+}} \varPhi^{-1}(\sqrt{\varPhi(3)}) \omega_{Y}} \approx \mathrm{e}^{\epsilon_{Y}-1.79 \omega_{Y}} $

(59)

Putting the values of $\epsilon_{t_{\mathrm{d}}} $ and ω_td derived in Eq.(34) and Eq.(35), where $ \epsilon_{t_{\mathrm{d}}}$ and ω_td represents location and scale parameters respectively, the minimum and maximum value can be calculated.

3 Results and Discussion

First of all, the drain current of the proposed model is verified and compared with the industry standard BSIM4 and TRM models^[17] for the entire range of the moderate inversion as shown in Fig. 2. It can be seen from the figure that the TRM model deviates at the end of moderate inversion or when the transistor moves towards the strong inversion region (The percentage error increases up to 20%). For the 22 nm technology node, the values of k_a, k₁ and k₂ are 0.532, 0.701 and 0.0037 respectively. The probability density function of the proposed model is compared with the results obtained from Fig. 3 which shows the PDF of the CMOS inverter gate with FO4 load. It can be seen from the figure that the proposed model has good agreement with Monte Carlo simulations. In order to evaluate the accuracy of the proposed model, the delay variability, maximum and minimum values of delay have been considered as the performance metrics. The delay variability is calculated by normalizing the standard deviation of the gate delay by its mean value. The maximum and minimum points are generally calculated to approximate the performance in extreme conditions. Figs. 4-6 show the delay variability, maximum and minimum values of the CMOS inverter with FO4 load at three different temperatures.

4 Conclusions

In this paper, an accurate mathematical log skew normal distributed gate delay model is presented for transistors operating specifically in the moderate inversion region. For worst-case analysis, the cumulative distribution function is used in conjunction with the position scale and shape parameter to determine the maximum and minimum value of delay. The proposed model has excellent agreement with Monte Carlo simulation and is better than the existing statistical delay models when compared variability in delay.