*Article* **Reliability Evaluation Based on Mathematical Degradation Model for Vacuum Packaged MEMS Sensor**

**Guizhen Du 1,2, Xianshan Dong 2,\*, Xinglong Huang 2,3, Wei Su <sup>2</sup> and Peng Zhang 1,\***


**Abstract:** Vacuum packaging is used extensively in MEMS sensors for improving performance. However, the vacuum in the MEMS chamber gradually degenerates over time, which adversely affects the long-term performance of the MEMS sensor. A mathematical model for vacuum degradation is presented in this article for evaluating the degradation of vacuum packaged MEMS sensors, and a temperature-accelerated test of MEMS gyroscope with different vacuums is performed. A mathematical degradation model is developed to fit the parameters of the degradation of Q-factor over time at three different temperatures. The results indicate that the outgassing rate at 85 ◦C is the highest, which is 0.0531 cm2/s; the outgassing rate at 105 ◦C is the lowest, which is 0.0109 cm2/s; and the outgassing rate at 125 ◦C is in the middle, which is 0.0373 cm2/s. Due to the different mechanisms by which gas was released, the rate of degradation did not follow this rule. It will also be possible to predict the long-term reliability of vacuum packaged MEMS sensors at room temperature based on this model.

**Keywords:** MEMS sensor; reliability evaluation; vacuum degradation; mathematical model

## **1. Introduction**

The micro-electro-mechanical system (MEMS) is a micro-intelligent system utilizing microelectronics, micromechanics, and other technologies to integrate sensors, actuators, and signal transmission units [1]. Micro-electro-mechanical sensors possess the characteristics of small volume, light weight, and low cost, and are on the verge of achieving passivity, miniaturization, and anti-interference capabilities [2]. As technology advances, MEMS sensors have become widely used in a variety of fields, including consumer electronics, aerospace, military equipment, biomedicine, and others [3–6]. It is common for MEMS sensors to have movable structures, such as gyroscopes, accelerometers, and filters that are based on resonant structures. Vacuum packaging can reduce the air damping of their movable structures, improving their performance. In spite of this, the vacuum in the MEMS chamber will gradually degenerate over time, affecting directly the performance of vacuum packaged MEMS sensors, while seriously affecting their long-term reliability. Vacuum degradation and failure are common problems associated with vacuum-packaged MEMS sensors. As a consequence, the study of the vacuum degradation and long-term reliability of vacuum packaged MEMS sensor is of great importance for the development and application of vacuum packaged MEMS sensors.

In the present day, some scholars have conducted in-depth research into the vacuum degradation of MEMS devices. Through vacuum reflow technology, Liang et al. [7] were able to ensure good sealing of the packaging cavity and placed getters within the cavity to maintain a stable vacuum inside the cavity. As a result of testing the resonant gyroscope

**Citation:** Du, G.; Dong, X.; Huang, X.; Su, W.; Zhang, P. Reliability Evaluation Based on Mathematical Degradation Model for Vacuum Packaged MEMS Sensor. *Micromachines* **2022**, *13*, 1713. https://doi.org/10.3390/mi13101713

Academic Editor: Wensheng Zhao

Received: 8 September 2022 Accepted: 3 October 2022 Published: 11 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

for nine months, Wu et al. [8] prepared multilayer getters, which were activated at 300 ◦C, and the getters showed effective performance. As reported by Tepolt et al. [9], the initial outgassing of the material in the packaging cavity, as well as the vacuum pressure required to stabilize the devices, was measured. In order to optimize the processing technology and minimize the exhaust amount of the MEMS devices cavity without serious degradation, the data was compared with the specific getter capacity data. A mathematical model has been developed by Lu et al. [10] that provides a theoretical basis for directly comparing the behavior of different gases entering the sealing cavity. According to the principles of thermodynamics and hydrodynamics, He et al. [11] determined the relationship between Q-factor, air pressure, and gas number of MEMS sensor cavities. For MEMS resonators, Candler proposed a single-chip packaging method [12]. Using thick polysilicon packaging, the MEMS resonator was packaged under a pressure of less than 1 Pa, and the reliability test was conducted under high temperatures. As a result, the majority of research focuses on improving the packaging process and developing packaging materials for MEMS devices. Some researchers devote a great deal of time to assessing the long-term reliability of vacuum-packaged sensors. Currently, there are few mature theoretical models of the effect of deflation of the internal materials of MEMS devices on the long-term reliability of these devices.

As part of this paper, an evaluation of the reliability of vacuum-packaged MEMS sensors is presented. In order to analyze the degradation failure of vacuum packaged MEMS sensors, a degradation failure model was developed. Additionally, the parameters of the model were obtained by fitting data obtained from MEMS sensors with varying vacuum degrees at multiple temperature points. According to the results, the model is well suited to the experimental data. Further, the model was used to predict the vacuum degradation of a vacuum-packaged MEMS sensor at room temperature based on the model obtained.

#### **2. Materials and Methods**

#### *2.1. Model of Vacuum Degradation*

In vacuum packaged MEMS sensors, vacuum degree is degraded mainly due to the release of gas inside the cavity [13], including gas molecules adsorbed on the inner wall of the cavity, gas molecules in the material inside the cavity, etc. High temperatures are the main source of sensitive stress [14]. It is generally true that vacuum-packaged MEMS sensors have resonant structures, which can be used as a measure of the level of vacuum inside the cavity of the package [15]. Hence, this paper first establishes a mathematical model for predicting the degradation of the Q-factor of vacuum packaged MEMS sensors when exposed to high temperatures. Here is a description of how the mathematical model is constructed.

Material outgassing from vacuum-packaged MEMS sensors is considered an unsteady process, and the diffusion coefficient in different directions is regarded as a constant. Based on Fick's second law:

$$\frac{\partial\_n}{\partial\_l} = D\nabla^2 n \tag{1}$$

In this formula, *n* represents the number of gas molecules per unit cross-sectional area, *t* represents the diffusion time, and *D* represents the diffusion coefficient. The Gilliland formula can be used to express the diffusion coefficient *D* as follows:

$$D = \frac{435.7\sqrt{2}\sqrt{T^3}}{2p\sqrt{\mu}\sqrt[3]{V}}\tag{2}$$

where *T* denotes the thermodynamic temperature, *p* is the total pressure, *µ* is the molecular weight, and *V* is the liquid gram molar volume of the gas at its normal boiling point.

After vacuum packaging, the MEMS sensor cannot maintain an absolute vacuum state in its cavity, so there will be a certain amount of gas in the initial state, which is directly related to the initial Q-factor of the samples. Considering this, assuming there are *N*<sup>0</sup> gas molecules in the sealed cavity at the beginning, and *N*tot is the total number of gas molecules released, the boundary condition for the number of molecules inside the cavity *N*(*t*) over time is:

$$N(t) = \begin{cases} N\_{0\prime}t = 0\\ N\_0 + N\_{\text{tot}\prime}t = \infty \end{cases} \tag{3}$$

Based on the above conditions, the Fourier series method may be used to calculate the diffusion equation as follows:

$$n(z,t) = \frac{N\_{\text{tot}}}{V} + \sum\_{i=1}^{\infty} \frac{2N\_{\text{tot}}}{V} \times e^{-\lambda\_i Dt} \times \cos\left(\frac{i\pi z}{l\_z}\right) \tag{4}$$

As a result, the number of free gas molecules can be expressed as follows:

$$\begin{split} \frac{1}{\mathbb{Q}(t)} = \mathbf{A} \cdot \mathbf{N}(t) &= \mathbf{A} \cdot \left[ \mathbf{N}\_0 + \mathbf{N}\_{\text{tot}}(t) \right] = \mathbf{A} \cdot \left[ \frac{\mathbf{A}}{\mathbb{Q}\_0} + \mathbf{N}\_{\text{tot}} \cdot \left( \mathbf{1} - \mathbf{e}^{-D\_T \cdot t} \right) \right] \\ &= \left( \frac{1}{\mathbb{Q}\_0} + \mathbf{A} \cdot \mathbf{N}\_{\text{tot}} \right) - \mathbf{A} \cdot \mathbf{N}\_{\text{tot}} \cdot \mathbf{e}^{-D\_T \cdot t} \end{split} \tag{5}$$

MEMS sensors that are vacuum packaged can be described by the following vacuum degradation model:

$$p(t) \propto \frac{1}{Q(t)} \propto N(t) = a - b \times \mathbf{e}^{-c \times t} \tag{6}$$

