hexagonal boron nitride heterostructures: Insights into the interfacial thermal transport properties

International Journal of Heat and Mass Transfer 151 (2020) 119395

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/hmt

Multilayer in-plane graphene/hexagonal boron nitride heterostructures: Insights into the interfacial thermal transport properties Ting Liang, Man Zhou, Ping Zhang∗, Peng Yuan, Daoguo Yang School of Mechanical and Electrical Engineering, Guilin University of Electronic Technology, No. 1 Jinji Road, Guilin, Guangxi 541004, China

a r t i c l e

i n f o

Article history: Received 27 September 2019 Revised 28 December 2019 Accepted 16 January 2020

Keywords: Multilayer-mixed heterostructures vdW interactions Interfacial thermal conductance Phonon transmission Coupling strength

a b s t r a c t Combining both vertical and in-plane two-dimensional (2D) heterostructures opens up the possibility to create an unprecedented architecture using 2D atomic layer building blocks. The thermal transport properties of such multilayer-mixed heterostructures, critical to various applications in nanoelectronics, however, have not been thoroughly explored. Herein, we construct two configurations of multilayer in-plane graphene/hexagonal boron nitride (Gr/h-BN) heterostructures (i.e., mixed heterostructures) via weak van der Waals (vdW) interactions and systematically investigate the dependence of their interfacial thermal conductance (ITC) on the number of layers using non-equilibrium molecular dynamics (NEMD) simulations. The computational results show that the ITC of two configurations of multilayer in-plane Gr/h-BN heterostructures (MIGHHs) decrease with increasing layer number n and both saturate at n = 3. Surprisingly, we find that the MIGHH is more advantageous to interfacial thermal transport than the monolayer in-plane Gr/h-BN heterostructure, which is in strong contrast to the commonly held notion that the multilayer structures of Gr and h-BN suppress the phonon transmission. The underlying physical mechanisms for these puzzling phenomena are probed through the stress concentration factor, overlap of phonon vibrational spectra and phonon participation ratio. In particular, by changing the stacking angle of MIGHH, a higher ITC can be obtained because of the thermal rectification behavior. Furthermore, we find that the ITC in MIGHH can be well-regulated by controlling the coupling strength between layers. Our findings here are of significance for understanding the interfacial thermal transport behaviors of MIGHH, and are expected to attract extensive interest in exploring its new physics and applications. © 2020 Elsevier Ltd. All rights reserved.

1. Introduction Van der Waals (vdW) heterostructures have attracted enormous interest in the past decade because of their novel device functionalities and superior electronic and thermal regulation properties [1–5]. Most of these captivating properties are strongly dependent on the different spatial configurations of heterostructures, typically including the vertical heterostructures [1,6,7] that vertically stacked multiple two-dimensional (2D) materials and the in-plane heterostructures [8–11] that seamlessly stitched together two different atomic monolayers through the covalent bonds. Because of the weak interactions of the vdW force in the vertical heterostructures, the heat flux through them is much lower than that through the in-plane heterostructures [4,12], making the

∗

Corresponding author. E-mail addresses: [email protected] (T. Liang), [email protected] (P. Zhang). https://doi.org/10.1016/j.ijheatmasstransfer.2020.119395 0017-9310/© 2020 Elsevier Ltd. All rights reserved.

in-plane heterostructures have a broader application prospect as thermal management materials. Still, for the vertical heterostructures, they can be considered as the basic materials for designing the next generation of high-performance wearable devices and the vdW heterostructure-enabled sensors, because of the diversity of their spatial configurations [13–15]. Because each configuration has outstanding advantages, it is conceivable that the combination of vertical and in-plane heterostructures will make it possible to create unprecedented architectures by using 2D atomic layer building blocks [16]. Graphene/hexagonal boron nitride (Gr/h-BN) heterostructure, which benefits from the small lattice mismatch between Gr and h-BN, has been successfully synthesized both vertical [17,18] and coplanar [2,19,20] configurations, paving the way for constructing the multilayer-mixed heterostructures. In parallel to the efforts on pursuing a more sophisticated synthesis method, the interfacial thermal transport properties of Gr/h-BN heterostructure, which play a pivotal role in the high-performance thermal interface materials and heterostructures devices [21–23], are urgently needed to

2

T. Liang, M. Zhou and P. Zhang et al. / International Journal of Heat and Mass Transfer 151 (2020) 119395

be understood. Accordingly, by using a transient heating technique, Zhang et al. [24] showed that the interfacial thermal resistance in the vertical Gr/h-BN heterostructure is affected by interatomic bond strength, heat flux direction, and functionalization. Relatively, Chen et al. [25,26] found that the in-plane Gr/h-BN heterostructure exhibits a negative differential thermal resistance behavior, which is caused by the phonon resonance effect and the lattice vibration mismatch. In addition, topological defects [4], doping and interface topography optimization [27] were employed to effectively improve the thermal energy transport across the in-plane Gr/h-BN heterostructure interface. These previous studies have enriched our understanding of the interfacial thermal transport properties of Gr/h-BN heterostructure. However, most of the Gr/h-BN heterostructure simulations are limited to a single vertical or inplane configuration without paying attention to their combinations, which may create many interesting thermal transport properties and contribute to the design of functional heterostructures with high-density multi-pixel capabilities. In the present study, the in-plane Gr/h-BN heterostructure is stacked vertically via weak vdW interactions, and two types of multilayer in-plane Gr/h-BN heterostructures (MIGHHs), i.e., mixed heterostructures, are constructed. It has been demonstrated that the thermal conductivity of multilayer Gr and h-BN decreases with the increasing number of layers [23,28–31], which prevents them from creating new applications in various emerging fields. Therefore, it is a meaningful topic to explore whether the change in the number of layers has the same effect on the thermal transport properties across the interface of multilayer-mixed heterostructures. To this end, by using a series of molecular dynamics (MD) simulations, the interfacial thermal conductance (ITC) of two types of MIGHHs is investigated. In addition, we further study the effect of coupling strength between different in-plane Gr/h-BN heterostructure layers on the ITC. To reveal the interfacial thermal transport mechanisms in MIGHH, the interfacial stress distributions, phonon density of states (PDOS) in the frequency domain and phonon participation ratio (PPR) reflecting the phonon characteristics are performed to analyze and compare in our research. 2. Models and methods To model the MIGHH, constructing the monolayer in-plane Gr/h-BN heterostructure is a top priority. Because previous theoretical [32] and experimental [32,33] studies have shown that the zigzag stitching edges are preferably formed in the in-plane Gr/hBN heterostructure, only the interface of the zigzag arrangement is considered here. Fig. 1a depicts two different structures of monolayer in-plane Gr/h-BN heterostructure, called C–NB and C–BN, respectively, depending on the bonding type of different atoms at the interface. It is known that although both C–N and C–B interactions are covalent in nature, the strength of the C–N bond is higher than that of C–B bond, leading to the advantage of the C– N interface for phonon transmission [4,34]. Therefore, when constructing the MIGHH model, we use the in-plane Gr/h-BN heterostructure of C–NB, while C–BN is only used for reference and comparison. In this work, two types of MIGHHs are constructed. One is to directly stack the monolayer in-plane Gr/h-BN heterostructure in order, named as multilayer parallel stacked heterostructure (MPSH), as shown in Fig. 1b and c. The other is to rotate the even-numbered in-plane Gr/h-BN heterostructure by 180 degrees on the basis of the MPSH configuration to form the multilayer staggered stacked heterostructure (MSSH), as shown in Fig. 1d and e. To construct the initial configuration of the MIGHH, the interlayer distance d (c.f. Fig. 1b) is temporarily set at 0.35 nm, and it is automatically optimized to a reasonable value after

