Multiphysics Analysis and Optimal Design of Compressible Micro-Interconnect for 2.5D/3D Heterogeneous Integration

: Compressible Micro-Interconnect (CMI) shows tremendous potential in 2.5D/3D heterogeneous integration due to its outstanding performance in integration, electrical isolation, and thermal management. In this work, an optimal design approach for CMIs is developed based on a coupling framework of multiphysics simulation and particle swarm optimization (PSO). In the framework, the mechanical simulation was conducted ﬁrstly to obtain the stress distributions as CMI switched from the initial state to the working position. The contact resistance between CMI and the top pad was modeled and quantitively analyzed. Subsequently, the PSO method was utilized to implement the structural optimization of CMI to improve the performance. Multiphysics simulations of both the original and optimized CMIs were carried out and compared. With the implementation of the optimized CMIs, the contact resistance dropped from 155.3 m Ω to 108.8 m Ω , which brought signiﬁcant improvement in both DC voltage drop and self-heating effect. The inﬂuence of the self-heating effect on the electrical performance of CMI is also discussed qualitatively.


Introduction
With the fast-approaching limit of Moore's law, three-dimensional (3D) heterogeneous integration technology was proposed as an alternative solution to improve the integration density and the system performance [1]. As the channel for signal transmission and power supply, 3D interconnect network plays a critical role. In past decades, solder bump was an indispensable component for systematic interconnection, and it could realize an electrical connection between different circuits or packaged integrated circuits (ICs). However, the conventional solder bump is no longer efficient to provide a reliable electrical connection between vertical stacked multi-dies, as it is challenging to achieve solder bump array with good electrical isolation at such a scaled dimension. Moreover, the rework and repairment are inconvenient or even nearly impossible for a permanently welded module of multi-die, thereby resulting in a costly implementation of a fault monitor [2]. Moreover, the stacked structure of 2.5D/3D integration leads to a significant increase in the power density, and the large temperature gradient along the vertical direction would lead to severe thermal crosstalk between layers since the solder bump array serves as an ideal heat-conducting route. The crucial thermal crosstalk would degrade the system performance dramatically.
To implement a more reliable interconnection for 2.5D/3D ICs and make the system convenient to disassemble, some mechanism-and pressure-contact-based compliant interconnect technologies have been developed. In [3], a land grid array (LGA) structure was proposed as a reliable vertical interconnect for electronic packaging. Some related research, including reliability assessments and contact resistance predictions, were carried out to accurately evaluate and predict the performance of LGA-like structures [4][5][6][7]. For

•
The contact plane of the reported compliant structure is usually horizontal without loading, i.e., only the edges of the plane contact with the pad, which causes small contact area and large contact resistance; • The interconnect resistance plays a critical role in IR drop. Due to the large contact resistance between the compliant structure and the pad, a large array should be arranged to meet the IR drop threshold, which would lead to an increased package size; • The self-heating effect caused by the large contact resistance at the interface would further enlarge the contact resistance and degrade the system power supply. At the same time, the rising temperature would increase the risk of noise mixing, thereby disturbing the purity of the supplied current; • Until now, the reported works mainly focused on the modeling, fabrication process, and parameter measurement of compliant interconnect structures. The mechanical simulation was employed to predict and evaluate the performance of designed structures [19,20]. However, due to the significant increase in power density in 3D IC, the influence of self-heating effect on the interconnect performance cannot be ignored. Multiphysics coupling simulation should be employed to precisely evaluate the performance of compliant interconnect structure.
In this study, a CMI structure with a sloping contact plane is proposed and optimized. The sloping contact would be squeezed to be horizontal to form a large contact area when the top dies are mounted, as shown in Figure 1. Optimal design of the proposed CMI was carried out based on a coupled framework of multiphysics simulation and particle swarm optimization (PSO) algorithm for the first time.

Modeling and Optimal Design Approach
In this section, detailed descriptions for the multiphysics coupling model, mathematical model for contact resistance, PSO approach, and optimal design approach are presented.