In the formula, *a* = 1/*Q*<sup>0</sup> + A × *N*tot, *b* = A × *N*tot, *c* = *DT*, and the physical meanings of model parameters *a*, *b*, and *c* are as follows:

Parameter *a* represents the final total number of gases inside the sample cavity, which determines the final Q-factor of the sample. Depending on parameter *b*, the degree of degradation of the sample is determined by the number of gases released inside the sample cavity. Parameter *c* characterizes the degradation speed of the sample, which is related to temperature.

With time, the vacuum degree of the vacuum-packaged MEMS sensor decays exponentially. By fitting the experimental data to the above-mentioned model, parameters *a*, *b*, and *c* can be obtained. The above analysis also reveals that 1/(*a*−*b*) is the initial value of the Q-factor of the sample.

## *2.2. Experiment of Reliability Evaluation*

## 2.2.1. Sample

As it is difficult to detect and activate the MEMS microstructure of the MEMS products electrically, we customized a batch of samples. It is typical to find MEMS gyroscopes in commercial applications. The structure of the device must be in a resonant state, and it is usually packaged in a vacuum. In order to investigate vacuum degradation of vacuum packaged MEMS devices, this paper selects the sensitive structure of a MEMS gyroscope as the subject of experimental research.

As shown in Figure 1b, MEMS gyroscopes have comb structures as their internal structure. During the operation of the gyroscope, the position of the combs will change, as will the capacitance between them. It is possible to determine the value based on the change in capacitance. As a result, the vacuum encapsulated gyroscope will have difficulty moving the comb, which will adversely affect the gyroscope's performance.

change in capacitance. As a result, the vacuum encapsulated gyroscope will have diffi-

As shown in Figure 1a, three types of device-level packaged MEMS gyroscopes were

culty moving the comb, which will adversely affect the gyroscope's performance.

customized with different vacuum degrees.

**Figure 1.** Schematic illustration of (**a**) test sample and (**b**) internal structure diagram. **Figure 1.** Schematic illustration of (**a**) test sample and (**b**) internal structure diagram.

Table 1 provides information regarding the vacuum packaging of the samples. There were a total of 18 samples, which were divided into three groups. Samples that were As shown in Figure 1a, three types of device-level packaged MEMS gyroscopes were customized with different vacuum degrees.

highly vacuum packaged were H1 and H6, medium vacuum packaged samples were M1 and M6, and low vacuum packaged samples were L1 and L6. **Table 1.** Design values of package air pressure and Q-factor for three kinds of samples. Table 1 provides information regarding the vacuum packaging of the samples. There were a total of 18 samples, which were divided into three groups. Samples that were highly vacuum packaged were H1 and H6, medium vacuum packaged samples were M1 and M6, and low vacuum packaged samples were L1 and L6.

**Sample Serial Number Package Air Pressure Q-Factor Table 1.** Design values of package air pressure and Q-factor for three kinds of samples.


#### Considering that the Q-factor can be used to characterize the vacuum of a MEMS 2.2.2. Test Method and System

sensors packaged in a vacuum.

cavity, the Q-factor of our sample was evaluated. As shown in Figure 2, we developed a Q-factor test method based on transient excitation [16] for the above-mentioned test samples with different levels of vacuum. Q value measurement is based on the oscillating principle [17]. Labview software, an NI PXI 4461 acquisition card, an analog circuit, and a MEMS sensor are included in the Considering that the Q-factor can be used to characterize the vacuum of a MEMS cavity, the Q-factor of our sample was evaluated. As shown in Figure 2, we developed a Q-factor test method based on transient excitation [16] for the above-mentioned test samples with different levels of vacuum. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 5 of 12

**Figure 2.** Q-factor test frequency response diagram and test system. **Figure 2.** Q-factor test frequency response diagram and test system.

2.2.3. Experiment The environment has a significant impact on the reliability of vacuum-packaged MEMS sensors. It has been reported that the Q-factor of vacuum packaged MEMS gyroscopes will degrade with time at room temperature [13]. Based on theory and analysis, it has been tentatively determined that temperature affects the outgassing of internal materials. At the beginning of the research, a vacuum packaged MEMS gyroscope was baked at three different high temperatures and the Q-factor was monitored over time. According to Figure 3, the Q-factor of the MEMS gyroscope dropped sharply at high temperatures. Accordingly, the analysis of the theoretical and experimental results indicates that high Q value measurement is based on the oscillating principle [17]. Labview software, an NI PXI 4461 acquisition card, an analog circuit, and a MEMS sensor are included in the test system. A part of the test process is as follows: by applying excitation signals to the resonant structure of the MEMS sensor using the acquisition card, Labview software causes the sensor to enter the steady state. After stopping the excitation, the sensor is in a state of damped vibration. As a final step, the acquisition card collects the transient response voltage signal and passes it to Labview for analysis and processing. Using Labview software, it is possible to determine the resonant frequency and Q value rapidly by using wave peak detection technology, the cycle average method, and exponential fitting algorithm.

temperature is the main stress contributing to the degree of vacuum degradation of MEMS

**Figure 3.** Degradation of Q-factor of vacuum packaged gyroscope with time at high temperature.

At different temperatures, we conducted accelerated degradation experiments on the samples using the high-temperature accelerated test method based on MIL-STD-883. As shown in Table 2, samples with high vacuum, medium vacuum, and low vacuum were

**Environmental Stress Samples Test Parameter** 

105 °C H3、H4、M3、M4、L3、L4 Q-factor

(1) Prior to the test, the Q-factor of the samples is measured at room temperature; (2) Incubator temperature is raised to the specified temperature, and samples are

85 °C H1、H2、M1、M2、L1、L2

125 °C H5、H6、M5、M6、L5、L6

In order to conduct a reliability test, the following steps must be followed:

placed in the incubator after the temperature reaches the specified temperature;

tested at three temperatures of 85 °C, 105 °C, and 125 °C.

**Table 2.** Test conditions and grouping.

#### 2.2.3. Experiment 2.2.3. Experiment The environment has a significant impact on the reliability of vacuum-packaged

The environment has a significant impact on the reliability of vacuum-packaged MEMS sensors. It has been reported that the Q-factor of vacuum packaged MEMS gyroscopes will degrade with time at room temperature [13]. Based on theory and analysis, it has been tentatively determined that temperature affects the outgassing of internal materials. MEMS sensors. It has been reported that the Q-factor of vacuum packaged MEMS gyroscopes will degrade with time at room temperature [13]. Based on theory and analysis, it has been tentatively determined that temperature affects the outgassing of internal materials.

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 5 of 12

**Figure 2.** Q-factor test frequency response diagram and test system.

At the beginning of the research, a vacuum packaged MEMS gyroscope was baked at three different high temperatures and the Q-factor was monitored over time. According to Figure 3, the Q-factor of the MEMS gyroscope dropped sharply at high temperatures. Accordingly, the analysis of the theoretical and experimental results indicates that high temperature is the main stress contributing to the degree of vacuum degradation of MEMS sensors packaged in a vacuum. At the beginning of the research, a vacuum packaged MEMS gyroscope was baked at three different high temperatures and the Q-factor was monitored over time. According to Figure 3, the Q-factor of the MEMS gyroscope dropped sharply at high temperatures. Accordingly, the analysis of the theoretical and experimental results indicates that high temperature is the main stress contributing to the degree of vacuum degradation of MEMS sensors packaged in a vacuum.

**Figure 3.** Degradation of Q-factor of vacuum packaged gyroscope with time at high temperature. **Figure 3.** Degradation of Q-factor of vacuum packaged gyroscope with time at high temperature.

At different temperatures, we conducted accelerated degradation experiments on the samples using the high-temperature accelerated test method based on MIL-STD-883. As shown in Table 2, samples with high vacuum, medium vacuum, and low vacuum were tested at three temperatures of 85 °C, 105 °C, and 125 °C. At different temperatures, we conducted accelerated degradation experiments on the samples using the high-temperature accelerated test method based on MIL-STD-883. As shown in Table 2, samples with high vacuum, medium vacuum, and low vacuum were tested at three temperatures of 85 ◦C, 105 ◦C, and 125 ◦C.

**Table 2.** Test conditions and grouping. **Table 2.** Test conditions and grouping.


In order to conduct a reliability test, the following steps must be followed: In order to conduct a reliability test, the following steps must be followed:

(1) Prior to the test, the Q-factor of the samples is measured at room temperature; (1) Prior to the test, the Q-factor of the samples is measured at room temperature;

(2) Incubator temperature is raised to the specified temperature, and samples are (2) Incubator temperature is raised to the specified temperature, and samples are

placed in the incubator after the temperature reaches the specified temperature; placed in the incubator after the temperature reaches the specified temperature;

(3) The samples are baked at a high temperature, and then taken out of the incubator after *t*<sup>1</sup> min;

(4) The Q-factor of the samples is determined after they have stood at room temperature for two hours;

(5) After the measurements have been completed, the samples are placed back into the incubator;

(6) Repeat steps (3) and (4) until the test is complete.

## **3. Results**

#### *3.1. Test Results*

In Tables 3–5, reliability test results of samples at 85 ◦C, 105 ◦C, and 125 ◦C under high-temperature stress are presented.

**Table 3.** Changes of Q-factor of different samples under high-temperature baking at 85 ◦C.


**Table 4.** Changes of Q-factor of different samples under high-temperature baking at 105 ◦C.


**Table 5.** Changes of Q-factor of different samples under high-temperature baking at 125 ◦C.


As a result of the reliability test, we used the mathematical model presented in this paper for fitting. As can be seen from Figure 4, the typical fitting curve reaches a fitting degree of 1.00, indicating that the mathematical model established in this paper is more capable of reflecting the changes in vacuum degree over time under a high-temperature environment. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 7 of 12

**Figure 4.** Q-factor degradation fitting curve of M4 sample at 105 °C. **Figure 4.** Q-factor degradation fitting curve of M4 sample at 105 ◦C.

Having analyzed the data obtained from the experiment, the fitting parameters for the test samples with different vacuum degrees at three temperature points are determined, as shown in Tables 6–8. Having analyzed the data obtained from the experiment, the fitting parameters for the test samples with different vacuum degrees at three temperature points are determined, as shown in Tables 6–8.

H1 14,281 2.44 × 10−4 1.73 × 10−4 0.0552 0.97 H2 19,390 1.91 × 10−4 1.37 × 10−4 0.0452 0.95 M1 7572 2.82 × 10−4 1.49 × 10−4 0.0576 0.98 M2 7822 2.76 × 10−4 1.47 × 10−4 0.0547 0.97 L1 2118 6.00 × 10−4 1.27 × 10−4 0.0544 0.98 L2 2018 6.12 × 10−4 1.15 × 10−4 0.0517 0.97

**Sample** *Q***<sup>0</sup>** *a b c* **R2** H3 13,048 4.96 × 10−4 4.19 × 10−4 0.0105 1.00 H4 16,388 4.18 × 10−4 3.57 × 10−4 0.0108 1.00 M3 7710 4.82 × 10−4 3.52 × 10−4 0.0118 1.00 M4 7744 4.93 × 10−4 3.64 × 10−4 0.0115 1.00 L3 2072 7.92 × 10−4 3.09 × 10−4 0.0102 1.00 L4 2061 8.06 × 10−4 3.21 × 10−4 0.0107 1.00

**Sample** *Q***<sup>0</sup>** *a b c* **R2** H5 12,504 7.10 × 10−4 5.73 × 10−4 0.0452 0.95 H6 19,571 6.82 × 10−4 5.88 × 10−4 0.0379 0.97 M5 8582 7.34 × 10−4 5.79 × 10−4 0.0354 0.98 M6 8405 6.95 × 10−4 5.38 × 10−4 0.0393 0.97 L5 1967 1.12 × 10−3 5.80 × 10−4 0.0325 0.98 L6 1933 1.13 × 10−3 5.76 × 10−4 0.0335 0.98

The curve fitting degree for 18 samples exceeds 1.00 or is close to it, which confirms the accuracy of our mathematical model. As a result, it can also be concluded that the established model can accurately predict the degree of vacuum degradation of MEMS

**Table 6.** Parameter values of Q-factor degradation fitting curve baking at 85 °C.

**Table 7.** Parameter values of Q-factor degradation fitting curve baking at 105 °C.

**Table 8.** Parameter values of Q-factor degradation fitting curve baking at 125 °C.

sensors that are vacuum packaged.


**Table 6.** Parameter values of Q-factor degradation fitting curve baking at 85 ◦C.

**Table 7.** Parameter values of Q-factor degradation fitting curve baking at 105 ◦C.


**Table 8.** Parameter values of Q-factor degradation fitting curve baking at 125 ◦C.


The curve fitting degree for 18 samples exceeds 1.00 or is close to it, which confirms the accuracy of our mathematical model. As a result, it can also be concluded that the established model can accurately predict the degree of vacuum degradation of MEMS sensors that are vacuum packaged.

Based on the physical meaning of the fitting parameters, this paper analyses and discusses the obtained parameter values.

Parameter *a* specifies the total number of gases inside the sample cavity and determines the final Q-factor of the sample. According to the parameter values above for the fitting curve parameters, there are 18 samples at three different temperature points, and there is no significant difference between the *a* values of the high vacuum and medium vacuum samples at the same temperature point. This is due to the fact that the amount of gas released from high and medium vacuum samples is significantly greater than the amount of gas released from the initial sample. Compared to the high and medium vacuum samples at the same temperature, the low vacuum sample has a higher *a* value. Since the initial gas number of the low vacuum samples is higher, it is comparable with the total amount of gas released.

In addition, parameter *b* determines the degree of sample degradation by characterizing the total number of gases released inside the sample cavity. At the same temperature, the outgassing amount represented by the *b* value is essentially the same, which means that it has nothing to do with the vacuum degree of the sample, but only with its temperature. The higher the temperature, the greater the amount of outgassing.

Again, *c* characterizes the degradation speed of the samples, that is, the outgassing rate inside the sample, which is related to temperature. In high, medium, and low vacuums, the outgassing rates of samples remain essentially the same, indicating this parameter is not related to the sample's vacuum degree, but only to its temperature.

#### *3.2. Data Analysis*

Parameters *b* and *c* of the six samples were relatively consistent at every temperature point, and the average value was taken as the value of parameter *b* and parameter *c* and summarized in Table 9.


**Table 9.** Test sample parameters *b* and *c* at three temperature points.

Parameter *b* characterizes the total amount of gas released inside the sample, and the total amount of released gas inside the sealed cavity is related to the temperature.

According to the relationship between the total amount of gas out of non-metallic materials and temperature [18]:

$$b = b\_0 \times \mathbf{e}^{-B \times \mathbf{T}} \tag{7}$$

The above equation can be transformed into a linear equation:

$$\text{In}b = \text{In}b\_0 - B \times T \tag{8}$$

The data of parameter *b* and temperature *T* at three temperature points are shown in Table 10.


Taking *T* as the abscissa and In*b* as the ordinate, the linear fitting curve is shown in Figure 5, and the fitting equation is *y* = 0.0350 × *x* − 21.3138. Taking *T* as the abscissa and In*b* as the ordinate, the linear fitting curve is shown in Figure 5, and the fitting equation is = 0.0350 × − 21.3138.

**Figure 5.** Linear fitting curve of *T* and ln*b.* **Figure 5.** Linear fitting curve of *T* and ln*b.*

The following parameters can be calculated from the equation of total gas output with temperature: The following parameters can be calculated from the equation of total gas output with temperature:

In a sealed cavity, the rate of gas release is determined by parameter *c*, and the rate of gas release is related to the temperature inside the sealed cavity. According to theory, the higher the temperature, the greater the rate of outgassing. Figure 6 illustrates the results of the tests conducted in this paper. The outgassing rate at 85 °C is the highest, which is 0.0531; the outgassing rate at 105 °C is the lowest, which is 0.0109; and the outgassing

As a result of the experiments conducted, this paper concludes that the outgassing rate of the gas inside the sample is equal to the desorption rate of the gas minus the adsorption rate of the gas. It is important to note that the maximum outgassing rate occurs

$$
\ln b\_0 = -21.3138; \; B = 0.0350\tag{9}
$$

<sup>0</sup> ln 21.3138 0.0350 *b B* =− = ; (9)

21.3138 0.0350 *<sup>T</sup> b e*− +× = (10)

rate at 125 °C is in the middle, which is 0.0373.

**Figure 6.** Outgassing rate at different temperature points.

Total outgassing volume and temperature are related as follows: Total outgassing volume and temperature are related as follows:

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 9 of 12

Figure 5, and the fitting equation is = 0.0350 × − 21.3138.

Taking *T* as the abscissa and In*b* as the ordinate, the linear fitting curve is shown in

$$b = e^{-21.3138 + 0.0350 \times T} \tag{10}$$

<sup>0</sup> ln 21.3138 0.0350 *b B* =− = ; (9)

In a sealed cavity, the rate of gas release is determined by parameter *c*, and the rate of gas release is related to the temperature inside the sealed cavity. According to theory, the higher the temperature, the greater the rate of outgassing. Figure 6 illustrates the results of the tests conducted in this paper. The outgassing rate at 85 ◦C is the highest, which is 0.0531; the outgassing rate at 105 ◦C is the lowest, which is 0.0109; and the outgassing rate at 125 ◦C is in the middle, which is 0.0373. In a sealed cavity, the rate of gas release is determined by parameter *c*, and the rate of gas release is related to the temperature inside the sealed cavity. According to theory, the higher the temperature, the greater the rate of outgassing. Figure 6 illustrates the results of the tests conducted in this paper. The outgassing rate at 85 °C is the highest, which is 0.0531; the outgassing rate at 105 °C is the lowest, which is 0.0109; and the outgassing rate at 125 °C is in the middle, which is 0.0373.

The following parameters can be calculated from the equation of total gas output

**Figure 5.** Linear fitting curve of *T* and ln*b.*

with temperature:

**Figure 6.** Outgassing rate at different temperature points. **Figure 6.** Outgassing rate at different temperature points.

As a result of the experiments conducted, this paper concludes that the outgassing rate of the gas inside the sample is equal to the desorption rate of the gas minus the adsorption rate of the gas. It is important to note that the maximum outgassing rate occurs As a result of the experiments conducted, this paper concludes that the outgassing rate of the gas inside the sample is equal to the desorption rate of the gas minus the adsorption rate of the gas. It is important to note that the maximum outgassing rate occurs at 85 ◦C as a result of the desorption of physically adsorbed gas atoms or molecules from the surface of the material, whereas a higher temperature is required for gas in the deeper part of the material. The collision frequency between gas molecules or atoms increases at 105 ◦C, which increases the adsorption rate. As a result, the outgassing rate decreases at 105 ◦C. Gas atoms and molecules are desorbed at 125 ◦C by physical adsorption and chemical adsorption, resulting in an increase in the rate at which they are outgassed.

## *3.3. Reliability Prediction*

Thus, we can obtain a mathematical model for describing the vacuum degradation of MEMS sensors as follows:

$$p(t) \propto \frac{1}{Q(t)} \propto N(t) = a + \left(e^{-21.3138 + 0.0350 \times T}\right) \times \left(e^{-c \times t}\right) \tag{11}$$

The majority of vacuum packaged MEMS gyroscopes and other MEMS sensors are used at room temperature, and the obtained mathematical model is used to predict the long-term reliability of MEMS gyroscopes under 25 ◦C vacuum degradation conditions.

The following parameters can be determined at a room temperature of 25 ◦C using the formula deduced above:

$$b = e^{-21.3138 + 0.0350 \times 298} = 0.00001876\tag{12}$$

Using the previous analysis, it can be determined that 1/(*a*−*b*) can obtain the initial value of the Q-factor of samples, and we also know the value of *b*. Therefore, parameter *a*, which represents the final gas volume of the cavity, determines the final value of the MEMS sensor. Since the value of *c* is uncertain, *t* = ∞ is taken.

Table 11 shows the changes in Q-factor during storage (T = 25 ◦C) with time for different initial Q-factor.


**Table 11.** Changes of different initial Q-factor with time during storage (T = 25 ◦C).

#### **4. Discussion**

According to the established mathematical degradation model of vacuum-packaged MEMS sensors, the long-term reliability of vacuum-packaged MEMS sensors at room temperatures (25 ◦C) was predicted. The data are shown in Table 11, the vacuum degree of samples with different initial Q-factor degrades continuously over time; whereas samples with small initial Q-factor do not exhibit obvious degradation trends, samples with larger initial Q-factors demonstrate more obvious degradation trends. Furthermore, the degradation degree of the Q-factor of all samples tends to be flat with the increasing storage time.

At room temperature, the gas molecules adsorbed on the surface material of the sensor cavity are continuously released, which leads to the decrease of the vacuum degree of the sensor, thus reducing the reliability of the sensor. In the process of high-temperature accelerated aging, physical and chemical desorption of gas will occur, especially the gas released by some polymers in the cavity, which has a great impact on the reduction of vacuum degree.

#### **5. Conclusions**

In this paper, a theoretical model for the prediction of vacuum degradation of vacuum packaged MEMS sensors is established. Furthermore, combination with the experiments, the MEMS gyroscopes with high vacuum degree, medium vacuum degree, and low vacuum degree are used as examples to study the vacuum degradation law at different high temperatures and get the parameters of the model. The reliability life prediction of the vacuum packaged MEMS sensor under a room temperature environment (25 ◦C) is also carried out. The results show that our theoretical model can achieve a good prediction of vacuum degradation of vacuum packaged MEMS sensors. This is of great significance for improving the reliability of vacuum-packaged MEMS sensors and promoting the commercial application of vacuum-packaged MEMS sensors. In the next step, we will conduct further research on the variation law of the outgassing rate in the cavity with temperature, such as using a lower temperature test to obtain the relationship between the outgassing rate and temperature and predict the variation of the vacuum degree with time.

**Author Contributions:** Conceptualization, G.D. and X.D.; methodology, G.D. and X.D.; software, X.H.; validation, G.D. and X.H.; formal analysis, G.D. and X.D.; investigation, X.D. and W.S.; resources, X.D. and W.S.; data curation, G.D.; writing—original draft preparation, G.D. and X.D.; writing review and editing, G.D. and X.D.; visualization, G.D. and P.Z.; supervision, P.Z. and X.D.; project administration, P.Z. and X.D.; funding acquisition, X.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Guangzhou Science and Technology Plan Project (Grant Number 202102080190), Industrial Technology Foundation Public Service Platform Project (Grant Number 2021-H021-1-1), and Ministry of Science and Technology Key R&D Program (Grant Number 2020YFB2008900).

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


**Wen-Sheng Zhao 1,\* , Rui Zhang <sup>2</sup> and Da-Wei Wang <sup>1</sup>**


**Abstract:** The advance of semiconductor technology not only enables integrated circuits with higher density and better performance but also increases their vulnerability to various aging mechanisms which occur from front-end to back-end. Analysis on the impact of aging mechanisms on circuits' reliability is crucial for the design of reliable and sustainable electronic systems at advanced technology nodes. As one of the most crucial back-end aging mechanisms, electromigration deserves research efforts. This paper introduces recent studies on physics-based modeling of electromigration aging of on-chip interconnects. At first, the background of electromigration is introduced. The conventional method and physics-based modeling for electromigration are described. Then studies on how electromigration affects powers grids and signal interconnects are discussed in detail. Some of them focus on the comprehensiveness of modeling methodology, while others aim at the strategies for improving computation accuracy and speed and the strategies for accelerating/decelerating aging. Considering the importance of electromigration for circuit reliability, this paper is dedicated to providing a review on physics-based modeling methodologies on electromigration and their applications for integrated circuits interconnects.

**Keywords:** integrated circuit interconnects; aging; reliability; electromigration; physics-based modeling

## **1. Introduction**

Although technology scaling enables integrated circuits (ICs) with higher density and better performance, it is still faced with serious vulnerability to various aging mechanisms appearing from front-end to back-end [1–10]. These aging mechanisms include Bias Temperature Instability (BTI), Hot Carrier Injection (HCI), Random Telegraph Noise (RTN), Gate-Oxide Breakdown (GOBD) at the front-end, Middle-of-line (MOL) time-dependent dielectric breakdown (TDDB), Back-end-of-line (BEOL) TDDB, and Electromigration (EM). BTI, HCI, RTN, and GOBD cause device parameter deviations. MOL and BEOL TDDBs cause *short circuit* between interconnects, while EM increases interconnects' resistance and it eventually results into *open circuit*. These phenomena lead to malfunction in circuits. From BTI to BEOL TDDB, there are numerous studies for their impact on device, circuit, and system performance and reliability [11–30], but they are not the topic discussed here. In this paper, our attention is focused on recent progress in physics-based modeling of EM in on-chip interconnects.

EM is the migration of interconnects' metal atoms after they obtain momentum from moving electrons. Since the interconnects are wrapped by barrier layer, the metal atoms' movement causes depletion regions and it, eventually, causes void nucleation and growth. EM induced degradation and failure is one of the most critical reliability issues for deeply scaled ICs. Such a degradation is expected to get even worse with further technology scaling. It is reported by ITRS-2015 that the operating current density has exceeded 1 MA/cm<sup>2</sup> and is rapidly approaching to 10 MA/cm<sup>2</sup> . Figure 1a shows the experiment and model

**Citation:** Zhao, W.-S.; Zhang, R.; Wang, D.-W. Recent Progress in Physics-Based Modeling of Electromigration in Integrated Circuit Interconnects. *Micromachines* **2022**, *13*, 883. https://doi.org/ 10.3390/mi13060883

Academic Editor: Aiqun Liu

Received: 16 May 2022 Accepted: 29 May 2022 Published: 31 May 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of lifetime scaling versus interconnect geometry, where an effective scaling model has been established by assuming the void is located at the cathode end of the interconnect which contains a single via with drift velocity dominated by interfacial diffusion [31]. The interconnects' EM lifetime is predicted to decrease by half for each new generation. Figure 1b shows evolution of *Jmax* (the maximum equivalent dc current density) and *JEM* (the current density for target EM lifetime) [31]. Both *Jmax* and *JEM* are limited by interconnect geometry. *Jmax* increases with the increase of operating frequency and the reduction of interconnect cross-section. Although there are several process options, such as usage of Cu alloys seed layer and the short length effect to overcome EM severity, EM is still an inescapable topic since the operating current density is entering *JEM* region, as shown in Figure 1b. It is necessary to apply novel models and methodologies to study the strategies which can help mitigate EM degradation by incorporating the innovations in material and process.

**Figure 1.** (**a**) Experiment and model of lifetime scaling versus interconnect geometry; (**b**) evolution of *Jmax* (the maximum equivalent dc current density) and *JEM* (the current density for target EM life-time). They are predicted by ITRS-2015 [31].

Depending on their functionality, on-chip interconnects can propagate signal within and between cells and deliver power to sub-circuits. Most of the previous EM studies have focused on the reliability of power grids and clock/signal nets. In traditional approach, the Blech limit and Black's equation are applied to investigate interconnects' EM reliability. The Blech limit or Blech product is [32]

$$(j \times L) \le (j \times L)\_{crit} = \frac{\Omega \sigma\_{crit}}{eZ\rho} \tag{1}$$

where *j* is the current density, *L* is the interconnect branch length, Ω is the atomic volume, *σcrit* is the critical stress for void nucleation, *e* is the electron charge, *eZ* is the effective charge of migrating atoms, and *ρ* is the metal resistivity.

The interconnect branches with (*j* × *L*) ≤ (*j* × *L*)*crit* are filtered out as EM immortal. The mean time to failure of remaining EM mortal branches is characterized with Black's equation [33]

$$MTTF = A j^{-n} \exp\{E\_a/k\_B T\} \tag{2}$$

where *A* and *n* are assumed to constants, *E<sup>a</sup>* is the activation energy, *k<sup>B</sup>* is the Boltzmann's constant, and *T* is the absolute temperature.

*A* and *n* are measured through accelerated test with a higher current density injected under higher temperature. Then the *MTTF* under normal use condition is extrapolated as a function of *MTTF* under accelerated test,

$$MTTF\_{\text{use}} = MTTF\_{\text{stress}} \left(\frac{j\_{\text{stress}}}{j\_{\text{use}}}\right)^{\text{n}} \exp\left\{\frac{E\_a}{k\_B} \left(\frac{1}{T\_{\text{use}}} - \frac{1}{T\_{\text{stress}}}\right)\right\} \tag{3}$$

where *jstress* and *Tstress* are current density and absolute temperature under accelerated test, while *juse* and *Tuse* are current density and absolute temperature under normal-use condition, respectively.

Based on *MTTFuse*, the probability of failure (at a specific stress time) of each branch is obtained from cumulative distribution function of its lifetime distribution (lognormal or weibull). Then the overall probability of failure of studied branches is computed by weakest-link statistics. Although this traditional approach is convenient to use, it does not provide accurate estimation on EM failure due to the following reasons: first, n is not really a constant under different current density especially when *jstress* is much larger than *juse*; second, the Belch limit is accurate only for single branch interconnect. It is inapplicable for a general interconnect tree with multiple segments/branches. As shown in Figure 2, since the segments are connected without barriers, it is necessary to consider the atom transportation between them [34]. Third, weakest-link statistics considers a single branch failure leads to failure of studied interconnect trees which is not true in state-of-the-art ICs, especially for power grids with a mesh structure [35]. Therefore, some physics-based EM compact models have been proposed to incorporate the atom movement between segments [36–44]. It is noted that there are already several excellent reviews on relevant topics in [34,45–48], but they are mainly on power grids. In this review article, we will summarize the studies on physics-based EM compact modeling and its applications on various on-chip interconnects and their most recent updates. The remainder of this paper is organized as follows. Section II introduces the physics-based three-phase compact EM model. The model should be imposed on interconnects with corresponding current density and boundary conditions. The applications of this compact EM model in current and future interconnects are discussed as well. Section III summaries recent progress in EM studies mainly on power grids. These studies cover from Black's equation-based simulations to physics-based simulations. Section IV shows the application of physicsbased modeling on interconnects in cache memory. The impact of physical dimensions and cache configurations on EM reliability is discussed. Section V concludes this paper.

**Figure 2.** An interconnect tree with multiple segments confined by diffusion barriers/liners in one-layer metallization [34].

#### **2. Physics-Based Three-Phase Compact Model**

EM is the phenomenon of metal atoms migration due to the momentum obtained from electrons when there is current flowing over interconnects. Since the momentum is also transferred between metal atoms, hydrostatic stress would appear in the interconnects confined by barrier material and capping layer. Although the stress gradient may prevent atoms to move, the maximum stress continues to increase until it reaches the saturation value or a critical stress value for void generation. A void is considered to be nucleated at the position where the maximum hydrostatic stress reaches critical stress value. Then the void is incubated and it grows under further stimulus by electrons' movement. A void will never appear in an interconnect if the saturation stress value is smaller than the critical value. Such an interconnect is EM immortal. Figure 3a shows a Cu dual damascene conductor structure [49]. The trench is lined with a Cobalt-based liner and the Cu is capped with either a dielectric or metallic layer. Figure 3b shows the voids possibly formed in an interconnect. The void formed under a via can cause *early failure* when the size is large enough to cover the via bottom and it leads to an open circuit status, while the void formed above a via causes *conventional failure* when the gradually growing void introduces an obvious resistance shift which causes circuit's functional failure with a significant performance deviation. Both *early* and *conventional failures* need to be checked during EM analysis. The EM process during void generation and growth can be described with a three-phase model which consists of nucleation phase, incubation phase, and growth phase. To calculate the time-dependent EM process, Korhonen's model has been proposed to capture the time-varying stress. The Korhonen's model together with appropriate initial condition and boundary conditions are able to accurately emulate hydrostatic stress evolution in interconnects. Thereafter, the time-dependent stress and atomic flux are combined to obtain resistance shift due to EM. The resistance shift can be inserted into circuit simulation to show its impact on performance. This section introduces Korhonen's model and the three-phase compact EM model.

**Figure 3.** (**a**) Cu dual damascene conductor structure, where the trench is lined with a Cobalt-based liner and the Cu is capped with either a dielectric or metallic layer. (**b**) Void formation in a Cu wire connected with a via-below (left) and a via-above (right) [49].

#### *2.1. Korhonen's Model*

In general, multi-physics 3-D simulation is able to present accurate estimation on hydrostatic stress and resistance shift in interconnects. However, since it requires vast computation resource and the simulation speed gets very slow when the studied interconnect structure (such as power grids for a microprocessor) is relatively large, 3-D numerical simulation is not always preferred for EM analysis. Fortunately, Korhonen's 1-D model offers a reasonable trade-off between accuracy and efficiency. Let us take a branch with a length of *L* as an example. The evolution of 1-D hydrostatic stress follows [50]

$$\frac{\partial \sigma}{\partial t} = \frac{\partial}{\partial \mathbf{x}} \left[ \frac{D\_d \mathcal{B}}{k\_B T} \left( \Omega \frac{\partial \sigma}{\partial \mathbf{x}} - e \mathbf{Z} \rho\_{\mathbb{C}u} j \right) \right] \tag{4}$$

$$J\_a = \frac{D\_a \mathbb{C}\_a \Omega}{k\_B T} \left[ \frac{\partial \sigma}{\partial x} - \frac{eZ\rho\_{\mathbb{C}u} j}{\Omega} \right] \tag{5}$$

where the effective atomic diffusivity *D<sup>a</sup>* = *D*0*exp*(−*Ea*/*kBT*). *D<sup>0</sup>* is a constant, *E<sup>a</sup>* is the effective activation energy, *k<sup>B</sup>* is Boltzmann's constant, *T* is absolute temperature, *B* is the effective bulk elasticity modulus, Ω is the atom lattice volume, *eZ* is the effective charge of migrating atoms, *ρCu* is copper resistivity, *j* is the current density over the studied interconnect, *C<sup>a</sup>* is the atom concentration, *J<sup>a</sup>* is the atomic flux, and *t* is the stress time.

