Revisiting Kinetic Monte Carlo Algorithms for Time-dependent Processes: from open-loop control to feedback control

This article begins by revisiting the theoretical foundation of the open-loop-control KMC method by providing a pedagogical proof that recapitulates various versions of these algorithms.Anderson (2007); Anderson, Ermentrout, and Thomas (2015); Anderson and Kurtz (2011); Chow and White (1996); Li et al. (2009) The theory is constructed by comparing the true probability of single-step trajectories and the simulated trajectory probability generated within one step of a modified KMC algorithm. This trajectory analysis not only provides an alternative proof for open-loop-control KMC methods described previouslyAnderson (2007); Anderson, Ermentrout, and Thomas (2015); Anderson and Kurtz (2011); Chow and White (1996); Li et al. (2009) but also provides the foundation from which a closed-loop-control KMC method can be derived.

Consider a continuous-time discrete-state Markov process that is subject to an external control protocol $\lambda(t)$ as an arbitrary function of time $t$. The system contains a discrete number of states $n_{s}$ and is initialized with probability distribution $\vec{p}(t=0)$. The dynamics of $\vec{p}(t)$ is characterized by the master equation

where matrix element $R_{ji}(\lambda(t))$ is the transient transition probability rates from $i$ to $j$ given the control parameter $\lambda$ at time $t$. The diagonal element satisfying $R_{ii}(t)=-\sum_{j\neq i}R_{ji}(t)$ is the negative transient escape rate from state $i$ at time $t$. The master equation evolution of $\vec{p}(t)$ can also be expressed in terms of the dynamics of $i$-th state’s probability:

For an open-loop controlled system, the control protocol $\lambda(t)$ is externally determined. In contrast, for a closed-loop controlled system (e.g., one with a series of feedback measurements), the specific form of $\lambda(t)$ may depend on the system’s state.

Each individual step of a KMC simulation is concerned with only generating the next stochastic event (i.e., evolving to state state $x_{k+1}$ at time $t_{k+1}$ from the previous state $x_{k}$ and time $t_{k}$). This event can be expressed as a single-step trajectory $X_{\rm single}=(x_{k},t_{k},x^{\prime},t^{\prime})$ given the previous step. The conditional probability for this single-step trajectory can be written as

which is also a joint probability distribution of the next event’s time $t^{\prime}$ and state $x^{\prime}$. This formula is true for any Markov dynamics regardless of whether the dynamical rates are constants or the rates change over time. We show below that, in the case of time-independent or time-dependent rates, the direct Gillespie Algorithm or an open-loop-control KMC algorithm, respectively, can correctly capture the dynamics.

Time-independent rates: In the case where the dynamics of the system are constant over time ¹¹1Here constant rate over time is used to describe the time period between two consecutive stochastic events. Naturally, when the system state changes after the next event, the possible transitions and their rates are naturally updated., the direct Gillespie AlgorithmGillespie (1977) is sufficient. From a probability theory perspective, the single-step transition probability of the next state $x^{\prime}$ and next time $t^{\prime}$ can be factored into the product of two statistically independent distributions and can be generated separately from their own probability distributions $P_{x}({x^{\prime}})$ and $P_{t}({t^{\prime}})$.

In this case, the direct KMC Gillespie (1977) is carried out by generating the next transition time $t^{\prime}$ according to $P_{t}({t^{\prime}})\propto e^{{R_{x_{k}x_{k}}\cdot({t^{\prime}}-t_{k})}}$ and generating the next state $x^{\prime}$ according to the state probability $P_{x}({x^{\prime}})\propto R_{{x^{\prime}}x_{k}}$.

Time-dependent rates: If the rates change over time, the direct KMC approachGillespie (1977) fails to generate the correct dynamics since the next event’s time and state are statistically correlated. In this case, we can still factorize the joint probability of $t^{\prime}$ and $x^{\prime}$, $P({x^{\prime}},{t^{\prime}})\equiv p[{x^{\prime}},{t^{\prime}}|x_{k},t_{k}]$ as the product of two probability distributions:

where the marginal distribution of the transition time $t^{\prime}$ (given the previous $x_{k}$ and $t_{k}$) is denoted by

and the conditional probability distribution of next state $x^{\prime}$ conditioned on the transition time being $t^{\prime}$ is

In this case, a modified KMC algorithm should be capable of faithfully generating the stochastic trajectory. In one simulation step of the modified KMC algorithm, one first generates the next transition’s time $t^{\prime}=t^{*}$ according to Eq. 6 and then generates the next state $x^{\prime}$ according to Eq. 7 given the proposed time $t^{\prime}=t^{*}$. One can verify that the generated pair of $x^{\prime}$ and $t^{\prime}$ follows the true single-step trajectory probability distribution (Eq. 3). This analysis is applicable to arbitrary time-dependent processes, where the system’s rates change over time either due to (1) open-loop control according to a scheduled protocolAnderson (2007); Anderson, Ermentrout, and Thomas (2015); Anderson and Kurtz (2011); Chow and White (1996); Li et al. (2009) or (2) closed-loop control with a series of feedback measurements.

In this section, we review the open-loop-control KMC methods previously described in Anderson (2007); Anderson, Ermentrout, and Thomas (2015); Anderson and Kurtz (2011); Chow and White (1996); Li et al. (2009). This review also serves as an alternative derivation and justification of these algorithms in the language of single-step trajectory generation introduced in Section II.1. The representation of the single-step trajectory probability as the product of a marginal probability and a conditional probability, Eq. 5, allows one to separately generate the next transition’s time $t^{\prime}$ and state $x^{\prime}$, even when the time-dependence of the system’s rate creates statistical correlation between $t^{\prime}$ and $x^{\prime}$.

Generating transition time $t^{\prime}$: The random generation of transition time $t^{\prime}>t_{k}$ according to the marginal probability density in Eq. 6 can be realized by first obtaining the cumulative density function

and then generating a random number $u$ following a uniform distribution between $0$ and $1$, and finally solving for $t^{*}$ according to the following equation:

Or equivalently, for a random number $v=1-u$ also following a uniform distribution between $0$ and $1$, one can solve $t^{*}$ as the root of the following equation

where

and $R_{xx}(t)$ is the negative escape rate from state $x$ at arbitrary time $t$.

Generating the next state $x^{\prime}$: Once the transition time $t^{*}$ is generated, the selection of the next state $x^{\prime}$ is straightforward – the new state’s probability is proportional to their transient transition rates evaluated at the chosen transition time $t^{*}$, as defined in Eq. 7. This resembles the state generation in the direct KMC method, with the only difference being the probability to jump to state $x^{\prime}$ here is proportional to the transient rate $R_{x^{\prime}x_{k}}(t^{*})$ evaluated at the chosen time $t^{*}$.

To summarize, one simulation step of the open-loop-control KMC algorithm can be enumerated in Section II.2. Furthermore, by including appropriate approximations (see Section II.3), one can reproduce various open-loop-control KMC algorithms previously described in Anderson (2007); Anderson, Ermentrout, and Thomas (2015); Anderson and Kurtz (2011); Chow and White (1996); Li et al. (2009).

Open-loop-control KMC

In practice, solving for the next event’s time $t^{*}$ by a given random number $v$ and the equation Eq. 10 for arbitrarily complex time-dependent rates may be challenging. We now briefly review a few ways to approximately solve for this integral equation. Various methods have been proposed in simulating specific biological processes.Anderson, Ermentrout, and Thomas (2015); Chow and White (1996); Ding et al. (2016); Li et al. (2009) For illustrative purposes, we demonstrate two approximations and sketch the corresponding algorithms for: (i) piece-wise linear function approximation and (ii) piece-wise step function approximation for the escape rate $-R_{xx}(t)$.

