Thermal rectification of graphene on substrates with inhomogeneous stiffness

Carbon 154 (2019) 81e89

Contents lists available at ScienceDirect

Carbon journal homepage: www.elsevier.com/locate/carbon

Thermal rectification of graphene on substrates with inhomogeneous stiffness Ning Wei a, **, Shanchen Li a, Yingyan Zhang b, Jige Chen c, Yang Chen a, Junhua Zhao a, * a

Jiangsu Key Laboratory of Advanced Food Manufacturing Equipment and Technology, Jiangnan University, 214122, Wuxi, China School of Computing, Engineering and Mathematics, Western Sydney University, Penrith, NSW, 2751, Australia c Shanghai Advanced Research Institute, Chinese Academy of Sciences, Shanghai, 201210, China b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 2 May 2019 Received in revised form 15 July 2019 Accepted 27 July 2019 Available online 27 July 2019

Thermal rectification is of great importance in various thermal management and thermal logic applications, where heat flux runs preferentially in one direction over the opposite one. This phenomenon is mostly reported in previous simulations and experiments by introducing asymmetry shape-tailoring, defects and doping, etc. However, these approaches break the structure of the heat conductor. In this paper, we report thermal rectification induced by modulating the substrate stiffness on the supported graphene, which preserves structural integrity of graphene, making it reusable. Our results show that the heat flux prefers to flow toward the stiffer part of a substrate. The maximum thermal rectification factor can reach ~48% by using the ideal substrate model prediction. Moreover, a wide range thermal rectification factor of 3.5%e22.7% can be achieved by using a sandwich model. These phenomena are explained by phonon vibration spectra and excitation simulations. Our results should be of great help for understanding the modulation of thermal rectification in supported material and shed light on the design of a new thermal logic gate and thermal management devices. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Heat transport properties of graphene have stimulated numerous experimental [1e3] and theoretical studies [4e7], since it possesses an ultrahigh yet easy-tuning thermal conductivity (TC). One can maximize (minimize) thermal conductance for cooling (thermoelectric) applications. Usually, the TC of graphene can be adjusted by doping [8], defects [9e11], tailoring [12e14], strain [15] and substrate [2,3,16e18], etc. When these tuning methods lead to the direction-dependent heat transport, the so-called thermal rectification (TR) occurs. It means that the TC value for the heat flux flowing in one direction is higher than that in the opposite one. The underlying mechanism is that the system allows heat flux in one direction unimpededly while it acts like an insulator when the temperature gradient is reversed [19]. TR phenomenon can be realized by an in-plane asymmetric design, e.g., graphene can be tailored to become asymmetrical

* Corresponding author. ** Corresponding author. E-mail addresses: [email protected] (N. Wei), junhua.zhao@163. com (J. Zhao). https://doi.org/10.1016/j.carbon.2019.07.088 0008-6223/© 2019 Elsevier Ltd. All rights reserved.

structures such as triangle [20], trapezoid [21], U-shape [22] and localized pores structures [11]. Triangular graphene was reported that the maximum TR factor h (h ¼ (khi - klo)/klo, where khi and klo are the high and low TC values, respectively) is around 60% at room temperature [20]. TR effect can also be manipulated by asymmetric doping [23], defects [11,21], chemical functionalization [24], boundary thermal contacts [25] and covalent bonds junctions [26]. For instance, a triangular shaped boron-doped graphene has a TR factor of 30% at 300 K [23]. For a defective monolayer graphene with nanopores on one side, its TR factor was found to be 26% [21]. Asymmetrically functionalized by silicon, the TR factor of graphene can be tuned in the range of 0e30% by changing Si/C ratio at room temperature [24]. The above-mentioned TR manipulation techniques involve the modification of graphene structure, which more or less deteriorate the intrinsic properties of graphene. To address this issue, herein we present a new manipulation method using substrate. Due to its low-dimensional structural instability, graphene is usually supported by a substrate. The presence of the substrate unavoidably affects the properties of graphene and other twodimensional materials as evidenced by the available researches [2,3,16e18,27,28]. For instance, TC of graphene was reduced by

82

N. Wei et al. / Carbon 154 (2019) 81e89

