Adsorption Energy Distribution Method

• Adsorption Energy Distribution Method is a framework that replaces a single adsorption energy with a spectrum of site energies, capturing surface heterogeneity.
• It integrates TPD inversion, isotherm analysis, and quantum-chemical sampling to derive detailed energy distributions for predicting adsorption kinetics and thermodynamics.
• Accurate parameter choices, such as the pre-exponential factor and energy grid, are critical to ensuring robust inversion and reliable modeling outcomes.

Searching arXiv for recent and directly relevant papers on the Adsorption Energy Distribution Method. The Adsorption Energy Distribution Method (AEDM) is a family of formalisms that replaces a single adsorption or binding energy with a distribution over site energies on a heterogeneous surface. In this formulation, adsorption heterogeneity is represented as $D(E)$, $\theta(E)$, $f(\epsilon)$, $g_i(E)$, or $p(E)$ depending on the application, and observables such as temperature-programmed desorption (TPD) traces, adsorption isotherms, diffusion coefficients, or astrochemical reaction rates are treated as averages over that distribution rather than as consequences of one “representative” site. Recent arXiv work has instantiated AEDM in four distinct but mathematically related ways: inversion of TPD spectra for SO$_2$ on water ice (Benoit et al., 2 Jun 2026), inversion of a single adsorption isotherm for CO$_2$ in CALF-20 (Gonçalves et al., 10 Jul 2025), quantum-chemical sampling of adsorption complexes on amorphous solid water (ASW) (Roy et al., 30 Aug 2025, Groyne et al., 25 Apr 2025), and occupation-weighted rate-equation frameworks for gas–ice astrochemistry (Furuya, 2024).

1. Core concept and mathematical structure

AEDM starts from the observation that amorphous, polycrystalline, and nanoporous substrates do not provide a single adsorption environment. Local coordination, dangling functional groups, hydrogen-bond topology, pore geometry, and adsorbate orientation create a spectrum of activation or binding energies. In the TPD formulation of Benoit et al., the adsorption energy distribution $D(E)$ represents the fraction of adsorbed molecules occupying sites with desorption activation energy $E$; in the CALF-20 formulation, $f(\epsilon)$ is a capacity density per unit energy; in astrochemical rate-equation models, $\theta(E)$0 is the normalized distribution of available sites for species $\theta(E)$1, while $\theta(E)$2 is the instantaneous occupation PDF over those sites (Benoit et al., 2 Jun 2026, Gonçalves et al., 10 Jul 2025, Furuya, 2024).

At the kinetic level, the common starting point is the Polanyi–Wigner law,

$\theta(E)$3

with linear heating $\theta(E)$4 giving

$\theta(E)$5

A measured TPD signal is then a superposition of site-resolved fluxes rather than the signature of a single barrier (Benoit et al., 2 Jun 2026).

At equilibrium, the same principle appears in integral isotherm models. In the CALF-20 paper, AEDM is written as

$\theta(E)$6

where the local occupancy obeys a Langmuir form,

$\theta(E)$7

The total saturation density is recovered by integrating the distribution:

$\theta(E)$8

In this sense, AEDM is not a single inversion algorithm but a distribution-first representation of heterogeneous adsorption. This suggests that its unifying content lies in the averaging over local energetics, whereas the inversion machinery depends on whether the data source is a desorption trace, an adsorption isotherm, or a set of electronic-structure calculations (Gonçalves et al., 10 Jul 2025).

A third mathematical form appears when the distribution itself is built directly from sampled adsorption complexes. In that case, the energy PDF is often modeled as Gaussian,

$\theta(E)$9

or, when multimodality is present, as a Gaussian mixture. This representation is used both for quantum-chemical BE distributions on ASW and for subsequent rate averaging in astrochemical models (Roy et al., 30 Aug 2025, Groyne et al., 25 Apr 2025).

2. TPD-based inversion on heterogeneous surfaces

For desorption experiments, AEDM is an inverse kinetic method. Benoit et al. measured TPD of SO$f(\epsilon)$0 deposited at 80 K on polycrystalline Au, compact amorphous solid water (c-ASW), and crystalline H$f(\epsilon)$1O films under ultra-high vacuum, with base pressure $f(\epsilon)$2 mbar, heating rate $f(\epsilon)$3 K/min $f(\epsilon)$4 K s$f(\epsilon)$5, and QMS detection calibrated by integrated areas (Benoit et al., 2 Jun 2026).

