Motion of water monomers reveals a kinetic barrier to ice nucleation on graphene

The interfacial behaviour of water remains a central question to fields as diverse as protein folding, friction and ice formation. While the properties of water at interfaces differ from those in the bulk, major gaps in our knowledge limit our understanding at the molecular level. Information concerning the microscopic motion of water comes mostly from computation and, on an atomic scale, is largely unexplored by experiment. Here, we provide a detailed insight into the behaviour of water monomers on a graphene surface. The motion displays remarkably strong signatures of cooperative behaviour due to repulsive forces between the monomers, enhancing the monomer lifetime ( ≈ 3 s at 125 K) in a free-gas phase that precedes the nucleation of ice islands and, in turn, provides the opportunity for our experiments to be performed. Our results give a molecular perspective on a kinetic barrier to ice nucleation, providing routes to understand and control the processes involved in ice formation.


Supplementary Methods
The growth and characterisation of the graphene layer on a Ni(111) surface has been published elsewhere 1 . In short, the nickel (111) single crystal used in the study was mounted onto a sample holder and can be heated radiatively by a filament on the back of the crystal, or cooled to 100 K using liquid nitrogen. Prior to the measurements, the Ni surface was cleaned by repeated cycles of Ar + sputtering and annealing at 870 K. A monolayer of graphene on Ni(111) was grown by dosing ethene (C2H4) while holding the crystal at 730 K for several hours.
The sample temperature was measured using a chromel-alumel thermocouple. While absolute temperatures can be determined to an accuracy of ±5 K, relative temperature values were determined with an accuracy of ±0.1 K, which was also confirmed by the reproducibility of the adsorption and dynamics measurements.
Water adsorption and desorption processes were studied during dosing with a precise water pressure control obtained by a motorised leak valve attached to the dosing supply. The leak valve itself was regulated by a feedback control system in order to maintain a constant pressure. Adsorption was measured at sample temperatures of 100, 110, 125, 130 and 150 K at a typical dosing pressure at the surface of (2 − 20) · 10 −9 mbar.
The helium reflectivity, at 100 K as shown in Figure 3 of the main article decreases continuously and remains low when the dose is stopped. Such a behaviour is typical of disordered structures forming, e.g. the growth of an amorphous layer. The formation of amorphous ice layers on surfaces, commonly referred to as amorphous solid water (ASW) has been observed since the 1960s 2 . For example, recent isothermal desorption measurements of water on highly oriented pyrolytic graphite (HOPG) at 100 K, showed a glass transition accompanied by a change in desorption rate and a growth of 3D water islands, rather than a wetting of the graphite surface 3 . At 110 K and 125 K, the helium reflectivity also decays during deposition, but reaches saturation. Based on the fact that exactly the same diffraction pattern is observed as from clean graphene, the deposition has been interpreted as forming separated islands of ice, leaving areas of bare graphene in between (see main text).
For sample temperatures above 120 K, when the applied pressure is reduced after dosing the helium reflectivity recovers very quickly. The system is thus in a dynamic equilibrium where small changes in the pressure immediately affect the coverage and hence the reflectivity. While with increasing overpressure the coverage increases, with increasing surface temperature the dynamic equilibrium is also reached faster. Within the available temperature range -where we could observe diffusion and where we are able to obtain a constant coverage by applying an overpressure -it was found that measurements at 125 K provided the best trade-off in order to clearly see dynamics and maintain a constant coverage. Dephasing rate measurements at 125 K are reported in the main text (Figure 1b), but further measurements at 130 K, indicate that within experimental uncertainties the measurements have the same variation with scattering momentum transfer, ∆K, and thus the same mechanism of motion at different temperatures. Finally, at even higher surface temperatures, i. e. around 150 K and above, negligible adsorption is observed, even when pressures of up to 1.5 · 10 −7 mbar were applied.

