The role of microstructure in the thermal fatigue of solder joints

Thermal fatigue is a common failure mode in electronic solder joints, yet the role of microstructure is incompletely understood. Here, we quantify the evolution of microstructure and damage in Sn-3Ag-0.5Cu joints throughout a ball grid array (BGA) package using EBSD mapping of localised subgrains, recrystallisation and heavily coarsened Ag3Sn. We then interpret the results with a multi-scale modelling approach that links from a continuum model at the package/board scale through to a crystal plasticity finite element model at the microstructure scale. We measure and explain the dependence of damage evolution on (i) the β-Sn crystal orientation(s) in single and multigrain joints, and (ii) the coefficient of thermal expansion (CTE) mismatch between tin grains in cyclic twinned multigrain joints. We further explore the relative importance of the solder microstructure versus the joint location in the array. The results provide a basis for designing optimum solder joint microstructures for thermal fatigue resistance.


INTRODUCTION
There is an ongoing need to improve the understanding and reliability of Pb-free solders, particularly for applications with high consequences of failure in the medical, automotive, aerospace and defence industries.Thermal fatigue resistance is an important design requirement that can be approached through electronic package design optimisation and by the design of improved materials, including new solder alloys and the optimisation of solder microstructures.
A range of past studies on accelerated thermal cycling of Sn-Ag-Cu solder joints 1-4 have shown that damage accumulation is sensitive to the solder microstructure including tin crystal orientation(s) and the intermetallic length scale, to the extent that failure can occur in a joint that is far from the location of maximum global shear strain 1 .A key factor is that tetragonal b-Sn is highly anisotropic, for example having a coefficient of thermal expansion (CTE) and Young's modulus that are approximately two and three times larger respectively along the c crystallographic direction than along the a=b directions 5,6 .
In the simplest case of single grain Sn-Ag-Cu solder joints, both experiments 1-3 and modelling 3,[7][8][9][10] have shown that the most severe thermal fatigue damage develops in single grain joints when the b-Sn c-axis lies in the plane of the printed circuit board (PCB).This has been attributed to two factors 1,2 : (1) this orientation maximises the in-plane CTE mismatch between the solder and package and, therefore, maximises in-plane shear during temperature cycling and (2) this orientation also minimises the b-Sn CTE component normal to the PCB which induces a tensile component of stress during heating.
The case of multigrain BGA joints is more complex and the optimum microstructure for thermal fatigue resistance is unclear from past work.Bieler et al. 1 found that some multigrain joints in a BGA package were not cracked when subjected to thermal cycling while, at the same time, unfavourably oriented single grain joints were cracked, and discussed this in terms of the majority orientation of the β-Sn grains.Arfaei et al. 11 reported longer thermal fatigue life for joints containing b-Sn grains with interlaced morphology compared with beachball morphology.Various authors [12][13][14][15] have noted the importance of CTE mismatch at grain boundaries in damage development in multigrain joints, but these have not quantitatively linked the different orientations within joints to the local CTE mismatch and subsequent damage.At the same time, microstructurally motivated models have been developed [16][17][18][19] to investigate microstructural effects on failure mechanisms and damage evolution in joints subject to thermal cycling and/or shear stress.However, these have not directly compared experimental measurements of microstructure and damage in BGA joints with mechanistic three-dimensional (3D) models of the joints.
Therefore, there is a need for coupled experimental and modelling studies on thermal cycling of BGA packages including complex multigrain solder microstructures to build mechanistic understanding on the role of microstructure in thermal fatigue damage.This will provide a basis for designing optimum solder microstructures for long-term thermal reliability in microelectronics that can be implemented in parallel with solder alloy design approaches (e.g. 20).
In this work, we build mechanistic understanding of how the Sn-3Ag-0.5Cu(wt%) solder microstructure affects thermal fatigue resistance in a ball grid array (BGA) package.This paper contains three areas of novelty in this direction that lead to insights into the role of microstructure in the thermal fatigue of BGA solder joints.(1) We develop an approach to processing electron backscatter diffraction (EBSD) maps that allows damage to be quantitatively compared across multiple complex Sn-Ag-Cu microstructure types, providing a versatile damage measure using the area fraction of the cross-section with misorientation larger than 8°.This approach is then applied to measure the damage accumulated in 20 time-zero joints and 84 thermally cycled joints from an 84CTBGA (thin chip array BGA) electronic test vehicle; (2) We develop a method to estimate the initial 3D microstructure of a BGA joint from a 2D EBSD map after mild thermal cycling; and (3) We perform a direct comparison between the experimentally-measured microstructure and damage in selected BGA joints and a 3D multi-scale thermal cycling model of the same joints, coupling a continuum finite element (FE) model at the scale of the whole 84CTBGA package and board, with a crystal plasticity finite element (CPFE) model of individual solder joints.The approach is overviewed in Figure 1.
Detailed methods are given at the end of the manuscript.

