In this work we propose a general strategy to calculate accurate He–surface interaction potentials. It extends the dispersionless density functional approach recently developed by Pernal et al [Phys. Rev. Lett. 103, 263201 (2009)] to adsorbate-surface interactions by including periodic boundary conditions. We also introduce a scheme to parametrize the dispersion interaction by calculating two- and three-body dispersion terms at coupled cluster singles and doubles and perturbative triples (CCSD(T)) level via the method of increments [H. Stoll, J. Chem. Phys. 97, 8449 (1992)]. The performance of the composite approach is tested on ^{4}He/graphene by determining the energies of the low-lying selective adsorption states, finding an excellent agreement with the best available theoretical data. Second, the capability of the approach to describe dispersionless correlation effects realistically is used to extract dispersion effects in time-dependent density functional simulations on the collision of ^{4}He droplets with a single graphene sheet. It is found that dispersion effects play a key role in the fast spreading of the ^{4}He nanodroplet, the evaporation-like process of helium atoms, and the formation of solid-like helium structures. These characteristics are expected to be quite general and highly relevant to explain experimental measurements with the newly developed helium droplet mediated deposition technique.

The ultra-low temperature helium droplet mediated synthesis and deposition technique of metal nanoparticles (NPs) on solid surfaces, originally developed by Vilesov's group,^{1–3} has attracted much interest over the last two years.^{4–7} This is due to both the exciting fundamental physics revealed via the technique, including earlier traces of quantum vorticity in superfluid ^{4}He droplets^{3} and the resulting potential applications in nanoscience and nanotechnology.^{5} Direct experimental evidences of quantum vorticity have just been reported through X-ray diffraction images of doped ^{4}He droplets by Gómez *et al.*^{8} It can be exploited to induce the formation of ultrathin wires of metal NPs,^{3,7} whereas the experimental set-up can be tailored to form metal core-shell morphologies^{5} and to produce NPs films also beyond the submonolayer regime.^{6} The microscopic understanding via first-principles simulations of the metal NPs grown inside the carrier droplet, subsequent deposition, and their interaction with either already deposited NPs or incoming ^{4}He droplets is relevant due to its potential to provide basic information and simplifications that can be then transferred for control purposes. These experiments frequently use amorphous carbon-based surfaces as the standard substrate to image the deposited metal NPs via transmission electron microscopy (TEM). As the simplest carbon-based surface, the work presented here focuses on a single graphene sheet (see Fig. 1). The relevance of graphene in nanomaterials research^{9} is very well-known due to its unique properties such as the high charge-carrier mobility, optical transparency, elasticity, and thermal conductivity. Very recent studies^{10,11} have considered the graphene application as a surface-enhanced Raman spectroscopy substrate upon deposition of plasmonic NPs. Owing to the higher conductivity and large surface area, it has been suggested that graphene can be used to improve the catalytic properties of supported metal NPs,^{12} which could be further enhanced through the one-dimensional quantum confinement effects induced by the helium droplets.

The collision of the ^{4}He droplet with the TiO_{2}(110) surface has recently been addressed using nuclear time-dependent density functional theory (TDDFT) and classical trajectory calculations.^{13} In contrast with the classical picture predicting the splashing of the helium droplet upon impact, TDDFT simulations revealed the droplet spreading and demonstrated the key role of quantum effects. Besides finite-temperature surface effects, the mobility of the deposited NPs and aggregation are expected to be mostly determined by the spreading dynamics of either the carrier or the incoming helium droplets. This process is naturally strongly influenced by the specific He–surface potential interaction so that its accurate description is necessary.

