Multilayer Structures of Graphene and Pt Nanoparticles - a Multiscale Computational Study

In order to calculate the Pt-graphene interaction using DFT we consider five-layer Pt(111) slabs adhering to monolayer graphene. Both Pt slabs and graphene sheet are oriented parallel to the $xy$ plane of a Cartesian coordinate systems. Periodic boundary conditions are used in all three spatial directions with a 10 Å vacuum gap between the periodic replicas in the $z$ direction. Two variants of this structure (graphene on Pt (111) surface) are constructed as shown in Figure 1; for model a, a larger supercell consisting of p $(7\times 7)$ graphene and p $(6\times 6)$ Pt slabs were used in the calculations to mimic an incoherent interface. The structure for model b consists of p $(2\times 2)$ graphene and three atoms per Pt layer and represents a coherent interface as considered in previous studies [36, 21]. For both models, the Pt lattice parameters were set to experimental values for bulk Pt.

All the calculations were carried out by the Quantum Espresso package version 6.1 [37], employing PBE-based projected augmented wave(PAW) potentials [38]. Previous studies have shown that the optB88-vdw [40] exchange-correlation functional gives a reasonably accurate prediction of both interlayer distance and binding energy for graphite [39] , so here we choose the optB88-vdw functional in all our DFT calculations. In our calculations, the kinetic energy cutoffs of 40 Ry and 450 Ry were used for wave function and charge density calculations respectively, and a Methfessel-Paxton smearing of 0.01 Ry was used for the electronic convergence. The convergence for self-consistency calculations was less than 0.0001 Ry. In the calculations of model a, only the gamma k-point was used due to the large system we chose, while an $8\times 8\times 1$ k-points grid was used in the calculations of model b.

Since our aim is to obtain data for parameterizing an interaction potential, rather than performing an accurate analysis of the interface structure and energy using DFT, we do not carry out any structural relaxation. Instead we simply separate the graphene sheet rigidly from the Pt(111) slab without geometry relaxation of either the Pt slab or the graphene. We evaluate the binding energy per carbon atom as

$\displaystyle E_{\rm b}=[E_{\rm tot}-(E_{\rm pt}+E_{\rm gr})]/n$

(1)

$E_{\rm tot}$ is the total energy of the graphene-Pt(111) system, $E_{\rm pt}$ represents the energy of the five-layer Pt(111) slab, and $E_{\rm gr}$ represents the energy of the free-standing graphene sheet. Results are depicted in Figure 1 for both models.

To obtain an estimate of the stress needed to shear the graphene-Pt(111) interface, we also evaluate the ’interface energy surface’ for model (b) which is obtained by sliding the graphene sheet rigidly on the Pt surface. Results are shown in Figure 2 in conjunction with energy profiles taken along three typical paths. It can be seen that the interface energy variations are quite small. Accordingly, the interface shear stresses (IFSS) – derivatives of the energy profiles – required to slide the graphene over the Pt surface are comparatively low, of the order of 20 MPa or less. The same is expected to hold for the typical shear stresses required to slide a Pt cluster on a graphene sheet.

We use the LAMMPS simulation package [41] to perform molecular mechanics (MM) and molecular dynamics (MD) simulations of binary systems consisting of platinum nanoparticles which are anchored between two layers of graphene sheets. For describing such a system we need to define Pt-Pt, Pt-C and C-C interaction potentials, where the latter have to account for both short and long range interactions between carbon atoms within the graphene sheets and between both sheets.