Microstructure variability
In both time-zero and thermally cycled joints, the β-Sn microstructures were single grain or multigrain with partial interlacing.Some multigrain joints contained three main orientations related by cyclic twinning as has been reported widely for Sn-3Ag-0.5Cujoints [21][22][23][24][25][26] .Other multigrain joints contained five main orientations related by double rings of cyclic twinning, similar to references [27][28][29][30] .The β-Sn morphology types are compared in time-zero versus thermally cycled joints in Figure 2. Figure 2 Supplementary Table 1 summarises the frequency of occurrence of β-Sn microstructure types for 20 time-zero joints and 84 thermally cycled joints, excluding recrystallised regions at the top and bottom.There is reasonable consistency in the percentage of each microstructure type in the time-zero and thermally cycled packages, although a relatively low number of timezero joints were analysed (20).Supplementary Table 1 also shows that >96% of multigrain joints contained grains all interrelated through a twin relationship and, therefore, originated from a single b-Sn nucleation event 21 .
From Supplementary Table 1 and Figure 1-3, we see highly variable β-Sn microstructures from joint to joint.The origin of this variability is: (i) the stochastic undercooling for β-Sn nucleation [21][22][23]30,31 , (ii) the variable β-Sn nucleation location, and (iii) the small number of grains developed in the small solder volume 1, 16,24 . The efect of different Sn nucleation locations can be seen clearly by comparing the joints in Figure 1(d) and Figure 3(a), where it can be inferred that nucleation occurred from the package (top) side in Figure 3(a,d) and from the PCB (bottom) side in Figure 1(d).
In many cases, the 3D grain structure could be reasonably estimated from 2D EBSD maps.
Our approach utilises the crystallographic geometry of the twin grain boundaries and is overviewed in Figure 3.In Figures 1-3, the multigrain joints contain some straight grain boundary segments that are clearest away from regions of interlacing.Prior work by Lehman et al. 21reported that such grain boundaries are {101} cyclic twin planes and result in "beachball" microstructures when the fraction of interlacing is low.Figure 3 (a) and (b) shows this using one joint from this study and Figure 3 (c) confirms {101} macroscopic beachball boundaries from 29 joints reanalysed from EBSD datasets in 29,32 .Note in Figure 3 (a) that the boundaries are not straight {101} planes at the micrometre scale but can be well-approximated as {101} boundaries at the scale of the joint.Combining the measured cyclic twinned orientations, the measured 2D vectors formed where the grain boundaries intersect the polishing plane, and the knowledge that the grain boundaries are {101} planes, enables the time-zero macroscopic 3D grain structure to be generated from measured 2D EBSD maps.This is demonstrated in Figure 3(d) and Figure 1(d) gives another example containing partial interlacing.Further details of the approach are given in the "Building 3D microstructure models from 2D EBSD maps" part of the Methods section.
Compared with the tin grain structure, there was much less variability in eutectic intermetallic size among joints before thermal cycling and among those after thermal cycling.This is overviewed in Supplementary Figure 1.There was also strong localised accelerated Ag3Sn coarsening near the top interface, as is common in thermally cycled Sn-3Ag-0.5Cujoints 12,25,33 due to coarsening with both a thermal [34][35][36] and strain-enhanced 12,37,38 component and further due to faster diffusion along recrystallised grain boundaries once recrystallisation has begun 12 .
Here, these accelerated coarsened Ag3Sn particles exceeded 2µm (Supplementary Figure 1(f)) and are considered further in the next section.

Microstructure evolution and damage metrics
The features of microstructure evolution in thermally cycled joints are overviewed in Figure 4 using one joint as an example.Equiaxed recrystallised grains (with MO>15°) have developed from the top-right corner and span ~80% of the way along the top region (Figure 4 (b)-(d)).
To the left of and below the recrystallised grains are subgrains with MO<15° (Figure 4 (c)-(d)) which appear as different shades of light blue in Figure 4 (c).There is a strong correlation between recrystallised grains and accelerated coarsened Ag3Sn particles.Comparing Figure 4 (a), (d) and (e), note that the large (~2µm) accelerated coarsened Ag3Sn particles are often located at recrystallised grain boundaries and the interior of recrystallised grains are, in many places, particle free zones.In contrast, subgrains (e.g.blue regions in (c) and regions within blue boundaries in (d)) contain numerous Ag3Sn particles that are only slightly coarser than those at the centre of the joint (e.g.see the top left region in Figure 4(d) and (e)).Comparing (f) and (d), the regions with anomalously wide Ag3Sn particle spacing (yellow triangles in (f)) match well with the regions of recrystallised grains in (d).A crack has propagated from the top-right along some of the grain boundaries between the recrystallised grains and the parent grain, Figure 4 (a).
The features of microstructure evolution in Figure 4 are similar to a large body of past work 12,22,[39][40][41][42] including by Bieler at al. 1,41 and Arfaei et al. 11,12,22 .However, the extent and detail of these features after 7580 cycles in the 84CTBGA varied strongly from joint to joint depending on the initial microstructure as well as the location in the array.To enable a quantitative comparison across all studied joints, we define two measures of damage as overviewed in the example in Figure 5 (a)-(d).First, we define the "cracking extent" as the ratio of the projected crack length to the top interface length in 2D sections (Figure 5 (a)).
Second, as a measure of microstructural damage, we define the area fraction of the crosssection with MO higher than 8°.The rationale for selecting MO>8° can be seen by comparing Figure 4 (h) with the rest of Figure 4.In Figure 4 (h), red has been assigned to regions with MO>15°, yellow to regions with 8≤MO≤15°, and blue to regions with MO<8°.The region with MO>15° correlates well with the regions of recrystallisation and accelerated Ag3Sn particle coarsening whereas the region with 8≤MO≤15° captures subgrains and regions of strong localised lattice curvature.The selection of a minimum value of 8° was to exclude small misorientations that are present even in time-zero joints (Figure 2).The MO>8° area is preferred here because lattice curvature and subgrain development are considered important aspects of damage accumulation in addition to recrystallisation.However, the conclusions of this paper are not significantly altered if MO>15° is used.When calculating the area fraction with MO>8°, we consider only the top quarter of the joint (see Figure 5 (d)) since this is where cracking and failure occurred.However, again, the results of this paper are not significantly altered if the whole joint area is used rather than the top quarter only.
Figure 5 (e) plots the extent of cracking vs. the area fraction with MO>8° for the 80 joints investigated in one thermally cycled 84CTBGA package.There is a high degree of scatter and a poor correlation.We note that measuring the extent of partial cracking in 2D (X-Y) crosssections has limitations.In prior work 43 , top-down (X-Z) sectioning of BGA packages has shown the complex shape of cracks as well as crack growth from multiple locations around the top circumference of the joint.It was demonstrated that the appearance of a partial crack in a 2D (X-Y) section would depend strongly on the location of the cross-section.
To proceed, we averaged the data for individual joints at equivalent locations in the package exploiting the 1/8 symmetry of the 84CTBGA as indicated in the inset of Figure 5 (f).
The correlation between the average extent of cracking and the area fraction with MO>8° is much stronger (Figure 5 (f)).We next selected individual joints for which the neutral point vector lies close to the sectioning plane as indicated in the inset of Figure 5 (g).For these individual joints (Figure 5 (g)), the general correlation is much clearer than when considering all individual joints in Figure 5 (e), albeit with some remaining scatter.Comparing Figure 5 (e) versus (g), the results are consistent with the extent of cracking depending on the crosssectioning plane and indicate that the greatest error (in crack extent) occurs for joints where the neutral point direction is near-perpendicular to the cross-section (e.g.M6,7 and A6,7).It is likely that there are true microstructure effects that determine how the extent of cracking depends on the evolving microstructure and therefore MO; however, this could not be studied here due to the limitations of measuring partial cracks in 2D sections.The MO>8° approach taken here is also a 2D method but is less prone to cross-sectioning effects since the MO features occupy a larger volume of the joint and EBSD mapping captures more of the 3D information than the (essentially 1D) extent of cracking.Therefore, for the remainder of this paper we consider MO>8° as the measure of damage.