Initial conditions and boundary conditions are necessary for a solution of hydrostatic stress. They are decided by the material difference in the coefficients of thermal expansion (*CTE*) and the difference between stress free annealing temperature and the circuit temperature. The position dependent (*x* ranges from *0* to *L*, *x* = *0* is at the branch's left end, and *x* = *L* is the branch's right end when it is assumed to be placed horizontally) initial value (at *t* = 0) is expressed as [51]

$$\sigma(\mathbf{x}, t=0) = B \left( a\_M - a\_{\text{Conf}f} \right) \left( T\_{\text{ZS}} - T(\mathbf{x}, t=0) \right) \tag{6}$$

where *TZS* is the stress free annealing temperature, *x* denotes the node position, *T*(*x*, *t* = 0) is the specific node temperature at *t* = 0, and *a<sup>M</sup>* and *aConf* are the *CTE* of the metal and confinement materials, respectively.

The studied branch may be connected to other branches at its left and right ends. The boundary conditions at the two ends depend on whether they are connected to other branches and whether the voids have appeared at their locations. At the left and right ends, the boundary conditions before and after void nucleation are given as [41,44]

$$\begin{cases} \sum\_{k=0}^{K\_l} J\_{a,k} w\_k h\_k = 0, \text{ for left node if it's nucleation phase} \\ \frac{\partial \sigma(\mathbf{x} = 0, t)}{\partial \mathbf{x}} = \frac{\sigma(\mathbf{x} = 0, t)}{\delta}, \text{ for left node if it's in inclusion and growth phases} \\ \sum\_{k=0}^{K\_r} J\_{a,k} w\_k h\_k = 0, \text{ for right node if it's in nucleation phase} \\ \frac{\partial \sigma(\mathbf{x} = L, t)}{\partial \mathbf{x}} = -\frac{\sigma(\mathbf{x} = L, t)}{\delta}, \text{ for right node if it's in inclusion and growth phases} \end{cases} (7)$$

where *K<sup>l</sup>* and *K<sup>r</sup>* are the total number of branches connected to the left and right nodes, respectively, *w<sup>k</sup>* and *h<sup>k</sup>* are the width and thickness of their *kth* connected branch, and *δ* is the thickness of the void surface. It is noted that the maximum hydrostatic stress (during void nucleation) in one branch appears at its boundary nodes and the void may not appear at the boundary nodes if the hydrostatic stress at their position never exceeds critical stress (*σth*), then the boundary condition should be always applied as the one in nucleation phase.

#### *2.2. Three-Phase Compact Model*

The three-phase compact model includes nucleation, incubation, and growth phases, as shown in Figure 4 [45]. When the interconnect branches are flown over by currents, the hydrostatic stress on them evolves with the initial value described in Section 2.1. In the nucleation phase, the resistance of a studied branch remains unchanged before the maximum stress exceeds a threshold (*σth*), while after it exceeds the threshold, the model enters the incubation phase where the void size is increasing but it is still smaller than a threshold size. The interconnect resistance is stays unchanged in nucleation and incubation phases. Later, the model enters growth phase if the time-dependent void size grows large enough to cover the interconnects' cross section. The resistance jump is caused by Joule heating. In the growth phase, the branch's time-dependent resistance shift can be expressed as [49]

