Development of a Simplistic Method to Simulate the Formation of Intermetallic Compounds in Diffusion Soldering Process

A simplistic simulation technique has been developed for computing the individual intermetallic compound (IMC) thickness which is formed in substrate-solder (Cu-Sn) systems during the diffusion soldering process in high-temperature power electronic applications. The method requires the time-dependent temperature profile for the soldering process and the growth rate parameters (e.g. concentration gradient, diffusion coefficient, activation energy, etc.) for the development of IMC layers as input. The method is suitable for predicting the thickness of an intermetallic phase layer during the diffusion soldering process. As such, it can be used in high-temperature power electronic application’s solder processing to enhance the reliability and lifetime of solder interconnections by allowing the control of the thickness of IMC layers. The method is demonstrated for IMC growth between pure copper as substrate and pure Sn as solder material. The growth behavior of the IMC layer is increased with increasing temperature over time according to the Arrhenius theory in the temperature range between 24°C to 260°C. To simulate the formation of IMC thickness in diffusion soldering interconnections, a simplistic way has been attempted using the popular commercial finite element simulation tool Comsol Multiphysics and scientific computing application ‘Matlab’. By means of transient thermal input, the diffusion-controlled intermetallic phase formation is simulated here. Few assumptions are taken care of this simulation process, for example, no convection, no reaction, solid-solid diffusion, no the pressure effect on the computational domain.


Introduction
Among all the interconnection technologies soldering still is the predominant technique for high-temperature power electronics applications. In the soldering process of Cu-Sn system, the Cu (Substrate) is diffused into liquid solder (Sn) which results in the nucleation of intermetallic compounds in the bulk of the solder material. During that time, molten solders (Sn) react with Substrate (Cu) films to form Cu 6 Sn 5 (ƞ) and Cu 3 Sn (ε) intermetallic compounds, after which the solder joints are established. These layers grow fairly rapidly during the diffusion soldering process. At the early stage of the reactions, the Cu 6 Sn 5 phase with a monoclinic structure formed first, followed by the appearance of the orthorhombic Cu 3 Sn phase. The morphology of Cu 6 Sn 5 IMC appears like scallops, while that of Cu 3 Sn is layer-type. The diffusivities of Cu in Cu 6 Sn 5 and Cu 3 Sn IMCs are much less than those in molten solders. It is reported that the Cu 6 Sn 5 IMC is the dominant phase during transient liquid phase soldering at temperatures below 260°C. Only after the solder is completely consumed, the Cu 3 Sn layer starts to grow thicker. [1] The growth is controlled by diffusion of the constituent atomic species to the reaction sites. While the presence of an intermetallic layer is evidence of a metallurgically sound solder joint, the resulting thickness of the layer is of particular importance to the integrity of a joint and, thus, to the reliability of the electronic device. In standard electronics, only thin layers of intermetallic compounds are realized for good solder connections because thick growth of intermetallic compounds leads to faster crack propagation in layer and failure of connection. But for high-temperature applications in power electronics, to avoid the problems with low-melting Sn-Cu based solders, a relatively more solder layer (Sn) between the chip and direct copper bonded (DCB) is to be transformed into intermetallic compounds as they have better thermal stability. In addition to this, intermetallic compounds are inherently hard and brittle and, therefore, have a low fracture toughness. The objective of this research work to develop a simplistic numerical method for calculating the thickness of an intermetallic layer formed in solder/substrate systems during the solder joint fabrication is introduced.

