Learning Localized Spatial Material Properties of Substrates in Ultra-Thin Packages Using Markov Chain Monte Carlo and Finite Element Analysis

Thinned silicon dies and thin substrates using thin core and coreless structures have enabled thin packages. For robust manufacturing and reliability of these parts, solving the warpage problem is key. While current finite element methodologies can provide some insights at the design stage, these simulations are only as accurate as the inputs such as the material properties and the stress-free temperatures. Electronic substrates are especially challenging to characterize and model as they are laminates consisting of a core with layers of resin and metal lines on either side. In this work, a hybrid approach using Markov Chain Monte Carlo (MCMC) and Finite Element Analysis (FEA) is used to learn the spatially varying properties of the substrate from Digital Image Correlation (DIC) measurements of the warpage. The analysis is carried out at room temperature and at an elevated temperature point. Image analysis on electrical artwork is also carried out to correlate the material properties to the substrate metal density. These results will be useful to package and substrate designers to understand how material properties vary over the substrate and how temperature and metal density affect material properties so that robust design for future packages to minimize warpage can be initiated by careful routing of metal lines depending on the locally desired properties of the stack.


I. INTRODUCTION
As consumers call for high performance mobile devices like laptops in thinner form factors, there is a push in the industry to make the electronic packages that go into them thinner. These thin packages need to be designed carefully as high warpage in these parts can lead to assembly and reliability issues like bump bridging and solder extrusion [1], [2]. Finite element analysis (FEA) would be a good tool for warpage prediction if accurate material properties were available. One study has shown that variation in package warpage measurements can be directly attributed to variation in the material properties of the bismaleimide triazine (BT) substrate [3]. The electronic substrate is inherently a composite material consisting of a core and alternating layers of The associate editor coordinating the review of this manuscript and approving it for publication was Jiajie Fan .
Ajinomoto Build-up Film (ABF) and copper metal lines on either side. The result is an anisotropic, spatially-varying, and temperature-dependent material that is cumbersome to characterize [3], [4], tedious to calculate analytically [5]- [7] and complex to model [8]. Several authors have demonstrated that FEA simulations with substrates modelled using trace mapping is the best method to model warpage accurately [9]- [13]. Trace mapping works by first importing the electrical artwork and then 'superimposing' it onto the elements of the metal layer in the FEA model. As a result, the volume fraction of the metal for each element is calculated and material properties are assigned accordingly. Regions of higher metal density would be assigned different material properties than regions of lower metal density. However, this method requires that the electrical artwork be available. During the design stage, this method is not feasible as the electrical artwork is often not even started. Clearly, an efficient method to determine VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ substrate material properties is required where electrical artwork is not a prerequisite. In one of the recent studies by Brandt et al., material properties of different layers and interfaces in the solar cell have been predicted with good accuracy using a methodology which combines Bayesian inference with a physics-based model [14]. A similar approach is attempted here to address this problem. The results achieved using the Bayesian inference are very promising as will be evident in the sections that follow. In an earlier work of the team, we proposed a hybrid approach of using Markov Chain Monte Carlo (MCMC) and FEA to learn the material properties of the electronic substrate at locations along the package diagonal, which showed that the simulation error reduced by as much as two orders of magnitude [15].
In this work, the approach is extended to learn the material properties across the whole substrate and at two different temperatures as well. Understanding how the material properties of the substrate vary spatially and result in the observed warpage measurements will add to our understanding of the warpage phenomena and help us design better packages. In addition, the learned material properties are then correlated with the metal densities at different regions over the substrate to gain a better understanding of how metal densities affect material properties and warpage.

