Finite element simulation of solid-liquid interdiffusion bonding process

—Solid–liquid interdiffusion (SLID) bonding ﬁnds a wide variety of potential applications toward die-attach, her-metic encapsulation of microelectromechanical systems (MEMS) devices and 3-D heterogeneous integration. Unlike soft soldering technique, the solidiﬁcation of intermetallic compound (IMC) formation in SLID bonding occurs during the process isother-mally, making it difﬁcult to predict and mitigate the sources of process-dependent thermomechanical stresses. Literature reports two dominant factors for the built-in stress in SLID bonds: volume shrinkage (due to IMC formation) and coefﬁcient of thermal expansion (CTE) mismatch. This work provides a detailed investigation of the Cu–Sn SLID bonding process by ﬁnite element (FE) simulations. Speciﬁcally, the FE simulation of the SLID bonding process is divided into three steps: ramp-up, hold-time, and ramp-down stages to understand the stresses formed due to each individual step. Plastic material properties for Cu as well as temperature-dependent material parameters for different entities are assigned. Process-dependent thermome-chanical stresses formed during the ramp-up and hold-time steps (IMC formation) were found not to be signiﬁcant. The hold-time step is governed by the reaction and diffusion kinetics, which determines the bond line quality including defects, such as voids. The ramp-down step is the dominant phase inﬂuencing the ﬁnal stress formations in the bonds. The results show an average of > 30% decrease in the stress levels in Cu 3 Sn layer (IMC) when the bonding temperature is brought down from 320 ◦ C to 200 ◦ C, thus demonstrating the importance of low-temperature SLID process.


Finite Element Simulation of Solid-Liquid
Interdiffusion Bonding Process: Understanding Process-Dependent Thermomechanical Stress Nikhilendu Tiwary , Member, IEEE, Vesa Vuorinen , Glenn Ross , Member, IEEE, and Mervi Paulasto-Kröckel , Member, IEEE Abstract-Solid-liquid interdiffusion (SLID) bonding finds a wide variety of potential applications toward die-attach, hermetic encapsulation of microelectromechanical systems (MEMS) devices and 3-D heterogeneous integration.Unlike soft soldering technique, the solidification of intermetallic compound (IMC) formation in SLID bonding occurs during the process isothermally, making it difficult to predict and mitigate the sources of process-dependent thermomechanical stresses.Literature reports two dominant factors for the built-in stress in SLID bonds: volume shrinkage (due to IMC formation) and coefficient of thermal expansion (CTE) mismatch.This work provides a detailed investigation of the Cu-Sn SLID bonding process by finite element (FE) simulations.Specifically, the FE simulation of the SLID bonding process is divided into three steps: rampup, hold-time, and ramp-down stages to understand the stresses formed due to each individual step.Plastic material properties for Cu as well as temperature-dependent material parameters for different entities are assigned.Process-dependent thermomechanical stresses formed during the ramp-up and hold-time steps (IMC formation) were found not to be significant.The hold-time step is governed by the reaction and diffusion kinetics, which determines the bond line quality including defects, such as voids.The ramp-down step is the dominant phase influencing the final stress formations in the bonds.The results show an average of >30% decrease in the stress levels in Cu 3 Sn layer (IMC) when the bonding temperature is brought down from 320 • C to 200 • C, thus demonstrating the importance of low-temperature SLID process.Index Terms-Cu-Sn, finite element (FE) simulations, heterogeneous integration (HI), solid-liquid interdiffusion (SLID) bonding, thermomechanical stress.