Regime	Kinetic assumption	Reported implementation
Multilayer	Zeroth order, $f(\epsilon)$6	Rising-edge $f(\epsilon)$7 vs $f(\epsilon)$8 regression above 90 K
Submonolayer	First order, $f(\epsilon)$9	Multi-bin forward model with constrained least squares

In the multilayer regime, the analysis uses

$g_i(E)$0

From the leading edges, the reported values are $g_i(E)$1 meV and $g_i(E)$2 s$g_i(E)$3 on Au, and $g_i(E)$4 meV and $g_i(E)$5 s$g_i(E)$6 on c-ASW. In the submonolayer regime, the spectrum is decomposed over a discretized energy axis,

$g_i(E)$7

with each bin obeying

$g_i(E)$8

The inversion determines the initial weights $g_i(E)$9 under non-negativity and mass-balance constraints, minimizing the sum of squared residuals between simulated and experimental fluxes (Benoit et al., 2 Jun 2026).

The parameterization is unusually explicit. Benoit et al. fixed $p(E)$0 s$p(E)$1 from transition-state theory at $p(E)$2 K, using $p(E)$3 amu, $p(E)$4 amu·Å$p(E)$5, $p(E)$6 amu·Å$p(E)$7, $p(E)$8 amu·Å$p(E)$9, $_2$0, and $_2$1 m$_2$2. The energy grid was $_2$3 meV with $_2$4 meV, and fitting was restricted to $_2$5 K in order to exclude volcano desorption from ASW crystallization and co-desorption with H$_2$6O (Benoit et al., 2 Jun 2026).

The extracted distributions distinguish Au from water ice. Au exhibits a single broad distribution peaking near $_2$7 meV, whereas water ice shows bimodal distributions consisting of a physisorbed population near $_2$8–450 meV and a more strongly bound population near $_2$9 meV. On c-ASW, the higher-energy shoulder becomes more prominent as coverage decreases; the mean increases from 424 meV at 1.00 ML to 456 meV at 0.07 ML. Crystalline H$_2$0O shows similar bimodality with only minor differences from c-ASW. The reported average binding energy of SO$_2$1 on water ice is $_2$2 meV, corresponding to $_2$3 kJ mol$_2$4 (Benoit et al., 2 Jun 2026).

Methodologically, this TPD implementation makes two points clear. First, full-spectrum forward modeling is not equivalent to peak-picking or Redhead-only analysis in the submonolayer regime. Second, the solution is sensitive to the prefactor: changing $_2$5 by an order of magnitude prevents convergence or produces poor fits. The paper therefore treats $_2$6 not as a nuisance parameter but as a central identifiability constraint (Benoit et al., 2 Jun 2026).

3. Isotherm-based AEDM and thermodynamic–kinetic parameter extraction

A second major AEDM lineage infers the energy distribution from a single adsorption isotherm rather than from desorption kinetics. In the CALF-20 study, the method derives from statistical mechanics and generalizes the Langmuir picture from a finite number of site classes to a continuous distribution over adsorption energies. The local equilibrium constant is written as

$_2$7

and the isotherm becomes

$_2$8

Discretization over $_2$9 bins converts the integral into a linear mixture of local Langmuir occupancies (Gonçalves et al., 10 Jul 2025).

The inversion is performed by an expectation–maximization (EM) update over the discretized energy axis. Initial guesses can be uniform or Gaussian. Convergence is assessed by a least-squares criterion on the difference between calculated and measured loadings. The paper emphasizes that the energy bounds must be wide enough and that limited pressure coverage can create artifacts; experimental datasets restricted to 1 bar are shown to generate spurious features and incorrect capacities. For CALF-20, a wide-pressure-range GCMC isotherm at 293.15 K yields a robust bimodal AED with site energies near 36.3 and 26.1 kJ/mol, which are corrected to $D(E)$0 kJ/mol and $D(E)$1 kJ/mol at 0 K (Gonçalves et al., 10 Jul 2025).

