# Ballistic thermophoresis of adsorbates on free-standing graphene

###### Abstract

The textbook thermophoretic force which acts on a body in a fluid is proportional to the local temperature gradient. The same is expected to hold for the macroscopic drift behavior of a diffusive cluster or molecule physisorbed on a solid surface. The question we explore here is whether that is still valid on a 2D membrane such as graphene at short sheet length. By means of a non-equilibrium molecular dynamics study of a test system – a gold nanocluster adsorbed on free-standing graphene clamped between two temperatures apart – we find a phoretic force which for submicron sheet lengths is parallel to, but basically independent of, the local gradient magnitude. This identifies a thermophoretic regime that is ballistic rather than diffusive, persisting up to and beyond a hundred nanometer sheet length. Analysis shows that the phoretic force is due to the flexural phonons, whose flow is known to be ballistic and distance-independent up to relatively long mean-free paths. Yet, ordinary harmonic phonons should only carry crystal momentum and, while impinging on the cluster, should not be able to impress real momentum. We show that graphene, and other membrane-like monolayers, support a specific anharmonic connection between the flexural corrugation and longitudinal phonons whose fast escape leaves behind a 2D-projected mass density increase endowing the flexural phonons, as they move with their group velocity, with real momentum, part of which is transmitted to the adsorbate through scattering. The resulting distance-independent ballistic thermophoretic force is not unlikely to possess practical applications.

Thermophoresis is the phenomenon by which a body immersed in a fluid endowed with a temperature gradient experiences a force and, independent of convection, drifts from hot to cold Duhr and Braun (2006). We address here the less common case of thermophoresis of a physisorbed nano-object caused by an in-plane temperature imbalance in the underlying solid substrate surface. Recent years have seen a surge of interest for methods to control nanoscale transport and manipulation, also in view of potential applications in nano-devices. The possibility to drive directional motion of adsorbates by means of thermal gradients is interesting and has been explored both theoretically Schoen et al. (2006) and experimentally Barreiro et al. (2008). By a similar principle the controlled directional motion on graphene was also explored by means of strain or wettability gradients Wang and Chen (2015); Huang et al. (2014); Liu and Xu (2015). Carbon systems such as graphene and carbon nanotubes (CNTs) are prime candidate substrates Barreiro et al. (2008); Schoen et al. (2007); Becton and Wang (2014); Schoen et al. (2006); Zambrano et al. (2008); Becton and Wang (2014); Shiomi and Maruyama (2009); Savin and Kivshar (2012) for these phoretic phenomena, owing to the their remarkable mechanical strength and thermal conductivity. Computational studies have highlighted the possibility to drive thermally gold nanoparticles, water clusters, graphene nanoflakes, C60 clusters, and small CNTs over graphene layers or inside CNTs. In spite of that, there appears to be so far insufficient intimate understanding of that phenomenon, besides the obvious consensus that the driving force stems from the spatial non-uniformity of the phonon population Barreiro et al. (2008); Becton and Wang (2014); Prasad and Bhattacharya (2016). Our scope will be to deepen that understanding.

To begin with, in a macroscopic system the heat current density is proportional to the temperature gradient (Fourier’s law) through the thermal conductivity :

(1) |

Existing discussions of adsorbate thermophoresis similarly assume that the adsorbate mass current will similarly be proportional to the temperature gradient

(2) |

where now also depends on adsorbate, substrate and temperature - but not on its gradient. According to Eqs. 1-2, systems with different substrate lengths but same gradient should yield the same heat current and the same adsorbate current – that is, and should not depend on . These macroscopic expectations will of course be necessarily borne out for sufficiently large systems where transport is diffusive.

However, for system sizes smaller than or comparable with the mean-free-paths (MFP) of the heat carriers (phonons in our case), the heat transport may instead turn from diffusive to ballistic, where phonons with MFP larger than carry the heat. In graphene, can extend to hundreds of nm Ghosh et al. (2008) or even more Lindsay et al. (2014), making Eq. 1 invalid at the nanoscale. Ballistic heat transport along with ballistic-diffusive crossover have been discussed experimentally Bae et al. (2013) and theoretically Fugallo et al. (2014); Barbarino et al. (2015). Collective phonon excitations have also been proposed Fugallo et al. (2014) with MFPs possibly larger than those of single phonons, which might render the pure diffusive regime of Eq. 1 only attainable with samples of size 0.1–1 mm Xu et al. (2014). To describe the small size regime one can phenomenologically introduce a non-local generalization of Eq. 1 by replacing the constant with , with the resulting convolution leading to the product of Fourier-trasformed quantities Allen (2016a). We note that will however in our case depend on system size, a point to which we shall return later.

The question which we address here is what will happen to thermophoresis in the small size and distance regime, in particular whether Eq. 2 would still be valid or not in that case. As we shall see, a ballistic regime emerges for thermophoresis too at short distances, where Eq. 2 breaks down, and a new understanding is necessary. This understanding will be mandatory for thermally induced transport of matter at the nanoscale.

As a specific test case we consider here the thermophoretic force felt by a gold cluster physisorbed on a graphene sheet of length suspended between two baths at temperatures apart. Graphene has an extremely large heat conductivity, amongst the highest of any known material, with measured values ranging from 2600 to 5300 Balandin (2011); Barbarino et al. (2015). Thermal conductivity of suspended graphene is known to be dominated by acoustical lattice vibrations (even if the electron contribution to the total heat conductivity, estimated initially to be as low as 1% Ghosh et al. (2008), might be underestimated especially in doped and in short samples Kim et al. (2016)). The acoustical lattice vibrations of free graphene are in-plane transverse (TA), in-plane longitudinal (LA) and out-of-plane flexural (ZA). The contribution of the latter, much lower in frequency and therefore much more populated, has been shown to dominate in suspended graphene Lindsay et al. (2010); Wang et al. (2010); Lindsay et al. (2014).

The thermophoretic motion of gold nanoclusters in carbon nanotubes (CNT) has been associated to collective motions of the carbon atoms of CNT by Schoen et al. Schoen et al. (2007), and the subnanometer motion of cargos adsorbed on suspended nanotubes has been observed in a thermal gradient Barreiro et al. (2008). Very recently it was shown via phonon wave packet molecular dynamics simulations in CNTs that thermophoretic motion is driven by the scattering of LA phonons between the external and internal CNTs Prasad and Bhattacharya (2016). In all theoretical studies so far the temperature gradient was nonetheless considered as the relevant quantity which control thermophoresis as in Eq. 2.