~80% and ~44% when it is placed on a SiO2 substrate [3] and Cu substrate [18], respectively. Xu et al. reported that TR effect can be realized in a half encased graphene model with a TR factor of 40% [29]. Manipulating TR effect of graphene by substrate can keep graphene intact in structure and make it reusable. Therefore, it is of great significance to explore the substrate-controlled TR. To this end, in the present work, we investigate the effect of stiffnessinhomogeneous substrates on TR effect of graphene. Inhomogeneous stiffness also induced mass directed transport toward the stiffer part of a substrate, known as durotaxis, has been found in living cells [30,31] and inorganic graphene [32]. Chang et al. [32] reported a mass directed transport phenomenon using a stiffness gradient substrate. In their theoretical model where graphene is placed on a substrate with controllable stiffness, they found that the mass prefers to move from softer region toward stiffer region of a substrate. Inspired by their work, herein an inhomogeneous substrate is used to explore the substrate stiffness effects on TR of graphene. The materials with inhomogeneous stiffness widely exist, e.g., functionalized graded materials [33] or composite materials with different elastic moduli [34], which could also be generated through aligned nanotubes, nanopillars or nanowires [35e38], etc. Therefore, it is necessary and important to quantitatively understand the effect of stiffness-inhomogeneous substrates on TR effect of graphene. In the present study, non-equilibrium molecular dynamics (NEMD) simulations are performed to investigate TR effect of graphene supported on substrates with homogeneous and inhomogeneous stiffness. We find that TR effect occurs when a graphene sheet is supported on an inhomogeneous stiffness substrate, and the TR factor can be further manipulated in a large range by the external out-of-plane strain. These findings are explained well by analyzing the vibrational density of states (VDOS) of carbon atoms in graphene. 2. Modeling and method A substrate with adjustable stiffness is shown in Fig. 1a. The atoms in the substrate are arranged in the same as the graphene lattice arrangement, and the interactions among atoms in the substrate are ignored to exclude the substrate in the heat transport because graphene is the best thermal conductor and the tunable substrate usually consists of polymer materials, which has three orders lower magnitude of the thermal conductivity (<1 W/mK for polyethylene) [39]. Each atom in the substrate is linked with a spring, which provides only normal support (out-of-plane direction) to graphene. The in-pane motion of each substrate atom is constrained to prevent possible thermal flux in the substrate, so that we can only focus on the effect of the stiffness gradient on the thermal conductance of the upper graphene layer. These restrictions of the substrate model have little influence on the thermal conductivity of graphene in comparison with that of regulating interfacial binding strength [2,3,16,40e42]. The interactions between substrate and graphene atoms are using 12-6 Lennard-Jones (LJ) potential V(r) ¼ 4ε[(s/r)12-(s/r)6], and CarboneCarbon parameters are employed with ε ¼ 0.00284 eV and s ¼ 0.34 nm [43], unless otherwise stated. The optimized Tersoff potential by Lindsay and Broido [44], with small modifications, we call modified optimized Tersoff potential, and is used to describe the interactions among the carbon atoms in graphene. In comparison with the original Tersoff potential [45], Lindsay and Broido optimized two parameters, one for the equilibrium bond angle and one for the attractive interaction strength. According to this optimized Tersoff potential, the equilibrium bond length in graphene is 1.4388 Å, which is larger than the experimental value of 0.142 nm [46]. As the only length-related

Fig. 1. (a) Schematic illustration of tunable substrate model where the atoms in substrate are arranged based on graphene lattice and linked to springs. (b) Graphene layer is supported by the substrate. (c) Schematic plot of reverse non-equilibrium molecular dynamics simulation. (A colour version of this figure can be viewed online.)

parameters in the Tersoff potential are l1 in the repulsive function (fR ¼ A exp(-l1r)) and l2 in the attractive function (fA ¼ B exp(l2r)), we can obtain the correct bond length by multiplying these two parameters by a factor of 1.4388/1.42. That is, we change l1 from 3.4879 Å1 to 3.5333 Å1 and change l2 from 2.2119 Å1 to 2.2407 Å1. These modifications only change the length scale of the potential in a global way. Both graphene sheet and substrate have a width of 2 nm (x direction) and a length of 40 nm (y direction) (see Fig. 1c). The thickness of graphene is assumed to be 0.142 nm [10,15] for the TC calculation. All molecular dynamics (MD) simulations are performed using the large-scale atomic/molecular massively parallel simulator (LAMMPS) package [47]. Periodic boundary conditions are applied on both x and y directions. The model is initially optimized using the Polak-Ribiered version of conjugated gradient algorithm, fol-Hoover thermostat [48,49] (coupling constant lowed by a 2 ns Nose of 0.1 ps) to reach the equilibrium state at 300 K (with a time step of 0.5 fs). Then the system is switched to the microcanonical ensemble (NVE) to conserve the energy. The TC is calculated by a reverse NEMD method [50], as illustrated in Fig. 1c, where the system is divided into 50 slabs evenly along the heat transport direction (y direction). The 1st and 26th slabs are used as the heat sink and the heat source, respectively. To generate a temperature gradient, the heat flux is imposed on the system by continuously exchanging the kinetic energies between the hottest atom in the heat sink slab and the coldest atom in the heat source slab. After 2 ns, the system reaches non-equilibrium steady state and the temperature of each slab is collected and averaged in the following 3 ns. The heat flux J can be obtained by calculating the exchanging amount of the kinetic energy between two selected slabs according to the following equation

N. Wei et al. / Carbon 154 (2019) 81e89

  mv2h  mv2c

P J¼

1 Nswap 2

tswap

83

(1)

where Nswap and tswap are the entire number and time of swaps, and m is the mass of atom, while vh and vc represent the velocities of the hottest and coldest atoms, respectively. The TC k is calculated by using the Fourier law as

k¼

J 2AvT=vL

(2)

where A is the cross sectional area of heat transfer (A is obtained by multiplying the width and the thickness of the model), and vT/vL denotes the temperature gradient along the heat flow direction. The factor 2 represents the fact that the heat flux propagates in two directions away from the heat source.