Integration under the two peaks gives $D(E)$2 mmol/g and $D(E)$3 mmol/g, with total saturation density $D(E)$4 mmol/g, close to the GCMC value $D(E)$5 mmol/g. Using an adsorbate vibration amplitude $D(E)$6 nm and

$D(E)$7

the inferred vibrational frequencies are $D(E)$8 THz and $D(E)$9 THz. The same parameter set is then used to predict isotherms between 293.15 and 393.15 K, Henry’s constants,

$E$0

and enthalpies of adsorption at infinite dilution (Gonçalves et al., 10 Jul 2025).

The CALF-20 paper extends AEDM beyond equilibrium thermodynamics by coupling it to transition-state-theory transport. For the dominant diffusion paths, it reports free-energy barriers of approximately 11 kJ/mol along [011] and 38 kJ/mol along [100], with the lower barrier linked to the difference in site energies, $E$1. The resulting self-diffusion coefficients are predicted as $E$2 m$E$3/s along [011] and $E$4 m$E$5/s along [100], qualitatively matching molecular dynamics trends and declining with pressure as the mean free path decreases from $E$6 nm toward $E$7 nm (Gonçalves et al., 10 Jul 2025).

In this framework, AEDM is not merely a descriptive fit to heterogeneity. It is a parameter-reduction scheme: one wide-range isotherm provides site energies, site capacities, and vibrational frequencies, after which thermodynamic extrapolation and approximate kinetic prediction become possible. The paper nevertheless restricts this claim to type-I isotherms and notes that framework rigidity, omitted barrier recrossing, and neglected temperature dependence of barriers limit the transport model (Gonçalves et al., 10 Jul 2025).

4. Quantum-chemical construction of adsorption-energy distributions on amorphous water ice

A third AEDM family constructs the distribution directly from electronic-structure calculations on many adsorption geometries. In the ASW-cluster study of alcohols, thiols, aldehydes, thioaldehydes, and their precursors, the binding energy is defined as a positive stabilization energy,

$E$8

with zero-point correction and BSSE removal giving

$E$9

Values are reported in Kelvin, i.e. divided by $f(\epsilon)$0. The surface model consists of nine distinct $f(\epsilon)$1 ASW clusters generated from crystalline ice by heating to 300 K and quenching to 10 K, following Shimonishi et al. (2018). Formaldehyde was used as the benchmark species, optimized on 75 distinct dangling-H sites across the nine clusters, while other species were sampled on ASW(1) over seven dangling-H sites or, for OH/SH donors, over both seven dangling-H and seven dangling-O sites (Roy et al., 30 Aug 2025).

For H$f(\epsilon)$2CO across the 75 sites, the +BSSE distribution is Gaussian with $f(\epsilon)$3 K and $f(\epsilon)$4 K. For low-temperature occupation, the study instead averages the highest d–H site in each cluster and obtains a representative +BSSE mean of 3696 K. Morphology dispersion is then transferred to other species through scaling factors: 1.172 for the mean and 2.385 for $f(\epsilon)$5 on +BSSE, with CO and CS using previously validated factors 1.177 for the mean and 0.721 for $f(\epsilon)$6. In subsequent modeling, each species is assigned a Gaussian $f(\epsilon)$7 truncated to $f(\epsilon)$8 (Roy et al., 30 Aug 2025).

The resulting distributions are chemically structured rather than generic. Representative mean $f(\epsilon)$9 values are CO, $\theta(E)$00 K; H$\theta(E)$01CO, $\theta(E)$02 K; CH$\theta(E)$03OH, $\theta(E)$04 K; CH$\theta(E)$05SH, $\theta(E)$06 K; CH$\theta(E)$07OH, $\theta(E)$08 K; and CH$\theta(E)$09SH, $\theta(E)$10 K. Oxygen-bearing species generally have higher means than sulfur analogues, as in CH$\theta(E)$11OH versus CH$\theta(E)$12SH and CH$\theta(E)$13OH versus CH$\theta(E)$14SH, while the study notes the exception CS $\theta(E)$15 CO, attributed to the larger dipole of CS and altered polarizability balance with the surface (Roy et al., 30 Aug 2025).

