Steric Interference in Bilayer Graphene with Point Dislocations

nanomaterials Article Steric Interference in Bilayer Graphene with Point Dislocations Francisco Arca 1, Juan Pedro Mendez 2, Michael Ortiz 2and Pilar Ariza 1,* 1Escuela Técnica Superior de Ingeniería, University of Seville, 41092 Seville, Spain 2 Division of Engineering and Applied Science, California Institute of Technology, Pasadena, CA 91125, USA *Correspondence: [email protected] Received: 21 June 2019; Accepted: 5 July 2019; Published: 14 July 2019   Abstract: We present evidence of strong steric interference in bilayer graphene containing offset point dislocations. Calculations are carried out with Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) using the Long-Range Carbon Bond-Order Potential (LCBOP) potential of Los et al.. We start by validating the potential in the harmonic response by comparing the predicted phonon dispersion curves to experimental data and other potentials. The requisite force constants are derived by linearization of the potential and are presented in full form. We then continue to validate the potential in applications involving the formation of dislocation dipoles and quadrupoles in monolayer configurations. Finally, we evaluate a number of dislocation quadrupole configurations in monolayer and bilayer graphene and document strong steric interactions due to out-of-plane displacements when the dislocations on the individual layers are sufficiently offset with respect to each other. Keywords: graphene; force constants model; distributed dislocations 1. Introduction Bilayer graphene was reported by Noselov et al. [ 1 ], who described graphene configurations containing one, two and three atomic layers. Interest in monolayer and bilayer graphene stems from their outstanding electronic and mechanical properties [ 2 – 4 ], including high thermal mobility, above 4000 W/mK [ 5 ], high electronic conductivity, above 15,000 cm 2 /Vs [ 6 ], low mass density, 0.77 mg/m 2 , and high breaking strength [ 7 ]. Owing to these exceptional properties, in conjunction with ever improving production techniques, graphene remains a material of interest for potential application in next-generation electronic devices [8,9]. The physical properties of graphene are strongly influenced by the presence of defects within the lattice. For instance, partial dislocations in graphene give rise to profound changes in transport properties [ 10 , 11 ] and interesting topological states [ 12 ]. For this reason, graphene defects, such as dislocations, vacancies, grain boundaries or doped graphene, have been studied using different interatomic potentials and several computational approaches, from ab initio methods to molecular dynamics. One of the most common graphene and carbon nanotube defect structures is the Stone-Wales [13] , which is composed of two pentagon-heptagon cell pairs resulting from the rotation of a single atomic bond through 90 ◦ . Meyer et al. [ 14 ] examined this kind of defect using transmission electron microscopy (TEM). Li et al. [ 15 ] calculated the formation energy and the activation barrier of the Stone-Wales defect using the density-functional theory. Xiao et al. [ 16 ] studied the deterioration of the mechanical properties of graphene due to the presence of this type of defects by using an atomistic based finite bond element model. It is found that the Stone-Wales defect modifies the chemical [ 17 ] and electronic [ 18 , 19 ] properties of pristine graphene lattices. 7-5 pair structures can also combine to form dislocation dipoles [ 20 ]. Lehtinen et al. [ 21 ] have applied high-energy electron irradiation to Nanomaterials 2019,9, 1012; doi:10.3390/nano9071012 www.mdpi.com/journal/nanomaterials Nanomaterials 2019,9, 1012 2 of 16 graphene lattices triggering vacancy-type defects, observing that for a certain number of vacancies the atoms locally reorganize into a dipole terminated by two edge dislocations. Warner et al. [ 22 ] have explained that dislocation dipoles can also be formed during the CVD growth, through the addition of surface adatoms or to Stone-Wales bond rotations. Jeong et al. [ 23 ] studied the stability of dislocation dipoles using density-functional theory. Carbon structures have been analyzed using a number of interatomic potentials ranging from harmonic potentials expressed in terms of force constants [ 24 , 25 ] to nonlinear potentials [ 26 , 27 ], more accurate but computationally costlier. The latter category includes reactive potentials, e.g., the Reactive Empirical Bond-Order (REBO) potential introduced by Brenner [ 28 ]. The addition of torsion and non-bonded interactions to this potential further resulted in the AIREBO potential [ 29 ]. Subsequently, Los et al. [ 30 ] developed the Long-Range Carbon Bond-Order Potential (LCBOP) and its second version, the LCBOPII potential [ 31 ], which accounts for interatomic interactions up to fourth neighbors. This latter potential is similar to AIREBO but presumed to be better suited to large displacements, extreme pressures and temperatures, liquid-solid phases and crystal lattices with topological defects. In this paper, we present evidence of strong steric interference in bi-layer graphene containing offset point dislocations resulting from the out-of-plane displacements of the individual layers. Calculations are carried on the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [ 32 ] of Sandia National Laboratories using the LCBOP [ 30 , 31 ] interatomic potential. We start by presenting an assessment of the fidelity of the LCBOP potential [ 30 , 31 ], cf. Section 2, in applications to graphene. In Section 2.1, we start by validating its harmonic response by comparing the predicted phonon dispersion curves to experimental data of Siebentritt et al. [ 33 ], Oshima et al. [ 34 ], Nicklow et al. [ 35 ] and Yanagisawa et al. [ 36 ] and those of other potentials. The requisite force constants are derived by linearization of the potential and are presented in full form in [ 37 ]. Subsequently, we turn to the anharmonic behavior of the LCBOP potential. Specifically, in Section 2.2 we evaluate the potential in configurations involving the formation of dislocation dipoles and quadrupoles in monolayer graphene and assess its ability to predict accurate and physically meaningful structures. With the LCBOP potential thus validated, we turn to the characterization of the core structure and energies of dislocation quadrupoles in bilayer graphene. Remarkably, we find that, if the dislocations in the individual layers are sufficiently offset, the out-of-plane displacements of the layers give rise to strong steric interactions which result in high energies. A summary of the main conclusions is presented in Section 4by way of closure. 2. Validation of the Theoretical Model For completeness and subsequent reference, we begin with a brief summary of the LCBOP potential [30,31]. The potential energy Ebof a system of Natoms is given by Eb=1 2 N ∑ i,jSdown sr,ij Vsr ij +Sup sr,ijVlr ij +1 Zmr i Sup mr,ijVmr ij , (1) where the short-range term Vsr ij gives the energy of the covalent bonds and Vlr ij =Vlr(rij) describes long-range interactions, with rij =|r r ri−r r rj| the distance between atoms i and j . Vmr ij is a term for bond breaking and formation based on ab initio calculations of the dissociation energy curves for single, double and triple bonds. It vanishes if the material is in a single phase, as for graphene. The S functions are switches expressed in terms of the Heaviside step function, Θ . Vsr is a Brenner type bond-order potential defined as Vsr ij =Vsr R,ij −BijVsr A,ij, (2) where Vsr R,ij and Vsr A,ij are repulsive and attractive radial pair potentials, Vsr R,ij =Asre−αri(3) Nanomaterials 2019,9, 1012 3 of 16 and Vsr A,ij =Bsr 1e−β1ri+Bsr 2e−β2ri, (4) respectively. The bond-order term Bij takes into account several many-body phenomena. It takes the form Bij =1 2bij +bji+Fconj ij +Tij, (5) where bij is a bond angle between first neighbors, Fij describes conjugate systems and Tij takes torsion into account. Long-range interactions are accounted for through the pair potential Vlr =hΘ(r0−r)Vlr 1(r) + Θ(r−r0)Vlr 2(r)iSdown lr (r), (6) where Vlr i(r)(i=1,2)are Morse functions, Vlr i(r) = eie−2λi(r−r0)−2e−λi(r−r0)+vi. (7) We note that all these functions are differentiable up to their first order and continuous up to their second order. Detailed expressions of the first and second derivatives of the potential up to fourth neighbors may be found in [ 37 ]. These derivatives are subsequently used in the calculation of force constants for the 18 neighbors of a reference atom. 2.1. Harmonic Response We begin by assessing the LCBOP potential as applied to graphene in the harmonic range. To this, we consider an infinite and perfect graphene lattice consisting of atoms, bonds and hexagonal cells [ 38 ], cf. Figure 1. We note that there are two types of atoms in the lattice, according to the collection of bonds connected to the atoms. Following standard notation [ 39 ], we label the atoms in the lattice with the double index (l l l , α) , where l l l= (l1 , l2 , l3) are integer lattice coordinates and α= 1,2 designates the type of atom. The lattice force constants are, then, Φ Φ Φl l l−m m m α β=∂2Etot ∂r r r(l l l,α)∂r r r(m m m,β), (8) where Etot is the total energy of the lattice and r r r(l α) are the spatial coordinates of atom (l l l , α) . The force constants depend on l l l−m m m only due to translation invariance. A systematic application of invariance under the symmetry group of graphene shows that the matrices Φ Φ Φ(l l l−m m m α β)are necessarily of the form Φ Φ Φi=−   αi0 0 0βi0 0 0 δi   i=1,3 , (9) Φ2=−   α2γ20 −γ2β20 0 0 δ2   , (10) Φ4=−   α4γ40 τ4β40 0 0 δ4   , (11) Nanomaterials 2019,9, 1012 4 of 16 where the subscripts represent the set of neighbors. Finally, an application of the discrete Fourier Transform furnishes the representation Φ Φ Φijl l l−m m m α β=1 (2π)2Zπ −πZπ −πb Φ Φ Φijθ θ θ α βe−iθ θ θ·(l l l−m m m)dθ1dθ2, (12) where the functions b Φ Φ Φij(θ θ θ α β)characterize the phonon dispersion curves of the material [39]. 121 2 3 1 Figure 1. 0-cells (atoms), 1-cells (atomic bonds) and 2-cells (hexagonal areas) in graphene. Table 1shows a comparison between the LCBOP constants, as computed from the expressions given in [ 37 ], and those put forth by Wirtz and Rubio [ 25 ], Tewary and Yang [ 27 ], Ariza et al. [ 38 , 40 ] and Mendez et al. [ 18 ]. As can be seen from the table, there are significant differences between the various potentials even in the harmonic range. Figure 2compares the phonon dispersion curves predicted by the LCBOP and the AIREBO potentials and the experimental data of Siebentritt et al. [ 33 ], Oshima et al. [ 34 ], Nicklow et al. [ 35 ] and Yanagisawa et al. [ 36 ]. As can be seen from the figure, the fidelity of the LCBOP and AIREBO phonon dispersion curves is comparable. Table 1. Parameters of the force constants models [ N/m ] by Wirtz et al. [ 25 ], Tewary et al. [ 27 ], Ariza et al. [38,40], Mendez et al. [18] and the present work. [25] [27] [38] [40] [18] Present Work α1399.0 409.7 364.0 527.7 497.2 423.6 β1135.7 145.0 247.0 68.1 173.7 144.3 δ1292.8 98.9 100.5 118.3 106.9 75.7 α2−79.6 −40.8 −30.8 5.8 −41.43 −6.5 β267.8 74.2 72.3 32.7 58.1 29.9 γ239.2 −9.1 −17.8 26.7 −3.0 −23.7 δ20.9 −8.2 −11.5 −16.9 −15.9 −8.8 α30.0 −33.2 0.0 −20.64 4.0 β30.0 50.1 0.0 34.51 −0.8 δ3−34.3 5.8 3.7 9.1 −0.8 α40.0 10.5 0.0 0.3 β40.0 5.0 0.0 0.0 γ40.0 2.2 0.0 0.1 τ40.0 −2.2 0.0 0.1 δ417.1 −5.2 −1.8 0.0 Nanomaterials 2019,9, 1012 5 of 16 0 500 1000 1500 T h i s w o r k A I R E B O P o t e n t i a l E m p i r i c a l v a l u e s G K M w ( c m - 1 ) G Figure 2. Comparison between the experimental point values of phonon dispersion obtained by Siebentritt et al. [ 33 ] (empty triangles), Oshima et al. [ 34 ] (empty squares), Nicklow et al. [ 35 ] (filled squares) and Yanagisawa et al. [ 36 ] (filled triangles); the phonon dispersion curves calculated using the AIREBO potential [40] and those calculated in the present work. 2.2. Lattice Defects We proceed to assess the fidelity of the LCBOP potential in a number of configurations including dipoles and quadrupoles in monolayer configurations. 2.2.1. Method of Analysis We introduce discrete dislocations in graphene by means of Mura’s theory of eigendeformations [ 41 ], as developed by Ariza et al. [ 42 , 43 ]. Specifically, the dislocations are introduced by slip on the three effective slip systems in graphene shown in Figure 3. Stable dipolar configurations are obtained when the unit slips occur across a chain of zig-zag bonds [ 43 ]. For an example, gliding along three consecutive zig-zag bonds leads to a 7-5-5-7 or Stone-Wales defect, cf. Figure 4. b1 b2 b3 Figure 3. Burgers vectors of graphene. Nanomaterials 2019,9, 1012 6 of 16 Figure 4. Unrelaxed Stone-Wales defect configuration generated by gliding along three consecutive zig-zag bonds. We seed calculations by solving for the equilibrium configuration of a prescribed distribution of lattice defects in the harmonic range using the force constants introduced in Section 2.1. Conveniently, such solutions can be obtained analytically in closed form in the Fourier domain. Next, in order to elucidate the fully-relaxed configurations of the defect, we use the harmonic atomic positions as initial conditions for a LAMMPS [ 32 ] calculation using the LCBOP [ 30 ]. For each defect configuration, we conduct two relaxation steps, as shown schematically in the flowchart in Figure 5. The first step consists of a molecular dynamics NVT relaxation at temperature T= 1K and imposing periodic boundary conditions. In the second step, a conjugate-gradient molecular-statics iteration is performed, giving the equilibrium configuration of the defects. In all cases, we have carried out studies in order to ensure convergence with respect to cell-size. For the harmonic calculations cell-size convergence is attained for 3360 atoms, while the anharmonic calculations require 18,720 atoms. Harmonic deformed lattice EH, uH Ei=0=EH ui=0=uHNVT MD relaxation CG energy minimization ΔE<ΔEmin Yes No Relaxed deformed lattice ENh, uNh ΔE=|Ei-Ei-1| ui=uH+u’ ui=uH+u’+u’’ (a) z x y σy=0 σx=0 Periodic Boundary Conditions σy=0 σx=0 Non-rigid PBC NVT relaxation CG minimization z x y T = 1K (b) Figure 5. ( a ) Flowchart of the nonlinear relaxation scheme and ( b ) outline of the boundary conditions for the nonlinear relaxation steps. 2.2.2. Dislocation Dipoles As a first validation case, we study a periodic dislocation dipole configuration of increasing length and compute the corresponding harmonic and fully-relaxed equilibrium configurations and energies. Nanomaterials 2019,9, 1012 7 of 16 We use the number n of gliding atomic bonds to measure the size of the dipoles. The linear size d , measured from the center of the rotating bonds in the undeformed lattice, follows as d=√3a(n− 1 )/ 2, where a=1.42 Å is the interatomic distance. We start by constraining the out-of-plane displacements. Under these conditions, the energy exhibits the expected logarithmic dependence with dipole size, Figure 6, in accordance with previous studies [ 38 ]. As expected, the anharmonic energies are much lower than the harmonic ones. Thus, the energy of the SW defect, or shortest dipole, decreases from 21.4 eV to 5.0 eV when the configuration is fully relaxed. This energy value is in keeping with—but lower than—the values 5.96 eV and 5.92 eV reported in [ 44 , 45 ]. It is also noteworthy that in the presence of defects the LCBOP is considerably softer than the semiempirical potential of Jain et al., which accounts for second-neighbor interactions only. 0 . 4 0 . 6 0 . 8 1 . 0 1 . 2 1 . 4 1 . 6 0 5 1 0 1 5 2 0 2 5 3 0 3 5 4 0 4 5 T h i s w o r k - h a r m o n i c T h i s w o r k - r e l a x e d 2 - D K . J a i n e t a l . 2 - D E n e r g y ( e V ) l o g (n) Figure 6. Dislocation energy as a function of dipole size, n , in 2-D with the force constants model, after relaxing harmonic configurations and presented by S.K. Jain et al. [44]. Remarkably, the fully relaxed 3-D dipole configurations exhibit two different stable modes: a symmetric, or S, mode in which the dislocation cores move in the same out-of-the plane direction; and an antisymmetric, or AS, mode in which the dislocation cores move in opposite out-of-plane directions, cf. [ 21 ], Figure 7. Since the stored energy per dislocation core of the Smode is approximately 10% higher than the corresponding AS mode, we regard the Smode as metastable, cf. [ 44 ]. In calculations, we prime the AS mode—and eschew the Smode—by applying an initial small out-of-plane displacement of the order of 0.2 Å to one atom in each core in opposite directions to each other. For the Stone-Wales defect, we compute a out-of-plane displacement range ( ∆z=zmax −zmin ) of 2.6 Å, comparable to the value of 2 Å reported by [46]. When out-of-plane displacements are allowed, the energy dependence on dipole size is also remarkably different than in the planar case, Figure 8. Thus, following an initial increase up to n= 13 approximately, the energy per unit periodic cell attains a constant value ostensibly independent of the dipole size. In particular, the logarithmic dependence characteristic of the planar solutions is lost. Correspondingly, the dislocation cores constrict when out-of-plane dislocations are allowed. The out-of-plane displacements thus effective localize the elastic field of the dislocations, which behave as non-interacting standing solitons. The results of Chen et al. [ 47 ] based on the REBO potential show similar behavior, although their maximum stored energy value is ∼ 2.5 eV higher, Figure 8. Nanomaterials 2019,9, 1012 8 of 16 This comparison shows that in the presence of defects the LCBOP is considerably softer than the REBO potential. 1.85 Å -2.28 Å (a) 3.29 Å -3.48 Å (b) Figure 7. Symmetric ( a ) and antisymmetric ( b ) stability modes in graphene. The color code indicates the out-of-plane position of the atoms. 0 1 0 2 0 3 0 4 0 5 0 4 6 8 1 0 1 2 T h i s w o r k 3 - D S . C h e n e t a l . 3 - D E n e r g y ( e V ) n Figure 8. Dislocation energy as a function of the size of the dipole for the 3-D fully relaxed configurations compared with the results presented by Chen et al. [47]. Nanomaterials 2019,9, 1012 9 of 16 2.2.3. Dislocation Quadrupole Dislocation quadrupolar arrangements are composed of two parallel dislocations of equal length and opposite Burgers vector, and thus encompassing four dislocation cores (Figure 9). Here, the periodic configurations are characterized by two parameters: the dipole length n and the dipole separation m , measured as the number of arm-chair bonds between the dipoles. The lowest energy state is always achieved when the out-of-plane displacements of the dislocation cores are the antisymmetric mode shown in Figure 10, corresponding to the periodic cell shown in Figure 9. n m Figure 9. Deformed configurations of periodic quadrupolar arrangement of discrete dislocations for the harmonic solution exhibiting pentagon– heptagon core structures (n=15, m=11). 6.58 Å -3.69 Å Figure 10. More stable mode of deformed configurations of periodic quadrupolar arrangement of dislocations in graphene. Proceeding as in the dipole calculations, we have computed the harmonic and nonharmonic atomic positions and stored energy per unit periodic cell of several quadrupolar configurations. We specifically aim to ascertain how the energy depends on the geometric of the quadrupoles, i.e., on n and m . The planar energies, Figure 11, exhibit the expected logarithmic or dipole-dipole elastic interactions, depending on geometry. Figure 12 shows the energy per unit periodic cell of the fully relaxed quadrupolar configurations. As in the case of dipoles, the out-of-plane displacements confer Nanomaterials 2019,9, 1012 16 of 16 35. Nicklow, R.; Wakabayashi, N.; Smith, H.G. Lattice Dynamics of Pyrolytic Graphite. Phys. Rev. B 1972 , 5, 4951–4962. [CrossRef] 36. Yanagisawa, H.; Tanaka, T.; Ishida, Y.; Matsue, M.; Rokuta, E.; Otani, S.; Oshima, C. Analysis of phonons in graphene sheets by means of HREELS measurement and ab initio calculation. Surf. Interface Anal. 2005 , 37, 133–136. [CrossRef] 37. Arca, F. Estudio de defectos topológicos en grafeno mediante un método no lineal basado en el potencial LCBOP. Ph.D. Thesis, University of Seville, Sevilla, Spain, 2019. 38. Ariza, M.P.; Ortiz, M. Discrete dislocations in graphene. J. Mech. Phys. Solids 2010,58, 710–734. [CrossRef] 39. Musgrave, M.J.P. Crystal Acoustics; Introduction to the Study of Elastic Waves and Vibrations in Crystals; Holden-Day: San Francisco, CA, USA, 1970. 40. Ariza, M.P.; Ventura, C.; Ortiz, M. Force constants model for graphene from Airebo potential. Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería 2011,27, 105–116. 41. Mura, T. Micromechanics of Defects in Solids; Kluwer Academic Publishers: Boston, MA, USA, 1987. 42. Ariza, M.P.; Ortiz, M. Discrete Crystal Elasticity and Discrete Dislocations in Crystals. Arch. Ration. Mech. Anal. 2005,178, 149–226. [CrossRef] 43. Ariza, M.P.; Serrano, R.; Mendez, J.P.; Ortiz, M. Stacking faults and partial dislocations in graphene. Philos. Mag. 2012,92, 2004–2021. [CrossRef] 44. Jain, S.; Barkema, G.; Mousseau, N.; Fang, C.; Huis, M. Strong Long-Range Relaxations of Structural Defects in Graphene Simulated Using a New Semiempirical Potential. J. Phys. Chem. C 2015 ,119, 9646–9655. [CrossRef] 45. Ma, J.; Alfè, D.; Michaelides, A.; Wang, E. Stone-Wales defects in graphene and other planar sp2 -bonded materials. Phys. Rev. B 2009,80, 033407. [CrossRef] 46. Jalali, S.; Jomehzadeh, E.; Pugno, N. Influence of out-of-plane defects on vibration analysis of graphene: Molecular Dynamics and Non-local Elasticity approaches. Superlattices Microstruct. 2016 ,91, 331–344. [CrossRef] 47. Chen, S.; Chrzan, D. Continuum theory of dislocations and buckling in graphene. Phys. Rev. 2011 ,84, 214103. [CrossRef] 48. Zhang, J.; Zhao, J. Mechanical properties of bilayer graphene with twist and grain boundaries. J. Appl. Phys. 2013,113, 043514. [CrossRef] 49. Jiang, B.Y.; Ni, G.X.; Addison, Z.; Shi, J.K.; Liu, X.; Zhao, S.Y.F.; Kim, P.; Mele, E.J.; Basov, D.N.; Fogler, M.M. Plasmon Reflections by Topological Electronic Boundaries in Bilayer Graphene. Nano Lett. 2017 , 17, 7080–7085. [CrossRef] [PubMed] 50. Yoo, H.; Engelke, R.; Carr, S.; Fang, S.; Zhang, K.; Cazeaux, P.; Sung, S.H.; Hovden, R.; Tsen, A.W.; Taniguchi, T.; et al. Atomic and electronic reconstruction at the van der Waals interface in twisted bilayer graphene. Nat. Mater. 2019,18, 448–453. [CrossRef] [PubMed] 51. Huang, S.; Kim, K.; Efimkin, D.K.; Lovorn, T.; Taniguchi, T.; Watanabe, K.; MacDonald, A.H.; Tutuc, E.; LeRoy, B.J. Topologically Protected Helical States in Minimally Twisted Bilayer Graphene. Phys. Rev. Lett. 2018,121, 037702. [CrossRef] c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).