# Unbiased backward Monte Carlo method for addressing stiff chemical kinetics

Abderrahim Sahim, Marion Carrier, Mouna El Hafi, Simon Eibner, Yaniss

Nyffenegger-pere, Stéphane Blanco, Richard Fournier

► To cite this version:

Abderrahim Sahim, Marion Carrier, Mouna El Hafi, Simon Eibner, Yaniss Nyffenegger-pere, et al. Unbiased backward Monte Carlo method for addressing stiff chemical kinetics. Chemical Engineering Journal, 2025, 520, pp.165009. 10.1016/j.cej.2025.165009. hal-05165554

HAL Id: hal-05165554

https://imt-mines-albi.hal.science/hal-05165554v1

Submitted on 15 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Unbiased backward Monte Carlo method for addressing stiff chemical kinetics

Abderrahim Sahim<sup>a,b,*</sup>, Marion Carrier<sup>a,b</sup>, Mouna El Hafi<sup>a,d</sup>, Simon Eibner<sup>a,b</sup>, Yaniss Nyffenegger-péré<sup>c,d</sup>, Stéphane Blanco<sup>d</sup>, Richard Fournier<sup>d</sup>

## ARTICLE INFO

### Keywords:

- Monte Carlo

- Stiffness

- Chemical kinetics

- Quasi-steady-state

- Partial equilibrium

## ABSTRACT

In this research, a new Monte Carlo-based approach is introduced to rethink kinetic models for reaction schemes involving coupled first-order reactions. The study examines two specific schemes: one with exclusively irreversible reactions and another with a reversible reaction, all within an isothermal regime. The kinetic model is reformulated under a system of coupled expectation-based formulations and subsequently solved via Monte Carlo simulation. The solutions are visualized in path space to demonstrate the double randomization technique, also known as the walk-on equation. This approach allows for a qualitative analysis of the algorithm's behavior under varying conditions. Two scenarios are analyzed: a non-stiff scenario and a stiff scenario, considering timescales related to quasi-steady-state (QSS) and partial equilibrium (PE). The new method provides results consistent with the analytical solution for the non-stiff case and the QSS stiff case. For the PE case, zero variance and classical approximations are employed to derive the optimal probability density functions, capturing the dynamics of fast reactions and eliminating rare events, which leads to good convergence of the method and comparable results to analytical solutions.

## 1. Introduction

Modeling the temporal multi-scale dynamics of chemical systems has been a longstanding challenge. Classical approaches from the early stages of research involved analytical approximations, notably the quasi-steady-state (QSS) assumption [1,2]. This assumption states that on the relevant time scale, certain intermediate species evolve far more rapidly than others, and their rates are approximately zero. Additionally, the partial equilibrium (PE) assumption was introduced, suggesting that certain fast reactions are constantly in equilibrium [3]. Several subsequent works attempted to demonstrate the validity of these approaches with a mathematical foundation from the perspective of singular perturbation theory [4–6]. However, Turany et al. [7] illustrated the limitations of the QSS approximation when studying propane pyrolysis. They showed how applying this approximation to some intermediate species can induce significant errors affecting the resolution of species where QSS was not applied. Farow and Edelson [8] used QSS to model a photochemical smog chamber, confirming that it introduced uncertainty that could dominate experimental uncertainties and affect the determination of reaction mechanisms and

rate constants. Another study demonstrated that QSS fails to capture oscillating dynamics in open enzyme-substrate reaction reduced systems [9]. A recent study [10] demonstrated that in the context of the Michaelis–Menten reaction mechanism, QSS is not broadly valid except when the initial substrate concentration is small. The partial equilibrium (PE) approximation has not been extensively examined, as it requires a somewhat more elaborate implementation procedure, but some studies have attempted to determine its limitations, particularly in the modeling of spatially heterogeneous systems [11,12] due to the oversimplification of the process.

With the advent of deterministic and stochastic numerical solution methods for modeling reactive systems, the multiscale problem has reappeared, particularly in the domain of deterministic integration methods under the name “stiffness” [13]. Since 1971, Gear [14] was the first to address stiffness using the Backward Differentiation Formula, a numerical method based on an implicit integration scheme. This method is considered the most popular method for integrating stiff problems, and it has led to the development of various solvers [15,16]. The use of implicit integration methods is due to their unconditional

<footnote><sup>a</sup> RAPSODEE - UMR CNRS 5302 - IMT Mines Albi - Campus Jarlard - 81013 Albi, CT Cedex 09, France</footnote>

<footnote><sup>b</sup> Institut national polytechnique de Toulouse, 6 allée Emile Monso -, BP 34038, 31029 Toulouse Cedex 4, France</footnote>

<footnote><sup>c</sup> Instituto de Astrofísica de Andalucía, Consejo Superior de Investigaciones Científicas (CSIC), Glorieta de la Astronomía s/n, 18008 Granada, Spain</footnote>

<footnote><sup>d</sup> LAPLACE - UMR CNRS 5213 - Université Paul Sabatier, Toulouse, France</footnote>

<footnote>* Correspondence to: IMT Mines Albi - Centre RAPSODEE (UMR CNRS 5302), Campus Jarlard, 81013 Albi Cedex 09, France.
E-mail address: abderrahim.sahim@mines-albi.fr (A. Sahim).</footnote>
## Nomenclature

### Greek letters

$\epsilon_i^{MC}$ relative error associated with the estimation of the mass fraction of species i (-)

$\epsilon_i$ relative standard deviation associated with the estimation of the mass fraction of species i (-)

$\pi_i$ probability density function associated with species i (s⁻¹)

$\sigma$ analytical standard deviation (%)

$\sigma_i^{MC}$ numerical standard deviation associated with the estimation of the mass fraction of species i (%)

$\tau$ random time (s)

### Sub-/Superscripts

$dim$ dimensionless

$i$ index (species)

$IS$ importance sampling

$KC$ kinetic control

$MC$ Monte Carlo

$obs$ observation

$PE$ partial equilibrium

$RE$ rare event

$ZV$ zero variance

$0$ initial

### Abbreviations

$DAEM$ Distributed Activation Energy Model

$kMC$ kinetic Monte Carlo

$ODE$ ordinary differential equation

$PE$ partial equilibrium

$QSS$ quasi-steady-state

### Roman letters

$\mathbb{E}$ mathematical expectation (-)

$\mathcal{T}$ random variable associated with time (s)

$W_{Y_i}$ random variable whose expectation is the mass fraction $Y_i$ (%)

$k$ reaction rate constant (s⁻¹)

$m_i$ mass of species i (kg)

$N$ number of samples (-)

$S_i$ source term in the governing equation of the mass fraction of species i (s⁻¹)

$t$ time (s)

$w_{RE}$ weight of rare events (%)

$w_{Y_i}$ realization of $W_{Y_i}$ (%)

$Y_i$ mass fraction of species i (%)

stability regardless of the choice of integration step, and their incorporation of future information into each time step calculation, allowing for larger steps and effective dampening of rapid oscillations [15]. The use of explicit methods to integrate stiff systems leads to instability if the chosen integration step is of the order of the slow dynamics, or very slow convergence if the integration step is of the order of fast dynamics [15].

Another approach has been proposed by Gillespie [17] who developed

kinetic Monte Carlo (kMC) to solve the master equation for modeling well-stirred chemically reacting systems, the problem of stiffness reappeared in a stochastic form. In stochastic modeling, time scales associated with QSS are well handled by kMC [18]. However, for time scales related to partial equilibrium (PE), kMC fails to capture slow reactions and spends the majority of time sampling rapid processes. Various proposals have emerged to address this stochastic stiffness, including the removal of fast variables from the master equation [19]. A hybrid analytical kMC method was proposed to address stiffness under PE by analytically calculating concentrations of species involved in fast reversible processes [20]. Another proposal to manage this stiffness is the tau-leaping method developed by Gillespie [21], and subsequently with its implicit implementation [22] to advance with large time steps, reducing simulation costs. Vlachos [23] developed a lumping method with an extension of the method to treat stiffness under partial equilibrium conditions with the work of Synder et al. [18], in a way that each fast reversible reaction can be treated as a single reaction called a "net event" to reduce time scale disparities with other reactions. Other studies have attempted to combine both QSS and PE approximations with stochastic modeling in order to mitigate timescale disparities [24,25].

The present work is part of a current trend involving the latest developments in the Monte Carlo method [26]. The most recent developments reflect a renewed and diversified interest in Monte Carlo methods within chemical engineering. Pahija et al. [27] revisited the fundamental principles of Monte Carlo methods – such as expectation estimation and the law of large numbers – and provided a comprehensive overview of their evolution into versatile tools for uncertainty quantification and experimental design. The review also presents a detailed survey of chemical engineering applications, with numerous references to kinetic Monte Carlo (KMC) methods, highlighting their growing role alongside deterministic approaches in modeling and decision-making. Sediawan et al. [28] further demonstrated their practical relevance through detailed applications in reactor design, process simulation, and risk analysis. Meanwhile, Zichi et al. [29] presented Reaction Network Monte Carlo (RNMC), an open-source framework for kinetic Monte Carlo simulations tailored to handle complex reaction networks, underscoring the increasing computational maturity and accessibility of MC-based tools. These recent contributions illustrate the methodological depth and software maturity that now support the integration of Monte Carlo methods in complex, multiscale chemical systems.

There is a multidisciplinary collective interest in processes that involve large ratios of spatial and temporal scales. As far as space is concerned, the use of applied mathematics and computer graphics libraries allows the computation times of Monte Carlo simulations to be insensitive to geometrical complexity [30–32]. The main scientific origin of all such geometric considerations started with the practical needs of the cinema industry up to present stage where all image synthesis is exclusively relying on Monte Carlo approaches. More recently, there has been even more developments concerning phenomenological complexity such as nested models, monolithic models transitioning from one phenomenology to another through path space formulations [26]. In chemical kinetics, the issue of coupling is crucial for one primary reason: most processes are significantly driven by thermal effects, and coupling with thermal heat transfer is needed. If we manage to address stiff kinetics with such a Monte Carlo approach then coupling with heat transfer will be straight forward using recent advances based on path-space formulations, defining a single monolithic model that couples conduction, convection, and radiation into statistical paths that can be simulated in complex geometries with the same ease as light paths for image synthesis [33]. As the very same path sampling algorithms also deal with temporal forcing involving scale ratios [30,32], we can envision kinetics coupled to thermal heat transfer with large scale ratios in both geometry and control (temporal forcing). However, to the best of our knowledge the Monte Carlo method encounter the same stiffness
problems as all other methods, and therefore these expected benefits are lost because enriching the small scale descriptions will activate the time scales of the stiffest phenomena (the smallest characteristic times of the interactive kinetics). We address this issue here, aiming to find a Monte Carlo approach that performs well, despite stiffness problems, and preserves at least partially the above envisioned benefits.

In this study, we develop a Monte Carlo method, different from kMC, to deal with stiff reactions models: the Integral Formulation Monte Carlo [34]. This approach reformulates alongside the solution of a differential equation into an integral formulation defined as an expected value. An estimator is associated to this solution alongside a statistical variance. Furthermore, the integral formulation is associated to a path space viewpoint [35] here adapted to the time domain. Consequently, the integral Monte Carlo estimates the solution at a given time by a backward propagation into the time domain, until it reaches the initial condition. It allows tracing back the different physical events that led to the solution of the model at a given time. This facilitates a better understanding of the physical processes driving the system evolution, thereby enhancing the analysis and the understanding. The proof of concept for the use of Integral Monte Carlo has already been demonstrated in a chemical kinetics case for studying lignocellulosic biomass pyrolysis [36] where a single chemical reaction is treated, but averaged over a global parameter set, that accounts for a large number of chemical reactions: the DAEM model. In addition, the model has been solved in an isothermal and non-isothermal case, i.e. when the rate constant (given by Arrhenius) depends on temperature change over time. In the present study, we will extend the treatment of one reaction to several coupled chemical reactions, and for the sake of simplicity, an isothermal configuration is considered. The methodological aspect of coupling several equations with Monte Carlo has already been the subject of many previous studies. If the solution of an equation is an expectation, the coupling between several expectations, so different equations, is achieved through a technique called double randomization, referred to “the walk in equation” in the literature [37,38]. In the present work, the coupling between several chemical reactions could be named "walk on reaction". One of the main original aspect of this work will therefore be the use of integral Monte Carlo to handle several coupled chemical reactions.

In this study, we will explore the applicability of the Monte Carlo Integral Formulation method to solve the reaction-rate equation, a deterministic kinetic model, for an isothermal batch reaction series. We will investigate cases involving successive irreversible first-order reactions and reversible reactions. The kinetic model results in a set of coupled first-order ordinary differential equations. We will consider different scenarios to model both non-stiff and stiff cases, accounting for timescale disparities associated with quasi-steady-state approximation and partial equilibrium approximation. Our aim is to demonstrate how stiffness manifests when treated by this Monte Carlo method and propose approaches based on statistical techniques such as variance reduction principles and importance sampling to address stiffness. To our knowledge, this is the first time such a temporal multi-scale problem will be addressed using this Monte Carlo approach.

## 2. Materials and methods

### 2.1. Irreversible case

A very common kinetic scheme involving a series of two irreversible first-order reactions and three species was first considered: a species A, which is solely consumed and transformed into an intermediate B, which subsequently reacts to give the final product, C, which is only formed as presented in Fig. 1

The system is assumed to be well-stirred and maintained in thermal equilibrium at a constant temperature. Given that the rate constant of each reaction is temperature-dependent under the thermal equilibrium

<figure><img src="image_4.png" /><figcaption>Fig. 1. Reaction scheme of the irreversible case.</figcaption></figure>

assumption, there will be consequently no dependence of the rate constants, $k_1$ and $k_2$, on time.

At the initial time, it is assumed that only the initial mass of species A is present in the system. The mass fraction for each species i, $Y_i(t)$, is defined as follows:

$$Y_i(t) = m_i(t)/m_0 \times 100 \quad (1)$$

where t is the time (s), $m_0$ is the initial mass of reactant A (kg), $m_i$ is the mass of species i (kg).

All reactions are assumed to be either unimolecular reactions or pseudo-unimolecular reactions, especially if the reaction is described in terms of pseudo (or lumped) species, which in turn represent complex materials or mixtures [39]. As a consequence of the previous assumptions, the mathematical model representing the time evolution of the mass fractions consists of a system of coupled first-order ordinary differential equations, formulated as follows:

$$\left\{
\begin{aligned}
\frac{dY_A(t)}{dt} &= -k_1 Y_A(t) \\
\frac{dY_B(t)}{dt} &= -k_2 Y_B(t) + k_1 Y_A(t) \\
\frac{dY_C(t)}{dt} &= k_2 Y_B(t)
\end{aligned}
\right.
\quad (2)$$

With initial conditions at the starting time $t = t_0$ are:

$$\begin{aligned}
Y_A(t_0) &= Y_{A,0} = 1 \\
Y_B(t_0) &= Y_{B,0} = 0 \\
Y_C(t_0) &= Y_{C,0} = 0
\end{aligned}
\quad (3)$$

**From integral formulation to expectation-based formulation.** In this framework, the objective is to estimate statistically the mass fraction of each species at a specific time t. The main question is then to express a given mass fraction, $Y_i$, at any specific time, t, as the expectation of a random variable $W_{Y_i}$.

The Integral Formulation Monte Carlo method relies on expressing the observable as an expectation, which is then evaluated with the arithmetic mean. In other words, it is a pointwise method – also called probe point – that can be used only to evaluate the mass fraction of species i at a given time: $Y_i(t)$. It differs from conventional discretization approaches, as they involve the evaluation of all mass fractions over the whole time range. This expression is mathematically defined as an expected value of a random variable. The integral Monte Carlo method relies on this expected value concept: if a quantity can be expressed as an expected value, it can be approximated by the arithmetic mean of numerous realizations sampled according to a distribution function, $\pi_i(\tau)$. Based on the Law of large numbers, the estimated value tends to the observable, expressed as an expected value.

Moreover, all statistical realizations are independent and tend towards a normal distribution according to the central limit theorem [40]. The Gaussian standard deviation represents the statistical error in our estimation. This aspect underscores a significant advantage of the integral Monte Carlo approach: it allows for the estimation of a quantity along with its corresponding error bar, which decreases as the square root of the number of realizations increases. Therefore, multiplying the realization number by four results in a standard deviation divided by two.

To switch from an integral formulation to an expectation-based formulation, we suggest the examination of the general analytical
solution form of a first-order ordinary differential equation presented in the following form:

$$ \frac{dY_i(t)}{dt} = -kY_i(t) + S_i(t) \qquad (4) $$

with $S_i(t)$ representing the source term, and $k$ denoting the sum of rate constants of the reactions according to which species $i$ is consumed, Further details about the explicit forms of $S_i(t)$ and $k$ are given in the Supplementary Material A.

The general analytical solution [41] for the above ordinary differential equation is expressed as follows:

$$ Y_i(t) = Y_i(t_0) \exp(-k(t-t_0)) + \int_{t_0}^{t} S_i(u) \exp(-k(t-u)) du \qquad (5) $$

Now, the first term represents the propagation in time of the initial condition with exponential decay, corresponding to the solution of the particular case of a homogeneous equation without a source contribution. In our context, this scenario specifically applies to the mass fraction of species A, where no other species contribute to its formation, resulting in the absence of a source term for A. The second term incorporates the influence of different sources over time, ultimately representing the sum of mass fractions of the reactants contributing to the formation of species $i$, each multiplied by the corresponding rate constants. This contribution is a cumulative effect of source terms between times $t_0$ and $t$, attenuated starting from the time $u$ when the source is encountered during the propagation through time until the time $t$. The overall behavior of the solution decays exponentially, primarily governed by the coefficient $k$, representing the consumption term. This insight arises from the mathematical formalism of the Green's function method [42–44], which suggests that the solution to a linear differential equation arises from the propagation of initial conditions and sources. Since the propagator is the exponential term, $\exp(-kt)$, which mainly governs the dynamics of the solution as shown in Eq. (5), the most intuitive approach is to choose a probability density of the same mathematical form as the propagator in order to reformulate the solution under an expectation-based formulation. We will explore the implications of this choice later on. To proceed step by step, let us first reformulate the solution in integral form:

$$ Y_i(t) = \int_{t-t_0}^{t} Y_i(t_0) k \exp(-k\tau) d\tau + \int_{t_0}^{t} S_i(u) \exp(-k(t-u)) du \qquad (6) $$

In the first term of Eq. (6), we already observe the appearance of the term $k \exp(-k\tau)$, which is defined between 0 and $+\infty$ as a probability density. However, since the integral here is defined only over $[t, +\infty[$, an additional step is necessary to extend the domain of integration. By performing a change of variable for the second term, setting $\tau = t - u$, we find:

$$ Y_i(t) = \int_{t-t_0}^{t} Y_i(t_0) k \exp(-k\tau) d\tau + \int_{0}^{t-t_0} S_i(t-\tau) \exp(-k\tau) du \qquad (7) $$

Now, starting from Eq. (7), we can extend the integration domain over the range $[0, +\infty[$ as follows:

$$ Y_i(t) = \int_{0}^{+\infty} k \exp(-k\tau) d\tau \times \begin{cases} [Y_i(t_0)] & \mathcal{H}(\tau > t - t_0) \\ + & \\ [S_i(t - \tau)/k] & \mathcal{H}(\tau < t - t_0) \end{cases} \qquad (8) $$

Notation $\mathcal{H}(x)$ stands for the Heaviside function, which takes a value of 1 for $x > 0$ and a value of 0 for $x < 0$.

We can define the random variable, $\mathcal{T}$, with probability density function, $\pi_i(\tau) = k \exp(-k\tau)$, which allows us to define the random variable, $\mathcal{W}_{Y_i}(\mathcal{T})$, as follows:

$$ \mathcal{W}_{Y_i}(\mathcal{T}) = \mathcal{H}(\mathcal{T} > t - t_0)Y_i(t_0) + \mathcal{H}(\mathcal{T} < t - t_0)S_i^*(t - \mathcal{T}) \qquad (9) $$

where $S_i^*(t - \mathcal{T}) = S_i(t - \mathcal{T})/k$. Indeed, we have constructed an expectation of a random variable, $\mathcal{W}_{Y_i}$, such as $Y_i(t) = \mathbb{E}(\mathcal{W}_{Y_i}(\mathcal{T}))$ with one realization denoted as $w_{Y_i}$, is obtained as follows:

* Sample a random time, $\tau$, following the probability density $k \exp(-k\tau)$;

* if the sampled time exceeds $t - t_0$, that could be interpreted as an observation time, it implies that no further chemical transformation of the species $i$ has occurred. Therefore, the only retained information is the initial condition, $Y_i(t_0)$; otherwise, retain the source term contribution at the time, $t - \tau$: $S_i^*(t - \tau)$.

The above approach represents a probabilization strategy of the Green's function formalism for the solution of a first-order chemical kinetic model. The practice of probabilizing a kinetic-type model has already been proposed in other physical fields to study coupled heat transfer modes [31,33] or gas kinetics [45], with a recent attempt in the realm of biomass pyrolysis kinetics by addressing the model of distributed activation energies (DAEM) using Monte Carlo [36].

**Expectation-based formulation for the mass fractions of the three species A, B and C.** Let us examine the general form of the expectation-based formulation provided in Eq. (8), and explore its specific form for each of the three species: A, B, and C.

For species A, the coefficient characterizing the consumption term is the rate constant $k_1$, and the source term is equal to zero since there are no species feeding into A (see Supplementary Material A). Therefore, the expectation-based formulation for the mass fraction of A is as follows:

$$ Y_{A}(t) = \mathbb{E}(\mathcal{W}_{Y_{A}}(\mathcal{T})) = \int_{0}^{+\infty} k_{1} \exp(-k_{1}\tau) d\tau \times \begin{cases} [Y_{A,0}] & \mathcal{H}(\tau > t - t_{0}) \\ + & \\ [0] & \mathcal{H}(\tau < t - t_{0}) \end{cases} \qquad (10) $$

Thus, a random variable, $\mathcal{W}_{Y_A}(t)$, has been successfully constructed, whose expected value is $Y_A(t)$. To obtain a realization of this random variable, the algorithm described in **Algorithm 1** simply needs to be followed.

**Algorithm 1 $\mathcal{W}_{Y_A}(\mathcal{T})$ sampling algorithm**

1: **Sample** a random time $\tau$ following the probability density $\pi_A(\tau) = k_1 \exp(-k_1 \tau)$.

2: **if** $\tau > t - t_0$:

3: **return** the initial condition $Y_{A,0}$.

4: **else**:

5: **return** 0.

Let us look at the dynamics of the mass fraction of species $i$ and try to understand how this dynamics translates algorithmically when attempting to obtain a realization of the random variable, $\mathcal{W}_{Y_A}$, at $t$. $Y_A$ represents a quantity that decreases over time following an exponential decay, $\exp(-k_1(t - t_0))$. At the initial time, the only contribution to the solution of the mass fraction of A will be its initial condition, and thus the initial condition is the sole outcome for the random variable, $\mathcal{W}_{Y_A}$, and the further we move from this time, the less the contribution will be from the initial condition and the more from the source term. As time tends to infinity, the only contribution that will remain is the source term.

For the species B, the coefficient characterizing the consumption term is the rate constant, $k_2$, and the source term is $k_1 Y_A(t)$ (see Supplementary Material A), therefore the expectation-based formulation is written as follows:
$$Y_B(t) = \mathbb{E}(W_{Y_B}(\tau)) = \int_0^{+\infty} k_2 \exp(-k_2 \tau) d\tau \\ \times \begin{cases} [Y_{B,0}] & \mathcal{H}(\tau > t - t_0) \\ + & \\ [k_1/k_2 Y_A(t - \tau)] & \mathcal{H}(t - t_0 > \tau) \end{cases} \quad (11)$$

The algorithmic interpretation is less straightforward for species B compared to species A because the source requires knowing the mass fraction, $Y_A$, at $t - \tau$. However, this mass fraction is not known beforehand, and we only have its expected formulation. We are then faced with two choices:

One option is to make a naive choice, where for each realization of the random variable, $W_{Y_B}$, we perform a Monte Carlo calculation with a large number of samples to statistically estimate $Y_A(t - \tau)$, and then estimate $Y_B$ with a large number of samples too. For example, if we set the number of samples to $N$ for the evaluation of $Y_A(t - \tau)$ and $M$ for the evaluation of $Y_B(t)$, we would need a total of $N \times M$ samples to estimate $Y_B$ at time $t$. However, this approach breaks down as the integration dimension increases each time, and it would be biased because the value of $Y_A$ would be approximated.

Alternatively, as shown in Algorithm 2, we can proceed with double randomization [46,47], as the expectation of the random variable, $W_{Y_B}$, depends linearly on the expectation of the random variable, $W_{Y_A}$. Instead of evaluating $Y_A$ at $t - \tau$, we will simply retain one realization of $W_{Y_A}$ at $t - \tau$. This allows us to evaluate $Y_B(t)$ without biasing the Monte Carlo calculation and with only $N$ samples, which is much more practical compared to the first option. The double randomization is a specific technique to Monte Carlo methods to deal with linearly coupled equations [37,38,48] as the case between the mass fraction of A and the mass fraction of B, it is a recursive procedure of coupling through nested expectations. Moreover, the double randomization is a very powerful concept, as it enables to think the question of the coupling locally, i.e, seen from the standpoint of a random walk backward in time, performed by random variables [33].

### Algorithm 2 $W_{Y_B}(\tau)$ sampling algorithm

1: **Sample** a random time, $\tau$, following the probability density, $\pi_B(\tau) = k_2 \exp(-k_2 \tau)$.

2: **if** $\tau > t - t_0$:

3: **return** the initial condition $Y_{B,0}$.

4: **else**:

5: **Sample** $W_{Y_A}(t - \tau)$ according to Algorithm 1

6: **return** $k_1/k_2 w_{Y_A}(t - \tau)$

Species C represents the final product, so it is solely formed by species B, thus there is no consumption term (see Supplementary Material A). Exploiting this particularity, we revisit Eq. (7), replacing $k$ with its null value and the source term, $S_C(t)$, with $k_2 Y_B(t)$, yielding:

$$Y_C(t) = \int_{t_0}^{t} k_2 Y_B(t - \tau) d\tau + Y_{C,0} \quad (12)$$

The Eq. (12) entails the direct integration of the kinetic equation governing the mass fraction of species C. To switch from this integral formulation to an expectation-based formulation, the approach involves the introduction of a defined probability density, $\pi_C$, function over the interval from $t_0$ to $t$ as follows:

$$Y_C(t) = \mathbb{E}(W_{Y_C}(\tau)) = \int_{t_0}^{t} \left[ \frac{k_2 Y_B(t - \tau)}{\pi_C(\tau)} + Y_{C,0} \right] \pi_C(\tau) d\tau \quad (13)$$

with, $\int_{t_0}^{t} \pi_C(\tau) d\tau = 1$.

It is essential, of course, to ensure that this probability density encapsulates the dynamics of the integrated quantity, which is none other than $Y_B(t)$. We will examine closely how to choose the probability density function, $\pi_C$; for now, we will settle for the most obvious choice: a uniform probability density between $t_0$ and $t$.

$$A \xrightarrow{k_1/k_{-1}} B \xrightarrow{k_2} C$$

<caption>Fig. 2. Reaction Scheme for the reversible case.</caption>

According to Algorithm 3, to obtain a realization of the random variable, $W_{Y_C}$, one must sample a random time, $\tau$, according to $\pi_C(\tau)$. Subsequently, employing the aforementioned double randomization technique, a realization of the random variable, $W_{Y_B}$, at time, $t - \tau$, is obtained by switching to the expectation-based formulation of the mass fraction of species B. The initial condition, $Y_{C,0}$, is accumulated alongside with the realizations of $W_{Y_C}$.

### Algorithm 3 $W_{Y_C}(\tau)$ sampling algorithm

1: **Sample** a random time $\tau$ following the probability density $\pi_C(\tau)$.

2: **Sample** $W_{Y_B}(t - \tau)$ according to Algorithm 2

3: **return** $\frac{k_2 w_{Y_B}(t - \tau)}{\pi_C(\tau)} + Y_{C,0}$

## 2.2. Reversible case

Let us consider now under the same conditions of the above irreversible reaction pathway a reversible reaction pathway with a reverse reaction between A and B. The reaction pathway of the reverse case is presented in Fig. 2.

Under the same assumptions as the irreversible case, and considering the same initial conditions as those specified in the Eq. (3), the mathematical model for the reaction pathway presented in Fig. 2 is written under the following system of coupled first order ODEs:

$$\left\{ \begin{aligned} \frac{dY_A(t)}{dt} &= -k_1 Y_A(t) + k_{-1} Y_B(t) \\ \frac{dY_B(t)}{dt} &= -(k_2 + k_{-1}) Y_B(t) + k_1 Y_A(t) \\ \frac{dY_C(t)}{dt} &= k_2 Y_B(t) \end{aligned} \right. \quad (14)$$

By employing the same switching procedure from the integral formulation of the general solution to the expectation-based formulation utilized for the irreversible case (see Eq. (8)), and considering the new changes in the reversible model (see Supplementary Material A for the specific forms of $k$ and $S_i$), the expectation-based formulations of the mass fractions of the three involved species are now presented as follows:

For species A:

$$Y_A(t) = \mathbb{E}(W_{Y_A}(\tau)) = \int_0^{+\infty} k_1 \exp(-k_1 \tau) d\tau \\ \times \begin{cases} [Y_{A,0}] & \mathcal{H}(\tau > t - t_0) \\ + & \\ [k_{-1}/k_1 Y_B(t - \tau)] & \mathcal{H}(\tau < t - t_0) \end{cases} \quad (15)$$

For species B:

$$Y_B(t) = \mathbb{E}(W_{Y_B}(\tau)) = \int_0^{+\infty} k_s \exp(-k_s \tau) d\tau \\ \times \begin{cases} [Y_{B,0}] & \mathcal{H}(\tau > t - t_0) \\ + & \\ [k_1/k_s Y_A(t - \tau)] & \mathcal{H}(t - t_0 > \tau) \end{cases} \quad (16)$$

with $k_s = k_2 + k_{-1}$.

And finally for species C:

$$Y_C(t) = \mathbb{E}(W_{Y_C}(\tau)) = \int_{t_0}^{t} \left[ \frac{k_2 Y_B(t - \tau)}{\pi_C(\tau)} + Y_{C,0} \right] \pi_C(\tau) d\tau \quad (17)$$
In the reversible case, the two way coupling between species A and B manifests in their respective expectation-based formulations, where the source of the mass fraction of species A at time, t, depends on the mass fraction of species B at time, t - τ, and vice versa for the mass fraction of B. Therefore, to obtain a realization of the random variable, $W_{Y_A}$, at an observation time, t, if the sampled time, τ, exceeds t - t₀, one must switch to the equation for the mass fraction of species B using the double randomization technique to evaluate $W_{Y_B}$ at new time, t - τ₁. This process may require knowledge of $W_{Y_A}$ obtained through another double randomization at another observation time, t - τ₁ - τ₂, and so on. This process unfolds backward in time. Consequently, the two-way coupling between the expectation-based formulations, of the mass fractions of species A and B, persists throughout the double randomization until reaching the initial condition $Y_{A,0}$ or $Y_{B,0}$ as shown in Algorithm 4 and Algorithm 5.

From a rigorous point of view, the above proposition is based on the estimation of potentially infinitely nested expectations, with an iterated kernel. As we are working on a kinetic-base model, the iterative kernels over the probabilized Fredholm in Eqs. (15) and (16) are contracting, and the solution appears as a development of a Neumann series that converges [33,40].

**Algorithm 4** $W_{Y_A}(\tau)$ sampling algorithm

1: **Sample** a random time $\tau$ following the probability density $\pi_A(\tau) = k_1 e^{-k_1 \tau}$.

2: **if** $\tau > t - t_0$:

3: **return** the initial condition $Y_{A,0}$.

4: **else**:

5: **Sample** $W_{Y_B}(t - \tau)$ according to Algorithm 5

6: **return** $k_{-1}/k_1 w_{Y_B}(t - \tau)$.

**Algorithm 5** $W_{Y_B}(\tau)$ sampling algorithm

1: **Sample** a random time $\tau$ following the probability density $\pi_B(\tau) = k_3 e^{-k_3 \tau}$.

2: **if** $\tau > t - t_0$:

3: **return** the initial condition $Y_{B,0}$.

4: **else**:

5: **Sample** $W_{Y_A}(t - \tau)$ according to Algorithm 4

6: **return** $k_1/k_3 w_{Y_A}(t - \tau)$.

For species C, the algorithmic procedure and interpretation remain the same as the irreversible case, utilizing the same Algorithm 3.

In practical terms, to estimate a mass fraction at an observation time, t, we introduce a mean estimator constructed as follows:

$$Y_i^{MC} = \sum_{j=1}^{N} W_{Y_{i,j}}/N \qquad (18)$$

where N is the number of samples for a Monte Carlo simulation, and $W_{Y_{i,j}}$ are independent and identically distributed random variables related to $W_{Y_i}(\tau)$.

Subsequently, this mass fraction, $Y_i(t)$, expressed as an expectation, $\mathbb{E}(W_{Y_i}(\tau))$, is approximated by the mean estimator, $Y_i^{MC}$, for a large number of samples, N, according to the law of large numbers:

$$Y_i(t) = \lim_{N \to +\infty} Y_i^{MC} \qquad (19)$$

Thus, there exists a relationship between the standard deviation of the mean estimator and the standard deviation of the random variable, $W_{Y_i}(\tau)$, such that:

$$\sigma(Y_i^{MC}) = \sigma_i^{MC} = \sigma(W_{Y_i}(\tau))/\sqrt{N} \qquad (20)$$

where $\sigma(W_{Y_i}(\tau)) = \sqrt{\mathbb{E}(W_{Y_i}(\tau)^2) - \mathbb{E}(W_{Y_i}(\tau))^2}$, is the standard deviation of the random variable $W_{Y_i}(\tau)$.

The standard deviation of the estimator, $\sigma_i^{MC}$, can also be determined through a standard deviation estimator defined according to the central limit theorem as:

$$\sigma_i^{MC} = \sqrt{\frac{1}{N-1}\left(\sum_{j=1}^{N} W_{Y_{i,j}}^2/N - \left(\sum_{j=1}^{N} W_{Y_{i,j}}\right)^2/N\right)} \qquad (21)$$

So, to estimate a mass fraction, $Y_i$, of a given species i at an observation time, t, the procedure outlined in Algorithm 6 is followed.

**Algorithm 6** Monte Carlo algorithm for estimating a mass fraction $Y_i$ at a given time t

1: $Som = 0$ and $Som_2 = 0$

2: **for** j = 0 : N **do**

3: $w_{Y_i} = 0$

4: **Sample** $W_{Y_i}(t)$ according to Algorithms 1, 2, 3, 4 or 5.

5: $Som = Som + w_{Y_i}$

6: $Som_2 = Som_2 + w_{Y_i}^2$

7: **end for**

8: $Y_i^{MC} = Som/N$

9: $\sigma_i^{MC} = \sqrt{(Som_2/N - (Som/N))^2/(N-1)}$

After demonstrating how to obtain a realization of a random variable, $W_{Y_i}$, following the Algorithms 1, 2, 3, 4 and 5, as well as how to estimate a mass fraction $Y_i$ using Monte Carlo according to Algorithm 6, let us now delve to a better interpretation of the different realizations outcomes of the random variables. In fact, the time evolution of the system and its changes over time (the conversion from A to B and from B to A or C) primarily stem from the non-zero initial conditions that have propagated over time. Thus, from a statistical standpoint, as soon as the observation time differs from the initial time, and since the algorithms are designed backward in time, it is crucial that the realizations of the random variables for a large number of iterations accurately reflect what is happening physically. In other words, they must effectively allow us to trace back to the origin of these evolutions of the system composed of A, B, and C, which are nothing but the non-zero initial conditions. This entails finding these initial conditions among the outcomes of the random variables employed to estimate the mass fractions of the various species. This is necessary to calculate a statistical mean that accurately represents the various contributions to the solution at a given observation time.

The main objective of this study is to assess the validity of our proposed expectation-based formulations of the presented kinetic models across various scenarios. These scenarios encompass situations with no time scale disparities between the different involved reactions, characterized as non-stiff, as well as stiff cases where there are disparities in time scales between the reactions.

In the stiff scenario, we delve into two distinct sub-cases. These sub-cases are characterized as follows: the quasi-steady-state (QSS) stiff situation [49], where the intermediate species B is formed slowly and consumed rapidly, and the partial equilibrium (PE) stiff situation [12, 50], where B is rapidly formed but consumed slowly.

The quasi-steady-state assumption assumes that, on the time scale of interest, the instantaneous rates of change for some intermediate species are approximately equal to zero. On the other hand, the partial equilibrium assumption implies that some fast reactions are always in equilibrium.

In our context, these approximations serve to describe specific arrangements of temporal scale disparities among various reactions. Their purpose is not to directly approximate the mass fraction of any involved species, but to derive approximate forms, allowing to determine more suitable density probabilities. We aim not to eliminate stiffness at these temporal scales, but rather to retain it and verify whether our method effectively captures the dynamics of the different involved species. Furthermore, we seek to understand how our algorithms behave under
<figure><img src="image_5.png" /><figcaption>Fig. 3. Graphical representation of backward-in-time random walks on the equation for estimating the mass fraction of C at $t_{obs}$ in the irreversible case. All paths start from species C at $t_{obs}$ (the orange dot in the figure). For all paths, random times are sampled according to the same probability density function $\pi_C$. At these new times, the mass fraction of species B needs to be evaluated. According to double randomization, a single realization is sampled from the B expectation-based model. Graphically, this results in a walk from C to B (dashed lines). Another random time is then sampled according to $\pi_B$. If the initial time is reached, the contribution of the initial condition $Y_{B,0}$ is taken into account in the Monte Carlo weight (represented by the red path ending in the red square). If the initial condition is not reached, the mass fraction of A needs to be evaluated. Following the same procedure from C to B, a walk from B to A is performed. Then, a random time is sampled according to $\pi_A$, which has two possible outcomes: a contribution from the initial condition $Y_{A,0}$ (the green path ending in the green square) or the null source term of $Y_A$ (the blue path ending in the blue star).</figcaption></figure>

extreme conditions, particularly when stiffness arises from coupling between different mass fractions whose dynamics are not of the same temporal order of magnitude.

To evaluate the accuracy of these expectation-based formulations across the various situations, we suggest employing the relative error criterion, defined as follows:

$$ \epsilon_i^{MC} = \begin{cases} \sigma_i^{MC}/Y_i^{MC}, & Y_i^{MC} \neq 0 \\ 0, & Y_i^{MC} = 0 \end{cases} \qquad (22) $$

This error can be either calculated numerically after each Monte Carlo run or determined analytically as:

$$ \epsilon_i^{MC} = \epsilon_i(t)/\sqrt{N} \qquad (23) $$

with,

$$ \epsilon_i(t) = \begin{cases} \sigma(W_{Y_i}(T))/\mathbb{E}(W_{Y_i}(T)), & \mathbb{E}(W_{Y_i}(T)) \neq 0 \\ 0, & \mathbb{E}(W_{Y_i}(T)) = 0 \end{cases} \qquad (24) $$

A well-converged Monte Carlo calculation is one where the relative error does not exceed 5%. Beyond this value, the error bars accompanying the Monte Carlo mean estimation become difficult to interpret physically [40].

This relative error will be calculated alongside the Monte Carlo estimation of the three mass fractions. It provides a convenient measure of statistical precision, expressing it as a fractional result relative to the estimated mean. Secondly, for worst-case scenarios, a more appropriate expectation-based formulation will be proposed utilizing importance sampling with variance reduction strategies [40] to determine the mathematical form of the optimal probability density.

In the irreversible case, the convergence of the Monte Carlo method will be assessed using the analytical forms of the relative error, as shown in Eq. (23). This is facilitated by the simplicity of the analytical expressions of the standard deviation and the expectation, allowing for the determination of the relative error. However, for the reversible case, the analytical forms of the standard deviation and the expectation are complex, thus we will rely solely on numerically determined relative errors after each Monte Carlo calculation following the Eq. (22).

### 3. Results and discussion

We will apply our Monte Carlo method for both irreversible and reversible scenarios, employing a varied selection of rate constant values. These values are chosen to reflect different timescale differences among the involved reactions. The maximum timescale difference considered is on the order of $10^8$ between the fastest and slowest reactions. Notably, beyond this value of timescale disparities, we observe similar patterns in the behavior of mass fractions and the convergence of our Monte Carlo calculations.

#### 3.1. Irreversible case

The evolution of mass fractions of the three species A, B, and C is computed using Monte Carlo simulation for a fixed number of samples, $N = 10^6$. Therefore, each result provided with an error bar relies on a Monte Carlo run and is estimated independently of the others, as was mentioned before. The results will be displayed and discussed for both non-stiff cases (where $k_1$ and $k_2$ are of similar magnitudes) and stiff cases (where $k_1$ and $k_2$ are not of similar magnitudes), and they will be compared with analytical solutions (see Supplementary Material B). Additionally, the relative errors will be plotted for the different mass fractions estimated by Monte Carlo considering the different cases. These relative errors are determined analytically according to Eq. (23).

##### 3.1.1. Non-stiff cases

In the non-stiff scenario, we investigate two distinct cases, each featuring a different analytical solution, resulting in varying forms of relative errors.

A graphical representation of the backward-in-time random walk in path space, for the irreversible and non-stiff case, is shown in Fig. 3. This illustration highlights the sampling process of the proposed methodology by depicting how the different random variables are sampled to reach the well-defined outcomes of the chemical system when tracing paths backward in time–namely, the initial mass fraction of species A, its null source term, or the initial mass fraction of species B.

When $k_1$ equals $k_2$, the Monte Carlo simulation demonstrates a qualitatively similar behavior to the analytical solutions for the mass fractions of the three species A, B and C, as depicted in Fig. 4. Similarly,
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Evolution of the mass fractions of the three species A, B, and C as a function of the dimensionless time, $t_{dim} = tk_1$, for the non-stiff where $k_1 = k_2$, given by Monte Carlo with the associated error bars (symbols) and the analytical solution (solid lines).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Evolution of the mass fractions of the three species A, B, and C as a function of the dimensionless time $t_{dim} = tk_2$, for the non-stiff where $k_1$ and $k_2$ are of the same order of magnitude (with $k_1 = 5s^{-1}$ and $k_2 = 1s^{-1}$), given by Monte Carlo with the associated error bars (symbols) and the analytical solution (solid lines).</figcaption></figure>

when $k_1$ is of the same order of magnitude as $k_2$ but differs from it, as shown in Fig. 5.

The relative error related to the estimation of the mass fraction species A depends solely on $k_1$. According to Figs. 4 and 5, the temporal evolution range of the mass fraction of A until reaching the steady state is approximately ~7/$k_1$. Consequently, within the relevant this temporal domain, the relative error peaks at around ~3.31%, occurring at $t = 7/k_1$ (see Supplementary Material C).

For species B and C, the temporal domain of interest lies between 0 and 7/$k_2$. Therefore, we plotted the relative errors $\epsilon_B^{MC}$ and $\epsilon_C^{MC}$ against a dimensionless time defined by $t_{dim} = k_2t$. Within this range of time, for the first non-stiff case, we observe that $\epsilon_B^{MC}$ does not exceed 0.5%, and $\epsilon_C^{MC}$ remains below 0.1%, and for the second non-stiff case, $\epsilon_B^{MC}$ does not exceed 4.6%, and the $\epsilon_C^{MC}$ remains below 0.14% as shown in Fig. 6. These results indicate that our Monte Carlo algorithm converges effectively within the relevant time domains, as the relative errors do not exceed 5%.

### 3.1.2. Stiff cases

For the stiff cases, we have two scenarios which will be treated separately. The case where species B is formed slowly and consumed rapidly ($k_1 \ll k_2$), which can be asymptotically approximated by the quasi-steady state (QSS) assumption, and the case corresponding to the rapid formation of species B and very slow consumption ($k_2 \ll k_1$), which can be asymptotically approximated by the partial equilibrium (PE) assumption.

#### 3.1.2.1. QSS situation.
In the scenario where species B is formed slowly and consumed rapidly ($k_1 \ll k_2$), corresponding to the QSS situation, only the extreme case where the ratio between $k_1$ and $k_2$ equals $10^{-8}$ was illustrated for the time evolution of the mass fractions of the three species as displayed in Fig. 7. According to this figure, Monte Carlo provides both qualitatively and quantitatively the same results as the analytical solutions for the mass fractions of all three species. It illustrates also, that in this QSS stiff case, the mass fraction of species B remains in quasi-steady state, nearly equal to its initial condition.

For species A, the corresponding relative error, within the time domain between 0 and 7/$k_1$ remains below 3.31% (see Supplementary
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Relative error calculated analytically according to Eq. (23) considering N = 10⁶ samples, as a function of dimensionless time, defined as $t_{dim} = tk_2$, for species B and C, in the two non-stiff cases where $k_1/k_2 = 1$ and $k_1/k_2 = 5$.</figcaption></figure>

Material C). For species B and C, the relative errors were plotted against the dimensionless time, $t_{dim} = tk_1$, following the rate constant of the slowest reaction, while considering different time scales disparities between the rate limiting step (A→B) and the fast reaction (B→C) as seen in Fig. 8. Remarkably, according to Fig. 8, the relative error of species B approaches the value of 0.01% for relatively small (one order of magnitude) disparities in time scales, and it remains close to 0 as the disparity in time scales becomes increasingly significant. Regarding the relative error of species C, it is notable that across all ratios, the relative error stabilizes around 0.1%.

In the QSS stiff case, our Monte Carlo method converges well, as the relative errors related to the estimations of the mass fractions of the three species do not exceed 5% in their corresponding transient temporal domains. This can be explained by the fact that, even though the disparities in time scales between the reactions (A→B) and (B→C) (the former being slow and the latter fast), the probability density function used to estimate the mass fraction of B is parametrized by the rate constant, $k_2$, which represents the dynamics of the fast reaction. Thus, it is possible to capture the dynamics of the slowest reaction when encountered using this probability density function, especially through double randomization that represents the coupling between the mass

fractions of A and B. In our scenario, the quantity we aim to estimate via Monte Carlo evolves rapidly compared to its source which dynamics are the dynamics of the rate limiting step. However, since B cannot vary in time without being generated by A, especially when the initial mass fraction of B is zero, the progression in time with realization of the random variable, $W_{Y_B}$, will be primarily governed by the rate limiting step dynamics, which is simply the dynamics of the consumption of A.

Our Monte Carlo algorithms are conceptualized to operate backward in time. To gain a deeper understanding of the algorithm's behavior under the quasi-steady-state (QSS) stiff case. Let us delve into the scenario of a single realization of the random variable, $W_{Y_A}$, according to Algorithm 2: Initially, a random time is sampled following $k_2 \exp(-k_2\tau)$. If the sampled time exceeds the observation time, $t-t_0$, we stop at the initial time and retain the initial condition of $Y_B$. However, if the sampled time is less than $t-t_0$, a double randomization process unfolds to determine a realization of the source term of the mass fraction, $Y_B$, corresponding to a realization of the random variable, $W_{Y_A}$, at $t - \tau$ up to a constant. As the probability density, $k_2 \exp(-k_2\tau)$, is parametrized by the rate constant, $k_2$, governing the fastest reaction, it leads to a slow random walk backward in time. Due to the slow consumption rate of species A, there is a higher likelihood of observing
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Evolution of the mass fractions of the three species A, B, and C as a function of the dimensionless time, $t_{\text{dim}} = tk_1$, for the QSS stiff case, characterized by an eight-order-of-magnitude separation in time scales between slow and fast reactions (with $k_1 = 1 \text{s}^{-1}$ and $k_2 = 10^8 \text{s}^{-1}$), given by Monte Carlo with the associated error bars (symbols) and the analytical solution (solid lines).</figcaption></figure>

it at $t - \tau$ while still capturing its consumption dynamics. By "capturing its consumption dynamics", we mean that among the possible outcomes of the random variable, $W_{Y_A}$, we will encounter the initial condition, $Y_{A,0}$, and not just the null source term (according to Green's formalism), as shown in Fig. 9. It is crucial to bear in mind that the coupling between $Y_A$ and $Y_B$ occurs through double randomization, which bases all reasoning on random variables rather than expectations. This makes understanding the algorithm's behavior significantly more complex, but feasible if we have an understanding based on averages of what may happen at any given time.

Therefore, under quasi-steady-state (QSS) stiff conditions, the employed density probabilities enable the construction of random variables that trace back to the non-zero initial condition, whose propagation over time is the origin of the dynamics of the mass fraction of B at a given time t.

A parallel can also be drawn between the integration step in deterministic solvers and the probability density in probabilistic solvers. An integration step that follows the fastest dynamics can also capture the slower dynamics without overlooking them during the resolution process [15]; the same principle applies to probability density here. However, this is not the case for the stiff PE scenario, as we will see in the following Section 3.1.2.2.

### 3.1.2.2. PE situation.
In the PE situation where species B is consumed slowly and formed rapidly ($k_1 \gg k_2$), only the extreme case where the ratio between $k_1$ and $k_2$ equals $10^8$ was illustrated for the time evolution of the mass fractions of the three species, as presented in Fig. 10. According to this figure, Monte Carlo provides both qualitatively and quantitatively the same results as the analytical solution for the mass fraction of species A. However, for species B and C, Monte Carlo simulations yield predominantly zeros with significant error bars, indicating the occurrence of rare events [51] with substantial weights. The rare events here are some outcomes among the possible ones that a realization of the random variables $W_{Y_B}$ or $W_{Y_C}$ can take.

For species A, the corresponding relative error within the time domain between 0 and $7/k_1$ remains below 3.31% (see Supplementary Material C). For species B and C, the relative errors were plotted this time against the dimensionless time, $t_{\text{dim}} = tk_2$, as $k_2$ is the rate constant of the rate limiting step, considering different orders of magnitude between $k_1$ and $k_2$ ranging from $10$ to $10^8$ as shown in Fig. 11. Analyzing the behavior of the relative error related to the mass fraction of species B across the time interval from 0 to $7/k_2$, as depicted on a logarithmic

scale (see Fig. 11), reveals that $\epsilon_B^{\text{MC}}$ remains below 20% for relatively small disparities in time (one to two orders of magnitude). Specifically, for the ratio of 10 between $k_1$ and $k_2$, $\epsilon_B^{\text{MC}}$ does not exceed 5% in the transient time domain. Over two orders of magnitude, $\epsilon_B^{\text{MC}}$ remains below 5% between 0 and $4/k_2$. In such cases, no special treatment is required; simply increasing the number of samples $N$ by a factor of 20 will reduce the relative error by $1/\sqrt{20} \approx 0.223$, resulting in a relative error of less than 5%. However, for larger time scale disparities, $\epsilon_B^{\text{MC}}$ exceeds by far 100%, indicating the occurrence of rare events with significant weights. For species C, no analysis can be drawn from the graph displayed in Fig. 11. We observe that for different ratios between $k_1$ and $k_2$, the relative errors do not exceed 0.15% in the transient domain. Nevertheless, in practice, the convergence of the Monte Carlo calculation for estimating the mass fraction of C depends on the convergence of the mass fraction estimation of B, since the random variable, $W_{Y_C}$, linearly depends on the random variable, $W_{Y_B}$. If rare events occur while estimating the mass fraction of B, they will also occur while estimating the mass fraction of C, but the analytical error does not allow this to be quantified, contrary to what is observed in Fig. 10.

We will analyze the behavior of our algorithm for different situations related to the observation time, $t$, at which we wish to estimate the mass fraction, $Y_B$. We want to emphasize that our analysis of the algorithm's behavior under partial equilibrium (PE) conditions relies on average trends, even though the behavior of random walks ultimately remain stochastic! For this purpose, we will define three time intervals: the first one $I_1 = [0, t_1]$, the second one $I_2 = [t_1, t_2]$, and the third one $I_3 = [t_2, +\infty[$, with $t_1 = 7/k_1$ and $t_2 = 7/k_2$.

For $t \in I_1$, i.e., when we are close to the starting time of the reaction, and within the interval where A is consumed to give B, the reaction scheme is under kinetic control [52] and far from equilibrium. The reaction scheme presented in Fig. 1 can be reduced to the following equivalent scheme presented in Fig. 12.

The rate limiting step, converting B to C has not yet come into play. Utilizing an exponential distribution, $\pi_B(\tau) = k_2 \exp(-k_2\tau)$, defined between 0 and $+\infty$ with a mean $1/k_2$ will result in sampling random times around $1/k_2$, thus falling outside the interval $I_1$. Consequently, the random variable, $W_{Y_B}$, has as an outcome only the initial condition, $Y_{B,0}$, as shown in Table 1, so it fails to capture the fast dynamics, which are the consumption of A over time, which leads to the formation of B. Therefore, this used probability density is inadequate for capturing the dynamics of the fast reaction (A→B). Even with a large number of
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Relative error calculated analytically according to Eq. (23) considering $N = 10^6$ samples, as a function of dimensionless time, defined as $t_{dim} = tk_1$, for species B and C, in the QSS stiff case. The analysis spans several orders of magnitude separation in time scales between slow and fast reactions, ranging from $k_1/k_2 = 10^{-1}$ to $k_1/k_2 = 10^{-8}$.</figcaption></figure>

samples $N = 10^6$, the outcome returned by Monte Carlo resembles the initial condition, $Y_{B,0}$, as if nothing occurred between A and B as seen in Fig. 10.

Likewise, for the other time intervals, $I_2$ and $I_3$, the observation times at which we want to estimate $Y_B$ are far from the time span during which A is consumed. For times within $I_2$, the random variable, $W_{Y_B}$, may exhibit two possible outcomes as the algorithm is written backward in time, walking randomly backward in time with a probability density parametrized by $k_2$. This walk in time may quickly revert to the initial condition of $Y_B$, or it will necessitate to continue by performing a double randomization to switch to the expectation-based formulation to estimate the mass fraction of A, in order to obtain a realization of the source term, $S_B$, at a previous time, $t-\tau$. The observation time, $t-\tau$, for the random variable, $W_{Y_A}$, will likely remain considerably distant from the time span during which the consumption of A occurs. Consequently, the outcomes of the random variable $W_{Y_A}$ will predominantly consist of the null source term of $Y_A$ as it is illustrated in Fig. 13. From a Green's perspective, only the contribution of the source term is considered, implying deterministically that A never exist, which is not the case.

According to Table 1, for an observation time $t = 1/k_2$, considering $10^6$ realizations, the random variable $W_{Y_B}$ is more likely to have as an

<caption>**Table 1**</caption>

The different outcomes of the random variable $W_{Y_B}$ for various observation times for $10^6$ realizations under PE conditions for the case $k_1/k_2 = 10^8$.

<table><thead><tr><th>Observation time</th><th>Y<sub>B,0</sub></th><th>S<sub>A</sub></th><th>Y<sub>A,0</sub></th></tr></thead><tbody><tr><td>t = 1/k<sub>1</sub></td><td>10<sup>6</sup></td><td>0</td><td>0</td></tr><tr><td>t = 1/k<sub>2</sub></td><td>367752</td><td>632248</td><td>0</td></tr><tr><td>t = 5/k<sub>2</sub></td><td>6701</td><td>993299</td><td>0</td></tr></tbody></table>

outcome the source term, $S_A$, followed by the initial condition of $Y_B$ and never the initial condition of $Y_A$. The longer the observation time, the more the null source term, $S_A$, is retained as an outcome, and the less likely it is to reach the initial condition of $Y_B$ as shown in Table 1. This is particularly the case for observation times in $I_3$, which are significantly distant compared to the time span in which A is consumed.

Thus, within algorithmic contexts, the rare event manifests as the outcome, $Y_{A,0}$, that is consistently disregard. However, under PE stiff conditions, the employed density probabilities do not enable the construction of random variables that trace back to the non-zero initial condition, whose propagation over time is the origin of the dynamics of the mass fraction of B at a given time t. This translates to an inability to
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Graphical representation of backward-in-time random walks on the equation for estimating the mass fraction of C at $t_{\text{obs}}$ in the irreversible QSS stiff case, where $k_2 \gg k_1$. All paths start from species C at $t_{\text{obs}}$ (the orange dot in the figure). The presented paths illustrate the behavior of the majority of paths under the QSS stiff case. Specifically, as $k_2 \gg k_1$, they predominantly end either at the null source term of $Y_A$ (blue star) or at the initial condition of $Y_A$ (red and green squares). They rarely, if ever, end at the initial mass fraction of $Y_B$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Evolution of the mass fractions of the three species A, B, and C as a function of the dimensionless time, $t_{\text{dim}} = tk_2$, for the PE stiff case, characterized by an eight-order-of-magnitude separation in time scales between slow and fast reactions (with $k_1 = 10^8$ s⁻¹ and $k_2 = 1$ s⁻¹), given by Monte Carlo with the associated error bars (symbols) and the analytical solution (solid lines).</figcaption></figure>

capture the dynamics of the kinetic regime wherein A undergoes rapid consumption to yield B.

Integral Formulation Monte Carlo method fails capturing the slow dynamics due to the reasons we elucidated earlier, which give rise to the occurrence of rare events. This limitation is also observed in other stochastic methods, such as kinetic Monte Carlo (kMC).

Indeed, kMC is a stochastic method where sampling is directly based on events, such as reactions or conversions from one species to another, each characterized by a probability in the form of $k_j^{\text{MC}} / \sum_j k_j^{\text{MC}}$, where $k_j^{\text{MC}}$ represents the rate constant of reaction j [17]. As a result, rare events manifest as infrequent occurrences of certain conversion processes or reactions.

In the literature, the same challenges arise when modeling homogeneous batch series reactions similar to the reaction schemes illustrated in Figs. 1 and 2 under Partial Equilibrium (PE) conditions, utilizing kinetic Monte Carlo (kMC) [18].

Essentially, the conventional kMC algorithm primarily samples the rapid reactions and seldom samples the slow reaction events under

PE conditions, hence encountering difficulty in capturing the slow dynamics.

Ideally, we should construct a probability density that enables us to accurately capture the consumption of A over time, ensuring the formation of B is adequately traced. To do so, we have to use the importance sampling technique [53] to effectively capture the consumption of A. Importance sampling requires that the new probability density follows the dynamics of A, what is called mathematically the target. In irreversible cases, where the expected value is known, the technique of zero-variance principle [54] can be used to search for the mathematical form of the ideal probability density function. The zero-variance principle allows finding the optimal probability density function that yields the exact solution in a single sampling by reducing or even suppressing the related variance [54]. Therefore, to obtain the exact solution in a single sampling, we must replace the mass fraction of A with its analytical solution.

We begin by restricting the integration domain from 0 to +∞ to the interval between $t_0$ and $t$. This transformation allows us to sample random times in the temporal domain of interest. The expectation-based
<figure><img src="image_13.png" /><figcaption>Fig. 11. Relative error calculated analytically according to Eq. (23) considering N = 10⁶ samples, as a function of dimensionless time, defined as $t_{dim} = tk_2$, for species B and C, in the PE stiff case. The analysis spans several orders of magnitude separation in time scales between slow and fast reactions, ranging from $k_1/k_2 = 10^1$ to $k_1/k_2 = 10^8$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Equivalent reaction scheme at kinetic control.</figcaption></figure>

formulation of species B becomes:

$$Y_B(t) = Y_{B,0} \exp(-k_2(t-t_0)) + \int_0^{t-t_0} [k_1 Y_A(t-\tau) \exp(-k_2\tau) d\tau] \quad (25)$$

By replacing $Y_A(t - \tau)$ with its analytical form and changing the integration bounds to $t_0$ and $t$, we find:

$$Y_B(t) = Y_{B,0} \exp(-k_2(t-t_0)) + \int_{t_0}^{t} [k_1 Y_{A,0} \exp(-k_1(t-\tau-t_0)) \exp(-k_2(\tau-t_0)) d\tau] \quad (26)$$

Now, we need to introduce a probability density defined between $t_0$ and $t$ that allows us to obtain the exact value of the mass fraction of B at time $t$ with no variance. This probability density is the integrand in Eq. (26) normalized between $t_0$ and $t$. Thus, the optimal probability

density function is:

$$\pi_B^{ZV}(\tau) = \frac{k_1 Y_{A,0} \exp(-k_1(t-\tau-t_0)) \exp(-k_2(\tau-t_0))}{\left( \int_{t_0}^{t} k_1 Y_{A,0} \exp(-k_1(t-\tau-t_0)) \exp(-k_2(\tau-t_0)) d\tau \right)} \quad (27)$$

Therefore,

$$\pi_B^{ZV}(\tau) = \frac{(k_1 - k_2) \exp(-(k_1 - k_2)\tau)}{\exp(-(k_1 - k_2)t) - \exp(-(k_1 - k_2)t_0)} \quad (28)$$

Substituting the new probability density into the expectation-based formulation of B, we find:

$$Y_B(t) = \int_{t_0}^{t} \pi_B^{ZV}(\tau) d\tau \left[ \frac{k_1 \exp(-k_2(\tau-t_0))}{\pi_B^{ZV}(\tau)} Y_A(t-t_0-\tau) + Y_{B,0} \exp(-k_2(t-t_0)) \right] \quad (29)$$

The algorithmic approach for determining the mass fraction of B is now a different **Algorithm 7**. Firstly, all randomly sampled times now
<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Graphical representation of backward-in-time random walks on the equation for estimating the mass fraction of C at $t_{obs}$ in the irreversible PE stiff case. All paths start from species C at $t_{obs}$ (the orange dot in the figure). The presented paths illustrate the behavior of the majority of paths under the PE stiff case. Specifically, as $k_1 \gg k_2$, they predominantly end either at the null source term of $Y_A$ (blue and green stars) or at the initial condition of $Y_B$ (red square). They rarely, if ever, end at the initial mass fraction of $Y_A$.</figcaption></figure>

fall within the interval from $t_0$ to $t$, preventing the possibility of sampling random times greater than the observation time. Consequently, the number of outcomes is now limited to two possibilities: the initial condition of $Y_A$ or it is related source term, which is equal to 0. This reduction in outcomes effectively reduces the variance. As a result, we always undergo a double randomization process to obtain a realization of $W_{Y_A}$ (see Algorithm 7).

Secondly, in this new algorithm, the initial condition accumulates along with the realization of the random variable, $W_{Y_B}$, while also considering its exponential decay (see Algorithm 7). Therefore, for each sampling of $\tau$, the accumulated initial condition is multiplied by $\exp(-k_2(t - t_0))$ to incorporate the exponential decay. Returning to a Green formalism interpretation, the probabilized contributions are solely the contributions of the source. By selecting an optimal probability density following the dynamics of this source, it becomes possible to reduce or even eliminate the variance.

**Algorithm 7** $W_{Y_B}(\tau)$ sampling algorithm

1: **Sample** a random time $\tau$ following the probability density $\pi_B^{ZV}(\tau)$.

2: **Sample** $W_{Y_A}(t - \tau)$ according to Algorithm 1.

3: **return** $\frac{k_1 \exp(-k_2(\tau - t_0))}{\pi_B^{ZV}(\tau)} w_{Y_A}(t - t_0 - \tau) + Y_{B,0} \exp(-k_2(t - t_0))$.

We have applied the importance sampling technique with the zero variance principle to search for the optimal density function for the expectation-based formulation of the mass fraction of C. Following the same zero variance procedure as described earlier, we have derived the optimal density function, expressed as follows:

$$ \pi_C^{ZV}(\tau) = \frac{\exp(-k_2\tau) - \exp(-k_1\tau)}{\frac{\exp(-k_2t_0) - \exp(-k_2t)}{k_2} - \frac{\exp(-k_1t_0) - \exp(-k_1t)}{k_1}} \quad (30) $$

For the expectation-based formulation of the mass fraction of species C, it is sufficient to replace the uniform distribution defined between $t_0$ and $t$ that we initially proposed in Eq. (13) with the new probability density $\pi_C^{ZV}$. We follow exactly the same Algorithm 3 to obtain a realization of the random variable, $W_{Y_C}$, but with a new sampling distribution this time.

Applying the new density function to estimate $Y_B$ and $Y_C$, we observe that for a significant time disparity, as illustrated in Fig. 14 (eight orders of magnitude), our Monte Carlo simulations yield both quantitative and qualitative results that are similar to the analytical solutions for the mass fractions of B and C, with minimal variances.

## 3.2. Reversible case

Now, we will explore the reversible case, maintaining the same sub-cases as for the irreversible case, except for the non-stiff scenario where we will focus on a situation where all non-zero values of the reaction rates ($k_1, k_{-1}$, and $k_2$) yield identical analytical solutions for the mass fractions of species A, B, and C.

The classification of cases and sub-cases will be based on a direct comparison between the rate constants $k_1$ and $k_s = k_{-1} + k_2$, as the probability densities proposed in Eqs. (15) and (16) are respectively parametrized by $k_1$ and $k_s$. We will continue to examine the time disparities by comparing these two constants, without paying any particular attention to $k_{-1}$ or $k_2$. Although this comparison may obscure several sub-cases, such as in a stiff scenario where $k_s$ is not of the same order as $k_1$ but not necessarily $k_2$ or $k_{-1}$, it is not a significant concern in our proposed method. What matters is solely the order of magnitude between $k_1$ and $k_s$, as ultimately the coupling is achieved through double randomization, and the realizations of the random variables used depend only on probability densities parametrized by $k_1$ and $k_s$. The same comparison between $k_1$ and $k_s$ has been conducted to assess the conditions of validity of both the Partial Equilibrium (PE) assumption and the Quasi-Steady-State (QSS) assumption within the context of the same reaction scheme presented in Fig. 2, as elaborated in [55,56]. Therefore, for a given order of magnitude between $k_1$ and $k_s$, if Monte Carlo adequately handles this temporal scale disparity, regardless of the sub-cases that may cause such a discrepancy, they will be properly addressed by our Monte Carlo method in the end. This holds true for the inverse case as well. The results will be displayed and discussed for both non-stiff case and stiff cases, and they will be compared with analytical solutions (see Supplementary Material B). Additionally, the relative errors will be plotted for the different mass fractions estimated by Monte Carlo considering the different cases. These relative errors are determined numerically according to Eq. (23).

### 3.2.1. Non-stiff case

In the non-stiff and reversible case, a graphical representation of the backward-in-time random walk in path space is provided in Fig. 15. This illustration highlights the sampling process of the proposed methodology by showing how the random variables are sampled to recover the well-defined outcomes of the chemical system when tracing backward in time–namely, the initial mass fraction of species A or species B. In contrast to the irreversible case, reversibility is reflected
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Evolution of the mass fractions of the three species A, B, and C as a function of the dimensionless time, $t_{dim} = tk_2$, for the PE stiff case after using the importance sampling and zero-variance principle, characterized by an eight-order-of-magnitude separation in time scales between slow and fast reactions (with $k_1 = 10^8$ s⁻¹ and $k_2 = 1$ s⁻¹), given by Monte Carlo with the associated error bars (symbols) and the analytical solution (solid lines).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Graphical representation of backward-in-time random walks on the equation for estimating the mass fraction of C at $t_{obs}$ in the reversible case. All paths start from species C at $t_{obs}$ (the orange dot in the figure). For all paths, random times are sampled according to the same probability density function $\pi_C$. At these new times, the mass fraction of species B needs to be evaluated. According to double randomization, a single realization is sampled from the B expectation-based model. Graphically, this results in a walk from C to B (dashed lines). Then, another random time is sampled according to $\pi_B$. If the initial time is reached, the contribution of the initial condition $Y_{B,0}$ is accounted for in the Monte Carlo weight (represented by the red path ending on the red square). If the initial condition is not reached, the mass fraction of A needs to be evaluated. Following the same procedure from C to B, a walk from B to A is performed. Then, a random time is sampled according to $\pi_A$, which has two possible outcomes: a contribution from the initial condition $Y_{A,0}$ (the blue path ending on the blue square) or returning to B due to reversibility (the green path ending on the green square), which explains the back-and-forth between A and B.</figcaption></figure>

here by the possibility of moving backward and forward between species A and B, capturing the bidirectional nature of the reaction.

The Monte Carlo simulation demonstrates a qualitatively similar behavior to the analytical solution for the mass fraction of the three species A, B and C, as depicted in Fig. 16.

The numerical relative error associated with the mass fraction of A increases until reaching 0.3% at $t = 19/k_1$, as illustrated in Fig. 17. The related numerical relative errors to the mass fractions of species B and C exhibit the same qualitative behavior. The two errors peak around an observation time close to the initial time, with a slight overshoot of the relative error associated with the mass fraction of C exceeding 5%. Then, both errors decrease near to 0 for the rest of the observation times.

These results indicate that our Monte Carlo algorithm converges effectively within the relevant time domains, as for the most observation times the relative errors do not exceed 5% for the three species A, B and C.

### 3.2.2. *Stiff cases*

As in the case of irreversible reactions, we will examine the stiffness by defining it through the intermediate species B standpoint. To do this, we will consider the two cases analyzed previously for irreversible reactions: the quasi-steady-state (QSS) situation, where ($k_1 \ll k_s$), and the partial-equilibrium (PE) situation, where ($k_1 \gg k_s$).

#### 3.2.2.1. QSS situation.
According to Fig. 18, the Monte Carlo yield for the mass fractions of the three species A, B and C are similar to the
<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Evolution of the mass fractions of the three species A, B, and C as a function of the dimensionless time, $t_{dim} = tk_1$, for the non-stiff where $k_1$ and $k_s$ are of the same order of magnitude (with $k_1 = 1 \text{ s}^{-1}$ and $k_s = 1 \text{ s}^{-1}$), given by Monte Carlo with the associated error bars (symbols) and the analytical solution (solid lines).</figcaption></figure>

<table><caption>**Table 2**
The different outcomes of the random variable $W_{Y_s}$ for various observation times for $10^6$ realizations under QSS conditions for the case $k_1/k_s = 10^{-8}$.</caption><thead><tr><th>Observation time</th><th>Y<sub>B,0</sub></th><th>Y<sub>A,0</sub></th></tr></thead><tbody><tr><td>t = 1/k<sub>s</sub></td><td>367750</td><td>632250</td></tr><tr><td>t = 10<sup>4</sup>/k<sub>s</sub></td><td>0</td><td>10<sup>6</sup></td></tr><tr><td>t = 1/k<sub>1</sub></td><td>0</td><td>10<sup>6</sup></td></tr></tbody></table>

analytical solutions for the case where an eight-order-of magnitude separation is considered between $k_1$ and $k_s$. We can clearly observe the quasi-stationary behavior of the mass fraction of B, with nearly negligible variations over time compared to its initial value.

The relative errors associated with the mass fraction estimations for the three species do not exceed the threshold of 5%, as shown in Fig. 19. This indicates that even in the case of a reversible reaction, Monte Carlo converges well when handling the QSS stiff case.

Let us imagine a backward random walk in time performed via double randomization, alternately by $W_{Y_B}$ according to $k_s \exp(-k_s \tau)$ and $W_{Y_A}$ according to $k_1 \exp(-k_1 \tau)$. Following Algorithm 5, the two outcomes of the random variable, $W_{Y_B}$, are $Y_{B,0}$ or $Y_{A,0}$ up to a constant factor. The first outcome, $Y_{B,0}$, is only reached by the random variable, $W_{Y_B}$, and the second outcome, $Y_{A,0}$, is only reached via $W_{Y_A}$ when performing double randomization. The dynamics of this walk are governed by the slowest process as its related probability density samples the greater random times, which is simply in this case $k_1 \exp(-k_1 \tau)$ (see Fig. 20), leading to a rapid return to the initial time. Thus, we can imagine that in the majority of cases, it is $W_{Y_A}$ that will first reach the outcome $Y_{A,0}$ to stop the random walk.

Certainly, the farther the observation time is from the initial time, the less likely it is to reach the initial condition of $Y_B$ as an outcome, thus leading to the occurrence of rare events for observation times with an order of magnitude greater than $1/k_s$. A simulation with $10^6$ realizations for various chosen observation times of the order of magnitude of the fastest reaction ($1/k_s$) and the slowest reaction ($1/k_1$), as well as an intermediate order of magnitude between the two ($10^4/k_s$), as shown in Table 2. The simulation precisely illustrates what we explained earlier. Moreover, the closer we get to the initial condition, the more likely we are to sample $Y_{B,0}$.

However, since $Y_{B,0}$ is equal to zero, it represents a non-significant rare event. Physically, the existence of species B does not arise from its initial condition, but rather from the consumption of A. This consumption is represented in terms of a random walk, with $Y_{A,0}$ being retained.

From the perspective of Green's formalism, it is the propagation in time of these initial conditions that yields B.

Mathematically, we can demonstrate how a null rare event does not influence the convergence of a Monte Carlo calculation. Let us take the definition of the standard deviation given by Eq. (21), and consider that out of $N$ realizations of the random variable, $W_{Y_i}$, which are denoted as $w_{Y_i}$, we have $M$ realizations which outcome is the rare event, where $M \ll N$ (because the occurrence of the rare event is rare!). We will denote the related weight to rare event as $w_{RE}$. We can then rewrite the standard deviation as:

$$ \sigma_i^{MC} = \left( \frac{1}{N-1} \left( \frac{1}{N} \left( (N-M)w_{Y_i}^2 + Mw_{RE}^2 \right) - \left( \frac{1}{N} \left( (N-M)w_{Y_i} + Mw_{RE} \right) \right)^2 \right) \right)^{0.5} \quad (31) $$

After simplification and considering a large number of realizations, $N$, we can approximate $N \approx N - M$, and we can neglect terms with $1/N^2$ compared to terms with $1/N$. Thus, we find:

$$ \sigma_i^{MC} \approx (M^2/N^2 w_{RE}^2 - 2M/N w_{Y_i} w_{RE})^{0.5} \quad (32) $$

So, if $w_{RE} \to 0$, then $\sigma_i^{MC} \to 0$. Therefore, a null rare event does not influence the Monte Carlo calculation, as from the standpoint of Green's formalism, this event does not contribute to the desired solution at a given time. However, if this event has a significant value, statistically, we will not be able to count its contribution, which will influence our Monte Carlo calculation and, consequently, increase its variance.

However, the analysis above remains valid only for the specific considered initial conditions, which involve no initial mass fraction of B. If this were not the case, the rare event would differ from zero, and from the perspective of A, QSS would result in partial equilibrium due to reversibility (A would be rapidly formed and slowly consumed). Nevertheless, to simplify the analysis, we confined ourselves to a scenario where the primary focus lies on B due to the chosen specific initial conditions. Consequently, for a configuration where the mass fraction of B differs from zero, we can apply the same analysis and solution developed for the mass fraction of B under PE conditions in the next Section 3.2.2.2 to also address the mass fraction of A (see Fig. 21).

### 3.2.2.2. PE situation.
Similarly to the irreversible case (see Section 3.1.2.2), Monte Carlo fails to capture the dynamics of the mass fraction of the three involved species, for the PE extreme situation,
<figure><img src="image_19.png" /><figcaption>**Fig. 17.** Relative error calculated numerically according to Eq. (22) considering $N = 10^6$ samples, as a function of dimensionless time, defined as $t_{dim} = tk_1$, for species A, B and C, in the non-stiff case where $k_1/k_s = 1$.</figcaption></figure>

considering a ratio of $10^8$ between $k_1$ and $k_s$, as shown in Fig. 21. For the three species A, B and C, according to Fig. 22, the proposed-expectation-based formulation converges well, with relative errors that do not exceed 5% for one to two orders of magnitude separation in timescales between the fast process (A→B) and the slow process (A ← B → C). Beyond these separations in timescales, the relative errors exceed 5%, reaching up to 100% or sometimes 0%, for the timescale disparity of $10^8$. However, a 0% error value does not necessarily indicate good convergence. When examining the results provided by Monte Carlo, where these errors are zero, these results also differ greatly from the analytical solution. This indicates the occurrence of rare events that have not been sampled with a sample size of $N = 10^6$. It is necessary to exceed by far this value to capture such events.

The behavior of the algorithm can be explained by the fact that the dynamics of the mass fraction of A is overlooked by the probability density $k_s \exp(-k_s \tau)$, used to estimate the mass fraction of B. Let us imagine a backward random walk in time performed via double randomization, alternately by $W_{Y_B}$ according to $k_s \exp(-k_s \tau)$ and $W_{Y_A}$ according to $k_1 \exp(-k_1 \tau)$. Following Algorithm 5, the two outcomes of the random variable $W_{Y_B}$ are $Y_{B,0}$ or $Y_{A,0}$ up to a constant factor. The

first outcome, $Y_{B,0}$, is only reached by the random variable $W_{Y_B}$, and the second outcome $Y_{A,0}$ is only reached via $W_{Y_A}$ when performing double randomization. The dynamics of this walk are governed by the slowest process as its related probability density samples the greater random times, which is simply in this case $k_s \exp(-k_s \tau)$, leading to a rapid return to the initial time. Thus, we can imagine that in the majority of cases, it is $W_{Y_B}$ that will first reach the outcome $Y_{B,0}$ to stop the random walk as shown in Fig. 23.

For observation times close to the initial time or under kinetic control, the scheme presented in Fig. 2 can also be reduced to the scheme presented in Fig. 12, since neither the reverse reaction nor the reaction converting B to C are involved yet. In light of this approximation, which facilitates our understanding, based on our analysis in the irreversible PE case (see Section 3.1.2.2), we can say that at the initial time, the only contribution to the solution of the mass fraction of B will be its initial condition as long as we remain under kinetic control. Consequently, the occurrence of rare events, which are the contribution of the initial condition of the mass fraction of A, as shown in Table 3 for an observation time equal to $1/k_1$ and for $10^6$ samples.
<figure><img src="image_20.png" /><figcaption>**Fig. 18.** Evolution of the mass fractions of the three species A, B, and C as a function of the dimensionless time $t_{dim} = tk_1$, for the QSS stiff case, characterized by an eight-order-of-magnitude separation in time scales between slow and fast reactions (with $k_1 = 1 \text{s}^{-1}$ and $k_s = 10^8 \text{s}^{-1}$), given by Monte Carlo with the associated error bars (symbols) and the analytical solution (solid lines).</figcaption></figure>

<table><caption>**Table 3**
The different outcomes of the random variable $W_{Y_s}$ for various observation times for $10^6$ realizations under PE conditions for the case $k_1/k_s = 10^8$.</caption><thead><tr><th>Observation time</th><th>Y<sub>B,0</sub></th><th>Y<sub>A,0</sub></th></tr></thead><tbody><tr><td>t = 1/k<sub>1</sub></td><td>99999</td><td>1</td></tr><tr><td>t = 1/k<sub>s</sub></td><td>10<sup>6</sup></td><td>0</td></tr><tr><td>t = 5/k<sub>s</sub></td><td>10<sup>6</sup></td><td>0</td></tr></tbody></table>

As we move away from kinetic control, i.e., for observation times far from the initial time, B will begin to be consumed, and thus the random variable $W_{Y_B}$ will not only have $Y_{B,0}$ as an outcome, but there will also be the contribution of the source term. This contribution is performed via double randomization. So, at this stage of the reaction, the backward random walks in time will become alternated between $W_{Y_B}$ and $W_{Y_A}$, which reflects well the reversibility of the conversion between A and B, but as we have already explained, this will be governed by $W_{Y_B}$, and thus we will always end up retaining only the initial condition of $Y_B$ in the end, as if A never existed. To illustrate this case, Table 3 shows the outcomes of the random variable, $W_{Y_B}$, for observation times in the kinetic control and far from it.

As a result, for the different situations related to the different observation times, it is almost never or rarely to attain the initial condition of $Y_A$ as an outcome, thereby failing to reach the non-zero initial condition.

In light of what has just been explained, we can understand why relative errors can sometimes even reach 100%. This value itself represents the extreme case of the occurrence of rare events, as an error reaches this value only when a single rare event occurs out of a certain number N of realizations of a given random variable $W_Y$. We can demonstrate this through the Eq. (32), as for our particular case, the event or outcome that is not rare is the initial condition of the mass fraction of B, which is itself zero, so the standard deviation reduces to:

$$ \sigma_i^{MC} \approx (w_{RE}^2/N^2)^{0.5} = Y_i^{MC} \qquad (33) $$

Hence the value of 100% for the relative error.

As in the irreversible case for the PE situation (see Section 3.1.2.2), to overcome the occurrence of rare events, we will have to construct a probability density that follows the dynamics of the mass fraction of A in this reversible case under PE conditions. Therefore, we are going to restrict our integration domain between $t_0$ and $t$ to avoid sampling random times far from the temporal domain of interest, and

consequently, the initial condition, $Y_{B,0}$, will be accumulated along the realization of the random variable, $W_{Y_B}$, as in the irreversible case.

Certainly, we have the analytical solution of A (see Supplementary Material B), but it is not very interesting to use it in this context as it will result in a probability density whose cumulative cannot be analytically inverted to sample random times, which will add a numerical bias. Additionally, we are in a perspective where we wish to propose a method that can work without knowing exactly the analytical solution but at least knowing the behavior of the function we wish to integrate through the approximations we can make in such conditions, specifically the conditions of partial equilibrium (PE). We wish to go through a physical understanding of the system's behavior and propose, in light of this understanding, an appropriate probability density.

For this purpose, we will examine how the mass fraction of A behaves under the kinetic control regime, i.e., close to the initial time, and the regime far from kinetic control, when we tend towards equilibrium or otherwise called thermodynamic control [52].

Under kinetic control, as we have already shown, the reaction scheme in the reversible case can be reduced to a reaction of A yielding B as shown in Fig. 12. Mathematically, the behavior of the mass fraction of A under this regime is similar to the behavior of the mass fraction of A in the irreversible case, and we can write:

$$ Y_A(t) \approx Y_{A,0} \exp(-k_1 t) \qquad (34) $$

So we can use the same probability density in Eq. (26), which we found using the zero-variance principle in the irreversible case under PE conditions, by simply replacing $k_2$ with $k_s$, which leads to finding:

$$ \pi_B^{KC}(\tau) = \frac{(k_1 - k_s) \exp((k_1 - k_s)\tau)}{\exp((k_1 - k_s)t) - \exp((k_1 - k_s)t_0)} \qquad (35) $$

Now, far from the kinetic control regime and as we tend towards thermodynamic control, the behavior of the reversible case becomes different from the irreversible case. As the reverse process and the reaction converting B to C come into play, a new probability density must be considered to capture the dynamics of the mass fraction of A in this regime. Let us delve into understanding how $Y_A$ behaves in this regime. At this point, A is slowly formed by B according to $k_{-1}$ and is rapidly consumed. Consequently, the changes in the mass fraction of A over time are minimal and can even be neglected, which is otherwise
<figure><img src="image_21.png" /><figcaption>Fig. 19. Relative error calculated numerically according to Eq. (22) considering N = 10⁶ samples, as a function of dimensionless time, defined as tdim = rk₁, for species A, B and C, in the QSS stiff case where. The analysis spans several orders of magnitude separation in time scales between slow and fast reactions, ranging from k₁/kₛ = 10⁻¹ to k₁/kₛ = 10⁻⁸.</figcaption></figure>

the manifestation of the partial equilibrium between A and B that has
been established. Mathematically, we can write:

$$
\frac{dY_A}{dt} \approx 0 \tag{36}
$$

Thus,

$$
-k_1 Y_A + k_{-1} Y_B \approx 0 \tag{37}
$$

Therefore,

$$
Y_A(t) \approx k_{-1}/k_1 Y_B(t) \tag{38}
$$

The Eq. (38) represents the consequence of the partial equilibrium approximation when we are far from kinetic control. By replacing this approximate form of the mass fraction of A in the kinetic equation of the mass fraction of B, we find:

$$
\frac{dY_B}{dt} \approx -k_2 Y_B \tag{39}
$$

This leads to:

$$
Y_B(t) \approx Y_{B,I} \exp(-k_2 t) \quad (40)
$$

with Y<sub>B,I</sub> being the mass fraction from the moment we are no longer
in kinetic control. As we want to determine only an approximate form
of the mass fractions, we can estimate that Y<sub>B,I</sub> ≈ Y<sub>A,0</sub> since from the
moment we are no longer in kinetic control, A is completely consumed
to give B and B has not yet started to be consumed.

From Eqs. (38) and (40), we have:

$$
Y_A(t) \approx k_{-1}/k_1 Y_{A,0} \exp(-k_2 t) \quad (41)
$$

At this stage, by replacing the source term of the expectation-based
formulation with the analytical approximate form of the mass fraction
of A proposed in Eq. (41), we can determine a density probability that
could accurately follow the dynamics of the mass fraction of A, far from
the kinetic control regime. By doing so, we find a density probability
under the following form:

$$
\pi_B^{PE}(\tau) = \frac{k_{-1} \exp(-k_2 t) \exp(-k_{-1} \tau)}{\exp(-k_2 t)(\exp(-k_{-1} t_0) - \exp(-k_{-1} t))} \quad (42)
$$

Therefore, it is of interest to construct a density probability that can
capture both, the behavior of the mass fraction of A at the early stages
<figure><img src="image_22.png" /><figcaption>Fig. 20. Graphical representation of backward-in-time random walks on the equation for estimating the mass fraction of C at $t_{obs}$ in the reversible QSS stiff case. All paths start from species C at $t_{obs}$ (the orange dot in the figure). The presented paths illustrate the behavior of the majority of paths under the QSS stiff case. Specifically, as, $k_x \gg k_1$, they predominantly end at the initial condition of $Y_A$ (red and blue squares). They rarely, if ever, end at the initial mass fraction of $Y_B$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 21. Evolution of the mass fractions of the three species A, B, and C as a function of the dimensionless time, $t_{dim} = tk_1$, for the PE stiff case, characterized by an eight-order-of-magnitude separation in time scales between slow and fast reactions (with $k_1 = 10^8$ s⁻¹ and $k_x = 1$ s⁻¹), given by Monte Carlo with the associated error bars (symbols) and the analytical solution (solid lines).</figcaption></figure>

of the reaction and the moment when the PE starts manifesting near the thermodynamically controlled regime. For this purpose, we propose a probability density under the following form:

$$ \pi_B^{IS}(\tau) = \frac{k_1 \exp(-k_1 t) \exp((k_1 - k_s)\tau) + k_{-1} \exp(-k_{-1} t) \exp(-k_{-1} \tau)}{F_N(t)} \quad (43) $$

with, $F_N(t) = \frac{k_1}{(k_1 - k_s)} \exp(-k_1 t) (\exp((k_1 - k_s)t) - \exp((k_1 - k_s)t_0)) + \exp(-k_2 t) (\exp(-k_{-1}t_0) - \exp(-k_{-1}t))$ is the integrand of the function $k_1 e^{-k_1 t} \exp((k_1 - k_s)\tau) + k_{-1} \exp(-k_{-1}t) \exp(-k_{-1} \tau)$ between $t_0$ and $t$.

We observe that the numerator of the probability density function proposed in Eq. (43) is composed of the sum of the functions that have been used to deduce the probability densities separately under the kinetic control regime and far from it. Then this sum is normalized by also the sum of the integrands of the probability densities given by Eq. (35) and (42) between $t$ and $t_0$. We can settle for this probability density, but its cumulative cannot be inverted analytically to sample random times.

However, by further developing the probability density given by Eq. (43), we can separate it into two distinct terms as follows:
Let us consider the first term:

$$ \frac{k_1 \exp(-k_1 t) \exp((k_1 - k_s)\tau)}{F_N(t)} \quad (44) $$

By dividing and multiplying by $\frac{k_1}{(k_1 - k_s)} \exp(-k_1 t) (\exp((k_1 - k_s)t) - \exp((k_1 - k_s)t_0))$, we obtain:

$$ \frac{k_1 \exp(-k_1 t) \exp((k_1 - k_s)\tau)}{\frac{k_1}{(k_1 - k_s)} \exp(-k_1 t) (\exp((k_1 - k_s)t) - \exp((k_1 - k_s)t_0)) \times \frac{k_1}{(k_1 - k_s)} \exp(-k_1 t) (\exp((k_1 - k_s)t) - \exp((k_1 - k_s)t_0))}{F_N(t)} $$

Similarly, for the second term:

$$ \frac{k_{-1} \exp(-k_{-1} t) \exp(-k_{-1} \tau)}{F_N(t)} \quad (45) $$
<figure><img src="image_24.png" /><figcaption>Fig. 22. Relative error calculated numerically according to Eq. (22) considering N = 10⁶ samples, as a function of dimensionless time, defined as $t_{dim} = tk_1$, for species A, B and C, in the PE stiff case where. The analysis spans several orders of magnitude separation in time scales between slow and fast reactions, ranging from $k_1/k_s = 10$ to $k_1/k_s = 10^8$.</figcaption></figure>

We perform the same development and obtain:

$$ \frac{k_{-1} \exp(-k_2 t) \exp(-k_{-1} \tau)}{\exp(-k_2 t) (\exp(-k_{-1} t_0) - \exp(-k_{-1} t))} \times \frac{\exp(-k_2 t) (\exp(-k_{-1} t_0) - \exp(-k_{-1} t))}{F_N(t)} $$

we can rewrite the probability density in Eq. (43) as:

$$ \pi_B^{IS}(\tau) = p(t)\pi_B^{KC}(\tau) + q(t)\pi_B^{PE}(\tau) \quad (46) $$

with:

$$ p(t) = \frac{\frac{k_1}{(k_1-k_s)} \exp(-k_1 t) (\exp((k_1-k_s)t) - \exp((k_1-k_s)t_0))}{F_N(t)} $$

and,

$$ q(t) = \frac{\exp(-k_2 t)(\exp(-k_{-1} t_0) - \exp(-k_{-1} t))}{F_N(t)} $$

We can observe that $p(t)$ and $q(t)$ are two complementary probabilities that depend on the observation time, $t$, rather than the random time, $\tau$. This set of probabilities accurately reflects the behavior of mass fraction of A. For instance, as $t$ approaches $t_0$, i.e., when near the initial time, $p(t) \to 1$ and $q(t) \to 0$, indicating that A is solely consumed following $k_1$ at these instances. As $t$ tends towards infinity, i.e., tending towards thermodynamic control, $p(t) \to 0$ and $q(t) \to 1$, which aptly reflects the assumptions formulated for this case arising from the partial equilibrium approximation.

Now by substituting the new probability density into the expectation-based formulation of the mass fraction of B, we find:

$$ Y_B(t) = \int_{t_0}^{t} \pi_B^{IS}(\tau) d\tau \left[ \frac{k_1 \exp(-k_s(\tau - t_0))}{\pi_B^{IS}(\tau)} Y_A(t - t_0 - \tau) + Y_{B,0} \exp(-k_s(t - t_0)) \right] \quad (47) $$
<figure><img src="image_25.png" /><figcaption>Fig. 23. Graphical representation of backward-in-time random walks on the equation for estimating the mass fraction of C at $t_{obs}$ in the reversible PE stiff case. All paths start from species C at $t_{obs}$ (the orange dot in the figure). The presented paths illustrate the behavior of the majority of paths under the QSS stiff case. Specifically, as $k_1 \gg k_s$, they predominantly end at the initial condition of $Y_B$ (red and blue squares). They rarely, if ever, end at the initial mass fraction of $Y_A$.</figcaption></figure>

By substituting $\pi_B^{IS}(\tau)$ with its expression given by Eq. (46), the previous equation becomes:

$$ Y_B(t) = p(t) \int_{t_0}^{t} \pi_B^{KC}(\tau) d\tau + q(t) \int_{t_0}^{t} \pi_B^{PE}(\tau) d\tau + \left[ \frac{k_1 \exp(-k_s(\tau - t_0))}{\pi_B^{IS}(\tau)} Y_A(t - t_0 - \tau) + Y_{B,0} \exp(-k_s(t - t_0)) \right] + \left[ \frac{k_1 \exp(-k_s(\tau - t_0))}{\pi_B^{IS}(\tau)} Y_A(t - t_0 - \tau) + Y_{B,0} \exp(-k_s(t - t_0)) \right] \quad (48) $$

The new expectation-based formulation proposed in Eq. (48) leads to a different algorithmic approach, as demonstrated in Algorithm 8. Here, the realization of the random variable, $W_{Y_B}$, begins with a Bernoulli test based on the observation time. The probability of retaining the weight related to $p(t)$, which determines random time sampling according to the probability density, $\pi_B^{KC}$, increases when closer to kinetic control. Conversely, when far from kinetic control and partial equilibrium between A and B is established, there is a higher chance of retaining the weight corresponding to the probability, $q(t)$. This leads to sampling random times according to the probability density, $\pi_B^{PE}$, constructed using the partial equilibrium approximation. Notably, in this new expectation-based formulation, the initial condition accumulates alongside the realization of the random variable, $W_{Y_B}$, considering its exponential decay as shown in Algorithm 8.

For species C, as the reaction converting B to C only occurs when partial equilibrium is established between A and B, we replaced $Y_B$ with its approximate analytical form given in Eq. (40). By applying the principle of zero variance, we found that the optimal probability density is:

$$ \pi_C^{PE}(\tau) = \frac{k_s \exp(-k_s \tau)}{\exp(-k_s t_0) - \exp(-k_s t)} \quad (49) $$

By applying the new formulation presented in Eq. (48) to estimate the mass fraction of B under partial equilibrium conditions, and incorporating the updated probability density from Eq. (49) for estimating the mass fraction of C, our Monte Carlo method yields this time results in close agreement with analytical solutions, according to Fig. 24. Furthermore, the relative errors associated with Monte Carlo estimations demonstrate satisfactory convergence, with errors in the estimation of

### Algorithm 8 $W_{Y_B}(\tau)$ sampling algorithm

1. **Sample** a random time $r$ uniformly in $[0,1]$.

2. **if** $r < q(t)$:

3. **Sample** a random time $\tau$ following the probability density $\pi_B^{KC}(\tau)$.

4. **Sample** $W_{Y_A}(t-\tau)$ according to Algorithm 4.

5. **return** $\frac{k_1 \exp(-k_s(\tau - t_0))}{\pi_B^{IS}(\tau)} w_{Y_A}(t - t_0 - \tau) + Y_{B,0} \exp(-k_s(t - t_0))$.

6. **else**:

7. **Sample** a random time $\tau$ following the probability density $\pi_B^{PE}(\tau)$.

8. **Sample** $W_{Y_A}(t-\tau)$ according to Algorithm 4.

9. **return** $\frac{k_1 \exp(-k_s(\tau - t_0))}{\pi_B^{IS}(\tau)} w_{Y_A}(t - t_0 - \tau) + Y_{B,0} \exp(-k_s(t - t_0))$.

the mass fractions of A, B, and C remaining within or only slightly exceeding 5% across the relevant transient domain, as displayed in Fig. 25.

### 3.3. Case study: model validation

To establish a proof of concept for the presented methodology and to demonstrate its modeling flexibility in managing stiffness within a complex reaction scheme combining both quasi-steady-state (QSS) and partial equilibrium (PE) stiff scenarios, the reaction scheme presented in Fig. 26 was considered.

All reactions are assumed to be either unimolecular reactions or pseudo-unimolecular reactions. As a consequence, the mathematical model representing the time evolution of the mass fractions consists of a system of coupled first-order ordinary differential equations, formulated as follows:

$$ \left\{ \begin{aligned} \frac{dY_A(t)}{dt} &= -k_1 Y_A(t) + k_{-1} Y_B(t) \\ \frac{dY_B(t)}{dt} &= -(k_2 + k_{-1}) Y_B(t) + k_1 Y_A(t) \\ \frac{dY_C(t)}{dt} &= -(k_3 + k_4) Y_C(t) + k_2 Y_B(t) + k_{-3} Y_D(t) \\ \frac{dY_D(t)}{dt} &= -k_{-3} Y_D(t) + k_3 Y_C(t) \\ \frac{dY_E(t)}{dt} &= k_4 Y_C(t) \end{aligned} \right. \quad (50) $$
<figure><img src="image_26.png" /><figcaption>Fig. 24. Evolution of the mass fractions of the three species A, B, and C as a function of the dimensionless time, $t_{dim} = tk_1$, for the PE stiff case, characterized by an eight-order-of-magnitude separation in time scales between slow and fast reactions (with $k_1 = 10^8 s^{-1}$ and $k_\tau = 1 s^{-1}$), given by Monte Carlo using the newly proposed density probabilities in Eqs. (48) and (49), to estimate respectively the mass fractions of B and C, with the associated error bars (symbols) and the analytical solution (solid lines).</figcaption></figure>

With initial conditions at the starting time, $t = t_0$, are:

$$ 
\begin{aligned}
Y_A(t_0) &= Y_{A,0} = 1 \\
Y_B(t_0) &= Y_{B,0} = 0 \\
Y_C(t_0) &= Y_{C,0} = 0 \\
Y_D(t_0) &= Y_{E,0} = 0 \\
Y_E(t_0) &= Y_{D,0} = 0
\end{aligned}
\qquad (51) $$

Given that $k_{-1}/k_1 = k_2/k_1 = k_{-3}/k_1 = 10^5$ and $k_1/k_3 = k_1/k_4 = 10^3$, the arbitrary choice of these kinetic constants was made to intentionally introduce different types of stiffness:

* The quasi-steady-state (QSS) stiff case for species B and D

* The partial equilibrium (PE) stiff case for species A and C

The purpose of considering several kinetic constants with different orders of magnitude was to demonstrate that our approach to handling stiffness is also valid for such scenarios.

Following the same procedure described in the main text, using importance sampling in combination with the approximations related to partial equilibrium ($dY_B/dt \approx 0$ and $dY_D/dt \approx 0$) we were able to derive an optimal probability density function under the form:

$$ \pi_C^{IS}(\tau) = p(t)\pi_C^{KC}(\tau) + q(t)\pi_C^{PE}(\tau) \qquad (52) $$

By setting $k_s = k_2 + k_{-1}$ and $k_\Sigma = k_3 + k_4$, the terms appearing in Eq. (52) have the following expressions:

$$ \pi_C^{KC}(\tau) = \frac{\left(\frac{k_1 k_2}{k_s} - k_\Sigma\right) \exp\left(\left(\frac{k_1 k_2}{k_s} - k_\Sigma\right) \tau\right)}{\exp\left(\left(\frac{k_1 k_2}{k_s} - k_\Sigma\right) t\right) - \exp\left(\left(\frac{k_1 k_2}{k_s} - k_\Sigma\right) t_0\right)} \qquad (53) $$

and

$$ \pi_C^{PE}(\tau) = \frac{k_3 \exp(-k_3 t) \exp(-k_4 \tau)}{\exp(-k_3 t) (\exp(-k_4 t_0) - \exp(-k_4 t))} \qquad (54) $$

The probabilities $p(t)$ and $q(t)$ are given by:

$$ p(t) = \frac{\left(\frac{k_1 k_s}{k_1 k_2 - k_s k_\Sigma}\right) \exp\left(-\frac{k_1 k_s}{k_s} t\right) \left(\exp\left(\left(\frac{k_1 k_s}{k_s} - k_\Sigma\right) t\right) - \exp\left(\left(\frac{k_1 k_s}{k_s} - k_\Sigma\right) t_0\right)\right)}{F_N(t)} $$

and

$$ q(t) = \frac{\exp(-k_3 t) (\exp(-k_4 t_0) - \exp(-k_4 t))}{F_N(t)} $$

where

$$ F_N(t) = \left( \frac{k_1 k_2}{k_1 k_2 - k_s k_\Sigma} \right) \exp\left(-\frac{k_1 k_2}{k_s} t\right) \times \left( \exp\left(\left(\frac{k_1 k_2}{k_s} - k_\Sigma\right) t\right) - \exp\left(\left(\frac{k_1 k_2}{k_s} - k_\Sigma\right) t_0\right) \right) + \exp(-k_3 t) (\exp(-k_4 t_0) - \exp(-k_4 t)) $$

By applying the new density probability presented in Eq. (52) to estimate the mass fraction of C under partial equilibrium conditions, our Monte Carlo method yields this time results in close agreement with the numerical BDF (backward differentiation formulas) [14] method, according to Fig. 27.

### 3.4. Computation time evaluation

All simulations were executed on a standard laptop (Lenovo 81LK) equipped with an Intel Core i5-8265U processor (quad-core, 1.60–3.90 GHz) and 8 GB of RAM.

As shown in Table 4, the CPU time remains modest for a relative error of 5%, although it increases substantially with both the stiffness and the complexity of the kinetic scheme. The time is multiplied by a factor of three when transitioning from a non-stiff to a stiff reversible scheme, and by an order of magnitude when moving from simple to moderately complex mechanisms (i.e., from 2–3 reactions to 6 reactions). This variation reflects both the increased number of elementary sampling steps and the fact that, in complex schemes, more reactions tend to operate under partial equilibrium. As a result, a larger number of mass fractions are estimated using importance sampling with newly defined probability densities over the interval $[t_0, t]$, which in turn reduces the occurrence of stopping conditions (see Algorithms 7 and 8).

These results indicate that, with the proposed integral formulations, the computational cost is not insensitive to model complexity. On the contrary, as the number of reactions and stiffness levels increase, so does the computational effort required. While memory usage remains negligible, the observed growth in CPU time suggests that naive extrapolation to large-scale reaction mechanisms (e.g., with hundreds of reactions) must be approached cautiously.

At this stage, the main objective was to validate and evaluate the algorithmic strategies on a controlled set of representative schemes, and not to demonstrate full scalability across large systems. A thorough benchmark on large-scale mechanisms (with 50–100 reactions or more) remains an important perspective for future work.
<figure><img src="image_27.png" /><figcaption>**Fig. 25.** Relative error calculated numerically according to Eq. (22) considering *N* = 10⁶ samples, as a function of dimensionless time, defined as *t*<sub>dim</sub> = *tk*₁, for species A, B and C, in the PE stiff case using the newly proposed density probabilities in Eq. (48) and (49), to estimate respectively the mass fractions of B and C. The analysis spans several orders of magnitude separation in time scales between slow and fast reactions, ranging from *k*₁/*k*<sub>*s*</sub> = 10 to *k*₁/*k*<sub>*t*</sub> = 10⁸.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Fig. 26.** Complex reaction scheme involving irreversible and reversible reactions.</figcaption></figure>

Although the reaction schemes considered in this study are relatively simple, the integral Monte Carlo method has already demonstrated strong scalability potential in other application domains [31, 57]. We can expect similar scalability in the field of chemical kinetics.

Indeed, the backward-in-time formulation, combined with a double-

randomization strategy, enables the estimation of the mass fraction of a target species at a specific observation time without requiring computation of the full system state. Moreover, unlike time-interval-based estimations typically used in deterministic approaches, the pointwise nature of Monte Carlo estimators proves particularly advantageous when focusing on selected species or discrete time outputs within large systems.

These two features – double randomization and pointwise estimation – naturally lead to a multi-level parallelism that can be efficiently exploited by modern computing architectures. Specifically, each Monte Carlo path can be sampled independently (parallelism at the path level), observations at multiple time points can be computed in parallel (parallelism at the temporal level), and the mass fractions of all species can be estimated concurrently (parallelism at the species level).

However, in the present work, parallelism was not activated, as the computational times remained reasonable for the considered examples. Nevertheless, this option becomes particularly relevant when scaling
<figure><img src="image_29.png" /><figcaption>**Fig. 27.** Evolution of the mass fractions of the species A, B, C, D and E as a function of the dimensionless time $t_{dim} = tk_4$, for the non-stiff where $k_3, k_4 \ll k_1 \ll k_2, k_{-1}, k_{-3}$, given by Monte Carlo with the associated error bars (symbols) and the numerical solution given by BDF method (solid lines).</figcaption></figure>

<caption>**Table 4**</caption>

Computational performance of the Integral Formulation Monte Carlo method for the different reaction schemes. The reported CPU time corresponds to the total computation time required to estimate all species mass fractions at the different observation times, with a target relative error of 5%.

<table><thead><tr><td>Scheme</td><td>Reaction number</td><td>Species number</td><td>CPU time (s)</td><td>Memory Allocation (%)</td></tr></thead><tbody><tr><td>Irreversible non-stiff</td><td>2</td><td>3</td><td>2</td><td>0.3</td></tr><tr><td>Irreversible stiff</td><td>2</td><td>3</td><td>4</td><td>0.3</td></tr><tr><td>Reversible non-stiff</td><td>3</td><td>3</td><td>2</td><td>0.3</td></tr><tr><td>Reversible stiff</td><td>3</td><td>3</td><td>6</td><td>0.3</td></tr><tr><td>Complex</td><td>6</td><td>5</td><td>23.5</td><td>0.3</td></tr></tbody></table>

up to mechanisms involving a significantly larger number of reactions and species. In addition, the method's inherent parallel structure makes it especially well-suited for deployment on multi-core and distributed architectures, particularly for the simulation of large-scale reaction networks in future work.

## 4. Conclusion

This work introduced a new backward Monte Carlo-based approach for modeling and simulating the chemical kinetics of coupled first-order reactions in relatively simple irreversible and reversible reaction schemes for pedagogical and clarity purposes. This was done under both non-stiff and stiff conditions, with the consideration of timescales related to the quasi-steady-state (QSS) and partial equilibrium (PE) situations.

The performance of the method was evaluated in these scenarios, revealing that stiff cases often involved the occurrence of rare events, especially under partial equilibrium conditions. These rare events led to erroneous values with significant error bars, highlighting the difficulty of convergence for our Monte Carlo method: i) in the irreversible case when estimating the mass fractions of species B and C and ii) in the reversible case for all species.

To address this numerical issue, importance sampling was used to target these rare events, either through the analytical solution of the target function using the zero variance technique or employing approximate forms derived from the partial equilibrium approximation. Optimal probability densities were identified through this approach. Error bars were reduced and results comparable to analytical solutions were achieved without approximating mass fractions for the three involved species, while still employing an unbiased estimator. The approach has been extended to address a more complex reaction pathway, incorporating various types of stiffness, including Quasi-Steady-State

(QSS) and Partial Equilibrium (PE), as well as successive and competitive reactions. By applying the proposed method for handling stiffness, our approach achieved results consistent with those obtained using Backward Differentiation Formulas (BDF).

Future work aims to extend this approach to a non-isothermal regime, incorporating the coupling of a thermal history. This concept has already been demonstrated using the DAEM model, in biomass pyrolysis modeling, for a prescribed thermal history and in radiation physics by coupling radiative transfer to the absorption coefficient model derived from a spectroscopic model, itself expressed as an expectation. The objective is to couple the kinetic model with a more complex thermal history that accurately describes the system's heating via the Arrhenius law and can be reformulated under an expectation-based formulation, using the double randomization technique, while considering the non-linearities of Arrhenius law.

Additionally, the method will be extended to more complex reaction schemes involving a large number of coupled reactions and chemical species. Another key objective is to address non-linear kinetics, where reaction orders deviate from unity. In this regard, recent developments in non-linear Monte Carlo methods – particularly those introduced for solving the Boltzmann equation in gas kinetic theory – offer promising pathways for tackling such challenges. These advances could provide a solid foundation for adapting Monte Carlo strategies to a broader class of chemically reactive systems with non-linear behavior.

### CRediT authorship contribution statement

**Abderrahim Sahim:** Writing – original draft, Visualization, Validation, Methodology, Conceptualization. **Marion Carrier:** Writing – review & editing, Supervision, Resources, Project administration, Methodology, Investigation, Funding acquisition, Conceptualization. **Mouna El Hafi:** Writing – review & editing, Supervision, Project
administration, Methodology, Funding acquisition, Conceptualization. **Simon Eibner:** Writing – review & editing, Validation, Supervision, Methodology, Conceptualization. **Yaniss Nyffenegger-péré:** Writing – review & editing, Validation, Methodology, Investigation, Conceptualization. **Stéphane Blanco:** Writing – review & editing, Visualization, Validation, Methodology, Investigation, Conceptualization. **Richard Fournier:** Writing – review & editing, Validation, Methodology, Investigation, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

The authors acknowledge both French scientific programs: the MOPGA (reference ANR-18-MPGA-0013), the MCMET project (ANR-23-CE6-0002) and the EUR Bio Eco (ANR-18-EURE-0021) managed by the National Research Agency and financially supported by the "Investissements d'Avenir" and also by the Region Occitanie (18016004) for the former.

## Appendix A. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.cej.2025.165009.

## Data availability

Data will be made available on request.

## References

[1] M. Bodenstein, Eine theorie der photochemischen reaktionsgeschwindigkeiten, Z. Phys. Chem. 85 (1) (1913) 329–397, http://dx.doi.org/10.1515/zpch-1913-8512.

[2] L. Michaelis, M.L. Menten, et al., Die kinetik der invertinwirkung, Biochem. Z. 49 (333–369) (1913) 352.

[3] D.L. Chapman, L.K. Underhill, Lv.—The interaction of chlorine and hydrogen. The influence of mass, J. Chem. Soc. Trans. 103 (1913) 496–508, http://dx.doi.org/10.1039/CT9130300496.

[4] J. Bowen, A. Acrivos, A. Oppenheim, Singular perturbation refinement to quasi-steady state approximation in chemical kinetics, Chem. Eng. Sci. 18 (3) (1963) 177–188, http://dx.doi.org/10.1016/0009-2509(63)85003-4.

[5] A. Goeke, S. Walcher, A constructive approach to quasi-steady state reductions, J. Math. Chem. 52 (2014) 2596–2626, http://dx.doi.org/10.1007/s10910-014-0402-5.

[6] D.A. Goussis, Model reduction: When singular perturbation analysis simplifies to partial equilibrium approximation, Combust. Flame 162 (4) (2015) 1009–1018, http://dx.doi.org/10.1016/j.combustflame.2014.09.022.

[7] T. Turanyi, A. Tomlin, M. Pilling, On the error of the quasi-steady-state approximation, J. Phys. Chem. 97 (1) (1993) 163–172, http://dx.doi.org/10.1021/j100103a028.

[8] L. Farrow, D. Edelson, The steady-state approximation: Fact or fiction? Int. J. Chem. Kinet. 6 (6) (1974) 787–800, http://dx.doi.org/10.1002/kin.550060603.

[9] E.H. Flach, S. Schnell, Use and abuse of the quasi-steady-state approximation, IEE Proc. Syst. Biol. 153 (4) (2006) 187–191, http://dx.doi.org/10.1049/ip-syb:20050104.

[10] J. Eilertsen, S. Schnell, S. Walcher, The unreasonable effectiveness of the total quasi-steady state approximation, and its limitations, J. Theoret. Biol. 583 (2024) 111770, http://dx.doi.org/10.1016/j.tbi.2024.111770.

[11] J. Goddard, Consequences of the partial-equilibrium approximation for chemical reaction and transport, Proc. R. Soc. Lond. A Math. Phys. Sci. 431 (1882) (1990) 271–284, http://dx.doi.org/10.1098/rspa.1990.0130.

[12] M. Rein, The partial-equilibrium approximation in reacting flows, Phys. Fluids A Fluid Dyn. 4 (5) (1992) 873–886, http://dx.doi.org/10.1063/1.858267.

[13] C. Curtiss, J. Hirschfelder, Integration of stiff equations, Proc. Natl. Acad. Sci. USA 38 (3) (1952) 235–243, http://dx.doi.org/10.1073/pnas.38.3.235.

[14] C.W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, in: Prentice-Hall Series in Automatic Computation, Prentice-Hall, 1971.

[15] G.D. Byrne, A.C. Hindmarsh, Stiff ODE solvers: A review of current and coming attractions, J. Comput. Phys. 70 (1) (1987) 1–62, http://dx.doi.org/10.1016/0021-9991(87)90001-5.

[16] A. Cuoci, Chapter 13 - numerical modeling of reacting systems with detailed kinetic mechanisms, in: T. Faravelli, F. Manenti, E. Ranzi (Eds.), Mathematical Modeling of Gas-Phase Complex Reaction Systems: Pyrolysis and Combustion, in: Comput. Aided Chem. Eng., vol. 45, Elsevier, 2019, pp. 675–721, http://dx.doi.org/10.1016/B978-0-444-64087-1.00013-9.

[17] D.T. Gillespie, A general method for numerically simulating the stochastic time evolution of coupled chemical reactions, J. Comput. Phys. 22 (4) (1976) 403–434, http://dx.doi.org/10.1016/0021-9991(76)90041-3.

[18] M. Snyder, A. Chatterjee, D.G. Vlachos, Net-event kinetic Monte Carlo for overcoming stiffness in spatially homogeneous and distributed systems, Comput. Chem. Eng. 29 (4) (2005) 701–712, http://dx.doi.org/10.1016/j.compchemeng.2004.09.016.

[19] E.L. Haseltine, J.B. Rawlings, Approximate simulation of coupled fast and slow reactions for stochastic chemical kinetics, J. Chem. Phys. 117 (15) (2002) 6959–6969, http://dx.doi.org/10.1063/1.1505860.

[20] J. He, H. Zhang, J. Chen, Y. Yang, Monte Carlo simulation of kinetics and chain length distributions in living free-radical polymerization, Macromolecules 30 (25) (1997) 8010–8018, http://dx.doi.org/10.1021/ma9614858.

[21] D.T. Gillespie, Approximate accelerated stochastic simulation of chemically reacting systems, J. Chem. Phys. 115 (4) (2001) 1716–1733, http://dx.doi.org/10.1063/1.1378322.

[22] M. Rathinam, L.R. Petzold, Y. Cao, D.T. Gillespie, Stiffness in stochastic chemically reacting systems: The implicit tau-leaping method, J. Chem. Phys. 119 (24) (2003) 12784–12794, http://dx.doi.org/10.1063/1.1627296.

[23] D.G. Vlachos, Stochastic modeling of chemical microreactors with detailed kinetics—induction times and ignitions of H2 in air, Chem. Eng. Sci. 53 (1) (1998) 157–168, http://dx.doi.org/10.1016/S0009-2509(97)00275-3.

[24] Y. Cao, D. Gillespie, L.R. Petzold, Multiscale stochastic simulation algorithm with stochastic partial equilibrium assumption for chemically reacting systems, J. Comput. Phys. 206 (2) (2005) 395–411, http://dx.doi.org/10.1016/j.jcp.2004.12.014.

[25] E.A. Mastny, E.L. Haseltine, J.B. Rawlings, Two classes of quasi-steady-state model reductions for stochastic kinetics, J. Chem. Phys. 127 (9) (2007) http://dx.doi.org/10.1063/1.2764480.

[26] M. Bati, S. Blanco, C. Coustet, V. Eymet, V. Forest, R. Fournier, J. Gautrais, N. Mellado, M. Paulin, B. Piaud, Coupling conduction, convection and radiative transfer in a single path-space: Application to infrared rendering, 42 (4) (2023) http://dx.doi.org/10.1145/3592121.

[27] E. Pahija, S. Hwangbo, T. Saulnier-Bellemare, G.S. Patience, Experimental methods in chemical engineering: Monte Carlo, Can. J. Chem. Eng. 102 (10) (2024) 3308–3321, http://dx.doi.org/10.1002/cjce.25374, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/cjce.25374. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/cjce.25374.

[28] W.B. Sediawan, Application of Monte Carlo Simulation in Chemical Engineering, UGM PRESS, 2024.

[29] L. Zichi, D. Barter, E. Sivonxay, E.W.C. Spotte-Smith, R.S. Mohanakrishnan, E.M. Chan, K.A. Persson, S.M. Blau, RNMC: kinetic Monte Carlo implementations for complex reaction networks, J. Open Source Softw. 9 (104) (2024) 7244, http://dx.doi.org/10.21105/joss.07244.

[30] N. Villefranque, F. Hourdin, L. d'Alençon, S. Blanco, O. Boucher, C. Caliot, C. Coustet, J. Dauchet, M. El Hafi, V. Eymet, et al., The "teapot in a city": A paradigm shift in urban climate modeling, Sci. Adv. 8 (27) (2022) eabp8934, http://dx.doi.org/10.1126/sciadv.abp8934.

[31] L. Ibarrart, S. Blanco, C. Caliot, J. Dauchet, S. Eibner, M. El-Hafi, O. Farges, V. Forest, R. Fournier, J. Gautrais, R. Konduru, L. Penazzi, J.-M. Tregan, T. Vourch, D. Yaacoub, Advection, diffusion and linear transport in a single path-sampling Monte-Carlo algorithm: getting insensitive to geometrical refinement, 2022, working paper or preprint.

[32] Y. Nyffenegger-Péré, R. Armante, M. Bati, S. Blanco, J.-L. Dufresne, M.E. Hafi, V. Eymet, V. Forest, R. Fournier, J. Gautrais, et al., Spectrally refined unbiased Monte Carlo estimate of the Earth's global radiative cooling, Proc. Natl. Acad. Sci. USA 121 (5) (2024) e2315492121, http://dx.doi.org/10.1073/pnas.2315492121.

[33] J.M. Tregan, J.L. Amestoy, M. Bati, J.-J. Bezian, S. Blanco, L. Brunel, C. Caliot, J. Charon, J.-F. Cornet, C. Coustet, L. d'Alençon, J. Dauchet, S. Dutour, S. Eibner, M. El Hafi, V. Eymet, O. Farges, V. Forest, R. Fournier, M. Galtier, V. Gattepaille, J. Gautrais, Z. He, F. Hourdin, L. Ibarrart, J.-L. Joly, P. Lapeyre, P. Lavielle, M.-H. Lecureux, J. Lluc, M. Miscevic, N. Mourtaday, Y. Nyffenegger-Péré, L. Pelissier, L. Penazzi, B. Piaud, C. Rodrigues-Viguier, G. Roques, M. Roger, T. Saez, G. Terrée, N. Villefranque, T. Vourc'h, D. Yaacoub, Coupling radiative, conductive and convective heat-transfers in a single Monte Carlo algorithm: A general theoretical framework for linear situations, PLoS One 18 (4) (2023) 1–54, http://dx.doi.org/10.1371/journal.pone.0283681.
[34] N. Metropolis, S. Ulam, The Monte Carlo method, J. Amer. Statist. Assoc. 44 (247) (1949) 335-341, http://dx.doi.org/10.1080/01621459.1949.10483310, PMID: 18139350.

[35] M. Kac, On some connections between probability theory and differential equations, in: Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability, vol. 1, University of California Press, 1951, pp. 1–10, http://dx.doi.org/10.1525/9780520370184-007.

[36] Y. Nyffenegger-Péré, M. Carrier, S. Eibner, M. El Hafi, M. Nasfi, A. Sahim, A Monte Carlo-based strategy to assess complex kinetics: Application of the null-reaction method to DAEM, J. Phys. Chem. A 127 (8) (2023) 1988–1997, http://dx.doi.org/10.1021/acs.jpca.2c06893, PMID: 36802603.

[37] I. Dimov, T. Dimov, T. Gurov, A new iterative Monte Carlo approach for inverse matrix problem, J. Comput. Appl. Math. 92 (1) (1998) 15–35, http://dx.doi.org/10.1016/S0377-0427(98)00043-0.

[38] I. Dimov, S. Maire, J.M. Sellier, A new walk on equations Monte Carlo method for solving systems of linear algebraic equations, Appl. Math. Model. 39 (15) (2015) 4494–4510, http://dx.doi.org/10.1016/j.apm.2014.12.018.

[39] A.K. Burnham, R.L. Braun, Global kinetic analysis of complex materials, Energy Fuels 13 (1) (1999) 1–22, http://dx.doi.org/10.1021/ef9800765.

[40] W.L. Dunn, J.K. Shultis, Exploring Monte Carlo Methods, Elsevier, 2022.

[41] W.A. Adkins, M.G. Davidson, First order differential equations, in: Ordinary Differential Equations, Springer New York, New York, NY, 2012, pp. 1–100, http://dx.doi.org/10.1007/978-1-4614-3618-8_1.

[42] I. Stakgold, M.J. Holst, Green's functions and boundary value problems, in: Pure and Applied Mathematics: A Wiley Series of Texts, Monographs and Tracts, Wiley, 2011.

[43] J. Mathews, R.L. Walker, Mathematical methods of physics, in: Addison-Wesley World Student Series, W.A. Benjamin, 1970.

[44] E.N. Economou, Green's functions in quantum physics, in: Springer Series in Solid-State Sciences, Springer Berlin Heidelberg, 2006.

[45] G. Terrée, M. El Hafi, S. Blanco, R. Fournier, J. Dauchet, J. Gautrais, Addressing the gas kinetics Boltzmann equation with branching-path statistics, Phys. Rev. E 105 (2022) 025305, http://dx.doi.org/10.1103/PhysRevE.105.025305.

[46] M. Bossy, N. Champagnat, H. Leman, S. Maire, L. Violeau, M. Yvinec, Monte Carlo methods for linear and non-linear Poisson-Boltzmann equation, ESAIM: Proc. 48 (2015) 420–446, http://dx.doi.org/10.1051/proc/201448020.

[47] S. Maire, G. Nguyen, Stochastic finite differences for elliptic diffusion equations in stratified domains, Math. Comput. Simul. 121 (2016) 146–165, http://dx.doi.org/10.1016/j.matcom.2015.09.008.

[48] J.H. Halton, Sequential Monte Carlo techniques for the solution of linear systems, J. Sci. Comput. 9 (2) (1994) 213–257, http://dx.doi.org/10.1007/BF01578388.

[49] L.A. Segel, M. Slemrod, The quasi-steady-state assumption: A case study in perturbation, SIAM Rev. 31 (3) (1989) 446–477, http://dx.doi.org/10.1137/1031091.

[50] J.D. Ramshaw, Partial chemical equilibrium in fluid dynamics, Phys. Fluids (1980) http://dx.doi.org/10.1063/1.863052.

[51] J.A. Bucklew, Introduction to Rare Event Simulation, vol. 5, Springer, 2004.

[52] R.B. Snadden, A new perspective on kinetic and thermodynamic control of reactions, J. Chem. Educ. 62 (8) (1985) 653, http://dx.doi.org/10.1021/ed062p653.

[53] M.-S. Oh, J.O. Berger, Adaptive importance sampling in monte carlo integration, J. Stat. Comput. Simul. 41 (3–4) (1992) 143–168, http://dx.doi.org/10.1080/00949659208810398.

[54] R. Assaraf, M. Caffarel, Zero-variance principle for Monte Carlo algorithms, Phys. Rev. Lett. 83 (1999) 4682–4685, http://dx.doi.org/10.1103/PhysRevLett.83.4682.

[55] C.W. Pyun, Steady-state and equilibrium approximations in chemical kinetics, J. Chem. Educ. 48 (3) (1971) 194, http://dx.doi.org/10.1021/ed048p194.

[56] V. Viossat, R.I. Ben-Aim, Validity of the quasi-stationary-state approximation in the case of two successive reversible first-order reactions, J. Chem. Educ. 75 (9) (1998) 1165, http://dx.doi.org/10.1021/ed075p1165.

[57] N. Villefranque, R. Fournier, F. Couvreux, S. Blanco, C. Cornet, V. Eymet, V. Forest, J.-M. Tregan, A path-tracing Monte Carlo library for 3-D radiative transfer in highly resolved cloudy atmospheres, J. Adv. Model. Earth Syst. 11 (8) (2019) 2449–2473, http://dx.doi.org/10.1029/2018MS001602.