3. Results and discussion 3.1. Effects of homogenous substrate stiffness First, we study the substrate stiffness effects on TC of graphene where the graphene sheet is supported on a homogeneous substrate. The schematic model is illustrated in Fig. 2a. As shown in Fig. 2b, TC of graphene decreases with increasing substrate stiffness. For example, when the substrate stiffness l increases from 1 nN/nm to 400 nN/nm, TC of graphene decreases from 571 W/mK to 517 W/mK. The reduction trend can be fitted well by an exponential function. It suggests that the maximum reduction of TC of graphene with infinite large substrate stiffness (which means that the atoms of the substrate are completely fixed) is ~12.8% (498 W/ mK) and there is a reduction of 15.9% when compared with the suspended graphene (592 W/mK). Our MD results of the thermal conductivity for the suspended graphene sheet are in good agreement with those of previous studies where the same or similar size of graphene was simulated [41,42]. The TC reduction is mainly attributed to the restriction of the out-of-plane movement of graphene and thus the out-of-plane phonon modes are suppressed [29]. The restriction effect increases with increasing substrate stiffness. For a given substrate stiffness of 320 nN/nm, the effect of ε (the LJ parameter between the carbon atom of graphene and the atom of the substrate) on the TC of graphene is studied, where the TC can be taken as a function of ε, as shown in Fig. 2c. The TC of graphene decreases with increasing ε (the LJ parameter between the carbon atom of graphene and the atom of the substrate) because the strengthening interactions for a large ε between graphene and the substrate further restrict the out-ofplane motion of graphene. For example, the TC drops by 15.1% when ε increases from 0.1 ε0 to 50.0 ε0 (ε0 ¼ 0.00284 eV). Larger ε indicates stronger interaction between graphene and the substrate, which hinder the out-of-plane motion of graphene and thus reduce its TC. When ε is set as an infinite large value, it implies that carbon atoms in graphene are fully bonded to the atoms of the substrate. In other words, the motion of carbon atoms in graphene is always synchronous with corresponding substrate atoms, which can be further understood as the carbon atoms in graphene being linked to springs directly, and the TC of graphene in this case is 425 W/mK. Then we extend the model to a more extreme condition with infinity of both l and ε. Under this condition, the out-of-plane motion of atoms in substrate is totally restrained and graphene can only move in the in-plane direction. The highly restrained graphene has a much reduced TC as 319 W/mK.

Fig. 2. (a) A graphene sheet is placed on a substrate with homogeneous stiffness. (b) TC of graphene sheet as a function of substrate stiffness. The data point can be fitted by an exponential function: k ¼ 498 þ 72.27e0.00335l. (c) TC k as a function of ε between graphene and substrate with a constant spring stiffness coefficient l ¼ 320 nN/nm, ε0 ¼ 0.00284 eV. (A colour version of this figure can be viewed online.)

3.2. Effects of inhomogenous substrate stiffness Since the substrate has a significant influence on TC of supported graphene, we designed a substrate with a stiffness gradient to explore its influence on TR effect of graphene. As shown in Fig. 3a1, a2, two types of stiffness gradient are designed, in which the high (low) stiffness coefficient l is set in the middle while the low (high) one is set in both ends of the substrate. The stiffness gradient g is defined as: g ¼ (lhard - lsoft)/(0.5L), where lhard and lsoft are the highest and lowest stiffness coefficients, respectively. L is the length of graphene (or substrate) along the y direction (heat transport direction). The heat flux flows from the hard (soft) region of substrate to the soft (hard) region is defined as: Jhs (Jsh), and the corresponding TC is denoted as khs (ksh). In view of the fact that TC changes more remarkably with l in the small stiffness region than that in the strong stiffness region, we change the stiffness gradient

84

N. Wei et al. / Carbon 154 (2019) 81e89

Fig. 3. (a1, a2) Graphene sheet adheres to substrate with inhomogeneous stiffness. The upper (lower) figure shows the substrate has the hardest (softest) stiffness part in the middle and softest (hardest) stiffness parts at both ends. ymid is the midpoint along the length of the graphene (y direction). (a3) Both k and h versus stiffness gradient g. The dot at g ¼ 0 represents the TC of graphene with homogeneous substrate (lsoft ¼ lhard ¼ 1.0 nN/nm). (a4) k as a function of ε between graphene and substrate with g ¼ 16 nN/nm2 (b1, b2) Within the schematic of sandwich model, the graphene sheet is stacked in the middle of two substrates with the same stiffness gradient. (b3) Both k and h versus stiffness gradient g. (b4) Both k and h versus △d/d0. d0 is the initial equilibrium distance between two substrates, △d is the variation distance from the equilibrium distance d0, and negative (positive) values represent compression (tension) in out-of-plane direction. (A colour version of this figure can be viewed online.)

g by keeping the minimum stiffness lsoft fixed as 1 nN/nm and changing the value of lhard. It is clearly seen from Fig. 3a3 that both khs and ksh decrease with increasing g, while ksh is higher than khs systematically. It means that heat flux prefers to flow from the soft region to the hard region, which also implies that heat flux prefers to flow from the high TC region to the low TC region. This directiondependent heat transport phenomenon, also called TR, is in good agreement with the findings reported in the previous studies [20]. TR factor h (h ¼ (ksh - khs)/khs) increases initially with g until it reaches a plateau (h ~10%) when the stiffness gradient g is greater than 8 nN/nm2. This plateau is attributed to the slight change of the TC at high l. For the case of g ¼ 20 nN/nm2, the maximum l is 400 nN/nm. After that, TC changes slowly with l. For instance, the TC of graphene decreases by only 3.7% (from 517 W/mK to 498 W/mK) when the stiffness coefficient l increases from 400 nN/nm to infinity. By contrast, TC of graphene decreases by 9.5% when l varies from 1 nN/nm to 400 nN/nm (see Fig. 2b). From Fig. 3a3, it is readily seen that the maximum TR factor h is ~11%, indicating that inhomogeneous support is not an efficient method in regulation of TR. In order to enhance the substrate effect for higher TR factor, we further manipulate the interfacial interaction strength between graphene and the substrate (see Fig. 3a4). Here, the substrate stiffness gradient is set as 16 nN/nm2 (corresponding to lsoft ¼ 1.0 nN/nm and lhard ¼ 320 nN/nm, respectively). Both khs and ksh decrease sharply with increasing ε, and TR factor h increases linearly to 22.9% when ε increases from 0.1 ε0 to 50.0 ε0. Obviously,