the system reaches a steady state, depending on the interactions strength of the vdW force between the layers. Meanwhile, the y-axis is set to the direction in which the heat flux is applied, while the x-axis is parallel to the heterostructure interface, as schematically shown in Fig. 1d. The dimensions of all the simulated systems are w = 62.12 and l = 215.50 A˚ in the x-axis (along the zigzag direction) and y-axis (along the armchair direction), respectively. All the MD simulations are performed using the large-scale atomic/molecular massively parallel simulator (LAMMPS) package [35]. Fixed boundary conditions are applied in the transport direction, whereas periodic boundary conditions are applied in the direction of w. To allow the interlayer distance to be adjusted freely, the free boundary conditions are applied in the transverse direction (z-direction). The optimized C–B–N parameters of the Tersoff potential given in ref. [36] are used to describe the atomic covalent interactions between C, N, and B atoms, because they have been proven to be of high quality in evaluating thermal transport properties [4,25,37]. Correspondingly, the Lennard–Jones (LJ) potential is used to build the weak vdW force between different in-plane Gr/h-BN heterostructure layers and is adopted as the following:

V ( r ) = 4χ ε

 12 σ r

−

 σ 6  r

,

(1)

where r is the distance between two atoms; ɛ and σ are the energy and distance constants, respectively. The parameter χ is used to control the coupling strength between different in-plane Gr/hBN heterostructure layers. The LJ parameters are listed in Table S1, which are calculated from a universal force field (UFF) model ˚ We have [38]. The cutoff distance of the LJ potential is set as 12 A. tested the effect of using different cutoff distances on simulation results in the Supporting Information (see Table S2). The ITC of the multilayer-mixed Gr/h-BN heterostructures is calculated by using the non-equilibrium molecular dynamics (NEMD) method, and the equations of atomic motion are integrated with a time step of 0.5 fs, ensuring good energy conservation. Each initial structure is first equilibrated in the constant atom number, volume, and temperature ensemble (NVT ensemble) with the Nosé–Hoover thermostat at 300 K for 2 ns (4 million time steps). Then, the equilibrated system is switched to an NVE (constant volume and no thermostat) ensemble. Meanwhile, we generate the non-equilibrium heat flux J by applying the Bussi–Donadio– Parrinello (BDP) thermostat (the corresponding parameter settings and reasons behind the selection of the BDP thermostat can be found in the Supporting Information) to the heat source and sink, with a higher temperature of 330 K and a lower temperature of 270 K, respectively, as schematically shown in Fig. 1d. When the steady state is achieved, the heat flux J can be calculated by continuously adding or removing energy in the thermostated regions at an energy exchanged rate dE/dt and then be expressed as follows:

J=

dE/dt , A

(2)

where A is the cross-sectional area of the simulation system in the direction perpendicular to heat transport. The thickness of monolayer in-plane Gr/h-BN heterostructure is chosen as 0.335 nm and that of MIGHH is the product of 0.335 nm and the layer number. The non-equilibrium heat flux is imposed by using the local thermal baths for 6 ns in total, and all the systems reach a steady state after 3 ns. Therefore, in the production stage, we sample the local temperatures and the accumulated energy exchange (see Fig. 2a) between the system and the thermal baths within the last 3 ns. Finally, the time-averaged temperature jump T (c.f. Fig. 2b) is obtained by linearly fitting these data, and the ITC G of

T. Liang, M. Zhou and P. Zhang et al. / International Journal of Heat and Mass Transfer 151 (2020) 119395

3

Fig. 1. The schematic models of monolayer and multilayer in-plane Gr/h-BN heterostructures. (a) Top view for the atomic model of monolayer in-plane Gr/h-BN heterostructure including C–NB and C–BN structures. (b) The atomic model of a multilayer parallel stacked heterostructure (MPSH, bilayer) viewed from the side. (c) Side view for the simplified model of MPSH (four-layer), where black and red lines represent Gr and h-BN, respectively. (d) Perspective view for the atomic model of multilayer staggered ˚ generating a heat stacked heterostructure (MSSH, bilayer). Heat source and sink are assigned to both ends of the simulation systems (excluding the fixed layers of 15.5 A), flux across the interface. (e) Side view for the simplified model of MSSH (four-layer). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 2. (a) The accumulated energy of the BDP thermostat (averaged over the heat source and sink) as a function of the time in steady state. The energy exchanged rate dE/dt is calculated as the slope of the linear fit (red dashed lines). (b) A typical steady-state temperature profile for the bilayer-mixed heterostructure system (MPSH). The temperature can be fitted by a linear function on each side of the interface of the simulation system, excluding the nonlinear region around the heat source or sink, and then the time-averaged temperature jump T can be extracted as the difference between the two linear functions at the interface. All simulation systems along the heat ˚ are divided into N = 40 blocks, where the number of blocks of the heat source and sink is set to 2. flux direction (excluding the fixed layers of 15.5 A)

4

