An Exact Recursion Relation Solution for the Steady-State Surface Temperature of a General Multilayer Structure

A recursion relation technique has been used in the past to determine the surface potential from the multilayer electrical Laplace equation. This has provided for a vastly simplified evaluation of the electrical spreading resistance and four-probe resistance. The isomorphism of the multilayer Laplace equation and the multilayer steady-state heat flow equation suggests the possibility of developing a recursion relation applicable to the multilayer thermal problem. This recursive technique is developed and is shown to provide the surface temperature of the multilayer steady-state beat flow equation. For the three-layer case, the thermal recursion relation readily yields the surface results, which are identical with those presented by Kokkas and the TXYZ thermal code. This recursive technique can be used with any number of layers while incurring only a small increase in computation time for each added layer. For the case of complete, uniform top surface coverage by a heat source, the technique gives rise to the generalized one-dimensional thermal resistance result. An example of the use of the new recursive method is provided by the preliminary calculations of the surface temperature of a buried oxide (SOI, SIMOX) structure containing several thicknesses of the surface silicon layers. This new technique should prove useful in the investigation and understanding of the steady-state thermal response of modem multilayer microelectronic structures.

I on the thermal analysis of multiple-layer structures [l], [2].This is based on the solution of the heat flow equation for a multilayer rectangular structure.While the interface boundary conditions for the continuity of the temperature and the heat flow are for a general multilayer structures, the equations have been solved in detail and evaluated for the case of up to a three-layer structure.The results of the three-layer steady-state heat flow analysis are exact for this situation.This three-layer solution has been numerically implemented in the TXYZ program [3], [4].This model and the code have found wide use in a variety of steady-state problems [5].In addition, the code has served as a benchmark for checking the numerical accuracy and consistency of finite-difference and finite-element programs [6].
Manuscript received April 6, 1994; revised September 1, 1994.This work was supported by the National Institute of Standards and Technology's Advanced Technology Program in support of the American Scaled-Electronics Consortium.This paper was presented at the 10th Semiconductor Thermal Management and Measurement Symposium, San Jose, CA, February 1-3, 1994.
The author is with the Semiconductor Electronics Division, National Institute of Standards and Technology, Gaithersburg, MD 20899 USA.
IEEE Log Number 9407173.
During this same twenty year time span, there have been a number of advances which have taken place in the theoretical analysis of the analogous electrical problem of two-probe spreading resistance.The starting point of the work on the electrical problem is the multilayer Laplace equation analysis of Schumann and Gardner [7], [SI.Originally, the numerical analysis required a mainframe computer to do the rather cumbersome matrix algebra which develops quickly once the calculations get past two or three layers.This limitation coupled with difficulties involved in the evaluation of the accompanying oscillatory integrals served to severly limit the applicability of the technique.However, a major advance came with the introduction of a recursion relation by Choo [9] for the construction of the solution of the N-layer problem from the solution of the ( N -1)-layer problem.The recursion relation had previously been introduced in geophysical literature [lo], [ 1 I].The book by Koefoed [ 101 provides a comprehensive discussion of Laplace's equation and its solution, the multilayer approach applied to geoelectric resistance measurements, and the derivation of the recursion relation.It should be noted that the mathematical description of spreading resistance and geoelectric resistance are exactly the same even though the length scales are vastly different.
The introduction of the recursion relation eliminated the cumbersome matrix manipulations and allowed for the relatively straightforward algebraic evaluation of the functions required for the multilayer analysis.The form of the recursion relation has been clarified by Berkowitz and Lux [ 121 and has provided one of the cornerstones of the routine calculation of the spreading resistance from the resistivity structure as well as the extraction of the resistivity structure from the spreading resistance.The final cornerstone was introduced by inventive quadrature techniques which allow for the rapid evaluation of the necessary integrals.Calculations which would have been prohibitive, if not impossible, twenty years ago can now be performed routinely on PC-based systems.The model, the numerical analysis, and the FORTRAN codes have been summarized recently [ 131.As the electrical problem and the thermal problem are similar from the mathematical viewpoint, the question arose as to the applicability of the recursion relation technique, with possible and appropriate modification, to the steadystate thermal problem.If this could be accomplished, then it would provide for an exact solution of the steady-state surface temperature for a multilayer structure with any number U.S. Government work not protected by U.S. Copyright of layers.The present work describes the modified recursion relation and its application to the thermal problem.The salient features of the recursion relation are presented along with a comparison of its results with those of the Kokkas three-layer solution.The use of the recursion relation for extension of the calculations beyond three layers is discussed with emphasis on the small increase in calculation time for each added layer.The computational efficiency of the recursion relation is also discussed.This new technique should prove extremely useful in the understanding of the thermal effects of multilayer structures.
As the similarity to the spreading resistance problem is part of the motivation, a brief discussion of the analysis involved in this problem is presented.This is followed by the corresponding analysis for the multilayer thermal problem.The connection of these two problems is then made along with a discussion of the need for a modified recursion relation for the thermal case.Its form is then determined and the thermal recursion relation is presented.Finally, the implications of this thermal recursion relation are discussed.

THE N-LAYER ELECTRICAL PROBLEM
For the case of the spreading resistance problem, the multilayer solution of the Laplace equation provides for the calculation of the resistance between the contacts on the top surface of a nonuniform conductivity (resistivity) structure.The fundamental assumption of the multilayer analysis is that Laplace's equation is satisfied in each 'layer' in the material.The problem is set up in cylindrical coordinates to emulate the symmetry of a single circular contact on the top surface of the material.The potential is assumed to be independent of the angular variable in this coordinate system.The Laplace equation may then be written as where V ( T , 2 ) is the potential, T is the radial coordinate, and z is the depth coordinate.This equation may be solved by means of separation of variables, with the result that a particular solution is where Ja(Xr) is the Bessel function and X2 is the separation of variables constant.The boundary condition on the T part of the solution is that V ( T , Z ) approaches zero as T tends to infinity.This is satisfied by the above for all values of A. Then, the general solution is an integral of the particular solution with a weighting factor and is of the form where the weighting functions B ( X ) and $(A) are determined from the z-dependent boundary conditions.The above is the general solution of Laplace's equation in cylindrical coordinates for a single layer.This provides the framework in which the multilayer solution may be addressed.
The depth-dependent portion of the multilayer geometry is presented in Fig. 1.Each layer is described by a thickness and electrical conductivity.Charge neutrality is assumed in each layer so that the potential in each layer satisfies Laplace's equation.The general one-layer solution given by (3) provides a convenient and useful basis for the required layer solution.Then, the solution of Laplace's equation in the ith layer in an N-layer structure may be written as The boundary conditions used to solve the system of equations (Le., determine {&(A)} and {?/12(A)}:i = l ; . .! N ) are provided by conditions on the top surface, the intermediate interfaces, and the bottom surface.On the top surface, current flow takes place only through the probe which is modeled as a circular plate of radius a. Then the top surface boundary condition is expressed as where is the electrical conductivity of the top layer, .J(T) is the probe current density [7], [8], and the derivative is evaluated at z = 0. On the bottom surface, the potential is assumed to be well behaved and, more specifically, is assumed to approach zero.This is usually expressed as Equations ( 5) and (6) provide 2 of the 2 N conditions required to solve the N-layer system of equations.The remaining 2(N -1) conditions are provided by requirements at the interfaces between the layers where the potentials and the current densities are assumed to be continuous [7], [8].These are expressed as and where the functions and their derivatives are to be evaluated at the interfacial boundaries.
For the case of an N-layer structure, the substitution of (4) into the boundary conditions given by ( 5)-( 8) gives rise to a set of 2 N equations in 2N unknowns ({O,(A)}, { l l j z ( X ) } , i = 1,. . . .N ) .The analytic solution of this system of equations requires the use of Cramer's rule of linear algebra.Clearly, this can become rather tedious especially since the expansion coefficients are functions of the continuous variable A. It is possible to show that the potential on the top surface of the N-layer structure may be written as [7], [SI where the kernel function A,v(X) = 1 + 20Av(X) depends upon the electrical conductivities and thicknesses of the 'layers' in the multilayer structure (through the solution of the above system of simultaneous equations).Schumann and Gardner were able to work out the system of equations for the cases up to I: = 3. Cases beyond three-layers presented a major stumbling block.
The development of the recursion relation by Koefoed [ 101 and its use by Choo et al. [9] provided a major breakthrough in the evaluation of (9) and effectively removed the numerical difficulty associated with matrix inversion.The philosophy and utility of a recursion relation is that the kernel function, AN(A), for an N-layer structure can be easily generated from the kernel of an ( N -1)-layer structure by means of an algebraic relation.This algebraic relation represents a reduction of the matrix algebra and the subsequent manipulations with determinants.The reader may wish to refer to Koefoed's book for a detailed discussion.In particular, if the kernel is known for the ( N -1)-layer case, then the kernel for the N-layer case is given by where 1 -exp( -2Ad,rq) and d,v is the thickness of layer N .In practice, the recursion relation is begun with the one-layer case from which the two-layer kernel is generated.Then the three-layer kemel is generated from the two-layer kernel.This sequence is repeated until the N-layer kemel is determined.Notice how the conductivity and thickness of the Nth layer enter in the recursion relation.Berkowitz and Lux [ 121 have shown that the recursion relation for the kernel function could be expressed as The recursion relation is most commonly used in this form.This is due, in part, to the hyperbolic tangent function which is bounded between (0, I ) making the numerical implementation relatively easy.This coupled with the quadrature techniques introduced by Berkowitz and Lux has moved the analysis into everyday applicability.

THE N-LAYER THERMAL PROBLEM
The previous section chronicles a problem and its inventive solution for the multilayer Laplace equation.In the present section, the multilayer steady-state heat flow equation is presented and discussed in a manner which points directly to the introduction of a recursion relation to solve analogous difficulties.Fig. 2 contains the basic geometric description of the multilayer thermal problem.The following analysis begins with the single-layer situation.The steady-state heat flow equation in Cartesian coordinates for a single layer is This problem is solved using separation of variables and Fourier expansion techniques.In this analysis, there is no heat transport by either convection or radiation; conduction is the only mechanism considered.The lateral boundary condition imposed upon the problem is that there is no heat transport out of lateral boundaries, i.e., Equation ( 15) is the general solution of the steady-state heat flow equation in Cartesian coordinates for a single layer.This provides the framework in which the multilayer problem may be addressed.For the case of a multiple-layer structure, the depth-dependent portion of the problem is pictured in Fig. 3.The temperature in each of the layers is assumed to satisfy the steady-state heat flow equation.The general one-layer solution provides a convenient and useful basis for the multilayer solution.Then, the Fourier coefficients for the solution of the steady-state heat flow equation in the ith layer of an N-layer structure may be written as where the expansion coefficients, (Y; and [I;, are determined from the r-dependent boundary conditions.Using the index notation similar to that of the electrical problem, i.e., the bottom layer has the index 1 and the top layer has the index N , the boundary conditions used to solve the set of equations in (1 8) may be expressed as follows.First, heat enters (or leaves) the structure through the portions of the top layer where power is applied.This is expressed as where the derivative is evaluated at z = 0, and where P ( x , y) is the power function.This is expressed as P ( x , y))= PoU(x, ly), where PO is the power density (watts/cm2), and V(x,y) describes the surface geometrical distribution of the power sources.Next, the bottom layer boundary condition is provided by the requirement that the temperature is continuous across the interface between the bottom layer and the heat sink and equal to the ambient which is taken as zero, where the function is evaluated at the interface.Finally, the remaining 2( N -1) conditions are provided by the requirement that the temperature and the normal component of the heat flow are continuous across the internal interfaces.These are expressed as where the functions and their derivatives are to be evaluated at the interfacial boundaries.
For the case of an N-layer structure, the substitution of (15) and (18) into the boundary conditions given by (19), (20), (21), and (22) gives rise to a set of 2 N equations in 2 N unknowns ({crZ(-y)}, { / j Z ( y ) } .i = 1, ..., N).The analytic solutionofthis system of equations, as the system of the multilayer Laplace equations, requires the use of the same kind of matrix algebra.Clearly, this can become rather tedious, especially since the expansion coefficients are functions of the continuous variable, y.It is possible to show that the temperature on the top surface of the N-layer structure may be written as [2], [3]   o i l 7 3 where U ( n , m) is Fourier cosine transform of the geometric portion of the input power function, The key elements in (23) are the Fourier coefficients TN (y) which depend upon the thermal conductivities and thicknesses of the layers in the multilayer structure (through the solution of the above system of simultaneous equations).
For cases up to 1: = 3, Kokkas [l], [2] was able to work out the system of equations for not only the surface temperature, but also the temperature inside the three-layer structure.The numerical implementation of this work is contained in the TXYZ code [3], [4].However, cases beyond three layers present a major difficulty.The analogous difficulty with the solution of the multilayer Laplace equation and its removal by means of the recursion relation provides the impetus for the introduction of a recursive scheme in the solution of multilayer thermal equation.
It is important to note that the factor of y has been explicitly written in the denominator of (23).This has been done to simplify the discussion of the recursion relation described in the next section.This recursion relation will be used to take into account the multilayer thicknesses and thermal conductivities.

IV. THE THERMAL RECURSION RELATION
The strong resemblances of (5) and (19), (7) and (21), and (8) and ( 22) lead to the possibility of expressing the Fourier expansion coefficient for the surface temperature of an N-layer structure in the form of (12) with the appropriate transcription from electrical conductivities to thermal conductivities, X to 7, A,v(X) to ~.v(y), and d~ to L N .Then the thermal recursion relation would take the form . (25) 6 ~r -1 As the electrical and thermal bottom layer boundary conditions are of different form, it is necessary to evaluate the corresponding T~ (7).The form of TI (y) may be obtained from (28) [2] or from (89) [ 3 ] .In either case, this may be done by setting z = 0 and then setting the thickness of the two top layers equal to zero.The result is that q(y) = tanh (yL1). (26) Equations ( 25) and ( 26) are the central results of the determination of the thermal recursion relation.They provide for the determination of the Fourier coefficients of the surface temperature for an N-layer structure from the set of thicknesses and thermal conductivities (Li, ~i , 1; = 1, N) by repeated use of the thermal recursion relation.This process begins with the one-layer Fourier coefficient as given by (26).This is substituted into (25) to generate the two-layer Fourier coefficient.This is substituted into (25) to determine the three-layer Fourier coefficient.This process is repeated until the Fourier coefficient for the desired number of layers is generated.Note that there is no restriction as to the number of layers.This means that the exact steady-state surface temperature can now be determined for a multilayer structure with an arbitrary number of layers.It is important to note that this exact steady-state surface temperature satisfies the heat flow equation as well as the necessary boundary conditions.It does not need verification with existing layered solutions.
The recursion relation technique should complement and possibly extend some of the other efforts [14]-[18] aimed at problems beyond three layers.

v. RELATION TO PREVIOUS SOLUTIONS
It is instructive to compare and connect the results of the thermal recursion relation with those presented in the literature.The principal one is based upon the case of a three-layer structure where the result of the thermal recursion relation is compared with that obtained from the Kokkas model and the TXYZ code.Using (25) and (26), it is straightforward to show that the Fourier coefficient for the three-layer structure is given by ( 27), shown at the bottom of the page.It is also straightforward to show that this is in exact formal agreement with the Kokkas and TXYZ expressions for z = 0 for the three-layer case.This does not serve as a limited verification of the recursion relation but rather shows that the two are in agreement, as expected.In addition to this exact formal agreement, numerical simulations have been performed with the TXYZ code and a preliminary code based upon the recursion relation.As expected, these are in numerical agreement.One interesting result is that the recursion relation code actually runs about 15% faster.This is due to the more compact nature of the recursion relation.
It is interesting to note in (27 is independent of the thermal conductivities and reduces to tanh(y(Ll + L2 + 153)) which is the result for a single layer of thickness equal to L1 + L2 + LJ.This line of reasoning can be extended to the N-layer equation.For the case where the thermal conductivities of the layers become equal, repeated use of the addition formula for the hyperbolic tangent shows that the form of the Fourier coefficient becomes tanh(y(L1 + . . .+ L 2 v ) ) as expected for a single layer of thickness equal to L1 + . . .+ L N .
The numerical implementation of the thermal recursion relation has been used in calculations of the surface temperature for structures of up to seven layers.The case of seven layers is used to provide a convenient range and does not represent a li mitation of the recursion relation.Any number of layers is possible and can be easily achieved on modem computing systems.These temperature calculations were performed for the purpose of testing the computational efficiency of the recursive scheme.The results are presented in Fig. 4 and were obtained on a VAX 11/785.They are shown for the purpose of illustration of the utility of the recursion relation.There it can be seen that there is only about a 5% increase in computation time for each layer added to the structure.This is strong evidence for the effectiveness of the algebraic nature of the recursion relation.
These test cases were also run for the situation where all the thermal conductivities were equal.In keeping with the above discussion, the results reduced to those for a single layer of thickness equal to L1+. ..+Ln; .This provided further support of the recursion relation and its numerical implementation.
The special case of uniform surface coverage by a single heat source for a three-layer structure has been presented in [3].There it was shown that only the y ---) 0 ( n = 0, rrL = 0) term needs to be considered.This same line of reasoning applies to the full N-layer case.Investigation of the recursion relation for small values of y shows for the case of uniform surface coverage by a single heat source that which is just the generalized one-dimensional thermal resistance result.

VI. APPLICATION TO BURIED OXIDE STRUCTURE
The thermal recursion relation should be especially useful in the understanding of the effects of multilayer structure thickness and thermal conductivities on the surface temperature.The possibility of using the surface temperature as a probe of the thermal conductivities of the structure is particularly intriguing.An example of the use of this technique is found in the case of the steady-state surface temperature of a buried oxide structure.Recent trends in semiconductor fabrication make use of the buried oxide as a way of electrically isolating the top active region.Probably the most popular version of this is found in SIMOX (Separation by IMplanted Oxygen) structures where the buried oxide is formed by high-dose implantation of oxygen followed by hightemperature annealing [19].The depth of the buried oxide is controlled by the energy of the incident oxygen ion beam.When processed correctly, the implantation damage should be removed and a buried layer of Si02 should be formed.Electrical isolation is certainly achieved, but the thermal conductivity of the buried oxide may impose a price of thermal isolation which may give rise to higher-than-desired operating temperatures in the device.These kinds of structures are prime candidates for the investigation of self-heating effects in multiple-layered structures [20], [21].
As an example of the recursion relation technique, preliminary calculations were performed on a series of buried oxide structures.These were done to ascertain the possibility of using surface temperature measurements as a probe of the thermal conductivity of the buried oxide.The basic structure was modeled as a chip (rectangular structure) with lateral dimensions of I .0275cm by 0.7737 cm.The vertical structure consisted of a surface layer of Si over a buried layer of Si02 over a substrate layer of Si.The buried oxide was 0.00001 cm thick and the substrate Si layer was 0.03 189 cm thick.Three thicknesses of the surface Si layer were used in the calculations.These included 0.00001 cm, 0.00002 cm, and 0.00003 cm, and are typical of SIMOX structures.The thermal conductivities of the surface and substrate Si layers were taken as the room temperature value of 1.55 and the thermal conductivities of the buried oxide were taken in the range from 0.001 to 0.015.A small resistive heat source on the top surface was modeled as a rectangular element of lateral dimension of 0.001005 cm by 0.0825 cm.The temperature (normalized by the power density, TIP,)) was calculated in the area of the heat source for the ranges of surface Si thickness and buried oxide thermal conductivities.This was done to ascertain: 1) the effects of the surface silicon thickness; and 2) the sensitivity of the calculated temperature on the buried oxide thermal conductivity.The results are presented in Fig. 5, where the calculated surface temperatures (normalized by PO) in the area of the small surface heat source are plotted against the inverse thermal conductivity of the buried oxide for the three surface silicon thicknesses.These results indicate the important dependence on the thickness of the surface silicon layer.In addition, the calculated temperature decreases with increasing thickness of the surface silicon layer.If the heat flow were one-dimensional and described by (28), then the calculated temperature should increase with increasing thickness of the surface silicon layer.The fact that the trend is in the opposite direction points to the two-dimensional heat flow caused by the more thermally conductive surface silicon layer.This necessitates careful interpretation of the temperature data on these structures.

VII. SUMMARY AND CONCLUSIONS
Using the strong mathematical similarity of the multilayer Laplace equation of spreading resistance and the multilayer steady-state heat flow equation, a recursion relation has been developed for the multilayer thermal problem.There are no restrictions as to the number of layers which can be used in the calculations.This makes the results truly multilayer.This recursion relation has been shown to produce the Kokkas-TXYZ results for the three-layer case.It also gives rise to a generalized one-dimensional thermal resistance result for the case of uniform surface coverage.A preliminary program based on this recursive scheme has been shown to provide very good computational speed.As an example, it has been applied to investigate the possible determination of the buried oxide thermal conductivity in SIMOX-type structures.This recursion relation technique is simple, elegant, and powerful and should be extremely useful in the investigation and understanding of the effects of layer thicknesses and thermal conductivities on the steady-state surface temperatures of multilayer structures.

ACKNOWLEDGMENT
The author wishes to thank Harry L. Berkowitz  His research was involved with statistical thermodynamics, the theory of fluids, quantum transport, pressure broadening of spectral lines, and nuclear spin relaxation.Since 1974, he has been a physicist in the Semiconductor Electronics Division at the National Institute of Standards and Technology.His original research was in the field of vibration analysis of ultrasonic bonding tools used in microelectronics.Over the years, he has been extensively involved in the modeling and analysis of electrical spreading resistance and four-probe resistance in the determination of carrier profiles in semiconductor device structures.His work helped to elucidate the need for the use of the Poisson equation for the phenomenon of carrier redistribution in shallow device structures.In the early 1980's.he became involved in the thermal analysis of semiconductor devices.His work in thermal modeling led to the TXYZ code, which has received wide distribution in both industry and academia.At this same time, he began work in the field of semiconductor process modeling, in general, and range and damage calculations for ion implantation, in particular.The ion implant calculations were instrumental in providing a sound basis for an ion implant standard needed for the Secondary Ion Mass Spectrometry atomic profiling technique.In 1990, seeing the need for a forum for the discussion of issues and exchange of information on ion implantation, he founded the Ion Implant Users Group.He has a continuing interest in atomic and electrical profiling, process modeling, ion implantation calculations, and the Ion Implant Users Group.He is actively involved in the investigation of the thermal conductivity of thin Si02 layers, the application of the multilayer thermal code to microelectronic structures, and the steady-state and transient thermal response of power device and other semiconductor structures.He has published scores of archival papers, reports and conference talks on the above topics.He is also responsible for the TXYZ and RESPAC software packages.

Fig. 1 .
Fig. 1.Schematic representation of the geometry used in the Schumann and Gardner multilayer Laplace equation analysis.In this figure, d , , U , , and I , are the thickness, electrical conductivity, and potential, respectively, in the rth layer.

Fig. 2 .Fig. 3 .
Fig. 2. Geometry of the multilayer structure for which the steady-state temperature is calculated.

Fig. 4 .
Fig. 4. Example of the small increase in computation time brought about by increasing the number of layers in the thermal structure.The slope represents about a 5% per layer increase.

Fig. 5 .
Fig. 5 .Example of the application of the thermal recursion relation to the calculation of the surface temperature of an SO1 (Silicon On Insulator) structure.The insert (upper left) depicts the structure where the layer thicknesses and heater dimensions are not to scale.The buried Si02 is 0.00001 cm thick.Curves A, B, and C are for the cases where the surface layers of Si are O.oooO1 cm, O.ooOo2 cm, and O.ooOo3 cm thick, respectively.Notice how the temperatures do not follow a one-dimensional interpretation in keeping with the two-dimensional heat flow in the more thermally conductive surface silicon layer.
previously of Fort Monmouth and Solid State Measurements for past collaborations which ultimately led to the present work.York, NY.He received the B.S. degree in chemistry from St. John's University in 1964, and the Ph.D. degree in chemical physics from the Massachusetts Institute of Technology on 1970.He has done Post-doctoral research at MIT and at the National Institute of Standards and Technology (previously the National Bureau of Standards) as an National Research Council Post-doctoral Fellow.