Phase Field Method
Phase-field method which solves the interfacial problems. At the interfacial region, it replaces the boundary conditions by a partial differential equation to evaluate a secondary field. Recently, the phase field formulations have been used to simulate soldering reactions. [2] In each of phases, the filed considers two different values (for example +1 and -1). At the interfacial region, it allows horizontal changes between both values, afterward, which can diffuse [3] with a finite width. It is one of the useful methods to simulate the evolution of microstructural.
In the phase field simulation, the microstructures are composed of different size of grains. The shape and common circulation of the grains is represented by the functions called phase field variables which are continuous in space and time. The below figure (1) Illustrates the interfacial region where the diffused interface is shown in figure (1-a) and the sharp interface is shown in figure (1-b) including the indication of different phases. The width h is also indicated in the figure which creates the difference between diffused and sharp interface. Within the grains, where the values of phase field variables have remained constant, which are relevant to the structure, positioning, and conformation of the grains. Grain size valuation or the interface positioning as a function of time is implicitly provided by the development of the phase field variables.
The temporal evolution is depicted by a number of differential conditions, which are possibly solved by different numerical methods.
In phase field model, the free energy function, composition constraint with the condition of local thermodynamic equilibrium, multiphase field function of time in term of relaxation equation and diffusion in respect of mass flux equation is to be defined [2]. The total free energy function is defined by the below equation (1).
Where fi is chemical free energy density in phase i, c i is the phase composition, ԑ ij is the gradient energy coefficient and ω ij represents the potential energy barrier between two phase ϕ i and ϕ j . In the phase-field modeling, the evolution of multiphase fields [4] as a function of time follows the following relaxation equation (2).
Where Np is the number of coexisting phases, M ij represents the mobility of the interfacial region. The derivative of total free energy F with respect to ϕ i and ϕ j is defined in the below equation (3).
For mass diffusion, the mass flux equation is set with respect to diffusion coefficient, free energy density and composition gradient for the multiphase system. The transformation of the flux equation into the diffusion equation is conceivable by using the spatial differentiation with diffusivity which fluctuates as a function of phase fields. The following equation (4) is given as the diffusion equation.

Monte Carlo Simulation Method
To solve the magnetic domain problem, Potts and Ising models have been used widely and Monte-Carlo grain growth simulation begins from that model. The Potts model permits multiple states namely Q states for individual particles in the entire system, on the other hand, Ising model comprises two spin states called up and down. For Modeling and simulation of material characteristics, for instance, the texture evolution and grain growth simulation, the Potts model has widely been used by several researchers.
Throughout the Monte-Carlo grain growth simulation, a 2D Monte-Carlo lattice is mapped by a continuum microstructure [5]. It can either be a rectangular or triangular lattice. [6] In each lattice location, an integer number between Q and 1 is assigned for initializing the lattice. In the system domain, Q denotes the total number of orientations. Two neighboring locations with dissimilar grain arrangement numbers are observed as being disconnected by a grain boundary and individual pair of different adjoining locations contributes a piece of grain boundary energy to the whole system. A collection of grain location having a similar arrangement number and bounded by grain boundaries are considered as a grain. The complete energy of the entire system is computed by the below equation (5) which is related to the grain boundary energy contributions through all the grain locations.
Where, summation of i is the overall Monte-Carlo locations in the system, the totaling of j is overall nearest neighbor locations of the site i, and δ ij is the Kronecker delta. This method simulates the grain development iteratively.

Multiscale Simulation Method
The combination of the Monte-Carol simulation method and the Finite Element method is called the Multiscale simulation which simulates the growth of intermetallic compounds by means of microstructural evaluation. The method is developed under the principle which says the stored energy of solder material is increased steadily throughout the each thermal cycle. The recrystallization is begun when an energy critical value is introduced. New grain growth and the nucleation are achieved by releasing the stored energy which consumes the strain hardened matrix of high dislocation density gradually. For modeling the macroscale in homogenous deformation, the finite element method is used in this particular simulation. On the other hand, modeling the mesoscale microstructural evolution Monte-Carlo simulation model is applied. A correlation between the Monte-Carlo simulation time and the real time is recognized while comparing with the in-situ experimental result. The effects of intermetallic phases (for example Cu 3 Sn and Cu 6 Sn 5 ) on recrystallization in the solder matrix are also been encompassed in the multiscale simulation. [4,7]

Dual Phase Lag Model
The dual phase lag model was originally developed for the diffusion medium to compute the formation of the interfacial compounds in thin films and metal matrix composites. The kinetics of interfacial phase formation in Multiscale simulation model and solder interconnections also applied the Dual phase lag model to demonstrating the intermetallic phase growth in solder joints, e.g. intermetallic growth in a pure Sn-Cu system at 170°C. Some researchers have further applied this model to simulate the growth of the interfacial intermetallic phase for the Sn-3.5Ag (solder) & Cu (substrate) and Sn-3.5Ag-0.7Cu (solder) & Cu (substrate) systems during the reflow soldering. The Dual Phase Lag model also be used to describe some of the singularities of the intermetallic phase formation in solder interconnections, however, it provides only a 1-D simulation model in nature. [8]