Cost-efficient determinations of adsorbate-surface interactions would presumable use periodic electronic structure codes and methods based on density functional theory with inclusion of dispersion corrections. During the last few years, numerous van der Waals (vdW)-corrected DFT methods (referred to as DFT+D) have been developed and implemented into standard periodic codes, with numerous applications to adsorbates on different surfaces (for a recent review see Ref. 14). For example, dispersion-corrected DFT energies can be obtained by adding interatomic vdW C_{6}/R^{6} terms using the parameterization DFT-D2 of Grimme.^{15} All the DFT+D approaches assume that the chosen DFT method accounts for dispersionless correlation effects realistically. The Perdew-Burke-Ernzerhof (PBE)^{16} functional has been commonly used in He-surface interaction problems. A localized molecular orbital decomposition analysis^{17} of the He–TiO_{2}(110) interaction has recently illustrated how the overestimation/underestimation of attractive/repulsive PBE-based interaction energy components can compensate for the shortcomings in the dispersion,^{18} disabling the further addition of vdW corrections. This study has also demonstrated the success of the dispersionless density functional approach (dlDF) by Pernal *et al.*^{19} This is a hybrid generalized gradient approximation (meta-GGA) functional which differs from the M05-2X functional^{20} in the number and values of DFT parameters which were optimized to reproduce benchmark dispersionless interaction energies of weakly bound dimers.^{19} The He–graphene problem is even more challenging due to its highly delocalized electronic structure. Table I collects the interaction energies with the He atom at the most stable adsorption site (i.e., the hollow site). To perform the periodic calculations, we used the computational setup reported for a study of water/graphene^{21} and the augmented polarized correlation-consistent triple-ζ basis of Dunning, Jr. and collaborators^{22} (aug-cc-pVTZ) for He. The interaction energies were counterpoise-corrected and the basis set quality was assessed by reproducing the same energies in plane-wave calculations. The maximum deviation was found to be 2 meV, corresponding to 5% of the total interaction energy. To get accurate dispersionless interaction energies, we have implemented the dlDF functional in a development version of the CRYSTAL14 code.^{23,24} As previously shown for He–TiO_{2}(110) at CCSD(T) level,^{18} the short-range *intramonomer* correlation contributions to the dispersionless interaction energy are expected to be repulsive owing to the correlation space truncation for each monomer (i.e., the adsorbate and the surface).^{25} By comparing dlDF and PBE interaction energies with the Hartree-Fock (HF) counterparts at the shortest He–surface distance, we notice that the dlDF correlation is repulsive as opposed to PBE. It follows that the PBE-D2 approach overestimates the attractive interaction. The dlDF interaction energies with inclusion of periodic conditions differ very little to those obtained using the dlDF MOLPRO implementation.^{18,26} However, the former are weakly attractive at medium range (around 4.0 Å). To include the vdW correction on top of the dispersionless interaction energies, Szalewicz and collaborators developed an effective pairwise interatomic functional^{19,27,28} named D_{as} with parameters optimized to reproduce symmetry-adapted perturbation theory-DFT (SAPT(DFT))^{29,30} dispersion energies on a training set data. Applying the dlDF+D_{as} *ansatz* to He–graphene, the interaction energy is given by

The sum in the second term (the D_{as} function) runs over as many graphene C atoms as necessary to get convergence and *f*_{n} are the damping functions of Tang and Toennies.^{31} The excellent performance of the dlDF+D_{as} approach to describe the He–TiO_{2}(110) interaction was first tested in Ref. 18 with the D_{as} functional parametrized using benchmark dispersion energies on a small He-cluster model. We propose here a complementary approach in parameterizing the D_{as} functional (referred to as incremental |${\rm D}_{as}^*$|$Das*$). It applies the method of increments originally developed by Stoll^{32} on surface cluster models of increasing size. Applying this method,^{18,21,25,32,33} the *intermonomer* correlation contribution to the He–surface interaction energy is written as a cumulant expansion in terms of localized orbital groups (LOGs) from the adsorbate (*A*) and the surface (*i*, *j*) which define *n*-body increments η (*n* denote the number of interacting LOGs within each increment),

