Role of Graphene Surface Ripples and Thermal Vibrations in

Jul 23, 2019 - (10,13,14) The variety of wheels for nanocars stimulated the first international ..... Namely, when surface ripples down (in the negati...
0 downloads 0 Views 3MB Size
Subscriber access provided by UNIV OF LOUISIANA

C: Physical Processes in Nanomaterials and Nanostructures

Role of Graphene Surface Ripples and Thermal Vibrations in Molecular Dynamics of C 60

Seyedeh Mahsa Mofidi, Hossein Nejat Pishkenari, Mohammad Reza Ejtehadi, and Alexey V. Akimov J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.9b03947 • Publication Date (Web): 23 Jul 2019 Downloaded from pubs.acs.org on August 3, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Role of Graphene Surface Ripples and Thermal Vibrations in Molecular Dynamics of C60 S. Mahsa Mofidi,1,2 Hossein Nejat Pishkenari,*3 Mohammad Reza Ejtehadi,4 Alexey V. Akimov*2 1. Institute for Nanoscience and Nanotechnology (INST), Sharif University of Technology, Tehran, Iran 2. Department of Chemistry, University at Buffalo, State University of New York 14260-3000, United States 3. Mechanical Engineering Department, Sharif University of Technology, Tehran, Iran 4. Department of Physics, Sharif University of Technology, Tehran, Iran

ABSTRACT: Nanocars are artificial molecular machines with chassis, axles, and wheels designed for nanoscale transport at materials’ surfaces. Understanding the dependence of surface dynamics of nanocars on the substrate’s physico-chemical properties is critical to design of the transport properties of such man-made nanoscale devices. Among the multitude of potential substrates for the nanotransporters, graphene exhibits intrinsic ripples on its surface, which may affect the surface dynamics of nanocars. In this work, we report our molecular dynamics study of motion of C60, a popular nanocar wheel, on the graphitic substrates with systematically-controllable surface ripples. We find that surface ripples increase the amplitude of fullerene fluctuation in the direction normal to surface, which leads to decrease of the desorption temperature from 650 K on double layer graphite system with less ripples to 550 K on single layer graphene with more ripples. The surface diffusion of C60 follows the rare hops mechanism for temperatures up to 30 K. It switches

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

to continuous semi-ballistic motion at 150 K. The surface ripples do not significantly affect the diffusion coefficients, but change the anomaly parameter, especially at low temperatures. The ripples exhibit no major effect on the rotational dynamics of C60, which is attributed to very small energy barriers for C60 rotation on graphitic surfaces.

1. Introduction Inspired by natural tiny machines, several artificial molecular machines for nanoscale transport have been synthesized recently, including nanocars,1,2 nanomotors,3 nanorotors,4–6 trimers,7 or unimolecular submersible nanomachines (USNs).8 Among these systems, nanocars are intended to operate on surfaces and are composed of chassis, axles, and wheels.9 Designed in analogy to a macroscopic automobile, nanocars are able to exhibit motion on solid surfaces by rolling of their wheels.10,11 The progress in synthetic methods has enabled the use of variable types of molecular moieties as nanowheels: C60,9 p-carborane,12 and adamantane.10,13,14 The variety of wheels for nanocars stimulated the first international nanocars race: activated by a 4-tips STM, the nanocars raced on gold and silver surfaces at 4 K.15 The interest to transport properties of nanocars as a function of their structure (e.g. the choice of the wheels) fuels further fundamental research with the goal to understand the dynamics of nanomachines on various surfaces at atomistic level. Activated by thermal energy, molecules adsorbed on surfaces may display diffusive motion.16 Although realizing the directional motion of nanocar molecule over long distance still remains a challenge, the insights into surface dynamics of single nanocars molecules revealed by the imaging techniques like STM17, AFM, or single-molecule fluorescence microscopy18 provide valuable information to address this challenge in the future. The experimental imaging studies are complemented by atomistic modeling that aims to reveal the role of nanocar’s wheels and establish

2 ACS Paragon Plus Environment

Page 2 of 32

Page 3 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

the relationship between the wheels’ chemical structure and transport characteristics of nanocars. For instance, using the rigid-body molecular dynamics (MD), Akimov et al.11 supported the hypothesis that nanocar motion of 3 and 4-wheeled C60-based nanocars on gold involved the thermally-activated rolling motion of fullerene wheels. Few years later they showed, that the charge transfer occurring at the C60/gold interface plays a critical role in the electric field-driven directional motion of nonpolar nanocars on gold surface.19 Teobaldi et al.20 confirmed that electric field causes the preferred rolling of wheels due to complex interaction with metal surface. The mobility of C60 fullerene,21,22 p-carborane23 as well as fullerene-based24,25 and p-carborane-based26 nanocars on gold surfaces was studied using MD simulations. The surface diffusion of most molecules can be described by the stick-slip or sliding mechanism.27 However, for fullerene and other spherically-symmetric molecules like p-carborane the translational diffusion is strongly coupled to the rolling motion, which greatly contributes to the motion of such systems.28,29 Despite the previous studies provided a great wealth of information about nanocar and nanocar’s wheels diffusion, they focused on relatively rigid surfaces like metal20,24 and silica.10 In contrast, graphene mono-layer is able to form surface ripples, which lead to up to a 1 nm out-of-plane, zdirection substrate deformations.30 Graphene’s deformation is governed by its mechanical properties such as number of layers, Young’s modulus, and interfacial energy.31 Thus, the graphene sheets are not flat, exhibit a wavy morphology and nearly periodic surface ripples.32 The thermal-fluctuation-induced ripples change with time dynamically as observed by scanning tunneling microscopy (STM).33 The influence of ripples’ wavelength and amplitude on graphene’s mechanical and chemical properties was discussed previously.31 Using density functional theory, Cui et al.34 showed that ripples in multilayered graphene increased its chemical activity.34 By changing the ripple or

3 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

wrinkle density, the properties of graphene can be tuned. Recent synthetic and engineering advances enable creation of multiscale wrinkles on 2D graphene structures.35 Such wrinkles can be tailored to yield desired mechanical stiffness and surface adhesion in graphene. In addition, graphene sheets can be folded along certain angles,36 leading to various kinds of surface curvatures. It is logical to expect that surface diffusion of molecules on graphene can be also affected by the ripples. With the unique electrical and structural properties37 and a wide range of potential applications, graphene may be considered the surface of choice for nanoscale molecular transporting applications. In this case, the wavy ripples may alter the motion of nanocars due to occasional decrease of the contact level and interlock effects to barricade molecule translation.38 As a result, the motion of single molecules and nanocars on carbonic structures has stimulated a new wave of interest. Ozmaian et al.39 studied the diffusive motion of C60 on five different graphyne structures. Also other works considered C60 on graphene40,41 and graphene nano-ribbons42 substrates and characterized the thermally-induced motion fullerene. However, the effect of surface ripples on the motion of C60 and fullerene-based nanocars has not been explored yet. The present study pursues the fundamental question of the role of structural effects on molecular dynamics. Our main objective is to examine the effect of graphene surface ripples on diffusion of C60 molecule. To answer this question, we consider three types of substrates with different ripple extents: a) single layer graphene (SLG), with maximally pronounced surface ripples; b) double layer graphene (DLG), in which the ripples are suppressed naturally due to layer-layer interactions in comparison to SLG; c) a hypothetical frozen SLG (FLG), which does not exhibit any ripples by the construction. Using classical MD calculations of C60 on these substrates at different temperatures and by examining the C60/substrate potential energy surfaces, we characterize the

