{'fraction_non_alphanumeric': 0.03627107858733694, 'fraction_numerical': 0.016179399252895913, 'mean_word_length': 4.688945498424617, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Anomalous Heat Conduction and Anomalous Diffusion in Low Dimensional Nanoscale Systems

Thermal transport is an important energy transfer process in nature. Phonon is the major energy carrier for heat in semiconductor and dielectric materials. In analogy to Ohm's law for electrical conductivity, Fourier's law is a fundamental rule of heat transfer in solids. It states that the thermal conductivity is independent of sample scale and geometry. Although Fourier's law has received great success in describing macroscopic thermal transport in the past two hundreds years, its validity in low dimensional systems is still an open question. Here we give a brief review of the recent developments in experimental, theoretical and numerical studies of heat transport in low dimensional systems, include lattice models, nanowires, nanotubes and graphenes. We will demonstrate that the phonon transports in low dimensional systems super-diffusively, which leads to a size dependent thermal conductivity. In other words, Fourier's law is breakdown in low dimensional structures.


INTRODUCTION
The conduction of heat is one fundamental energy transport mechanisms in nature. In addition to electron and photon, phonon-the heat pulse through lattice, also carries and processes energy, which has broad applications for heat control/management in the real world. Actually, for non-metallic materials, phonon is the most dominated heat energy carrier. Its contribution to heat conduction is much larger than those from electrons and photons. Traditionally, the phenomenon of thermal transport is believed to follow the Fourier's law of heat conduction where J is the heat flux in the system, ∇T is the gradient of temperature. The conductivity κ is a geometry-independent coefficient which mainly depends on the composition and structure of the material and the temperature. This empirical law has received great success in describing macroscopic thermal transport and is widely accepted as a general truth.
In the past decades, low dimensional nano scale systems have been extensively studied due to their promising potential applications for future electronic, optoelectronic, and phononic/thermal devices. In addition to the electrical and optical properties, the thermal properties of nanoscale materials are also important but less studied compared with electronic and optical properties. The thermal properties in nanostructures differ significantly from their bulk counterpart because the phonon characteristic lengths are comparable to the characteristic length of the nanostructures. However, it is still an open and much debated question whether the Fourier's law is valid in low dimensional systems. A rigorous proof for this empirical law from microscopic Hamiltonian dynamics is still absent even though it is already two-hundred-year old. Here "invalid" means the thermal conductivity κ will depend on the system size.
Such size dependent thermal conductivity has been observed in theoretical models, such as the harmonic chains [1,2], the FPU-β model [3,4] and the hard point gas model [5][6][7]. The discovery of this anomalous behavior in general low dimensional models has then inspired enormous research studies. It is found that in these models, the thermal conductivity typically diverges with system size: for 1D models and for 2D models in thermodynamic limit, which predicts extraordinarily high thermal conductivity for low dimensional systems with finite but large sizes. Therefore, not only it is a fundamental demand for the development of statistical physics to understand normal and anomalous heat conduction in low dimensional systems, but it is also of great interest from the application point of view, since the achievement of modern nano fabrication technology allows one to access and utilize one dimensional (1D) and two dimensional (2D) structures with sizes in the range of few nanometers up to few hundred nanometers.
In recent years, much effort has been devoted to study the anomalous heat conduction in low dimensional nanostructures by numerical simulations such as [8] and [9]. These studies provide evidence that low dimensional nanostructures are very promising platforms to verify fundamental thermal transport theories. Inspired by these theoretical and numerical studies, it is exciting that length dependent thermal conductivity has been observed in carbon nanotubes and boron-nitride nanotubes [10]. Unfortunately, an experimental investigation of the quantitative size dependent thermal conductivity in nanoscale is exceptionally challenging. The difficulty lies primarily with the technique associated with measuring temperature distribution and heat flux of low dimensional systems, with the added complexity of accurately varying the size of the investigated object. Due to these challenges in experimental measurement, the experimental data need to be complemented by theoretical study on a microscopic footing.
Generally, heat transport can be studied in the aspects of two different yet deeply related phenomena. One is to study the heat flux in a non-equilibrium steady state to get the heat conductivity as we normally do in simulation or in experiments, in which temperature biases are applied to the system. The other approach is to study the heat diffusion in a nonequilibrium transient process in which the system is firstly excited away from equilibrium and then relax to equilibrium. Both approach should reveal some properties of heat conduction of the system. Moreover, the connection between thermal conductivity and heat diffusion, especially for the anomalous cases, is still an open question and attract wide research interests.
As heat conduction in low dimensional systems is a fast growing area, combination of experimental, theoretical and numerical investigations are indispensable to speed up its development. In this articles, we would like to give a review on the recent development of heat transport in low dimensional systems, from both experimental and theoretical point of view. The rest of the article is organized as follows: section introduces the experiments on quasi-1D nanostructures and 2D graphene that display the anomalous size dependent thermal conductivity. Section discusses the numerical studies on thermal conduction of nanotubes, nanowires and graphenes, which would shed light on understanding the phenomenon of size dependent thermal conductivity from a theoretical point of view. Section is devoted to the general theories on anomalous heat conduction and anomalous energy diffusion in low dimensional systems. Finally, in section and , we present a few conclusions and a short outlook.
There exists a large amount of literatures and studies on different aspects of thermal property in low dimensional systems. For a comprehensive review on thermal property of nano materials, please refer to these articles Refs. [11][12][13]. There are also a few reviews talk about anomalous heat transport in low dimensional systems from the viewpoint of fundamental statistical physics, such as Ref. [14] and [15]. Due to the limit of space, we only address the most fundamental physical issue, i.e., the violation of Fourier's law in this article.

Anomalous Heat Conduction in Quasi 1D Nanostructures
When the dimensions of materials shrink to the nanoscale, thermal transport properties can be very different from those of their counterpart bulk materials. Recent experimental in-vestigations in (quasi) 1D structures have revealed two opposite intriguing properties: unusual high thermal conductivities and significantly suppressed thermal conductivities. The unusual high thermal conductivities are attributed to the unique crystalline structures of the 1D materials, such as carbon nanotubes [16], and stretched polymer nanofibres [17], whereas the significantly suppressed thermal conductivities are due to an increased phonon-boundary scattering, which has been observed in various nanowires, such as Si [18][19][20], Si/SiGe superlattice [21], ZnO [22], Bi [23], etc. More recently, it has been observed that the thermal conductivity can be further suppressed by the coherent phonon resonance effect beyond the phonon-boundary scattering limit in Ge-Si core-shell ultrathin nanowires [24].
Although significant experimental efforts have been devoted to study the diameter-dependence of the thermal conductivity to understand the phonon-boundary scattering mechanisms in the nanostructures, the study of length dependence of the thermal conductivity to address the validity of Fourier's Law in 1D nanostructure remains rare. The primary technical challenges lie in the problem of suspicious thermal contact resistance and the difficulty in guaranteeing the diameter exactly the same when varying the length of the samples in different measurements. In 2008, Chang et al. developed a sequential multiprobe method that can establish the deviation from Fourier's Law behavior in one-dimensional nanostructures, without suffering from the problem of suspicious contact resistance [10]. A suspended microelectromechanical system (MEMS) device was used to measure thermal conductivity of carbon nanotubes (CNTs) and boron-nitride nanotubes (BNNTs) as a function of length. Unlike the traditional method using a movable probe as a local probe, Chang et al. deposited a series of thermal contacts to vary the length along the same sample, thus ensuring the diameter of the sample unchanged. Fig.1(a) shows experimental procedure, and the inset in Fig.  1(a) shows the nanotube bonded to the test fixture by Pt-C composite using electron-beam-induced deposition. The right contact attaches the nanotube to the top of a preformed vertical "rib" on the right thermal pad; this elevation ensures that the nanotube is fully suspended between the two contacts. In order to vary the length of the suspended segment, a series of additional thermal contacts, labeled as 1 to 5 in Fig. 1(a), was deposited inward of the original right-hand contact using Pt-C composite.
The thermal conductance of the nanotube as a function of length was measured by the method previously described by L. Shi et al. [25]. In brief, the Pt loop on one of the suspended thermal pads acts as a heater , while the other loop acts as a sensor ( Fig. ??a). Heat loss from the heater is through the low-thermal-conductance leads suspending the thermal pad, and through the sample bridging to the sensor. By measuring the temperature rises at the heater and the sensor, taking into account the heat loss through its suspending leads, the thermal conductance of the sample can be calculated. From the analysis of nanotube thermal conductance versus sam- 1. (a) upper: SEM image of a thermal conductivity test fixture with a BNNT after five sequences of Pt-C deposition. The numbers denote the nth deposition. The inset shows the SEM image after the first Pt-C deposition. The arrow denotes the preformed 'rib' for suspending the BNNT. Lower: Two circuit models for analyzing the data taking into account of the contact resistance. Rs(n) and Rc(n)denote the sample resistance and the contact resistance at each deposition, respectively.(b) Normalized thermal resistance vs normalized sample length for CNT sample (solid black circles), best fit assuming β=0.6 (open blue stars), and best fit assuming Fourier's law (open red circles). (c) Normalized thermal resistance vs normalized sample length for BNNT sample (solid black diamonds), best fit assuming β=0.4 (open blue stars), and best fit assuming Fourier's law (open red circles). For further details see in [10]. ple length, taking account into the finite contact resistance between the nanotube and the thermal pads ( Fig.1(a)), it is clear that thermal conductivity of nanotube does not follow the Fourier's law, regardless of whether L ≫ λ (λ is phonon mean free path), as shown in Fig.1(b) and Fig.1(c). Chang et al.. found that the thermal conductivity diverges with tube length as κ ∝ L β (L is the length of the nanotube; β=0 corresponds to Fourier's Law). The value of β ranges from 0.6 to 0.8 for CNTs, whereas β ranges from 0.4 to 0.6 for BNNTs.
The different length behavior might be due to the difference in isotropic disorder between BNNTs and CNTs. It is worth pointing out that the observed β on multi-walled CNTs differs from the theoretical predictions on single-walled CNTs [8]. The possible physical origin for this discrepancy might be the inter-shell phonon scattering in multi-walled CNTs.
In addition to the experimental observation of the violation of the Fourier's Law in both BNNTs and CNTs, the anomalous heat conduction was also experimentally observed in polymer nanofibres by Shen et al. [17]. Bulk polymers are generally regarded as thermal insulators because they have very low thermal conductivity on the order of 0.1 W/mK. Inspired by the theoretical work which suggests that individual polyethylene chains can have extremely high thermal conductivity [26], Shen et al. fabricated high-quality ultra-drawn polyethylene nanofibres with diameters of 50-500 nm and lengths up to tens of millimeters. Shen et al. used a sensitive bi-material AFM cantilever in combination with a micro thermocouple to measure the thermal conductivities of the individual nanofibres, as shown in Fig. 2a. It was found that the nanofibres have an extraordinary high thermal conductivity which increases with increasing draw ratios (Fig. 2b). The highest measured thermal conductivity ( 104 W/mK) is about three times higher than that of micrometer-sized fibres and 300 times that of bulk polyethylene ( 0.35 W/mK). A value of 104 W/mK is higher than many metals, including Pt, Fe and Ni. The high thermal conductivity was attributed to the restructuring of the polymer chains by stretching, which improves the fibre quality toward an ideal single crystalline fibre. Very recently, the similar anomalous high thermal conductivity is also observed in spider silk [27].

