**4. The Absorption, Reflection, and Transmission of the Top Silicon Optical Window** *4.1. Modeling*

The TMOS cap is fabricated by two layers: a silicon wafer of the order of 100 µm and a much thinner silicon dioxide layer. The optical wave transmission of this Fabry–Perot cavity can be analytically calculated with transfer matrices of the electric fields [23].

Figure 3 exhibits the absorption, reflection, and transmission of the top optical window.

**Figure 3.** The modeled absorption, reflection, and transmission of the optical window as a function of the wavelength. The silicon wafer size is ~100 µm and the oxide layer is of the order of 0.1 µm.

#### *4.2. LUMERICAL Simulation*

A FDTD model has been created in LUMERICAL as described in Figure 4 in order to simulate the TMOS optical window.

**Figure 4.** The simulation model and setup. The simulated area is outlined by the orange rectangle. The lateral dimensions are in the order of 0.1 µm. PML B.C.: perfectly matched layer boundary conditions.

As can be seen in Figure 4, on the sides of the FDTD region periodic boundary conditions were applied and on the top and bottom of the FDTD region a perfectly matched layer (PML) boundary conditions were applied. A radiation source with the wavelength range of 5–20 µm was set above the TMOS cap. The refractive indices of the materials were set by the software using Palik data [20,24].

The simulation results are shown in Figure 5.

**Figure 5.** Simulations results for the transmission of the TMOS optical window as a function of the wavelength. PSA: partial spectral averaging (δ: delta, see Section 4.3).

**Figure 6.** Simulation and modeling of transmission results as a function of the wavelength.

The measured transmission results, compared with the simulation, as shown in Figure 7, are very different and do not show the resonance lines. A good fit between measurements and simulation or modeling requires partial spectral averaging, and δ (delta) estimation, as explained below.

**Figure 7.** The simulated optical transmittance of the TMOS optical window for different silicon dioxide thickness as a function of the wavelength. The simulated results are in good agreement with measured results.

#### *4.3. Comparison with Measurements*

The optical window that was simulated was a Fabry–Perot cavity with two dielectric slab devices. One of the slabs was much wider (>100 µm) than the typical wavelength used for the measurements. Therefore, inside this slab, we expected to find a standing wave. During measurement, the transmitted power through the device was measured for a range of wavelengths. The simulation calculates the transmission for a specific single wavelength, whereas measurement is experimentally conducted around a certain bandwidth, since it is difficult to produce a monochromatic source. Therefore, the calculated simulation should be averaged by a convolution with a Lorentzian weighting function, which is a function mainly used to characterize narrow spectrum lines, such as emission by atoms, laser radiation, etc. [20].

The Lorentzian weighting function is defined by:

$$\left| h(\omega, \omega') \right|^2 = \frac{\delta}{\left( \omega - \omega' \right) + \left( \pi \delta \right)^2} \tag{1}$$

where *δ* is the FWHM (full width at half maximum) and *ω* is the Lorentzian center frequency. The value of *δ* can be evaluated by differentiating the relation λ = c/f with respect to f, resulting in |df| = f/λ dλ, where dλ is the measurement error of the wavelength. The value of df can be used for the evaluation of the FWHM. As described in Figure 5, the value of the FWHM used in the simulation was *δ* = 0.6 (THz), which corresponds to a value of dλ = 200 (nm). This value yields a match between the measurements and simulation when using that value.

#### **5. Optimized Absorption of the TMOS Sensor with an Impedance Matching Layer with the Right Thickness and Location**

To gain an understanding of the electromagnetic absorption of the TMOS sensor, the classical propagation of waves should be applied. The dielectric layers of the TMOS released transistor (see Figure 1c) are comprised of silicon dioxide and silicon nitride. The silicon oxide bond absorbs at ~9.3 µm, while the silicon nitride bond absorbs at ~11.5–12 µm [25]. The physical description becomes complicated since there are several built-in Fabry–Perot resonators and the LUMERICAL simulations become very useful here. Furthermore, it is well established from bolometers that an impedance matching layer with the correct thickness and placement allows optical absorption with 90% efficiency along large arrays. The concept of impedance matching layer for an ideal λ/4 optical cavity is described in Appendix A [26]. The transmission line description, discussed in Appendix B, provides an intuitive description of an ideal optical cavity.

It is apparent that in order to reduce reflection and enhance absorption, an impedance matching layer is mandatory. This layer should provide the required impedance at the frequency of the incident radiation. It should be based on a standard metallization available in the CMOS FAB. Luckily, Ti and TiN may provide this layer. At DC, the specific resistivity of these layers is too low, but at the IR frequency of interest (around 10<sup>13</sup> Hz) the values increase by a factor of ~3 because of the plasma effect [25].

In the TMOS sensor case, the modeling is more complicated than that of Appendix B. We assume that a thin "impedance matching layer" is sandwiched between two optical materials (the interlevel dielectric material of the TMOS). When the skin depth into the metal is larger than its thickness, the metallic film can be considered as a regular optical layer (that may be described by the Drude model [25]). In this case, the Salisbury screen becomes a dielectrically-coated Salisbury sheet (DSS) [26–29] (see an example in Figure 8):

**Figure 8.** Schematic of an optical system consisting of an optical material (1) with n<sup>1</sup> , thin metal layer (M), and optical material (2) with n<sup>2</sup> .