4 ACS Paragon Plus Environment

Page 4 of 32

Page 5 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

rotational and translational diffusion coefficients of fullerene, the activation energy barriers, and the diffusion anomaly parameter of these systems.

2. Computational Methodology To study the motion of a single-molecule C60 fullerene on graphene substrates, we utilize a combination of the potential energy surface (PES) analysis and the classical MD simulations. To be able to systematically vary the magnitude of graphene surface ripples and elucidate their effects, we define three types of substrates (Figure 1): A) This system consists of a single C60 molecule and a single layer of graphene (SLG). All of graphene’s atoms are allowed to move, which induces the surface ripple formation. This system is designed as the one with the maximally pronounced ripples, the effect facilitated by the absence of interactions with nearby layers. B) This system consists of a single C60 molecule and a double layer of graphene (DLG). The bottom layer of this DLG structure is kept fixed, whereas the top (closest to C60) is not constrained, similar to the setup A. The fixed atoms are also excluded from the bath (thermostat) calculations, but still, contribute to potential energy calculations. This design is intended to minimize the rippling of the top SLG, which is suppressed naturally by the layer-layer interactions. C) This system consists of a single C60 molecule and a single layer of graphene, which is kept frozen (FLG, for frozen layer graphene). This design corresponds to a hypothetic case of an ideally flat graphene surface, without ripples and without thermal motions.

5 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a)

(b)

(c) Figure 1. Representation of the three types of C60/substrate systems used in the present work. Fixed atoms are demonstrated in gray color. (a) System type A (C60/SLG), where the long-range ripples of graphene substrate can be observed. (b) System type B (C60/DLG), the bottom layer is kept fixed and ripples decreased. (c) System type C (C60/FLG), all substrate atoms are kept fixed. In all systems, the fullerene molecule is treated at the all-atomic level, with all internal vibrations enabled. The substrates are modeled as a square sheet of 12×12 nm2 constructed with 5744 carbon atoms in the system A and C, and 11488 atoms in the DGL of system B. The graphene sheet is positioned at a 𝑧 = 0 plane. For the case of the system B, the 𝑧 = 0 level corresponds to the upper layer. The C60 molecule is initially placed above the substrate, near the equilibrium distance, which

6 ACS Paragon Plus Environment

Page 6 of 32

Page 7 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

corresponds to 𝑧 ≈ 6.4 Å for the C60 center of mass (COM). Periodic boundary conditions are applied in x and y directions to allow unlimited diffusion of the C60 over the surfaces. Molecular interactions are described using classical force fields as implemented in the LAMMPS software.43 Tersoff potential is used to describe covalent bonds in graphene and C60 systems. Lennard-Jones 6-12 (LJ6-12) potential is used to describe the non-bonded interaction between each atom of graphene with each atom of C60: 𝑈𝑙𝐽 = 4𝜀

𝜎 12

𝜎 6

𝑟

𝑟

[( ) ― ( ) ]

(1)

𝑟 < 𝑟𝑐𝑢𝑡,

The parameters of  =3.4 Å and =2.41 meV are used for carbon van-der-Waals (vdW) radii and depth of C-C interaction potential, respectively. These values were used by Rafii-Tabar44 for sp2 C-C nonbonding interactions. The cut-off radius, 𝑟𝑐𝑢𝑡 of 12 Å is utilized in this work to reduce computational expenses. This selection is motivated by the commonly used criterion 𝑟𝑐𝑢𝑡 > 2.5𝜎. Using the vdW potentials, we also construct the PES corresponding to the interaction of rigid C60 and FLG substrate. The PES provides a quantitative description of potential modes of C60 motion on the substrate surface and has been widely utilized in the previous studies.22,23,42 In this work, we focus on lateral and vertical translational coordinates as well as pivoting coordinate of C60 with respect to the graphene lattice. The computed PES (Figure 2a) is periodic in both x and y directions with the energy barrier separating symmetrically-equivalent sites to be 1.4 meV, which corresponds to thermal energy at 𝑇 = 16 𝐾. The PES profile along the vertical shift coordinate suggests the equilibrium distance between C60 COM and graphene plane is about 6.4 Å with C60 hexagon facing down the surface. The desorption energy is on the order of 650 meV, which corresponds to thermal energy at ca. 𝑇 = 7500 𝐾 (Figure 2b). The desorption energy inferred from our PES analysis agrees well with the estimates of Savin et al.42 and Neek-Amal et al.,45 whose estimates are 730 and 760 meV, respectively.

7 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 32

400

C60/Graphene interaction 200 0 -200 -400 -600 -800 5.5

6

(a)

6.5

7

7.5

8

8.5

9

(b)

Figure 2. The C60/graphene PES along 3 coordinates: (a) the PES slice along the lateral displacement coordinates x and y with z value fixed at equilibrium value of z = 6.4 Å and the C60 oriented with its hexagon facing down the graphene surface; (b) the vertical displacement along zdirection with x and y coordinates fixed at center of the graphene sheet of 12×12 nm2 just above a graphene carbon and with the C60 hexagon facing down the graphene surface. MD trajectories are computed for each system for 40 ns using the velocity Verlet integration scheme with the nuclear integration time step of 1 fs. Nose-Hoover (NVT ensemble simulations) thermostat is used to control the temperature. The thermostat inertia parameter that determines the rate of heat exchange between the system and thermostat (Tdamp) is set to 100 fs, which is a typical value used in molecular simulations.46 For each system, the simulations are performed for temperatures ranging from 1K to 1000K. To compute the desorption temperatures, we utilize the following procedure. As Figure 2b suggests, the equilibrium distance between the C60 COM and outermost graphene surface is about 6.4 Å. The interaction energy becomes negligible when this distance is greater than 15 Å. Therefore, we treat the molecules that move away from the surface by 15 Å or more and never return during the simulation (40 ns) as desorbed. When we observe such “desorption” events, we

8 ACS Paragon Plus Environment

Page 9 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

reduce the temperatures by 25 K to examine nearby temperature values and test whether the desorption still occurs. If the desorption process is still observed, we lower our estimate of the desorption energy and repeat the process. Otherwise, we consider the obtained value as the resulting desorption temperature with 25 K error bar. To quantify the translational motion of C60, we compute mean square displacement (MSD) for every type of simulation (as defined by the system type and MD conditions). MSD quantifies the mobility of molecule due to random-walk-like motion. The MSD is used to compute the surface (2D) diffusion coefficient:47 𝑀𝑆𝐷(𝑡) = 〈(𝑥(𝑡) ― 𝑥(0))2 + (𝑦(𝑡) ― 𝑦(0))2〉 = 4𝐷𝑡𝛼.