A more explicitly method-testing variant appears in the ONIOM-based ASW slab study of NH$\theta(E)$16, CO, and CH$\theta(E)$17. That work uses a 2000-water heat–quench TIP4P/2005 slab, hemispherical cuts, and ONIOM(B3LYP-D3(BJ)/6-311+G**:GFN2-xtb) with an 8 Å high-level model radius and a converged real hemisphere radius of 16 Å. Sampling is twofold: 100 regularly spaced spatial sites on a 10×10 grid and multiple starting orientations per site. The required orientation counts are species dependent: three for NH$\theta(E)$18, two for CO, and one for CH$\theta(E)$19. NH$\theta(E)$20 is best fit by a two-component Gaussian mixture, with a low-energy component spanning roughly 15–30 kJ/mol and carrying 16.3% of the probability mass, and a main component spanning approximately 35–63 kJ/mol; CO and CH$\theta(E)$21 are adequately described by single Gaussians centered near $\theta(E)$22 and $\theta(E)$23 kJ/mol, respectively (Groyne et al., 25 Apr 2025).

This computational lineage shows that “distribution” does not imply a single canonical shape. Depending on the adsorbate and on the sampling protocol, the inferred AED can be a simple Gaussian, a truncated Gaussian, an empirical histogram, or a Gaussian mixture. It also establishes that orientation sampling is not a secondary detail: in NH$\theta(E)$24, the coexistence of donor-like and acceptor-like motifs is precisely what generates the observed bimodality (Groyne et al., 25 Apr 2025).

5. Embedding AEDM into astrochemical rate equations

AEDM has also been reformulated as a way to modify reaction-network kinetics without introducing explicit site classes into the ODE system. In the rate-equation framework for gas–ice astrochemistry, the normalized distribution of available sites for species $\theta(E)$25 is $\theta(E)$26, the occupied fraction of those sites is $\theta(E)$27, and the total surface coverage is

$\theta(E)$28

The occupation PDF is defined as

$\theta(E)$29

and is approximated by a Fermi–Dirac-like quasi-steady form,

$\theta(E)$30

with $\theta(E)$31 determined by balancing site population and depletion processes (Furuya, 2024).

The resulting rate coefficients are occupation-weighted averages rather than single Arrhenius factors. For thermal desorption,

$\theta(E)$32

and for hopping,

$\theta(E)$33

The directional diffusion barrier is taken as

$\theta(E)$34

which enforces microscopic reversibility. Representative values used in that work are $\theta(E)$35 s$\theta(E)$36, default $\theta(E)$37, with $\theta(E)$38 for H and H$\theta(E)$39 and $\theta(E)$40 for CO (Furuya, 2024).

The principal algorithmic claim is that surface heterogeneity can be included without increasing the number of ODEs. The computational cost is moved to evaluating energy-grid sums, scaling as $\theta(E)$41 per species per RHS call. The paper reports that $\theta(E)$42 is typically sufficient for $\theta(E)$43, $\theta(E)$44 for $\theta(E)$45, and that 1 Myr evolutions remain practical with a stiff solver such as LSODES (Furuya, 2024).

The physical impact is substantial. In the hydrogen-only tests, the distribution-aware method reproduces the site-resolved exact RE solution while conventional single-energy RE fails at $\theta(E)$46–16 K because it misses deep-site trapping. In the small H/O/CO network at 10 K, O$\theta(E)$47 and CO$\theta(E)$48 are enhanced by orders of magnitude through shallow-site diffusion of O and CO. In the large 709 gas + 302 ice network, differences for major ices are within a factor of a few at 10 K, but at 16 K the distribution model gives H$\theta(E)$49O ice more abundant than CO$\theta(E)$50, reversing the conventional ordering, and at 20 K yields H$\theta(E)$51CO and CH$\theta(E)$52OH higher by roughly an order of magnitude (Furuya, 2024).

The ASW-distribution study using the Rokko code with the REPDF extension applies the same averaging idea more directly to desorption and diffusion:

$\theta(E)$53

$\theta(E)$54