Our present non-equilibrium molecular dynamics (NEMD) study between = 475K and = 325K shows that in the nanoscale regime below graphene sheet lengths 150 nm the accurately measured thermophoretic force acting on the gold clusters is not at all proportional to the temperature gradient. Instead, the force is found to depend basically only on the absolute temperature difference between heat source and sink, independently of the sheet size between them. That unmistakeable evidence of ballistic thermophoresis is, we further establish, associated with the known ballistic heat transport, caused chiefly by ZA flexural vibrations. The phoretic force arises due to scattering of these phonons on the adsorbed cluster which is immersed in the phonon flow. Momentum is tangibly transfered from these phonons to the cluster. That is surprising at first, since harmonic phonons are not supposed carry real momentum, that can only be carried by a moving mass.

As it turns out, a specific anharmonic process of flexural phonons in a flexible but nearly inextensible 2D membrane concentrates extra 2D-projected mass when it is corrugated, as the ZA phonons do. These phonons then carry real momentum as they move with the ballistic ZA phonon phase velocity. Some of that momentum is picked up by the adsorbed object which as a result is phoretically pushed from hot to cold regions. The overall ballistic thermophoresis of clusters proposed here, yet to be verified experimentally, appears to bear some similarity with that which enables fast diffusion of water clusters on graphene ripples Ma et al. (2016).

Our paper is organized as follows. First we describe the model system, gold clusters on a graphene sheet, and the simulation approach. The approximate temperature profile of suspended graphene clamped between two unlike temperatures is discussed, and some of its features related to the phonon flux as described by McKelvey-Shockley-type theory Maassen and Lundstrom (2015). Next, a gold cluster is adsorbed on graphene and initially shown to be freely diffusing in thermal equilibrium. After that, a left-to-right temperature difference = is turned on in the graphene sheet, and phoretic motion of the cluster is readily observed in the simulation. To measure accurately the phoretic force, the actual cluster motion is subsequently harnessed by a harmonic spring, acting as a dynamometer. The phoretic force so obtained and its dependence on is examined and found independent of up to at least 150 nm, indicating ballistic thermophoresis.

In order to gauge the ballistic flux of ZA phonons which appears as the driving agent of thermophoresis, the frequency-selected energy transmission spectrum of monochromatic flexural waves is examined. The ability of a ZA phonon mode to carry physical momentum is shown to occur as a result of the mass-carrying mechanism also associated with an anharmonically entangled LA mode, a mechanism best understood by viewing graphene as a nearly inextensible membrane. Finally, the thermophoretic force resulting from scattering on the adsorbed cluster of this “ZA+LA” complex is demonstrated. The gradient-independent character of this nanoscale phoretic force invites a short final discussion.

## I System and methods

Because the physical results to be reached in this work are entangled with technical aspects of the simulation, we find it best to describe the system and methods first thing here, rather than letting the reader wonder about them until later chapters.

Graphene is described by a C–C Tersoff potential, reparametrized to better reproduce the experimental phonon spectrum Lindsay and Broido (2010). A gold cluster ( = 459) with internal fcc structure and truncated-octahedral shape was described with Au–Au interactions of the embedded atom method (EAM) type, modified through a smooth cutoff Johnson (1988). The cluster 36-atom (111) facet is physisorbed on graphene. Interaction between gold and graphene atoms is assumed to be of Lennard-Jones type with = 22 meV and = 2.74 Å, as parametrized by Lewis et al. Lewis et al. (2000), a choice meant to reproduce the gold-graphite corrugation, rather than the adhesion energy.

We simulate a suspended graphene sheet of -length and -width 6.5 nm, with periodic boundary conditions along . The first C-atom row ( = 0) and the last one ( = ) are frozen, thus clamping the graphene sheet. A left-right temperature difference is introduced by coupling the first and the last 40 mobile atomic rows ( 4.5 nm) of the graphene sheet to two Langevin thermostats at temperatures and , respectively (see Fig. 1). We typically use 475 K, 325 K, and a Langevin damping coefficient of 10 ps strictly limited to the two thermostated regions, leaving the largest middle part of the graphene sheet (and the cluster when present) unthermostated. All initial equilibrium and subsequent NEMD simulations are conducted using our home-developed code. Two main approximations are the neglect of quantum effects, and of the finite-size acoustical phonon gaps at = 0.

Quantum effects are absent in our entirely classical simulations. For that very reason (besides practical ones, including experimental accessibility) we choose to work at 400 K, a temperature that compromises between three constraints: 1) it is well above the temperature where the quantum effects of ZA flexural phonons (the dominant thermophoresis agent) become irrelevant Bonini et al. (2012); 2) it is 1/5 of the LA (and TA) Debye temperatures Pop et al. (2012), so that their specific quantum effect, although not irrelevant, are at least not dominant. LA and TA modes will anyway turn out not to contribute to thermophoresis; 3) it is still below the higher temperature regimes where some phonon mean free paths get anharmonically shorter than the system size Chen et al. (2012).

The finite-size acoustical phonon gaps at = 0 on the other hand would be a problem if their magnitude came close to our working temperature = 400 K. Even for the smallest size considered, = 30 nm, however, the LA phonon gap /2L is only meV, or K (for 22 km/s); the ZA gap C/4 is for graphene meV or mK. At K both finite-size gaps are therefore irrelevant.

The adsorbate-free graphene sheet is simulated first. For an approximate evaluation of the local temperature we subdivide the graphene sheet in slices of 0.5 nm along the gradient direction , averaging the steady-state atomic temperature – obtained by the equipartition theorem, = – over the atoms in the slice and over a simulation time of at least 15 ns, see Fig. 2.

The cluster is then deposited near the center of the graphene sheet, following the adiabatic procedure described in Supporting Information Panizon et al. (), a protocol which also permits the direct calculation of the (temperature-dependent) adsorbtion free energy in thermal equilibrium. In these equilibrium conditions, the cluster undergoes thermal diffusion, both positional and angular Guerra et al. (2010). Once the temperature difference = - is turned on, the cluster is observed to drift from hot to cold, as shown in Fig. 3.

The approximate uniformity of the cluster center-of-mass (CM) motion indicates that the cluster-graphene friction is viscous Guerra et al. (2010). Given an average thermophoretic force and assuming a viscous friction coefficient , the equation of motion for the average CM velocity v is