(2)

Here, 𝑥 and 𝑦 are coordinates of the C60 center of mass. 𝐷 and 𝛼 are the 2D diffusion coefficient and the diffusion anomaly parameter, respectively. For normal diffusion, the latter parameter is close to unity, 𝛼 = 1, whereas for super- and sub-diffusion regimes 𝛼 > 1, and 𝛼 < 1, respectively. The angle brackets indicate the ensemble averaging, which is performed in the following manner. The initially obtained 40 ns trajectory is split into 60 intervals of 660 ps each. These intervals are regarded as 60 independent (sub)-trajectories, each started with distinct initial conditions, sampled from NVT ensemble (by the NVT MD of the original long trajectory). Finally, the averaging over these 60 sub-trajectories is used to compute the ensemble-averaged MSD values for the duration of up to 660 ps. The approach follows closely the recipe of Ernst and Kohler.48

3. Results and Discussion 3.1. Vertical motion As the PES analysis suggests, the diffusion of C60 on the ideally flat graphene surface is dictated by the barriers for lateral displacement. However, the displacement of C60 in the vertical direction may facilitate this process, since the lateral diffusion barriers would decrease when the C60 is

9 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 32

moved away from the surface. The latter type of motion is temperature activated. Figure 3 demonstrates the variation of z-component of C60 COM along representative trajectories, which corresponds to the vertical height of the molecule on top of the substrates in two different temperatures.

(a)

(b)

Figure 3. Variation of z coordinate of C60 center of mass during 500 ps for 3 system types at: (a) 30 K; (b) 300 K. As it is expected, at low temperatures (e.g. 30 K, Figure 3a), the C60 is mainly located 6.4 Å away from the surface, consistent with the predictions from PES (Figure 2b). Already at this temperature, one may observe the signs of surface rippling as manifested in larger C60 height fluctuation for the system of type A with flexible SLG as opposed to two other systems with more constraints on GL motion. The surface rippling effect in the system of type A is enhanced with temperature. Figure 3b illustrates that the z value of C60 COM varies from 5 to 7.5 Å range at high temperature. However, one may sight that these “turning points” on the PES profile would not be consistent with each other on the basis of energy: as Figure 2b suggests, the interaction energy at 𝑧 = 7.5 Å would be on the order of -400 meV, whereas for 𝑧 = 5.0 Å, the system would be in a highly repulsive region of the potential, with the energy much more positive than even 400 meV.

10 ACS Paragon Plus Environment

Page 11 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

This inconsistency can be explained if one considers that the z coordinate of the surface itself changes due to rippling. Namely, when surface ripples down (in the negative z-direction), the C60 molecule can dip together with it, leading to z coordinates of COM down to 5.0 Å. Yet, the C60graphene equilibrium distance may be notably larger than 5.0 A. To demonstrate that the rippling we observe in C60/graphene systems is an intrinsic property of graphene, we conduct the MD simulations of pure graphene, without C60. The calculations are repeated for a range of temperatures. Figure 4 illustrates the COM z component fluctuation, 𝜎𝑍 =

〈(𝑧(𝑡) ― 〈𝑧〉)2〉, as a function of temperature. The angle brackets here indicate the ensemble and trajectory-averaging. We observe a steady increase of the fluctuation magnitude with temperature. The fluctuation magnitude saturates to 1.0 Å at high temperatures, but reaches a substantial magnitude of ca. 0.6 Å already at room temperature. To show it is independent of the interatomic potential, we repeat these calculations for two distinct reactive potentials (only covalent C-C) parts: Airebo49 and Tersoff.50 In both cases, the dependencies appear similar and match with each other even quantitatively.

Figure 4. SLG ripple amplitude as a function of temperature. The results for two distinct bonding potentials Airebo and Tersoff, are shown. Inset: Variation of ripple amplitude during simulation

11 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

time for several temperatures: 1 K, 10 K, 100 K and 500 K as simulated using the Tersoff potential (the diagrams are shifted for better perception). The C60 COM z coordinate oscillation can lead to more probable desorption of molecule from substrate. The affinity for C60 to desorb from the substrate is examined based on the outcomes of the 40 ns MD simulations for each type of system in the range of temperatures from 1 K to 1000 K. For each combination of temperature and system, the simulations are repeated 3 times. Each repetition corresponds to a distinct initial velocities initialization (different seed numbers). Table 1 summarizes the lowest temperatures at which the C60 desorption occurs. It should be noted that any molecule tracked individually can, in principle, desorb from the surface at almost any temperature, given infinite time to evolve. It is, however, important to look at the statistics of such desorption events. In determining the desorption temperature, we look at the fraction of the desorbed molecules at every temperature. The desorption temperatures we report below correspond to temperatures when the fraction of desorbed trajectories changes rapidly from 0 to 1. We examine the dependence of such transition temperature on the MD trajectory length used in the analysis. We find that 10 and even 20 ns trajectories yield a wide range of temperatures at which the transition is observed, but both the 30 ns sub-trajectory and full 40 ns trajectory yield comparable values. We therefore believe 40 ns trajectories used in our calculations is sufficient to obtain accurate estimates of the upper limit of the desorption temperatures. Note that because of the 40 ns trajectory limit, there the values computed with longer runs may be slightly lower than the upper bounds we report here. Counterintuitively, the graphene surface rippling decreases the desorption temperature. For system A with the most flexible substrate, it is only 550 K. For the systems B and C, with the decreased degree of substrate’s flexibility, the temperature rises to 650 and 850 K, respectively. Experimentally, Ulbricht et al.51 demonstrated that the C60 desorption

12 ACS Paragon Plus Environment

Page 12 of 32

Page 13 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

from few-layered graphite occurs at 580 K, which falls in the 550 – 650 K temperature range dictated by SLG vs. semi-flexible DLG dynamics and interactions. The desorption temperature for the fully flat FLG with the rippling disabled (system C), is 850 K, which is notably larger than the experimental value. This indicates that the inclusion of the rippling effects is critical for obtaining accurate description of thermodynamics of C60 interaction with graphitic surfaces. Table 1. Estimated C60 desorption temperatures in each of 3 types of studied systems. Experiment51, Simulated system

C60 /SLG (A)

C60 /DLG (B)

C60 /FLG (C)

550±25

650±25

850±25

C60/graphite

Desorption 580

temperature, K

The counterintuitiveness of the relationship between the desorption energy and substrate’s flexibility originates from static thinking that “more bent” surface of the flexible graphene layer (e.g. SLG) may provide a better possibility for C60-SLG interactions, alike a “lock-key” or “hostguest” type. However, because the surfaces are being bent dynamically, such energeticallyfavorable configurations exist only transiently, quickly turning into energetically-unfavorable ones (e.g. with the C60 “sitting” on top of a hill formed in SLG). The potential energy profile (Figure 2b) fluctuates (the position of its minimum shifts back and forth) in time due to surface rippling, as opposed to a “constant” profile of an ideally flat substrate. As a result, the configurations that correspond to a more negative potential energy values at one times can find themselves at more positive values of potential energy at other times, due to occasional rippling of the substrate towards the C60 molecule that minimizes the vertical distance between the molecule and the surface. Hence, on average the potential energy of the system is more positive in the system with