As discussed above, the definition of "impedance matching layer" in this case is much more complicated: What should be the thickness and the location of this layer in the cross section of Figure 1c? Theoretically, Fabry–Perot or transmission line modeling may be performed, but LUMERICAL simulations are very effective and useful here.

At this point, it is important to remind the reader that the accuracy of the LUMERICAL simulation is very much dependent upon the mesh dimensions. Accuracy increases if the mesh is less than one tenth of the wavelength. Furthermore, a simple "sanity check" should be performed to validate the results. The sum of absorption, reflection, and transmission should be 1. If it is higher than 1, the mesh should be redefined and be reduced in dimensions.

We present simulation results below that demonstrate how the thickness of the layer (Section 5.1) and location (Section 5.2) affect absorption. At certain bandpass regions, high absorption of 90% may be obtain (Section 5.3).

#### *5.1. The Effect of the Thickness of the TiN*

To begin with, we wanted to examine the influence of the thickness of the TiN layer on the absorbance. We simulated the model for varying thickness from 10 nm up to 30 nm, all positioned in the same height of 2.16 µm from the bottom of the Box. The absorption, reflection, and transmission of this simulation can be seen in Figure 9.

**Figure 9.** Absorption (A), reflection (R), and transmission (T) of the transistor with varying thicknesses of the TiN layer. Simulation is performed at 9.5 µm wavelength.

It can be seen from Figure 9 that the absorption in the transistor is more significant as the layer thickness decreases and can be improved from 36% up to 52%.

#### *5.2. The Effect of the Location of the TiN along the Upper Silicon Oxide*

After optimizing the thickness of the TiN layer, we continued with optimizing the absorbance of the sensor by changing the location of the TiN along the upper silicon oxide, which starts after the last metal layer at the height in the order of ~1.5 µm (see Figure 10).

**Figure 10.** Schematic of the TiN location in the TMOS sensor.

The absorption, reflection, and transmission were simulated for different locations along the silicone oxide for 10-nm and 20-nm layers of TiN.

It can be seen in Figure 11 that the absorbance is higher at a lower TiN placement. As shown in Section 5.1, the higher absorption is reached by the 10–20 nm TiN layer and is slightly higher for 20 nm layer of TiN.

**Figure 11.** Absorption (A), reflection (R), and transmission (T) of the transistor with varying locations of the TiN layer along the upper silicon oxide for 10-nm and 20-nm layer thickness. This simulation is performed at a 9.5-µm wavelength.

#### *5.3. Simulation of Absorption at an Optical Bandpass*

In practice, the TMOS sensor is provided with a bandpass filter that is determined by the use case. We present the simulation below within a bandpass of 6–13 µm. The simulations were performed with default continuous wave normalization, where in this state the power monitors are normalized by the Fourier transform of the source pulse. In other words, it returns the impulse response of the system. For most applications this is the best choice, but in some applications there may be critical mismatch. For this paper, there is a good agreement between the field's results yield by the normalized power over the source spectrum to the impulse response of single frequency simulations, as can be seen in Figure 12.

**Figure 12.** Absorption of the TMOS transistor with a TiN layer thickness of 20 nm located at 1.58 µm (see Figure 10) as a function of the wavelength.

For the optimized thickness and location of the TiN layer, an absorption of about 90% can be achieved. The simulation results for various TiN layer thickness and location can be seen in Figure 13.

**Figure 13.** Absorption of the TMOS transistor for various thicknesses (T) and locations (L) as a function of the wavelength.

Figure 14 shows the simulation results for additional SiN layer with thickness of 0.1 µm for various locations. Though the maximum absorption decreases by about 10% around 10 µm for most of the bandpass spectrum, the absorption increases by about 20%.

**Figure 14.** Absorption of the TMOS transistor with additional SiN layer with a thickness of 0.1 µm for various locations (L) as a function of the wavelength.

#### **6. Summary**

The goal of this study was to provide physical insight and an intuitive approach to the layers and mechanisms that play primary roles regarding 3D, vacuum-packaged, nano-machined TMOS electromagnetic absorption.

The main parameters that determine the overall absorption efficiency of the vacuumpackaged TMOS sensor are (i) the *transmission* of the optical window and (ii) the *absorption* of the released TMOS sensor, which should be as high as possible. The latter is significantly increased by the impedance matching layer (a dielectric-coated Salisbury sheet) [27–29].

As can be seen in Figure 15, the degrading effect of the oxide, which is part of the optical window, upon the transmission of the electromagnetic radiation in the bandpass of interest (>5 µm) is clearly observed. The average transmission is ~55% while at the resonance wavelengths (around 9.3 µm) it is significantly reduced. By integrating an antireflection (AR) coating on the optical window, its transmission may be significantly increased.

**Figure 15.** The effect of the thickness of the oxide on the transmission of the optical window. The simulation is achieved with PSA averaging and δ = 0.6 (THz).

The absorption of the TMOS sensor may be optimized with the LUMERICAL simulations. By optimizing the thickness of the *impedance matching layer* of a thin TiN layer and its location along the upper inter-level dielectric of the TMOS, a significant absorption of 90% may be achieved. Adding more than one impedance matching layer may also enhance absorption [27] (chapter 9). Physical insight may be gained by applying transmission lines modeling and Smith chart modeling, but this is beyond the scope of this paper [27–29].

