Finite-size effect of the thermal conductivity in one dimensional chain

Finite-size effect of the thermal conductivity in one dimensional chain is theoretically investigated with a phenomenological model. Boundary temperature dropping and the spectrum-dependent linear temperature profile within the chain are given self-consistently at high temperatures. A general formula of thermal conductivity is presented and a good agreement with the molecular dynamics simulations is obtained. It is found that while d 1 κ / d 1 L decreases with increasing of 1/L monotonically, the well known linear dependence of thermal resistance 1/κ on the reciprocal length 1/L of systems is, however, valid only in two limit cases of 1 / L → ∞ and/or 1 / L → 0 .


Introduction
Thermal conductivity of systems is one of the most basic and important property in solid physics [1]. By designing and constructing modern spintronic and semiconductor devices into micro-and/or nano-scales, a finite-size effect known as the Casimir limit [2] stem from the sizes and boundaries emergences. In addition to traditional interests, bearing the potential applications in thermoelectrics [3][4][5], the finite-size effect is one of the most hot topics in the study of thermal conductivity . In the early study of molecular dynamics [30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45] (which has been used widely for decades), a linear dependence [41,42,45] of the thermal resistance 1/κ on the reciprocal length of the system 1/L was proposed and extrapolated to infinite systems [45,46]. Such relation is supported by some experimental observations [47,48] and is used comprehensively. However, Sellan et al [17] reported that the relation is hold only if the size of the sample is larger than a critical length, and they expanded 1/κ with Taylor's series of 1/L. Nevertheless, recent studies show [28,29] that neither the linear dependence nor the Taylor's expansion (up to the second power) can be extrapolated to all the scenarios. The problem is still not fully appreciated and more investigations are needed.
Qualitatively, it is easy to see that the rapid nonlinear temperature drops [22,29,30,46] in boundary layers are chief reasons for the decrease of thermal conductivities with decreasing system size. As a result, a smaller temperature gradient in the middle regions, a smaller heat current, and a smaller thermal conductivity than ideal cases emergence. The smaller the sample size, the more significant of the effect is observed. What can be determined is that the effect relate to the comparison of the mean free path l and the sample size L, but the quantitative knowledge of how the temperature drops change with the relation (l and L) is still not very clear. Another essential issue is that the mean free paths of phonons are spectrum dependent. In the traditional models (such as Normal Model of Vibrations of a lattice, relaxation time approximation and Callway Model) [49][50][51][52][53][54][55] and current molecular dynamics studies of thermal conductivity, a uniform local temperature profile (where a local thermal equilibrium is assumed) either predefined or derived is comprehensively used. However, we must note that with the decrease of the sample size, some phonons will not follow but separate from the local thermal equilibrium, for their relaxation times is much longer than their transit time along the thermal non-equilibrium chain, thus a differentiation of temperature profiles for phonons with different frequencies can be deduced. For the non-equilibrium and mixed property, models for limiting cases (such as diffusion model and ballistic model) cannot be used and a spectral theory is needed in the present general study.
Theoretically, the passing rate of spectral phonons can always be counted in the view of quasi-particles, by examining and comparing the passing rates of spectral phonons in non-equilibrium systems with the ones of equivalent systems with equilibrium temperatures, the effective local temperatures of the spectral phonons can be obtained [56]. Aimed at the two issues mentioned above, in this paper, a phenomenological model that is not based on local thermal equilibrium has been developed, and the size effects of the temperature profiles and the thermal conductivities in one dimensional lattice have been studied. The contents are arranged as follows: in the second section, the phenomenological model is suggested; in section 3, the effective local temperature of spectral phonons has been studied, temperature function has been established, the self-consistent numerical calculations have been carried out at high temperatures, and spectrum dependent linear temperature profiles with temperature drops in boundaries has been founded. In section 4, the thermal conductivity has been studied. Based on the effective local temperatures, general formula of thermal conductivities has been established and discussed in general and limiting cases, discussions and fittings of our result with molecular dynamics (EMD and NEMD) have been carried, a good agreement is identified. Finally, a conclusion is summarized.