Here we illustrate the algorithms by focusing on one step of the simulation. Let us denote the previous state as $x_{k}=x$ and the previous time $t_{k}=t_{\rm last}$; then the algorithm’s task is to generate the next event’s time $t^{*}$ and the next state $x^{\prime}$. As illustrated in Fig. 1a, the negative escape rate from the state $x$, $R_{xx}(t)$ (shown as black curve) can be approximated by (i) a consecutive set of straight line segments (top panel), or (ii) approximated by a set of step-wise step functions (bottom panel).

(i) Piece-wise linear escape rates. Let us denote each window $w_{i}$ by a time interval $(\tau_{i},\tau_{i+1})$ and also set the beginning of the first window at $\tau_{1}=t_{\rm last}$. The negative escape rate from state $x$ at each window’s left edge is denoted by $R_{xx}^{\tau_{i}}=-\sum_{x^{\prime}\neq x}R_{x^{\prime}x}(\tau_{i})$. Under this assumption, the numerical integration involved in Eqs. 10 and 11 can be considered as:

for $\tau_{i}<t^{*}<\tau_{i+1}$. Here the integral for each window starting from $\tau_{1}=t_{\rm last}$ is defined as

Here, $A_{i}$ is the integrated negative escape rate within the window $w_{i}$ defined from $\tau_{i}$ to $\tau_{i+1}$. Successive windows $w_{i}$ starting from $i=1$ will be evaluated until the $j$-th window where the transition time is chosen, which can be identified by $\sum_{i=1}^{j-1}A_{i}<-\ln v<\sum_{i=1}^{j}A_{i}$. Given Eqs. 12 and 13, we can solve for the $t^{*}$ according to Eq. 10 in sequence, starting from $i=1$ and continuing until the appropriate $t^{*}$ can be drawn from the $j-$th window according to:²²2In practical implementation of the numerical code, when $q_{i}\approx q_{i-1}$, the quadratic equation reduces to the linear, which can be solved from Eq. 15.

(ii) Piece-wise step-function escape rates. When the time dependent escape rate from state $x$ is considered as a piece-wise step function, where $R_{xx}(t)=R_{xx}^{\tau_{i}}$ for each window $t\in(\tau_{i},\tau_{i+1})$, we can determine the next event’s time $t^{*}$ by solving for a linear equation

where $v$ is a random number uniformly distributed between $0$ and $1$ and the integrated negative escape rate for each window is denoted by

By solving for the above equations starting from the first window, one can find the next event’s time as $t^{*}\in(\tau_{k},\tau_{k+1})$ occurring at $k$-th window. We demonstrate this method in a 1D lattice diffusion problem as shown in the Appendix.

Consider a feedback control scheme where a scheduled series of measurements will be taken at a set of observation time points $\tau_{1},\cdots,\tau_{i},\tau_{i+1},\cdots\tau_{i+k},\cdots$ (see Fig. 2a). At each observation time $\tau_{j}$, the measurement outcome is determined by the transient state of the stochastic system: $m_{j}=x(\tau_{j})$.³³3Here the measurement of the system may not be as detailed as the specific state that the system is in. It is possible that a subset of states will give a same measurement outcome (i.e., coarse-grained measurement). The method can be applied to both cases. Immediately, the measurement outcome updates the system’s kinetic rates via a rapid feedback control mechanism. Let us denote the control signal as an explicit function of time $\lambda(t)$, which can be represented by a piece-wise function where each piece is defined for the time between two consecutive measurements and is determined by all previous measurement outcomes. The actual protocol can be piece-wise constant (see illustrative example in Sections II.5 and III.1), or in more general cases, could be a solution of an ODE parameterized by the measurement outcomes at each feedback measurement. The general algorithm proposed here is applicable to arbitrary functional forms of the control protocol.

