Qiang Yu

Tadahiro Shibutani

Masaki Shiratori

Tomio Matsuzaki

Tsubasa Matsumoto

MODELING WITH STRUCTURE OF RESINS IN ELECTRONIC COMPONENTS

In recent years, interfacial fracture becomes one of the most important problems in the assessment of reliability of electronics packaging. Especially, underfill resin is used with solder joints in flip chip packaging for preventing the thermal fatigue fracture in solder joints. In general, the interfacial strength has been evaluated on the basis of interfacial fracture mechanics concept. However, as the size of devices decrease, it is difficult to evaluate the interfacial strength quantitatively. Most of researches in the interfacial fracture were conducted on the basis of the assumption of the perfectly bonding condition though the interface has the micro-scale structure and the bonding is often imperfect. In this study, the mechanical model of the interfacial structure of resin in electronic components was proposed. Bimaterial model with the imperfect bonding condition was examined by using a finite element analysis (FEA). Stress field in the vicinity of interface depends on the interfacial structure with the imperfect bonding. In the front of interfacial crack tip, the behavior of process zone is affected by interfacial structure. However, the instability of fracture for macroscopic crack which means the fracture toughness is governed by the stress intensity factor based on the fracture mechanics concept.

INTRODUCTION

In recent years everyday life has become more and more convenience through significant progress in computerization. On the other hand, there is great stress put upon the technology that supports these new technologies. There is no end to demanding smaller, lighter, more high tech and faster electronic devices, and the modern field of science is struggling with many problems that need to be overcome. These problems include expanding the variation of forms and materials used in the production, adaptation to a wider range of environment where these technologies may be used, making device forms more complicated, consideration of environmental issues (i.e. reduction of lead usage, recycling, long life expectancy designing), and foundation of a method for reliability designing. One of the significant methods to improve the strength of chip assemblies these days is the dispensing of resin between the spaces between the silicon chip and the substrate, where the circuits come in action through solder balls 1),2) . The most critical failures in chip assembles are the disconnections of the circuit, and the joints of different materials, such as the chips and the solder ball joints, are the most susceptible to failure. These failures occur when stress concentrates at the weak points in the structure often caused by shocks and differences in the coefficient of thermal expansion (CTE) in different materials. Resin materials are used for the purpose of redistributing the stress and reducing the effect of differences in CTE.

Since electronic devices consist of various materials, there are many interfaces: metal/ceramic, metal/resin, ceramic/resin and so on. Therefore, interfacial cracking affects the reliability of devices. Interfacial strength has been often evaluated as the fracture toughness on the basis of fracture mechanics concept 3),4) . Perfect bonding was assumed and stress intensity factor is dominant parameter of interfacial failure. However, as the size of devices decreases, interfacial fracture behavior becomes more complicated. Figure 1 shows an example of interfacial cracking between resins. These materials are bonded partly and interfacial structure can be seen at the crack tip. As the devices size becomes smaller, the effect of interfacial structure appears. Additionally, in the interfacial facture mechanics, the failure process of interface has been ambiguous yet. Therefore, interfacial structure-based evaluation should be constructed. In this study, the simple partly bonded interface model was proposed. FEA for the proposed model was compared with the bimaterial interface model. The effect of interfacial structure was extracted. Considering the failure criteria of the interfacial structure, the behavior of interfacial fracture was examined.

INTERFACIAL STRUCTURE BASED MODEL

Figure 2 shows the bimaterial model with interfacial structure. Bonding parts of which width and height are a and h, are arrayed periodically with the interval s. FEA analysis was conducted under two boundary conditions: uniaxial tensile and thermal loading. In the former cases as shown in Fig. 2, the bottom and right sides were fixed. Uniform tensile stress was applied on the upper of resin. In the latter case, temperature changed from 150 C to 20 C. No stress is applied at 150 C. Five models were prepared where sizes of bonding structure, a and h range as shown in Table I. 4-node plain strain element was used and their numbers of elements and nodes are shown in Table I. Fine mesh was constructed near the interfacial zone. The solver is Marc 2005. Materials are assumed to be isotropic elastic body and each material property is listed in Table II. Additionally, by considering the fracture property of interfacial structure, the failure behavior at the crack tip was investigated. In this case, the interfacial structure and resin is assumed to be elasto-plastic body. In order to determine the criteria of fracture, the maximum plastic strain was employed. Using the subroutine "uactive", elements which reach the criteria are deactivated. along the interface, respectively. Each plotted data is taken at the middle of interfacial structure. The stress curves in bimaterial model with perfect bonding are also shown by solid lines. When s is small, the curve is close to the solid line. Since separated area is small, the structure is almost that in the perfect bonding. Other curves are higher than the solid line. When interfacial height h is low, the stress σ y increases. The mismatch of deformation concentrates near the interface and the strain in the normal direction depends on the height of interfacial structure. On the other hands, in the case of τ xy , when interfacial height h is low, the stress decreases. It suggests that the interfacial structure affects the mix mode condition at the edge of interface. The effect of interface's height appears in the middle of model. Though the horizontal mismatch of displacement is generated at the edge of interface, the vertical mismatch affects the shear stress in the middle of model.

RESULT AND DISCUSSION

Uni-axial tensile

Residual stress

Figure 4 (a) and (b) shows the distribution of the residual stresses. When a=5 and s=5, (the fraction of bonded area is 50%), the residual stress, σ y increases. Not only the bonded area decreases, and but also the mismatch of deformation is changed due to the rigidity of interfacial structure. It implies that the decrease in the bonding area causes that the effect of interfacial structure appears. When the height of interfacial structure, h becomes high, the stress increases. The effective zone of localized deformation due to the mismatch depends on the bonding structure. When h is low, not only the interfacial structure but also bulk materials can be deformed by the mismatch of CTE. On the other hand, when h is high, the deformation by mismatch concentrates on the bonding area.

Comparing with the uniaxial tensile state, the effect of interface height on the shear stress is not so big. CTE mismatch does not depend on the interfacial structure and subject each interfacial structure. Therefore, the width and interval of bonding part mainly affect the behavior of shear stress.

Fracture of interface

Figure 5 and 6 show the failure process, near the crack tip under uniform tensile stress state. The fracture in the interfacial structure occurred in the front of the tip. As the applied stress increase, failure process proceeds. Finally, instability of fracture takes place. Figure.7 and Figure.8 are the stress distribution from the tip at the instability of fracture. When the stress singularity field near the crack tip reaches critical value, the instability begins. It implies that macroscopic stress field governs the interfacial fracture. On the other side, stress depends on the interfacial structure in r < 100 µm. Therefore, the process zone is affected by the interfacial model. In the case of electronics components, the size of failure is less than that of process zone. It suggests the importance of the interfacial structure model concept proposed in this paper.

CONCLUSION

1) Interfacial stress behavior in the state of elasticity of the proposed two dimension model was able to be confirmed.

2) The stress field along the interface depends on the interfacial structure. At the edge of interface, mixed mode state is affected by the bonding structure. 3) In the front of the crack, the process zone is affected by the interfacial structure. On the other hand, when the crack length is long, the fracture toughness obeys on the stress singularity field along the interface. 4) In electronics components, the fracture behavior of process zone affects the reliability of devices. FEA results in this study suggest the importance of the interfacial structure.

Figure 3

 3 Figure 3 (a) and (b) are distributions of σ y and τ xy

REFERENCES 1 )

 1 Hirohata, K., Kawamura, N., Mukai, M., Kawakami, T., Aoki, H., Takahashi, K., "Mechanical fatigue test method for chip/underfill delamination in flip-chip packages", IEEE Trans. on Electronics Packaging Manufacturing, Vol. 25, No. 3, pp. 217-222, 2002. 2) Dutta, I, Gopinath, A, Marshall, C, "Underfill constraint effects during thermomechanical cycling of flip-chip solder joints", Journal of Electronic Materials, Vol. 31, No. 4, 2002. 3) Timothy Ferguson and Jianmin Qu, "Effect of Moisture on the Interfacial Adhesion of the Underfill/Solder Mask Interface", Journal of Electronic Packaging, Vol. 124, No. 2, pp. 106-110, 2002. 4) Zhong, Z.W. , Wong, K.W. , Shi, X.Q., "Interfacial behavior of a flip-chip structure under thermal testing", IEEE Trans. on Electronics PackagingManufacturing, Vol. 27, No. 1, pp. 43-48, 2004.

Fig. 1

 1 Fig.1 Interfacial structure between resins in the front of crack tip with imperfect bonding.

Fig. 3 .

 3 Fig.3.(a) Distributions of σ y when uniform tensile stress was applied on the upper of resin (h1=h2=500 w=1000 µm)

Fig. 3 .

 3 Fig.3. (b) Distribution of τ xy when uniform tensile stress was applied on the upper of resin (h1=h2=500µm w=1000µm)

Fig

  Fig.4. (a) Distribution of σ y (20C from 150C h1=h2=500 w=1000µm)

Fig

  Fig.4. (b) Distribution of τ xy (20C from 150C h1=h2=500µm w=1000µm)

FigFig. 7 .Fig. 8 .

 78 Fig.5. Crack growth from the crack tip

Table II

 II

						Complete Interface
		4.00E+01				a=5 s=5 h=1
						a=5 s=3 h=1
		3.50E+01				a=5 s=1 h=1
		3.00E+01				a=5 s=5 h=0.5
						a=5 s=5 h=2
	σy(MPa)	2.00E+01 2.50E+01
		1.50E+01
		1.00E+01
		5.00E+00
		0.00E+00
		0	100	200	300 x(µm)	400	500	600
								Material properties
								Young's modulus (GPa) Young's modulus (GPa)	Poisson's ratio (-) Poisson's ratio (-)	CTE (*10 -6 /K) CTE (*10 -6 /K)	Initial yield stress (MPa) Initial yield stress (MPa)
								resin resin	10 10	0.3 0.3	19 19	60 60
								interface interface	10 10	0.3 0.3	19 19	60 60
								silicon silicon	130 130	0.26 0.26	3.6 3.6	--


A Augustin
email: adam.augustin@freescale.com

T Hauck

B Maj

J Czernohorsky

E B Rudnyi

J G Korvink
email: korvink@imtek.de

Model Reduction for Power Electronics Systems with Multiple Heat Sources

The paper demonstrates the model order reduction procedures applied to semiconductor devices with multiple heat sources. The approach is demonstrated for a device with nine heat sources where some of them are permanently active and other work under switching conditions. For the order reduction the software package MOR for ANSYS is used, which is based on the Krylov subspace method via the Arnoldi algorithm.