I. INTRODUCTION
H ETEROGENEOUS integration (HI) is receiving consid- erable attention for stacking components with different functionalities [logic, memory, application-specified integrated circuit (ASIC), microelectromechanical systems (MEMS), etc.] to continue the trend of faster, efficient, and cost-effective electronic products [1], [2].Small volume, high-density interconnects are of utmost importance for bonding and electrically connecting different components for the realization of HI [3], [4].The interconnect size and density should conform at every hierarchical level in a 3-D stack, resulting in a wide range of available interconnect size/pitch spanning from flip-chip solder bumps (≈100 μm) to μbumps (≈30 μm) to the hybrid bonding (<1 μm).Moreover, as the size of the solder bumps reduces, the bump almost entirely gets transformed into intermetallic compounds (IMCs) [5].
Solid-liquid interdiffusion (SLID) bonding is a technology, which relies on complete formation of IMCs in the bond line.The melting point of the resulting IMCs is higher than the processing/bonding temperature allowing subsequent processes at higher temperature.Thus, owing to its superior thermal stability, more stacks could be piled up making this process attractive for 3-D ICs [6].Moreover, SLID bonds are metallic bonds, which provides very good hermeticity, and thus, it also has great applications toward hermetic sealing of the MEMS devices with cap wafers [7].Therefore, SLID could offer potential solution 1) to fill the gap between μbump technology and hybrid bonding and 2) hermetic encapsulation of MEMS devices.Furthermore, it also finds potential applications toward die attach for high-power semiconductor devices [8], [9].Although there are different SLID techniques reported, such as Cu-Sn, Au-Sn, Ag-Sn, Ni-Sn-Cu, Cu-In, Cu-Sn-In, and so on [6], [10], the most popular and widely researched is the Cu-Sn system.The main reasons for that are the low cost compared to Au-Sn system, easy processing steps, lower resistivity of the Cu 3 Sn IMC, and easy integration with the Cu interconnects and Cu through silicon vias (TSVs) in the system [11].Thus, all the above aspects make SLID bonds very attractive and aligned with the HI approach.
Employing the SLID technique for 3-D HI is also accompanied with the formation of process-dependent thermomechan-ical stress, due to thermal mismatch of different layers and components at the bonding interfaces.These stresses impact the reliability of the bonds and can cause silicon cracking [12], [13], change in carrier mobilities of transistors [14], and on the timing and circuit performance [15].Thus, it is required to understand the origin of SLID process-dependent thermomechanical stress in detail.
Limited studies have been done so far focusing on the process-dependent residual stress formation during the SLID bonding process.Ladani [16] performed finite element (FE) simulations on a package with Cu-Sn SLID bonds incorporating silicon dies, underfill, die attach, substrate, and molding compound.The stresses ascribed to coefficient of thermal expansion (CTE) mismatch were studied due to the thermal cycle loading condition from −55 • C to 125 • C. It was found that the die and substrate thickness were the dominant factors influencing the stress levels at the interfaces and on Cu interconnects.Similarly, Xiong et al. [17] performed FE simulation of a 3-D IC package with Cu 3 Sn interconnects and Sn-3.9Ag-0.6Cusolders for thermal cycling from 218 K to 398 K to study the stress and their fatigue life.Contrary to Ladani [16], they found that the solder joint array design and interconnect materials influence the stresses, whereas the chip and substrate thickness have less effect.Apart from that, Cu-Sn SLID seal rings are simulated with defects to understand the crack initiation and propagation mechanisms through the bond line because of thermal cycling [18].However, in all the above studies [16]- [18], process-dependent residual stresses in Cu-Sn SLID bonds were not considered, but the geometrical and structural parameters of the package were considered.Ivankovic et al. [19] studied the impact of stress, formed due to Cu-Sn SLID bumps with no-flow underfill (NUF) material, on the performance of field effect transistors (FETs) on thinned silicon dies.However, only the cooling step from the bonding temperature of 250 • C to room temperature was considered in the model.In another work [20], W-filled TSV in a thinned silicon die was bonded to a thick silicon die via Cu-Sn SLID process.Although the entire manufacturing process was considered in the FE simulations by incorporating both heating and cooling steps, the impact of volume shrinkage during the formation of IMCs was not clarified.Taklo et al. [21] performed FE simulations for a Cu-Sn SLID bond frame to estimate the processdependent residual stresses.According to their model, stress generated due to volume shrinkage of IMCs far exceeds the stress due to CTE mismatch, and thus, stress due to the later was ignored in their study.They also did not consider the heating and cooling step in their models.However, stresses formed due to CTE mismatch could not be ignored, and entire processing step of Cu-Sn SLID bonding process should be analyzed in detail to have a comprehensive understanding of the sources of stress.
Since the stresses have a profound effect on the functionality of devices and systems reliability [13]- [15], [18], it is important to understand the stress builtup during each process step in the SLID bonding process.Furthermore, parallel information could be derived from literature regarding the dominant mechanism of stress builtup: volume shrinkage of IMCs and CTE mismatch.In this work, detailed FE simulations are performed considering entire SLID bonding manufacturing steps to understand their related effect on the formation of process-dependent thermomechanical residual stress, unlike previous studies.This work elucidates differing perspectives on volume shrinkage due to IMC formation and CTE mismatch on the built-in stress in SLID bonds.More specifically, Cu-Sn SLID bonding process is employed due to the easy availability of material parameters of the related IMC, i.e., Cu 3 Sn.Furthermore, the IMC layer in SLID bonds often consists of defects, such as voids.The IMCs are strong and hard material, and due to its low fracture ductility, they are sensitive to defects, unlike metals.In this work, voids are also incorporated in the IMC layer to understand their role impacting the reliability of the bonds.In the following, a brief overview of Cu-Sn SLID bonding process is presented, and then, the modeling methodology and results are presented in detail.