Moreover, the properties of thin layers are different from the bulk materials and depend on the deposition technology. Appendix C reports the values used in this study. Better simulations may be achieved if the actual TMOS main dielectric layers and refractive index are measured. The values of n(λ) and k(λ) play a role in determining the exact values of absorption. For example, the refractive index for the main TMOS material SiO<sup>2</sup> used in the simulation can be seen in Figure 16.

**Figure 16.** The refractive index of SiO<sup>2</sup> as a function of the TMOS bandpass wavelength where (**n**) is the real part and (**k**) is the imaginary part.

In summary, the LUMERICAL software package is an excellent tool; however, like all software, it requires understanding of the limitations of the tool and the underlying approximation of results. A designer that is also an expert in transmission line theory and Smith charts may achieve very high absorption of the order of 90% within the bandpass of interest.

**Author Contributions:** Conceptualization, Y.N. and G.G.; methodology, Y.N., G.C. and M.A.; software, G.C., M.A. and S.B.-L.; validation, G.C., M.A. and S.B.-L.; formal analysis, G.C. and M.A.; investigation, G.C. and M.A.; resources, Y.N.; data and M.A.; writing—review and editing, Y.N.; visualization, G.C. and M.A.; supervision, Y.N. and G.G.; project administration, Y.N.; funding acquisition, Y.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by TODOS technologies, grant number 2024790.

**Acknowledgments:** We are grateful for the technical support provided by Greg Baethge, ANSYS LUMERICAL. The authors highly appreciate the contributions of the research team at Technion as well as at STM involved with the TMOS project. The funding by TODOS technologies is gratefully acknowledged.

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

#### **Appendix A. Impedance Matching Layer: Wave Absorption with a Salisbury Sheet**

It is well established that if a sheet with a thin layer of 377 Ω/ (ohm per square), known as a Salisbury sheet [26], is placed at λ/4 in front of a perfectly conducting sheet, the incident wave in air will be completely absorbed in the sheet and dissipated as heat, without any reflection, as shown in Figure A1.

**Figure A1.** Ideal optical cavity based on Salisbury screen.

The 377 Ω/ (ohm per square) is the intrinsic impedance of air or vacuum and the sheet of resistive film (the Salisbury Sheet) is known as "impedance matching layer".

In this arrangement, the impedance presented to the incident wave at the sheet is the impedance of the sheet backed by an infinite impedance, since at λ/4 the perfectly conducting layer at the bottom is transformed to an infinite impedance. As a result, the incident wave is totally absorbed; however, there is a standing wave and energy circulation between the resistive and conducting sheets and this energy is dissipated as heat in the resistive sheet.

#### **Appendix B. Modeling** *the Ideal Optical Cavity* **with Transmission Lines**

The easiest way to model the ideal optical cavity is by using transmission line theory. Optical waves travel through space as through an infinite transmission line. A pixel with an electrical impedance Z can be modeled as a parallel load between two halves of this line (see Figure A2). The refection coefficient is given by:

$$
\Gamma = \frac{Z\_0 - (Z||Z\_0)}{Z\_0 + (Z||Z\_0)} \tag{A1}
$$

For *Z* = *Z*<sup>0</sup> the reflection is 1/3.

**Figure A2.** Schematic model of transmission line with pixel with impedance Z. (note Z<sup>0</sup> is the impedance of the air or vacuum).

To maximize absorption, we need to minimize reflection. This is achieved with λ/4 optical cavity. We can achieve a perfect matching by using a "λ/4 transformer" terminated by a perfect conductor (*Z*<sup>L</sup> = 0), which transforms the impedance seen by the waves to an infinite impedance, in parallel to the impedance of the pixel (see Figure A3).

**Figure A3.** Schematic model of transmission line with pixel with impedance *Z* and λ/4 transformer for perfect matching.

Accordingly, the reflection coefficient is:

$$\Gamma = \frac{Z\_0 - (Z||Z\_0)}{Z\_0 + (Z||Z\_0)} = \frac{Z\_0 - (Z||\infty)}{Z\_0 + (Z||\infty)} = \frac{Z\_0 - Z}{Z\_0 + Z} \tag{A2}$$

Hence, if *Z* = *Z*0, the reflection is canceled. Since both, the optical cavity and the reflector, do not absorb energy, all the energy is absorbed by the pixel—the only dissipative element.

## **Appendix C. Optical and Material Parameters of the Main Layers Used in the Simulations**

**Table A1.** Refractive index of TiN, SiN, and SiO<sup>2</sup> at λ = 9.5 µm.


Note: It is well known that the physical properties of thin layers are different from the properties of bulk materials and depend on the deposition technology.

#### **References**


**Pao-Hsiung Wang, Yu-Wei Huang and Kuo-Ning Chiang \***

Department of Power Mechanical Engineering, National Tsing Hua University, Hsinchu 300, Taiwan; s101033812@gmail.com (P.-H.W.); s106033852@m106.nthu.edu.tw (Y.-W.H.)

**\*** Correspondence: knchiang@pme.nthu.edu.tw

**Abstract:** The development of fan-out packaging technology for fine-pitch and high-pin-count applications is a hot topic in semiconductor research. To reduce the package footprint and improve system performance, many applications have adopted packaging-on-packaging (PoP) architecture. Given its inherent characteristics, glass is a good material for high-speed transmission applications. Therefore, this study proposes a fan-out wafer-level packaging (FO-WLP) with glass substrate-type PoP. The reliability life of the proposed FO-WLP was evaluated under thermal cycling conditions through finite element simulations and empirical calculations. Considering the simulation processing time and consistency with the experimentally obtained mean time to failure (MTTF) of the packaging, both two- and three-dimensional finite element models were developed with appropriate mechanical theories, and were verified to have similar MTTFs. Next, the FO-WLP structure was optimized by simulating various design parameters. The coefficient of thermal expansion of the glass substrate exerted the strongest effect on the reliability life under thermal cycling loading. In addition, the upper and lower pad thicknesses and the buffer layer thickness significantly affected the reliability life of both the FO-WLP and the FO-WLP-type PoP.