I. INTRODUCTION

Modern semiconductor technologies allow the development of power electronics systems with manifold functionality integrated on one single die. Due to this everlasting process, the devices include more and more heat sources, which have to be considered in thermal simulations. For complex package and heat source configurations as well as for different applications, accurate thermal models are needed.

As module configuration is close to a multilayered block structure, one of the simulation approaches is based on an analytical solution of the heat transfer equation [START_REF] Augustin | Two-dimensional temperature fields in ASICs for multiple heat sources of arbitrary power dissipation derived from conformal mapping theory[END_REF]- [START_REF] Khatir | Real-time computation of thermal constraints in multichip power electronic devices[END_REF]. Obviously, real geometry is different from the idealized geometry required by the analytical solution and the finite element method seems to be the only valid alternative for detailed thermal simulations in the general case (see for example [START_REF] Ammous | Choosing a thermal model for electrothermal simulation of power semiconductor devices[END_REF]). Unfortunately, the biggest disadvantage of this approach is high simulation time because of high dimensional models with hundreds thousands degrees of freedom. One possible way to overcome this problem is to use mathematical methods in order to reduce the system order with a minimum loss of accuracy.

A formal model reduction approach [START_REF] Antoulas | Approximation of large-scale dynamical systems: An overview[END_REF] (overview from engineering viewpoint in [START_REF] Rudnyi | Review: Automatic Model Reduction for Transient Simulation of MEMS-based Devices[END_REF]) allows us to take a high dimensional finite element model and generate its lowdimensional approximation. As such, it is an ideal candidate for the goal above. Several research groups have already documented its successful application to a thermal problem [START_REF] Bechtold | Automatic generation of compact electro-thermal models for semiconductor devices[END_REF]- [START_REF] Yang | Extraction of heat-transfer macromodels for MEMS devices[END_REF].

However, the formal model reduction approach also puts some limits for the original thermal problem. For example, it happens that its application is not straightforward for cooling of electronic components as considered by the DELPHI and PROFIT consortium [START_REF] Rosten | The world of thermal characterization according to DELPHI -Part I: Background to DELPHI[END_REF]- [START_REF] Lasance | Two benchmarks to facilitate the study of compact thermal modeling phenomena[END_REF]. The main reason is that the use of film coefficients is limited for forced convection within a volume with complex geometrical constraints. In this case, one cannot avoid a joint thermal-CFD simulation and compact thermal models should possess special properties.

In our case, the situation is different. 85 % of heat is removed by conduction and there is no forced convection. The remaining 15 % of heat is removed by natural convection. Additionally, film coefficients have been validated by means of experimental monitoring of temperature. As a result, formal model reduction can be employed for the thermal problem without any further modification. Finally, we should stress that dynamic behavior must be preserved in our case.

The structure of this paper is as follows. In section II a short introduction into the model order reduction theory will be presented together with the reduction process for devices with multiple heat sources. In section III the necessary preparatory steps for calculations with reduced models will be explained including the consideration of some permanent active heat sources. In section IV the presented simulation results will be discussed. Finally, in section V some concluding remarks will be made.

II. MODEL REDUCTION PROCESS

Using the finite element simulation software ANSYS, a geometrical and physical model is generated. The model should concern all thermally relevant components to represent the real device as good as possible.

A. General Theory

After the discretization, the finite element model is expressed as follows:

E dT (t) dt + K T (t) = B u(t) y(t) = C T (t) (1)
where y(t) is the vector of unknown temperatures at the nodes, and T (t) is the state vector. E and K are the heat capacity and heat conductivity system matrices, B is the input matrix, and C is the output matrix. The vector u(t) comprises n input functions for the n heat sources. The output matrix specifies particular linear combinations of temperatures which are of interest to an engineer. The system matrices as well as the load vectors have been read from ANSYS binary *.full and *.emat files. Model reduction is performed in the Laplace domain, that is, for the transfer function of eq. ( 1)

H(s) = C (sE + K) -1 B ( 2
)
and it is based on an assumption that there exists a lowdimensional subspace V that accurately enough captures the dynamics of the state vector T (t)

T (t) = Vz(t) . (3)
Obviously, z(t) indicates here the state vector of the reduced system. To generate the low-dimensional subspace V we expand eq. ( 2) around the Laplace variable s 0 = 0 and then we use implicit moment matching via the Krylov subspace method. Finally, for the reduced system we obtain

V T E V dz(t) dt + V T K V z(t) = V T B u(t) y(t) = C V z(t) . (4)

B. Reduction Process for Multiple Heat Sources

In order to generate the reduced model with multiple heat source devices the right strategy has to be chosen. For example, for applications under switching conditions in kHz range, model generation for each particular switching event is not practicable, since for the description of only few second response, hundreds of load steps are needed. Fortunately, the input function in eq. ( 1) does not take part during model reduction process and can be passed to the reduced model in eq. ( 4) without changes. This means that provided we have split the load to the product of B and u, for each particular heat source only one load step has to be generated. Switching conditions can be implemented afterwards by applying pulsed input functions to the reduced model. This procedure allows a very comfortable control on the input function characteristics without any additional model order reduction runs.

As a demonstration, we have chosen a device with nine heat sources. In Fig. 1 the finite element model of the investigated device is shown. Since the printed circuit board with the heat sink is much bigger than the package only the much more interesting section around the investigated devices (together with the package) is shown. The device is assembled in a package with an exposed lead frame. This allows a direct connection between the package and the heat sink, which is protruded from the bottom side to the top surface of the board. Fig. 2 indicates schematically the die in a top view with all nine investigated heat sources. Five of them are permanently active, the other four sources are dynamically switched on and off. In principle, only five load steps could be generated, as all the permanent sources could be joined with each other. However, we have generated nine separate load steps for each heat source to If the power dissipation is set to 1 W during the load generation, the effective power dissipation in the reduced model can be controlled by the amplitude of the input function. Obviously, using this method the input functions for five heat sources are constant functions f i (t) = P i with i = 1, . . . , 5 and the effectively dissipated power P i . The input functions for the dynamic heat sources are pulse functions as shown schematically in Fig. 3 with f j (t) = P j • f pulse,j (t), j = 6, . . . , 9 and the amplitudes P j . In order to draw a distinction between the input functions f pulse,j (t) in Fig. 3 both functions are scaled. Finally, the input matrix B is a m × 9 matrix, with m as the order of the reduced model.

Because of the special operation mode of MOR for ANSYS [16] all boundary conditions are taken from the first load step file, therefore these conditions must be set only for the first load step besides the definition of power dissipation for the first Two heat sources are using the function (solid) with 0.4 s period and the other two use the 0.8 s period function (dashed). The amplitudes are adjusted later to appropriate values heat source. For the remaining eight load steps, only a particular source (two until nine) was active in each load generation step. Finally, nine load step files are created by ANSYS: file1.full (convection and only source one is on with P = 1 W), file2.full (only source two is on with P = 1 W), file3.full (only source three is on with P = 1 W), and so on until nine. As was mentioned before, the dynamic state of the four transient heat sources will be implemented later, simply by the use of the pulsed input function as shown in Fig. 3. In order to begin the reduction process all output positions (nodes in the finite element model) must be set, for which the temperature will be monitored, as well. Using all generated files and the information about the outputs, the reduction process can be performed with MOR for ANSYS.

III. PREPARATORY STEPS

At the beginning of the simulation process, on one hand, the decision about the order of the reduced model must be made as well as, on the other hand, the heat generation of the permanent heat sources must be considered.

A. Model Order Estimation

Implicit moment matching does not have global error estimates. In order to determine the dimension of the reduced model, we suggest to use a local error estimate at some frequency that is roughly inversely proportional to the rise time. Fig. 4 displays logarithm of the relative error for three frequencies, 0.1 Hz, 1 Hz and 10 Hz. We define the relative error as

H f ull (f ) -H red (f ) H f ull (f ) (5)
with H as transfer functions of the reduced and full system.

As the expansion point was zero, then the higher the frequency the slower the convergence. The error of 1 % can be approximated by first fifty generalized states at 0.1 Hz. At 1 Hz the error of 1 % can be reached at the dimension of 100, first. Based on these results, we chose the order of

B. Consideration of Permanent Heat Sources

To perform transient calculations with permanent heat sources one have to shift the pulse inputs (necessary for the dynamic sources) by a sufficient time in order to reach the steady state caused by all permanent sources as shown in Fig. 5. As soon as this state is reached, the pulsed heat sources can be activated. In terms of simplicity, we start here with 2000 s to make sure that the steady state is reached. In reality the steady state is reached much faster.

T [ • C] t [s]

C. Model Verification by Measurements

Since the investigated device contains several build-in temperature sensors distributed over the entire die surface, we have used seven of them to verify the finite element model. First, all sensors have been calibrated separately. Using this calibration approach we estimated an accuracy of 3 % for the steady state. Using the calibrated model, transient measurements and simulations with the reduced model were done.

IV. RESULTS

First, we show here the improvements in time consumption using the reduced model for transient simulations.

In ANSYS for a transient simulation with 15 s duration time and a pulse period of 0.4 s, at least 38 load steps are needed. Therefore, transient simulation with very bad accuracy due to high numerical errors during integration in time needs at least 874 s of runtime. The entire calculation with the reduced model, including the model reduction process and a time vector with 4100 time steps, needs 158 s of runtime, from which 150 s take for the model reduction and 8 s for transient integration. Note that integration of the full model (with 104346 nodes) in ANSYS for 4100 time steps would take about 26 hours.

For the presentation of the results we chose two sensors/outputs. As can be seen in Fig. 2 one build-in sensor is placed near the permanent active heat source and the other in the vicinity of a dynamic source. In Fig. 6 sensor shows very high gradients, caused by the vicinity to the dynamic heat sources. The measurements at this sensor are indicated by dots. Since the distance of the other sensor to the dynamic sources is much larger, the gradients here are much lower. On both solid curves, the switching times of the particular sources can be seen very good. However, the left sensor shows these switching activities much more pronounced. The measurements (dots and dashed line)

show a small difference that can be explained by both measurement errors and approximations made during the finite element modeling of the original high dimensional model. Nevertheless, a very good correlation with the simulated results can be observed.

V. CONCLUSIONS

In this paper, we have shown the application of the model order reduction method for a power device with multiple heat sources. The treatment of two different heat sources (permanently active and dynamic sources) was explained. Using complex measurements, the temperature was monitored with several build-in temperature sensors on the die, in order to verify the finite element model. All simulation results, for steady state and transient conditions, show very good correlations with measurements.