Modeling and Optimal Design Approach
In this section, detailed descriptions for the multiphysics coupling model, mathematical model for contact resistance, PSO approach, and optimal design approach are presented.

Mathematic Model for Multiphysics
In multiphysics simulations, governing equations for different physical process are coupled together with the assistance of a coupling model. In this work, the multiphysics coupling process involved elastic mechanics, heat conduction, and electrical processes. The coupling model is usually composed of a field-source direct coupling process and a field-dependent nonlinear parameters indirect coupling process. The current continuity equation, Fourier thermal conduction equation, and equilibrium equations were utilized to describe the electrical, thermal, and mechanical processes, respectively. With the assistance of corresponding boundary models, the governing equations of the multiphysics coupling process can be solved numerically with the finite element method.
The coupling process between electrical, thermal, and mechanical analysis is shown in Figure 2. The electric simulation was carried out with Dirichlet boundary conditions assigned to the top boundary of the top pad and the bottom boundary of the CMI. Then, thermal simulation was conducted with the Joule heating achieved from the electric simulation. The calculated temperature distribution is adopted in the mechanical simulation to update the pressure distribution. Based on the updated pressure, a new contact resistance was achieved. The updated contact resistance and temperature distribution were fed back to an electric simulation to start a new iteration. The coupling simulation continued until the convergence criteria had arrived. In this section, detailed descriptions for the multiphysics coupling model, mathematical model for contact resistance, PSO approach, and optimal design approach are presented.

Mathematic Model for Multiphysics
In multiphysics simulations, governing equations for different physical process are coupled together with the assistance of a coupling model. In this work, the multiphysics coupling process involved elastic mechanics, heat conduction, and electrical processes. The coupling model is usually composed of a field-source direct coupling process and a field-dependent nonlinear parameters indirect coupling process. The current continuity equation, Fourier thermal conduction equation, and equilibrium equations were utilized to describe the electrical, thermal, and mechanical processes, respectively. With the assistance of corresponding boundary models, the governing equations of the multiphysics coupling process can be solved numerically with the finite element method.
The coupling process between electrical, thermal, and mechanical analysis is shown in Figure 2. The electric simulation was carried out with Dirichlet boundary conditions assigned to the top boundary of the top pad and the bottom boundary of the CMI. Then, thermal simulation was conducted with the Joule heating achieved from the electric simulation. The calculated temperature distribution is adopted in the mechanical simulation to update the pressure distribution. Based on the updated pressure, a new contact resistance was achieved. The updated contact resistance and temperature distribution were fed back to an electric simulation to start a new iteration. The coupling simulation continued until the convergence criteria had arrived.

Calculation of Contact Resistance
As one of the most important power supply network metrics, the DC IR drop plays a significant role in the reliability of an electronic system. The contact resistance between CMI and the top pad is the major source of the CMI resistance. To implement the performance analysis and optimization of the proposed CMI structure, an analytical model of contact resistance was employed.
According to the Theory of Contact [21], the contact surface is usually not ideally smooth, and the actual contact area is far smaller than the geometric area. Current flows through the contact spots on the rough surface when two surfaces are squeezed with each other, resulting in the so-called constriction resistance [22]. Besides, to accurately describe the contact resistance, complex physical factors, such as geometrical parameters, material properties, the degree of surface cleanliness, the ambient humidity, and temperature, should be considered. In this study, the contact resistance model proposed in [23] was utilized, i.e., where H c is the microhardness of the softer material, σ contact is the harmonic mean of the contacting surface conductivities, m asp is the average slope, and P denotes the pressure between two contact objects. To calculate the contact resistance, numerical mechanical simulation was carried out to obtain the pressure distribution on the interface of the deformed CMI, and then the pressure on each mesh grid was input to Equation (2) to achieve the contact resistance. The contact resistance of the CMI can be treated as the parallel connection of the mesh grid contact resistances, i.e., where n denotes the grid number and R ci is the contact resistance of the ith grid.