**Keywords:** fan-out wafer-level package; finite element; glass substrate; reliability life; packaging-on-packaging

#### **1. Introduction**

As electronics technology progresses, 3D integrated systems are becoming increasingly important in the realization of lightweight devices with higher performance and better miniaturization. These systems not only effectively reduce the package footprint and weight, but also improve system performance by reducing the system circuit communication length. The fan-out package structure has good electrical and thermal performance and is used in many applications for system integration. Moreover, packaging-on-packaging (PoP) technology allows packages to be stacked three-dimensionally, thus achieving highdensity integration and improving chip-to-chip performance (e.g., in applications with high-frequency data exchange between application processes and memory) [1–3].

Glass is an insulating material, and hence its electrical characteristics are more favorable than those of silicon. In addition, the thickness of a glass substrate can be modified without additional thinning and polishing processes, and through glass via (TGV) technology does not require an additional barrier layer, greatly reducing the production cost. Moreover, the coefficient of thermal expansion (CTE) of glass can be optimized to reduce warpage [4,5] and improve the reliability life of stacked fan-out wafer-level packaging (FO-WLP). Hence, a fan-out structure with a glass substrate is favorable for high-frequency applications [4,6]. However, before mass production, packaging must pass reliability life testing under thermal cycling loading; in the JEDEC(Joint Electron Device Engineering Council) standard (JESD22-A104D), the thermal range is −40 to 125 ◦C. Finite element analysis is widely used to optimize package structures [7–12]. In this study, we proposed an FO-WLP with a glass substrate architecture and used it as a PoP. We fabricated FO-WLP test samples, subjected them to onboard thermal cycling testing (OBTCT), and verified

**Citation:** Wang, P.-H.; Huang, Y.-W.; Chiang, K.-N. Reliability Evaluation of Fan-Out Type 3D Packaging-On-Packaging. *Micromachines* **2021**, *12*, 295. https://doi.org/10.3390/mi12030295

Academic Editor: Seonho Seok

Received: 5 February 2021 Accepted: 5 March 2021 Published: 10 March 2021

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

**Copyright:** © 2021 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/).

the results against those of the finite element models described below. In the fabricated FO-WLP structure, the corner solder joint is the critical failure point because during heating and cooling, stress and strain accumulate in these joints due to mismatches in the CTE of the package and the printed circuit board (PCB), eventually leading to failure.

We further established two-dimensional (2D) and three-dimensional (3D) models and verified their consistency in terms of the mean-time-to-failure (MTTF). To reduce computing time, the 2D model was used to optimize the reliability life by varying the upper and lower pad diameters, the CTE of the glass substrate, and the thickness of the buffer layer. Parametric studies revealed that the reliability life of the optimized FO-WLP and PoP exceeds 1000 cycles, satisfying JEDEC condition G.

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

#### *2.1. Shape Prediction of the Reflowed Solder Joint*

The geometry profile of the solder joints strongly affects the reliability life of a package. Therefore, before analyzing the reliability of the solder joints, the shape of the solder balls must be accurately described. In this study, Surface Evolver [13–15] software based on the energy method was used to describe and predict the solder joint shape.

When a liquid reaches static equilibrium, its total energy tends to be the lowest and its surface area the smallest. The energy of a liquid mainly comprises surface tension energy, gravitational energy, and external energy. From the total energy (Equation (1)), the restoring force in the direction of gravity can be calculated, and the shape and height of the solder ball can be estimated:

$$\delta \mathcal{E}\_{\text{total}} = T\_{\mathcal{S}} \iiint\_{\mathcal{S}} \left( \vec{div} \vec{h} - \vec{n} \cdot \vec{D} \vec{h} \cdot \vec{n} \right) dA + \rho g \iint\_{\mathcal{S}} \left( \vec{div} \vec{\mathcal{Z}} \frac{\rightarrow}{2} \vec{k} \right) \vec{h} - \text{curl} \left( \vec{\hat{h}} \times \frac{z}{2} \vec{k} \right) \cdot \text{d}A - P \iint\_{\mathcal{S}} \vec{h} \cdot \text{d}A \tag{1}$$

If Equation (1) is differentiated once, the restoring force of the solder ball can be expressed as follows:

$$F\_r = \frac{\partial E\_{total}}{\partial H} = \frac{\partial E\_{surface\ tension} + \partial E\_{gravity} + \partial E\_{external\ force}}{\partial H} \tag{2}$$

where *Etotal* is the total energy related to the height of the solder ball *H*, *T<sup>s</sup>* is the surface tension of the solder ball, *P* is the pressure caused by the external force, *A* is the surface area of a single element on the solder ball surface, *z* is the height of the solder ball of a single element surface parallel to the direction of gravity, \* *k* represents the unit vector in the direction of gravity, \* *n* is the unit vector along the positive direction of the element surface, *Fr* is the restoring force, *ρ* is the density of the solder ball, *g* is the acceleration due to gravity, \*