The key to correctly simulate the feedback controlled system is to recognize that, when carrying out each step of a KMC simulation, one only needs to predict the next transition event given the present state and the transition rates set by the control protocol. Now consider a system placed at state $x$ and time $t_{\rm last}$ set by the previous step of simulation. Without losing generality, let us consider the current time as $\tau_{i-1}<t_{\rm last}<\tau_{i}$, so that the current system has gone through the $(i-1)-$th measurement and the very next feedback measurement (the $i-$th measurement) will occur at a future time $\tau_{i}$.

To carry out the present simulation step and predict the next transition time $t^{*}$ and next state $x^{\prime}$, there are only two possible scenarios to consider, which are separated by a milestone time – the immediate next measurement time, $\tau_{i}$. The two scenarios are the following: (i) the next event takes place at a time before the milestone time ($t^{*}<\tau_{i}$, Fig. 2b[i]) or (ii) the next event takes place at any time after the milestone time ($t^{*}>\tau_{i}$, Fig. 2b[ii]). The milestone time $\tau_{i}$ is highlighted by the color scheme transition from red to yellow in Fig. 2b[ii]. This color scheme is also used in Eq. 17 and Eq. 18.

In the first case, one does not need to consider the control protocol after the next measurement and only needs to consider the protocol before $\tau_{i}$. For this case, the control protocol for any $t<\tau_{i}$ is fully determined without any ambiguity, given the historical trajectory of the simulation. This pre-milestone part of control protocol before $\tau_{i}$ can be generally represented by the following function and is fully determined:

where $\lambda_{j}$ denotes the control protocol for time window $\tau_{j}<t\leq\tau_{j+1}$, and its functional form can be impacted by the all previous measurement results $m_{1},m_{2},\cdots,m_{j}$ through any prescribed feedback mechanisms. Here, $j$ must be smaller than $i$.

In the second case, we need to treat the system with a step-wise protocol where the system still evolves according to the pre-milestone protocol (Eq. 17) up to the milestone time $\tau_{i}$, and then the protocol is updated to a future post-milestone protocol. It may appear that the post-milestone protocols (for $t>\tau_{i}$) cannot be fully determined due to the ambiguity of the future system states. However, as we have argued, each step of the KMC simulation involves the generation of a single step trajectory that starts at $(x,t_{\rm last})$, remains at state $x$ and only jumps to an unknown state $x^{\prime}$ at a future unknown transition time $t^{*}$. Thus, for any unknown transition time $t^{*}>\tau_{i}$ larger than the milestone time, each futuristic feedback measurement according to the trajectory under consideration will return a measurement outcome based on the system state $x$. In other words, the post-milestone protocol for any future $t\geq\tau_{i}$ is denoted by:

which is determined without ambiguity because all future measurements must detect the system at the current state $x$ (i.e., $m_{i},m_{i+1},m_{i+2},\cdots=x$). In this argument, all futuristic protocols are then determined by the combination of the historical measurement outcomes and the asserted futuristic measurement outcomes, $m_{1},m_{2},\cdots,m_{i-2},{\color[rgb]{0.65,0,0}m_{i-1}},{\color[rgb]{0.85,0.34,0.17}x,x,x,\cdots}$.

Closed-loop-control KMC

In summary, through our single-step trajectory probability argument introduced in Section II.1, we are able to assert a fully determined protocol for a feedback control system. This determined protocol combines Eq. 17, the pre-milestone protocol $\lambda(t)$ for $t<\tau_{i}$, which is fully determined by historical feedback measurements, and Eq. 18, the futuristic post-milestone protocol $\tilde{\lambda}(t)$ for $t>\tau_{i}$, which appears to be unknown due to the future feedback measurement outcomes. As a result, even for closed-loop controlled systems, one can utilize an open-loop-control KMC algorithmAnderson (2007); Anderson, Ermentrout, and Thomas (2015); Anderson and Kurtz (2011); Chow and White (1996); Li et al. (2009) reviewed in Section II.2 to generate the next transition event according to the correct probability distribution Eqs. 5, 6 and 7. This general algorithm is shown in Section II.4. For illustrative purpose, the following demonstrates a simple application of the proposed feedback-control KMC algorithm.