PSO Approach
The particle swarm optimization algorithm proposed by Eberhart and Kennedy was employed to optimize the CMI as it is much more efficient in implementation and convergence than wolf swarm and ant swarm optimization (ASO) algorithms [24]. Subsequently, the detailed principle of the PSO approach is introduced.
The swarm intelligence optimization algorithm, inspired by the social behavior mechanism of animal population, is a distributed solution strategy. As one of them, particle swarm optimization algorithm is an approach imitating the cooperation and information sharing behavior of birds in foraging. To implement the PSO approach, the number of particles and swarm's dimensions, i.e., the number of target variables and their range, should first be defined. Then, the initial position and the velocity for all particles within the pre-defined range are randomly set. Starting from the initial guess, the optimization process iterates until the pre-defined threshold values are achieved. During the iteration process, the calculated individual optimal result (p best ) and global optimal result (g best ) are retained. A cost function or fitness function is defined to assess the quality of a certain set of parameters. Here, the fitness function of contact resistance is chosen to evaluate the output. In each iteration, the particles' positions and velocities are updated with the returned fitness function value. Assuming where w is the inertia factor, c 1 and c 2 are the learning factors, r 1 and r 2 are any random numbers between zero and one, and best denote the optimal position of the ith particle and best fitness function value in the first (m − 1)th iteration. The second and third items in the right-hand side of Equation (3) are the self-cognition item and group cognition item, respectively. The whole process of optimization would stop when the convergence condition is achieved. The implementation of the PSO approach is presented in Algorithm 1.

Optimal Design Approach
In this subsection, the implementation of optimization approach is presented. Based on the COMSOL platform, the PSO algorithm was programmed, and the interaction between the PSO program and CMI modeling program was set up. Figure 3 illustrates the detailed flowchart of the co-simulation process where MATLAB codes serve as a tool for data processing and COMSOL serves as the tool for multiphysics simulations. The detailed implementation steps are as follows:

Optimal Design Approach
In this subsection, the implementation of optimization approach is presented. Based on the COMSOL platform, the PSO algorithm was programmed, and the interaction between the PSO program and CMI modeling program was set up. Figure 3 illustrates the detailed flowchart of the co-simulation process where MATLAB codes serve as a tool for data processing and COMSOL serves as the tool for multiphysics simulations. The detailed implementation steps are as follows:  Initialization: The optimized parameters include the width, the slope of the contact plane, and the fillet radium of CMI. The ranges of the optimized parameters are defined first, and the ranges of velocities are defined to be ten percent of the corresponding parameters. After that, the initial guess for all particles' positions and velocities are allocated randomly.
Parameters transferring: In this step, all particles information is transferred to the program of modeling and the corresponding CMI is built.
Modeling and simulation: The mechanical characteristics of CMI are simulated and the pressure distribution along the contact interface is fed back to the optimization program.
Fitness function computing: Based on the feedback data, the value of the fitness function is updated to assess the quality of the related parameters.
Renew the particles' information: For all particles, if the achieved contact resistance is smaller than previous one, the position and p best are updated. The global optimal result g best is also updated if better parameters are achieved.
Judgement of convergence: All particles proceed to the next iteration until the convergence condition or maximum iteration is reached. The best result is chosen to be the candidate.

Modeling and Initial Simulation
Generally, the contact plane of the conventional CMI is in a horizontal state without loading. Once the upper die or core is mounted, the CMI is deformed and only a sole edge is attached to the pad, which leads to a small area of electrical contact and, consequently, a large contact resistance. In this study, a CMI based on the structure in [17] is proposed, and the significant improvement is that the contact plane would stay sloping in the free-standing state and remain in a horizontal state when subjected to the gravity of the upper layers. In this section, the geometrical model is built and simulated in a mechanical environment.