Choosing the Equation of the Physics
During intermetallic formation, inter-diffusion which is driven by the atomic concentration gradient can be described by the second Fick's law. In COMSOL Multiphysics, the Transport of Diluted Species module is used to compute the concentration field. Transport and reactions of the species dissolved in a gas, liquid, or solid can be computed in this module. As described above, the reaction is assumed to be zero. So the transport of diluted species module is arranged according to the assumption mentioned above.
The COMSOL model with Transport of diluted species module follows the generic diffusion equation which has the same structure as the heat conduction equation and is governed by the following mass transfer equation (6). As the convection is assumed to be zero, the equation (6) is prepared according to the assumption.
Where G is the reaction rate. For the simplicity of the model, the reaction rate will be assumed to be zero here (G=0), then the equation (6) would be in below form of equation (7).
Equation (12) is the governing transport equation for the diffusion soldering process. Where D is the growth rate constant which is equivalent to diffusivity if intermetallic compound formation is diffusion controlled. Here in this study, it has been assumed that the TLP soldering is diffusion controlled. According to the Arrhenius equation (8) as below, the diffusivity is determined by the activation energy Q (KJ/mol), the diffusion coefficient Do (m2/s), and the absolute temperature T (K).

Model Setup
The steps are taken to set up the simulation model is given below: (a) Geometry setup and Material selection (b) Meshing (c) Boundary conditions (d) Running the simulation.

Geometry Setup
To reduce the computational complexity and time, a simple 2D geometry was considered as a part of preparing the simulation model to see the concentration gradient due to the diffusional variation. In the geometry, the top and bottom domain are selected as tin (Sn) and copper (Cu) respectively. The length/width of both Sn and Cu is selected as 5 mm, on the other hand, the height/thickness of both Sn and Cu is 40 µm. Material data plays an important role in finite element simulation. [9] Pure Tin (Sn) and Copper (Cu) was selected to simulate the intermetallic phase formation in this research work. In the real application some other soldering materials are used but in most of the tin (Sn) based solder alloy where the rest of the composition percentage is less than 3%.

Meshing
In order to analyze the concentration gradients, computational domains are splitted into smaller subdomains. The governing equations are then discretized and solved inside each of these subdomains. For this work, finite element method (FEM) was used to solve the approximate version of the system of equations where quadrilaterals mesh was chosen for the computational domain.
To mesh the model, the user-defined mesh parameter option was selected in COMSOL and meshed according to the maximum element size. The maximum element size determines the minimum number of nodes for the model. The more important factor for the transient model is the maximum time stepping size for the model. As a time-dependent COMSOL model, the maximum time stepping size determines the smoothness and accuracy of the data. There are three types of meshing with different element size according to the computational domain was tested. All are given below with the number of elements, element size and element ratio according to the X and Y axis of the computational domain. All three types of meshing are illustrated graphically with geometrical setup in the below figure 3. Meshing parameter such as the number of mesh element, element size, and ratio are explicitly shown in table 1 which are examined later.  After testing the different meshing, type 3 was implemented in the current simulation model. Type 3 mesh gives better convergence in the result because it has optimum element size in the computational domain.