Single grain joints
For all 15 single grain joints studied after thermal cycling, Figure 6 (a) plots the area fraction with MO>8° versus the angle between the ⟨001⟩Sn and the Y-axis.There is a clear increase in damage as the angle increases, which is consistent with prior experimental studies that used other measures of damage 1-3 .The result in Figure 6 (a) can be partially understood from Figure

Multigrain joints
A key finding in this work is that many multigrain joints were more damaged than the worstoriented single grain joints, both in terms of MO>8° and in terms of cracking extent.EBSD analysis of five thermally cycled solder joints is shown in Figure 7  For single grain joints, the extent of damage generally increases with increasing in-plane ΔCTE to the package which is the same result as in Figure 6 (a), although plotted as a CTE mismatch.For multigrain joints with a single orientation at the top in Figure 7 (d), the MO>8° data are slightly higher (although with a high degree of overlap) than the single grain data at the same CTE mismatch values.
For joints with two or more β-Sn orientations at the top, the MO>8° data are significantly higher than for single grain joints in many cases.

Multi-scale modelling on multi grain joints
To understand the effects of microstructure while accounting for the boundary conditions at different locations in the array, the multiscale model was applied to the joints at locations M09, C03, M07 and C07 as overviewed in could be inferred from the presence of some small cyclic twinned domains near the edge of the recrystallized grains in the EBSD maps although the extent of the initial interlacing could not be determined accurately.
Geometrically necessary dislocation (GND) distributions were calculated and used to derive simulated MO distributions (based on the methodology in ref 3,44 ) at the end of 10 thermal cycles in Figure 9 (e).Experimental MO distributions are given in Figure 9 (f).A different colour range is used since the experimental MO is after 7580 thermal cycles when significant recrystallisation has occurred whereas the simulated MO is after 10 thermal cycles.
We see that regions with high simulated MO in Figure 9 (e) correspond reasonably well with experimental regions of subgrains and recrystallisation in Figure 9 (f).For example, in single grain joint M09, both the simulated MO hotspots and the measured regions of highest MO occur in a band adjacent to the top interface.In the multi-grain joints, the simulations contain localised MO within interlaced regions which is also apparent in the experiments.This is particularly clear in C03 where the initial interlacing was on the bottom side and the simulated MO in this region reasonably matches the experimentally measured extent of high MO spreading from the bottom towards the centre of the joint.In the three multi-grain joints, there is a simulated MO hotspot at one or both top corners and the experimentally high MO spreads laterally from these regions adjacent to the top interface.In C07, the simulated MO contains an MO hotspot towards the bottom right of the joint and the experimental MO includes a region of high MO on the bottom side, slightly to the left of the simulated hotspot.
The predicted SED developed within these solder joints is shown in Figure 9 (d) at the end of 10 thermal cycles.Compared with the single grain joint M09, the SED in the multigrain joints is more strongly localised at hotspots (peaks) along the top circumference in C03, M07 and C07, at sites associated with the intersections between βSn-βSn grain boundaries where stress concentrates 45 .To investigate the peaks in SED, the perimeter distribution of the SED along the path A-A', indicated in Figure 9 (d), is compared among the four joints in Figure 10 (a).Note that the maximum on the SED colour scale in Figure 9 (d) is G = 2 J m -2 whereas the peaks in SED in Figure 10 (a) are significantly higher.Thus, the hotspots in Figure 9 (d) are useful for confirming that peaks are located at the intersections of grain boundaries but cannot be used to compare the magnitude of the peaks.In Figure 10 (a), we see that the peak SED developed in joint C07 is considerably larger than peaks in M07 and C03, consistent with the experimental damage observations in Figure 9 (a) and (f).Also notice in Figure 10 (a) that the worst-oriented single grain joint (M09, see Figure 6) develops relatively low peaks in SED compared with the multi-grain solder joints, consistent with the experimentally measured damage in Figure 9 and Figure 7(d)-(e).