13 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

a more flexible substrate. Consequently, the desorption energy decreases, leading to the observed C60 detachment at lower temperatures.

3.2. Lateral diffusion We next focus on lateral diffusion of C60 on graphitic substrates. The PES analysis (Figure 2a) reveals that molecule tends to slide in a horizontal direction such that the projection of COM passes the mid-point of a carbon-carbon bond of the lattice hexagons, since energy barrier through this route is the lowest. The estimated energy barrier between minimum and saddle points on PES for such translation is 1.4 meV (equivalent to the thermal energy of ca. 16 K), which is consistent with the value for sliding motion of C60 on graphene nanoribbons of 1.5 meV reported by Savin et al.42 Figure 5 panels a-c show representative trajectories of C60 molecule moving on SLG (system A) at temperatures ranging from 1 K to 500 K, with the upper limit selected to be below the desorption temperature threshold. As expected, the diffusivity of C60 increases with temperature. The diffusivity is quantified by the MSD vs. time (Figure 5, panels d-f). A close examination of the results suggests two notable features. First, for temperatures below 20 K, the C60 molecule does not show any considerable motion, with only a few accidental jumps to adjacent sites occurring. In addition, the

𝑑𝑀𝑆𝐷 𝑑𝑡

slope (diffusion constant) increases fast when transitioning from 𝑇 = 20 𝐾 to

𝑇 = 30 𝐾. This behavior is consistent with the PES-based estimates of the diffusion activation energy of 1.4 meV (ca. 16 K). Second, at higher temperatures, corresponding to 150 K, the MSD exhibits a power-law dependence on time, approaching the nearly quadratic (𝛼 = 1.8) limit at higher temperatures. It is due to the high kinetic energy of the molecule which overcomes the energy barrier of molecule-surface interactions easily. This type of motion can be attributed to

14 ACS Paragon Plus Environment

Page 14 of 32

Page 15 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

semi-ballistic motion and Levy flights (𝛼 > 1) in Eq. 2, for which the effect of surface structure diminishes.52 400

T=1 K T=5 K T=10 K T=20 K T=30 K T=35 K

200 0 -200

3000

T=35 K T=50 K T=60 K T=75 K T=100 K T=150 K

2000 1000 0

0

200

400

600

-2000 -2000

-1000

(a) 1500

500

0

0

1000

2000

104

4

2

400

600

(d)

0

0

1

400

10000

T=150 K T=200 K T=250 K T=300 K T=400 K T=500 K

2

200

5000

105

3

T=35 K T=50 K T=60 K T=75 K T=100 K T=150 K

3

200

0

(c)

1

0

-10000 -5000

(b)

T=1 K T=5 K T=10 K T=20 K T=30 K T=35 K

1000

0 -5000

-1000

-400 -200

T=150 K T=200 K T=250 K T=300 K T=400 K T=500 K

5000

600

(e)

0

0

200

400

600

(f)

Figure 5. (a-c) Representative trajectories of C60 motion on SLG at different temperatures. The length of each trajectory is 40 ns. (d-f) the mean square displacement of fullerene molecule for system A at different temperatures, corresponding to panels (a-c). To quantify the observed changes in the surface diffusion regime, we compute the diffusion coefficients, 𝐷, and the diffusion anomaly parameters, 𝛼, using the long-time limit of Eq. 2. It should be emphasized that sufficiently long trajectories must be used in the analysis. As Figure 6 suggests, the short-time

𝑑 𝑙𝑜𝑔(𝑀𝑆𝐷) 𝑑 log (𝑡)

slopes, which are directly related to the α parameter, are similar

for all temperatures used. However, they differ from each other for the trajectories longer than ca. 100 ps and don’t change later. At least, this is not observed on the nanosecond timescale of the

15 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 32

simulation. Only these data (beyond the first 100 ps) are used to compute the diffusion coefficients and diffusion regime parameters reported in this work. Most importantly, Figure 6 illustrates that the slopes change from relatively low numbers around 1.0 or below for temperatures below 𝑇 = 35 ― 50 𝐾 to larger numbers (e.g. 1.0-1.8) for higher temperatures, manifesting the change of the diffusion regime. 6 4 2 0

T=5 K T=10 K T=20 K T=30 K T=60 K T=100 K T=300 K T=500 K

-2 -4 -1

0

1

2

3

Figure 6. Log-log diagram of MSD vs. time for the C60/SLG system (A) at different temperatures. The slope is practically independent of temperature in the short simulation time (below 10 ps) limit, but is temperature-dependent in the long simulation time limit. Two dash lines indicate slope values of 1 (lower) and 2 (upper) at long simulation times. The diffusion coefficient and diffusion regime parameter computed for all three types of systems are shown in Figure 7. The diffusion coefficients are not significantly affected by the type of system – whether the ripples of SLG are allowed or not (Figure 7a). From the Arrhenius plot (Figure 7b), we estimate the activation energy of diffusion to be 1.6 meV at low temperatures and 39.2 meV at high temperatures. The low-temperature activation energy is consistent with the value of 1.4 meV obtained from the PES analysis (Figure 2a) and it composed mainly of the internal

16 ACS Paragon Plus Environment

Page 17 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

energy contribution. The increase of the activation barrier with the temperature likely reflects the entropic contribution to the Helmholtz free energy that defines the activation energy barriers. The activation energies are found to be practically independent of the surface model. Therefore, the entropic factors that contribute to the activation energies of C60 diffusion may be attributed to the activation of the C60, degrees of freedom (DOFs), in particular its rotational DOFs. This explanation is also consistent with the energy barrier of ca. 40 meV for the rolling motion activation we infer from the PES analysis (vide infra). The obtained diffusion coefficient is about 50 Å2/ps at room temperature. This values is comparable to the room-temperature diffusion coefficient on the order of 30 Å2/ps reported for C60 on a graphene sheet.41 The magnitude of the diffusion coefficient of C60 on graphene is certainly on a larger side in comparison to those of other systems. For instance, the room-temperature diffusion coefficient of gold nanoparticles on graphite surfaces on the order of 0.1 Å2/ps.53 The room-temperature diffusion coefficient of C60 on gold surfaces is estimated to be on the order of 10-2 Å2/ps. In the latter two cases, the difference in electronegativities of gold and carbon enable metal-to-carbon charge transfer. The increased Coulombic interactions lower the diffusion coefficients dramatically in comparison to the situation when no charge transfer is present. Indeed, in previous computational study, Akimov et al.20 demonstrated that the inclusion of charge transfer effects in the description of C60/gold interactions can lower the diffusion coefficients by as much as 4 orders of magnitude. Extending the analysis of the range of diffusion coefficients, Claytor et al.54 reported diffusion coefficient for carborane-based nanocars on glass surfaces on the order of 10-9 Å2/ps. In this situation, the carborane/silica interactions are dominated by strong hydrogen bonds55 and there are 4 wheels per molecule, which further lowers the diffusion coefficients well below what we compute for C60/graphene. Considering the above points, the values on the order

