Kapitza thermal resistance in linear and nonlinear chain models: isotopic defect
11 Kapitza thermal resistance in linear and nonlinear chain models: isotopic defect.
O.V.Gendelman * , Jithu Paul Faculty of Mechanical Engineering, Technion – Israel Institute of Technology, Haifa, 3200003, Israel * - contacting author, [email protected] Abstract
Kapitza resistance in the chain models with internal defects is considered. For the case of the linear chain, the exact analytic solution for the boundary resistance is derived for arbitrary linear time-independent conservative inclusion or defect. A simple case of isolated isotopic defects is explored in more detail. Contrary to the bulk conductivity, in the linear chain, the Kapitza resistance is finite and well-defined. At the same time, the resistance strongly depends on the parameters of the thermostats, and therefore cannot be considered as a local property of the defect. Asymptotic scaling behavior of the heat flux in the case of very heavy defect is explored and compared to the nonlinear counterparts; similarities in the scaling behavior are revealed. For the lightweight isotopic defect in the linear chain, one encounters a typical dip of the temperature profile, related to weak excitation of the localized mode in the attenuation zone. If the nonlinear interactions are included, this dip can still appear at a relatively short time scale, with subsequent elimination due to the nonlinear interactions. This observation implies that even in the nonlinear chains, the linear dynamics can predict the main features of the short-time thermal profile if the temperature is low enough. Introduction
The history of thermal boundary resistance starts from the famous work of Kapitza on temperature discontinuity at the interface between a metal and liquid Helium II[1]. Multiple and versatile experimental results on the boundary thermal resistance are available in the literature[2–5]. The first attempt to offer a theoretical explanation based on the phonon transmission coefficient, particularly at low-temperature experiments, was made by Khalantikov in 1952 [6] and by Little in 1959 [7], and it is generally known as the acoustic mismatch model (AMM). The well-known limitations of AMM are that, first; it fails to predict the Kapitza resistance when the contacting structures are similar. Then, it considers only long-wavelength phonons since the theory has been developed for low-temperature experiments. The first drawback was somewhat corrected by including a multicomponent phonon distribution functions on each side of the interface 8]. Additional developments included account of the phonon scattering at and near the interface[9]. Several other works were based on transport theory approach, but complexity of the calculations limited its practical use[10–12]. Sluckin developed a simple model for the chain of billiard particles and the results were surprisingly matching with experiments[13]. There were attempts to extend Khalatnikov’s theory to harmonic systems [14,15]; in[15] the results were compared with molecular dynamics simulations. In 1987, Swartz [16] proposed diffuse mismatch model (DMM) to explain the higher temperature behavior of thermal boundary resistance where the AMM becomes ineffective. This theory considers other extremes compared to the AMM, that is, all the phonons are considered as diffusely scattered at the interface, and any acoustic correlation at the interface is destroyed. Harmonic crystals are the simplest model to study thermal properties of dielectrics. It is well-known that the dynamics of homogeneous harmonic chain in modal coordinates reduces to an assembly of freely moving phonons, therefore the thermal conductivity in this model is anomalous. The modern era studies of the heat transport in harmonic chain started from the famous work of Rieder-Lebowitz-Lieb (RLL) [17] that presented the exact solution for one-dimensional homogeneous chain coupled to Langevin thermostats. It was demonstrated that the heat flux is proportional to the temperature difference between the boundaries. Besides, the exponential boundary layer near the thermostats has been observed. Then the studies progressed into more complicated problems that involved various scattering mechanisms by inserting various impurities like defect mass or defective coupling, or the phonon-phonon scattering caused by the nonlinearity [18]. In[19], the problem of infinite disordered harmonic chain has been analyzed. It was demonstrated that the heat flux is proportional to the boundary temperature difference, provided that the spectral measure exhibits an absolutely continuous part. For the special case of alternating-mass-chain, the temperature profile oscillates through the chain as follows from the exact solution in [20]. Somewhat surprisingly, such temperature oscillations also revealed themselves in the alternate-mass FPU model [21]. Recent paper [22], explored numerically the Kapitza resistance in a variety of benchmark models with isolated isotopic defect. Current work revisits this problem, and presents exact analytic solution for the Kapitza resistance in the linear chain with isolated defect (interface) in the limit of infinite chain length. This exact solution allows comprehensive study of the boundary resistance in the linear chain, including the near-field around the defect. Besides, the heat flux exhibits interesting asymptotic behavior in the limit of large defect mass. This behavior is compared to numeric data concerning similar asymptotic limit for β-FPU and Lennard-Jones (LJ) chains [23, 24]. If the mass of the inclusion is less than the mass of other chain particles, the thermal profile exhibits an interesting "dip", indicating a "cold" point in the chain. It is a well-known artifact of the linear model; similar behavior is observed also at larger scales [23,24]. It is demonstrated below that account of the nonlinearity removes this anomaly. More exactly, when the temperature gradient is imposed, the dip is initially formed also in the nonlinear chains. Then the dip is gradually destroyed, ending up with monotonous thermal profile. Possible role of mobile breathers [29, 30] in destruction of the dip is discussed. The general outline of this paper is as follows. Section 2 presents the exact solution for the Kapitza resistance in the linear chain containing linear time-independent (LTI) inclusion. In Section 3, the expressions for special case-of the isotopic defect are derived and asymptotic limit for the heat flux in the case of very heavy defect is compared to similar phenomena in
FPU and Lennard-Jones models. Section 4 addresses in detail the thermal "dip" at light defect in the linear model and its connection to localized mode in the chain spectrum. The treatment is extended to the FPU chain, and disappearance of the dip is explored numerically. Then, the thermal near-field of the defect in linear model is explored and exponential convergence to the average temperature is established. The exponent is related to the defect mass. Section 5 is devoted to discussion and concluding remarks. Kapitza resistance in linear chain with inclusion: general treatment.
Let us consider a heat transport in a linear chain with linear time-independent (LTI) inclusion, with thermostats attached to the terminal particles at the right and at the left. General sketch of the model system is presented in Fig. 1. FIG. 1. Sketch of the model system All masses and stiffnesses within the chain fragments beyond the inclusion are set to unity without affecting the generality. At this stage, the structure of the inclusion, or defect, in the chain is left generic. The only assumption is that the motion of all particles that belong to the inclusion is described by linear time-independent equations. Beyond the inclusion and apart from the terminal particles, the chain is described by common equations of motion: n n n n u u u u (1) Traveling waves in this chain obey a well-known dispersion relation: q q (2) Here ( ) q denotes the frequency of the traveling wave, q is the wave vector. The system is excited through thermostats attached to the terminal particles. For this sake, we assume that both these particles are subject to linear viscous damping with coefficient and excited by stochastic forces. Appropriate equations of motion are written as follows: N N N NN N N N u u u u tu u u u t (3) At this stage, the forcing functions in (3) will also stay generic, under assumption of zero mean and lack of correlation: ( ) 0, ( ) ( ) 0 t t t (4) Moreover, to simplify the expressions, we set ( ) ( ), ( ) 0 t t t . The care of original problem may be taken by appropriate superposition. For sufficiently long chains, one can admit that the heat flux is realized only by waves that belong to the propagation zone of the chain. Thus, for fragments of the chain far from the terminal particles, the displacements of the particles can be expressed as linear combinations of the waves travelling to the left and to the right, with the frequencies belonging to the propagation zone: exp( )( ( ) exp( ) ( ) exp( )) (Fragment 1)exp( )( ( ) exp( ) ( ) exp( )) (Fragment 2) n Pn P u i t iqn iqn du i t iqn iqn d (5) In the following, - and q - dependence of all functions will be suppressed, where it will not cause the misunderstanding. It is easy to derive that the temperatures and the heat fluxes through the chain fragments (T , T , J and J respectively) are expressed as follows: , ; , ; cos 2 grP P grP Pgr T d J v dT d J v dd qv dq (6) Here gr v is the group velocity in the chain. Note that the inclusion can, in principle, include the damping elements, and therefore the heat fluxes in fragments 1 and 2 are not necessarily equal. Equations (3) lead to the following boundary conditions for Fourier components of the heat flux in the propagation zone: exp( ) exp( )exp( ) exp( ) 0exp( ) iqN c iqN ciqN c iqN cc iq i (7) Here Ξ denotes the Fourier transform of the forcing function. To close System (7) we recall that the inclusion is assumed linear and time-independent. Therefore, the amplitudes of incoming and outgoing waves in the chain fragments are related through certain transfer-matrix, defined as , g gg g G G (8) Then, System (7) is rewritten in the form ( exp( ) exp( ) ) ( exp( ) exp( ) )exp( ) exp( ) 0 g iqN c g iqN c g iqN c g iqN ciqN c iqN c (9) System (9) is solved straightforwardly: exp( ) exp( ),exp(2 ) ( ) exp( 2 ) iqN c iqN cD DD g iqN c g g c c g iqN c (10) To achieve certain simplification, it is instructive to explore the properties of transfer matrix G in more depth. Let us consider the reciprocal system, in which the forcing is applied to the right end of the chain. Then, System (9) is substituted by the following equations: ( exp( ) exp( ) ) ( exp( ) exp( ) ) 0exp( ) exp( ) g iqN c g iqN c g iqN c g iqN ciqN c iqN c (11) According to Raleigh reciprocity theorem [25], the displacement of the rightmost particle according to System (9) is equal to the displacement of the leftmost particle according to System 11: exp( ) exp( ) exp( ) exp( ); iqN iqN iqN iqN G (12) From (10-12), after some simple algebra, one obtains:
11 22 12 21 det 1 g g g g G (13) This conclusion is valid, even if the inclusion is not symmetric and contains the damping elements, i.e. is not conservative. Further simplification is achieved if one assumes that the inclusion is conservative. In this case, due to the energy conservation, in the stationary regime the heat fluxes through Fragments 1 and 2 are equal. Then, for any forcing configuration the relationship (14) must hold identically. Therefore, the elements of the transfer matrix obey the following relationships:
1; 1; g g g g g g g g (15) It is easy to derive that the most generic transfer matrix that respects conditions (13) and (15) may be represented with the help of the following convenient parametrization: cosh exp( ) sinh exp( ) ; [0, ), , [ , )sinh exp( ) cosh exp( ) x i x i xx i x i G (16) For further derivations, we'll assume that the inclusion does not contain the damping elements, and therefore the transfer matrix is described by parametrization (16). It is interesting to note that from equations (15) one can derive that det 1 G , but cannot obtain stronger condition (13). With account of solution (10) and parametrization (16), expressions (6) for the temperatures and the heat flux in the case of conservative inclusion are presented as follows: ; (cosh sinh )( ) ; sinh 2 ( exp ( 2 ) exp ( 2 )); cosh exp (2 ) 2 sinh sin cosh exp ( 2 grPPP c cJ v dDx x c cT dD x c c i qN c c i qNc cT dDD x i qN c i x c c x i ) qN c (17) In the limit of very long chains, N one can further simplify expressions (17) by averaging over rapidly oscillating phase qN in exponents. Details of this straightforward but somewhat awkward procedure are presented in Appendix A. Finally, the heat flux and the average temperatures at both sides of the inclusion are presented in the following general form: ; cosh 2 1 2 tanh sin(cosh sinh )( ) 8 sinh sin ; cosh 2 1 2 tanh sincosh 2 grPP c cJ v dx c c c c xx x c c c c xT dx c c c c c c xc cT x c c c c c c ; 1 2 tanh sin P dx (18) At this stage we adopt, in addition, that the forcing function in equation (3) corresponds to white Gaussian noise, i.e. the self-correlations obey the well-known relations ( ) ( ) 2 ( ) t t T t t (19) Then one sets , 0, TT T T . For the trivial case of the homogeneous chain the transfer matrix is unit, therefore x , and from equations (19) the following well-known results [17] are easily obtained: / 21 2 1 44 h T T TTJ (20) Further treatment requires one to specify the considered inclusion or defect. In the next Section we consider the case of a single isotopic defect in the chain. Kapitza resistance in the case of a single isotopic defect.
Linear chain
Let us consider that the only inhomogeneity in the chain is the single isotopic defect – a particle with mass m at the site n . The chain (beyond the thermostats) is therefore described by the following equations: (1 ( 1)) 2 0 no n n n n m u u u u (21) Then, Fragment 1 corresponds to all particles with n<0, and Fragment 2 – to all particles with n>0. Substituting expansions (5) into (21), one obtains the following relationships: exp( ) exp( ) (2 exp( )) (2 exp( )) iq iq m iq m iq (22) Therefore, the transfer matrix is expressed as follows: iq iq iq iq e e m e m em i q mi q m i q m G (23) The values necessary for evaluation of integrals (18) are expressed as: ( 1) ( 1) 1sinh , sgn4sin 4 2 sin m m mx q q (24) Substituting (24) into integrals (18), one obtains the following expressions: T z z z dzJ z m zT z m z dzT z m z zT m T z dzz m z zT z dzT z m z zT m T z dzz m z z (25) From (25), Kapitza resistance can be calculated as K T TR J . (26) Integrals in (25) can be explicitly evaluated since the denominators contain a cubic polynomial with respect to z . The details of somewhat awkward calculation are given in Appendix B. Dependency of the Kapitza resistance (26) on the defect mass is plotted in Fig. 2, and on the coupling friction - in Fig. 3. FIG. 2. (Color online) Analytical results of Kapitza resistance plotted by varying the defect mass.
1, 1 T FIG. 3. (Color online) Analytical results of Kapitza resistance plotted by varying the coupling friction.
1, 1.2
T m
Fig. 2 indicates that, contrary to the bulk conductivity, the Kapitza resistance is well-defined and finite even in the linear chain. However, from Fig. 3 one learns that the resistance is not local property of the defect – it strongly, and even non-monotonously depends on the friction coefficient in the thermostats. Therefore, one requires nonlinear models to get in line with basic physical intuition concerning the locality of the boundary resistance [22]
Asymptotic behavior of the boundary resistance.
It is instructive to explore asymptotic behavior of the heat flux and the Kapitza resistance in the special cases of weak and heavy isotopic defects. The case of weak defect is defined as m . In this limit, expressions (25, 26) yield 0
21 2 2
20 1 m TT T ; J T ;
40 ( 1) K R m (27) One observes that the resistance analytically depends on the mass mismatch (cf. Fig 2). The case of heavy isotopic defect is defined as
1, 1 1 m Denoting ˆ 8( 1) z m z , a simple expansion yield the following expression for the heat flux: ˆˆ3 1 6 T dz TJ m z m (28) It is somewhat surprising that the heat flux in the basic approximation does not depend on the friction in this limit case. For comparison, we explored the similar limit in FPU , and Lennard-Jones chains with nearest-neighbor interaction. The FPU potential and Lennard-Jones are given by,
12 4 uuV u ;
14 2 uV u (29) where , and are constants. As one can expect, the nonlinearity has significant effect. Still, it is possible to observe that the data for both nonlinear models collapse on the curves somewhat similar to expression (28) (see Fig. 4 and Fig. 5). FIG. 4. (Color online) Heat flux variation with the chain temperature and defect mass for
FPU chain at
1, 1 1 m . Here N
1 FIG. 5. (Color online) Heat flux variation with the chain temperature and defect mass for Lennard-Jones chain at
1, 1, 1 1 m m .Here N One can conjecture that in the case of very heavy inclusion only long-wave phonons can significantly contribute to the heat flux through the defect, and thus the asymptotic behavior for linear and weakly nonlinear regimes is somewhat similar. The exponents from Figures 4 and 5 require more detailed exploration. Temperature distribution at and near the defect. Thermal "dip" in linear and nonlinear models.
In previous Sections global properties of the heat flux and temperature drop were addressed. Here we consider the details of temperature distribution at the defect site and nearby. According to (5, 6) the local temperature at each chain site beyond the defect is expressed as: exp exp exp exp nn T inq inq dT inq inq d (30) After expanding the inner part of the integral, one obtains exp exp exp 2 exp 2exp exp exp 2 exp 2 inq inq inq inqinq inq inq inq (31) Straightforward simplification of (30) yields the following set of equations 2 n T z m zz m z zz m z m z nqT z m z zm z z z z nqz m z z m z z , if 11 4 128 11 4 64 1 11 1 128 1 2 1 4 4 sin6 1 11 1 4 64 1 sin 21 1 1 4 64 1 1 dz mz m zz m z zz m z m z nqT z m z zm z z z z nqz m z z m z z
22 22 22 2 2 6 21 22 2 2 4 20 122 22 2 2 2 2 6 20 , if 11 41 4 64 1 1 , if 11 1 4 64 1 sin 21 1 ( 1 4 64 1 ) 1 n dz mzz m z zT dz mm z z z z nqz m z z m z zT
22 22 22 2 2 6 21 22 2 2 4 20 122 22 2 2 2 2 6 2 zz m z zT dz mm z z z z nqz m z z m z z (32) Here m zz m z and q z (33) At the defect site itself, the temperature is expressed as 3
22 21 2 22 2 2 6 20 22 2 2 4 2 222 22 2 2 2 2 6 222 21 2 22 2 2 6 20 22 2 2 4 2 222 2 2 20 n T zmT z m z zm m z z z z dz mz m z z m z zzmT z m z zm m z z z zz m z z , if 164 1 1 dz mm z z (34) For the case of very small mass detuning, m , one obtains the following expression for the temperature at the defect:.
21 20 2 2 2221 20 2 220 2 n m m zmT dT z mz m z zmm m zmT dz mz m z zm (35) These expressions simplify to , if 1, if 12 n mT mT T m (36) This peculiar non-analytic dependence of the defect temperature on m calls for closer exploration. In Fig. 6, T , T and n T are presented versus the defect mass m . 4 FIG. 6. (Color online) T , T and n T versus m (
1, 1 T ). T and T are symmetric to both sides of m , whereas n T is not and it is less than T (thermal dip) for m . It is clear from Fig. 6 that for the case m (lighter isotopic impurity), the temperature of the defect is less than T , therefore the temperature profile exhibits a thermal "dip" at the light defect. Such dips, apparently inconsistent with the Second Law of thermodynamics, appear due to lack of modal interactions in the linear models [23,24]. To explain the appearance of the dip, one notes that for m , the spectrum of the chain with defect contains a localized mode in the attenuation zone. For chain long enough, the thermostats cannot excite this localized mode, thus causing the dip. To illustrate this point, we simulate the modal energy distribution [26] for rather short chain with N and light isotopic defect placed at the middle or close to one of the thermostats. The results are presented in Fig. 7. FIG. 7. (Color online) The temperature profile (above) and the energy spectrum (below). The isotopic defect is placed at the middle of the chain at n (left) and near the hot thermostat at n (right),
51, 1.1, 0.9, 1, 0.5
N T T m . 5 In Fig. 7, when the defect is at the middle of the chain, the localized mode, is almost not excited. As the defect is closer to the thermostat, a peak replaces the dip. In this case, the localized mode is strongly excited. This latter case is not described by the previous analytic treatment, since the transition to the infinite chain length is irrelevant here. The aforementioned anomalies of the thermal profile are caused by peculiarities of the frequency spectrum and by lack of interaction between the modes in the linear system. Thus, one can expect that the nonlinearity will remove, or, at least, significantly modify the temperature distribution at the defect. To check this assumption, we simulated the FPU chain [27] with the light isotopic defect. In this case both thermostats are present, the simulation parameters are specified in the figure caption. Time dynamics of the thermal profile in this system is presented in Fig. 8. FIG. 8. (Color online) Destruction of the dip in the temperature profile in the FPU chain. Each plot shows temperature profile after time steps (1), time steps (2), time steps (3), and time steps ( 4) respectively,
51, 1 0.1, 1, 0.5, =0.1
N T m . Fig. 8 demonstrates that at relatively short time scale the dip is observed even in the nonlinear chain. As the simulation evolves, the dip gets gradually destroyed and the defect particle temperature gets lifted close to the average chain temperature. Simple explanation is that as the nonlinearity is relatively small, the initial thermal profile is established through approximately linear dynamics. Nonlinearity reveals itself at higher time scale, related to the temperature. For lower temperatures, the process of the dip elimination is expected to take more time. In the FPU chain one can discuss two possible basic mechanisms of the dip elimination - resonance of the renormalized waves and excitation of the defect by mobile discrete breathers, existing in the frequency range of the localized mode [28,29]. Detailed 6 analysis of these mechanism is beyond the scope of this paper. We report on some numeric observations related to the dip elimination in two different settings – in the chain with only one thermostat (equilibrium simulation), and with temperature gradient and nonzero heat flux (non-equilibrium). In both settings, initially the dip exists, and is gradually eliminated; however, the dynamics of this elimination demonstrates substantial differences. The easy way to observe the breather is to observe the spatiotemporal profile for the frequency range of interest. Such picture is obtained by filtering out all unnecessary frequencies from the time history of the system. We consider a one-dimensional chain of length N with light isotopic defect at the middle. As it was mentioned above, two systems are considered: the non-equilibrium system with N , fixed boundary and thermostats at each boundary, and the equilibrium system with N , periodic boundary conditions and one thermostat in the closed chain. The Langevin thermostat with a coupling friction , temperature T is considered, according to (4). To improve the visibility, rather light defect with m is taken. The system is run initially with (linear, see (29)) for steps, and subsequently the nonlinearity is switched on with . Immediately from the first step, for all particles, the particle-position history is recorded for steps. Then, Fourier transform is taken for these data and all frequency components other than in the range of linear localized mode are filtered out. Here, the linear frequency of the localized mode is equal to 6.4 and the filtering range is . After filtering, the inverse Fourier transform is taken. In Fig. 9, the interaction between a moving breather and the defect particle is clearly observed for the non-equilibrium simulation. The energy of the defect particle increases after the interaction with the breather. In the case of single thermostat (equilibrium simulation) the breather propagation has not been observed. 7 FIG. 9. (Color online) Breather starts from the thermostat and excites the defect particle in the non-equilibrium system. Frequency filtering range , m T . To quantify the dynamics of the dip elimination, the normalized temperature of the defect particle is defined as follows: max Ave DefAve Def
T TT T T (37) Ave T is the average temperature of the system and Def T is the temperature of the defect particle, as functions of time. First, the non-equilibrium case is considered. To suppress the fluctuations, the average of 1800 realizations is taken. The excitation profile is studied for 20 different cases of the defect mass. The latter were chosen to yield particular linear localized frequencies, selected between loc , with 0.1 increment. The time series for * T in the case of the non-equilibrium simulation are presented in Fig. 10 and for the equilibrium simulation - in Fig. 11. In general, these time series have three regions. The first region is almost horizontal, where the nonlinearity almost does not reveal itself. Then one observes the excitation region where the defect particle receives energy as the time progresses. In the third region, the defect particle reaches a steady state at which T starts fluctuating around positive and negative values. One can observe that in the non-equilibrium case the temperature dip is removed at shorter time scale. It is possible to conjecture that the process is expedited by interaction of the defect with mobile breathers created at the thermostats. For the equilibrium case, the excitation profile shows linear trend with the evolution time for all simulated values of defect masses (cf. Fig. 11) and temperatures (cf. Fig. 12). 8 FIG. 10. (Color online) Normalized defect temperature for the non-equilibrium simulations, for 20 different defect masses, m T FIG. 11. (Color online) Normalized defect temperature for the equilibrium simulations, for 20 different defect masses, m T FIG. 12. (Color online) Normalized defect temperature for the equilibrium simulations in temperature interval T . As in Fig. 11, the excitation profile follows linear trend. m It is well known [17] that in linear chains the temperature profile near the thermostats exponentially achieves the average values. This transient occurs due to the waves irradiated by the thermostat, with frequencies belonging to the attenuation zone of the chain. In the present case, in the vicinity of the isolated defect, one also observes the transient to the constant temperature (see, e.g. Fig. 7), despite the fact that in the case of very long chain all frequencies active in this region belong to the propagation zone. 9 To explain this peculiarity of the temperature profile around the defect site, one defines ;Δ Δ n n n n
T T T T T T , (38) Then, with account of (32), it is easy to obtain
22 2 2 5 12 22 2 2 6 21 22 2 2 4 20 122 22 2 2 2 2 6 222 2 2 5 1100 n z m z m z nqz m z zT dz mm z z z z nqz m z z m z zz m z mT zT
21 4 64 1 1 , if 11 1 4 64 1 sin 21 1 1 4 64 1 11 1 4 64 1 sin 2 1 1 n nqz m z z dz mm z z z z nqz m z z m z zm z z z z nqT z m zT , if 11 4 64 1 11 1 4 64 1 sin 2 , if 11 1 1 4 64 1 1 dz mz m z zm z z z z nqT dz mz m z z m z z (39) We consider the case
0, 1 n m , with trivial generalization for other cases. To proceed with the analytic evaluation, one assumes m and obtains
222 2 210 22 2 n m z n zz m z zT dzm z n zz mT z (40) 0 We denote m and sin 2 qz , and obtain.
22 2 20 2 2 22 02 sin 2 cos 21 sin sin2 2 2 sin cos2 2 sin 21 sin sin2 21 cos cos 21 cos 4 1 sin sin 21 cos n q nqq qT dqq q nqq qq nqqT dqq nqqT (41) Here . Then, equation (41) is transferred to the complex plane by substituting exp i , and then reduced to the form (2Re 2 1 8 1 Im 2 1 n n nn nn iTT d (42) The poles in (42) are and the one in which inside the unit circle is , then the expression simplifies to, nn n nn n TT (43) The result means that the near - field of the defect follows staggering exponential decay to the average value. Interestingly, the exponent is governed solely by the mass of defect. In Fig. 1 13, the exponential decay to the average temperature is illustrated by the numerical integration of (37). FIG. 13. (Color online) The absolute value of n T versus particle number, semi-logarithmic scale.
1, 1.2, 1.6
T m . Concluding remarks.
For linear chain with arbitrary local conservative LTI defect, we presented the explicit exact expressions for the heat flux and the temperature drop. Therefore, the Kapitza resistance in this sort of models is well-defined. However, the resistance is shown to depend on particular parameters of the Langevin thermostats; therefore, it is not local property of the defect, as one would expect in physical reality. Other salient feature of the linear model is the thermal dip at the light defect site, related to insufficient excitation of the localized mode by the thermostat, and by lack of the intermodal interactions. When the nonlinearity is switched on, the dip gradually disappears, thus repairing the apparent "violation" of the Second Law. Time dynamics of the dip removal is different in the equilibrium and non-equilibrium simulations. We conjecture that the difference is related to creation of propagating breathers in the non-equilibrium system, that can interact with the defect and further excite it. Intrinsic dynamics of this process of the dip removal requires further exploration. It is important to mention that the dip indeed appears in the temperature profile even for the nonlinear system, although for relatively short time. This time scale is governed by the system size and the defect mass, as the linear substructure is non-dimensional and normalized. The time scale necessary to reveal the nonlinearity strongly depends on the temperature. Thus, if the temperature is low enough, the initial response of the system to the imposed thermal 2 gradient will be governed by the linear dynamics. In other terms, the linear dynamics may be useful to understand the short-time response of more realistic models, although the steady-state regime will require full account of the nonlinearity.
Acknowledgment
The authors are very grateful to Israel Science Foundation (grant 1696/17) for financial support.
References [1] P. L. Kapitza, J. Phys. USSR , 181 (1941). [2] G. L. Pollack, Kapitza Resistance , Rev. Mod. Phys. , 48 (1969). [3] N. S. Snyder, Heat Transport through Helium II: Kapitza Conductance , Cryogenics , 89 (1970). [4] L. J. Challis, Kapitza Resistance and Acoustic Transmission across Boundaries at High Frequencies , J. Phys. C Solid State Phys. , 481 (1974). [5] E. T. Swartz and R. O. Pohl, Thermal Boundary Resistance , Rev. Mod. Phys. , 605 (1989). [6] I. M. Khalatnikov, An Introduction To The Theory Of Superfluidity (CRC Press, 2018). [7] W. A. Little,
The Transport of Heat between Dissimilar Solids at Low Temperatures , Can. J. Phys. , 334 (1959). [8] S. Simons, On the Thermal Contact Resistance between Insulators , J. Phys. C Solid State Phys. , 4048 (1974). [9] R. E. Peterson and A. C. Anderson, The Kapitza Thermal Boundary Resistance , J. Low Temp. Phys. , 639 (1973). [10] P. Erdös and S. B. Haley, Low-Temperature Thermal Conductivity of Impure Insulators , Phys. Rev. , 951 (1969). [11] H. Budd and J. Vannimenus,
Thermal Boundary Resistance , Phys. Rev. Lett. , 1637 (1971). [12] W. M. Saslow, Kapitza Conductance, Temperature Gradients, and Solutions to the Boltzmann Equation , Phys. Rev. B , 2544 (1975). [13] T. J. Sluckin, A Simple Model for the Kapitza Conductance , Phys. Lett. A , 390 (1975). [14] Ch. Steinbrüchel, The Scattering of Phonons of Arbitrary Wavelength at a Solid-Solid Interface: Model Calculation and Applications , Z. Für Phys. B Condens. Matter Quanta , 293 (1976). [15] M. E. Lumpkin, W. M. Saslow, and W. M. Visscher, One-Dimensional Kapitza Conductance: Comparison of the Phonon Mismatch Theory with Computer Experiments , Phys. Rev. B , 4295 (1978). [16] E. T. Swartz and R. O. Pohl, Thermal Resistance at Interfaces , Appl. Phys. Lett. , 2200 (1987). [17] Z. Rieder, J. L. Lebowitz, and E. Lieb, Properties of a Harmonic Crystal in a Stationary Nonequilibrium State , J. Math. Phys. , 1073 (1967). [18] A. Dhar and D. Roy, Heat Transport in Harmonic Lattices , J. Stat. Phys. , 801 (2006). [19] A. Casher and J. L. Lebowitz,
Heat Flow in Regular and Disordered Harmonic Chains , J. Math. Phys. , 1701 (1971). [20] V. Kannan, A. Dhar, and J. L. Lebowitz, Nonequilibrium Stationary State of a Harmonic Crystal with Alternating Masses , Phys. Rev. E , 041118 (2012). [21] T. Mai, A. Dhar, and O. Narayan, Equilibration and Universal Heat Conduction in Fermi-Pasta-Ulam Chains , Phys. Rev. Lett. , 184301 (2007). 3 [22] J. Paul and O. V. Gendelman, Kapitza Resistance in Basic Chain Models with Isolated Defects , Phys. Lett. A 126220 (2020). [23] X. Cao and D. He,
Interfacial Thermal Conduction and Negative Temperature Jump in One-Dimensional Lattices , Phys. Rev. E , 032135 (2015). [24] Y. Liu and D. He, Anomalous Interfacial Temperature Profile Induced by Phonon Localization , Phys. Rev. E , 062119 (2017). [25] J. W. S. Rayleigh, Treatise on Sound Vol II (London: Macmillan), 1878). [26] L. Meirovitch,
Fundamentals of Vibrations (McGraw-Hill, Boston, Mass, 2001). [27] E. Fermi, P. Pasta, S. Ulam, and M. Tsingou, Studies of the Nonlinear Problems, Los Alamos Scientific Lab., N. Mex., 1955. [28] B. Gershgorin, Y. V. Lvov, and D. Cai,
Renormalized Waves and Discrete Breathers in β -Fermi-Pasta-Ulam Chains , Phys. Rev. Lett. , 264302 (2005). [29] N. Li, B. Li, and S. Flach, Energy Carriers in the Fermi-Pasta-Ulam β Lattice: Solitons or Phonons? , Phys. Rev. Lett. , 054102 (2010).
APPENDIX A: AVERAGING IN EXPRESSIONS (17)
First, we denote qN . Averaging in the limit N is performed with the help of the following expressions:
22 2 20 2*2 * * *224 * *2 1exp( ) 1 2 1 2 i ii iN i ii izz i dxe e c i x c c xe e cD xe e c i x c c xe e ci zdzz z z z zz zzx c (A1) In (A1) z are the roots of polynomial ( ) cosh 2 sinh sin cosh i i P z xe c z i x c c z xe c : sinh sin cosh sinh sincosh exp( ) i x x xcz c x i (A2) It is easy to demonstrate that both roots lie inside or at the boundary of the unit circle: c qz c q (A3) Therefore, only these two roots contribute to the integral in (A1), two others lie outside the unit circle. Thus, one obtains: N z zD x c z z z z z z z z z zc cx c c c c c c x
4 (A4) The other auxiliary expression is evaluated in a similar way:
22 2 2 20 2*2 * * *224 * *2 1exp( ) 1 2 1 2
12 (cosh 2 sinh sin cosh )(cosh 2 sinh sin cosh )( )( )( 1)( 1)2 cosh iqN ii ii iN i ii izz i e e dxe e c i x c c xe e cD xe e c i x c c xe e ci dzz z z z zz zzx c * *1 2 1 22 4 4 2 22 2 21 2 2* 1 22 2 4 4 2 22 2 2 2 cosh 1 ( ) 2 1 2 tanh sin2 exp( ) tanh sincosh 2 1 2 tanh sin z z z zcx z z c c c c xc ic c i xx c c c c c c x (A5) APPENDIX B: EXACT VALUES OF INTEGRALS IN EQ. (25)
Let us consider the following part from (25) containing polynomial in the denominator,
11 4 64 1
I z m z (B1) Put z
22 2 3 22 2 2 2
I m m m m (B2) If , , be the roots of the polynomial in the denominator of (B2), then,
164 1
I m (B3) Let us write the depressed cubic equation by introducing the variable r ,
12 1 r m (B4) 5
223 2 22 2 32 22 4 2 2
11 38 1 12 164 1 12 12 1 96 1 64 1
I r r m mm m m m (B5) If we define,
222 2 p m m ;
32 22 4 2 2
12 12 1 96 1 64 1 q m m m , (B6) then one root of the cubic equation can be written as trigonometric and hyperbolic solutions. , >0 an1 3 32 sinh arcsin , if 4 27 03 3 212 1 31 32 cosh arccos if 4 27 03 3 212 1 1 3 32 cos arccos , if 3 d 43 212 >0 and 1 p qh p q pp pm q qp h p q pq p pm p q pp pm q (B7) Then the remaining roots can be found as follows,
164 1
I m b c (B8) b m ;
164 1 c m ; b b c ; b b c Now, it is easy to do the partial fraction following the exact integration.
11 164 1 1
I m (B9) Putting (B9) back to (25) yields 6 z z zzz z zTJ dzzm z z zz (B10) J TmTm (B11) Similarly we can calculate T and T
81 3 4 8 181 3 4 82 16 181 3 4 8 1
T TT (B12) 7
81 3 4 8 181 3 4 82 16 181 3 4 8 1
T TT