In a first step, we examine different sets of interaction potentials to establish which set of potentials is most appropriate for the problem at hand. The first potential we examined is the reactive force field (ReaxFF) which has been developed by Sanz-Navarro et al. [42, 43]. ReaxFF is based on a bond-order model in conjunction with a charge-equilibration scheme[44]. This potential is in principle able to describe all the above mentioned interactions in the Pt-C system, however, it is optimized to correctly represent reactive processes rather than mechanical properties. The second potential we considered is a Brenner force field (BrennerFF) parametrized by Albe et al. [45] to describe short-range interatomic interactions in the Pt-C systems, with a Lennard-Jones potential added to mimic the weak van der Waals interaction between the carbon atoms of the two graphene sheets [46]. A further possibility is to combine potentials of different type that are separately optimized to describe the energetics and mechanical properties of covalently bonded carbon, and of platinum. A similar approach was used by [50, 51] for simulating Ni-graphene interactions: These authors use a standard AIREBO potential [46, 47] for describing covalent C-C interactions in conjunction with an EAM potential for Ni, and a Lennard-Jones potential for Ni-carbon interactions at the Ni-graphene interface. We emphasize that such addition of structurally dissimilar potentials is, in general, problematic. However, in the present context the approach may be feasible because the components of the system (metal, graphene) maintain their structural integrity, and the van der Waals type interaction between these components can be considered as additive to their internal (metallic, covalent) interactions. Accordingly we consider a Hybrid combination comprising the standard AIREBO potential [46, 47] for covalent C-C interactions and the EAM potential of Zhou and coworkers [48] for Pt-Pt interactions. Pt-C interactions at the interface are described using either a LJ potential taken from the literature [49]; this combination is denoted as HybridFF-LJ. Alternatively we consider a Morse potential of the form $E_{\rm M}=D(\exp[-2\alpha(r-r_{\rm 0})]-2\exp[-\alpha(r-r_{\rm 0})])$ where $r$ is the Pt-C distance, $r_{0}$ an equilibrium bond distance, $D$ is the well depth, and $\alpha$ controls the stiffness of the potential (the smaller $\alpha$ is, the smaller the attractive/repulsive forces). The combination of AIREBO (C-C) and EAM (Pt-Pt) with this Morse potential, whose parameters we determine by matching with our own DFT calculations, is denoted as HybridFF-M.

The performances of the different potentials or combinations of potentials are compared in view of two main criteria. First, we ask how well the potentials reproduce theoretical predictions for the mechanical properties of the system components. Second, we study how well they describe the adhesion between Pt and graphene. Regarding mechanical properties, the data in Table 1 show that all potentials correctly represent the lattice constant and cohesive energy of Pt, however, the ReaxFF produces isotropic elastic behavior for Pt, completely disregarding the significant cubic anisotropy of the material. The other potentials yield acceptable results regarding the elastic properties of Pt. Lattice properties of monolayer graphene as well as the interlayer distance $d_{\rm 0}$ and interlayer adhesion energy $\gamma_{\rm GG}$ of bilayer graphene are reported in Table 2.

Concerning the binding between Pt and graphene, we use our DFT results for the interaction between Pt(111) and graphene as a reference. We consider the same configuration as in the DFT simulations, i.e., we separate the graphene sheet rigidly from the Pt(111) surface, which allows us to directly compare the MD and DFT energies. First, we observe that the BrennerFF completely fails to reproduce the interfacial bonding between graphene and Pt(111) – the interaction reduces to a short-range repulsion and the binding energy is thus zero – which renders this potential unsuitable. The ReaxFF significantly over-estimates the bonding between Pt and graphene while it under-estimates the interfacial separation. Therefore, and also in view of its poor performance in modeling the elastic properties of Pt, this potential is also ruled out, which leaves us with the two hybrid combinations of potentials (HybridFF-LJ and HybridFF-M).

The two hybrid potentials offer the advantage that the parameters of the LJ or Morse potential, which describes Pt-C interactions at the interface between graphene sheets and Pt, can be fitted to correctly reproduce the interface separation and binding energy. At variance with previous studies [49], we find that the LJ potential has clear deficiencies as it is too stiff and thus significantly over-estimates the forces acting across the Pt-graphene interface (see Figure 3 (c,d)). The Morse potential is the only one that can be fitted very well to the DFT data because the well depth, well width and equilibrium distance parameters can be optimized to reproduce the quantum mechanical interaction energy curve from the density functional calculations. The resulting parameters for the Pt-C Morse potential are $D$ = 0.0071 eV, $r_{\rm 0}$ = 4.18 Å and $\alpha$ = 1.05 $\mathrm{\AA}^{\rm-1}$.