17 ACS Paragon Plus Environment

The Journal of Physical Chemistry

of 50 Å2/ps for C60/graphene diffusion coefficients are not that unreasonable – no charge transfer or polar/hydrogen bonds are present in such systems. The van der Waals interactions are relatively week. We further note that the molecular size is not necessarily a direct determinant of its diffusion coefficient, especially in the case of C60. If an increase of the molecule’s size is associated with its vertical (normal to surface) elongation, the atoms that are more distant from the surface contribute only little to the short-ranged van der Waals adsorbate-substrate interaction energies. In addition, the symmetry of the system along a particular degree of freedom may play a role. Adding more atoms to a system increases the overall van der Waals interactions and decreases the likelihood of the molecule’s desorption. On the contrary, symmetric organization of atoms (as in C60) lowers the effective diffusion energy barriers. The rolling motion of C60 proceeds in such a way that the number of broken vdW contacts on one side of C60 is approximately equal to the number of the formed contact of the same type. Such “compensation” effect minimizes the overall energy variation, resulting in smaller activation barrier for diffusion. The “compensation” mechanism was earlier discussed by Akimov et al.20 in the context of dynamical charge-transfer. Thus, we consider the symmetry of C60 being another important factors contributing to its large diffusion coefficient on graphene. 250

10 C60/SLG C60/DLG C60/FLG

200

D (Å2/ps)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 32

150

2 C60/SLG C60/DLG C60/FLG

5

1.5

0

1

-5

0.5

100 50 0

0

200

400

600

-10

0

0.05

0.1

0.15

0.2

0

C60/SLG C60/DLG C60/FLG 0

100

200

Temperature (K)

(a)

(b)

18 ACS Paragon Plus Environment

(c)

300

400

500

Page 19 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 7. Quantification of the C60 mobility on three substrates: (system A – blue, system B – dark red, system C - yellow): (a) Diffusion coefficient as a function of temperature; (b) Arrhenius plot to compute activation energies; (c) diffusion anomaly parameter, 𝛼. Unlike the diffusion coefficient, the regime of diffusion is notably affected by the type of substrate at low temperature, 𝑇 < 30 𝐾 (Figure 7c). In this temperature range, only the system A, with the most flexible SLG substrate and relatively more pronounced rippling, shows normal diffusion (𝛼~1), whereas the systems B and C with the rippling inhibited or totally excluded, 𝛼 < 1 characteristic to a sub-diffusion regime. We conclude that ripples have more pronounced effect at lower temperature and help the molecule move on the substrate. When temperature increases, the signs of super-diffusion (𝛼 > 1) start appearing, although the effect is insensitive to the type of substrate. At this temperatures range, thermal fluctuations of substrate atoms exceed the effect of ripples.

Figure 8. The evolution of x-coordinate of C60 COM for several temperatures. One can observe notable Levy flights in x-direction for at 300 K.

19 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The tendency for super-diffusive motion increases with temperature, with the 𝛼 parameter approaching 1.8. A closer examination of the C60 trajectories indicates the appearance of Levy flights already at room temperature (Figure 8). According to such types of motion, the molecule may travel unidirectionally over a significant distance in a relatively short period of time (e.g. see the dynamics at around 25 ns in Figure 8, red curve). Such a motion is essentially a semi-ballistic (1 < 𝛼 < 2) and caused by a sudden deposition of kinetic energy into a translational mode. To revert such a motion of the molecule, it needs to “collide” with enough thermostat quasi-particles to revert its momentum. Because each collision may change C60’s momentum in a random direction, a series of effective collisions are needed to stop the Levy flight. Because the Levy flights occur within a finite time interval, the 𝛼 parameter reduces on average in comparison to the value attributable to a purely ballistic motion.

3.3. Rotational motion Finally, we characterize the rotational dynamics of C60 on the substrates. We compute potential energy profiles for fullerene rotation around x (horizontal) and z (vertical) axes when it is positioned on top of distinct high-symmetry points of the planar substrate (Figure 9). In this study, we consider five such points, as illustrated in Figure 9a. These points correspond to the projections of C60 COM onto the surface. For each axis/projection point combination, we consider the initial positioning of C60 facing the surface with its hexagonal or pentagonal side. The distance between the C60 COM of the molecule and substrate is adjusted for every rotational angle, to be sure we sample minimal energy pathway along the coordinates studied. Figure 9b shows the potential energy profiles for C60 rotation around the horizontal axis for various positions of C60 (different colors). The figure illustrates that rotational energy barriers vary depending on both the orientation

20 ACS Paragon Plus Environment

Page 20 of 32

Page 21 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

and position of C60. For the horizontal rotation (rolling) axis, the maximal value reaches 40 meV (equivalent temperature of 460 K) – and corresponds to a situation when C60 is rolls fixed on top of a C-C bond (point 3, green line). For such location of the C60 COM, the potential energy profile contains several minima. Two of them correspond to the penta-down configuration (angle 120°), when C60 faces the SLG surface with one of its pentagons (e.g. as shown in the inset in panel 9e) and the hexa-down configuration (angle 160°), when C60 faces the SLG surface with one of its hexagons (e.g. as shown in the inset in panel 9f). These energy minima are separated by the highest-energy barrier along this coordinate at angle ca. 138° (as also better detailed in Figure 9c). The largest energy difference is between this maximum (ca. -0.61 eV) and the hexa-down configuration (ca. -0.65 eV). Correspondingly, the energy barrier for rotation hexapenta is about 40 meV. Analogous analysis can be performed for other positions of the C60 COM, leading to different numbers of the maximal energy variation. As one can see from Figure 9 (panels b-d), the largest energy variation corresponds to hexapenta transition for all COM positions tested, but only the magnitude varies. In particular, this number reduces down to 30 meV (equivalent temperature of 350 K) for the on-top (point 2) configuration, to 20 meV (equivalent temperature of 230 K) when C60 COM projects to the middle of graphene hexagon (point 1). Rolling of C60 also involves transition between two hexagon faces. We find that such transitions generally require less energy (Figure 9d) than the hexapenta transitions at equivalent COM positions. For instance, the barrier for hexahexa transition on top of the point 1 (blue line) is only 10 meV, whereas the corresponding value for the hexapenta transition is about twice as large. Thus, the rolling motion of C60 is facilitated when it is in a hexagon-down orientation and resides over the center of hexagons of SLG surface.

21 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 9. Analysis of energetics of rotation of C60 on graphene: (a) Positions of C60 center of mass projected onto graphene plane; (b) potential energy profiles for C60 horizontal rotation on SLG when C60 is fixed at each of the five points shown in (a); (c, d) zooms into potential energy profiles shown in panel (b). The range of rotational angles corresponds to C60 orientations with pentagon (c) or hexagon (d) facing down the SLG; (e, f) potential energy variation due to vertical rotation around z-axis when C60 is oriented with its pentagon (e) or hexagon (f) facing down the SLG. The

22 ACS Paragon Plus Environment

Page 22 of 32

Page 23 of 32