Under dark-cloud conditions with $\theta(E)$55, $\theta(E)$56, $\theta(E)$57, and $\theta(E)$58, 15, and 20 K, the paper reports significant differences in predicted ice-phase abundances between RE and REPDF. In that setup, methanol is lower in REPDF than in RE, a trend explicitly attributed to omission of an H-atom distribution; the same paper notes that Furuya (2024) found the opposite tendency when H’s distribution was included (Roy et al., 30 Aug 2025).

6. Sensitivities, limitations, and interpretive issues

AEDM is best viewed as a controlled approximation to heterogeneous adsorption rather than as a direct readout of microscopic site energies. Several papers identify distinct sources of non-uniqueness. In TPD inversion, the pre-exponential factor is critical: for SO$\theta(E)$59 on water ice, changing $\theta(E)$60 by one order of magnitude prevents convergence or produces poor fits. The reported widths $\theta(E)$61 therefore reflect both genuine site heterogeneity and deviations from an ideal first-order law. Volcano desorption, co-desorption, and possible thermal reactivity near 100–120 K are additional sources of spectral distortion, which is why Benoit et al. restrict inversion to $\theta(E)$62 K and treat the $\theta(E)$63 meV shoulder cautiously as strong SO$\theta(E)$64–H$\theta(E)$65O interactions, reaction products, or processes near water-desorption onset (Benoit et al., 2 Jun 2026).

In isotherm-based AEDM, the main sensitivity lies in data coverage and model closure. The CALF-20 study shows that narrow pressure windows generate spurious peaks and incorrect capacities, and its predictive thermodynamics are explicitly tied to type-I behavior. Its diffusion model further assumes one-dimensional TST, rigid framework dynamics, $\theta(E)$66 for the Bennett–Chandler correction, and temperature-independent barriers. These are simplifying assumptions rather than generic consequences of AEDM (Gonçalves et al., 10 Jul 2025).

In quantum-chemical AEDM, the dominant issues are finite sampling and transferable morphology corrections. The ASW-cluster study derives morphology scaling from H$\theta(E)$67CO and propagates it to chemically distinct species; the paper presents this as a pragmatic approximation, not as a universally validated mapping. The ONIOM slab study similarly shows that inadequate orientation sampling can qualitatively miss multimodality, especially for NH$\theta(E)$68, and that very small real-system sizes can generate artifacts such as the outlier NH$\theta(E)$69/site-2 binding energy above 80 kJ/mol at $\theta(E)$70 Å. Its recommended surrogate $\theta(E)$71 is a practical expedient for large datasets, not a replacement for explicit frequencies in all systems (Roy et al., 30 Aug 2025, Groyne et al., 25 Apr 2025).

In rate-equation embeddings, the quasi-steady occupation PDF assumes that site occupation adapts instantaneously to current physical conditions and that species do not compete explicitly for the same sites. The framework also neglects site topology beyond a maximally connected random-site picture and does not by itself resolve stochasticity when the mean population per grain is below unity. These approximations explain why the method preserves ODE dimensionality; they are also the reason it cannot reproduce every correlation captured by explicit multi-site RE or kinetic Monte Carlo (Furuya, 2024).

A recurring misconception is that AEDM always yields a single Gaussian “distribution of adsorption energies.” The recent literature does not support that simplification. Depending on the observable and on the system, the inferred distribution can be bimodal, truncated Gaussian, empirical histogram, or Gaussian mixture. Another misconception is that the method necessarily provides unique microscopic site assignments. In practice, AEDM usually provides an effective energetic decomposition whose interpretability depends on prefactor choices, sampling completeness, and the extent to which neglected physics—reactivity, entrapment, lateral interactions, framework flexibility, or co-adsorption—can bias the inferred shape (Benoit et al., 2 Jun 2026, Gonçalves et al., 10 Jul 2025, Groyne et al., 25 Apr 2025).

Across these variants, AEDM has converged on a common methodological role: it supplies a compact representation of heterogeneous adsorption that can be inferred from desorption traces, equilibrium uptake, or quantum chemistry and then propagated into thermodynamic, kinetic, and astrochemical models. What differs from one implementation to another is not the core idea of energy-distribution averaging, but the inversion assumptions that connect that idea to a specific dataset or simulation target.