The correlation term associated with the dispersion is mainly determined by the two-body increments η_{Ai}. These terms are defined as the non-additive part of the correlation energy ε associated with the simultaneous correlation of electrons from two LOGs centered at the He atom (*A*) and the surface (*i*), respectively, η_{Ai} = ε_{Ai} − ε_{A} − ε_{i}. Modeling the graphene surface with coronene and using a Foster-Boys localization procedure,^{34} the six equivalent LOGs are formed from (see also Fig. 1) (1) the 1*s* He orbital; (2) the six σ C−C bond orbitals located at the central benzene-like ring; (3) the twelve σ orbitals of the outer rings; (4) and (5) the six σ + π C=C bond orbitals of the outer rings oriented either radially or tangentially from the midpoint; (6) and the twelve dangling C−H σ bond orbitals. One-body increments characterize the *intramonomer* correlation contributions to the interaction energy. These contributions agreed very well with those obtained as the difference between dlDF and HF interaction energies on He/TiO_{2}(110).^{18} We have assumed that the dlDF approach is capable of providing accurate *intramonomer* correlation contributions also on He/graphene and, therefore, the corresponding one-body increments have not been calculated. The *C*_{n} coefficients and the damping D_{as} parameters β are obtained through the global fitting of the total *intermonomer* correlation contribution including all the two-body increments and their effective reduction when the most relevant three-body counterparts are accounted for. Since all the incremental contributions are summed, the rings inequivalency produced by the localization procedure has no effects.

He on top of the graphene . | Energy (meV) . | |||
---|---|---|---|---|

“hollow” site . | 2.4 Å . | 3.2 Å . | 4.0 Å . | 6.0 Å . |

|${\rm E}_{{\rm int}}^{{\rm PBE}}$|$E int PBE $ | 38.8 | −5.01 | −3.23 | −0.02 |

|${\rm E}_{{\rm int}}^{{\rm PBE-D2}}$|$E int PBE \u2212D2$ | −9.03 | −21.4 | −9.95 | −1.27 |

|${\rm E}_{{\rm int}}^{{\rm HF}}$|$E int HF $ | 109.0 | 7.94 | 0.48 | 0.0 |

|${\rm E}_{{\rm int}}^{{\rm dlDF}}$|$E int dlDF $ | 133.7 | 9.78 | −0.16 | −0.02 |

|${\rm E}_{{\rm int}}^{{\rm dlDF}}$|$E int dlDF $(coronene) | 135.9 | 8.99 | 0.06 | 0.0 |

Two-body increments | ||||

|$\sum _{i}\eta _{A\sigma _{i}}$|$\u2211i\eta A\sigma i$(He/C−C inner ring) | −23.6 | −4.98 | −1.36 | −0.12 |

|$\sum _{i}\eta _{A\sigma _{i}}$|$\u2211i\eta A\sigma i$(He/C−C outer ring) | −4.21 | −1.82 | −0.82 | −0.13 |

|$\sum _{i}\eta _{A(\sigma +\pi )_{i}}$|$\u2211i\eta A(\sigma +\pi )i$(He/C=C radial) | −59.5 | −16.9 | −5.37 | −0.57 |

|$\sum _{i}\eta _{A(\sigma +\pi )_{i}}$|$\u2211i\eta A(\sigma +\pi )i$(He/C=C tangential) | −10.8 | −3.99 | −1.69 | −0.27 |

|$\sum _{i}\eta _{A\sigma _{i}}$|$\u2211i\eta A\sigma i$(He/C−H) | −1.63 | −0.88 | −0.47 | −0.09 |

∑_{i}η_{Ai}(total) | −99.7 | −28.6 | −9.71 | −1.18 |

Three-body increments | ||||

|$\sum _{i < j}\eta _{A(\sigma +\pi )_{i}(\sigma +\pi )_{j}}$|$\u2211i<j\eta A(\sigma +\pi )i(\sigma +\pi )j$ | 3.84 | 1.15 | 0.37 | 0.01 |