Anomalous Heat Conduction in Graphene
Although significant progress has been made for 1D systems, the study of heat conduction in 2D systems is still in its infancy due to the lack of true 2D materials, lack of 2D phonon transport theory and short of computational power . The discovery of graphene has changed this scenario [28]. The first experiment, based on Raman spectroscope measurement ( Fig. 3), was carried out by UC riverside group in 2008 [29][30][31]. This experiment was performed on suspended single layer graphene (SLG) exfoliated from HOPG. The room temperature thermal conductivity κ was measured to reach a high value of ∼4800-5300 W/mK, exceeds that of bulk graphite and carbon nanotube [16]. Such a high thermal conductivity is related to the long phonon mean free path of graphene, which is calculated to be around 750 nm at room temperature [29].
The thermal conductivity in fully supported graphene is however much smaller. The measurement of SiO 2 -supported graphene revealed a value of ∼600 W/mK at room temperature, which is due to backside scatterings and the flexural phonons leakage into the substrate [32]. By introducing an exact numerical solution for Boltzman's transport equation,  2. (a) Schematic of experimental set-up used to measure the thermal properties of a single ultra-drawn nanofibre. The thermal sensor is a silicon nitride AFM cantilever coated with a 70-nm gold film. A laser beam (wavelength, 650 nm; output power, 3 mW) is focused on the tip of the cantilever and reflected onto a bi-cell photodiode. The nanofibre drawn from the AFM cantilever is loosely suspended between a micro thermocouple and the AFM cantilever. (b) Thermal conductivities of three samples vs their corresponding draw ratios. For further details see in Ref. [17]. Seol et al. claimed that the thermal conductivity of suspended graphene should be around 5 times larger, which is ∼ 3000 W/mK at room temperature and is consistent with the results from UC riverside group.
Several following Raman based measurements, performed by other independent groups, confirmed the high thermal conductivity both in exfoliated and CVD graphene with corbino geometry [33][34][35][36]. Cai et al. reported that the thermal conductivity of suspended CVD graphene exceeds ∼2500W/mK near 350K and decreases to around ∼1400 W/mK at about 500K [33]. Lee et al. found that the value in suspended exfoliated graphene is ∼1800W/mK near 325K, and decreases to ∼700W/mK at around 500K [35]. Another group repeated the similar experiment on exfoliated graphene and found a value of ∼600W/mK [36]. The differences can be explained by the difference in the actual temperature in graphene under laser heating, optical absorption, laser spot size, and the well-know uncertainty in the Raman based measurements [37]. Despite the variation in the measured values, the thermal conductivity in both suspended and supported graphene is believed to be much larger than that in silicon and most of the metals.
Apart from Raman based measurement, typical thermal- bridge method, using the MEMS membranes mentioned above, was also used to measure the thermal conductance of graphene by the group from National University of Singapore [38,39]. Fig. 4a and 4b show scanning electron microscope (SEM) images of the suspended thermal-bridge used for the measurements of suspended and supported graphene, respectively. For details of the fabrication process and measurement accuracy, please refer to L. Shi et al. [25]. Suspended CVD SLG was measured with length of 300 nm by Xu et al. using the thermal-bridge method [39]. σ/A (thermal conductance per unit cross section area) was found to reach a high value of around 1.8×10 5 T 1.53 W/m 2 K below 100K (Fig. 5). Mingo and Broido [40] have predicted that the thermal conductance will approach the ballistic limit in clean devices, with calculated σ/A to reach 6×10 5 T 1.5 W/m 2 K. The experimentally observed values are within 30% of the predicted ballistic thermal conductance in graphene, demonstrate that the thermal transport in suspended micron graphene is indeed in the ballistic regime. Furthermore, the temperature dependence of κ can be well fitted with κ =bT α (α ∼ 1.5). This observed temperature dependence of κ(T ) and the high value in σ/A clearly identify the dominant phonons contribution to the thermal conduction in clean graphene. Mingo et al. have argued that the ZA phonons carry most of the heat in SLG; at low temperatures, the out-of-plane ZA modes are predicted to lead to a T 1.5 behavior of the thermal conductivity [41].
In another experiment, the room temperature thermal conductivity was measured to be around 1000W/mK in suspended SLG with 5 µm in length and 1.6 µm in width [42], which is 40% to 50% larger than that obtained using the same measurement technique in suspended exfoliated bilayer graphene of similar sample size [43]. This is expected due to the lack of interlayer coupling in SLG [31]. While the Raman based measurements in CVD/exfoliated SLG show an average value which is much larger [29,33,34], a direct comparison is challenging raised from the uncertainty of Raman measurements [37], and, more likely, the differences in sample geometry.
Efforts have been made to study the size effect of thermal conductivity in SLG. Chen et al. studied the thermal conductivity in suspended CVD SLG of corbino geometry with different diameters ranging from 2.9 µm to 9.7 µm using Raman method [34]. The measured thermal conductivity randomly varies from ∼ 2600 W/mK to ∼ 3100 W/mK with diameter changes, showing no sample sized dependence. The authors attributed this to a relatively large measurement uncertainty as well as grain boundaries, defect or possible polymer residues.
Alternatively, it was found experimentally that the thermal conductivity of multi-layer graphene changes with sample size [38]. In this experiment, both the suspended and supported multi-layer graphene were measured by using the suspended thermal-bridge configuration mentioned above. A room temperature thermal conductivity of 1250 W/mK was obtained for a 5 µm-long supported flake (Fig. 6). Experimen- FIG. 6. Thermal conductivity of multi-layer graphene vs. temperature. The lengths of samples S1, S2 and S3 are 5 µm (three layers), 2 µm (five layers) and 1 µm (three layers), respectively. All the samples have the same width of 5 µm. Sample S1 and S3 are supported by a thin SiNx membrane and S2 is suspended. For further details see Ref. [38].
tal results show that thermal conductivity depends strongly on the graphene size and is lowered by 85% when the length of the flake reduces from 5 µm to 1 µm (S1 and S3 in the same Figure). This length dependence is consistent with the calculation in SLG that thermal conductivity is predicted to be length dependent even when graphene flake is as long as 100 µm, due to the emergence of the long wave phonons when sample size increases [44]. More interestingly, theories had envisioned the length dependent behavior in 2D systems before the thermal conductivity in graphene was carried out [15,45]. A detailed discussion on the size dependence of the thermal conductivity in a 2D Fermi-Pasta-Ulam lattice was carried out and a ∼ logL behavior is proposed in a wide range of L and W/L ratio [46].