We therefore use in the following the HybridFF-M combination of potentials: The AIREBO potential [46, 47] for covalent C-C interactions, the EAM potential of Zhou and coworkers [48] for Pt-Pt, and the DFT parameterized Morse potential for C-Pt interactions. Van der Waals interactions between different graphene sheets are described by a standard LJ potential.

In nanoscale clusters, the configuration of Pt atoms differs from the planar slab used in our reference DFT calculation because many atoms have reduced coordination number and/or local environments of reduced symmetry. To verify that the HybridFF-M combination of potentials correctly describes the adhesion of small Pt clusters to graphene, we consider the extreme case of a very small $\mathrm{Pt}_{\rm 13}$ nanocluster adhering to a pristine graphene sheet. DFT calculations for this system were performed by Fampiou and Ramasubramanian [26] who report an adhesion energy of -0.84 eV for a minimum energy configuration of non icosahedral symmetry. Repeating the calculation with the same structure file but using the HybridFF-M potential for energy evaluation gives an adhesion energy of -0.94 eV, i.e., the deviation from the DFT result is only 11%. (We note that calculations using HybridFF-LJ, Tersoff and ReaxFF yield discrepancies of over 300 % of the DFT adhesion energy.) This gives us confidence that the HybridFF-M combination used in the following MM and MD simulations adequately captures the adhesion of larger Pt clusters on graphene.

For simulating the energy and geometrical properties of graphene-encapsulated Pt nanoparticles, we consider two types of particles: On the one hand, we use Pt particles constructed with cylindrical geometry (radius $R_{\rm p}$, particle height $2h$) to produce a prescribed aspect ratio of $q_{\rm p}=h/R_{\rm p}$. On the other hand, we consider annealed particles where an initially cylindrical particle is, before encapsulation, heated above the melting point of Pt and slowly cooled down to near 0 K using a Nosé-Hoover-thermostat in order to minimize the particle energy. Such annealed particles generally exhibit a spheroidal particle shape. In order to encapsulate the particles, we start from two graphene sheets distanced significantly beyond the cut-off radius in order to avoid van der Waals interaction between the layers. One platinum particle is positioned on the lower graphene sheet within the cut-off radius of the Pt-C interaction (Figure 4: step 1). To prevent edge effects, periodic boundary conditions are used in all directions. To form a system where the platinum particle is acting as a spacer between the two graphene sheets, the upper graphene sheet is rigidly moved right above the platinum particle (Figure 4: step 2). Then, the outer areas of the two sheets are brought into contact with each other with 3.4 $\mathrm{\AA}$ distance between them. A minimization step brings the system into mechanical equilibrium, where the platinum particle functions as a mechanical obstacle that prevents the two graphene sheets from adhering to each other over a ’detachment area’ surrounding the particle (Figure 4: step 3). In all simulations, we determine the total energy of the relaxed structure as well as the elastic energy of the graphene sheet and the energies of graphene-graphene, Pt-Pt, and graphene-Pt interaction. In addition, we determine the detachment area $A_{\rm d}$ as the area over which the spacing between the graphene layers exceeds the value for pure bilayer graphene by more than 1 Å, and for single encapsulated particles we define the detachment radius via $A_{\rm d}=:\pi R_{\rm d}^{2}$.

