Thermal conductivity of graphene nanoribbons under shear deformation: A molecular dynamics simulation

Tensile strain and compress strain can greatly affect the thermal conductivity of graphene nanoribbons (GNRs). However, the effect of GNRs under shear strain, which is also one of the main strain effect, has not been studied systematically yet. In this work, we employ reverse nonequilibrium molecular dynamics (RNEMD) to the systematical study of the thermal conductivity of GNRs (with model size of 4 nm × 15 nm) under the shear strain. Our studies show that the thermal conductivity of GNRs is not sensitive to the shear strain, and the thermal conductivity decreases only 12–16% before the pristine structure is broken. Furthermore, the phonon frequency and the change of the micro-structure of GNRs, such as band angel and bond length, are analyzed to explore the tendency of thermal conductivity. The results show that the main influence of shear strain is on the in-plane phonon density of states (PDOS), whose G band (higher frequency peaks) moved to the low frequency, thus the thermal conductivity is decreased. The unique thermal properties of GNRs under shear strains suggest their great potentials for graphene nanodevices and great potentials in the thermal managements and thermoelectric applications.

Graphene, a single sheet 2D material of carbon atoms, has recently attracted significant attention due to its extraordinary mechanical and electrical properties. The attractive properties of graphene are being explored for series applications including nanomechanical, nanoelectronics,etc. It has the especially excellent applicable potentials since it was experimentally accessible in 2004 [1][2][3][4][5][6] . Graphene is regarded as an ideal heat transfer material as it is endowed with extremely high in-plane thermal conductivity 7 . The thermal conductivity of graphene can be adapted according to different requirement of different devices, such as on one hand the better radiating properties for the devices when the thermal conductivity is increased and on the other hand the more effective thermalelectric properties of the devices when the thermal conductivity is decreased. However, the properties of graphene can be influenced by defects and tailored geometry shapes (size and asymmetry), strain etc. It should be pointed out that the control of graphene with defects and doping is irreversible and unrecoverable due to that its pristine structure will break, in contrast the properties of graphene under strain can be recycled. The mechanical properties of graphene nanoribbons (GNRs) under tension and shear have been explored by experiments and atomistic simulation methods [8][9][10][11][12][13][14] . The study of Jiang et al. observed that the Young's modulus is closely linked with the temperature in the certain region 15 . The Young's modulus has been shown depending on the concentration of defects to a certain extent, and the thermal conductivity is explored to be much more sensitive than the mechanical modulus to the defects 16 etc. Recently, a new structure of graphene shows the extreme reduction of the thermal conductivity by tailoring sizes in graphene nanoribbon kirigami, whose thermal conductivity is about two orders of magnitude lower than that of the corresponding GNR 23 . Furthermore, the thermal conductivity of graphene under uniaxial tensile has been verified to appear a decrease of 45% and 60% in zigzag and armchair direction 24 , respectively. Therefore, it is critical to clear the effect of strain on the thermal conductivity of such structural materials.
Strain effects on thermal conductivity of the low dimensional (two dimensional and one dimensional) materials have been studied for last at least two decades [25][26][27] . It's reported that the thermal conductivity of the low dimensional nanostructures as the silicon and diamond nanowires and thin films is decreased continuously from compressive strain to tensile strain 27 . The thermal conductivity of GNRs is also found to appear a tendency of decreasing under tensile strain by classical molecular dynamics (MD) simulations 24,28 . Hence, the strain effects of graphene can play a key role in the continuous tunability and applicability of its thermal conductivity property at nanoscale, and the dissipation of thermal conductivity is a obstacle for the applications of thermal management, such as in the complementary metal-oxide-semiconductor technology and chemical vapor deposition technique. Up to now, the thermal conductivity of graphene under shear deformation has not been investigated yet. From a practical point of view, good thermal managements of GNRs have significantly potential applications of future GNR-based thermal nanodevices, which can greatly improve performances of the nanosized devices due to heat dissipations. Meanwhile, graphene is a thin membrane structure, it is also important to understand the wrinkling behavior under shear deformation. The wrinkling behavior of aluminized Kapton membranes and graphene has been verified to show an essential capability to keep its one-atom layer nanostructures stable 29 . Furthermore, graphene will form corrugation under small shear strain, and failure with bonds break when shear deformation is large enough.
In this work, we have performed computer simulations based on reverse nonequilibrium molecular dynamics (RNEMD) 30 to investigate shear strain effects on thermal conductivity of graphene.