it is more effective to regulate TR factor by tuning the strength of the interactions between graphene and substrate, rather than the substrate stiffness. 3.3. A sandwich model To enhance the TR factor more effectively, a sandwich model (bisubstrates model) is proposed in Fig. 3b1, b2, where the graphene sheet is stacked in the middle of two inhomogeneous substrates. The same stiffness gradient g is set as the single substrate model (Fig. 3a1, a2). Compared to the single substrate model, the sandwich model shows that h monotonously increases initially and then reaches a plateau (h ~14%) when g is greater than 12 nN/nm2. The maximum TR factor of the sandwich model (14.4%) is slightly higher than that of the single substrate model (~11%). Based on the linear relationship between h and ε (see Fig. 3a4), the TR factor h is assumed as 12.7% at ε ¼ 2.0 ε0 (g ¼ 16 nN/nm2), which is still smaller than that of sandwich model (14.4%) at ε ¼ ε0 (g ¼ 16 nN/ nm2). Therefore, the sandwich model cannot be simply treated as a single substrate model with double interfacial interaction strength. It is evident that the sandwich model structure also makes contributions to restrict the out-of-plane vibration of middle graphene and thus leads to a higher TR effect. Next we use strain engineering to further tune the TC and TR of substrate-supported graphene [51,52]. In this simulation, tensile and compressive strains are applied on the sandwich model along the out-of-plane direction, as shown in Fig. 3b4 (inset). The out-of-

N. Wei et al. / Carbon 154 (2019) 81e89

Fig. 4. VDOS versus △d. (a) VDOS in in-plane direction. (b) VDOS in out-of-plane direction. (A colour version of this figure can be viewed online.)

plane strain can be experimentally conducted by compressing the sandwich model using water as the pressure transmission medium [53]. The variation of interlayer distance between two substrates △d is used to demonstrate the strain effect quantitatively. The initial distance d0 is 0.67 nm based on the assumption of the equilibrium interlayer distance between graphene and substrate as 0.335 nm. Based on the above finding in Fig. 3b3, i.e., h converges to a constant for g > 12 nN/nm2, a stiffness gradient parameter of 16 nN/nm2 is fixed here to study the strain effect (corresponding to lsoft ¼ 1.0 nN/nm and lhard ¼ 320 nN/nm). The process of the strain engineering can be divided into three distinct stages based on the status of graphene, as illustrated in Fig. 3b4. In the first stage (0.35 < △d/d0 < 0 or 0.23 nm < △d < 0), graphene is under compression along the out-of-plane direction. TC of graphene decreases sharply under squeezing process (reducing interlayer distance between two substrates d) at both heat flowing directions (from hard to soft region or

Fig. 5. (a) Schematic of ideal spring restricted model, all carbon atoms in graphene are linked to springs directly. (b) k as a function of spring stiffness coefficient l. The corresponding exponential fitting function: k ¼ 358.89 þ 221.76e0.00378l. (A colour version of this figure can be viewed online.)

85

Fig. 6. VDOS versus spring stiffness coefficient l. (a) VDOS in in-plane direction. (b) VDOS in out-of-plane direction. The values in figures are spring stiffness coefficient l (nN/nm). (A colour version of this figure can be viewed online.)

from soft to hard region). The main reason is that the compression suppresses the phonon modes in the out-of-plane direction and thus reduces TC of graphene. When subjected to out-of-plane compression, initially h increases sharply from 14.4% (strain-free state at △d ¼ 0) to 22.7% with △d/d0 ¼ 0.06 (△d ¼ 0.04 nm). Thereafter, h decreases quickly with decreasing △d. At the largest compression state considered here, i.e., △d/d0 ¼ 0.35 (△d ¼ 0.23 nm), h is reduced to 3.5%. In the second stage (0 <△d/d0 < 0.3 or 0 <△d < 0.2 nm), graphene is under out-of-plane tension. At the small tensile strain (0 <△d/d0 < 0.12 or 0 <△d < 0.08 nm), the graphene sheet is subjected to a small attractive force by van der Waals interactions from the substrate, which leads to a small enhancement in the TC due to the less restriction of the substrate induced by the tensile strain. After that, the TC of graphene at both heat flowing

Fig. 7. (a, b) Schematic of ideal spring restricted model. (c) Both TC k and TR factor h versus stiffness gradient g. (A colour version of this figure can be viewed online.)

86

N. Wei et al. / Carbon 154 (2019) 81e89