A. Cu-Sn SLID Bonding Process
In Cu-Sn SLID bonding process, Sn (low melting point metal) is deposited on Cu (high melting point metal) on both the top and bottom substrates [Fig.1(a)].The substrates could be anything depending on the application such as Si, interposers, direct bonded copper, and so on.The two substrates are then brought together in contact.The application of force and temperature greater than the melting point of Sn [Fig.1(b)] results in the melting of Sn and subsequent formation of IMC [Fig.1(c)].More specifically, when Sn melts, Cu starts to dissolve in the melt and when liquid Sn becomes supersaturated, the IMC nucleates.The growth of the IMC and the consumption of liquid then proceed through solid-state diffusion of Cu through the IMC layer.The melting point of the IMCs is much higher than the melting point of Sn, which is 232 • C. For the Cu-Sn system, depending on the typical bonding temperature (≈300 • C), the IMC phases that could form are Cu 6 Sn 5 : M.P.-415 • C and Cu 3 Sn: M.P.-676 • C [6].Out of those, Cu 3 Sn is the most stable phase with superior thermal, mechanical stability, and electromigration resistant [5].In this work, only Cu 3 Sn intermetallic is studied by incorporating it into the bond line.Fig. 1(a)-(c) shows the block diagram of the Cu-Sn SLID bonding process.

B. Modeling Approach
SLID bonding despite sharing some similarities with soft soldering process, such as the presence of liquid phase and IMCs formation, is fundamentally different [22].One of the main differences is the reversibility of the soldering process where remelt is possible.In contrast to that, SLID technology is not reversible due to the complete formation of IMCs in the bond line.To have a comprehensive understanding of the key parameters impacting stress formation, it is necessary to study the entire SLID processing steps in detail.Therefore, the SLID bonding process is divided in three steps, and FE simulations for each of those were carried out to understand the formation of stress.The different steps in the SLID bonding process after the wafers are brought into contact are: 1) ramp-up-when the temperature is ramped up to the bonding temperature above the melting point of Sn; 2) hold-time-when the temperature is held constant at the bonding temperature with the application of force for the desired bonding time, resulting in isothermal solidification of IMCs with complete consumption of liquid phase; and 3) ramp-down-when the temperature is ramped down from the bonding temperature to the room temperature.Fig. 1(d) shows the schematic of the three steps with the points under consideration for modeling and explained in detail in the following paragraphs.
1) Ramp-Up: When the temperature is ramped up, stresses are formed due to different CTE mismatches of the layers.Previous work has reported Cu undergoing plastic deformation in the case of thin films used in electronic interconnects [23], which implies assigning plastic material properties for Cu.Moreover, Sn melts and then dissolution of Cu into Sn commences toward the formation of IMCs.However, the thickness fraction of the IMC formed during the ramp up step is rather low [24].Therefore, the IMC layer is not incorporated at this stage, but the melting of Sn is accounted for in the model.
2) Hold-Time: In this step, the temperature is kept constant during the bonding time with the application of the bonding force.This step is governed by the formation of IMCs.As the temperature is constant in this step, there is no CTE mismatch among the layers.The bonding force serves two purposes: ensuring proper contact between the bumps by overcoming the topology difference across substrate and breaking the native oxides of the Sn layer.For Cu and Sn thicknesses of 6 μm and 2 μm, respectively, 1-h hold time at 320 • C and a bonding force of 14.4 MPa results in full Cu 3 Sn formation in the bond line [18].When the bonding temperature is reduced, the hold time could be increased to have full Cu 3 Sn formed in the bond line.Thus, the hold time depends on the bonding temperature, thickness ratio (Cu and Sn), and the desired IMC (Cu 6 Sn 5 , Cu 6 Sn 5 + Cu 3 Sn, and Cu 3 Sn).Therefore, in this step, time constraints are not included and the IMC of interest, i.e., Cu 3 Sn is directly incorporated in the model.Simulating the IMC formation is complicated as it is governed by diffusion kinetics and thermodynamical constraints [24].Moreover, when Cu 3 Sn is formed, there is also a volume shrinkage.So, in this step, the volume shrinkage of the Cu 3 Sn layer is incorporated in the model.Finally, the bonding pressure is included to estimate its impact on the stress formation.
3) Ramp-Down: This is the step after the hold-time, when the bonding is accomplished with complete formation of the desired IMC in the bond line, i.e., Cu 3 Sn.The heating is terminated with the release of bonding force, and then, the temperature is ramped down to the room temperature.As the temperature is reduced to the room temperature, different CTE mismatches among the layers will result in the formation of stress.Here again, plastic deformation in Cu cannot be ignored, and so, the elastoplastic material model is assigned to it.In this work, the time-dependent studies and the related effect, such as creep, are not considered.

