Thermal simulations for integrated circuits

Miniaturizing components has always been a primary goal in the semiconductor industry because it cuts production cost and lets companies build smaller computers and other devices. Miniaturization, however, has increased dissipated power per unit area and made it a key limiting factor in integrated circuit performance. Temperature increase becomes relevant for relatively small-cross-sections wires, where it may affect normal semiconductor behavior. Besides, since the generation of heat is proportional to the frequency of operation for switching circuits, fast computers have larger heat generation than slow ones, an undesired effect for chips manufacturers. This article summaries physical concepts that describe the generation and conduction of heat in an integrated circuit, and presents numerical methods that model heat transfer from a macroscopic point of view.
Generation and transfer of heat.
Fourier's law.
At macroscopic level, Fourier's law states a relation between the transmitted heat per unit time per unit area and the gradient of temperature:
Where formula_2 is the thermal conductivity, [W·m−1 K−1].
Joule heating.
Electronic systems work based on current and voltage signals. Current is the flow of charged particles through the material and these particles (electrons or holes), interact with the lattice of the crystal losing its energy which is released in form of heat. Joule Heating is a predominant mechanism for heat generation in integrated circuits and is an undesired effect in most of the cases. For an ohmic material, it has the form:
Where formula_4 is the current density in [A·m−2], formula_5 is the specific electric resistivity in [formula_6·m] and formula_7 is the generated heat per unit volume in [W·m−3].
Heat-transfer equation.
The governing equation of the physics of the heat transfer problem relates the flux of heat in space, its variation in time and the generation of power by the following expression:
Where formula_2 is the thermal conductivity, formula_10 is the density of the medium, formula_11 is the specific heat, formula_12, the thermal diffusivity and formula_13 is the rate of heat generation per unit volume. Heat diffuses from the source following the above equation and solution in an homogeneous medium follows a Gaussian distribution.
Techniques to solve heat equation.
Kirchhoff transformation.
To get rid of the temperature dependence of formula_14, Kirchhoff transformation can be performed 
where formula_16 and formula_17 is the heat sink temperature. When applying this transformation, the heat equation becomes:
where formula_19 is called the diffusivity, which also depends on the temperature. To completely linearize the equation, a second transformation is employed:
formula_20
yielding the expression:
formula_21
Simple, direct application of this equation requires approximation. Additional terms arising in the transformed Laplacian are dropped, leaving the Laplacian in its conventional form.
Analytical solutions.
Although analytical solutions can only be found for specific and simple cases, they give a good insight to deal with more complex situations. Analytical solutions for regular subsystems can also be combined to provide detailed descriptions of complex structures. In Prof. Batty's work, a Fourier series expansion to the temperature in the Laplace domain is introduced to find the solution to the linearized heat equation.
Example.
This procedure can be applied to a simple but nontrivial case: an homogeneous cube die made out of GaAs, L=300 um. The goal is to find the temperature distribution on the top surface. The top surface is discretized into smaller squares with index i=1...N. One of them is considered to be the source.
Taking the Laplace transform to the heat equation:
where formula_23
Function formula_24 is expanded in terms of cosine functions for the formula_25 and formula_26 variables and in terms of hyperbolic cosines and sines for formula_27 variable. Next, by applying adiabatic boundary conditions at the lateral walls and fix temperature at the bottom (heat sink temperature), thermal impedance matrix equation is derived:
Where the index formula_29 accounts for the power sources, while the index formula_30 refers to each small area.
For more details about the derivation, please see Prof. Batty's paper.
The below figure shows the steady state temperature distribution of this analytical method for a cubic die, with dimensions 300 um. A constant power source of 0.3W is applied over a central surface of dimension 0.1L x 0.1L. As expected, the distribution decays as it approaches to the boundaries, its maximum is located at the center and almost reaches 400K
Numerical solutions.
Numerical solutions use a mesh of the structure to perform the simulation. The most popular methods are: Finite difference time-domain (FDTD) method, Finite element method (FEM) and method of moments (MoM).
The finite-difference time-domain (FDTD) method is a robust and popular technique that consists in solving differential equations numerically as well as certain boundary conditions defined by the problem. This is done by discretizing the space and time, and using finite differencing formulas, thus the partial differential equations that describe the physics of the problem can be solved numerically by computer programs.
The FEM is also a numerical scheme employed to solve engineering and mathematical problems described by differential equations as well as boundary conditions. It discretizes the space into smaller elements for which basis functions are assigned to their nodes or edges. Basis functions are linear or higher order polynomials. Applying the differential equation and the boundary conditions of the problem to the basis functions, a system of equations is formulated using either the Ritz or Galerkin method. Finally, a direct or iterative method is employed to solve the system of linear equations. For the thermal case, FEM method is more suitable due to the nonlinearity nature of the thermal properties.
Example.
The previous example can be solved with a numerical method. For this case, the cube can by discretized into rectangular elements. Its basis functions can be chosen to be a first order approximation (linear):
formula_31
formula_32
formula_33
where formula_34. If formula_35, then formula_36.
Using this basis functions and after applying Galerkin's method to the heat transfer equation, a matrix equation is obtained:
formula_37
where,
This expressions can be evaluated by using a simple FEM code. For more details, please see. The figure below shows the temperature distribution for the numerical solution case. This solution shows very good agreement with the analytical case, its peak also reaches 390 K at the center. The apparent lack of smoothness of the distribution comes from the first order approximation of the basis functions and this can be solved by using higher order basis functions. Also, better results might be obtained by employing a denser mesh of the structure; however, for very dense meshes the computation time increases a lot, making the simulation non-practical.
The next figure shows a comparison of the peak temperature as a function of time for both methods. The system reaches steady state in approximately formula_41.
Model order reduction.
The numerical methods such as FEM or FDM derive a matrix equation as shown in the previous section. To solve this equation faster, a method called Model order reduction can be employed to find an approximation of lower order. This method is based on the fact that a high-dimensional state vector belongs to a low-dimensional subspace .
Figure below shows the concept of the MOR approximation: finding matrix V, the dimension of the system can be reduced to solve a simplified system. 
Therefore, the original system of equation:
becomes:
Whose order is much lower than the original making the computation much less expensive. Once the solution is obtained, the original vector is found by taking the product with V.
Conclusion.
The generation of heat is mainly produced by joule heating, this undesired effect has limited the performance of integrated circuits. In the preset article heat conduction was described and analytical and numerical methods to solve a heat transfer problem were presented. Using these methods, steady state temperature distribution was computed as well as the peak temperature as a function of time for a cubic die. For an input power of formula_44 (or formula_45) applied over a single surface source on the top of a cubic die a peak increment of temperature in the order of 100 K was computed. Such increase in temperature can affect the behavior of surrounding semiconductor devices. Important parameters like mobility change drastically. That is why the heat dissipation is a relevant issue and must be considered for circuit designing.