Here, we demonstrate the proposed new algorithm with a general yet simple scenario. In this example, the rates within each measurement window remain constant and are instantaneously updated based on each measurement. For more complicated feedback control regimes (e.g., the control protocol evolves deterministically in time according to differential equations parameterized by historical measurement outcomes), one can simply extend this algorithm in accordance with general feedback controlled KMC algorithm described in Section II.4.

We here describe a model algorithm where the control parameter $\lambda(t)$ between two consecutive measurement events are always maintained at a constant value. In this case, both the determined protocol $\lambda_{i-1}(t)=\lambda_{a},~{}t_{\rm last}<t<\tau_{i}$ and the projected protocols $\tilde{\lambda}(t)=\lambda_{b},~{}t>\tau_{i}$ are both piece-wise constants. Therefore, the escape rate from state $x$ is a piece-wise constant function with only two pieces (see Fig. 2b):

In this case, the simulation chooses the next event’s occurrence time by solving for the solution of the piece-wise function:

where $v$ is a uniformly distributed random number between $0$ and $1$, and the $r.h.s$ of the equation is defined by Eq. 11. There are two types of possible simulation outcomes of the single-step trajectories. At possibility $\#1$ (top panel of Fig. 2b), when the next event occurs before the next measurement, $t_{\rm{last}}<t^{*}<\tau_{i}$, the dynamical rates remains unchanged throughout the time, and the rates was determined by the last measurement outcome $m_{i-1}=x(\tau_{i-1})$ (see red shaded regions in Fig. 2b). At possibility $\#2$ (bottom panel of Fig. 2b), when the next event occurs after $k$ more future measurements $\tau_{i+k-1}<t^{*}<\tau_{i+k}$, for $k=1,2,\cdots$, the dynamics are initially governed by the old rates $\hat{R}(\lambda_{a})$ before the $i$-th measurement, and then by a new rate $\hat{R}(\lambda_{b})$ after $\tau_{i}$ (see yellow shaded regions in Fig. 2b). We demonstrate this method in a two-state refrigerator model in Section III.1.

To demonstrate the new feedback controlled KMC algorithm proposed in Section II.4, we utilize an example feedback control system whose behavior can be solved analytically. Feedback control systems have been extensively studies in various stochastic systemsBechhoefer (2021); Annby-Andersson et al. (2022); Horowitz and Jacobs (2014); Horowitz and Parrondo (2011a, b); Horowitz and Vaikuntanathan (2010); Bhattacharyya and Jarzynski (2022); Sagawa and Ueda (2012a, 2010, b); Parrondo, Horowitz, and Sagawa (2015). Here we choose to simulate the stochastic trajectories of a 2-state refrigerator inspired by referenceMandal, Quan, and Jarzynski (2013).

Consider a 2-level system coupled to a heat bath. The ground state level’s energy is $E_{0}=0$ and the excited level’s energy is controlled externally by a binary toggle switch. At control $A$ the excited state energy is $E_{1}^{A}=1$; at control condition $B$, $E_{1}^{B}=\epsilon>1$. The feedback control is facilitated by a periodic sequence of instantaneous measurements to the system’s state at times $\tau_{i}=i/\nu$, where $\nu$ is the measurement frequency. If the measurement detects the system at the ground state “0”, the control is switched to $B$, and if excited state “1”, control is $A$(see Fig. 3a).

