# A Qualitative Study of the Disorder Effect on the Phonon Transport in a Two-Dimensional Graphene/h-BN Heterostructure Yinong Liu, et al.

May 6, 2022

## Introduction

With the development of low-dimensional materials and the increasing demand for microelectronic devices, effective management of nanoscale heat transport has gradually become an urgent problem to be solved. Traditionally, the thermal conductivity can manipulate by introducing additional defects (Chen et al., 2010; Hao et al., 2011; Zhang et al., 2011; Ding et al., 2015; Feng et al., 2015; Nai et al., 2015), impurities (Chen et al., 2009; Chen et al., 2012), nanoparticles (Maldovan, 2013; Wang et al., 2013; Lin et al., 2016; Mendoza and Chen, 2016; Lu et al., 2021), and ion-intercalation (Qian et al., 2016) into the pristine materials. Additionally, the underlying physical mechanism can be well explained by solving the Boltzmann transport equation (BTE), in which phonons are regarded as incoherent particles (particle nature of phonons).

Recently, periodic nanostructures (Hopkins et al., 2011; Zen et al., 2014; Alaie et al., 2015; Yang et al., 2015; Hu et al., 2016; Ma et al., 2016; Sledzinska et al., 2019; Vasileiadis et al., 2021; Nomura et al., 2022) constructed based on the phonon wave interference (wave nature of phonons), which can modify the phonon dispersion and further reduce the group velocity (Swinteck et al., 2013; Latour et al., 2014; Xiong et al., 2016; Ma et al., 2018), have attracted widespread attention from researchers. Subsequently, the researchers additionally introduced randomness into periodic nanostructures. The results showed that a certain degree of randomness could significantly suppress the thermal conductivity (Hu et al., 2018; Hu et al., 2019). However, the mentioned results above were obtained by calculating a limited number of random structures, mainly because it is impossible to perform global optimization with the current computer resource for all possible configurations (a vast design space).

However, material informatics (machine learning), which applies informatics principles to materials science and engineering to improve the understanding, use, selection, development, and discovery of materials, seems to have brought dawn to the solution of large-scale design space problems. Such as, in some recent studies, Wan et al. used a convolutional neural network method to effectively design the structure of porous with the lowest thermal conductivity (Wan et al., 2020). Wei et al. obtained the result that the thermal conductivity of disordered structures is unexpectedly more significant than that of the periodic system through the efficient genetic algorithm search (Wei et al., 2020). Hu et al. realized the ultimate impedance of coherent heat conduction in van der Waals graphene-MoS2 heterostructures based on Bayesian optimization (Hu et al., 2021). Nevertheless, it is worth noting that the studies mentioned above were all performed in a relatively small design space, even if machine learning (Bayesian optimization) can reduce the number of calculations to a few percent (Ju et al., 2017; Hu et al., 2020). Therefore, to effectively broaden the application scope of machine learning and expand the design space, exploring the typical characteristics of target properties is becoming an urgent task and research hotspot.

In this work, we use non-equilibrium molecular dynamics (NEMD) simulations to investigate the effect of disorder on the phonon transport in a two-dimensional graphene/hexagonal boron nitride (G/h-BN) heterostructure. First, the signature for coherent phonon transport is observed in G/h-BN heterostructure by varying the period length. Furthermore, we introduce five types of disorders to identify the heterostructure’s common feature with low thermal conductivity. The underlying mechanism is further uncovered by performing wave packet simulation. The present findings help to reveal the underlying physical mechanism of phonon transport in two-dimensional heterostructures, and will be instructive for the future synthesis of heterostructures with unique thermal properties.

## Simulation Method

In our simulations, the G/h-BN heterostructure (Figures 1D,E) is composed of the unit cells constructed by distributing h-BN periodically/disorderly into the pristine graphene (Figures 1A,B). Meanwhile, the edge of h-BN is guaranteed to be zigzag, mainly considering that different defect type have different effects on the thermal transport of graphene (Feng et al., 2015). The period length (Lp) and replacement ratio (Rer) are used to characterize the unit cells. Figures 1A,B shows two representative unit cells with different Lp. The replacement ratio is NB/NG, where NB and NG are the numbers of atoms in h-BN and the total number of atoms in pristine graphene, respectively. Here, the replacement ratio fixes at 10%.

FIGURE 1. Schematic figures of G/h-BN lateral heterostructure. (A)(B) Schematic diagrams of the unit cells with

and

. The replacement rate is fixed at 10%. (C) Schematic diagram of the introduction of disorder. (D)(E) Schematic diagram of periodic structure and TU-DS.