Supplementary Note 1: Details about coverage calibration and the monomer lifetime
All measurements were performed at the same coverage of 0.07 ML, which we define using a particular value of reflectivity, I/I0 = 0.25. The reflectivity is adjusted, at each temperature, by varying the overpressure of water vapour from the capillary doser (see Methods in the main text). We estimate the value of the coverage when I/I0 = 0.25 using the measured dynamical data ( Figure 5, main paper) and, in particular, the shape and position of the features arising from the pairwise inter-adsorbate forces. The method provides an absolute measure of coverage since the prominent minimum in the data (|∆K| ≈ 0.8Å −1 ) corresponds to a quasi-hexagonal overlayer with a spacing that defines the coverage. Varying the coverage in the kinetic Monte-Carlo calculations leads to the red curve in Figure 5 of the main paper, which corresponds to a coverage of (0.07 ± 0.02) ML (where 1 ML is defined as one water molecule per adsorption site). A second method of coverage calibration is less direct but it serves to confirm the internal consistency in our measurements and their interpretation. Here, we use the dosing rate to calculate the coverage, assuming a constant sticking coefficient of unity (see Refs. [4][5][6] , with monolayer coverage corresponding to n = 0.115 molecules/Å 2 , 7 which is close to the density of an ice I h overlayer 5 ). The rate of impingement is calculated from the measured partial-pressure of water in the chamber and the known enhancement generated by the microcapillary array 8 . Using that approach we can show that upon water adsorption at 100 K, the reflectivity follows a model for random adsorption (stick and sit, see ref. 9 ), with remarkable agreement over two-orders of magnitude.
Hence it confirms our assumption that the water monomers are static at 100 K. Also, the method allows us to estimate a scattering cross section for isolated water monomers of Σ = (120 ± 20)Å 2 , consistent with the coverage calibration derived from the dynamic measurements. Here, most of the quoted tolerance arises from uncertainties in the doser calibration and doser position relative to the sample. The determined cross section is in good agreement with other water cross sections in the literature on similar systems (e.g. Σ = 130Å 2 in Ref 10 ).
The fraction of the surface covered by ice islands cannot be estimated from the scattering data since the diffraction patterns with and without water adsorption are essentially identical (Figure 2d, main paper). It follows that the area of the surface covered by ice is extremely small. We can estimate the separation of ice islands from the r.m.s. distance travelled by an H2O monomer using known conditions for dosing and the measured diffusion constant. At 125 K the impinging flux of molecules was 2.5 · 10 17 m −2 s −1 and the corresponding rate of adsorption is 0.022 monolayer per second. Taking the coverage of water monomers, at equilibrium, to be 0.07 ML we obtain 3.4 s as the average lifetime of an H2O monomer, after adsorption and before being bound to an island. From the dynamics measurements in the main part of the manuscript we know that the hopping rate at 125 K is 1.5 · 10 10 s −1 . The mean jump length during diffusion is 3.3Å and it follows that a monomer travels about 70 µm before sticking to an island. The island separation must therefore be on a similar scale, which leads to the approximate scale-bar indicated in Figure 2c (main paper).

Supplementary Note 2: Desorption measurements
Several groups have conducted thermal desorption spectroscopy (TDS) measurements of water on the (0001) basal plane of graphite. Consistently, a single desorption peak was observed that corresponds to a desorption energy in the range of 0.4 to 0.5 eV which is close to the sublimation enthalpy of ice, 0.49 eV 5,6,11 . The observed desorption energy does not change with coverage, indicating the formation of separated islands on the graphene surface 11 .
On the surfaces of graphene/Ni(111) and of graphene/Ir(111), TDS spectra reveal pseudo-zeroth order desorption and desorption energies of (356±23) meV in the first case, and (585±31) meV in the latter case, respectively, were found 12 . Smith et al. report a desorption energy of about 490 meV for graphene/Pt(111), however they noted that the desorption of the water monolayer is complicated by de-wetting upon desorption and only multi-layer water films show zero-order desorption 7 .
We have also conducted thermal desorption spectroscopy while monitoring the m/z = 18 peak on a mass spectrometer and simultaneously measuring the helium reflectivity. A single desorption peak with a maximum at 163 K coincides with a rapid recovery of the specular signal. The Redhead equation can be applied, in order to estimate the desorption energy E d . Using ν = 9 · 10 −14 s −1 according to Ulbricht et al. 11 for the peak maximum at (163 ± 5) K at a heating rate β = 0.22 K · s −1 , we obtain a desorption energy of E d = (520 ± 20) meV.
Furthermore, we can use the recovery of the helium reflectivity to determine the desorption energy. We exposed the graphene surface to a water overpressure of 2 · 10 −8 mbar and waited until the system was in equilibrium, before turning off the exposure and monitoring the reflectivity recovery. From this we calculated the corresponding surface coverage as a function of time. The surface coverage first rises during exposure and then decays exponentially after exposure has been turned off. The initial desorption rate, which is identical to the exponential decay rate, exhibits an activated temperature dependence. The desorption energy can then be determined from the slope in an Arrhenius plot and gives a value of E d = (510 ± 10) meV, in good agreement with the conventional TDS method.
As mentioned in the main text, it suggest that water molecules tend to desorb from the surface of water islands into the vacuum rather than as individual molecules which are adsorbed on the bare graphene surface -with the latter being more likely to diffuse and bind to an island. These findings are further supported by the diffraction measurements stated in the main text and the mentioned de-wetting upon desorption as observed for water on graphene/Pt(111) 7 .