II. METHODOLOGY
The Bayesian inference approach uses each new evidence (experimental observation) to update the probability of a hypothesis. For this problem, the warpage model parameters are represented as a random vector (x) with a prior distribution → p 0 (x). Using Bayes' theorem, the posterior distribution of parameters given the observation (or likelihood), where ϒ is the observed data. Therefore, the posterior considers the observed data and refines our hypothesis about the distribution of the parameters. MCMC is a subset of this approach where evidence is collected by sampling the probabilistic space randomly [16], [17]. In this work, FEA simulations were carried out using random vectors of parameters generated by MCMC random sampling. The likelihood, L(ϒ|x) is obtained from where SS is the sum of squared error between the measurements and the FEA model output, ν is the standard deviation of the measurement errors and n is the number of parameters. The MCMC analysis is carried out on MATLAB using a toolbox developed by Haario et al. [18]. A short introduction to MCMC can be found in our previous work, where this approach is used to learn the material properties of the stiffener and different regions of the substrate [15].
As the objective of this work is to learn the spatially varying material properties of the substrate at two different temperatures, experimental measurements of warpage over the whole substrate and at two temperatures are required. In the next section the warpage measurements from digital image correlation (DIC) that are used in the MCMC analysis as the experimental observations (likelihood) will be presented.

III. DIGITAL IMAGE CORRELATION
DIC is a deformation measurement tool. A typical DIC setup consists of the camera to image the sample at certain time intervals, the speckled sample under external force or temperature loading and a computer for image correlation, speckle tracking and deformation measurement. The 2D version of this tool consists of one camera and can measure the inplane deformation only. On the other hand, a 3D setup of this system consists of two or more cameras. As each camera captures different views of the sample, image processing and tracking can be used to measure in-plane and out-of-plane deformation at once, enabling curved surfaces to be measured as well. A speckle pattern on the sample is typically required for accurate measurements [19]. 3D DIC measurements were carried out on three package samples at two temperatures sequentially → 25 • C and 60 • C. The temperature ramp rate used was about 40 • C/min to simulate the fast ramp rate of a typical reflow profile. The package sample has a body size of less than 40 × 40 mm 2 with a total thickness of less than 1.5mm. The sample consists of a silicon die assembled onto an electronic substrate with a stiffener ring to control the warpage. In turn, the electronic substrate is a laminate consisting of the core and alternating layers of metal line and ABF. As we are interested in measuring the warpage of the whole package, the speckle pattern was applied on the back of the substrate. A sample speckle pattern can be found in our earlier publication [8]. The contour plot of a typical DIC result is shown in Fig. 1. Out-of-plane deformation of these parts was measured over four lines -two diagonal lines (AC and BD), one horizontal (FH) and one vertical line (EG) as shown in Fig. 1, in order to fully characterize the warpage of the part.
Line plots of the warpage measurements are shown in Figs. 2-4. Fig. 2 and 3 shows the warpage of all three samples (referred to as S1, S2 and S3), along each line, at room temperature and the elevated test temperature of 60 • C respectively. The mean of the measurements from the three samples is represented by the dashed black line in both figures. From Fig. 2, in general, there is good agreement between the measurements from the three samples at room temperature, with the largest variability observed along the horizontal line F-H. Fig. 3 indicates that greater variability between the samples is observed at the higher temperature. This could be the result of different amount of residual stresses accumulated in each package during the assembly process, resulting in different amounts of stress relief and hence, different warpages.
For the MCMC analysis, the observation (likelihood) data used is the mean of the DIC measurements from the three  samples. This is shown in Fig. 4. It is clearly evident that warpage is higher at room temperature and decreases as temperature increases. Similar warpage shapes and trends are observed along all four lines.

IV. FEM AND MCMC MODEL DETAILS A. FINITE ELEMENT ANALYSIS (FEA)
Typical FEA analysis for package warpage assumes that the package is symmetric about the two axes. Using symmetrical boundary conditions, the model size can be reduced by four times. However, in this case, there can be no assumption of symmetry as the objective is to determine material properties over the entire substrate. Therefore, a finite element model of the whole package was built on ANSYS v19. The model consists of the silicon die, substrate, underfill and stiffener, as shown in Fig. 5(a). The loading condition used was a single step temperature ramp down from the stress-free temperature of 150 • C to the temperature of interest (25 • C or 60 • C). The stress-free temperature was determined from warpage measurements, where the magnitude of the warpage is the smallest.  In order to learn the spatially-varying material properties of the substrate, the substrate was divided into nine parts as shown in Fig. 5(b). The model was parameterized such that each of the nine sections of the substrate can be assigned its own unique in-plane coefficient of thermal expansion (CTE) value. This approach of assigning each section its own unique CTE will allow the MCMC analysis to learn the material properties of the substrate that reflect the variation of copper density across the substrate.
A mesh convergence study is carried out to optimize the run-time of the warpage model to about 15 sec without compromising on the warpage simulation accuracy. The model is post-processed to extract warpage at the same locations as in the DIC measurements i.e. along the two diagonals and the central vertical and horizontal lines and output these values to a text file. Warpage contours from the simulations and  the lines along which the warpage is extracted are shown in Fig. 6. Fig. 7 shows the initial warpage plots from simulation when assuming uniform CTE of 14 ppm/ • C for the whole substrate. As expected, there is a large discrepancy between the measured and simulated results. The total squared sum difference between measured and simulated results for the four curves is 9.9 x10 −2 .