|${\rm E}_{{\rm int}}^{{\rm tot}}$|$E int tot $ | 37.8 | −17.7 | −9.50 | −1.19 |

He on top of the graphene . | Energy (meV) . | |||
---|---|---|---|---|

“hollow” site . | 2.4 Å . | 3.2 Å . | 4.0 Å . | 6.0 Å . |

|${\rm E}_{{\rm int}}^{{\rm PBE}}$|$E int PBE $ | 38.8 | −5.01 | −3.23 | −0.02 |

|${\rm E}_{{\rm int}}^{{\rm PBE-D2}}$|$E int PBE \u2212D2$ | −9.03 | −21.4 | −9.95 | −1.27 |

|${\rm E}_{{\rm int}}^{{\rm HF}}$|$E int HF $ | 109.0 | 7.94 | 0.48 | 0.0 |

|${\rm E}_{{\rm int}}^{{\rm dlDF}}$|$E int dlDF $ | 133.7 | 9.78 | −0.16 | −0.02 |

|${\rm E}_{{\rm int}}^{{\rm dlDF}}$|$E int dlDF $(coronene) | 135.9 | 8.99 | 0.06 | 0.0 |

Two-body increments | ||||

|$\sum _{i}\eta _{A\sigma _{i}}$|$\u2211i\eta A\sigma i$(He/C−C inner ring) | −23.6 | −4.98 | −1.36 | −0.12 |

|$\sum _{i}\eta _{A\sigma _{i}}$|$\u2211i\eta A\sigma i$(He/C−C outer ring) | −4.21 | −1.82 | −0.82 | −0.13 |

|$\sum _{i}\eta _{A(\sigma +\pi )_{i}}$|$\u2211i\eta A(\sigma +\pi )i$(He/C=C radial) | −59.5 | −16.9 | −5.37 | −0.57 |

|$\sum _{i}\eta _{A(\sigma +\pi )_{i}}$|$\u2211i\eta A(\sigma +\pi )i$(He/C=C tangential) | −10.8 | −3.99 | −1.69 | −0.27 |

|$\sum _{i}\eta _{A\sigma _{i}}$|$\u2211i\eta A\sigma i$(He/C−H) | −1.63 | −0.88 | −0.47 | −0.09 |

∑_{i}η_{Ai}(total) | −99.7 | −28.6 | −9.71 | −1.18 |

Three-body increments | ||||

|$\sum _{i < j}\eta _{A(\sigma +\pi )_{i}(\sigma +\pi )_{j}}$|$\u2211i<j\eta A(\sigma +\pi )i(\sigma +\pi )j$ | 3.84 | 1.15 | 0.37 | 0.01 |

|${\rm E}_{{\rm int}}^{{\rm tot}}$|$E int tot $ | 37.8 | −17.7 | −9.50 | −1.19 |

The correlated individual increments η are listed in Table I. They were calculated at CCSD(T) level of theory with the MOLPRO code,^{26} using cc-pVQZ and aug-cc-pVQZ basis sets for C_{24}H_{12} and He, respectively. The vdW contribution to the interaction is clearly dominated by the attractive two-body η_{Ai} increments which involves the π localized orbitals closest to the He atom. Notice also that their values decay slowly as the He–graphene distance increases. The repulsive three-body terms including the He orbital and two nearest neighbors π orbital groups represent a minor fraction (below 6%) and are short-ranged. Other three-body terms (not shown) contribute with less than 2%.

As a stringent test for the accuracy of the composite approach, we have calculated the nuclear bound-state energies. The He–surface interaction was averaged by considering the three adsorption sites shown in Fig. 1 (see the supplementary material^{36}). The resulting vibrational energies were compared with the best available theoretical data^{35} (see Table II and the supplementary material^{36}). We can clearly notice from Table II the improvement that the |${\rm D}_{as}^{*}$|$Das*$ scheme brings over the original D_{as} parametrization.^{28} The results become systematically better upon increasing the cluster model used so that for coronene |${\rm dlDF+D}^{*}_{as}$|$ dlDF +Das*$ bound-state energies differ from the theoretical reference values by less than 0.3 meV. Notice that the negligible role of the coronene dangling bonds on the new |${\rm D}^{*}_{as}$|$Das*$ parametrization marks the convergence with respect to the cluster model size. The nuclear energies agreed also very well to those obtained by using the SAPT(DFT) method for the parametrization (i.e., the periodic dlDF + D_{as} scheme) as well as the experimental data on graphite.^{36}