insets in panels (e) and (f) illustrate the view on the C60 molecule from SLG surface. They correspond to pentagon-down and hexagon-down orientations, respectively. The analysis of the PES along the vertical rotation axis suggests that pivoting motion of C60 should be easily activated even at very low temperatures. With the pentagon-down orientation of C60 (Figure 9e), the pivoting motion is essentially barrier-less for high-symmetry configurations (points 1, 2, and 3) or has a very small activation energy, on the order of 5 meV for slightlyasymmetric points (4 and 5). Such energy corresponds to temperature of 60 K. The hexagon-down orientation of C60 (Figure 9f) shows qualitative differences: the largest pivoting energy barrier (also on the order of 5 meV) is observed to the placement of C60 on top of the hexagon center (point 1), whereas for all other placements the barriers are smaller, but not negligible (as for some pentagon-down configurations). We conclude that the pivoting motion of C60 dominates its surface rolling motion, especially when C60 is oriented with its pentagon face down the substrate and it’s COM on the top of several high-symmetry points. Meanwhile, the hexagon-down orientation is more stable in relation to pentagon-down and molecule tend to have hexagon-down orientation at low temperatures. So at extremely low temperature the molecule may orient in hexagon-down on the surface with dominant pivoting motion. As Figure 10 shows, vertical rotation is the only major type of motion at low temperatures. At 10 K 300

Rot

z

200

2000 200 0 -200

0

20

40

Time (ns)

100 0 -100 0

10

20

30

40

Rotation (rad)

Rotx

At 300K C60 Angular velocity (rad/ns)

Rotation (rad)

400

C60 Angular velocity (rad/ns)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Rotx

1500

Rot

z

1000

5000 0 -5000

0

20

40

Time (ns)

500 0 -500 0

Time (ns)

10

20

Time (ns)

23 ACS Paragon Plus Environment

30

40

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a)

Page 24 of 32

(b)

Figure 6. Horizontal and vertical angular velocity of C60 on SLG relating to (a) a low and (b) a high temperature. Insets: C60 rotation at 10 K and 300 K around horizontal and vertical axis for 3 systems. At 10 K, C60 can rotate around vertical axis. The MD simulation has been performed for all three systems and repeated 3 times each. The rotational dynamics shown in Figure 10 corresponds to the SLG system. The results for other systems reveal no significant differences in rotational dynamics. The rotational motion of C60 is characterized by the temporal evolution of angular velocity, 𝜔, components along the MD trajectories. We focus only on the components that correspond to horizontal (rolling,𝜔𝑥) and vertical (pivoting,𝜔𝑧) rotation of the molecule around its center of mass. The insets in Figure 10 𝑡

show the corresponding integrated cumulative rotation angle 𝜙(𝑡) = ∫0𝜔(𝑡′)𝑑𝑡′ that characterize the overall rotation degree. We observe that although both 𝜔𝑥 and 𝜔𝑧 are non-negligible and fluctuate randomly, their integrals are different: the cumulative rotation angle 𝜔𝑧 exhibits an increase of its magnitude, suggesting the overall rotation of the C60 molecule around z axis. The degree of rotation as measured by the integral is strongly temperature-dependent. Of course, the spontaneous unidirectional motion is prohibited thermodynamically, so it is expected the spinning reversal will occur, but may require longer times than in our calculations. In contrast to the observed spinning around z axis, rolling (rotation around x axis) is significantly hindered, especially at low temperatures, these observations are also consistent with our PES analysis that suggests much higher activation barrier for the C60 rolling and almost no barrier for pivoting. We observe relatively similar trends of vertical or horizontal rotations in all three systems. So, the ripples exhibit no major effect on the rotational dynamics of C60. This can be attributed to very

24 ACS Paragon Plus Environment

Page 25 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

small energy barriers for C60 rotation on graphitic surfaces, which are demonstrated on energy diagrams of figure 9 (c-f) relating to PES analysis. In passing, we note that a typical ripple size in SLG may be as large as 2.5 nm at room temperature. Ripples of such large spatial extent are therefore unlikely to notably affect the rotational dynamics of C60, which is quite localized. They have a negligible effect on surface diffusion at high temperatures, when kinetic energy is sufficient to surpass any “barriers” created by ripples. At the same time, the ripples become more important at lower temperatures since at those conditions the energy barriers for C60 to move from one valley to another become comparable to the kinetic energy available. All these effects are indeed confirmed by our calculations. Importantly, one can extend this principle to formulate new rules for surface design in nanotransporting applications. In particular, it is reasonable to expect that small-scale defects such as highly-localized ripples, vacancies or adatoms, would have higher impact both on rotational motion and high-temperature surface diffusion of C60 than large-scale, smoothly-varying ripple defects.

4. Conclusions We demonstrate that the presence of graphene ripples significantly affects the variation of the C60/substrate vertical distance and consequently decreases the desorption temperature from 850 K in a hypothetic system with the ideally flat substrate to 550 K in a more realistic SLG system with ripples. We observe two regimes of diffusive motion of the C60 molecule. First, at temperatures below 30 K, the molecule exhibits rare jumps to the adjacent sites. Second, at temperature above 150 K, the molecule undergoes a semi-continuous motion with long Levy flights. In the middle temperature range, 20 – 150 K, a mixture of the two types of motion is observed. The motion

25 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

regime change is also reflected by the diffusion anomaly parameter, 𝛼, which increases up to 1.8 at the highest temperature before desorption. We estimate the activation energy for C60 surface diffusion is ca. 1.6 meV at low temperatures (< 20-30 K) and ca. 39.2 meV at higher temperatures. We observe that for most part of the temperature range examined, ripples have no major effects on the diffusion coefficient and the anomaly parameter because of the dominance of thermal fluctuations over the topology of the potential energy surface. At low temperatures, the presence of ripples increases the diffusion anomaly parameter, making the C60 to follow the normal diffusion as opposed to the sub-diffusion in ideally flat systems. We find that pivoting motion of C60 around the axis normal to the substrate occurs even at temperatures as low as 10 K. However, the rolling motion (around the axis parallel to the substrate’s surface) is not activated up until higher temperatures. We also demonstrate that the rotational energy barriers depend on both the orientation and position of C60. We observed that ripples have no significant effect on rotation. The computational protocols used in this work, the key input and output files, as well as important structural data are available online at https://github.com/AkimovLab/Project_C60_graphene_ripple. The repository also provides the digital equivalents of some figures shown in the manuscript.

Author Information Corresponding Authors * Email: [email protected] Twitter: @AkimovLab * Email: [email protected] Phone: +98 21 6616 5543 Notes The authors declare no competing financial interests. Acknowledgements

26 ACS Paragon Plus Environment

Page 26 of 32

Page 27 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Support of computations is provided by the Center for Computational Research at the University at Buffalo and by the high-performance computing facilities at the Information and Communication Technology Center of Sharif University of Technology. S.M.M. acknowledges the Ministry of Science, Research, and Technology of Iran for financial support. Also would like to thank Alireza Nemati, Seyyed Mohammad Lavasani and Brendan Smith for helping with the data analysis. A.V.A. acknowledges the financial support from the University at Buffalo, The State University of New York startup package. References