and *V* is the volume of the solder ball. Furthermore, *h* is the perturbation equation:

$$\stackrel{\rightarrow}{h} = \left[ \left( z\_{top} - z \right) / \left( z\_{top} - z\_{base} - H \right) \right] \stackrel{\rightarrow}{k} \tag{3}$$

where *ztop* and *zbase*, respectively, represent the upper and lower boundary conditions of the solder ball (Figure 1). By applying a slight upward or downward interference on the pad, the restoring force of the solder ball in the direction of gravity can be determined. When the restoring force of the solder ball equals the gravity exerted on the solder ball, the molten solder ball achieves static equilibrium. Then, the height and geometry of the solder ball under static equilibrium can be determined.

**Figure 1.** Geometry of the solder joints in the reflow process.

#### *2.2. Life Prediction of the Solder Joints*

In electronic packages subjected to accelerated thermal cycling tests, CTE mismatch between the package and the PCB can cause excessive stress and strain to accumulate in the solder joint with the largest distance from the neutral point (DNP); this may cause the first failure of the packaging. The Coffin–Manson strain-based empirical model is widely used to estimate the fatigue life of solder joints, and the empirical equation [16–18] is as follows:

$$N\_f = \mathbb{C}\left(\mathfrak{e}\_{eq}^{in}\right)^{-\eta} \tag{4}$$

where *N<sup>f</sup>* is the mean time to failure (MTTF), and *C* and *η* are material constants, usually obtained experimentally. In our model and for SAC305 solder material, these are 0.235 and 1.75, respectively [10–12]. The incremental equivalent inelastic strain in each temperature cycle, *ε in eq*, is defined as follows:

$$
\Delta \varepsilon\_{\epsilon \underline{\eta}}^{\dot{m}} = \frac{\sqrt{2}}{3} \left[ \left( \Delta \varepsilon\_x^{\dot{m}} - \Delta \varepsilon\_y^{\dot{m}} \right)^2 + \left( \Delta \varepsilon\_y^{\dot{m}} - \Delta \varepsilon\_z^{\dot{m}} \right)^2 + \left( \Delta \varepsilon\_z^{\dot{m}} - \Delta \varepsilon\_x^{\dot{m}} \right)^2 + \frac{3}{2} \left( \Delta \gamma\_{xy}^{\dot{m}2} + \Delta \gamma\_{yz}^{\dot{m}2} + \Delta \gamma\_{zx}^{\dot{m}2} \right) \right]^{\frac{1}{2}} \tag{5}
$$

where ∆*ε in x* , ∆*ε in y* , ∆*ε in z* , ∆*γ in xy*, ∆*γ in yz*, and ∆*γ in zx* are the incremental inelastic strains in the *x, y*, and *z* directions and the incremental inelastic shear strain in the *xy, yz*, and *zx* directions, respectively.

In this study, we used finite element simulation and the Coffin–Mason empirical model to evaluate the reliability life of the solder joint with the largest DNP.

## **3. Test Vehicle Structure and Thermal Cycling**

#### *3.1. Structure of the Test Vehicle*

The target stackable FO-WLP was a cavity-down chip mounted on a glass interposer using TGV and redistribution lines. These chips can be stacked and molded to form PoPtype packaging using through molding via (TMV). The test vehicle [19] was a simplified FO-WLP (a 10 mm × 10 mm × 0.1 mm chip) assembled on a glass substrate of size 14 mm × 14 mm × 0.1 mm using a die attach film (Nitto, EM700). The chip and substrate were covered by a molding compound of size 14 mm × 14 mm × 0.2 mm to form the package. Next, the package was subjected to mechanical debonding (Figure 2). The underside of the glass substrate contained 432 solder joints arranged in a peripheral layout and connected in series. The diameter and pitch of the solder ball were 250 mm and 400 µm, respectively. The test board dimensions were 77 mm × 132 mm, following the

JEDEC (JESD22-B111) design rule. The dimensions of all of the components are listed in Table 1. In the test board, a daisy-chain structure was used to check the electrical resistance of the chained solder joints (Figure 3). The test vehicle was deemed to have failed if its daisy-chain resistance was infinite.

**Figure 2.** (**a**) Fan-out wafer-level package (FO-WLP) after debonding; (**b**) schematic of the fabricated FO-WLP onboard.


**Table 1.** Dimensions of the components in a fan-out wafer-level package (FO-WLP).

## *3.2. Thermal Cycling and Weibull Distribution*

For reliability testing, we subjected 16 test vehicles to onboard thermal cycling. The test conditions follow the JEDEC standard (condition G): temperature range, −40 to 125 ◦C; ramp rate, 16.5 ◦C/min; dwell time, 10 min. The Weibull distribution of the test results (Figure 4) revealed that the mean time to failure was 249 cycles when the cumulative failure percentage equaled 63.2%. Moreover, all 16 samples failed at the top of the outermost solder joint.

**Figure 4.** Weibull distribution of the glass interposer fan-out package.

#### **4. Finite Element Analysis of FO-WLP**

We established a 2D diagonal semi-symmetric plane strain model and a 3D quarter model to study the thermomechanical behavior of the solder joints in the FO-WLP structure. As the 2D model had fewer nodes than the 3D model, it required less simulation time. However, because plane strain was assumed in the 2D model, it cannot truly represent the actual state of the test vehicle (e.g., the semi-spherical nature of the solder ball). Therefore, a 3D model closer in geometry to the actual vehicle was established to verify the robustness of the 2D model.