directions (khs and ksh) reach a plateau. On the other hand, h drops monotonously in this region. In the third stage with larger interlayer distance (△d/d0 > 1.0 or △d > 0.67 nm), the graphene sheet is adhered to only one substrate while the other one moves to a far distance and exerts negligible effect on the graphene. TC and TR factor are thus insensitive to △d. At this stage, the sandwich model is equivalent to the single substrate model. The strain-dependent TC can be explained by VDOS of phonons of carbon atoms in graphene. The VDOS are calculated by taking the Fourier transform of the atomic velocity auto-correlation functions as:

1 PðuÞ ¼ pffiffiffiffiffiffi 2p

∞ ð

* iut

e 0

N X

+ vj ðtÞvj ð0Þ dt

(3)

j¼1

where P(u) is the VDOS at frequency u and vj(t)vj(0) is the velocity

auto-correlation functions. Fig. 4 shows that the in-plane vibration has little influence under out-of-plane compression/tension, while the out-of-plane vibration is significantly suppressed. As a result, TC of graphene is reduced significantly, because the majority of TC of graphene is contributed by the phonon along the out-of-plane direction [54]. 3.4. An ideal spring restricted model To further explore the limitation of restriction effect on the TC of graphene, an ideal spring restricted model is proposed (see Fig. 5a). In this ideal model, no substrate is present while the carbon atoms in graphene are linked to springs directly. The out-of-plane motion is restrained by the spring for the large regulation range. The present ideal model is helpful to focus only on the underlying physics of the thermal rectification of graphene. Although a real substrate with a stiffness gradient is available, it is difficult to extract the

Fig. 8. Thermal energy transport function of the out-of-plane (a) and in-plane (b) excitations from the soft to hard regions, and vice versa. (c) The average transport distance of thermal energy d1 of the out-of-plane excitations, and d2 of the in-plane excitations. (A colour version of this figure can be viewed online.)

N. Wei et al. / Carbon 154 (2019) 81e89

effects from other factors such as mass density. The relationship between the spring stiffness coefficient l and k of graphene is plotted in Fig. 5b. The exponential fitting results show that the TC drops by 38.2% (from 581 W/mK to 359 W/mK) with increasing l from 1.0 nN/nm to infinity. The minimum value of TC (359 W/mK) obtained by exponential fitting is in a reasonable agreement with the TC of graphene (319 W/mK with l ¼ infinity, i.e., all atoms of graphene are completely fixed in the out-of-plane direction). The VDOS of graphene supported by different stiffness spring are shown in Fig. 6. Similar to the sandwich model, the VDOS are affected explicitly along the out-of-plane direction by comparison with the VDOS along the in-plane direction. The VDOS along the out-of-plane direction are significantly suppressed and thus lead to the reduction of TC.

87

We further explore the stiffness gradient effect on TR of graphene. Both TC and TR factors with various g are shown in Fig. 7. It is found that the value of khs (ksh) sharply decreases by 63.6% (49.5%) when g increases from 2 nN/nm2 to 10 nN/nm2. In particular, the reduction is more significant than that of uniform l, which has the reduction factor of 38.2%. This phenomenon is mainly attributed to the stiffness gradient which induces further phonon scattering in thermal transport. Similar to the mass gradient, the mismatch of phonon vibration spectra enhances the reduction of TC [55,56]. Fig. 7 shows that h reaches its maximum value (48.3%) at g ¼ 12 nN/nm2, and then it drops with increasing g. When g is over 12 nN/nm2, both ksh and khs show the tendency of convergence and the difference between them gradually narrows with g. Due to the fact that the thermal conductivity of graphene is mainly

Fig. 9. Thermal energy transport function, dE, as a function of heat transport direction (a1-a3) from soft to hard region of a substrate and (b1-b3) from hard to soft region of a substrate, respectively, under excitation strength cS ¼ 5 km/s at simulation time 1.0 ps. The leading wave profile, which leaves the fluctuating wave packets far behind, represents a shockwave induced by pulse. (A colour version of this figure can be viewed online.)

88

N. Wei et al. / Carbon 154 (2019) 81e89

contributed by out-of-plane phonon (over 70%) [57], the effect of the substrate on the thermal conductivity becomes more obviously with increasing length for the involvement of more long wavelength phonons. Thus the thermal rectification will increase with increasing sample length until the thermal conductivity of the supported graphene is converged. To confirm the microscopic origin of TR effect from the out-ofplane phonons, we study the thermal energy transport process of the out-of-plane and in-plane phonons through the excitation method [27,28]. The thermal excitation method provides a kinetic wave packet picture of the excited phonons so that we can fully understand the thermal conduction behavior. This method is widely used to study thermal conduction in low-dimensional materials such as black phosphorene [27,28], graphene [58], carbon nanotubes [59,60], interfacial water [61,62], anharmonic lattice [63,64] and gases [65], etc. The thermal energy transport function dE is defined as,

dE(x, y, t) ¼ E(x, y, t) - E0(x, y, t)

(4)

where E represents the energy of each carbon atoms in the nano sheet with the exerted heat pulse, and E0 represents the energy of the carbon atoms in the same nano sheet with the exact same initial condition other than the heat pulse. The snapshots of dE provide the propagating profiles of the heat pulse. Fig. 8a illustrates the thermal energy transport behavior of the out-of-plane excitation from the soft to the hard regions, and vice versa. It shows that more energy is lagging behind if the excitation is from the hard to the soft regions. Similarly, Fig. 8b illustrates the thermal energy transport behavior of the in-plane excitation from the soft to the hard regions, and vice versa. The two wave packets exhibit similar profile in the in-plane excitations (see Fig. 9). To numerically quantify the thermal energy transport efficiency of the out-of-plane and in-plane phonons, we consider the average transport distance of thermal energy d1 in the out-of-plane excitations, and d2 in the in-plane excitations as [27,28],