(3) |

Assuming the thermophoretic force to be sufficiently large to dominate fluctuations, the force can be extracted from the initial acceleration of the cluster from a rest condition with 0 (see Supporting Information Panizon et al. ()). Alternatively, once both the heat flow and the adsorbate drift reach the steady state regime, the mean velocity = can be extracted Barreiro et al. (2008); Schoen et al. (2006, 2007); Zambrano et al. (2008); Rurali and Hernandez (2010).

Neither method is very precise, however. Past work Guerra et al. (2010) showed that is, in the diffusive regime, the time averaged result of a sequence of cluster jumps, rotations, and pinning events. The cluster needs a thermal fluctuation in order to misalign relative to the graphene lattice, and then it will positionally diffuse only when misaligned, events which occurs stochastically. Such erratic behavior, interesting as it is, requires very long simulation times, complicating the estimate of , and thus of in our case. Some improvement can be obtained by artificially locking the cluster orientation relative to the underlying graphene lattice into a poorly commensurate state, e.g. at = 30, so as to obtain an artificially small but well-defined and a large and smooth steady state speed suitable to an accurate extraction of . The clear difference between the two phoretic regimes is visible in Fig. 3. Different values would lead to a different and . Ignoring a possible angular dependence of momentum pickup rate by the cluster from the sheet, their product should be relatively independent of as desired.

In all cases however, the extraction of from the cluster trajectory is indirect, and thus tricky. It is therefore more convenient and immediate to measure directly the thermophoretic force from the NEMD simulation, by harnessing the mobile object, as also done by others Becton and Wang (2014). We artificially tie a spring to the cluster CM by means of a potential term , where is the cluster CM coordinate relative to the center of the sheet, and 0 is a spring constant, strong enough to restrict the cluster motion to the most representative central (unthermostated) part of the system. Typical trajectories of the harnessed cluster c.o.m. are shown in Fig 4. Experimentally, this kind of harnessed geometry might possibly be reproduced by an AFM tip, playing the role of our cluster, kept in place by a soft cantilever spring.

Owing to the modest cluster displacement permitted by the spring, can be obtained with the desired accuracy within a much shorter simulation time, = , where is the average cluster CM displacement along the gradient direction. To further reduce the error and accelerate the averaging, the cluster alignment angle can be kept fixed by constraining = 30. Before systematically doing that, we checked that the angular constraint does not influence the calculated thermophoretic force (see SI for more info Panizon et al. ()).

The statistical error affecting our results can be estimated both by the fluctuations of itself, = , and additionally checked by the fluctuations of the force along the direction perpendicular to the thermal gradient, = . Since the position of the cluster has large autocorrelation times due to the slow dynamics, a block-analysis has been performed. Our overall error in this procedure is estimated to be 0.5–1.0 pN, compared to 2.0–10 pN. In order to enhance the resolution on the calculated average force we set the spring constant = 0.001 meV/Å 16 N/m, a value much smaller than the typical 1 N/m of AFM cantilevers. However, we note that for a stiffer spring the same resolution could be achieved by just increasing the simulation time.

## Ii Temperature profile and heat flux

Let us begin with the results of a NEMD simulation of the free graphene sheet with the temperature imbalance, but without the adsorbed cluster. Fig. 2 shows the typical average steady-state temperature profile. The thermal gradient that will become relevant later is the slope, obtained as a linear fit, in the central region. The steep temperature jumps near the two thermostated regions are a common feature found in all NEMD simulations Zhou et al. (2009); Cao and Qu (2012); Becton and Wang (2014). A reasonable explanation generally offered for the jumps is the same as the Kapitza resistance jump Pollack (1969). The constrained left and right borders are mechanically different from the inside, and a travelling wave through that interface is partially reflected, causing a thermal resistance similar to that at the interface between two different materials Rurali et al. (2014).

In reality, the temperature jumps and their magnitude contain a much more interesting underlying source. Recent work by Maassen et al. Maassen and Lundstrom (2015) has shown that near-border jumps in the effective temperature profile are present even for “ideal” contacts – where no reflection occurs – when at least some of the thermal current is ballistic, in the following sense.

In standard Boltzmann theory of thermal transport, the net heat current is the difference between forward and backward currents, = - . Following McKelvey McKelvey et al. (1961) and Shockley Shockley (1962), the scattering between forward and backward currents is controlled by a parameter, , roughly reflecting the phonon mean free paths. In this respect, the evolution of the heat currents inside the sample depends on the scattering probability between the two opposite currents. Considering only phonons with a given energy ,

(4) |

The temperatures of the two thermostats enter as boundary condition in the values of the two currents. With an assumed ideal nature of the contacts, the right-flowing current through the left contact is in equilibrium with the left thermostat at , and similarly the left-flowing current through the right contact is in equilibrium at temperature , then

(5) | |||

(6) |

where is the distribution of modes of the thermal conductor at energy , is Planck’s constant, is the Bose-Einstein function. The coupled equations 5-6 for the currents can be solved to yield a value for the thermal gradient which is not, in general, equal to , thus implying the two temperature jumps at the two thermostated regions with ballistic transport variables. The ballistic thermal resistance equals in fact , where is the temperature jump between thermostated-unthermostated regions, assumed to be the same at left or right contacts since the thermostat efficiency is fixed (given by the Langevin damping rate and by the extent of the thermostated area). Note that this resistance is present even for ideal, non-reflective contacts.

Summing up, the border temperature jumps represent the prime evidence for at least some ballistic heat transport. If all the heat current was transported ballistically, with zero dissipation from hot to cold, then both temperature jumps would be (in a symmetric case) , and = 0 in between. If on the contrary the heat flux was completely diffusive, then the jumps would disappear, Fourier’s law would be obeyed throughout, and = everywhere. The temperature profile and the jumps in our simulated graphene are intermediate. The jumps are sizable, indicating a large amount of ballistic transport, besides some back-reflection. The jumps are clearly smaller than , so that the average gradient is nonzero, corresponding to some amount of dissipation taking place in the middle region. As can be seen in the Supporting Information, the border jumps decrease consistently for increasing , so that they will vanish for . With the new nanothermometric techniques Menges et al. (2016) it should be possible to verify this behavior.