DISCUSSION
While the CPFE approach does not explicitly capture subgrain formation, recrystallisation and crack initiation, all these processes are directly driven by peaks in stored energy density (SED).
For example, subgrains arise from thermally activated dislocation rearrangement and first develop in regions of locally high SED.During continued cycling, the same regions produce more dislocations that migrate to nearby subgrain boundaries, increase the subgrain boundary angle, and ultimately lead to continuous recrystallisation 41 .Similarly, a locally high SED is the driving force required to nucleate discontinuous recrystallisation 46,47 as well as to initiate a fatigue crack.Indeed, a similar SED criterion has been shown to predict fatigue crack nucleation sites in other ductile alloys including steels 48 , nickel 49,50 and zirconium alloys 51 .
Quantitative comparison between the peak SED predicted with the CPFE model and the experimentally measured damage is shown in Figure 10  (area fraction with MO>8°) and the simulated peaks in SED, consistent with local peaks in SED being the driving force for recrystallisation 47,52,53 and fatigue crack nucleation [48][49][50] .
Comparing the simulations and experiments in Figures 6, 9 and 10 highlights the key role of the magnitude of local peaks in SED in determining the extent of damage accumulation (area fraction of subgrains and recrystallisation in the region where cracking occurs) in this study.
In Figures 9 and 10, we see that multi-grain joints develop higher peaks in SED than single grain joints due to stress concentration near grain boundaries where there is mismatch in CTE and elasticity due to the anisotropy of b-Sn.The negative effects of grain boundaries are most significant when the boundaries intersect the top interface where they accumulate the highest SED and when there is a large intrinsic CTE mismatch between these grains.Experimentally, in single grain joints, the MO>8° increased as the angle between ⟨001⟩Sn and the Y-axis increased.Increasing this angle increases the in-plane mismatch in coefficient of thermal expansion (CTE) between the tin grain and the package (Figure 7d).Many multigrain joints had significantly higher damage (MO>8°) than the worst-oriented single grain joints because there is additional CTE mismatch between grains which causes SED localisation near grain boundaries.The extent of damage increased as the CTE mismatch between grains in the substate plane increased (Figure 7e).
The findings of this study reinforce the importance of the solder microstructure in determining the development of damage in solder joints.Our findings suggest that the optimum microstructure for thermal fatigue resistance is a single Sn grain with ⟨001⟩Sn perpendicular to the substrate.In multi-grain joints, avoiding grain boundaries intersecting the top interface and minimising the CTE mismatch between tin grains may also be useful.For our findings to be usable by industry, it will be necessary to develop industrial scale approaches to control the microstructure formed during the assembly/soldering process and combine this with alloy development.Some methods have been published to control tin orientations 24,54,55 or promote the formation of single grain joints 24 in Sn-Ag-Cu alloys, and certain alloying additions are also known to promote single grain joints 56 .It will be interesting in future work to optimise these microstructure control methods and measure the extent to which these approaches can improve the thermal fatigue performance of BGA and other packages.

Electronic test vehicle and thermal cycling
The test vehicle involved daisy-chained 84 input/output (IO) thin core chip array (84CTBGA) packages with an electrolytic Ni/Au pad finish.The printed circuit board (PCB) had high temperature Panasonic R-1755V laminate and a surface finish of organic solderability preservative (OSP), with 16 sites for the 84CTBGA.The 84CTBGA packages had 300 µm Sn-3Ag-0.5Cu(wt%) balls attached and were assembled on the PCB using Sn-3Ag-0.5Cu(wt%) type 4 no-clean solder paste applied with a 130 µm thick stencil with 305 x 305 µm square apertures.The test vehicles were then reflowed in a 14-temperature zone convection oven in a N2 atmosphere with a nominal peak temperature of 240 ºC measured on the board adjacent to the solder joints.The test vehicle, assembly and soldering procedures are the same as ref. 20.
A photograph of an assembled package is given in Figure 1(a).A pin diagram of the 84CTBGA package is shown in Figure 2.
Two PCBs with a total of 32 84CTBGA packages were thermally cycled in accordance with the IPC-9701A guidelines 57 .The thermal profile involved ramping from 0/100°C at 10°C per minute with 10-minute dwells at the hot and cold extremes of the cycle.The solder joints were monitored in situ using an event detector set at a resistance limit of 1000 ohms.The characteristic life (N63.2) from a 2-parameter Weibull analysis of 32 packages in this thermal cycling profile was 7286 cycles.The strategy in this work was to apply a fixed number of thermal cycles to all joints in selected BGA packages and then to study how the accumulated damage depends on the solder microstructure and the location in the array.The thermally cycled packages studied here had undergone 7580 cycles.