Methods
The open source Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) 31 is employed in our simulation. The adaptive intermolecular reactive empirical bond-order (AIREBO) potential function is employed to calculate interaction between carbon atoms in graphene, which includes both covalent bonding (stretching, bending, torsion) expressed in terms of bond orders, as well as van der Waals interactions between atoms at a larger distance: where the E ij REBO term has the same functional form as the hydrocarbon REBO potential 32 . The coefficients for E ij REBO are the same as in Brenner's potential. The E ij LJ term considers the long-range interactions (2 < r < cut off) using a form similar to the standard 12-6 Lennard Jones potential, and the E k jil TORSION term is an explicit 4-body potential that describes various dihedral angle preferences in configurations of carbon atoms.
The thermal of 300 K is maintained as the background temperature by a Nosé-Hoover thermostat 33,34 . The system reaches equilibrium state in 200 picosecond (ps), and then a 0.5 femtosecond (fs) is imposed as the time step with a heat flux. The shear strains are applied according to Fig. 1(a,b). The different boundary conditions may result in different thermal conductivity. In this work, periodic boundary condition is used in all three directions..
To compute the thermal conductivity properties of sheared graphene, the size of the model is initially established in 4 nm × 15 nm. The constant pressure and constant temperature (NPT) simulation is used until the energy of the system is fully minimized. The deformation-control method 13 is employed by applying the shear strain rate of 0.001 ps −1 to the bulk graphene in both zigzag graphene nanoribbons (ZGNR) and armchair graphene nanoribbons (AGNR) cases, separately. The constant volume and constant temperature (NVT) simulation is used in every 1000 time steps, which followed applied a strain increment. The time step, 0.1 fs, is applied by velocity-Verlet time stepping scheme. The shear modulus, G, can be gotten by the following calculation: where V 0 is the volume of graphene, E is the strain energy, and γ xy is the shear strain. The shear strain expresses as γ xy = Δ x/L, in which x and L are the transverse displacement and the initial length of graphene, respectively. The linear portion of the shear stress and strain curve is used, xy xy to verify the calculation of shear modulus. The τ xy is the shear stress at 0.005 strain. The result here shows that the linear portion in Eq. (3) is consistent with the result obtained based on the strain energy method described as Eq. (2). Here, we compute the thermal conductivity by using RNEMD simulation. The main idea of this method is that the model is imposed a heat flux J as Fig. 2 to form a temperature gradient. The advantage is that the convergence of the temperature gradient can be achieved quickly in a nonequilibrium steady system with a known quantity of heat flux. The heat flux J transfers energy continuously from the "hot" slab, located in the middle of the structure, Scientific RepoRts | 7:41398 | DOI: 10.1038/srep41398 to the "cold" slab, located at the beginning and end of the structure. The detail of the setup description is shown in Fig. 2. The heat flux J is from the heat bath to the system, and it is given by where k B is Boltzmann's constant, N j is the number of the particles in the ith slab, p j is the momentum of atom j, m j and v j is the mass and velocity of the jth carbon atom, respectively. The temperature profile is obtained by time-averaged with averaging in a time interval of 100 ps. The time-averaged temperature profile, as a function of atomic position along the nanoribbon's axis, is used to compute the temperature gradient ∂ ∂ T x . Then, the thermal conductivity of the system is calculated by using Fourier law: where A is the cross-sectional area perpendicular to the direction of heat transfer, defined by the width of GNR times the thickness of the GNR. The thickness of graphene is the carbon-carbon bond length, and the value is assumed to be 0.142 nm 28,35 here. Here, the ensemble average is substituted with time averaging. In this MD simulation, the time step is 0.5 fs and the velocities are recorded every two steps. The averaging is performed over 10 5 steps after the system reaches equilibrium for 200 ps.