In view of the above, the nonlocality of heat conductivity , very appropriately introduced by Allen Allen (2016b) to account for the nonlinear temperature distribution in a thermal hetero-contact turns out in our case to be connected with the ballistic fraction of the heat flow: the range of varying from infinity to zero (i.e., locality) from the ballistic regime at small to the diffusive limit at very large .

### ii.1 Cluster adsorption

Next, we adsorb on graphene the close-packed facet of a truncated octahedral cluster. Prior to studying its thermophoretic motion, it is of interest to consider the adsorbtion thermodynamics and to characterize its thermal diffusion in full equilibrium, at zero temperature gradient. Following early experimental studies Bardotti et al. (1996), previous simulations of the diffusion of such gold clusters on graphene or graphite found a strong correlation between the rotational and translation degrees of freedom of the cluster Luedtke and Landman (1999); Lewis et al. (2000); Guerra et al. (2010). At some angular orientations of the cluster the gold-graphene interface shows a strong interlocking, generally corresponding to minima of the zero-temperature adhesion energy (a more detailed discussion is reported in the SI Panizon et al. ()). At a majority of other orientations there is no interlocking, the adhesion is worse, and the cluster is much more mobile. In this ”incommensurate” state the translational barrier drops dramatically, allowing positional diffusion to occur. The overall diffusive motion of the cluster consists therefore of an alternation of locked states at specific angles where only vibrations take place, and of diffusive states where the rotationally depinned cluster executes fast translations Luedtke and Landman (1999); Guerra et al. (2010).

A second feature of the adsorbed cluster to be understood before turning on thermal gradients is the cluster-graphene interface free energy , measuring the cluster adhesion as a function of temperature, in full thermal equilibrium. Unlike bulk free energies, the (negative) interface free energy may upon heating either rise (commonest) or drop (rare, but possible). In either case, should happen to depend strongly enough upon , that dependence would in itself provide a source of thermophoretic force, once a spatial temperature gradient was introduced. We extract from thermodynamic integration of simulation data, by the standard method of gradually removing the graphene-cluster interaction, passing from a fully interacting system to a purely free-standing cluster (more information in the SI Panizon et al. ()). The adsorption free energy , a negative quantity, is found here to weaken in magnitude with increasing temperature. Due to this dependence of on temperature, a thermodynamic force is calculated as =- = - 0.29 pN for = 1 K/nm, a result that includes the small additional detachment of the cluster at higher temperatures. As it turns out, this is 15 to 40 times smaller than the value necessary to account for the actual phoresis, to be described below. Even without any pretense to accuracy, this is far too small. Moreover, this contribution to the thermophoretic force scales by construction as , which as we shall see is not compatible to our results. The origin of the main thermphoretic force must therefore be different.

## Iii Cluster thermophoretic motion and force

By turning on the left-right temperature difference across graphene, we simulate thermophoresis. Firstly, as a purely visual and extreme illustration, we present in Fig. 5 a few snapshots where a thermal corrugation front, starting from = 700 K, = 0 K, hits the cluster from the hot region, carrying it along as if swept along by a ”tidal” wave. Unrealistic as this sketchy situation is, it does provide an initial suggestion that the ZA flexural vibrations are responsible for the thermophoretic motion.

More realistic simulations carried out by the standard protocol of Section II confirm the thermophoretic cluster drift, the X-distance as a function of time shown in Fig. 3. Based on that drift we estimate, using Eq. 3, the thermophoretic force obtained from the evolution of the cluster velocity, with = 150 K and with = 30 constrained so as to improve lubricity, as explained earlier. By fitting the initial acceleration of the cluster (see Supporting Information Panizon et al. ()) we obtain estimates of in the range 6.79.3 pN. This direct approach is rather hard to converge for small sizes and not ideal for quantitative and accurate purposes, and as anticipated we proceed in the following to measure directly the average force by harnessing the cluster motion by a harmonic spring . Measuring, as a in a dynamometer, the mean displacement from the rest position at the midpoint of the sheet, the average thermophoretic force is extracted as . After choosing a convenient value which roughly optimizes the time needed to determine the cluster displacement we carry out a systematic series of harnessed cluster simulations for increasing graphene sheet size , two different temperature differences , and two different graphene orientations, obtaining the spring-constrained thermophoretic forces reported in Table 1.

(nm) | (K/nm) | (K) | (K) | (pN) |

30 | 0.42 | 38 | 40 | 2.1 0.5 |

40 | 0.39 | 37 | 40 | 1.7 0.6 |

50 | 0.33 | 37 | 40 | 1.9 0.7 |

60 | 0.29 | 37 | 40 | 3.6 0.6 |

70 | 0.28 | 37 | 40 | 2.9 0.8 |

30 | 1.88 | 148 | 150 | 9.7 0.5 |

40 | 1.47 | 147 | 150 | 9.2 0.9 |

50 | 1.26 | 147 | 150 | 8.9 0.5 |

60 | 1.10 | 147 | 150 | 10.4 0.6 |

70 | 1.04 | 147 | 150 | 10.1 0.6 |

150 | 0.54 | 147 | 150 | 8.8 0.6 |

30 | 1.95 | 147 | 150 | 9.0 0.4 |

40 | 1.71 | 147 | 150 | 9.8 0.6 |

50 | 1.27 | 147 | 150 | 9.3 0.5 |

60 | 1.15 | 147 | 150 | 8.7 0.4 |

70 | 1.03 | 147 | 150 | 7.9 0.5 |

The thermophoretic force displays, within tolerable errors, a linear dependence on the overall temperature difference but, remarkably, no dependence on the local gradient as would be expected from Eq. 2, which is valid in the diffusive regime. This is a striking hallmark of ballistic behaviour. Even at our largest graphene sheet size = 150 nm the expected decline of is still insignificant.

We conclude that up to at least the sizes considered in this work the thermophoretic force on absorbed nano-object on graphene remains largely ballistic. In turn, this matches our previous observation that heat flux has a large ballistic component as indicated by the graphene temperature profile. The dependence of the force on generally invoked in literature even for nanoscale sizes is therefore not borne out here.

The close parallel with ballistic heat transport of graphene suggests that the ZA flexural vibrations are the actors responsible for the thermophoretic force. To check that, we repeat simulations by artificially freezing the out-of-plane graphene motion, only allowing in-plane LA and TA vibrations. In that case we obtain, for = 50 nm and = 150 K, a thermophoretic force of 0.50.7 pN, negligible in comparison with the force found with unconstrained graphene, and more comparable to the error.

