Equilibration of sinusoidal modulation of temperature in linear and nonlinear chains
Elena A. Korznikova, Vitaly A. Kuzkin, Anton M. Krivtsov, Daxing Xiong, Vakhid A. Gani, Aleksey A. Kudreyko, Sergey V. Dmitriev
aa r X i v : . [ n li n . PS ] J a n Equilibration of sinusoidal modulation of temperature in linear and nonlinear chains
Elena A. Korznikova , , ∗ Vitaly A. Kuzkin , , † Anton M. Krivtsov , , ‡ Daxing Xiong , § Vakhid A. Gani , , ¶ Aleksey A. Kudreyko , ∗∗ and Sergey V. Dmitriev , †† Institute of Molecule and Crystal Physics, Ufa Federal ResearchCentre of the Russian Academy of Sciences, Ufa 450054, Russia Ufa State Aviation Technical University, Ufa 450008, Russia Peter the Great Saint Petersburg Polytechnical University, Saint Petersburg 195251, Russia Institute for Problems in Mechanical Engineering, RAS, Saint Petersburg 199178, Russia MinJiang Collaborative Center for Theoretical Physics,Department of Physics and Electronic Information Engineering,Minjiang University, Fuzhou, Fujian 350108, China National Research Nuclear University MEPhI (Moscow Engineering Physics Institute), Moscow 115409, Russia Institute for Theoretical and Experimental Physics of NationalResearch Centre “Kurchatov Institute”, Moscow 117218, Russia Department of Medical Physics and Informatics,Bashkir State Medical University, Ufa 450008, Russia Institute of Mathematics with Computing Centre,Ufa Federal Research Centre of RAS, Ufa 450008, Russia
The equilibration of sinusoidally modulated distribution of the kinetic temperature is analyzedin the β -Fermi-Pasta-Ulam-Tsingou chain with different degrees of nonlinearity and for differentwavelengths of temperature modulation. Two different types of initial conditions are used to showthat either one gives the same result as the number of realizations increases and that the initialconditions that are closer to the state of thermal equilibrium give faster convergence. The kineticsof temperature equilibration is monitored and compared to the analytical solution available for thelinear chain in the continuum limit. The transition from ballistic to diffusive thermal conductivitywith an increase in the degree of anharmonicity is shown. In the ballistic case, the energy equi-libration has an oscillatory character with an amplitude decreasing in time, and in the diffusivecase, it is monotonous in time. For smaller wavelength of temperature modulation, the oscillatorycharacter of temperature equilibration remains for a larger degree of anharmonicity. For a givenwavelength of temperature modulation, there is such a value of the anharmonicity parameter atwhich the temperature equilibration occurs most rapidly. PACS numbers: 05.45.Yv, 63.20.-eKeywords: harmonic chain, nonlinear chain, temperature, thermal equilibration
I. INTRODUCTION
The miniaturization of electronic and microelectrome-chanical systems has led to the need to control heat fluxesat the micro- and nanolevels, where the laws establishedfor macrosystems turned out to be inaccurate; for ex-ample, the Fourier law of thermal conductivity is notvalid [1–6]. On the other hand, the possibility of creat-ing nanostructured metamaterials with desired proper-ties endows researchers and technologists with new leversof heat flow control. Heat in materials is transferredmainly by phonons, and phononics is a rapidly growingfield of knowledge at the intersection of physics, mate-rials science, and nanotechnology, studying phonon en- ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] § Electronic address: [email protected] ¶ Electronic address: [email protected] ∗∗ Electronic address: [email protected] †† Electronic address: [email protected] ergy transport and its applications [7, 8]. Major achieve-ments in this field are the development of thermal tran-sistors [9, 10], thermal diodes [11–15], and thermal logicdevices [16–18]. Such developments require a better the-oretical and experimental understanding of anomalousthermal transport in miniature low-dimensional systems.The Fourier law of thermal conductivity states thatfor macroscopic bodies, the heat flux is proportional tothe temperature gradient with a proportionality coeffi-cient κ , referred to as thermal conductivity. In a num-ber of theoretical works [19–31], it was shown that inone-dimensional structures, Fourier’s law does not work,in the sense that κ depends not only on the material,but also on the dimensions of the conductor, in particu-lar, on its length L according to the power law κ ∼ L α ,where the exponent is in the range 0 ≤ α ≤
1. The case α = 1 corresponds to ballistic thermal conductivity, andfor α = 0 we have normal, diffusive thermal conductivityobeying the Fourier law. If 0 < α <
1, we have anoma-lous thermal conductivity. Note that defect-free linearsystems of any complexity always demonstrate ballisticheat propagation [32–40]. One might expect that in sys-tems with weak anharmonicity, the linear theory can bevery helpful [41].It was shown that the ballistic propagation of heat inlinear discrete systems can be described with high accu-racy by continuum equations [38, 41]. In particular, equi-libration of sinusoidal modulation of temperature wasdescribed analytically [38], and it was shown that tem-perature oscillations in the α -Fermi-Pasta-Ulam-Tsingou( α -FPUT) chain lead to the excitation of mechanical vi-brations. This physical phenomenon is referred to as theballistic resonance [41].For the description of anomalous thermal conductivity,the nonlocal fractional-type diffusion equation [20] andnonlinear fluctuating hydrodynamics [42, 43] have beendeveloped.The equilibration of sinusoidal modulation of temper-ature was studied earlier in anharmonic chains [44, 45].The authors found that the equilibration of the short-wavelength modulation of temperature occurs throughoscillations, while the long-wavelength modulation isequilibrated monotonically. Here we give explanationsfor those observations, demonstrating that the oscillatoryregime of temperature equilibration can be described bythe linear theory and that the oscillations of the long-wavelength modulations are suppressed by a weaker an-harmonicity as compared to short-wavelength modula-tions.Here the β -FPUT chain with symmetric anharmonic-ity is considered [46] because it is free of thermal expan-sion and the mechanical and thermal oscillations are notcoupled [41].The main goal of this study is to analyze the transi-tion from ballistic [38, 41] to diffusive [20, 44, 45] energyequilibration as the degree of anharmonicity increases.We describe the β -FPUT model and simulation setupin Sec. II, numerical results are presented in Sec. III, andconclusions are drawn in Sec. IV. II. THE MODEL AND SIMULATION SETUPA. The Fermi-Pasta-Ulam-Tsingou chain
We consider the β -FPUT chain of particles [46] havingmass m , whose dynamics is defined by the Hamiltonian H = K + P = X n K n + X n P n , (1)which is the sum of the kinetic ( K ) and potential ( P )energies of the chain with the kinetic, potential, and totalenergies of individual particles being equal to K n = m ˙ u n , (2) P n = k (cid:0) u n − u n − (cid:1) + β (cid:0) u n − u n − (cid:1) + k (cid:0) u n +1 − u n (cid:1) + β (cid:0) u n +1 − u n (cid:1) , (3) H n = K n + P n , (4) respectively. Here u n ( t ) is the displacement of the n thparticle from its equilibrium position, which is an un-known function of time t , and ˙ u n ≡ du n /dt is its velocity.The particles are coupled to their nearest neighbors bythe potential which includes the quadratic term with theharmonic force constant k and the quartic term with theanharmonic force constant β .The equations of motion that stem from Eqs. (1)–(3)are m ¨ u n = k (cid:0) u n − − u n + u n +1 (cid:1) − β (cid:0) u n − u n − (cid:1) + β (cid:0) u n +1 − u n (cid:1) . (5)Without loss of generality, we set m =1, k = 1, andtake different values for β . The lattice spacing is h =1. The chain of N particles with the periodic boundaryconditions ( u n = u n + N ) is considered.In the case of small amplitude vibrations, one can ne-glect the nonlinear term by setting β = 0 and find thesolutions of the linearized equation (5) in the form ofnormal modes u n ∼ exp[ i (2 πqn/N − ω q t )] with the wavenumber q = 0 , , ..., N/ ω q . These modesobey the following dispersion relation: ω q = 2 r km sin πqN . (6)The considered chain supports the small-amplitude run-ning waves (phonons) with frequencies within the bandfrom ω min = 0 for q = 0 to ω max = 2 p k/m for q = N/ T = ¯ K = 1 N X n m h ˙ u n i , (7)is used, where h·i denotes the mathematical expectation. B. Two types of initial conditions
Since we deal with temperature and take a relativelysmall number of particles, the physical picture emergesas a result of averaging over many realizations.We aim to set the initial conditions that create, afteraveraging over many realizations, the initial distributionof the total energy over the particles of the form h H n i = H b + ǫ h (cid:16) πnN (cid:17)i , (8)where H b is the background level of the total energy andthe term with a multiplier ǫ ≪ H b adds the desired si-nusoidal modulation.The initial total energy distribution (8) can be achievedin many different ways. We will compare the results fortwo types of initial conditions. The small sinusoidal ad-dition, in both cases, will be taken in the form of kineticenergy as described below. On the other hand, the back-ground energy, which constitutes the main part of thetotal energy of the system, will be introduced either askinetic energy or it will be nearly equally shared betweenthe kinetic and potential forms. Initial conditions of the first type.
We assign randomvelocities to the particles with zero initial displacementsso that the initial total energy includes only kinetic en-ergy, while the potential energy is zero. To do so, in eachrealization, we set the initial velocities of the particlessuch that˙ u n = ρ n r m r H b + ǫ h (cid:16) πnN (cid:17)i , (9)where ρ n is a random variable with the standard normaldistribution having zero expectation and unit variance.For generation of random variable ρ n , the standard nor-mal distribution given by the probability density function P ( x ) = exp( − x / / √ π is used. Then the desired spa-tial distribution of the total energy at t = 0 is achieved, h H n i = h K n i = m h ˙ u n i = H b + ǫ h (cid:16) πnN (cid:17)i . (10) Initial conditions of the second type.
In this case, thebackground total energy H b is obtained by summing upall N/ u n = b N/ cos( πn )+ N/ − X q =1 b q cos (cid:16) qπnN ± ω q t + δ q (cid:17) , (11)with ω q given by Eq. (6), with random phase shifts δ q uniformly distributed in the domain (0 , π ), and with theamplitudes b q chosen such that each harmonic has sametotal energy equal to H b / ( N/ ω q t term is taken with equal probability inorder to have equal contribution to the energy from thewaves running to the right and to the left.For running waves in the linear lattice ( β = 0), thekinetic and potential energies are exactly equal, whilefor β >
0, there appears a deviation from equality due tothe effect of nonlinearity.The sinusoidal modulation of the energy distributionalong the chain is achieved for each realization by increas-ing the particle kinetic energies by∆ K n = ǫ h (cid:16) πnN (cid:17)i , (12)so that each particle having velocity ˙ u n gets the velocityincrement, ∆ ˙ u n = ± r ˙ u n + 2 m ∆ K n − ˙ u n , (13)where the upper (lower) sign is for positive (negative) ˙ u n . Summary . In the initial conditions of the first type, thetotal energy of the chain at t = 0 is in the form of kineticenergy with the potential energy being exactly zero. In the initial conditions of the second type, the backgroundenergy at t = 0 is almost equally shared between thekinetic and potential forms (exactly equally in the linearcase) and a small amount of the kinetic energy is addedto achieve the desired sinusoidal modulation of the totalenergy distribution. C. Analytical solution for the linear chain
Heat transfer in the linear chain can be described witha high accuracy by the continuum equation [38]¨ T + 1 t ˙ T = c T ′′ , (14)where T ( x, t ) is the temperature field and c = h r km , (15)is the sound velocity.For the initial conditions T = T b + ∆ T sin( λx ) , ˙ T = 0 , (16)where λ = 2 πL , (17)and L is the modulation wavelength, the solution toEq. (14) reads [38] T = T b + A ( t ) sin( λx ) , A ( t ) = ∆ T J ( ωt ) , (18)where ω = λc, (19)and J is the Bessel function of the first kind.Initially, the left half of the chain, 0 ≤ x ≤ π/λ , has anaveraged temperature greater than the right half, π/λ There are two ways to address the effect of anhar-monicity on the dynamics of the chain. The first ap-proach is to fix the anharmonicity parameter β and studysystems with different levels of energies (different temper-atures). The second approach is to fix the energy of thesystem and to change β . In the present study, the secondmethod is used.The equations of motion given by Eq. (5) are integratednumerically using the symplectic, sixth-order St¨ormermethod with the time step τ = 10 − p m/k . The accu-racy of integration is controlled by monitoring the totalenergy of the system which, in our simulations, is con-served with the relative error not exceeding 10 − withinthe whole numerical run.In the simulations, we take H b = 1 . ǫ = 0 . β including the case of the linear chain with β = 0. The ini-tial conditions are taken in one of two forms, as describedin Sec. II B.As mentioned above, initially the left half of the chainhas an averaged temperature T L greater than the righthalf, T R . In order to study the kinetics of temperatureequilibration, the difference between these temperatures, δT = T L − T R = 2 N N/ − X n =0 h K n i − N N − X n = N/ h K n i , (24)is calculated as a function of time and compared to theprediction of the linear theory, given by Eq. (23). Notethat Eq. (24) is the discrete version of Eq. (20).We also carry out simulations for chains of differentlengths N . We emphasize that in most calculations, thetemperature modulation wavelength is equal to the chainlength N and, only at the end of Sec. III C, we considerthe case of a chain including several temperature modu-lation wavelengths.The period of the Bessel function J ( ξ ) graduallyshortens, rapidly approaching the value of 2 π . Sup-pose we want to analyze the dynamics of the chain forabout I periods of the Bessel function, in other words,for ξ ≤ πI . Then, in view of Eq. (23), the simulationtime is about t s = IN r mk . (25)Recall that in our simulations, we set m = 1 and k = 1. III. SIMULATION RESULTSA. Characteristics of the initial conditions The initial conditions of the first type, given by Eq. (9),and of the second type, given by Eqs. (11) and (13), are < H n > n (a) < H n > n (b) < H n > n (c) Figure 1: Values of particle energies h H n i in the chain of N = 1024 particles for the initial conditions of the first typeaveraged over (a) 10 , (b) 10 , and (c) 10 realizations. Thered line shows the desired energy distribution (8) for H b = 1and ǫ = 0 . 1. Here the chain is linear, β = 0. < H n > n (a) < H n > n (b) < H n > n (c) Figure 2: The same as in Fig. 1, but for the initial conditionsof the second type. stochastic, and it is instructive to estimate their statisti-cal characteristics for chains of different lengths N .First, in Figs. 1 and 2, we show the initial distributionof total energy in the linear chain ( β = 0) of N = 1024particles for the initial conditions of the first and secondtypes, respectively, for H b = 1 and ǫ = 0 . 1. The energiesof the particles are averaged over (a) 10 , (b) 10 , and(c) 10 realizations. The red line in each panel shows thedesired distribution of energy, given by Eq. (8).A comparison of the results presented in Figs. 1 and 2shows that the initial conditions of the second type areless stochastic than the initial conditions of the first type.This is understandable because in the initial conditions ofthe second type, the background energy H b is obtainedby summation of phonon modes followed by correctionof the particles’ velocities to achieve the sinusoidal dis-tribution of energy, while in the first type of the initial 100 1000 100001E-30.010.11 s N Figure 3: Standard deviation of the difference between thetemperatures of the left and right halves of the chain, δT , at t = 0 calculated for 10 realizations of the initial conditionsof the first type (dots) and the second type (circles) for thechains of different length N . The dashed line shows the slopeof − / 2. Results for β = 0 (linear case). conditions, all particles have random initial velocities.In this work, energy equilibration in the chain will bemonitored by observing the time evolution of an inte-gral characteristic, namely, the difference between thetemperatures of the left and right halves of the chain, δT = T L − T R . Let us analyze the statistical characteris-tics of δT at t = 0. According to Eq. (23), at t = 0 oneshould have the mean value of δT = 2 ǫ/π . We generatesets of 10 initial conditions for different N and calculatethe standard deviation s of δT for β = 0 (linear case).In Fig. 3, the standard deviation of δT is presented as afunction of the number of particles in the chain N for theinitial conditions of the first type (dots) and the secondtype (circles). The log-log plot shows that s = a/ √ N (the offset dashed line has the slope − / 2) and a = 1 . a = 0 . s is smallerthan for the initial conditions of the first type.Recall that when the initial conditions of the first typeare used, the potential energy of the chain at t = 0 iszero and the total energy is equal to the kinetic energy.According to the exact result obtained for linear chains,the energy exchange between the kinetic and potentialparts follows the Bessel function [39]. In Fig. 4, we plotthe time evolution of the kinetic and potential energiesin the linear chain ( β = 0) of N = 2 = 32768 parti-cles for the initial conditions of the first [Fig. 4(a)] andsecond [Fig. 4(b)] types. The result is the average over100 realizations. In the case of the first (second) initialconditions at t = 0, one has K − P = H b + ǫ = 1 . K − P = ǫ = 0 . P K , P t (a) KPK (b) K , P t Figure 4: Kinetic ( K ) and potential ( P ) energies of the chainas functions of time. Initial conditions of the (a) first and (b)second type are used. The presented result is the averagedone over 100 realizations. The linear chain ( β = 0) includes N = 2 = 32768 particles. at t = 0 are much closer than in the case of the first typeof initial conditions.It should be pointed out that the period of oscillationsof the kinetic and potential energies is about π/ ≈ . N ; see Eq. (25). For example, for achain of N = 2 = 1024 particles, the period of temper-ature equilibration is three orders of magnitude longerthan the oscillation period of the kinetic and potentialenergies.The main idea of this work is to analyze the effectof weak anharmonicity on the equilibration of sinusoidaltemperature modulation, and thus it is important toquantify the proximity of the chain to the harmonic case.It is well known that in a harmonic chain, the kinetic en-ergy is exactly half of the total energy, but this is notthe case in the presence of anharmonicity. Let us use thetotal to kinetic energy ratio, H/K , as the measure of de-viation from the harmonic case, which can also be used asthe measure of heat capacity [48–50]. In Fig. 5, the H/K ratio is shown as a function of the anharmonicity pa-rameter β for the case of total energy per particle equalto H/N = 1, with the number of particles N = 2 =32768. The range β ≤ . β < . 1, the deviation of the H/K ratio from the value corresponding to the harmonic limit( β = 0) is within 5%.The presented analysis allows us to draw the follow-ing conclusions about the two types of the initial condi-tions. (i) The initial conditions of the second type are lessstochastic (see Fig. 3), which means that they should givea better convergence of the results with increasing num-ber of realizations. (ii) Simulations with longer chains H / K Figure 5: Total to kinetic energy ratio as a function of theanharmonicity parameter β = 0 for the case of total energyper particle equal to H/N = 1, with the number of particles N = 2 = 32768. should give a better convergence of the results with in-creasing number of realizations since s decreases withincreasing N ; see Fig. 3. (iii) The initial conditions ofthe second type are closer to thermal equilibrium withcloser values of the kinetic and potential energies in thesystem at t = 0; see Fig. 4. (iv) For β < . 1, the totalto kinetic energy ratio, H/K , differs by no more than5% from 2, which is the value for the harmonic case; seeFig. 5. This means that for β < . 1, we have the regimeof weak anharmonicity. B. Temperature equilibration in the linear chain For the linear chain ( β = 0), the analytical result givenby Eq. (23) is available, according to which the temper-ature difference of the left and right halves of the chainoscillates with a reduction in the time amplitude follow-ing the Bessel function of the first kind. In Fig. 6, thenormalized difference δT = T L − T R as a function of thenormalized time is plotted for N = 16384 [Figs. 6(a) and6(b)] and N = 32768 [Figs. 6(c) and 6(d)]. In Figs. 6(a)and 6(c), the initial conditions of the first type are used,while in Figs. 6(b) and 6(d), those of the second type areused. Black curves show the numerical results and lightgreen curves stand for the analytical prediction given byEq. (23). In all cases, the result of averaging over 50realizations is shown.It can be seen from Fig. 6 with the naked eye that theresults for the longer chain shown in Figs. 6(c) and 6(d)are closer to the theoretical prediction than the resultsfor the shorter chain presented in Figs. 6(a) and 6(b). Inorder to quantify the difference between the numericaland theoretical curves, we calculate the area S betweenthem within the interval 0 ≤ (2 π/N ) t ≤ 30. The resultis (a) S = 0 . S = 0 . S = 0 . S = 0 . S is considerably larger than for (c) and (d),and this is in line with the result presented in Fig. 3,showing that the standard deviation s of δT reduces with T / (a) First type IC, N=16384 T / (b) Second type IC, N=16384 T / (c) First type IC, N=32768 T / (2 /N)t (d) Second type IC, N=32768 Figure 6: Results for β = 0 (linear case). The normalizeddifference between the averaged temperatures of the left andright halves of the chain, δT , is plotted as a function of nor-malized time. The numerical result is shown by a black line,while the theoretical prediction given by Eq. (23) (the Besselfunction of the first kind) is shown by the light green line. Ini-tial conditions (a),(c) of the first type and (b),(d) of the sec-ond type are used. In (a),(b), the chain length is N = 16384,and in (c),(d), it is N = 32768. In all cases, the result ofaveraging over 50 realizations is shown. increasing chain length. It can also be seen that the area S is slightly smaller when the initial conditions of thesecond type are used. This again agrees with Fig. 3,showing that s is smaller for the second type of initialconditions.From the results presented in Fig. 6, one can see thatthe number of realizations needed to achieve certain ac-curacy strongly depends on the chain length N . This isso because we monitor the time evolution of the integralparameter δT , and for a longer chain, this parameter isestimated with a higher accuracy for each particular re-alization. In the following, if the chain length is halved,the number of realizations is at least doubled to get ap-proximately the same accuracy.According to our simulations, for increasing numberof realizations, the kinetics of temperature equilibrationconverges not only for the harmonic chain, but also for β > T / ( ) (2 /N) t Figure 7: Normalized δT as a function of normalized time fordifferent values of the nonlinearity parameter β , as specifiedfor each curve. The initial conditions of the second type areused. The chain length is N = 32768. All curves are theresult of averaging over 50 realizations. *=0.09430.120.140.080.060.040.02 T / (2 /N)t =0.0 Figure 8: A set of curves showing normalized δT as a functionof time for different values of the nonlinearity parameter β inthe vicinity of the first minimum. For β ∗ = 0 . δT curve reaches zero value (shown by thesolid line). The initial conditions of the first type are used.The chain length is N = 512. All curves are the result ofaveraging over 2 . × realizations. It is assumed that for β < β ∗ , the linear theory is still informative and, for larger β , the anharmonicity plays an essential role. (14) is capable of describing the heat flux in linear chains. C. Effect of anharmonicity The effect of the anharmonicity parameter β on equi-libration of the sinusoidal modulation of temperature ispresented in Fig. 7. The time evolution of normalized δT is presented and the value of β is indicated for eachcurve. The chain length is N = 32768. The results are 100 1000 10000 1000000.010.1 N * N * Figure 9: The characteristic value of the anharmonicity pa-rameter β ∗ as a function of the wavelength of temperaturemodulation N . The dashed line shows the slope − / 2. Thetable shows the numerical values of β ∗ for different N . averaged over 50 realizations.It can be seen from Fig. 7 that with increasing β , theamplitude of oscillations of the curves decreases and, at β = 0 . 08, one has an almost monotonically decreasingcurve. The oscillation frequency increases with increasing β . It can also be noted that the fastest temperature equi-libration is observed for β = 0 . 02 because, in this case,the energy difference δT vanishes faster. For smaller β ,the amplitude of oscillations decays slower, and for larger β , the relaxation becomes slower, as compared to the caseof β = 0 . 02. We note the analogy with the damped os-cillator, which relaxes to the equilibrium position veryslowly if the viscosity is very high (overdamped motion)or when it is very small (slowly decaying oscillations), sothat there exists a value of the viscosity parameter withthe fastest relaxation [51].Our next step is to estimate the accuracy of the lineartheory for different wavelengths of the temperature mod-ulation, which, in our simulations, is equal to the chainlengths N . We need to define a characteristic value of theanharmonicity parameter β ∗ such that for β < β ∗ , onecan rely on the linear theory, but for large β , the effectof anharmonicity must be taken into account.As one can see from Fig. 7, with increasing β , theminima of the oscillating curve go up. Let us focus on thefirst minimum of the curve and specify β ∗ as the valueof β when the first minimum of δT reaches zero. Thisdefinition of β ∗ is further illustrated in Fig. 8. Here aset of curves is plotted for different values of β in thevicinity of the first minimum and interpolation betweenthe minimal points allows one to find β ∗ = 0 . δT vanishes (shown by the solidline). Note that for β = β ∗ , the oscillations of δT intime are still pronounced and not yet suppressed by theanharmonicity. This means that for β < β ∗ , the lineartheory is valid in the sense that it explains the oscillatorycharacter of temperature equilibration in the chain. Theresults presented in Fig. 8 are obtained for N = 512,initial conditions of the first type are used, and averagingover 2 . × realizations is performed for each value of β . Note that a larger number of realizations should betaken for short-wavelength temperature modulation.Figure 9 shows how the characteristic value of the an-harmonicity parameter β ∗ depends on the wavelengthof the temperature modulation N . The dashed line inthe log-log plot has the slope − / 2. This means that β ∗ ∼ N − / for large N , and for smaller N the decreaseof β ∗ with increasing N is somewhat faster. The pre-sented result says that the linear theory describes short-wavelength modulations of temperature in the chain withlarger anharmonicity, while for long-wavelength modula-tions the effect of anharmonicity is stronger.So far, the temperature modulation wavelength wasequal to the chain length N . Let us check the effect ofthe chain length taking the modulation wavelength equalto 512 in the chain of 32768 particles, thus having 64 tem-perature modulation periods in the chain. In this case,the temperatures T L and T R are calculated as the sumsover 64 half periods of temperature modulation. Simu-lations with the initial conditions of the first type haveshown that for the long chain, β ∗ = 0 . β ∗ = 0 . β ∗ is nearly the same forthe chains that include one and 64 periods of temperaturemodulation. IV. CONCLUSIONS In the present study, the effect of anharmonicity onthe equilibration of sinusoidal modulation of tempera-ture in the β -FPUT chain was analyzed. The results fordifferent values of the anharmonicity parameter β andfor different wavelengths of temperature modulation N are obtained numerically and compared to the analyticalsolution available for the linear case ( β = 0). Applicabil-ity of the linear theory to a weakly nonlinear chain wasassessed for different wavelengths of temperature modu-lation N . Initial conditions of two types were used: (i) at t = 0, the energy of the system is in the form of kineticenergy with zero potential energy and (ii) the other, ma-jor part of the energy is initially shared between kineticand potential energies.Our main findings can be summarized as follows:(1) For the linear chain ( β = 0), the numerical re-sults averaged over increasing number of realiza-tions converged to the analytical solution given byEq. (18). This solution predicts that equilibrationof sinusoidal modulation of temperature demon-strates oscillations with decrease in time amplitude,following the Bessel function of the first kind. Thiswas true for the initial conditions of both types, though convergence with increasing number of re-alizations was faster for the initial conditions withnearly equal kinetic and potential energies. Con-vergence was also faster for larger wavelength oftemperature modulation N ; see Sec. III A. The ki-netics of temperature equilibration, for increasingnumber of realizations, converges not only for theharmonic chain, but also for β > β , oscillations of temperature de-cay slowly, and for larger β , the monotonic decayis slow; see Fig. 7.(3) Linear theory remains informative for weakly an-harmonic chains when β < β ∗ , with β ∗ defined asshown in Fig. 8. As can be seen from Fig. 9, β ∗ decreases with increasing temperature modulationwavelength N . This means that temperature mod-ulation with short wavelength is less affected by theanharmonicity or, in other words, linear theory re-mains valid for larger values of β , as compared tothe long-wavelength temperature modulation.Overall, our results have confirmed that (i) the con-tinuum equation (14) derived in [38] accurately describesthe temperature flow in linear chains, (ii) linear theory re-mains informative for weakly anharmonic chains, and (iii)short-wavelength modulations of temperature are less af-fected by the anharmonicity and linear theory remainsvalid for larger values of β , as compared to the long-wavelength modulations of temperature.In this regard, the results presented in previousworks [18, 44] have found their explanation. Oscillationsof the short-wavelength sinusoidal temperature modula-tion, observed by the authors of those works, can be wellexplained by the linear theory [38]. The oscillations werenot observed by the authors for long-wavelength tem-perature modulation because, in this case, the effect ofanharmonicity is much stronger. The oscillations of long-wavelength temperature modulation can be observed forsmaller values of the anharmonicity parameter. Acknowledgments