d¼

∬ ðy  y0 ÞdEdxdy ∬ dEdxdy

(5)

where x0 and y0 represents the initial excitation position at the sheet edges, and dE represents the respective excited thermal energy. From Fig. 8c, it is shown that the average thermal energy transport speed of the out-of-plane phonons from the soft to the hard region is faster than in the opposite direction, where d1 (softhard) ¼ 6 nm/ps and d1 (hard-soft) ¼ 2 nm/ps. The out-of-plane phonons exhibit 3-times higher transport speed along the thermal rectification preferred direction. On the other hand, identical thermal energy transport speed is observed for the in-plane phonons, where d1 (soft-hard) ¼ d1 (hard-soft) ¼ 12 nm/ps. The inplane phonons take no preference in the thermal transport process, further explaining that thermal rectification effect originates from the out-of-plane phonons rather than the in-plane phonons. 4. Conclusion We have proposed an undamaged thermal rectification method using rational substrate stiffness design in supported graphene. We find that TR phenomenon occurs when a graphene sheet is supported on an inhomogeneous stiffness substrate and heat flux prefers to flow toward the stiffer part of a substrate. The maximum thermal rectification factor can reach ~48% by using the ideal substrate model prediction. Moreover, a sandwich model can manipulate the thermal rectification factor in a wide range of 3.5%e 22.7%. The orientation bias of heat transfer is attributed to

suppression of out-of-plane phonon transport and it is well explained by diversity of energy transport in graphene using excitation simulation. Our results should be of great help towards understanding the modulation of thermal rectification in supported material and shed light on the design of a new thermal logic gate and thermal management devices. Financial interest The authors declare no competing financial interests. Acknowledgements We gratefully acknowledge support by the National Natural Science Foundation of China (Grant Nos. 11572140, 11302084, 11972102 11502217 and 11405245), HPC of NWAFU, the Youth Training Project of Northwest A&F University (No. Z109021600) the Natural Science Foundation of Jiangsu Province (Grant No. BK20180031) and the Fundamental Research Funds for the Central Universities (No. 2452015054). Natural Science Foundation of ShaanXi (Grant. 2018JM1008), China Postdoctoral Science Foundation (No. 2015M570854 and 2016T90949), the 111 project (Grant No. B18027). References [1] S. Ghosh, I. Calizo, D. Teweldebrhan, E.P. Pokatilov, D.L. Nika, A.A. Balandin, et al., Extremely high thermal conductivity of graphene: prospects for thermal management applications in nanoelectronic circuits, Appl. Phys. Lett. 92 (15) (2008) 151911. [2] W.W. Cai, A.L. Moore, Y.W. Zhu, X.S. Li, S.S. Chen, L. Shi, et al., Thermal transport in suspended and supported monolayer graphene grown by chemical vapor deposition, Nano Lett. 10 (5) (2010) 1645e1651. [3] J.H. Seol, I. Jo, A.L. Moore, L. Lindsay, Z.H. Aitken, M.T. Pettes, et al., Twodimensional phonon transport in supported graphene, Science 328 (5975) (2010) 213e216. [4] D.L. Nika, E.P. Pokatilov, A.S. Askerov, A.A. Balandin, Phonon thermal conduction in graphene: role of umklapp and edge roughness scattering, Phys. Rev. B 79 (15) (2009) 155413. [5] L. Lindsay, D.A. Broido, N. Mingo, Flexural phonons and thermal transport in graphene, Phys. Rev. B 82 (11) (2010) 115427. [6] W.J. Evans, L. Hu, P. Keblinski, Thermal conductivity of graphene ribbons from equilibrium molecular dynamics: effect of ribbon width, edge roughness, and hydrogen termination, Appl. Phys. Lett. 96 (20) (2010) 203112. [7] H.J. Zhang, G. Lee, A.F. Fonseca, T.L. Borders, K. Cho, Isotope effect on the thermal conductivity of graphene, J. Nanomater. 2010 (1) (2010) 537657. [8] J.W. Jiang, J.H. Lan, J.S. Wang, B.W. Li, Isotopic effects on the thermal conductivity of graphene nanoribbons: localization mechanism, J. Appl. Phys. 107 (5) (2010), 054314. [9] F. Hao, D.N. Fang, Z.P. Xu, Mechanical and thermal transport properties of graphene with defects, Appl. Phys. Lett. 99 (4) (2011), 041901. [10] T.Y. Ng, J.J. Yeo, Z.S. Liu, A molecular dynamics study of the thermal conductivity of graphene nanoribbons containing dispersed StoneeThrowereWales defects, Carbon 50 (13) (2012) 4887e4893. [11] Y. Wang, S. Chen, X. Ruan, Tunable thermal rectification in graphene nanoribbons through defect engineering: a molecular dynamics study, Appl. Phys. Lett. 100 (16) (2012) 163101. [12] M. Yarifard, J. Davoodi, H. Rafii-Tabar, In-plane thermal conductivity of graphene nanomesh: a molecular dynamics study, Comput. Mater. Sci. 111 (2016) 247e251. [13] L. Hu, D. Maroudas, Thermal transport properties of graphene nanomeshes, J. Appl. Phys. 116 (18) (2014) 184304. [14] N. Wei, Y. Chen, K. Cai, J.H. Zhao, H.Q. Wang, J.C. Zheng, Thermal conductivity of graphene kirigami: ultralow and strain robustness, Carbon 104 (3) (2016) 203e213. [15] N. Wei, L.Q. Xu, H.Q. Wang, J.C. Zheng, Strain engineering of thermal conductivity in graphene sheets and nanoribbons: a demonstration of magic flexibility, Nanotechnology 22 (10) (2011) 105705. [16] J. Chen, G. Zhang, B. Li, Substrate coupling suppresses size dependence of thermal conductivity in supported graphene, Nanoscale 5 (2) (2013) 532e536. [17] Z.X. Guo, J.W. Ding, X.G. Gong, Substrate effects on the thermal conductivity of epitaxial graphene nanoribbons, Phys. Rev. B 85 (23) (2012) 235429. [18] L. Chen, S. Kumar, Thermal transport in graphene supported on copper, J. Appl. Phys. 112 (4) (2012), 043502. [19] B.W. Li, L. Wang, G. Casati, Thermal diode: rectification of heat flux, Phys. Rev. Lett. 93 (18) (2004) 184301.