(1)

Jin, T.; García-López, V.; Chiang, P.-T.; Kuwahara, S.; Tour, J. M.; Wang, G. Enhancing Photostability of Fluorescent Dye-Attached Molecular Machines at Air– Glass Interface Using Cyclooctatetraene. J. Phys. Chem. C 2019, 123, 3011–3018. (2) Simpson, G. J.; García-López, V.; Petermeier, P.; Grill, L.; Tour, J. M. How to Build and Race a Fast Nanocar. Nat. Nanotechnol. 2017, 12, 604–606. (3) Kudernac, T.; Ruangsupapichat, N.; Parschau, M.; Maciá, B.; Katsonis, N.; Harutyunyan, S. R.; Ernst, K.-H.; Feringa, B. L. Electrically Driven Directional Motion of a Four-Wheeled Molecule on a Metal Surface. Nature 2011, 479, 208– 211. (4) Leigh, D. A.; Wong, J. K.; Dehez, F.; Zerbetto, F. Unidirectional Rotation in a Mechanically Interlocked Molecular Rotor. Nature 2003, 424, 174–179. (5) Sanada, K.; Ube, H.; Shionoya, M. Rotational Control of a Dirhodium-Centered Supramolecular Four-Gear System by Ligand Exchange. J. Am. Chem. Soc. 2016, 138, 2945–2948. (6) Tierney, H. L.; Baber, A. E.; Sykes, E. C. H.; Akimov, A.; Kolomeisky, A. B. Dynamics of Thioether Molecular Rotors: Effects of Surface Interactions and Chain Flexibility. J. Phys. Chem. C 2009, 113, 10913–10920. (7) Morin, J.-F.; Sasaki, T.; Shirai, Y.; Guerrero, J. M.; Tour, J. M. Synthetic Routes toward Carborane-Wheeled Nanocars. J. Org. Chem. 2007, 72, 9481–9490. (8) García-López, V.; Chiang, P.-T.; Chen, F.; Ruan, G.; Martí, A. A.; Kolomeisky, A. B.; Wang, G.; Tour, J. M. Unimolecular Submersible Nanomachines. Synthesis, Actuation, and Monitoring. Nano Lett. 2015, 15, 8229–8239. (9) Shirai, Y.; Osgood, A. J.; Zhao, Y.; Kelly, K. F.; Tour, J. M. Directional Control in Thermally Driven Single-Molecule Nanocars. Nano Lett. 2005, 5, 2330–2334. (10) Jin, T.; García-López, V.; Kuwahara, S.; Chiang, P.-T.; Tour, J. M.; Wang, G. Diffusion of Nanocars on an Air–Glass Interface. J. Phys. Chem. C 2018, 122, 19025–19036. 27 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(11) Akimov, A. V.; Nemukhin, A. V.; Moskovsky, A. A.; Kolomeisky, A. B.; Tour, J. M. Molecular Dynamics of Surface-Moving Thermally Driven Nanocars. J. Chem. Theory Comput. 2008, 4, 652–656. (12) Vives, G.; Kang, J.; Kelly, K. F.; Tour, J. M. Molecular Machinery: Synthesis of a “Nanodragster.” Org. Lett. 2009, 11, 5602–5605. (13) Chen, F.; García-López, V.; Jin, T.; Neupane, B.; Chu, P.-L. E.; Tour, J.; Wang, G. Moving Kinetics of Nanocars with Hydrophobic Wheels on Solid Surfaces at Ambient Conditions. J. Phys. Chem. C 2016, 120, 10887–10894. (14) García-López, V.; Alemany, L. B.; Chiang, P.-T.; Sun, J.; Chu, P.-L.; Martí, A. A.; Tour, J. M. Synthesis of Light-Driven Motorized Nanocars for Linear Trajectories and Their Detailed NMR Structural Determination. Tetrahedron 2017, 73, 4864– 4873. (15) Ariga, K.; Mori, T.; Nakanishi, W. Nano Trek Beyond: Driving Nanocars/Molecular Machines at Interfaces. Chem. – Asian J. 2018, 13, 1266–1278. (16) Wasio, N. A.; Murphy, C. J.; Patel, D. A.; Wei, D.; Sholl, D. S.; Sykes, E. C. H. Towards the Directional Transport of Molecules on Surfaces. Tetrahedron 2017, 73, 4858–4863. (17) Kelly, K. F. Probing Molecular Machines on Surfaces: The Nanocar and Beyond. Microsc. Microanal. 2008, 14, 956–957. (18) Vives, G.; Tour, J. M. Synthesis of Single-Molecule Nanocars. Acc. Chem. Res. 2009, 42, 473–487. (19) Akimov, A. V.; Kolomeisky, A. B. Unidirectional Rolling Motion of Nanocars Induced by Electric Field. J. Phys. Chem. C 2012, 116, 22595–22601. (20) Akimov, A. V.; Williams, C.; Kolomeisky, A. B. Charge Transfer and Chemisorption of Fullerene Molecules on Metal Surfaces: Application to Dynamics of Nanocars. J. Phys. Chem. C 2012, 116, 13816–13826. (21) Teobaldi, G.; Zerbetto, F. C60 on Gold: Adsorption, Motion, and Viscosity. Small 2007, 3, 1694–1698. (22) Nejat Pishkenari, H.; Nemati, A.; Meghdari, A.; Sohrabpour, S. A Close Look at the Motion of C60 on Gold. Curr. Appl. Phys. 2015, 15, 1402–1411. (23) Hosseini Lavasani, S. M.; Nejat Pishkenari, H.; Meghdari, A. Mechanism of 1, 12Dicarba-Closo-Dodecaborane Mobility on Gold Substrate as a Nanocar Wheel. J. Phys. Chem. C 2016, 120, 14048–14058. (24) Nemati, A.; Pishkenari, H. N.; Meghdari, A.; Sohrabpour, S. Directing the Diffusive Motion of Fullerene-Based Nanocars Using Nonplanar Gold Surfaces. Phys. Chem. Chem. Phys. 2018, 20, 332–344. (25) Nemati, A.; Meghdari, A.; Pishkenari, H. N.; Sohrabpour, S. Investigation into Thermally Activated Migration of Fullerene-Based Nanocars. Sci. Iran. 2018, 25, 1835–1848. (26) Hosseini Lavasani, S. M.; Nejat Pishkenari, H.; Meghdari, A. How Chassis Structure and Substrate Crystalline Direction Affect the Mobility of Thermally Driven PCarborane-Wheeled Nanocars. J. Phys. Chem. C 2019, 123, 4805–4824. 28 ACS Paragon Plus Environment

Page 28 of 32