T. Liang, M. Zhou and P. Zhang et al. / International Journal of Heat and Mass Transfer 151 (2020) 119395

Fig. 3. (a) ITC and (b) thermal rectification (TR) ratio of different heterostructure systems as a function of the number of layers. The red and black dotted lines in (a) and (b) are guides for the eyes. Heat flux and temperature jumps of (c) MPSH and (d) MSSH as a function of the number of layers. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

the steady system can be calculated according to Fourier’s law as follows:

G=

J

T

.

(3)

It is worth mentioning that to avoid the interface dislocations of different in-plane Gr/h-BN heterostructure layers during the equilibrium process, resulting in multiple temperature jumps when heat flux is applied, we have frozen some of the atoms at both ends of the system (15.5 A˚ in total) before the equilibrium, as visualized in Fig. 1d. Because these frozen atoms do not participate in thermal transport, the effective dimensions of all simulated sys˚ To prevent the significant boundtems are w × l = 62.12 × 200 A. ary temperature jump between the temperature-controlled parts (i.e., thermostated regions) and the rest parts, the thermostated dimensions should not be too small [39,40]. This temperature jump is generally considered to be the result of thermal boundary resistance, and it results in both larger calculation error and longer simulation time. In our simulations, the length of the heat source and ˚ ensuring a good temperature gradient curve, as sink is set to 10 A, shown in Fig. 2b. For each simulation system, we performed five independent runs and calculated the error bounds in terms of the standard error. The PDOS is a simple method and powerful tool for characterizing phonon activities in materials [41], and it is computed from the Fourier transform [42] of the velocity autocorrelation function

(VACF) of all the atoms as follows:

PDOS(ω ) =



∞ −∞

eiωt VACF(t )dt,

(4)

where PDOS(ω) is the total PDOS at the vibrational frequency ω and VACF(t) is given by the following:

VACF(t ) =

N 1 vi (0 )vi (t ). N

(5)

i=1

vi (t) is the velocity vector for particle i at time t, N is the number of atoms in the system, and the space-average is the average over the particles. The ensemble average (denoted by ) in Eq. (5) is realized by time-averaging over a period of 15 ps, with the sample velocities extracted from the simulation every 5 fs. Because VACF is normalized and dimensionless, the unit of PDOS is 1/THz. For polarized PDOS, the correlation function is calculated as v(x )v(x ) + v(y )v(y ) and v(z)v(z) for the in-plane and out-ofplane PDOS, respectively. To quantitatively the PDOS mismatch of MIGHH on the two sides of the interface, the overlap area factor S is defined as follows:



S=

∞ 0

min {PGr (ω ), Ph−BN (ω )}dω,

(6)

where PGr and Ph−BN are the two PDOS of Gr and h-BN, respectively.

T. Liang, M. Zhou and P. Zhang et al. / International Journal of Heat and Mass Transfer 151 (2020) 119395

5

Fig. 4. Snapshots of von Mises stress distribution in different heterostructure systems. (a) Von Mises stress distribution in the whole configuration of C–NB. The average stress for the atoms in fixed layers is relatively high because these atoms have not been relaxed. (b) Von Mises stress distribution in the C–BN near the interface. The snapshot is the representative section of C–BN because the whole configuration is too long to be viewed in detail. (c) Top and (d) side views for von Mises stress distribution in MPSH (bilayer).

Calculating the PPR is another effective method of providing insight into the phonon activities, especially for quantitatively describing the phonon localization effect. Without lattice dynamic calculation, the PPR at any temperature can be calculated directly from the MD simulations, thereby implicitly including the all-order of anharmonic scatterings [41], and is defined as [43]:

1 PPR(ω ) = N

 i



PDOSi (ω )2 i

2

PDOSi (ω )4

,

(7)

where N is the total number of atoms, and PDOSi (ω) is the singleparticle PDOS obtained by calculating a single-particle VACF. 3. Results and discussion 3.1. ITC of multilayer in-plane Gr/h-BN heterostructures First, we focus on the layer effect on the ITC of the MIGHH by varying the layer number n from 2 to 6. From Fig. 2 and Eq. 3, we obtained the values of ITC of different heterostructure systems at room temperature (300 K), which are shown in Fig. 3a. The results show that the ITC of C–NB is larger than the C–BN, with the values of 4.35 ± 0.12 and 3.54 ± 0.07 GW m−2 K−1 , respectively, which are consistent with some previous studies [4,44,45]. However, the results are about two to three times smaller than those of other studies [27,46–48], or even smaller. The discrepancy stems from the differences in simulation dimensions, boundary conditions, bonding types, especially the thermostatting methods that generate heat flux [40,49]. Fig. 3a shows that the ITC of MSSH

and MPSH decrease with increasing layer number n and both saturate at n = 3. Nevertheless, both the saturation values are still slightly larger than the ITC value of monolayer in-plane Gr/h-BN heterostructure, indicating that the MIGHH is more conducive to interfacial thermal transport. As mentioned, the multilayer structures of Gr and h-BN suppress the phonon thermal transport because the out-of-plane flexural phonon modes are restricted by the interlayer vdW interactions [28–30]. Therefore, the behavior of MIGHH to enhance the ITC is very confusing, and its mechanisms are analyzed later. Another striking feature observed here is that the ITC value saturated by MSSH is larger than that saturated by MPSH, suggesting that a better interfacial thermal transport effect can be obtained by changing the stacking angle of MPSH. This result is reasonable because it has been proven in many studies [25,37,50] that the magnitude of heat flux from h-BN to Gr domain is higher than that in the opposite direction, indicating a clear thermal rectification (TR) performance in the monolayer in-plane Gr/h-BN heterostructure. In support of this, we further calculate the TR ratio in different heterostructure systems via the following:

TR(% ) =

JH − JL × 100%, JL

(8)

where JH and JL are the high and low heat currents obtained by swapping the temperature bias of the system. When the temperature bias is inverted, the temperature difference between thermal reservoirs remains the same. Here, the JH and JL are the magnitudes of the heat flux in the directions from the right to left region and the left to right region (see Fig. 1b and d), respectively. Note

6

T. Liang, M. Zhou and P. Zhang et al. / International Journal of Heat and Mass Transfer 151 (2020) 119395

Fig. 5. (a) The variation of von Mises stress at different locations in MPSH (bilayer) along the y-direction. (b) The dependence of stress concentration factor Vs on the number of layers in different heterostructure systems.