Our study’s NEMD simulations are implemented using the LAMMPS package (Plimpton, 1995). The optimized Tersoff potential function describes the covalent-bonds interaction in graphene and h-BN(Sevik et al., 2011; Kınacı et al., 2012). The time step is set as 0.5 fs The fixed and periodic boundary conductions are adopted along the length and width direction in our simulations, respectively. The width is fixed at 7.38 nm, which is enough to eliminate the numerical size effect caused by the periodic boundary condition (Hu et al., 2017). Two Langevin thermostats with temperatures of 310 and 290 K are applied at both ends of the simulation system (red and blue region in Figure 1D) to establish a temperature gradient. The system first runs 50 ps in the NPT ensemble and then relaxes for 2 ns in the NVE ensemble to reach the steady-state. The cumulative energy ΔE added/subtracted to the heat source/sink region and the temperature are recorded for another 5 ns. The energy change per unit time (ΔE/Δt) was obtained by linearly fitting the raw data of the accumulated energy ΔE, which was used to calculate the heat flux

$J=ΔEΔt⋅S$

. Here S is the cross-sectional area of G/h-BN heterostructure. We use Fourier’s Law to calculate the thermal conductivity (κ), κ =—J/

$∇T$

, where

$J$

and

$∇T$

are, respectively, the heat flux and the temperature gradient.

## Results and Discussion

### Thermal Conductivity of G/h-BN Heterostructure

We first study the effective thermal conductivity of G/h-BN heterostructure with different Lp in a finite-size system. Here, the system length fixes to 76.7 nm. Such treatments have been extensively employed in literature (Chang et al., 2008; Liu et al., 2014; Xu et al., 2014) when the thermal transport behavior is not necessarily diffusive. As shown in Figure 2, the thermal conductivity of the periodic structure first decreases and then increases as the increase of

$Lp$

. The competition mechanism between coherent and incoherent phonons in heat conduction leads to this nonlinear trend (Wu and Han, 2022). At smaller

$Lp$

, the heat transport mainly depends on coherent phonons. In contrast, the heat transport mainly contributes from the incoherent phonons at larger

$Lp$

. As

$Lp$

increases, the transition from coherent phonon to incoherent phonon-dominated thermal transport is observed, and a similar trend was found in previous studies (Hu et al., 2018; Hu et al., 2019; Roy Chowdhury et al., 2020). The detailed discussion of the relevant physical mechanisms can be found in the Refs (Latour et al., 2014; Lee et al., 2017; Hu et al., 2018; Xie et al., 2018; Felix and Pereira, 2020). In addition, we further construct five traditional unguided disorder structures (TU-DS) by randomly introducing disorders at each

$Lp$

to explore the disorder’s effect on thermal conductivity. Figure 1C shows the seven center positions of h-BN (denoted 1–7) that can be selected when constructing the TU-DS. The results show that different TU-DS heterostructures exhibit different thermal conductivities, which originate from different degrees of phonon localization caused by different disorders. In addition, introducing the disorder at smaller

$Lp$

can greatly reduce the thermal conductivity (Figure 2). However, as

$Lp$

increases, the effect of disorder on the thermal conductivity is becoming not apparent anymore. For instance, as shown in the inset, when

, the thermal conductivity of the TU-DS heterostructures is reduced to 37% of the corresponding periodic structure. Inversely, when the disorder is introduced into the heterostructure with

$Lp=$

, the thermal conductivity is reduced by only 5.4%. Therefore, to gain further insight into how disorder affects coherent phonon transport and thus reduces thermal conductivity, we performed a qualitative analysis of several different disorder types, with

fixed

in the following simulations.

FIGURE 2. The effect of

$Lp$

on the thermal conductivity of periodic structure and TU-DS at 300 K. The inset shows the relative difference between the thermal conductivity of the TU-DS and the corresponding periodic structure with different

$Lp$

. The relative difference is defined as

$(κperiodic−κTU−DS)/κperiodic×100%$

. The error bars are the standard deviations of five independent simulations with different initial conditions.

### Qualitative Analysis of Disorder Effects

To classify different types of disorder, we first introduced a descriptor, shift distance (SD), to characterize the strength of the disorder. The SD is defined as the deviation length from the current h-BN center (dots of different colors in Figures 3A,C,E) to the original center (denoted as “O” in Figures 3A,C,E). For example, as shown in Figure 3A, the distance between the red point and the original center O is defined as

$a$

(SD =

). In addition to SD, we further propose another descriptor, the direction of disorder, to complement the description of the disorder type. We randomly generated 20 disordered structures (changing the disorder direction and SD) for each specific disorder type in the following qualitative analysis simulations.