This system falls into the category of the step-wise constant feedback control scenario (discussed in II.5). In this case, the system evolves either according to the transition rates set by control condition $A$ or by condition $B$. At each instantaneous time, the transition rates takes the Arrhenius form $R_{ij}=e^{B_{ij}-E_{j}}$, where we have taken the unit such that the inverse temperature $\beta=1/(k_{B}T)=1$. In our demonstration, we set the barrier for the transition to $B_{ij}=2$, and set the energies according to $E_{1}^{A}=1$, and $E_{1}^{B}=\epsilon=1.5$. Under this feedback control, it was shown that the system behaves as a refrigerator, constantly drawing energy from the heat bath and performs work against the external control. This refrigerator does not require energy expenditure but rather operates at the cost of information obtained via the measurements. For a comprehensive overview of the definition of heat and work in the stochastic systems, refer to the book Sekimoto . For the connection between information and thermodynamic entropy, a few representative works can be found in Penrose (2005); Bennett (1973); Landauer (1961); Mandal and Jarzynski (2012); Zurek (1989); Deffner (2013); Strasberg et al. (2013); Sagawa and Ueda (2012a, 2010, b); Parrondo, Horowitz, and Sagawa (2015); Lu, Mandal, and Jarzynski (2014).

At each stochastic transition from the ground state to the excited state, the system draws from the heat bath energy $Q_{\rm sys}=1$ (for control $A$) or $Q_{\rm sys}=\epsilon$ (for control $B$). At the transition from excited state to the ground state, the system’s heat gain is $Q_{\rm sys}=-1$ (for control $A$) and $Q_{\rm sys}=-\epsilon$ (for control $B$). At any schedule measurement causing a control switch from $A$ to $B$ (the system is at the ground state), there is no work or heat exchange ($W=0$, $Q=0$). At the scheduled measurement causing a control switch from $B$ to $A$, the system performs positive work to the controller $W_{\rm ext}=\epsilon-1$, and there is no heat exchange $Q=0$ within the instantaneous measurement time. By counting the frequencies of the stochastic transitions of system ($0$ to $1$ or $1$ to $0$) under different control conditions $A$ and $B$ respectively, one can calculate the average heat rate (the average heat drawn by the system per unit time).

By implementing the feedback controlled KMC algorithm, we are able to obtain an ensemble of 100 statistically independent stochastic trajectories from such a feedback controlled refrigerator. Each trajectory is of time length 500 and we exclude the initial relaxation period of length 100 from statistics. From the trajectory, one can compute the bath’s entropy change rate, which is equal to the total entropy change rate of the system and the bath:

where the net transition current $J$ is defined by the difference between the detailed currents:

Furthermore, since each measurement is binary, one can calculate the information rate simply by multiplying the information per measurement with the measurement frequency: $k_{B}\dot{I}=\nu\ln 2$. Finally, we also computed the efficiency of the refrigerator as the entropy decrease rate over the information gain rate:

which characterizes the system’s ability to decrease entropy at the expenditure of information change. Notice that according to the generalized second law of thermodynamics $\dot{S}_{tot}+k_{B}\dot{I}\geq 0$,Sagawa and Ueda (2012a, b, 2010); Parrondo, Horowitz, and Sagawa (2015), the efficiency must be lower than or equal to 1. The simulation results for a range of measurement frequencies are shown in Fig. 3b-c.

One can verify the proposed feedback controlled KMC by comparing its result to the analytical solution of the system. In the limit of the infinitely high frequency measurements, $\nu\gg 1$, one can obtain an effective constant rate matrix of the system and analytically obtain its average heat rate and entropy change rate. Due to infinitely frequent measurements, the system’s transition rate from state $1$ to $0$ is always equal to $R^{\rm eff}_{01}=R^{A}_{01}=e^{-1}$ and the system’s transition rate from state $0$ to $1$ is always equal to $R^{\rm eff}_{10}=R^{B}_{10}=e^{-2}$. As a result, the effective rate matrix under infinite frequency feedback control is

The steady state probability for this system can be written as $(p^{\rm ss}_{0},p^{\rm ss}_{1})=(\frac{e}{1+e},\frac{1}{1+e})$. This allows us to compute the detailed transition current from 0 to 1 as $\frac{e}{1+e}e^{-2}$, and current from 1 to 0 as $\frac{1}{1+e}e^{-1}$. Due to the frequent measurement, each transition from 0 to 1 must occur under control $B$, which takes energy $\epsilon$ from the bath, and each transition from 1 to 0 must occurs under control $A$, which deposits energy $1$ into the bath. As a result, for a system at the steady state (no system entropy change), the total entropy change of the system and the bath is purely the bath entropy change with the rate