that the heat flux used to calculate the ITC in this study is from left to right, and the TR ratio is calculated only in this subsection. Fig. 3b illustrates the TR ratios for the different heterostructure systems. Definitely, in the monolayer in-plane Gr/h-BN heterostructure, heat flux preferentially flows from h-BN to Gr domain, demonstrating a pronounced TR behavior. In particular, we collected the average values of heat flux J and temperature jumps T in different heterostructure systems (heat flux from left to right), which are visualized in Fig. 3c and d. By visualizing their magnitude changes, we can better understand the underlying mechanisms of ITC changes in different systems. When the applied heat flux direction is from left to right, the heat flux of evennumbered layers in MSSH flows from h-BN to Gr, resulting in the magnitude of heat flux in MSSH being larger than that in MPSH, as shown in Fig. 3c and d. This also directly leads to the heat flux in MSSH surges at n = 2. Hence, the reason for contributing to this striking feature is as follows: because of the difference in the structural arrangement between MPSH and MSSH (see Fig. 1b–e), the TR behavior leads to a larger heat flux in the MSSH configuration, further leading to the improvement of ITC. Furthermore, the TR in MPSH decreases first and then increases as n increases and finally converges, and the convergence value is smaller than the TR in the monolayer heterostructure. The main reason for this behavior is that the heat flux in MPSH is suppressed as n increases, as depicted in Fig. 3c. Besides, we find that the TR of MSSH changes non-monotonically (inverted V-shaped) with the increase of n, and the TR effect disappears (i.e., TR ratio can be considered statistically zero) when n is even. This is because when n = even, the structure of the MSSH is morphologically identical at both ends of the interface, and at n = odd, the TR in MSSH only comes from the role of a monolayer in-plane Gr/h-BN heterostructure. In other words, the TR in MSSH is physically rooted from the morphological asymmetry. Thus, in Fig. 3b, it is intuitive to expect a lower TR in the MSSH than MPSH, and for the even-layered MSSH, the TR is zero because of the morphological symmetry. The presented results indicate that it is difficult to achieve both high ITC and excellent TR performance in MIGHH at the same time. The following section discusses the physical origin of MIGHH’s ITC saturation at n = 3 and the reason why the saturation values are larger than the ITC value of the monolayer in-plane Gr/h-BN heterostructure. The interfacial thermal transport properties can be

characterized laterally by the stress field, which has been demonstrated in many studies [4,44,46,48]. Accordingly, we calculate the initial stress distribution in different heterostructure systems after they are fully relaxed. The von Mises stress on each atom is calculated as follows:

σm =

1 2

(σ11 − σ22 )2 + (σ11 − σ33 )2 + (σ33 − σ11 )2

+6



2 2 2 + σ23 + σ31 σ12

 12

(9)

where σ 11 , σ 22 , and σ 33 are the normal stress along the x, y, and z directions, respectively; σ 12 , σ 23 , and σ 31 are the shear stress. From Fig. 4, a clear stress concentration is observed near the interface in different heterostructure systems because the lattice mismatch between Gr and h-BN and is more severe in the Gr domain because of its higher in-plane stiffness [51]. In addition, the average stress near the interface of C–NB is lower than that of C– BN, and their heat flux magnitudes are 113.80 (±0.67) × 109 and 104.77 (±1.59) × 109 W m−2 (also see Fig. 3c and d), respectively. In the structure of C–NB and C–BN, the thermal resistance of other parts except the interface is completely the same. This means that the mitigating/severe interfacial stress level of the C–NB/C–BN intensifies/constrains vibrations of atoms near the interface, thus facilitating/obstructing the transfer of heat flux across the interface. Therefore, the analysis of interfacial stress levels can help us gain qualitative insight into the thermal transport across the interface of different heterostructure systems. Unfortunately, when comparing the monolayer in-plane Gr/h-BN heterostructure with the MPSH (bilayer) (Fig. 4a and c, respectively), the results show that the difference in stress distribution near the interface between the two is not obvious enough. Thus, to further analyze the contribution of the interfacial stress field, the stress concentration factor Vs is calculated by the following:

σA , (10) (σB + σC )/2 where σ A , σ B and σ c are the averaged stress extracted from the

Vs =

stress concentration field near the interface (Region A in Fig. 4a) and far from the stress concentration field (Regions B and C in Fig. 4a), respectively. Fig. 5a intuitively shows the numerical trend of stress distribution in MPSH (bilayer) along the y-direction. In Fig. 5b, the Vs

T. Liang, M. Zhou and P. Zhang et al. / International Journal of Heat and Mass Transfer 151 (2020) 119395

7

Fig. 6. (a)–(c) Out-of-plane and in-plane PDOS in different heterostructure systems. (d) A schematic diagram of the regions selected to calculate PDOS. To avoid the interference from the fixed layers, we chose the regions on both sides of the interface and away from the fixed layers. (e) The dependences of overlap of out-of-plane (So ) and in-plane (Si ) PDOS on the number of layers in different heterostructure systems.

of C–NB is lower than that of C–BN, which is consistent with the larger ITC in C–NB. In addition, the Vs of MIGHH dramatically decreases at n = 2, and the ITC of MIGHH is the maximum at n = 2 (see Fig. 3a). These self-consistent results are because of the reduction in the stress concentration factor Vs caused by the effect of the vdW force to intensify the vibration of atoms near the interface, thereby promoting the thermal transport across the interface. Nevertheless, as the number of layers (n > 2) increases, the declining trends of Vs in MPSH and MSSH slow down and even slightly increase. The possible reason for this anomaly is that the atomic vibration near the interface responds rapidly to the interactions of the vdW force at the initial stage, but the response slows down as the n increases. Eventually, the Vs in MPSH and MSSH both show a downward trend. In Fig. 5b, the Vs values in MSSH are larger than that in MPSH, which is caused by the structural deformation of MSSH after relaxation. For a more detailed explanation, see the Supporting Information (Table S3).