$$
\Delta R = L\_v(t) \left( \frac{\rho\_{linear}}{(W + 2H + 2t\_{linear})t\_{linear}} - \frac{\rho\_{cu}}{WH} \right) \tag{8}
$$

where *W* and *H* are the width and height of the studied branch, respectively, *ρliner* and *tliner* are the resistivity and thickness of the liner. *Lv*(*t*) is the time-dependent void length, and its growth can be calculated with ∆*Lv*(*t*) *= Ja*(*t*)·Ω·∆*t*.

**Figure 4.** Resistance change over time under the three-phase EM model [45].

The three-phase model can be validated by experimental data. Figure 5 shows resistance trace of interconnects in experimental measurements [52]. Obviously, most of the traces follow resistance shift behavior described by the three-phase EM model.

**Figure 5.** Resistance trace of interconnects under EM testing [52].

#### *2.3. Model Applications*

The Korhonen's model describes time-dependent evolution of hydrostatic stress in specific branch. Since the interconnects are confined by diffusion barriers/liners in onelayer metallization, the whole interconnect network can be divided into many individual interconnect trees which are similar as the one shown in Figure 2. For interconnect tree, their hydrostatic stress evolution can be evaluated independently. It should be noted that the current density over the branches may vary under the appearance of resistance shift of some branches. For each branch in a specific tree, the Korhonen's model can be applied

with appropriate initial condition and boundary conditions to accurately compute the hydrostatic stress evolution. Based on this, the time-dependent resistance shift of the EM mortal branch is obtained conveniently. The model has been applied to check the reliability of state-of-the-art Cu interconnects. More details are given in the following sections. With the continuous scaling, the Cu interconnects in sub-10 nm technology node suffer from high resistance due to serious surface scattering of the electrons flowing over them. It is found in [53,54] that reducing the linewidth to 10 nm results in a drop of *jmax* to below 1 MA/cm<sup>2</sup> and scaling linewidth from 25 nm to 10.5 nm leads to a 90% drop of *jfail* i.e. the current density that induces failure at 10 y. Under such a scenario, interconnects based on Ru and Co are potential replacements with better reliability than Cu because of their lower resistivity and higher EM activation energies [55–60]. The barrierless Ru interconnect together with an integration scheme have been identified to be more EM reliable than Cu [61]. And it is found that full Ru vias have no risk of voiding after long thermal storage at high temperature [62]. With respect to Co interconnects, the first estimation of an effective *<sup>D</sup>0Z\** (~1.72 <sup>×</sup> <sup>10</sup>−10) for Co is performed [63]. It is two orders of magnitude lower than Cu. While Co vias may be EM immortal, the lack of a barrier may induce diffusion along the Co/dielectric interfaces and Co/Cu intermixing [64]. The physics-based compact EM model introduced in this section is also applicable to interconnects based on new materials such as Co and Ru only if the diffusion coefficients and microstructures are extracted from experiment result [65].

## **3. Modeling of EM Impact on Power Grids**

In the most previous studies of EM modeling, the Black's equation and Blech limit are applied to analyze the reliability of signal interconnects and power grids [66–73]. In [66–68], by incorporating the Joule heating effect, Gracieli Posser et al. have developed approaches for modeling and efficient characterization of cell-internal EM and have simulated EM effects on different metal layers at different wire lengths. On the one hand, they found the cell-internal EM reliability can be optimized with layout modification and constraints on output pin position. On the other hand, it is concluded that larger metal layers have smaller EM effects and, consequently, a higher EM lifetime for the wires. Palkesh Jain et al. proposed a SoC-level logic-IP-internal EM verification methodology which provides onthe-fly retargeting capability for reliability constraints [69,70]. The proposed approach is demonstrated on a 28-nm design. Meanwhile, they presented a fast and stochastic analysis methodology to overcome the lifetime under-estimation by conventional methodologies based on weakest-link assumption for EM assessment of click grids and power grids [71,72]. Vidya et al. also applied Black's equation and Blech limit to estimate the self-heating impact on EM reliability of FinFET and GAAFET designs [73]. In order to overcome the reliability under-estimation due to the traditional series model for EM checking and the pessimistic assumptions about the chip workload and the corresponding supply currents, Mohammad Fawaz and Sandeep Chatterjee et al. proposed a framework for EM checking that allows users to specify conditions-of-use type constraints which help capture realistic chip workload and which includes the use of a novel mesh model for EM prediction in the grid, instead of the traditional series model [74–76]. They developed a framework to estimate the change in statistics of interconnects as their effective-EM current varies and developed a novel vector less mesh model technique to estimate the average minimum time-to-failure of a power grid under uncertain workload. Their results indicate that the series model causes pessimistic estimation of power grid MTF and reliability by a factor of 3–4 [75,76].

Although the novel methodologies/frameworks based on Black's equation can estimate EM reliability, since the atom flow within segment trees is ignored there, they are not able to provide suitable and accurate results as physics-based methods. Sandeep Chatterjee et al. found that the power grid's lifetime estimated by their physics-based approach is on average 2.75× longer than those based on Black's model [44]. Xin Huang et al. verified that the lifetime of IBMPG2 predicted by the traditional approach with the series

and mesh model is 7.82 y and 10.67 y while the lifetime predicted by a physics-based EM model is 15.66 y [37]. It means the frameworks based on conventional estimation method is too conservative, therefore, the physics-based EM model is preferred for designs with space limited. In this section, our attention is focused on the studies which applied physics-based model to investigate the EM impact on power grids. This section consists of three subsections. The first subsection covers modeling and simulation methodologies proposed to study power grids' reliability. The second subsection focuses on the strategies adopted to improve EM evaluation speed while ensuring good accuracy and the optimization methodologies for better EM reliability. The third subsection is mainly on techniques for EM acceleration and deceleration.

#### *3.1. Modeling and Simulation Methodologies*

There are many studies on modeling and simulation methodologies for EM reliability of power grids with applying physics-based models. Some of them are based on numerical simulations while the others are based on analytical solutions. At an early time, Vivek Mishra et al. modeled the impact of EM in Cu interconnects on power grid integrity with using probability analysis [36,38–40,43]. Figure 6 shows the CDF plots for IR drop of a studied power grid for different circuit lifetime [40]. For tlife = 5 y, the studied power grid remains functional because all samples' IR drop are under a 10% threshold, however, it has a worst resistance degradation of 48% which is much more than a typical 10%~20% resistance increase criteria used by circuit designers. It verifies that power grids have inherent resilience to EM failures. They also studied circuit delay variability due to interconnects' resistance shift under AC EM [38]. It shows that even non-catastrophic EM on critical paths can cause serious performance degradation which ultimately result into circuit malfunction. As shown in Figure 7, the impact of EM on absolute delay shifts increases under technology scaling [38]. It is mainly because the higher current density in smaller interconnects exacerbate EM degradation. Under a specific technology node, EM effect becomes more obvious because EM void size and number increase with stress time. Meanwhile, they developed methods to evaluate transient stress evolution in interconnects, and presented simple and practical criteria for EM mortality checking. It is demonstrated that the number of EM mortal interconnects highly depend on lifetime target and reliability expectations [39]. It is also observed in [43] that power grids' EM degradation is impacted by configuration of via arrays which connect interconnects at different layers.

**Figure 6.** CDF plots for IR drop of a studied power grid for different circuit lifetimes, tlife [40].

**Figure 7.** Absolute delay shifts due to various aging mechanisms for advanced technology nodes at circuit operation time of 5, 10, 15, and 20 y [38].

In order to accurately model EM degradation in power grids, Xin Huang et al. proposed a new physics-based assessment method [37,42]. This method incorporates power grids' redundancy by assuming the circuits get failed only when the IR drop reaches a threshold value. The hydrostatic stress in each interconnect tree is evaluated with considering atom transportation between connected branches. The experimental result not only verifies that the result obtained from Black's equation is too pessimistic, it but also shows that IBM P/Gs' EM failure is more likely to happen at the places with large initial stress value and is more likely to happen at longer time when the void volume saturation phenomenon is taken account. The time-dependent IR drop can be captured by placing the EM induced resistance shift into the P/G circuit model. Figure 8 shows the time-dependent voltage drop of the first failed node and maximum voltage drop in IBMPGNEW1 [42]. The proposed method is also applied to study the impact of cross-layout temperature and thermal stress distributions on full-chip EM assessment.

**Figure 8.** Time-dependent voltage drop of the first failed node and maximum voltage drop in IBMPGNEW1 [42].

Since most of the previous studies are based on the uniform temperature assumption, Xin Huang et al. implemented a flow of EM assessment for multi-layer P/G in a 32 nm test chip [77,78]. The cross-layout temperature variation due to devices' power consumption and interconnects' Joule heating are characterized and incorporated into EM assessment. It is found that uniform temperature assumption causes inaccurate prediction on TTF and the thermal stress variation results into better evaluation on EM reliability. It means the on-chip temperature variation is needed to get more reasonable EM assessment results. With applying the same physics-based EM model, Kai He et al. proposed a lightweight on-chip aging sensor [79]. The interconnects in this sensor are designed to have detectable EM failure at specific time. A number of parallel interconnects are used in the sensor to mitigate inherent variations. The EM-based aging sensor can provide more accurate prediction of chip usage time and offers simpler circuit implementation and smaller area footprints than the ring-oscillator based sensor. Chase Cook et al. has applied finite difference method to solve 1-D EM problem in multi-branch interconnects [80]. The new method can easily accommodate non-uniformly distributed residual stress and timedependent temperature and current during circuit operation. The numerical results match well with that of COMSOL which is based on finite element method. Taeyoung Kim et al. presented an approach for system-level EM reliability management for multi/many core microprocessors [81]. They proposed a task migration method to balance EM resource consumption by all the cores. It treats TTF as a resource to consume during task execution and uses task migration to balance TTF consumption across the cores. It equalizes the probability of failure of each core to maximize the lifetime of multi-core system. The simulation results show a balanced TTF consumption by the cores and the system's EM reliability has been maximized.

Since the current in P/G is unidirectional, the interconnects' EM immortality can be determined by checking steady-state stress distribution. If the maximum stress is larger than critical value, the void is nucleated in a mortal interconnect, otherwise, the interconnect is immortal. However, since there is not closed form for steady state stress in multi-branch interconnect trees, numerical solutions generally need long time simulation. New techniques for convenient immortality checking are necessary for EM study. In [82], Zeyu Sun et al. proposed a new parameter called Critical EM voltage (*V*Crit,EM) to evaluate EM immortality at steady-state stress in multi-branch interconnect tree. The *V*Crit,EM is an extension of Belch limit concept. The difference is that Blech limit is for single branch while *V*Crit,EM is applicable to multi-branch tree. Since this voltage-based EM (VBEM) method overcomes the problem of current-density-based criteria, it can handle the impact of interconnect tree structure on EM-induced stress with easy implementation. With this method, the EM voltage at the ground node or cathode node of a tree is determined with the total area of branches in the tree, the total area of the branches connected to each node in the tree, and the nodal voltage at each node. The EM immortality of studied tree depends on whether the EM voltage is smaller than the *V*Crit,EM. This criterion is applicable for immortality checking in void nucleation phase. The VBEM analysis not only agrees with results from finite difference method at steady state but also matches well with COMSOL and XSim results, as shown in Figure 9. Since the VBEM analysis assumes that current density is evenly distributed in one branch, the impact on current crowding is investigated. Comparison with COMSOL result shows that current crowding effect is unobvious if the length of a branch is much greater than its width. Void saturation volume is another important issue for EM immortality checking. Since previous saturation volume model only works for single branch, Zeyu Sun et al. derived a void saturation volume model for multi-branch tree, as shown in Figure 10a [83]. The model follows atom conservation at steady state of void growth phases. The overall void saturation volume in a tree at steady state is the sum of saturation volume contributed by each branch. The tree is considered as EM immortal if the overall void volume is smaller than the critical size. This criterion is applicable for immortality checking in void incubation phase. Transient analysis is necessary on the tree only when it is EM mortal. The authors proposed a new EM

immortality checking flow which considers the checking criteria in both void nucleation phase and void incubation phase. The algorithm is given in Figure 10b. The new flow reduced conservativeness of existing EM assessment methods, and it helps the designer quickly identify the new type of immortal branches which have void nucleated but with a size smaller than the critical value.

**Figure 9.** Steady-state EM stress comparisons for each straight-line 3-terminal interconnect case [82].

**Figure 10.** (**a**) EM immortality check algorithm flow, and (**b**) simulation framework for *EMSpice* simulator [83,84].

Based on their EM immortality check algorithm shown in Figure 10a and the existing transient EM analysis theories, Zeyu Sun et al. developed a new full-chip EM simulator called *EMSpice* to evaluate EM reliability of P/Gs [84]. Figure 10b shows the simulation flow of *EMSpice*. It starts from P/G layout information from Synopsys IC Compiler. In the first step, it disregards immortal trees by considering the immortality criteria in nucleation and nucleation phases. Then, the mortal trees are applied with a FDTD solver to extract time-dependent hydrostatic stress in both nucleation and post-voiding phases. Since the EM-induced resistance shift cause current variation in P/G, the EM analysis is interacting with IR drop analysis of a whole P/G at each time step to ensure the comprehensiveness of *EMSpice*. *EMSpice* simulator can reduce the over-conservation in EM assessment. It predicts the failed tree number 76.7% less than the Black's method and 66.7% less than another full-chip EM analysis method.

EM postvoiding analysis attracts researchers' attention because it is hard to handle the interactions between current density, hydrostatic stress, and temperature when the target is to investigate detailed void growth. Hengyang Zhao et al. proposed a multi-physics finite-element-method-based (FEM-based) analysis method for void growth simulation in confined copper interconnect [85,86]. The method considers time-varying interactions between hydrostatic stress in the confined interconnects structure illustrated in Figure 11a, the current density and Joule heating induced temperature rise. The interactions are realized by solving a set of coupled partial differential equations, including the Korhonen's equation, the phase field equation, the Laplace equation, and the heat diffusion equation. The FEM-based EM solver is capable to predict unique transient resistance change for copper interconnects. Figure 11b shows the time-dependent total resistance, hotspot temperature, resistance at the hotspot, and void size. It verifies the statement in Section 2.2 that the resistance jump between nucleation phase and growth phase is due to Joule heating. Meanwhile, the lifetime distribution from this EM solver can provide a higher fitting accuracy on the current density exponent parameter (*n*) described in Black's equation in the previous study, as shown in Figure 11c.