Phenomenology model of effective local temperatures
Consider the one-dimensional chain shown in figure 1. Where T 1 and T 2 are the temperatures of the heat reservoirs on the left and right end of the sample respectively, with T 1 >T 2 . Please note that the heat reservoirs here are ideal, phonons in the reservoirs are set identical to the ones in equivalent equilibrium systems with equal temperatures. Based on the effective local temperatures which are spectrum dependent (a similar spectrum theory has been studied and used in the magnon-phonon systems [57,58]), by surveying the propagation of spectrum phonons, the spectrum dependent temperature profiles and heat transfers can be studied. With these considerations, we propose the following phenomenological model: (1) The effective local temperature T(x, ω) in the system without sources and/or sinks, is defined by the passing rate of ω phonons at position x, )and N ps (x, ω) are the forward, reverse and total passing rate of ω phonons respectively at position x, and N , ps  w ( ) is the passing rate of the ω phonons with equilibrium temperature  in the equivalent one-dimensional system. That is to say, if the passing rates (ω phonons in the sample somewhere, and in equivalent equilibrium system) are equal, the temperatures should be equal. And vice versa, if The observing temperature here is defined by the energy passing rate of full spectrum phonons, (2) Annihilation properties are depicted by the lifetimes and the mean free paths of the phonons. At time t, the existing probability of a ω phonon that exists at zero time satisfies [56]  where t w ( )is the lifetime of the ω phonons. The mean free path relate to the lifetime as follows [17,29,38] l where v(ω) is the propagation velocity of ω phonons. An exact knowledge of how the mean free path depend on ω is extremely hard [14,19]. We will see in section 3 that mean free paths of spectrum phonons influence the effective local temperature profoundly, correspondingly, mean free paths of spectrum phonons can be determined by their temperature profiles, which may be used to verify the validity of the assumptions proposed here.
(3) The probabilities of a creating phonon along the forward direction and the reverse direction are equal, and the creation rate of ω phonons is determined as follows, where g(x, ω) is the creation rate of ω phonons at x position, and g ,  w ( ) is the creation rate of ω phonons in the equivalent system with equilibrium temperature T x, ).
(4) The heat current of ω phonons is determined by the difference of the forward and reverse passing rates, can be written as follows ps ps where I(x, ω) is the heat current of ω phonons at x position, I + (x, ω)/I − (x, ω) is the forward/reverse heat current due to the ω phonons at x position (forward/reverse propagation phonons corresponding to positive/negative heat current), and ÿ is the Plank constant.

Effective local temperature profiles of the thermal non-equilibrium chain
In the equivalent one dimensional chain with equilibrium temperature , the passing rates of the phonons are equal everywhere, and satisfy: ) is the total number of ω phonons, and f k , ] is the Bose-Einstein distribution function, factor 2 comes from the double degeneracy of ω phonons in the one dimensional chain. The total creation rate should be equal to the negative total annihilation rate as the system is in equilibrium, according to equation (3) ) is the total creation rate of ω phonons. The creation rate on unite length is Considering the non-equilibrium one dimensional chain in figure 2, the temperature profile of w phonons can be investigated as follows. The spectral phonons passing through x position come from everywhere of the chain, and the passing rate of ω phonons can be written as follows , , , d , 9 )is the probability that a ω phonon created at position x¢ and passing through position x. According to the model assumption in item (3), using equations (3) and (4), )can be expressed as follows Boundaries are the concentrate regions of phonon creating and annihilating, so boundaries have an important influence on equation (9). Assuming that the reservoirs and the one dimensional lattice are full coupled, the created phonons at the left/right end of the chain per unit time should be equal to the half passing rate of the phonons in left/right reservoir, with N ps (T 1 , ω)/2 at −L/2 and N ps (T 2 , ω)/2 at L/2. Thus equation (9) can be expressed as where x 1 and x 2 are the coordinates of the left and right end of the one dimensional chain respectively, with x 2 −x 1 =L and x 1 =−L/2, x 2 =L/2. The creation rate of ω phonons on unit length at position x¢ should be equal to the creation rate in the equivalent system with equilibrium temperature T x , ), combing equations (7), (8), (10), and (11), the passing rate can be expressed as