B. MARKOV CHAIN MONTE CARLO (MCMC)
For this analysis, the equations governing the MCMC are similar to those in our previous work [10]. The only difference is that we are trying to fit to the measurements over four discrete lines, instead of just one. In order to accomplish this, the sum of squared error for each line is added together.
The relationship between the displacement y i,l and position p i,l is given by the FEM model in Eqn 3. Subscript i denotes Here, x = (x 1 , . . . , x n p ) is a vector of n p parameters of the model. Given that the measurement of displacement data along each line Υ = z i,l = y i,l + ς i,l , i = 1, . . . , N , l = 1, . . . , 4 and errors, ς i,l , are available, the total sum squared error is now given by:

A. NINE-PARAMETER MODEL: IN-PLANE CTE AT DIFFERENT LOCATIONS ACROSS THE SUBSTRATE AT ROOM TEMPERATURE
A nine-parameter MCMC is executed to determine the spatially varying material properties in the substrate. This analysis is run for 2000 iterations. The total time taken for this analysis is about 9 hours. The convergence plot is shown in Fig. 8. Note that all the nine CTE values seem to have converged. The corresponding probability density plots are shown in Fig. 9. The best estimate for CTE1 through CTE9 are determined to be 16.0, 14.8, 12.0, 10.7, 8.2, 11.2, 12.9, 15.1 and 14.9 ppm/ • C respectively by averaging the latter values of the MCMC trials. These values are represented graphically in Fig. 10. Note that there is quite a wide variation in the CTE in the nine regions with variations of up to 50%. Similar to earlier work, it is worth noticing that the CTE increases from the center to the corner of the substrate.
Displacement plots of the experiment and simulation are shown in Fig. 11 using the optimum values of CTE obtained above. Also shown on this plot are the 5% and 95% quantiles for the 2000 iterations. The total squared sum of the difference between the experimental and simulated displacement curves across the four lines has now decreased from 50166 VOLUME 8, 2020  9.9 × 10 −2 in the previous section to 3.2 × 10 −3 , much more than an order of magnitude. This implies that an extremely good fit is achieved with the learnt material properties.

B. NINE-PARAMETER MODEL: IN-PLANE CTE AT DIFFERENT LOCATIONS ACROSS THE SUBSTRATE AT ELEVATED TEMPERATURE
The MCMC analysis is repeated to determine the material properties, this time at 60 • C. The analysis is run for 1500 iterations and the total time taken for this analysis is about 6 hours. The convergence plot in Fig. 12 shows that the nine CTE values have converged satisfactorily. The corresponding probability density plots are shown in Fig. 13.
The best estimate for CTE1 through CTE9 is determined to be 16.    Just as in the results for the room temperature, the CTE is again found to increase from the center to the corner of the substrate. This trend will be investigated further in the next VOLUME 8, 2020  section when we correlate the learnt material properties to the metal layer density. Displacement plots of the experiment and simulation are shown in Fig. 15 with the 5% and 95% quantiles for the 1500 iterations. The total squared sum of the difference between the experimental and simulated curves across the four lines is 3.7 × 10 −3 . This is in the same range as that for the room temperature analysis. Again, an extremely good correlation between the simulated and experimental curves is achieved with the learnt material properties.

