Thermal Transport in Two-Dimensional Heterostructures

Heterostructures based on two-dimensional (2D) materials have attracted intense attention in recent decades due to their unusual and tunable physics/chemical properties, which can be converted into promising engineering applications ranging from electronics, photonics, and phononics to energy recovery. A fundamental understanding of thermal transport in 2D heterostructures is crucial importance for developing micro-nano devices based on them. In this review, we summarized the recent advances of thermal transport in 2D heterostructures. Firstly, we introduced diverse theoretical approaches and experimental techniques for thermal transport in low-dimensional materials. Then we briefly reviewed the thermal properties of various 2D single-phase materials beyond graphene such as hexagonal boron nitride (h-BN), phosphorene, transition metal dichalcogenides (TMDs) and borophene, and emphatically discussed various influencing factors including structural defects, mechanical strain, and substrate interactions. Moreover, we highlighted thermal conduction control in tailored nanosystems—2D heterostructures and presented the associated underlying physical mechanisms, especially interface-modulated phonon dynamics. Finally, we outline their significant applications in advanced thermal management and thermoelectrics conversion, and discuss a number of open problems on thermal transport in 2D heterostructures.


INTRODUCTION
The excellent physicochemical properties of graphene (superior thermal stability, high conductivity, favorable mechanical strength, broadband absorption) have inspired great interest not only in its fundamental and applied aspects, but also in exploring other two-dimensional (2D) materials. Subsequently, numerous post-graphene 2D materials ranging from hexagonal boron nitride (h-BN) , silicene (Vogt et al., 2012), transition-metal dichalcogenide (TMDC) like MoS 2 (Smithe et al., 2018) and MoSe 2 , to phospherene (Chen et al., 2018c) and borophene (Ranjan et al., 2019) have been synthesized, and present diversified electronic, thermal and optical properties, which could be complementary to those of graphene. For example, compared to the zero band-gap of graphene, MoS 2 exhibit direct band gap, making it an attractive material for nanoelectronic applications like field-effect transistors with a high on-off ratio and low power consumption . Another example, stanene, a single-layer buckled honeycomb structure of Tin atom, is demonstrated to be near-room-temperature quantum anomalous Hall effect  and ultra-low thermal conductivity (Cherukara et al., 2016), which is particularly suited for thermoelectric devices. Furthermore, 2D organic carbon nitride materials as graphene derivatives exhibited high photocatalytic activity owing to their interface interaction, and can be used as the next-generation highly effective photocatalyst (Zhao et al., 2018a;Zhao et al., 2021).
Recently, research on heterostructures composed of distinct 2D materials has captured much attention because their emergence could bring novel exciting physical properties beyond either of the individual component monolayer structures and provide additional degrees of freedom for device engineering. Normally, heterostructures are divided into two classes: in-plane heterostructures where two 2D materials are stitched seamlessly into a shared plane with an atomically sharp 1D interface; van der Waals (vdW) heterostructures where distinct 2D materials are stacked vertically with non-bonded vdW forces between adjacent layers. Benefited from the progresses in the microtechnology, high-quality graphene/ h-BN in-plane heterostructures can be synthesized via chemical vapor deposition (CVD) method, and their compositional and structural diversities could translate into greater freedom for tuning the physical properties such as the presence of a metal-insulator transition (Ci et al., 2010;Liu et al., 2014d). Based on the same preparation method, the integration of MoS 2 , MoSe 2 or WSe 2 monolayers into a heterostructure has provided the possibility to produce strong localized photoluminescence and in-plane p-n junctions, showing perspectives in phototransistors and field-effect transistor (Huang et al., 2014;Li et al., 2015). The advantages of inplane heterostructures mainly contain simpler band alignment, more distinct phase separation and strong interfacial action . Compared to in-plane heterostructures which demands the matched lattice structures of constituent materials, vdW heterostructures intrinsically create an ultraclean interface regardless of the typical interfacial latticematching constraints, which is crucial for tunneling devices that suffer from interfacial defects and dislocations. Additionally, vdW heterostructures offer new mechanisms for manipulating the properties of emergent device via the stacking sequence, the relative rotation between adjacent layers as well as interlayer spacing. A few pioneering works have reported prototypical field-effect-transistor with graphene-based (Bertolazzi et al., 2013;Georgiou et al., 2013) or phosphorus vdW heterostructures (Kang et al., 2016), photodiodes based on MoS 2 /WSe 2 stacks (Roy et al., 2015), and high-efficiency photovoltaic cells based on graphene/WS 2 vertical junctions (Iannaccone et al., 2018). On the other hand, the non-bonded vdW interactions of adjacent layers can preserve the electronic properties of each layer, and help resolve a critical challenge that the poor electronic quality of 2D layered materials deposited on SiO 2 substrates (Yankowitz et al., 2019). For instance, the quality of graphene encapsulated in h-BN is dramatically improved, and possesses high mobility and ballistic transport over tens of micrometers, which approaches those of suspended graphene (Wang et al., 2013a). It is reported h-BN/phosphorus/h-BN sandwiched heterostructure exists high field-effect mobility and high on-off ratios at room temperature, and the embed phosphorus remains intact and preserves high quality (Chen et al., 2015b). In a word, the physicochemical properties of different heterostructures are closely linked to material composition and structural configuration, representing the significance of structure-property relations.
Thermal transport in 2D heterostructures is not only of a fundamental issue in physics but also crucial to engineering applications. On the one hand, heat dissipation arising from charge carrier scattering processes has to be considered in the development of such nanostructures for electronic and photoelectric devices because the generation of waste Joule heat often leads to the formation of high-temperature hot spots in devices (Wang et al., 2017d), which severely affect the stability and performance of devices. Hence, the choice of 2D building blocks with high thermal conductivity helps to alleviate Joule heat dissipation, and a deeper understanding of the underlying mechanisms for heat conduction across the heterogeneous interface would aid in the design of functional devices. On the other hand, low thermal conductivity materials are necessary for the application in thermoelectric devices converting waste heat into electrical energy. The introduction of 2D heterostructures probably causes increasing Seebeck coefficient and lowering thermal conductivity based on carrierenergy filtering or quantum confinement effects. Compared to dissipating heat and thermoelectric generator, the design of functional thermal devices such as rectifier, transistor, and logical gate to handle heat signals requires controllable thermal transport. In fact, the structural diversities combined with structural defect and external stress field can realize the rich diversity of thermal properties in 2D heterostructures, and enables them to be used in different thermal fields. From another perspective, the interfaces of heterostructures make the thermal transport mechanisms deviate from that in single 2D materials and conventional bulk materials, either through inducing a contact thermal resistance or new phonon modes. For example, it is reported that the thermal conductivity of graphene and graphene nanoribbons (GNRs) monotonously decreases with the tensile strain due to phonon softening (Wei et al., 2011;Shen et al., 2014;Guo et al., 2020), while the load of tensile strain of graphene/h-BN in-plane heterostructures significantly enhance the interfacial thermal conductance owing to improving the alignment of out-of-plane phonon modes near the interface (Ong et al., 2016). Further research demonstrates that the thermal conductivity of graphene/h-BN superlattices first increases and then decreases with increasing the tensile strain, which stems from the competition between in-plane phonons softening and flexural phonons stiffening (Zhu and Ertekin, 2015). Furthermore, some abnormal thermal phenomena such as thermal rectification (TR) (Yuan et al., 2015;Zhang et al., 2018), negative differential thermal resistance (NDTR) , thermal bandgap (Zhu and Ertekin, 2014b) and thermophoresis (Mandelli et al., 2017) have been reported in 2D heterostructures, and sparked heated debate.
As we know, phonons are the main heat carriers in semiconductors like TMDCs and silicene, and insulators like h-BN. Specifically, the ultra-high thermal conductivity of graphene dominated by phonons has raised the exciting prospect for heat removal in electronic devices and peculiar phonon transport properties, as reviewed by (Nika and Balandin, 2012;Xu et al., 2014c;Al and Farías, 2016;Xu et al., 2016b;Nika and Balandin, 2017) which has also laid the solid basis for other emerging 2D materials. Phonon transport in novel 2D materials and their heterostuctures probably differs from that in graphene because of crystal structures and atom species, suggesting that the existing knowledge of graphene couldn't be transplanted to them directly. Fortunately, sustaining advances of theory and experiment offers microscopic details from the atomistic level to the device scale. To date, the thermal properties of 2D materials beyond graphene have been discussed in many previous review articles. Regarding to the progress of theoretical approaches and experimental techniques in 2D materials, please refer to (Cahill et al., 2014;Wang et al., 2017d;Zeng et al., 2018b;Fu et al., 2019) For the phonon transport in phosphorene, borophenes, TMDs and other 2D materials, please refer to (Gu et al., 2018;Qin and Hu, 2018;Chen and Chen, 2020a;Li et al., 2020a;Zhao et al., 2020b) For a detailed survey of thermoelectricity applications of 2D materials, please refer to Zhang and Zhang, 2017a;Li et al., 2020b). In this review, we not only summarizes the existing results of phonon transport in pristine 2D materials and their heterostuctures, but also gives a physical picture of various phonon scattering mechanisms induced by interior imperfections or external field, to pave the way for better thermal management in practical devices. In Theoretical Approaches and Experimental Measurements, we introduced diverse theoretical approaches and experimental techniques for investigating heat conduction in low-dimensional nanostructures. In Thermal Transport in Some Typical 2D Materials, we briefly reviewed the unique thermal properties of several typical 2D materials. In Thermal Transport in 2D Heterostructures, we highlighted the manipulation of heat conduction in tailored nanosystems-2D heterostructures and present the associated underlying transport mechanisms, especially interface-modulated phonon dynamics. Finally, we show the potential applicability of 2D heterostructures in advanced thermal and thermoelectrics devices, and give the conclusions and a brief outlook.

Theoretical Methods
Currently, different theoretical method have been developed to study the thermal transport properties of 2D materials, principally including the Boltzmann transport equation (BTE), atomistic Green's functions (AGF) and molecular dynamics (MD) simulations. Note that the latter two are especially suitable for describing phonon transmission/scattering across an interface in 2D heterostructures.

