Machine learning exploration of binding energy distributions of H2O at astrochemically relevant dust grain surfaces

Binding energies (BEs) of adsorbates on interstellar dust grains critically control adsorption, desorption, diffusion, and surface reactivity, and therefore strongly influence astrochemical models of star- and planet-forming regions. While recent computational studies increasingly report full distributions of BEs rather than single representative values, these distributions are typically derived for either bare grain surfaces or thick water-ice mantles. In this work, we bridge these regimes by systematically investigating the BE distributions of water on partially and fully ice-covered dust grain surfaces. We employ machine-learning interatomic potentials (MLIPs) based on graph neural networks to model water adsorption on graphene and on the Mg-terminated (010) surface of forsterite, representing carbonaceous and silicate grains, respectively. The models enable extensive sampling of adsorption sites on water clusters, monolayers, and bilayers generated under both crystalline (thermally processed) and amorphous (low-temperature) growth conditions. At submonolayer coverage, the chemical nature of the underlying grain strongly affects both ice morphology and binding energies, with Mg-O interactions on silicate surfaces producing particularly deep binding sites. From monolayer coverage onward, adsorption on both substrates is dominated by hydrogen bonding within the ice, reducing the influence of the grain material. Across all coverages, amorphous ice structures systematically shift the BE distributions toward stronger binding compared to crystalline ice, introducing highly stable defect and pocket sites. These results demonstrate that BE distributions in the submonolayer to few-layer ice regime are broad and highly surface dependent, and they provide physically motivated input for next-generation astrochemical models incorporating surface heterogeneity.

💡 Research Summary

This paper addresses a central challenge in astrochemical modeling: the accurate representation of binding energies (BEs) for water molecules adsorbed on interstellar dust grains. Traditional astrochemical networks typically employ a single BE value per surface, an approach that neglects the intrinsic heterogeneity of realistic grain–ice systems. The authors bridge the gap between bare grain surfaces and thick water‑ice mantles by systematically exploring BE distributions from the sub‑monolayer regime up to a few ice layers, using both crystalline (thermally processed) and amorphous (low‑temperature) ice morphologies.

Methodology
Two prototypical grain materials are selected: (i) graphene, representing carbonaceous grains, and (ii) the Mg‑terminated (010) facet of forsterite (Mg₂SiO₄), representing silicate grains. High‑level density‑functional theory (DFT) calculations (PBE‑D3/D4) generate reference energies and forces for a diverse set of water/graphene and water/forsterite configurations. These data are used to train message‑passing neural network potentials (PaiNN), a class of graph‑neural‑network interatomic potentials that can predict energies and forces with near‑DFT accuracy while being orders of magnitude faster. Training incorporates both global structure‑optimization (GOFEE) searches (up to 12 water molecules on graphene, 65 on forsterite) and ab‑initio molecular dynamics (AIMD) at various temperatures, ensuring coverage of stable minima, high‑energy “defect” structures, and non‑equilibrium configurations. The final datasets contain roughly 5 700 structures for each substrate, after careful curation to remove unphysical events (e.g., Mg extraction from the surface).

Surface Generation
With the trained MLIPs, realistic ice‑covered grain models are generated. Two growth protocols are employed: (a) low‑temperature (10 K) NVT molecular dynamics to produce amorphous ice, and (b) global optimization to yield crystalline ice. For each substrate three coverage regimes are constructed: (1) clusters (partial coverage, exposing most of the underlying grain), (2) monolayers (nearly full coverage), and (3) bilayers (almost complete coverage). Water molecules are placed both on top of and beneath the ice, allowing the study of “sandwiched” configurations that may occur on porous grains.

Binding‑Energy Analysis
Binding energies are computed for every sampled adsorption site. The key findings are:

1. Sub‑monolayer regime: The chemical nature of the grain dominates the BE distribution. On the Mg‑terminated forsterite surface, strong Mg–O interactions create deep potential wells, yielding average BEs around –1.2 eV (≈ 30 meV deeper than literature values). Graphene, governed by weaker van‑der‑Waals forces, shows a broader, shallower distribution centered near –0.8 eV.

2. From monolayer onward: Hydrogen‑bond networks within the ice become the primary binding mechanism. Consequently, the substrate influence diminishes; both graphene and forsterite converge to average BEs of ≈ –0.6 eV, reflecting the intrinsic strength of water–water hydrogen bonds.

3. Crystalline vs. amorphous ice: Amorphous ice consistently shifts the BE distribution toward stronger binding. Defect and pocket sites in the disordered network generate high‑energy tails, with some sites reaching –1.5 eV—substantially deeper than any crystalline counterpart. This effect is most pronounced in the sub‑monolayer and cluster regimes, where the underlying grain still contributes to the local environment.

4. Distribution shape: The BE histograms are markedly non‑Gaussian, exhibiting asymmetric tails and multiple peaks corresponding to distinct binding motifs (e.g., Mg‑O coordination, surface dangling OH, buried water molecules). The breadth of the distributions (standard deviations up to 0.3 eV) underscores the inadequacy of a single‑value description.

Implications for Astrochemical Modeling
The authors discuss how these detailed BE distributions can be incorporated into kinetic Monte Carlo (KMC) simulations, where each lattice site can be assigned its own BE, naturally reproducing heterogeneous diffusion and desorption rates. For rate‑equation models, the paper suggests two pathways: (i) use the full distribution to compute temperature‑dependent effective rates via statistical averaging, or (ii) adopt a multi‑well approach where the distribution is discretized into a few representative BE classes. The strong enhancement of binding in amorphous ice implies that water (and other hydrogen‑bonding species) can remain trapped on grain surfaces at higher astrophysical temperatures than previously assumed, potentially affecting the timing of ice sublimation in protostellar envelopes and protoplanetary disks.

Conclusions and Outlook
The study demonstrates that (a) graph‑neural‑network MLIPs enable exhaustive, DFT‑quality sampling of water adsorption on realistic grain–ice systems; (b) substrate chemistry matters most at low coverage, while hydrogen‑bonding dominates once a continuous ice layer forms; and (c) amorphous ice introduces high‑energy binding sites that could significantly alter surface chemistry and desorption behavior in cold interstellar environments. The authors propose extending the methodology to mixed ices (e.g., H₂O/CO₂), other adsorbates (NH₃, CH₃OH), and direct comparison with temperature‑programmed desorption experiments, thereby providing a robust, physics‑based foundation for next‑generation astrochemical networks that explicitly account for surface heterogeneity.