MM simulations of single embedded Pt nanoparticles of varying radius but constant aspect ratio indicate that the detachment radius surrounding a particle is proportional to the particle radius, i.e. $R_{\rm d}=4.2R_{\rm p}$ (Figure 5b). The energy of the graphene bilayer is raised by embedding a Pt nanoparticle. This energy increase (’GG excess energy’) is due to the energy of adhesion which must be provided in order to detach the two graphene sheets over the detachment area (detachment energy: $E_{\rm d}$), diminished by the adhesion between the two graphene sheets and the nanoparticle (attachment energy: $E_{\rm a}$), and due to the elastic distortion of the sheets (elastic energy: $E_{\rm el}$). All three energy contributions, which we evaluate from the MM/MD simulations, are found to increase in mutual proportion as the Pt particle radius increases. A simplified continuum model formulated in section 3 allows to understand this behavior.

In order to establish the dependence of the GG excess energy on the distance between Pt particles, we perform calculations on nanoparticle pairs as shown in Figure 6. To control the inter-particle separation distance $d$, in these MM simulations we keep the coordinates of the Pt atoms fixed during the final relaxation step. If the detachment areas surrounding two Pt nanoparticles overlap, this leads to a reduction in GG excess energy since the overall detachment area is reduced and also the elastic deformation of the graphene sheet decreases. The energy reduction becomes more pronounced as the distance between the two particles decreases, which leads to an attractive configurational force acting on the particles. Curves showing the GG excess energy of two-nanoparticle systems ($E^{\rm 2p}$) with different nanoparticle size vs the nanoparticle separation $d$ are shown in Figure 7. A data collapse can be obtained when the GG excess energy is scaled by its value for a single particle ($E^{\rm 1p}$), and the inter-particle separation is measured in units of the detachment radius $R_{\rm d}$ of the graphene sheet around that particle. The energy vs. distance curve can be well approximated by a linear dependency, $E^{\rm 2p}=E^{\rm 1p}(1+d/d_{\rm c})$ where $d\approx 2.4R_{\rm d}$ is the critical distance where the particles first interact with each other. Above that distance, we are dealing with two separate single particles (see Figure 7 inset).

To obtain semi-analytical relations between the size of the Pt nanoparticles encapsulated between the graphene sheets, the detachment between the two sheets, and the energy of the system, we use a simplified continuum model as illustrated in Figure 8, top. The particle is idealized as a cylindrical body with a radius of $R_{\rm p}$ and a height of $2h$, with its axis perpendicular to the adhering graphene sheets. Accordingly, a cylindrical coordinate system is used with its origin in the particle center. The $z$-axis runs along the particle axis, and the radial coordinate $r$ and angular coordinate $\theta$ denote the position parallel to the the plane $z=0$, which corresponds to the plane of adhesion between the two graphene sheets. The deformation of the graphene sheets is described by the displacements $u_{r},u_{\theta{}},u_{z}$. Pt-graphene binding length ($h_{\rm b}^{*}$) and graphene-graphene inter-layer distance ($2h_{\rm b}$) are assumed to be equal. A matching atomistic model is provided by a graphene sheet indented by a cylindrical indenter of radius $R_{\rm p}$. In this model all atoms within $R_{\rm p}$ from the center of the bilayer sheet are displaced rigidly in the respective up- and downward direction by a distance $h$ while the rest of the sheet is allowed to relax (Figure 8, bottom).

The graphene sheets are modeled as 2D elastic membranes with a bending stiffness that is negligible compared to their in-plane stiffness. The latter is characterized by the planar elastic modulus $E_{\rm 2D}$ and Poisson number $\nu$. The radius of the detached area surrounding the particles is denoted as $R_{\rm d}>R_{\rm p}$. Because of the mirror symmetry of the problem with respect to the plane $z=0$, it is sufficient to consider a single membrane deformed by the half particle, hence, the problem is similar to the problem of indentation of a 3D membrane of finite thickness as described by Begley [55], whose treatment we adapt to a 2D membrane sheet. We thus consider a membrane that is displaced, along the radius $r=R_{\rm p}$, by $u_{z}=h$ whereas for $r\geq R_{\rm d}$ the displacement is fixed at $u_{z}=0$. A solution of the elastic problem with these boundary conditions is sought for $R_{\rm p}\leq r\leq R_{\rm d}$. The approximate solution of this elastic boundary value problem is discussed in supplementary information. Here we only summarize the main results:

• •

  The displacement profile of the membrane is approximately given by

  $$u_{z}(r)=h\frac{q_{\rm d}^{a}-(r/R_{\rm p})^{a}}{q_{\rm d}^{a}-1}\quad,q_{\rm d}=\frac{R_{\rm d}}{R_{\rm p}},\quad a\approx 2/3$$

  (2)

• •

  The elastic energy of the system is given by

  $$E_{\rm el}=\pi R_{\rm p}^{2}\frac{2E_{2D}}{27}\frac{q_{\rm p}^{4}}{(q_{\rm d}^{2/3}-1)^{3}}$$

  (3)

The parameter $q_{\rm p}$ is the aspect ratio of the sandwiched particle. To determine the parameter $q_{\rm d}$ and hence the radius of the detached area surrounding the particle, we use a virtual work argument. The virtual work incurred upon expanding the detached area consists of the work released by the change in elastic energy and the work of adhesion that must be expended to detach the adhering graphene sheets. This detachment energy is simply given by ${\rm d}E_{\rm d}=2\pi(\gamma_{\rm GG}/2)R_{\rm d}{\rm d}R_{\rm d}$ where $\gamma_{\rm GG}$ is the adhesion energy per unit area of the graphene bilayer, and $\gamma_{\rm GG}/2$ the adhesion energy per sheet. We thus set ${\rm d}W={\rm d}E_{\rm el}+{\rm d}E_{\rm d}=0$ which leads to the equation

$$\frac{{\rm d}E_{\rm el}}{dR_{\rm d}}+\pi\gamma_{\rm GG}R_{\rm d}=0$$

(4)

From this virtual work balance we derive an implicit relation connecting $q_{\rm p}$ and $q_{\rm d}$:

$$\frac{q_{\rm p}}{q_{\rm d}^{1/3}(q_{\rm d}^{2/3}-1)}=\left(\frac{27\gamma_{\rm GG}}{4E_{2D}}\right)^{1/4}=0.283$$

(5)

where we used the adhesion energy $\gamma_{\rm GG}=0.288$ J/m² for bilayer graphene and a planar elastic modulus of $E_{2D}=300$ N/m. Taking the detachment ratio for a given particle aspect ratio from this relation, we can write the detachment energy as a function of particle radius

$$E_{\rm d}=\pi R_{\rm d}^{2}\gamma_{\rm GG}=\pi\gamma_{\rm GG}q_{\rm d}^{2}R_{\rm p}^{2}$$

(6)

From this relation and the detachment energies calculated by MM simulation in Figure 5, we find $q_{\rm d}\approx 4.0$ (red datapoints and red line in Figure 5a). This is consistent with direct determination of the detachment area as a function of nanoparticle radius, which gives $q_{\rm d}\approx 4.2$ (Figure 5b). By using Eq. (5) elastic energy of one graphene sheet can be obtained as

$$E_{\rm el}=\pi R_{\rm d}^{2}\frac{\gamma_{\rm GG}}{2}(1-q_{\rm d}^{-2/3})=\pi R_{\rm p}^{2}\gamma_{\rm GG}\frac{q_{\rm d}^{2}}{2}(1-q_{\rm d}^{-2/3})$$

(7)

which, with $q_{\rm d}=4.0$, gives us the black lines in Figure 5a and 5c. Finally, geometrical similarity requires that, also in case of non cylindrical particles, the effective area of contact between the graphene sheet and the Pt nanoparticle is proportional to the square of the particle radius. This gives us for the Pt-graphene adhesion energy

$$E_{\rm a}=\pi R_{\rm p}^{2}\gamma_{\rm PG}q_{\rm c}^{2}$$

(8)