After verifying the finite element simulations against the experimental results, we added a TMV component to the FO-WLP structure to make it stackable. The stacked PoP architecture was evaluated by simulation under the JEDEC standard (condition G) testing condition.

#### *4.1. Material Parameters*

The test vehicle was composed of an Si chip, a die attach film, a glass substrate, a stress buffer layer (SBL), copper, SAC305 solder balls, a PCB, and a molding compound. We used three types of glass substrates whose CTE ranged from 3.17 to 9.8 ppm/◦C and whose modulus ranged from 64 to 74 GPa, with little difference in Poisson's ratio. The material parameters are listed in Table 2. The temperature-dependent characteristics of all of the materials were linear, except for the SAC305 solder balls. Figure 5 shows the stress–strain curve of the solder balls at different temperatures [20].

**Table 2.** Summary of the material properties of the components in the FO-WLP. CTE, coefficient of thermal expansion.


**Figure 5.** Stress–strain curves of Sn–3Ag–0.5Cu solder balls at different temperatures.

To model the creep behavior [7] of the solder balls, we adopted the Garofalo–Arrhenius hyperbolic sine model (Equation (6)), which has been widely used to simulate the thermal cycle load of solder joints:

$$\frac{d\varepsilon}{dt} = A \left( \sinh(B\sigma)^n \exp\left(\frac{-Q}{RT}\right) \right) \tag{6}$$

where *ε* is the strain, *σ* is the stress, *A* and *B* are the material constants, *R* is the gas constant, *T* is the absolute temperature, *Q* is the activation energy, and *n* is the stress index. Table 3 lists the values of these constants as used in this study [13].

**Table 3.** Constants in the creep equation [13].


#### *4.2. 2D Plane Strain Model*

For the finite element analysis of the 2D model, we assumed plane strain and used PLANE42 and PLANE182 (ANSYS, 2020R2). Surface Evolver was used to define the shape of the solder joints. Figure 6 depicts the top view of the test vehicle, where the solder joints located in the periphery array format are marked in white. The finite element model was a 2D half model built along the diagonal of the test vehicle (Figure 7). Its position is the yellow line in Figure 6. Figure 8 illustrates the finite element model of the stacked FO-WLP (PoP) structure. The upper and lower stacks of the package were connected through soldering and TMV. The TMV was 0.25 mm in diameter and located atop the solder balls.

**Figure 6.** Top view of the test vehicle.

**Figure 7.** Two-dimensional (2D) finite element model of the FO-WLP.

**Figure 8.** 2D finite element model of the stacked FO-WLP (i.e., packaging-on-packaging (PoP)).

## *4.3. 3D Quarter Symmetry Model*

Unlike a 3D finite element model, a 2D plane strain model cannot accurately represent semi-spherical-type features. However, a 3D model with a fine mesh contains a large number of elements and necessitates excessive simulation computing time, making parametric study infeasible. Two approaches can be used to overcome these drawbacks. The first is to build a 3D model to verify the 2D plane strain model, which can then be used to evaluate the packaging structure, and the second is to reduce the number of elements in the 3D model. In this work, we used multipoint constraint (MPC) technology to significantly reduce the number of elements. Panels (a) and (c) in Figure 9 illustrate 1/4 models of a simple WLP with full fine meshes, and panels (b) and (d) illustrate their MPC-reduced counterparts. This simple 1/4 WLP model will not take too much simulation time and was selected in this research to verify that the MPC model can obtain stress/strain results similar to the full fine mesh model.

**Figure 9.** Three-dimensional (3D) finite element model: (**a**) Full fine-mesh model; (**b**) multipoint constraint (MPC) model; (**c**) cross-section view of the full fine-mesh model; (**d**) cross-section view of the MPC model.

After verifying the MPC model against the full model, the MPC model was used to simulate the test vehicle and the stacked PoP architecture. Figure 10 shows the 3D quarter symmetry finite element model.

**Figure 10.** Three-dimensional MPC-reduced finite element models for the (**a**) FO-WLP and (**b**) FO-WLP PoP.

#### **5. Results and Discussion**

#### *5.1. 2D Diagonal Plane Strain Model*

After OBTCT, we inspected the failed test vehicle and found that the crack was located atop the solder joint with the largest DNP. Consistent with the experiment wherein we used glass substrate A, the simulation also showed that the maximum equivalent inelastic strain/stress was atop the outermost solder joint (Figure 11). In addition, the 248 life cycles simulated by Equation (4) are in good agreement with the experiment results of 249 mean cycles to failure (Figure 4). Thus, the failure prediction of the 2D model is in good agreement with the experiment data.

**Figure 11.** Failure point of solder joint after onboard thermal cycling testing (OBTCT): (**a**) Cross-section view; (**b**) equivalent inelastic strain distribution.

Next, we conducted a parametric study to optimize the reliability life under different combinations of upper and lower pad diameters. As evident in Figure 12, an upper pad larger than a lower pad resulted in better reliability. The Young's modulus of silicon material is stronger than that of PCB, so a larger upper pad size is required to reduce the equivalent inelastic strain. A lower strain can have better reliability life, and the strain of the solder ball is related to the pad size, contact angle, and standoff height of the solder ball. In the simulation, when the upper pad diameter was 250 µm and the lower pad diameter was 180 µm, the best reliability life cycles could be achieved. That is, the pad on the interposer side should be larger than that on the PCB side. The optimal reliability life of 368 cycles was achieved at an upper pad/lower pad diameters ratio of approximately 1:0.72.