Supplementary Note 3
The density functional theory (DFT) approach has been applied a number of times to the adsorption of water on graphene. DFT calculations generally agree that the potential energy surface is rather flat and that the binding energy depends more on the orientation than on the position of the adsorbent. Most calculations predict a preferential water adsorption with the hydrogen atoms pointing downwards (see Figure 4b in the main paper and Supplementary Figure 1a). We obtained an adsorption energy, E ads , of about 250 meV for the global minimum, which is located at the centre of the graphene hexagonal unit cell. The value is in very good agreement with the 183 meV obtained by Li et al. 13 , while other DFT calculations give slightly smaller adsorption energies 12,14-17 . We also find the water molecule is adsorbed with the orientation usually predicted, i.e. that it has the two OH bonds pointing towards the surface, so that the plane of the molecule is perpendicular to the surface plane itself. A general agreement on an adsorption distance of about 3.3Å can be observed between our results and previous reports 13,15,18,19 .
In addition to calculations on bare graphene, we performed a set of vdW DFT calculations including a Ni substrate (modelled as a five-layers nickel slab with a ( √ 7 × √ 7) surface unit cell) to verify that the Ni substrate does not affect significantly the nature of water interaction with graphene. We found that the water to graphene distance and the orientation of the water monomers are very close to the calculations on the pristine graphene without the Ni substrate, although the absolute adsorption energy with the substrate is slightly different, ≈ 230 meV compared to 250 meV for suspended graphene. Including the substrate necessarily reduces the size of the surface unit cell and, by omitting the Ni, it is possible to increase the unit cell significantly (a 9-fold increase in the number of carbon atoms). The larger supercell allows us to model adsorption at significantly lower coverages, much close to those in the experiment. For these reasons, we used a larger graphene supercell without substrate to investigate the dynamics of water motion at low coverage, as presented in the manuscript.
We estimate the magnitude of the energy barrier to dimer formation by performing a series of DFT calculations for two water monomers (in the same orientation as the isolated monomer) with varying distance. The resulting barrier is about 90 meV, while the binding energy of the dimer is approximately 200 meV. Thus, the barrier to dimer formation is significant and we can conclude that, once a dimer forms, it will rarely dissociate. The dimer exhibits then a total adsorption energy of 595 − 696 meV (depending if the substrate is frozen or not). These results support also our observation of hysteresis during isobaric cooling/heating and they are consistent with previous calculations of of H2O clusters adsorbed on graphite where the association energy (including re-orientation, binding and adsorption) with the cluster is in the range of 450 − 500 meV, while the binding energy of a monomer to the graphene surface is much lower 15 In the measured ∆K-range of this study, the dephasing rate α is in the order of 10 ns −1 at 125 K. This temperature corresponds to a mean kinetic energy in the order of 10 meV. Since the adsorption energy of an H2O molecule in an ice cluster is predicted to be in the order of 500 meV 15 , while for the adsorption energy of a molecule on the graphene surface, values in the order of 100 − 200 meV have been calculated, one would expect to observe the diffusion of H2O on graphene, rather than on the surface of an ice cluster. Together with the adsorption and diffraction results this is another evidence that we are seeing the motion of single water molecules on graphene.