The analyses of the stress concentration factor show that as the number of layers increases, the heat flux through the interface should increase gradually, resulting in the gradual increase of ITC of multilayer heterostructures. However, this is inconsistent with what we have observed, so it is not sufficient to focus only on the stress distribution at the interface. Next, we further explore the intrinsic mechanism of the abnormal changes of ITC with the increase of the number of layers by understanding the phonon activities in depth. According to Eqs. (4)–(6), the popular PDOS and its overlap representing the phonon coupling degree are obtained by changing the speed random number and performing three independent calculations. Because of the complexity of the MSSH configuration, it is difficult to separate the PDOS of Gr and h-BN, so we only calculate the PDOS of the MPSH as a typical representative. Fig. 6a–c clearly show that the out-of-plane PDOS of Gr and hBN in MPSH have different degrees of decline in the low-frequency region (0–10 THz) compared with that in the monolayer in-plane

8

T. Liang, M. Zhou and P. Zhang et al. / International Journal of Heat and Mass Transfer 151 (2020) 119395

Fig. 7. (a) Phonon participation ratio of monolayer Gr/h-BN heterostructure and pure Gr and h-BN configurations. (b) and (c) Comparison of phonon participation ratio in MIGHH. The 0.4 value of participation rate is indicated by a black horizontal dotted line. (d) The average participation ratio in different heterostructure systems as a function of the number of layers. Similar to the PDOS calculation, the average phonon participation ratio is obtained by performing three independent calculations.

Gr/h-BN heterostructure. In addition, as the number of layers increases, the decreasing region increases, directly leading to the reduction in out-of-plane overlap So (see in Fig. 6e). Actually, in the MIGHH, the out-of-plane flexural phonon modes of Gr and h-BN are restricted by the interlayer vdW interactions [28–30], causing reduced PDOS of Gr and h-BN in the low-frequency region and thus resulting in the reduction of So . From the perspective of the phonon coupling theory, the reduced So means the weakening of phonon coupling between Gr and h-BN [25,26,37], leading to the deterioration of thermal transport performance in MIGHH. In comparison with the out-of-plane PDOS resulting in a large drop in overlap, the in-plane PDOS, which has been shown in many studies [26,37,52,53] to play a minor role in phonon coupling, does not many differences, as seen in Fig. 6a–c. Although, as shown Fig. 6e, the overlap (Si ) of in-plane PDOS in different heterostructure systems decreases monotonically as the n increases, the decreased amplitude is much smaller than that of So . This proves that the out-of-plane acoustic phonons play a dominant role in the thermal transport of MIGHH. Through the analysis of PDOS, it is known that both So and Si show a monotonic downward trend, indicating that the phonon coupling between Gr and h-BN in MIGHH gradually weakens (the

weakening of out-of-plane phonon coupling is the main responsibility), leading to the deterioration of interfacial thermal transport performance in MIGHH. This is consistent with the result that the thermal conductivity of MIGHH’s parent materials (e.g., multilayer Gr and multilayer h-BN) decreases as the number of layers increases. Combined with the previous analyses of stress concentration at the interface, we conclude that as the number of layers increases, the interfacial thermal transport in MIGHH is promoted but also suppressed. To further capture the phonon activities (especially for the degree of phonon localization) and their effects on the ITC in different heterostructure systems, we calculate the PPR using Eq. (7). It measures the fraction of atoms participating in a given mode and can effectively indicate the localized modes with O (1/N) and delocalized modes with O (1) [54]. Therefore, the PPR is widely used to reveal the phonon confinement, and TR behavior in different Gr/hBN heterostructure systems [25,37,50]. In our analyses, PPR < 0.4 is taken as the criterion to ensure the delocalization of most phonon modes. Fig. 7a shows an overall reduced PPR in monolayer Gr/hBN heterostructure (compared with the pure Gr and h-BN) because of the mismatch of the intrinsic phonon modes of Gr and h-BN. Despite this, most of the phonon modes in the monolayer

T. Liang, M. Zhou and P. Zhang et al. / International Journal of Heat and Mass Transfer 151 (2020) 119395

9

Fig. 8. The typical spatial distribution of the phonon localized modes in different heterostructure systems. The intensity of the localization is depicted according to the color bar.

Gr/h-BN heterostructure have PPR greater than 0.4, showing a characteristic of delocalized mode. However, Fig. 7b–d show that the PPR of MPSH and MSSH decreases with increased n, whether in low-frequency phonons (0–10 and 20–25 THz) or high-frequency phonons (50–60 THz). In other words, as n increases, more and more phonon modes in MPSH and MSSH have PPR less than 0.4, showing a characteristic of localized mode. The contribution to thermal transport mainly comes from delocalized modes rather than localized modes. Therefore, these analyses show that increasing the number of layers induces a greater degree of phonon localization in MIGHH, causing interfacial phonon thermal transport to be blocked. In fact, the deterioration of the thermal transport in the parent materials of MIGHH with increased n is the direct cause of the obstruction of interfacial phonon transport in MIGHH. This result is consistent with the conclusion of PDOS analyses. To get a better physical picture about the phonon localization, we further provide the information about the spatial distribution of a specific mode ( ∈ PPR < 0.4) by [55]:

φi α ,  =

1 N

 PDOSiα dω .   1≤i≤N  PDOSiα d ω

(11)

A larger value of φ iα ,  indicates stronger localization of modes  on the ith atom. Both the in-plane and out-of-plane phonon modes (α = x, y, and z) are considered here. Fig. 8 shows the spatial distribution of the phonon localized modes in different heterostructure systems at 300 K. Obviously, the intensity of the phonon localized modes on the Gr domain is greater than that on the h-BN domain. In other words, the Gr domain is the main bottleneck of interfacial phonon thermal transport channels, which has been confirmed in the studies of thermal rectification [25,37]. Moreover, Fig. 8 shows that as n increases, the intensity of the phonon localized modes in MPSH and MSSH gradually increases, i.e., more phonon modes are converted from delocalization to localization both in Gr and h-BN domains, thus causing a gradual blockage of phonon transmission channels. Furthermore, it can be seen from Fig. 8 that the localized modes in the Gr and h-bn domains are both somewhat far from the interface, indicating that the channel obstructing phonon transport mainly comes from the region far from the interface. When comparing MPSH and MSSH with the same number of layers (see a3 and b3 in Fig. 8), the results show that phonon localization is more severe in MSSH

10

T. Liang, M. Zhou and P. Zhang et al. / International Journal of Heat and Mass Transfer 151 (2020) 119395

Fig. 9. (a) The ITC, (b) stress concentration Vs , (c) overlap of out-of-plane and in-plane PDOS, and (d) average participation ratio of MPSH (bilayer) as a function of the coupling strength χ . The black horizontal dotted line in (a) indicates that ITC is a constant within the range of the error bars.