Microstructure characterisation
20 joints were studied in the as-soldered (time-zero) condition and 84 other joints were studied after thermal cycling.One thermally cycled package was given a near-complete analysis in which 80 out of 84 joints were characterised to study the combined effects of microstructure and location in the array on damage evolution.This involved grinding into the first row of joints, characterising the microstructures of each joint, grinding into the next row of joints, and so on through the package.
Packages were mounted vertically on the side of row M in Struers VersoCit-2 acrylic coldmounting resin, ground down to 4000 grit SiC paper before being polished with OPS on a MD/DP-Nap pad.Backscattered electron (BSE) images were taken using a Zeiss Sigma300 field emission gun scanning electron microscope (FEG-SEM).To measure the eutectic Ag3Sn and Cu6Sn5, a fine imaging pixel size of 19.49 nm was used for time-zero samples and 38.98 nm for thermally cycled samples.To capture the entire joint at this resolution, the SEM stage was automated to scan across the joint, collecting 14×10 images for each time-zero joint and 7×5 images for each thermally cycled joint.The collected images were then stitched into a single image in ImageJ.Fine Ag3Sn and Cu6Sn5 particles were segmented using the in-house MATLAB image processing algorithms published in 36 .The eutectic Ag3Sn particle spacing was found from a Delaunay triangulation in MATLAB, using the centroids of all detected Ag3Sn particles as the vertices of the triangulation.Particle spacing maps were generated by filling the triangles with a colour scale that was proportional to the mean side length of each triangle.
Crystal orientation maps were collected by a Bruker e-Flash electron backscatter diffraction (EBSD) detector with the ESPRIT 2.1 software.All EBSD datasets were postanalysed using the MTEX toolbox in MATLAB.A data-processing pipeline was developed to characterise the misorientation maps for different types of β-Sn microstructures: single grain, cyclic twinned beachball and interlaced.For single grain joints, the reference orientation was taken as the orientation with the highest density in the MTEX orientation distribution function (ODF) near the centre of the joint and the misorientation of each pixel was calculated relative to this reference orientation.For joints with multiple β-Sn grains, the pipeline included: 1) select all major β-Sn orientations from the ODF near the joint centre as the reference orientations, 2) generate one misorientation map for each reference β-Sn orientation, 3) take the minimum misorientation value at each pixel and form a new map.The method developed here enabled misorientation maps to be quantified in complicated β-Sn microstructures and allowed numerous joints to be analysed within a single framework.
Building 3D microstructure models from 2D EBSD maps 3D beachball microstructures were generated from 2D EBSD maps by measuring the crystallographic orientations of the β-Sn grains (from the ODF near the joint centre) and the location of the grain boundary vectors in the 2D section and using the following 3D information/assumptions. (i) A beachball is a cyclic twin with six twin segments made up of three β-Sn orientations related through 60 ° rotations around a common ⟨010⟩ axis 21 .(ii) The macroscopic grain boundaries are {101} planes (Figure 3).(iii) The β-Sn nucleation point is on the Cu6Sn5 layer, either the top or the bottom layer 23 .(iv) Measured cross-sections such as Figure 3 are vertical (X-Y) sections through the centre of the joint.An example of this construction is given in Figure 3(d).For joints with interlacing, an interlaced volume was generated by intersecting the joint volume with a spherical interlaced volume whose radius and location were tuned to reproduce a digital cross section similar to the EBSD cross section.The interlaced volume was then filled by Voronoi tetrahedra with sizes mimicking those in the EBSD map, and with the measured orientations randomly assigned among the interlaced grains.

Macro-scale modelling at the package/board scale
The properties used for the individual materials in the PCB and package excluding the solder alloy are given in 3 (Table 1 in that paper).For the Sn-3.0Ag-0.5Cumaterial model, two different macro-scale continuum models were compared to investigate the effect of the model assumptions on the cyclic inelastic deformation (displacements) of solder joints.Approach 1 combined separate models for kinematic hardening plasticity (Chaboche formulation) and for creep (Garofalo-Arrhenius formulation).Approach 2 used the Anand model where the sum of plastic strain and creep strain are considered in a unified inelastic strain.
In approach 1, the non-linear kinematic hardening rule allows the modelling of cyclic hardening and can capture the Bauschinger effect.The nonlinear kinematic hardening model of Sn-3.0Ag-0.5Cudeployed in this work is a rate-independent version of the kinematic hardening model for the back-stress tensor  proposed by Chaboche 58,59 and obtained by superimposing three evolving kinematic back-stress tensors  !,  = 1,2,3.
The evolution of each back-stress model is defined with the kinematic hardening rule, where  ! and  ! are model material parameters, ̇% & is the plastic strain rate, and ̅ %& is the magnitude of the plastic strain rate.The Chaboche model parameters for Sn-3.0Ag-0.5Cuare detailed in Supplementary Table 2.
The creep behaviour part of approach 1 was modelled with the Garofalo-Arrhenius constitutive law, where ̇' ( is the steady-state creep strain rate,  is the applied stress,  is the temperature in Kelvin,  and B are material constants,  is the gas constant,  is an activation energy, and  is a stress exponent.The Garofalo-Arrhenius creep model reported in 60 was deployed with temperature-dependant values of the material constant  to capture the creep behaviour of Sn-3.0Ag-0.5Cu.The creep constitutive law parameters are detailed in Supplementary Table 3. For approach 2, the Anand model constitutive law was used for the inelastic strain rate of Sn-3.0Ag-0.5Cuand the respective model constants were the same as those reported in 3 (Table 2 in that paper).
The finite element analysis (FEA) revealed that the displacement predictions for the solder joints in approach 1 (separate kinematic hardening plasticity and pure creep constitutive laws) were very closely matched to the equivalent simulation run with approach 2 (the Anand model with unified inelastic strain rate).Supplementary Figure 2 and Supplementary Figure 3 show a comparative analysis of the displacement results for solder joint M09 from the two macroscale model approaches, showing that nodal displacements calculated from the two approaches are very similar, and the extreme point displacements in the solder joint differ by at most 3.7% and generally by <1%.Statistical analysis of the relative difference in predicted displacements for all solder joint mesh nodes, predicted with the two different material constitutive law approaches for Sn-3Ag-0.5Cu,show that the mean (average) value of the difference is only 0.8% and the respective standard deviation is 0.9%.This high similarity suggests that, in this study, kinematic hardening does not have a significant effect on the deformation of solder joints, and that displacement is dominated by creep.Several reasons for this are identified: (1) Sn-