. | ε_{0}
. | ε_{1}
. | ε_{2}
. | ε_{3}
. | ε_{4}
. | ε_{5}
. |
---|---|---|---|---|---|---|

D_{as}^{a} | −14.69 | −8.28 | −3.99 | −1.52 | −0.40 | −0.05 |

|${\rm D}^{*}_{as}$|$Das*$(C_{6}H_{6}) | −14.00 | −7.85 | −3.77 | −1.45 | −0.39 | −0.05 |

|${\rm D}^{*}_{as}$|$Das*$(C_{24}H_{12}) | −12.88 | −6.95 | −3.16 | −1.12 | −0.27 | −0.03 |

|${\rm D}^{*}_{as}$|$Das*$(C_{24}H_{12})^{b} | −12.92 | −6.96 | −3.16 | −1.11 | −0.26 | −0.03 |

Theory^{c} | −12.63 | −6.68 | −2.93 | −1.03 | −0.24 | −0.04 |

. | ε_{0}
. | ε_{1}
. | ε_{2}
. | ε_{3}
. | ε_{4}
. | ε_{5}
. |
---|---|---|---|---|---|---|

D_{as}^{a} | −14.69 | −8.28 | −3.99 | −1.52 | −0.40 | −0.05 |

|${\rm D}^{*}_{as}$|$Das*$(C_{6}H_{6}) | −14.00 | −7.85 | −3.77 | −1.45 | −0.39 | −0.05 |

|${\rm D}^{*}_{as}$|$Das*$(C_{24}H_{12}) | −12.88 | −6.95 | −3.16 | −1.12 | −0.27 | −0.03 |

|${\rm D}^{*}_{as}$|$Das*$(C_{24}H_{12})^{b} | −12.92 | −6.96 | −3.16 | −1.11 | −0.26 | −0.03 |

Theory^{c} | −12.63 | −6.68 | −2.93 | −1.03 | −0.24 | −0.04 |

By construction, the composite approach is designed to describe accurately not only the total but also the dispersionless He–surface interaction energy, allowing to extract truly dispersion effects in the nuclear dynamics. For this purpose, we have applied the TDDFT method using both the dlDF and dlDF+incremental |${\rm D}_{as}^{*}$|$Das*$ laterally averaged He–graphene potentials which have well-depths of 0.5 meV (at 4.3 Å) and 16.4 meV (at 3.3 Å). The details of the method when applied to the (zero-temperature) collision dynamics of ^{4}He droplets on surfaces can be found in Ref. 13. Briefly, the droplet is described by a complex effective wavefunction Ψ(**r**, *t*) such that ρ(**r**, *t*) = |Ψ(**r**, *t*)|^{2}, with the number of ^{4}He atoms fixed to 300 in this work. The helium wave-packet follows the 3D time-dependent equation,

where |${\cal E}_{{\rm He}}[\rho ]$|$E He [\rho ]$ is a modification^{37} of the Orsay-Trento density functional^{38} for the helium nuclei which is capable of describing very structured helium configurations.^{39} The term *V*_{ext}(**r**) denotes the He–graphene potential and Λ(**r**) is a damping function avoiding the wave-packet reflection on the box boundaries. The initial wave-function is given by

