Simulation of TSV Protrusion in 3DIC Integration by Directly Loading on Coarse-Grained Phase-Field Crystal Model

: As thermal management in 3DIC integration becomes increasingly important in advanced semiconductor node processes, novel experimental and modeling approaches are in great demand to reveal the critical material issues involving multiscale microstructures that govern the behavior of through-silicon-via (TSV) protrusion. Here, a coarse-grained phase-ﬁeld crystal model properly coupled with mechanics through the atomic density ﬁeld is used to simulate the formation of polycrystalline structures and protrusion of nano-TSVs from the atomic scale. TSVs with different grain structures are directly loaded, and protrusion/intrusion proﬁles are obtained along with displacement, stress, and strain ﬁelds. Thermodynamic driving forces from external loadings and the mismatch of Young’s modulus between adjoining grains as well as detailed displacement and strain distributions are ascribed to control the complex deformation in TSVs. TSVs with sizes up to around 30 nm and an aspect ratio of 4 are successfully investigated, and a further increase in the size and aspect ratio to cover the micrometer range is feasible, which lays down a solid basis toward a multiscale material database for simulation inputs to the design of TSV-based 3DIC integration and relevant electronic design automation (EDA) tools.


Introduction
As the main driving forces nowadays for electronics, big data, AI, IoTs, and their convergence have pushed the semiconductor industry to face unprecedented challenges in the past decade. With the semiconductor advanced node progressing to 3/2 nm and beyond, advanced packaging, i.e., 2.5D and 3D packaging, enables more chips in a system and makes spaces between chips much closer, thus achieving a even higher density of integration [1,2]. Interconnects for 3D ICs, e.g., through-silicon-vias (TSVs), serving as connections through the vertical direction among different functional chips or ICs in 3D packages, play an increasingly important role in electronic devices. The TSV-based technologies can be used at different levels, e.g., from transistors to chips and to board, and they can empower devices to exhibit desirable features such as lower delay, lower consumption, and higher performance. Along with parallel development in size miniaturization, manufacturing processes of TSVs have been improved to accommodate the scaling of semiconductors and satisfy the demands for high performances at the same time. Novel processes of 3D packaging technologies utilizing TSVs have been invented, such as CoWoS from TSMC [3], Foveros from Intel [4], and X-stacking from Changjiang Storage [5].
Despite the above-mentioned advantages, integration schemes based on TSVs have also brought about reliability concerns. Higher integration density often means difficulties in heat dissipation, leading to thermal stresses and consequently plastic deformation in TSVs [6,7]. The plastic deformation drives atom diffusion, dislocation motion, voids, and crack formation, which may obstruct the electrical connection in devices. TSV protrusion, in particular, destroys the redistribution layer (RDL) and leads to chip failures [8,9]. In addition, the material factor becomes critical and can no longer be ignored in the miniaturized interconnects in 3D ICs [10,11]. The complexity of microstructure coupled with thermal stresses makes the behavior of protrusion or intrusion of TSVs unpredictable [12,13]. The plastic deformation in polycrystalline TSV strutures is determined by the defect motion, including the motion of vacancies, dislocations, and grain boundaries (GBs). Both a highly heterogeneous distribution of defects in TSVs, and a variety of defect motion mechanisms, e.g., vacancy diffusion, dislocation slip, climb, and grain boundary motion, make contributions to the complexity and unpredictable nature of TSV protrusion or intrusion. Therefore, 3D microelectronic packaging technologies in the near future have to be upgraded to keep pace with the post-Moore era by taking the material factor into consideration.
To find the root cause of TSV protrusion, experimental methods such as Micro-Raman Spectroscopy, X-ray diffraction used in conjunction with a cross-section TEM, or a synchrotron radiation source have been used to identify the stress states in TSVs. In addition, protrusion profiles were found to show different surface characteristics. The influence of annealing conditions on the protrusion profiles was recently examined by Zhao [14,15], in which the protrusion profiles were classified into global and local protrusion according to the variation of grain structures after annealing at different temperatures. Due to the impact of stresses from the adjoining TSVs, the deformation behavior of TSVs is also influenced by the pitch distance. Jalilvand [16] reported that the protrusion morphology is featured by a more annular shape in small-pitch TSVs and more granular shape in large-pitch TSVs. Jalilvand [16] also found that the overall stress in small-pitch TSVs is larger than that in large-pitch counterparts.
However, the overall picture on the atomic scale microstructure as well as the stress and strain distributions along the depth direction of the whole TSV remain unclear. As a result, different simulation methods have been developed to study the stress and strain distributions in TSVs and to understand the connection between the TSV protrusion morphology and microstructure. The traditional finite element method (FEM)-based mechanical analyses, defining TSVs as homogeneous block materials and adopting constitutive laws of bulk materials, can provide the distribution of residual stresses [17,18]. In Liu's work [19,20], the phase-field crystal (PFC) model, based on the classical density functional theory, was used to simulate the microstructures in the TSVs at the atomic scale. To study the characteristics of protrusion morphologies, a model based on convolutional neural networks (CNNs) was developed by Jalilvand [21] to discriminate different classes of protrusion profile. The complex amplitude phase field crystal model, hereinafter referred to as the APFC model, was developed by a complex amplitude expansion based on the PFC model by Goldenfeld [22,23]. Furthermore, Skaugen [24,25] derived and studied the crystal plastic behavior in the PFC model from the perspective of mechanics, making the stresses inside the crystal linkable to the atomic density field.
A method capable of directly loading on the TSV with microstructure simulated by the APFC model has been developed in this study based on the principles outlined by Skaugen [24] and Salvalaglio [26]. This method enables, first, an increase in the aspect ratio and size of the simulated TSVs, and second, direct mechanical loading on the boundaries of TSVs. Finally, the stress and strain states in the TSVs and the resulting protrusion characteristics along with the complex amplitudes and atomic density field can also be obtained simultaneously. Therefore, it is expected that the model and results presented in this study are beneficial to further understand the mechanisms behind TSV protrusions and intrusion, and pave a way, at the same time, for scale-bridging modeling and property predictions from atomic to nano and from nano to micron scales for 3D TSV-based interconnects.