C. CORRELATION BETWEEN LEARNT MATERIAL PROPERTIES AND METAL LAYER DENSITIES
Though learning material properties of the substrate is found to be greatly useful in modeling the current design of the substrate and warpage accurately, the learnt material properties will not help in future robust design unless they can be correlated to the metal layer densities. Unfortunately, the software that generates the electrical artwork, Cadence , is only able to output the metal densities of the whole layer and Significantly lower SS of 3.7 × 10 −3 is obtained between these two results. The 5% and 95% quantiles for the 1500 simulation runs are also plotted (dashed black lines).
not subsections of the layer. In order to get this information, a code in Matlab is written to carry out image analysis on the electrical artwork of each of the eight layers in the substrate and compute the metal densities based on counting of the white and black pixels. The metal density results from the code are then validated using similar data obtained from Cadence . There is an error of about 6% in the copper density calculated by the Matlab code, as compared to the Cadence report for the first metal layer. This is possibly caused by the reduced resolution of the narrow trace lines in the imported image. Other layers with fewer lines have their areas calculated with an error of less than 3%.
With the validated code, image analysis is carried out to determine the copper density in each of the 8 metal layers, in each of the 3 by 3 cells (Fig. 5(b)) in the substrate, totaling 72 different regions. The metal density of each cell is determined by averaging the metal densities of the eight layers at that cell location. The plot of the extracted localized CTE determined from the MCMC analysis against the metal density determined for each cell from the Matlab image analysis is shown in Fig. 16. Clearly, as the copper density in the cell increases, the corresponding effective CTE decreases. This is because the metal layer consists of copper, whose CTE is 17 ppm/ • C, and ABF material of higher CTE (CTE > 20 ppm/ • C). Therefore, with a higher metal content, the volume average of CTE is expected to be lower. A similar trend is observed at higher temperature as well. However, as indicated by the large scatter of the datapoints in the plot, the average copper density is not the only factor governing the CTE at each cell. Other factors such as copper layer pattern and the resulting anisotropy of the sample, residual stresses in the substrate and edge effects, to mention a few, also contribute to the local CTE of the sample and are not accounted for explicitly in this work. The substrate consists of three materials → core, ABF and copper metal. At the central cell of the substrate, the highest metal density and lowest CTE is found. This is because the core is significantly thicker than the ABF and metal layers. At the substrate center, the low CTE of the core material (CTE < 10 ppm/ • C) easily dominates the material property of the cell. The outer cells in the periphery have higher CTEs, as in those regions, the accumulated effect of the ABF and copper layers plays a larger role in the deformation.

VI. CONCLUSION AND RECOMMENDATIONS
In this study, a hybrid approach of FEA and MCMC is proposed and used to calculate the spatially varying material properties (specifically, CTE) of the electronic substrate at two different temperatures. When the optimized local values of CTE from the MCMC routine are used in the FEA, excellent agreement with experimental measurements using DIC is obtained, indicating that the CTEs calculated using the methodology are indeed very accurate. The metal densities for each cell in the substrate were determined through image analysis of electrical artwork. It was found that the CTE is inversely proportional to the metal density. These results will be very useful for package designers to design future packages with better spatial planning so as to ensure much lower warpage. Though accurate, the hybrid approach of FEA and MCMC is still quite time-consuming, especially for larger FEA models.
In the future, we will attempt to replace the FEA model with a neural network as demonstrated by several authors [20], [21]. The neural network is a computational tool that learns the weights and biases between the input and output layers and then uses the learned parameters to make efficient predictions for other values in the parameter design space [22], [23], resulting in a significant speed up of simulation time. This would enable the number of MCMC iterations to increase well beyond 2000 in a fraction of the time.
In addition, this would also enable us to increase the number of parameters included in the analysis to better understand the factors (copper pattern, material property changes in constituent materials, residual stresses, etc.) that affect CTE and hence the warpage of the package.
The hybrid approach of using the MCMC and neural networks is the first step in realizing an inverse design model for electronic substrates. The next step would be to develop a model to link substrate characteristics like the number of layers and the material properties of the constituent materials to the material properties learned from the MCMC. With this framework in place, designers would be able to input a warpage profile and the model would return the necessary substrate characteristics (and perhaps the stack design, material choice and metal pattern density itself) that are required to achieve the desired warpage contour. ACKNOWLEDGMENT AMD, the AMD Arrow logo, and combinations thereof are trademarks of Advanced Micro Devices, Inc. Other product names used in this publication are for identification purposes only and may be trademarks of their respective companies. 2019 Advanced Micro Devices, Inc. All rights reserved.