Geometric Modeling
For convenience, the geometric model of an isolated CMI was built on the platform of COMSOL, as shown in Figure 4. The part of cube above the CMI can be regarded as the bottom metal pad of the upper layers. The contact pressure of the contact surface can be extracted by the interaction between the pad and the CMI, and it plays a key role in the process of establishing the mathematical model of contact resistance. The materials used in the construction of the CMI were as follows. The pad was made of copper. For the CMI, which is the most critical part and plays the role of mechanical support and electrical connection, the material usage is particularly important. The maximum stress should be limited to the range of material parameters from the perspective of mechanical characteristics. It is known that nickel-tungsten (NiW) is widely used in MEMS devices due to its excellent yield strength and other unique properties. NiW can achieve a high yield strength of 1.93 GPa [25], whereas it is only 136 MPa for copper. Thus, it can be expected that the implementation of NiW can enhance the ability of resisting plastic deformation during the process of bending in comparison with copper. The properties of materials utilized in this work are listed in Table 1.

Mechanical Simulation
The solid mechanical module of COMSOL Multiphysics software was utilized to conduct the mechanical simulation of the CMI. The mechanical simulation of the CMI was conducted based on the initially established model. As CMIs are usually applied to the space of fixed height, the displacement boundary condition was set up to be −12.68 μm, as shown in Figure 5a. Additionally, the bottom of the anchor was defined as a fixed constraint, meaning there was no displacement in any direction. Based on the initially established model in Figure 4, the mechanical simulation of the CMI was conducted. The de-

Mechanical Simulation
The solid mechanical module of COMSOL Multiphysics software was utilized to conduct the mechanical simulation of the CMI. The mechanical simulation of the CMI was conducted based on the initially established model. As CMIs are usually applied to the space of fixed height, the displacement boundary condition was set up to be −12.68 µm, as shown in Figure 5a. Additionally, the bottom of the anchor was defined as a fixed constraint, meaning there was no displacement in any direction. Based on the initially established model in Figure 4, the mechanical simulation of the CMI was conducted. The deformed model is shown in Figure 5a, and the vertical displacement is presented in Figure 5b.