Molecular Dynamics Simulations
MD simulations with classical potentials are a powerful technique to handle many-body problems at the atomistic level. Here, allorder anharmonic effect is inherently contained in atomistic interactions, differing from that in BTE and AGF. At present, the MD method is widely adopted to study phonon transport problems such as structural effect and various influences like doping, defect, strain and substate. Note that the accuracy of classical MD calculations depends on the quality of the interatomic potentials. Recently, intensive efforts have been conducted to construct reliable potential function from the calculations of first principles (Rowe et al., 2018). In practice, equilibrium MD (EMD) and nonequilibrium MD (NEMD) simulations are adopted to extract the thermal properties of 2D materials. In EMD, the calculation of thermal conductivity is based on the Green-Kubo formula or Einstein Relation ration (Sellan et al., 2010;Sevik et al., 2011). For NEMD, both sides of sample are connected to a hot and cold reservoir, respectively. After sufficient run-time, the stationary-state heat flux and temperature are obtained to calculate thermal conductivity in terms of Fourier's law of heat conduction, Obviously, the obtained thermal conductivity using the two methods does not provide microscopic perspective into thermal conduction in crystals. For this, a spectral energy density (SED) technique on the base of harmonic lattice dynamics and EMD calculations has been developed to extract mode-resolved thermal conductivity (Thomas et al., 2010;Qiu and Ruan, 2012a). However, the contribution of some long wavelengths phonons can't be contained in the simulation results because of size effect Gu et al., 2018). As a result, the calculated thermal conductivity is sensitive to the size of simulation domain. Note that the NEMD method has also been applied to study interface-modulated phonon dynamics such as graphene/MoS 2 , graphene/ phosphorene in-plane heterostructures  and graphene/h-BN (Xu et al., 2014a) MoS 2 /WSe 2 , phosphorene/graphene vdW heterostructures .
The presence of phonon-interface scattering will produce thermal resistance, also called as Kapitza thermal resistance. The phenomenon could be partly explained by the acoustic mismatch model (AMM) and the diffuse mismatch model (DMM), which is good for some simple material (Kakodkar and Feser, 2017). While for 2D heterostructrues, both kinds of theoretical models are likely to fail because of the dimensionality mismatch and limited thickness of actual system. In addition to this, the models couldn't give a detailed picture for interfacial thermal transport as AGF calculation. For this, a spectral decomposition of interfacial thermal conductance method based on NEMD simulations is developed by (Sääskilahti et al., 2015) and extended by (Zhou and Hu, 2017c), offering a detailed description of inelastic effects. Under nonequilibrium steady state, the spectrally resolved heat flux between two atoms i and j is given by where t s is the simulation time, v α i (ω) is the discrete Fourier transforms of atomic velocity v iα (t) and K αβ ij is the two-order force constant matrix calculated by the finite displacement. "L" and "R" denotes the left and right sides along the interface. The phonon transmission function can be obained in terms of the spectral heat flux where k B is the Boltzmann parameter, and ΔT is the applied temperature bias. At low temperatures, inelastic effects caused by the anharmonicity were negligible and the spectral conductance matched the results obtained from the harmonic Green's function method. Under this case, this transmission function is simply equal to M(ω), the number of propagating modes. The anharmonic phonon effect incorporated in NEMD render T(ω) dependent on the sample length L, which can be phenomenologically taken into account through the relation (Sääskilahti et al., 2015).
where Λ(ω) is the effective phonon mean free path (MFP). Additionally, the phonon wave-packet (PWP) method based on wavelet transforms and MD simulation has also developed to calculate the phonon transmission coefficient. A detailed derivation process is referred to the articles by (Schelling and Phillpot, 2003). The advantage of the PWP method is that the phonon scattering process could be observed directly from the simulation, not just the transmission result as the AGF. The basic idea of PWP is to form a phonon wave packet so that it is localized in both real space and reciprocal space. To produce a wave packet at wave vector q, the initial displacement from equilibrium of the ith atom in the nth unit cell along the a direction is where A m is the wave amplitude, m i is the mass of the ith atom, ε ]iα (q) is the polarization vector, z 0 and z n are the positions of the packet center. It is worth noting that because of the very small amplitude of the incoming wave packets, only elastic scattering is captured, thus the frequencies of incoming and outgoing phonons are always the same. Inelastic scattering is very difficult to pick out the signal of wave-packet against the background of the thermal vibrations.

Atomistic Green's Functions
As we know, phonons have wave nature. The AGF method can rigorously include wave effects on a discrete atomic lattice. Actually, this method was initially developed to handle quantum electron transport in nanostructures (Fan and Chen, 2010;Chen et al., 2018a;Deng et al., 2018;Zeng et al., 2018a;Chen et al., 2019a;Zhou et al., 2019;Fan et al., 2020). By making a few careful substitutions, the approach has been extended to study phonon transport in a wide variety of nanostructures (Xu et al., 2009;Peng and Chen, 2014;Zhou et al., 2017b;He et al., 2019) especially suitable for low-dimensional hetrostuctures such as Si/Ge (Tian et al., 2012b), graphene/h-BN (Peng et al., 2017), MoS 2 /metal (Yan et al., 2016b) interfaces and others (Ma et al., 2016). Here, we simply introduce this method from the numerical simulation point of view. The phonon Hamiltonian of entire system is consisting of three parts, namely Left (L), Center (C) and Right (R) regions, the matrix form can be written as where u α is a column vector consisting of all the atomic displacement in region, H n is nonlinear part of the interaction, V CR is the coupling matrix between the central region and right lead.
Here K is the force constant matrix obtained by empirical force field or the first principle calculations. Within the harmonic approximation (H n 0), the retarded/advanced Green's function for central region is given by Where I C is an identity matrix, and the Π L(R) (ω) H CL(CR) g r(a) L(R) H LC(RC) is the self-energy of left and right leads. The retarded/advanced surface Green's functions for the leads is defined as The surface Green's function usually solved recursively or using the decimation technique. According to Caroli formula, the phonon transmission coefficient as phonon frequency is represented as In terms of the Landauer formula under ballistic limit, the phonon thermal conductance is obtained from: Here f is the Bose-Einstein distribution, ω max is the maximum phonon frequency. Note that the AGF method couldn't directly provide the mode-resolved thermal condcutance. Even so, it is superior to the elastic wave continuum model, which is normally used at low temperatures (Peng et al., 2010;Xie et al., 2013;Xie et al., 2014b;Zhang et al., 2014b). Recently, a convenience technique developed by (Ong et al., 2016) has been widely used for the investigation of interface thermal transport, where the phonon reflection, absorption and transmission could be obtained, as described next. The first step is defined a non-Hermitian matrix that acts on the eigenstates of left or right lead. Here, take the left lead for example: Diagonalizing F r(a) L one can obtain the eigenvector matrix U r L ( ± ) and the eigenvalue matrix Λ r L ( ± ) that consist of the normalized eigenstates and eigenvalue. These matrices are needs for calculation of individual phonon transmission. Next, the transmission matrix is computed by: Once these matrices are known, the transmission coefficient of individual phonon mode on the left or right leads is given by: By summing over the subscript m and n, the total phonon transmission in Eq. 9 can be obtained. If there is significant nonlinear effect in the central scattering region, the third-order force constant can be consierded as a perturbation of harmonic system. The three-phonon scattering self-energy s Π is computed by using the following relations: where H(x) represents the Hilbert transform, V (3) ijk z 3 V/zu i zu j zu k is the third-order force constant. Finally, the effective transmission function is computed by: Boltzmann Transport Equation From the particle-based picture perspective, phonon transport is governed by the BTE, where phonons obey the Bose-Einstein distribution function n 0 λ at thermal equilibrium but deviating from n 0 λ under a temperature gradient. Here λ indicates the phonon mode of wave vector q and phonon polarization s. The deviating distribution function can be obtained from the BTE: where v λ is the group velocity of phonon mode. Assuming a small enough ∇T in most practical situations, n λ can be expanded to in first order so that n λ n 0 λ + Δn λ , where Δn λ depends on ∇T linearly. If only the scattering mechanism of two-and threephonon processes are considered, Eq. 1 is written as: where F λ τ 0 λ (υ λ + Δ λ ) and τ 0 λ is the relaxation time of mode λ, which corresponds to the relaxation time approximation (RTA).
where Γ + λλ′λ ″ denotes absorption processes that the combined energy of two incident phonons leads to one new phonon (ω λ′ + ω λ″ ω λ ), and Γ − λλ′λ ″ describes emission processes in which the energy of one incident phonon is split among two phonons (ω λ′ + ω λ ω λ″ ). A detailed descriptions of Γ and Δ λ is available in . Eq. 2 is numerically solved by iterating from a zero-th order approximation start of F λ τ 0 λ υ λ . The iterative procedure has a great influence on diamond, graphene and carbon nantubes (CNTs) in which the normal processes play a significant role in the phonon-phonon scattering (Lindsay et al., 2009;Lindsay et al., 2010;Lee et al., 2015a;Lee and Lindsay, 2017). However, for Si, Ge and amorphous materials with strong Umklapp scattering (Ward and Broido, 2010;Zhou et al., 2020a) iterating to convergence typically has little effect on the obtained thermal conductivity. Apart from phonon-phonon scattering, there may be various other phonon scattering processes, including phonon-defect (τ λ,pd ), phonon-impurity (τ λ,pd ), and phonon-boundary (τ λ,pb ). Based on Matthiessen's rule where various scattering mechanisms are presumed to be independent of each other, the total relaxation time is given by τ λ,t − 1 τ λ − 1 + τ λ,pd − 1 + τ λ,pi − 1 + τ λ,pb − 1 . By the above solution, the lattice thermal conductivity κ can be obtained here C λ Zωzn 0 /zT is the specific heat of the mode λ. At present, the BTE in combination with first-principles calculations becomes popular to predict the thermal conductivity of various 2D materials, showing reasonable agreement with experimental measurements (Lindsay and Broido, 2011;Fugallo et al., 2014;Gu and Yang, 2014;Lindsay et al., 2014). It is writhing pointed out that the over-simplified of phonon scattering processes and a finite number of k points causes obvious ambiguity in the calculation of thermal conductivity. Another point worth noting is that the coupling effect of various scattering mechanisms play important roles in phonon transport.