Results and Discussions
The shear strain effect on the mechanical structure of graphene. Shear stress-strain curve. We study the mechanical properties of graphene along the zigzag and armchair direction under the shear strain as shown in Fig. 1. The strain at fracture along the zigzag and armchair direction is 30% and 27%, and the broken shear stress is 53.8 GPa and 53.1 GPa, respectively. The difference of anisotropic property under shear strain is less than that under uniaxial tensile. Under uniaxial tensile, the broken strain in zigzag direction of graphene is 150% 36 of the strain in armchair direction, and broken stress in zigzag direction of graphene is 167% 36 of the broken stress in armchair direction, respectively.
The change of the corrugation size under shear strain in macro view. The change of amplitude and wavelength of the wrinkle graphene under different shear strain are studied. As illustrated in Fig. 3(a), shear strains cause helical wrinkles in the graphene structure. The relation of wrinkle size and shear strain is studied next. Figure 3(c,d) show the amplitude of the wave is increased as the shear strain increases, but the wavelength is decreased as the shear strain increases, probably since the difference in the bending stiffness between zigzag and armchair direction is minor. Hence, the wrinkle structure in zigzag and armchair direction under shear strain are nearly identical. The micro structure change of graphene under shear strain. The micro structure change is commonly reported based on the bond-angle and bond-length see Fig. 1. Distinct from uniaxial tensile, all the bond-angles and bond-lengths are changed independently with the shear strain, due to the broken axisymmetry. The bond-angle and bond-length due to shear strain relation are shown in Fig. 4. The carbon-carbon bonds are elongated and angles are changed under shear stress. In order to analyze the results quantitatively, we calculated the carbon-carbon bond lengths and angles with every 0.01 strain step and the obtained values are averaged overall. We can clearly see that the bonds along the shear direction, bond type C for ZGNR and bond type A for AGNR (see Fig. 1),undergo more deformation.
Variations of bond length and bond angle of graphene are studied under shear strain in contrast with the results under uniaxial tensile. Figure 4 clearly shows the effect of shear stress to both bond-angle and bond-length are less than that of uniaxial tensile. The carbon-carbon bonds and their angles change for different values of shear stress (approximately 1.40 Å and 120° for bond length and angle, respectively without strains) is illustrated in Fig. 4. From shear stress 0.05 to 0.25 Gpa, the bond length of ZGNR and AGNR increase, in the zigzag direction, the bond type C grows faster than the bond type A, while the bond type B almost never changes; in the armchair direction, the bond type A grows faster than the bond type C. The bond type C remains nearly unaltered. The ZGNR and AGNR all have one type angle increases, but the other two types angle decrease, in the ZGNR, the angle three increase, the other angle all decrease, in the AGNR, the angle two increase, the other angle all decrease. The different changes of bond types are the reason of the different changes of angle types. These changes can also affect the thermal conductivity, and the atomic structure of graphene sheet under shear strain provides the physical insight for the mechanism of thermal transport.
Thermal conductivity of graphene under shear strain. Thermal conductivity. Now, we study the relative thermal conductivity of graphene under shear strain. Figure 5 shows the changing process of thermal conductivity of graphene with the increase of the shear strain. With increasing shear strain coverage from 0% to 20% or 25%, the thermal conductivity increases firstly and decreases finally in both zigzag and armchair cases. At 5% coverage of the shear strain, the thermal conductivity in zigzag direction reaches the peak value, then it decreases subsequently. Different from zigzag direction, the thermal conductivity in armchair direction increases from 0% to 15% of the shear strain, and achieves its peak value at 15% of the shear strain, then it drop down rapidly. The thermal conductivity barely decreases in both zigzag and armchair direction. Compared with thermal conductivity which decrease 45% and 60% 24 in zigzag and armchair direction under uniaxial tensile, respectively,  Fig. 1(a,b). Here, (a) is the shear bond-length in armchair and zigzag direction. (b) is the shear bond-angle in armchair and zigzag direction.
Scientific RepoRts | 7:41398 | DOI: 10.1038/srep41398 the thermal conductivity under shear strain decreases only 12% and 16% (see Fig. 5) in zigzag and armchair direction, respectively. Thus, we conclude that the thermal conductivity of graphene is not sensitive w.r.t. shear strain. Hence, graphene nanodevices should remain stable under shear strain, but should not be used under tensile strain.
Phonon density of states (PDOS). The strain effect on thermal transport properties is considered in this work by investigating the variation of the phonon spectra (ZGNR and AGNR cases) under different strains (see Fig. 5). The underlying mechanism of the thermal conductivity can be explained by the phonon spectra or the phonon density of states (PDOS) for the graphene under shear strain. The phonon spectrum function G(ω) is computed from the velocity autocorrelation function 37 (VACF) by using Fourier transform as follows: where v j (t) and v j (0) is the velocity of jth particle at time t and time 0, respectively, and ω denotes the vibrational wavenumber. The G band (higher frequency peaks) has a red shift toward the lower frequency region (see the arrow in Fig. 6). It can be seen that the G band is softened by the shear strain, and then the G band slows down the phonon group velocities. Additionally, the shear strain can produce the corrugation of graphene, which modifies the local vibrational characteristics and induces an acoustic mismatch. These appearances scatter phonons and reduce the phonons mean free path (MFP) l. Hence, the thermal conductivity is decreased according to the classical lattice thermal transport theory: where m is the phonon mode value occupied at a specific temperature; C m , v m , l are the specific heat, group velocity and MFP of the phonon, respectively. It's well known that the acoustic phonons are the main heat carriers in graphene near room temperature. Moreover, it was shown both theoretically and experimentally that transport properties of phonons are substantially different in graphene, which may lead to the unusual thermal conduction in graphene [38][39][40][41] . Additionally, the contribution of long-wavelength phonons is usually underestimated. We further analyze the distribution of PDOS under different shear strain. It can be seen from Fig. 6, there is no clear difference of out-plane PDOS under shear strain, but the G-peak of the in-plane PDOS under shear strain shows a clear red-shift to the low frequency, which is similar to the decrease of thermal conductivity of graphene under uniaxial tensile. The tensile of the bond-length and the change of the bond-angle under shear strain causes phonon softening, followed by a decrease in the velocity of the phonon group, resulting in the decrease of the thermal conductivity. Furthermore, the change of phonon red-shift under shear strain is far less the value of uniaxial tensile, thus the decrease of the thermal conductivity is minor.