Depending on the frequency of input functions, the order of the reduced model can be reduced down to 100. However, for faster switching conditions multiple expansion points should be used, in order to keep the dimension of the reduced model applicable.

Fig. 1 .Fig. 2 .

 12 Fig. 1. ANSYS model of a power device placed on PCB with a heat sink. To have a view on the die, the enclosing mold compound has been removed

Fig. 3 .

 3 Fig. 3. Schematical view of the input functions for transient simulations. Two heat sources are using the function (solid) with 0.4 s period and the other two use the 0.8 s period function (dashed). The amplitudes are adjusted later to appropriate values

2 Fig. 4 .

 24 Fig. 4. Local error estimation for the reduced model depending on the frequency and model order

Fig. 5 .

 5 Fig. 5. Temperature rise on one output showing the consideration of the permanent active heat sources. To be certain for a steady state and in terms of simplicity, all pulse functions for the dynamic heat sources are shifted by 2000 s

Fig. 6 .

 6 Fig. 6. Temperature rise caused by four dynamic heat source with two different input functions as well as by five permanent active sources at T amb = 25 • C. Solid lines indicate the simulation results of the reduced model with order 100. Dots (left sensor in Fig. 2) and the dashed line show transient thermal measurements. Time scale is shifted to the begin of the pulsed input functions


Dirk Schweitzer

THERMAL TRANSIENT MULTISOURCE SIMULATION USING CUBIC SPLINE INTERPOLATION OF ZTH FUNCTIONS

This paper presents a very straightforward method to compute the transient thermal response to arbitrary power dissipation profiles in electronic devices with multiple heat sources. Using cubic spline interpolation of simulated or measured unit power step response curves (Zthfunctions), additional errors due to model reduction can be avoided. No effort has to be spent on the generation of compact models. The simple analytic form of the interpolating splines can be exploited to evaluate the convolution integral of the Zth-functions with arbitrary power profiles at low computational costs. An implementation of the algorithm in a spreadsheet program (EXCEL) is demonstrated. The results are in very good agreement with temperature profiles computed by transient Finite Element simulation but can be obtained in a fraction of the time.

INTRODUCTION

The trend towards higher integration and higher power densities in semiconductor devices has also driven the needs of methods to predict the transient temperature of devices with multiple heat sources. Since the computational effort for direct thermal simulation of transient switching processes with detailed (e.g. finite element) models is huge, several methods of model reduction have evolved which drastically reduce computing time.

If material properties and boundary conditions are independent of temperature, a single heat source is thermally fully characterized by its unit power step response or Zth function Z th (t). I.e. the temperature response T(t) to any power dissipation profile P(t) can be computed by the convolution integral

∫ ⋅ - + = t th d P t Z T t T 0 0 τ τ τ ) ( ) ( ) ( & , (1)
T 0 being the ambient temperature [START_REF] Gerstenmair | A New Procedure for the Calculation of the Temperature Development in Electronic Systems[END_REF]. For real materials with temperature dependent density, specific heat, and thermal conductivity, equation [START_REF] Gerstenmair | A New Procedure for the Calculation of the Temperature Development in Electronic Systems[END_REF] holds not exactly true, but the errors can be neglected as long as the temperature difference T(t) -T o is not too large. For N heat sources the mutual heating has to be taken into account, extending eq.

(1) to τ τ τ

d P t Z T t T N j t j ij th i ∑ ∫ = ⋅ - + = 1 0 0 ) ( ) ( ) ( & , (2)
for the temperature T i (t) of the i-th heat source [START_REF] Gerstenmair | Calculation of the Temperature Development in Electronic Systems by Convolution Integrals[END_REF]. Z th ij (t) is a transfer function, which describes the temperature rise at the location of source i, when only source j is heated. Again the superposition principle strictly speaking holds true only for temperature independent material properties but is a good approximation in most cases.

Based on these principles, any reduced model which correctly reproduces the thermal (transfer-) Zth curves can be used to simulate the temperature development at the locations of the heat sources of a device. In most cases thermally equivalent RC networks, which had been fitted to simulated or measured power step responses, have been used for this purpose [START_REF] Gerstenmair | Calculation of the Temperature Development in Electronic Systems by Convolution Integrals[END_REF] [3] [START_REF] Hauck | Thermal RC-Network Approach To Analyze Multichip Power Packages[END_REF]. However the process of model reduction introduces new errors into the simulation in addition to the errors inherent to the underlying Zth functions. An obvious idea is therefore to use directly the Zth values Z th ij (t i ) sampled at times t i and interpolate at intermediate points of time.

INTERPOLATION OF THE ZTH FUNCTIONS

In the following it is assumed that a thermal unit power step response Z th (t) has been sampled at m points of time t i which are equidistant on a logarithmic timescale. The interpolation can either be done on a linear timescale, i.e. for points

{ } ) , ( i th i Z t , i = 1 .. m (3)
or on a logarithmic timescale, i.e. for points

{ } ) , ( i th i Z z , i = 1.. m with i i t z ln = (4)
For a Foster RC network (figure 1), the exact value of the Zth function of its junction node can be calculated by the equation

) ( ) ( ∑ = - - = n i C R t i th th i th i th e R t Z 1 1 (5)
for any time t. Using eq. [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF] easily perform some numerical experiments to compare the deviations between interpolated and exact Z th curves for different interpolation methods.

Figure 1: Foster RC-network

Infineon's TLE 6220 in a P-DSO-20 exposed pad package, an automotive device for engine control with four power dissipating channels on a single chip, will be used to demonstrate the methods presented herein. More details on the device model follow below. At the moment the RC parameters fitted for this device (table 1) just serve as an example. Using eq. ( 5), a data set containing m=51 points, equidistant on a logarithmic timescale, has been generated (figure 2). [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF], this data set has been interpolated both on a linear timescale (interpolating between points given by eq. ( 3)) and on a logarithmic timescale (interpolating between points given by eq. ( 4)). The resulting relative error between interpolated and analytic thermal impedance,

i R th i [K/W]

C
) ( ) ( ) ( . t Z t Z t Z err rel th th ip th - = , ( 6
)
is shown in figure 3. The maximum interpolation error is 0.03% for interpolation on a linear timescale and 0.01% for interpolation on a logarithmic timescale over the whole time range from 1.0e-6 s to 1000 s. These values depend of course on the shape of the interpolated Zth curve as well as on the distance of the interpolation points. But for a reasonable number of points, the interpolation error can be neglected. The errors that occur when fitting a RC-network to a unit power step response are typically in the range of a few percent (and often much larger for the transfer functions).

NUMERICAL CONVOLUTION

Cubic spline interpolation not only saves one from the trouble of fitting RC networks but can also be exploited to efficiently compute the convolution integrals in eq. ( 1) and (2). Any power dissipation profile P(t) can be approximated by a sequence of trapezoidal pulses as shown in figure 4. Using the superposition principle, the temperature response to the whole power profile is just the sum of the temperature response curves to all trapezoidal pulses. For a single pulse, starting at time t 1 with power P 1 and ending at time t 2 with power P 2 , partial integration of the convolution integral in eq.( 1) yields:

             ≥ ⋅ - - + - - - < ≤ ⋅ - - + - < = ∆ ∫ ∫ - - - 2 1 2 1 2 2 2 1 1 2 1 0 1 2 1 2 1 1 1 , ) ( ) ( ) ( , ) ( ) ( , 0 ) (
) ( ) ( T t T t T - = ∆ .
The integrals in eq. ( 2) can be evaluated accordingly. I.e. to compute the temperature response T(t) for a trapezoidal pulse, not only the Z th (t) has to be evaluated, but also the integral:

∫ ≡ t th th d Z t IZ 0 ) ( ) ( τ τ (8)
Both can be easily computed since the (piecewise) interpolating splines are just third degree polynomials which can be integrated analytically: 9) is the i-th interpolating spline between Z th values given at times t i and t i+1 and its integral if the interpolation is done on a linear timescale. For cubic spline interpolation on a logarithmic timescale, the interpolating spline is