We conclude that the flexural ZA vibrations of graphene are the main source of thermophoretic force on the adsorbed cluster. This conclusion now invites a more detailed analysis.

### iii.1 Scattering of monochromatic flexural wave packets

We need to qualify the role of phonons in the observed ballistic thermophoresis. The connection between the thermophoretic force and the atomic vibrations which transport the heat is direct at the nanoscale sheet size, where phonons, whose mean-free path is larger, remain sufficiently well defined. The temperature difference gives rise to a net imbalance of phonon population between the two sides, and . In coaxial nanotubes, Prasad et al. Prasad and Bhattacharya (2016) showed how a phonon wave packet traveling in the outer (longer) tube scatters at the edges of the inner (shorter) tube, exchanging energy and momentum. This “push” from a single wave is felt identically whether it comes from the hot reservoir or from the cold reservoir. However, since the phonons that constitute the wave packets from the two opposite regions have different populations a net phoretic force arises. In the ballistic transport regime, where Eq. 1 is violated and dominates, the latter is a function only of , and so is the thermophoretic force resulting from this unbalanced population. In our case, the ballistic propagation of flexural phonons is identified as the source of the thermophoretic force. Harmonic phonons however should carry crystal momentum but no physical momentum: how can then the flexural phonons cause a net force?

We address this issue by simulating the time evolution of nearly monochromatic flexural phonon wave-packets. A phonon is injected in a free graphene sheet by applying an external vertical force on the left side of the non-thermostated region, shaking along a single column of carbon atoms with a force = for a certain number of cycles . The shaking produced two symmetric flexural wave packets, a left-traveling one which enters the left thermostat and gets dissipated, and another right-travelling across the graphene sheet, finally absorbed by the right thermostat. By measuring the heat absorbed by the right thermostat, alternatively without and with the adsorbed cluster, we calculate the variation of transmitted energy caused by the cluster. For each flexural frequency the heat transmission coefficient = is extracted in this manner. To allow for a computationally viable and sufficiently precise calculation we study running waves of just a few ( = 3–5) oscillations, so as to avoid interference with the backscattered wave, identifying the transmitted and reflected parts, finally absorbed by the opposite thermostats. Although ideally one could in the same way measure the flux of momentum besides energy, that is in practice substantially more difficult, and not really necessary. First, by injecting a small number of oscillations the resulting force on the cluster is often insufficient to overcome the static friction, and the cluster just oscillates around its average position. The magnitude of the oscillation is related to the momentum transfer between the wave packet and the cluster, but the net momentum transfer = is zero, absorbed in this regime by the lattice underneath. Secondly the flux of energy and its variation are well defined even in the absence of adsorbates, and can always be measured with good precision, unlike the total momentum flux, which as we will show in the following, has multiple components and cannot be measured with the same precision. We will therefore stick to the energy transmission, whose analysis is nonetheless quite informative.

The phonon-resolved energy transmission of graphene sheets of tranverse size = 9 nm and lengths in the range = 53–88 nm is shown in Fig. 6. At frequencies larger than 500 GHz the transmittance increases smoothly from 0.80 to 0.97. The loss of transmission can be attributed to all processes that produce backscattering.

In the low frequency region two transmittance dips are visible, approximately at = 49 GHz and = 175 GHz. These two features are connected with resonances between cluster-induced local vibrational modes and the incoming wave. The first mode, , coincides with the natural frequency of out-of-plane vibration of the cluster as a whole on the graphene substrate. This frequency depends on the Au–C interaction, and on the cluster size. Assuming approximately , a cluster mass and effective spring constant scaling proportional to and , respectively, one obtains . Using the value obtained for = 459 we obtain = GHz. The validity of this estimate is well verified for gold truncated octahedra clusters up to = 5635, where it yields = 33 GHz, close to the 31 GHz value obtained in simulation. This cluster -mode has the same symmetry as the ZA flexural modes of graphene, and is therefore linearly coupled to that continuum. The coupled problem actually provides an example of Fano resonance Fano (1961). The transmittance at a Fano resonance at = is of the general form

(7) |

where = , being the width of the resonance and an asymmetry parameter, ranging from infinity (when has a symmetric peak at = 0) to zero (when has a symmetric dip at = 0). The ZA phonon transmittance near can be fitted with 0.0 and 5 GHz. To understand the second resonance, we recall that adsorption of the cluster on graphene is accompanied at = 0 by a local “dimple” of the graphene sheet under the cluster contacting face Guerra et al. (2016). We find that the dip frequency corresponds to a wavelength which is about twice the size of this dimple under the cluster. Thus, when the wavelength matches the size of the cluster contact facet the flexural phonon can more effectively set it into oscillation, and the incresed scattering maximum is reflected in a transmittance minimum. This kind of process may be connected to the recent observation that the fast diffusion of water drops on graphene is dominated by the graphene thermal ripples which have wavelengths large enough to “contain” the drops Ma et al. (2016).

At both dips in the transmission spectra, therefore, the resonant coupling between the incoming phonon and the cluster enhances the scattering process, contributing to the overall increase of the cluster-induced thermal resistance. Interesting as it is, quantitative analysis however shows that the contribution of these resonances to the thermal resistance, and as a consequence to the thermophoretic force, is minor, compared to the remaining integrated continuum contribution over the whole frequency spectrum, which dominates. By integrating the two on-resonance parts, we find that resonances only contribute 1.5% of the total cluster-induced reduction of the transmitted energy.

In conclusion, all thermally excited flexural wavelengths, especially the longest, low-frequency ones, contribute to the transmitted heat flux. The same occurs to the momentum transfer to the cluster, which we describe next.

### iii.2 Physical momentum of flexural traveling waves

The energy-reflecting (and momentum-reflecting) nature of the cluster is first of all visible in the temperature profile shown in Fig. 2: the cluster acts as a localized thermal resistance, resulting in a minute drop of 2 K in the temperature profile. The cluster acts as a scatterer for the incoming phonon packets, and picks up momentum that way.

Here, as anticipated, we face a problem. In a completely harmonic setting phonons carry crystal momentum but not physical momentum Ashcroft and Mermin (1976). At the harmonic level, there should therefore be zero net momentum exchange with the cluster even in presence of scattering, and zero ballistic thermophoretic force. This is not the case, and anharmonicity must be involved in the thermophoretic force.