C. Governing Equations
The thermal strain for a restricted layer is given by where ε th is the thermal strain, α is the CTE, and T and T ref are the steady state and reference temperature, respectively.Based on the above equation, the corresponding stress in the layer can then be expressed as where σ th is the thermal stress and E is the Young's modulus.
In our model, Cu layers were assigned the elastoplastic material model.The plastic strain rate given by the associative flow rule where the plastic strain rate .
pl is related to the partial derivative of plastic potential Q with respect to the stress S is given as where λ is the scalar plastic multiplier.In the associative flow rule, the plastic potential is the same as the yield function, Q = F.The yield function F was provided as von Mises given by where σ mises and σ ys are the von Mises stress and yield stress of the material, respectively.
A linear isotropic hardening model was employed, which is governed by the following equations: where σ ys0 is the initial yield stress and E Tiso is the isotropic tangent modulus of the material.Kinematic hardening was not employed in the model.

D. Material Parameters
The material parameters for all the materials were taken as the default values available in the COMSOL library.The parameters for which the values are not present in the library were gathered from literature.The material parameters employed are shown in Table I.The plastic material domain is only assigned to Cu, whereas linear elastic material domain is assigned for all other materials.The assumption of linear elastic for Si and Cu 3 Sn is valid because of their high yield strength: Si-7000 MPa [25] and Cu 3 Sn-1900 MPa [21].In the case of Sn, it could well be defined as linear elastic until its melting point.The yield stress values of electroplated Cu are reported to be in the range of 233 ± 15 MPa [26].It has also been found to depend on the processing conditions [27], structure, texture, grain size, and temperature [28], [29].The yield stress for electroplated Cu has been found to decrease from ≈250 MPa at room temperature to <100 MPa at 250 • C [29].Moreover, a tangent modulus of 1 GPa has been employed previously for Cu in FE modeling [30].In this work, a temperature independent initial yield stress (σ ys0 ) of 250 MPa and an isotropic tangent modulus (T iso ) of 2 GPa were used for Cu to account for its plastic material properties.The higher values were employed even though there would be some overestimation of stress in the study.The corresponding stress-strain plot with the above values is shown in Fig. 2 For Cu 3 Sn, Young's modulus and Poisson's ratio were assigned as 108.3 GPa and 0.299, respectively [31].