The APFC Model
The APFC model is utilized to simulate TSV protrusion. For the purpose of completeness, the predecessor of the APFC model, i.e., the PFC model, is briefly introduced first. The advantage of the PFC model lies in that it can naturally incorporate the elastic and plastic physics from the atomic scale to the continuum field. An order parameter, i.e., the density field ψ, is introduced to describe the crystalline solid state in the system, and ψ can be expressed as follows [26][27][28]: where A j s are the complex amplitudes of plane waves, q j s are the reciprocal lattice vectors describing the crystal structure, r is the position vector, and c.c. is the acronym for complex conjugate. ψ 0 is the average density field and set to zero in this study. The free energy of the system can be well defined as the functional of ψ and expressed by the following equation [26,[29][30][31]: Since ψ is a conservative field, its kinetics can be obtained by solving the following Cahn-Hilliard type of equation [26]: where M is the mobility of the interface, and δF/δψ is the first variation of the free energy functional F with respect to ψ. By the complex amplitude expansion of the density field, the free energy functional can be formulated directly in terms of A j and A * j . Then, the free energy functional F ψ in Equation (2) is replaced by F A in the APFC model [22,23,32]: where g j = ∇ 2 + 2i q j · ∇, and Φ = 2 ∑ N j=1 |A j | 2 . ν, ∆B 0 , and B x 0 are parameters and set to 1/3, 0.42, and 0.98, respectively. Polynomial f S ({A j }, {A * j }) is related to A j and A * j , and in our study, it reads f S = −2∆T(A 1 A 2 A 3 + c.c.), where ∆T is set to 1.2. The evolution of amplitudes A j can also be obtained by solving the governing equation [22,23,26]: A set of the fourth-order partial differential equations (PDEs) are obtained, and the details on order reduction are referred to Salvalaglio [31,32]. After solving the equations, the time-dependent evolution of the density field ψ can be reconstructed by using Equation (1). In this study, crystals with a triangular lattice are considered, which can be described by: with k 0 = 1.

Initial Condition
The initial conditions (I.C.s) for the complex amplitudes A j in different nuclei placed in the TSV are set by the following equation [26,32]: where θ controls the grain orientation and φ 0 = . To obtain a diversity of polycrystalline structures in the TSV, different combinations of numbers, radii, and locations of nuclei are generated as the ICs for solidification. Then, after solidification, TSVs filled with different grain structures are relaxed for a certain time to equilibrate the structures before external loadings are applied to study the protrusion behavior.

Deformation and Boundary Condition
In addition to solving Equation (5), it requires computing an additional smooth distortion u δ to fulfill the condition of mechanical equilibrium, i.e., ∇ · σ = 0. The stress field σ can be decomposed into two parts σ ψ and σ δ . σ ψ ij is derived from the density field ψ, and then it can be expressed by A j [24][25][26]: To find σ, the Airy function χ is introduced, which is given by [26]: where κ = λ/[2(λ + µ)], with λ and µ representing the two Lamé coefficients and ij representing the Levi-Civita symbol. Then, σ can be calculated by [26]: Once σ δ is obtained, the corresponding strain field can be determined by [26]: According to [26], the deformation field u can be determined through a Helmhotz decomposition into curl-and divergence-free parts : Once the strain field ε is determined, ϕ and α can be computed by solving the following Possion equation [26]: Then, the deformation field can be obtained by Equation (12).
To simulate the protrusion behavior of TSVs, an external loading T along the left and right boundaries is directly applied through boundary condition (B.C.) settings. According to the relationship between T and the stress field as shown in Figure 1, we derive: where T x = Tcosβ and T y = Tsinβ.

Model Implementation
The simulations reported in this study are performed by integrating COMSOL Multiphysics (hereinafter referred to as COMSOL) with MATLAB. All the PDEs are solved using the FEM and implemented in COMSOL, and postprocessing and visualization of the results are conducted in MATLAB. In order to match the geometry of real TSVs after the dry reactive ion etching (DRIE) or RIE process [33,34], a trapezoidal geometry is adopted in this study, as shown in Figure 1. Note that all the variables in the model are dimensionless. Enforcing an external loading to the TSV is implemented in COMSOL by applying constraints and weak contributions to the left and right boundaries of the TSV. The constraint related to T x is directly implemented as a "constraint". The principles behind the weak contribution can be briefly introduced as follows. According to the variational principle [35], (16) where λ is the Lagrangian multiplier.
where g represents the constraint that needs to be applied to the boundaries. In this case, the constraint related to T y needs to be applied by a "weak contribution": For a better understanding of the present work, an overview of the APFC model implementation with cross-referenced mathematical equations is shown in Figure 2.

Microstructure Formation in TSV
Simulation of solidification in the TSV requires solving six equations related to A j and A * j in Equation (5). Nuclei with a triangular lattice are placed inside the TSV and initialized according to Equation (7). The atom density fields of the initial circular grains, with radii setting to the same value of 15, can be obtained by the relationship between the atom density field ψ and the complex amplitudes A j and A * j according to Equation (1). The snapshots of the microstructures in the TSV at t = 0, 50, and 180 during solidification are shown in Figure 3. At the beginning, the TSV contains only six nuclei with different orientations, and the remaining area is filled with supercooled liquid, as shown in Figure 3a. As solidification proceeds, the liquid in the TSV is gradually consumed, and finally, the inner cavities are filled with crystalline solids when all the liquid is consumed up. To obtain an explicit picture on the grain structure, the defects atoms, which are identified as the atoms missing one or more of the first nearest neighboring atoms, are outlined by red closed curves, as shown in Figure 3d. Figure 4 shows the number of atoms in the growing nuclei with different orientation θ at t = 0, 10, and 20 during solidification. Comparing the growth rates of those grains with θ ranging from −40 • to 40 • , it is found that the growth rates are not constants but symmetric about the line θ = 0. Furthermore, the growth of the nuclei also depends on the parameter ∆T in the APFC model, which relates to the degree of undercooling. If ∆T is smaller than 1.2, the nuclei inside the TSV will not grow but be absorbed back into the surrounding liquid and vanish. Figure 3. Snapshots of solidified microstructure in the TSV by placing six initial nuclei with different orientations at a time of (a) t = 0; (b) t = 50; and (c) t = 180; with (d) highlighting defect atoms in (c) by red closed curves. Figure 5 shows two different grain structures, obtained by setting I.C.s with the same number, size, and location of the nuclei, but the misorientations between adjacent nuclei are set to ≤10 • and ≥15 • for Figure 5a,b, respectively. In Figure 5a, there exists only a few defect atoms, which cannot form continuous GBs. However, in Figure 5b, it can be observed that GBs are formed between adjacent grains after solidification, featuring grains with different orientations properly wrapped by the GBs. This difference in microstructure is simply caused by the above-mentioned orientation setting in the ICs. An explanation can be given that grains with different orientations compete to grow with their neighboring grains during solidification, forming a complex nework of GBs. Once the grains come in contact and join with each other, their orientations gradually converge to the same, resulting in dislocation annihilation. When misorientations between grains are smaller, almost a single crystal can form upon solidification completes, featuring a few defect atoms scattered in the TSVs. However, with misorientation becoming larger, coordination and orientation convergence between neighboring grains become more difficult and thus result in the formation of a polycrystalline structure in the TSV.  Figure 5 shows two different grain structures, obtained by setting I.C.s with the same number, size, and location of the nuclei, but the misorientations between adjacent nuclei are set to ≤10 • and ≥15 • for Figure 5a,b, respectively. In Figure 5a, there exists only a few defect atoms, which cannot form continuous GBs. However, in Figure 5b, it can be observed that GBs are formed between adjacent grains after solidification, featuring grains with different orientations properly wrapped by the GBs. This difference in microstructure is simply caused by the above-mentioned orientation setting in the ICs. An explanation can be given that grains with different orientations compete to grow with their neighboring grains during solidification, forming a complex nework of GBs. Once the grains come in contact and join with each other, their orientations gradually converge to the same, resulting in dislocation annihilation. When misorientations between grains are smaller, almost a single crystal can form upon solidification completes, featuring a few defect atoms scattered in the TSVs. However, with misorientation becoming larger, coordination and orientation convergence between neighboring grains become more difficult and thus result in the formation of a polycrystalline structure in the TSV. Figure 6 shows the microstructures at t = 0, 100, and 320 during solidification in a TSV with a height around 30 nm and an aspect ratio of 4, which is three times as large as the size studied in Liu's work [19]. Increasing the size and aspect ratio of the simulated TSV is meaningful for the design of 3DIC, where a multiscale material database is required for simulation inputs [36]. In this context, the APFC model presented in this study can play a role as the size of TSV-based interconnects is envisioned to span the whole spectrum from nanometers to micrometers for near future applications.  Figure 6 shows the microstructures at t = 0, 100, and 320 during solidification in a TSV with a height around 30 nm and an aspect ratio of 4, which is three times as large as the size studied in Liu's work [19]. Increasing the size and aspect ratio of the simulated TSV is meaningful for the design of 3DIC, where a multiscale material database is required for simulation inputs [36]. In this context, the APFC model presented in this study can play a role as the size of TSV-based interconnects is envisioned to span the whole spectrum from nanometers to micrometers for near future applications.  After solidification completes, a stress field σ ψ induced by the defects inside the TSV can be calculated according to Equation (8). Figure 7b plots σ ψ yy present in the microstructure shown in Figure 7a. It can be clearly observed that the stress concentration occurs near defect atoms as well as at the walls of the TSV. Noted that the magnitude of σ ψ is around 0.25 and substantially smaller compared with σ δ caused by external loadings, as discussed in the following sections.
After solidification completes, a stress field σ ψ induced by the defects inside the TSV can be calculated according to Equation (8). Figure 7b plots σ ψ yy present in the microstructure shown in Figure 7a. It can be clearly observed that the stress concentration occurs near defect atoms as well as at the walls of the TSV. Noted that the magnitude of σ ψ is around 0.25 and substantially smaller compared with σ δ caused by external loadings, as discussed in the following sections.

Effect of Microstructure
As illustrated in Figure 1, a symmetric external loading is applied as a B.C. to simulate TSV protrusion, while the displacement u y at the TSV bottom is fixed to be zero. Consisting of polycrystalline structures, TSVs exhibit anisotropy in Young's modulus. Moreover, the atomic-scale structure and its evolution, in particular, the evolution of the grain structure and the motion of defect atoms, make the mechanisms behind the mechanical deformation of TSV difficult to understand. In the following, the Young's moduli of the grains are set according to their crystallographic orientation. Otherwise, a uniform Young's modulus across the whole TSV makes the effect of grain structure on the protrusion profile of TSV negligible.
Sketched by the gray curves, the protrusion profiles of TSVs with initial microstructures shown in Figure 5 and subjected to the same external loading with a magnitude of 8 × 10 7 and β = π/3 are plotted in Figure 8a,b. It can be clearly observed that the protrusion profiles of the two TSVs are different with the average protrusion height in Figure 5a being greater. Note that essentially there is only one grain with a few scattered defect clusters present in Figure 5a, while a polycrystalline structure formed in Figure 5b. In terms of deformation easiness, the TSV shown in Figure 5b contains more grains and GBs, and thus, its deformation is less easy. Figure 9 shows the initial microstruture in a larger size TSV together with its protrusion profile after applying the external loading. It can be observed that dislocations are always hindered by GBs, which means increasing the number of grains leads to less deformation. Meanwhile, grain interiors are usually more susceptible to deformation than GBs, so once GB motion is switched on, protrusion in large grains is greater than in small grains.

Effect of Microstructure
As illustrated in Figure 1, a symmetric external loading is applied as a B.C. to simulate TSV protrusion, while the displacement u y at the TSV bottom is fixed to be zero. Consisting of polycrystalline structures, TSVs exhibit anisotropy in Young's modulus. Moreover, the atomic-scale structure and its evolution, in particular, the evolution of the grain structure and the motion of defect atoms, make the mechanisms behind the mechanical deformation of TSV difficult to understand. In the following, the Young's moduli of the grains are set according to their crystallographic orientation. Otherwise, a uniform Young's modulus across the whole TSV makes the effect of grain structure on the protrusion profile of TSV negligible.
Sketched by the gray curves, the protrusion profiles of TSVs with initial microstructures shown in Figure 5 and subjected to the same external loading with a magnitude of 8 × 10 7 and β = π/3 are plotted in Figure 8a,b. It can be clearly observed that the protrusion profiles of the two TSVs are different with the average protrusion height in Figure 5a being greater. Note that essentially there is only one grain with a few scattered defect clusters present in Figure 5a, while a polycrystalline structure formed in Figure 5b. In terms of deformation easiness, the TSV shown in Figure 5b contains more grains and GBs, and thus, its deformation is less easy. Figure 9 shows the initial microstruture in a larger size TSV together with its protrusion profile after applying the external loading. It can be observed that dislocations are always hindered by GBs, which means increasing the number of grains leads to less deformation. Meanwhile, grain interiors are usually more susceptible to deformation than GBs, so once GB motion is switched on, protrusion in large grains is greater than in small grains.  Figure 5a,b, respectively, subjected to the same external loading with T = 4.8 × 10 7 and β = 4π/12.
In general, the migration of GBs during polycrystalline deformation is controlled by the gradients of thermodynamic driving forces [37]. For the case studied here, the same loading is applied to the TSV, and the gradient of driving force can be provided by the difference in Young's modulus across two neighboring grains. To demonstrate this effect from the gradient of Young's modulus, a simple grain structure as shown in Figure 10a is considered here. For the convenience of setting Young's modulus, grains with an initial rectangular geometry are evolved from t = 0 to t = 10 according to Equation (5) to obtain a regular grain structure free of dislocations in the TSV, such that the protrusion behavior is simply governed by the gradient of Young's modulus. Five scenarios, denoted as case n with n ranging from 1 to 5, are investigated, with the Young's modulus of grain A in the top layer for case n assigned to 1.2(n + 2) × 10 13 , while those of grains B, C, and D  Figure 5a,b, respectively, subjected to the same external loading with T = 4.8 × 10 7 and β = 4π/12.
In general, the migration of GBs during polycrystalline deformation is controlled by the gradients of thermodynamic driving forces [37]. For the case studied here, the same loading is applied to the TSV, and the gradient of driving force can be provided by the difference in Young's modulus across two neighboring grains. To demonstrate this effect from the gradient of Young's modulus, a simple grain structure as shown in Figure 10a is considered here. For the convenience of setting Young's modulus, grains with an initial rectangular geometry are evolved from t = 0 to t = 10 according to Equation (5) to obtain a regular grain structure free of dislocations in the TSV, such that the protrusion behavior is simply governed by the gradient of Young's modulus. Five scenarios, denoted as case n with n ranging from 1 to 5, are investigated, with the Young's modulus of grain A in the top layer for case n assigned to 1.2(n + 2) × 10 13 , while those of grains B, C, and D are set to 1.2, 3.0, and 1.8 × 10 13 , respectively. Although sharing the same grain structure, different protrusion profiles result for the five cases studied here. To reveal further details on this difference, protrusion morphologies and displacement fields are plotted together and shown in Figure 10b-d for cases 1, 3, and 5, respectively. The results clearly suggest that the protrusion is greater when the mismatch of the Young's modulus between the two grains in the top layer is larger. The displacement field of atoms in the TSV can be useful to understand this result. When the Young's modulus mismatch is 2.4 × 10 13 and 3.6 × 10 13 , atoms are prone to rotation in the middle of TSV, and only a small proportion of the atoms in the top join the movement along the vertical direction, which makes the protrusion smaller. As the Young's modulus mismatch turns larger to 4.8 × 10 13 , the rotation in the middle of the TSV come to a stagnation, while a small proportion of the atoms in the top protrude or intrude more compared with Figure 10b. When the Young's modulus of grain A is up to four and five times as great as that of grain B, more atoms in the top move upward and thus lead to more protrusion. It can be explained that with a larger Young's modulus, grain A can resist deformation better than other grains in the TSV, thus first preventing the movement of atoms, and second causing stress concentration in the adjoining grains with smaller Young's moduli, i.e., grain B. In addition, the induced stress concentration can be considered as the driving force for boundary migration. Therefore, the driving force becomes larger with the mismatch between grains A and B turning larger and results in larger protrusion. Regardless of the number of grains in the top region, the stress concentration induced by a larger mismatch of Young's modulus can be regraded as a larger driving force, impelling a higher protrusion.
Not only the mismatch of Young's modulus amongst grains but also the number of grain layers could be correlated to the protrusion profile of TSV. With the same magnitude of external loading, i.e., 4.8 × 10 7 in this case, more protrusion occurs in TSVs with more grain layers. An increase in the number of grain layers, assuming a constant number of grains present in each layer, means that the number of grains inside the TSV becomes larger and the specific surface area of GBs increases, thus making plastic deformation under the same external loading more difficult. The same principle can be applied to study the effect of the number of grains in the top layer on protrusion. When the number of grains in the top layer increases, greater driving force is required to initiate the migration of GBs.
(a) (b) Figure 9. The microstructure (a) and the protrusion profile together with the displacement field (b) of a lager size TSV subjected to an external loading with t = 8 × 10 7 and β = 4π/12.
(a) (b) Figure 9. The microstructure (a) and the protrusion profile together with the displacement field (b) of a lager size TSV subjected to an external loading with t = 8 × 10 7 and β = 4π/12.
(c) (d) Figure 10. The microstructure (a) and protrusion profiles together with displacement fields of TSVs with different settings on the mismatch of Young's modulus for (b) case 1; (c) case 3; and (d) case 5.

Effect of External Loading
The protrusion profiles of TSVs are governed by the atomic-scale microstructure as well as the stress states caused by external loadings. Therefore, in this section, the effects of stress states on TSV protrusion are investigated. Selected protrusion profiles together with displacement fields for the TSV with microstructure shown in Figure 5b under external loadings with the same magnitude 6.4 × 10 7 but along different directions with β systematically varying from 0 to 23π/12 with an interval of π/12 are shown in Figure  11. The relationship between the protrusion height and the loading direction β is plotted in Figure 12, in which four characteristic regions can be clearly observed and referred to as region I, II, III, and IV, respectively. Despite a symmetric nature of the applied loadings, asymmetric displacement and strain fields result because of the presence of a polycrystalline structure in the TSV. Loadings along different directions result in different distributions of strain concentration inside the TSV, which in turn control the movement of the atoms. It is found that the normal strain ε yy generated at the top of the TSV is comparable in four regions as shown in Figure 13, and such that the shear strain needs to be carefully examined. In region I, the positive (counterclockwise) and negative (clockwise) shear strains disperse as fine spots due to the underlying atomic-scale structure in the top region of the TSV, and positive and negative shear strains are counter-balanced so that the shear strains exert little influence on the movement of the atoms in this region. Therefore, the net effect of an external loading in region I is only from ε yy , which drives the atoms in the top region to move upwards and results in protrusion. In regions II and IV, the positive and negative shear strains are distributed on the right and left parts of the TSV, causing the atoms in the top region to rotate, thus leading to less protrusion. In region III, the pattern of shear strain is different from the other regions. It is also distributed at the top as in region I, but the shear strain is almost counterclockwise, which drives the top of the TSV to also undergo a larger protrusion but with a different protrusion profile compared with region I. In general, a relatively smaller protrusion occurs in regions II and IV, but a larger protrusion occurs in regions I and III. Although a higher protrusion can occur in both regions I and III, different characteristics of the protrusion profiles can be observed, i.e., protrusion covering the entire top end in region I, while both protrusion and intrusion appear in region III.

Effect of External Loading
The protrusion profiles of TSVs are governed by the atomic-scale microstructure as well as the stress states caused by external loadings. Therefore, in this section, the effects of stress states on TSV protrusion are investigated. Selected protrusion profiles together with displacement fields for the TSV with microstructure shown in Figure 5b under external loadings with the same magnitude 6.4 × 10 7 but along different directions with β systematically varying from 0 to 23π/12 with an interval of π/12 are shown in Figure 11. The relationship between the protrusion height and the loading direction β is plotted in Figure 12, in which four characteristic regions can be clearly observed and referred to as region I, II, III, and IV, respectively. Despite a symmetric nature of the applied loadings, asymmetric displacement and strain fields result because of the presence of a polycrystalline structure in the TSV. Loadings along different directions result in different distributions of strain concentration inside the TSV, which in turn control the movement of the atoms. It is found that the normal strain ε yy generated at the top of the TSV is comparable in four regions as shown in Figure 13, and such that the shear strain needs to be carefully examined. In region I, the positive (counterclockwise) and negative (clockwise) shear strains disperse as fine spots due to the underlying atomic-scale structure in the top region of the TSV, and positive and negative shear strains are counter-balanced so that the shear strains exert little influence on the movement of the atoms in this region. Therefore, the net effect of an external loading in region I is only from ε yy , which drives the atoms in the top region to move upwards and results in protrusion. In regions II and IV, the positive and negative shear strains are distributed on the right and left parts of the TSV, causing the atoms in the top region to rotate, thus leading to less protrusion. In region III, the pattern of shear strain is different from the other regions. It is also distributed at the top as in region I, but the shear strain is almost counterclockwise, which drives the top of the TSV to also undergo a larger protrusion but with a different protrusion profile compared with region I. In general, a relatively smaller protrusion occurs in regions II and IV, but a larger protrusion occurs in regions I and III. Although a higher protrusion can occur in both regions I and III, different characteristics of the protrusion profiles can be observed, i.e., protrusion covering the entire top end in region I, while both protrusion and intrusion appear in region III.     Figure 13. The strain distribution in the TSV with microstructure shown in Figure 5b and β varying from 3π/12 to 21π/12. (a) +ε yy , β = 3π/12; (b) −ε yy , β = 3π/12; (c) +ε xy , β = 3π/12; (d) −ε xy , β = 3π/12; (e) +ε yy , β = 9π/12; (f) −ε yy , β = 9π/12; (g) +ε xy , β = 9π/12; (h) −ε xy , β = 9π/12; (i) +ε yy , β = 15π/12; (j) −ε yy , β = 15π/12; (k) +ε xy , β = 15π/12; (l) −ε xy , β = 15π/12; (m) +ε yy , β = 21π/12; (n) −ε yy , β = 21π/12; (o) +ε xy , β = 21π/12; (p) −ε xy , β = 21π/12.
With β = 2π/12, referred to as Model 1, and β = 6π/12, referred to as Model 2, a systematic change in the magnitude of the external loading applied to the TSV with the microstructure shown in Figure 5b can be studied by conducting a parametric sweep in COMSOL. Thus, a parameter p is introduced and set to change from 0 to 2 with an interval of 0.2, and T can be expressed as 8 × 10 7 × p. The protrusion height versus p for Models 1 and 2 is plotted in Figure 14, and the protrusion profiles caused by different values of p for model 2 are shown in Figure 15. It can be seen that maximum and minimum protrusion heights increase almost linearly with p, and the slopes of the curves vary in the two models. Moreover, it is found that the protrusion height becomes more and intrusion depth becomes less as the magnitude of the external loading increases but the original deformation trend, i.e., the protrusion and intrusion profile, remains similar. The external loading applied to TSV can be regarded as a driving force for deformation. A greater loading provides a larger driving force and results in more plastic deformation. interval of 0.2, and T can be expressed as 8 × 10 7 × p. The protrusion height versus p for Models 1 and 2 is plotted in Figure 14, and the protrusion profiles caused by different values of p for model 2 are shown in Figure 15. It can be seen that maximum and minimum protrusion heights increase almost linearly with p, and the slopes of the curves vary in the two models. Moreover, it is found that the protrusion height becomes more and intrusion depth becomes less as the magnitude of the external loading increases but the original deformation trend, i.e., the protrusion and intrusion profile, remains similar. The external loading applied to TSV can be regarded as a driving force for deformation. A greater loading provides a larger driving force and results in more plastic deformation.

Conclusions
Thermal management of 3DIC and related reliability concerns on TSV protrusion present serious challenges today. Although materials issues have been thoroughly investigated for the advancement of IC technologies, co-design of materials, processes, IC, and packaging is critical to realize the full potential of 3DICs. In this study, an APFC model coupled with mechanics is successfully applied to study the formation of grain structures in TSVs and the subsequent deformation behavior under external loadings from the atomic scale to a few tens of nanometers. TSV protrusion and intrusion have been observed, and the following conclusions can be drawn : (1) Providing an atomic-scale resolution, the APFC model implemented in this study has successfully extended the size and aspect ratio of simulated TSVs to around 30 nm and 4, respectively. Further increase of the size and aspect ratio to cover the micrometers range is possible with high-performance computing; (2) The mismatch of Young's modulus between grains in TSVs and external loadings can both be regarded as thermodynamic driving forces that govern the TSV protrusion and intrusion. A more homogenous grain structure can reduce the part of thermodynamic driving force caused by the mismatch of Young's modulus.
(3) The stress fields σ ψ , before loading, and σ δ , after loading, along with the displacement and strain fields present in TSVs can be obtained and used to understand the complex behavior of TSV protrusion and intrusion; (4) The protrusion height is found to vary nonlinearly with the direction of the external loading but vary almost linearly with the magnitude of the external loading.

Conclusions
Thermal management of 3DIC and related reliability concerns on TSV protrusion present serious challenges today. Although materials issues have been thoroughly investigated for the advancement of IC technologies, co-design of materials, processes, IC, and packaging is critical to realize the full potential of 3DICs. In this study, an APFC model coupled with mechanics is successfully applied to study the formation of grain structures in TSVs and the subsequent deformation behavior under external loadings from the atomic scale to a few tens of nanometers. TSV protrusion and intrusion have been observed, and the following conclusions can be drawn : (1) Providing an atomic-scale resolution, the APFC model implemented in this study has successfully extended the size and aspect ratio of simulated TSVs to around 30 nm and 4, respectively. Further increase of the size and aspect ratio to cover the micrometers range is possible with high-performance computing; (2) The mismatch of Young's modulus between grains in TSVs and external loadings can both be regarded as thermodynamic driving forces that govern the TSV protrusion and intrusion. A more homogenous grain structure can reduce the part of thermodynamic driving force caused by the mismatch of Young's modulus.
(3) The stress fields σ ψ , before loading, and σ δ , after loading, along with the displacement and strain fields present in TSVs can be obtained and used to understand the complex behavior of TSV protrusion and intrusion; (4) The protrusion height is found to vary nonlinearly with the direction of the external loading but vary almost linearly with the magnitude of the external loading.

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