Experimental Measurement
The measurement of thermal conductivity or interfacial thermal conductance (ITC) in 2D nanostructures is helpful to grasp the law of phonon transport. At present, the optothermal Raman and microbridge techniques are the most common method for measuring the thermal properties of 2D materials. In addition the 3ω and time-domain thermoreflectance (TDTR) techniques have been adopted to obtain ITC. Because of space limitation, we briefly introduce some measurement results for 2D materials and discuss the measurement errors.
With the aid of optothermal Raman technique, the thermal conductivity measurement of single-layer graphene was reported by (Balandin et al., 2008), and the obtained value ranges from 2,000 to 5,000 Wm −1 K −1 at room temperature. Analogously, the thermal conductivity of bilayer graphene was also identified , coinciding with the theoretical predictions (Nika et al., 2014). Furthermore (Malekpour et al., 2016), reported that the room-temperature thermal conductivity of defective graphene reduces from 1800 to 400 Wm −1 K −1 when the defect density induced by electron beam irradiation slightly increases from 0.0005 to 0.005%.  measured the thermal conductivity of graphene with different 13 C isotope concentrations, and revealed that 50% concentration results in a 50% reduction in thermal conductivity. Except graphene, this technique has also applied to black phosphorous  MoS 2 and MoSe 2 (Zhang et al., 2015c), In a word, although the uncertainty of the optical absorption coefficient, weak temperature-dependence of Raman peak shift and uncertainty of contact thermal resistance will affect the measurement accuracy, the optothermal Raman measurement is still a simple and effective method to measure the thermal conductivity of 2D structures.
The microbridge technique was initially used to measure the thermal conductivity of 1D nanostructures (Li et al., 2003), and then was applied in 2D materials (Pettes et al., 2011). For example, (Xu et al., 2014b) experimentally found the logarithmic dependence of thermal conductivity on the sample length L (κ ∼ log L) in suspended graphene. Moreover, (Jang et al., 2013) reported that there is positive correlation between thermal conductivity and the number of layers in few-layer graphene, and the value of ten-layer graphene is over 1,000 Wm −1 K −1 . On the other hand, the optothermal Raman measurement indicated that the thermal conductivity of 4-layer MoS 2 is 44-50 W m −1 K −1 at 300 K. (Jo et al., 2014). Analogously, (Lee et al., 2015b) verified that the thermal conductivity anisotropic ratio approaches two at 300 K. Note that the presence of contact resistance, polymer residue and inconsistency of sample qualities limit the precision of measurement.
The ITC measurement of in-plane heterostructures is quite challenging compared to that of vdW interfaces. There exist far fewer measurements on the ITC. Usually, the measured ITC of graphene/substrate interface is in the range from 20 to 100 MWm −2 K −1 at 300 K. Taking this value as reference, the measured ITC MoS 2 /substrate interface is considerably smaller. Using the TDTR-based approach, (Lee et al., 2015b) reported that the ITC of MoS 2 /metal thin-films interface is ∼26 MW m −2 K −1 . Moreover, the optothermal Raman measurement indicated that the ITC of MoS 2 /SiO 2 interface is about 1.94 MW m −2 K −1 (Taube et al., 2015). Analogously,  reported that the ITC between WSe 2 /SiO 2 interface ranges from 10 to 32 MW m −2 K −1 via analyzing the temperaturedependence of the low-frequency shear mode. ) measured the ITC of vertically stacked graphene/h-BN interface and found that the value reaches 7.4 MW/m 2 K at room temperature, which stems from the remaining organic residues at the interface.

THERMAL TRANSPORT IN SOME TYPICAL 2D MATERIALS Graphene
It is well known that single-layer graphene possesses ultra-high room-temperature thermal conductivity. The results of thermal conductivity measurements for graphene have been discussed in Molecular Dynamics Simulations. Note that the measured values are closely related to experimental approaches, external condition (strain, substrate), sample size (length, layerthickness) and quality (like point-defect, isotopic abundance, grain size). Hence, theoretical simulations could provide deeper insight into concerning the thermal property of graphene. Using the PBT method combined in conjunction with the forceconstant model,  claimed that the flexural acoustic (ZA) phonons contribute 70-80% of heat conduction in graphene, which is ascribed to the existence of selection rules in phonon scattering processes. The phenomenon has also been confirmed by others research groups based on similar theoretical approach (Fugallo et al., 2014;Lindsay et al., 2014). Using the BTE method and Raman technique, (Seol et al., 2010) reported that the thermal transport contribution from ZA modes for graphene is as large as 77%, and this ratio is decreased by the substrate or organic residues (Pettes et al., 2011). Recently, using the iterative BTE method, (Lee et al., 2015a) found that there is an obvious hydrodynamic phonon transport in graphene over a wide temperature range due to the presence of strong N scattering processes. Subsequently, (Machida et al., 2020) experimentally showed an extremely high thermal conductivity in the very thin graphite samples (8.5-micrometer-thick) at higher temperature, which extends the hydrodynamic regime from cryogenic to room temperatures. That is, there are three phonon transport behaviors, namely ballistic, diffusive and hydrodynamic transport. Furthermore, (Al and Farías, 2016) outlined the phonon dynamics in supported graphene, and revealed the physical picture of phonon-substrate scattering. (Nika and Balandin, 2017) minutely discussed the mode-resolved thermal conductivity in graphene. In addition, our previous work highlighted the manipulation of thermal transport in graphene from the viewpoint of the wave/particle nature of phonons .  summarized experimental and theoretical results about the thermoelectric properties of graphene and its derivatives, and described how the poor thermoelectric performance of graphene achieves high thermoelectric figures by tuning on a wide range both the thermal and electronic transport based on different strategies of nanostructuring and bandgap engineering. Those studies on thermal transport of graphene have stimulated great interest in exploring other emerging 2D materials, as shown in Figure 1.

Hexagonal Boron Nitride
Monolayer h-BN is also called "white graphene" because of their similar honeycomb structure. Unlike zero band-gap graphene, h-BN with 5.8 eV band-gap has been viewed as the top ranked candidate for the dielectric material in electronic devices or the new nonmetal catalyst in the field of photocatalysis (Zhao et al., 2019). It was reported that the electron mobility of graphene electronic devices supported on few-layer h-BN far outweighs that supported on silicon dioxide (SiO 2 ) because of the atomically smooth surface of h-BN (Serrano et al., 2007). In addition, h-BN also holds high chemical stability, strong mechanical strength and good thermal property. The basal-plane thermal conductivity of high-quality bulk h-BN sample could reach up to 390 W m −1 K −1 at 300 K (Geick et al., 1966), which shows great potential in current-generation dielectric materials. For sing-layer h-BN, the BTE calculations indicated that the system's thermal conductivity is about 600 Wm −1 K −1 (Lindsay and Broido, 2011), but the EMD results showed that the value is 400 Wm −1 K −1 (Sevik et al., 2011).
Experimentally, because the intensity of Raman peaks in h-BN is very weak, the optothermal Raman method may not be applicable for h-BN.  found only 5% optical absorption for nine-layer h-BN in optothermal Raman measurment, and the obtained thermal conductivity was 227-280 m −1 K −1 . Using the microbridge method,  reported that the thermal conductivity of an 11-layer sample reaches ∼360 W m −1 K −1 at 300 K, approaching the value of bulk h-BN (Lindsay and Broido, 2011). Actually, the accurate measurement of thermal properties of 2D heterostructures still has a long way to go.

Silicene
Silicene possesses a buckled honeycomb atomic arrangement of Si atoms. There is a spin-orbit coupling of 1.55 meV owing to the breaking of mirror symmetry in this bucked structure. In addition, silicene can be easily integrated into electronic devices (Liu et al., 2013a). From the view of the heat dissipation of electronic devices, the influence of buckled structure on phonon transport is well worth exploring. Based on the BTE method, (Gu and Yang, 2015) demonstrated the thermal conductivity of silicene is only a fifth of that of bulk silicon (140 W m −1 K −1 ). Moreover, it was also found the contribution of ZA modes to the thermal conductivity of silicene is about 7.5%. Similarly,  found the relaxation times of acoustic phonon modes with long-wavelength are now close to zero at q → 0, accounting for the low thermal conductivity for silicene.
Based on different MD simulations, the calculated thermal conductivity for silicene ranged from 5 to 50 W m −1 K −1 , which could be ascribed to the adopted empirical potentials including the Tersoff potential, embedded-atom method (MEAM) potential, optimized Stillinger-Weber potential Liu et al., 2014b). Compared with planar 2D materials, the structure buckling has important effect on the phonon transport in silicene via providing additional scattering or breaking the symmetry selection rule. As a result, there is an unusual thermal response to the applied tensile strain in silicene: a small strain can enhance the thermal conductivity while further increase will cause the decrease of thermal conductivity. That is differentiated from the monotone decreasing with tensile strain for graphene (Hu et al., 2013). This can be explained as follows: for a small strain, the bucked configuration becomes flat because of bond rotation, rendering the enhancement of in-plane stiffness and facilitating phonon transport; under large strain, the Si-Si bonds are significantly elongated and the in-plane stiffness weakens, lowering the thermal conductivity. However, the experimental study of thermal transport in silicene hasn't been reported by now.

Transition Metal Dichalcogenides
Unlike graphene and h-BN with single atomic layer, a transition metal layer is sandwiched between two chalcogenide atomic layers for TMDs. In general, these TMDs possess a finite Frontiers in Materials | www.frontiersin.org December 2020 | Volume 7 | Article 578791 bandgap, and are suitable for many electronic and photonic applications (Zhao et al., 2020a). Owing to the difference of lattice structure between TMDs and graphene, TMDs demonstrates different thermal transport properties (Zhao et al., 2020b). Clearly, the strong covalent bonds between the carbon atoms result in a high thermal conductivity of graphene. By comparison, the transition metal-chalcogen bond is weak because of its special sandwich structure, raising the particularity of phonon properties. Besides the crystal structure, the atomic masses are another factor affecting the thermal conductivity of TMDs. A large mass difference between transition metal and chalcogen atoms often induces substantial phonon-gap (Jiang, 2014). For instance, the phonon gap in WS 2 achieves 3.3 THz, and some three-phonon processes (like acoustic + acoustic → optical) are strongly suppressed owing to the energy conservation requirement; while the phonon gap is only 1.4 THz in MoS 2 , far less than the frequency range (0∼6.9 THz) of acoustic branches (Gu and Yang, 2014). As a result, long phonon relaxation time is observed in WS 2 , and the thermal conductivity of WS 2 (142 Wm −1 K −1 ) is significantly greater than that of MoS 2 (103 Wm −1 K −1 ). MoS 2 and MoSe 2 as the most typical TMDs have been experimentally fabricated (Keyshar et al., 2017), and their thermal transport properties have also been widely explored. (Wei et al., 2014) claimed that the thermal conductivity of MoS 2 is about 26 Wm −1 K −1 at room temperature using BTE method under the RTA assumption. Analogously, (Gu and Yang, 2014) reported the thermal conductivity of MoSe 2 is 54 Wm −1 K −1 under the iterative solution, which is over 20% than that under the RTA assumption. Furthermore, the relative contribution of different phonon polarization branches to thermal conductivity resembles the case in silicene rather than graphene because MoS 2 is not an individual atomic plane. (Cai et al., 2014) predicted a room-temperature thermal conductivity of 23.2W Wm −1 K −1 for MoS 2 by solving the AGF. Based on NEMD simulations, (Hong et al., 2016b) showed that the thermal conductivities of MoS 2 and MoSe 2 are 110.43 and 43.88 Wm −1 K −1 , respectively. Lately, 2D MXenes, namely hexagonal transition metal carbides/nitrides, are highly attractive because of their promising technological applications such as energy storage, thermoelectric conversion and flexible devices (Ma et al., 2014;Sarikurt et al., 2018). Clearly, the performance of aforementioned devices strongly depends on the thermal transport properties of constituent material. As we know, low thermal conductivity and excellent conductor are desired in thermoelectric devices, while micro-nano electronic devices require the material with high thermal conductivity. Some previous studies demonstrated that the heat dissipation is a major challenge for MXene-based devices owing to their large aspect ratios and heterogeneous structures. For example,  found the ITC of Ti 3 C 2 T z /SiO 2 interface is the range of 10-27 MW m −2 K −1 under different annealing temperatures, and a similar range of ITC was obtained in Ti 3 C 2 /SiO 2 interface with or without AlO x encapsulation . On the other hand, the in-plane thermal conductivity of MXenes has also been studied. For instance, (Guo et al., 2018) demonstrated that the thermal transport performance in luorine-terminated Ti 2 CT z exceeds that in oxygen-terminated one, but is comparable to that of bare Ti 2 C using BTE method. Analogously, (Gholivand et al., 2019) reported that compared to oxygen-terminated Ti 3 C 2 T z , the fluorine-terminated structure has remarkable advantage in heat dissipation. In another theoretical study, (Sarikurt et al., 2018) investigated the influences of transition metals and surface functionalization on the thermal conductivity of M 2 CO 2 , and the obtained results showed that their thermal conductivities range from 15-60 Wm −1 K −1 .
Experimentally, the optothermal Raman method was used to measure the room-temperature thermal conductivities of 1-layer and 11-layer MoS 2 , and the obtained values were ∼34.5 Wm −1 K −1  and ∼52 Wm −1 K −1 (Sahoo et al., 2013), respectively, which exceeds the results calculated from the NEMD method . Meanwhile, the measured value for 4-layer MoS 2 was in the range of 44∼50 Wm −1 K −1 by thermal bridge method (Jo et al., 2014). In addition, the temperature-dependent thermal conductivity of substrate-supported MoS 2 was extracted from the optothermal Raman method, and the measured room-temperature thermal conductivity was 62.2 Wm −1 K −1 , but reduced to 7.45 W Wm −1 K −1 at 450 K (Taube et al., 2015). Apart from MoS 2 , (Peimyoo et al., 2015) used the optothermal Raman method to determine the thermal conductivity of monolayer and bilayer CVD-grown WS 2 . The measured values are 32 and 53 Wm −1 K −1 for monolayer and bilayer WS 2 at 300 K, respectively. Analogously, (Zhang et al., 2015c) reported that the room-temperature thermal conductivity of monolayer MoSe 2 is 59 ± 18 Wm −1 K −1 . The similar phonon and Raman scattering in different kinds of TMDs have been discussed earlier. It is believed that other than phonon frequency due to the different mass between Mo and W atom, the common features of 1 L, 2 L and bulk structures in both MoS 2 and WS 2 should be similar. ) demonstrated a low thermal conductivity of 9 W m −1 K −1 for TaSe 2 with a thickness of 45 nm. The experimental discrepancies were primarily due to the sample quality like surface disorders, and the measurement error. Most recently, (Liu and Li, 2018c) experimentally observed that the thermal conductivity of Ti 3 C 2 T z to be 55.88 Wm −1 K −1 , Note that neither explicit configuration nor surface termination atoms are clarified in this study.

