While the approach holds promise, dynamic catalysis presents unique modeling challenges. For example, time-varying parameters introduce complexity into surface coverage dynamics, which are typically modeled using coupled systems of ordinary differential equations (ODEs) that must be solved numerically over multiple periods to reach cyclic steady state. In our previous work2, we used a simultaneous simulation approach to directly compute cyclic steady state behavior, as well as gradient-free and gradient-based optimization algorithms to design waveforms that could lead to the higher rate enhancement within bounds. In this work, we revisit the dynamic catalysis problem fundamentally by turning to analytical methods to obtain results and physical insights into the behavior of dynamic catalysis systems.
We develop and validate a multiple time-scale (or multiscale, for short) analytical framework for predicting surface coverage behavior under high-frequency periodic modulation. This approach decomposes the system response into averaged and oscillatory components, yielding closed-form expressions for surface coverage fluctuations. We compare these results with numerical simulations for the same system and assess the conditions under which the analytical model is accurate. Prior to this work, Foley and Razdan3 explored similar goals through different analytical techniques, mainly by describing two distinct mechanisms or regimes, which they call quasi-static and stepwise. These methods were implemented in a subsequent work that investigated the relationship between forcing frequency and turnover frequency response4.
In contrast to our previous work based on numerical methods, here we introduce analytical approaches for solving, initially, simplified dynamic catalysis systems. The proposed methods are divided into two steps: calculation of the time-averaged surface coverage of adsorbates, and retrieval of the limit cycle, i.e., the periodic fluctuation around the average. We initially focus on reactive systems that lead to linear systems of ODEs, under the assumption of constant gas-phase concentrations. In this abstract, we show a condensed derivation for a simple system, for which results can be directly compared with the ones obtained by Foley and Razdan3.
Assuming a reactive system with two irreversible elementary steps as follows:
A(g) + * → I* → B(g) + *
Since I* is the only adsorbate and the constant concentrations in the the gas-phase molecules are equal to unity, a single ODE describing the evolution of the system over time can be derived:
dθI/dt = − k2 (t) θI + k1 (t) (1−θI ), (1)
in which k1 is the rate constant of the first reaction, and k2 is the rate constant of the second reaction. By solving analytically in the long-time limit, where a periodic steady-state is attained, , and averaging, we reach
⟨θI ⟩ = ⟨k1 ⟩ / (⟨k1 ⟩ + ⟨k2 ⟩, (2)
from which we obtain the average of the oscillations for the surface coverage fraction by I in the cyclic steady state. Next, to solve for the limit cycle, we apply a multiscale method, separating the system into slow and fast time scales. The slow time variable s captures the long-term average behavior, while the fast scale τ captures rapid oscillations within each period. Once again, we show here a condensed and summarized version of the derivation.
Assuming that the rate constants follow square-wave profiles of the form k1 = A1f1(t) + offset1, k2 = A2f2(t) + offset2, where A is the amplitude of each square wave for k, offset is the value that shifts the wave vertically and f are the piecewise square functions that describe the behavior of the rate constants. We will assume s = A1t, and δ = A2/A1. For this case, the value of both the offsets are equal to zero. The ODE for the slow time scale becomes:
dθI /ds = − δ f2 θI + f1 (1 − θI) (3)
We then assume that θI has the multiple-scale expansion θI = θI,0 (s,τ) + ε θI,1 (s,τ) + ⋯, which depends on both s and τ. By applying the expansion to (3) and grouping terms by powers of ε, we find that the leading-order terms arise at O(ε0), while the next correction appears at O(ε1). Grouping the leading order terms, we obtain an ODE that describes the sum of both contributions. By averaging such ODE over the fast time scale we retrieve an averaged version of (3). Subtracting that version from the ODE for both contributions finally leads to the differential equation for the fast time scale oscillation:
dθI,1/dτ = − δ (f2 (τ) − ⟨f2 ⟩) ⟨θI⟩ + (f1 (τ) − ⟨f1 ⟩) (1− ⟨θI⟩) (4)
This equation is solved by integrating both sides. Since f₁ and f₂ are piecewise functions, the solution is also piecewise and depends on the shape of the waveform. The fluctuations are given by closed-form integrals over the waveform, and constants of integration are determined by enforcing zero-mean cyclic behavior. Hence, the final output of this framework is an explicit analytical prediction of the periodic surface response to square-wave forcing. When employing our proposed analytical method to this system, we were able to obtain similar results to those obtained by Foley and Razdan3.
In this abstract, we omit the derivation for systems of coupled linear ODEs, but the approach is similar. For the time-averaged behavior, we set the ODEs equal to zero, average each term, and solve the resulting linear system using standard linear algebra. As an example, in a more complex system (A + * ⇌ A* ⇌ B* ⇌ B + *), with 0.2 eV amplitude and 1.4 eV offset for the BE oscillation, the analytical and simulated average coverages differ by less than 0.5%. The analytical limit cycle solution also matches the simulated result in shape and magnitude.
Overall, we observe a good agreement between the analytical and numerical results under the key assumptions of the multiscale method: high frequency relative to intrinsic kinetics, and small amplitude for the oscillation responses, such that surface coverages remain close to their average values. As amplitudes increase or frequencies decrease, the analytical solution loses accuracy, specifically when coverage dynamics no longer remain close to the average. This provides a quantitative criterion for the applicability of the multiscale methods.
This work presents an analytical framework for predicting surface coverage behavior in dynamic catalysis under high-frequency periodic forcing. By separating slow and fast dynamics through multiscale analysis, we obtain closed-form expressions for both average coverage and limit cycle fluctuations. Our results are validated by the numerical solutions, under conditions that can be encountered in dynamic catalytic systems design. The method is particularly valuable in regimes where simulations are computationally expensive but the assumption of high frequency holds. Future work will further explore complex mechanisms, including those that lead to non-linear ODEs, and different waveform shapes. Ultimately, this approach bridges the gap between simulation and design, offering fast, interpretable tools for dynamic/programmable catalysis systems.
1- Ardagh, M. A.; Abdelrahman, O. A.; Dauenhauer, P. J. Principles of Dynamic Heterogeneous Catalysis: Surface Resonance and Turnover Frequency Response. ACS Catalysis 2019, 9, 6929–6937.
2- Tedesco, C. C.; Kitchin, J. R.; Laird, C. D. Cyclic Steady-State Simulation and Waveform Design for Dynamic/Programmable Catalysis. The Journal of Physical Chemistry C 2024, 128, 8993–9002
3- Foley, B. L.; Razdan, N. K. Clarifying mechanisms and kinetics of programmable catalysis. iScience 2024, 27, 109543.
4- Canavan, J. R., Hopkins, J. A. Hopkins, Foley B. L., Abdelrahman, O. A., Dauenhauer, P. J. Catalytic Resonance Theory: Turnover Efficiency and the Resonance Frequency. ACS Catalysis 2024, 15, 653−663.