At high temperatures, Bose-Einstein function can be approximated as
, thus the distribution equation of the effective local temperatures is obtained, The formula shows that the effective local temperatures of spectral phonons are related to their mean free paths: The temperature profiles of spectrum phonons with different mean free paths are different and a spectrum dependent property can be expected. In order to investigate how the mean free paths influence on the temperatures, by equation (13), temperatures of spectrum phonons with different mean free paths (L/100, L/3, L/2, L, 2L, 3L, 1000L) have been calculated self-consistently with numerical method, where the temperatures of the reservoirs are set as 300 K and 270 K respectively and the results are shown in figure 3. Two points can be identified easily: (1) effective local temperatures varies linearly along the chain. (2) Temperature drops are observed in boundaries, and the drop increase with l(ω)/L; while l(ω)?L, such as l(ω)/L=1000, the drop is near (T 1 −T 2 )/2=15 K; while l (ω)=L, the drop tends to vanish. Using the linear property, the relationship between the temperature drop and the mean free path of the spectral phonons can be obtained. According to appendix A, the temperature gradient of ω phonons is and the drop is In general cases, the effective local temperatures is In the region of l(ω)/L = 1, the drop vanishes and the finite-size effect can be ignored, we have which is the familiar result. As the chain is much shorter than the mean free path of the phonons, L = l , equation (16) can be written as That's to say, while l(ω)/L ? 1 , effective local temperatures of the phonons are all the same, equal to (T 1 +T 2 )/2, which is the ballistic case in literatures. Equation (16) stems from the model assumptions proposed in this work, the validity is to be verified by experiments and theoretical calculations. In principle, if the motion of local molecules can be obtained by experimental measurements and/or theoretical calculations (such as MD Simulations) exactly, by spectral analysis, effective local temperature of spectrum phonons can be obtained. Thus equation (16) can be verified, meanwhile, mean free paths of spectrum phonons can be derived by equation (16) with known temperatures. The discussions and the numerical results show that when the size of the sample is comparable to the mean free paths of the phonons, finite-size effect must be considered and the temperature distributions are spectrum dependent: phonons with longer mean free paths possesses smaller temperature gradient and larger temperature drop in boundaries. The emergence of the drop is consistent with the assumptions we have adopted: at the high temperature boundary, ω phonons come from the high temperature reservoir and low temperature zones, the passing rate is less than the passing rate in high temperature heat reservoir (a similar analysis at low temperature boundary holds also). Thus, the effective local temperature, by our definition, drops appear. Please note that the temperature profiles obtained here are only for the phonons whose relaxation times are much longer than their transit time so that the temperature relaxation processes can be neglected. However, according to the derivations and discussions in B, the influence of phonon-phonon relaxation on temperature distribution is estimated to be small, thus equation (16) can be used for full spectrum phonons.

Thermal conductivity in the one dimensional chain
According to equation (6), the heat current of spectral phonons can be expressed as x Substituting equations (8), (10) and (20) to (19), the heat current can be written as follows Integrating and simplifying the function, the heat current of ω phonons can be obtained (please see appendix C) where c(ω)=C(ω)/L is the heat capacity of the ω phonons per unit length. The total heat current is the summation of the heat currents due to different frequencies where T T T 2 1 D =is the temperature difference of the heat reservoirs.
Thermal conductivity is defined by the Fourier's law [22,29], First the formula will be studied by investigating the two limit cases of L ? l(ω) and L = l(ω). While the length of the chain is much longer compared with all the mean free paths of the phonons with l(ω)/L = 1, equation (24) can be approximated as, , and i is the energy index. The following approximation is obtained immediately, where α 1 =1/κ 0 , and β 1 =2λ 1 /κ 0 . In this case 1/κ dependent on 1/L linearly and the slope is β 1 , which is consistent with MD simulations [29] and experimental results [23]. While in the case of L = l(ω) , the thermal conductivity can be approximated as In the condition of L l w  ( ), by equation (24), a linear dependence of 1/κ on 1/L with slope β 2 and intercept α 2 is also obtained.
The calculation shows that 1/κ dependent on 1/L linearly under the condition of L l w  ( )and/or L l w  ( ), the relation has been identified and extrapolated to general cases [28,29,[39][40][41][42][43][44][45][46]. However, in our study, the slops (as well as the intercepts) of the two limit cases are not equal, with ). Reasoning that the inequality of the slopes (β 1 and β 2 ), the dependence (1/κ versus 1/L) must involve nonlinear terms in general cases and which has been founded by recent studies [17,22,28,29]. To general cases, the derivative d d To solve the equation, l=l(ω) must be known. However, we do not plan to deal with the sophisticated question in this study, but to discuss the general properties of equation (27). Firstly, it shows that d d , the thermal conductivity κ increases with L, which is consistent with the experimental and theoretical results. Secondly, as has been shown before, in the limiting cases of L ? l (ω) and/or L = l(ω), a linear dependence between 1/κ and 1/L is valid. Thirdly, d d 0 must be true. Thus a bounded property can also be obtained, which is consistent with the report that fittings between the Taylor's expansion and the MD (NEMD) simulations are not very well [29]. If Taylor's series with terms greater than the first power, like where α, β, γ and ν are constant coefficients, δ2 is the highest power index. As 1/L goes to infinity (L 0  ), κ≈L δ /ν. According to the definition of thermal conductivity in equation (23), the heat current is , which is impossible and is conflict with the bounded property proposed here also. Thus the expansion of Taylor's series such as formula above can not be extrapolated to all the regions. To illustrate the properties intuitively, sketches have been drawn to depict the dependence (1/κ on 1/L and κ on L) qualitatively in figure 4.
To verify the results, in quasi-one dimensional SiNW [29], equation (24) and the MD (EMD and NEMD) simulations were fitted. By considering the reciprocal processes, and neglecting the influence of defects and higher order processes, we follow the literatures [49,51] (where a quadratic attenuation between τ(ω) and ω is proposed), with l l B 1 where l T max and l L max are the maximal mean free path of transverse and longitude phonons respectively, and 1/B T and 1/B L are the decay factors of the quadratic attenuation between l and ω. By Using the reported effective velocity v g =7.5 km s −1 [29] in the Q1D case, and choosing l , nm. Figure 5 shows a good agreement between MD (NEMD) simulations and our result, and l l , , are in their reasonable regions [29]. However, to the validity of parameter values of l T max , l L max , B L and B T , further studies are needed for the dependence between l and ω used here are not confirmed.
At high temperatures, equation (2) can be simplified as where  is the number of energy levels. The temperature drop in boundaries can be expressed as With the parameters used in figure 5, the observing local temperatures in the Q1D SiNW have been calculated, the results are shown in figure 6. The calculation may be used as a criterion to justify ours result. For this purpose, we need to measure the effective temperature profiles (and the drops in boundaries) in the chains (with different sizes), and compare it with our predictions in figure 6, a verdict (the validity of our theory) can thus be obtained if the quadratic attenuation (l l w = ( )) mentioned before is true.