FIGURE 3. Qualitative analysis of the effect of disorder on thermal conductivity at 300 K (A)(B) Schematic diagrams of

$Ya$

-DS and

$Ca$

-DS and their corresponding thermal conductivity. (C)(D) Schematic diagrams of

$Ca$

-DS and

$C2a$

-DS and their corresponding thermal conductivity. (E)(F) Schematic diagrams of

$C3a$

-DS and

$X3a$

-DS and their corresponding thermal conductivity.

We first construct the disordered structures by shifting the h-BN with fixed SD =

in the y-direction (

$Ya$

-DS) and cross direction (

$Ca$

-DS), respectively. The possible centers of the h-BN in

$Ya$

-DS and

$Ca$

-DS are denoted as black and red dots, respectively. (Figure 3A). The corresponding thermal conductivity results (Figure 3B) show that the average thermal conductivity of

$Ca$

-DS (

$124.77$

W/m-K) is lower than that of the

$Ya$

-DS

W/m-K). It inspires us that the introduction of disorder in the cross-direction can more effectively hinder phonon transport, thereby further reducing the thermal conductivity of the heterostructures. Therefore, we fix the disorder direction as the cross-direction and increase the SD to

$2a$

(

$C2a$

-DS). The possible centers are denoted as purple points in Figure 3C. The thermal conductivity results show that the

$κ¯C2a−DS$

(

$118.78$

W/m-K) are further reduced compared to the

$κ¯Ca−DS$

(Figure 3D). This result suggests that the larger SD appears more promising to block phonon transport. Along this design route, we introduce two additional disorder types at x (

$X3a$

-DS) and cross-direction (

$C3a$

-DS) with a fixed SD =

$3a$

. As we expected, the

$κ¯C3a−DS$

(

$116.16$

W/m-K) is significantly lower than the

$κ¯X3a−DS$

(

$135.46$

W/m-K), which again confirms that larger SD and introducing disorder at cross-direction can effectively reduce the thermal conductivity. In Figures 3B,D,F, we plot the horizontal dashed line to denote the thermal conductivity of the periodic structure (

$κ¯periodic=339.42$

W/m-K). Our qualitative analysis establishes a rationale between the disorder characteristics (direction and SD) and the thermal conductivity of heterostructure with different disorder types. It manifests that the heterostructure with a larger

$SD$

and disorder in the cross-direction leads to a further reduction in thermal conductivity. Our study focuses on manipulating phonon thermal transport in heterostructures through different disorder types, which provides two characteristics for considering the introduction of the disorder in structural design, bringing new insights into further narrowing the design space.

### Structural Analysis

To understand the underlying physical mechanism of the above results, we conducted phonon wave packet simulations in pristine graphene, periodic structure, TU-DS, and

$C3a$

-DS, which can directly record the wave characteristics of phonon transport behavior (phonon transmission and reflection) in the coordinate space. As shown in Figures 4A–D, we launched the same phonon wave packet from the pristine graphene side with a fixed frequency of 6.11 THz (Hu et al., 2019) and recorded its propagation along the length direction. The length ratio between the pristine graphene and the four targeted structures mentioned above is 1:1. Periodic boundary conditions are used for all directions. Before the wave packet simulation performing, each structure was relaxed with an NPT ensemble at 0 K. The initial wave packet was formed via linear combinations of normal vibration modes as follows:

$μn,α=Re{Aεjα,λeik0(xn−x0)e−(xn−x0)2/η2},(1)$

where

denotes

$αth$

displacement component of the

$n$

th atom in the simulated structure.

$A$

is the amplitude of the wave packet,

$η$

is the spatial width of the wave packet, and

represents the αth eigenvector component of the eigenmode λ for the jth atom in the pristine graphene. To initialize velocities, we added time dependence to Eq. 1 and differentiated it as

$vn,α=Re{−iωλμn,α}.(2)$

FIGURE 4. Phonon wave-packet simulations in different structures. Snapshots of displacement for a TA wave packet in (A) pristine graphene, (B) periodic structure, (C)TU-DS, and (D)

$C3a$

-DS, at 0 K. The black dashed line at 1,240 nm represents the position of the interface. (E) The kinetic energy in the region of 1,240 nm–1,340 nm in length versus time for pristine graphene, periodic structure, TU-DS, and

$C3a$

-DS.

We chose the transverse acoustic (TA) phonon mode of pristine graphene as the excited wave packet with the wavevector

$k0⇀=0.1*b1⇀+0.1*b2⇀$

, where

$b1⇀$

and

$b2⇀$

are two reciprocal basis vectors of the graphene unit cell.

$xn$

