Theoretical and Experimental Investigation of Warpage Evolution of Flip Chip Package on Packaging during Fabrication

This study attempts to investigate the warpage behavior of a flip chip package-on-package (FCPoP) assembly during fabrication process. A process simulation framework that integrates thermal and mechanical finite element analysis (FEA), effective modeling and ANSYS element death-birth technique is introduced for effectively predicting the process-induced warpage. The mechanical FEA takes into account the viscoelastic behavior and cure shrinkage of the epoxy molding compound. In order to enhance the computational and modeling efficiency and retain the prediction accuracy at the same time, this study proposes a novel effective approach that combines the trace mapping method, rule of mixture and FEA to estimate the effective orthotropic elastic properties of the coreless substrate and core interposer. The study begins with experimental measurement of the temperature-dependent elastic and viscoelastic properties of the components in the assembly, followed by the prediction of the effective elastic properties of the orthotropic interposer and substrate. The predicted effective results are compared against the results of the ROM/analytical estimate and the FEA-based effective approach. Moreover, the warpages obtained from the proposed process simulation framework are validated by the in-line measurement data, and good agreement is presented. Finally, key factors that may influence process-induced warpage are examined via parametric analysis.


Introduction
In recent years, there has been explosive and continuous growth in the consumer market for various smart products and Internet of Things (IoT) products, as well as the developing requirements of 5G communication, artificial intelligence (AI), and autonomous vehicles. Advanced packaging technology like flip chip packaging [1], wafer level packaging [2], and flip chip chip-scale packaging (FCCSP) [3,4] was introduced to achieve high I/O density, excellent electrical performance and miniaturization, and thus is commonly used in high-end smart chips in recent years. However, the physical limitations of Moore's law [5] make it difficult for electronic packaging to continuously shrink and functionally improve, prompting the development of new packaging technologies. Among the many solutions, the system-in-package (SiP) for heterogeneous integration is the current alternative, and is one of the most feasible methods for "More than Moore" or even "Beyond CMOS".
SiP technology may have a two-dimensional (2D) planar configuration, a threedimensional (3D) vertical stacking configuration, or an integrated (hybrid) configuration. 3D packaging technology can be categorized into package stacking, like package-onpackage (PoP) and package-in-package (PiP), wire-bonding [6], and through silicon via (TSV)-based 3D IC stacking [7]. In addition to high I/O quantity and miniaturization, further requirements of multi-functionality have aroused the development of the flip chip HBM for heterogeneous integration, the CCSBs were utilized to connect the bottom FCCSP and the core interposer, and the EMC was used to protect the solder balls. To minimize the package profile, a three-layer 100 µm embedded trace substrate (ETS) was applied, as schematically depicted in Figure 2, which mainly included two solder mask (SM) protective layers, two prepreg (PP) dielectric layers, and three metal (Cu) layers, with the circuitries filled with either SM or PP material. The chip was 9.36 mm in length, 8.76 mm in width, and 70 µm thick. The chip was connected to the coreless substrate using 2500 Cu pillar bumps, which were 40 µm in length, 70 µm in width and 58 µm in height. The gap between the chip and coreless substrate was filled with an underfill via the capillary action. The top layer was stacked on a 90 µm thick two-layer core interposer. In the FCPoP assembly, there were a total of 550 CCSBs with a diameter of 190 µm. Finally, the EMC was filled between the substrate and interposer to form an FCPoP assembly with a length and width of 14 mm and a thickness of 500 µm. Figure 3 describes the main fabrication process steps; i.e., the die bonding process (steps 1-3), underfill cure process (steps 3-6), interposer bonding process (steps 6-9), mold cure process (steps 9-12) and temperature elevation process (steps [12][13], and also the corresponding process temperatures.

Structure and Fabrication Process of FCPoP
The research vehicle was an FCPoP assembly, as shown in Figure 1, that was primarily composed of an FCCSP package, an EMC, a core interposer, and Cu core solder balls (CCSBs). The main structure of the FCCSP package included a silicon chip, an underfill, Cu pillar bumps, and a coreless substrate. The FCCSP package usually is used for highend processors. In addition, the core interposer was applied to facilitate the connection with the HBM for heterogeneous integration, the CCSBs were utilized to connect the bottom FCCSP and the core interposer, and the EMC was used to protect the solder balls. To minimize the package profile, a three-layer 100 μm embedded trace substrate (ETS) was applied, as schematically depicted in Figure 2, which mainly included two solder mask (SM) protective layers, two prepreg (PP) dielectric layers, and three metal (Cu) layers, with the circuitries filled with either SM or PP material. The chip was 9.36 mm in length, 8.76 mm in width, and 70 μm thick. The chip was connected to the coreless substrate using 2500 Cu pillar bumps, which were 40 μm in length, 70 μm in width and 58 μm in height. The gap between the chip and coreless substrate was filled with an underfill via the capillary action. The top layer was stacked on a 90 μm thick two-layer core interposer. In the FCPoP assembly, there were a total of 550 CCSBs with a diameter of 190 μm. Finally, the EMC was filled between the substrate and interposer to form an FCPoP assembly with a length and width of 14 mm and a thickness of 500 μm. Figure 3 describes the main fabrication process steps; i.e., the die bonding process (steps 1-3), underfill cure process (steps 3-6), interposer bonding process (steps 6-9), mold cure process (steps 9-12) and temperature elevation process (steps [12][13], and also the corresponding process temperatures.

Structure and Fabrication Process of FCPoP
The research vehicle was an FCPoP assembly, as shown in Figure 1, that was primarily composed of an FCCSP package, an EMC, a core interposer, and Cu core solder balls (CCSBs). The main structure of the FCCSP package included a silicon chip, an underfill, Cu pillar bumps, and a coreless substrate. The FCCSP package usually is used for highend processors. In addition, the core interposer was applied to facilitate the connection with the HBM for heterogeneous integration, the CCSBs were utilized to connect the bottom FCCSP and the core interposer, and the EMC was used to protect the solder balls. To minimize the package profile, a three-layer 100 μm embedded trace substrate (ETS) was applied, as schematically depicted in Figure 2, which mainly included two solder mask (SM) protective layers, two prepreg (PP) dielectric layers, and three metal (Cu) layers, with the circuitries filled with either SM or PP material. The chip was 9.36 mm in length, 8.76 mm in width, and 70 μm thick. The chip was connected to the coreless substrate using 2500 Cu pillar bumps, which were 40 μm in length, 70 μm in width and 58 μm in height. The gap between the chip and coreless substrate was filled with an underfill via the capillary action. The top layer was stacked on a 90 μm thick two-layer core interposer. In the FCPoP assembly, there were a total of 550 CCSBs with a diameter of 190 μm. Finally, the EMC was filled between the substrate and interposer to form an FCPoP assembly with a length and width of 14 mm and a thickness of 500 μm. Figure 3 describes the main fabrication process steps; i.e., the die bonding process (steps 1-3), underfill cure process (steps 3-6), interposer bonding process (steps 6-9), mold cure process (steps 9-12) and temperature elevation process (steps [12][13], and also the corresponding process temperatures.

Linear Viscoelasticity
The properties of viscoelastic materials can be divided into elastic and viscoelastic

Linear Viscoelasticity
The properties of viscoelastic materials can be divided into elastic and viscoelastic parts. The elastic part can react immediately and obey Hooke's law when subjected to a fixed load. Conversely, the viscoelastic part gradually increases the strain (creep) or reduces the stress (relaxation), and eventually reaches stability. Polymer materials typically show viscoelastic relaxation behavior [13][14][15][16]. The stress relaxation effect is usually dominated by chemical phenomena at high temperatures for a long period of time. To describe the viscoelastic relaxation behavior, the generalized Maxwell model is most commonly used. It is composed of several Maxwell elements and an independent spring combined in parallel. The time-dependent viscoelastic stress relaxation modulus can be expressed by a Prony series mathematical representation [11]: where n is the number of Maxwell elements, E k denotes the modulus of each Maxwell element, E ∞ represents the relaxed modulus, ς k stands for the relaxation time, and t is the time. Furthermore, the modulus of each Maxwell element can be described as: where c k denotes the weighting factor and E 0 represents the unrelaxed modulus, expressed as Combining Equations (1)-(3) yields: The properties of polymer materials, in either a glass or a rubbery state, show a strong relationship with temperature. The glassy state refers to the polymer material at a temperature higher than the glass transition temperature (T g ). In contrast, the rubbery state refers to the polymer material at a temperature lower than T g . The Young's modulus and CTE have a large variation during the phase transition state. The time-temperature superposition (TTS) principle is often applied to depict the time-temperature dependence of the linear viscoelastic behavior. The TTS principle illustrates that the relaxation curve of the material's modulus and time (or frequency) at a certain temperature is analogous to the relaxation curve of the adjacent temperature. The relaxation curve at each temperature, except the reference temperature, is translated in the logarithmic time domain to form the relaxation curve of the material at the reference temperature. The value of this translation highly depends on the temperature, the reference temperature and the properties of the polymer materials. The TTS principle can be simply expressed as follows: where τ stands for the reduced time when T 0 < T, which can be derived below: In Equation (6), κ T is the shift factor of temperature. This parameter can be approximated by the well-known Williams-Landel-Ferry (WLF) model [17] as: where T ref is the reference temperature, and a 1 and a 2 represent the coefficients of the curve fit. Basically, they are highly dependent on the materials and the reference temperature.

Linear Elastic Mechanics
Polymer materials, such as EMCs, will exhibit volumetric changes or chemical shrinkage during the mold cure process due to chemical reactions. In essence, the extent of volumetric change highly depends on the material's cure state in an isothermal isobaric ensemble. Consider that a polymer cube has a length of one-unit side before curing. The volumetric change ( V) of the unit cube after full curing is written as: where ∆l is the variation of the unit side length. The corresponding strain due to the volumetric change ε vc can be expressed as: The total strain of the polymer materials is the sum of the volumetric change-inducted strain (ε vc ), elastic strain (ε e ), and thermal strain (ε thermal ). According to Hooke's law, the relationship between the elastic strain and stress can be written as follows: where {w} is the nodal displacement vector and [B] represents the strain-displacement matrix. To effectively capture the process-induced warpage behavior of the FCPoP assembly relies on a dependable and accurate thermo-mechanical characterization of the core interposer and coreless substrate, consisting of several Cu circuit layers with multi-material and multi-scale structures and complex geometric features. Accurately and fully modeling them presents great challenges due to the requiring extensive, tedious effort and cost required. To ease the modeling challenges, the core interposer and coreless substrate were approximated as an equivalent homogeneous medium and their effective elastic properties were evaluated by using an effective approach that made use of the powerful electronic computer-aided design (ECAD) trace mapping (TM) method together with the rule-of-mixture (ROM) technique and finite element analysis (FEA).

Numerical Modeling and Material Characterization
The ECAD TM method enables a more efficient and accurate representation of tiny, delicate, complex Cu traces, pads and vias surrounded by the PP dialectic material and the SM protective material on the coreless substrate and core interposer [18]. A flowchart of the ECAD TM method is shown in Figure 4. It is noted that the Cu circuit layers consisted of not only Cu traces, pads and vias, but also SM or PP dielectric material. First of all, based on the ECAD model, a uniform regular background mesh with a significant number of very fine first-order brick elements was established on each Cu circuit layer of the coreless substrate and core interposer. Then, the spatially non-uniformly distributed Cu circuitries in the ECAD model were mapped onto the background mesh to obtain a high resolution (HR) Cartesian Cu circuit map and a finite element (FE) model. Based on the volume ratio of Cu and neighboring materials, such as PP or SM on each brick element of the background mesh, the effective isotropic elastic properties were calculated using an ROM technique, by which a complete material property map of the Cu circuit layer was derived ( Figure 5). The benefits of the ECAD TM method were a very uniform regular mesh, flexible mesh density control and close match to the geometry of the Cu circuitries.
sisted of not only Cu traces, pads and vias, but also SM or PP dielectric material. First of all, based on the ECAD model, a uniform regular background mesh with a significant number of very fine first-order brick elements was established on each Cu circuit layer of the coreless substrate and core interposer. Then, the spatially non-uniformly distributed Cu circuitries in the ECAD model were mapped onto the background mesh to obtain a high resolution (HR) Cartesian Cu circuit map and a finite element (FE) model. Based on the volume ratio of Cu and neighboring materials, such as PP or SM on each brick element of the background mesh, the effective isotropic elastic properties were calculated using an ROM technique, by which a complete material property map of the Cu circuit layer was derived ( Figure 5). The benefits of the ECAD TM method were a very uniform regular mesh, flexible mesh density control and close match to the geometry of the Cu circuitries. As soon as the material property map and 3D FE model of the Cu circuit layer were created, FEA was applied to calculate its effective orthotropic elastic properties. For modeling simplification, the orthotropic Cu circuit layer could be approximated as a transversely isotropic material by averaging the effective in-plane elastic properties and further as an isotropic material by averaging the effective in-plane and out-of-plane elastic properties. The 3D FE models of the mapped Cu circuit layer, PP dielectric and SM layers could be combined together to form an integrated 3D FE model of the substrate and interposer. Finally, the effective orthotropic elastic properties of the substrate and interposer as a whole could be derived using FEA. It is worth mentioning that the final modeling step; sisted of not only Cu traces, pads and vias, but also SM or PP dielectric material. First of all, based on the ECAD model, a uniform regular background mesh with a significant number of very fine first-order brick elements was established on each Cu circuit layer of the coreless substrate and core interposer. Then, the spatially non-uniformly distributed Cu circuitries in the ECAD model were mapped onto the background mesh to obtain a high resolution (HR) Cartesian Cu circuit map and a finite element (FE) model. Based on the volume ratio of Cu and neighboring materials, such as PP or SM on each brick element of the background mesh, the effective isotropic elastic properties were calculated using an ROM technique, by which a complete material property map of the Cu circuit layer was derived ( Figure 5). The benefits of the ECAD TM method were a very uniform regular mesh, flexible mesh density control and close match to the geometry of the Cu circuitries. As soon as the material property map and 3D FE model of the Cu circuit layer were created, FEA was applied to calculate its effective orthotropic elastic properties. For modeling simplification, the orthotropic Cu circuit layer could be approximated as a transversely isotropic material by averaging the effective in-plane elastic properties and further as an isotropic material by averaging the effective in-plane and out-of-plane elastic properties. The 3D FE models of the mapped Cu circuit layer, PP dielectric and SM layers could be combined together to form an integrated 3D FE model of the substrate and interposer. Finally, the effective orthotropic elastic properties of the substrate and interposer as a whole could be derived using FEA. It is worth mentioning that the final modeling step; As soon as the material property map and 3D FE model of the Cu circuit layer were created, FEA was applied to calculate its effective orthotropic elastic properties. For modeling simplification, the orthotropic Cu circuit layer could be approximated as a transversely isotropic material by averaging the effective in-plane elastic properties and further as an isotropic material by averaging the effective in-plane and out-of-plane elastic properties. The 3D FE models of the mapped Cu circuit layer, PP dielectric and SM layers could be combined together to form an integrated 3D FE model of the substrate and interposer. Finally, the effective orthotropic elastic properties of the substrate and interposer as a whole could be derived using FEA. It is worth mentioning that the final modeling step; i.e., approximation of the substrate and interposer as a homogeneous equivalent continuum, may not be indispensable, since FE modeling of the substrate and interposer could be directly carried out using the integrated 3D FE model. This effective approach is hereinafter termed the TM/FEA effective method.

The ROM/Analytical Estimate
The effective in-plane and out-of-plane CTEs of the Cu circuit layers could be also assessed using an analytical estimate integrated with an ROM method (It is alternatively termed the ROM/analytical estimate). Specifically, the effective in-plane CTE α x,y of the Cu circuit layers was evaluated according to the literature [19] using an energy approach, and the effective out-of-plane CTE α z also was derived based on the work of [19], and also as presented in [7]: where E 1 (E 2 ), α 1 (α 2 ), and ξ 1 (ξ 2 ) are the Young's modulus, CTE and volume fraction of the Cu (SM or PP), respectively. In Equation (12), the effective Poisson's ratio υ can be simply approximated using ROM as: where ξ 1 and ξ 2 are the Poisson's ratio of the Cu and SM or PP respectively. Likewise, the effective in-plane elastic modulus E x,y and out-of-plane elastic modulus E z of the Cu circuit layers can be also estimated as: FEA using a detailed fine mesh model could be directly applied to derive the effective orthotropic elastic properties of the Cu circuit layers. This method can be very effective in accurately grasping the crucial parameters affecting the effective properties but require a very tedious, time-consuming and complex procedure to model and simulate the material models [20]. The method is briefly termed the FEA-based effective approach [7]. The underlying idea behind this approach is that the elastic responses of the homogeneous equivalent continuum should be consistent with those of the original continuum.
The effective CTEs of the Cu circuit layers could be simply calculated based on the strength of the materials, where δ i is the thermal deformation, α i (i = x, y, z) stands for the effective CTE in the i-th direction, ∆T denotes the temperature increment, and L i represents the side length of the Cu circuit layers in the i-th direction.
In accordance with the generalized Hooke's law, the stress-strain relationship of an orthotropic material is expressed as: where ε(ε x , ε y , ε z ) and σ(σ x , σ y , σ z ) are the normal strain and stress, respectively, γ(γ xy , γ yz , γ xz ) and τ(τ xy , τ yz , τ xz ) represent the shear strain and stress, respectively, and υ(υ xy , υ yz , υ xz ) denotes the Poisson's ratio. In total, there are nine independent effective elastic constants to be determined for an orthotropic elastic material, which are E x , E y , E z , υ xy , υ yz , υ xz , G xy , G yz , and G xz . These constants can be simply derived based on Equations (17)-(22) through FEAs with a set of different loading and boundary conditions. The rest of the effective elastic constants υ yx , υ zy , υ zx can be readily derived from the fact that the compliance matrix is symmetric.

Process Modeling
In this study, a process simulation framework that incorporated the proposed TM/FEA effective method, thermal and mechanical FEAs and the element death and birth method in ANSYS was introduced. Due to symmetry, a quarter of the FCPoP assembly was simulated via the proposed process simulation framework. Figure 6 reveals the constructed 3D FE model of the FCPoP assembly, which comprised 193,401 nodes and 185,211 solid elements. It consisted of the main components of the FCPoP assembly, including a coreless interposer, a core substrate, an EMC, a silicon chip, solder balls, Cu pillar bumps, and an underfill. The displacement boundary conditions are set to simulate the symmetry boundary condition, where the out-of-plane displacement of the nodes on the symmetry planes were constrained. In addition, the bottom node on the intersecting line of these two symmetry planes was fixed in the z-direction to prevent rigid body motion. All the materials in the assembly were assumed to be either linearly elastic and isotropic or orthotropic except the EMC, which was assumed to be linearly viscoelastic. It was noteworthy that the temperature dependence of these materials and the effects of curing shrinkage of the EMC were also taken into account in this investigation. The temperature-dependent elastic properties of the components in the assembly are characterized using a thermal-mechanical analyzer (TMA) (TA Instruments, New Castle, DE, USA) and a dynamic mechanical analyzer (DMA) (TA Instruments, New Castle, DE, USA), as shown in Figure 7. The stress-free temperature of the EMC was defined as its cure temperature. The curing process of the EMC involved two processes: in-mold cure (IMC) and post-mold cure (PMC). Typically, the IMC process applies a lower temperature and a shorter curing time All the materials in the assembly were assumed to be either linearly elastic and isotropic or orthotropic except the EMC, which was assumed to be linearly viscoelastic. It was noteworthy that the temperature dependence of these materials and the effects of curing shrinkage of the EMC were also taken into account in this investigation. The temperature-dependent elastic properties of the components in the assembly are characterized using a thermal-mechanical analyzer (TMA) (TA Instruments, New Castle, DE, USA) and a dynamic mechanical analyzer (DMA) (TA Instruments, New Castle, DE, USA), as shown in Figure 7. All the materials in the assembly were assumed to be either linearly elastic and isotropic or orthotropic except the EMC, which was assumed to be linearly viscoelastic. It was noteworthy that the temperature dependence of these materials and the effects of curing shrinkage of the EMC were also taken into account in this investigation. The temperature-dependent elastic properties of the components in the assembly are characterized using a thermal-mechanical analyzer (TMA) (TA Instruments, New Castle, DE, USA) and a dynamic mechanical analyzer (DMA) (TA Instruments, New Castle, DE, USA), as shown in Figure 7. The stress-free temperature of the EMC was defined as its cure temperature. The curing process of the EMC involved two processes: in-mold cure (IMC) and post-mold cure (PMC). Typically, the IMC process applies a lower temperature and a shorter curing time to increase the stiffness of the EMC. Subsequently, a PMC process with a higher temper- The stress-free temperature of the EMC was defined as its cure temperature. The curing process of the EMC involved two processes: in-mold cure (IMC) and post-mold cure (PMC). Typically, the IMC process applies a lower temperature and a shorter curing time to increase the stiffness of the EMC. Subsequently, a PMC process with a higher temperature and a longer duration was utilized to completely cure the EMC. However, during the curing, the EMC would experience a volumetric (chemical) expansion or contraction. The measured volumetric change data provided by the manufacturer was applied. It was reported that the EMC during the curing process from the gel point (where the stiffness of the EMC is nearly developed) to a full cure state would cause a 0.09% volumetric shrinkage. The process modeling for the warpage prediction of the FCPoP assembly during fabrication closely adhered to the process steps shown in Figure 3.

Characterization of EMC Viscoelastic Properties
In this work, a DMA measurement system was applied to conduct the stress relaxation experiments in the frequency domain through a three-point bending mode. The storage moduli of the EMC under a 0.5% applied strain over a wide frequency scan ranging from 0.1 Hz to 100 Hz and a broad isothermal temperature range of 25-260 • C with 5 • C increment were derived, and some of the results are shown in Figure 8a. These stress relaxation storage moduli were further approximated by the Ninomiya-Ferry method [21] as: where ω = 1/t, E represents the stress relaxation modulus, E is the storage modulus, and E denotes the loss modulus. It is clear that the stress relaxation storage modulus would have a strong temperature correlation at temperatures neighboring the T g of the EMC, which was around 100 • C. In addition, the stress relaxation storage modulus near the T g showed a great time dependence, and at temperatures lower than 50 • C and higher than 200 • C exhibited trivial time and temperature correlations. shrinkage. The process modeling for the warpage prediction of the FCPoP assembly during fabrication closely adhered to the process steps shown in Figure 3.

Characterization of EMC Viscoelastic Properties
In this work, a DMA measurement system was applied to conduct the stress relaxation experiments in the frequency domain through a three-point bending mode. The storage moduli of the EMC under a 0.5% applied strain over a wide frequency scan ranging from 0.1 Hz to 100 Hz and a broad isothermal temperature range of 25 °C -260 °C with 5 °C increment were derived, and some of the results are shown in Figure 8a. These stress relaxation storage moduli were further approximated by the Ninomiya-Ferry method [21] as: where 1/ ω t , E represents the stress relaxation modulus, E is the storage modulus, and E denotes the loss modulus. It is clear that the stress relaxation storage modulus would have a strong temperature correlation at temperatures neighboring the Tg of the EMC, which was around 100 °C. In addition, the stress relaxation storage modulus near the Tg showed a great time dependence, and at temperatures lower than 50 °C and higher than 200 °C exhibited trivial time and temperature correlations. Based on the TTS principle, a single master curve could be constructed by shifting these frequency-dependent storage moduli at different temperatures along the time axis, as shown in Figure 8b, where the reference temperature was set to the Tg of the EMC. The master curve could be well fitted by a Prony series equation using 22 Prony elements; the fitted weighting coefficients k c and the relaxation times k ς are listed in Table 1. The corresponding temperature shift factors in logarithmic scale are shown in Figure 9, as a function of temperature. These shift factors were further fitted to the curve using the WLF model and the curve fitting result is also shown in Figure 9. Clearly, there also was a very good fit to these shift factors with the fitted constant values a1 = 208.9 and a2 = 1092.0. Based on the TTS principle, a single master curve could be constructed by shifting these frequency-dependent storage moduli at different temperatures along the time axis, as shown in Figure 8b, where the reference temperature was set to the T g of the EMC. The master curve could be well fitted by a Prony series equation using 22 Prony elements; the fitted weighting coefficients c k and the relaxation times ς k are listed in Table 1. The corresponding temperature shift factors in logarithmic scale are shown in Figure 9, as a function of temperature. These shift factors were further fitted to the curve using the WLF model and the curve fitting result is also shown in Figure 9. Clearly, there also was a very good fit to these shift factors with the fitted constant values a 1 = 208.9 and a 2 = 1092.0.

Verification of the Effective Models
To demonstrate the feasibility of the TM technique, a fraction of a Cu circuit layer was considered as a test vehicle. The fractional Cu circuit layer was modeled using both detailed FE modeling and the TM technique, and the results are presented in Figure 10. Noticeably, there was a high agreement between them, indicating that the TM technique could not only robustly but also precisely distinguish the Cu circuitries from the PP dielectric or SM material.
The effectiveness of the proposed TM/FEA effective method was verified by comparing the predicted effective orthotropic elastic properties of the fractional Cu circuit layer displayed in Figure 10 with the ROM/analytical estimate and with the FEA-based effective approach. The latter was considered a benchmark model. It was worth mentioning that for simplification, the orthotropic elastic material is simplified as a transversely isotropic elastic material by averaging the effective in-plane elastic moduli ( x E and y E ) and CTEs  Tables 2 and 3. These two tables illustrate that the effective elastic moduli and CTEs showed a strong temperature dependence, where an elevated temperature would considerably lessen the effective elastic moduli but enlarge the CTEs. Moreover, Table 2 shows that the calculated effective in-plane and out-of-plane elastic moduli by the proposed TM/FEA effective method were much more consistent with those of the FEA-based effective method over the temperature range of 25 °C 260 °C, as compared to the ROM/analytical estimate. On the other hand, the calculated effective in-plane elastic moduli by the ROM/analytical estimate deviated considerably from those of the other two effective approaches across the temperature range. Furthermore, a similar result could be also found for the predicted effective CTEs, as shown in Table 3. Similar to the effective elastic moduli, there is a very pronounced difference in the effective in-plane CTEs between the ROM/analytical estimate and the other two effective approaches.

Verification of the Effective Models
To demonstrate the feasibility of the TM technique, a fraction of a Cu circuit layer was considered as a test vehicle. The fractional Cu circuit layer was modeled using both detailed FE modeling and the TM technique, and the results are presented in Figure 10. Noticeably, there was a high agreement between them, indicating that the TM technique could not only robustly but also precisely distinguish the Cu circuitries from the PP dielectric or SM material.   The effectiveness of the proposed TM/FEA effective method was verified by comparing the predicted effective orthotropic elastic properties of the fractional Cu circuit layer displayed in Figure 10 with the ROM/analytical estimate and with the FEA-based effective approach. The latter was considered a benchmark model. It was worth mentioning that for simplification, the orthotropic elastic material is simplified as a transversely isotropic elastic material by averaging the effective in-plane elastic moduli (E x and E y ) and CTEs (α x and α y ) to be E x,y and α x,y , respectively. The calculated effective properties are shown in Tables 2 and 3. These two tables illustrate that the effective elastic moduli and CTEs showed a strong temperature dependence, where an elevated temperature would considerably lessen the effective elastic moduli but enlarge the CTEs. Moreover, Table 2 shows that the calculated effective in-plane and out-of-plane elastic moduli by the proposed TM/FEA effective method were much more consistent with those of the FEA-based effective method over the temperature range of 25-260 • C, as compared to the ROM/analytical estimate. On the other hand, the calculated effective in-plane elastic moduli by the ROM/analytical estimate deviated considerably from those of the other two effective approaches across the temperature range. Furthermore, a similar result could be also found for the predicted effective CTEs, as shown in Table 3. Similar to the effective elastic moduli, there is a very pronounced difference in the effective in-plane CTEs between the ROM/analytical estimate and the other two effective approaches. The orthotropic Cu circuit layer was further approximated as a transversely isotropic material by simply averaging the effective in-plane elastic moduli (E x and E y ) and CTEs (α x and α y ) as E x,y and α x,y , respectively. The transversely isotropic material could be further simplified as an isotropic material through an average of the effective in-plane and out-of-plane elastic moduli (E x,y ,E z ) and CTEs (α x,y , α z ) as E x,y,z and α x,y,z , respectively. The feasibility of the three constitutive models, i.e., orthotropic, transversely isotropic and isotropic, was further examined through FEA of the fractional Cu circuit layer shown in Figure 10 when subjected to a temperature load from 25 • C to 260 • C, and the calculated thermal deformations in the x-, yand zdirections are displayed in Table 4. For comparison, the detailed FEA results, serving as benchmark data, are also listed in the table, which shows that the effective orthotropic model presented the best consistency with the detailed FEA, followed by the transversely isotropic and isotropic models. This result matched with mechanical intuition. Accordingly, the effective orthotropic constitutive model was used in the subsequent warpage process simulation.

Thermal Analysis of the Interposer Bonding Process
The temperature distribution of the FCPoP assembly during the interposer bonding process may not be uniform across the assembly due to the uneven applied process thermal loading, which would cause a more excessive deformation. Thus, prior to conducting the warpage process simulation, the temperature distribution of the FCPoP assembly in natural convection during the interposer bonding process was characterized using a 3D transient heat conduction FEA. The natural convective heat transfer model proposed in [22] and the standard radiative heat transfer model [23] were applied to depict the natural convective and radiative surface heat transfer, respectively. According to the process condition, a preheat temperature was first set on the top and bottom surfaces of the assembly for 4 s, which were 185 • C and 145 • C, respectively, followed by a temperature increase up to 245 • C in 10 s on the top surface. The ambient temperature is 25 • C. The thermal analysis result at the end of the process is demonstrated in Figure 11. It is important to note that the Cu core solder balls were arranged in two to three rows around the periphery of the substrate and interposer. Evidently, the heat was conducted from the top coreless substrate to the bottom core interposer mainly by way of these Cu core solder balls; as a result, the part of the substrate and interposer adjacent to these periphery Cu core solder balls experienced a higher temperature. In addition, there was a significant temperature non-uniformity and gradient across the assembly. The characterized temperature distribution was imposed as a thermal load in the warpage process simulation for a better prediction accuracy.

Thermal Analysis of the Interposer Bonding Process
The temperature distribution of the FCPoP assembly during the interposer bonding process may not be uniform across the assembly due to the uneven applied process thermal loading, which would cause a more excessive deformation. Thus, prior to conducting the warpage process simulation, the temperature distribution of the FCPoP assembly in natural convection during the interposer bonding process was characterized using a 3D transient heat conduction FEA. The natural convective heat transfer model proposed in [22] and the standard radiative heat transfer model [23] were applied to depict the natural convective and radiative surface heat transfer, respectively. According to the process condition, a preheat temperature was first set on the top and bottom surfaces of the assembly for 4 s, which were 185 °C and 145 °C, respectively, followed by a temperature increase up to 245 °C in 10 s on the top surface. The ambient temperature is 25 °C. The thermal analysis result at the end of the process is demonstrated in Figure 11. It is important to note that the Cu core solder balls were arranged in two to three rows around the periphery of the substrate and interposer. Evidently, the heat was conducted from the top coreless substrate to the bottom core interposer mainly by way of these Cu core solder balls; as a result, the part of the substrate and interposer adjacent to these periphery Cu core solder balls experienced a higher temperature. In addition, there was a significant temperature non-uniformity and gradient across the assembly. The characterized temperature distribution was imposed as a thermal load in the warpage process simulation for a better prediction accuracy.

Warpage Process Simulation
The calculated temperature-dependent effective orthotropic elastic properties of the coreless substrate and core interposer using the TM/FEA effective approach are demonstrated in Tables 5 and 6, respectively. It was interesting to see that the core interposer was much stiffer than the coreless substrate across the temperature range; On the contrary, the CTEs of the coreless substrate tended to slightly greater than those of the core interposer. The warpage evolution of the FCPoP assembly during the fabrication process is shown in Figure 12a. It can be clearly observed in the figure that the process-induced warpage extensively varied with the process steps, and also showed a significant increase after the die process bonding, underfill curing, and interposer bonding processes. In addition, the

Warpage Process Simulation
The calculated temperature-dependent effective orthotropic elastic properties of the coreless substrate and core interposer using the TM/FEA effective approach are demonstrated in Tables 5 and 6, respectively. It was interesting to see that the core interposer was much stiffer than the coreless substrate across the temperature range; On the contrary, the CTEs of the coreless substrate tended to slightly greater than those of the core interposer. The warpage evolution of the FCPoP assembly during the fabrication process is shown in Figure 12a. It can be clearly observed in the figure that the process-induced warpage extensively varied with the process steps, and also showed a significant increase after the die process bonding, underfill curing, and interposer bonding processes. In addition, the maximum warpage occurred after the interposer bonding process; i.e., at around 653.7 µm, rather than after the fabrication process, i.e., at about 82.6 µm. The reason that the interposer bonding process created the maximum warpage was the lack of the EMC helping to resist the shear force caused by the global CTE mismatch between the substrate and interposer.     Figure 12b illustrates the simulated warpages during the increasing temperature process. For comparison, the in-line warpage measurement data are also listed in Figure 12, presented as an average value with a standard deviation (SD) (error bar). It was evident that the simulated warpages over the temperature range of 25-260 • C closely followed the measurement data. The fair difference in warpage between the measurement and simulation could be attributed to the uncertainty in the measured temperature-dependent Young's modulus and CTE of the EMC. Additionally, the increased temperature reduced the warpage probably due to the softening of the EMC when exposed to temperatures greater than the T g . The simulated and measured warpage contour plots of the FCPoP assembly at step 12 (after the mold cure process) are shown in Figure 13. Once again, there were very consistent results between the simulation and measurement. Table 7 illustrates the simulated and average measured residual warpages at 30 • C and 260 • C together with the smallest and largest values of the measured data shown in the bracket. It is clear that these simulation data fell in the respective ranges of the measured data, and in addition, the results of the process simulation were very comparable to those of the measurement, where the maximum warpage difference between them was only around 5%. Moreover, the residual warpage at 30 • C is nearly double that found at 260 • C, and due to this, the residual warpage at 260 • C was not considered in the subsequent parametric analysis.

Effects of Component CTEs
The influences of the effective CTEs of the EMC, core interposer and coreless substrate on the warpage of the FCPoP assembly at 30 °C were addressed. The parametric results of the effect of the EMC CTE are presented in Figure 14a. In the parametric analysis, the effective CTE of the EMC nominally varied from −10% to +10%. The figure shows that the EMC CTE had a minor impact on the residual warpages due to the relatively rigid substrate and interposer. Specifically, an increase in the EMC CTE somewhat decreased the residual warpage. This could be due to an increased EMC CTE reducing its local CTE mismatch with the core interposer and coreless substrate, thereby leading to a lessened residual warpage.
The effects of the effective CTEs of the core interposer and coreless substrate are investigated, and the parametric results are also displayed in Figure 14a. Note that the three effective CTEs (αx, αy, αz) shown in Tables 5 and 6 simultaneously underwent a 10% variation from the original value. The figure demonstrates that the effective CTEs had a significant impact on the residual warpage. Specifically, an increase in the effective interposer CTEs dramatically reduced the warpage, whereas there was a totally opposite trend  The influences of the effective CTEs of the EMC, core interposer and coreless substrate on the warpage of the FCPoP assembly at 30 • C were addressed. The parametric results of the effect of the EMC CTE are presented in Figure 14a. In the parametric analysis, the effective CTE of the EMC nominally varied from −10% to +10%. The figure shows that the EMC CTE had a minor impact on the residual warpages due to the relatively rigid substrate and interposer. Specifically, an increase in the EMC CTE somewhat decreased the residual warpage. This could be due to an increased EMC CTE reducing its local CTE mismatch with the core interposer and coreless substrate, thereby leading to a lessened residual warpage. for the effective substrate CTEs. This was principally due to the effective CTEs of the interposer being smaller than those of the substrate. This suggests that the increase in the interposer CTEs reduced the CTE mismatch with the substrate, thereby leading to a reduced residual warpage. Likewise, the result of the effects of the effective CTEs of the coreless substrate can be also explained in the same way.

Effect of Component Orthotropic Elastic Properties
The effects of the effective orthotropic elastic properties of the core interposer and coreless substrate on the warpage at 30 °C were considered. Similar to the parametric analysis of the effective CTEs, there was a 10% variation in these 9 independent effective elastic property data ( The effects of the effective CTEs of the core interposer and coreless substrate are investigated, and the parametric results are also displayed in Figure 14a. Note that the three effective CTEs (α x , α y , α z ) shown in Tables 5 and 6 simultaneously underwent a ±10% variation from the original value. The figure demonstrates that the effective CTEs had a significant impact on the residual warpage. Specifically, an increase in the effective interposer CTEs dramatically reduced the warpage, whereas there was a totally opposite trend for the effective substrate CTEs. This was principally due to the effective CTEs of the interposer being smaller than those of the substrate. This suggests that the increase in the interposer CTEs reduced the CTE mismatch with the substrate, thereby leading to a reduced residual warpage. Likewise, the result of the effects of the effective CTEs of the coreless substrate can be also explained in the same way.

Effect of Component Orthotropic Elastic Properties
The effects of the effective orthotropic elastic properties of the core interposer and coreless substrate on the warpage at 30 • C were considered. Similar to the parametric analysis of the effective CTEs, there was a ±10% variation in these 9 independent effective elastic property data (E x , E y , E z , υ xy , υ yz , υ xz , G xy , G yz and G xz ) shown in Tables 5 and 6. The parametric results are presented in Figure 14b. They indicated that increased effective elastic properties of the interposer and decreased effective elastic properties of the substrate would amplify the residual warpage. This was because the core interposer is stiffer than the coreless substrate due to its possessing greater effective elastic and shear moduli. The growth of the effective elastic shear and elastic moduli of the interposer tended to result in a more excessive shear force resulting from the CTE mismatch between the interposer and substrate, thereby causing a greater warpage. On the other hand, the structural rigidity of the FCPoP assembly increased with the increase of the effective elastic and shear moduli of the substrate, which thus led to a reduced warpage. The results totally differed from the effects of the effective CTEs of the interposer and substrate.

Effect of Component Thickness
The dependence of the warpage at 30 • C on the thickness of the core interposer, coreless substrate and EMC was examined. Similarly, the thickness variation is also ±10% from their original value. Figure 14c illustrates the parametric results, where the residual warpage would increase both with an increasing interposer thickness and with a decreasing substrate thickness. The results were very consistent with the effects of the elastic properties of the interposer and substrate, and the explanation for this is the same as stated in the pervious section. In the parametric analysis of the thickness effect of the EMC, parametrizing the EMC thickness would, in the meantime, change the height of the CCSBs. Before conducting the parametric study of the influence of the EMC thickness, the height effect of the CCSBs in the assembly without an EMC was first explored. The parametric study, which is not presented here due to limited space, suggested that the CCSBs' height had little impact for the residual warpage. As a result, the parametric results on the thickness effect of the EMC was barely affected by the CCSBs' height. The dependence of the residual warpage on the EMC thickness is also presented in Figure 14c, which shows that the residual warpage significantly decreased with the EMC thickness. This can be attributed to the structural stiffness of the assembly substantially increasing with the EMC thickness, thus resulting in a less residual warpage.

Conclusions
This research successfully conducted an effective and robust prediction of the warpage performance of an FCPoP assembly during the fabrication process through the proposed process simulation framework. In this framework, the temperature-dependence of the elastic properties of the components, as well as the viscoelastic behavior and chemical shrinkage of the EMC were taken into account in this investigation. The temperaturedependent elastic properties and viscoelastic properties were experimentally characterized. In order to improve the computational and modeling efficiency while also preserving good prediction accuracy, a novel effective approach; i.e., the TM/FEA effective method, was introduced to assess the effective elastic properties of the orthotropic coreless substrate