Supplementary Discussion
Water diffusion on other substrates Supplementary Table 1 compares our activation energies and diffusion constants (top row) with water measurements on other substrates in the literature.
In addition, the diffusion of water on graphene has been recently studied by means of molecular dynamics (MD) simulations. Tocci et al. predict a substantially lower macroscopic friction coefficient in comparison to adsorption on a hexagonal boron nitride surface 23   stant between 2 · 10 −7 m 2 /s and 8.6 · 10 −7 m 2 /s depending on the size of the droplet (at 298 K). Both values are way beyond the diffusion constant found in our experiments, yet they are considering the motion of water clusters and droplets at much higher temperatures (room temperature) rather than the diffusion of monomers. The diffusion coefficient for single water molecules on graphene has been estimated to be 6 · 10 −9 m 2 /s at a temperature of 100 K by Ma et al. 14 which is somewhat closer to the conditions in our experiments. Indeed their value is closer to our result but still one order of magnitude larger. However, all calculations mentioned above were performed on free standing graphene while our measurements are on graphene/Ni(111) where the motion of the ripples which gives rise to the ultra-fast diffusion 25 is suppressed 1 .
Compared to the diffusion in (bulk) amorphous ice on the other hand, where for the translational motion D0 is in the range of (0.5 − 5) · 10 −17 m 2 /s 26 , or the diffusion of ASW at the liquid/ice interface with D ≈ 10 −21 m 2 /s at 125 K 27 , the diffusion of water monomers on graphene with D = 4.1 · 10 −10 m 2 /s at 125 K is incomparably faster. We also note that the mobility on graphene is higher than translational diffusion of water in nanoconfinement 28 .

Adsorption on graphene and bulk graphite
The water dynamics may be affected by the metal substrate 7 or the presence of the bulk, in the case of graphite 29 . Most previous studies are performed in a different temperature or coverage regime so detailed comparisons are problematic. For example, on graphite 29 there is evidence that ice with a thickness of hundreds of monolayers coexists with regions of bare graphite, similar to our observations. In terms of the adsorption energies, vdW DFT predicts that for supported graphene, about 30% of the vdW interactions between the water and the substrate are transmitted through graphene 30 . Moreover, for graphene on metal substrates where a Moiré superstructures with a periodic height variation of the graphene layer forms, it has been reported that the regions closest to the metal substrate act as nucleation centres 31 . The desorption energies on the other hand are all within a similar energy region, both for water on metal supported graphene as well as for graphite 7,29,32 , however, desorption results are typically complicated by de-wetting upon desorption as further discussed in Supplementary Note 2: Desorption measurements Repulsive forces in adsorbate structure and dynamics On the one hand, inter-adsorbate repulsion occurs widely at surfaces and as mentioned in the main text these can limit the adsorbate density as well as defining the adsorbate structure. A typical example where it is commonly assumed that pairwise repulsion exists are (oriented) CO molecules 33,34 . On the other hand, processes such as thin-film growth and nucleation require attractive interactions and attractive forces are therefore often assumed to dominate the process. Our observations show that repulsion can control the surface motion, creating in the case of water molecules adsorbed on graphene a kinetic barrier to nucleation. However, the specific signature of repulsive interactions in surface diffusion remains a subject which is typically addressed via indirect approaches 35 while direct experimental evidence about repulsive interactions in adsorbate dynamics has only been addressed in very few cases such as for the diffusion of Na 36 . For example, the repulsion between CO molecules in the high-coverage regime based on information from structural organisation and adsorption sites 33,34 , could not be confirmed by surface diffusion measurements in the low-to medium coverage regime 37 .