Indeed, according to standard theory, a phonon is associated to a physical momentum (different from crystal momentum ) due to anharmonicity. If = , is the anharmonic Gruneisen parameter, where is the volume, then = Lewis (1969); Lee and Lee (2006). For the flexural phonons of graphene, the Gruneisen parameter is indeed large. Alas, it is negative Mounet and Marzari (2005), which at first is puzzling.

While of course contained in the full derivations that can be found in the literature Michel et al. (2015) the physical reason underlying the negative Gruneisen parameter of graphene, or any other membrane, is actually simple and transparent. Thermally excited flexural phonons increase the total length of a clamped membrane causing a state of tension, therefore pulling the boundaries inwards - the contrary of normal bulk materials where temperature causes expansion. This connection between the negative Gruneisen parameter and the negative thermal expansion coefficient of graphene is well established Schelling and Keblinski (2003).

If applied blindly to our case, standard theory would thus suggest that flexural phonons carry a negative physical momentum, at odds with the positive thermophoretic force observed in simulation. The solution of this puzzle is, as it turns out, quite interesting. We find that a localized wave packet of flexural phonons in graphene conveys in reality a net positive physical momentum, part of which is communicated to the adsorbate. This does not violate, as it would seem at first, the Gruneisen requirement because of the anharmonic association of the ZA phonon with a LA phonon of double frequency, longer wavelength and much larger group velocity. The combined total momentum of the ZA-LA phonon pair is negative as in Gruneisen theory. However the two wave packets separate and the flexural phonon, the only one that scatters the adsorbate, carries a positive physical momentum.

To establish that, we apply a vertical, time-dependent force to the row of graphene carbon atoms at position 20 nm. This localized force creates two opposite-moving traveling flexural packets, whose oscillating parts are respectively of the form = and = for and . Consider now the forward-traveling phonon (the backward one behaves in exactly the same manner). Anharmonically associated with this harmonic ZA phonon there is an increase of the C–C bond lengths along , since the ZA corrugation extends the effective length of clamped graphene by (assuming ).

This forced bond length modulation gives rise to a LA phonon wavepacket of double frequency = , which travels away from with velocity and a small wavevector such that . LA phonons are orders of magnitude faster than ZA phonons, and in extremely short times (t 1 ps) the LA phonon carries away all longitudinal perturbations, leaving the remaining flexurally corrugated graphene with equilibrium C–C bond-lengths. The combination of corrugation and equilibrium bond-length produces a 2D projected density excess in the x-y plane , and therefore a positive physical momentum density associated the slow flexural phonon = .

Figs. 7 and 10 show as an example the case of = 80 GHz and a z-shaking force = 30 pN applied at 20 nm of the suspended graphene sheet (no adsorbed cluster for now). Upon increasing the intensity it is verified that the ZA corrugation scales linearly with , accompanied by an LA excitation amplitude which scales like , as expected for its anharmonic origin. In the first picoseconds after switching on of the force we monitor the projected bond length, defined as the distance between carbon atoms projected on the x-y plane, and the true C-C bond-length.

In Fig. 7 the in-plane projected – bond lengths are reported as a function of position and time. The propagating flexural ZA phonon involves a projected bond length which is shorter than the static equilibrium value, indicating that graphene is slightly “over-dense” under the slowly moving mode of velocity 0.9 km/s , close to the theoretical value = 1.06 km/s Pop et al. (2012). The reason for that higher density is clarified in Fig. 8, where the true bond lengths are shown. At very short times there are compensating “under-dense” modulations (bonds longer than the equilibrium value) which move away from the excitation point. Their apparent speed is 22 km/s, close to the experimental speed of sound for LA phonons 21 km/s Lindsay and Broido (2010). Note the scale difference of the two pictures: while the true bond length modulation is of 0.03%, the projected length shrinks by as much as 2%; and this for all ZA frequencies . The explanation for this mismatch lies in the fact that the modulations in the ZA and LA phonons are “diluted” over the relative wavelengths, which are very different: indeed and 7.9 nm. To better compare the two contribution we can estimate the momentum of the two modulations produced in a cycle of the vertical force. The ratio, for the simulation discussed above, can be calculated as . This result indicates that the total momentum is indeed negative as Gruneisen’s theory indicates (the factor 2 in the formula above comes from the double frequency).

Simulations thus confirm that the longitudinal modulation, although anharmonically created together with the flexural mode, moves away and separates very quickly. The remaining flexurally corrugated graphene has an increased mass density per unit projected area and its associated physical momentum is positive. This mechanism bears a resemblance with the picture by Bassett et al. Bassett and Pryce (1966) who investigated the physical momentum of localized running waves for 1D systems; to the best of our knowledge no further investigation of this mechanism was done following that work.

Consider now the interaction of these artificially excited phonons with the gold cluster adsorbed in the middle of the graphene sheet. The cluster produces intense scattering of the flexural wave packets, as visible from the projected bond lengths (Fig. 9). In this way, it picks up positive physical momentum, which explains the thermophoretic force. At the same time, there is no visible cluster-related scattering of the fast-moving longitudinal phonons. The LA cross section on the cluster is very small as can be seen comparing the evolution of the true bond length in simulations with and without the adsorbed cluster, Fig. 10 and Fig. 8: indeed the presence of the cluster has little effect until the ZA wave packet hits the cluster and additional anharmonic effects take place. As a result, the negative LA momentum is not transferred to the cluster.

### iii.3 Phoretic force by a flexural phonon

We are now in a position to present the resulting formulation for the phoretic force caused by a single injected flexural phonon in graphene, valid also for a more general membrane-like 2D material. While (predominantly) propagating from hot to cold, each flexural phonon en route impinges on the adsorbed cluster, which picks up momentum by parly reflecting it backwards. This can be described as follows.

The excess density embedded in a traveling ZA wave moves with its phase velocity. This excess projected density is

(8) |

where is the projection along of the average carbon-carbon bonds length whose equilibrium length is . The physical momentum density per unit area is therefore:

(9) |

It should be stressed again that this physical -momentum is not directly produced by the external force which generates the flexural phonon, a force that acts only along . The momentum is an anharmonic effect, proportional to the square of the ZA phonon amplitude , and the result of the anharmonic pairing of the ZA and an LA phonon of very different wavevector and velocity.