By matching this with MM simulation data, we find the geometrical factor of $q_{\rm c}=0.87$ for annealed particles (blue data points and blue line in Figure 5a). Where $\gamma_{\rm PG}=0.476$ J/m² represents Pt-graphene adhesion energy per unit area. For cylindrical particles $q_{\rm c}=1$.

The elastic energy and the detachment radius are shown in Figure 9 as functions of the nanoparticle aspect ratio. By comparing the continuum results with the results of our simulations of annealed particles we can show that, for those particles, the elastic energy amounts to about $0.62E_{\rm d}$, hence $q_{\rm d}\approx 4$ according to the mechanical model. This is in good agreement with the observed proportionality $R_{\rm d}\approx 4.2R_{\rm p}$ found in the MM simulations in Figure 5c and Figure 9. These values correspond to an effective aspect ratio of $q_{\rm p}\approx 0.67$.

As noted in the MD section, if two particles get close and the respective detachment radii overlap, they experience an attractive mutual interaction. A semi-quantitative model can be formulated by approximating the conformation of the graphene sheets around a pair of interacting particles as two detached semicircles of radius $R_{d}$ connected by a detached ribbon of constant width $w\approx 2R_{\rm d}$ and area $2R_{\rm d}d$ (see Figure 6). The ribbon forms as soon as the particles approach below a critical distance $d_{\rm c}$, and the interaction energy then changes linearly with distance $d$.

Since the elastic distortion of the graphene sheet over the two detached semicircles is very similar to that around a single particle, we approximate the detachment energy in these regions by that of a single particle. On the other hand, across the detached ribbon the elastic energy is significantly reduced, and we therefore approximate the energy per unit area of the ribbon by $\gamma_{\rm GG}$. The total GG excess energy of the interacting two-particle system is then

$$E^{\rm 2p}(d,R_{\rm d})=1.6\pi\gamma_{\rm GG}R_{\rm d}^{2}+2\gamma_{\rm GG}dR_{\rm d}.$$

(9)

The detached ribbon forms as soon as this energy falls below the GG excess energy of two non interacting particles, $2E^{\rm 1p}=2\times 1.6\pi\gamma_{\rm GG}R_{\rm d}^{2}$. Based on this argument the critical separation follows via $1.6\pi R_{\rm d}=2d_{\rm c}$ as $d_{\rm c}\approx 2.5R_{\rm d}$ which is in good agreement with the MD data shown in Figure 7. The GG excess energy can then be written as

$$E^{\rm 2p}(d)\approx 2\gamma_{\rm GG}R_{\rm d}^{2}\frac{d_{\rm c}+d}{R_{\rm d}}$$

(10)

from which the interaction force of the two particles follows as

$$F^{\rm 2p}=\frac{\partial E^{\rm 2p}(d)}{\partial d}\approx 2\gamma_{\rm GG}R_{\rm d}$$

(11)

or, expressed in terms of $R_{\rm p}\approx R_{\rm d}/4$

$$F^{\rm 2p}=\frac{\partial E^{\rm 2p}(d)}{\partial d}\approx 8\gamma_{\rm GG}R_{\rm p}$$

(12)

The interaction force is, hence, approximately proportional to the Pt nanoparticle radius. Two nanoparticles are stable if this force is unable to overcome the shear resistance at the Pt-Graphene interface. The shear force required to slide a Pt nanoparticle inside the graphene bilayer can be estimated as $F_{\rm c}\approx 2\pi R_{\rm p}^{2}\tau_{\rm c}$ where $\tau_{\rm c}$ is the critical shear stress for interface sliding. Equating this to the inter-particle force we find a critical radius below which the particles will slide towards each other:

$$R_{\rm p,c}=\frac{4\gamma_{\rm GG}}{\pi\tau_{\rm c}}$$

(13)

With a typical critical shear stress of $\tau_{\rm c}=15$ MPa (Figure 3) we find a critical nanoparticle radius $R_{\rm p}\approx 24$nm. Particles below this radius, if closer than $d_{\rm c}\approx 10R_{\rm p}$, are likely to spontaneously aggregate.