Phosphorene
Phosphorene as a 2D puckered honeycomb structure possesses high carrier mobility and large direct band-gap (∼1.5 eV) (Liu et al., 2014c), indicating abundant nanoelectronic applications. In addition, phosphorene has shown promising application in thermoelectrics devices (Fei et al., 2014). Actually, either electrical or thermoelectrical devices strong depend on the thermal transport properties of phosphorene. Hence, the investigation of the thermal conductivity of such 2D material is necessary.
Theoretically, the thermal conductivity of phosphorene has been explored using different methods including EMD (Xu et al., 2015), BTE under RAT (Zhu et al., 2014a) or iterative rule . For example,  showed that the simulated thermal conductivity along zigzag (ZZ) and armchair Frontiers in Materials | www.frontiersin.org December 2020 | Volume 7 | Article 578791 (AC) directions are respectively 15.33 Wm −1 K −1 and 4.59 Wm −1 K −1 using the iterative BTE method. The low thermal conductivity could be understood from the softening transverse optical (TO) phonons resulted from the long-ranged interaction in phosphorene structure. Combined the spectral energy density (SED) technique with EMD simulations, the ZZ and AC -direction thermal conductivity were predicted to be 152.7 and 33 Wm −1 K −1 in phosphorene, respectively (Xu et al., 2015). (Zhu et al., 2014a) used the RTA to predict the thermal conductivity of black phosphorene of 84 Wm −1 K −1 (24 Wm −1 K −1 ) in the ZZ (AC) direction (armchair) at 300 K.
Although the results reported are discrepant from different research groups, the thermal transport anisotropy behaviors have been clearly demonstrated for all techniques. To understand the phenomenon, the thermal conductivity contribution from ZA modes has also been calculated, which corresponds to 16% (8%) for the case of the ZZ (AC) direction . That is analogous to silicene (7.5%) but different from graphene (75%), because the hinge-like structure of phosphorene broken the selection rule in three phonon processes. Recently, (Zhao et al., 2018b) thought that the large anisotropic thermal transport stems from phonon group velocity rather than relaxation time. To be specific, the speed of sound velocity v 2 s is positively correlated with the Young's modulus E, and the group velocity of low-frequency phonons could be treated as the speed of sound. Meanwhile, the three-point bending method was employed to measure the Young's modulus of black phosphorene along crystallographic directions. It was found that the thermal conductivity anisotropy ratio is closely related to the Young's modulus one. Moreover, further enhancement of the anisotropy ratio can be realized by introducing phononic crystal (Xu and Zhang, 2016a) or loading tensile strain (Zhu et al., 2018). The latter regulation has also been experimentally proved, where Raman sensitive peak A 1 g , B 2g and A 2 g modes have opposite responses when the tensile strain was respectively applied to ZZ-and AC-directions.
Considering the limitations of synthesis technique, only few experiments were performed to extract the thermal conductivity of different thicknesses phosphorene films. Based on the measurement of micro-Raman spectroscopy,  demonstrated that the thermal conductivity of a 15 nm thick sample is ∼40 Wm −1 K −1 (20 Wm −1 K −1 ) along AC (ZZ) direction and decreases with the decrease of film thick. That is, the thermal conductivity anisotropic level has negative relation with film thickness, as shown in Figure 1. Another example, (Jang et al., 2015) adopted the TDTR method to investigate the thermal transport properties of silicon-supported thin films, and claimed that the thermal conductivity anisotropic ratio is about 2.5. Using the thermal bridge method, (Lee et al., 2015b) showed that when the thickness of phosphorous flakes reduces from 300 to 50 nm, the ZZ (AC) direction thermal conductivity drops to 12 from 27 Wm −1 K −1 (5 from 15 Wm −1 K −1 ). The first-principle calculation revealed the anisotropic thermal transport stems from orientation-dependent phonon dispersion and phonon-phonon scattering based on firstprinciple calculation. By comparison, the out-of-plane thermal transport is independent of the thickness of phosphorene film. Interestingly, (Su and Zhang, 2015) performed the optothermal Raman method to evaluate the substrate effect on the thermal conductivity phosphorene film, and the measured results indicated that the introduction of SiO 2 /Si substrate can enhance the thermal conductivity by more than one fold. The novel phenomenon seems to contradict the generally held notion held notion that the substrates have negative effect on the thermal conductivity of 2D materials.