The phoretic force produced by partial back-reflection by the the adsorbate of contact y-width of a wavepacket of wavevector , phase velocity , amplitude , and transmission coefficient , is

(10) |

where the factor 2 accounts for backwards reflection of the scattered phonon (1D behavior assumed for simplicity).

To check this result in a specific case we carry out a specific simulation of the suspended graphene sheet of = 85 nm with the adsorbed cluster at 50 nm by applying a z-shaking force = 50 pN applied at 20 nm and frequency = 80 GHz. Here we can compare the above prediction with the phoretic force as measured by the initial acceleration of the free cluster.

From an independent, cluster-free simulation with the same z-shaking force = 50 pN we extract - after reaching steady state - an average = 1.5%. Assuming the momentum transmission coefficient to be well approximated by the energy transmission coefficient extracted earlier on, = 0.28 = , we obtain a phoretic force of 20 pN, in reasonably good agreement with the actual force acting on the cluster of 25 2 pN, directly obtained in simulation. In principle, by repeating the procedure for all ZA modes and integrating over the thermal distribution with its left-right imbalance one will approximate the total thermophoretic force. That however would demand a massive effort which we do not undertake. Since the single phonon force calculation works, there is no reason to doubt that the overall integrated thermophoretic force would also work.

## Iv Conclusions

An externally imposed temperature difference between the extremes of a vacuum suspended graphene sheet is predicted to push thermophoretically an adsorbed cluster from hot to cold. For submicron nanoscale sheet sizes the theoretical thermophoretic force is proportional to but independent of sheet length and thus independent of the thermal gradient , a key evidence of ballistic phoresis. The main agents of this effect are the flexural phonons, whose long mean-free path and ability to transport heat ballistically are well known. Besides heat, flexural phonons also carry physical momentum, owing to their membrane-like increase of projected density. In turn, that is associated with the special anharmonic coupling to a longitudinal phonon that takes the compensating density decrease away from the scattering region. The flexural phonon flux flowing from hot to cold cedes some of its physical momentum by scattering onto the adsorbed cluster, which is ballistically trasported. Previously known examples of ballistic phoresis include the Knudsen force exerted by gas particles on tips, described by Passian Passian et al. (2002) and observed by Gotsmann Gotsmann and Dürig (2005). The phonon-induced ballistic thermophoresis described here is specific to submicron sheet sizes and will eventually disappear as flexural phonons evolve from ballistic to diffusive once the sheet length is large enough, a crossover which our simulations do not yet detect at . Below that crossover size, the possibility to realize a distance-independent ballistic force is remarkable, and not unlikely to find practical applications. As an example the ballistic character of thermophoresis on graphene could allow long-range, non-contact action by a moving heat source such as a hot cantilever Menges et al. (2012) on adsorbed clusters or molecules. Other current uses of the thermophoretic effect at micro- and nano-scale, i.e. particle separation Lervik and Bresme (2014) or evaluation of molecular interactions Jerabek-Willemsen et al. (2011), could also benefit from the ballistic regime. Sensitive tools such as so-called “pendulum” Atomic Force Microscopes Gysin et al. (2011); Langer et al. (2014) or Quartz Crystal Microbalances Krim and Widom (1988); Bruschi and Mistura (2001) should be able to detect this effect, and to characterize the expected ballistic-diffusive crossover.

## Acknowledgments

Study conducted under ERC Advanced Grant 320796 MODPHYSFRICT, also partly sponsored by European COST Action MP1303.

## References