Mechanical Simulation
The solid mechanical module of COMSOL Multiphysics software was utilized to conduct the mechanical simulation of the CMI. The mechanical simulation of the CMI was conducted based on the initially established model. As CMIs are usually applied to the space of fixed height, the displacement boundary condition was set up to be −12.68 μm, as shown in Figure 5a. Additionally, the bottom of the anchor was defined as a fixed constraint, meaning there was no displacement in any direction. Based on the initially established model in Figure 4, the mechanical simulation of the CMI was conducted. The deformed model is shown in Figure 5a, and the vertical displacement is presented in Figure  5b. The models of two materials, copper and gold-plated NiW, were simulated and are compared in Figure 6. The maximum stress of the model with copper is about 1.01 GPa, which exceeds the yield strength of copper and can lead to a heavy plastic deformation or even cracks. However, with the implementation of gold-plated NiW, the maximum stress is about 1.66 GPa, which is still within the range of yield strength of NiW. It is worth noting that the oxidation and corrosion would shorten the lifetime of ambient-exposed CMI. To resolve this issue, the method of gold coating proposed in [11] was applied in The models of two materials, copper and gold-plated NiW, were simulated and are compared in Figure 6. The maximum stress of the model with copper is about 1.01 GPa, which exceeds the yield strength of copper and can lead to a heavy plastic deformation or even cracks. However, with the implementation of gold-plated NiW, the maximum stress is about 1.66 GPa, which is still within the range of yield strength of NiW. It is worth noting that the oxidation and corrosion would shorten the lifetime of ambient-exposed CMI. To resolve this issue, the method of gold coating proposed in [11] was applied in this study. Although the thin gold layer, about 0.3 µm in thickness, cannot have an apparent effect in terms of mechanical property, its low resistivity would provide a relatively lower contact resistance than the CMI without the passivation layer. this study. Although the thin gold layer, about 0.3 μm in thickness, cannot have an apparent effect in terms of mechanical property, its low resistivity would provide a relatively lower contact resistance than the CMI without the passivation layer.  Figure 7 shows the contact pressure distribution on the interface of the NiW-made CMI. It is evident that the contact pressure distribution is non-uniform where the largest pressure appears around the edges ➀ and ➂, and the pressure near edge ➃ is the smallest. This is due to the inevitable deformation of the contact plane. As the CMI squeezed, the large deformation in the arm led to the concavity around the edge ➃ and warped deformation near the edges ➀ and ➂. As the material parameters are usually temperature-dependent, temperature would influence the electrical property of the CMI. For example, large contact resistance between CMI and the top pad would induce self-heating, and then the temperature rise in the CMI structure would result in pressure reduction, i.e., increased contact resistance. Consequently, finding the optimal geometric parameters of CMI is of importance. In this study, an intelligent approach was employed to optimize the CMI geometry to achieve minimum contact resistance.  Figure 7 shows the contact pressure distribution on the interface of the NiW-made CMI. It is evident that the contact pressure distribution is non-uniform where the largest pressure appears around the edges 1 and 3 , and the pressure near edge 4 is the smallest. This is due to the inevitable deformation of the contact plane. As the CMI squeezed, the large deformation in the arm led to the concavity around the edge 4 and warped deformation near the edges 1 and 3 . this study. Although the thin gold layer, about 0.3 μm in thickness, cannot have an apparent effect in terms of mechanical property, its low resistivity would provide a relatively lower contact resistance than the CMI without the passivation layer.  Figure 7 shows the contact pressure distribution on the interface of the NiW-made CMI. It is evident that the contact pressure distribution is non-uniform where the largest pressure appears around the edges ➀ and ➂, and the pressure near edge ➃ is the smallest. This is due to the inevitable deformation of the contact plane. As the CMI squeezed, the large deformation in the arm led to the concavity around the edge ➃ and warped deformation near the edges ➀ and ➂.  As the material parameters are usually temperature-dependent, temperature would influence the electrical property of the CMI. For example, large contact resistance between CMI and the top pad would induce self-heating, and then the temperature rise in the CMI structure would result in pressure reduction, i.e., increased contact resistance. Consequently, finding the optimal geometric parameters of CMI is of importance. In this study, an intelligent approach was employed to optimize the CMI geometry to achieve minimum contact resistance. As the material parameters are usually temperature-dependent, temperature would influence the electrical property of the CMI. For example, large contact resistance between CMI and the top pad would induce self-heating, and then the temperature rise in the CMI structure would result in pressure reduction, i.e., increased contact resistance. Consequently, finding the optimal geometric parameters of CMI is of importance. In this study, an intelligent approach was employed to optimize the CMI geometry to achieve minimum contact resistance.

Parameters Optimization and Results
In this section, the modeling of CMI was carried out based on the COMSOL Multiphysics software, and through the application programming interface (API) between COMSOL and MATLAB, the modeling and simulation were parameterized, and programming controlled. As the geometric structure of CMI is determined by multiple parameters, a multi-objective optimization algorithm with the ability of global searching was considered.
As the CMI is usually used in cases where the height of the deformed CMI is fixed, the deformation displacement was set as a constant to ensure a reasonable comparison. Due to the shunt effect of grids' equivalent resistances, a wider contact area is more likely to form a smaller contact resistance. Consequently, the excessive optimized results which led to a small contact area were discarded during the optimization process. Figure 8 shows the convergence curve of the optimization. It can be observed that the convergence was achieved at the 39th iteration step, and in total, 57 iterations took about 9.1 h. The geometric parameters of the initial and optimized CMIs are listed in Table 2. The simulated pressure distributions along the contact boundaries of initial and optimized CMIs are compared in Figure 9. The maximum pressures in the initial and optimized CMIs were 31 MPa and 40.5 MPa, respectively. The improvement in the contact pressure could reduce contact resistance. The contact resistances of initial and optimized CMIs were calculated as 155.3 mΩ to 108.8 mΩ, respectively, i.e., a 29.9% drop in contact resistance was achieved. The improvement in contact resistance can benefit the DC IR drop and reduce the CMI array size to achieve high integration density.
Electronics 2021, 10, x FOR PEER REVIEW mized CMIs are compared in Figure 9. The maximum pressures in the initial an mized CMIs were 31 MPa and 40.5 MPa, respectively. The improvement in the pressure could reduce contact resistance. The contact resistances of initial and op CMIs were calculated as 155.3 mΩ to 108.8 mΩ, respectively, i.e., a 29.9% drop in resistance was achieved. The improvement in contact resistance can benefit the drop and reduce the CMI array size to achieve high integration density.

