Universal heat superdiffusion in one-dimensional nonequilibrium lattices

PHYSICAL REVIEW RESEARCH 6, 043211 (2024) Universal heat superdiffusion in one-dimensional nonequilibrium lattices Yuanyang Ren ,1,2Yang Wang ,3Kai Wu ,1and David Cubero 4,* 1State Key Laboratory of Electrical Insulation and Power Equipment, Xi’an Jiaotong University, No. 28 Xianning West Road, Xi’an 710049 Shaanxi, China 2Electric Power Research Institute of Yunnan Power Grid Company Ltd., China Southern Power Grid, Kunming 650217, China 3School of Electronics and Information, Xi’an Polytechnic University, Xi’an 710048, China 4Departamento de Física Aplicada I, EPS, Universidad de Sevilla, Calle Virgen de África 7, 41011 Sevilla, Spain (Received 22 January 2024; accepted 30 October 2024; published 26 November 2024) We propose a generalized Fourier’s law, based on recent mode-coupling calculations, for one-dimensional nonlinear lattices in nonequilibrium simulations, where the boundaries are kept at different temperatures. Though not useful for the typical textbook numerical examples, we are able to demonstrate numerically the theory for an anharmonic lattice with force field parameters similar to those of realistic low-dimensional structures such as covalently bonded molecular chains or carbon nanotubes, specifically showing that the asymptotic heat conductivity is described by the mode-coupling theory in good approximation, without any adjustable parameter. DOI: 10.1103/PhysRevResearch.6.043211 I. INTRODUCTION Standard Fourier’s law, Jx=−κT/L, relating the heat flux Jxwith the (small) temperature gradient in a system of size L, is only guaranteed for materials where heat transport is diffusive. In pure harmonic systems, it is well known [1] that transport is ballistic, and thus the thermal conductivity κdiverges as the system size L→∞. It came as a surprise when it was shown, in 1955 by Fermi-Pasta-Ulem in a celebrated numerical study [2], that simple one-dimensional (1D) anharmonic lattices are also not able to drive the system to thermal equilibrium, for this kind of divergent behavior was only expected from harmonic lattices, where phonon scattering is absent. The Fermi-Pasta-Ulem (FPU) lattice was one of the first examples [3] of a nonlinear lattice displaying superdiffusive behavior, defined as a thermal conductivity κthat diverges as κ∼ALα, where Ais a nonuniversal factor, and 0<α<1(α=0 corresponds to normal materials satisfying standard Fourier’s law and α=1 to ballistic transport). Recent mode-coupling calculations [4,5], based on local hydrodynamics, have predicted for 1D anharmonic lattices the exponents α=1/3 and α=1/2, with explicit expressions for nonuniversal coefficients for the case α=1/3, and placing it under the class of the Kardar-Parisi-Zhang (KPZ) equation. They have been validated numerically by looking at the spatiotemporal correlation functions [6–8]—thus boundaryless simulations. *Contact author: [email protected] Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Nonequilibrium simulations, where a temperature gradient is maintained by using two thermostats at the boundaries, is a natural setup to study heat transport. They were the first choice [3] to show superdiffusive behavior and can be directly applied to simulate realistically low-dimensional structures such as carbon [9,10], or boron nitride nanotubes [11], where superdiffusive behavior has been experimentally demonstrated [12]. The validity of the 1D mode-coupling calculations in this kind of bounded setup has so far remained elusive. Here we demonstrate that they are able to quantitatively predict the asymptotic anomalous heat conductivity in nonequilibrium simulations of 1D lattices with realistic force field parameters, typical in the modeling of covalently bonded low-dimensional structures. A generalized Fourier’s law is derived and validated numerically in the simulations for this purpose. It must be stressed that this kind of success is only observed in the stated force field range, as, outside of that scope, as discussed below, the generalized Fourier’s law and the associated mode-coupling predictions are not seen to be validated in the nonequilibrium setup. II. MODEL We start by considering the 1D lattice as a model of a molecular chain in three-dimensional (3D) space, made of carbon atoms in a straight line, with each carbon separated from its neighbor by an equilibrium distance of a=1.286Å, and with the many-body Tersoff force field [13]. The 1D lattice consists of Nparticles, with position xj=−L/2+aj+qj, j=1,...,N, momentum pj,massm=12 g/mol, and governed by the Hamiltonian H= N  j=1p2 j 2m+V(qj+1−qj),(1) 2643-1564/2024/6(4)/043211(6) 043211-1 Published by the American Physical Society REN, WANG, WU, AND CUBERO PHYSICAL REVIEW RESEARCH 6, 043211 (2024) with the prototypical FPU potential V(x)=u2x2/2+ u3x3/3! +u4x4/4!. The potential parameters, u2=72.54 eV/Å2,u3=−376.9eV/Å3,and u4=1241.8eV/Å4,(2) were obtained by fitting the molecular chain potential energy when one-half of the chain was displaced a distance xfrom its equilibrium position. The nonequilibrium molecular dynamics (NEMD) setup was formed by adding Nbath =400 particles at each end of the chain [14]. These particles, which prolong the chain in both directions, are of the same kind as in the chain bulk, but subject to Langevin thermostats which keep the temperature at the left and right ends to Tc=280 and Th=320 K, respectively. The dissipative forces associated with the Langevin thermostats have a friction coefficient 1/τbath, which was chosen as τbath =0.5ps. In addition, to avoid the chain translational movement, we added, at the chain’s ends, two fixed particles with no equilibrium deviations, q−Nbath =qN+Nbath+1=0. The average length of the part of the chain that is not thermostatted is then given by L=a(N−1). Naturally, for one-dimensional systems at rest, the flux Jxis defined as the energy per unit of time that crosses one particle, which can be easily computed in the simulations as the energy rate transferred to the left j=1 or right j=Nparticle. III. FIRST RESULTS AND THE VELOCITY EXCHANGE MECHANISM Figure 1shows the simulation results (symbols), with the solid lines displaying fits to the data assuming the power law JxN∝Nα. The 1D Hamiltonian model (blue diamonds) is seen to exhibit a purely ballistic behavior (α=1) in the size range that is shown. Clearly, the harmonic term dominates at this range, though, strikingly, it does not in the carbon chain in 3D space (black crosses), which does not show a clear power-law behavior in the figure. Here, phonon scattering, responsible for that deviation, is mediated by the two additional dimensions, which provides a local mechanism for thermalization. In order to capture this scattering mechanism, we have added, to the 1D lattices, stochastic collisions between neighboring particles, in the form of random exchanges of the two momenta at neighboring sites. Figure 1shows a good agreement between the 1D simulation with random velocity exchanges at a rate γ=2.2ps −1(black diamonds) and the 3D results (black crosses). It could be tempting to conclude that the two additional dimensions just provide a simple randomization mechanism, mimicked by the random velocity exchanges, so that anharmonic effects are not needed to account for the anomalous behavior observed in the unperturbed lattice. On the contrary, triangle-up points in Fig. 1show a distinctly different power law, with α=1/2, for perturbed lattices without the cubic term (u3=0), one that is purely harmonic (u4=0) and another in which the quartic term is increased fourfold. The latter displays identical behavior to that of the perturbed harmonic lattice. 1000 10000 N 100 102 104 Jx N (eV/ps) Jx N ∝ Ν0.3333 Jx N ∝ Ν Jx N ∝ Ν0.5 FIG. 1. Anomalous conductivity for several chains in nonequilibrium simulations. The (black) crosses correspond to a carbon straight chain in 3D space, while all other curves are from 1D lattices. The (blue) diamonds are from the anharmonic lattice with the potential parameters fitted to the carbon chain, given by Eq. (2). The (black) diamonds are from the same anharmonic lattice, but with random velocity exchanges between neighboring particles at a rate γ=2.2ps −1. The triangles also correspond to lattices with the same rate of velocity exchanges, but in the case of the (magenta) filled triangle-up points, the potential is harmonic (u3=u4=0); in the (red) empty triangle-up points, the potential has no cubic term and the quartic term is amplified (u3=0, u4is increased fourfold); and in the (orange) filled triangle-down points, the potential is purely quartic (being amplified eightfold). The solid lines are fits to the simulation data assuming the asymptotic behavior indicated. This α=1/2 exponent is not unexpected, for it has been rigorously proven for harmonic lattices under stochastic collisions that conserve linear momentum [15,16]. Furthermore, mode-coupling calculations [4,5] also predict an exponent α=1/2 in unperturbed anharmonic lattices when the potential is symmetric (u3=0 in our system) and the pressure is zero [5]. In other cases, the mode-coupling prediction for 1D systems is an asymptotic power law with exponent α=1/3. In Fig. 1, this behavior is observed for the perturbed anharmonic lattice without quadratic and cubic potential terms (triangle-down points). In the size range shown in Fig. 1, it is not clear if both the 3D carbon chain and its 1D model have reached an asymptotic power law with an intermediate exponent or, rather, are in an intermediate transient region, having not yet reached the final asymptotic state predicted by the mode-coupling theory for N→∞. The prediction α=1/3 relies on local hydrodynamics and, therefore, it should also be of application to the 3D carbon chain. Larger simulation sizes should then reveal it, after the shorter lived modes, which are blurring the behavior for the size range shown in Fig. 1, die out. But observing the universal exponent 1/3 for larger sizes is not enough to validate the mode-coupling explanation because the nonuniversal prefactor, which can be predicted by the modecoupling calculation, should also be observed. Indeed, a few considerations are in place before comparing with the modecoupling prediction. 043211-2 UNIVERSAL HEAT SUPERDIFFUSION IN ONE- … PHYSICAL REVIEW RESEARCH 6, 043211 (2024) IV. THE NEMD SOLUTION In the NEMD setup, the temperature profile T(x)= m˙q2 j/2, with x=−L/2+aj, reaches a steady state, which is smooth in the relevant interval −L/2<x⩽L/2. Due to the boundary conditions T(−L/2) =Tcand T(L/2) =Th, the steady profile is obviously not space periodic. However, the temperature gradient dT/dx = [T(x+a)−T(a)]/awas observed to be periodic at the boundaries, dT/dx(−L/2+)=dT/dx(L/2−), in all our simulations, which allows us to consider the gradient as a periodic function, with period L. This fact, after going into Fourier space, Jk=L/2 −L/2dxexp(ikx)Jx,T k= L/2 −L/2dxexp(ikx)dT/dx,k=2πn/L, with ninteger, allows us to write a generalized anomalous definition of the thermal conductivity, Jk=−κkT k.(3) However, for this definition to be as useful as in diffusive materials—where κkis a constant—it is desirable that the transport coefficient κkdepends only on the chain properties and k, and not on the specific details of the boundaries. This is expected for small temperature gradients if local hydrodynamics is at play. Linearization of the hydrodynamics equations to leading order in k, taking Eq. (3) into account, yields (see Appendix A) κk=λk cp k2,(4) where cpis the specific heat per particle at constant pressure (see Appendix A) and λkis the eigenvalue associated [17]to the damping of the thermal mode. Mode-coupling calculations [4,5] yield explicit expressions for λkleading to κk=Aα|k|−α in the asymptotic limit L→∞, where Aαis a nonuniversal coefficient, and the exponent can take the value α=1/3 (general case) or 1/2 [particular case of symmetric potential V(x) and zero pressure]. At the steady state, energy conservation dictates the flux Jxto be a constant throughout the interval −L/2<x<L/2, which implies that Jk=0 for all k= 0, and Jx=J0 L=−κ0 T(L/2) −T(−L/2) L.(5) Assuming a nonzero conductivity, we arrive to T k=0for k= 0, which implies dT dx (x)=constant in −L/2<x<L/2.(6) Note that the linear gradient prediction of Eq. (6) is a direct consequence of the generalized Fourier law (3). In all our NEMD simulations, the temperature profile is observed to be linear in good approximation, as shown in Fig. 2. There are still small visible deviations from linearity, which are more pronounced near the boundaries. This is not unlike in heat diffusive materials, which in NEMD setups also shows temperature slips at the boundaries, i.e., the so-called Kapitza or boundary resistances. The failure of Fourier law at these regions is sometimes attributed to the failure of the hydrodynamics description because of the mean free path not being 280 290 300 310 320 T (K) -0.4 -0.2 0 0.2 0.4 x / L 0.995 1 1.005 temp. ratio N=10000 N=20000 N=40000 N=60000 FIG. 2. Typical temperature profile in the explored region, indicating linearity in good approximation. The data correspond to a FPU 1D lattice perturbed with random velocity exchanges at rate γ=8ps −1. Top panel: The results for N=40 000 (diamonds) and its linear fit (red solid line)—which yielded a correlation coefficient larger than 0.999. Bottom panel: The ratio of the local temperature to its linear fit value plotted for several system sizes, highlighting the deviations from linearity—note they are less than 0.5%, except very near the boundaries. sufficiently small compared to the hydrodynamics variation length. However, while in diffusive materials the size of the boundary layer is approximately independent of the system size, and thus negligible as L→∞,Fig.2reveals that this is not the case for superdiffusive materials, with the boundary layer scaling instead with the size L. In fact, the bottom panel of Fig. 2suggests that the temperature profile takes a scaling form independent of the system size. Note, however, these deviations from linearity are of the order of 0.5% or less. Let us remark that in the parameter range of the simulations, which is taken after realistic models of covalent bonded low-dimensional structures at room temperature, the pressure is not zero, but negligible. This conclusion is arrived at after noticing that the mode-coupling prediction is almost the same (within 0.6%) if we consider p=0 in the calculation. This is not the case for the typical parameter range considered in the literature [18,19], where all coefficients u2,u3, and u4 are taken of the order of unity. There, the mode-coupling calculation indicates a non-negligible pressure and larger appreciable deviations from the linear temperature profile are observed. Further increasing the pressure—by reducing the distance between the outermost (fixed) particles—produces a more distorted nonlinear profile, in agreement with previous observations [20], suggesting unavoidable boundary effects. It can be readily checked that for 0 <α<1, the Fourier transform gk=L/2 −L/2dxexp(ikx)g(x) of the function g(x)=|x|α−1 2(α)cos(πα/2) (7) yields gk=|k|−αin the asymptotic limit L→∞. Then, by taking κk=Aαgkfor a finite system of size L,wearriveat κ0=AαLα 2α(α)αcos(πα/2).(8) 043211-3 REN, WANG, WU, AND CUBERO PHYSICAL REVIEW RESEARCH 6, 043211 (2024) The expression (8) is also valid in the limiting cases α=0 (normal diffusion), where g(x)=δ(x) and κ0=A0; and α= 1 (ballistic limit), where κ0=ˆ A1L, with the constant ˆ A1= limα→1Aα/[2 cos(πα/2)]. For a harmonic chain (u3=u4= 0) with Nbath =1, and without the fixed particles at the ends, it has been found [1] analytically ˆ A1=u2kBτbath 21+μ 2(1 −1+4/μ),(9) where μ=mu2τ2 bath, i.e., fully dependent on boundary details such as τbath. In the general superdiffusive case α=1/3, mode coupling predicts [4,5]λk=λh|k|5/3, where the heat eigenvalue λh depends on the potential parameters, i.e., temperature and pressure, thus determining the nonuniversal coefficient as A1/3=λhcp.(10) Note than unlike the ballistic case, the mode-coupling prediction of Eq. (10) is fairly independent of the boundary details, pretty much like the diffusive case α=0. It thus remains to be seen if it can be actually observed in a practical application—a central issue we focus on in the remainder of the paper. Computer limitations prevent us from simulating the 1D lattice γ=2.2ps −1of Fig. 1with a large enough size so that the expected asymptotic limit α=1/3 is observed. We can overcome this limitation by observing that the mode-coupling prediction is the same regardless of the specific value of γ, since the random velocity exchanges do not change the local thermodynamics—which is the base of the mode-coupling calculations. On the other hand, increasing γcan change other transient modes. Intuitively, this increase should introduce more phonon scattering, thus lowering the flux and, with luck, revealing the final mode-coupling asymptotic behavior in a manageable size range. That is precisely what is observed, as illustrated by Fig. 3. Increasing the rate γproduces a curve that is increasingly converging towards a clear power law with exponent α=1/3, with a prefactor that is quite near the full mode-coupling prediction, i.e., the solid (red) line in Fig. 3. Note that the nonuniversal factor is calculated ab initio using Eq. (10), that is, with no fitting parameters. As a second example, we consider a FPU 1D lattice with the following parameters: m=12 g/mol, u2=9.4eV/Å2,u3=−44.5eV/Å3,and u4=137.5eV/Å4,(11) and a particle separation a=1.25 Å. This 1D lattice is associated to carbon nanotubes (CNTs), which can be simulated with 3D molecular dynamics codes, such as LAMMPS [21], using Tersoff potentials parametrized to recover the phonon spectrum of carbon nanostructures [22]. In order to speed up the universal superdiffusive behavior expected for very large sizes, such as we have demonstrated above for the first example, we have implemented in the code LAMMPS the same velocity exchange mechanism for the CNTs as in the 1D lattice. The CNT(5,5) can be viewed as 10 linear (zigzag) carbon chains along the tube symmetry axis, covalently joined in the perpendicular directions to form the nanotube. For the comparison, we regard the 1D lattice as a model of one of 0246810 γ (ps-1) 101 1000 10000 N 100 101 Jx N / ΔT (eV/ps K) γ=0.9 ps-1 γ=2.2 ps-1 γ=3.5 ps-1 γ=6 ps-1 γ=8 ps-1 ab initio FIG. 3. Speeding up the universal superdiffusive behavior. The velocity exchange rate γis demonstrated as a convenient control parameter to reveal the asymptotic universal 1/3 tail. The diamonds correspond to the anharmonic 1D lattice with various values of velocity exchanges rates. The solid (red) line is the mode-coupling 1/3 prediction, with the nonuniversal factor calculated using Eq. (10) with no fitting parameters. Inset: Convergence towards a plateau value for N=40 000 as the rate γis increased, which is in very good agreement with the mode-coupling prediction. those chains and approximate the nanotube’s flux as 10 times that of a single chain, ignoring the extra scattering introduced by the lateral bonding. As shown in Fig. 4, the 1D lattice (11) captures very well the behavior of the carbon nanotube as a function of the velocity exchanges rate γ, providing a very similar thermal conductivity per chain, except for very low γ, which is expected because for γ=0 the 1D lattice has 01234 γ (ps-1) 0 2 4 6 Jx N / ΔT (eV/ps K) FIG. 4. A 1D chain and a carbon nanotube. Thermal conductivity as a function of the exchange rate γfor a fixed system size L=1µm (N=7990). The crosses are the results of the carbon nanotube CNT(5,5), with the flux divided by 10 to account for the 10 linear chains, covalently joined transversally, which forms the nanotube. The diamonds correspond to the FPU 1D chain defined in Eq. (11). The solid red line is the ab initio prediction of the 1D chain given by the mode coupling in Eq. (10). 043211-4 UNIVERSAL HEAT SUPERDIFFUSION IN ONE- … PHYSICAL REVIEW RESEARCH 6, 043211 (2024) no inherent thermalization mechanism, in contrast to the 3D carbon nanotube. Note, however, that the 1D lattice is not expected to mimic all the heat transport properties of the nanotube; it should be regarded as a 1D lattice with the same mode-coupling behavior as the CNT. Finally, let us notice that, again, the mode-coupling prediction of Eq. (10) is in very good agreement with the simulation data. V. CONCLUSIONS In conclusion, by using a simple thermalization mechanism in 1D lattices to speed up convergence into a manageable system size, we have shown that a generalized Fourier law, based on mode-coupling theory, is able to describe in good approximation superdiffusive heat transport in a standard nonequilibrium setup, without any adjustable parameter. Unlike normal diffusive systems, the validity of the approximation is observed to be limited to a range of force field parameters typical of covalently bonded molecular structures. However, this is within the scope of many practical applications of low-dimensional structures, in which heat transport is of crucial importance. ACKNOWLEDGMENTS Financial support from the Ministerio de Ciencia e Innovación of Spain, Grant No. PID2019-105316GB-I00, is acknowledged (D.C.), and from the Yunnan Lancang-Mekong National Joint Laboratory for Electric Power Technology, Grant No. 202203AP140001 (Y.R.). APPENDIX A: DERIVATION OF EQUATION (4) Our starting point is the detailed calculations by Spohn in Ref. [5]. Here, we follow his notation so that we can express the temperature gradient in terms of β=1/(kBT)as ∂T ∂x=− 1 kBβ2 ∂β ∂x,(A1) and write the variations of βin terms of elongation and internal energy per particle e(both defined in Eq. (2.10) of [5]) as ∂β ∂x=∂β∂ ∂x+∂eβ∂e ∂x.(A2) We are now ready to lift the linearized Euler equations (Eqs. (3.6) and (3.7) of [5]) to include the generalized thermal conductivity, defined by Eq. (3) of the main paper, yielding the matrix A(Eq. (3.7) of [5]) replaced, after going into Fourier space, by A=⎛ ⎜ ⎝ 0−10 ∂p0∂ep −k2κk kBβ2∂βpk2κk kBβ2∂eβ ⎞ ⎟ ⎠.(A3) The eigenvalue of (A3) associated with the heat mode is given, in the lowest order of k,by λh=k2κk kBβ2 ∂β∂ep−∂eβ∂p ∂p−p∂ep.(A4) 0500 1000 1500 Nbath 20000 25000 30000 35000 40000 Jx N / ΔT FIG. 5. Flux dependence with the size of the thermostatted regions, Nbath, for the anharmonic 1D lattice of size N=10 000 perturbed with a velocity exchange rate of γ=8ps −1.Thered solid line is the mode-coupling prediction. The error bars were estimated from the difference between the fluxes at the left and right boundaries. Next, we consider the map (, e)→(p,β)→(, e), and the four equations that follow from considering ∂/∂ =1, ∂/∂e=0, ∂e/∂e=1, ∂e/∂ =0, for example, 1=∂e/∂e=∂βe∂eβ+∂pe∂ep,(A5) and so on. From these equations, we find ∂βe=− ∂p ∂β∂ep−∂eβ∂p,(A6) ∂β=∂ep ∂β∂ep−∂eβ∂p.(A7) -0.4 -0.2 0 0.2 0.4 x / L 280 290 300 310 320 T (K) N=3000 y = (284.79±0.06) + (27.7±0.1) * x y = (281.3±0.2) + (56±3) * x y = (260.6±1.9)+ (54±2) * x γ=2.2 ps-1 FIG. 6. Typical temperature profile illustrating the periodicity of the temperature profile at the boundaries in the NEMD setup. The system is the first 1D lattice example, with N=3000, under random velocity exchanges with rate γ=2.2ps −1. The solid blue line is the linear fit in the bulk. The solid red line and the solid orange lines at the ends are a linear fit in the range that is plotted. The calculated slopes at the ends are within the uncertainties. 043211-5 REN, WANG, WU, AND CUBERO PHYSICAL REVIEW RESEARCH 6, 043211 (2024) Using (A6) and (A7), the specific heat per particle at constant pressure can be written as cp=−kBβ2(∂βe+p∂β)=kBβ2(∂p−p∂ep) ∂β∂ep−∂eβ∂p.(A8) Combining Eqs. (A8) and (A4) yields Eq. (4)ofthemain manuscript. APPENDIX B: THE HEAT EIGENVALUE In the general case α=1/3, the heat eigenvalue is given by Eq. (4.12) of [5], though the expression λh=λ−2/3 sc−1/3G0 112ac,(B1) given in Ref. [7], is more appropriate with our definition of Fourier transform, with the sound mode eigenvalue given by λs=2√2G1 11.(B2) Note, however, that the numeric value of the coefficient acis given by [8] ‘ac=2√3(1/3)∞ −∞ dxfKPZ(x)2=3.6175.... (B3) The values of the sound speed c, i.e., |G1 11|,G0 11, and cp,are easily obtained from the simulation pressure and average temperature (Tc+Th)/2 following the formulas given in [5]. Note that the fixed particles at the boundaries in the nonequilibrium setup make the average elongation =0. APPENDIX C: INSENSITIVITY TO Nbath The system in the nonequilibrium setup is sensitive to the size of the thermostatted regions at the boundaries, determined by Nbath, only if they are small. Figure 5shows the results for a typical system. For large enough Nbath, the simulation results become independent of the value of Nbath. APPENDIX D: PERIODICITY OF THE TEMPERATURE GRADIENT As reported in the manuscript, we have observed in all our simulations that the temperature gradient is periodic at the boundaries, i.e., it shows the same value, within the uncertainties, at both ends, as illustrated in Fig. 6. [1] Z. Rieder, J. L. Lebowitz, and E. Lieb, Properties of a harmonic crystal in a stationary nonequilibrium state, J. Math. Phys. 8, 1073 (1967). [2] E. Fermi, J. Pasta, and S. Ulam, Studies of nonlinear problems, Vol. II Los Alamos Report number LA–1940 (1955). [3] S. Lepri, R. Livi, and A. Politi, Heat conduction in chains of nonlinear oscillators, Phys. Rev. Lett. 78, 1896 (1997). [4] Henk van Beijeren, Exact results for anomalous transport in one-dimensional Hamiltonian systems, Phys. Rev. Lett. 108, 180601 (2012). [5] Herbert Spohn, Nonlinear fluctuating hydrodynamics for anharmonic chains, J. Stat. Phys. 154, 1191 (2014). [6] C.B. Mendl and H. Spohn, Dynamic correlators of Fermi-Pasta- Ulam chains and nonlinear fluctuating hydrodynamics, Phys. Rev. Lett. 111, 230601 (2013). [7] S. G. Das, A. Dhar, K. Saito, C.B. Mendl, and H. Spohn, Numerical test of hydrodynamic fluctuation theory in the Fermi-Pasta-Ulam chain, Phys.Rev.E90, 012124 (2014). [8] J. Schmidt, G.M. Schütz, and H. van Beijeren, A lattice gas model for generic one-dimensional Hamiltonian systems, J. Stat. Phys. 183, 8 (2021). [9] C. W. Chang, D. Okawa, H. Garcia, A. Majumdar, and A. Zettl, Breakdown of Fourier’s law in nanotube thermal conductors, Phys. Rev. Lett. 101, 075903 (2008). [10] U. Ray and D. T. Limmer, Heat current fluctuations and anomalous transport in low-dimensional carbon lattices, Phys. Rev. B 100, 241409(R) (2019). [11] Y. Ren, Y. Wu, B. Xiao, K. Wu, and D. Cubero, Heat transport and surface functionalization in nanocomposites of boron nitride nanotubes and polyethylene, Phys. Chem. Chem. Phys. 23, 9604 (2021). [12] V. Lee, C.-H. Wu, Z.-X. Lou, W.-L. Lee, and C.-W. Chang, Divergent and ultrahigh thermal conductivity in millimeter-long nanotubes, Phys. Rev. Lett. 118, 135901 (2017). [13] A. Kinaci, J. B. Haskins, C. Sevik, and T. Çagin, Thermal conductivity of BN-C nanostructures, Phys.Rev.B86, 115410 (2012). [14] The specific value of Nbath =400 used in the simulations was large enough to guarantee insensitivity to this parameter; see Appendix C. [15] G. Basile, C. Bernardin, and S. Olla, Momentum conserving model with anomalous thermal conductivity in low-dimensional systems, Phys. Rev. Lett. 96, 204303 (2006). [16] C. Bernardin and G. Stoltz, Anomalous diffusion for a class of systems with two conserved quantities, Nonlinearity 25, 1099 (2012). [17] L. Kadanoff and P. Martin, Hydrodynamic equations and correlation functions, Ann. Phys. 24, 419 (1963). [18] A. Dhar, Heat transport in low-dimensional systems, Adv. Phys. 57, 457 (2008). [19] G. Benenti, S. Lepri, and R. Livi, Anomalous heat transport in classical many-body systems: Overview and perspectives, Front. Phys. 8, 292 (2020). [20] T. Mai, A. Dhar, and O. Narayan, Equilibration and universal heat conduction in Fermi-Pasta-Ulam chains, Phys.Rev.Lett. 98, 184301 (2007). [21] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 117, 1 (1995). [22] L. Lindsay and D. A. Broido, Optimized Tersoff and Brenner empirical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene, Phys. Rev. B81, 205441 (2010). 043211-6