is the coordinate of the

$n$

th atom in the structure, and

$x0$

determines the center of the wave packet.

After the initial wave packet is excited, the wave packet continues to run for 50 ps at 0 K under the NVE system. We limited the wave packet simulation to 0 K, mainly because the anharmonic phonons-phonons interaction can be ignored at this temperature, which is beneficial for us to focus only on the wave nature of the phonon transport process. The method of introducing low-temperature conditions in the wave packet simulation has been widely used in previous research works (Wei et al., 2012; Chen et al., 2015; Zhang et al., 2017). Figures 4A–D show the wave packets propagation process in the pristine graphene, periodic structure, TU-DS, and

$C3a$

-DS, respectively. As expected, the phonon can be ballistically transported at a constant speed (group velocity) and maintain its shape in the pristine graphene since the anharmonic phonons-phonons interaction is negligible (Figure 4A).

As shown in Figure 4B, when the wave packet enters the periodic structure, the original wave is divided into a transmitted and a reflected wave after colliding with the interface. In addition, it is worth noting that most waves can transmit through the periodic structure and continue to propagate, exhibiting the wave propagation characteristics in the coordinate space. Interestingly, when the same wave packet enters the TU-DS (Figure 4C), the phonons in the reflected part increase significantly. Furthermore, the transmitted wave gradually dissipates during the propagation process in the TU-DS system, reflecting a solid signal that impedes the phonons’ transport. Moreover, as shown in Figure 4D, the wave packet shows a similar transport behavior when propagating in the

$C3a$

-DS, except that it exhibits a much stronger dissipation. The wave packet simulation results correspond to our thermal conductivity results above, that is, the

$C3a$

-DS has the smallest thermal conductivity. We show that the G/h-BN heterostructure relying on local-specific disorder can impede heat transport at the nanoscale through wave packet simulations.

To further quantify the difference, we also monitored the lattice vibrational kinetic energy of the target region near the interface (1,240–1340 nm) in the above four different structures, which can directly reveal the propagation capability and localization degree of phonons. As shown in Figure 4E, all the energy of the wave-packet can go into the target region in pristine graphene, and it then exits completely (energy decays to zero after 27 ps), reflecting the propagating nature of the phonon. In contrast, when the h-BN structures are introduced into pristine graphene to construct the heterostructures, less energy could enter the target region, and the entering energy decrease sequentially in the periodic, TU-DS, and

$C3a$

-DS heterostructure. This decreasing trend is consistent with the case of thermal conductivity. In addition, it is worth noting that the power in the target region cannot completely decay to zero within the same simulation time, which further confirms the localized nature of the phonon modes in these heterostructures. Moreover, in the

$C3a$

-DS, the most energy stays in the target region, reflecting the most substantial phonons localization. In the disordered G/h-BN heterostructures, the non-periodic distribution of the h-BN induces randomness and thus leads to the localization of the phonons. In most cases, localized strength is difficult to discern. However, we effectively identify the localization difference with the lattice vibrational kinetic energy method, which is a significant feature of this paper. It is worth noting that since the essence of thermal conductivity decrease is the localization of coherent phonons caused by disorder, this effect will not be limited to graphene system, but also applies to other 2D materials and thin films.

## Conclusion

In summary, we have investigated the different types of the disorder effect on the thermal conductivity of the G/h-BN heterostructure via non-equilibrium molecular dynamics simulations The study found that when hexagonal boron nitrides are disorderly distributed in the cross–direction with a larger shift distance in graphene, the thermal conductivity could be reduced by 37%, which is much larger than the random disorder introduction. The wave packet simulation revealed that a larger shift distance in the cross-direction would induce a much more robust phonon localization. Moreover, our research proposes two useful descriptors (direction and shift distance) to effectively reduce the design space of the random structures and gives a more in-depth physical understanding of the disorder effect on phonon transport.

## Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

## Author Contributions

WR, SH and CZ conceived the idea and supervised the entire project. YL performed the MD simulations, supported by SH. YL, WR, SH, and CZ analyzed the data and prepared the manuscript. All authors participated in the discussion.

## Funding

This research was funded in parts by the National Natural Science Foundation of China (Grant Nos. 12105242, 12004242, 11864020, and 61904072), by Yunnan Fundamental Research Project (Grant Nos. 202001AU070047, 202201AT070161, 2019FI002, 202101AS070018, 202001AU070025, and 202101BE070001-049), and by Shanghai Rising-Star Program (No. 21QA1403300), the Yunnan Ten Thousand Talents Plan Young and Elite Talents Project, and the Yunnan Province Computational Physics and Applied Science and Technology Innovation Team.

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.