Crystal plasticity finite element model at the microstructure scale
For the CPFE model, the 3D microstructure models (built from 2D EBSD maps) and the temperature-dependent, anisotropic material properties of Sn-3.0Ag-0.5Cu in Supplementary Table 4 61-64 were used.These material properties were extracted from a mechanics-based homogenization of β-Sn and embedded IMC phases, where the material properties were obtained through a combination of crystal plasticity modelling and micro-pillar experimental tests 61 .The homogenization also incorporates the IMC size effect into the hardening rate.
In the CPFE model, the slip rate ̇( !) along the i th slip system is determined by temperature-dependent thermal activation events where pinned dislocations jump over obstacles (such as lattice defects, inclusions etc.) in both forward and backward directions, which is given by: where ρssdm is the mobile statistically stored dislocation density, b the magnitude of the Burgers vector, vD the dislocation jumping frequency,  the activation energy, DV activation volume, k the Boltzmann constant, and  (!) and  ' (!) are the resolved shear stress and critical resolved shear stress on the i th slip system, respectively.
The geometrically necessary dislocation (GND) density term ρGND is established from the plastic strain gradients that accommodate lattice curvature, which are calculated from Nye's dislocation tensor Λ that originates from the curl of plastic deformation gradient F P : A microstructure-sensitive, Griffith-Stroh type energy-density driving force is employed to investigate fatigue damage that is given by: where G is the stored energy density (SED) and  is the fraction of the total plastic energy stored in the dislocation structure, taken as 0.05 65 . ++, is the statistically stored dislocation density and  2", is the geometrically necessary dislocation density.
The statistically stored dislocation density accumulation rate is related to the effective plastic strain rate ̇ and current accumulated plastic strain p according to: where  $ and  .control the slip system hardening and recovery rates respectively and depend on the average Ag3Sn and Cu6Sn5 IMC radius as given in Supplementary Table 4.
and the CTE mismatch analysis and contributed to the discussion; S.S. performed and discussed the macro-scale FE analysis; Y.L.X.performed and discussed the CPFE analysis; C.M.G. and F.P.E.D. developed and supervised the project and contributed to the analysis and discussion.All authors wrote parts of the paper.

Competing Interests Statement:
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
figure (IPF) maps with respect to the Y-direction) for time-zero joints and different thermally cycled joints respectively.The main β-Sn orientations are plotted as unit cell wireframes from the measured Euler angles and have been translated to highlight the cyclic twin relationships and the shared orientation in double rings.(b) and (d) are β-Sn misorientation (MO) maps, where the colour scale has been selected to clearly show small misorientations in the range 0-5° as changing shades of blue.Comparing the MO maps in Figure 2(b) and (d), note that the blue/white features are similar in the time-zero joints and the thermally cycled joints, both in terms of the microstructure revealed by the blue/white features and by the magnitude of the MO (the shade of blue).In contrast, there are regions with localised high MO near the top (the package side), and sometimes the bottom, in joints after thermal cycling (Figure 2(d)) which are absent in all time-zero joints (Figure 2(b)).These are regions of recrystallised grains, subgrains and local lattice curvature that are considered in more detail in the next section.From Figure 2(b) and (d), we see that strain and damage localise at the top and bottom of joints with less in central regions.Therefore, it is reasonable to use the major β-Sn orientation(s) from near the centre of a thermally cycled joint to infer the microstructure of that joint at time zero.

6
(b) which plots the mean CTE of b-Sn in the X-Z plane as a function of the angle between the ⟨001⟩Sn and the Y-axis.Also shown is the mean CTE of the package in the X-Z plane which was measured to be 11 × 10 -6 K -1 in20 .The mean in-plane CTE of b-Sn increases as the angle between the ⟨001⟩Sn and the Y-axis increases, and the CTE mismatch with the package increases in a similar manner.However, this is only one factor.A more complete analysis needs to consider the CTE mismatch in other directions, the anisotropic elasticity and plasticity as well as the different boundary conditions at each position in the array of joints.To do this, the multi-scale model was applied to the six single grain joints at locations A6, A7, A8, A9, M8 and M9 using the measured orientations for each joint.The distribution of stored energy density (SED) developed in three of the solder joints at the end of 10 thermal cycles is shown on the left of Figure6 (c).The localization of SED occurs at hotspots along the top interface and particularly along the top circumference.Experimental MO maps in the X-Y cross-section are shown on the right of Figure 6 (c), where regions of localised MO are in reasonable agreement with the simulated regions of high SED.The peak in SED was found by averaging the highest 10% of SED values (i.e. the maximum 20 elemental values) along the top interface of each joint and these peak SED values are plotted against the angle between the ⟨001⟩Sn and the Y-axis in Figure 6 (d).The experimental MO>8° of the same six joints is given as a second axis.The correlation in Figure 6 (d) shows that the peak SED calculated with the CPFE model after 10 thermal cycles is a good indicator of the area fraction of microstructure occupied by subgrains and recrystallised grains after 7580 thermal cycles in the experiments.
Figure 8 (a)-(b).In Figure 8 (a), the in-plane CTE of the b-Sn grain is the blue ellipse on a grey background and the dashed ellipse is the in-plane CTE of the 84CTBGA package from 20 .