Boundary Condition
Initially, 100% Cu and 0% Sn was selected in the left computational domain. On the other hand, 100% Sn and 0% Cu was selected in the right computational domain. No flux boundary condition was implemented at the edge of the domain. One of the important input is at the interface between two media. At the boundary, two conditions are specified. A condition that links the dependent variable in the two regions and a condition that links the flux of the dependent variable in each region. For this simulation, it was assumed that the continuity of the flux across the boundary and the dependent variable.
Flux boundary conditions are interesting in a number of contexts ranging from multiscale simulations to simulations of concentration gradients in dilute species at microscale systems. Present research work introduced the local boundary conditions that allow the fluxes of momentum and energy to be simultaneously specified either as free parameters or in order to simulate an external parameter. [10] Such boundary conditions may be required in a wide variety of contexts where the diffusional variation in an inter-diffusion system needs to be simulated on the atomic level. Here, flux boundary conditions are interesting in particular for the design of concentration simulation of an inter-diffusion scheme where heat transfer and diffusion simulations are coupled by flux exchange algorithms. For these schemes to be conservative the flux leaving one system must enter the other. Parameter preparation inside COMSOL is one of the important parts of the running simulation. In below Table 3, important parameters are listed with their names, value, unit and the relevant description. [11,12] Table 2 represents the parameter used in the current simulation process. Those important parameters are imported into the COMSOL computational module which are the necessary constants and variables to simulate the concentration gradient of Cu-Sn system.

Importing Temperature Profile as a Global Function
After setting all the parameters, COMSOL imports MATLAB generated temperature profile as a global function. An interpolated function is used in current time-dependent study. An interpolated function can create a smooth function which allows the dependent variable is plotted against the independent variable. [13] Matlab script generates transient thermal input which I call it here as temperature profile. In the profile, Cu-Sn system is heated from initial temperature to a peak temperature with a certain heating gradient. Then it is heated at constant peak temperature for some time which is called as handling time (t h ). After that, the system is cooled down to the room temperature with a cooling gradient. In this research work, the initial/ room temperature (T i ) is selected as 24°C. Peak temperature (T p ) is selected as 260°C. The heating (∆ H ) and cooling (∆ C ) gradients vary from 0.5 K/s. The handling time (t h ) vary from 1.5 minutes to 30 minutes. The temperature profile is shown in the below figure 4.

Setting the Diffusion Coefficient as Variable
As the soldering process is depending on diffusion between solder and substrate, the concentration changes in the Cu-Sn system also depends on the diffusion coefficient. As the diffusion coefficient is temperature dependent [14], it changes with temperature in every time steps. So, both Sn and Cu were taken as a global variable which changes in every time step during the simulation.

Simulation Process Overview
A MATLAB script was written to show the output of Cu-Sn from COMSOL in terms of concentration gradients and reforms the result in Intermetallic phase thickness. Figure 5 shows the flow chart of visualizing individual intermetallic phase recognition and formation during the simulation. The total simulation process is divided into six stages which is graphically represented in figure 6. In the beginning, a MATLAB script generates the user-defined temperature profile. After defining the inputs, COMSOL import it as a global variable where temperature changes over time. After importing the temperature inputs into COMSOL, necessary other variables and constants are defined. Subsequently, COMSOL simulates the concentration gradient using finite element method. For simulating the concentration gradients, transport of dilute species (tds) module was used in COMSOL which represents the mass transportation equation due to diffusion. The simulation output contains the information of concentration in every location of the computational domain over time. Finally, a MATLAB script was used to show the intermetallic phase formation. The MATLAB script reads the simulation file and based on the composition of two important intermetallic phases (Cu 6 Sn 5 and Cu 3 Sn), an atomic concentration range was arranged which defines the phase formation of two intermetallic compounds as well as pure Cu and Sn. An algorithm is defined in MATLAB to show the concentration change of every phase in every time steps which is shown in figure 5. In the below figure 6, the overall simulation process is drawn step by step where all the applied tools are described.