Borophenes
Borophene is quite unique because there are rich hexagonal holes distributions in the triangular lattice such as the structure of borophene is mixed by triangle lattice and hexagonal lattice. At present borophene sheets with different phases and borophene nanoribbons have been experimentally synthesized on various substrates using the molecular beam epitaxy (MBE) approaches. For example, (Mannix et al., 2015) fabricated borophene sheets on the Ag (111) surface and found the planar islands with single atoms thickness occur at substrate temperature 570 and 650 K. Combined with theoretical calculations, the two phases of borophene can be identified as β 12 and χ 3 sheets. Moreover, (Li et al., 2018d) chose Al (111) substrate to grow boron because more charge can be transferred from Al (111) to boron atoms, which help enhance the stability of honeycomb borophene. The experimental results are very exciting, in which the perfect honeycomb-shape islands with single atom thickness form at substrate temperature 500 K and the lattice constant (0.29 nm) agrees well with the calculated value of freestanding honeycomb borophene. In others words, the type and temperature of substrate and is important for the formation of borophene morphology like 3D clusters or planar borophene. The successful preparation of several borophene structures provides opportunities for further exploring their electronic properties and device applications.
It is well known that light atomic mass and strong interatomic interaction are the crucial factors for the high thermal conductance of a nanostructure. Compared to carbon atom, the atomic mass of boron is light, and the Young's module of borophene is close to that of graphene , indicating that borophene probably possesses superior thermal conductivity as high as that of graphene. Using AGF method, (Zhou et al., 2017a) found the phonon thermal conductance of δ 6 borophene along the AC direction is as high as 7.8 nWK −1 at 300 K n the ballistic transport region, which is almost twice that of graphene. Moreover, using the ReaxFF potential,  predicted the thermal conductivity of δ 6 borophene to be about 147 and 75.9 Wm −1 K −1 along AC and ZZ directions, while it is about 233.3 and 102.5 Wm −1 K −1 using the SW potential. Here, the low thermal conductivity for borophene was explained by the following two reasons. On the one hand, small extension of acoustic dispersions and out-of-plane symmetry breaking leaded to stronger anharmonic scattering . On the other hand, the N-processes have very little contribution to thermal conductivity. All of that causes short phonon lifetime, and lowers thermal conductivity.  (Liu et al., 2013b;Gao et al., 2013;Li et al., 2015). Nevertheless, because of lattice matching constraints, only certain combinations of materials can enable continued heteroepitaxial growth. On the other hand, the presence of interfaces in heterostructures also induces some novel properties that are distinct from the parent materials. In particular, the lattice mismatch-induced strain would change the force environment of atoms, which cause the deformation of 2D crystalline structures and the variation of phonon spectra. At present, the ITC measurement of in-plane heterostructures is quite challenging compared to that of vdW interfaces. Hence, theoretical simulations provide a convenient way to understand interfacial thermal transport of in-plane heterostructures.
Benefitting from the mall lattice mismatch between grpahene and h-BN, their in-plane heterostructures as a typical case have been successfully synthesized (Gao et al., 2013). The simulated results of EMD method showed that the thermal conductivity along the vertical-interface direction is limited by the less conductive component (h-BN) and is insensitive to the details of the interface (Song and Medhekar, 2013), while that along the parallel-interface direction stabilises a constant which approaches the average of the parent materials, and gradually declines with increasing BN composition. The vibrational eigen-mode analysis reveals that increasing the BN composition would cause intermediate-frequency phonons to be localised in graphene domain . (Hong et al., 2016c) studied the effect of sample length and ambient temperature on the ITC in graphen/h-BN in-plane heterostructures via NEMD method, and showed that the ITC varies from 1.9 to 4.5 GW/m 2 K as the sample length in the range from 20∼100 nm and rising temperature could obviously enhance the ITC. Based on the same method and analogous models, (Liu et al., 2019a) found that increasing the number of atomic row can elevate (lower) the ITC in the zigzag (armchair) interface topography. To understand the mechanism of phonon transmission/scattering at interface, (Feng et al., 2017) calculated the temperatures of phonons under the framework of NEMD simulations. Taking graphene/h-BN planar interfaces as example, the temperature jump of the ZA mode is only 10∼20% of the overall temperature jump because there is considerable overlap between the acoustic branches of graphene and h-BN, while that of the ZO mode is as high as 220% of the overall one, stemming from the negligible coupling at the interface, as shown in Figure 2. Moreover, when interface mixing is introduced, the contribution ratio of ZA mode to the overall temperature jump increases from 19 to 41% and that of the in-plane modes decrease, indicating that a certain number of heat transfer channels shifts from the ZA mode to the TA/LA modes. Using AGF method, (Ong et al., 2016) demonstrated that alternating the bond orientation and type at the interface significantly affect the scattering and transmission of higher frequency modes. In the case of the zigzag-bond interface, the bonding type is able to dramatically alter the preferred orientation of the transverse optical (TO) phonons. During the crystal preparation process, a series of defects like single-vacancy (SV), stone−wales (SW) ineluctably exist in heterostructures. (Sutter et al., 2011) experimentally prepared graphene/h-BN in-plane heterostructures via a sequential CVD growth method, and pointed out that the interface of prepared heterostrucutres will be intermixing when the growth temperature becomes high. Based on NEMD simulations,  showed that the topological defect near the graphene/ h-BN interface could promote the interfacial thermal transport, which may be attributable to the localization of the stress fields. Analogously,  systematically studied the influences of SV and SW interfacial defects on the ITC in the in-plane heterostructures. The results indicated that the ITC quickly decreases with SV concentration but slowly with SW defects, and ultimately remain unchanged for large concentration. The phonon spectra and stress distribution analysis revealed that the regression of in-plane modes induced by SV defects results in the reduced thermal conductivity. Moreover, (Gao and Xu, 2017) reported that the existence of SW interfacial defects can realize the regulation of phonon transport direction, being equivalent to applying mechanical tensile strain. More recently, (Ni et al., 2019) reported the dependence of thermal conductivity on interface composition diffusion in graphene/h-BN in-plane heterostructures, and demonstrated that the system's thermal conductivity is closely related to the interface composition diffusion length because of the Anderson localization of phonons. (Liang et al., 2020) found that the ITC decreases with the increase of layer number n and saturates at n 3 in multilayer in-plane graphen/h-BN heterostructures. Even so, the asymptotic value still surpasses that in monolayer heterostructures.
Graphene/h-BN superlattices as special hetrostructures were forecasted to own numerous excellent properties like good thermoelectric performance (Yokomizo and Nakamura, 2013) and opto-electronics properties (Liu et al., 2013b). Using AGF method, (Jiang et al., 2011b) reported that there is a minimum in the period length l P dependence of the thermal conductance of graphene/h-BN superlattice. In general, the thermal transport characteristic of a superlattice indicates that thermal conductivity decreases first, and then increases with the increase of period length. In principle, to the left of the minimum value, heat conduction mostly relies on the coherent phonon transport. In turn, it is dominated by coherent phonons induced by interfacial modulation, which has been verified experimentally in 3D superlattices (Ravichandran et al., 2014) and theoretically in other 2D superlattices (Liu et al., 2014b). It is worth noting that the thermal conductivity of superlattices often decreases remarkably compared to that of parent materials, suggesting promising potential thermoelectric applications. Furthermore, (Zhu and Ertekin, 2014b) proposed that the initial reduction of thermal conductivity in graphene/h-BN superlattices is principally accounted for the effect of harmonic wave like phonon band gap and flat phonon branches. Subsequently, weakening interface scattering is conducive to thermal transport. Moreover, they also found that the critical period length is insensitive to the sample length. In addition, (Da Silva et al., 2016) carried out normal mode decomposition (NMD) approach to compute the thermal conductivity of graphene/h-BN superlattices with short l P , and showed that the phonons with frequency range 0∼25 THz accounts for over 90% of the total thermal conductivity regardless of the period length. Moreover, the reducing thermal conductivity with increasing l P and was attributed to the decrease of phonon group velocities rather than relaxation times. Recently, (Chen et al., 2016c) found that the critical period length decreases at higher temperatures in the superlattices, suggesting the suppression of phonon wave effect. Interestingly, for certain period lengths, the system's thermal conductivity at 200 K is very close to that at 300 K, contrary to previous studies that high temperature normally leads to the decrease of thermal conductivity, as shown in Figure 3. The reason for this unusual phenomenon was identified as the emergence of strong phonon wave interference under lower temperature. Thereafter, (Wang et al., 2017c) observed that the temperature profile behaves wave-like fluctuation in pristine/nitrogenated holey graphene superlattice at 100 K, which is also believed to be the result of the phonon wave interference. (Mu et al., 2015) reported that the existence of smooth interface will disrupts the coherence of phonons in supperlattice structure, and then lowers its thermal conductivity. Furthermore, the phonon transport properties of other 2D superlattices including silicene/ germanene superlattice  and MoS 2 /MoSe 2 superlattice  have also been studied using different simulation methods. These previous studies provide valuable insights on the phonon transport in in-plane heterostructures.
Except for graphene/h-BN in-plane heterostructures, some previous research also focused on the interfacial thermal transport across other 2D heterojunctions. Combined the firstprinciples with NEMD methods,  indicated that the ITC of graphene/MoS 2 in-plane heterostructure is higher than covalently bonded graphene/metal interfaces. Here, every Mo-C bond could act as a heat transfer channel. (Liu et al., 2014a) explored the dependence of ITC on temperature and strain in graphene/silicene in-plane heterostructure, and concluded that high temperature and small strain is beneficial to improve the ITC. In this work, it is worth noting that when the heat flux imposed into the system is lower than 42 GW m −2 , the value ITC almost varies with increasing the heat flux, corresponding to the Frontiers in Materials | www.frontiersin.org December 2020 | Volume 7 | Article 578791 Fourier heat conduction; while it exceeds the critical value, a lowfrequency kinetic wave is stimulated, constructing an extra channel for the non-Fourier heat transfer. (Qin et al., 2019) found that the ITC of MoS 2 /WSe 2 in-plane heterostructure is about 10% of that of graphene/h-BN heterojunction because of strong phonon scatterings at the interface, which could be useful for thermoelectric applications.  claimed that the ITC of phosphorene-based coplanar heterostructures doesn't present obvious anisotropy as the thermal conductivity of phosphorene because the interfacial anharmonicity results in a "mixing" effect for phonon transport. Moreover, the dependence of ITC on strain is very sensitive to the interface morphology, which can be understood by the analysis of stress distribution and phonon spectra overlap.

van der Waals Heterostructures
Benefitting from the progress of ingenious experimental methods, the synthesis of vdW heterostructures via self-assembly technology has achieved an unprecedented degree of control and accuracy, and provides new insights to explore structurefunction relationship (Frisenda et al., 2018) and to assemble complex devices (Abidi et al., 2018;Ding et al., 2018b) which are able to be tuned by stacking sequence, orientation, interlayer coupling, and external field. For example, graphene/hBN vdW heterostructures were prepared via an all-dry poly (methyl methacrylate) transfer process using polyvinylalcohol as the water-soluble sacrificial layer, which can avoid the polymer residual problem and prevent contamination from chemical solutions (Tien et al., 2016). Another example is that monolayer WS 2 and MoS 2 can be directly deposited on h-BN by CVD, which maintain their intrinsic bandgap of 2.01 and 1.85 eV (Fu et al., 2016), indicating a high quality of heterostructures with uniform distribution and a clean interface being developed. Thermal properties of vdW heterostructures is of great interest because it is highly relevant to the functionality, performance and reliability of such heterostructures-enabled devices. For bilayer graphene/h-BN heterostructures, there are periodic moire patterns on graphene layer. Three typical stacking configurations including AA, AB, AB′ stacking have been verified in experimental results. Moreover, different moire patterns may also affect the phonon transport in such heterostructures. In terms of the analysis of phonon dispersion, (Jung et al., 2012) found that the AB stacking is superior to the AA stacking for graphene/h-BN vdW heterostructures.  carried out the NEMD simulations to calculate the thermal conductivity of h-BN-supported graphene, and demonstrated that the system's thermal conductivity is reduced only 23% compared to the suspended case. That is because the force environment of graphene is barely affected by h-BN substrate. Additionally, they also discovered that the thermal transport properties strongly depend on the stacking ways of vdW heterostructures. Based on AGF method, (Yan et al., 2016a) showed that when the interlayer spacing between grapheme and h-BN layers is small, the perpendicular heat conduction is determined by in-plane acoustic modes. Another important finding in this work is that the ITC could be enhanced by more than 50% from C-N matched to C-B matched interface. Using the same method, (Wang et al., 2013b) found that the thermal conductance of graphene coupled with h-BN greatly reduces at low temperatures but can be still preserved as high as 96% that of graphene at room temperature. This could be explained by the fact that the h-BN coupling affects the out-of-plane vibrations, linearizing the flexural phonon modes and suppressing the transmittance. However, it has a minor effect on the in-plane vibrations, leaving the TA and LA phonons unchanged. Furthermore, the effect of different layers of h-BN (N) on the thermal conductance of graphene has also been studied, and the results showed that the presence of multi-layers h-BN (N ≥ 2) hardly affect the heat conduction performance of graphene, indicating that the design of the graphene/h-BN vdW heterostuctures is robust from the standpoint of heat dissipation. Analogously, (Pak and Hwang, 2016) predicted that the thermal conductivity of graphene changes very little after introducing h-BN substrate. The reason was that the strong interlayer coupling could minimize the internal stress fluctuations of graphene, as shown in Figure 4. Moreover, the ITC at gaphene/h-BN interface was computed to be 4.5 M Wm −2 K −1 . Acoording to lattice dynamics theory, (Zou and Cao, 2017) thought that the phonon group velocity of graphene supported by h-BN have minor changes in comparison with suspend graphene because of weak interlayer coupling. Meanwhile, the lifetimes of ZA phonons had significant reduction because of the invalidity of the selection rule. Using a fast transient technique, (Zhang et al., 2015a) observed that hydrogenating graphene can improve the ITC of grapheme/h-BN vdW heterostructures by more than 70%. That is because the hydrogenation could enhance the coupling degree between in-plane and out-of-plane phonons, and widens phonon channels. Furthermore, the reported ITC monotonically decreased with temperature and interatomic bond strengths.  used the Raman spectroscopy technique to measure the ITC of graphene/h-BN vdW heterostructures and the obtained ITC was 7.4 M WK −1 m −2 at room temperature, which exceeds the value of graphene/SiC interface. Apart from graphene/h-BN vdW heterostructure, the ITC of other vdW heterostructures have also been investigated. Based on MD simulations, (Pei et al., 2017) showed that constructing the graphene/phosphorene/graphene sandwiched structures can elevate the thermal stability and thermal conductivity of phosphorene.  found that the in-plane thermal transport in graphene/phosphorene vdW heterostructures exists obvious anisotropy and increasing the interlayer coupling could enhance the thermal conductivity of It is worth pointing out that for some specific l P , the temperature profile clearly behaves wave-like, suggesting that the number on the curves don't represent the real temperature originating from random atom motion. Reproduced from (Chen et al., 2016c) with permission from AIP Publishing LLC.
Frontiers in Materials | www.frontiersin.org December 2020 | Volume 7 | Article 578791 phosphorene layer because of phonons stiffening. With the aid of NEMD simulations,  found defect, hydrogenation and compressive strain can increase the ITC of graphene/phosphorene/grapheme multilayer structures. By BTE method, the room-temperature thermal conductivity of MoS 2 / MoSe 2 bilayer heterostructure was predicted to be 25.4 Wm −1 K −1 , which is intermediate between MoSe 2 and MoS 2 . Moreover, it was also found that increasing temperature will decrease the in-plane thermal conductivity, but increase the ITC. (Gao et al., 2016) showed that the interlayer coupling can reduce the in-plane thermal conductivity of graphene layer in graphene/ MoS2 bilayer heterostructure, but barely affects the thermal transport in MoS 2 layer. Similar phenomenon has been reported in the h-BN/silicene bilayer heterostructure by (Cai, et al., 2016). More importantly, they also found that the electronic and phonon transport are dominated by the silicone and h-BN layers, repectively, which realizes the decupling of phonon and electron transport.  summarized systematically the ITC of different grapheme-based heterostructurs, and the obtained results range from 2.8 to 25 M WK −1 m −2 at temperature range from 100 to 700 K. In a word, the ITC can be modulated by temperature, contact pressure, hydrogenation and vacancy defects. Experimentally, the ITC measurement of vdW heterostructures has been reported in previous studies.
Using the Raman optothermal technique, (Liu et al., 2017c) observed that the ITC of MoS 2 /h-BN vdW heterostructure (∼17 MW m −2 K −1 ) is a third of that of graphene/h-BN one (∼52 MW m −2 K −1 ). That is because the lighter mass of carbon atoms leads to larger frequency range that can be available for phonon transmission from the graphene to h-BN region. Another example, (Kim et al., 2018) used the space-resolved Raman spectroscopy technique to show that the doping level regulated by high electric field is beneficial to the spatial heat dissipation of graphene/h-BN devices, particularly close to the temperature of graphene breakdown.