ATOMISTIC SIMULATIONS OF HEAT CONDUCTION IN NANO MATERIALS
It is obvious that the systematic applications of nanomaterials based devices will be greatly accelerated by a detailed understanding of their material property. For nanoscale materials, the fundamental question one may ask is whether Fourier's law is still valid. The question is not trivial, since nanoscale structures is of finite number of atoms and far from the thermodynamic limit. Due to the challenge in nanoscale experimental measurement of temperature and heat flux, atomistic simulations have an important contribution to the development of this area. In this section, we focus on the simulation investigations about the thermal conduction and diffusion in low dimensional nano materials. Because a large variety of studies on thermal transport of nanoscale materials have been done in the past decade, here we only addressed the most fundamental aspects of validity of Fourier's law in nano materials, in particular, in nanotubes, nanowires and graphene.

Anomalous Heat Conductivity and Energy Diffusion In Nanotubes and Nanowires
Carbon nanotube (CNT) is one of the promising nanoscale materials come to the spotlight of research after it was discovered in the 1990s [47]. In addition to electronic and optical properties, thermal property of CNTs has attracted more and more interests. It was found experimentally that at room temperature the thermal conductivity of a single CNT is about 3000 W/mK [16]. Recently, thermal contact resistance of carbon nanotube junctions have been explored experimentally [48]. Both electron and phonon can be heat carrier. Yamamoto et al. has demonstrated that even for the metallic nanotubes, the electrons give limited contribution to thermal conductivity of single-walled CNTs (SWNTs) at low temperature, and with increase of temperature this part decreases quickly [49]. It is found that the phonon mean free path of CNT can exceed the characteristic length of the structural ripples. This makes CNT an ideal phonon waveguide which can have superior phonon transport property to its counterpart of photon in optical fibers [50]. Thus SWNT is ideally suited for molecular dynamics (MD) investigations of the phonon thermal conduction law in low dimensional system. By using MD simulations, it was found that the thermal conductivity of SWNT diverges with the length of system as κ ∼ L β , where the exponent β depends on the temperature and SWNT diameters, and the value of β is between 0.12 and 0.4 [51] and between 0.11 and 0.32 [52]. As discussed in above section, this length dependent thermal conductivity has been verified experimentally, although the absolute value of β varies.
To understand the physical mechanism for the length dependent thermal conductivity in SWNT observed both experimentally and theoretically, Zhang and Li investigated the vibration energy diffusion along SWNT [8]. In the calculation the Tersoff empirical potential is used to derive the force form, which has been widely used in the study of heat conduction in nanostructured carbon systems. They first thermalize the CNT to a temperature T, then a packet of energy (heat pulse) is excited in the middle of the tube to study how it spreads along the system. To suppress statistical fluctuations, an average over 1000 realizations is performed. Fig. 7(a) shows the representatives of pulse propagation along the lattice at 2K. The pulse diffusion profiles at room temperature are very similar to those at low temperature. In the SWNT, one single peak expands with time, and it is reflected back from the two ends after the wave front reaches the tube ends, no peak splits are observed. Indeed, there exist three time periods in the diffusion process. At the start period, the pulse spreads between the two nearby layers, which is denoted as "initial period." Then, the pulse diffuses along the tube, this is "diffusion period"; finally, the wave front hits the boundary and is reflected back and transport in opposite direction, which is "reflection period." The snapshots of representative pulse propagation along the CNT at room temperature is shown in Fig. 7(b). It is clear that in CNT, one peak expands with time, this pulse spread can be described quantitatively as below: is the position of energy pulse at t = 0. The averaged energy profile spreads as: where . denotes the ensemble average over different realizations.
In Fig. 7(c) and (d) we show σ 2 (t) versus time in double logarithmic scale, so the slope of the curve gives the value of α. It is clearly seen that the slope is two at 2K. This means that phonon energy transports ballistically at low temperature. This can be understood from the Taylor expansion of the atomic interaction potential by keeping up to the secondorder term. At low temperature, the vibrations of atoms are very small, the Tersoff potential can be approximated by a harmonic one, thus corresponding to a ballistic transport. However for SWNT at room temperature, energy transports superdiffusively with α ≈ 1.2. This is slower than ballistic transport ( α = 2 ) but faster than normal diffusion ( α = 1 ). With temperature increases, the anharmonic terms appear due to excitation of the transverse vibrational mode. According to a theory proposed by Wang and Li [53], it is known that for a quasi-1D lattice model, the interaction between the transverse modes and the horizontal modes will result in a superdiffusive phonon transport. Moreover, Li and Wang [54] demonstrate a connection that the anomalous energy diffusion induces the divergent thermal conductivity. Combine the observed superdiffusion with α ≈ 1.2 and their theory, the thermal conductivity diverges with tube length for the tube length up to few micrometers, namely, the thermal conduction of SWNT does not obey Fourier's law. This provides the physical mechanism for the experimental and theoretical observed length dependent thermal conductivity of CNTs.
Like in CNTs, the length dependent thermal conductivity is also an important physical question to SiNWs. By using NEMD simulation, Yang et al. have studied the length dependence of thermal conductivity of SiNWs [9]. In order to establish a temperature gradient along the SiNW, the atoms close to the two ends are put into heat bathes with high and low temperature. Nosé-Hoover and Langevin heat bathes [74] are applied to ensure the conclusions are independent of heat bath. The dependence of thermal conductivity on the length of SiNW is shown in Fig. 8. Both types of heat baths give rise to the same results. It is obvious that the thermal conductivity increases with the length as, κ ∝ L β , even when the wire length is as long as 1.1µm. Traditionally, the phonon mean free path (λ) is a characteristic length scale beyond which phonons scatter and lose their phase coherence. In three-dimensional systems, it gives Fourier's law when system scale L ≫ λ . By using the phonon relaxation time (∼ 10ps) in SiNWs, and the group velocity of phonon as 6400m/s, the mean free path λ is about 60nm [9]. The maximum SiNW length (1.1µm) in this study is obviously much longer than the phonon mean free path of SiNW. This demonstrates that in SiNWs, Fourier's law is broken even when the length is obviously longer than the traditionally mean free path.
More interestingly, it is found that the length dependence of thermal conductivity is different in different length regimes. At room temperature, when SiNW length is less than about 60nm, the thermal conductivity increases with the length linearly ( β ≈ 1 ). For the NW with length larger than 60nm, the divergent exponent β reduces to about 0.27. This critical length (60nm) is just the value of mean free path, then this length dependent divergent exponent can be fully understood as below. There is weak interaction among phonons when the length of SiNW is shorter than mean free path, thus the phonons transport ballistically, like in the harmonic lattice. However, when the length of SiNW is longer than the mean free path, phonon-phonon scattering dominates the phonon transport and the phonon cannot flow ballistically.
In addition, the diverged exponent β also depends on temperature. At 1000K, β is only about 0.15 when NW length is longer than 60nm. The decrease of β can be understood from the temperature dependent phonon interaction. In SiNWs at high temperature, the displacement of atoms increases, which induces more phonon-phonon interaction, and reduction in β.
The divergence of thermal conductivity can be understood from the phonon density of states spectra (PDOS) as shown in Fig. 9. In the left column of Fig. 9, the higher energy part (short wavelength) in PDOS spectra is not sensitive to the wire length. The peak at about 16 THz in PDOS of bulk silicon appeares even when the NW length is only 4 nm. However, as shown in right column of Fig. 9, the lower frequency part in PDOS spectra is much different from that of bulk Si. In contrast to the continuous spectra of bulk silicon, there are many discrete peaks in PDOS spectra of NW. For a short SiNW, the low energy phonon density is very low, thus the thermal conductivity is low. With the increase of length, more and more (long wavelength) phonons are excited, which leads to the increase of thermal conductivity. However, due to the size confinement effect, the energy density of acoustic phonon in PDOS spectra of SiNW is much smaller than that of bulk Si, and leads to the low thermal conductivity of SiNW compared to that of bulk Si.
So far we show the length dependent thermal conductivity of SiNWs. Now let us turn to the energy diffusion process in SiNWs. The heat energy diffusion in SiNW was firstly simulated by Yang et al. in 2010 [9]. Followed the approach applied in diffusion in CNT by Zhang and Li [8], they first thermalize the system to an equilibrium state with temperature T 0 , then atoms in the middle layer are given a much higher temperature T 1 . The evolution of the energy profile along the chain is then recorded afterwards. Fig. 10 shows σ 2 (t) versus time in double logarithmic scale, so that the slope of the curve gives the value of α. For the NWs of length of 140nm, we obtain α = 1.15 and 1.07 at 300K and 1000K, respectively. Different theory describes the physical connection between energy diffusion and thermal conductivity [54,75]. For instance, in normal diffusion, the phonon transports diffusively, it corresponds to a size-independent thermal conductivity. This is what we have in bulk material and the Fourier's law is valid. However, in the ballistic transport, the thermal conductivity of the system is infinite when the system goes to thermodynamic limit. This is the case for one-dimensional harmonic lattice. In another case, if α < 1, which we call sub-diffusion case, corresponds to β < 0, namely, the thermal conductivity of the system goes to zero. Thus the system is an insulator. In the super diffusion regime, it predicts that the thermal conductivity increases as the length of the system increases. This is what we observed in SiNWs.
Thus the anomalous heat diffusion is responsible for the length dependent thermal conductivity in SiNWs. Combine with the results in CNTs, all these reports provide strong evidence that Fourier's law of heat conduction is not valid in one dimensional and quasi one dimentional nanostructures. Moreover, these theoretical studies demonstrate that nanowire and nanotube are promising platforms to verify phonon transport mechanisms.