Sandeep Chatterjee et al. also proposed a methodology to check P/G EM by using physics-based model [41,44]. They firstly listed the detailed steps to extend the physical models for EM in branches to compute the hydrostatic stress evolution in multi-branch interconnect trees. The boundary condition for branches under different scenarios are provided as well. Then filtering and predictor-based schemes are designed to speed up the overall EM assessment. It enables the statistical computation on IBM benchmarks finish in 2.3 hrs. Therefore, the proposed method has potential to be applied on large-scale circuits. The simulation results verified the inaccuracy of Black's model and explained the importance of *early failures* for EM assessment. As shown in Figure 12, exclusion of early failures leads to optimistic evaluation on voltage drop increase and on MTF of P/Gs [44]. In [87], a finite difference method-based EM analysis methodology is applied to 3-D IC test structure. Its comparison with finite element analysis and experiment measurement demonstrates that EM in 3D IC structures can be suitably evaluated with finite difference simulation. In [88], the authors presented a systematic approach to resize the grid metal lines to achieve a design target lifetime at the minimal extra cost of metal area. With the help of this approach, on a grid with 1.2 million nodes, the authors can increase its MTF from 10.5 y to 12.2 y under a cost of 0.02% extra metal area by scaling only 14 interconnect trees. Current density variation is an important factor for EM evaluation. To ensure the comprehensiveness of their EM modeling, Adam Issa et al. has investigated EM checking by using stochastic effective current model [89]. It is observed that current variations bring us worse EM reliability. Based on his EM modeling experience, Farid N. Najm derived the equivalent circuits for EM under different model phases [90]. It is shown that the dynamic behavior of stress and flux in metal line is identical to dynamic behavior of voltage and current in RC circuit, thus EM assessment can be executed by simulating its equivalent RC circuit. It has potential to drastically improve EM assessment capability on large circuits. In [91], an industry-level physics-based tool for EM assessment in commercial-grade PDNs was introduced. As shown in Figure 13, after analysis the tool can highlight voided metal line segments with a voiding probability. The tool's accuracy has been validated with the experimental data [92]. There is good fit between lifetime statistics derived from measurement (*MTTFEXP* = 62,305 s, ∆EXP = 14,012 s) and simulation (*MTTFSIM* = 60,344 s, ∆SIM = 12,613 s).

**Figure 11.** (**a**) 3D illustration of up-stream interconnect structure and simulated physical systems; (**b**) simulation of Joule heating effect. Rtotal: total wire resistance. T: temperature at the hotspot. Rhs, in ohms: copper resistance at the hotspot. Void, in cubic micrometer: simulated void size; (**c**) extracted current density exponent compared to experiment data and previous postvoiding EM analysis work [85,86].

**Figure 12.** Impact of early failures on (**a**) voltage drop of a sample grid and (**b**) estimated mesh MTF for ibmpg2 [44].

**Figure 13.** Highlighted wire with high voiding probability on the layout of a metal layer [91].

Meanwhile, Houman Zahedmanesh et al. investigated the EM limits of Cu nanointerconnects by using a novel hybrid physics-based model [93]. The modeling framework incorporates variations of materials, dimensions, interfaces, and operation conditions. It considers void dynamics and resistance shift by using a local cellular automation module with a resistive network. The simulation result only shows that the nucleation phase gets more significant under a narrower linewidth, it but also predicts complex R-shift signatures which match well with experimental data. Sarath Mohanachandran Nair et al. proposed a variation-aware physics-based EM modelling which is experimentally calibrated [94]. The model can be used to explore the impact of material and dimension on design space and to study failure time variation at various operating conditions. Then the model was extended to handle both *early* and *conventional failures* [95]. In [96,97], the system-level simulation on EM under 3 nm technology node was performed. It is shown that Ru rails reduced IR-drop penalty by a factor of ~0.6 than the Cu rails. Although EM voids appear in multiple PDN segments, the EM induced IR-drop always stay below 3.3% without any failed operation of standard cells.

Except for the studies on EM reliability of full power grids, there are other novel works on EM by using physics-based model. In [98–100], the authors explored an approach to enhance the TSV grid reliability. The main idea is to allow the nonfaulty TSVs to be temporarily deactivated so that it can take advantage of EM recovery property. To achieve this goal, a reconfigurable routing network for a (4:2) TSV group was adopted, as depicted in Figure 14. Depending on a recovery schedule, all TSVs cab operate under active mode, recovery mode. The recovery-aware proactive repair approach helps improve EM lifetime of the entire TSV grid by up to 12 times relative to conventional reactive method without an extra area cost.

**Figure 14.** The reconfigurable routing network for a (4:2) TSV group consisting of four f-TSVs and two s-TSVs [100].

#### *3.2. Fast EM Assessment and Optimization*

Although FDM and FEM methods can be applied for transient EM analysis with good accuracy, they are not always preferred when the studied interconnect structures are too large or the computation resource is limited. To overcome this problem, a number of analytical and semi-analytical solutions and speed up techniques have been proposed. At the very beginning, Valeriy Sukharev and Xin Huang et al. derived the analytical equation for transient hydrostatic stress in a confined metal wire [101–103]. The analytical equation captures the impact of time-dependent current density and the stress recovery effect, as shown in Figure 15. The calculation results not only show that in the case of high frequency currents with periods much smaller than the characteristic time of stress evolution, the pulse duty factor decides stress buildup, it but also demonstrates that temperature oscillation can cause notable resistance increase in a short metal line with preexisted voids. On the other hand, the stress recovery effect is applicable to improve on-chip interconnect lifetime by properly managing driving powers at run time.

Later, Hai-Bao Chen et al. developed a first principle based analytical solution for hydrostatic stress evolution in 3 specific interconnect trees [104,105]. It solves stress evolution in multi-branch tree by de-coupling individual branches with suitable boundary conditions which account for interactions between adjacent branches. The solution is based on Laplace transformation technique. Since time-varying temperature and current density and branch length have non-negligible impact on stress evolution, they incorporated these factors into their analytical solution under the same method [106,107]. However, since these analytical solutions only work for specific tree structures, Xiaoyi Wang et al. applied eigenfunction technique for stress evolution in multi-branch trees, and they have extended it to EM analysis in full-chip P/Gs [108–110]. This method handles different current densities and nonuniform thermal distribution. Since this method does not require discretization, except

for its excellent scalability for large-scaled interconnect trees, it brings 10–100 times of computation speed improvement over FDM method. In order to improve EM analysis speed, Liang Chen et al. proposed accelerated separation of variables (ASOV) method which offers improvements over the existing plain SOV-based method [111]. It exhibits 3–5 times of speedup on a number of multi-branch interconnects benchmarks. Furthermore, the SOV-based approach was adopted to obtain a semi-analytical stress transient analysis method which considers the impact of temperature gradient in the studied trees [112]. This method is about an order of magnitude faster than COMSOL with 10× less memory footprint and negligible error loss. It shows that the impact of Joule heating on EM process is significant. Mohammadamir Kavousi et al. studied EM immortality check with considering Joule heating for multi-branch trees [113]. They improved EM stress analysis speed with the benefit of Krylov subspace-based reduction technique which reduces the size of system matrices [114]. Their analysis on interconnect with up to 1000 branches for both void nucleation and growth phases has been accelerated by 28 times. Last, but not least, Mohammad Abdullah Al Shohel et al. found a linear-time approach for immortality check on general tree/mesh interconnects and developed an analytical model of transient stress based on boundary reflections [115,116]. With respect to speedup techniques different from analytical solution, Sandeep Chatterjee et al. presented a fast and scalable methodology for P/G EM verification [117]. The PDE system was converted to a succession of homogeneous linear time invariant (LTI) system. Then the LTI system was solved with an optimized backward differentiation formulas (BDF) solver. Under further help from preconditioned conjugate gradient and parallel programming, this method gives around 23 times of speedup. In [118], a Krylov subspace-based method was proposed for fast stress evolution under finite difference method. After discretization, the original system matrices are reduced so that they can be simulated more efficiently in time domain. This optimized method brings 1–2 orders of magnitude speedup over ordinary finite difference time domain methods.

**Figure 15.** Stress evolution caused by periodic unipolar pulse current densities at cathode end of the metal line under varying temperature [103].

Recently, machine learning techniques have been applied to speed up studies on P/G performance and reliability. In [119], a generative adversarial networks-based (GAN-based) tool (called EM-GAN) was built to do fast analysis on transient stress in multi-branch trees. This work was inspired by the image synthesis feature of GAN. As shown in Figure 16a,

the GAN's inputs include P/G topology, current density distribution at a given aging time. Its output is EM stress distribution. This tool can achieve high prediction accuracy with an average error of 6.6%, and it exhibits 8.3 times speedup over analytic EM solver. In order to achieve even higher accuracy, another graph convolution network-based (GCN-based) tool (called EMGraph) was developed for transient EM stress estimation. Its basic framework is shown in Figure 16b. It shows less than 1.5% average error compared with training data and orders of magnitude faster than COMSOL. Moreover, EMGraph surpasses EM-GAN with 4 times higher prediction accuracy and 14 times faster speed.

**Figure 16.** (**a**) EM-GAN framework for stress estimation and (**b**) framework of EMGraph with multilayer perceptron network [119,120].

In [121], a conditional generative adversarial networks-based (CGAN-based) framework (called *GridNet*) was developed, as shown in Figure 17, to accelerate the incremental full-chip EM-induced IR drop analysis and the optimization for IR drop violation fixing. GridNet provides accurate prediction on IR drop as compared with the ground truth obtained from EMSpice. Since GridNet also provides sensitivity information of node voltage with respect to branch resistance, it expediates localized IR drop violation fixing for P/G design.

tilayer perceptron network [119,120].

(**b**) **Figure 16.** (**a**) EM‐GAN framework for stress estimation and (**b**) framework of EMGraph with mul‐

In [121], a conditional generative adversarial networks‐based (CGAN‐based) frame‐ work (called *GridNet*) was developed, as shown in Figure 17, to accelerate the incremental

GridNet provides accurate prediction on IR drop as compared with the ground truth

**Figure 17.** (**a**) CGAN architecture for EM‐induced voltage prediction and (**b**) comparison of infer‐ ence results from GridNet to EMSpice [121]. **Figure 17.** (**a**) CGAN architecture for EM-induced voltage prediction and (**b**) comparison of inference results from GridNet to EMSpice [121].

Successful P/G design targets at a good enough EM lifetime with a reasonable cost on area. P/G design optimization is an important process to get a balance between EM reliability and area cost. Han Zhou et al. proposed P/G optimization techniques based on fast EM nucleation phase immortality check method for multi-branch interconnect trees [122,123]. They, firstly, verified that the issue can be formulated as a sequence of linear programming problem. Then they proposed an aging-aware optimization method which improves mortal wires' lifetime by adding reservoir branches and allows some interconnects to age/breakdown then just optimizes EM reliability of remaining branches. This strategy ensures the optimization operate effectively. Numerical results demonstrated that the new method can effectively reduce P/G area while ensuring immortality or target lifetime of all the wires. Later, the EM incubation phase immortality check method is introduced into the optimization framework as well [124,125]. The updated method can fix IR drop violation due to EM in minutes for P/Gs from ARM core designs.

#### *3.3. EM Acceleration*

In order to effectively validate the reliability of dual damascene interconnect trees under specific process and structure, it is necessary to detect EM failures in a relatively short testing time. Since the traditional acceleration techniques mainly focus on high temperature and current density, they lead to higher probability of failure due to not only EM but also the other aging mechanisms such as BTI and HCI, which leads to less distinguishability between the mechanisms. In [126–128], several techniques have been proposed to accelerate and decelerate EM failure with the help of reservoir branches and sink branches. Figure 18 shows an active branch with a reservoir branch enabled/disabled and the time-dependent hydrostatic stress at cathode. Obviously, the reservoir branch highly extends void nucleation time at the cathode which indicates much better EM reliability. The impact of reservoir could be disabled by applying a current to it in a reversed direction. And the current density decides how much the main branch's lifetime is suppressed.

**Figure 18.** (**a**) An active interconnect branch with a reservoir branch enabled and disabled and (**b**) Hydrostatic stress at cathode under various scenarios [128].

In [127,128], a sink branch was proposed to accelerate EM lifetime of the main branch. As shown in Figure 19, the void nucleation can be accelerated by an active sink branch because the atomic flux gets improved when the currents over main branch and sink are in the same direction [128]. Then the authors proposed a hybrid structure with both sink and reservoir applied to accelerate EM failure. The proposed structures can achieve desired very short TTF under acceleration mode while the main branch itself has 10+ y lifetime under normal usage mode. These novel structures together with temperature control can further accelerate EM testing about 10<sup>5</sup> time under the 150 ◦C temperature limit.

**Figure 19.** (**a**) An active interconnect branch with a sink branch enabled and disabled and (**b**) hydrostatic stress at cathode with the sink branch enabled and disabled [128].

## **4. Modeling of EM Impact on Interconnects in Cache Memory**

EM not only affects P/G but also worsen the reliability of bitlines which is frequently stressed by currents during read/write operations of cache memory. During read/write operations, the unbalanced currents flowing over interconnects cause voids which ultimately lead to operation failure by introducing large delay. In [129–131], the authors designed a methodology for SRAM EM reliability assessment with considering process variations. The equivalent current distributions over bitlines are calculated with an AFDbased current conversion scheme. As shown in Figure 20, the equivalent current is much different from average value of pulsed DC which indicates different conclusion on EM reliability. The current distribution is combined with process variations including threshold voltage variation, gate length variation, and bitline edge roughness, to evaluate bitlines' reliability by using statistical modelling methodology. The authors adjusted bitline width to find an optimal value to minimize the total probability of failure for an SRAM array and to maximize its yield. The experiment results indicate that a 22 nm technology-based 256 rows × 128 columns SRAM array suffer from serious EM issue if the bitline width is chosen as <sup>1</sup> ⁄<sup>2</sup> metal pitch, as shown in Figure 21. And a tradeoff between functional failure and EM failure can be reached for a 46 nm width bitline when the edge roughness is incorporated.

**Figure 20.** Typical current distribution comparison between time-average value and converted value for all segments within an SRAM bitline [131].

The current distributions over bitlines are decided by activity of cells attached to them. The impact of cell activity on EM reliability of an SRAM cell array was studied in [132]. However, since none of the previous prior works have taken SRAM workload in realistic usage scenarios to evaluate interconnects' reliability, a simulator called CacheEM was proposed for reliability analysis on cache memory aging due EM by considering the realistic application scenarios of cache in an ARM microprocessor [49]. CacheEM includes five parts: microprocessor emulation, memory cell array activity extraction, computation of current in long interconnects, evaluation on time-dependent hydrostatic stress and the resistance shift of interconnects, and characterization of the interconnect EM lifetime distribution of a cache memory.

Figure 22a shows a simplified schematic of an SRAM cell array which is handled in CacheEM with currents corresponding to specific operations (read 0 and1 and write 0 and1) marked. Figure 22b,c show interconnect-array in SRAM cell array and the representative interconnects corresponding to a column of cells [49]. *LStart* and *LEnd* are the segments from the array to the pre-charge and write drivers, respectively. In order to perform timedependent EM assessment on these interconnects, the current density distributions over them need to be determined in the first step. To obtain the current density distributions, it is required to know the number of load and store to each cell, which means the cells' activity in cache memory needs to be extracted. On the other hand, the activity of cells in the cache memory highly depends on the application of the microprocessor where the cache is configured. It needs a simulation flow for CacheEM from a top–down perspective as described in Figure 23. In the first step, several testbenches on the microprocessor which is configured with specific cache memory settings are executed. The trace of the target cache memory is recorded. In the second step, the recorded memory trace is fed into a cache simulator to extract the cells' activity. Then, the variations in effective atomic diffusivity and threshold stress are incorporated into CacheEM with a Monte Carlo simulation to obtain the cache EM lifetime distribution. For each sample, the current flowing in each branch is derived individually with considering the cells' activity and the currents during each operation (read 0/1 and write 0/1).