where ρ_{0}(**r**) = |Ψ_{0}(**r**)|^{2} is obtained by minimizing |${\cal E}_{{\rm He}}[\rho ]$|$E He [\rho ]$ without including the He-graphene interaction. Following the experimental setup,^{2} the boost |${\bf k}_{0}=-1.26\;\hat{e}_{z}$|$k0=\u22121.26e\u0302z$ Å^{−1} provided to the droplet with a collective initial velocity towards the surface plane of 200 m/s.

Figure 2 shows the evolution of the wave-packet during the first 27 ps (multimedia view). The dynamical simulations start with the droplet mass center at *z*_{0} = 27.4 Å from the surface. During the first picoseconds, the droplet is accelerated towards the substrate. The inclusion of the dispersion in the He–graphene potential causes an earlier and much more pronounced compression of the droplet. After approximately 9 ps, the droplet reaches the graphene surface and pressure density waves start to propagate backwards. The spreading of the droplet on the graphene sheet starts at about 11 ps and it is markedly influenced by the spontaneous symmetry breaking and the appearance of high density fluctuations like the snowball observed in ^{4}He droplets with a high attractive impurity inside^{39} (Multimedia view). Simultaneously, the density waves propagating backwards reach the droplet surface, starting the evaporation-like process which is recognized in the helium densities at *t* = 18 ps in Fig. 2. Although the droplet evaporation is more evident in the normal direction, it shows up along different angles to the surface plane. With the dispersion-accounting potential, the spreading along the graphene sheet and the evaporation-like process continue till the end of the simulation at *t* = 27 ps, when the droplet has reached the box boundaries and the calculation has been necessarily finished.

Owing to the extremely weak He–graphene attractive interaction without including the dispersion (even weaker than the He–He interaction), the spreading is not complete and a thick helium layer on the graphene sheet remains standing. This is reminiscent of the partial wetting of ^{4}He droplets when spread on the Cs surface.^{40} On the contrary, a liquid ^{4}He film on graphene has been predicted to be metastable with respect to a commensurate solid.^{41} The effects of the He–graphene dispersion interaction can be discerned already in the density contours at *t* =9 ps presented in Fig. 2 as the marked bending of density waves towards the middle of the surface plane. They are even more evident at *t* =18 ps (see Fig. 2) when solid-like helium spots with very high densities appear. These solid-like helium structures, however, are not stable. They are annihilated at impact with the graphene sheet, releasing energy into the droplet and contributing to the further evaporation of helium atoms.

In concluding, our dynamical calculations have shown the key role played by dispersion effects in the fast spreading of a ^{4}He droplet on a graphene sheet, including the evaporation of helium atoms along different angular orientation from the surface plane and the appearance of solid-like helium structures. The fast flow of helium density along the surface plane is expected to promote both the mobility of deposited metal NPs laterally to the surface plane and their subsequent aggregation. Overall, the spreading process is more complex and richer than previously shown using the PBE approach for the He–TiO_{2}(110) interaction^{13,42} (even including a long-range dispersion term). It highlights the importance of using accurate He-surface potentials to capture the microscopic details of the dynamical process. The excellent performance and simplicity of the periodic dlDF+incremental |${\rm D}^{*}_{as}$|$Das*$ and periodic dlDF + D_{as} approaches make them suitable in first-principles simulations of helium-mediated deposition processes.

We thank Andrey Vilesov and Massimiliano Bartolomei for very useful discussions and suggestions. This work has been supported by Grants Nos. CCG08-CSIC/ESP-3680 from CSIC-CM, FIS2011-29596-C02-01 and FIS2011-28617-C02-01 from DGI, Spain (FEDER), and 2009SGR1289 from Generatitat de Catalunya. The aid of COST Action CM1002 “COnvergent Distributed Environment for Computational Spectroscopy (CODECS)” is also acknowledged. The Cesga Super-Computer Center (Galicia) and the Centro Téchnico de Informática (CTI, CSIC) are acknowledged for providing computer time.

## REFERENCES

*ab initio*programs, 2012, see http://www.molpro.net.

^{4}He/graphene and

^{4}He/graphite as well as multimedia view.