where the detailed steady-state currents are

From Fig. 3 one can verify that the entropy change rate at high frequency limit obtained from our proposed feedback controlled KMC algorithm matches well with the analytical result in Eq. 26

This work proposed a new algorithm to simulate systems under feedback control. Here we illustrate that for systems with time-dependent control (both open-loop and closed-loop), the direct implementation of the original Gillespie algorithmGillespie (2001) can lead to erroneous results, where the system’s dynamics appears to be immobilized due to frequent changes of the external control parameter or frequent feedback controls. We name this erroneous results the artificial Zeno effect. This effect appears to resemble the quantum Zeno effectFischer, Gutiérrez-Medina, and Raizen (2001) but is fundamentally different – the former is an artificial numerical error, and the latter is a true physical phenomenon.

In the quantum Zeno effectFischer, Gutiérrez-Medina, and Raizen (2001) frequent projective measurements may reset the system’s wave function into an eigenstate and, in the limit of infinitely high measurement frequency, the system’s wave function is frozen at the eigenstate and the system appears to be non-evolving. In contrast, classical stochastic systems, where the measurement does not directly impact the system’s state⁴⁴4Under feedback control measurements, the system’s kinetic rates can be updated after the measurement, but the system’s state is not directly altered by the action of measurement., should not have any Zeno effect. The artificial Zeno effect is a numerical error only when one naively implements the direct KMC algorithm in a piece-wise manner for time-dependent systems.

For illustrative purpose we first demonstrate the erroneous artificial Zeno effect in the feedback control system discussed in Fig. 3. In the Appendix, this discussion has been extended to an open-loop controlled time-dependent systems (Fig. 4).

For the feedback controlled system, one naive way to implement the original Gillespie algorithmGillespie (2001) is to reset the KMC simulation at each measurement time point. Specifically, one may perform a constant-rate KMC simulation for each measurement time window. Consider a system at state $x$ at time $<\tau_{i-1}<t<\tau_{i}$. If the direct KMC approach proposes the next transition’s time to be $t^{*}>\tau_{i}$, one simply evolves the system time to $t^{\prime}=\tau_{i}$ while maintaining the system state to be $x$ and then starts another step of the direct KMC simulation, with the initial state $x$ and initial time $\tau_{i}$. This approach appears to be valid, since the predicted time $t^{*}$ is indeed longer than the time of the next measurement, and the system should not change until after the next measurement. However, this direct implementation of the KMC method leads to wrong results that can become obvious by considering the infinitely frequent measurement limit, where the time till the next measurement approaches $0$ and the chance that the direct KMC proposes a $t^{*}$ before the next measurement becomes zero. Thus the simulation is frozen to state $x$. This artificial Zeno effect becomes more prominent for higher frequency measurements.

To illustrate this mistake, we implement the above direct KMC simulation for the Maxwell’s refrigerator discussed in Section III.1 and show that the artificial Zeno effect significantly underestimates the transition events of the systems dynamics, resulting in an underestimation of the refrigerator’s entropy decreasing rate and thermodynamic efficiency (see Fig. 3bc).

For a feedback controlled system where one measures and acquires system’s states at a scheduled sequence of times, plainly restarting the clock of the KMC simulation from the measured state to simulate the dynamics of the system will lead to an artificial Zeno effect. The closed-loop-control algorithm demonstrated here avoids this pitfall and preserves the correct dynamics. In the Appendix, we show that similar mistakes may occur in open-loop-control systems only if one overlooks previously introduced algorithmsAnderson (2007); Anderson, Ermentrout, and Thomas (2015); Anderson and Kurtz (2011); Chow and White (1996); Li et al. (2009)