**Figure 21.** (**a**) EM probability of failure, (**b**) functional probability of failure, and (**c**) yield comparison for 256 rows × 128 columns 22-nm SRAM array w/and w/o considering edge roughness [131].

Figure 24 shows the representative operation distributions on the L1 I-Cache cells in an ARMv8 core while it is running the *sjeng*, *specrand*, and *patricia* benchmarks [49]. The cells' activity and currents during read/write operations under user-defined parameters, including gate length, bit-line unit capacitance, row number of the cells, supply voltage, and temperature are combined together to calculate the effective current density distributions over each BL, BLB, VDD, and GND interconnects. Then they are provided to FDM-based EM solver to evaluate time-dependent stress evolution. The EM lifetime distribution is measured under predefined threshold value. Figure 25 shows the hydrostatic stress distribution of VDDs, GNDs, BLs, and BLBs in randomly selected caches [49]. The current density in the interconnect which suffers the most serious EM stress is plotted on the basis of the left *y*-axis in each subplot. In BL and BLB, the voids are most likely formed in the middle of the lines. The exact void position depends on the ratio between read and write operations. The voids most likely appear at the ends of VDDs and GNDs. However, the voids most likely appear near the via for GNDs, but far from the vias for the VDDs due to the difference of current direction. Stress is highest when current flows into the vias, and lowest when current flows away from the vias. Under a two-port implementation and bi-directional flow, the stress is highest around the middle of the VDDs and at the ends of the GNDs. After the maximum stress reaches *σth*, the voids are generated, and then after the void size becomes larger than the threshold size, the interconnect resistance increases. An interconnect's lifetime is obtained when the resistance shift reaches a pre-defined threshold value.

**Figure 22.** (**a**) Simplified schematic of SRAM cell array with currents relevant with various operation marked. The left and right parts show the currents relevant with read 0 and 1 and write 0 and 1, respectively. The currents for read 0 and write 0 in red, while the currents for read 1 and write 1 are in blue, (**b**) interconnect array in SRAM cell array, with periodic cells shown functional probability of failure, and (**c**) representative interconnects corresponding to a column of cells [49].

**Figure 23.** Simulation flow of CacheEM [49].

*Micromachines* **2022**, *13*, 883 23 of 31

**Figure 24.** (**a**) Read 0, (**b**)read 1, (**c**) write 0, and (**d**) write 1 distributions of L1 I‐Cache cells in ARMv8 core which has run sjeng, specrand, and *patricia* benchmarks [49]. **Figure 24.** (**a**) Read 0, (**b**) read 1, (**c**) write 0, and (**d**) write 1 distributions of L1 I-Cache cells in ARMv8 core which has run sjeng, specrand, and *patricia* benchmarks [49].

(**a**) (**b**) (**c**) (**d**) Figure 26a shows the lifetime distribution of BL, BLB, and GND interconnects which suffer from the most stress in 100 samples of an I-Cache. With the help of the two-ports connection, GND has a much better EM reliability than BL and BLB. BL is the most vulnerable interconnect in the I-Cache. Figure 26b shows the time-dependent resistance shift of the corresponding BLs in 100 cache samples [49]. BL is not always more vulnerable than BLB. Their relative vulnerability might change with the cache configuration. CacheEM was applied to investigate the impact of configuration parameters on performance and EM degradation of an I-Cache. The configuration parameters include cache line size, replacement policy, set associativity, latency, and cache size. Cache line size has options of 8B, 16B, and 32B. The I-Cache hit rate is 94.43%, 93.65%, and 93.75%, when the line size is 8B, 16B, and 32B, respectively. 32B brings better lifetime than the other two options because the interconnect length is much shorter while the current density is smaller. Replacement policy, associativity, and latency do not have obvious impact on I-Cache performance and reliability. Cache size has a non-negligible impact on the hit rate and the lifetime distribution of the studied I-Cache. The hit rate is 89.16% under 8 KB, 93.65% under 16 KB, 97.39% under 32 KB, and 99.19% 64 KB, respectively. 8 KB cache has better lifetime than the other three options because the interconnect length is shorter than the other cases and the current density is smaller. For the remaining three cases, EM reliability improves with the increase of the cache size.

**Figure 25.** The hydrostatic stress distribution before a void appears in (**a**) BLs, (**b**) BLBs, (**c**) VDDs in a one port implementation, (**d**) GNDs in a one port implementation, (**e**) VDDs in a two-port implementation, and (**f**) GNDs in a two-port implementation [49].

**Figure 26.** (**a**) Lifetime distribution of the BL, BLB, and GND interconnects which suffer from the most serious stress in samples of an I-Cache, and (**b**) the time-dependent resistance shift of the corresponding BLs [49].

## **5. Conclusions**

Recent progress in physics-based modeling of EM aging of on-chip interconnects has been reviewed in this paper. Due to the over-conservative EM assessment on multi-branch interconnect tree by conventional Black's equation, physics-based modeling methodologies which consider the atom flow between adjacent branches are required to provide more accurate EM assessment. These modeling methodologies can be achieved with analytical equations or numerical computations. Voltage-based immortality check can simplify EM assessment by eliminating the immortal trees which do not need time-dependent evaluation on stress and resistance. The speedup techniques, such as separation of variable method and Krylov subspace-based reduction technique, are applicable for fast assessment on large-scaled circuits. Machine learning algorithms like GAN also help improve assessment efficiency. With respect to interconnects in cache memory, their EM reliability depends on physical dimensions, cache configurations, and workloads under specific benchmarks.

Since the integration density in SoC is always getting higher to achieve better system performance, the number of interconnect trees to be studied becomes tremendous. It is a challenging task to evaluate the full EM reliability of on-chip interconnects with good accuracy as well as fast speed. Further application of machine learning techniques deserves more attention to accelerate the evaluation speed and to improve the models' prediction capability. On the other hand, since Cu is predicted to be replaced by novel materials such as Co and Ru to ensure interconnects' reliability at sub-10 nm technology nodes, more efforts need to be made to novel modeling methodologies which suitably estimate the aging of new material-based interconnects and networks.

**Author Contributions:** Writing—Original Draft Preparation, R.Z. and W.-S.Z.; Writing—Review and Editing, W.-S.Z. and D.-W.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Zhejiang Provincial Natural Science Foundation under Grants LXR22F040001 and LD22F040003, and the Natural Science Foundation of China (NSFC) under Grants 61874038 and 62101170.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Dongzhe Pan <sup>1</sup> , Bin You 1,\*, Xuan Wen <sup>2</sup> and Xungen Li <sup>1</sup>**


**Abstract:** This article presents a novel wideband bandpass filter based on the integration of a substrate integrated waveguide (SIW) and a spoof surface plasmon polariton (SSPP). An SIW cavity with periodic arrays of meander-slot units is etched on the top metallic layer to achieve the characteristics of a multi-order filter with good performance. The passbands can be flexibly selected by varying the geometric parameters of the SIW and SSPP to adjust the lower and upper sidebands independently. Using a redistribution layer (RDL) process, a novel 3D capacitive interconnection called a throughdielectric capacitor (TDC) is proposed and collaboratively designed with an interdigital capacitor to achieve capacitive source-load cross-coupling. The proposed filter has a center frequency of 60 GHz, with a wide 3-dB fractional bandwidth of about 45.8%. The improved simulated sideband suppression has a 30 dB rejection at 40 GHz and 75.4 GHz, corresponding to a 30-dB rectangular coefficient of 1.28.

**Keywords:** substrate integrated waveguide (SIW); spoof surface plasmon polaritons (SSPPs); integrated passive device; through-dielectric capacitor (TDC)

## **1. Introduction**

With the availability of 60 GHz unlicensed frequency bands and millimeter-wave spectrum for fast data transmission, wireless systems are increasingly operating at higher frequencies. Meanwhile, filters are key front-end components in communication systems, and higher transmission rates must be supported by wider bandwidths and higher frequencies. High-quality and compact BPFs are essential for the improvement of the overall performance of an RF receiver [1]. Acoustic-wave resonators (AWR) such as surface acoustic wave (SAW) [2] and film bulk acoustic resonators (FBAR) [3] exhibit remarkable features in the field of integrated circuit design, with excellent quality factor (Q) and compact size. However, the confined electromechanical coupling coefficient of piezoelectric material limits their bandwidth and center frequency expansion, and starting from C-band, the performance of SAW/FBAR decreases sharply in broadband and high frequency [4–6], making it unsuitable for high frequency bands.

The substrate integrated waveguide (SIW) is popular, with high-pass characteristics, and it exhibits many advantages, such as low insertion loss, high quality factor, and ease of integration. However, its applications in mmW frequency are limited by the harmonics in stopband and bulk size. Frequency selectivity can be improved by introducing transmission zeros (TZs) into the SIW filter. The filter in [7] employs a source-load coupling to improve the slope of the sideband. The conventional folded coupling configuration in [8] introduces a TZ on each side of the passband. An SIW filter with higher-order mode resonators is proposed in order to achieve modal bypass coupling [9]. Multi-order filters with couplings require complicated coupling structures that are larger in size and create some inevitable

**Citation:** Pan, D.; You, B.; Wen, X.; Li, X. Wideband Substrate Integrated Waveguide Chip Filter Using Spoof Surface Plasmon Polariton. *Micromachines* **2022**, *13*, 1195. https://doi.org/10.3390/ mi13081195

Academic Editor: Stefano Mariani

Received: 13 June 2022 Accepted: 27 July 2022 Published: 28 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

losses. Miniaturization is not possible because cascaded multiple resonators prohibit full miniaturization. The use of 3D-stacked structures to achieve miniaturization has also become a popular method in recent years [10–16]; however, it introduces additional losses.

To avoid the deficiencies mentioned above, a new design structure to implement a quasi-elliptic SIW filter is proposed. Coplanar waveguides (CPW) [17] and complementary split-ring resonators (CSRR) [18–20] have been embedded in the SIW cavities to achieve miniaturization, but they result in bad insertion loss. The spoof surface plasmon polariton (SSPP) waveguides that support the surface wave transmission on thin, planar corrugated metals [21,22] have attracted the interest of researchers due to their exceptional properties, such as high confinement [23], low loss [24], potential for solving severe on-chip signal integrity and interference issues [25], as well as the possibility for minimizing the circuit area [26,27]. By embedding the SSPP directly in the SIW, it has been proven to further minimize the filter and improve the transmission performance [28,29]. However, only the upper sideband is improved, due to the characteristics of the SSPPs.

This article proposes a novel SIW-based bandpass filter with periodic arrays of meander-slot etching on the top metallic layer. In [30], an advanced RDL (redistribution layer) process was adopted to design passive components, and it proved to be feasible. In this article, we design the filter using a silicon substrate with 4 RDLs. The dispersion and transmission characteristics are numerically studied by simulation. Due to the distinctive structure of the SSPPs, the asymptotic frequency is lower than that of the original groove structure; only one SIW cavity can achieve the characteristic of multi-order filter, and the size of BPF is greatly diminished. The low-cutoff and high-cutoff frequencies can be flexibly adjusted with the variation of size of the SIW and SSPPs, respectively. With the RDL process, an improved TDC structure is proposed and analyzed. After obtaining the equivalent lump model, a TDC is designed as a part of the coupling circuit to achieve capacitive source-load cross-coupling, and a left transmission zero (TZ) appears successfully. The location of the TZ id can be selected by tuning the TDC and interdigital capacitor. The filter simulation results show that the sideband suppression and stopband suppression are significantly improved, TZ is generated in 37.92 GHz with the rejection of −45 dB, and the 30-dB rectangular coefficient is 1.28.

#### **2. Design of the Proposed BPF**

This work adopts a redistribution layer (RDL) process whose stack-up is illustrated in Figure 1a. All circuit structures are designed in a four-layer dielectric system interconnected by through-dielectric vias (TDV). The material of the dielectric layer is polyimide (PI-HD4100) with a relative permittivity of 3.2. As shown in Figure 1b, two dielectric layers (P3, P4) are utilized as the SIW cavity. The device layer (M4) is employed to etch the SSPPs grooves and M2 is used as a ground plane. The M1 layer is used to design interdigital capacitors. The low-cutoff frequency is determined by the size of the SIW cavity, and the high-cutoff is finally acquired from the dispersion curve of the SSPPs. All electromagnetic (EM) simulations were carried out using the finite-element method (FEM) of the HFSS 3D simulator.

The vertical view of the proposed bandpass filter is shown in Figure 2., An array of periodic curved slots is etched on the top surface metal of the SIW cavity along the direction of electromagnetic wave propagation. The top layer comprises three parts: the GCPW input (region I), the SIW–SSPP transition (region II), and the periodic array part for the propagation of the SSPP. The GCPW transmission line with lower loss and better heat dissipation is more appropriate for mmW integrated circuits; the filter is excited by two 50-Ω GCPW transmission lines with a pair of quarter wavelength coupling slots. Three curved slot cells of different lengths (*Ls*1, *Ls*2, *Ls*3), which are gradually increased in length, constitute region II. The increase is optimized to achieve smooth transition and mode matching between the SIW and SSPP, and then the SSPP mode is effectively excited. The dimensions of the BPF are calculated below.

P2

P3

P4

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 3 of 12

M1

M2

M3

M4

via

RDL (5.4μm )

Device Layer (5.4μm )

(**a**) (**b**)

The dimensions of the BPF are calculated below.

**Figure 1.** (**a**) Stack-up of the RDL process. (**b**) Configuration of the proposed filter (without Si layer). **Figure 1.** (**a**) Stack-up of the RDL process. (**b**) Configuration of the proposed filter (without Si layer).

mode matching between the SIW and SSPP, and then the SSPP mode is effectively excited.

**Figure 2. Figure 2.** Vertical Vertical view of the proposed SIW BPF based on SSPP. view of the proposed SIW BPF based on SSPP.

*W*g1 *W*g2

*<sup>L</sup> <sup>W</sup>*SIW <sup>g</sup> The TE101 mode resonant frequency of the SIW cavity is adopted as the low-cutoff frequency, and the dimensions are determined as in [31]:

$$f\_0 = \frac{c\_0}{2\sqrt{\varepsilon\_r}} \sqrt{\frac{1}{W\_{SIW}^2} + \frac{1}{L\_{SIW}^2}}\tag{1}$$

Vias

SSPPs groove

where *WSIW* and *LSIW* are the width and length of the equivalent rectangular waveguide, calculated as:

$$\mathcal{W}\_{SIW} = w - 1.08 \frac{d^2}{p} + 0.1 \frac{d^2}{w} \tag{2}$$

$$L\_{SIVW} = l - 1.08\frac{d^2}{p} + 0.1\frac{d^2}{l} \tag{3}$$

calculated as:

where *w* and *l* represent the width and length of the SIW cavity, respectively, *d* is the diameter of Cu via, and *p* is the center-to-center pitch between the adjacent via holes. *c*<sup>0</sup> is the light velocity in vacuum, and *ε<sup>r</sup>* is the relative permittivity. In order to make the resonance frequency *f*<sup>0</sup> = 47 GHz, *WSIW* = 2 mm and *LSIW* = 3.8 mm is finally determined by the above equations. is the light velocity in vacuum, and is the relative permittivity. In order to make the resonance frequency <sup>0</sup> = 47 GHz, = 2 mm and = 3.8 mm is finally determined by the above equations. The schematic diagram of the unit cell of the SSPP is illustrated in Figure 3. According to [32], the dispersion curve for the SSPP mode propagated in the metallic grove array can