Conclusion
To solve the challenging non-equilibrium problem, we proposed a spectrum dependent model, in which the temperature is defined by the passing rate of the spectrum phonons. With the model, finite-size effects of effective local temperatures and thermal conductivity have been studied. By taking into account of the difference of the mean free paths, spectrum dependent temperature profiles are founded. Two important properties: boundary temperature dropping and spectrum dependence are presented. As a result, thermal conductivities of spectral phonons are different as well, and the size effect appears naturally. The general formula of thermal conductivity equation (24) gives a pretty good agreement with the MD simulations and experimental measurements. This agreement suggests that the spectrum dependent mechanism plays a key role in the thermal non-equilibrium transport. The proposed mechanism can be used in seeking good (and/or poor) thermal conductors at macro/ nano scales, and has potential applications in thermoelectric devices.  Considering the linear distribution of effective local temperatures on the one-dimensional lattice, equation (13) can be expressed as Integrating the formulas, we get Sum over the two terms in equation (30) T L  T  L  T  T  T T  2  ,  2  , , , To estimate the influence of phonon relaxation, we list the estimations as follows: (1) In the limiting cases of L l L l and   for all the phonons, the conclusion of linear temperature profile is valid obviously, because the temperature profiles are all the same for all the phonons, thus the influence of relaxation (on temperature distribution) can be neglected.
(2) Usually, the relaxation time is much longer than the lifetime. At high temperatures studied in this paper, the Bose-Einstein function can be approximated as The relaxation time can be estimated as follows: for the phonons are bosons, the creation and annihilation rate of ω phonons must be proportional to f (ω)+1 and f (ω) respectively, can be expressed as P T f T , , and P(T, ω)f (T, ω) in a system with equilibrium temperature T, satisfies stand at high temperatures, τ(ω) is the lifetime of ω phonons. Supposing that ω phonons are stimulated into temperature T ω =T+ΔT for some reasons. The relaxation function is where t ¢ w is the relaxation time for ω phonons, and the following formula stand also where C≈k B is the thermal capacity of ω phonons, with the relaxation time is much larger than the lifetime of ω phonons.
To phonons whose mean free path is near/larger the size of the sample (l∼L or l ? L), the creation/ annihilation process is less than several times, the influence of the relaxation can be neglected, and the linear temperature profiles are hold for these phonons.
To phonons whose mean free path is much shorter than the size of the sample (l = L), the creation/ annihilation process are very frequent along the chain. According to figure 3 in this paper, we can see that the temperature profiles very close to the line connecting (−L/2, T 1 ) and (L/2, T 2 ). Considering that the creation/ annihilation processes are very rare for large mean free path phonons(l∼L or l ? L), we can deduce that the creation/annihilation processes of short mean free path phonons are mainly involved themselves (l = L). Thus a local thermal equilibrium between the phonons (l = L) can be distinguished. For the temperature profiles of these phonons (l = L) are very close to the line connecting (−L/2, T 1 ) and (L/2, T 2 ), the local thermal equilibrium temperature profile is very close to the line (connecting (−L/2, T 1 ) and (L/2, T 2 )). Thus a linear profile for these phonons is also holds. For the temperature profiles of these phonons only slightly deviate their un-relax lines, our results can be promoted to the phonons that their mean free paths are much shorter than the sample size, without considering the relaxation between the phonons.
In conclusion, in the condition of high temperatures, the linear temperature profiles for most phonons are hold and the conclusion can be used to study the thermal conductivity of small size samples.
the addition of the two terms of the integration and the heat flows injected by the reservoirs can be calculated as the heat capacity formula have been used here, 1 2 comes from the double-degenerate of the ω phonons in the one-dimensional lattice, satisfying For the other part of equation (21), we have We get the formula of the heat flows in the one-dimensional lattice  where  is the number of the energy levels, the terms in the square bracket can be deduced as which is less than zero (the only exception is ∀i, j satisfies l j =l i , the second derivative of 1/κ to 1/L equal to zero).