Computational Study on Thermal Conductivity of Graphenes
Besides to the one-dimensional nano materials, Graphenes [28,76] have attracted immense interests recently, mostly because of their unusual one-atom-thick structure. Superior thermal conductivity (as high as 5000 W m − 1 K − 1) has been observed in graphene [30], which has raised the exciting prospect of using them for thermal devices.
Based on graphene, narrow (sub-10 nm) graphene nanoribbon (GNR) can be obtained either by cutting mechanically exfoliated graphenes, by patterning epitaxially growth graphenes or by unzipping carbon nanotubes (CNTs) [28,77,78]. GNRs are materials with distinctive electronic and thermal transport properties and are candidates for field effect transistors of future carbon-based nanoelectronics. There are rich physical phenomena about thermal property of GNRs. For instance, intrinsic anisotropy originate from different boundary condition at ribbon edges is observed in narrow GNRs, results to the room temperature thermal conductance of zigzag GNRs is about thirty percent larger than that of armchair GNRs [79,80]. The effects of edge, roughness, and hydrogen termination have been investigated by different simulation methods [81][82][83][84][85][86]. And isotopic doping can reduce remarkably the thermal conductivity of graphene [87]. In addition to the isotope randomly doping, the impact of disorder defect on thermal conductance [88,89] and thermoelectric figure of merit of GNRs [90] has been investigated by using density functional theory calculations combined with the nonequilibrium Green's function method. It is found that the figure of merit can be remarkable enhanced five times by randomly introducing vacancies to the graphene. And remarkable thermal conductivity enhancement has been found in graphene nanoribbons under homogeneous uniaxial stretching due to a lot of dispersive phonon modes are converged to the low frequency region [91]. And thermal conductivity of GNRs is very sensitive to tensile strain [92]. More interesting, the presence of a substrate [93,94] and layer-layer interaction [95][96][97][98] lower the thermal conductivity. The excellent thermal property of GNRs makes them ideal candidate for nanoscale phononic devices [95,[99][100][101].
Similar to the cases in CNT and SiNW, thermal conductivity of GNRs also depends on their sizes. By using equilibrium molecular dynamics simulations, the effect of GNR width on thermal conductivity has been studied [81]. They found that the thermal conductivity increases with width, but for width larger than 5nm the thermal conductivity becomes relatively size independent. The larger conductivity values at very small sizes is explained from the limited number of phonons in the system. In their simulation, the width range is from 2nm to 12nm. The similar width dependence of thermal conductivity was also observed by non-equilibrium molecular dynamics simulations [102]. Moreover, thermal conductivity also depends on GNR length. In their work, for GNR with width of 2nm, its thermal conductivity increases from about 200 W/mK to 400W/mK when length increases from 10nm to 60nm. Based on the extrapolation procedure, they predicted the thermal conductivity of GNRs with a length of 2 micrometer will be about 3000W/mK, which is close to the experimental reports. However, as the maximum length in their simulation [102] is only about 60nm, it is not long enough to provide detailed information about the divergent exponent. Further simulation studies on much larger GNRs will be greatly helpful to explore whether Fourier's law is still validity in twodimensional one-atom-thick systems (See also [103]).

HEAT TRANSPORT IN LOW DIMENSIONAL MODELS
The anomalous heat transport phenomena observed in the above mentioned low dimensional nanostructures can be partially understood by general theories on low dimensional models, which generally demonstrate that the phenomenological Fourier's law may not be valid. Actually, ever since Fourier proposed this law on heat conduction about 200 year ago, enormous number of studies has been conducted to fully understand it from fundamental physics. But a rigorous derivation from microscopic Hamiltonian dynamics is still absent. However, from these studies, we are already able to draw the surprising conclusion that heat conduction in low dimensional systems may not follow the Fourier's law. In this section, we will briefly review the fundamental studies on heat transport in low dimensional systems: from theories to numerical simulations on the toy models. Especially, we introduce the works that study anomalous heat conduction from the aspect of anomalous energy diffusion.
In order to understand the realistic models from a rigorous perspective, it is better to begin with the simplest models which involve only the basic but important ingredients, such as linearity and disorder, that will possibly lead to the expected normal or anomalous transport behaviors.

Low Dimensional Harmonic Lattices
The ordered harmonic lattice is the simplest lattice model that we can start with. From the early Boltzmann transport equation, it is not surprising that this linear model does not have well defined transport coefficients due to the lack of interaction between phonon modes. Therefore, heat transport in harmonic models should not obey Fourier's law. Actually, a clearer clue was found after Debye extended Boltzmann kinetic theory of ideal gas and established the expression of heat conductivity κ = cvl, where c is the heat capacity, v the phonon velocity, and l the phonon mean free path, because the non-interacting phonons should have infinite mean free paths.
The first explicit results on heat transport in classic harmonic models were given by Rieder, Lebowitz and Lieb [1]. They studied a harmonic chain connected to stochastic Langevin heat baths. The lattice has Hamiltonian H = N l=1 The particles at the ends are subjected to additional forces F 1,N = η 1,N − γv 1,N , where the noise terms η 1,N are independent Gaussian random processes with mean zero and variances η 1,N (t)η 1,N (0) = 2γk B T 1,N δ(t). Rieder et al. proved that, at large N , the heat flux will saturate to [1] where ν = mk γ 2 . This heat flux is independent of the length N for large N , which demonstrates the expected ballistic transport. Another important result on harmonic model is that the non-equilibrium temperature profile is flat at the center part with temperature T = (T 1 + T N )/2. Temperature jumps occur at the boundaries. In this sense, temperature gradient cannot be established in harmonic lattices. Another work by Nakazawa extended Rieder et al.'s model by introducing on-site harmonic potential to all sites and a similar lengthindependent heat flux was found [2]. In the same work [2], Nakazawa also proved that high dimensional harmonic lattices can be reduced to a 1D problem. Therefore the conclusions are similar.
Due to the absence of scattering mechanism in harmonic models, heat transport in these lattices is ballistic. Therefore, it is necessary to introduce scattering of phonons to get the desired diffusive transport. Basically, two approaches are commonly adopted: by introducing disorder or nonlinearity to the system. Both of these two methods can cause phonon-phonon interaction. However, as we will see later, neither method works effectively for diffusive transport in low dimensional lattices.
For the first approach, we consider the harmonic lattices with disorder. Disorder can be generally introduced by randomly assigning varies particle masses or spring constants, or both. These two types of disorder do not change the underlying physics very much. The former with random masses is just like isotopic doping, which is more interesting and relevant in applications. Therefore, it is studied intensively. The presence of disorder will generally cause localization of normal modes, which could be partially understood as an analogue to the Anderson localization of electrons. In the Anderson tightbinding model, all the eigenstates of electrons are localized in one and two dimensional cases, therefore the systems are electric insulators. However, in the phonon case, the picture is much more complicated.
For one dimensional disordered harmonic chains, Allen and Ford first observed that the thermal conductivity of an infinite disordered harmonic chain is finite [104], unlike the Anderson model of electrons in which the conductivity will exponentially decay to zero. Furthermore, the eigenstates of electrons in the strongly disordered cases are all exponentially localized. However, in disordered media, Matsuda and Ishii [105] showed that only the high frequency phonons has this property, while phonons with low frequencies ω < ω d can be considered as delocalized phonons, since the localization lengths are greater than the system length. The characteristic frequency ω d depends on the variance of mass distribution ∆m 2 as ω d ∼ ( km N ∆m 2 ) 1/2 , where k, m and N are the spring constant, the average mass and the lattice length, respectively. Therefore, low frequency phonons will contribute to heat transport, which prevents the disordered harmonic lattices from being a thermal insulator. Matsuda and Ishii also studied the thermal conductivity of disordered chains connected to white noise heat baths (Langevin baths) and baths modeled by ordered semi-infinite harmonic chains (Rubin baths), respectively [105]. They found that for both models (a): disordered harmonic chain connected to Langevin baths with fixed boundary conditions; and model (b): disordered harmonic chain connected to Rubin baths with free boundary conditions, their thermal conductivity will diverge with the system length as κ ∼ N β with β = 1/2. However, Casher and Lebowitz later rigorously proved that the correct exponent for model (a) should be β = −1/2 [106], while for model (b) the result β = 1/2 is supported by numerical simulations [107] and a rigorous proof was provided by Verheggen later [108].
Altogether, the thermal transport properties of disordered harmonic chains would probably depend on the both the boundary conditions and the heat baths. To clarify this point, Dhar restudied the thermal transport in disordered harmonic chains connected to various baths modeled by generalized Langevin equations [109]. Using a Langevin equations and Green's function (LEGF) formalism, Dhar found that in disordered harmonic chains, the exponent β depends on the spectral properties of the baths. The previous model (a) and (b) are only two special cases of the generalized Langevin heat baths, which have β = −1/2 and 1/2, respectively. A special choice of heat bath spectrum can even lead to normal heat transport with β = 1, obeying Fourier's law [109]. A more detailed study on the effects of the boundary conditions on the transport properties were carried out by Roy and Dhar later using the same LEGF formalism [110]. The results are a little contradictory to those of [109], which states that the different exponents β are indeed dependent on the different boundary conditions, instead of the spectral properties of the baths. Clear numerical evidence also supports this conclusion as shown in Fig. (11), i.e., for disordered harmonic chains with free boundary conditions β = 1/2 and with fixed boundary conditions β = −1/2, no matter whether Langevin baths or Rubin baths are used. Ref. [110] also showed that, in the low frequency regime, the transmission coefficient of the disordered chains T (ω) can be approximated by the transmission of the ordered chains T ordered (ω). Therefore, heat flux in a disordered harmonic chain can be estimated as where ω d is the low frequency boundary ω d ∼ ( km N ∆m 2 ) 1/2 . Considering the fact that the transmission coefficients for low frequency modes scale as T ordered (ω) ∼ 1 for free boundary conditions and T ordered (ω) ∼ ω 2 for fixed boundary conditions, respectively, we can obtain J free ∝ (T L − T R ) ω d ∝ (T L − T R )/N 1/2 for free boundary conditions and J fixed ∝ (T L − T R ) ω 3 d ∝ (T L − T R )/N 3/2 for fixed boundary conditions, respectively.
Till now, we only discussed the disordered harmonic chain without onsite pinning potential (except the case of fixed boundary conditions, the particles at the ends are pinned). The absence of pinning in this model makes the system translational invariant in a coarse grained manner. As a consequence, the phonon modes are not fully localized and the low frequency ones can transport heat. When onsite potentials are involved, the picture is a little different. In their work [111], Dhar and Lebowitz proved that with onsite potentials, all the phonon modes are localized. Consequently, the heat flux J decays exponentially with the chain length and the system is a thermal insulator in the thermal dynamic limit. This could be understood because the onsite potentials will typically break the translational invariance of the chain and eliminate the low frequency phonon modes in this structure. The remaining phonon modes have high frequencies and are localized.
Compared with the one dimensional lattices, the studies on two dimensional disordered ones are much less due to the mathematical difficulties. In a renormalization group study by John et al. [112], it was found that the low frequency modes in 2D disordered harmonic lattices are not localized either as in 1D cases. The localization length diverges as e 1/ω 2 in the low frequency regime. Therefore, modes with frequency ω < ω d = log −1/2 (N ) can be considered as extended, which will contribute to heat transport.
In another study by Lee and Dhar [113], 2D disordered harmonic lattices connected to two types of stochastic heat baths are investigated. The first type is modeled by uncorrelated Gaussian processes and the other by exponentially correlated Gaussian processes. In their work, disorder is introduced as isotopic doping. A mass M = 2 is randomly assigned to half of the particles, while the rest have mass m = 1. The authors simulated square lattices with widths up to L = 256 and found a power law divergent heat conductivity in both two types of models. For uncorrelated baths the exponent β ≈ 0.41 while for correlated baths β ≈ 0.49. In the same work [113], Lee and Dhar also studied a special case of correlated disorder in which the lattice is only disordered in the conducting direction. It was found that this model can be transformed into an effective one dimensional problem thus being mathematically tractable. Analytical results showed that the exponent β = −1/2 which was also verified by numerical simulations using square samples with widths up to L = 512. In another work [114], Yang considered 2D lattices with bond-missing defects connected to Nosé-Hoover heat baths. The results showed that when the defect density is large enough, temperature gradient can build up and the heat conductivity is finite, while for the case of small defect density, the conductivity will diverge logarithmically.