E. Model Schematic and Boundary Conditions
The FE simulations were carried out in COMSOL Multiphysics utilizing the solid mechanics physics interface.The target in this work is to evaluate the process-dependent local stresses on the individual SLID bonds and not global stresses at the wafer level.Therefore, it is not necessary to simulate the entire silicon wafer area with all micro bumps incorporated, which will be very complicated and computationally impractical to solve.The stresses formed in the silicon region also do not extend much beyond few tens of micrometers from the bonds [21].Again, it is not necessary to incorporate the entire thickness of silicon wafer in the model.The individual SLID bonds are a couple of micrometers thick, so keeping low thickness of Si substrate would ensure lower aspect ratios in the model and thus fewer mesh elements.Similarly, no adhesion layer (typically TiW or Cr located between Cu and Si) is included in the model, which would otherwise make the model computationally expensive.This is because the thickness of adhesion layers is in the range of few tens of nanometers.
Two separate schematics were used in the modeling process.One for the ramp-up step where the model consists of Si, Cu, and Sn layers [Fig.3(a)] and the other for the hold-time and ramp-down step, where it is assumed that full IMCs are formed [Fig.3(b)].It consists of both the top and bottom Si and Cu layers with Cu 3 Sn layer sandwiched between them.As explained above to keep the model simple, just one bump is considered, and the thickness of Si layers is assigned as 50 μm.This ensures minimum computational costs with reduced impact on the results.
A subnode of thermal expansion and plasticity was added under the linear elastic material node.The plasticity domain was assigned to Cu layer.A rigid motion suppression boundary condition boundary condition (BC) was applied to the bottom Si.It suppresses the rigid body motions in 3-D space with minimal constraints and so does not impose any additional stresses in the model.Since the motivation in this work is to analyze the thermomechanical stresses developed due to the CTE mismatch between the layers depending on the processing conditions, we have taken this BC as the default.An auxiliary sweep study was carried out by parametrizing the temperature, starting from the room temperature to the bonding temperature for ramp-up step and from the bonding temperature to the room temperature for ramp-down step.In both the cases, the initial stresses at the start temperature was zero.In this work, the stresses evaluated are the von Mises stress, which is a convenient way of characterizing a tridimensional stress state of a body with the yield limit equal to the uniaxial stress state of the material.Three bonding temperatures were studied, 200 • C, 250 • C, and 320 • C. The room temperature was the default value 20 • C.

A. Ramp-Up
The ramp-up simulations were carried out for three temperature values viz.200 • C, 250 • C, and 320 • C. Fig. 4(a) shows the schematic of the diagonal cut plane across the model, where the results are evaluated.A diagonal cut plane captures the stress profile at all regions of interest, i.e., corners and at halfway locations across the layers.Cut lines are shown at two locations: 2 μm from the corner (red) and at the middle position (black).The sharp corners are generally the stress concentration regions due to singularity at those points in FE analysis.To remove any artifacts due to those effects, the cut line plot (red) is evaluated 2 μm away from the corner of the bump.However, a major point to note from this step is that Cu goes into plastic deformation at this step, which is evident from the stress value in the Cu layer which is ≈250 MPa equal to its yield strength which was assigned.Fig. 4(h) shows the maximum equivalent plastic strain generated in the Cu layer at different temperatures.These plastic strains were found to be concentrated at the edges and corners of Cu layer near to the Si interface.It is interesting to note the two different linear plastic strain profiles with temperature: one below the melting point of Sn while the other above it.When Sn melts, it results in some strain relaxation in the Cu layer.This could also be seen from the red cut line plots in the Cu layer of Fig. 4(e) and (g), where the stress in Cu layer is more relaxed while approaching Sn layer as compared to Fig. 4(c).Thus, the maximum von Mises stress generated at this step is in the Cu layer limited by its yield strength.