Figure 8 (
Figure 8 (b) shows another (orange) single grain example with a larger substrate-plane CTE ellipse in the grey substrate plane because the ⟨001⟩Sn lies nearly parallel with the substrate.Now consider beachball microstructures consisting of a ring of three twinned β-Sn orientations such as that in Figure 8 (c) and (d). Figure 8 (e) plots the triple CTE ellipsoids corresponding to the three cyclic twinned grains oriented as in Figure 8 (c).Figure 8 (f) thenshows the digitally reorientated 3D triple ellipsoids with a translucent X-Z substrate plane, and its 2D section in the grey substrate plane plus the dashed ellipse of the in-plane CTE of the 84CTBGA package.This visualises and quantifies the in-plane CTE mismatch among the three β-Sn grains in the substrate plane and between these β-Sn grains and the package.The CTE mismatch among the three tin orientations in the substrate plane depends on the orientation of the tricrystal with respect to the substrate.The same is true of double ring cyclic twins where there are five interrelated β-Sn orientations, and a similar procedure of 2D sectioning the centre of five ellipsoids in the substrate plane allows the in-plane CTE to be calculated.This is overviewed in Figure8(g) which plots the calculated maximum difference in the mean in-plane CTE (the max ΔCTE) between β-Sn grains in the substrate plane (i.e. the difference in the mean CTE between the largest and smallest β-Sn ellipses on the grey substrate plane in Figure8 (f)) for different single ring and double ring cyclic twinned β-Sn orientations.Our calculations assume that all β-Sn orientations contact the top package.The plots were generated by applying a set of 50,000 random orientations to single ring and double ring cyclic twins, calculating the maximum ΔCTE in the substrate plane for each, plotting the orientations as data points in Euler space, and colouring each data point by the calculated maximum ΔCTE between β-Sn grains in the substrate plane.The 50,000 maximum ΔCTE values are plotted as histograms in Figure8(h) for both single and double ring microstructure types, revealing various features: (i) there is overlap between the distributions such that some single ring microstructures have a higher maximum ΔCTE value than some double ring microstructures

Figure 7 (
e) plots the extent of damage (area fraction with MO>8°) versus the maximum ΔCTE between β-Sn grains.The data points are coloured by the number of β-Sn orientations at the top, e.g. the dark red points are joints with five β-Sn orientations (double ring cyclic twins) at the top.The extent of damage generally increases with increasing maximum in-plane ΔCTE between the β-Sn grains.We also see that joints with more β-Sn orientations at the top were more likely to have a high maximum inplane ΔCTE, consistent with the probability distributions in Figure8 (g,h), and therefore accumulated more damage.Comparing Figure7 (d)and (e), note that the vertical axis is the same and spans from 0 (no MO>8° in the top quarter) to 1 (100% of the area in the top quarter having MO>8°).The most damaged joints are those with both multiple β-Sn grains at the top and a high CTE mismatch between those grains.The least damaged joints are those with a single β-Sn grain at the top and a small CTE mismatch to the package.Note that Figure7(d) and (e) do not account for any effects of the location in the array (e.g.distance to the neutral point) except for excluding the inner ring, and yet there is a reasonable trend between the CTE mismatch created by the b-Sn orientations and the area fraction with MO>8°.This highlights the strength of the b-Sn microstructure effect in determining damage accumulation in Sn-3Ag-0.5Cusolder joints in this package.At the same time, some location effects are evident.Outer corner joints are marked with a square around their datapoint and two of the three corner joints in Figure7(e) are somewhat anomalous.It is possible that their high degree of damage resulted in minor assolidified β-Sn grain orientations being missed and it is also likely that the different boundary conditions at these locations played a role.

Figure 9 .
These joints contain a single b-Sn grain in M09, cyclic twinned b-Sn with a single ring in C03 and M07, and a cyclic twinned double ring in C07 (Figure 9 (a)).C03 grew from the bottom (PCB) side with apparently one β-Sn orientation at the top in the 2D section whereas M07 and C07 grew from the top side with three orientations at the top.The generated digital microstructures are shown in Figure 9 (b) and (c), revealing that C03 has two β-Sn orientations at the top in 3D.The interlacing in C03, M07 and C07 (b) for the six single grain joints and three multigrain joints.The vertical axis is the average SED in the highest 10% values (i.e. the maximum 20 elemental values) along the top interface of each joint.Figure 10 (b) and Figure 6 (d) plot the same simulated SED data for the six single grain joints but with different SED axis limits.In Figure 10 (b), all joints show a strong correlation between the measured damage

Figure 1
Figure 1 overviews the parallel experiments and modelling approach.The multi-scale model couples (i) a macro-scale continuum model of the BGA package, solder joints and printed circuit board which considers the solder as an isotropic viscoplastic continuum with (ii) a dislocation based, temperature-dependant CPFE model of the individual solder joints, which

3 .
0Ag-0.5Cu has a high homologous temperature over the investigated mild temperature cycle 0°C to 100°C (T/Tm = 0.54 -0.74) which makes creep behaviour strongly dominant over kinematic hardening plasticity, (2) the long temperature cycle duration of 40 minutes, with 10minute dwells at both the low (0°C) and the high (100°C) temperature extreme, is associated with low strain rates that favour creep deformation and stress relaxation throughout the cycle, (3) solder joint movements at the component assembly level are kinematically constrained by the movement of the printed circuit board and the BGA package.While these conclusions have been drawn for the problem in this paper, solder joints subjected to temperature cycles with lower-temperature regimes and fast cycling rates could promote kinematic hardening behaviour over creep.Due to the similar results from the two macro-scale model approaches here, only the Anand model transient displacement boundary conditions were passed over to the microstructure level CPFE model.

Figure 2 .
Figure 2. Typical solder joint microstructures in the 84CTBGA before and after thermal cycling.a,b Time-zero.c,d after 7580 thermal cycles from 0/100°C.a,c EBSD orientation maps (IPF-Y) and Sn unit cell wireframes translated to highlight the cyclic twins.b,d EBSD misorientation maps with reference to the major Sn orientation(s).Pin diagrams indicate the locations of each joint in the array and labels of selected joints are marked under EBSD maps.

Figure 3 .
Figure 3. 3D microstructure models from 2D EBSD maps.a EBSD β-Sn orientation map with IPF-X colouring of a typical beachball microstructure with three twinned β-Sn grains whose 3D boundaries intersect the EBSD cross section forming near-linear 2D trace vectors A and B. b A stereogram plotting the 2D trace vectors A and B from (a) with respect to the β-Sn orientations on either side of Sn-Sn boundaries.The four datapoints (A1, A2, B2 and B3) are repeated 8 times according to the 4/mmm point group symmetry of β-Sn.c Summary of 82 IPF-trace vectors measured from 29 solder joints showing that macroscopic beachball grain boundaries are approximately {101}Sn planes.d A 3D digital model of the microstructure assuming that (a) was polished to the centre of the joint and that {101}Sn beachball boundaries emanate from a point where the common ⟨010⟩ twin axis touches the intermetallic reaction layer (the top side here).The 3D viewing direction in (d) is close to the ⟨010⟩ twin axis.

Figure 4 .
Figure 4. Typical features of microstructure and damage evolution due to thermal cycling using joint G08 as an example.a BSE-SEM image near the top of the joint, b EBSD orientation map with IPF-X colouring, c misorientation map with reference to the as-solidified major Sn orientations, d grain boundary segments categorized by misorientation: 5-15° = subgrains; >15° = recrystallised grain boundaries, e IMC layer and particle segmentation, f-g Ag3Sn particle spacing maps generated with a Delaunay triangulation of all measured Ag3Sn particle centroids.h Misorientation map using threshold ranges of 15>MO>8° and MO>15°.

Figure 5 .
Figure 5. Damage metrics.a Definition of 2D cracking extent measured using a BSE-SEM image.b-d Definition of microstructural damage extent: b) IPF-X map, c) misorientation (MO) map, d) filtered MO maps after removing MO>8° and MO>15°.e Plot of cracking extent vs. area fraction of pixels with MO>8° from the top (1/4) region of all 80 investigated joints in a 84CTBGA.f The data from (e) after averaging symmetry-related positions in the array as indicated in the inserts in (e-f).g Data from individual joints in the cyan regions of the insert.