MANIPULATION OF THERMAL TRANSPORT IN 2D MATERIALS AND THEIR HETEROSTRUCTRES
The manipulation of thermal transport in 2D heterostructres is crucial for the design of micro-nano devices based on them. However, since thermal transport is carried by a large population of phonons spanning a frequency spectrum and existing diverse interactions, controlling phonon transport is actually a complex task. Basically, the phonon scattering processes could be enhanced by introducing structural defects like isotopes, vacancy or dislocation, substrates, chemisorption, and strain. Note that using the mentioned-above controlling method, only the phonons with short wavelengths gets scattered efficiently by above tuning methods, while some long-wavelength phonons are nearly not influenced. Recently, a new thinking for manipulating phonon transport on the base of the wave nature of phonons like interference or local resonance has drawn considerable attention. Besides, the pseudospins and topological effects of phonons as new quantum degrees of freedom might be utilized to manipulate phonons (Gao et al., 2018;Liu et al., 2019b;Xiong et al., 2018).

From the Perspective of Phonon-Particle Picture
During the preparation of 2D materials, point-defect, isotopes and grain boundaries, inevitably resides in samples. The presence of imperfections normally introduces phonon scattering and thus causes decreasing thermal conductivity. The effect of randomly distributed isotopes was studied in h-BN via BTE calculations (Lindsay and Broido, 2011), demonstrating that the reduced thermal conductivity is the result of high-frequency phonon modes scattering induced by mass difference. Using AGF method,  found that the presence of isotope clustering could strengthen the cross section scattering. Analogously, (Morooka et al., 2008) demonstrated two interesting phonon transport behavior in GNRs with SWdefects: the first is that there is edge heat flux and the second is the appearance of circulating heat flux near SW defect. Subsequently, (Jiang et al., 2011a) showed that the probability of phonon transmission across the Si-doped region in graphene is very low, causing a significant reduction of thermal conductivity. (Huang et al., 2013) reported the influence of extended linedefects on the thermal conductance of GNRs, and showed that the modulation of the relative orientation between line-defect and transport direction can substantially change the system's thermal conductance. With aid of the NEMD simulations, (Ding et al., 2015) indicated the thermal conductivity of MoS 2 is reduced by over 60% for a 0.5% mono-Mo vacancy concentration. The analysis of SDE revealed that the reduction of thermal conductivity in defective MoS 2 is mainly determined by the phonon relaxation time. Analogously, (Zhou and Chen, 2015) found that the presence of nanoholes in graphyne nanoribbons will result in the intensive localization of phonons. In particular, a thought of "nanodomains in 2D alloy" has been proposed (Gu and Yang, 2016), because alloying can strongly scatter the highfrequency phonons, as shown in Figure 5. Moreover,  found that an unexpected enhancement of the ITC is observed when vacancy defects are introduced at the interface layers. The defect-induced enhancement in thermal transport is attributed to the increased friction and thus enhanced phonon coupling at the interfaces. Using the suspended thermal bridge method, (Aiyiti et al., 2018) experimentally verified that the introduction of defect is an effective approach to lower the thermal conductivity of few layer MoS 2 . Mechanical strain as a regular method has been widely applied to modulate the physical properties of 2D materials. In general, the thermal conductivities of 2D planar nanostructures like graphene or graphyne monotonically decrease with increasing the tensile strain because of phonon softening, which lowers the phonon group velocity (Wei et al., 2011). By contrast, in the buckled 2D materials like silicene or phosphorene, the dependence of thermal conductivity on the tensile strain shows a transition from the initial increase to subsequent reduction variation. For instance, based on MD simulations,  demonstrated that the loading of tensile strain can increase the relaxation times of ZA phonon modes in graphene. Using BTE method under RTA, (Kuang et al., 2016) reported that the applied tensile strain can either reduce or increase the thermal conductivity of graphene, which depends on the sample length. Apart from graphene,  demonstrated that the thermal conductivity of MoS 2 monotonically decreases with the increase of tensile strain.  performed the BTE simulations to study the controllable thermal conductivity in monolayer silicene by mechanical strain. Here, there was an up-and-down relation between the obtained thermal conductivity and the tensile strain, because a small tensile strain is able to partially recover the selection rule that only exists in 2D planar nanostructure. While the tensile strain becomes large, the appearance of phonon softening lowers the phonon group velocity. Moreover, using NEMD method, (Zhu and Ertekin, 2015) reported that the thermal conductivity of graphene/h-BN superlattice initially increases to a peak value, after which it then decreases with further strain. The non-monotonic curve was considered as a competition between in-plane softening and flexural stiffening of phonons.  showed that as a compressive strain of 0.18 is applied in the graphene/black phosphorus vdW heterostructures, the ITC increase more than 15-fold compared to the strain-free value. The significant enhancement of ITC induced by compressive strain could be understood from the following two aspects. The first is that the applied normal pressure on the multilayer structures clearly strengths the interlayer coupling. The second is that the compressive strain also enhances the shear interaction between the interlayers and thus improves the thermal transport contribution from in-plane phonons. Experimentally, (Guo et al., 2020) reported the tensile strain dependence of the thermal conductivity of monolayer graphene using optothermal Raman method. The measured results demonstrated that under a 0.12% biaxial tensile strain, the system's thermal conductivity drops approximately by 20% at 350 K, which verifies the previous theoretical literatures prediction (Varshney et al., 2018). Similar phenomenon has also been reported in polycrystalline multilayer graphene based on the microbridge measurement (Zeng et al., 2020c). Additionly, using TDTR method, (Meng et al., 2019) showed that a 9% crossplane compressive strain can cause the 7-fold increase of the cross-plane thermal conductivity of multilayer MoS 2 , which is available for improving MoS 2 -based device performance and stability, as shown in Figure 6. The analysis of lattice dynamics and phonon spectroscopy revealed that the drastic increase of cross-plane thermal conductivity with pressure is mainly the result of strengthened interlayer forces and enhanced group velocity of longitudinal acoustic phonons, while the saturation of thermal conductivity above 15 GPa is associated with the combined effects from increasing group velocity and reduced phonon lifetimes. It is conceivable that the observed phenomena should occur in most vdW heterosturctues. The third common approach to modulate the thermal conductivity of 2D heterostructres is chemical functionalization. Using SED technology, (Pei et al., 2011;Kim et al., 2012) found that the hydrogenation coverage dependence of thermal conductivity in functionalized graphene shows a U-shaped variation, similar to the composition dependence of thermal conductivity in Si/Ge alloy (Xie et al., 2020). (Mu et al., 2014) reported that a 5% oxidation coverage leads to a 90% reduction of thermal conductivity in graphene samples. It can be explained from the following two aspects: one is that the existence of the oxygen atoms adds phonon-impurity scattering; the other is that the reflection symmetry of pristine graphene plane is destructed owing to oxygen coverage, which enhances the scattering of flexural modes. In some cases, hydrogenation is beneficial for the thermal transport in 2D materials. According to the BTE calculated results,  demonstrated that the hydrogenations can enhance the phonon transport in pentagraphene (up to 76% increase), which ascribes to the wakening anharmonicity of C-C bond. Using EMD method,  found that when the coverage of the adsorbed groups exceeds 0.1, the thermal conductivity of functioned h-BN remains constant (about 100 Wm −1 K −1 ), which are still sufficiently thermally conductive as fillers of composites. (Hong et al., 2016a) showed that the ITR of graphene/phosphorene bilayer heterostructure increases monotonically with the hydrogen coverage in graphene layer. The maximum addition of 84.5% is observed at 12% hydrogen coverage. The enhanced thermal transport can be understood that an extra H-P heat dissipation channel is created in addition to that between C-P atoms and the FIGURE 6 | A typical atomic configuration for the Mo 1−x W x S 2 alloy embedded with triangular WS 2 nanodomains (B) The accumulative thermal conductivity of Mo 0.5 W 0.5 S 2 with and without nanodomains embedded as a function phonon frequency. (C) The phonon nanodomain scattering rates of LA phonons from the Γpoint to the M point and the phonon nanodomain scattering rates of long-wavelength LA phonons and short-wavelength LA phonons. Reproduced from (Gu and Yang, 2016) with permission from American Physical Society.
Frontiers in Materials | www.frontiersin.org December 2020 | Volume 7 | Article 578791 absorbed H atoms on graphene behave as scattering centers, thereby enhancing the graphene's lateral to flexural direction phonon coupling which indirectly strengthens the thermal transmission from graphene to phosphorene.  found that the in-plane thermal conductivity of phosphorene layer in graphene/phosphorene vdW heterostrucutres is significantly improved by hydrogenating the graphene layer.
When 2D materials or related heterostructres are integrated for device applications, their thermal transport properties have been drastically altered on account of the phonon-substrate scattering. The heat generated in device geometry may flow in the sample, and is dissipated directly into the substrate through their interfaces. (Sadeghi et al., 2013) experimentally observed that when the number of layers in SiO 2 -supported few-layer graphene achieves 34, the system's thermal conductivity almost equated with the natural graphite, because there are longwavelength out-of-plane phonons and partially diffusive interface-phonon scattering in thick sample. In general, the substrate will cause the enhancement of optical phonon scattering in 2D materials (Berciaud et al., 2009). In terms of SED analysis, (Qiu et al., 2012b) argued that the substrate effect is responsible for the reduction of the phonon relaxation time of ZA modes. Interestingly, (Su and Zhang, 2015) experimentally demonstrated that the SiO 2 substrate could double the thermal conductivity of phosphorene film, which differs from the case in SiO 2 -supported CNTs (Ong et al., 2011). Based on EMD simulations,  revealed that the degree of thermal transport anisotropy can be significantly enhanced in supported phosphorene. That was because the substrate effect significantly suppresses the ZA modes that dominate the thermal transport along armchair direction.  found that the ITC of graphene/h-BN in-plane heterostructrue is enhanced as it couples with a SiO 2 substrate, which ascribes to the increase of thermal transport channels. Similar simulated results have been reported in asymmetric graphene/h-BN inplane heterosturctures (Sandonas et al., 2017).