B. Hold-Time
In this step, the wafers are already in contact after the rampup stage, and the bonding temperature is held constant for the desired bonding time.Thus, the CTE mismatch among the layers is ignored, and this step is governed by the formation of IMCs.Formation of IMCs results in volume shrinkage.In literature, there is a large scatter in the data regarding volume shrinkage when IMCs are formed.Different percentage values for volume shrinkage such as 2% [24], 5% [32], 10% [33], and 12.4% [21] are reported when Sn transforms to Cu 6 Sn 5 .Similarly, 2% [24], 7.5% [21], 8.2% [33], and 8.5% [32] are reported for volume shrinkage when Cu 3 Sn is formed.The density of the materials is temperature dependent, which could explain few variations in the reported volume shrinkage values.The volume shrinkage for Cu 3 Sn could be estimated from the following equation [34]: Considering the density at 320 • C for Sn, Cu, and Cu 3 Sn to be 7.2251 [Fig.8(A2)], 8.96, and 8.8313 g/cm 3 [Fig.9(A4)], respectively, and mass of Sn, Cu, and Cu 3 Sn to be 118.71,63.546, and 309.348 g, respectively; the V will be ≈7.1%.This is close to the other values reported in [21], [32], and [33].
The other aspect is how this volume contraction is incorporated in the model.Taklo et al. [21] incorporated ≈17% of volume shrinkage in their model in the IMC layers and found the stresses to be extremely high (≈GPa), even exceeding the levels resulting from CTE mismatch.Fig. 5(a) shows  the cross-sectional SEM micrograph of a Cu-Sn SLID bond manufactured in our lab with Cu 3 Sn IMC formed and typical Sn squeeze out.The squeeze out is apparent in the two axial directions (x-and y-axes).Incorporating the volume shrinkage percentage by resolving it in terms of initial strain in the axial directions inside the IMC layer results in extremely high levels of stresses (≈GPa), which is unrealistic.Moreover, the growth of IMCs in the perpendicular direction (x-and y-axes) of the interdiffusion direction (z-axis) is much slower [35].Therefore, the volume shrinkage could be assumed to take place along the thickness, i.e., z-axis.Since the wafers are also free to move in the z-axis during the process, this assumption is realistic.Li et al. [36] calculated volume shrinkage in micro Ni/Sn/Ni joints based on the height measurements of the sandwiched structure.The volume shrinkage calculations based on the height measurements (along the z-axis) including the void fraction has been shown to match well with the molar volume calculations [36].
Therefore, the volume shrinkage condition is applied uniaxially in the model along the z-axis for the bump.The total volume shrinkage percentage then transforms entirely to the percentage strain along the z-axis.Thus, an initial strain condition of 7.1% was incorporated along the z-axis in the Cu 3 Sn layer, and the model was solved under two BCs: 1) without and 2) with bonding force.In the first case, rigid motion suppression BC was applied, and in the later, bonding pressure of 20 MPa was applied on the top Si layer, while the bottom Si layer was given a fixed constraint BC.  and at the middle (black) demonstrate that the maximum stress generated in the bump decreases from the corners (<140 MPa) to the middle (<40 MPa) of the bonding area, Fig. 5(f).Once the bonding pressure is released, the stresses due to that would also be released, so the stresses formed due to the bonding pressure are insignificant.
Therefore, in this step, the stresses formed due to the actual process parameters are insignificant.Indeed, this step is purely dominated by IMC growth, which are governed by diffusion kinetics model [24].Thus, the hold-time step decides the quality of bonds and the bond-line but does not have significant effect on the residual thermomechanical stress formations in the IMC layer based on the physical process parameters.