Low Dimensional Anharmonic Systems
It is clear now that disorder alone cannot effectively make the phonon transport in the harmonic chains diffusively. As an alternative approach to introduce phonon-phonon scattering, the effect of nonlinear interaction on the properties of heat transport in low dimensional systems has also been intensively studied, especially by numerical simulations. Because of the intrinsic non-integrability of the system Hamiltonian, analytical results are rare.
Basically, there are three different approaches to study heat transport analytically: the mode coupling theory (MCT) [4,115], the renormalization group theory [45] and the Peierls-Boltzmann kinetic theory [116]. All are based on Green-Kubo formula [117,118] where C(t) = J (t)J (0) /L is the time correlation of the total heat flux J . Strictly speaking, Green-Kubo formula can only be applied to infinite systems, as indicated by the limit taken lim L→∞ , which should be performed before lim t→∞ strictly. However, as we already know, in many low dimensional systems, the heat conductivity in thermal dynamic limit is divergent. The focus of such anomalous situations is how the conductivity diverges with the system size. To utilize the Green-Kubo formula to study the length dependence of heat conductivity in these cases, a usual procedure is to truncate the upper limit of the integration in Eq. (8) to t c ∼ L/v s . This is commonly believed as that the sound waves propagate to the boundaries at a finite speed v s , which will lead to a fast decay of correlation C(t) at time ∼ L/v s [14,15]. Using this treatment, the previously mentioned three theories, i.e. mode coupling theory, renormalization group theory and Peierls-Boltzmann equation kinetic theory, are all aimed to calculate the long time tail of the correlation function C(t). A basic result shared by all these theories is that the systems with nonlinear interaction remain anomalous except those subjected to onsite potentials. For 1D cases without onsite potential, these theories predict that the conductivity of the system would diverge as a power law κ ∼ N β . However, the values of β differ from one to another. For 2D cases, the mode coupling theory and the renormalization group theory predict logarithmic divergence.
The mode coupling theory approach was first used by Lepri et al. [4,115]. The basic idea is that the divergence of heat conductivity is due to the slow relaxation of phonon modes with long wavelengths. A later development of this theory was done by Delfini et al. [119]. By constructing a correlation function G(k, t) = Q * (k, t)Q(k, 0) of the normal mode coordinates Q(k) = 1 √ N N n=1 x n exp(−ikn), the authors were able to get the exact evolution equations for G(k, t) using the Mori-Zwanzig projection approach [115]. Using some approximations, these equations were solved self-consistently for the Fermi-Pasta-Ulam (FPU) chains with interaction potential U (x) = k 2 x 2 /2 + k 3 x 3 /3 + k 4 x 4 /4, which gives G(k, t) = A(k, t)e iω(k)t + c.c. for a small wavenumber k. In this formula, ω(k) is the temperature dependent dispersion relation which is obtained using harmonic approximation ω(k) = 2| sin(k/2)|, and A(k, t) has the form (9) Finally, the heat flux correlation function C(t) was obtained from C(t) ∝ k ( dω(k) dk ) 2 G 2 (k, t). It was shown that C(t) ∼ t −2/3 for the case k 3 = 0 and C(t) ∼ t −1/2 for k 3 = 0 but k 4 = 0. Therefore, by inserting C(t) into the Green-Kubo formula and adopting the cut-off time t ∝ N , we obtain β = 1/3 and β = 1/2 for each case, respectively. The mode coupling theory can also provide predictions for higher dimensional systems. It was shown that [14] for two dimensional systems, C(t) ∼ t −1 , so that the heat conductivity would diverge logarithmically with system size. While for three dimensional systems, C(t) ∼ t −3/2 , with which we obtain normal heat conduction.
In a later work, Wang and Li investigated the effect of transverse degrees of freedom in one dimensional chains using this mode coupling theory [53]. The system under investigation has Hamiltonian where (r i , r i ) is the canonical coordinates, and φ l is the angle between vectors r i+1 − r i and r i−1 − r i . This model can be regarded as a simplification of more realistic polymer chains. The model is still one dimensional but the motions of particles are two dimensional. Both numerical simulations and mode coupling analyses suggest that in the presence of interaction with transverse modes, β = 1/3. When the coupling between transverse modes and longitudinal modes is weak, β = 2/5. The second approach was proposed by Narayan and Ramaswamy which uses hydrodynamic equations and renormalization group theory [45]. It was argued that in an interacting 1D system with large system length, the thermal fluctuation will wipe out the long-range order, and consequently the system will behave like fluid. Therefore it is expected that this approach can be used generally to describe heat transport in 1D phenomenon. By assuming that the only conserved quantities in the system are the total number of particles, the total momentum and the total energy, one can obtain three hydrodynamic equations describing the evolution of the particle density field and the velocity field, with addition terms describing the thermal noise. Narayan and Ramaswamy then solved these equations using linear response approximation and finally obtained the thermal flux correlation by considering the symmetries of the system. It was demonstrated that the flux correlation function C(t) ∼ t −2/3 asymptotically for the 1D case. Therefore, one obtains β = 1/3 by the cut-off time reasoning. For 2D systems, the authors showed that the conductivity would diverge logarithmically with system sizes [45].
In another work by Mai and Narayan [120], the authors claimed that the renormalization group theory analyses can be applied to lattice models as well, even for stiff one dimensional oscillator chains including the FPU model. The conductivity of 1D lattices should also diverges in the same way as fluid, i.e. diverges in a power-law with exponent β = 1/3. However, in a consequent work by Hurtado [121], the breakdown of hydrodynamics in a simple one dimensional fluid was reported.
The estimation of the exponent β in anomalous heat transport using Peierls-Boltzmann equation is first carried out by Pereverzev [116]. The Boltzmann equation is originally used to describe the the phase space density evolution of kinetic gases. From it we have already known that the heat conductivity can be written as κ = c v l, where c, v and l are the specific heat, sound velocity and mean free path of the gases, respectively. Peierls developed this theory to describe the phonon transport in solids and found a similar formula for the conductivity κ ∼ c k v 2 k τ k dk, where k specifies the phonon modes and τ k is the phonon relaxation time. Using this Peierls-Boltzmann approach, Pereverzev studied the FPU-β mode [116]. Considering that anomalous heat conduction is mainly due to the ballistic-like behavior of phonons with small wavenumber k. the dependence of τ k on the wavenumber k in this regime is the main interest. After some approximation, Pereverzev was able to obtain τ k ∼ k −5/3 for small k. It was also found that heat flux correlation function can be calculated by Because v k is almost constant for small k, one obtains C(t) ∼ t −3/5 and β = 2/5. This result was confirmed later by Lukkarinen and Spohn using a more rigorous Peierls-Boltzmann approach [122]. For momentum non-conserving cases, a finite heat conductivity was obtained by Aoki et al. [123] using this approach.
Compared to the few analytical results on low dimensional heat conduction, numerical results are much more abundance.
Here we only briefly highlight the results on some well studied lattice oscillator models.
Due to the landmark work by Fermi, Pasta and Ulam, the FPU model has received great attraction in understanding nonlinear statistical mechanics. In the FPU-β model, interaction between nearest particles has the form U (x) = k 2 x 2 /2 + k 4 x 4 /4. This model is one of the simplest lattice models which can be used to investigate the effect of nonlinearity on heat transport.
The first numerical study on heat conduction of FPU-β model was carried out by Kaburaki and Machida [124]. In their study, the ends of FPU-β lattices were connected to Boltzmann heat baths at different temperatures. Detailed calculations showed that linear temperature profiles can be established in the steady state. Moreover, it was found that the conductivity tends to saturate when the chain length increases. However, it is clear now that neither of these findings reveals the true heat transport behaviors of FPU-β lattices.
The first indication of anomalous heat conduction violating Fourier's law in FPU-β lattices was presented by Lepri et al. [3]. FPU-β lattices connected to Nosé-Hoover heat baths with lengths up to N = 400 was studied. It was found that the heat conductivity diverges as κ ∼ N β with exponent β = 0.55 ± 0.55. Another finding in this work is that the temperature profile is non-linear even when the temperature bias is very small, a result contradictory to the prediction from the Fourier's law. In a subsequent work [4], the divergent exponent was improved to β = 0.37 with length up to N = 2048. These pioneering works in one dimensional lattices have inspired enormous subsequent studies on anomalous heat conduction in low dimensional systems.
It should be mentioned that the system lengths in the studies above are too small to draw any quantitative conclusions. Actually, in non-equilibrium simulations, the heat conductance will be greatly affected by the contact resistances between the lattice and heat baths. Therefore, it is necessary to go to a very large system length to make the contact resistances negligible. In a recent work [125] In the rightmost part of the figure, it can be seen that the running slope increases as the system length is increased. Solid and dashed lines with slope 2/5 and 1/3, respectively, are drawn for reference. For more details please see Ref. [126]. this into consideration. They studied the FPU-β chain with length up to N = 65536 using two different types of heat bath: the stochastic Langevin heat baths and the deterministic Nosé-Hoover heat baths. It was found that when the system size is small, the calculated thermal conductances differ a lot when different heat baths are used. Meanwhile, when the system size is large enough, the discrepancy between the calculated heat conductances is negligible, which ensures us that the boundary effects could be omitted. In this situation, it was then found that the exponent β decreases to a final result β = 0.333 ± 0.004. The authors claimed that this result is consistent with the renormalization group theory [45] which supports its universality.
However, in a quite recent work, Wang et al. [126] restudied the FPU-β model and the purely quartic model with an even longer length up to N = 131072. The interaction in purely quartic model is U (x) = 1 4 k 4 x 4 . This model is the high temperature limit of the FPU-β model. They found that the slope d ln κ/d ln N is not monotonically decreasing as shown by Mai et al. [125]. The slope starts to increase when N is close to 10 5 and finally reaches β = 2/5 instead of 1/3 (see Fig. 12). In order to remove the boundary effects, the authors also performed equilibrium simulations without heat baths to calculate the heat flux correlation function C(t) in the Green-Kubo formula. Again, the slope d ln C(t)/d ln N increases when t is very large. The heat flux correlation finally saturates at C(t) ∼ t −3/5 , which supports their nonequilibrium results β = 2/5 and the Peierls-Boltzmann theory [116].
For one-dimensional lattices with momentum conserving interparticle interactions, all theories predict a divergent heat conductivity κ ∝ L β with universal β. However, in a most recent study on FPU-β lattices with nearest-neighbor (NN) and next-nearest-neighbor (NNN) coupling [127], Xiong et al. found that the exponent β strongly depends on the ratio of the NNN coupling to the NN coupling γ. The Hamiltonian of the FIG. 13. Dependence of β on γ in the model with Hamiltonian Eq. (11). For more details please see Ref. [127]. system under investigation is where V (x) is of FPU-β type 1 2 x 2 + 1 4 x 4 . In their study, the authors used the reverse nonequilibrium molecular dynamics method (RNEMD) [128] to simulate the system. Firstly, this method was tested under condition γ = 0 and length L = 2496 which showed that the temperature gradient was well established and β = 0.325 ± 0.002. Consequent simulations using different γ revealed that the exponent β varies continuously with γ. Starting from β = 0.325 when γ = 0, β keeps decreasing until reaching its minimum 0.25 at γ = 0.25. After this point, β increases to 0.35 at γ = 1 as shown in Fig. 13. Therefore, it was suggested that a universal exponent β does not exist. However, it should be mentioned the result β at γ = 0 is quite different from Wang et al.'s which was calculated using both non-equilibrium method with Langevin heat baths and equilibrium method with Green-Kubo formula with length up to 131072 [126].
For momentum non-conserving models, it has been mentioned that all theories predict a finite heat conductivity. Numerical simulations also support this conclusion. Actually, before the three theoretical predictions were announced, Hu et al. had realized that momentum non-conservation should be a necessary ingredient to obtain the Fourier's law from their work on the Frenkel-Kontorva (FK) model with a periodic onsite potential V (x) ∝ cos(ax) in addition to the harmonic interpartical interaction [129], Two other studies, one by Aoki et al. [130], and the other by Hu et al. [131], on φ 4 model with onsite potential V (x) = k 4 x 4 /4 further confirms this conclusion.
What is the disorder effect to anharmonic lattices? This problem was first investigated by Payton et al. [132]. They performed non-equilibrium simulations on disordered anharmonic lattice connected to stochastic baths. It was found that the anharmonicity could greatly enhance the heat current. But due to the limited computer facility at that time, they were not able to verify the Fourier's law. The first systematic study of the disorder effect in anharmonic lattice was taken by Li et al. [133]. In their work, Li et al. studied the FPU chain with mass disorder connected to Nosé-Hoover heat baths. Their results showed that at low temperatures, the system obeys the Fourier's law. While at high temperatures, the conductivity diverges with exponent β = 0.43. However, in a later study by Dhar and Saito [134], the authors found that the exponent β is still 1/3 at low temperatures if one goes to much larger system size. Therefore, there is no crossover from diffusive to superdiffusive in disordered FPU chain. Disorder plays less important role in β in anharmonic chains.
For disordered anharmonic chains with pinning, Dhar and Lebowitz studied the general Hamiltonian systems [111] in which x 0 = x N +1 = 0. It was found that when anharmonicity is absent, which corresponds to the disordered harmonic chain with harmonic onsite potential, the heat flux decays exponentially with length and the system is a thermal insulator. However, introduction of a small will lead to the J ∼ 1/N dependence, resulting in a diffusive heat transport arises.
Regarding the two dimensional oscillator lattices, only few numerical studies are carried out to investigate the size dependent heat conductivity. But the good news is that predictions from theories do not contradict one another. Both the renormalization group theory and the mode coupling theory predict a logarithmic divergence for the heat conductivity.
For numerical investigations, Lippi and Livi studied two dimensional lattices with size N x × N y [135]. Heat was transported along x direction. The Hamiltonian of the system has the form where p ij , x ij are the canonical coordinates. The interaction was taken to be the FPU-β type U (x) = x 2 /2 + k 4 x 4 /4 and the Lennard-Jones type U (x) = A/x 12 − B/x 6 + B 2 /4A, respectively. The parameters in the Lennard-Jones potential were chosen so that its Taylor expansion at the minimum coincide with the FPU-β potential. Using non-equilibrium method with Nosé-Hoover heat baths, Lippi and Livi firstly demonstrated that with increasing N y and fixed N x , the current saturates at a small ratio N y /N x . Subsequent calculations were performed using fixed ratio N y /N x = 1/2. It was found that the conductivity will diverge logarithmically with N x for both models. Additional equilibrium simulations showed that the heat flux correlation functions decay as t −1 . Using Green-Kubo formula, we again obtain κ ∼ ln(N x ). This result is consistent with predictions from theories.
There also exists numeric works that contradicts the prediction of logarithmic divergence. Yang and Grassberger studied a 2D Hamiltonian system involving FPU type interaction [136]. The Hamiltonian of the system is similar to Eq. (13) but with scalar displacements. The results showed that when the ratio N x /N y is small which corresponds to the 2D behavior, the conductivity diverges exponentially with β = 0.22 ± 0.03. When N y = 1 which is actually a 1D system, β = 0.37±0.01. In another recent work, Shiba and Ito also studied the same system as Eq. (13) with the FPU-type interaction [137]. However, a power law divergence was obtained by using lattices with sizes up to 384 × 768. It is found the conductivity also diverges in a power-law with β ≈ 0.268, even it was assured that the boundary effect is negligible in that size.

