- Email: [email protected]
Numerical simulation of low temperature thermal conductance of corrugated nanofibers
Physica E 44 (2012) 1189–1195
Contents lists available at SciVerse ScienceDirect
Physica E journal homepage: www.elsevier.com/locate/physe
Numerical simulation of low temperature thermal conductance of corrugated nanofibers ¨ a, L. Baskin b T. Puurtinen a,n, P. Neittaanmaki a b
Department of Mathematical Information Technology, University of Jyv¨ askyl¨ a, Finland Department of Mathematics, St. Petersburg State University for Telecommunications, Russia
a r t i c l e i n f o
abstract
Article history: Received 29 August 2011 Received in revised form 16 January 2012 Accepted 17 January 2012 Available online 25 January 2012
In this work we numerically investigated the decrease of low temperature thermal conductance of cylindrical Silicon nanofibers with periodical and almost periodical corrugation on the boundary. We used 1D and 3D finite element models to study scattering of low frequency acoustic phonons from the corrugation. By assuming weak deformation and high enough cross-sectional symmetry it was found that 1D models yield a sufficiently good approximation and could be used in the calculation of thermal conductance for long nanofibers. Landauer 1D formalism was applied for the calculation of the thermal conductance of several long axially symmetric and elliptical fibers. Without significantly hindering the structural strength of the fiber we found boundary shapes which decrease the thermal conductance by each mainly contributing phonon down to 10% of the universal quantum on temperature range 0.1–1 K. We concluded that the similar drop in total thermal conductance could be obtained by combining these designs in series. & 2012 Elsevier B.V. All rights reserved.
1. Introduction It is known that the process of heat conductivity renders an important influence on the stable operation of nanoscale electronic devices. For this reason the mechanism of heat transfer in nanostructures and nanoelements has come presently to notice. Heat conductivity, conditioned by acoustic phonons, is investigated in a number of nanostructures, such as superlattices [1–4], thin films [5–7], nanowires [8–15] and nanotubes [16]. Heat conductivity in quantum wires began to be studied especially intensively after the theoretical prediction of existence of universal quantum of heat conductivity [17,18] and its experimental discovery at very low temperatures [19]. Based on Green functions formalism Santamore and Cross investigated the influence of surface inhomogeneity on universal heat conductivity in a dielectric quantum wire at sub-Kelvin temperatures [20,21]. In Ref. [22] a contribution of longitudinal acoustic waves to heat conductivity was calculated taking into account the influence of the form of contact between a source of heat and the quantum wire. It was assumed that a contact had a catenoidal form. Heat conductivity was analyzed in structures with a split (Y-branch structures) [23], a T-shaped quantum wire [24] and a doublen Corresponding author. Present address: Nanoscience Center, Department of ¨ ¨ PL 35, 40014 Finland. Tel.: þ 358 50 367 5000; Physics, University of Jyvaskyl a, fax: þ 358 14 260 4756. E-mail address: tuomas.a.puurtinen@jyu.fi (T. Puurtinen).
1386-9477/$ - see front matter & 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.physe.2012.01.009
bend quantum waveguide [25]. In Refs. [26,27] low temperature heat conductivity was taken into consideration in nanowires with structural defects such as absent lattice atoms (voids) and compressions (clamped materials). The influence of boundary conditions on defects was analyzed in Ref. [27]. In a number of applications it is necessary to substantially (in several times) lower the heat conductivity. This necessity arises, for example, in the creation of thermo-insulating suspenders for the sensors of infrared radiation (bolometers) operating at ultralow temperatures. Sensitivity of such instruments is inversely proportional to heat conductivity. The analysis of the results of the mentioned works shows that thermal insulation can be enhanced with non-homogeneity of fiber border and with the designing of the forms of contact in a proper way. Noticeable changes can also take place in the case of the presence of defects in the volume of fiber due to the internal reflection of acoustic waves. However, it takes place only at an experimentally unrealizable condition where a surface of a defect is a hard wall. In addition, to obtain a complete picture of sub-Kelvin heat conductivity in such wires requires a systematic analysis of all three types of elastic waves present in low-energy regime. In Ref. [28] the process of transport of three types of the elastic waves (longitudinal, transversal and torsional) was investigated in a cylindrical fiber with deformed boundary. It was shown that the small periodic deformation of the boundary of the fiber, which only slightly alters the mechanical strength, forms conditions for the reflection of waves at resonant frequencies. In order to effectively
1190
T. Puurtinen et al. / Physica E 44 (2012) 1189–1195
lower the thermal conductivity of fiber at sub-Kelvin temperatures it is necessary to decrease transmission probability of waves in a wide range of frequencies ð0:2kB T=_ o o o5kB T=_Þ. To that end, it was suggested to monotonically change the period of border deformation (so-called ‘‘chirping’’) of the fiber. Because the mentioned three types of acoustic waves have different speeds of propagation, it is necessary, in principle, to create three areas of deformation, each designed for different intervals of wavelengths and for the reflection of different types of waves. In Ref. [28] the calculations of transmission probabilities of acoustic waves in chirp-deformed corrugation were executed in one-dimensional approximation. The purpose of this research is to check the validity of such approximation by comparing 1D and 3D models, and to analyze the influence of boundary design for thermal conductance in a nanofiber depending on temperature.
For a cylindrical fiber, there are four independent modes for which oj,min ¼ 0: a longitudinal mode, a torsional mode and two flexural modes. For all other modes oj,min ¼ ncj =D, where cj is the characteristic propagation speed of elastic wave j, D is the diameter of the rod, n ¼ 1; 2, . . .. At temperatures T o1 K and D o 100 nm only the mentioned four modes give the main contribution to O, so O ¼ 4Ouniv . It is known that at low temperatures an unstructured surface roughness of a wire can diminish quantum heat conductance due to the constructive interference of scattered waves [20]. The roughness structured in a proper way can sharply decrease the transmission probability pj ðoÞ. It is possible to design the surface of a wire in such a way that pj ðoÞ would be small for all phonons in a wide frequency band centered near o kB T=_, so that it would induce a substantial drop in total thermal conductance.
2. Physical background
3. 3D model of scattering in an elastic wire
In an elastic dielectric nanofiber heat energy is transferred by ! elastic waves. Such a wave with wavevector k and frequency ! oð k Þ can be represented as a flow of virtual particles known as phonons. Every phonon possesses energy E ¼ _o and crystal ! ! momentum p ¼ _ k . At temperatures T r1 K the phonon mean free path exceeds 100 mm (for defect-free fiber) and phonons propagate without collisions through the fiber. In other words, phonon flow is within the ballistic regime. Let reservoirs R1 and R2 with respective temperatures T1 and T 2 4T 1 be connected by a dielectric fiber G, operating in ballistic regime. Then the thermal energy passing through G from R2 to R1 in a unit time is [17] Z @Q 1 X ¼ _o½nðo,T 2 Þnðo,T 1 Þpj ðoÞ do @t 2p j
Let P ¼ fðx,y,zÞ A R3 : ðx,yÞ A Y, 1o z o þ 1g be an infinite cylinder whose cross-section Y R2 is a two-dimensional bounded and connected domain with smooth boundary. Assume that an elastic wire G of homogeneous and isotropic material coincides with P for 9z9 4L, and let
where nðo,T i Þ ¼ ðexpð_o=ðkB T i ÞÞ1Þ1
ð1Þ
is the Bose–Einstein distribution function; pj ðoÞ is the probability for a phonon of mode j to pass through G from R2 to R1; kB is the Boltzmann constant. Integration is over the whole spectrum of phonon frequencies. In the case 9T 2 T 1 9 5 T 2 ,T 1 we have Z @Q T 2 T 1 X @nðo,TÞ ¼ p j ð oÞ d o _o @t @T 2p j The value
O¼
1 @Q T 2 T 1 @t
is called the thermal conductance. Taking into account Eq. (1) we obtain the 1D Landauer formula for phonons [17]: Z _2 X oj,max o2 expðb_oÞ O¼ pj ðoÞ do ð2Þ 2 ½expðb_oÞ12 2pkB T j oj,min where b ¼ ðkB TÞ1 ; T T 1 T 2 and oj,min ðoj,max Þ is the minimal (maximal) phonon frequency for mode j. If for mode j the values pj ðoÞ ¼ 1, oj,min ¼ 0 and oj,max ¼ 1 the corresponding contribution in O is
Ouniv ¼
GðLÞ ¼ fðx,y,zÞ A G : 9z9 r Lg be a truncated part of it (see Fig. 1). For computational purposes parameter L4 0 is chosen large enough so that there are sufficiently long straight cylinders (with cross-section Y) contained in G(L). Propagation of 3D phonons in wire G can be investigated by the homogeneous ‘‘frequency domain’’ boundary value problem of elasticity theory: ! L u ðx,y,zÞ ¼ 0, ðx,y,zÞ A G ! n sð! u Þðx,y,zÞ ¼ 0, ðx,y,zÞ A @G ! where u ¼ ðu1 ,u2 ,u2 Þ stands for displacements and
ð3Þ
! ! ! ! L u ¼ mD u þ ðl þ mÞr div u þ ro2 u
ð4Þ
Lame´’s constants l and m satisfy m 4 0 and 3l þ 2m 4 0, and density r together with the frequency o is supposed to be ! constant. Vector n is the outward unit normal on @G. Finally, sðÞ is the stress tensor with Cartesian components
sð! u Þ ¼ ðlekk dij þ 2meij Þij where eij ¼ 12 ð@i uj þ @j ui Þ, i,j ¼ 1; 2,3. By setting the stress tensor to zero one considers a free boundary problem, i.e. no outer mechanical forces affect the boundary. A solution Y to problem (3) is called a continuous spectrum eigenfunction (CSE) corresponding to the frequency o 4 0, if Y satisfies the following conditions: Y is a bounded function in G : 9Yðx,y,zÞ9 rC, ðx,y,zÞ A G, C constant; and Y does not converge to zero at infinity: 9Yðx,y,zÞ9Q0 as 9ðx,y,zÞ9-1. The CSEs play the role of phonon wave functions. It is known that the total number NðoÞ of linearly independent CSEs below eigenfrequency o increases in steps as o increases, so that it is a piece-wise constant function. The step values oi of N are usually called ‘‘threshold’’ values. We assume here to consider only frequencies
p k2B T 6 _
For ballistic regime of thermal conductivity this universal value is independent of the properties of a fiber material and its geometry. If _oj,min bkB T, then the contribution of mode j in thermal conductivity is exponentially small and can be neglected.
Fig. 1. Illustration of a cylindrical fiber located on the z-axis with axially symmetric perturbation near origin. G(L) is used as a computational domain for scattering of elastic waves.
T. Puurtinen et al. / Physica E 44 (2012) 1189–1195
below the first threshold. It is known that for problem (4) the number of linearly independent CSEs is eight. In the space of CSEs in G corresponding to a frequency o, a basis fY 1 ð, oÞ, . . . ,Y 8 ð, oÞg can be chosen such that 8 X
Y j ð, oÞ ¼ wjþ ð, oÞ þ
g9z9 sjk ðoÞw Þ k ð, oÞ þ Oðe
ð5Þ
k¼1
where ‘’ represents a point ðx,y,zÞ in G, j ¼ 1, . . . ,8 and g is a sufficiently small positive number. The existence of such solutions was discussed e.g. in Ref. [30]. Functions wj7 :¼ wj expðaj ikzÞ ZR ðbj zÞ are called incoming ( þ) and outgoing ( ) waves. They can be explicitly described by choosing correct signs aj ¼ 71 and bj ¼ 7 1. However, depending on the cross-sectional shape Y, functions wj :¼ wj ðx,y,zÞ do not usually possess an analytical form. We can choose the numbering in such a way that wj7 with j ¼ 1, . . . ,4 ðj ¼ 5, . . . ,8Þ have a support on negative (positive) z-half space, respectively. If Y is a circle, the analytical expression exists and functions wj can be identified as flexural (with two polarizations), torsional and longitudinal waves. Here ZR is any smooth cutoff function with the properties ZR ðz ZRÞ ¼ 1, and ZR ðz r REÞ ¼ 0 for 0 o E o R. Asymptotic behavior of Yj is independent of the choice of ZR . Finally, the matrix SðoÞ ¼ ðsjk ðoÞÞ8j,k ¼ 1 is unitary and it is known as the scattering matrix. Let us examine the solutions Yj in Eq. (5) for j ¼ 1, . . . ,4. Solution Yj describes the phonon incoming from 1 as the wave wjþ , which scatters from the perturbed part of the wire G located near the origin into eight outgoing waves w k . With such a numbering convention the sum pj ðoÞ :¼
8 X
2
9sjk ðoÞ9
k¼5
is called the transmission coefficient. It can be interpreted as the transmission probability for the phonon Yj through the wire G P 2 from 1 to þ1. Finally, 1pj ðoÞ ¼ 4k ¼ 1 9sjk ðoÞ9 is the reflection coefficient.
Non-zero real functions Bx(z) and By(z) are the elliptical semiaxes. We also assume that the cross-section is a circle of constant radius outside of the corrugated part GðlÞ :¼ G \ f9z9 r lg (see Fig. 1), i.e. Bx ðzÞ ¼ By ðzÞ for 9z9 Z l. We can now use three simple 1D ordinary differential equation models to approximate various 3D waves. The equations can be derived e.g. from physical basis [29]. For longitudinal modes we use the equation: 1 @ @uðzÞ r AðzÞ ¼ o2 uðzÞ ð6Þ AðzÞ @z @z E and for torsional modes the equation: " 3 # Bx ðzÞB3y ðzÞ @uðzÞ @ r Bx ðzÞBy ðzÞðB2x ðzÞ þ B2y ðzÞÞo2 uðzÞ ¼ @z B2x ðzÞ þ B2y ðzÞ @z 4m
In various situations a 3D system can be approximately described by means of 1D model. Such a situation arises in our case for wave propagation in a fiber when the fiber length as well as the wavelength is much greater than the diameter of the fiber. In this case the equation describing a 3D model can be approximately replaced by equations corresponding to various 1D models. Such approaches allow getting description of the process which is adequate over limited range of parameters. The limits of the range will be clarified with a reasonable accuracy by numerical experiments. The low temperature thermal conductance of elastic dielectric fibers is mainly determined by waves which have wavelengths l about 0:121 mm at sub-Kelvin temperatures. Thus, for nanofibers with diameter 10–20 nm the passage is possible from 3D to 1D model. The transition from 3D to 1D is done separately for each mode. It is known (see e.g. Ref. [29]) that in a cylinder each of the elastic modes j with oj,min ¼ 0 is mainly one-dimensional in their nature, i.e. each mode has major displacement component in only one direction in some coordinate system. In cylindrical fiber, flexural modes have major displacement in the x, y-directions, longitudinal in the z-direction and torsional in the y-direction in cylindrical coordinate system. We henceforth assume that the cross-section of 3D fiber G is a circle or an ellipse, so that G can be defined by G ¼ fðx,y,zÞ A R3 : x2 =B2x ðzÞ þ y2 =B2y ðzÞ r1g
ð7Þ
To simulate flexural phonon modes with polarization in the x-direction we use the equation: 2 @2 4r 3 @ uðzÞ Bx ðzÞBy ðzÞo2 uðzÞ B ðzÞB ðzÞ ð8Þ ¼ x y E @z2 @z2 Here r is a density, E ¼ mð3l þ2mÞ=ðl þ mÞ is Young modulus and o is the angular frequency of the wave. In Eq. (6) AðzÞ ¼ pBx ðzÞBy ðzÞ is the area of cross-section of 3D domain at z. For y-polarized flexural waves indices x and y in the above equation are interchanged.
5. Computational method A finite part G(L), L4 l, of fiber G is taken as a computational domain. In order to solve the complex scattering amplitudes sjk of þ outgoing waves w k for incoming wave wj in Eq. (5) we simplify the calculation by further assuming that the perturbation on G(L) should have a cross-sectional symmetry in such a way that there is no mixing of normal mode types after scattering. In other words, ansatz Y j ¼ wj expðikzÞZR ðzÞ þ sjj wj expðikzÞZR ðzÞ þsj,j þ 4 wj expðikzÞZR ðzÞ þ Oðeg9z9 Þ
4. 1D models of scattering in an elastic wire
1191
in GðLÞ
ð9Þ
is valid in the domain G(L). Here wj :¼ wj ðx,y,zÞ and ZR is any smooth cut-off function with the properties ZR ðz Z RÞ ¼ 1 and ZR ðz r REÞ ¼ 0 as in Section 3. Number k4 0 is the wavenumber and it depends on the angular frequency o and the type of the phonon j. The exact calculation of frequency relation oj ðkÞ for circular cylinders was conducted in Ref. [29]. However, we decided to use 1D approximations of this relation also for 3D model, because at low frequencies such an approximation was found to provide a sufficient numerical accuracy. Inverse relations kj ðoÞ are sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffi rffiffiffiffi rffiffiffiffi r r o 4r ð10Þ , ktors ðoÞ ¼ o klong ðoÞ ¼ o , kflex ðoÞ ¼ E R E m For flexural waves the constant R :¼ Bx ð 7 LÞ ¼ By ð 7 LÞ is the radius of the cylindrical fiber at the ends. We use the following convention to number the phonon modes: j ¼ 1; 2 are the two flexural modes, j¼3 is the torsional mode and j¼4 is the longitudinal mode. We now explain the choice of boundary conditions on the cylindrical ends Gð 7 Þ :¼ fðx,y,zÞ : z ¼ 7 Lg \ GðLÞ of the domain. They will be used to numerically solve the amplitudes sjj and sj,j þ 4 needed in the calculation of transmission probabilities. Let us first take a look at boundary Gð þ Þ on positive z-axis where there is no incoming wave present: Y j wj expðikzÞ. We therefore require ð@z þ ikÞY j ¼ 0
on Gð þ Þ
1192
T. Puurtinen et al. / Physica E 44 (2012) 1189–1195
Now at the boundary on negative z-axis we must allow both incoming and outgoing waves of the same type to exist, so that Y j wj expðikzÞ wj expðikzÞ. We therefore put ð@z þ ikÞY j ¼ 2ikwj expðikzÞ
on G
ðÞ 2
Transmission and reflection probabilities pj :¼ 9sj,j þ 4 9 and 2 pj,ref :¼ 9sjj 9 can be calculated from the elastic field Yj which is a solution to frequency domain problem: LY j ¼ 0 in GðLÞ ! n sðY j Þ ¼ 0 on @GðLÞ\ðGð þ Þ [ GðÞ Þ on Gð þ Þ
ð@z þ ikÞY j ¼ 0
ð@z þ ikÞY j ¼ 2ikwj expðikzÞ
on GðÞ
ð11Þ
where L can be any of the operators defined by Eqs. (4), (6), (7) or (8) and k ¼ kj ðoÞ is the inverse 1D frequency equation of the corresponding wave type given by Eq. (10). Note that in 1D case the domain G(L) reduces to ½L,L and the free boundary condition in problem (11) is not needed. Transmission coefficient pj ðoÞ is calculated by evaluating the ! ! field value Y j ð z Þ on a point z ¼ ð½x,y,zÞ near the boundary Gð þ Þ . By taking squared absolute values from both sides of Eq. (9) and ! 2 by normalizing with 9wj ð z Þ9 : 2
pj ðoÞ ¼ 9sj,j þ 4 9
! 2 9Y j ð z Þ9 2 ! 9wj ð z Þ expðikzÞ9 2
2
Generally pj,ref ¼ 9sjj 9 ¼ 19sj,j þ 4 9 holds, but in order to observe possible numerical error we calculate the reflection coefficient ! separately. We therefore choose a point z near GðÞ and calculate ! from Y j ð z Þ wj ðx,y,zÞ expðikzÞ þ sjj wj ðx,y,zÞ expðikzÞ 2
pj,ref ¼ 9sjj 9
2 ! ! 9Y j ð z Þwj ð z Þ expðikzÞ9 2 ! 9wj ð z Þ expðikzÞ9
Even in the simplest fiber geometries functions wj can have quite complicated analytical forms (see e.g. Ref. [29]), so in order to further simplify the analysis we use approximations for them as well. We found the following forms sufficient (only small discrepancy with numerically solved eigenmodes) at low frequencies and small amplitudes. For longitudinal wave wlong ðx,y,zÞ ð0; 0, EÞ, for x-polarized flexural wave wflex ðx,y,zÞ ¼ ðE,0,iEkxÞ and for torsional wave wtors ðx,y,zÞ ¼ ððx cosðEÞy sinðEÞxÞ,ðx sinðEÞ þ y cosðEÞyÞ,0Þ where E 4 0 is small. We also remark that for different modes ! different evaluation points z ¼ ðx,y,zÞ must be chosen. In 3D model, longitudinal and flexural modes are evaluated on the z-axis while torsional modes are evaluated near the cylinder boundary where displacement is large. All three 1D models and 3D model were solved with finite element method (FEM). Implementation was done in Comsol v3.5 and Matlab 2008b using coefficient form PDE application mode for 1D models and structural mechanics application mode for 3D Lame´ system. For each elastic mode we put suitable boundary conditions on the cylindrical ends to transform problem (11) ! ! ! into a standard form used in Comsol: n cr u þq u ¼ g, where cr ¼ sðÞ is the stress tensor. To solve the arising FEM problem we used swept mesh generation in Comsol, which sweeps a 2D mesh from GðÞ to Gð þ Þ by constant steps. This method provided numerically reliable solutions with ease, as much finer mesh was needed if free mesh generation was used instead. Especially with flexural 3D modes the stability of solutions was heavily dependent on the mesh used. In the simulations second-order elements were used and the number of degrees of freedom was usually around 60,000 depending on the fiber geometry. For the linear equation we used direct linear solver Spooles v2.2.
6. Periodic structures of finite length It is known that infinite periodic structures can possess spectral gaps in frequency spectrum which prevent propagation of waves with such frequencies. It was numerically demonstrated in e.g. Ref. [28] that as the number of periods increases the transmission probability curve converges to a characteristic function of spectral gaps of similar infinite periodical structure. Therefore an effective phonon reflection over a thin frequency band can be established already with a small number (5–20) of periods. We used a simple harmonic periodical structure to test the suitability of 1D models for calculation of transmission probabilities. Assume that the radius of circular fiber along the z-axis is defined by n o a BðzÞ ¼ R 1 ½1 þ signðl9z9Þ½1 þ cosðLðzÞÞ ð12Þ 4 where R 40 is the radius of fiber, 0 r a o 1 is perturbation strength and l, 0 ol o L, is the half-length of perturbation. For a strictly periodical perturbation function LðzÞ ¼ pððz þ lÞn=l þ 1Þ is linear. Number n stands for the number of perturbations and 2l/n then becomes the length of one period in the perturbation. Total length of the fiber is 2L and total length of the perturbed part 2l (see Fig. 2). In Fig. 3 there are transmission probabilities of 1D and 3D torsional and longitudinal waves in a cylindrical Silicon fiber with a finite periodical deformation. It can be seen that 1D and 3D models give a very good agreement in the results with only small increasing error as the frequency o increases. This error can be explained as a numerical error arising from a fact that wavelength l approaches the element length in our FEM-system and it can be negated by increasing the mesh size. In Fig. 4 we show the transmission probability of flexural waves from 1D and 3D models in a Silicon fiber with finite harmonic perturbation on boundary. It is again seen that at low frequencies 1D model gives a fairly good approximation but the error slightly increases as o increases. Differences can be again explained as an error arising from diminishing of wavelength but also with the fact that 1D approximation of kflex ðoÞ is not as good for flexural waves as it is for torsional and longitudinal waves. Even that the 1D equations cannot model scattering into different mode types, such effects (if present) would induce noticeable changes in 3D transmission probabilities so that there would be significant discrepancy in 1D and 3D transmission probability curves. Therefore the above results predict that 1D models (6), (7) and (8) can provide sufficiently reliable transmission probability curves for the calculation of thermal conductance OðTÞ in high cross-sectional symmetry geometries with tens or hundreds of periods of perturbation.
Fig. 2. Illustration of a typical function B, which was used to define the crosssectional radius of cylindrical fiber G with radius R ¼ 5 nm and total length 2L ¼ 240 nm. There are eight periods of harmonic perturbation on the boundary with wavelength 20 nm and amplitude 0:251R.
1
1
0.8
0.8
Ω/Ωuniv, ptors, plong
ptors, plong
T. Puurtinen et al. / Physica E 44 (2012) 1189–1195
0.6
0.4
0.2
1193
0.6
0.4
0.2
0 0
100
200
300
400
500
600
700
0
800
0
ω [108 s−1]
2
4
6
8
10
T [K], X [hω/(2πkBT)]
Fig. 3. Transmission probability of longitudinal and torsional waves in cylindrical Silicon fiber with finite periodical perturbation on the boundary. Black lines correspond to torsional and grey lines to longitudinal waves. Dashed lines are calculated from 1D models and solid line from 3D model. Parameters for geometry: radius R ¼ 10 nm, number of periods n¼ 4, period length 2l=n ¼ 400 nm, perturbation strength a ¼ 0:251 and total fiber length 2:24 mm.
Fig. 5. Thermal conductance by torsional and longitudinal 1D waves in Silicon fiber with chirp. Solid lines correspond to thermal conductance normalized by the universal quantum and dashed lines correspond to scaled transmission probabilities ptors ðXÞ (black) and plong ðXÞ (dark grey) where X ¼ _o=ðkB TÞ and T ¼ 1 K. Light grey dash-dotted line is the scaled phonon distribution D(X). Parameters for geometry: radius R ¼ 5 nm, total fiber length 107 mm, number of periods n ¼80, first period length l1 ¼ 30 nm, period ratio C¼ 1.0664, chirp ratio C¼ 160 and perturbation strength a ¼ 0:251. Inset picture shows the periodicity cell.
1
1 0.8
Ω/Ωuniv, pflex
pflex
0.8 0.6
0.4
0.6
0.4
0.2
0.2 0 0
100
200
300
400 8 −1
500
600
700
0
ω [10 s ]
0
Fig. 4. Transmission probability of flexural waves in Silicon fiber with periodical perturbation. Dashed line corresponds to 1D model and solid line to 3D model. Parameters for geometry: radius R ¼ 10 nm, number of periods n¼ 4, period length 2l=n ¼ 100 nm, perturbation strength a ¼ 0:5 and total fiber length 960 nm.
7. Almost periodic structures To broaden the reflection band caused by periodic perturbation we change the periodicity slowly. This structure is called a chirp. The underlying idea is that chirp looks locally strictly periodical and each of these locally periodical parts contributes a thin reflection band at its resonant frequency. We now explain the mathematical formulation of such structure. Let the perturbation on the boundary consist of n periods of harmonic deformation with lengths lj , j ¼ 1, . . . ,n. Assume that lj þ 1 =lj ¼: d 4 1 is constant for all j and call it a period ratio. Now the radius of the fiber is defined by Eq. (12), but instead of linear LðzÞ we put
LðzÞ ¼ p 1 þ 2 logd 1 þ
ðzlÞðd1Þ l1
ð13Þ
2
4
6
8
10
T [K], X [hω/(2πkBT)] Fig. 6. Thermal conductance by flexural 1D waves in Silicon fiber with elliptical chirp. Solid line corresponds to thermal conductance normalized by the universal quantum and dashed line corresponds to scaled transmission probability pflex ðXÞ, where X ¼ _o=ðkB TÞ and T ¼ 1 K. Light grey dash-dotted line is the scaled phonon distribution D(X). Parameters for geometry: radius R ¼ 5 nm, total fiber length 4:3 mm, number of periods n¼50, first period length l1 ¼ 20 nm, period ratio d¼ 1.04, chirp ratio C¼ 6.8 and perturbation strength a ¼ 0:7. Inset picture shows the periodicity cell.
where l1 is the length of the first period (when looking from GðÞ ) and l¼
n1 1X j l1 d 2j¼0
is the haft-length of the chirp. We also mark ln =l1 ¼: C and call it the chirp ratio. In Fig. 5 we show the thermal conductance by torsional and longitudinal 1D waves (solid lines) through a long chirp with 80 periods of perturbation. According to simulations thermal conductance by these waves can decrease down to 10% of the
1194
T. Puurtinen et al. / Physica E 44 (2012) 1189–1195
The appearance of multiple gaps can be utilized in reflecting torsional and longitudinal waves as well. When the boundary of the cylindrical perturbation is defined as a superposition of harmonic functions of different period length (as in Fourier expansion), each of the present functions contributes to reflection at their respective frequencies. If the boundary is chosen to locally obey triangle wave shape (the chirp is built of truncated cones) we can see from Fig. 7 that thermal conductance is reduced on wider temperature range when compared to similar harmonic perturbation with matching fundamental component of Fourier expansion of each period.
1
Ω/Ωuniv
0.8
0.6
0.4
8. Conclusions
0.2
0 0
0.5
1
1.5
2
2.5
3
3.5
4
T [K] Fig. 7. Thermal conductance by torsional 1D waves in Silicon fibers with conical (triangular) axially symmetric chirp and similar harmonic chirp. Black line corresponds to thermal conductance of harmonic chirp and grey line to conical chirp. Parameters for geometry: radius R ¼ 5 nm, total fiber length 28:5 mm, number of periods n ¼30, first period length l1 ¼ 100 nm, period ratio d ¼ 1.0826, chirp ratio C¼ 10 and perturbation strength a ¼ 0:251.
universal quantum Ouniv on temperature range 0.1–1 K. In the figure there are also scaled transmission probability curves pj(X) for torsional and longitudinal waves, where X ¼ _o=ðkB TÞ and T ¼ 1 K is constant and, for comparison, the phonon distribution DðXÞ ¼
X 2 expðXÞ expðXÞ1
from the 1D Landauer formula (2). It can be seen that transmission probability curves are close to zero on 0:1 o X o 5 where there is most of the weight of the distribution D(X). We also found that harmonic perturbation with strength a ¼ 0:251 is not very effective in reflecting flexural waves. To reach 80–90% drop in thermal conductance for temperature range 0.1–1 K we had to use 350 periods in chirp with d ¼1.005, l1 ¼ 15 nm and R ¼ 5 nm. To reflect flexural waves efficiently one should be able to use stronger perturbation without deteriorating the tensile strength of the fiber. One possible solution to this is to use elliptical deformation where cross-sectional area is kept constant. Let cðzÞ ¼ 12 ð1 þsignðsin zÞÞ and let SðzÞ ¼ cðzÞð1 þ a sin zÞ þ cðz þ pÞ
1 1 þ a sinðz þ pÞ
where a is the strength of the perturbation. Periodic function SðÞ has a property SðzÞ ¼ Sðz þ pÞ1 . We now define Bx ðzÞ ¼
R 1signðl9z9Þ þð1 þ signðl9z9ÞÞSð2pLðz þ lÞÞ 2
where function LðÞ is given by Eq. (13). Finally by setting By ðzÞ ¼ RðBx ðzÞ=RÞ1 the cross-sectional area of an ellipse with semi-axes Bx and By is A :¼ pBx ðzÞBy ðzÞ ¼ pR2 . In Fig. 6 we show the thermal conductance by flexural 1D waves (solid line) through a circular fiber with elliptically deformed chirp. It can be seen that normalized thermal conductance goes down to (0.05–0.1)Ouniv on temperature range 0.1–1 K. Scaled transmission probability curve does not approach unity when X o 10, because elliptical (as well as circular) deformation opens multiple gaps in frequency spectrum of flexural waves, which are located at higher frequencies. The effect of higher gaps on hindering the wave propagation is reduced due to the fact that gap width is lower.
We investigated the possibility of decreasing the thermal conductance of Silicon nanofibers by designing the shape of the boundary to reflect propagating acoustic phonons. In cylindrical nanofibers of defect-free dielectric material only the four lowest energy acoustic phonons (longitudinal, torsional and two flexural waves) mainly contribute to thermal flow in temperatures 0.1–1 K. Propagation and scattering of these phonons can be studied within elasticity theory. We numerically simulated with 1D and 3D finite element models propagation of elastic waves in cylindrical nanofibers with periodical and almost periodical deformation on the boundary by calculating the transmission probability curves. We compared the curves obtained from 1D and 3D models for various simple corrugations and suggested designs where transition from 3D to 1D is possible. Finally, thermal conductance by each of the four low frequency phonons was calculated with 1D Landauer formula for long nanofibers. It was found that by assuming high-enough cross-sectional symmetry (to neglect phonon scattering into separate modes) and weak enough perturbation 1D approximation was reliable and the differences in transmission probability curves resulted into negligible discrepancies in thermal conductance curves. In order to effectively decrease thermal conductance of the fiber without deteriorating the structural strength we designed a long locally periodical corrugation (chirp) for each phonon. For torsional and longitudinal modes we used circular cross-section where the radius is varied locally periodically along the corrugation. For flexural modes such construction was found inefficient. To reflect flexural modes we designed an elliptical corrugation where the cross-sectional area is kept constant. It was found that with these designs the thermal conductance could be reduced down to 10% of the universal quantum on wide temperature range 0.1–1 K individually for all phonons. We conclude that a similar reduction in total thermal conductance could be obtained by combining such corrugations in series. However, with harmonic and elliptical fiber designs in series the total length of the corrugation can approach the reference phonon mean free path 100 mm. Even greater decreases in thermal conductance could be obtained using longer corrugations. The total fiber length can be reduced not only by designing the corrugation for narrower temperature range, but also by choosing a different period shape for chirp, which could simplify the possible experimental measurements in future. To this end, we finally compared thermal conductances by torsional waves in similar axially symmetric harmonic and conical periodical chirp structures. It was found that the temperature zone with reduced thermal conductance could be slightly expanded using a conical chirp design. Advantage of such a structure is that it may be easier to fabricate by etching processes. References [1] S.M. Lee, D.G. Cahill, R. Venkatasubramanian, Applied Physics Letters 70 (1997) 2957. [2] G. Chen, Physical Review B 57 (1998) 14958.
T. Puurtinen et al. / Physica E 44 (2012) 1189–1195
[3] W.E. Bies, R.J. Radtke, H. Ehrenreich, Journal of Applied Physics 88 (2000) 1498. [4] M.V. Simkin, G.D. Mahan, Physical Review Letters 84 (2000) 927. [5] X.Y. Yu, G. Chen, A. Verma, J.S. Smith, Applied Physics Letters 67 (1995) 3554. [6] A. Balandin, K.L. Wang, Physical Review B 58 (1998) 1544. [7] J. Philip, et al., Journal of Applied Physics 93 (2003) 2164. [8] T.S. Tighe, J.M. Worlock, M.L. Roukes, Applied Physics Letters 70 (1997) 2687. [9] A. Balandin, Physics of Low-Dimensional Structures 1–2 (2000) 1. [10] G. Chen, International Journal of Thermal Sciences 39 (2000) 471. [11] G. Chen, T. Zeng, Microstructure Thermophysical Engineering 5 (2001) 71. [12] B.A. Glavin, Physical Review Letters 86 (2001) 4318. [13] S. Volz, D. Lemonnier, J.B. Saulnier, Microstructure Thermophysical Engineering 5 (2001) 191. [14] J. Zou, A. Balandin, Journal of Applied Physics 89 (2001) 2932. [15] W.-Q. Huang, et al., Physics Letters A 336 (2005) 245. [16] P. Kim, L. Shi, A. Majumdar, P.L. McEuen, Physical Review Letters 87 (2001) 215502.
[17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]
1195
L.G.C. Rego, G. Kirczenow, Physical Review Letters 81 (1998) 232. M.P. Blencowe, Physical Review B 59 (1999) 4992. K. Schwab, E.A. Henriksen, J.M. Worlock, M.L. Roukes, Nature 404 (2000) 974. D.H. Santamore, M.C. Cross, Physical Review Letters 87 (2001) 115502. D.H. Santamore, M.C. Cross, Physical Review B 63 (2001) 184306. X. Lu, J.H. Chu, W.Z. Shen, Journal of Applied Physics 93 (2003) 1219. W.-X. Li, et al., Applied Physics Letters 85 (2004) 822. W.-X. Li, et al., Journal of Physics: Condensed Matter 16 (2004) 5049. W.-Q. Huang, et al., Journal of Applied Physics 98 (2005) 93524. K.-Q. Chen, et al., Physical Review B 72 (2005) 45422. L.-Y. Nie, et al., Physics Letters A 359 (2006) 234. T. Puurtinen, Numerical Simulation of Low Temperature Thermal Conduc¨ tance of Corrugated Nanofibers, Jyvaskyl a¨ Studies in Computing, vol. 124, ¨ ¨ 2010. University of Jyvaskyl a, [29] K.F. Graff, Wave Motion in Elastic Solids, Dover, New York, 1991. [30] S.A. Nazarov, B.A. Plamenevsky, Elliptic Problems in Domains with Piece-wise Smooth Boundaries, Walter de Gruyter, Berlin, 1994.
Contents lists available at SciVerse ScienceDirect
Physica E journal homepage: www.elsevier.com/locate/physe
Numerical simulation of low temperature thermal conductance of corrugated nanofibers ¨ a, L. Baskin b T. Puurtinen a,n, P. Neittaanmaki a b
Department of Mathematical Information Technology, University of Jyv¨ askyl¨ a, Finland Department of Mathematics, St. Petersburg State University for Telecommunications, Russia
a r t i c l e i n f o
abstract
Article history: Received 29 August 2011 Received in revised form 16 January 2012 Accepted 17 January 2012 Available online 25 January 2012
In this work we numerically investigated the decrease of low temperature thermal conductance of cylindrical Silicon nanofibers with periodical and almost periodical corrugation on the boundary. We used 1D and 3D finite element models to study scattering of low frequency acoustic phonons from the corrugation. By assuming weak deformation and high enough cross-sectional symmetry it was found that 1D models yield a sufficiently good approximation and could be used in the calculation of thermal conductance for long nanofibers. Landauer 1D formalism was applied for the calculation of the thermal conductance of several long axially symmetric and elliptical fibers. Without significantly hindering the structural strength of the fiber we found boundary shapes which decrease the thermal conductance by each mainly contributing phonon down to 10% of the universal quantum on temperature range 0.1–1 K. We concluded that the similar drop in total thermal conductance could be obtained by combining these designs in series. & 2012 Elsevier B.V. All rights reserved.
1. Introduction It is known that the process of heat conductivity renders an important influence on the stable operation of nanoscale electronic devices. For this reason the mechanism of heat transfer in nanostructures and nanoelements has come presently to notice. Heat conductivity, conditioned by acoustic phonons, is investigated in a number of nanostructures, such as superlattices [1–4], thin films [5–7], nanowires [8–15] and nanotubes [16]. Heat conductivity in quantum wires began to be studied especially intensively after the theoretical prediction of existence of universal quantum of heat conductivity [17,18] and its experimental discovery at very low temperatures [19]. Based on Green functions formalism Santamore and Cross investigated the influence of surface inhomogeneity on universal heat conductivity in a dielectric quantum wire at sub-Kelvin temperatures [20,21]. In Ref. [22] a contribution of longitudinal acoustic waves to heat conductivity was calculated taking into account the influence of the form of contact between a source of heat and the quantum wire. It was assumed that a contact had a catenoidal form. Heat conductivity was analyzed in structures with a split (Y-branch structures) [23], a T-shaped quantum wire [24] and a doublen Corresponding author. Present address: Nanoscience Center, Department of ¨ ¨ PL 35, 40014 Finland. Tel.: þ 358 50 367 5000; Physics, University of Jyvaskyl a, fax: þ 358 14 260 4756. E-mail address: tuomas.a.puurtinen@jyu.fi (T. Puurtinen).
1386-9477/$ - see front matter & 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.physe.2012.01.009
bend quantum waveguide [25]. In Refs. [26,27] low temperature heat conductivity was taken into consideration in nanowires with structural defects such as absent lattice atoms (voids) and compressions (clamped materials). The influence of boundary conditions on defects was analyzed in Ref. [27]. In a number of applications it is necessary to substantially (in several times) lower the heat conductivity. This necessity arises, for example, in the creation of thermo-insulating suspenders for the sensors of infrared radiation (bolometers) operating at ultralow temperatures. Sensitivity of such instruments is inversely proportional to heat conductivity. The analysis of the results of the mentioned works shows that thermal insulation can be enhanced with non-homogeneity of fiber border and with the designing of the forms of contact in a proper way. Noticeable changes can also take place in the case of the presence of defects in the volume of fiber due to the internal reflection of acoustic waves. However, it takes place only at an experimentally unrealizable condition where a surface of a defect is a hard wall. In addition, to obtain a complete picture of sub-Kelvin heat conductivity in such wires requires a systematic analysis of all three types of elastic waves present in low-energy regime. In Ref. [28] the process of transport of three types of the elastic waves (longitudinal, transversal and torsional) was investigated in a cylindrical fiber with deformed boundary. It was shown that the small periodic deformation of the boundary of the fiber, which only slightly alters the mechanical strength, forms conditions for the reflection of waves at resonant frequencies. In order to effectively
1190
T. Puurtinen et al. / Physica E 44 (2012) 1189–1195
lower the thermal conductivity of fiber at sub-Kelvin temperatures it is necessary to decrease transmission probability of waves in a wide range of frequencies ð0:2kB T=_ o o o5kB T=_Þ. To that end, it was suggested to monotonically change the period of border deformation (so-called ‘‘chirping’’) of the fiber. Because the mentioned three types of acoustic waves have different speeds of propagation, it is necessary, in principle, to create three areas of deformation, each designed for different intervals of wavelengths and for the reflection of different types of waves. In Ref. [28] the calculations of transmission probabilities of acoustic waves in chirp-deformed corrugation were executed in one-dimensional approximation. The purpose of this research is to check the validity of such approximation by comparing 1D and 3D models, and to analyze the influence of boundary design for thermal conductance in a nanofiber depending on temperature.
For a cylindrical fiber, there are four independent modes for which oj,min ¼ 0: a longitudinal mode, a torsional mode and two flexural modes. For all other modes oj,min ¼ ncj =D, where cj is the characteristic propagation speed of elastic wave j, D is the diameter of the rod, n ¼ 1; 2, . . .. At temperatures T o1 K and D o 100 nm only the mentioned four modes give the main contribution to O, so O ¼ 4Ouniv . It is known that at low temperatures an unstructured surface roughness of a wire can diminish quantum heat conductance due to the constructive interference of scattered waves [20]. The roughness structured in a proper way can sharply decrease the transmission probability pj ðoÞ. It is possible to design the surface of a wire in such a way that pj ðoÞ would be small for all phonons in a wide frequency band centered near o kB T=_, so that it would induce a substantial drop in total thermal conductance.
2. Physical background
3. 3D model of scattering in an elastic wire
In an elastic dielectric nanofiber heat energy is transferred by ! elastic waves. Such a wave with wavevector k and frequency ! oð k Þ can be represented as a flow of virtual particles known as phonons. Every phonon possesses energy E ¼ _o and crystal ! ! momentum p ¼ _ k . At temperatures T r1 K the phonon mean free path exceeds 100 mm (for defect-free fiber) and phonons propagate without collisions through the fiber. In other words, phonon flow is within the ballistic regime. Let reservoirs R1 and R2 with respective temperatures T1 and T 2 4T 1 be connected by a dielectric fiber G, operating in ballistic regime. Then the thermal energy passing through G from R2 to R1 in a unit time is [17] Z @Q 1 X ¼ _o½nðo,T 2 Þnðo,T 1 Þpj ðoÞ do @t 2p j
Let P ¼ fðx,y,zÞ A R3 : ðx,yÞ A Y, 1o z o þ 1g be an infinite cylinder whose cross-section Y R2 is a two-dimensional bounded and connected domain with smooth boundary. Assume that an elastic wire G of homogeneous and isotropic material coincides with P for 9z9 4L, and let
where nðo,T i Þ ¼ ðexpð_o=ðkB T i ÞÞ1Þ1
ð1Þ
is the Bose–Einstein distribution function; pj ðoÞ is the probability for a phonon of mode j to pass through G from R2 to R1; kB is the Boltzmann constant. Integration is over the whole spectrum of phonon frequencies. In the case 9T 2 T 1 9 5 T 2 ,T 1 we have Z @Q T 2 T 1 X @nðo,TÞ ¼ p j ð oÞ d o _o @t @T 2p j The value
O¼
1 @Q T 2 T 1 @t
is called the thermal conductance. Taking into account Eq. (1) we obtain the 1D Landauer formula for phonons [17]: Z _2 X oj,max o2 expðb_oÞ O¼ pj ðoÞ do ð2Þ 2 ½expðb_oÞ12 2pkB T j oj,min where b ¼ ðkB TÞ1 ; T T 1 T 2 and oj,min ðoj,max Þ is the minimal (maximal) phonon frequency for mode j. If for mode j the values pj ðoÞ ¼ 1, oj,min ¼ 0 and oj,max ¼ 1 the corresponding contribution in O is
Ouniv ¼
GðLÞ ¼ fðx,y,zÞ A G : 9z9 r Lg be a truncated part of it (see Fig. 1). For computational purposes parameter L4 0 is chosen large enough so that there are sufficiently long straight cylinders (with cross-section Y) contained in G(L). Propagation of 3D phonons in wire G can be investigated by the homogeneous ‘‘frequency domain’’ boundary value problem of elasticity theory: ! L u ðx,y,zÞ ¼ 0, ðx,y,zÞ A G ! n sð! u Þðx,y,zÞ ¼ 0, ðx,y,zÞ A @G ! where u ¼ ðu1 ,u2 ,u2 Þ stands for displacements and
ð3Þ
! ! ! ! L u ¼ mD u þ ðl þ mÞr div u þ ro2 u
ð4Þ
Lame´’s constants l and m satisfy m 4 0 and 3l þ 2m 4 0, and density r together with the frequency o is supposed to be ! constant. Vector n is the outward unit normal on @G. Finally, sðÞ is the stress tensor with Cartesian components
sð! u Þ ¼ ðlekk dij þ 2meij Þij where eij ¼ 12 ð@i uj þ @j ui Þ, i,j ¼ 1; 2,3. By setting the stress tensor to zero one considers a free boundary problem, i.e. no outer mechanical forces affect the boundary. A solution Y to problem (3) is called a continuous spectrum eigenfunction (CSE) corresponding to the frequency o 4 0, if Y satisfies the following conditions: Y is a bounded function in G : 9Yðx,y,zÞ9 rC, ðx,y,zÞ A G, C constant; and Y does not converge to zero at infinity: 9Yðx,y,zÞ9Q0 as 9ðx,y,zÞ9-1. The CSEs play the role of phonon wave functions. It is known that the total number NðoÞ of linearly independent CSEs below eigenfrequency o increases in steps as o increases, so that it is a piece-wise constant function. The step values oi of N are usually called ‘‘threshold’’ values. We assume here to consider only frequencies
p k2B T 6 _
For ballistic regime of thermal conductivity this universal value is independent of the properties of a fiber material and its geometry. If _oj,min bkB T, then the contribution of mode j in thermal conductivity is exponentially small and can be neglected.
Fig. 1. Illustration of a cylindrical fiber located on the z-axis with axially symmetric perturbation near origin. G(L) is used as a computational domain for scattering of elastic waves.
T. Puurtinen et al. / Physica E 44 (2012) 1189–1195
below the first threshold. It is known that for problem (4) the number of linearly independent CSEs is eight. In the space of CSEs in G corresponding to a frequency o, a basis fY 1 ð, oÞ, . . . ,Y 8 ð, oÞg can be chosen such that 8 X
Y j ð, oÞ ¼ wjþ ð, oÞ þ
g9z9 sjk ðoÞw Þ k ð, oÞ þ Oðe
ð5Þ
k¼1
where ‘’ represents a point ðx,y,zÞ in G, j ¼ 1, . . . ,8 and g is a sufficiently small positive number. The existence of such solutions was discussed e.g. in Ref. [30]. Functions wj7 :¼ wj expðaj ikzÞ ZR ðbj zÞ are called incoming ( þ) and outgoing ( ) waves. They can be explicitly described by choosing correct signs aj ¼ 71 and bj ¼ 7 1. However, depending on the cross-sectional shape Y, functions wj :¼ wj ðx,y,zÞ do not usually possess an analytical form. We can choose the numbering in such a way that wj7 with j ¼ 1, . . . ,4 ðj ¼ 5, . . . ,8Þ have a support on negative (positive) z-half space, respectively. If Y is a circle, the analytical expression exists and functions wj can be identified as flexural (with two polarizations), torsional and longitudinal waves. Here ZR is any smooth cutoff function with the properties ZR ðz ZRÞ ¼ 1, and ZR ðz r REÞ ¼ 0 for 0 o E o R. Asymptotic behavior of Yj is independent of the choice of ZR . Finally, the matrix SðoÞ ¼ ðsjk ðoÞÞ8j,k ¼ 1 is unitary and it is known as the scattering matrix. Let us examine the solutions Yj in Eq. (5) for j ¼ 1, . . . ,4. Solution Yj describes the phonon incoming from 1 as the wave wjþ , which scatters from the perturbed part of the wire G located near the origin into eight outgoing waves w k . With such a numbering convention the sum pj ðoÞ :¼
8 X
2
9sjk ðoÞ9
k¼5
is called the transmission coefficient. It can be interpreted as the transmission probability for the phonon Yj through the wire G P 2 from 1 to þ1. Finally, 1pj ðoÞ ¼ 4k ¼ 1 9sjk ðoÞ9 is the reflection coefficient.
Non-zero real functions Bx(z) and By(z) are the elliptical semiaxes. We also assume that the cross-section is a circle of constant radius outside of the corrugated part GðlÞ :¼ G \ f9z9 r lg (see Fig. 1), i.e. Bx ðzÞ ¼ By ðzÞ for 9z9 Z l. We can now use three simple 1D ordinary differential equation models to approximate various 3D waves. The equations can be derived e.g. from physical basis [29]. For longitudinal modes we use the equation: 1 @ @uðzÞ r AðzÞ ¼ o2 uðzÞ ð6Þ AðzÞ @z @z E and for torsional modes the equation: " 3 # Bx ðzÞB3y ðzÞ @uðzÞ @ r Bx ðzÞBy ðzÞðB2x ðzÞ þ B2y ðzÞÞo2 uðzÞ ¼ @z B2x ðzÞ þ B2y ðzÞ @z 4m
In various situations a 3D system can be approximately described by means of 1D model. Such a situation arises in our case for wave propagation in a fiber when the fiber length as well as the wavelength is much greater than the diameter of the fiber. In this case the equation describing a 3D model can be approximately replaced by equations corresponding to various 1D models. Such approaches allow getting description of the process which is adequate over limited range of parameters. The limits of the range will be clarified with a reasonable accuracy by numerical experiments. The low temperature thermal conductance of elastic dielectric fibers is mainly determined by waves which have wavelengths l about 0:121 mm at sub-Kelvin temperatures. Thus, for nanofibers with diameter 10–20 nm the passage is possible from 3D to 1D model. The transition from 3D to 1D is done separately for each mode. It is known (see e.g. Ref. [29]) that in a cylinder each of the elastic modes j with oj,min ¼ 0 is mainly one-dimensional in their nature, i.e. each mode has major displacement component in only one direction in some coordinate system. In cylindrical fiber, flexural modes have major displacement in the x, y-directions, longitudinal in the z-direction and torsional in the y-direction in cylindrical coordinate system. We henceforth assume that the cross-section of 3D fiber G is a circle or an ellipse, so that G can be defined by G ¼ fðx,y,zÞ A R3 : x2 =B2x ðzÞ þ y2 =B2y ðzÞ r1g
ð7Þ
To simulate flexural phonon modes with polarization in the x-direction we use the equation: 2 @2 4r 3 @ uðzÞ Bx ðzÞBy ðzÞo2 uðzÞ B ðzÞB ðzÞ ð8Þ ¼ x y E @z2 @z2 Here r is a density, E ¼ mð3l þ2mÞ=ðl þ mÞ is Young modulus and o is the angular frequency of the wave. In Eq. (6) AðzÞ ¼ pBx ðzÞBy ðzÞ is the area of cross-section of 3D domain at z. For y-polarized flexural waves indices x and y in the above equation are interchanged.
5. Computational method A finite part G(L), L4 l, of fiber G is taken as a computational domain. In order to solve the complex scattering amplitudes sjk of þ outgoing waves w k for incoming wave wj in Eq. (5) we simplify the calculation by further assuming that the perturbation on G(L) should have a cross-sectional symmetry in such a way that there is no mixing of normal mode types after scattering. In other words, ansatz Y j ¼ wj expðikzÞZR ðzÞ þ sjj wj expðikzÞZR ðzÞ þsj,j þ 4 wj expðikzÞZR ðzÞ þ Oðeg9z9 Þ
4. 1D models of scattering in an elastic wire
1191
in GðLÞ
ð9Þ
is valid in the domain G(L). Here wj :¼ wj ðx,y,zÞ and ZR is any smooth cut-off function with the properties ZR ðz Z RÞ ¼ 1 and ZR ðz r REÞ ¼ 0 as in Section 3. Number k4 0 is the wavenumber and it depends on the angular frequency o and the type of the phonon j. The exact calculation of frequency relation oj ðkÞ for circular cylinders was conducted in Ref. [29]. However, we decided to use 1D approximations of this relation also for 3D model, because at low frequencies such an approximation was found to provide a sufficient numerical accuracy. Inverse relations kj ðoÞ are sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffi rffiffiffiffi rffiffiffiffi r r o 4r ð10Þ , ktors ðoÞ ¼ o klong ðoÞ ¼ o , kflex ðoÞ ¼ E R E m For flexural waves the constant R :¼ Bx ð 7 LÞ ¼ By ð 7 LÞ is the radius of the cylindrical fiber at the ends. We use the following convention to number the phonon modes: j ¼ 1; 2 are the two flexural modes, j¼3 is the torsional mode and j¼4 is the longitudinal mode. We now explain the choice of boundary conditions on the cylindrical ends Gð 7 Þ :¼ fðx,y,zÞ : z ¼ 7 Lg \ GðLÞ of the domain. They will be used to numerically solve the amplitudes sjj and sj,j þ 4 needed in the calculation of transmission probabilities. Let us first take a look at boundary Gð þ Þ on positive z-axis where there is no incoming wave present: Y j wj expðikzÞ. We therefore require ð@z þ ikÞY j ¼ 0
on Gð þ Þ
1192
T. Puurtinen et al. / Physica E 44 (2012) 1189–1195
Now at the boundary on negative z-axis we must allow both incoming and outgoing waves of the same type to exist, so that Y j wj expðikzÞ wj expðikzÞ. We therefore put ð@z þ ikÞY j ¼ 2ikwj expðikzÞ
on G
ðÞ 2
Transmission and reflection probabilities pj :¼ 9sj,j þ 4 9 and 2 pj,ref :¼ 9sjj 9 can be calculated from the elastic field Yj which is a solution to frequency domain problem: LY j ¼ 0 in GðLÞ ! n sðY j Þ ¼ 0 on @GðLÞ\ðGð þ Þ [ GðÞ Þ on Gð þ Þ
ð@z þ ikÞY j ¼ 0
ð@z þ ikÞY j ¼ 2ikwj expðikzÞ
on GðÞ
ð11Þ
where L can be any of the operators defined by Eqs. (4), (6), (7) or (8) and k ¼ kj ðoÞ is the inverse 1D frequency equation of the corresponding wave type given by Eq. (10). Note that in 1D case the domain G(L) reduces to ½L,L and the free boundary condition in problem (11) is not needed. Transmission coefficient pj ðoÞ is calculated by evaluating the ! ! field value Y j ð z Þ on a point z ¼ ð½x,y,zÞ near the boundary Gð þ Þ . By taking squared absolute values from both sides of Eq. (9) and ! 2 by normalizing with 9wj ð z Þ9 : 2
pj ðoÞ ¼ 9sj,j þ 4 9
! 2 9Y j ð z Þ9 2 ! 9wj ð z Þ expðikzÞ9 2
2
Generally pj,ref ¼ 9sjj 9 ¼ 19sj,j þ 4 9 holds, but in order to observe possible numerical error we calculate the reflection coefficient ! separately. We therefore choose a point z near GðÞ and calculate ! from Y j ð z Þ wj ðx,y,zÞ expðikzÞ þ sjj wj ðx,y,zÞ expðikzÞ 2
pj,ref ¼ 9sjj 9
2 ! ! 9Y j ð z Þwj ð z Þ expðikzÞ9 2 ! 9wj ð z Þ expðikzÞ9
Even in the simplest fiber geometries functions wj can have quite complicated analytical forms (see e.g. Ref. [29]), so in order to further simplify the analysis we use approximations for them as well. We found the following forms sufficient (only small discrepancy with numerically solved eigenmodes) at low frequencies and small amplitudes. For longitudinal wave wlong ðx,y,zÞ ð0; 0, EÞ, for x-polarized flexural wave wflex ðx,y,zÞ ¼ ðE,0,iEkxÞ and for torsional wave wtors ðx,y,zÞ ¼ ððx cosðEÞy sinðEÞxÞ,ðx sinðEÞ þ y cosðEÞyÞ,0Þ where E 4 0 is small. We also remark that for different modes ! different evaluation points z ¼ ðx,y,zÞ must be chosen. In 3D model, longitudinal and flexural modes are evaluated on the z-axis while torsional modes are evaluated near the cylinder boundary where displacement is large. All three 1D models and 3D model were solved with finite element method (FEM). Implementation was done in Comsol v3.5 and Matlab 2008b using coefficient form PDE application mode for 1D models and structural mechanics application mode for 3D Lame´ system. For each elastic mode we put suitable boundary conditions on the cylindrical ends to transform problem (11) ! ! ! into a standard form used in Comsol: n cr u þq u ¼ g, where cr ¼ sðÞ is the stress tensor. To solve the arising FEM problem we used swept mesh generation in Comsol, which sweeps a 2D mesh from GðÞ to Gð þ Þ by constant steps. This method provided numerically reliable solutions with ease, as much finer mesh was needed if free mesh generation was used instead. Especially with flexural 3D modes the stability of solutions was heavily dependent on the mesh used. In the simulations second-order elements were used and the number of degrees of freedom was usually around 60,000 depending on the fiber geometry. For the linear equation we used direct linear solver Spooles v2.2.
6. Periodic structures of finite length It is known that infinite periodic structures can possess spectral gaps in frequency spectrum which prevent propagation of waves with such frequencies. It was numerically demonstrated in e.g. Ref. [28] that as the number of periods increases the transmission probability curve converges to a characteristic function of spectral gaps of similar infinite periodical structure. Therefore an effective phonon reflection over a thin frequency band can be established already with a small number (5–20) of periods. We used a simple harmonic periodical structure to test the suitability of 1D models for calculation of transmission probabilities. Assume that the radius of circular fiber along the z-axis is defined by n o a BðzÞ ¼ R 1 ½1 þ signðl9z9Þ½1 þ cosðLðzÞÞ ð12Þ 4 where R 40 is the radius of fiber, 0 r a o 1 is perturbation strength and l, 0 ol o L, is the half-length of perturbation. For a strictly periodical perturbation function LðzÞ ¼ pððz þ lÞn=l þ 1Þ is linear. Number n stands for the number of perturbations and 2l/n then becomes the length of one period in the perturbation. Total length of the fiber is 2L and total length of the perturbed part 2l (see Fig. 2). In Fig. 3 there are transmission probabilities of 1D and 3D torsional and longitudinal waves in a cylindrical Silicon fiber with a finite periodical deformation. It can be seen that 1D and 3D models give a very good agreement in the results with only small increasing error as the frequency o increases. This error can be explained as a numerical error arising from a fact that wavelength l approaches the element length in our FEM-system and it can be negated by increasing the mesh size. In Fig. 4 we show the transmission probability of flexural waves from 1D and 3D models in a Silicon fiber with finite harmonic perturbation on boundary. It is again seen that at low frequencies 1D model gives a fairly good approximation but the error slightly increases as o increases. Differences can be again explained as an error arising from diminishing of wavelength but also with the fact that 1D approximation of kflex ðoÞ is not as good for flexural waves as it is for torsional and longitudinal waves. Even that the 1D equations cannot model scattering into different mode types, such effects (if present) would induce noticeable changes in 3D transmission probabilities so that there would be significant discrepancy in 1D and 3D transmission probability curves. Therefore the above results predict that 1D models (6), (7) and (8) can provide sufficiently reliable transmission probability curves for the calculation of thermal conductance OðTÞ in high cross-sectional symmetry geometries with tens or hundreds of periods of perturbation.
Fig. 2. Illustration of a typical function B, which was used to define the crosssectional radius of cylindrical fiber G with radius R ¼ 5 nm and total length 2L ¼ 240 nm. There are eight periods of harmonic perturbation on the boundary with wavelength 20 nm and amplitude 0:251R.
1
1
0.8
0.8
Ω/Ωuniv, ptors, plong
ptors, plong
T. Puurtinen et al. / Physica E 44 (2012) 1189–1195
0.6
0.4
0.2
1193
0.6
0.4
0.2
0 0
100
200
300
400
500
600
700
0
800
0
ω [108 s−1]
2
4
6
8
10
T [K], X [hω/(2πkBT)]
Fig. 3. Transmission probability of longitudinal and torsional waves in cylindrical Silicon fiber with finite periodical perturbation on the boundary. Black lines correspond to torsional and grey lines to longitudinal waves. Dashed lines are calculated from 1D models and solid line from 3D model. Parameters for geometry: radius R ¼ 10 nm, number of periods n¼ 4, period length 2l=n ¼ 400 nm, perturbation strength a ¼ 0:251 and total fiber length 2:24 mm.
Fig. 5. Thermal conductance by torsional and longitudinal 1D waves in Silicon fiber with chirp. Solid lines correspond to thermal conductance normalized by the universal quantum and dashed lines correspond to scaled transmission probabilities ptors ðXÞ (black) and plong ðXÞ (dark grey) where X ¼ _o=ðkB TÞ and T ¼ 1 K. Light grey dash-dotted line is the scaled phonon distribution D(X). Parameters for geometry: radius R ¼ 5 nm, total fiber length 107 mm, number of periods n ¼80, first period length l1 ¼ 30 nm, period ratio C¼ 1.0664, chirp ratio C¼ 160 and perturbation strength a ¼ 0:251. Inset picture shows the periodicity cell.
1
1 0.8
Ω/Ωuniv, pflex
pflex
0.8 0.6
0.4
0.6
0.4
0.2
0.2 0 0
100
200
300
400 8 −1
500
600
700
0
ω [10 s ]
0
Fig. 4. Transmission probability of flexural waves in Silicon fiber with periodical perturbation. Dashed line corresponds to 1D model and solid line to 3D model. Parameters for geometry: radius R ¼ 10 nm, number of periods n¼ 4, period length 2l=n ¼ 100 nm, perturbation strength a ¼ 0:5 and total fiber length 960 nm.
7. Almost periodic structures To broaden the reflection band caused by periodic perturbation we change the periodicity slowly. This structure is called a chirp. The underlying idea is that chirp looks locally strictly periodical and each of these locally periodical parts contributes a thin reflection band at its resonant frequency. We now explain the mathematical formulation of such structure. Let the perturbation on the boundary consist of n periods of harmonic deformation with lengths lj , j ¼ 1, . . . ,n. Assume that lj þ 1 =lj ¼: d 4 1 is constant for all j and call it a period ratio. Now the radius of the fiber is defined by Eq. (12), but instead of linear LðzÞ we put
LðzÞ ¼ p 1 þ 2 logd 1 þ
ðzlÞðd1Þ l1
ð13Þ
2
4
6
8
10
T [K], X [hω/(2πkBT)] Fig. 6. Thermal conductance by flexural 1D waves in Silicon fiber with elliptical chirp. Solid line corresponds to thermal conductance normalized by the universal quantum and dashed line corresponds to scaled transmission probability pflex ðXÞ, where X ¼ _o=ðkB TÞ and T ¼ 1 K. Light grey dash-dotted line is the scaled phonon distribution D(X). Parameters for geometry: radius R ¼ 5 nm, total fiber length 4:3 mm, number of periods n¼50, first period length l1 ¼ 20 nm, period ratio d¼ 1.04, chirp ratio C¼ 6.8 and perturbation strength a ¼ 0:7. Inset picture shows the periodicity cell.
where l1 is the length of the first period (when looking from GðÞ ) and l¼
n1 1X j l1 d 2j¼0
is the haft-length of the chirp. We also mark ln =l1 ¼: C and call it the chirp ratio. In Fig. 5 we show the thermal conductance by torsional and longitudinal 1D waves (solid lines) through a long chirp with 80 periods of perturbation. According to simulations thermal conductance by these waves can decrease down to 10% of the
1194
T. Puurtinen et al. / Physica E 44 (2012) 1189–1195
The appearance of multiple gaps can be utilized in reflecting torsional and longitudinal waves as well. When the boundary of the cylindrical perturbation is defined as a superposition of harmonic functions of different period length (as in Fourier expansion), each of the present functions contributes to reflection at their respective frequencies. If the boundary is chosen to locally obey triangle wave shape (the chirp is built of truncated cones) we can see from Fig. 7 that thermal conductance is reduced on wider temperature range when compared to similar harmonic perturbation with matching fundamental component of Fourier expansion of each period.
1
Ω/Ωuniv
0.8
0.6
0.4
8. Conclusions
0.2
0 0
0.5
1
1.5
2
2.5
3
3.5
4
T [K] Fig. 7. Thermal conductance by torsional 1D waves in Silicon fibers with conical (triangular) axially symmetric chirp and similar harmonic chirp. Black line corresponds to thermal conductance of harmonic chirp and grey line to conical chirp. Parameters for geometry: radius R ¼ 5 nm, total fiber length 28:5 mm, number of periods n ¼30, first period length l1 ¼ 100 nm, period ratio d ¼ 1.0826, chirp ratio C¼ 10 and perturbation strength a ¼ 0:251.
universal quantum Ouniv on temperature range 0.1–1 K. In the figure there are also scaled transmission probability curves pj(X) for torsional and longitudinal waves, where X ¼ _o=ðkB TÞ and T ¼ 1 K is constant and, for comparison, the phonon distribution DðXÞ ¼
X 2 expðXÞ expðXÞ1
from the 1D Landauer formula (2). It can be seen that transmission probability curves are close to zero on 0:1 o X o 5 where there is most of the weight of the distribution D(X). We also found that harmonic perturbation with strength a ¼ 0:251 is not very effective in reflecting flexural waves. To reach 80–90% drop in thermal conductance for temperature range 0.1–1 K we had to use 350 periods in chirp with d ¼1.005, l1 ¼ 15 nm and R ¼ 5 nm. To reflect flexural waves efficiently one should be able to use stronger perturbation without deteriorating the tensile strength of the fiber. One possible solution to this is to use elliptical deformation where cross-sectional area is kept constant. Let cðzÞ ¼ 12 ð1 þsignðsin zÞÞ and let SðzÞ ¼ cðzÞð1 þ a sin zÞ þ cðz þ pÞ
1 1 þ a sinðz þ pÞ
where a is the strength of the perturbation. Periodic function SðÞ has a property SðzÞ ¼ Sðz þ pÞ1 . We now define Bx ðzÞ ¼
R 1signðl9z9Þ þð1 þ signðl9z9ÞÞSð2pLðz þ lÞÞ 2
where function LðÞ is given by Eq. (13). Finally by setting By ðzÞ ¼ RðBx ðzÞ=RÞ1 the cross-sectional area of an ellipse with semi-axes Bx and By is A :¼ pBx ðzÞBy ðzÞ ¼ pR2 . In Fig. 6 we show the thermal conductance by flexural 1D waves (solid line) through a circular fiber with elliptically deformed chirp. It can be seen that normalized thermal conductance goes down to (0.05–0.1)Ouniv on temperature range 0.1–1 K. Scaled transmission probability curve does not approach unity when X o 10, because elliptical (as well as circular) deformation opens multiple gaps in frequency spectrum of flexural waves, which are located at higher frequencies. The effect of higher gaps on hindering the wave propagation is reduced due to the fact that gap width is lower.
We investigated the possibility of decreasing the thermal conductance of Silicon nanofibers by designing the shape of the boundary to reflect propagating acoustic phonons. In cylindrical nanofibers of defect-free dielectric material only the four lowest energy acoustic phonons (longitudinal, torsional and two flexural waves) mainly contribute to thermal flow in temperatures 0.1–1 K. Propagation and scattering of these phonons can be studied within elasticity theory. We numerically simulated with 1D and 3D finite element models propagation of elastic waves in cylindrical nanofibers with periodical and almost periodical deformation on the boundary by calculating the transmission probability curves. We compared the curves obtained from 1D and 3D models for various simple corrugations and suggested designs where transition from 3D to 1D is possible. Finally, thermal conductance by each of the four low frequency phonons was calculated with 1D Landauer formula for long nanofibers. It was found that by assuming high-enough cross-sectional symmetry (to neglect phonon scattering into separate modes) and weak enough perturbation 1D approximation was reliable and the differences in transmission probability curves resulted into negligible discrepancies in thermal conductance curves. In order to effectively decrease thermal conductance of the fiber without deteriorating the structural strength we designed a long locally periodical corrugation (chirp) for each phonon. For torsional and longitudinal modes we used circular cross-section where the radius is varied locally periodically along the corrugation. For flexural modes such construction was found inefficient. To reflect flexural modes we designed an elliptical corrugation where the cross-sectional area is kept constant. It was found that with these designs the thermal conductance could be reduced down to 10% of the universal quantum on wide temperature range 0.1–1 K individually for all phonons. We conclude that a similar reduction in total thermal conductance could be obtained by combining such corrugations in series. However, with harmonic and elliptical fiber designs in series the total length of the corrugation can approach the reference phonon mean free path 100 mm. Even greater decreases in thermal conductance could be obtained using longer corrugations. The total fiber length can be reduced not only by designing the corrugation for narrower temperature range, but also by choosing a different period shape for chirp, which could simplify the possible experimental measurements in future. To this end, we finally compared thermal conductances by torsional waves in similar axially symmetric harmonic and conical periodical chirp structures. It was found that the temperature zone with reduced thermal conductance could be slightly expanded using a conical chirp design. Advantage of such a structure is that it may be easier to fabricate by etching processes. References [1] S.M. Lee, D.G. Cahill, R. Venkatasubramanian, Applied Physics Letters 70 (1997) 2957. [2] G. Chen, Physical Review B 57 (1998) 14958.
T. Puurtinen et al. / Physica E 44 (2012) 1189–1195
[3] W.E. Bies, R.J. Radtke, H. Ehrenreich, Journal of Applied Physics 88 (2000) 1498. [4] M.V. Simkin, G.D. Mahan, Physical Review Letters 84 (2000) 927. [5] X.Y. Yu, G. Chen, A. Verma, J.S. Smith, Applied Physics Letters 67 (1995) 3554. [6] A. Balandin, K.L. Wang, Physical Review B 58 (1998) 1544. [7] J. Philip, et al., Journal of Applied Physics 93 (2003) 2164. [8] T.S. Tighe, J.M. Worlock, M.L. Roukes, Applied Physics Letters 70 (1997) 2687. [9] A. Balandin, Physics of Low-Dimensional Structures 1–2 (2000) 1. [10] G. Chen, International Journal of Thermal Sciences 39 (2000) 471. [11] G. Chen, T. Zeng, Microstructure Thermophysical Engineering 5 (2001) 71. [12] B.A. Glavin, Physical Review Letters 86 (2001) 4318. [13] S. Volz, D. Lemonnier, J.B. Saulnier, Microstructure Thermophysical Engineering 5 (2001) 191. [14] J. Zou, A. Balandin, Journal of Applied Physics 89 (2001) 2932. [15] W.-Q. Huang, et al., Physics Letters A 336 (2005) 245. [16] P. Kim, L. Shi, A. Majumdar, P.L. McEuen, Physical Review Letters 87 (2001) 215502.
[17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]
1195
L.G.C. Rego, G. Kirczenow, Physical Review Letters 81 (1998) 232. M.P. Blencowe, Physical Review B 59 (1999) 4992. K. Schwab, E.A. Henriksen, J.M. Worlock, M.L. Roukes, Nature 404 (2000) 974. D.H. Santamore, M.C. Cross, Physical Review Letters 87 (2001) 115502. D.H. Santamore, M.C. Cross, Physical Review B 63 (2001) 184306. X. Lu, J.H. Chu, W.Z. Shen, Journal of Applied Physics 93 (2003) 1219. W.-X. Li, et al., Applied Physics Letters 85 (2004) 822. W.-X. Li, et al., Journal of Physics: Condensed Matter 16 (2004) 5049. W.-Q. Huang, et al., Journal of Applied Physics 98 (2005) 93524. K.-Q. Chen, et al., Physical Review B 72 (2005) 45422. L.-Y. Nie, et al., Physics Letters A 359 (2006) 234. T. Puurtinen, Numerical Simulation of Low Temperature Thermal Conduc¨ tance of Corrugated Nanofibers, Jyvaskyl a¨ Studies in Computing, vol. 124, ¨ ¨ 2010. University of Jyvaskyl a, [29] K.F. Graff, Wave Motion in Elastic Solids, Dover, New York, 1991. [30] S.A. Nazarov, B.A. Plamenevsky, Elliptic Problems in Domains with Piece-wise Smooth Boundaries, Walter de Gruyter, Berlin, 1994.