The TE101 mode resonant frequency of the SIW cavity is adopted as the low-cutoff

1 2

+ 0.1 2 

+ 0.1 2 

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 4 of 12

frequency, and the dimensions are determined as in [31]:

<sup>0</sup> =

0 2√ √ 1 <sup>2</sup> +

= − 1.08

= − 1.08

where and are the width and length of the equivalent rectangular waveguide,

where and represent the width and length of the SIW cavity, respectively, is the diameter of Cu via, and is the center-to-center pitch between the adjacent via holes. <sup>0</sup>

 2 

 2 

The schematic diagram of the unit cell of the SSPP is illustrated in Figure 3. According to [32], the dispersion curve for the SSPP mode propagated in the metallic grove array can be expressed as be expressed as = 0√1 + <sup>2</sup> 2(0) (4)

$$k = k\_0 \sqrt{1 + \frac{W^2}{C^2} \tan^2(k\_0 L)}\tag{4}$$

(1)

(2)

(3)

where *k*<sup>0</sup> = 2*π*/*λ* refers to the propagation constant in free space, *C* is the width of the unit cell, *W* and *L* represent the width and depth of the grooves, respectively. unit cell, and represent the width and depth of the grooves, respectively.

**Figure 3.** The vertical view of (**a**) rectangular groove unit cell, (**b**) meander–slot unit cell, (**c**) SIW with meander-slot unit cell. **Figure 3.** The vertical view of (**a**) rectangular groove unit cell, (**b**) meander–slot unit cell, (**c**) SIW with meander-slot unit cell.

When the width of the groove and the width of the unit cell is fixed, the depth of the groove is the main factor affecting the dispersion. Meanwhile, to keep the occupied area of the groove unchanged, a meander-slot structure is proposed. Figure 3a and b show a traditional groove unit and a meander-slot unit with the same height (). Figure 3c shows the unit cell consisting of an SIW and meander-slot units. The orange and white parts represent the metal surface on the top layer and the slot line, respectively. The geometric parameters are set as shown in Table 1. When the width of the groove *W* and the width of the unit cell *C* is fixed, the depth of the groove *L* is the main factor affecting the dispersion. Meanwhile, to keep the occupied area of the groove unchanged, a meander-slot structure is proposed. Figure 3a,b show a traditional groove unit and a meander-slot unit with the same height (*L*). Figure 3c shows the unit cell consisting of an SIW and meander-slot units. The orange and white parts represent the metal surface on the top layer and the slot line, respectively. The geometric parameters are set as shown in Table 1.

**Table 1.** The geometric parameters of the SSPPs unit cell.


To visualize the relationship between the dimension and frequency of the SSPP, the dispersion curves are calculated using the eigenmode solver of commercial electromagnetic software. In the simulations, the dispersion relation was obtained by calculating the eigenfrequency of the SSPPs unit. Figure 4 shows the dispersion curves of fundamental SSPP modes, *k* represents propagation constant, and *k* is swept from 0◦ to 180◦ between the period boundaries in the propagation direction. It is clear that the two dispersion curves have similar frequency trends, increasing as *k* increases. It is obvious that these two dispersion curves have similar frequency variation trends, but the asymptotic frequency of the groove unit is 97 GHz when the asymptotic frequency of meander-slot unit is 73.7 GHz. This means that adopting a meander-slot SSPP in the same occupied area can reduce the available asymptotic frequency, which is effective in minimizing the device. Figure 4

demonstrates that the SIW can be combined with SSPPs to achieve bandpass characteristics, starting from the cutoff frequency, and its dispersion curve behaves like the SIW in the low frequency range and like SSPPs in the high frequency range. Figure 4 demonstrates that the SIW can be combined with SSPPs to achieve bandpass characteristics, starting from the cutoff frequency, and its dispersion curve behaves like the SIW in the low frequency range and like SSPPs in the high frequency range.

*C L L<sup>s</sup> W WSIW* 0.4 mm 0.39 mm 1 mm 0.075 mm 2 mm

To visualize the relationship between the dimension and frequency of the SSPP, the dispersion curves are calculated using the eigenmode solver of commercial electromagnetic software. In the simulations, the dispersion relation was obtained by calculating the eigenfrequency of the SSPPs unit. Figure 4 shows the dispersion curves of fundamental SSPP modes, represents propagation constant, and is swept from 0° to 180° between the period boundaries in the propagation direction. It is clear that the two dispersion curves have similar frequency trends, increasing as increases. It is obvious that these two dispersion curves have similar frequency variation trends, but the asymptotic frequency of the groove unit is 97 GHz when the asymptotic frequency of meander-slot unit is 73.7 GHz. This means that adopting a meander-slot SSPP in the same occupied area can reduce the available asymptotic frequency, which is effective in minimizing the device.

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 5 of 12

**Table 1.** The geometric parameters of the SSPPs unit cell.

**Figure 4.** Dispersion curve of the unit cell. **Figure 4.** Dispersion curve of the unit cell.

To confirm the passband alteration characteristics of the proposed filter, a parameter inquiry was conducted. The simulated results of transmission coefficients are demonstrated in Figure 5 with different geometric parameters. It can be inferred that the highcutoff frequency of the proposed filter is determined by the SSPPs length , and the lowcutoff frequency is determined by the SIW width . Figure 5a shows that when the value of increases gradually, the upper sideband frequency of the filter will move to the left, while the lower sideband remains constant; the bandwidth becomes narrower accordingly. It can be attributed to the gradual increase of the propagation constant and momentum, and the decrease of the asymptotic frequency, with the increase of the geometric length of the curved slot. In addition, as shown in Figure 5b, the decreases and the lower sideband shifts to the right because the resonant frequency of the SIW is determined by the SIW width . Therefore, the bandwidth of the filter can be flexibly controlled by adjusting the size of the SIW cavity and the SSPP slot. Finally, the dimensions of the BPF are determined by the low-cutoff frequency (47 GHz) and highcutoff frequency (74 GHz), as shown in Table 2. To confirm the passband alteration characteristics of the proposed filter, a parameter inquiry was conducted. The simulated results of transmission coefficients are demonstrated in Figure 5 with different geometric parameters. It can be inferred that the high-cutoff frequency of the proposed filter is determined by the SSPPs length *L<sup>s</sup>* , and the low-cutoff frequency is determined by the SIW width *WSIW*. Figure 5a shows that when the value of *L<sup>S</sup>* increases gradually, the upper sideband frequency of the filter will move to the left, while the lower sideband remains constant; the bandwidth becomes narrower accordingly. It can be attributed to the gradual increase of the propagation constant and momentum, and the decrease of the asymptotic frequency, with the increase of the geometric length *L<sup>s</sup>* of the curved slot. In addition, as shown in Figure 5b, the *WSIW* decreases and the lower sideband shifts to the right because the resonant frequency of the SIW is determined by the SIW width *WSIW*. Therefore, the bandwidth of the filter can be flexibly controlled by adjusting the size of the SIW cavity and the SSPP slot. Finally, the dimensions of the BPF are determined by the low-cutoff frequency (47 GHz) and high-cutoff frequency (74 GHz), as shown in Table 2. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 6 of 12

**Figure 5.** Simulated results of the proposed bandpass filter with different (**a**) SSPPs length . (**b**) SIW width . **Figure 5.** Simulated results of the proposed bandpass filter with different (**a**) SSPPs length *LS*. (**b**) SIW width *WSIW*.

The SIW filter using SSPP without source-load coupling was simulated, and the re-

With the multi-layer stacked structure, 3D interconnections can be applied. To further improve the performance of the filter, source-load cross-coupling was employed. The sectional drawing of the proposed filter with source-load cross-coupling is shown in Figure 7a. With the adoption of the TDC, signals can be transmitted vertically through different layers. There are two microstrip lines and an interdigital capacitor in the M1 layer. The interdigital capacitor is more suitable for applications where low values of capacitance are required. Letting the finger width (*X* = 0.025 mm) equal the slot width to achieve maximum capacitance density, the expression [33] for estimation of capacitance of the in-

> 

where n is the number of fingers and C is in pF. <sup>1</sup> = 0.75 and <sup>2</sup> = 0.175 are deter-

[( − 3)<sup>1</sup> + <sup>2</sup>

refer to in [33], where *T* is the height of the substrate ( = 0.01

] (5)

*LSIW WSIW L<sup>g</sup> Wg1 Wg2 Ls1 Ls2 Ls3* 4 mm 2 mm 0.76 mm 0.055 mm 0.2 mm 0.26 mm 0.49 mm 0.78 mm

**<sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> <sup>90</sup> −80 −70 −60 −50 −40 −30 −20 −10**

**Frequency (GHz)**

**Figure 6.** Simulation results of the BPF based on SIW and SSPP.

GHz to 76.5 GHz; and the corresponding rectangular coefficient is 1.4.

 S11 S21

= ( + 1)

**Table 2.** The dimensions of the BPF.

**0**

terdigital capacitor can be given by

mined by the value of

**S-parameter (dB)**

**Table 2.** The dimensions of the BPF. **Table 2.** The dimensions of the BPF.

SIW width .

**0**

**S21 (dB)**

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 6 of 12

**<sup>20</sup> <sup>25</sup> <sup>30</sup> <sup>35</sup> <sup>40</sup> <sup>45</sup> <sup>50</sup> <sup>55</sup> <sup>60</sup> <sup>65</sup> <sup>70</sup> <sup>75</sup> <sup>80</sup> <sup>85</sup> <sup>90</sup> <sup>95</sup> −60 −40 −20**

*L***<sup>S</sup> = 1.9** *L***<sup>S</sup>**

*L***<sup>S</sup>**

*L***<sup>S</sup> = 2** 

 **= 2.1**

 **= 1.8**

**Frequency (GHz)**)


(**a**) (**b**)

**Figure 5.** Simulated results of the proposed bandpass filter with different (**a**) SSPPs length

**S21 (dB)**

**<sup>20</sup> <sup>25</sup> <sup>30</sup> <sup>35</sup> <sup>40</sup> <sup>45</sup> <sup>50</sup> <sup>55</sup> <sup>60</sup> <sup>65</sup> <sup>70</sup> <sup>75</sup> <sup>80</sup> <sup>85</sup> <sup>90</sup> <sup>95</sup> −60 −40 −20**

*WSIW* **= 2** *WSIW* **= 1.9** 

*WSIW* **= 2.1**

*WSIW* **= 1.8**

. (**b**)

**Frequency (GHz)**)

**0**

The SIW filter using SSPP without source-load coupling was simulated, and the result is shown in Figure 6. Good bandpass features and high-efficiency propagation are obtained. The BPF operates from 46.13 GHz to 74 GHz; the 30-dB bandwidth is from 37.5 GHz to 76.5 GHz; and the corresponding rectangular coefficient is 1.4. The SIW filter using SSPP without source-load coupling was simulated, and the result is shown in Figure 6. Good bandpass features and high-efficiency propagation are obtained. The BPF operates from 46.13 GHz to 74 GHz; the 30-dB bandwidth is from 37.5 GHz to 76.5 GHz; and the corresponding rectangular coefficient is 1.4.

**Figure 6.** Simulation results of the BPF based on SIW and SSPP. **Figure 6.** Simulation results of the BPF based on SIW and SSPP.

With the multi-layer stacked structure, 3D interconnections can be applied. To further improve the performance of the filter, source-load cross-coupling was employed. The sectional drawing of the proposed filter with source-load cross-coupling is shown in Figure 7a. With the adoption of the TDC, signals can be transmitted vertically through different layers. There are two microstrip lines and an interdigital capacitor in the M1 layer. The interdigital capacitor is more suitable for applications where low values of capacitance are required. Letting the finger width (*X* = 0.025 mm) equal the slot width to achieve maximum capacitance density, the expression [33] for estimation of capacitance of the interdigital capacitor can be given by With the multi-layer stacked structure, 3D interconnections can be applied. To further improve the performance of the filter, source-load cross-coupling was employed. The sectional drawing of the proposed filter with source-load cross-coupling is shown in Figure 7a. With the adoption of the TDC, signals can be transmitted vertically through different layers. There are two microstrip lines and an interdigital capacitor in the M1 layer. The interdigital capacitor is more suitable for applications where low values of capacitance are required. Letting the finger width (*X* = 0.025 mm) equal the slot width to achieve maximum capacitance density, the expression [33] for estimation of capacitance of the interdigital capacitor can be given by

$$\mathbf{C}\_{l} = (\varepsilon\_{r} + 1) \frac{L\_{l}}{\mathcal{W}\_{l}} [(n - 3)A\_{1} + A\_{2}] \tag{5}$$

where n is the number of fingers and C is in pF. <sup>1</sup> = 0.75 and <sup>2</sup> = 0.175 are determined by the value of refer to in [33], where *T* is the height of the substrate ( = 0.01 where *n* is the number of fingers and *C* is in pF. *A*<sup>1</sup> = 0.75 and *A*<sup>2</sup> = 0.175 are determined by the value of *<sup>T</sup> X* refer to in [33], where *T* is the height of the substrate (*T* = 0.01 mm). As shown in Figure 7b, the complete coupling circuit consists of two TDCs, two microstrips and an interdigital capacitor which are lumped models. The collaborative design method of the TDC and the coupling circuit can be used to minimize the influence of the interconnection structure. The series connection of the interdigital capacitor and TDC can also reduce the circuit capacitance and increase the adjustment range of the interdigital capacitor. The admittance of the equivalent coupling circuit with the TDC in Figure 7b can be expressed as

$$Z = Z\_t + R\_m + j(\omega L\_m - \frac{1}{\omega \mathbb{C}\_i}) \tag{6}$$

where *Z<sup>t</sup>* is the impedance of the TDC structure, which is analyzed below. 7b can be expressed as

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 7 of 12

where

Gap

coupling circuit (**left**).

TDV

7b can be expressed as

where

mm). As shown in Figure 7b, the complete coupling circuit consists of two TDCs, two microstrips and an interdigital capacitor which are lumped models. The collaborative design method of the TDC and the coupling circuit can be used to minimize the influence of the interconnection structure. The series connection of the interdigital capacitor and TDC can also reduce the circuit capacitance and increase the adjustment range of the interdigital capacitor. The admittance of the equivalent coupling circuit with the TDC in Figure

= + + ( −

is the impedance of the TDC structure, which is analyzed below.

mm). As shown in Figure 7b, the complete coupling circuit consists of two TDCs, two microstrips and an interdigital capacitor which are lumped models. The collaborative design method of the TDC and the coupling circuit can be used to minimize the influence of the interconnection structure. The series connection of the interdigital capacitor and TDC can also reduce the circuit capacitance and increase the adjustment range of the interdigital capacitor. The admittance of the equivalent coupling circuit with the TDC in Figure

= + + ( −

is the impedance of the TDC structure, which is analyzed below.

*Li*

1 

1 

) (6)

*T*

) (6)

**Figure 7.** (**a**) Cross-section schematic view of the filter; (**b**) Equivalent lump model of this proposed coupling circuit (**left**). **Figure 7.** (**a**) Cross-section schematic view of the filter; (**b**) Equivalent lump model of this proposed coupling circuit (**left**). Figure 8a shows the structure and the main parasitic components of the TDC in the

Figure 8a shows the structure and the main parasitic components of the TDC in the designed filter. Between the two metal layers, two polyimide bonding layers act as capacitive coupling media, while one is filled with copper through the dielectric via. The equivalent lump model of the TDC can be referred to [34], as shown in Figure 8b. Figure 8a shows the structure and the main parasitic components of the TDC in the designed filter. Between the two metal layers, two polyimide bonding layers act as capacitive coupling media, while one is filled with copper through the dielectric via. The equivalent lump model of the TDC can be referred to [34], as shown in Figure 8b. designed filter. Between the two metal layers, two polyimide bonding layers act as capacitive coupling media, while one is filled with copper through the dielectric via. The equivalent lump model of the TDC can be referred to [34], as shown in Figure 8b.

The impedance *Z<sup>t</sup>* can be approximately expressed by

$$Z\_t = \frac{2}{Y\_1} + j\omega L + R \tag{7}$$

$$Y\_1 = G + j\omega \mathcal{C} \tag{8}$$

where *Y*<sup>1</sup> is the admittance of the coupling medium, *R* and *L* represent via resistance and inductance, respectively, when *C* and *G* represent the parasitic capacitance and conductance, respectively.

To accurately verify the compatibility and analyze the electrical performance of the TDC, Ansys Q3D was used to extract all the parasitic components of TDC. Substituting where <sup>1</sup>

ance, respectively.

all the parasitic components into the proposed lump model then enables the transmission characteristics to be obtained by the ADS. The comparison of transmission characteristics between the TDC simulation result by HFSS and the proposed lumped model is shown in Figure 9. Two experimental cases with different heights of through dielectric via (*ht*) are used in simulation. The lumped models of the two cases both match well with the HFSS simulations. The insertion loss S<sup>21</sup> of TDC increases with frequency below 10 GHz, which implies that the coupling behavior of the TDC is mainly capacitive coupling, and slightly decreases in high frequency. The capacitive TDC effect is dominant in the lower frequency range, and the inductive TDV effect becomes dominant in the higher frequency range [34]; the inductance of the TSV channel starts to affect the insertion loss over 10 GHz. The above analysis shows that the proposed TDC can be regarded as a CGRL lumped model, collaboratively designed with the coupling circuit. In addition, it will not seriously influence the circuit during the working frequency. all the parasitic components into the proposed lump model then enables the transmission characteristics to be obtained by the ADS. The comparison of transmission characteristics between the TDC simulation result by HFSS and the proposed lumped model is shown in Figure 9. Two experimental cases with different heights of through dielectric via (ℎ ) are used in simulation. The lumped models of the two cases both match well with the HFSS simulations. The insertion loss S<sup>21</sup> of TDC increases with frequency below 10 GHz, which implies that the coupling behavior of the TDC is mainly capacitive coupling, and slightly decreases in high frequency. The capacitive TDC effect is dominant in the lower frequency range, and the inductive TDV effect becomes dominant in the higher frequency range [34]; the inductance of the TSV channel starts to affect the insertion loss over 10 GHz. The above analysis shows that the proposed TDC can be regarded as a CGRL lumped model, collaboratively designed with the coupling circuit. In addition, it will not seriously influence the circuit during the working frequency.

is the admittance of the coupling medium, *R* and *L* represent via resistance and

inductance, respectively, when *C* and *G* represent the parasitic capacitance and conduct-

To accurately verify the compatibility and analyze the electrical performance of the TDC, Ansys Q3D was used to extract all the parasitic components of TDC. Substituting

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 8 of 12

The impedance can be approximately expressed by

 = 2 1

+ + (7)

<sup>1</sup> = + (8)

**Figure 9.** Comparison of S parameters between TDC and model with different ℎ : (**a**) ℎ = 0.024 mm; (**b**) ℎ = 0.08 mm. **Figure 9.** Comparison of S parameters between TDC and model with different *h<sup>t</sup>* : (**a**) *h<sup>t</sup>* = 0.024 mm; (**b**) *h<sup>t</sup>* = 0.08 mm.

To further verify the lumped model of the whole coupling circuit, the S-parameter is extracted from Figure 6 and simulated with the above coupling circuit model, and the result is compared with the HFSS simulation of the proposed filter, as shown in Figure 10. The frequency responses of the equivalent circuit and filter simulations both exhibit a transmission zero at 37.6 GHz, and their simulation results are in good agreement. To further verify the lumped model of the whole coupling circuit, the S-parameter is extracted from Figure 6 and simulated with the above coupling circuit model, and the result is compared with the HFSS simulation of the proposed filter, as shown in Figure 10. The frequency responses of the equivalent circuit and filter simulations both exhibit a transmission zero at 37.6 GHz, and their simulation results are in good agreement. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 9 of 12

**Figure 10.** Comparison of equivalent circuit and filter simulations. **Figure 10.** Comparison of equivalent circuit and filter simulations.

According to the above analysis and verification, the proposed TDC structure can be applied to transmit signals vertically through different layers, and it is possible to design a circuit with a lumped model. With the adoption of the previous lumped model, capaci-According to the above analysis and verification, the proposed TDC structure can be applied to transmit signals vertically through different layers, and it is possible to design a circuit with a lumped model. With the adoption of the previous lumped model, capacitive source-load cross-coupling is investigated.

capacitance of the interdigital capacitor, the location of the TZ is moved, as shown in Fig-

performance of the filter, a transmission zero is finally generated at 37.92 GHz, and the height of the through dielectric via ℎ = 0.0246 mm and the radius of the via = 0.04

The simulated results of the S-parameters with reflection coefficients (S11) and transmission coefficients (S21) are demonstrated in Figure 12. An SIW cavity is employed to provide multi-transmission poles in order to achieve a passband. The compact BPF operates at 46.1–73.7 GHz with a wide 3dB FBW of 45.8%. Good bandpass features and highefficiency propagation are obtained; the minimum insertion loss is 1.08 dB, and the return loss is better than −10 dB in the passband. Due to the steep upper sideband, the 30-dB bandwidth is from 40 GHz to 75.4 GHz, and the corresponding rectangular coefficient is 1.28 when the BPF without TZ is 1.4. Meanwhile, the stopband rejection is better than 30

It is notable that the proposed BPF has the merits of competitive wideband performance, broad stopband, and compact size. The performance of the proposed filter is compared

with other work in Table 3; all the data are based on the simulation results.

). The final core size of the BPF is 2 mm × 4.4 mm (0.74 × 1.63).

. By adjusting the

.

The capacitance of the interdigital capacitor will change over

**<sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> <sup>90</sup> −80 −70 −60 −50 −40 −30 −20 −10**

*L*i=0.54mm *L*i=0.5mm *L*i=0.46mm

**Figure 11.** Simulated results with different interdigital capacitor length

*L*i=0.58mm

**Frequency (GHz)**

tive source-load cross-coupling is investigated.

mm are finally determined.

**3. Results and Analysis**

dB up to 125 GHz (2

**0**

**S-parameter (dB)**

The capacitance of the interdigital capacitor will change over *L<sup>i</sup>* . By adjusting the capacitance of the interdigital capacitor, the location of the TZ is moved, as shown in Figure 11, and the suppression deteriorates as TZ approaches the passband. To optimize the performance of the filter, a transmission zero is finally generated at 37.92 GHz, and the height of the through dielectric via *h<sup>t</sup>* = 0.0246 mm and the radius of the via *r<sup>t</sup>* = 0.04 mm are finally determined. capacitance of the interdigital capacitor, the location of the TZ is moved, as shown in Figure 11, and the suppression deteriorates as TZ approaches the passband. To optimize the performance of the filter, a transmission zero is finally generated at 37.92 GHz, and the height of the through dielectric via ℎ = 0.0246 mm and the radius of the via = 0.04 mm are finally determined.

The capacitance of the interdigital capacitor will change over

According to the above analysis and verification, the proposed TDC structure can be applied to transmit signals vertically through different layers, and it is possible to design a circuit with a lumped model. With the adoption of the previous lumped model, capaci-

. By adjusting the

.

 

**FBW (%)**

), the location of

**Figure 11.** Simulated results with different interdigital capacitor length **Figure 11.** Simulated results with different interdigital capacitor length *L<sup>i</sup>* .

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 9 of 12

**<sup>20</sup> <sup>30</sup> <sup>40</sup> <sup>50</sup> <sup>60</sup> <sup>70</sup> <sup>80</sup> <sup>90</sup> −80 −70 −60 −50 −40 −30 −20 −10**

<sup>S</sup><sup>11</sup> <sup>S</sup><sup>21</sup>

 HFSS simulation Equivalent circuit TZ

**Frequency (GHz)**

tive source-load cross-coupling is investigated.

**Figure 10.** Comparison of equivalent circuit and filter simulations.

**0**

**S-parameter (dB)**

#### **3. Results and Analysis**

**3. Results and Analysis** The simulated results of the S-parameters with reflection coefficients (S11) and transmission coefficients (S21) are demonstrated in Figure 12. An SIW cavity is employed to provide multi-transmission poles in order to achieve a passband. The compact BPF operates at 46.1–73.7 GHz with a wide 3dB FBW of 45.8%. Good bandpass features and highefficiency propagation are obtained; the minimum insertion loss is 1.08 dB, and the return loss is better than −10 dB in the passband. Due to the steep upper sideband, the 30-dB bandwidth is from 40 GHz to 75.4 GHz, and the corresponding rectangular coefficient is 1.28 when the BPF without TZ is 1.4. Meanwhile, the stopband rejection is better than 30 dB up to 125 GHz (2 ). The final core size of the BPF is 2 mm × 4.4 mm (0.74 × 1.63). It is notable that the proposed BPF has the merits of competitive wideband performance, The simulated results of the S-parameters with reflection coefficients (S11) and transmission coefficients (S21) are demonstrated in Figure 12. An SIW cavity is employed to provide multi-transmission poles in order to achieve a passband. The compact BPF operates at 46.1–73.7 GHz with a wide 3 dB FBW of 45.8%. Good bandpass features and high-efficiency propagation are obtained; the minimum insertion loss is 1.08 dB, and the return loss is better than −10 dB in the passband. Due to the steep upper sideband, the 30-dB bandwidth is from 40 GHz to 75.4 GHz, and the corresponding rectangular coefficient is 1.28 when the BPF without TZ is 1.4. Meanwhile, the stopband rejection is better than 30 dB up to 125 GHz (2 *fc*). The final core size of the BPF is 2 mm × 4.4 mm (0.74*λ<sup>g</sup>* × 1.63*λg*). It is notable that the proposed BPF has the merits of competitive wideband performance, broad stopband, and compact size. The performance of the proposed filter is compared with other work in Table 3; all the data are based on the simulation results. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 10 of 12

on the simulation results).

**4. Conclusions**

of the metal via (ℎ

**Ref.**

**Figure 12.** Simulation result of the bandpass filter with source-load coupling. **Figure 12.** Simulation result of the bandpass filter with source-load coupling.

**Table 3.** Performance comparisons of BPFs operating above millimeter-wave. (All the data are based

In this article, an SIW-based 46.1 GHz to 73.7 GHz bandpass filter is proposed, which is designed on the RDL process. The simulation results show low insertion loss, good frequency selectivity, and wideband harmonic suppression. The etching of periodic arrays of meander slot units on the top metallic layer enables the substrate integrated waveguide (SIW) to attain bandpass characteristics with high-efficiency and strongly confined microwave SSPP transmission. In addition, its asymptotic frequency can be significantly reduced compared to the conventional groove SSPP. This means that the propagation of the structure can occupy a smaller area with the benefit of lower cost, especially for the process of integrating a passive device, and the leakage loss will decrease as the gap area reduces. The simulated results show that the bandwidth of the proposed filter can be flexibly elected by tuning the geometric parameters of the SIW and SSPPs. To improve the performance of the upper sideband, a novel 3D capacitive interconnection is proposed and investigated, the lump model of the TDC is obtained, and a collaborative design with interdigital capacitance is adopted to achieve a transmission zero. By adjusting the height

) and the capacitance of the interdigital capacitor (

the TZ can be selected. The proposed filter based on one SIW resonator is the same as the multi-order filter with coupling, realizing the miniaturization with good performance.

**Author Contributions:** Conceptualization, D.P. and B.Y.; Data curation, D.P., X.W. and X.L.; Formal analysis, D.P. and X.W.; Funding acquisition, B.Y.; Methodology, D.P., X.W. and X.L.; Validation, B.Y. and X.L.; Writing—original draft, D.P.; Writing—review & editing, D.P., B.Y. and X.W. All au-

**Funding:** This research was funded by the Natural Science Foundation of Zhejiang Province under Grant LZ22F010006 and National Natural Science Foundation of China under Grant 61671195.

thors have read and agreed to the published version of the manuscript.

**Conflicts of Interest:** The authors declare no conflicts of interest.

[7]-2 9.1 0.218 × 0.218 0.84 2.5 19.8 [9] 93 2.31 × 1.57 4 1.59 3.5 [29] 236.5 1.16 × 0.37 2 1.86 12 This work 60 1.63 × 0.74 1.08 1.28 45.8

**(GHz) Size (** × **) MIN. IL (dB)**


**Table 3.** Performance comparisons of BPFs operating above millimeter-wave. (All the data are based on the simulation results).

## **4. Conclusions**

In this article, an SIW-based 46.1 GHz to 73.7 GHz bandpass filter is proposed, which is designed on the RDL process. The simulation results show low insertion loss, good frequency selectivity, and wideband harmonic suppression. The etching of periodic arrays of meander slot units on the top metallic layer enables the substrate integrated waveguide (SIW) to attain bandpass characteristics with high-efficiency and strongly confined microwave SSPP transmission. In addition, its asymptotic frequency can be significantly reduced compared to the conventional groove SSPP. This means that the propagation of the structure can occupy a smaller area with the benefit of lower cost, especially for the process of integrating a passive device, and the leakage loss will decrease as the gap area reduces. The simulated results show that the bandwidth of the proposed filter can be flexibly elected by tuning the geometric parameters of the SIW and SSPPs. To improve the performance of the upper sideband, a novel 3D capacitive interconnection is proposed and investigated, the lump model of the TDC is obtained, and a collaborative design with interdigital capacitance is adopted to achieve a transmission zero. By adjusting the height of the metal via (*ht*) and the capacitance of the interdigital capacitor (*C<sup>i</sup>* ), the location of the TZ can be selected. The proposed filter based on one SIW resonator is the same as the multi-order filter with coupling, realizing the miniaturization with good performance.

**Author Contributions:** Conceptualization, D.P. and B.Y.; Data curation, D.P., X.W. and X.L.; Formal analysis, D.P. and X.W.; Funding acquisition, B.Y.; Methodology, D.P., X.W. and X.L.; Validation, B.Y. and X.L.; Writing—original draft, D.P.; Writing—review & editing, D.P., B.Y. and X.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Natural Science Foundation of Zhejiang Province under Grant LZ22F010006 and National Natural Science Foundation of China under Grant 61671195.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