**Figure 12.** Simulated reliability life under different combinations of pad thicknesses.

To meet condition G of the JEDEC standard, the thermal cycling life of the product must exceed 1000 cycles. To meet this standard with the optimized pad combination (upper pad 250 µm and lower pad 180 µm), we further optimized the SBL thickness and CTE of the glass substrate (Figure 13). Specifically, we simulated SBL thicknesses of 5, 10, and 20 µm and both substrate glasses A and B. Glasses A and B differ little in terms of Poisson's ratio, but differ substantially in CTE (3.17 and 8.37 ppm/◦C, respectively). Under the same conditions, the reliability life was higher for glass B than for glass A by 600 cycles on average, meaning that the reliability life is most sensitive to glass CTE. Moreover, a thicker SBL reduced solder stress/strain and thus favorably affected the thermal cycling reliability. Overall, the reliability life improved from 248 cycles to 1432 cycles with the following optimized combination: upper pad diameter, 250 µm; lower pad diameter, 180 µm; SBL thickness, 20 µm; glass CTE, 8.35 ppm/◦C.

**Figure 13.** Simulated reliability life for stress buffer layers of different thicknesses.

#### *5.2. 3D Quarter Symmetry Model*

The 3D FO-WLP PoP simulation model is much larger than the nonstacked FO-WLP, which makes 2D simulation the only feasible option for the parametric study of a stacked FO-WLP. To verify whether a 2D model can be applied to stacked packaging, a 3D model must first be established. To reduce the computing time, the number of elements and nodes must be decreased. First, we established a small WLP model (four solder balls per quarter model; Figure 9) both as a full fine-mesh model (86,656 elements; Figure 9a) and an MPC model (32,648 elements; Figure 9b). The corresponding computing times were 5 h 43 min and 2 h 13 min. We applied the same thermal loading as in the earlier analyses (−40 to 125 ◦C). Strain was found to be concentrated in a corner of the ball and the upper edge (Figure 14). After one to three thermal cycles, the strain increment in the two models differed by less than 1% (Figure 15). These simulation results demonstrate the feasibility of applying the MPC approach to simulate FO-WLP PoP with 112 solder balls per quarter model, as shown in Figure 10.

Next, we applied the design optimized in Section 5.1 (upper pad diameter, 250 µm; lower pad diameter, 180 µm; SBL thickness, 20 µm; glass CTE, 8.35 ppm/◦C) to the 3D quarter symmetry model. Figure 16 shows the equivalent inelastic strain distribution in the FO-WLP. Strain accumulated in the solder ball with the maximum DNP, and the cycling life predicted by the 3D model differed from the 2D model prediction by only 1.3%.

**Figure 14.** Equivalent inelastic strain distribution: (**a**) Full fine-mesh model; (**b**) MPC model; (**c**) maximum strain in the full model (top view); (**d**) maximum strain in the MPC model (top view).

**Figure 15.** Strain increment in the full and MPC models.

**Figure 16.** Equivalent inelastic strain distribution of the FO-WLP: (**a**) All solder joints; (**b**) location of the maximum strain; (**c**) top view of the bottom solder joint with the maximum distance from the neutral point (DNP); (**d**) bottom view of the solder joint with the maximum DNP.

We repeated the above simulation for the FO-WLP PoP structure. Again, we found that strain concentrated in the bottom layer of the solder ball with the maximum DNP, replicating the trend of failure at a corner location. Moreover, the maximum strain was in the upper pad (Figure 17). The upper and lower parts of the solder balls in the bottom layer were connected to the glass substrate and the PCB, respectively. Therefore, these solder balls could withstand greater stress than could solder balls in the upper layer. However, this double-layer structure reduced the cycling life to 974 cycles. In the PoP structure, the percentage of glass in this package increased; therefore, the equivalent inelastic strain due to CTE mismatch also increased, and a higher strain may lead to earlier failure.

*Micromachines* **2021**, *12*, x FOR PEER REVIEW 12 of 14

**Figure 17.** Equivalent inelastic strain distribution in FO-WLP PoP: (**a**) All solder joints; (**b**) location of the maximum strain; (**c**) top view of the solder joint with the maximum DNP; (**d**) bottom view of the solder joint with the maximum DNP. Herein, the 2D and 3D models exhibited good consistency, with only 5% difference in reliability life. Therefore, we used the verified 2D models in the parametric analysis to

*Model*

Herein, the 2D and 3D models exhibited good consistency, with only 5% difference in reliability life. Therefore, we used the verified 2D models in the parametric analysis to optimize the design parameters in order to extend the reliability life of the FO-WLP PoP structure to more than 1000 cycles. optimize the design parameters in order to extend the reliability life of the FO‐WLP PoP structure to more than 1000 cycles. *5.3. Parametric Analysis of FO‐WLP PoP Structure Using the 2D Plane Strain Finite Element*

the percentage of glass in this package increased; therefore, the equivalent inelastic strain due to CTE mismatch also increased, and a higher strain may lead to earlier failure.

#### *5.3. Parametric Analysis of FO-WLP PoP Structure Using the 2D Plane Strain Finite Element Model* In addition to the already optimized pad combination, three design parameters were considered in our parametric study: SBL thickness, chip thickness, and CTE of the glass