including Gr and h-BN domains, which may be attributed to the disrupted geometric symmetry along the z-direction in the MSSH configuration. Through the analyses of Vs , overlap of PDOS and PPR, we obtain the intrinsic mechanisms of the abnormal changes of ITC. When the number of layers increases, the competition between the enhanced vibration at the interface (confirmed by the Vs analyses) and the suppressed phonon interfacial transport (confirmed by the overlap of PDOS and PPR analyses) results in these abnormal phenomena. When n = 2, the enhancement of atomic vibration at the interface is dominant, so the ITC of MPSH and MSSH is the largest. When n = 3, suppressed interfacial phonon transport (mainly comes from the region far from the interface) dominates, so the ITC decreases slightly. As the value of n increases, the competition between the two has reached a dynamic equilibrium; thus, the ITC is finally saturated. In addition, because of the enhanced effect of atomic vibration, the saturation values are larger than the ITC value of the monolayer in-plane Gr/h-BN heterostructure, providing enlightenment for the application of MIGHH as a thermal interface material. 3.2. Effect of the coupling strength between layers We know that the vdW force coupling strength between layers in the vertical heterostructures can greatly affect the thermal

transport properties, which has been revealed in many studies [7,47,52]. Coincidentally, in MIGHH, the vibration of the atoms near the interface is stimulated by the vdW force, so it is necessary to explore the effect of coupling strength on the ITC of MIGHH. To this end, we adjust the coupling strength by changing the parameter χ from 1 to 10. For convenience, we only consider here the effect of different coupling strength on the ITC of MPSH (bilayer), and the conclusions obtained are also true in the MSSH configuration. From Fig. 9a, it can be seen that as the coupling strength χ increases, the ITC of MPSH (bilayer) first remains almost constant from χ = 1 to χ = 4, and then rises to the peak and converges at χ = 6. Obviously, the ITC of MPSH (also for MSSH) can be enhanced by controlling the coupling strength χ . To gain insights into the underlying mechanisms for this phenomenon, the stress concentration factor Vs , overlap of PDOS, and PPR analyses are also performed here. Figure S2 also visualizes the changes in heat flux and temperature jumps with different coupling strengths. Moreover, the trend of Vs decreasing nonlinearly with the χ increasing is shown in Fig. 9b. There is an abrupt point (i.e., χ = 4) for the change of Vs , indicating that when χ > 4, the atomic vibration near the interface is sharply intensified. Meanwhile, we calculate the overlap of out-of-plane (So ) and in-plane (Si ) PDOS of MPSH (bilayer), as shown in Fig. 9c. Obviously, as the coupling strength increases, the Si , which plays a minor role, increases monotonically,

T. Liang, M. Zhou and P. Zhang et al. / International Journal of Heat and Mass Transfer 151 (2020) 119395

while the So , which plays a major role, almost decreases gradually. In addition, the decline in phonon average participation ratio (see Fig. 9d) is in good agreement with the trend of So . The detailed information on the spatial distribution of the phonon localized modes in MPSH (bilayer) is also presented here, as illustrated in Fig. S3. The results show that as the χ increases, the phonon localization in the h-BN domain is progressively more severe than that in the Gr domain, indicating that the effect of enhanced coupling strength χ on the phonon modes in h-BN is more significant. The increasingly severe phonon localization is located relatively far from the interface in both the Gr and h-BN domains. These selfharmonized results show that as χ increases from 1 to 4, the enhanced atomic vibration at the interface and the suppressed interfacial phonon transport are in a dynamic equilibrium, so the ITC of MPSH remains a constant. As χ continues to increase, the atomic vibration of the interface sharply intensifies and becomes dominant, leading to the increased ITC, and then the competition returns to equilibrium again; thus, the ITC finally reaches saturation.

11

sualization, Investigation, Data curation. Daoguo Yang: Resources, Project administration. Acknowledgments The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors acknowledge the financial support provided by National Natural Science Foundation of China (Project No. 51506033), Guangxi Natural Science Foundation (Grant Nos. 2017JJA160108 and 2019JJG160011), and Guangxi Colleges and Universities Program of Innovative Research Team and Outstanding Talent. Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.ijheatmasstransfer. 2020.119395.

4. Conclusions To summarize, in this paper, we construct two types of multilayer in-plane Gr/h-BN heterostructures (including MPSH and MSSH configurations) and systematically investigate the effect of the number of layers on their ITC by performing NEMD simulations. It is interesting to find that the ITC of MPSH and MSSH decrease with increasing layer number n and both saturate at n = 3. Nevertheless, both the saturation values are still larger than the ITC value of monolayer in-plane Gr/h-BN heterostructure, suggesting that the MIGHH is more conducive to interfacial thermal transport. Based on the analyses of stress concentration factor, overlap of PDOS, and phonon participation ratio, we conclude that the enhanced interfacial thermal transport originates from the strengthening of atomic vibration near the interface (confirmed by the mitigation of stress concentration); the suppression of interfacial thermal transport results from the weakening of phonon coupling between Gr and h-BN and the increase of phonon localized modes, and ultimately, the competition between the two leads to these phenomena. Moreover, by comparing the ITC in MPSH and MSSH, the results show that a higher ITC can be obtained by changing the stacking angle of MPSH, because of the thermal rectification behavior. Finally, we explore the effect of coupling strength on the ITC of MPSH (bilayer) and find that the coupling strength can be effective in enhancing the ITC in MIGHH. We anticipate that our findings not only can shed some light on exploring the thermal transport properties of other in-plane heterostructures, such as phosphorene/graphene [9], MoS2 /graphene [10], and silicene/graphene [11], but also provide a new design idea for the applications of multilayer in-plane Gr/h-BN heterostructures in thermal management and thermoelectric devices. Declaration of Competing Interest The authors declared that they have no conflicts of interest to this work. We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted. CRediT authorship contribution statement Ting Liang: tion, Writing Data curation, administration,

Conceptualization, Methodology, Software, Validaoriginal draft. Man Zhou: Investigation, Software, Writing - review & editing. Ping Zhang: Project Supervision, Funding acquisition. Peng Yuan: Vi-