From the Perspective of Phonon-Wave Picture
Considering that phonons are waves of the atomic lattice, it is natural to think phonons interference, which could be used to regulate thermal transport (Chen et al., 2005;Xie et al., 2018a;Xie et al., 2018b;Zeng et al., 2020a). Actually, in periodic superlattice, phonons would be reflected by the periodic interfere as that the photons in photonic crystals. To produce the interference patterns, it should satisfy the following two conditions. First, the interface roughness is much less than the phonon wavelength, in which the phase of phonons could be preserved after undergoing interface scattering. However, the wavelengths of most phonons at room temperature are rather short in superlattice. Hence, atomically smooth interface are required for distinct phonon interference. Second, the spacing of two adjacent interfaces (namely period length) should also be relatively short so that propagating phonons will not be scattered before being reflected at the interfaces. For this, the coherent phonons consisting of the constructive interference of incident and reflected phonons could preserve their phases across multiple interfaces. That is, when the period length is sufficiently long, phonons are difficult to produce interference patterns, and so the phonon transport is diffusive and corresponding thermal conductivity monotonically increases with the number of scattering interfaces. From the kinetic theory perspective, the phonon interference might generate a phonon-gap, and decline the system's thermal conductivity to some extent. More importantly, the interference could flatten the phonon dispersion in a large frequency range, which affects the group velocity, as discussed in In-Plane Heterostructures. Note that the strict condition of phonon interference actually limits the applications of one-dimensional superlattice.
Lately, a metamaterial constructed by arrays of pillars on the nanofilms/nanowires has attracted more and more attention. The presence of nanopillars causes the formation of various resonant modes in the whole frequency range. If the resonant modes have the same polarization and frequency with the propagating modes from the base material, both kinds of phonons will hybridize due to the band anticrossing effect, and the phonon group velocity is lowered. The advantage of this design is that the long-wavelength phonons can be effectively controlled. For instance, using BTE method, (Davis and Hussein, 2014) calculated the thermal conductivity of nanopillared silicon thin films, and found that the obtained value is less than half of that in pristine nanowires. This is attributed to the hybridizations taking place at both subwavelength and superwavelength frequencies. Moreover, based on the NEMD simulations, (Xiong et al., 2016) reported that introducing resonant structures to alloying structures could achieve ultra-low thermal conductivity (0.9 Wm −1 K −1 ). It can be understood that the local resonance effect impedes the transport of low-frequency phonons and alloying can effectively scatter high-frequency phonons. Subsequent studies also confirmed this conclusion Nomura et al., 2018). Apart from silicon-based nanostrucutres, (Ma et al., 2018) reported a reduction of half thermal conductivity in nanopillared GNRs compared to pristine nanoribbons. It was also found that isotopic doping in nanopillars is conducive to phonon transport, because doping could weaken the hybridization degree. By designing the pristine/branched GNR junctions,  proposed a high-performance thermal rectifier where its TR ratio is as high as 400% under a small temperature bias. The underlying mechanism was found to be that when the longitudinal phonons exists obvious local resonance for the negative heat flux. Experimentally, (Anufriev et al., 2017) observed a reduced thermal conductivity in silicon films by attaching aluminum pillars. However, they argued that the reduction of thermal conductivity originates from interface roughness induced by the metal deposition rather than phonon local resonance. In other words, the validity of phonon local resonance lowering thermal conductivity is still not clear. Furthermore, the design 2D hole-based phononic crystal is another feasible method to manipulate phonon transport. (Feng and Ruan, 2016) claimed that the thermal conductivity of graphene phononic crystal (GPC) is three orders of magnitude lower than the pristine one and decreases exponentially with increasing porosity. The Frontiers in Materials | www.frontiersin.org December 2020 | Volume 7 | Article 578791 vibrational eigen-mode analysis revealed that the phonon localization and phonon back scattering around the nanopores are the major factor responsible for the ultra-low thermal conductivity of GPC. Similar results have also been reported in some other studies . Analogously, (Cui et al., 2018a) reported that the thermal conductivity of silicene phononic crystal is not sensitive to temperature but strongly depends on the porosity ratio, as that in 3D phononic crystal (Hopkins et al., 2011). There were three reasons for that, namely increasing phonon scatterings, inducing phonon localization and phonon coherent effects.

PHONON-DRIVEN EMERGING APPLICATIONS OF 2D HETEROSTRUCTRES Thermal Management
From the point of view of informatics, phonons as heat carriers can be also used for data storage and processing, namely called phononics. Various thermal devices (corresponding to the electronic counterparts) like rectifier, transistor, and logic gate have been conceptualized Jia et al., 2020), making the design of thermal logic circuits possible. Then, some of them are gradually realized via designing different sophisticated nanojunctions. To be specific, thermal rectifier as basic components that enables heat flow preferentially in one direction has been explored both theoretically and experimentally. For example, (Chang et al., 2006) experimentally observed the TR phenomenon in gradual massloading BN nanotubes opened the active regulation of heat flux in low-dimensional materials (Chang et al., 2006). After that, numerous architecture models like graded core-shell nanowires (Liu et al., 2014e), asymmetric nanobibbons Sandonas et al., 2017; and phase change materials (Kang et al., 2018) have been proposed to use as nanoscale thermal rectifiers. For example, the TR ratio of the single-layer Y-shaped GNRs was predicted to be about 45% under a normalized temperature bias (Δ) of 0.5 and increases with the number of layer . With the aid of AGF method, (Ouyang et al., 2010) investigated the ballistic phonon transport in three-terminal GNRs and showed that the TR effect rises quickly with increasing the structural asymmetry. Moreover, the corner shape of nanojunctions was found to play an important role in the rectification effect, and the TR ratio could approach 200% via the parameters optimization.  confirmed that the TR effect in asymmetric GNRs gradually weakens with increasing the sample width. The mechanism origin was that phonon lateral confinement results in the mismatch of phonon spectra, phonon edge localization and the temperature and space dependence of thermal conductivity, all of which colludes to produce asymmetric phonon transport. (Tian et al., 2012a) experimentally designed a novel solid-state thermal rectifier based on the triangular-shaped graphene oxide. The measured TR ratio was up to 121% under a temperature bias of 12 K and increased with the increasing of the angle.
Apart from geometric asymmetry and mass asymmetry, the introduction of functionalized or defected engineering is also able to realize efficient TR effect. (Pal and Puri, 2014) proposed asymmetrically polymer-functionalized carbon nanotubes (CNTs) to obtain high TR ratio. When the applied temperature bias was 40 K, the TR ratio could be as high as 204% in this complex structure. From this work, it can be summarized that the structural asymmetry and the strong temperature dependence of vibration density of states are two necessary conductions to obtain high TR ratio. (Melis et al., 2015) reported that the TR ratio is about 54% in pristine/hydrogenated graphene heterojuction with triangular morphology. (Hu et al., 2017) proposed the thought of improving TR performance by tendering multi thermal rectifiers. (Nobakht et al., 2018) studied TR in defect-engineered graphene with asymmetric holing arrangements. It was found that the TR ratio increases as porosity increases, and can reach as high as 78% at room temperature with triangular hole arrangement with porosity of 75%, enabling effective and precise TR. Using microbridge method,  experimentally demonstrated the measured TR ratio is about 26% in asymmetric defect graphene, which is ascribed to temperature and space dependence of thermal conductivity. It is worth pointing out that the aforementioned good TR effect commonly requires large temperature bias and delicate nanostructures namely complicated synthetic procedures, which limits the practical applications. Recently, phase change materials show promising potential in high-performance thermal rectifier. Utilizing the crystalline polymer nanofiber/cross-linked polymer nanofiber heterojunctions, (Zhang and Luo, 2015b) constructed a new thermal rectifier and its workflow was that the polymer nanofiber portion undergoes a phase transition and its thermal conductivity dramatically alters under temperature gradient. As a result, a 20 K temperature bias stimulated ∼100% TR ratio. Detailed spectral analysis revealed that the large TR ratio was ascribed to by the vibrational mismatch of sulfur and carbon atoms. (Kang et al., 2018) developed an analytical framework upon a material architecture for actualizing one type of nonlinear thermal transport, and applied the analytical framework to a junction comprising two phase-change materials (polyethylene and vanadium dioxide). The results showed that such heterojunction has a maximal theoretical rectification factor of approximately 140%. However, associated problems include slow phase change and narrow temperature range, thus hindering practical applications.
Actually, constructing nanojunctions is an effective method to implement efficient thermal rectifier. Base on NEMD simulations,  demonstrated that the TR ratio is as high as 1,200% in graphene/CNTs junction under 300 K temperature bias, approaching the actual application, which is considered as the formation of standing wave along negative heat flux direction. (Chen et al., 2019b) reported a high TR ratio of 334% in carbon/h-BN heteronanotubes can achieve 334%. Similar phenomenon has also been found in gaphene/h-BN in-plane heterostructures (Chen et al., 2016d). Note that a moderate C/BN ratio and parallel arrangement help improve the TR ratio in CBN hybrids . (Dong et al., Frontiers in Materials | www.frontiersin.org December 2020 | Volume 7 | Article 578791 2019) reported that the TR ratio in Au/CNTs heterojunctions decreases with the number of the molecular bridges but increases with the length of CNTs, reaching 375% under 60 K temperature bias. (Sandonas et al., 2017) indicated that the TR ratio in supported graphene/h-BN in-plane heterojunctions increases from 8.8 to 79% by lowering the substrate temperature.
Recently, (Wei et al., 2019) proposed a new thermal rectifier base on multilayer graphene-based with interlayer gradient functionalization. The TR mechanism was that the vocational of each graphene layer exists obvious difference.  reported a distinct TR effect in asymmetric graphene/ h-BN vdW heterostructure. Moreover, the TR ratio increased rapidly with the interlayer coupling level, as shown in Figure 7A. Detailed spectral analyses revealed that for strong interlayer coupling, the low-frequency phonons presents the characteristic of localized modes along negative heat flux direction, while Δ > 0, medium-frequency phonons become high-efficient heat conduction channels, as illustrated in Figures 7B,C. Apart from TR effect, another crucial thermal transport phenomenon is negative differential thermal resistance (NDTR) in which heat flux is inversely correlated with thermal bias. Actually, NDTR is the foundation not only for thermal switching, but also for thermal memory. By now, the phenomenon has been widely investigated in non-linear lattice models (Chan et al., 2014). Some obtained results showed that the existence of interface effect is a requirement for generating NDTR. Subsequently, (Ai et al., 2012) found that the phenomenon of NDTR only exists in single-layer GNRs with appropriate width. For multi-layer GNRs, NDTR regime becomes gradually smaller and even disappears with increasing the number of layers. By designing sandwiched nanofluids,  demonstrated that the heat flux is suppressed and even prohibited when temperature bias becomes sufficiently large, because nanofluids are adsorbed on the solid surface under negative temperature bias, and the free fluid molecules is decreased. Furthermore,  studied the nonlinear thermal transport in graphene/h-BN in-plane heterostructure and showed that when Δ surpasses a threshold, NDTR phenomenon appears. It could be understood that out-of-plane kinetic wave is stimulated in graphene domain under large temperature bias, which results in the out-of-plane vibration mismatch between graphene and h-BN regions, impeding interfacial thermal transfer, as illustrated in Figure 8. It was also found that the regime of NDTR becomes small as the sample length increases, but barely changes with the sample width. That is, the wave-dominated thermal transport mechanism becomes dominant when the sample length is small or the applied temperature bias is large. To realize logic calculations uisng phonons, (Pal and FIGURE 7 | (A) Schematic of thermal conductivity measurement with a diamond anvil cell integrated with a picosecond TRTR, (B) Extracted cross-plane thermal conductivity (dominated by phonons) as a function of pressure. The red curve is included only as a guide to the eye, (C) Pressure-dependent coherent longitudinal acoustic phonons (LAP) frequencies extracted from coherent phonon spectroscopy (CPS) measurements (black circles) and LAP group velocities from first-principles calculations (red squares), (D) Pressure-dependent LAP lifetimes extracted from CPS measurements (black circles) and from first-principles calculations (red squares), as well as lifetimes of A 1g optical phonons extracted from Raman measurements (green triangles). Reproduced from (Meng et al., 2019)  Puri, 2015) designed a thermal AND gate consisting of asymmetric GRN and the test of a truth table demonstrated that with either inputs at "0" state, the output value is in the range 0∼0.27, which is comparable to silicon-based logic gate (0∼0.12). Note that the running speed of the thermal device reaches 100 ps owing to ultrahigh phonon group velocity in graphene. (Liu et al., 2015) experimentally demonstrated a thermal modulator consisting of a partially clamped singlelayer graphene, and found the high/low heat flux ratio can remarkably is as high as 150% by regulating the external pressure. The reason for the decrease of heat flux was that the phonon spectra mismatch of the unclamped and clamped graphene sections leads to strong phonon-interface scattering. Another experiment, (Sood et al., 2018) prepared a thermal transistor based on multi-layer MoS 2 driven by reversible electro-chemical intercalation of Li ions. The thermal device is in the following: during the charging the cell, Li ions continue to enter the multilayer MoS 2 and some 2H phase MoS 2 layers are converted into the 1T phase, namely 2H-1T phase interface, which can lowers phonon group velocity and enhances phonon scattering rates. Moreover, lithiation may induce disordered stacking among multilayer MoS 2 and hinders the interlayer thermal transport, causing the decrease of ITC from 15 to 1.6 MW m −2 K −1 . For the delithiation step, the cross-plane thermal conductance increases back to the pre-lithiation value. The obtained on/off ratio is about 10, opening the possibility of exactly controlling heat flux, as shown in Figure 9.