Anomalous Energy Diffusion in Low Dimensions
The nonlinearity of the interaction in the lattices prohibits us to analytically solve most of the models. Instead, a lot of quasi-one-dimensional gas channel systems consists of noninteracting particles have been proposed to study heat transport in low dimensions. These works helped a lot in understanding the basic ingredients required by the Fourier's law. Moreover, they helped to understand heat transport from another aspects of energy diffusion.
Alonso et al. first studied the heat conduction in Lorentz gas channels [138]. The model consists of a quasi-1D billiard with periodically distributed semicircular scatterers as shown in Fig. 14(a). Particles carrying heat move along the billiard and exchange energy with Maxwell thermal reservoirs at the ends. In this model, no particle can move between the two reservoirs without being scattered by the semicircles. Due to the exponential separation of trajectories, this system has a positive Lyanpunov exponent, i.e., it is chaotic. Alonso et al. verified that the heat conduction obeys Fourier's law which was explained by deterministic particle diffusion with energy dependent diffusivity D(E) ∼ E 1/2 .
Later, in order to resolve the role of dynamic chaos on determining the heat transport behavior, Li et al. studied the Ehrenfest gas channels with isosceles right angles replacing the semicircles in the Lorentz gas channels [139]. Their configurations are shown in Figs. 14(a) and 14(b). The linear separation of nearby trajectories in these models make them different from the Lorentz gas channels in underlying dynamics. It was found that when the scatterers are periodic, the heat conductivity diverges with system length as κ ∼ N 0.81 . Normal heat conduction κ ∼ N 0 can only be reached when either position disorder (Fig. 14(b)) or size disorder (Fig. 14(c)) is introduced. These results imply that chaos is not necessary for normal heat conduction and Fourier's law. Moreover, the diffusive behavior of the particles are also investigated using the mean square displacement. Numerical simulations showed that for the disordered case, the mean square displacement increase linearly with time, while the periodic case is superdiffusive ∆x(t) 2 ∼ t α with α = 1.672. Subsequently, two additional works on quasi-1D billiard model but with different scatterers are studied by Alonso et al. [140] and by Li et al. [141], respectively. In the first one, the scatterers have polygonal shape (Fig. 14d). By fixing φ 1 = ( √ 5 − 1)π/8 and adjusting φ 2 = π/q with q = 3, 4, 5, 6, 7, 8, 9, Alonso et al. found that the time dependence of the MSD ∆x(t) 2 ∼ t β is superdiffusive for q = 3 (β = 1.30), subdiffusive for q = 4 (β=0.86) and normal for q = 5, 6, 7, 8, 9 (β ≈ 1). Numerical simulations also show that for q = 5, 6, 7, 8, 9, the heat conductivity of the system converges to a non-zero constant in thermal dynamical limit, while it diverges for q = 3 and decays to zero for q = 4 [140]. In the second work [141], the scatterers are triangles (Fig. 14). It was found that the transport behaviors of this model depend sensitively on whether the angles, measured in unit of π, are rational numbers or irrational numbers. For the rational cases, the conductivity diverge as power law with exponent β = 0.22 and the MSD α = 1.178. While for the irrational case, the conductivity saturates to 0.225 and the diffusivity to 0.15.
All of these works in the billiard models lead us to think about the intrinsic relation between normal (anomalous) particle diffusion and normal (anomalous) heat conduction. Consequently, there were two independent works trying to relate the diffusion exponent α to β, one by Li and Wang [54] and the other by Denisov et al. [75]. Both the studies concentrated on the gas models. Particles transport between the ends of dynamic channels and exchange energy with the heat baths connected to the ends. In both models, particles do not interact just as the billiard models described above. The first one [54] used a general continuous time random walk description. By studying the mean first passage time of the particles, Li and Wang finally obtained the relation β = 2 − 2/α.
In the second study [75], Denisov et al. focused on the spe-cial Lévy type random walk, which generally assumes that a particle waits at each point for a waiting time following distribution ψ w (t) ∼ t −γw−1 (γ w > 1) and then make a jump with distance following distribution ψ f (x, t). For Lévy walk [142], ψ f (t) ∼ t −γt−1 and ψ f (x, t) ∼ ψ f (t)δ(|x| − vt). This means that during a single flight, the particles travel at a constant speed v. Depending on γ t , the particles diffuse according to x 2 (t) ∼ t α with α = 2 for 0 < γ f < 1; α = 3 − γ f for 1 < γ f < 2 and α = 1 for γ f > 2. For the superdiffusive cases 1 < γ f < 2, the authors obtained the mean first passage time τ ∝ L γ f . The heat flux contributed from a single particle then scales as J ∝ L −γ f = L α−3 . Considering that the temperature gradient ∇T ∝ L −1 and the total number of particles in the channel N ∝ L, it was obtained that the heat conduc- However, due to the divergence of mean square displacement, |x| 2 ∼ t α was used instead of the mean square displacement For the subdiffusive cases, the conductivity vanishes because the mean passage time diverges.
Although the billiard gas models helped us clarify a few puzzles on the conditions for diffusive heat transport, we should admit that due to the lacking of interactions and local thermal equilibrium as pointed out by A. Dhar and D. Dhar [143], conclusions drawn from them could not be directly transferred to oscillating lattices. The properties of the billiard gas models are dominated by their dynamic properties, or even simpler, geometric properties. while for lattices, thermodynamic properties are more important and of interest. Moreover, heat flux in lattices takes place by phonon transport without net particle flow, which also make the lattice systems different from the billiard gases in underlying physics.
For energy diffusion in lattices, Cipriani et al. studied a one dimensional diatomic hard-point model [144]. The model consists of a chain of hard-point particles with alternating masses, m 2i = m and m 2i+1 = rm, lying on a line segment with length L. The particles move in one dimension and elastically collide with other nearest neighbors when they meet at the same point. Therefore, the order of the particles are conserved and the system is more like a lattice model compared to gas models. In their simulation, the average kinetic energy is chosen to be m i v 2 i /2 = 1 and m = 1. It is clear that when r = 1, the collisions only lead to exchange of velocities and the dynamics is integrable. When r deviates from 1, the system is non-integrable but still remains non-chaotic because the evolution equations is linear: if the velocities of two successive particles before collision are v i,j , then after collision they will change according to a linear form where plus (minus) sign in ± corresponds to i even (odd) and j = i ± 1.
After introducing an infinitesimal perturbation δv i (t = 0) to the i'th particle, the spread of the energy perturbation at later time can be characterized by δ (2) This quantity can be regarded as the energy distribution in the energy diffusion process. It is observed by numerical simulation that the profile of δ (2) (i, t) (See Fig. 15) satisfies the power-law ansatz δ (2) (i, t) ≈ t −γ δ (2) (i/t γ ) very well. A best fit for i = 0 gives γ = 0.606 ± 0.008. The evolution of the profile δ (2) (i, t) was then compared with Lévy walk with waiting time distribution ψ(t) = t −µ−1 [142]. The scaling law of Lévy walk is [145] P (x, t) = 1 t 1/µ P ( for the central part. Therefore, µ = 1/γ ≈ 5/3. Numerical simulations showed a good agreement between this two models except that the ballistic wavefronts are slightly broadened (Fig. 15). Since for Lévy walk, the mean square displacement x 2 (t) ∝ t α with α = 3 − µ, therefore α should be 4/3 approximately. Numerical calculation of mean square displacement of their diatomic gas model showed that β ≈ 1.35, which supported the connection of their model with Lévy walk.
The heat conduction of this model had been separately studied earlier by Hatano [5], Dhar [6] and by Grassberger et al. [7]. The first work [5] used r = 1.22 and found that β = 0.33 ∼ 0.37. The second one [6] observed a small divergence of the heat conductivity with β ≈ 0.17 when r = 1.22 using length N up to 1281. The last one [7] used a very efficient event driven algorithm and simulated chains with length up to 16383. It was found that for very large N , the divergent exponent tends to β = 0.32 +0.03 −0.01 for all mass ratio between 1 and ∼ 5, which was claimed to be consistent with the prediction of renormalization group theory, β = 1/3, by Narayan and Ramaswamy [45].
The hard-point gas model cannot mimic the real lattice models because the interactions between particles are present only when they collide. For interacting lattices connected by springs, another important systematic study on energy diffusion were carried out by Zhao [146]. In this work, Zhao showed a heuristic approach to investigate the energy distribution in a diffusion process, denoted as ρ(i, t), using equilibrium energy spatiotempral correlations, i.e.
This quantity is in equivalence with δ (2) (i, t) in Cipriani et al.'s work [144]. Using this method, Zhao studied three well known lattices showing different behaviors of heat conduction: the Toda lattice which is ballistic, the FPU-β lattice which is superdiffusive and the φ 4 lattice which is diffusive. It was found that for all these models, the energy correlation functions can well capture the characteristics of the corresponding energy diffusion (Fig. 16). Specifically, for Toda lattice, the energy distribution is dominated by two clear ballistic wave front; for φ 4 lattice, the distribution is Gaussian-like; for FPU-β lattice, it is a combination of the two above, two seemingly ballistic peak and a Gaussian-like central part. Quantitative results for the mean square displacement x 2 (t) = i i 2 ρ(i, t) were calculated, which shows power law dependence on time x 2 (t) ∼ t α . The exponent α ≈ 2 for Toda lattice, α ≈ 1.4 for FPU-β lattice and α ≈ 1 for φ 4 lattice, which coincide with their ballistic, superdiffusive and diffusive behaviors, respectively. Moreover, it was argued that this supports the β = α − 1 relation because the exponent β for FPU-β model is 2/5 as predicted by the Peierls-Boltzmann theory.
In a recent work, Zaburdaev et al. restudied the energy diffusion in one dimensional systems using the phenomenological Lévy walk approach [147]. They first generalized the Lévy walk model to involve random fluctuation in the movement of particles. Specifically, compared to the standard Lévy walk, it was additionally assumed that during a single flight between two successive collisions, a particle's position evolves according to Langevin equationẋ(t) = v 0 + ξ(t), where ξ(t) is a Gaussian random process with correlation ξ(t)ξ(s) = D v δ(t − s). When the randomness is absent (D v = 0), this model reduces to the standard Lévy walk model with moving speed v 0 . Analytic analysis showed that for their generalized Lévy walk model, the central part of the particle density profile P (x, t) follows the scaling of standard Lévy walk, i.e., Eq. (15), while the two ballistic humps follow another scaling wherex = x − v 0 t. Using this model, Zaburdaev et al. restudied the energy diffusion in hard-point particles (HPG) model [5][6][7]144] and the FPU-β model using Eq. (16) by Zhao [146]. For the hard-point particles model, Zaburdaev et al. calculated the evolution of an infinitesimal perturbation ρ(i, t) following Ref. [144]. It was found that the scaling relations Eq. (15) with γ = 5/3 and Eq. (17) are fulfilled separately for the central part and the ballistic part, respectively (Fig. 17). For the FPU-β model, the authors performed a direct comparison between the energy spatiotempral correlation functions e(i, t) as in Eq. (16) and the probability density for the generalized Lévy walk model. The comparison was made at time t = 1000 and 2000 respectively (Fig. 17). A good agreement was observed by choosing the average energy per particle ǫ = 1, v 0 = 1.384 and D v = 0.49.
Although there are a few numerical studies showing that energy diffusion in some one dimensional systems can be described by (generalized) Lévy walk, it can be argued that a fundamental theory connecting these two phenomena does not exist in general. All existing works are based on phenomenological approaches. It is not known either for which kind of systems, the Lévy walk description is applicable.
However, regarding the connection between (anomalous) heat conduction and (anomalous) energy diffusion, in a most recent work, Liu et al. have rigorously proved an equality that can be generally applied to all homogeneous one dimensional systems near thermal equilibrium [148]. In their work, Liu et al. considered a general one dimensional Hamiltonian system which was initially (t = −∞) at thermal equilibrium with temperature T . Then by assuming local thermal equilibrium, a small perturbation with perturbative Hamiltonian H ext = − 1 T ǫ(x, t)∆T (x)dx was imposed on the system from t = −∞, where ǫ(x, t) is the local energy density. This form of perturbation will drive the system away from equilibrium and finally reach a steady state with temperature profile T (x) = T + ∆T (x), provided that local thermal equilibrium is satisfied. After removing the perturbation at t = 0, the system will start to relax to equilibrium at temperature T again. This relaxation process is interpreted as energy diffusion. In their work, Liu et al. first analytically proved that in the diffusion process (t ≥ 0), the normalized distribution of the energy profile can be written as where C ǫǫ = ∆ǫ(x, t)∆ǫ(0, 0) is the autocorrelation function of energy density fluctuation. Consequently, it was proved that the mean square displacement of energy diffusion should satisfy a very general equality independent of the system Hamiltonian and the initial temperature profile ∆T (x), which reads where C(t) is again the heat flux autocorrelation function appearing in the Green-Kubo formula Eq. (8). Because in their derivation there is no any explicit assumption except the local thermal equilibrium, this equality should be universal, even for isotropic systems in higher dimensions.
The authors also showed a few corollaries of this equality for different diffusive behaviors. When the system is normal diffusive, i.e. x 2 (t) ∼ 2Dt in large t limit, by inserting Eq. (19) into the Green-Kubo formula, we will get κ = c V D. For the superdiffusive case x 2 (t) ∼ t α with 1 < α < 2, the Green-Kubo formula with a cut-off time t c ∝ L will lead to κ ∼ L α−1 , which means β = α − 1. While for subdiffusive cases (0 < α < 1), the heat conductivity decays to 0 as κ ∼ 1/L 1−α [148].