Supplementary Note 4: Details about single particle and collective motion
The trajectories of the molecules versus time resulting from the KMC simulations can be used to calculate the intermediate scattering function (ISF) which is also obtained in the experiment. From the KMC simulation both the coherent and the incoherent ISF can be calculated. The subtle difference between the coherent and incoherent ISF is the averaging procedure. While the coherent ISF is obtained by averaging over all particles, the incoherent ISF is obtained by first calculating the ISF of a single particle followed by averaging over all particles. Details on how to obtain both the coherent and the incoherent ISF can be found elsewhere 38 .
The ISFs obtained from the simulation are then analysed in the same way as the experimental data: The ISF is fitted with a single exponential decay which allows to determine the dephasing rate α(∆K) from the simulation in analogy to the curve determined from the experiments. The trajectories from the KMC simulation can be used to calculate both the coherent and the incoherent ISF.
On the other hand, He spin-echo is a coherent scattering method, hence the measurements provide the coherent ISF. As shown for neutron scattering 39 as well as for X-ray photocorrelation spectroscopy 40 one can de-correlate the effect of adsorbate interactions, i.e. obtain the corresponding incoherent αinc(∆K) from the measured coherent α coh (∆K) 38 . Following the derivation of Sinha and Ross 41 , where the interaction forces are considered as a mean field, the scattering function Ss of the non-interacting system becomes: where c is the concentration of dynamic adsorbates, and the quasi-elastic broadening Γs follows the well established Chudley-Eliott lineshape (the analytical model as in the main text) 38 .
The only region where this approach does not apply is in the vicinity of the substrate diffraction peaks. Here the structure factor of the substrate becomes important while at the same time the uncertainty of the quasi-elastic amplitude becomes large. As a consequence the blue dots in Figure 1b of the main text show an offset for ∆K close to zero and around the diffraction peaks. Finally, here we note again that only the implementation of repulsive interactions in the KMC simulation can reproduce the peak-and-dip structure as evident in the experimental data. Supplementary Figure 2a shows both the dephasing rate α (top panel) and the corresponding static structure factor S(∆K) (lower panel) as extracted for the KMC simulations along the ΓKazimuth. We see that only for repulsive forces (B > 0, equation (4) in the main text) there appears a clear peak in S(∆K) at the same position, where α(∆K) shows a dip, as illustrated by the dash-dotted vertical line.
Supplementary Figure 2: a Comparison between the dephasing rate α (upper panel) and the corresponding static structure factor S(∆K) (lower panel) along the ΓK-azimuth, as extracted from the KMC simulations at the same conditions (including temperature and water coverage) as for the experiment. Only repulsive interactions (B > 0) give rise to a clear peak in S(∆K) at the same position where α(∆K) shows a dip as illustrated by vertical dash-dotted line. b Static structure factor S(∆K) along the ΓM-azimuth.
The repulsive forces between the adsorbates give rise to a deviation of the dephasing rate α(∆K) with respect to the analytical expression (equation (1) in the main text). Adsorbates repelling each other prefer a long-range quasi-hexagonal structure leading to a preferred, coverage dependent average distance (see Supplementary Note 1) between the adsorbates and reduced mobility on these length scales 36,42,43 . At the same time, when adsorbates approach each other their mobility increases compared to the nonrepelling case. The result is a peak at lower ∆K followed by a dip feature, termed as "de Gennes narrowing" as illustrated by the red line in the top panel of Supplementary Figure 2a. We see that the grey line from the KMC simulation without repulsive interactions follows the same sinusoidal curve as for the analyti-cal expression while the red line -illustrating the case for interadsorbate repulsion -exhibits a peak appearing at low ∆K values due to the increased mobility at certain length scales, followed by a dip (vertical dash-dotted line) occuring at the length scale of the quasi-hexagonal arrangement. The location of this dip corresponds to a peak in the static structure factor 44 , as seen in the lower panel of Supplementary Figure 2a and Supplementary Figure 2b. The peak in S(∆K) appearing at ∆K ≈ 0.8Å −1 , designated by the green vertical line, corresponds to an intermolecular distance of about 8Å in real space. The preferred average distance during the diffusion of water monomers on graphene in the low-coverage regime is thus much larger compared to any spacing observed in X-ray diffraction from amorphous or polycrystalline bulk ice. In the latter case the first peak in S(∆K) appears at 1.7Å −1 and hence at much smaller intermolecular distances in real space 26 , as expected. However, it is important to note that the dynamics measurements in the low coverage regime do not provide us with information about the actual structure of the adsorbate layer. Structural information can instead be taken from diffraction measurements such as Figure 2d in the main text -which in the current scenario are conclusive about the formation of a disorderd film in the high coverage regime. Moreover, the mentioned static structure factor in the case of X-ray diffraction is an isotropically averaged S(∆K) due to being from a polycrystalline/amorphous sample 26 meaning that the radial distribution function can be obtained in a straightforward Fourier transform. S(∆K) for 2D diffusion in a HeSE experiment, on the other hand, will be direction dependent. It is determined by the azimuthal orientation of the crystal -given by ∆K -and provides a measure of the projection along that direction 42 as seen for the ΓK and ΓM azimuth in Supplementary  Figure 2. Finally, we note that the current observation of long-range repulsive interactions does not exclude the possibility of short-range attractive interactions and it is the implementation in the KMC that reproduces the feature in the experimental data. Short-range interactions may rather occur within a length scale that corresponds to intra-cell diffusion 45 while the discrete grid in terms of the KMC simulations allows just for interactions at the the intercell diffusion length-scale to be taken care of.