Thermoelectric Conversion
Thermoelectric materials enable heat energy to directly convert into electric power. The efficiency of thermoelectric conversion is determined by a dimensionless parameter, namely figure-ofmerit (S 2 σT/(κ e + κ p )), where S is Seebeck coefficient, σ is electrical conductance, T is absolute temperature, κ e and κ p is electron and phonon thermal conductance, respectively. It is not difficult to find that enhancing power factor (S 2 σ) and reducing thermal conductivity are both effective strategies to promote thermoelectric performance. However, according to the Wiedeman-Franz law, S, σ and κ e depend on each other. In contrast, inhibiting phonon transport is more viable realizing high conversion efficiency. In the last few decades, a huge amount of effort has been conducted to reduce the thermal conductivity of thermoelectric materials via nanostructure engineering . Reducing dimensionality as an important strategy for enhancing thermoelectric performance has been proposed by (Heremans et al., 2013), and its working mechanisms are as follows: 1) the size-quantization in low-dimensional The schematic diagram illustrates the mode-matching between region A and B. Clearly, the mode with frequency ω A in left side is basically able to match a specific mode with ω B in right side, but the reverse is not true. (C) The spatial distribution of a low-frequency vibrational modes at 4.3 THz and 6.8 THz. Reproduced from  with permission from American Chemical Society.
Frontiers in Materials | www.frontiersin.org December 2020 | Volume 7 | Article 578791 nanostrucutres facilitates increasing the Seebeck coefficient; 2) the phonon-boundary scatterings can reduces the system's thermal conductivity. In fact, it has reported that the electronic transmission exists obvious fluctuation within the first conductance band for asymmetric GNRs, and consequently the Seebeck coefficient is enhanced (Xie et al., 2012;Cao et al., 2019). (Yokomizo and Nakamura, 2013) demonstrated that the value of S in graphene/h-BN superlattices is twenty times larger than that in graphene. (Mazzamuto et al., 2011) reported that the electron transport exits resonant tunneling behavior in mixed armchair/zigzag GNRs, which corresponds to a ZT value of 1 at room temperature. On the other hand, (Tran et al., 2015) found the phonon conductance is significantly inhibited in GNRs with BNNRs branches, but the electron transmission is equivalent to that in pristine GNRs, similar to that in vdW graphene junction (Hung Nguyen et al., 2014). Moreover, (Cui et al., 2018b) reported that the ZT value of graphyne nanoribbons with edge disorder can achieve 2.5 because of the strong phonon-boundary scatterings. Along the way, the MoS2/WSe 2 (Jia et al., 2019) and graphene/ MoS 2 vdW heterostructures (Sadeghi et al., 2016) have been proposed to construct thermoelectronic devices with excellent performance. Lately, the electronic devices based on organic single-molecules have been shown excellent electronic properties such as quantum tunneling and magnetic change Deng et al., 2017;Xie et al., 2017;Zhang et al., 2019) and some of them can be also used as the competitive thermoelectric materials Wu et al., 2020;Zeng et al., 2020b). Actually, sometimes it should balance the power output and thermoelectric conversion efficiency. That is, blindly reducing the system's thermal conductivity is not always better (Narducci, 2011). As shown in Figure 10C,  reported an 80% reduction for the thermal conductivity of polycrystalline graphene/h-BN in-plane heterostructures. Meanwhile, the electrical conductivity also declined rapidly. As a result, the calculated ZT value was only 0.01. Experimentally,  measured the electron/phonon transport in the twisted bilayer graphene, and demonstrated that the cross-plane thermoelectric transport is driven by the scattering of electrons and breathing phonon modes that represent a unique "phonon drag" effect across atomic distances, as illustrated in Figure  10A. Moreover, it was also found that the value of S in gold/ graphene/WSe 2 vdW heterostructrues approaches 100 μV/K (Rosul et al., 2019). To sum up, 2D heterostructures have great potential for the applications of thermoelectric conversion. However, theoretical calculations still have large limitations due to the limitation of computation cost and the handing question of van der Waals force FIGURE 9 | (A) Heat flux J vs. temperature bias Δ for the graphene/h-BN heterojunction. Here, T L 500 K and TR ranges from 480 K to 100 K. (B) The temperature bias dependence of matching coefficients S between the in-plane/out-of-plane power spectra of graphene and h-BN regions. The dependence of (C) NDTR and (D) out-of-plane spectra vs. the heterojunction width. Reproduced from  with permission from AIP Publishing LLC.
Frontiers in Materials | www.frontiersin.org December 2020 | Volume 7 | Article 578791 in the DFT calculation. In particular, the electron-phonon coupling was not considered in the calculation process. In fact, some previous studies have indicated that the electron-phonon coupling can severely affect the electron transport properties such as spin filter and magnetoresistances Li and Chen, 2017). Again,  found that the presence of bipolar effect in SnP 3 could lower the system thermal conductivity and help enhance the ZT value (∼3.4). Hence, new theoretical method including multiparticle coupling is necessary to predict the thermoelectric performance of nanostructurs.

CONCLUSION AND OUTLOOK
In this review, we have introduced various experimental and theoretical methods for the thermal transport in low-dimensional materials. Then, the thermal transport properties of 2D materials and their in-plane/vdW heterostructures were discussed. Furthermore, we have summarized different thoughts to manipulate thermal conductivity. In practical terms, the introduction of isotopes, dislocation, substrates, chemisorption, and strain can strength diffuse phonon scattering (particle-like), which is valid with a limited range of high-frequency phonons.
Meanwhile, the wave nature of phonons like interference or local resonance can regulate low-frequency phonons. The diversity of thermal properties in 2D heterostructures indicates important potential applications in thermoelectric, thermal dissipation and thermal devices fields, which needs further investigations. Advanced fabrication techniques have made multistage circuit integration of 2D heterostructures possible, and resolving the thermal dissipation across the interfaces also looms ahead when designing nanodevices. Considering that heat transfer in semiconductor device is carried by a large population of phonons spanning a frequency spectrum and existing diverse interactions, controlling phonon transport is actually a complex task in these nanostructures and some important questions are still not clear. First, the phonon transport in nanostructures under coupling multi-field effect like electron-phonon or magneton-phonon couplings is in the primary stage of development. Some already obtained results are based on either theoretical models or based on multiple hypotheses, and haven't been verified by experiments. Second, because of the multi-scale character of phonon transport, the accurate control of phonon transport in specific frequency range remains surprisingly difficult. Moreover, the anomalous non- Fourier transport caused by the excitation of the kinetic wave needs further confirmation in experiment. Third, on the experiment side, the accurate measurement of thermal properties of 2D heterostructures still has a long way to go.

AUTHOR CONTRIBUTIONS
X-KC, Y-JZ and K-QC conceived the manuscript. All authors reviewed and discussed the manuscript.