C. Ramp-Down
This is the step where the bonding is complete with the full formation of Cu 3 Sn in the bond line.Subsequently, the bonding pressure is released, and the temperature is ramped down to the room temperature.Therefore, stresses due to CTE mismatch will be evident.Since, the Cu 3 Sn formation is often associated with defects, such as voids; we have performed this study considering two conditions: with and without voids.
1) Without Voids: Fig. 6 The stress in the Si layer also increases with the bonding temperature.Moreover, the stress in Si layer at Si-Cu interface increases beyond 250 MPa for the bonding temperature of 320 • C.These stresses in the Si layer are also expected to have some effect on the adhesion layer, which we have not included in our model.Thus, a poor adhesion layer quality could easily become a source of adhesion-related failures if it is not able to adjust to the resultant stress level.Fig. 6(g) shows the maximum equivalent plastic strain plot as a function of temperature.Higher plastic strains are generated in this step (≈0.023 max) as compared to the ramp-up step [Fig.4(h)] (≈0.010 max), implying this to be the dominant step in the residual stress formations in the bonds.
2) Voids at the IMC: The SLID bonds are often accompanied by voids impacting their reliability.Several independent and coupled phenomena are responsible for the formation of voids, such as Kirkendall effect, volume shrinkage, electroplating parameters, impurities in the electroplating solutions, crystallinity of the electroplated films, bonding pressure, stress relaxation, and so on [36]- [41].These voids could form because of the above factors during the hold-time step or during the ramp-down step.The voids are observed to be randomly distributed across the IMCs or mostly at Cu/Cu 3 Sn interface [37], [42].When the bonding pressure is low, voids have also been found to be concentrated at the midplane of the IMC layer [40].A large percentage of void formation could significantly impact the reliability of the bonds.Efforts are there to understand the origin of voids and suppress it by using high-purity Cu electroplating solutions [39] and increasing bonding pressure [40].Nevertheless, its formation could not be ignored, and the ramp-down simulation was repeated by incorporating voids in the model.The size of the voids could vary from few nanometers to few micrometers [43].Moreover, small voids could also coalesce together to form larger voids with random shapes [44].To keep the model simple, the size of the voids is taken as a sphere of radius 0.5 μm.Fig. 7 shows the results along the diagonal cut plane of the model, with the voids incorporated at the middle of the bonding interface for the bonding temperature of 320 • C and 200 • C. The stress concentration across the voids could be observed from the zoomed-in image and appears to be significantly reduced for the lower bonding temperature.Increasing the mesh elements by refining it increases the stress values asymptotically at the void boundaries.This is due to the limitations of the FE techniques, which create singularity at those points while evaluating stress.Moreover, crack initiation and propagation spots due to the stress concentration across the voids because of thermal cycling are well known and have been reported in the past [18].
Herein, stresses were also evaluated across cut lines 1 μm away from the voids toward the left.The two cut lines: red (left void) and black (middle void) are at 1 μm away from the voids (toward left), and a slight notch near the void position is apparent in the cut line plots [Fig.7(b) and (d)].These cut line plots were tested with different mesh sizes and were found to be consistent.The large stresses builtup in the intermetallic layer also tends to get compensated by the formation of voids, which are reported in the literature [41].The results were also consistent with the position of the voids across the IMC layer and just the notch (stress relaxation) shifts with the void position in the stress plots.Thus, in stationary case, the formation of voids supplements some stress relaxation in IMCs.In the dynamic case, e.g., thermal cycling, the voids are also the probable failure spots in the IMC layer posing significant reliability issues.
Therefore, the ramp-down simulation results demonstrate that the cooling is the major step, which will have the most effect on residual stress builtup in the bonds.Higher the bonding temperature, higher the stresses would form on the bonds which is quite intuitive.On the other hand, a lower bonding temperature in the case of Cu-Sn would then lead to a larger hold-time for complete Cu 3 Sn formation in the bond line.Therefore, to have an optimal stress level in the bonds, there is a tradeoff between the processing time and the temperature.

IV. CONCLUSION
In this work, detail FE simulations of the SLID bonding process are performed by dividing the entire bonding process in individual steps and simulating each of those.The stresses formed due to the actual process parameters in the ramp-up and hold-time step is not significant.The ramp-down step was found to be the one largely influencing the overall stress level in the bonds.Thus, the major contribution to the thermomechanical stresses in SLID bonds is due to CTE mismatch and not volume shrinkage during the formation of IMCs.The results show an average of >30% decrease in the von-Mises stress level in Cu 3 Sn bonds when the bonding temperature is brought down from 320 • C to 200 • C.This brings the importance of low-temperature SLID (LT-SLID) process, wherein the bonding temperature is reduced with the incorporation of other lower melting temperature material like In (melting point: 156.6 • C) [10].The lower bonding temperature reduces the thermomechanical stresses on the bonds, which, in turn, reduces the adverse impact of the residual stresses on the system reliability with different stacked modules.This FE simulation approach could be repeated on any SLID process concerning different metals and with different combinations of Si and cap materials to have a qualitative and quantitative understanding of the stress levels.The requirement will be on the material properties of the metallurgies involved.Furthermore, the methodology presented could also be employed in investigating the time-dependent studies, for which the requirement will be on time-dependent material parameters.APPENDIX See Figs. 8 and 9.

Fig. 1 .
Fig. 1.Schematic of the Cu-Sn SLID bonding process.(a) Cu-Sn stacks deposited on top and bottom Si wafer, (b) bonding process, (c) bonded pair, and (d) modeling methodology: different steps in SLID bonding process.

Fig. 2 .
Fig. 2. (a) Stress-strain plot for Cu fed in the model and (b) E versus T for Sn, E is set to 10 −5 for T > 232 • C.
(a).Various temperature-dependent material properties from the COMSOL library are employed to mimic the actual situation as close as possible.To mimic the liquid Sn during the rampup step, the temperature-dependent Young's modulus plot was modified by assigning it to ≈0 (10 −5 to avoid singularity) at the melting point of Sn, i.e., 232 • C. Fig. 2(b) shows the corresponding plot of E versus T for Sn.Constant extrapolation was employed for the data beyond the temperature range and is shown in all the temperature-dependent plots with red dashed line.All the other temperature-dependent material parameter plots for Sn and Cu 3 Sn are shown in Appendix.