Multiphysics Coupling Simulation and Results
In this section, the multiphysics coupling simulation was carried out for both original and optimized CMIs. Furthermore, a cell of 3 × 3 CMI array was simulated to evaluate the thermal performance for an array structure. To conduct the multiphysics simulations, the corresponding settings for individual physics and coupling module should be set up, including the discretization of geometric model, setting up boundary conditions, and assigning a coupling path. For electrical simulation, a current source with 1 A amplitude was assigned to the top pad while the bottom surface of the anchor was connected to the ground. The fixed temperature boundary condition with a value of 293.15 K and a heat transfer boundary condition with 7 × 10 6 W/m 2 ·K heat flux were assigned to the bottom surface of the anchor and top pad for thermal simulation, respectively. Rest boundaries

Multiphysics Coupling Simulation and Results
In this section, the multiphysics coupling simulation was carried out for both original and optimized CMIs. Furthermore, a cell of 3 × 3 CMI array was simulated to evaluate the thermal performance for an array structure. To conduct the multiphysics simulations, the corresponding settings for individual physics and coupling module should be set up, including the discretization of geometric model, setting up boundary conditions, and assigning a coupling path. For electrical simulation, a current source with 1 A amplitude was assigned to the top pad while the bottom surface of the anchor was connected to the ground. The fixed temperature boundary condition with a value of 293.15 K and a heat transfer boundary condition with 7 × 10 6 W/m 2 ·K heat flux were assigned to the bottom surface of the anchor and top pad for thermal simulation, respectively. Rest boundaries were designated to be natural convection boundary conditions with 20 W/m 2 ·K heat flux. The setting for the elastic mechanical simulation was kept the same as that in Section 3. For the coupling module, the Joule heat achieved by the electrical simulation was set to be the power source for thermal simulation, and the temperature distribution in the simulated region was updated and coupled back to electrical and mechanical simulations through temperature-dependent parameters. Figure 10 presents the simulated potential distributions at the contact interface of both the initial and optimized CMIs. It can be observed that the potential distributions are consistent with those of pressure (see Figure 9). As the pressure distribution is non-uniform, the contact resistance is varying with the position. When the current flows through the contact interface, the potential difference takes place. As shown in Figure 10, the potential difference of optimized CMI was smaller than that of the initial one, which shows that the optimized CMI has a more uniform distribution of potential. The average DC IR drop of the initial and optimized CMIs were 11.7 mV and 9.3 mV, respectively, indicating that the optimized CMI performed better in power integrity. Figure 11 presents the simulated temperature distributions. The temperature variations were about 36.6 • C and 26.8 • C and the average temperatures are about 48.8 • C and 40.5 • C, respectively. About an 8.3 • C decrease in average temperature was achieved after the optimization. were designated to be natural convection boundary conditions with 20 W/m 2 ·K heat flux. The setting for the elastic mechanical simulation was kept the same as that in Section 3. For the coupling module, the Joule heat achieved by the electrical simulation was set to be the power source for thermal simulation, and the temperature distribution in the simulated region was updated and coupled back to electrical and mechanical simulations through temperature-dependent parameters. Figure 10 presents the simulated potential distributions at the contact interface of both the initial and optimized CMIs. It can be observed that the potential distributions are consistent with those of pressure (see Figure 9). As the pressure distribution is non-uniform, the contact resistance is varying with the position. When the current flows through the contact interface, the potential difference takes place. As shown in Figure 10, the potential difference of optimized CMI was smaller than that of the initial one, which shows that the optimized CMI has a more uniform distribution of potential. The average DC IR drop of the initial and optimized CMIs were 11.7 mV and 9.3 mV, respectively, indicating that the optimized CMI performed better in power integrity. Figure 11 presents the simulated temperature distributions. The temperature variations were about 36.6 °C and 26.8 °C and the average temperatures are about 48.8 °C and 40.5 °C, respectively. About an 8.3 °C decrease in average temperature was achieved after the optimization.  Due to the large contact resistance at the interface, a single CMI unit cannot satisfy the requirement of power integrity, and an array of CMIs should be employed to provide a wider channel for current and reliable ability of mechanical supporting. As such, performance of the 3 × 3 CMI array shown in Figure 12 was simulated and analyzed. The CMIs were arranged in an array with a gap of 250 μm × 250 μm, lying in the center of the IC chip. Except for the fixed distance between the silicon interposer and IC chip, other settings of mechanical simulation were the same as those utilized in CMI unit simulation. In the electrical simulation, a current source with the amplitude of 1A was applied to the bottom surface of the IC and the anchor of CMI was connected to the ground. Thermal Due to the large contact resistance at the interface, a single CMI unit cannot satisfy the requirement of power integrity, and an array of CMIs should be employed to provide a wider channel for current and reliable ability of mechanical supporting. As such, performance of the 3 × 3 CMI array shown in Figure 12 was simulated and analyzed. The CMIs were arranged in an array with a gap of 250 µm × 250 µm, lying in the center of the IC chip. Except for the fixed distance between the silicon interposer and IC chip, other settings of mechanical simulation were the same as those utilized in CMI unit simulation. In the electrical simulation, a current source with the amplitude of 1A was applied to the bottom surface of the IC and the anchor of CMI was connected to the ground. Thermal simulation was carried out with fixed temperature boundary conditions (T 0 = 20 • C) and air convection conditions (15 W/m 2 ·K) assigned to the bottom of the interposer and top of the IC, respectively. Considering the limited computing resource and computational efficiency, the details of the internals of the IC chip were not modeled and an isopyknic cube with equivalent material parameters was utilized. In order to include the power dissipation of the IC chip in our simulation, a heat source with 3 × 10 9 W/m 3 power density was applied to the IC region.
Electronics 2021, 10, x FOR PEER REVIEW 1 simulation was carried out with fixed temperature boundary conditions (T0 = 20 °C air convection conditions (15 W/m 2 ·K) assigned to the bottom of the interposer and the IC, respectively. Considering the limited computing resource and computationa ciency, the details of the internals of the IC chip were not modeled and an isopyknic with equivalent material parameters was utilized. In order to include the power di tion of the IC chip in our simulation, a heat source with 3 × 10 9 W/m 3 power densit applied to the IC region.

IC chip
Interposer/Substrate CMI Array Cell Figure 12. The schematic of a 3 × 3 CMI array.
The simulated temperature distribution of the CMI array is presented in Figure  is evident that the average temperature of the CMI in the array was 38 °C which is than the 40.5 °C achieved in the simulation of a single CMI unit. This is mainly be the current density flowing through the CMI unit in the array is only 1/9 A. As the mechanical parameters of the materials are usually temperature-depen the temperature variation in the CMI would lead to a change in its mechanical pro and further influence the electric property. The Young's modulus, a paramete measures the non-deformability of materials, decreases when temperature rises. Fro microscopic point of view, this is because the thermal movement of atoms is inten and the bonding force between atoms is weakened when the temperature rises. In The simulated temperature distribution of the CMI array is presented in Figure 13. It is evident that the average temperature of the CMI in the array was 38 • C which is lower than the 40.5 • C achieved in the simulation of a single CMI unit. This is mainly because the current density flowing through the CMI unit in the array is only 1/9 A. simulation was carried out with fixed temperature boundary conditions (T0 = 20 °C) and air convection conditions (15 W/m 2 ·K) assigned to the bottom of the interposer and top of the IC, respectively. Considering the limited computing resource and computational efficiency, the details of the internals of the IC chip were not modeled and an isopyknic cube with equivalent material parameters was utilized. In order to include the power dissipation of the IC chip in our simulation, a heat source with 3 × 10 9 W/m 3 power density was applied to the IC region.

IC chip
Interposer/Substrate CMI Array Cell Figure 12. The schematic of a 3 × 3 CMI array.
The simulated temperature distribution of the CMI array is presented in Figure 13. It is evident that the average temperature of the CMI in the array was 38 °C which is lower than the 40.5 °C achieved in the simulation of a single CMI unit. This is mainly because the current density flowing through the CMI unit in the array is only 1/9 A. As the mechanical parameters of the materials are usually temperature-dependent, the temperature variation in the CMI would lead to a change in its mechanical property, and further influence the electric property. The Young's modulus, a parameter that measures the non-deformability of materials, decreases when temperature rises. From the microscopic point of view, this is because the thermal movement of atoms is intensified and the bonding force between atoms is weakened when the temperature rises. In other words, the material is softer than before and CMI has a tendency for further bending, as As the mechanical parameters of the materials are usually temperature-dependent, the temperature variation in the CMI would lead to a change in its mechanical property, and further influence the electric property. The Young's modulus, a parameter that measures the non-deformability of materials, decreases when temperature rises. From the microscopic point of view, this is because the thermal movement of atoms is intensified and the bonding force between atoms is weakened when the temperature rises. In other words, the material is softer than before and CMI has a tendency for further bending, as shown in Figure 14a. Transient multiphysics coupling simulation was carried out to capture transient electrothermal effect on the performance of the optimized CMI. In the simulation, a step current pulse with 1A amplitude and 100 µs rising time was applied to the top surface of the metal pad. The simulated time-varying contact resistance is presented in Figure 14b, where the contact resistance increases from the initial value (108.8 mΩ) to a steady value (109.27 mΩ) in 1ms. Table 3 shows the stress, contact pressure, and contact resistance of cases with and without the consideration of temperature influence in steady multiphysics simulations. It is evident that the contact pressure at the interface becomes smaller when the temperature effect is included in the simulation, and the drop in pressure leads to a slight increase in the contact resistance.

Conclusions
In this study, a model of CMI was proposed and optimized based on the coupled framework of the PSO approach and multiphysics simulations for the first time. A contact resistance model, which is the function of contact pressure, was proposed to quantitatively analyze the electrical performance of the CMI. Based on the simulated data, the intelligent optimization algorithm, PSO, was employed to conduct the optimization. A simulation environment based on MATLAB and COMSOL was set up to support the optimization process. After the optimization, a 29.9% reduction in the contact resistance and 2.4 mV decrease in DC IR drop were achieved. Multiphysics simulations were conducted for both initial and optimized CMIs, and an 8.3 °C reduction in temperature was achieved. Finally, the influence of the self-heating effect on the mechanical property and contact resistance was studied.

Conclusions
In this study, a model of CMI was proposed and optimized based on the coupled framework of the PSO approach and multiphysics simulations for the first time. A contact resistance model, which is the function of contact pressure, was proposed to quantitatively analyze the electrical performance of the CMI. Based on the simulated data, the intelligent optimization algorithm, PSO, was employed to conduct the optimization. A simulation environment based on MATLAB and COMSOL was set up to support the optimization process. After the optimization, a 29.9% reduction in the contact resistance and 2.4 mV decrease in DC IR drop were achieved. Multiphysics simulations were conducted for both initial and optimized CMIs, and an 8.3 • C reduction in temperature was achieved. Finally, the influence of the self-heating effect on the mechanical property and contact resistance was studied.

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