References [1] A.K. Geim, I.V. Grigorieva, Van der Waals heterostructures, Nature 499 (7459) (2013) 419–425. [2] Z. Liu, L.L. Ma, G. Shi, W. Zhou, Y.J. Gong, S.D. Lei, et al., In-plane heterostructures of graphene and hexagonal boron nitride with controlled domain sizes, Nat. Nanotechnol. 8 (2) (2013) 119–124. [3] J.E. Barrios-Vargas, B. Mortazavi, A.W. Cummings, R. Martinez-Gordillo, M. Pruneda, L. Colombo, et al., Electrical and Thermal transport in coplanar polycrystalline graphene-hBN heterostructures, Nano Lett. 17 (3) (2017) 1660–1664. [4] X. Liu, G. Zhang, Y. Zhang, Topological defects at the graphene/h-BN interface abnormally enhance its thermal conductance, Nano Lett. 16 (8) (2016) 4954–4959. [5] Y. Gao, Q. Liu, B. Xu, Lattice mismatch dominant yet mechanically tunable thermal conductivity in bilayer heterostructures, Acs Nano 10 (5) (2016) 5431–5439. [6] L. Britnell, R.V. Gorbachev, R. Jalil, B.D. Belle, F. Schedin, A. Mishchenko, et al., Field-effect tunneling transistor based on vertical graphene heterostructures, Science 335 (6071) (2012) 947–950. [7] B. Liu, J.A. Baimova, C.D. Reddy, A.W. Law, S.V. Dmitriev, H. Wu, et al., Interfacial thermal conductance of a silicene/graphene bilayer heterostructure and the effect of hydrogenation, Acs Appl. Mater. Inter. 6 (20) (2014) 18180–18188. [8] L. Liu, J. Park, D.A. Siegel, K.F. McCarty, K.W. Clark, W. Deng, et al., Heteroepitaxial growth of two-dimensional hexagonal boron nitride templated by graphene edges, Science 343 (6167) (2014) 163–167. [9] X.J. Liu, J.F. Gao, G. Zhang, Y.W. Zhang, Design of phosphorene/graphene heterojunctions for high and tunable interfacial thermal conductance, Nanoscale 10 (42) (2018) 19854–19862. [10] X.J. Liu, J.F. Gao, G. Zhang, Y.W. Zhang, MoS2-graphene in-plane contact for high interfacial thermal conduction, Nano Res. 10 (9) (2017) 2944–2953. [11] B. Liu, J.A. Baimova, C.D. Reddy, S.V. Dmitriev, W.K. Law, X.Q. Feng, et al., Interface thermal conductance and rectification in hybrid graphene/silicene monolayer, Carbon 79 (2014) 236–244. [12] K. Tielrooij, N.C.H. Hesp, A. Principi, M.B. Lundeberg, E.A.A. Pogna, L. Banszerus, et al., Out-of-plane heat transfer in van der Waals stacks through electron-hyperbolic phonon coupling, Nat. Nanotechnol. 13 (1) (2018) 41–46. [13] D. Jariwala, T.J. Marks, M.C. Hersam, Mixed-dimensional van derWaals heterostructures, Nat. Mater. 16 (2) (2016) 170–181. [14] Y. Gao, B. Xu, On the generalized thermal conductance characterizations of mixed one-dimensional–two-dimensional van der Waals heterostructures and their implication for pressure sensors, Acs Appl. Mater. Inter. 10 (16) (2018) 14221–14229. [15] Y. Gao, B. Xu, van der Waals graphene Kirigami heterostructure for strain-controlled thermal transparency, Acs Nano 12 (11) (2018) 11254–11262. [16] Y.J. Gong, J.H. Lin, X.L. Wang, G. Shi, S.D. Lei, Z. Lin, et al., Vertical and in– plane heterostructures from WS2/MoS2 monolayers, Nat. Mater. 13 (12) (2014) 1135–1142. [17] S. Dai, Q. Ma, M.K. Liu, T. Andersen, Z. Fei, M.D. Goldflam, et al., Graphene on hexagonal boron nitride as a tunable hyperbolic metamaterial, Nat. Nanotechnol. 10 (8) (2015) 682–686. [18] C. Dean, A.F. Young, L. Wang, I. Meric, G.H. Lee, K. Watanabe, et al., Graphene based heterostructures, Solid State Commun. 152 (15) (2012) 1275–1282. [19] M.P. Levendorf, C. Kim, L. Brown, P.Y. Huang, R.W. Havener, D.A. Muller, et al., Graphene and boron nitride lateral heterostructures for atomically thin circuitry, Nature 488 (7413) (2012) 627–632. [20] L.J. Ci, L. Song, C.H. Jin, D. Jariwala, D.X. Wu, Y. Li, et al., Atomic layers of hybridized boron nitride and graphene domains, Nat. Mater. 9 (5) (2010) 430–435.

12

T. Liang, M. Zhou and P. Zhang et al. / International Journal of Heat and Mass Transfer 151 (2020) 119395