1 0 2 1 2 1 3 2 3 1 4 3 4 1 0 1 2 2 3 3 , ) ( ) ( ) ( ) ( ) ( ) ( + < ≤ - + - + + - + - = + + + = ∫ i i i i i i i i i i t t
+ = + i i th i ip th Z t Z . Eq. (
1 0 1 2 2 3 3 + < ≤ = + + + = i i i i i i ip th z z z t z b z b z b z b z Z i , ln , ) ( ~ (10) with i i th i ip th Z z Z = ) ( ~ and 1 ) ( ~1 + = + i i th i ip th Z z Z .
The coefficients a ij (b ij ) are computed such as to ensure that the interpolating function is smooth in the first derivative and continuous in the second derivative (see [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing[END_REF] for a detailed description of the algorithm). In the latter case (eq. ( 10)) the transformation t → z yields

∫ ∫ ∫ = = z z ip th t t t t ip th ip th i i i i i i d e Z d Z d Z ζ ζ τ τ τ τ ζ ) ( ) (ln ) ( (11)
Which can be further evaluated using the elementary integrals

∫ - = ) 1 (x e dx e x x x ∫ + - = ) 2 2 ( 2 2 x x e dx e x x x (12) ∫ - + - = ) 6 6 3 ( 3 3 x x x e dx e x x x . to ) ( ) ( ) ( ~0 3 2 2 3 3 0 1 2 2 3 3 1 i i i i i i i z i i i i z z z ip th c z c z c z c e c z c z c z c e d e Z i i + + + - + + + = ∫ ζ ζ ζ (13) with 3 2 1 0 0 3 2 1 1 3 2 2 3 3 6 2 6 2 3 i i i i i i i i i i i i i i b b b b c b b b c b b c b c - + - = + - = - = = (14)
The coefficients a ij or b ij and c ij should be stored for subsequent evaluations of Z th and IZ th .

PERIODIC POWER PULSES

If the temperature response to a sequence of periodic power pulses is to be computed for a time t much larger than the pulse period t p (t >> t p ), the contributions of all preceding pulses within the decay time t decay (i.e. the time it takes until the rise in temperature due to a single power pulse has dropped to a negligible value) must be taken into account. The ratio t decay / t p and therefore the number of pulses to be included can be quite large, meaning a considerable computational effort.

But the difference between the temperature response to a periodic power excitation and to a constant power excitation with the corresponding average power vanishes quickly after the power has been switched off (figure 5). In the long run the decay curve depends merely on the total energy dissipated during the pulse(s), not on the

Dirk Schweitzer Thermal Transient Multisource Simulation Using Cubic Spline Interpolation of Zth functions.

©TIMA Editions/THERMINIC 2006 -page-ISBN: 2-916187-04-9 specific shape of the power profile P(t). Therefore the computing time can be reduced drastically if the first n pulses of the periodic sequence are replaced by a single constant pulse of length n•t p and average power P avg = E Pulse / t p , E Pulse being the energy per pulse. In practice, n is chosen such that the temperature response ∆T n (t) to the n-th pulse is < ε (small) at the time t of interest and ≥ ε for all subsequent pulses > n (figure 6).

APPLICATION TO AN AUTOMOTIVE ENGINE CONTROL DEVICE

The Infineon device TLE 6220 in a P-DSO-12 package controls the ignition of a four-cylinder engine. The power is dissipated periodically in four regions on the chip, named channel 1-4 (figure 7). The exposed pad of the package is soldered to a PCB, the bottom side of which being fixed at a temperature of 125°C.

The methods described herein shall be used to compute for this device the transient temperature response to a sequence of application specific pulses. Due to the symmetry of the device, all necessary (transfer) Zthfunctions can be obtained by one transient simulation, heating only channel 1 and monitoring the temperature rise at the center of channels 1-4. Using the macro programming language of the Finite Element software ANSYS (APDL), the simulated step responses are written directly to a file which can be imported by a simulator without further post-processing. The standardized XML file format has been chosen for this purpose since it can be easily extended and XML parser exist for almost all programming languages.

The Zth functions have been computed for 74 equidistant points on a logarithmic timescale, ranging in between 0.5 µs and 1000 s (steady-state). Figure 8 shows the power profile which is to be applied to the channels: after a switch-on phase of 0.5 ms during which the power rises linearly to 1.0 W, the power remains constant until the channel is switched off after 4 ms. The energy stored in a connected load inductance is released within 0.3 ms, causing a peak power dissipation of 60 W in the device. This power profile is subsequently applied to all channels each shifted by 4 ms. . Using these Zthfunctions, an EXCEL simulator which implements the algorithm described above can compute the same temperature response within seconds (figure 9). In addition, changes to the power profile can be analyzed immediately, which would otherwise require a new timeconsuming FE simulation.

Figure 9: EXCEL Thermal Transient Simulator

In Figure 10 the temperature curves obtained by FE simulation and by the EXCEL simulator are so close, that differences are not distinguishable. The absolute error plot shows that the maximum difference < 0.3°C occurs during the short high power pulses. Discrepancies arise from numerical errors and non-linearities of the FE model (thermal conductivity of silicon).

CONCLUSION

Using interpolated Zth functions for the evaluation of the convolution integral in eq. ( 2), simulated or measured Zth curves can be used directly for the calculation of the temperature response to arbitrary power profiles in multi- source devices. This can be several heat-sources on a single chip as well as multi-chip devices or even several devices on a PCB. An implementation of the algorithm in EXCEL has demonstrated that the temperature curves for the selected source nodes of the model can be computed much faster than by a pure Finite Element simulation, even with the interpreted VBA programming language of EXCEL. This enables a fast analysis of different power profiles as well as the simulation of fast switching processes over a long period of time. A compiled code would perform even better, allowing also the simulation of temperature dependent power sources. The approximation error of the cubic spline fit is very small; deviations are mainly caused by non-linearities of the FE model.

Figure 2 :Figure 3 :

 23 Figure 2: Model curve to be interpolated

Figure 4 :

 4 Figure 4: Approximation of arbitrary power profiles

Figure 5 :Figure 6 :

 56 Figure 5: Temperature response curves for a periodic sequence of 5 pulses and for the corresponding constant pulse with average power P avg .

Figure 7 :

 7 Figure 7: Finite Element model of the TLE 6220

Figure 8 :

 8 Figure 8: Power profileIn a direct transient FE simulation 544 points of time were necessary to resolve the temperature response in a time range from 0 to 20 ms. The simulation took about 11 h on a 2.8 GHz CPU compared to 90 min for the computation of the Zth-functions (74 points). Using these Zthfunctions, an EXCEL simulator which implements the algorithm described above can compute the same temperature response within seconds (figure9). In addition, changes to the power profile can be analyzed immediately, which would otherwise require a new timeconsuming FE simulation.

Figure 10 :

 10 Figure 10: Comparison FE simulation -EXCEL simulator

  as a model function, one can Dirk Schweitzer Thermal Transient Multisource Simulation Using Cubic Spline Interpolation of Zth functions.

©TIMA Editions/THERMINIC 2006

-page-ISBN: 2-916187-04-9

th i [Ws/K]

	1	0.07746	0.000273
	2	0.48958	0.000365
	3	1.55159	0.000772
	4	2.70800	0.003253
	5	1.45388	0.052680
	6	2.51305	0.274200
	7	1.05932	4.496140
	8	0.11046	266.4490

Table 1 : Parameters of the test model Using cubic splines

 1


Peter E Raad

Pavel L Komarov

Mihai G Burzo

A COUPLED THERMOREFLECTANCE THERMOGRAPHY EXPERIMENTAL SYSTEM AND ULTRA-FAST ADAPTIVE COMPUTATIONAL ENGINE FOR THE COMPLETE THERMAL CHARACTERIZATION OF THREE-DIMENSIONAL ELECTRONIC DEVICES: VALIDATION

This work builds on the previous introduction [1] of a coupled experimental-computational system devised to fully characterize the thermal behavior of complex 3D submicron electronic devices. The new system replaces the laser-based surface temperature scanning approach with a CCD camera-based approach. As before, the thermo-reflectance thermography system is used to noninvasively measure with submicron resolution the 2D surface temperature field of an activated device. The measured temperature field is then used as input for an ultra-fast inverse computational solution to fully characterize the thermal behavior of the complex threedimensional device. For the purposes of this investigation, basic micro-heater devices were built, activated, and measured. In order to quantitatively validate the coupled experimental-computational system, the system was used to extract geometric features of a known device, thus assessing the system's ability to combine measured experimental results and computations to fully characterize complex 3D electronic devices.

INTRODUCTION

Faster and more powerful devices mean hotter devices, which can lead to a decrease in performance and reliability. Thus, understanding and determining the thermal behavior of modern electronics has become a key issue in their design. As a result, there is a critical demand for methods that can be used to determine the temperature of features at the submicron level, particularly when most important features are physically inaccessible [START_REF] Goodson | Short-time-scale Thermal Mapping of Microdevices using a Scanning Thermoreflectance Technique[END_REF][START_REF] Grauby | High Resolution Photothermal Imaging of High Frequency Using Visible Charge Couple Device Camera Associated with Multichannel Lock-in Scheme[END_REF][START_REF] Tessier | Quantitative Thermoreflectance Imaging: Calibration Method and Validation on a Dedicated Integrated Circuit[END_REF][START_REF] Grauby | Qualitative Temperature Variation Imaging by Thermoreflectance and SThM Techniques[END_REF]. Computational approaches can provide insight into the internal thermal behavior of such complex devices, but can be limited by the inherent necessity of modeling the heat sources, which in the case of self-heating microelectronic devices, are the result of electrical fields whose exact shapes and locations are difficult to specify with reasonable certainty. Moreover, such devices can actually experience irreversible changes in thermophysical properties and/or geometries that cannot be otherwise predicted from theory or monitored. Experimental approaches can also be helpful in determining thermal behavior, but require either physical access or a visual path to the region of interest. Contact methods, for example, present the difficulties of having to access features of interest with an external probe, or in the case of embedded features, fabricate a measuring probe into the device, and then having to isolate and exclude the influence of the probe itself. Non-contact methods, on the other hand, can provide surface temperature profiles, but in and of themselves cannot impart information on internal behavior. In other words, these methods provide a two-dimensional perspective on what otherwise is, in the case of stacked complex devices, an intricate three-dimensional thermal behavior.

We show herein that by combining an experimental method capable of mapping the surface temperature of a complex device with high spatial and temporal accuracy together with a computational engine capable of rapidly and accurately resolving the geometric and material complexities of a full three-dimensional microelectronic device, it becomes possible to use the independent information from the experimental measurements to mitigate the lack of knowledge in the source model parameters, which directly affect the usefulness of the computational results. This article builds on the previous introduction by the authors of a proof of concept of a coupled computational-experimental approach that uses a measured two-dimensional surface temperature mapping to help obtain a fully three-dimensional thermal characterization of an active micro-device [START_REF] Raad | Coupling Surface Temperature Scanning and Ultra-Fast Adaptive Computing to Thermally Fully Characterize Complex Three-Dimensional Electronic Devices[END_REF].

Peter E. Raad, Pavel L. Komarov, and

METHODOLOGY

The overall approach combines computational and experimental methods previously developed by the authors. The transient two-dimensional surface temperature is measured by the use of the thermoreflectance thermography system [START_REF] Komarov | Thermal Characterization of Pulse-Activated Microelectronic Devices by Thermoreflectance-Based Surface Temperature Scanning[END_REF], while the three-dimensional thermal behavior of multi-layered integrated circuits (with embedded features) is inferred by solving the inverse heat transfer problem with the selfadaptive ultra-fast numerical technique [START_REF] Raad | System and Method for Predicting the Behavior of a Component[END_REF][START_REF] Wilson | A Transient Self-Adaptive Technique for Modeling Thermal Problems with Large Variations in Physical Scales[END_REF]. To minimize the number of uncertainty sources that the inverse method must deal with, the thermo-physical properties of the various thin layers can be measured independently with the Transient Thermo-Reflectance technique [START_REF] Burzo | A Study of the Effect of Surface Metallization on Thermal Conductivity Measurements by the Transient Thermo-Reflectance Method[END_REF][START_REF] Komarov | Performance Analysis of the Transient Thermo-Reflectance Method for Measuring the Thermal Conductivity of Single Layer Materials[END_REF].

Thermo-Reflectance Thermography System

The experimental temperature mapping system is based on the thermoreflectance (TR) method, where the change in the surface temperature is measured by detecting the change in the reflectivity of the sample. The measurement methodology requires two steps. First, the thermoreflectance coefficient must be determined for each of the surface materials in the measurement area. Second, the changes in the surface reflectivity as a function of changes in temperature are measured over the area of interest, with submicron spatial resolution. The resulting reflectivity data is combined to obtain a temperature field over that area of interest. The thermoreflectance coefficient, C TR , varies as a function of the material under test, material layering, and the wavelength of the probing light [START_REF] Ujihara | Reflectivity of Metals at high Temperatures[END_REF]. Thus, in order to maximize the signal to noise ratio, it is important to use a light source whose wavelength produces the maximum value of C TR for the exposed material layer.

A schematic of the TR thermography (TRTG) system is shown in Fig. 1. The probing light reflects from the heated surface back along the optical path to the sensitive element of a CCD camera (512×512 pixels). The intensity of the reflected light depends on the reflectivity (temperature) of the sample's surface. The frames containing the change in surface reflectivity induced by the temperature variations of the DUT are acquired and scaled according to the calibrated data. The details of the TRTG method are presented in a companion article.

The calibration approach consists of determining the relationship between the changes in reflectance and surface temperature. The change in reflectance is measured by a differential scheme involving two identical PDs in order to minimize the influence of fluctuations in the energy output of the probing laser. The sample temperature is controlled by a thermoelectric (TE) element and measured with a thermocouple. The calibration must be performed for each of the materials on the surface of each device where a temperature mapping is carried out.

Computational Engine

The numerical engine is capable of simulating the transient thermal behavior of active multi-layered devices whose dimensions vary over several orders of magnitude and where the thermophysical properties of the materials used may not be isotropic. The thermal modeling engine is used for solving the required heat transfer problem of the corresponding physical device. The measured surface temperature field (outlined white square in Fig. 2) is then used as input signature for an optimization scheme that varies control parameters (e.g., source power, length) until the RMS error between the computed solution and the input signature is minimized.

The novel approach begins by solving the corresponding steady-state problem by the use of a grid nesting technique. Since the physical dimensions of the various materials used in modeling high performance electronic devices vary greatly, a uniform mesh that resolves all of the details in three dimensions results in a prohibitively large computational grid. A common method for dealing with dimensional variation is to skew the mesh and concentrate more grid points in areas where higher resolutions are needed. The shortcoming of using a biased-mesh approach to resolve the geometry is that the problem geometry, and not the temperature gradients, will end up dictating the meshing. The meshing strategy used in the development of the numerical engine [START_REF] Raad | System and Method for Predicting the Behavior of a Component[END_REF] was set on ensuring that it is (i) automatic and adaptive, (ii) independent of user expertise, and (iii) independent of materials (including air), geometry features, embedded vias, and heat source locations. The approach makes it possible to start with the full 3D geometry of a real device in all of its complexity, and then uses a physics-based automatic error predictor to focus the entire available computational power on only those regions that require further refinement in order to achieve the level of acceptable error prescribed by the analyst [START_REF] Wilson | A Transient Self-Adaptive Technique for Modeling Thermal Problems with Large Variations in Physical Scales[END_REF]. The power of the method is that it uses effective thermal properties that are consistent with the local grid spacing at the particular grid level in use. As a result, dealing with air, embedded vias, and ultra-thin multi-layered structures requires no special treatment.

Sample Characterization: Geometry and Thermal Properties

The fidelity of a computational solution depends directly on the accuracy of the material properties and geometric characteristics of the features making up the device of interest. Since it has been shown [START_REF] Burzo | Thermal Transport Properties of Gold-Covered Thin-Film Silicon Dioxide[END_REF] that the thermal properties of thin-films vary from those of bulk materials, it is necessary for the numerical simulation to use the real values of the thermal conductivity for all of the materials making up the system under study. The previously described TTR measurement system [START_REF] Burzo | A Study of the Effect of Surface Metallization on Thermal Conductivity Measurements by the Transient Thermo-Reflectance Method[END_REF] can be used to determine any unknown properties of thin-film materials and their interface resistances.

In order to carefully characterize the geometry of a device of interest, an ellipsometer was used to measure the thicknesses of transparent layers as well as to confirm the optical properties of surface metals. In addition, a profiler was used to measure the thicknesses of various opaque and transparent layers making up the device.

To prove the concept and validate the described method, we constructed basic reference aluminum microresistor devices buried in a layer of silicon dioxide, and investigated the accuracy with which the coupled experimental-computational approach is capable of determining key geometric features of a known device. In the micro-resistor device, shown schematically in Fig. 2, an aluminum (Al) strip heater is sandwiched between top and bottom oxide layers in the vertical direction, and between two Al activation pads in the horizontal direction. When activated at known electrical power levels, the known heat source will generate a threedimensional temperature field throughout the device. Since the approach calls for measuring the temperature signature on the surface of the device, a gold (Au) layer was deposited over the anticipated area of interest, which includes the heat source and large surrounding regions. Gold was chosen since when used in conjunction with an LED light source at 485 nm, it will maximize the C TR value in the TRTG technique.

The simple construction of this micro-resistor device makes it possible to specify and measure all essential heat transfer problem parameters. Specifically, (i) the geometry of the different layers can be controlled in the fabrication process and later measured for confirmation; (ii) the oxide layers and Al strip provide a Joule heat source with known uniform power distribution; and (iii) the large pads make it possible to use a four-wire scheme to simultaneously activate the strip heater and measure its electrical power. For the purposes of this article, only a single device is reported on, which has a width of 14 µm and a length of 200 µm. All other pertinent geometric parameters are provided in Fig. 2.

RESULTS AND DISCUSSION

By using the experimental system and combining it with the numerical approach presented above it becomes possible to solve the inverse conduction problem associated with a complex, multi-layered, deep submicron electronic device in order to infer the thermal behavior of the embedded features that cannot be otherwise accessed. The inverse solution is obtained by varying key parameters that define the heat transfer problem under consideration. For the test micro-resistor in this work, these parameters include the size and location of the heater strip, its power and distribution, the thicknesses of the top and bottom layers, and their thermal properties. Obviously, a solution that comprehends all of these variables would be impractical. Furthermore, many of these parameters can be either directly measured (e.g., thermal properties, applied power, thicknesses of layers) and/or have smaller influences on the final temperature Since heat flows primarily toward the substrate, the thickness of the bottom oxide layer is expected to strongly affect the temperature distribution because of the resistance that this layer presents. Therefore this parameter was chosen for its sensitivity on the final result. A second parameter was chosen to be the length of the heat source because (i) it cannot be directly measured and (ii) it is expected to be longer than the length of the heater strip itself because of the end effects at the junction between the aluminum heater strip and its activation pads. In fact, one would expect that the end effects would extend beyond the strip by a distance that scales with the heater's width.

Figure 3 compares the experimental data along the mid-plane of the heater strip to the corresponding numerically computed temperature distributions for different values of the bottom oxide thickness and for a fixed heat source length of 200µm. As expected, thicker layers of bottom oxide result in higher strip temperatures. It appears that the numerical temperature distribution for the thickness of 1,500Å fits nicely the measured data in the center of the heater, X∈[400,600], but misses inside the pad areas. The discrepancy is due to the end effect previously mentioned.

To examine the end effects, the bottom oxide thickness was held fixed at 2000Å, and the heat source length was varied as depicted by the results in Fig. 4. Indeed, extending the heat source length beyond the ends of the strip heater pulls the temperature curves outwardly toward the experimental data at the lower temperatures. However, the agreement between the numerical and experimental distributions becomes worse toward the center of the domain. It should be pointed out that simply extending the rectangular heat source, as done here, only captures the spherical nature of the heat distribution at the pad-heater junction to a first order approximation. In order to more precisely simulate the end effects, one would have to introduce a more sophisticated heat power distribution model in this region.

Nevertheless, taken together, the results of Figs. 3 and4 point clearly to the importance of both physical parameters and hence the need to optimize over both of them simultaneously. The optimization method used in this work is a variant of the "steepest descent" method [START_REF] Snyman | Practical Mathematical Optimization[END_REF]. The search begins with a numerical solution of the heat transfer problem at nominal values of the two parameters being considered, i.e., bottom oxide thickness (h) and heat source length (L). This numerical solution is compared with the measured signature field (outlined white square in Fig. 2) on the gold pad at every common location to compute an RMS error. Then, additional trial solutions and associated errors are obtained at neighboring pair values of h and L in the twodimensional parameter space. Evaluation of the errors at the nominal and eight neighboring pairs provides the direction for modifying h and L in order to decrease the error. This process continues until the changes in both h and L are acceptably small.

For the specific device under consideration, Fig. 5 shows the path taken in the (h, L) parameter space to converge onto the final values of h and L which yield a numerical solution that matches the measured surface signature to within the minimum RMS error. Of particular interest here is that the process of reaching the final result required 57 solutions of the full three-dimensional heat transfer problem, each of which was converged to less than 1% numerical error at the specified h and L values. It is obvious that such an approach would be impractical To compare the experimental and numerical surface temperature fields, a surface slice is extracted from the full 3D solution that corresponds to the size, location, and resolution of the experimental area. Figure 6(a) shows the surface temperature slice at stage 9, which represents the "optimal" solution that is closest to the experimental signature shown in Fig. 6(b). While the agreement is very favorable, it is clear that the end effects give the experimental contours a more "rectangular" shape on the heater edges. As discussed above, a more sophisticated model of the power distribution at the pad-heater junction would be required to further improve the agreement, which is beyond the scope of this validation article.

The ultimate benefit of the coupled experimental-numerical system is that it takes an inherently 2D experimental approach and provides a full 3D thermal characterization of the complete device. The optimization-based coupling ensures that the 3D solution is consistent with the 2D experimental signature within the chosen heat transfer model as defined by the combination of known and unknown (i.e., optimizable) parameters. Figure 7 provides an example of the optimal 3D numerical solution. For clarity, only two slices are shown and the spatial domain is restricted to the region surrounding the heater, even though the computational domain is 1000×500×1000 µm. The horizontal slice is on the surface of the device and the vertical slice cuts across the mid-plane of the heater.

CONCLUSIONS

This work builds on the previous introduction [START_REF] Raad | Coupling Surface Temperature Scanning and Ultra-Fast Adaptive Computing to Thermally Fully Characterize Complex Three-Dimensional Electronic Devices[END_REF] of a coupled experimental-computational system that makes it possible to extend an experimentally-obtained 2D surface field to a full 3D characterization of the thermal behavior of complex submicron electronic devices. The new system improves on the previous one by replacing the laser-based, point-by-point, surface temperature scanning approach with a CCD camera-based approach. The aim of this work was to validate the coupled system by applying it to a simple test micro-resistor device with known geometric, material, and thermal characteristics. The investigation focused on the two parameters that most influence the heat transfer solution within the device, which are the thickness of the bottom oxide (h) and the length of the heat source (L). The obtained optimal value of h was found to be within 4% of the physical thickness as measured with an ellipsometer during the fabrication process. The optimal heat source length was found to extend into the side pads by a distance that is approximately equal to half the heater's width, which is qualitatively consistent with the physics of Joule heating in the sudden expansion areas connecting the strip and the pads. The results herein have provided the desired validation for the coupled experimentalcomputational system as well as demonstrating the power of the method in providing full 3D temperature field that agrees very well with the measured 2D surface temperature signature.

For existing devices, the highly resolved and accurate 3D temperature field would provide the ability to detect hot spots, diagnose performance, and assess reliability. In the design and manufacturing of new devices, this new tool has the potential to provide a rapid approach for analyzing the thermal behavior of complex stacked structures, to identify regions of excessive heat densities, and ultimately to contribute to improved thermal designs, better device reliability, and shorter design cycle time.

Fig. 1

 1 Fig. 1 CCD-based TRTG scanning system

Fig. 4

 4 Fig. 3 Influence of the bottom oxide thickness

Fig. 6 Fig. 7

 67 Fig. 5 Convergence of optimization path toward final solution at h = 1,413 Å and L = 212.5 µm

  Mihai G. Burzo A Coupled Thermoreflectance Thermography System and Ultra-Fast Computational Engine...

©TIMA Editions/THERMINIC 2006

-page-ISBN: 2-916187-04-9

  Peter E. Raad, Pavel L. Komarov, and Mihai G. Burzo A Coupled Thermoreflectance Thermography System and Ultra-Fast Computational Engine...

	Probing Hole 400x300 µm	Al Heater Pad 500 µm wide
	Au Test Pad
	0.8x0.8 mm
	Al Heater, 14x200 µm
	Al Heater and Pads	Bottom SiO2 Layer
	1000 Å	1360 Å
	Au Test Pad, 2000 Å
	Top SiO2 Layer, 3200 Å
	Si Substrate

Fig. 2 Geometry of test micro-resistor ©TIMA Editions/THERMINIC 2006 -page-ISBN: 2-916187-04-9 distribution in the device (e.g., layers of SiO 2 and Au above the heater).

ACKNOWLEDGEMENTS

We wish to thank Mr. Jay Kirk of the EE Department at SMU for his help in fabricating the devices for this work.


Tetsuya Baba

Kazuko Ishikawa

Takashi Yagi

Naoyuki Taketoshi

N Taketoshi Measurements Of

Thermophysical Prop

Kimiaki Tamano

Tetsuro Ohtsuka

Hiroshi Watanabe

Yuzo Shigesato

MEASUREMENTS OF THERMOPHYSICAL PROPERTY OF THIN FILMS BY LIGHT PULSE HEATING THERMOREFLECTANCE METHODS

whether they are published or not. The documents may come

INTRODUCTION

Phase-change optical disk media (DVD Rewritable), high-density integrated circuit and flat panel display are comprised of several nanometers to several 100 nanometers thick of thin films. To know how the heat flows and what distribution of temperature is caused when they are used, information of the thermophysical properties of thin films and the boundary thermal resistance between thin films is required [START_REF] Baba | General needs on nanoscale thermal metrology and the Japanese program on this subject[END_REF][START_REF] Cahill | Thermometry and thermal transport in micro/nanoscale solid devices and structures[END_REF].

Picosecond thermoreflectance method was developed to measure thermal diffusivity of subnanometer thick thin films by Paddock and Eesley [START_REF] Paddock | Transient thermoreflectance from thin metal films[END_REF]. The optical reflection intensity of the temperature detection light is detected by photodiode. Since reflectivity of material surface changes dependent on the surface temperature, the change of specimen front face temperature can be observed by the change of reflected light amplitude. This temperature measurement method with the temperature change of such a reflectivity is called as thermoreflectance method [START_REF] Rosencwaig | Detection of thermal waves through optical reflectance[END_REF]. Thermal diffusivity of submicrometer thin films perpendicular to the surface was calculated from the cooling rate of the surface temperature and the penetration depth of the heating light.

National Metrology Institute of JAPAN, AIST has succeeded in developing the thermoreflectance methods by picosecond / nanosecond light pulse heating [START_REF] Taketoshi | Thermal diffusivity measurement of metal thin films with a picosecond thermoreflectance technique[END_REF][START_REF] Taketoshi | Picosecond thermoreflectance measurements of thermal diffusion in film/substrate two layer systems[END_REF][START_REF] Baba | Thermophysical property measurement by light pulse heating[END_REF] and realized to measure the thermal diffusivity of metallic thin films from several 10 nm to several micrometers thick on transparent substrate in thickness direction under the configuration of rear face heating / front face detection picosecond thermoreflectance method [START_REF] Baba | General needs on nanoscale thermal metrology and the Japanese program on this subject[END_REF][START_REF] Taketoshi | Observation of heat diffusion across submicrometer thin films using a picosecond thermoreflectance technique[END_REF][START_REF] Taketoshi | Development of a thermal diffusivity measurement system for metal thin films using a picosecond thermoreflectance technique[END_REF][START_REF] Taketoshi | A new picosecond thermoreflectance technique for thermal diffusivity measurements of nanoscale metal thin films[END_REF][START_REF] Taketoshi | Electric delay technique in a picosecond thermoreflectance method for thermophysical property measurements of thin films[END_REF].

Since the geometrical configuration of this method is the same as the laser flash method which is the standard measurement technique for the thermal diffusivity of bulk materials [START_REF] Baba | Improvement of the laser flash method to reduce uncertainty in thermal diffusivity measurements[END_REF][START_REF] Akoshima | Laser flash thermal diffusivity measurements of isotropic graphite and glass-like carbon[END_REF], thermal diffusivity value can be calculated reliably from heat diffusion time across well-defined length of the film thickness under one-dimensional heat flow [START_REF] Taketoshi | Observation of heat diffusion across submicrometer thin films using a picosecond thermoreflectance technique[END_REF][START_REF] Taketoshi | Development of a thermal diffusivity measurement system for metal thin films using a picosecond thermoreflectance technique[END_REF].

HIGH SPEEDLIGHT PULSE HEATING THERMOREFLECTANCE METHODS

Front face heating / front face detection

It was not easy by conventional measurement technique to determine the thermal conductivity and the thermal diffusivity of thickness direction of thin films of less than one micrometer thick. In order to solve this problem, picosecond thermoreflectance method was developed to observe the temperature changes of thin film front face by heat diffusion to the inside [START_REF] Paddock | Transient thermoreflectance from thin metal films[END_REF]. The optical reflection intensity of the temperature detection light is detected by photodiode based on the thermoreflectance method.

In this thin film thermal diffusivity measurement system by a picosecond thermoreflectance method, the laser beam emitted from a mode lock titanium sapphire laser and is divided into transmitted beam and reflected beam by a quartz plate. About 90 % is used for pulse heating and the other about 10 % is used for temperature detection to measure the temperature changes of the thin film front face. Light travels 0.3 mm in one picosecond. By adjusting distance to a specimen after it was divided, the time difference that a heating light and detection light arrive at the specimen front face can be controlled. The response time of the thermoreflectance method is much faster than thermocouples, resistance temperature sensor or radiation thermometers. According to the pump probe method, ultra fast thermometry is possible only limited by time duration of the pulses. On the other hand, it is a weak point that the sensibility of temperature detection is low.

Figure 1 Heating area in front face heating / front face detection picosecond thermoreflectance method Figure 2 shows the result which observed the change of the front face temperature by the picosecond thermoreflectance method about three kinds of aluminium thin films of different thickness synthesized on a glass substrate [START_REF] Taketoshi | Picosecond thermoreflectance measurements of thermal diffusion in film/substrate two layer systems[END_REF].

For the specimen of 500 nm thick, the heat has not arrived at the substrate within 120 ps after pulse heating in and the front face temperature change represented by red line corresponds to internal heat diffusion of the aluminium thin film.

On the other hand, for the specimen of 100 nm thick, the temperature change speed decreases around 30 ps after pulse heating as show by black line and deviates from the temperature change of the 500 nm thick specimen.

Because the thermal diffusivity of glass substrate is much smaller than the thermal diffusivity of aluminium thin film, the heat diffusion to the substrate is suppressed when the heat arrives at the interface between the thin film and the substrate.

As shown by blue line in figure 2, for the specimen of 50 nm thick, the temperature change only for the inside of the thin film cannot be observed because of heat diffusion to the substrate just after pulse heating. Area of a diameter of several 10 µm on thin film front face is heated by the picosecond laser beam and the same position is irradiated by the probe beam. Then, the history of front face temperature is observed by the conventional thermoreflectance method.

In this method, the thermal diffusivity can be calculated from the cooling rate after pulse heating. However, it is not easy to make quantitative and reproducible measurement because the cooling rate changes sensitively dependent on the condition of thin film front face.

Rear face heating / front face detection

National Metrology Institute of JAPAN, AIST has developed rear face heating / front face detection picosecond light heating thermoreflectance methods which are evolution of the conventional laser flash method and the picosecond thermoreflectance method as shown in figure3 [START_REF] Taketoshi | Observation of heat diffusion across submicrometer thin films using a picosecond thermoreflectance technique[END_REF][START_REF] Taketoshi | Development of a thermal diffusivity measurement system for metal thin films using a picosecond thermoreflectance technique[END_REF][START_REF] Taketoshi | A new picosecond thermoreflectance technique for thermal diffusivity measurements of nanoscale metal thin films[END_REF][START_REF] Taketoshi | Electric delay technique in a picosecond thermoreflectance method for thermophysical property measurements of thin films[END_REF]. Figure 4 shows the block diagram of a measurement system. This configuration is essentially equivalent to the laser flash method which is the standard measurement method to measure the thermal diffusivity of bulk materials. The thermal diffusivity of the thin films can be calculated with small uncertainty from the thickness of a thin film and the heat diffusion time across a thin film. The temperature detection beam is focused to a spot diameter of about 50 µm at the specimen front face just opposite to the heating light focus position. The reflected light intensity of the temperature detection light is in proportion to the change of specimen front face temperature and the change of the reflected light intensity is detected by photodiode. The temperature coefficient of reflectivity for normal metal such as aluminium is small with 10 -5 /K order. Since the transient temperature rise of the specimen front face after picosecond pulse heating is smaller than a few degrees, the thermoreflectance signal is much smaller than the offset level of the reflected light. Such a small signal can be measured by lock-in detection at modulation frequency of heating light by an acoustic optic modulator.

Delay time (ps)

Figure 5 shows the temperature history curves of an aluminium single-layered thin film of 100 nm thick and a molybdenum single-layered thin film of 100 nm measured by the picosecond light pulse heating thermoreflectance method [START_REF] Taketoshi | Picosecond thermoreflectance measurements of thermal diffusion in film/substrate two layer systems[END_REF]. Both films were synthesized on a Pyrex glass substrate by magnetron DC sputtering method. Here, film thicknesses are nominal values.

Since thermoreflectance signals are similar to those observed by the laser flash method for bulk specimen suggest, the heat energy transport of these metallic thin films of about 100 nm thick at room temperature in time scale of several 100 ps can be expressed by the classic thermal diffusion equation .

Nanosecond light pulse heating thermoreflectance method

Initially, the electrical delay method was developed to expand the observation time of the picosecond light pulse heating thermoreflectance method longer for measurements of thicker films [START_REF] Taketoshi | Electric delay technique in a picosecond thermoreflectance method for thermophysical property measurements of thin films[END_REF]. Since the pulse duration can be longer than picosecond and the repetition period of pulses is flexible, nanosecond pulse laser can be used for light pulse heating thermoreflectance method [START_REF] Baba | General needs on nanoscale thermal metrology and the Japanese program on this subject[END_REF][START_REF] Baba | Development of a high speed laser flash system to measure thermophysical properties of thin films, Nanosecond optical pulse heating thermoreflectance method[END_REF]. For the pump pulse, pulse duration is 2 ns and pulse interval is 20 µs with intensity modulation of 1 kHz by an acoustic optical modulator. Typical size and shape of the specimen are a disk of 10 mm in diameter or a square of 10 mm on each side. The specimen is irradiated from bottom to the rear face of the specimen by the heating beam and the reflected light of the probe beam is detected by a photodiode. The thermoreflectance signals are detected using a lock-in amplifier. In this system it is possible to measure the thermal diffusivity of thin films with thickness up to several micrometers.

MEASREMENTS OF THIN FILMS FOR INDUSTRIAL USE

Transparent conductive films

The thermal diffusivity of tin doped indium oxide (ITO) thin films with the thickness of 200 nm and 300 nm have been measured using the nanosecond light pulse heating thermoreflectance method [START_REF] Tamano | Thermal diffusivity measurements of transparent conductive thin films using a nanosecond thermoreflectance technique[END_REF]. As shown in figure 6, ITO films between molybdenum (Mo) films of 70 nm thick were prepared on fused silica substrate by DC magnetron sputtering using ITO and Mo multi-targets because the wavelengths of pulse lasers used in this study are 1064 nm as pump beam and 830 nm as probe beam where ITO is transparent. The thermal diffusivity measurements of three-layered films were carried out using the nanosecond light pulse heating thermoreflectance system. During deposition of Mo films, total gas pressure was maintained at 1.0 Pa. On the other hands, ITO layers were deposited under various total gas pressure between 0.5 and 3.0 Pa. Such a Mo/ITO/Mo layered structure was fabricated without exposure to the atmosphere between each deposition. Analysis of heat diffusion across three layer thin films on substrate must be considered in order to calculate the thermal diffusivity of the thin film between metal thin films. It is also necessary to know the boundary thermal resistance between the layers as well as the thermophysical properties of each layer. The measured temperature history curve for the three-layered films were analyzed by the response function method which is general technique to analyze heat diffusion across a multilayer films developed by AIST [START_REF] Baba | General needs on nanoscale thermal metrology and the Japanese program on this subject[END_REF][START_REF] Baba | Thermophysical property measurement by light pulse heating[END_REF][START_REF] Baba: Jpn | Introduction of a response function method for analysis of functionally gradient materials[END_REF][START_REF] Baba | Analysis of heat diffusion in multi-layer thin films by a response function method[END_REF][START_REF] Baba | Analysis of heat diffusion across multi-layer thin films and boundary thermal resistances by the response function method general solution for onedimensional heat diffusion after impulse heating[END_REF]. The areal heat diffusion time is defined as the area surrounded by the horizontal line at the height of the maximum temperature rise and the transient temperature response curve at the rear face after pulse heating as shown in figure 7 [START_REF] Baba | General needs on nanoscale thermal metrology and the Japanese program on this subject[END_REF]. Thermal diffusivity value of the thin unknown thin film and the boundary thermal resistance can be calculated from the observed areal heat diffusion time based on the response function method assuming thickness of each film, thermal diffusivity and volume heat capacity of molybdenum thin films are known. The thermal diffusivity of the ITO films decreased as the total gas pressure increased as shown in figure 8 [START_REF] Tamano | Thermal diffusivity measurements of transparent conductive thin films using a nanosecond thermoreflectance technique[END_REF].

Hard coating films

Titanium Nitride (TiN) thin films are used as hard coating films. TiN thin films with the thickness from 200 nm to 800 nm were synthesized on glass substrate by reactive RF magnetron sputtering using Ti target and Ar-N2 plasma. Thermal diffusivity measurements of these films were carried out using the nanosecond light pulse heating thermoreflectance system [START_REF] Ohtsuka | Thermal diffusivity measurement of TiN x thin films using a nanosecond thermoreflectance technique[END_REF]. As shown in figure 9, the thermal diffusivity of TiN film changes dependent on the ratio of N over Ti which is controlled by N 2 gas flow ratio in Ar-N 2 . There is good linear relation between the thermal diffusivity and the electric conductivity as shown in figure 10. The thermal conductivity was calculated using the specific heat capacity and the density of bulk TiN. Then, Lorenz number of TiN thin film was calculated assuming Wiedemann-Franz law. As electrical conductivity increases, Lorenz number converges to Sommerfeld value based on the free electron model. This result means that the major heat carriers are phonons for lower electrical conductive TiN and electrons for higher electrical conductive TiN [START_REF] Ohtsuka | Thermal diffusivity measurement of TiN x thin films using a nanosecond thermoreflectance technique[END_REF].

Multilayered films of the next generation phase change optical disk

Dependence of the thermal diffusivity of Ge 2 Sb 2 Te 5 thin films, used for optical recording media, on phase change has been investigated using the nanosecond light pulse heating thermoreflectance system [START_REF] Watanabe | Thermal diffusivity measurements of Ge 2 Sb 2 Te 5 thin films using a nanosecond thermoreflectance technique[END_REF]. Two kinds of thin films are prepared; one is Ge 2 Sb 2 Te 5 single layered thin film and the other is Mo/Ge 2 Sb 2 Te 5 /Mo three-layered thin film, where 70 nm thick Mo films act as the absorbing layer for the heating beam and the reflection layer for the temperature detection beam. Films, with which the thickness of Ge 2 Sb 2 Te 5 is from 100 nm to 400 nm, were deposited on non-alkali glass substrates by RF magnetron sputtering from Ge 2 Sb 2 Te 5 and Mo targets. It is observed by X-ray diffraction study that the as-deposited films show amorphous structure and it is transformed to FCC crystal phase by heat treatment of 5 min at 573 K. The thermal diffusivity of Ge 2 Sb 2 Te 5 amorphous thin film is determined to be 2.5×10 -7 m 2 /s and the thermal diffusivity of the crystal film is 4.8×10 -7 m 2 /s from the temperature history curves shown in figure 11 observed by the nanosecond light pulse heating thermoreflectance system [START_REF] Watanabe | Thermal diffusivity measurements of Ge 2 Sb 2 Te 5 thin films using a nanosecond thermoreflectance technique[END_REF].

SUMMARY

Thermoreflectance methods by picosecond light pulse heating and by nanosecond light pulse heating have been developed under the same configuration as the laser flash method which is the standard method to measure the thermal diffusivity of bulk materials. Thermal diffusivity values across thin films were measured with small uncertainty. These high speed light pulse heating thermoreflectance methods can be observed the heat diffusion time across well-defined length of the film thickness under one-dimensional heat flow. Using these light pulse heating methods, the thermal diffusivity of each layer of multilayered thin films and the The thermophysical properties of the transparent conductive films for flat panel display and the hard coating films were measured with the high speed light pulse heating thermoreflectance methods. The boundary thermal resistance between thin films and the thermal diffusivity of each layer of multilayer films constitute next generation phase-change optical disc were also measured with these methods. Reliable thermal design can be realized by heat conduction simulation using these reliable thermal diffusivity values of thin films and boundary thermal resistance values between thin films measured with the high speed light pulse heating thermoreflectance methods.

Figure 2

 2 Figure 2 Thermoreflectance signals of three kinds of different thickness of aluminium thin film that synthesized on a glass substrate

Figure 3

 3 Figure 3 Principle of rear face heating / front face detection picosecond light heating thermoreflectance method

Figure 4

 4 Figure 4 Block diagram of the thermoreflectance method thin film thermophysical property measuring system by picosecond pulsed laser heating

Figure 5

 5 Figure 5 Thermoreflectance signals of aluminium thin film and molybdenum thin film of 100 nm thick by rear face heating / front face detection picosecond thermoreflectance method

Figure 7 Figure 8 Figure 6

 786 Figure 7 Area surrounded by the maximum temperature rise line and the temperature response at the specimen rear face after the pulse heating

Figure 9 Figure 10

 910 Figure 9 Thermoreflectance signals of TiN thin films of different ratio of N over Ti

  Tetsuya Baba, Kazuko Ishikawa, Takashi Yagi, Naoyuki Taketoshi MEASUREMENTS OF THERMOPHYSICAL PROPERTY OF THIN FILMS ©TIMA Editions/THERMINIC 2006 -page-ISBN: 2-916187-04-9 boundary thermal resistance between the layers can be determined from the observed transient temperature curves based on the response function method.

Figure 11

 11 Figure 11 Thermoreflectance signals of Mo/Ge 2 Sb 2 Te 5 / Mo three layered thin films in FCC structure and amorphous structure

©TIMA Editions/THERMINIC 2006-page-ISBN: 2-916187-04-9


Pierre-Olivier Chapuis

Kumar Sourabh

Sebastian Saha

Volz

QUANTITATIVE 3 ΟΜΕGΑ SCANNING THERMAL MICROSCOPY : MODELING THE TIP AC/DC COUPLING AND THE SAMPLE HEAT CONDUCTION

A way to increase the Scanning Thermal Microscope (SThM) sensitivity in the harmonic 3ω mode is to heat the probe with an AC current sufficiently high to generate a coupling between the AC and the DC signals. We detail in this paper how to properly take into account this coupling with a Wollaston-probe SThM. We also show how to link correctly the thermal conductivity to the thermal conductance measured by the SThM.

INTRODUCTION

The continuous decrease of sizes in the material science applications and in the microelectronics industry has led to the development of instruments that aim at analyze heat transfer at micro and submicroscales. Thermal characterization at nanoscales has recieved significant amount of efforts since a decade. While the optical methods have reached their maturity and remain limited by diffraction to half a micron, contact probes were considered as the only means to reach 10-100nm scale resolutions. The Scanning Thermal Microscope (SThM) based on an Atomic Force Microscope mounted with thermal probe was invented in 1986 by Williams and Wickramasinghe [START_REF] Williams | Scanning Thermal Profiler[END_REF] to provide the topography of electrically insulating materials. Since then, various techniques have been proposed with spatial resolutions varying between 5 microns and a few tens of nanometers. However, quantitative measurements are still affected by two main artefacts, i.e. (i) the topography-thermal coupling and (ii) the thermal resistance between tip and sample. We previously showed that the artefact (i) is amplified by heat conduction through air due to poor spatial resolution [START_REF] Lefèvre | 3-omega Scanning Thermal Microscope[END_REF]. The thermal resistance between the tip and the sample (ii) is due to conduction in air, through the solid-solid contact and through the water meniscus between tip and sample. The mechanical properties, the roughness and the hydrophilicity of the sample are involved in those interactions.

We will focus on the Wollaston wire probe [START_REF] Dinwiddie | Thermal conductivity contrast imaging with a scanning thermal microscope[END_REF][3][START_REF] Lefèvre | Modélisation et élaboration des métrologie de microscopie thermique à sonde locale résistive[END_REF]. A way to increase the SThM sensitivity in the harmonic (3ω) mode [START_REF] Lefèvre | 3-omega Scanning Thermal Microscope[END_REF] is to heat the sample with an AC current sufficiently high to generate a coupling between the AC and the DC signals. We detail in this paper how to properly take into account this coupling. We also give the correct factor that takesnto account the constriction of the heat flux lines in the sample.

PROBE MODEL

The probe is a Wollaston wire consisting of a platinum core 5 microns in diameter and a silver coating 70 microns in diameter. The silver coating is etched to uncover the platinum wire over a 2Ls=200 microns length. This tip was studied in several of our previous works [3][START_REF] Lefèvre | Modélisation et élaboration des métrologie de microscopie thermique à sonde locale résistive[END_REF][START_REF] Lefèvre | 3-omega Scanning Thermal Microscope[END_REF]. A modulated (AC) electrical current is used to Joule heat the wire and the second harmonic of the temperature is measured with a Lock In Aplifier at the third harmonic of the voltage (see Reference [START_REF] Lefèvre | 3-omega Scanning Thermal Microscope[END_REF]).

The previous papers did not take the coupling between the constant (DC) and the modulated (2ω) temperature into account. The AC-DC coupling is solved in this section when the tip is in contact and out of contact with the sample. An estimation of the Biot number shows that heat conduction can be considered as monodimensional in the wire [START_REF] Lefèvre | Modélisation et élaboration des métrologie de microscopie thermique à sonde locale résistive[END_REF]. The size of the tip is micrometric, so we can use the classical Fourier conduction equation in the tip: where x denotes the abscissa in the probe (x=0 is at the silver-platinum interface and x=Ls is the middle of the tip). as and λ are the probe thermal diffusivity and conductivity. S and p correspond to the section and the

©TIMA Editions/THERMINIC 2006 -page- ISBN: 2-916187-04-9 2 2 2 2 ( ) 1 o s ph T T T T I a t x S S ρ λ λ - ∂ ∂ = - + ∂ ∂ 2 3 2 o o R I V ω ω α θ = P.
-O. Chapuis, S. K. Saha, and S. Volz, Quantitative 3ω Scanning Thermal Microscopy : Modeling the tip AC/DC coupling and the sample heat conduction perimeter of the wire probe. ρ is the tip electrical resistivity and h represents the heat transfer coefficient between the wire and the ambient. When the current is modulated:

Taking into account the dependance on temperature of the probe resistivity, it turns out that where α is the temperature coefficient of the wire. The first equation can then be written as follows:

if we consider a time origin so that the phase is zero. We consider that the local temperature in the tip can be written as the sum of a DC term and a AC term:

where Note that the AC term is complex : a phase lag between the current and the intensity is possible. Thus we can decompose the equation for θ into its DC component:

and into its 2ω harmonic component:

The The silver end of the probe is a thermal sink for the harmonic component. For the DC component, the flux at the silver end is taken into account through the high conductance of one of the legs of the silver coating [START_REF] Lefèvre | Modélisation et élaboration des métrologie de microscopie thermique à sonde locale résistive[END_REF].

The other end can be considered insulated due to symmetry when the tip is out of contact with the sample. The boundary conditions are described as follows. At x=0:

At x=Ls :

When the tip is in contact, the boundary conditions are modified by the heat flux from the tip to the sample. We assume that this flux is confined in the volume around the tip/sample contact. The heat transfer through air is neglected on the remaining surface of the tip. We argue that the hotter part of the tip is in the vicinity of the sample. This was confirmed by simulations in a previous paper [START_REF] Lefèvre | Nanoscale Heat Transfer at Contact Between a Hot Tip and a Substrate[END_REF]. When the pressure decreases, the mean free path of the air molecules increases, for instance 10 -2 bar corresponds to a 10 microns mean free path. But the surface of the tip legs have a lower contribution because 1) (i) the temperature is lower on these legs than in the middle of the tip and (ii) also note that the moleculemolecule scattering is still lower in the heat transfer between the very end of the tip and the sample than the one between the tip legs and the contact area. We finally write at x=Ls:

©TIMA Editions/THERMINIC 2006 -page- ISBN: 2-916187-04-9 2 2 2 2 ( ) (1 )(1 cos(2 )) 1 2 o o o s ph T T I T t T T a t x S S ρ αδ ω λ λ - + + ∂ ∂ = - + ∂ ∂ 2 2 2 2 2 2 2 2 2 ( 1) 2 o o D C d ph I j a dx S S ω ω ω ω θ θ ρ ω θ αθ αθ λ λ = - + + + 2 2 1 2 2 o o I ph m S S ρ α λ λ = - 2 2 2 o o I J S ρ λ =
= + + + + 2 2 2 2 2 2 o o s I ph m j S S a ρ α ω λ λ = - + 1 2 2 1 2 2 2 2 1 2 2 1 2 2 (
Jc A m m Jc B m m J J m C m α α α = - - = - - + = 2 0, 0 D C d d dx dx ω θ θ = = 2 , 0 D C A g D C d S G dx ω θ λ θ θ = = cos( ) o I I I t ω ϕ = + (1 ) o T ρ ρ αδ = +
where S is the contact surface and the prime on Geq denotes that it is half of the total conductance since the equations are solved only for one half of the tip. We surely do not reproduce the real heat flux lines in the tipsample contact area. But the size of this area is small in comparison to the one of the tip. The average temperature will not be significantly modified. Note also that Geq,AC is not known yet. We assume that Geq,AC=Geq,DC=Geq because the thermal characteristic time is smaller than the reverse of the frequency f=3ω/2π. For heat transfer through air, the characteristic time is the diffusion time through air over 100 microns at atmospheric pressure, which is of the order of 10 -4 s. The heat transfer through the meniscus and the mechanical contact are faster because the system sizes are of the order of nanometers. This justifies the equivalence of the tipsample heat transfers in the DC mode and in the AC mode in our range of frequencies.

Applying the boundary conditions we have: and with G'eq=0 for the out of contact case. At this step, we know θ(x). What is measured in the 3ω mode is the mean temperature of the tip. It is found by integration over the length of the tip.

RESULTS

We present on the following figure the results showing that the coupling is responsible for the heating of the probe of a temperature on the order of 5 K for an input current of 50 mA at ambient pressure (here h=3500 Wm -1 K -1 ). If the coupling was not taken into account, a bad exploitation of the results could lead to an error of one order of magnitude in the conductance calculation. One should note as seen on Figure 2 that the coupling does not affect the phase lag : the two curves (with and without coupling) are completely superimposed.

HEAT FLUX IN THE SAMPLE

We have considered in the first part the heat flux at the contact by using a the thermal conductance Geq. Its physical meaning corresponds to the contributions of the contact, Gc, and the sample, Gs, conductances:

The contact conductance Gc takes into account conduction in air (Ga) and in the water meniscus (Gw) and through the solid-solid mechanical contact (Gss). We stress in the following that a summation rule is applicable Gc=Ga+Gw+Gss The sample conductance Gs corresponds to the effect of constriction of the heat flux lines in the semi-infinite sample. If we consider that the heat flux generated by the tip in the sample arises from a hot disk on the top of the sample (see Figure 3), the geometrical factor taking into account the constriction of the heat flux lines is 3π 2 /8 and we have then :

This relation is retrieved by performing a finite element simulation of a hot disk on the top of the semi-infinite volume. The geometrical factor which is independent of the thermal conductivity of the material matches with the analytical one within 10%. Note that this is the first time that this geometrical factor is correctly taken into account : there was previously [3][START_REF] Lefèvre | Modélisation et élaboration des métrologie de microscopie thermique à sonde locale résistive[END_REF][START_REF] Lefèvre | 3-omega Scanning Thermal Microscope[END_REF] a confusion with the geometrical factor of the lines arising from a semihemisphere (2π).

CONCLUSIONS

We presented here a new step on the road to quantitative Scanning Thermal Microscopy by showing how to take into account the AC/DC coupling at high input currents, which are necessary if one wants to increase the sensitivity of the instrument. We also showed how to link correctly the thermal conductivity and the heat flux flowing from the tip to the sample by introducing the correct geometrical factor.

  DC and AC equations are coupled through the dependence of the harmonic temperature on the DC temperature. The here-called DC equation is the same as the one obtained in the DC regime if the Joule term is divided by a factor of 2. The general solution of the DC equation is where and The general solution of the AC equation is: where and The constants c1,c2,c3 and c4 are obtained from the boundary conditions. The symmetry of the problem allows us to solve the problem for only one half of the probe tip.

  O. Chapuis, S. K. Saha, and S. Volz, Quantitative 3ω Scanning Thermal Microscopy : Modeling the tip AC/DC coupling and the sample heat conduction

Figure 1 .

 1 Figure 1. Harmonic temperature amplitude calculated as a function of the frequency of the input current : cases with and without AC/DC coupling.

Figure 2 .

 2 Figure 2. Phase lag of the harmonic temperature calculated as a function of the frequency of the input current : there is no difference with or without coupling

Figure 3 .

 3 Figure 3. Heat flux lines arising from the hot disk of radius b on the surface.