Figure 6 .
Figure 6.Influence of tin grain orientation on the experimental damage and simulated stored energy density (SED) for single grain joints.a Measured area fraction with MO>8° versus the angle between Sn⟨001⟩ and Y for all 15 single grain joints studied.b Calculated connection between in-plane CTE and the angle between Sn⟨001⟩ and Y, compared with the in-plane CTE of the package.c Simulated 3D SED distribution and measured misorientation distribution in the 2D cross-section of three joints.d Correlation between simulated peak SED and experimental area fraction with MO>8°.

Figure 7 .
Figure 7. EBSD and damage analysis of typical microstructures.a,b Five examples spanning the range of microstructures observed after thermal cycling: a) Sn IPF orientation maps with MO>8° removed and measured unit cell wireframes, b) misorientation maps defined with reference to the original major Sn orientations.c Damaged area fraction values for different types of microstructure.d Area fraction with MO>8° increases with in-plane CTE mismatch to package when single orientation at top (based on experimental 2D sections).e Area fraction with MO>8° increases with increasing maximum in-plane ∆CTE between Sn grains when multiple original orientations at top.

Figure 8 .
Figure 8. Calculated connections between tin orientations and CTE mismatch.a,b Two single-grain examples: Sn IPF orientation map, unit cell, 3D CTE ellipsoid and 2D CTE in the substrate plane.c A beachball joint coloured by angle between Sn⟨001⟩ and Y plus cyclic Sn unit cells, d the geometric model at different viewing angles, e corresponding 3D CTE ellipsoids in the EBSD mapping coordinates, f 3D CTE triple ellipsoids at arbitrary viewing angle and 2D CTE section in the substrate (x-z) plane.g Maximum in-plane CTE mismatch in the substrate plane between tin grains in single and double ring cyclic twins calculated from 50,000 random orientations and plotted in Euler space, h histograms of the max in-plane CTE mismatch from (g).On all in-plane 2D CTE plots (with grey background), the dashed ellipse is the in-plane CTE of the package.

Figure 9 .
Figure 9. Direct comparison between experimentally-measured microstructure and damage and the multi-scale thermal cycling model of the same joints.a Experimental EBSD IPF orientation maps after thermal cycling.b,c 3D geometrical models before thermal cycling, input into the CPFE model.Yellow dashed lines in (c) are the sectioned planes in (b) which are at the same location as in (a).d Simulated stored energy density distribution after 10 cycles.e Simulated misorientation distribution after the 10 cycles.f Experimental misorientation distribution after 7580 thermal cycles.(a) and (f) are experimental results.(b)-(e) are CPFE modelling.

Figure 10 .
Figure 10.Calculated stored energy in selected joints.a Plots of stored energy density G along the perimeter of the top interface, i.e. path A-A' in Figure 9(d).b Representative stored energy G* that is averaged within top 10% range around the peak positions in (a).