N. Wei et al. / Carbon 154 (2019) 81e89 [20] J.N. Hu, X.L. Ruan, Y.P. Chen, Thermal conductivity and thermal rectification in graphene nanoribbons: a molecular dynamics study, Nano Lett. 9 (7) (2009) 2730e2735. [21] H. Wang, S. Hu, K. Takahashi, X. Zhang, H. Takamatsu, J. Chen, Experimental study of thermal rectification in suspended monolayer graphene, Nat. Commun. 8 (2017) 15843. [22] J.G. Cheh, H. Zhao, Thermal rectification in asymmetric U-shaped graphene flakes, J. Stat. Mech. 6 (2012) P06011, 2012. [23] Q. Liang, Y. Wei, Molecular dynamics study on the thermal conductivity and thermal rectification in graphene with geometric variations of doped boron, Physica B 437 (2014) 36e40. [24] K.P. Yuan, M.M. Sun, Z.L. Wang, D.W. Tang, Tunable thermal rectification in silicon-functionalized graphene nanoribbons by molecular dynamics simulation, Int. J. Therm. Sci. 98 (2015) 24e31. [25] X. Yang, J. Xu, S. Wu, D. Yu, B. Cao, The effect of structural asymmetry on thermal rectification in nanostructures, J. Phys. Condens. Matter 30 (43) (2018) 435305. [26] X. Yang, D. Yu, B. Cao, A.C. To, Ultrahigh thermal rectification in pillared graphene structure with carbon nanotube-graphene intramolecular junctions, ACS Appl. Mater. Interfaces 9 (1) (2016) 29e35. [27] J.G. Chen, S.D. Chen, Y. Gao, Anisotropy enhancement of thermal energy transport in supported black phosphorene, J. Phys. Chem. Lett. 7 (13) (2016) 2518e2523. [28] J.G. Chen, S.D. Chen, Y. Gao, Supersonic thermal excitation-induced shock wave in black phosphorene, Phys. Rev. B 95 (13) (2017) 134301. [29] W. Xu, G. Zhang, B.W. Li, Interfacial thermal resistance and thermal rectification between suspended and encased single layer graphene, J. Appl. Phys. 116 (13) (2014) 134303. [30] C.-M. Lo, H.-B. Wang, M. Dembo, Y.-l. Wang, Cell movement is guided by the rigidity of the substrate, Biophys. J. 79 (2000) 144e152. [31] D.E. Discher, P. Janmey, Y.L. Wang, Tissue cells feel and respond to the stiffness of their substrate, Science 310 (5751) (2005) 1139e1143. [32] T.C. Chang, H.W. Zhang, Z. Guo, X.M. Guo, H.J. Gao, Nanoscale directional motion towards regions of stiffness, Phys. Rev. Lett. 114 (1) (2015), 015504. [33] S.C. Hernandez, C.J.C. Bennett, C.E. Junkermeier, S.D. Tsoi, F.J. Bezares, R. Stine, et al., Chemical gradients on graphene to drive droplet motion, ACS Nano 7 (6) (2013) 4746e4755. [34] X.F. Wei, T. Zhang, T.F. Luo, Thermal energy transport across hardesoft interfaces, ACS Energy Lett 2 (10) (2017) 2283e2292. [35] Q. Cao, S.J. Han, G.S. Tulevski, Y. Zhu, D.D. Lu, W. Haensch, Arrays of singlewalled carbon nanotubes with full surface coverage for high-performance electronics, Nat. Nanotechnol. 8 (3) (2013) 180e186. [36] N. Liebing, S. Serrano-Guisan, K. Rott, G. Reiss, J. Langer, B. Ocker, et al., Tunneling magnetothermopower in magnetic tunnel junction nanopillars, Phys. Rev. Lett. 107 (17) (2011) 177201. [37] M. Dietzel, S.M. Troian, Formation of nanopillar arrays in ultrathin viscous films: the critical role of thermocapillary stresses, Phys. Rev. Lett. 103 (7) (2009), 074501. [38] J. Ji, J. Zhao, W. Guo , Novel nonlinear coarse-grained potentials of carbon nanotubes, J. Mech. Phys. Solids 128 (2019) 79e104 . [39] J.H. Zhao, J.W. Jiang, N. Wei, Y.C. Zhang, T. Rabczuk, Thermal conductivity dependence on chain length in amorphous polymers, J. Appl. Phys. 113 (18) (2013) 184304. [40] 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) 396e400. [41] Z. Wei, Z. Ni, K. Bi, M. Chen, Y. Chen, In-plane lattice thermal conductivities of multilayer graphene films, Carbon 49 (2011) 2653e2658. [42] T. Li, Z. Tang, Z. Huang, J. Yu, Interfacial thermal resistance of 2D and 1D

