Thermal transport in long-range interacting Fermi-Pasta-Ulam chains

Studies of thermal transport in long-range (LR)interacting systems are currently particularly challenging. The main difficulties lie in the choice of boundary conditions and the definition of heat current when driving systems in an out-of-equilibrium state by the usual thermal reservoirs. Here, by employing a reverse type of thermal baths that can overcome such difficulties, we reveal the intrinsic features of thermal transport underlying a LR interacting Fermi-Pasta-Ulam chain. We find that under an appropriate range value of LR exponent $\sigma =2$, while a \emph{nonballistic} power-law length ($L$) divergence of thermal conductivity $\kappa$, i.e., $\kappa \sim L^{\alpha}$ still persists, its scaling exponent $\alpha \simeq 0.7$ can be much larger than the usual predictions in short-range interacting systems. The underlying mechanism is related to the system's new heat diffusion process, weaker nonintegrability and peculiar dynamics of traveling discrete breathers. Our results shed light on searching for low-dimensional materials supporting higher thermal conductivity by involving appropriate LR interactions.

which implicitly states that the heat current J is proportional to the temperature gradient ∇T with κ the thermal conductivity being constant for a bulk material, is not valid. Instead, κ shows a sublinear power-law divergence as increasing system size L, i.e., κ ∼ L α (0 < α < 1). This issue is in part motivated by the recent carbon nanotubes (CNTs) technology [4,5], and also by the desire to understand the microscopic origin of non-Fourier's thermal transport. Additionally, recently how to efficiently control heat is a fascinating topic, which has led a new emerging field of phononics [6][7][8][9][10].
Indeed, continuous efforts of theoretical studies on systems of classical coupled oscillators with nearest-neighbor (NN) interactions, such as the Fermi-Pasta-Ulam (FPU) models, and experimental investigations on 1D CNTs have confirmed such anomalous effects in a quite good level. The theoretical predictions of α are varying from 0.2 to 0.5 [11,12], which has been experimentally corroborated in 1D single-walled CNTs [4,5]. However, the same experiments performed on 1D multi-walled CNTs [4] showed that α can be even higher (α = 0.6-0.8). As the real systems are usually plagued with many effects, such as defects, isotopic disorders, and impurities that would greatly lower down the divergence, this poses a fundamental question: Is there any theoretical model free of such effects can support a high divergent exponent comparative to experimental observations?
In this Letter we show that a new class of α 0.7 that falls within α = 0.6-0.8 can achieve in a theoretical model of long-range (LR) interacting FPU chain. This result provides a more solid theoretical ground of sufficiently high thermal conductivity appearing in 1D ma-terials. It also suggests the need to extend the study of 1D thermal transport beyond the short-range interacting systems. In principle this indicates a new direction for further exploring thermal transport in 1D systems.
Before proceed, we would like to stress that the study of thermodynamics in LR interacting systems is certainly far from trivial and of its practical interest. Theoretically, LR interactions are ubiquitous in many physical systems, ranging from self-gravitating systems to nanoscale systems, and to quantum systems [13][14][15][16]. These LR interacting systems can display many peculiar features, such as ensemble inquivalence, the lack of additivity, and anomalous diffusion of energy. Practically, the present technology allows to fabricate materials with LR interactions. The Coulomb crystals [17], Ising pyrochlore magnets [18,19], and permalloy nanomagnets [20] are some of the noticeable examples.
Our model is a FPU type system of N particles including LR interactions with Hamiltonian where x i is the ith particle's displacement from its equilibrium position; p i is its momentum; |j − i| σ represents the interaction strength of the ith particle with its |j−i|th neighborings (σ denotes the range value of LR exponent). Its thermal conduction properties have been studied by the usual method [21], i.e., directly coupling two thermal reservoirs to two ends of the system [22]. However, when one performs studies like that, the feature of nonadditivity of the system will cause the central bulk particles and the two thermalized ends implicitly identified. This results that eventually the whole system may display properties that do not correspond to the original one. Due to this, it has been pointed out that one needs to distinguish between heat flows towards or from the reservoirs and those within the system [23]. View-arXiv:1906.11086v1 [cond-mat.stat-mech] 26 Jun 2019 ing this, here we use an alternative approach that called "reverse nonequilibrium molecular dynamics simulation (RNEMD)" method [24,25] to perform our study. This method can get rid of such boundary effects and wherein the heat current is naturally defined. Note that our model only considers the LR interactions in the quartic anharmonic term [21], which differs from the model in [23], but this does not violate our general conclusion [26]. In addition we are not going to include the Kac scaling factorÑ = 1 This factor was designed to restore the system's extensivity as increasing system size, but it does not help improve the system's nonadditivity [16]. It only constructs an "artificial" extensive system, but our cost for thermal transport is that the phonon's group velocity should depend on system size, which is an unwanted effect. In fact, when looking at dynamical aspects, for a correspondence between both treatments, the only difference is that time should beÑ 1/2 scaled [27].
We shall focus on the case of σ = 2 and also present the result of σ = 8 for comparison. The latter might correspond to the FPU model with NN interactions, although this equivalence requires σ → ∞. The former is of particular interest since wherein ballistic transport (α 1) [21] has been conjecture. This is attributed to the strong finite-size effects and it has been pointed out that to gain convincing results is extremely hard [23]. But anyway, this implies new physics. Another reason for σ = 2 of interest is that the system can support travelling discrete breathers (DBs) without tail under zero temperature [28]. Then one might ask: what will happen of these moving excitations at finite-temperature systems and how they would affect transport?
The RNEMD method produces temperature gradient in an unusual way. Unlike the traditional approach [22] to directly induce ∇T , it imposes the heat current by frequently exchanging the kinetic energy (or momentum) between particles. While the nonequilibrium stationary state reaches, ∇T will establish. This "reversion" makes such method an ideal candidate for studying thermal transport in the LR interacting systems.
We implement the method in such a way: First, periodic boundary conditions are used and the chain forms like a ring. The ring is then decomposed into M = 32 equaling slabs (each contains n = N/M particles). We give each slab a serial number and label the cold one slab 1 and accordingly, the hot one slab M/2 + 1. This labeling allows us to interchange the momentum of the hottest particle in slab 1 with that of the coldest particle in slab M/2 + 1 at a frequency f exc = 0.1. Such interchanges cause a redistribution of kinetic energy of the system with an energy difference: ∆E = where the subscripts h and c refer to the hottest and coldest particles whose momenta are exchanged, and the sum runs over all exchange events in time t. As a con- sequence, the relaxation of energy difference will drive two heat currents to flow from hot to cold slabs along the two semiring sides bridging them (each side has an effective length L = N/2 − n). After the stationary state eventually reaches, the long-time averaged current across each side is defined by accordingly the time averaged kinetic temperature of each slab is T k = 1 nk B nk i=n(k−1)+1 p 2 i , where k B is the Boltzmann constant (set to be unit) and the sum runs over all n particles in slab k. The heat conductivity κ can then be obtained by κ = − J /∇T according to Fourier's law, with ∇T being evaluated over the slabs between the cold and hot ones.
We start our calculations with several fully thermalized systems under T = 0.5. These systems are evolved by velocity-Verlet algorithm [29] with a time step 0.01, that guarantees energy conservation with a relative accuracy of O(10 −5 ). As for LR interacting systems the calculations for force at each time step that demands O(N 2 ) operations is extremely huge, we adopt an algorithm based on the Fast Fourier Transform [30] to accelerate our computations. With this, a transient stage of time 10 6 for system to reach the stationary state, is discarded, and the next evolution of time 10 6 is performed for average. Figures 1(a) and (b) depict two typical temperature profiles for σ = 2 and 8. In both cases a well-behaved temperature gradient is identified. The difference is that ∇T for σ = 2 is much smaller than that for σ = 8. Despite this difference, both results are obviously not the case of flat temperature profile of ballistic transport shown in integrable systems [31,32]. Therefore, here the integrable dynamics is excluded. κ ∼ L α has been observed. This is the case for σ = 2 for the all considered range of L. While for σ = 8 an asymptotic behavior only achieves for a large L as the crossover to the NN interaction model. Indeed, the best fitting gives α = 0.34 ± 0.01, which is very close to the prediction of α = 1 3 in FPU chains with NN interaction [33] and within the recent predictions of two universality classes of α [34,35]. Remarkably, for σ = 2, we have obtained an enhanced κ (at least one order of magnitude compared to σ = 8 for the same L) and a quite large α ( 0.71). We emphasize that this α just falls within the experimental observations of α = 0.6-0.8 in 1D multi-walled CNTs [4], and thus it can be expected, as L approaches macroscopic scales, an extremely high thermal conductivity can achieve, which undoubtedly suggests potential applications. From the theoretical point of view, this indicates new dynamical mechanisms for thermal transport.
We understand such mechanisms from the following three aspects. First, it is rooted in a new type of heat diffusion process. This exhibits in the propagation of heat fluctuations following a new type of density function with a new scaling. To characterize such process, we employ the spatiotemproal correlation function [36,37] Here, due to the translational invariance, the correlation depends only on the relative distance m; · represents the spatiotemporal average; l labels a coarsegrained bin's number similar to that adopted in the RNEMD method (here each bin has n = 8 particles). In the definition of Q l (t), g l (t), |j−k| σ and the sum runs over all particles within bin l), and F l (t)( F ≡ 0) correspond to the particle number, energy and pressure den-  [38] with a slowly relaxed central peak together with two side peaks moving ballistically [see Fig. 2(b)], like what we have observed in the short-range interacting models. In contrast, for σ = 2, ρ Q (m, t) indeed shows a new shape [see Fig. 2(a)]. Remarkably, the central peak in a long time now turns to a platform, indicating a much faster relaxation. We have further checked that the short time's ρ Q (m, t) behaviors similarly to the sort of ballistic transport [39], and thus our long-time scaled dynamics by excluding the Kac scaling does help reveal the real regime of transport.
We figure out this new transport regime by studying This scaling formula is based on the Lévy walk theory [38]. As shown, it is only well satisfied for σ = 8 (especially the central parts), and for σ = 2 [see Fig. 2(c)], to see an excellent collapse requires a much longer time. However, this does not mind our using such a formula. In fact, another relation α = 2−γ [38] based on the same theory connecting γ to α just gives an excellent estimation α 2 − 1.29 = 0.71, in agreement with our above thermal conduction calculation.
Second, the new exponent is related to the system's weaker nonintegrability. The nonintegrability is featured by the maximal Lyapunov exponent λ max (> 0) (see Fig. 3), obtained from the standard Benettin-Galgani-Strelcyn technique [40]. As comparison, another completely integrable Toda chain [32,41] with λ max = 0 is also demonstrated in Fig. 3. As shown, λ max for σ = 2 just lies in between the results for σ = 8 and Toda chain, and this seems not changed as further increasing system size. It thus indicates a weaker nonintegrability of the system compared to the short-range interacting FPU model. Indeed, a nonmonotonic variation of λ max on σ (see the inset) also confirms this [21], but the integrable dynamics is certainly ruled out. Therefore, a weaker nonintegrability seems to provide a mechanism to raise the divergent exponent. We finally conjecture that this weaker nonintegrability can make the system support a new type of excitations, i.e., the travelling DBs [28], and it is these moving DBs together with their relatively weak interactions leading to the high divergence. To verify this conjecture is of interest but usually greatly challenging since it is hard to catch out these moving excitations at equilibrium states due to their mobility. Viewing this we choose to study the spatiotemporal evolutions of local energy densities E i (t) for two-time scales (t = 100 and 1000) to visualize this process (see Fig. 4). The evolutions are obtained by considering a short chain with N = 200 for facilitating the computation. The chain is first thermalized to T = 0.5, then the thermal baths are removed and the results are recorded and displayed by a suitable time step [∆t = 1(10) for t = 100(1000)]. As indicated in the upper panel of Fig. 4, both σ = 2 and 8's short-time scaled dynamics exhibit transport similar to the ballistic regime. This explains the ballistic scaling that we have observed in a short time. In contrast, for a relatively long-time scale (see the lower panel of Fig. 4), the ballistic regime for σ = 8 disappears, suggesting strong interactions between heat carriers. But this is apparently not the case for σ = 2: the signature of the localized excitations is still recognized, but probably due to their weak interactions, their identification now becomes weaker. Such evidence is in good accord with our conjecture.
To summarize, we have presented a striking result that a high length-divergent exponent α 0.71 of thermal conductivity, that beyond the existing theoretical predictions and falls within the experimental observations of α = 0.6-0.8, can occur in a theoretical model of LR interacting FPU chain under an appropriate range value of σ = 2. This finding is of fundamental importance as it provides a promising ground for the highly sufficient thermal conductivity observed in 1D materials, thus pointing towards new manipulations of heat in practice [42,43]. It also opens up a new avenue for further exploring thermal transport since such a new type of divergence should be supported by new mechanisms. We have revealed that the new exponent is related to the system's more rapid heat propagation, weaker nonintegrable chaotic dynamics, and also the new microscopic dynamics of travelling DBs. The heat propagation shows a new type of density with its scaling exponent γ 1.29 that can be connected to α by the formula α = 2 − γ from Lévy walk theory. This suggests that the Lévy walk model, with appropriate variations, can still be useful for understanding the transport in LR interacting systems. The system's weaker nonintegrability has been confirmed and it indicates that although the chaotic dynamics is not an ingredient necessary for the invalidity of Fourier's law [22], the change of its strength does sufficiently influence the system's thermal conduction, thus paving new way to use nonlinearity to control transport. More-interestingly, the weaker nonintegrability can result in weaker dynamics of interactions of travelling DBs at thermal equilibrium and this seems responsible for the high divergence. This last evidence, together with the above, would undoubtedly encourage more studies of transport from the underlying microscopic dynamics.