Fig. 4 (
b), (d), and (f) shows the von Mises stress plots for the three temperatures under investigation along the diagonal cut plane of the model.The corresponding line graphs at the cut locations are shown in Fig. 4(c), (e), and (g) with Si, Cu, and Sn regions marked.As 200 • C is lower than the melting point of Sn (232 • C), some stresses in Sn layer are evident from Fig. 4(b) and (c).The effect of the low CTE mismatch between Cu and Sn layers and the high CTE mismatch between Cu and Si layer can be clearly seen from the corresponding black cut line plot through the middle location of the bump [Fig.4(c)].Moreover, the stress at the edge (red) in Cu layer gets relaxed as it approaches the Sn layer.On the other hand, Fig. 4(d) and (f) shows the von Mises stress plots for ramp-up temperature of 250 • C and 320 • C greater than the Sn melting point.As is evident, the stresses in the Sn layer are zero.This could also be more properly visualized from the corresponding line graphs at the two marked locations [Fig.4(e)-(g)].Thus, our model mimics the actual melting of Sn during the bonding process, and once it is liquid, the stresses in the Sn layer are absent.

Fig. 4 .
Fig. 4. (a) Schematic showing the diagonal cut plane and the cut lines.von Mises stress plots for the ramp-up step along the diagonal cut plane for (b) 200 • C, (d) 250 • C, and (f) 320 • C, the corresponding cut line graphs for (c) 200 • C, (e) 250 • C, and (g) 320 • C, and (h) plot of plastic strain with temperature.

Fig. 5 .
Fig. 5. (a) Typical Cu-Sn SLID bump with Sn squeeze out.The dimensions of the μbump are ≈100 μm × 100 μm × 10 μm.(b) Schematic showing the diagonal cut plane and the cut lines.von Mises stress plots for the holdtime step with volume contraction of 7.1% along the z-axis without applying bonding pressure, (c) along the diagonal cut plane, (d) corresponding cut line graphs, and with 20 MPa bonding pressure, (e) along the diagonal cut plane, and (f) corresponding cut line graphs.

Fig. 5 (
b) shows the schematic of the diagonal cut plane along which the results are evaluated.The red and black cut lines are marked at locations 2 μm from the corner and at the middle location.Fig.5(c) and (d) shows the results along the diagonal cut plane for the case when no bonding pressure BC is applied, and the resultant stresses (≈10 −7 MPa) are negligible.On the other hand, when the bonding pressure of 20 MPa is applied, stresses ≈240 MPa are seen at the corner of the bumps [Fig.5(e)].These stress levels are below the yield stress value for Cu (250 MPa), which was fed in the model.Since Cu already yielded during the ramp-up stage, it will not have any dominant effect on the bump corners.The cut line plots at the location of 2 μm from the corner (red)
(a), (c), and (e) shows the ramp down simulation results for the three bonding temperature values: 200 • C, 250 • C, and 320 • C, along the diagonal cut plane of the model.The corresponding cut line graphs at the two locations [shown in Fig. 5(b)], 2 μm from the corner and at the middle, are shown in red and black colors, respectively [Fig.6(b), (d), and (f)].Here again, Cu goes into yielding due to plastic deformations with the stress value equal to its yield strength, i.e., 250 MPa.The stress in the Cu 3 Sn layer increases with the bonding temperature.At 200 • C, the stress in the Cu 3 Sn layer is near to the yield strength of Cu ≈ 256 MPa and it increases to ≈383 MPa at 320 • C. Thus, there is ≈33% decrease in the stresses generated in Cu 3 Sn layer at the bonding temperature of 200 • C as compared to the bonding temperature at 320 • C [Fig.6(b) and (f)].

Fig. 7 .
Fig. 7. von Mises stress plots for the ramp-down step with voids along the diagonal cut plane for (a) 320 • C and (c) 200 • C and the corresponding cut line graphs for (b) 320 • C and (d) 200 • C.

TABLE I MATERIAL
PROPERTIES INCORPORATED IN THE MODEL