- Email: [email protected]
Phonon transport in silicon–silicon and silicon–diamond thin films: Consideration of thermal boundary resistance at interface
Physica B 406 (2011) 2186–2195
Contents lists available at ScienceDirect
Physica B journal homepage: www.elsevier.com/locate/physb
Phonon transport in silicon–silicon and silicon–diamond thin films: Consideration of thermal boundary resistance at interface S. Bin Mansoor, B.S. Yilbas n ME Department, King Fahd University of Petroleum and Minerals, Saudi Arabia
a r t i c l e i n f o
a b s t r a c t
Article history: Received 14 February 2011 Accepted 13 March 2011 Available online 30 March 2011
Phonon transport across the silicon–silicon and silicon–diamond dielectric films is examined. To simulate the phonon transport in the films EPRT is accounted and numerical solutions are obtained with the consideration of 1 K difference across the films prior to the initiation of the phonon transport. The diffuse mismatch model is introduced to account for the thermal boundary resistance at the interface of the films. To validate the code developed, the steady and transient cases for phonon transport in silicon film are carried out. It is found that the equilibrium temperature of phonons attains higher value at the interface of silicon films in silicon–silicon films than that corresponding to silicon– diamond films. The predictions of phonon temperature distribution in the silicon film agree well with its counterpart reported in the previous studies. & 2011 Elsevier B.V. All rights reserved.
Keywords: Phonon transport Silicon Diamond EPRT
1. Introduction Thin solid films of dielectric materials are used in electronic components. The heat transfer from the electronic device is important, since it has an effect on the performance of the device. Moreover, heat conduction enhancement is possible from the barrier layers between the electrically conducting regions. In this case, to transfer heat from the electronic device, high thermal conductivity dielectric films such as diamond film is desirable. Heat transfer in dielectric films is governed by the phonon transport through atomic or crystal vibrations. The phonon transport can be analyzed through radiative transfer, which is formulated based on the Boltzmann Transport Theory. The resulting equation of phonon radiative transfer (EPRT) governs the energy transport in the dielectric film. Although the fundamentals of EPRT are established and being well presented in the literature [1–7], the transient solution of energy transport within two adjacent thin dielectric films is interesting in electronic cooling applications. Considerable research studies were carried out to examine the phonon transport in dielectric films. Microscale heat conduction in dielectric films was studied by Majumdar [1]. He showed that hyperbolic heat equation could be derived from the equation of phonon radiative transfer (EPRT) only in the acoustically thick film limit. Joshi and Majumdar [2] introduced transient ballistic and diffusive phonon heat transport in thin films. They indicated that the solution of the EPRT for diamond thin film not only produced
n
Corresponding author. Tel.: þ96638604481; fax: þ 96638602949. E-mail address: [email protected] (B.S. Yilbas).
0921-4526/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.physb.2011.03.028
wall temperature jumps under ballistic transport, but showed markedly different transient response from that of the Fourier law as well as the hyperbolic heat equation even for predominantly diffusive transport. Chen [3] investigated nonequilibrium heat conduction in the vicinity of nanoparticles. He indicated that the solutions approached the prediction of the Fourier law when the particle radius became much larger than the heat-carrier mean free path of the host medium. However, the heat transfer rate from the particle became significantly less and the particle temperature rise was much larger than the prediction of the Fourier conduction theory when particle radius was smaller than the mean free path. Zeng and Chen [4] studied phonon heat conduction in thin films including the thermal boundary resistance and the internal heat generation. They showed that when the free surface behaved as a black phonon emitter, the thermal boundary resistance (TBR) for thin diamond film with internal heat generation was the same as that of without the internal heat generation. However, when the free surface was adiabatic and the film thickness increased, the TBR increased and approached the value of the corresponding black surface. Prasher and Phelan [5] examined a scattering-mediated acoustic mismatch model for the prediction of thermal boundary resistance. They introduced a scattering-mediated acoustic mismatch model exploiting the analogy between phonon and radiative transport by developing a damped wave model to describe the phonon transport. A new numerical scheme was introduced by Pilon and Katika [6] for simulating the multi-dimensional transient and the steady-state microscale energy transport. They used a fixed computational grid and indicated that the modified method of characteristics was unconditionally stable, accurate and compatible with other numerical schemes. The electron–phonon coupling in
S. Bin Mansoor, B.S. Yilbas / Physica B 406 (2011) 2186–2195
Nomenclature
k
x t
v C T I
thermal conductivity (W m 1 K 1) mean free path phonon speed volumetric specific heat capacity (J m 3 K 1) temperature (K) phonon intensity
L
y
m
distance along the film thickness time (s) azimuthal angle cosine of the azimuthal angle
thermal conductance of metal–non-metal interfaces was investigated by Majumdar and Reddy [7]. They indicated that the predictions of the total thermal conductance including electron–phonon and phonon–phonon interfacial conductances showed reasonable agreement with the experimental data at the interface of TiN/MgO. The energy transport and nanostructuring of dielectrics by femtosecond laser pulse trains were investigated by Jiang and Tsai [8]. They indicated that repeatable structures at the desired small nanoscales could be achieved in dielectrics by using the femtosecond pulse trains. Heat conduction in microstructured materials was studied by Miyazaki et al. [9]. They showed that the effective temperature distributions were very different from those of the conventional heat diffusion equation because of the ballistic phonon transport at microscale. The ballistic phonon transport in dielectric membranes was studied by Kuhn and Maasilta [10]. They indicated that for a fixed membrane wall temperature, the radiated power first decreased with decreasing membrane thickness. However, further decrease in membrane thickness resulted in increasing radiated power. Submicron heat transfer model in silicon accounting for phonon dispersion and polarization was investigated by Narumanchi et al. [11]. They incorporated transverse acoustic and longitudinal acoustic phonons as well as optical phonons in the analysis. In the present study, transient phonon energy transport between two adjacent dielectric films, namely silicon–silicon and silicon– diamond, is carried out using the equation of phonon radiative transfer (EPRT). The thermal boundary resistance at the interface of the two films is incorporated in the analysis. The transient variation of temperature in each film is computed. To validate the code developed, temperature variation in a single silicon film is computed after incorporating the initial and boundary conditions reported in Refs. [1,6].
2. Mathematical analysis for equation of phonon radiative transport In thin dielectric films, the film thickness is comparable to the mean-free-path of the phonons and it is necessary to solve the Boltzmann Transport Equation (BTE) to determine the transport properties. The BTE is derived from the Liouville Equation, which governs the time evolution of the N particle distribution function of a collection of particles in 6N dimensional phase space. The Liouville equation can be modified to yield an equation for f (1)(x,y,z,px,py,pzt) (one particle distribution function) [12]: @f p _ rp f ¼ Collision Integral þ :rr f þ p: @t m
ð1Þ
The collision integral is a complicated expression accounting for the affects of collisions between the particles. It may be replaced by a simpler expression by assuming the relaxation time approximation [12]. Therefore, Eq. (1) may then be written as: @f p ðf fo Þ _ rp f ¼ þ :rr f þ p: @t m t
ð2Þ
In the above expression fo is the equilibrium distribution function. For a phonon ‘gas’ fo is the Bose–Einstein distribution function.
2187
A thin dielectric film has a thickness, which is very small as compared to its other two dimensions. In consideration, we may use the one-dimensional BTE to describe the dynamics of the ‘‘dilute phonon gas’’. The independent variable is the one-dimensional distance along the thickness of the thin film. Eq. (2) simplifies to: @f @f @f @f @f ðf fo Þ þvx þ p_ x þ p_ y þ p_ z ¼ @t @x @px @py @pz t
ð3Þ
where f is the single-particle distribution function of the phonons in the phase space. It should be noted that even though f is a function of a single distance variable, it is a function of all the three momentum variables, since the phonons travel in all directions. The dimensions of the distribution function remain same m 3(kg m/s) 3 and a unit length is assumed for the other two directions in the film. The atomic/ molecular lattice vibrates at different frequencies. Although these frequencies are quantized (only certain frequencies are allowed), the difference between two allowable frequencies is negligible. Therefore, the frequency spectrum can be assumed to be continuous. All the phonons travel with a constant speed v, which is the speed of sound in the dielectric; however, their momentums are not the same and are related to the frequency of lattice vibrations by the deBroglie formula [12]. For a single phonon, this can be written as: p¼
h
l
¼
hn ¼ _k v
ð4Þ
k is the wave number. In three dimensions the three components of the momentum are related to the wave vector k. These are given as: px ¼ _kx ,
py ¼ _ky ,
pz ¼ _kz ,
or
p ¼ _k
ð5Þ
^ n, ^ is the unit vector in the direction of propagawhere k ¼ ð2p=lÞn; tion of the wave. The momentum of a phonon could be different, but its speed is the same as those of the other phonons. The momentum and the wave vector may be used interchangeably. It is assumed that the phonons are not under the influence of any attractive/repulsive force and therefore during their motion, i.e.: p_ x ¼ p_ y ¼ p_ z ¼ 0. Therefore, Eq. (3) yields: @f @f ðf fo Þ þvx ¼ @t @x t
ð6Þ
where f (x,px,py,pz,t). We may also write equivalently, f (x,kx,ky,kz,t). The phonons are pictured as moving between the two surfaces of the thin film, across the thickness, which is exactly the same as the photons when they are transferring thermal energy between two surfaces across a vacuum or across a participating medium (a gas which transmits, absorbs and scatters the thermal radiation). Therefore, it is suggested to derive an Equation of Radiative Transfer for phonons just as there exists one for the photons, which is called the equation of phonon radiative transfer (EPRT). Moreover, the spectral intensity is defined as the amount of energy transferred per unit time per unit perpendicular area per unit solid angle per unit wavelength interval, which is Il ¼Il(x,y,z,y,f,t) with unit Wm 2 O 1 m 1. Similarly, the spectral intensity for the phonons can also be defined. Note that the dependence on (y,f) is actually dependence on the solid angle, which in turn implies dependence on the direction. Similarly, dependence on l implies a dependence on the wave number or the frequency, which in turn implies dependence on the momentum and
2188
S. Bin Mansoor, B.S. Yilbas / Physica B 406 (2011) 2186–2195
Substituting Eq. (10) into the right-hand-side Eq. (16) yields: Z oD o Z oD Z @qx Io do 1 ¼2 d o Io dm ð17Þ @x vt vt 1 0 0
the energy of the phonons. The phonon spectral intensity can be related to the phonon distribution function [13], but instead of using the wave vector components, one can use the frequency and the solid angle to specify momentum (energy) and direction, respectively. These two are completely equivalent representations. It is important to note that when the distribution function f is a function of the momentums and its dimension is m 3(kg m/s) 3. However, when it is treated as a function of the frequency and solid angle, then, its dimensions become m 3(rad/s) 1 sr 1. For one-dimensional heat transfer situation, the heat flux in the x-direction can be written as [14]: ZZZ qx ðx,tÞ ¼ vx f _oDðoÞsin f do dy df ð7Þ
Under a phonon radiative equilibrium, @qx/q¼0 and the above equation simplifies to, ! Z 1 Z 1 Z oD do o o 2Io Io dm ¼ 0 or Io ¼ 12 Io dm ð18Þ vt 0 1 1
where D(o) is the density of states of phonons with frequency between o and o þdo, which is:
Note that in general, t ¼ t(o,T). Eq. (19) is the EPRT for the transient heating situation.
DðoÞ ¼
1 o2 V 2p2 v3
ð8Þ
The phonon spectral intensity is defined as [1]: Io ðy, f,x, o,tÞ ¼ vðy, fÞf ðy, f,x, o,tÞ_oDðoÞ so that the heat flux in the x-direction is: ZZZ mIo ðy, f,x, o,tÞsin f dy df do qx ðx,tÞ ¼
ð9Þ
ð10Þ
where, m ¼ nx/n and n is the velocity magnitude. Combining Eqs. (6),(9) and (10), and multiplying Eq. (6) with v_oDðoÞ, it yields:
Substituting Eq. (18) into Eq. (12) yields, R 1 1 Io dmIo 1 @Io @Io þm ¼ 2 1 v @t @x vt
ð19Þ
2.2. Initial and boundary conditions Initially, it is considered that the two films are at the constant ambient temperature of To ¼300 K. This means that the intensity in each film is constant everywhere and corresponds to the temperature To. The intensity may be calculated from the following formula: Ii ðxi , mi ,0Þ ¼
Ci vi To 4p
ð20Þ
where, Ci is the volumetric specific heat capacity of the ith layer,
ni is the phonon speed in the ith layer.
o @I @I ðIo Io Þ þvx ¼ @t @x t
ð11Þ
o ¼ Io ðo,TÞ. T¼T(x) is the temperature. where Io ðy, f,x, o,tÞ and Io Dividing the above equation by the velocity magnitude v of the phonons yields: o 1 @Io @Io ðIo Io Þ þm ¼ v @t vt @x
ð12Þ
Eq. (12) is known as the equation of phonon radiative transfer. In the case of thin films, the equation of phonon radiative transfer in each film is given as: R1 Ii dmi Ii 1 @Ii @I þ mi i ¼ 1 ð13Þ vi @t @xi Li where, Ii(Ii(xi,mi,t)) Ii is the phonon intensity in the ith layer, xi is the coordinate variable along the film thickness in the ith layer, mi is the cosine of the azimuthal angle in the ith layer, t is the time variable, and Li(Li ¼ nt) is the mean free path of the phonons in the ith layer. The heat flux in the x-direction may also be written as: ZZZ mIo ðx, o, O,tÞdo d2 O ð14Þ qx ðx,tÞ ¼
2.2.1. Left boundary condition The left boundary condition corresponds to the front surface of film 1 as depicted in Fig. 1. At all times t40 it is considered that the front surface of the film 1 is maintained at a constant temperature of Tleft ¼301 K. This gives the following boundary condition: I1 ð0, m1 ,tÞ ¼
C1 v1 Tleft 4p
ð21Þ
2.2.2. Right boundary condition The right boundary condition corresponds to the back surface of film 2 as depicted in Fig. 1. At all times t40 it is considered that the back surface of the layer 2 is maintained at a constant temperature of Tback ¼300 K. This gives the following boundary condition: I2 ðL2 , m2 ,tÞ ¼
C2 v2 Tback 4p
I1+
ð22Þ
I1-
I2+
I2-
2.1. Steady-state heat transfer For steady-state heat transfer, we set @Io/@t ¼0 in Eq. (12) and integrating the resulting equation first with respect to m and then with respect to o yields: # # Z oD "Z 1 Z oD "Z 1 o @I ðIo Io Þ m o dm do ¼ dm do ð15Þ vt @x 0 0 1 1
First Film
Second Film
or @ @x
"Z 0
oD
Z
1
!
#
mIo dm do ¼ 1
Z oD Z 0
1 o
Io dm do 1 vt
Z Z oD do 1 Io dm vt 1 0 ð16Þ
Tfront
Interface
Tback
Fig. 1. Schematic view of two thin films and direction of phonon radiation in the films.
S. Bin Mansoor, B.S. Yilbas / Physica B 406 (2011) 2186–2195
2.2.3. Interface condition The interface condition corresponds to the interface between the first and second films as depicted in Fig. 1. The intensity at the interface for the 1st layer is due to the phonons of the 1st film that are reflected back plus the phonons of the 2nd film that are transmitted. ( ) Z 1 Z 0 I1 ðL1 , m1 ,tÞ ¼ 2 R12 I1þ ðL1 , m1 ,tÞm1 dm1 T21 I2 ð0, m2 ,tÞm2 dm2 0
1
ð23Þ The intensity at the interface for the 2nd film is due to the phonons of the 2nd film that are reflected back plus the phonons of the 1st film that are transmitted. ( ) Z 1 Z 0 þ þ I2 ð0, m2 ,tÞ ¼ 2 T12 I1 ðL1 , m1 ,tÞm1 dm1 R21 I2 ð0, m2 ,tÞm2 dm2 0
1
ð24Þ In the above formula T12 is the transmittance for the phonons of the 1st film traveling towards the 2nd film and R12 is the reflectance for these phonons. Similarly, T21 is the transmittance for the phonons of the 2nd film traveling towards the 1st film and R21 is the reflectance for these phonons. These are given below: T12 ¼
C2 v2 C1 v1 : T21 ¼ : R12 ¼ 1T12 C1 v1 þ C2 v2 C1 v1 þ C2 v2 : R21 ¼ 1T21
the films. This resistance can be measured experimentally and its magnitude may be found in the literature for different interfaces. Theoretically, this resistance may be explained on the basis of the assumption that phonons are the dominant thermal energy carriers in a dielectric. Those phonons in the 1st film that are traveling towards the 2nd film get intercepted at the interface. A fraction of these phonons are transmitted in to the 2nd film and the remaining fraction is reflected back in to the 1st film. The same scenario applies to those phonons in the 2nd film that are traveling towards the 1st film. This is the main reason of the thermal boundary resistance. Here we are assuming that the two surfaces are perfectly planar and that there is a perfect contact between the two surfaces. If this is not the case, then, additional thermal resistance due to the roughness of the surfaces and imperfect contact between them may be present. There are two main theoretical models that describe the thermal boundary resistance. These are the Acoustic Mismatch Model and the Diffuse Mismatch Model [16]. In the present study, we consider the diffuse mismatch model. The transmittance and reflectance are calculated based on the (i) energy conservation at the interface and (ii) the diffuse mismatch model. Energy conservation gives: T12 C1 v1 ¼ T21 C2 v2
T12 ¼ 1T21
2.3.2. Equivalent equilibrium phonon temperature (defined anywhere in the films) It represents the average energy of all phonons around a local point and it is equivalent to the equilibrium temperature of phonons when they redistribute adiabatically to an equilibrium state [15]. Z 1 Z 1 2p 2p T1 ¼ I1 dm1 : T2 ¼ I2 dm2 C1 v1 1 C2 v2 1
T12 ¼
C2 v2 , C1 v1 þC2 v2
2.3.4. Heat flux (defined anywhere in the films) 1
1
00
I1 ðx1 , m1 ,tÞm1 dm1 : q2 ¼ 2p
Z
1
I2 ðx2 , m2 ,tÞm2 dm2 1
T21 ¼
C1 v1 C1 v1 þ C2 v2
ð28Þ
The reflectances are given as: R12 ¼ 1T12 ,
R21 ¼ 1T21
ð29Þ
The thermal boundary resistance is defined as: TBR ¼
DT12
ð30Þ
00
q12
In Eq. (30), DT12 ¼T1–T2 is the equivalent equilibrium phonon 00 temperature at the interface for layers 1 and 2, respectively. q12 is the heat flux at the interface. Using the definition of the equivalent equilibrium phonon temperature at the interface and using the expression for the heat flux, the TBR may be expressed as: TBR ¼
2.3.3. Gray-body temperature (defined only at a diffusely scattering interface of the films) It is assumed that the radiative equilibrium over all frequencies takes place at the interface of the films [1]. Z 0 Z 1 4p 4p 4p 4p þ Tgb1 ¼ I1 dm1 ¼ I1 : Tgb2 ¼ I þ dm2 ¼ I C1 v1 1 C1 v1 C2 v2 0 2 C2 v2 2
Z
ð27Þ
Simultaneous solution of the above two equations yields the transmissions:
2.3.1. Emitted phonon temperature (defined at the interface of the films) Emitted phonon temperature is defined due to the emitted phonon radiation at the interface of the films. Z 1 Z 1 4p 4p Te1 ¼ I1þ dm1 : Te2 ¼ I dm2 C1 v1 0 C2 v2 0 2
00
ð26Þ
The diffuse mismatch model gives: ð25Þ
2.3. Definition of temperatures and heat flux
q1 ¼ 2p
2189
DT12 q00 12
¼
2 2 ¼ T12 C1 v1 T21 C2 v2
ð31Þ
3. Numerical method The Equation of Phonon Radiative Transfer in each film with the associated initial and boundary conditions are solved numerically using the finite-difference method. These are integro-differential equations and may be solved using an iterative technique. Since the equations are time-dependent, therefore, the implicit scheme is used to solve them. Consider the EPRT in the first film: R 1 1 2 1 I1 d 1 I1
m
2.4. Thermal boundary resistance
1 @I1 @I1 þ m1 ¼ v1 @t @x1
At the interface of two dielectric films or a metal and a dielectric film there exists a resistance to heat flow, which is called the thermal boundary resistance. Due to this resistance, there is a small but finite temperature drop across the interface of
The variable m1 varies between 1 r m1 r1. To facilitate the solution of the EPRT it is customary to break the intensity into parts. When 1 r m1 r 0, I1 ¼ I1 ðx1 , m1 ,tÞ and when 0 r m1 r1, I1 ¼ I1þ ðx1 , m1 ,tÞ. Therefore, we may write separate equations for
L1
2190
S. Bin Mansoor, B.S. Yilbas / Physica B 406 (2011) 2186–2195
the two intensities. These are given below: R R1 þ 0 þ 1 @I1þ 1 I1 dm1 þ 0 I1 dm1 I1 2 1 @I1þ þ m1 ¼ , v1 @t @x1 L1
0 r m1 r 1
and @I 1 @I1 þ m1 1 ¼ v1 @t @x1
1 2
R 0
m
1 I1 d 1 þ
R1 0
I1þ dm1 I1
L1
Layer 1 Layer 2
,
and 1 2
R 0
m
1 I2 d 2 þ
R1 0
I2þ dm2 I2
L2
,
1r m2 r 0
For the intensity I2þ , the initial condition is I2þ ðx2 , m2 ,0Þ ¼ C2 v2 To = 4p and the boundaryn condition comes from the interface condiR1 R0 tion I2þ ð0, m2 ,tÞ ¼ 2 T12 0 I1þ ðL1 , m1 ,tÞm1 dm1 R21 1 I2 ð0, m2 ,tÞm2 dm2 g. For the intensity I2 , the initial condition is I2 ðx2 , m2 ,0Þ ¼ C2 v2 To =4p and the boundary condition is I2 ðL2 , m2 ,tÞ ¼ C2 v2 Tright =4p. Starting from the initial intensity distributions in each film, at each time-step two equations are solved in each film. This solution gives the approximate intensity distribution in each film. From this solution the interface temperature is calculated and the next iteration is performed. After convergence the computations proceed to the next time-step. Once the simulation is complete, the equivalent equilibrium phonon temperature may be evaluated at each time-step and at each location inside the two films. Table 1 gives the properties of the silicon and the diamond films used in the simulations. Table 2 gives the geometric configuration of the films and initial and boundary temperatures.
Table 1 Properties of silicon and diamond films used in the simulations.
Silicon Diamond
Length (mm)
To (K)
Tleft (K)
Tright (K)
0.1 0.1
300 300
301 –
– 300
1r m1 r 0
To solve the above equations an implicit method is used for the time discretisation and discrete ordinates method is used for the discretisation of the angular term. The time-derivative in each equation is approximated by means of the first-order backwarddifference. The partial derivative with respect to the coordinate variable x1 is descritised using a backward-difference for the first equation in which 0r m1 r1 and is descritised using a forwarddifference for the second equation in which 1r m1 r0. This procedure is important to obtain stable discretisations. The integrals are calculated by first discretising the cosine of the azimuthal angle m1 into 40 divisions and then using Simpson’s rule to evaluate the integrals numerically for each time-step and each coordinate location. It is important to keep in mind that the intensities I1þ and I1 are two separate variables and each needs an initial condition and a single boundary condition. For the intensity I1þ , the initial condition is I1þ ðx1 , m1 ,0Þ ¼ C1 v1 To =4p and the boundary condition is I1þ ð0, m1 ,tÞ ¼ C1 v1 Tleft =4p. For the intensity I1 , the initial condition is I1 ðx1 , m1 ,0Þ ¼ C1 v1 To =4p and the boundary ncondition R1 comes from the interface condition I1 ðL1 , m1 ,tÞ ¼ 2 R12 0 I1þ R0 ðL1 , m1 ,tÞm1 dm1 T21 1 I2 ð0, m2 ,tÞm2 dm2 g. Exactly similar considerations apply to the intensity calculations in the 2nd film. The EPRT is first written in terms of the two intensities as given below: R R1 þ 0 þ 1 @I2þ 1 I2 dm2 þ 0 I2 dm2 I2 2 1 @I2þ þ m2 ¼ , 0 r m2 r1 v2 @t @x2 L2
@I 1 @I2 þ m2 2 ¼ v2 @t @x2
Table 2 Geometric configurations of the film and initial and boundary temperatures of the films.
Thermal conductivity (W m 1 K 1)
Volumetric specific heat (J m 3 K 1)
Phonon speed (m s 1)
Mean free path (nm)
157 3320
1.653 106 1.815 106
8430 12,288
33.8 447
Table 3 Transmittance and reflectance phonon radiation determined from diffusive mismatch model. Transmittance
Si/Si Si/Di
1-2 0.5 0.6154
Reflectance 2-1 0.5 0.3846
1-2 0.5 0.3846
2-1 0.5 0.6154
Table 3 gives the transmittance and reflectance of phonon radiation at the film surfaces.
4. Results and discussion The transient analysis of phonon transport in dielectric thin films is carried out and temperature variation along the films is predicted. Silicon–silicon and silicon–diamond thin films are incorporated in the simulations. The thermal boundary resistance at the interface of the films is accounted in the analysis. To validate the code developed for EPRT, the steady state transport across a single silicon film is simulated in line with the previous study [1]. The simulation results for different film thickness are shown in Fig. 2. The findings are exactly the same as the previous results [1]. Moreover, the validation study is extended to include the transient phonon transport across the silicon film in accordance with the early study [6]. Dimensionless temperature variation along the film thickness is shown in Fig. 3 for different dimensionless time. The findings are the same as those reported in the early study [6], i.e. dimensionless temperature slip at diffused film surfaces and dimensionless temperature distributions for the previously reported results. Fig. 4 shows the dimensionless temperature distribution along the silicon–silicon films for three dimensionless heating durations. It should be noted that the dimensionless temperature is set at 1 at time t*Z0 and radiative phonon transport is initiated from the front surface of the silicon film (x*¼0). In addition, initially both films are considered to be in temperature equilibrium. The emitted, equilibrium and gray body temperatures are shown at the interface. The temperature slip occurs at the interface because of the transmission and the reflectance of the radiative phonon energy at the interface. Due to the transient nature of the phonon transport, temperature decay in the first film of silicon is sharp while it is gradual in the adjacent silicon film. In addition, the temperature slip at the interface increases with progressing time. This is attributed to the phonon scattering in the film and the phonon energy reaching at the interface, which is small during the early periods. As the time progresses, the radiative phonon energy reaching at the interface becomes high, which in turn increases the equilibrium temperature at the interface. Moreover, as the heating period progresses further, the steady state is reached and temperature variation along the film thickness becomes linear. This is true for both films. It should be noted that temperature in the film corresponds to the energy of all phonons at a local point in the film. Furthermore, although the local thermodynamic equilibrium condition is not satisfied in terms of a
S. Bin Mansoor, B.S. Yilbas / Physica B 406 (2011) 2186–2195
2191
Fig. 4. Dimensionless temperature with dimensionless thickness along silicon– silicon. thin films for different heating periods. Dimensionless Temperature is defined as (Tequilibrium–Tback)/(Tfront–Tback).
Fig. 2. Dimensionless temperature with dimensionless film thickness for steady state condition similar to that of Ref. [1]. Dimensionless temperature is defined as (Tequilibrium–Tback)/(Tfront–Tback).
Fig. 3. Dimensionless temperature with dimensionless film thickness for transient condition similar to that of Ref. [6]. Dimensionless temperature is defined as (Tequilibrium–Tback)/(Tfront–Tback).
thermal equilibrium state, the equilibrium temperature is considered for the phonons when they are redistributed adiabatically to an equilibrium state. Consequently, the equivalent equilibrium temperatures on both sides of the interface across the films are not the
same due to the reflection and transmission of phonons at the interface. Fig. 5 shows temporal variation of temperatures at both surfaces of the interface. The emitted phonon temperature represents the phonons that are emitted at each side of the interface at this temperature. Therefore, the emitted phonon temperatures are different than the equilibrium temperature of phonons. It should be noted that the equilibrium temperature is more appropriate for the predictions of thermal boundary resistance across the interface and represents the measurable temperature at both sides of interface [15]. Temperature at the interface of the front film (silicon) increases with progressing time. This is because of the reflected phonons from the interface towards the front film. Since the phonon transport is not diffusive in the film, the difference between the emitted temperature and the gray body temperature becomes large with progressing time. It should be noted that due to non-diffusive phonon transport in the film, temperature slip occurs at the boundary and the Fourier theory of heating fails to govern the phonon transport in the film. The temperature slip is associated with the phonons reaching at the interface, which is small in the early heating period. Moreover, the temperature slip is also noted at the interface of the second film (silicon), provided that the equilibrium temperature and the gray body temperatures are higher than of the emitted phonon temperature at the interface. This is associated with reflected phonons from the outer boundary of the second film towards the interface, since the emitted phonon temperature is related to phonons reaching at the interface when emitted from the free surface of the layers. In addition, the thermal boundary resistance lowers the equilibrium phonon temperature at the interface of the second layer. Since the gray body temperature is determined from the consideration of the diffusive boundary [1], it remains the same at both sides of the interface due to energy balance. Fig. 6 shows temporal variation of heat flux due to phonon transport at the front surface of the film and the back surface of the second film. It should be noted that the front surface of the first layer is considered to be at a higher temperature than the back surface of the second layer prior to heating is initiated (t*o0). Consequently, heat flux at the front surface of the first film is higher in the early heating period and, first, reduces sharply tr100 ps and decays gradually to reach a steady value. The opposite is true for the
2192
S. Bin Mansoor, B.S. Yilbas / Physica B 406 (2011) 2186–2195
Fig. 5. Dimensionless temperatures with time at the silicon–silicon interface. Dimensionless temperature is defined as (Tequilibrium–Tback)/(Tfront–Tback).
second layer. The rapid decay of heat flux at the surface of the first film is because of the high rate of phonon transport from surface region towards the bulk of the film during the early heating period. However, temperature rise in the film with progressing time (Fig. 4) lowers the amount of radiative phonon transport from this region to the bulk. Once, the steady temperature (Fig. 4) is achieved within both films, the radiative phonon transport within the film becomes steady, i.e. the heat flux remains the same. Fig. 7 shows temperature distribution in the silicon–diamond films for three heating periods. The first film is silicon and the second film is diamond. Temperature behavior in the silicon film is similar to that shown in Fig. 4, provided that interface temperature in the silicon film is lower than that corresponding to the silicon–silicon films. This is attributed to (i) phonon
transmittance and reflectance at the interface as given in Table 3 and (ii) the thermal boundary resistance determined from the diffuse mismatch model, which is different than that of the silicon films. Moreover, the large phonon mean free path of diamond lowers the phonon scattering in this film. Consequently, the transmitted phonons reflect from the back surface of the diamond film towards the interface. The attainment of low equilibrium phonon temperature at the silicon film interface is associated with the large phonon transmission from the silicon film to the diamond film. Fig. 8 shows temporal behavior of the emitted phonon temperature, the equilibrium phonon temperature and the gray body temperature at the interface of the silicon and diamond films. Temperature behavior is similar to those shown in Fig. 5.
S. Bin Mansoor, B.S. Yilbas / Physica B 406 (2011) 2186–2195
2193
Fig. 6. Heat flux with time at the front and back surfaces of silicon–silicon films.
However, the maximum temperatures at interface of both sides are lower than that corresponding to the silicon–silicon films. This is attributed to small scattering of phonons in the diamond film due to large phonon mean free path of diamond, and transmittance and reflectance at the interface. Moreover, temporal variation of heat flux is shown in Fig. 9 for the front surface of silicon film and the back surface of the diamond film. The flux at the back surface of diamond film is much higher than that of the front surface of silicon film. This is more pronounced as the heating period progresses. The attainment of high heat flux is associated with the phonon reaching at the back surface of the diamond film, which is large due to the large mean phonon mean free path of diamond. In addition, phonon transmission from the silicon interface contributes to the attainment of the high heat flux at the back surface of diamond.
5. Conclusion The phonon transport across the silicon–silicon films and silicon–diamond films is considered. The steady state analysis is carried out to validate the phonon temperature distribution in the silicon layer. The transient analysis is carried out after assuming the presence of temperature difference of 1 K across the films. The EPRT is solved numerically to obtain phonon transport in the films. The diffuse mismatch model is incorporated to formulate and determine the thermal boundary resistance across the interface of the films. It is found that phonon temperature variation across the silicon film predicted from the present study for steady and transient situations agree well with the previously reported results. In silicon–silicon films, temperature slip occurs at the interface, despite the fact that the phonon
2194
S. Bin Mansoor, B.S. Yilbas / Physica B 406 (2011) 2186–2195
Fig. 7. Dimensionless temperature with dimensionless thickness along silicon– diamond thin films for different heating periods. Dimensionless temperature is defined as (Tequilibrium–Tback)/(Tfront–Tback).
mean free path remains the same in both films. This is associated with the transmittance and reflectance of phonons at the interface of the films. The rise of the equilibrium, emitted and gray body temperatures is higher at the interface of the silicon films during the heating period 36rt r72 ps. This is because of the small mean free path of silicon; in which case, phonons reach the interface at higher rates. Temperature slip is higher at the interface of silicon–diamond films, which is associated with the thermal boundary resistance and the large phonon mean free path of diamond. In this case, the transmittance at interface increases and less scattering in the diamond film allows transmitted phonons to reflect from the back surface of the diamond to interface at a large amount due to small scattering. The heat flux increases at the front surface of silicon film and at the back surface of the diamond film with progressing time. This is more pronounced at the back surface of the diamond film due to the large number of phonon reaching at the back surface of diamond film because of less scattering of transmitted phonons in the diamond film.
Fig. 8. Dimensionless temperatures with time at the silicon–diamond interface. Dimensionless temperature is defined as (Tequilibrium–Tback)/(Tfront–Tback).
S. Bin Mansoor, B.S. Yilbas / Physica B 406 (2011) 2186–2195
2195
Fig. 9. Heat flux with time at the front and back surfaces of silicon–diamond films.
Acknowledgments The authors acknowledge the support of Center of Excellence for Scientific Research Collaboration with MIT and King Fahd University of Petroleum and Minerals, Dhahran, Saudi Arabia for this work. References
[6] [7] [8] [9] [10] [11] [12] [13]
[1] [2] [3] [4] [5]
A. Majumdar, ASME J. Heat Transfer 115 (1993) 7. A.A. Joshi, A. Majumdar, J. Appl. Phys. 74 (1) (1993) 31. G. Chen, ASME J. Heat Transfer 118 (1996) 539. T. Zeng, G. Chen, ASME J. Heat Transfer 123 (2001) 340. R.S. Prasher, P.E. Phelan, ASME J. Heat Transfer 123 (2001) 105.
[14] [15] [16]
L. Pilon, K.M. Katika, ASME J. Heat Transfer 126 (2004) 735. A. Majumdar, P. Reddy, Appl. Phys. Lett. 84 (23) (2004) 4768. L. Jiang, H. Tsai, ASME J. Heat Transfer 128 (2006) 926. K. Miyazaki, T. Arashi, D. Makino, H. Tsukamoto, IEEE Trans. Components Packag. Technol. 29 (2) (2006) 247. T. Kuhn, I.J. Maasilta, Nucl. Instrum. Methods Phys. Res. A 559 (2006) 724. S.V.J. Narumanchi, J.Y. Murthy, C.H. Amon, ASME J. Heat Transfer 126 (6) (2004) 946. C. Kittel, Introduction to Solid State Physics, sixth ed., Wiley, New York, 1986. W.G. Vincenti, C.H. Kruger, Introduction to Physical Gas Dynamics, Robert Krieger Publications, New York, 1977. M.N. Ozisik, Radiative Transfer and Interactions with Conduction and Convection, Wiley, New York, 1973. G. Chen, Phys. Rev. B 57 (23) (1998) 14958. E.T. Swartz, R.O. Pohl, Rev. Mod. Phys. 61 (3) (1989) 605.
Contents lists available at ScienceDirect
Physica B journal homepage: www.elsevier.com/locate/physb
Phonon transport in silicon–silicon and silicon–diamond thin films: Consideration of thermal boundary resistance at interface S. Bin Mansoor, B.S. Yilbas n ME Department, King Fahd University of Petroleum and Minerals, Saudi Arabia
a r t i c l e i n f o
a b s t r a c t
Article history: Received 14 February 2011 Accepted 13 March 2011 Available online 30 March 2011
Phonon transport across the silicon–silicon and silicon–diamond dielectric films is examined. To simulate the phonon transport in the films EPRT is accounted and numerical solutions are obtained with the consideration of 1 K difference across the films prior to the initiation of the phonon transport. The diffuse mismatch model is introduced to account for the thermal boundary resistance at the interface of the films. To validate the code developed, the steady and transient cases for phonon transport in silicon film are carried out. It is found that the equilibrium temperature of phonons attains higher value at the interface of silicon films in silicon–silicon films than that corresponding to silicon– diamond films. The predictions of phonon temperature distribution in the silicon film agree well with its counterpart reported in the previous studies. & 2011 Elsevier B.V. All rights reserved.
Keywords: Phonon transport Silicon Diamond EPRT
1. Introduction Thin solid films of dielectric materials are used in electronic components. The heat transfer from the electronic device is important, since it has an effect on the performance of the device. Moreover, heat conduction enhancement is possible from the barrier layers between the electrically conducting regions. In this case, to transfer heat from the electronic device, high thermal conductivity dielectric films such as diamond film is desirable. Heat transfer in dielectric films is governed by the phonon transport through atomic or crystal vibrations. The phonon transport can be analyzed through radiative transfer, which is formulated based on the Boltzmann Transport Theory. The resulting equation of phonon radiative transfer (EPRT) governs the energy transport in the dielectric film. Although the fundamentals of EPRT are established and being well presented in the literature [1–7], the transient solution of energy transport within two adjacent thin dielectric films is interesting in electronic cooling applications. Considerable research studies were carried out to examine the phonon transport in dielectric films. Microscale heat conduction in dielectric films was studied by Majumdar [1]. He showed that hyperbolic heat equation could be derived from the equation of phonon radiative transfer (EPRT) only in the acoustically thick film limit. Joshi and Majumdar [2] introduced transient ballistic and diffusive phonon heat transport in thin films. They indicated that the solution of the EPRT for diamond thin film not only produced
n
Corresponding author. Tel.: þ96638604481; fax: þ 96638602949. E-mail address: [email protected] (B.S. Yilbas).
0921-4526/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.physb.2011.03.028
wall temperature jumps under ballistic transport, but showed markedly different transient response from that of the Fourier law as well as the hyperbolic heat equation even for predominantly diffusive transport. Chen [3] investigated nonequilibrium heat conduction in the vicinity of nanoparticles. He indicated that the solutions approached the prediction of the Fourier law when the particle radius became much larger than the heat-carrier mean free path of the host medium. However, the heat transfer rate from the particle became significantly less and the particle temperature rise was much larger than the prediction of the Fourier conduction theory when particle radius was smaller than the mean free path. Zeng and Chen [4] studied phonon heat conduction in thin films including the thermal boundary resistance and the internal heat generation. They showed that when the free surface behaved as a black phonon emitter, the thermal boundary resistance (TBR) for thin diamond film with internal heat generation was the same as that of without the internal heat generation. However, when the free surface was adiabatic and the film thickness increased, the TBR increased and approached the value of the corresponding black surface. Prasher and Phelan [5] examined a scattering-mediated acoustic mismatch model for the prediction of thermal boundary resistance. They introduced a scattering-mediated acoustic mismatch model exploiting the analogy between phonon and radiative transport by developing a damped wave model to describe the phonon transport. A new numerical scheme was introduced by Pilon and Katika [6] for simulating the multi-dimensional transient and the steady-state microscale energy transport. They used a fixed computational grid and indicated that the modified method of characteristics was unconditionally stable, accurate and compatible with other numerical schemes. The electron–phonon coupling in
S. Bin Mansoor, B.S. Yilbas / Physica B 406 (2011) 2186–2195
Nomenclature
k
x t
v C T I
thermal conductivity (W m 1 K 1) mean free path phonon speed volumetric specific heat capacity (J m 3 K 1) temperature (K) phonon intensity
L
y
m
distance along the film thickness time (s) azimuthal angle cosine of the azimuthal angle
thermal conductance of metal–non-metal interfaces was investigated by Majumdar and Reddy [7]. They indicated that the predictions of the total thermal conductance including electron–phonon and phonon–phonon interfacial conductances showed reasonable agreement with the experimental data at the interface of TiN/MgO. The energy transport and nanostructuring of dielectrics by femtosecond laser pulse trains were investigated by Jiang and Tsai [8]. They indicated that repeatable structures at the desired small nanoscales could be achieved in dielectrics by using the femtosecond pulse trains. Heat conduction in microstructured materials was studied by Miyazaki et al. [9]. They showed that the effective temperature distributions were very different from those of the conventional heat diffusion equation because of the ballistic phonon transport at microscale. The ballistic phonon transport in dielectric membranes was studied by Kuhn and Maasilta [10]. They indicated that for a fixed membrane wall temperature, the radiated power first decreased with decreasing membrane thickness. However, further decrease in membrane thickness resulted in increasing radiated power. Submicron heat transfer model in silicon accounting for phonon dispersion and polarization was investigated by Narumanchi et al. [11]. They incorporated transverse acoustic and longitudinal acoustic phonons as well as optical phonons in the analysis. In the present study, transient phonon energy transport between two adjacent dielectric films, namely silicon–silicon and silicon– diamond, is carried out using the equation of phonon radiative transfer (EPRT). The thermal boundary resistance at the interface of the two films is incorporated in the analysis. The transient variation of temperature in each film is computed. To validate the code developed, temperature variation in a single silicon film is computed after incorporating the initial and boundary conditions reported in Refs. [1,6].
2. Mathematical analysis for equation of phonon radiative transport In thin dielectric films, the film thickness is comparable to the mean-free-path of the phonons and it is necessary to solve the Boltzmann Transport Equation (BTE) to determine the transport properties. The BTE is derived from the Liouville Equation, which governs the time evolution of the N particle distribution function of a collection of particles in 6N dimensional phase space. The Liouville equation can be modified to yield an equation for f (1)(x,y,z,px,py,pzt) (one particle distribution function) [12]: @f p _ rp f ¼ Collision Integral þ :rr f þ p: @t m
ð1Þ
The collision integral is a complicated expression accounting for the affects of collisions between the particles. It may be replaced by a simpler expression by assuming the relaxation time approximation [12]. Therefore, Eq. (1) may then be written as: @f p ðf fo Þ _ rp f ¼ þ :rr f þ p: @t m t
ð2Þ
In the above expression fo is the equilibrium distribution function. For a phonon ‘gas’ fo is the Bose–Einstein distribution function.
2187
A thin dielectric film has a thickness, which is very small as compared to its other two dimensions. In consideration, we may use the one-dimensional BTE to describe the dynamics of the ‘‘dilute phonon gas’’. The independent variable is the one-dimensional distance along the thickness of the thin film. Eq. (2) simplifies to: @f @f @f @f @f ðf fo Þ þvx þ p_ x þ p_ y þ p_ z ¼ @t @x @px @py @pz t
ð3Þ
where f is the single-particle distribution function of the phonons in the phase space. It should be noted that even though f is a function of a single distance variable, it is a function of all the three momentum variables, since the phonons travel in all directions. The dimensions of the distribution function remain same m 3(kg m/s) 3 and a unit length is assumed for the other two directions in the film. The atomic/ molecular lattice vibrates at different frequencies. Although these frequencies are quantized (only certain frequencies are allowed), the difference between two allowable frequencies is negligible. Therefore, the frequency spectrum can be assumed to be continuous. All the phonons travel with a constant speed v, which is the speed of sound in the dielectric; however, their momentums are not the same and are related to the frequency of lattice vibrations by the deBroglie formula [12]. For a single phonon, this can be written as: p¼
h
l
¼
hn ¼ _k v
ð4Þ
k is the wave number. In three dimensions the three components of the momentum are related to the wave vector k. These are given as: px ¼ _kx ,
py ¼ _ky ,
pz ¼ _kz ,
or
p ¼ _k
ð5Þ
^ n, ^ is the unit vector in the direction of propagawhere k ¼ ð2p=lÞn; tion of the wave. The momentum of a phonon could be different, but its speed is the same as those of the other phonons. The momentum and the wave vector may be used interchangeably. It is assumed that the phonons are not under the influence of any attractive/repulsive force and therefore during their motion, i.e.: p_ x ¼ p_ y ¼ p_ z ¼ 0. Therefore, Eq. (3) yields: @f @f ðf fo Þ þvx ¼ @t @x t
ð6Þ
where f (x,px,py,pz,t). We may also write equivalently, f (x,kx,ky,kz,t). The phonons are pictured as moving between the two surfaces of the thin film, across the thickness, which is exactly the same as the photons when they are transferring thermal energy between two surfaces across a vacuum or across a participating medium (a gas which transmits, absorbs and scatters the thermal radiation). Therefore, it is suggested to derive an Equation of Radiative Transfer for phonons just as there exists one for the photons, which is called the equation of phonon radiative transfer (EPRT). Moreover, the spectral intensity is defined as the amount of energy transferred per unit time per unit perpendicular area per unit solid angle per unit wavelength interval, which is Il ¼Il(x,y,z,y,f,t) with unit Wm 2 O 1 m 1. Similarly, the spectral intensity for the phonons can also be defined. Note that the dependence on (y,f) is actually dependence on the solid angle, which in turn implies dependence on the direction. Similarly, dependence on l implies a dependence on the wave number or the frequency, which in turn implies dependence on the momentum and
2188
S. Bin Mansoor, B.S. Yilbas / Physica B 406 (2011) 2186–2195
Substituting Eq. (10) into the right-hand-side Eq. (16) yields: Z oD o Z oD Z @qx Io do 1 ¼2 d o Io dm ð17Þ @x vt vt 1 0 0
the energy of the phonons. The phonon spectral intensity can be related to the phonon distribution function [13], but instead of using the wave vector components, one can use the frequency and the solid angle to specify momentum (energy) and direction, respectively. These two are completely equivalent representations. It is important to note that when the distribution function f is a function of the momentums and its dimension is m 3(kg m/s) 3. However, when it is treated as a function of the frequency and solid angle, then, its dimensions become m 3(rad/s) 1 sr 1. For one-dimensional heat transfer situation, the heat flux in the x-direction can be written as [14]: ZZZ qx ðx,tÞ ¼ vx f _oDðoÞsin f do dy df ð7Þ
Under a phonon radiative equilibrium, @qx/q¼0 and the above equation simplifies to, ! Z 1 Z 1 Z oD do o o 2Io Io dm ¼ 0 or Io ¼ 12 Io dm ð18Þ vt 0 1 1
where D(o) is the density of states of phonons with frequency between o and o þdo, which is:
Note that in general, t ¼ t(o,T). Eq. (19) is the EPRT for the transient heating situation.
DðoÞ ¼
1 o2 V 2p2 v3
ð8Þ
The phonon spectral intensity is defined as [1]: Io ðy, f,x, o,tÞ ¼ vðy, fÞf ðy, f,x, o,tÞ_oDðoÞ so that the heat flux in the x-direction is: ZZZ mIo ðy, f,x, o,tÞsin f dy df do qx ðx,tÞ ¼
ð9Þ
ð10Þ
where, m ¼ nx/n and n is the velocity magnitude. Combining Eqs. (6),(9) and (10), and multiplying Eq. (6) with v_oDðoÞ, it yields:
Substituting Eq. (18) into Eq. (12) yields, R 1 1 Io dmIo 1 @Io @Io þm ¼ 2 1 v @t @x vt
ð19Þ
2.2. Initial and boundary conditions Initially, it is considered that the two films are at the constant ambient temperature of To ¼300 K. This means that the intensity in each film is constant everywhere and corresponds to the temperature To. The intensity may be calculated from the following formula: Ii ðxi , mi ,0Þ ¼
Ci vi To 4p
ð20Þ
where, Ci is the volumetric specific heat capacity of the ith layer,
ni is the phonon speed in the ith layer.
o @I @I ðIo Io Þ þvx ¼ @t @x t
ð11Þ
o ¼ Io ðo,TÞ. T¼T(x) is the temperature. where Io ðy, f,x, o,tÞ and Io Dividing the above equation by the velocity magnitude v of the phonons yields: o 1 @Io @Io ðIo Io Þ þm ¼ v @t vt @x
ð12Þ
Eq. (12) is known as the equation of phonon radiative transfer. In the case of thin films, the equation of phonon radiative transfer in each film is given as: R1 Ii dmi Ii 1 @Ii @I þ mi i ¼ 1 ð13Þ vi @t @xi Li where, Ii(Ii(xi,mi,t)) Ii is the phonon intensity in the ith layer, xi is the coordinate variable along the film thickness in the ith layer, mi is the cosine of the azimuthal angle in the ith layer, t is the time variable, and Li(Li ¼ nt) is the mean free path of the phonons in the ith layer. The heat flux in the x-direction may also be written as: ZZZ mIo ðx, o, O,tÞdo d2 O ð14Þ qx ðx,tÞ ¼
2.2.1. Left boundary condition The left boundary condition corresponds to the front surface of film 1 as depicted in Fig. 1. At all times t40 it is considered that the front surface of the film 1 is maintained at a constant temperature of Tleft ¼301 K. This gives the following boundary condition: I1 ð0, m1 ,tÞ ¼
C1 v1 Tleft 4p
ð21Þ
2.2.2. Right boundary condition The right boundary condition corresponds to the back surface of film 2 as depicted in Fig. 1. At all times t40 it is considered that the back surface of the layer 2 is maintained at a constant temperature of Tback ¼300 K. This gives the following boundary condition: I2 ðL2 , m2 ,tÞ ¼
C2 v2 Tback 4p
I1+
ð22Þ
I1-
I2+
I2-
2.1. Steady-state heat transfer For steady-state heat transfer, we set @Io/@t ¼0 in Eq. (12) and integrating the resulting equation first with respect to m and then with respect to o yields: # # Z oD "Z 1 Z oD "Z 1 o @I ðIo Io Þ m o dm do ¼ dm do ð15Þ vt @x 0 0 1 1
First Film
Second Film
or @ @x
"Z 0
oD
Z
1
!
#
mIo dm do ¼ 1
Z oD Z 0
1 o
Io dm do 1 vt
Z Z oD do 1 Io dm vt 1 0 ð16Þ
Tfront
Interface
Tback
Fig. 1. Schematic view of two thin films and direction of phonon radiation in the films.
S. Bin Mansoor, B.S. Yilbas / Physica B 406 (2011) 2186–2195
2.2.3. Interface condition The interface condition corresponds to the interface between the first and second films as depicted in Fig. 1. The intensity at the interface for the 1st layer is due to the phonons of the 1st film that are reflected back plus the phonons of the 2nd film that are transmitted. ( ) Z 1 Z 0 I1 ðL1 , m1 ,tÞ ¼ 2 R12 I1þ ðL1 , m1 ,tÞm1 dm1 T21 I2 ð0, m2 ,tÞm2 dm2 0
1
ð23Þ The intensity at the interface for the 2nd film is due to the phonons of the 2nd film that are reflected back plus the phonons of the 1st film that are transmitted. ( ) Z 1 Z 0 þ þ I2 ð0, m2 ,tÞ ¼ 2 T12 I1 ðL1 , m1 ,tÞm1 dm1 R21 I2 ð0, m2 ,tÞm2 dm2 0
1
ð24Þ In the above formula T12 is the transmittance for the phonons of the 1st film traveling towards the 2nd film and R12 is the reflectance for these phonons. Similarly, T21 is the transmittance for the phonons of the 2nd film traveling towards the 1st film and R21 is the reflectance for these phonons. These are given below: T12 ¼
C2 v2 C1 v1 : T21 ¼ : R12 ¼ 1T12 C1 v1 þ C2 v2 C1 v1 þ C2 v2 : R21 ¼ 1T21
the films. This resistance can be measured experimentally and its magnitude may be found in the literature for different interfaces. Theoretically, this resistance may be explained on the basis of the assumption that phonons are the dominant thermal energy carriers in a dielectric. Those phonons in the 1st film that are traveling towards the 2nd film get intercepted at the interface. A fraction of these phonons are transmitted in to the 2nd film and the remaining fraction is reflected back in to the 1st film. The same scenario applies to those phonons in the 2nd film that are traveling towards the 1st film. This is the main reason of the thermal boundary resistance. Here we are assuming that the two surfaces are perfectly planar and that there is a perfect contact between the two surfaces. If this is not the case, then, additional thermal resistance due to the roughness of the surfaces and imperfect contact between them may be present. There are two main theoretical models that describe the thermal boundary resistance. These are the Acoustic Mismatch Model and the Diffuse Mismatch Model [16]. In the present study, we consider the diffuse mismatch model. The transmittance and reflectance are calculated based on the (i) energy conservation at the interface and (ii) the diffuse mismatch model. Energy conservation gives: T12 C1 v1 ¼ T21 C2 v2
T12 ¼ 1T21
2.3.2. Equivalent equilibrium phonon temperature (defined anywhere in the films) It represents the average energy of all phonons around a local point and it is equivalent to the equilibrium temperature of phonons when they redistribute adiabatically to an equilibrium state [15]. Z 1 Z 1 2p 2p T1 ¼ I1 dm1 : T2 ¼ I2 dm2 C1 v1 1 C2 v2 1
T12 ¼
C2 v2 , C1 v1 þC2 v2
2.3.4. Heat flux (defined anywhere in the films) 1
1
00
I1 ðx1 , m1 ,tÞm1 dm1 : q2 ¼ 2p
Z
1
I2 ðx2 , m2 ,tÞm2 dm2 1
T21 ¼
C1 v1 C1 v1 þ C2 v2
ð28Þ
The reflectances are given as: R12 ¼ 1T12 ,
R21 ¼ 1T21
ð29Þ
The thermal boundary resistance is defined as: TBR ¼
DT12
ð30Þ
00
q12
In Eq. (30), DT12 ¼T1–T2 is the equivalent equilibrium phonon 00 temperature at the interface for layers 1 and 2, respectively. q12 is the heat flux at the interface. Using the definition of the equivalent equilibrium phonon temperature at the interface and using the expression for the heat flux, the TBR may be expressed as: TBR ¼
2.3.3. Gray-body temperature (defined only at a diffusely scattering interface of the films) It is assumed that the radiative equilibrium over all frequencies takes place at the interface of the films [1]. Z 0 Z 1 4p 4p 4p 4p þ Tgb1 ¼ I1 dm1 ¼ I1 : Tgb2 ¼ I þ dm2 ¼ I C1 v1 1 C1 v1 C2 v2 0 2 C2 v2 2
Z
ð27Þ
Simultaneous solution of the above two equations yields the transmissions:
2.3.1. Emitted phonon temperature (defined at the interface of the films) Emitted phonon temperature is defined due to the emitted phonon radiation at the interface of the films. Z 1 Z 1 4p 4p Te1 ¼ I1þ dm1 : Te2 ¼ I dm2 C1 v1 0 C2 v2 0 2
00
ð26Þ
The diffuse mismatch model gives: ð25Þ
2.3. Definition of temperatures and heat flux
q1 ¼ 2p
2189
DT12 q00 12
¼
2 2 ¼ T12 C1 v1 T21 C2 v2
ð31Þ
3. Numerical method The Equation of Phonon Radiative Transfer in each film with the associated initial and boundary conditions are solved numerically using the finite-difference method. These are integro-differential equations and may be solved using an iterative technique. Since the equations are time-dependent, therefore, the implicit scheme is used to solve them. Consider the EPRT in the first film: R 1 1 2 1 I1 d 1 I1
m
2.4. Thermal boundary resistance
1 @I1 @I1 þ m1 ¼ v1 @t @x1
At the interface of two dielectric films or a metal and a dielectric film there exists a resistance to heat flow, which is called the thermal boundary resistance. Due to this resistance, there is a small but finite temperature drop across the interface of
The variable m1 varies between 1 r m1 r1. To facilitate the solution of the EPRT it is customary to break the intensity into parts. When 1 r m1 r 0, I1 ¼ I1 ðx1 , m1 ,tÞ and when 0 r m1 r1, I1 ¼ I1þ ðx1 , m1 ,tÞ. Therefore, we may write separate equations for
L1
2190
S. Bin Mansoor, B.S. Yilbas / Physica B 406 (2011) 2186–2195
the two intensities. These are given below: R R1 þ 0 þ 1 @I1þ 1 I1 dm1 þ 0 I1 dm1 I1 2 1 @I1þ þ m1 ¼ , v1 @t @x1 L1
0 r m1 r 1
and @I 1 @I1 þ m1 1 ¼ v1 @t @x1
1 2
R 0
m
1 I1 d 1 þ
R1 0
I1þ dm1 I1
L1
Layer 1 Layer 2
,
and 1 2
R 0
m
1 I2 d 2 þ
R1 0
I2þ dm2 I2
L2
,
1r m2 r 0
For the intensity I2þ , the initial condition is I2þ ðx2 , m2 ,0Þ ¼ C2 v2 To = 4p and the boundaryn condition comes from the interface condiR1 R0 tion I2þ ð0, m2 ,tÞ ¼ 2 T12 0 I1þ ðL1 , m1 ,tÞm1 dm1 R21 1 I2 ð0, m2 ,tÞm2 dm2 g. For the intensity I2 , the initial condition is I2 ðx2 , m2 ,0Þ ¼ C2 v2 To =4p and the boundary condition is I2 ðL2 , m2 ,tÞ ¼ C2 v2 Tright =4p. Starting from the initial intensity distributions in each film, at each time-step two equations are solved in each film. This solution gives the approximate intensity distribution in each film. From this solution the interface temperature is calculated and the next iteration is performed. After convergence the computations proceed to the next time-step. Once the simulation is complete, the equivalent equilibrium phonon temperature may be evaluated at each time-step and at each location inside the two films. Table 1 gives the properties of the silicon and the diamond films used in the simulations. Table 2 gives the geometric configuration of the films and initial and boundary temperatures.
Table 1 Properties of silicon and diamond films used in the simulations.
Silicon Diamond
Length (mm)
To (K)
Tleft (K)
Tright (K)
0.1 0.1
300 300
301 –
– 300
1r m1 r 0
To solve the above equations an implicit method is used for the time discretisation and discrete ordinates method is used for the discretisation of the angular term. The time-derivative in each equation is approximated by means of the first-order backwarddifference. The partial derivative with respect to the coordinate variable x1 is descritised using a backward-difference for the first equation in which 0r m1 r1 and is descritised using a forwarddifference for the second equation in which 1r m1 r0. This procedure is important to obtain stable discretisations. The integrals are calculated by first discretising the cosine of the azimuthal angle m1 into 40 divisions and then using Simpson’s rule to evaluate the integrals numerically for each time-step and each coordinate location. It is important to keep in mind that the intensities I1þ and I1 are two separate variables and each needs an initial condition and a single boundary condition. For the intensity I1þ , the initial condition is I1þ ðx1 , m1 ,0Þ ¼ C1 v1 To =4p and the boundary condition is I1þ ð0, m1 ,tÞ ¼ C1 v1 Tleft =4p. For the intensity I1 , the initial condition is I1 ðx1 , m1 ,0Þ ¼ C1 v1 To =4p and the boundary ncondition R1 comes from the interface condition I1 ðL1 , m1 ,tÞ ¼ 2 R12 0 I1þ R0 ðL1 , m1 ,tÞm1 dm1 T21 1 I2 ð0, m2 ,tÞm2 dm2 g. Exactly similar considerations apply to the intensity calculations in the 2nd film. The EPRT is first written in terms of the two intensities as given below: R R1 þ 0 þ 1 @I2þ 1 I2 dm2 þ 0 I2 dm2 I2 2 1 @I2þ þ m2 ¼ , 0 r m2 r1 v2 @t @x2 L2
@I 1 @I2 þ m2 2 ¼ v2 @t @x2
Table 2 Geometric configurations of the film and initial and boundary temperatures of the films.
Thermal conductivity (W m 1 K 1)
Volumetric specific heat (J m 3 K 1)
Phonon speed (m s 1)
Mean free path (nm)
157 3320
1.653 106 1.815 106
8430 12,288
33.8 447
Table 3 Transmittance and reflectance phonon radiation determined from diffusive mismatch model. Transmittance
Si/Si Si/Di
1-2 0.5 0.6154
Reflectance 2-1 0.5 0.3846
1-2 0.5 0.3846
2-1 0.5 0.6154
Table 3 gives the transmittance and reflectance of phonon radiation at the film surfaces.
4. Results and discussion The transient analysis of phonon transport in dielectric thin films is carried out and temperature variation along the films is predicted. Silicon–silicon and silicon–diamond thin films are incorporated in the simulations. The thermal boundary resistance at the interface of the films is accounted in the analysis. To validate the code developed for EPRT, the steady state transport across a single silicon film is simulated in line with the previous study [1]. The simulation results for different film thickness are shown in Fig. 2. The findings are exactly the same as the previous results [1]. Moreover, the validation study is extended to include the transient phonon transport across the silicon film in accordance with the early study [6]. Dimensionless temperature variation along the film thickness is shown in Fig. 3 for different dimensionless time. The findings are the same as those reported in the early study [6], i.e. dimensionless temperature slip at diffused film surfaces and dimensionless temperature distributions for the previously reported results. Fig. 4 shows the dimensionless temperature distribution along the silicon–silicon films for three dimensionless heating durations. It should be noted that the dimensionless temperature is set at 1 at time t*Z0 and radiative phonon transport is initiated from the front surface of the silicon film (x*¼0). In addition, initially both films are considered to be in temperature equilibrium. The emitted, equilibrium and gray body temperatures are shown at the interface. The temperature slip occurs at the interface because of the transmission and the reflectance of the radiative phonon energy at the interface. Due to the transient nature of the phonon transport, temperature decay in the first film of silicon is sharp while it is gradual in the adjacent silicon film. In addition, the temperature slip at the interface increases with progressing time. This is attributed to the phonon scattering in the film and the phonon energy reaching at the interface, which is small during the early periods. As the time progresses, the radiative phonon energy reaching at the interface becomes high, which in turn increases the equilibrium temperature at the interface. Moreover, as the heating period progresses further, the steady state is reached and temperature variation along the film thickness becomes linear. This is true for both films. It should be noted that temperature in the film corresponds to the energy of all phonons at a local point in the film. Furthermore, although the local thermodynamic equilibrium condition is not satisfied in terms of a
S. Bin Mansoor, B.S. Yilbas / Physica B 406 (2011) 2186–2195
2191
Fig. 4. Dimensionless temperature with dimensionless thickness along silicon– silicon. thin films for different heating periods. Dimensionless Temperature is defined as (Tequilibrium–Tback)/(Tfront–Tback).
Fig. 2. Dimensionless temperature with dimensionless film thickness for steady state condition similar to that of Ref. [1]. Dimensionless temperature is defined as (Tequilibrium–Tback)/(Tfront–Tback).
Fig. 3. Dimensionless temperature with dimensionless film thickness for transient condition similar to that of Ref. [6]. Dimensionless temperature is defined as (Tequilibrium–Tback)/(Tfront–Tback).
thermal equilibrium state, the equilibrium temperature is considered for the phonons when they are redistributed adiabatically to an equilibrium state. Consequently, the equivalent equilibrium temperatures on both sides of the interface across the films are not the
same due to the reflection and transmission of phonons at the interface. Fig. 5 shows temporal variation of temperatures at both surfaces of the interface. The emitted phonon temperature represents the phonons that are emitted at each side of the interface at this temperature. Therefore, the emitted phonon temperatures are different than the equilibrium temperature of phonons. It should be noted that the equilibrium temperature is more appropriate for the predictions of thermal boundary resistance across the interface and represents the measurable temperature at both sides of interface [15]. Temperature at the interface of the front film (silicon) increases with progressing time. This is because of the reflected phonons from the interface towards the front film. Since the phonon transport is not diffusive in the film, the difference between the emitted temperature and the gray body temperature becomes large with progressing time. It should be noted that due to non-diffusive phonon transport in the film, temperature slip occurs at the boundary and the Fourier theory of heating fails to govern the phonon transport in the film. The temperature slip is associated with the phonons reaching at the interface, which is small in the early heating period. Moreover, the temperature slip is also noted at the interface of the second film (silicon), provided that the equilibrium temperature and the gray body temperatures are higher than of the emitted phonon temperature at the interface. This is associated with reflected phonons from the outer boundary of the second film towards the interface, since the emitted phonon temperature is related to phonons reaching at the interface when emitted from the free surface of the layers. In addition, the thermal boundary resistance lowers the equilibrium phonon temperature at the interface of the second layer. Since the gray body temperature is determined from the consideration of the diffusive boundary [1], it remains the same at both sides of the interface due to energy balance. Fig. 6 shows temporal variation of heat flux due to phonon transport at the front surface of the film and the back surface of the second film. It should be noted that the front surface of the first layer is considered to be at a higher temperature than the back surface of the second layer prior to heating is initiated (t*o0). Consequently, heat flux at the front surface of the first film is higher in the early heating period and, first, reduces sharply tr100 ps and decays gradually to reach a steady value. The opposite is true for the
2192
S. Bin Mansoor, B.S. Yilbas / Physica B 406 (2011) 2186–2195
Fig. 5. Dimensionless temperatures with time at the silicon–silicon interface. Dimensionless temperature is defined as (Tequilibrium–Tback)/(Tfront–Tback).
second layer. The rapid decay of heat flux at the surface of the first film is because of the high rate of phonon transport from surface region towards the bulk of the film during the early heating period. However, temperature rise in the film with progressing time (Fig. 4) lowers the amount of radiative phonon transport from this region to the bulk. Once, the steady temperature (Fig. 4) is achieved within both films, the radiative phonon transport within the film becomes steady, i.e. the heat flux remains the same. Fig. 7 shows temperature distribution in the silicon–diamond films for three heating periods. The first film is silicon and the second film is diamond. Temperature behavior in the silicon film is similar to that shown in Fig. 4, provided that interface temperature in the silicon film is lower than that corresponding to the silicon–silicon films. This is attributed to (i) phonon
transmittance and reflectance at the interface as given in Table 3 and (ii) the thermal boundary resistance determined from the diffuse mismatch model, which is different than that of the silicon films. Moreover, the large phonon mean free path of diamond lowers the phonon scattering in this film. Consequently, the transmitted phonons reflect from the back surface of the diamond film towards the interface. The attainment of low equilibrium phonon temperature at the silicon film interface is associated with the large phonon transmission from the silicon film to the diamond film. Fig. 8 shows temporal behavior of the emitted phonon temperature, the equilibrium phonon temperature and the gray body temperature at the interface of the silicon and diamond films. Temperature behavior is similar to those shown in Fig. 5.
S. Bin Mansoor, B.S. Yilbas / Physica B 406 (2011) 2186–2195
2193
Fig. 6. Heat flux with time at the front and back surfaces of silicon–silicon films.
However, the maximum temperatures at interface of both sides are lower than that corresponding to the silicon–silicon films. This is attributed to small scattering of phonons in the diamond film due to large phonon mean free path of diamond, and transmittance and reflectance at the interface. Moreover, temporal variation of heat flux is shown in Fig. 9 for the front surface of silicon film and the back surface of the diamond film. The flux at the back surface of diamond film is much higher than that of the front surface of silicon film. This is more pronounced as the heating period progresses. The attainment of high heat flux is associated with the phonon reaching at the back surface of the diamond film, which is large due to the large mean phonon mean free path of diamond. In addition, phonon transmission from the silicon interface contributes to the attainment of the high heat flux at the back surface of diamond.
5. Conclusion The phonon transport across the silicon–silicon films and silicon–diamond films is considered. The steady state analysis is carried out to validate the phonon temperature distribution in the silicon layer. The transient analysis is carried out after assuming the presence of temperature difference of 1 K across the films. The EPRT is solved numerically to obtain phonon transport in the films. The diffuse mismatch model is incorporated to formulate and determine the thermal boundary resistance across the interface of the films. It is found that phonon temperature variation across the silicon film predicted from the present study for steady and transient situations agree well with the previously reported results. In silicon–silicon films, temperature slip occurs at the interface, despite the fact that the phonon
2194
S. Bin Mansoor, B.S. Yilbas / Physica B 406 (2011) 2186–2195
Fig. 7. Dimensionless temperature with dimensionless thickness along silicon– diamond thin films for different heating periods. Dimensionless temperature is defined as (Tequilibrium–Tback)/(Tfront–Tback).
mean free path remains the same in both films. This is associated with the transmittance and reflectance of phonons at the interface of the films. The rise of the equilibrium, emitted and gray body temperatures is higher at the interface of the silicon films during the heating period 36rt r72 ps. This is because of the small mean free path of silicon; in which case, phonons reach the interface at higher rates. Temperature slip is higher at the interface of silicon–diamond films, which is associated with the thermal boundary resistance and the large phonon mean free path of diamond. In this case, the transmittance at interface increases and less scattering in the diamond film allows transmitted phonons to reflect from the back surface of the diamond to interface at a large amount due to small scattering. The heat flux increases at the front surface of silicon film and at the back surface of the diamond film with progressing time. This is more pronounced at the back surface of the diamond film due to the large number of phonon reaching at the back surface of diamond film because of less scattering of transmitted phonons in the diamond film.
Fig. 8. Dimensionless temperatures with time at the silicon–diamond interface. Dimensionless temperature is defined as (Tequilibrium–Tback)/(Tfront–Tback).
S. Bin Mansoor, B.S. Yilbas / Physica B 406 (2011) 2186–2195
2195
Fig. 9. Heat flux with time at the front and back surfaces of silicon–diamond films.
Acknowledgments The authors acknowledge the support of Center of Excellence for Scientific Research Collaboration with MIT and King Fahd University of Petroleum and Minerals, Dhahran, Saudi Arabia for this work. References
[6] [7] [8] [9] [10] [11] [12] [13]
[1] [2] [3] [4] [5]
A. Majumdar, ASME J. Heat Transfer 115 (1993) 7. A.A. Joshi, A. Majumdar, J. Appl. Phys. 74 (1) (1993) 31. G. Chen, ASME J. Heat Transfer 118 (1996) 539. T. Zeng, G. Chen, ASME J. Heat Transfer 123 (2001) 340. R.S. Prasher, P.E. Phelan, ASME J. Heat Transfer 123 (2001) 105.
[14] [15] [16]
L. Pilon, K.M. Katika, ASME J. Heat Transfer 126 (2004) 735. A. Majumdar, P. Reddy, Appl. Phys. Lett. 84 (23) (2004) 4768. L. Jiang, H. Tsai, ASME J. Heat Transfer 128 (2006) 926. K. Miyazaki, T. Arashi, D. Makino, H. Tsukamoto, IEEE Trans. Components Packag. Technol. 29 (2) (2006) 247. T. Kuhn, I.J. Maasilta, Nucl. Instrum. Methods Phys. Res. A 559 (2006) 724. S.V.J. Narumanchi, J.Y. Murthy, C.H. Amon, ASME J. Heat Transfer 126 (6) (2004) 946. C. Kittel, Introduction to Solid State Physics, sixth ed., Wiley, New York, 1986. W.G. Vincenti, C.H. Kruger, Introduction to Physical Gas Dynamics, Robert Krieger Publications, New York, 1977. M.N. Ozisik, Radiative Transfer and Interactions with Conduction and Convection, Wiley, New York, 1973. G. Chen, Phys. Rev. B 57 (23) (1998) 14958. E.T. Swartz, R.O. Pohl, Rev. Mod. Phys. 61 (3) (1989) 605.