- Duhr and Braun (2006) S. Duhr and D. Braun, Proceedings of the National Academy of Sciences 103, 19678 (2006).
- Schoen et al. (2006) P. A. Schoen, J. H. Walther, S. Arcidiacono, D. Poulikakos, and P. Koumoutsakos, Nano letters 6, 1910 (2006).
- Barreiro et al. (2008) A. Barreiro, R. Rurali, E. R. Hernandez, J. Moser, T. Pichler, L. Forro, and A. Bachtold, Science 320, 775 (2008).
- Wang and Chen (2015) C. Wang and S. Chen, Scientific reports 5 (2015).
- Huang et al. (2014) Y. Huang, S. Zhu, and T. Li, Extreme Mechanics Letters 1, 83 (2014).
- Liu and Xu (2015) Q. Liu and B. Xu, Langmuir 31, 9070 (2015).
- Schoen et al. (2007) P. A. Schoen, J. H. Walther, D. Poulikakos, and P. Koumoutsakos, Applied Physics Letters 90, 253116 (2007).
- Becton and Wang (2014) M. Becton and X. Wang, Journal of chemical theory and computation 10, 722 (2014).
- Zambrano et al. (2008) H. A. Zambrano, J. H. Walther, P. Koumoutsakos, and I. F. Sbalzarini, Nano letters 9, 66 (2008).
- Shiomi and Maruyama (2009) J. Shiomi and S. Maruyama, Nanotechnology 20, 055708 (2009).
- Savin and Kivshar (2012) A. V. Savin and Y. S. Kivshar, Scientific reports 2 (2012).
- Prasad and Bhattacharya (2016) M. V. Prasad and B. Bhattacharya, Nano letters 16, 2174 (2016).
- Ghosh et al. (2008) S. Ghosh, I. Calizo, D. Teweldebrhan, E. Pokatilov, D. Nika, A. Balandin, W. Bao, F. Miao, and C. N. Lau, Applied Physics Letters 92, 151911 (2008).
- Lindsay et al. (2014) L. Lindsay, W. Li, J. Carrete, N. Mingo, D. Broido, and T. Reinecke, Physical Review B 89, 155426 (2014).
- Bae et al. (2013) M.-H. Bae, Z. Li, Z. Aksamija, P. N. Martin, F. Xiong, Z.-Y. Ong, I. Knezevic, and E. Pop, Nature communications 4, 1734 (2013).
- Fugallo et al. (2014) G. Fugallo, A. Cepellotti, L. Paulatto, M. Lazzeri, N. Marzari, and F. Mauri, Nano letters 14, 6109 (2014).
- Barbarino et al. (2015) G. Barbarino, C. Melis, and L. Colombo, Physical Review B 91, 035416 (2015).
- Xu et al. (2014) X. Xu, L. F. Pereira, Y. Wang, J. Wu, K. Zhang, X. Zhao, S. Bae, C. T. Bui, R. Xie, J. T. Thong, et al., Nature communications 5 (2014).
- Allen (2016a) P. B. Allen, arXiv preprint arXiv:1612.01173 (2016a).
- Balandin (2011) A. A. Balandin, Nature materials 10, 569 (2011).
- Kim et al. (2016) T. Y. Kim, C.-H. Park, and N. Marzari, Nano letters 16, 2439 (2016).
- Lindsay et al. (2010) L. Lindsay, D. Broido, and N. Mingo, Physical Review B 82, 115427 (2010).
- Wang et al. (2010) Z. Wang, R. Xie, C. T. Bui, D. Liu, X. Ni, B. Li, and J. T. Thong, Nano letters 11, 113 (2010).
- Ma et al. (2016) M. Ma, G. Tocci, A. Michaelides, and G. Aeppli, Nature materials 15, 66 (2016).
- Maassen and Lundstrom (2015) J. Maassen and M. Lundstrom, Journal of Applied Physics 117, 035104 (2015).
- Lindsay and Broido (2010) L. Lindsay and D. Broido, Physical Review B 81, 205441 (2010).
- Johnson (1988) R. Johnson, Physical Review B 37, 3924 (1988).
- Lewis et al. (2000) L. J. Lewis, P. Jensen, N. Combe, and J.-L. Barrat, Physical Review B 61, 16084 (2000).
- Bonini et al. (2012) N. Bonini, J. Garg, and N. Marzari, Nano letters 12, 2673 (2012).
- Pop et al. (2012) E. Pop, V. Varshney, and A. K. Roy, MRS bulletin 37, 1273 (2012).
- Chen et al. (2012) S. Chen, Q. Wu, C. Mishra, J. Kang, H. Zhang, K. Cho, W. Cai, A. A. Balandin, and R. S. Ruoff, Nature materials 11, 203 (2012).
- (32) E. Panizon, R. Guerra, and E. Tosatti (????).
- Guerra et al. (2010) R. Guerra, U. Tartaglino, A. Vanossi, and E. Tosatti, Nature materials 9, 634 (2010).
- Rurali and Hernandez (2010) R. Rurali and E. Hernandez, Chemical Physics Letters 497, 62 (2010).
- Zhou et al. (2009) X. Zhou, S. Aubry, R. Jones, A. Greenstein, and P. Schelling, Physical Review B 79, 115201 (2009).
- Cao and Qu (2012) A. Cao and J. Qu, Journal of Applied Physics 111, 053529 (2012).
- Pollack (1969) G. L. Pollack, Reviews of Modern Physics 41, 48 (1969).
- Rurali et al. (2014) R. Rurali, X. Cartoixà, and L. Colombo, Phys. Rev. B 90, 041408 (2014), URL http://link.aps.org/doi/10.1103/PhysRevB.90.041408.
- McKelvey et al. (1961) J. McKelvey, R. Longini, and T. Brody, Physical Review 123, 51 (1961).
- Shockley (1962) W. Shockley, Physical Review 125, 1570 (1962).
- Menges et al. (2016) F. Menges, P. Mensch, H. Schmid, H. Riel, A. Stemmer, and B. Gotsmann, Nature communications 7 (2016).
- Allen (2016b) P. B. Allen, arXiv preprint arXiv:1612.01173 (2016b).
- Bardotti et al. (1996) L. Bardotti, P. Jensen, A. Hoareau, M. Treilleux, B. Cabaud, A. Perez, and F. C. S. Aires, Surface science 367, 276 (1996).
- Luedtke and Landman (1999) W. Luedtke and U. Landman, Physical review letters 82, 3835 (1999).
- Fano (1961) U. Fano, Physical Review 124, 1866 (1961).
- Guerra et al. (2016) R. Guerra, E. Tosatti, and A. Vanossi, Nanoscale 8, 11108 (2016).
- Ashcroft and Mermin (1976) N. W. Ashcroft and N. D. Mermin, Solid state physics (Holt, Rinehart and Winston, 1976).
- Lewis (1969) M. Lewis, Physics Letters A 29, 644 (1969).
- Lee and Lee (2006) Y. C. Lee and W. Z. Lee, Phys. Rev. B 74, 172303 (2006), URL http://link.aps.org/doi/10.1103/PhysRevB.74.172303.
- Mounet and Marzari (2005) N. Mounet and N. Marzari, Physical Review B 71, 205214 (2005).
- Michel et al. (2015) K. H. Michel, S. Costamagna, and F. M. Peeters, Phys. Rev. B 91, 134302 (2015), URL https://link.aps.org/doi/10.1103/PhysRevB.91.134302.
- Schelling and Keblinski (2003) P. K. Schelling and P. Keblinski, Phys. Rev. B 68, 035425 (2003), URL https://link.aps.org/doi/10.1103/PhysRevB.68.035425.
- Bassett and Pryce (1966) I. Bassett and M. Pryce, Physical Review 150, 640 (1966).
- Passian et al. (2002) A. Passian, A. Wig, F. Meriaudeau, T. Ferrell, and T. Thundat, Journal of applied physics 92, 6326 (2002).
- Gotsmann and Dürig (2005) B. Gotsmann and U. Dürig, Applied Physics Letters 87, 194102 (2005).
- Menges et al. (2012) F. Menges, H. Riel, A. Stemmer, and B. Gotsmann, Nano letters 12, 596 (2012).
- Lervik and Bresme (2014) A. Lervik and F. Bresme, Physical Chemistry Chemical Physics 16, 13279 (2014).
- Jerabek-Willemsen et al. (2011) M. Jerabek-Willemsen, C. J. Wienken, D. Braun, P. Baaske, and S. Duhr, Assay and drug development technologies 9, 342 (2011).
- Gysin et al. (2011) U. Gysin, S. Rast, M. Kisiel, C. Werle, and E. Meyer, Review of scientific instruments 82, 023705 (2011).
- Langer et al. (2014) M. Langer, M. Kisiel, R. Pawlak, F. Pellegrini, G. E. Santoro, R. Buzio, A. Gerbi, G. Balakrishnan, A. Baratoff, E. Tosatti, et al., Nature materials 13, 173 (2014).
- Krim and Widom (1988) J. Krim and A. Widom, Physical Review B 38, 12184 (1988).
- Bruschi and Mistura (2001) L. Bruschi and G. Mistura, Physical Review B 63, 235411 (2001).