CONCLUSIONS
In this Colloquium we took the reader on a tour presenting the state of the art of the topic about anomalous thermal conduction and heat diffusion in low dimensional nanoscale structures, an emerging research direction which is expected to shade light for renewable energy and thermal management for future nano devices. Particularly, we surveyed and explained physical mechanisms that respond to the anomalous thermal conductivity.
We firstly reviewed recent experiments on thermal conductivity of nanotube, nanowire and graphene. Then, we provided a review on the computational investigations of thermal conduction in nanascale, concentrating on atomistic simulations. At last, various theoretical models to explain the anomalous thermal conductivity and the connection between thermal conduction and heat diffusion were overviewed. Both numerical and experimental studies in nanostructures ranging from nanotubes [8,10,51] nanowires [9], to polyethylene nanofibers [17,26,149] have shown that heat conduction in low-dimensional systems does not obey Fourier's law even though the system length is much longer than the phonon mean free path. In these structures phonons transports superdiffusively-a process faster than random walk but slower than ballistic motion, which leads to a length dependent thermal conductivity. Thus nano structured materials are promising platforms to testify fundamental phonon transport theories.
Moreover, given the fact that thermal management devices generally can be considered for applications that require power ranging from milli-watts up to several thousand watts, the general range of applications where they are indispensable is stupendous. Compared with decades ago, we have a better and clearer understanding of the thermal conduction and heat diffusion in nanoscale materials.