In addition to the already optimized pad combination, three design parameters were considered in our parametric study: SBL thickness, chip thickness, and CTE of the glass substrate. In this parametric analysis, SBL thicknesses of 20, 25, and 30 µm, chip thicknesses of 90, 100, and 150 µm, and glass substrate CTEs of 8.35 and 9.8 ppm/◦C were evaluated, yielding 16 combinations (Table 4). SBLs thicker than 30 µm were not evaluated due to the manufacturing difficulty. The results showed that the thicker the SBL, the higher the reliability and the higher the released stress/strain concentration. Furthermore, a lower chip thickness yielded a higher molding compound volume, which in turn increased the CTE of the whole package and thus narrowed the CTE mismatch between the package and PCB. Most importantly, glass CTE was found to have the strongest decreasing effect on the CTE mismatch between the glass substrate and the PCB. In summary, the reliability life of FO-WLP PoP can be increased by optimizing the SBL thickness and glass CTE (Figure 18). Of the 16 simulated combinations, combination 13—SBL thickness, 30 µm; chip thickness, 90 µm; glass substrate CTE, 9.8 ppm/◦C—yielded the longest life cycle of 1427 cycles (Table 4). substrate. In this parametric analysis, SBL thicknesses of 20, 25, and 30 μm, chip thick‐ nesses of 90, 100, and 150 μm, and glass substrate CTEs of 8.35 and 9.8 ppm/°C were eval‐ uated, yielding 16 combinations (Table 4). SBLs thicker than 30 μm were not evaluated due to the manufacturing difficulty. The results showed that the thicker the SBL, the higher the reliability and the higher the released stress/strain concentration. Furthermore, a lower chip thickness yielded a higher molding compound volume, which in turn in‐ creased the CTE of the whole package and thus narrowed the CTE mismatch between the package and PCB. Most importantly, glass CTE was found to have the strongest decreas‐ ing effect on the CTE mismatch between the glass substrate and the PCB. In summary, the reliability life of FO‐WLP PoP can be increased by optimizing the SBL thickness and glass CTE (Figure 18). Of the 16 simulated combinations, combination 13—SBL thickness, 30 μm; chip thickness, 90 μm; glass substrate CTE, 9.8 ppm/°C—yielded the longest life cycle of 1427 cycles (Table 4). 

**Figure 18.** Life predictions of the FO‐WLP PoP under different (**a**) SBL thicknesses, (**b**) chip thicknesses, and (**c**) CTEs. **Figure 18.** Life predictions of the FO-WLP PoP under different (**a**) SBL thicknesses, (**b**) chip thicknesses, and (**c**) CTEs.

**Table 4.** Combination of the design parameters in the parametric analysis of the FO‐WLP PoP.

**Commented [M2]:** We recommend five or more digits, separated by commas between every three digits both in the text and in the picture, for four digits, commas are not needed, please revise.


**Table 4.** Combination of the design parameters in the parametric analysis of the FO-WLP PoP.

#### **6. Conclusions**

The simulation of 3D models is infeasible because of the high processing power and time requirements. In this study, to optimize the reliability life of FO-WLP and FO-WLP-type PoP under thermal cycling loading, 2D and 3D MPC technique finite element simulation models were established and experimentally verified to have similar accuracy. The reliability life of the verified 2D models was then optimized through parametric analysis. The CTE of the glass substrate was found to exert the strongest effect on reliability life. In addition, the upper and lower pad diameters and buffer layer thickness significantly affected the reliability life of the FO-WLP and FO-WLP-type PoP. Specifically, the reliability life of the solder joints was highest when the upper pad was larger than the lower pad (1:0.72). Moreover, a thicker buffer layer released a greater stress/strain concentration and thus positively affected the reliability life. However, SBLs thicker than 30 µm were not evaluated in this study due to the manufacturing difficulty.

**Author Contributions:** P.-H.W. and Y.-W.H. studied the reliability life prediction and structural optimization methods of the PoP architecture. K.-N.C. is a professor who led and discussed the content of this research. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors would like to thank the Ministry of Science and Technology (Project MOST 108-2221-E-007-077-MY3) of Taiwan for supporting this research.

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

#### **References**


**Qian-Yo Lee <sup>1</sup> , Ming-Xuan Lee <sup>2</sup> and Yen-Chun Lee 3,\***


**Abstract:** Integrated devices incorporating MEMS (microelectromechanical systems) with IC (integrated circuit) components have been becoming increasingly important in the era of IoT (Internet of Things). In this study, a hybrid fuzzy MCDM (multi-criteria decision making) model was proposed to effectively evaluate alternative technologies that incorporate MEMS with IC components. This model, composed of the fuzzy AHP (analytic hierarchy process) and fuzzy VIKOR (VIseKriterijumska Optimizacija I Kompromisno Resenje) methods, solves the decision problem of how best to rank MEMS and IC integration technologies in a fuzzy environment. The six important criteria and the major five alternative technologies associated with our research themes were explored through literature review and expert investigations. The priority weights of criteria were derived using fuzzy AHP. After that, fuzzy VIKOR was deployed to rank alternatives. The empirical results show that development schedule and manufacturing capability are the two most important criteria and 3D (three-dimensional) SiP (system-in-package) and monolithic SoC (system-on-chip) are the top two favored technologies. The proposed fuzzy decision model could serve as a reference for the future strategic evaluation and selection of MEMS and IC integration technologies.

**Keywords:** technology evaluation; MEMS and IC integration; MCDM; fuzzy AHP; fuzzy VIKOR