Page 29 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(27) Rosei, F.; Schunack, M.; Naitoh, Y.; Jiang, P.; Gourdon, A.; Laegsgaard, E.; Stensgaard, I.; Joachim, C.; Besenbacher, F. Properties of Large Organic Molecules on Metal Surfaces. Prog. Surf. Sci. 2003, 71, 95–146. (28) Yamachika, R.; Grobis, M.; Wachowiak, A.; Crommie, M. F. Controlled Atomic Doping of a Single C60 Molecule. Science 2004, 304, 281–284. (29) Shirai, Y.; Minami, K.; Nakanishi, W.; Yonamine, Y.; Joachim, C.; Ariga, K. Driving Nanocars and Nanomachines at Interfaces: From Concept of Nanoarchitectonics to Actual Use in World Wide Race and Hand Operation. Jpn. J. Appl. Phys. 2016, 55, 1102A2 1-15. (30) Li, Y.; Liu, X.; Chen, C.; Duchamp, J.; Huang, R.; Chung, T.-F.; Young, M.; Chalal, T.; Chen, Y. P.; Heflin, J. R.; et al. Differences in Self-Assembly of Spherical C60 and Planar PTCDA on Rippled Graphene Surfaces. Carbon 2019, 145, 549–555. (31) Deng, S.; Berry, V. Wrinkled, Rippled and Crumpled Graphene: An Overview of Formation Mechanism, Electronic Properties, and Applications. Mater. Today 2016, 19, 197–212. (32) Bao, W.; Miao, F.; Chen, Z.; Zhang, H.; Jang, W.; Dames, C.; Lau, C. N. Controlled Ripple Texturing of Suspended Graphene and Ultrathin Graphite Membranes. Nat. Nanotechnol. 2009, 4, 562–566. (33) Xu, P.; Neek-Amal, M.; Barber, S. D.; Schoelz, J. K.; Ackerman, M. L.; Thibado, P. M.; Sadeghi, A.; Peeters, F. M. Unusual Ultra-Low-Frequency Fluctuations in Freestanding Graphene. Nat. Commun. 2014, 5, 3720–3726. (34) Cui, T. T.; Li, J. C.; Gao, W.; Jiang, Q. Geometric and Electronic Structure of Multilayered Graphene: Synergy of the Nondirective Ripples and the Number of Layers. Phys. Chem. Chem. Phys. 2018, 20, 2230–2237. (35) Lee, W.-K.; Kang, J.; Chen, K.-S.; Engel, C. J.; Jung, W.-B.; Rhee, D.; Hersam, M. C.; Odom, T. W. Multiscale, Hierarchical Patterning of Graphene by Conformal Wrinkling. Nano Lett. 2016, 16, 7121–7127. (36) Wang, B.; Huang, M.; Kim, N. Y.; Cunning, B. V.; Huang, Y.; Qu, D.; Chen, X.; Jin, S.; Biswal, M.; Zhang, X.; et al. Controlled Folding of Single Crystal Graphene. Nano Lett. 2017, 17, 1467–1473. (37) Meyer, J. C.; Geim, A. K.; Katsnelson, M. I.; Novoselov, K. S.; Booth, T. J.; Roth, S. The Structure of Suspended Graphene Sheets. Nature 2007, 446, 60–63. (38) Li, Z.; Liu, Z.; Sun, H.; Gao, C. Superstructured Assembly of Nanocarbons: Fullerenes, Nanotubes, and Graphene. Chem. Rev. 2015, 115, 7046–7117. (39) Ozmaian, M.; Fathizadeh, A.; Jalalvand, M.; Ejtehadi, M. R.; Allaei, S. M. V. Diffusion and Self-Assembly of C60 Molecules on Monolayer Graphyne Sheets. Sci. Rep. 2016, 6, 21911–21919. (40) Lohrasebi, A.; Neek-Amal, M.; Ejtehadi, M. R. Directed Motion of C 60 on a Graphene Sheet Subjected to a Temperature Gradient. Phys. Rev. E 2011, 83, 0426011–0426014. (41) Jafary-Zadeh, M.; Reddy, C. D.; Sorkin, V.; Zhang, Y.-W. Kinetic Nanofriction: A Mechanism Transition from Quasi-Continuous to Ballistic-like Brownian Regime. Nanoscale Res. Lett. 2012, 7, 148–155. 29 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(42) Savin, A. V.; Kivshar, Y. S. Transport of Fullerene Molecules along Graphene Nanoribbons. Sci. Rep. 2012, 2, 10121–10128. (43) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1–19. (44) Rafii-Tabar, H. Computational Modelling of Thermo-Mechanical and Transport Properties of Carbon Nanotubes. Phys. Rep. 2004, 390, 235–452. (45) Neek-Amal, M.; Abedpour, N.; Rasuli, S. N.; Naji, A.; Ejtehadi, M. R. Diffusive Motion of C 60 on a Graphene Sheet. Phys. Rev. E 2010, 82, 0516051–0516057. (46) Yu, T.-Q.; Alejandre, J.; López-Rendón, R.; Martyna, G. J.; Tuckerman, M. E. Measure-Preserving Integrators for Molecular Dynamics in the Isothermal–Isobaric Ensemble Derived from the Liouville Operator. Chem. Phys. 2010, 370, 294–305. (47) Michalet, X. Mean Square Displacement Analysis of Single-Particle Trajectories with Localization Error: Brownian Motion in an Isotropic Medium. Phys. Rev. E 2010, 82, 04191401–04191413. (48) Ernst, D.; Köhler, J. Measuring a Diffusion Coefficient by Single-Particle Tracking: Statistical Analysis of Experimental Mean Squared Displacement Curves. Phys. Chem. Chem. Phys. 2013, 15, 845–849. (49) Stuart, S. J.; Tutein, A. B.; Harrison, J. A. A Reactive Potential for Hydrocarbons with Intermolecular Interactions. J. Chem. Phys. 2000, 112, 6472–6486. (50) Tersoff, J. Modeling Solid-State Chemistry: Interatomic Potentials for Multicomponent Systems. Phys. Rev. B 1989, 39, 5566–5568. (51) Ulbricht, H.; Moos, G.; Hertel, T. Interaction of C 60 with Carbon Nanotubes and Graphite. Phys. Rev. Lett. 2003, 90, 0955011–0955014. (52) Luedtke, W. D.; Landman, U. Slip Diffusion and Levy Flights of an Adsorbed Gold Nanocluster. Phys. Rev. Lett. 1999, 82, 3835–3838. (53) Lewis, L. J.; Jensen, P.; Combe, N.; Barrat, J.-L. Diffusion of Gold Nanoclusters on Graphite. Phys. Rev. B 2000, 61, 16084–16090. (54) Claytor, K.; Khatua, S.; Guerrero, J. M.; Tcherniak, A.; Tour, J. M.; Link, S. Accurately Determining Single Molecule Trajectories of Molecular Motion on Surfaces. J. Chem. Phys. 2009, 130, 164710. (55) Kupchenko, I. V.; Moskovsky, A. A.; Nemukhin, A. V.; Kolomeisky, A. B. On the Mechanism of Carborane Diffusion on a Hydrated Silica Surface. J. Phys. Chem. C 2011, 115, 108–111.

30 ACS Paragon Plus Environment

Page 30 of 32

Page 31 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

TOC graphic

31 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

80x41mm (220 x 220 DPI)

ACS Paragon Plus Environment

Page 32 of 32