Conclusion
In this paper, the thermal conductivity of graphene under shear strain is systematically studied by MD simulations. The results show that, in contrast to the dramatic decrease of thermal conductivity of graphene under uniaxial tensile, the thermal conductivity of graphene is not sensitive to the shear strain, and the thermal conductivity decreases only 12-16%. Studying the thermal conductivity of graphene under shear strain, we analyze the relationship of wrinkle structure of graphene sheet under shear strain. The wrinkle evolve when the shear strain is around 5-10%, but the thermal conductivity barely changes. There is evidence that the thermal conductivity of graphene does not converge even when its size reaches 9 μm 42 . When the size of simulation model is smaller, considering the size effect, the wrinkle deformation could be restricted and the shape change would mainly on the atomic structure. Therefore, the thermal conductivity of graphene is more sensitive with decreasing the size of graphene. When a bigger simulation model size is employed, in contrast, the strain will express more as the wrinkle form. We found that the wrinkle has a relative small effect on the thermal conductivity of graphene in our previous research. Thus, we think that the thermal conductivity of graphene is much less sensitive to the shear strain with the size increasing. Furthermore, according to the PDOS results, we find that there is no clear difference of out-plane PDOS under shear strain, but the in-plane PDOS shows that G-peak under shear strain moves to the low frequency and shows red-shift. Since the change of phonon red-shift under shear strain is far beyond the value of uniaxial tensile, the decrease of the thermal conductivity is minor. In summary, we come to a conclusion that the thermal conductivity of graphene is not sensitive to the shear strain. Therefore, graphene nanodevices should be kept under shear strain to get relative stable thermal conductivity, and they should be avoided to be under tensile strain.