OUTLOOK AND CHALLENGES
However, there are still many open questions. On the experimental side, experimental demonstration of low-dimensional phonon transport theory depends on how to measure its thermal conductivity accurately. To this end, one needs to make the sample smaller, for example towards a dozens of nanometers or even a few nanometers only. Difficulties in integrating nanomaterials with suspended structures suitable for thermal measurements have resulted in rare experimental studies of size-dependent thermal conductivity in 1D and 2D systems. And another challenge is to get rid of the contact thermal resistance. The unknown thermal contact resistance remains the primary technical challenge in testing the Fourier's law. More efforts are needed to explore a brand new method of measuring thermal contact resistance directly.
Because of the complexities inherent in thermal contact resistance measurement, this task needs the help from theoretical predictions. However, so far a truly comprehensive theory for this effect is lacking. The current approaches for thermal transport across an interface, such as the acoustic mismatch (AMM) theory and the diffusive mismatch (DMM) theory, are based on the assumption that phonon transport via a combination of either ballistic or diffusive transport on either side of the interface. Both schemes offer limited accuracy for nanoscale interfacial resistance predictions, due to as we have shown above, phonons transport super-diffusively in low-dimensional systems. Therefore, it is necessary to establish an improved theory describing thermal transport across the interface by taking into account the anomalous thermal transport characteristics of nanostructures. To provide quantitative prediction from theory and simulation, how to set up a transport theory by incorporating nonlinearity in the quantum regime is still a challenge. The nonequilibrium Green's (NEGF) function method [150] certainly serves as an elegant mathematical framework. However, when including the phonon-phonon interactions, the NEGF presents a cumbersome task. Thus further systematical investigations combine experimental and theoretical efforts on fundamental mechanisms will be helpful to advance the field.