[43] [44]

[45] [46] [47] [48] [49] [50]

[51]

[52]

[53]

[54] [55]

[56]

[57]

[58] [59] [60] [61]

[62] [63] [64]

[65]

89

carbon/hexagonal boron nitride van der Waals heterostructures, Carbon 105 (2016) 566e571. S.J. Stuart, A.B. Tutein, J.A. Harrison, A reactive potential for hydrocarbons with intermolecular interactions, J. Chem. Phys. 112 (14) (2000) 6472e6486. L. Lindsay, D.A. Broido, Optimized tersoff and brenner empirical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene, Phys. Rev. B 81 (20) (2010) 262e265. J. Tersoff, Modeling solid-state chemistry: interatomic potentials for multicomponent systems, Phys. Rev. B 39 (8) (1989) 5566e5568. Y. Baskin, L. Meyer, Lattice constants of graphite at low temperatures, Phys. Rev. 100 (2) (1955), 544-544. S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117 (1) (1995) 1e19. , A unified formulation of the constant temperature molecular dyS. Nose namics methods, J. Chem. Phys. 81 (1) (1984) 511e519. W.G. Hoover, Canonical dynamics: equilibrium phase-space distributions, Phys. Rev. A 31 (3) (1985) 1695e1697. F. Müller-Plathe, A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity, J. Chem. Phys. 106 (14) (1997) 6082e6085. Y. Chen, Y.Y. Zhang, K. Cai, J.W. Jiang, J.C. Zheng, J.H. Zhao, et al., Interfacial thermal conductance in graphene/black phosphorus heterogeneous structures, Carbon 117 (2017) 399e410. J.J. Jiang, W.C. Fu, J.G. Chen, H. Zhao, Anharmonicity induced thermal modulation in stressed graphene, Sci. China Phys. Mech. Astron. 60 (7) (2017), 070512. L.G.P. Martins, M.J.S. Matos, A.R. Paschoal, P.T.C. Freire, N.F. Andrade, A.L. Aguiar, et al., Raman evidence for pressure-induced formation of diamondene, Nat. Commun. 8 (1) (2018) 96. M. Hu, X.L. Zhang, D. Poulikakos, Anomalous thermal response of silicene to uniaxial stretching, Phys. Rev. B 87 (19) (2013) 195417. Nuo Yang, Nianbei Li, Lei Wang, B. Li, Thermal rectification and negative differential thermal resistance in lattices with mass gradient, Phys. Rev. B 76 (2) (2007) 20301. M. Alaghemandi, F. Leroy, E. Algaer, M.C. Bohm, F. Muller-Plathe, Thermal rectification in mass-graded nanotubes: a model approach in the framework of reverse non-equilibrium molecular dynamics simulations, Nanotechnology 21 (7) (2010), 075704. Z.Y. Fan, L.F.C. Pereira, P. Hirvonen, M.M. Ervasti, K.R. Elder, D. Donadio, et al., Thermal conductivity decomposition in two-dimensional materials: application to graphene, Phys. Rev. B 95 (14) (2017) 144309. J. Chen, W. Qi, M. Zhang, H. Zhao, Kinetic behavior of subsonic solitary wave in graphene nanoribbon, J. Stat. Mech. Theory Exp. 2015 (2015) P06007. L. Chen, S. Kumar, Thermal transport in double-wall carbon nanotubes using heat pulse, J. Appl. Phys. 110 (7) (2011), 074305. E. Gonzalez Noya, D. Srivastava, M. Menon, Heat-pulse rectification in carbon nanotube Y junctions, Phys. Rev. B 79 (11) (2009) 115432. J. Cheh, Y. Gao, C. Wang, H. Zhao, H. Fang, Ice or water: thermal properties of monolayer water adsorbed on a substrate, J. Stat. Mech. Theory Exp. 2013 (2013) P06009. J. Chen, C. Wang, N. Wei, R. Wan, Y. Gao, 3D flexible water channel: stretchability of nanoscale water bridge, Nanoscale 8 (10) (2016) 5676e5681. A. Dhar, Heat transport in low-dimensional systems, Adv. Phys. 57 (2008) 457e537. E. Arevalo, F.G. Mertens, Y. Gaididei, A.R. Bishop, Thermal diffusion of supersonic solitons in an anharmonic chain of atoms, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 67 (2003), 016610. J. Yang, Y. Zhang, J. Wang, H. Zhao, Energy transfer process in gas models of Lennard-Jones interactions, Phys. Rev. E 83 (2011), 052104.