Evolution of Intermetallic Phases
The direct output from the simulation is the atomic concentration of Cu-Sn system with respect to the location and time. In the below figure shows four individual phase formation over time. In the beginning, there was 100% Sn in the right side and 100% Cu on the left side. But over time two more phase forms due to diffusion. It is commonly known that Cu and Sn react with each other at the interface [15] to form an intermetallic compound of Cu 6 Sn 5 and Cu 3 Sn. The intermetallic Cu 6 Sn 5 forms very fast and grows with scallop-like morphology because Cu atoms diffuse inside Sn through fast atomic transport phenomena. Atomic concentration profile can be plotted along the perpendicular direction of IMC interface at any time points by this way the determination of total IMC thickness is possible.
In the early stages of the diffusion process, Cu atoms dissolved into the molten Sn and form Cu 6 Sn 5 intermetallic compounds at the room temperature because Cu atoms diffuse faster in Sn which is shown in Figure 7. After that Cu 6 Sn 5 reacts with Cu and form Cu 3 Sn IMC. [6] Figure 7 shows how the different phases are evolved over time. The reaction kinetics are given below.
6Cu + 5Sn = Cu 6 Sn 5 Cu 6 Sn 5 + 9Cu = Cu 3 Sn  Table 3 describes all the parameter of the simulation study. In this case study both the heating (∆ H ) and cooling (∆ C ) gradients are kept as constant where the value of both gradients remains 2.5 K/s. The peak temperature is selected as 260°C in this particular study. The initial/room temperature is selected as 24°C. A variant is added here which is handling time (t h ) that varies from 1.5 minutes to 30 minutes. The maximum time of the simulation or total time of a study is calculated as 33 minutes at handling time of 30 minutes and the peak temperature of 260°C. The dimension of the geometry (computation domain) is kept as same (Cu=40µm, Sn=40µm) in the entire case study of this simulation work. All the assumptions are the same as described at the beginning. The intermetallic phase evolution over time is graphically represented in figure 8 and 9.   . The light red color (left side) represents the substrate (Cu) and blue color (right side) represents the pure solder material (Sn). In between two of them, another two colors are grown from bottom to top and those are two intermetallic phases. The bigger percentage or the dark red color (on the side of Sn) is Cu 6 Sn 5 and lower percentage or green color (on the side of Cu) is Cu 3 Sn. This graphic is prepared using Matlab. Basically, both of the figures have been captured from real-time animation during the simulation while data was exchanging simultaneously between COMSOL and Matlab. Figure 10 shows that intermetallic phase thickness of Cu 6 Sn 5 and Cu 3 Sn at the peak temperature of T p = 260°C having different handling time while the heating and cooling gradients are constant (∆ H = ∆ C = 2.5 K/s). The handling time varies from 1.5 minutes to 30 minutes where the maximum intermetallic thickness of Cu 6 Sn 5 and Cu 3 Sn observed as 9.9 µm and 1.6 µm at handling time of t h = 30 minutes. From figure 8, it is shown that the intermetallic compound growth of Cu 6 Sn 5 is higher than Cu 3 Sn which is actually seen in experiments. [16] Normally Cu 6 Sn 5 start to form in the room temperature when Sn (solder) and Cu (substrate) react together. On the other hand, Cu 3 Sn need more time to form even more temperature. So after a certain temperature Cu 3 Sn start to form and grow more with increasing handling time.

Conclusion
In this research work, a simplistic simulation method was attempted to develop by which we can predict the intermetallic compounds formation in diffusion soldering process [17] in the solder-substrate interface as well as the thickness of those phases with a specific time-dependent temperature profile. For the simplicity of the design of the simulation method, solder and substrate material is chosen as pure tin (Sn) and copper (Cu) in this study. Two particular intermetallic compounds such as Cu6Sn5 and Cu3Sn were described here which is grown during the diffusion soldering process of Cu-Sn system. The total simulation process was divided into three sub-process. At the first step, a MATLAB script was used to generate the time-dependent on the temperature profile. After that popular CAE tool called COMSOL Multiphysics simulated the concentration gradient of the system during the diffusion process between Cu and Sn where the temperature profile was the input function in the COMSOL global definition area. The diffusion coefficient was also taken as a variable where the depending variable was temperature. Finally, a MATLAB script was written to illustrate the intermetallic phase formation over the time period based on the result of concentration gradients from COMSOL. For simulating the concentration gradients, transport of dilute species (tds) module was used inside the COMSOL module which represents the mass transportation equation due to diffusion. The maximum thickness of the intermetallic phases was observed at peak temperature, Tp= 260°C, heating and cooling gradient, ∆ H = ∆ C = 2.5 K/s and the handling time, t h = 30 minutes where the thickness of Cu 6 Sn 5 and Cu 3 Sn were calculated as 9.90 µm and 1.60 µm respectively.