[21] P. Zhang, P. Yuan, X. Jiang, S.P. Zhai, J.H. Zeng, Y.Q. Xian, et al., A theoretical review on interfacial thermal transport at the nanoscale, Small 14 (2) (2018) 1702769. [22] J. Wang, X. Mu, X. Wang, N. Wang, F. Ma, W. Liang, et al., The thermal and thermoelectric properties of in-plane C–BN hybrid structures and graphene/h-BN van der Waals heterostructures, Mater. Today Phys. 5 (2018) 29–57. [23] Y. Wang, N. Xu, D. Li, J. Zhu, Thermal properties of two dimensional layered materials, Adv. Funct. Mater. 27 (19) (2017) 1604134. [24] J. Zhang, Y. Hong, Y. Yue, Thermal transport across graphene and single layer hexagonal boron nitride, J. Appl. Phys. 117 (13) (2015) 134307. [25] X. Chen, Z. Xie, W. Zhou, L. Tang, K. Chen, Thermal rectification and negative differential thermal resistance behaviors in graphene/hexagonal boron nitride heterojunction, Carbon 100 (2016) 492–500. [26] X. Chen, J. Liu, Z. Peng, D. Du, K. Chen, A wave-dominated heat transport mechanism for negative differential thermal resistance in graphene/hexagonal boron nitride heterostructures, Appl. Phys. Lett. 110 (9) (2017) 91907. [27] F. Liu, R. Zou, N. Hu, H. Ning, C. Yan, Y. Liu, et al., Enhancement of thermal energy transport across the graphene/h-BN heterostructure interface, Nanoscale 11 (9) (2019) 4067–4072. [28] D.A. Broido, N. Mingo, L. Lindsay, Flexural phonons and thermal transport in multilayer graphene and graphite, Phys. Rev. B 83 (23) (2011) 235428. [29] Z. Wei, Z. Ni, K. Bi, M. Chen, Y. Chen, In-plane lattice thermal conductivities of multilayer graphene films, Carbon 49 (8) (2011) 2653–2658. [30] D.A. Broido, L. Lindsay, Theory of thermal transport in multilayer hexagonal boron nitride and nanotubes, Phys. Rev. B 85 (3) (2012) 35436. [31] I. Jo, M.T. Pettes, J. Kim, K. Watanabe, T. Taniguchi, Z. Yao, et al., Thermal conductivity and phonon transport in suspended few-layer hexagonal boron nitride, Nano Lett. 13 (2) (2013) 550–554. [32] Y. Gao, Y. Zhang, P. Chen, Y. Li, M. Liu, T. Gao, et al., Toward single-layer uniform hexagonal boron nitride–graphene patchworks with zigzag linking edges, Nano Lett. 13 (7) (2013) 3439–3443. [33] M. Liu, Y. Li, P. Chen, J. Sun, D. Ma, Q. Li, et al., Quasi-freestanding monolayer heterostructure of graphene and hexagonal boron nitride on Ir(111) with a zigzag boundary, Nano Lett. 14 (11) (2014) 6342–6347. [34] G. Zhang, Y. Zhang, Z. Ong, Controlling the thermal conductance of graphene/h-BN lateral interface with strain and structure engineering, Phys. Rev. B 93 (7) (2016) 75406. [35] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117 (1) (1995) 1–19. [36] J.B. Haskins, C. Sevik, T. Çag˘ ın, A. Kınacı, Thermal conductivity of BN-C nanostructures, Phys. Rev. B 86 (11) (2012) 115410. [37] Y. Gao, B. Xu, Controllable interface junction, in-plane heterostructures capable of mechanically mediating on-demand asymmetry of thermal transports, Acs Appl. Mater. Inter. 9 (39) (2017) 34506–34517. [38] A.K. Rappé, C.J. Casewit, K.S. Colwell, W.A. Goddard III, W.M. Skiff, UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations, J. Am. Chem. Soc. 114 (25) (1992) 10024–10035.

[39] J. Chen, J. Wang, B. Li, J. Jiang, Edge states induce boundary temperature jump in molecular dynamics simulation of heat conduction, Phys. Rev. B 80 (5) (2009) 52301. [40] J. Chen, G. Zhang, B. Li, Molecular dynamics simulations of heat conduction in nanostructures: effect of heat bath, J. Phys. Soc. Jpn. 79 (7) (2010) 74604. [41] H. Bao, J. Chen, X. Gu, B. Cao, A review of simulation methods in micro/nanoscale heat conduction, ES Energy Environ. 1 (2018) 16–55. [42] A. PASKIN, J.M. DICKEY, Computer simulation of the lattice dynamics of solids, Phys. Rev. 188 (3) (1969) 1407–1418. [43] G.C. Loh, E.H.T. Teo, B.K. Tay, Phonon localization around vacancies in graphene nanoribbons, Diam. Relat. Mater. 23 (2012) 88–92. [44] Y. Zhang, Q. Pei, C. Wang, C. Yang, Y. Zhang, Interfacial thermal conductance and thermal rectification of hexagonal BCnN/graphene in-plane heterojunctions, J. Phys. Chem. C 122 (39) (2018) 22783–22789. [45] S. Lu, A.J.H. McGaughey, Thermal conductance of graphene/hexagonal boron nitride heterostructures, J. Appl. Phys. 121 (11) (2017) 115103. [46] M. Li, B. Zheng, K. Duan, Y. Zhang, Z. Huang, H. Zhou, Effect of defects on the thermal transport across the graphene/hexagonal boron nitride interface, J. Phys. Chem. C 122 (26) (2018) 14945–14953. [47] T. Li, Z. Tang, Z. Huang, J. Yu, Substrate effects on the thermal performance of in-plane graphene/hexagonal boron nitride heterostructures, Carbon 130 (2018) 396–400. [48] T. Liang, P. Zhang, P. Yuan, M. Zhou, S. Zhai, Interfacial Thermal conductance and thermal rectification across in-plane graphene/h-BN heterostructures with different bonding types, in: Proceedings of the ASME 2019 6th International Conference on Micro/Nanoscale Heat and Mass Transfer, 2019. [49] Z. Li, S. Xiong, C. Sievers, Y. Hu, Z. Fan, N. Wei, et al., Influence of thermostatting on nonequilibrium molecular dynamics simulations of heat conduction in solids, J. Chem. Phys. 151 (23) (2019) 234105. [50] L. Medrano Sandonas, G. Cuba-Supanta, R. Gutierrez, A. Dianat, C.V. Landauro, G. Cuniberti, Enhancement of thermal transport properties of asymmetric Graphene/hBN nanoribbon heterojunctions by substrate engineering, Carbon 124 (2017) 642–650. [51] S. Ciraci, M. Topsakal, Elastic and plastic deformation of graphene, silicene, and boron nitride honeycomb nanoribbons under uniaxial tension: a first-principles density-functional theory study, Phys. Rev. B 81 (2) (2010) 24107. [52] T. Liang, P. Zhang, P. Yuan, S. Zhai, In-plane thermal transport in black phosphorene/graphene layered heterostructures: a molecular dynamics study, Phys. Chem. Chem. Phys. 20 (32) (2018) 21151–21162. [53] Y. Chen, Y. Zhang, K. Cai, J. Jiang, J. Zheng, J. Zhao, et al., Interfacial thermal conductance in graphene/black phosphorus heterogeneous structures, Carbon 117 (2017) 399–410. [54] J. Chen, G. Zhang, B. Li, Remarkable reduction of thermal conductivity in silicon nanotubes, Nano Lett. 10 (10) (2010) 3978–3983. [55] Y. Wang, A. Vallabhaneni, J. Hu, B. Qiu, Y.P. Chen, X. Ruan, Phonon lateral confinement enables thermal rectification in asymmetric single-material nanostructures, Nano Lett. 14 (2) (2014) 592–596.