Destructive Impact of Molecular Noise on Nanoscale Electrochemical Oscillators

TL;DR

This study reveals that molecular noise significantly disrupts electrochemical oscillations at nanoscale electrodes, affecting their stability, parameter regions, and the relationship between correlation time and period variance.

Contribution

It demonstrates how stochastic electron-transfer events cause potential fluctuations, fundamentally altering the behavior of nanoscale electrochemical oscillators compared to chemical ones.

Findings

Oscillations shift in phase space with decreasing system size.

Minimum molecules for sustained oscillations are over 10 times higher than in chemical systems.

Correlation time and period variance relationship is only valid in limited conditions.

Abstract

We study the loss of coherence of electrochemical oscillations on meso- and nanosized electrodes with numeric simulations of the electrochemical master equation for a prototypical electrochemical oscillator, the hydrogen peroxide reduction on Pt electrodes in the presence of halides. On nanoelectrodes, the electrode potential changes whenever a stochastic electron-transfer event takes place. Electrochemical reaction rate coefficients depend exponentially on the electrode potential and become thus fluctuating quantities as well. Therefore, also the transition rates between system states become time-dependent which constitutes a fundamental difference to purely chemical nanoscale oscillators. Three implications are demonstrated: (a) oscillations and steady states shift in phase space with decreasing system size, thereby also decreasing considerably the oscillating parameter regions; (b) the minimal number of molecules necessary to support correlated oscillations is more than 10 times as large as for nanoscale chemical oscillators; (c) the relation between correlation time and variance of the period of the oscillations predicted for chemical oscillators in the weak noise limit is only fulfilled in a very restricted parameter range for the electrochemical nano-oscillator.

Tables

Table 1. Table 1 : Preexponential factors of the transition rates.

reaction	$W_{\rho }^0$
(9)	$\frac{k_1^0}{\Omega }{C}_{\text{ce}HO}^b\left(\Omega -N_{\text{ce}H}-N_{\text{ce}HO}-N_{\text{ce}X}\right)$
	$\left(\Omega -N_{\text{ce}H}-N_{\text{ce}HO}-N_{\text{ce}X}-1\right)$
(10)	$k_2^0{C}_{\text{ce}H}^b{N}_{\text{ce}HO}$
(11)	$\frac{k_3^0}{\Omega }{N}_{\text{ce}HO}\left(N_{\text{ce}HO}-1\right)$
(12)$\to$	$k_4^0{C}_{\text{ce}H}^b\left(\Omega -N_{\text{ce}H}-N_{\text{ce}HO}-N_{\text{ce}X}\right)$
(12)$\leftarrow$	$k_5^0{N}_{\text{ce}H}$
(13)$\to$	$k_6^0{C}_{\text{ce}H}^b\left(\Omega -N_{\text{ce}{H}^{ot}}\right)$
(13)$\leftarrow$	$k_7^0{N}_{\text{ce}{H}^{ot}}$
(14)	$\frac{k_8^0}{\Omega }{N}_{\text{ce}{H}^{ot}}\left(N_{\text{ce}{H}^{ot}}-1\right)$
(15)$\to$	$k_9^0{C}_{\text{ce}X}^s\left(\Omega -N_{\text{ce}H}-N_{\text{ce}HO}-N_{\text{ce}X}\right)$
(15)$\leftarrow$	$k_{10}^0{N}_{\text{ce}X}$

Table 2. (a) Time independent values of the reaction rate constants k ~ ρ 0 superscript subscript ~ 𝑘 𝜌 0 \widetilde{k}_{\rho}^{0} and of the exponential coefficients c ρ subscript 𝑐 𝜌 c_{\rho} .

${\tilde{k}}_{\rho }^0$	Value		$c_{\rho }$	Value $\left[V^{-1}\right]$
${\tilde{k}}_1^0$	$4.0\cdot {10}^{-2}$	$\mathrm{cm}\cdot s^{-1}$	$c_1$	$0$
${\tilde{k}}_2^0$	$1.0\cdot {10}^{-5}$	$\mathrm{cm}\cdot s^{-1}$	$c_2$	$-19.33$
${\tilde{k}}_3^0$	$1.0\cdot {10}^{-8}$	$\mathrm{mol}\cdot s^{-1}\cdot {\mathrm{cm}}^{-2}$	$c_3$	$19.33$
${\tilde{k}}_4^0$	$1.0\cdot {10}^{-2}$	$\mathrm{cm}\cdot s^{-1}$	$c_4$	$-11.60$
${\tilde{k}}_5^0$	$1.0\cdot {10}^{-5}$	$\mathrm{mol}\cdot s^{-1}\cdot {\mathrm{cm}}^{-2}$	$c_5$	$11.60$
${\tilde{k}}_6^0$	$5.0\cdot {10}^{-3}$	$\mathrm{cm}\cdot s^{-1}$	$c_6$	$-19.33$
${\tilde{k}}_7^0$	$5.0\cdot {10}^{-6}$	$\mathrm{mol}\cdot s^{-1}\cdot {\mathrm{cm}}^{-2}$	$c_7$	$-19.33$
${\tilde{k}}_8^0$	$5.0\cdot {10}^{-6}$	$\mathrm{mol}\cdot s^{-1}\cdot {\mathrm{cm}}^{-2}$	$c_8$	$0$
${\tilde{k}}_9^0$	$5.0\cdot {10}^{-5}$	$\mathrm{cm}\cdot s^{-1}$	$c_9$	$19.33$
${\tilde{k}}_{10}^0$	$5.0\cdot {10}^{-8}$	$\mathrm{mol}\cdot s^{-1}\cdot {\mathrm{cm}}^{-2}$	$c_{10}$	$-19.33$

Table 3. (a) Time independent values of the reaction rate constants k ~ ρ 0 superscript subscript ~ 𝑘 𝜌 0 \widetilde{k}_{\rho}^{0} and of the exponential coefficients c ρ subscript 𝑐 𝜌 c_{\rho} .

${\tilde{k}}_{\rho }^0$	Value		$c_{\rho }$	Value $\left[V^{-1}\right]$
${\tilde{k}}_1^0$	$4.0\cdot {10}^{-2}$	$\mathrm{cm}\cdot s^{-1}$	$c_1$	$0$
${\tilde{k}}_2^0$	$1.0\cdot {10}^{-5}$	$\mathrm{cm}\cdot s^{-1}$	$c_2$	$-19.33$
${\tilde{k}}_3^0$	$1.0\cdot {10}^{-8}$	$\mathrm{mol}\cdot s^{-1}\cdot {\mathrm{cm}}^{-2}$	$c_3$	$19.33$
${\tilde{k}}_4^0$	$1.0\cdot {10}^{-2}$	$\mathrm{cm}\cdot s^{-1}$	$c_4$	$-11.60$
${\tilde{k}}_5^0$	$1.0\cdot {10}^{-5}$	$\mathrm{mol}\cdot s^{-1}\cdot {\mathrm{cm}}^{-2}$	$c_5$	$11.60$
${\tilde{k}}_6^0$	$5.0\cdot {10}^{-3}$	$\mathrm{cm}\cdot s^{-1}$	$c_6$	$-19.33$
${\tilde{k}}_7^0$	$5.0\cdot {10}^{-6}$	$\mathrm{mol}\cdot s^{-1}\cdot {\mathrm{cm}}^{-2}$	$c_7$	$-19.33$
${\tilde{k}}_8^0$	$5.0\cdot {10}^{-6}$	$\mathrm{mol}\cdot s^{-1}\cdot {\mathrm{cm}}^{-2}$	$c_8$	$0$
${\tilde{k}}_9^0$	$5.0\cdot {10}^{-5}$	$\mathrm{cm}\cdot s^{-1}$	$c_9$	$19.33$
${\tilde{k}}_{10}^0$	$5.0\cdot {10}^{-8}$	$\mathrm{mol}\cdot s^{-1}\cdot {\mathrm{cm}}^{-2}$	$c_{10}$	$-19.33$

Table 4. (b) Values of the equilibrium redox potentials E ρ 0 superscript subscript 𝐸 𝜌 0 E_{\rho}^{0} and further constants entering the transition rates W ρ subscript 𝑊 𝜌 W_{\rho} .

$E_{\rho }^0$	Value	Parameter	Value
$E_2^0$	$0.8V$	$C_{\text{ce}HO}^b$	$0.7\cdot {10}^{-3}\mathrm{mol}\cdot {\mathrm{cm}}^{-3}$
$E_3^0$	$0.4V$	$C_{\text{ce}H}^b$	$0.3\cdot {10}^{-3}\mathrm{mol}\cdot {\mathrm{cm}}^{-3}$
$E_4^0$	$-0.15V$	$C_{\text{ce}X}^b$	$5\cdot {10}^{-4}\mathrm{mol}\cdot {\mathrm{cm}}^{-3}$
$E_5^0$	$-0.15V$	$N_s$	$2.2\cdot {10}^{-9}\mathrm{mol}\cdot {\mathrm{cm}}^{-2}$
$E_6^0$	$-0.32V$	$C_{dl}$	$2.0\cdot {10}^{-5}F\cdot {\mathrm{cm}}^{-2}$
$E_7^0$	$-0.32V$	$T$	$300K$
$E_9^0$	$-0.23V$	$\alpha$	0.5
$E_{10}^0$	$-0.23V$

Table 5. (c) Initial values of the dynamical variables used in the simulations of the \ce ​ H ​ 2 ​ O ​ 2 \ce 𝐻 2 𝑂 2 \ce{H2O2} reduction system.

Variable	$t=0$
$\varphi$	$-0.09$ $V$
$N_{\text{ce}H}$	$0.1\Omega$
$N_{\text{ce}OH}$	$0.1\Omega$
$N_{\text{ce}{H}^{ot}}$	$0.4\Omega$
$N_{\text{ce}X}$	$0.4\Omega$

Equations

$$k_{\rho }(\varphi )={\tilde{k}}_{\rho }^0{\mathrm{e}}^{c_{\rho }(\varphi -E_{\rho }^0)}=k_{\rho }^0{\mathrm{e}}^{c_{\rho }\varphi }$$

$$c_{\rho }=\frac{({\beta }_{\rho }-\alpha )|n_{\rho }|F}{RT}$$

$$\frac{\mathrm{d}P(\vec{N},t)}{\mathrm{d}t}=\sum _{\rho =1}^r\left[W_{\rho }(\vec{N}-{\nu }_{\rho },k_{\rho })P(\vec{N}-{\nu }_{\rho },t)-W_{\rho }(\vec{N},k_{\rho })P(\vec{N},t)\right]$$

$$c_{dl}\frac{\mathrm{d}\varphi }{\mathrm{d}t}=-i_F+\frac{U-\varphi }{RA}$$

$${\varphi }_{j+1}={\varphi }_j-\frac{n_{\rho }\mathrm{e}}{c_{dl}A}+\frac{U-\varphi }{c_{dl}RA}{\tau }_{\rho }$$

$${\tau }_{\rho }=\frac{R{c}_{dl}A}{c_{\rho }(U-{\varphi }_j)}\ln \left[1+\frac{c_{\rho }(U-{\varphi }_j)}{W_{\rho }^{0}R{c}_{dl}A{\mathrm{e}}^{c_{\rho }{\varphi }_j}}\ln \left(\frac{1}{r_{\rho }}\right)\right]$$

$${\tau }_{\rho }=\frac{1}{W_{\rho }^0}\ln \left(\frac{1}{r_{\rho }}\right)$$

$$W_{\rho }(\vec{N})=W_{\rho }^0{\mathrm{e}}^{c_{\rho }\varphi }$$

$$\text{ce}$$

$$\text{ce}(Pt-OH)+H++e-$$

$$\text{ce}2(Pt-OH)$$

$$\text{ce}Pt+H++e-$$

$$\text{ce}2(Pt-H^{ot})$$

$$\text{ce}Pt+X-$$

$$\Omega =A{n}_s{N}_A$$

$$\mathrm{SN}(x)=2\phi (x)\Phi (\alpha x)$$

$$\mathrm{SN}(x)=\frac{1}{\omega \sqrt{2\pi }}\exp \left[\frac{(x-\chi {)}^2}{2{\omega }^2}\right]\cdot \left[1+\mathrm{erf}\left(\alpha \frac{x-\chi }{\omega }\right)\right]$$

$$\mathrm{E}(x)=\chi +\omega \delta \sqrt{2/\pi }$$

$$\mathrm{var}(x)={\omega }^2(1-2{\delta }^2/\pi )$$

$$\zeta (x)=\frac{4-\pi }{2}\frac{\mathrm{E}(x{)}^3}{\mathrm{var}(x{)}^{3/2}}$$

$$\ln \frac{\tau }{⟨T⟩}=-\ln \frac{{\sigma }^2}{{⟨T⟩}^2}+\ln \frac{1}{2{\pi }^2}$$

$$\frac{\mathrm{d}{p}_{\rho }(\vec{N},t)}{\mathrm{d}t}=-W_{\rho }(\vec{N}+{\nu }_{\rho }|\vec{N})p_{\rho }(\vec{N},t)$$

$${\varphi }_{dl}(t+\Delta t)-{\varphi }_{dl}(t)=\frac{U-{\varphi }_{dl}(t)}{RA{c}_{dl}}\Delta t$$

$$P_{\rho }({\tau }_{\rho })=\exp \left(-{\int }_{t_n}^{t_n+{\tau }_{\rho }}{W}_{\rho }\left(\vec{N}+{\nu }_{\rho }|\vec{N}\right)\mathrm{d}t\right)=\exp \left(-W_{\rho }^0{\int }_{t_n}^{t_n+{\tau }_{\rho }}\exp \left(c_{\rho }{\varphi }_{dl}\right)\mathrm{d}t\right)$$

$$P_{\rho }({\tau }_{\rho })=\exp \left[-W_{\rho }^0{\int }_{t_n}^{t_n+{\tau }_{\rho }}\exp \left(c_{\rho }{\varphi }_{dl}(t_n)+c_{\rho }\frac{U-{\varphi }_{dl}(t)}{RA{c}_{dl}}(t-t_n)\right)\mathrm{d}t\right]$$

$$P_{\rho }({\tau }_{\rho })=\exp \left[-\frac{W_{\rho }^{0}RA{c}_{dl}\exp \left(c_{\rho }{\varphi }_{dl,j}\right)}{c_{\rho }(U-{\varphi }_{dl,j})}\left(\exp \left(c_{\rho }\frac{U-{\varphi }_{dl,j}}{RA{c}_{dl}}{\tau }_{\rho }\right)-1\right)\right]$$

$${\tau }_{\rho }=\frac{R{c}_{dl}A}{c_{\rho }(U-{\varphi }_{dl,j})}\ln \left[1+\frac{c_{\rho }(U-{\varphi }_{dl,j})}{W_{\rho }^{0}R{c}_{dl}A\exp \left(c_{\rho }{\varphi }_{dl,j}\right)}\ln \left(\frac{1}{r_{\rho }}\right)\right]$$

Full text

1 Introduction

Self-sustained oscillations in chemical systems are a prominent manifestation of self-organization under non-equilibrium conditions Nicolis:1995 . They occur in laboratory systems as well as in nature, the circadian clock governing pivotal functions of our body being just one prototype example Goldbeter2002 . The small volume of single cells, in which oscillations might take place, as well as nano-sized catalytic surfaces that support oscillatory surface reactions Kruse:2010 ; Barroo:2014 naturally lead to the question of the impact of molecular noise on the oscillatory behavior, or, in other words, of the minimum number of molecules necessary to support correlated oscillations. This question has been studied theoretically by Gaspard Gaspard2002 . For the weak noise limit, he derived an expression of the correlation time in terms of the extensivity parameter, the oscillation period, and the derivative of the period with respect to the Hamilton-Jacobi pseudoenergy. For a typical chemical oscillator this analysis predicts the lower limit of the number of molecules to be between 10 and 100. This result could be validated in recent experimental studies on the oscillatory $\mathrm{NO_{2}+H_{2}}$ reaction on nano-sized Pt tips Barroo:2015 . Here, the product of the correlation time of the oscillations and the variance of the period of the oscillations was shown to be constant.

Electrochemical oscillators differ from isothermal chemical oscillators because, in general, the electrostatic potential of the electrode is an oscillating quantity, in addition to the chemical variables Krischer:1999 ; Krischer:2002 . For nano- or mesoscale electrodes that are resistively coupled to some conducting support each stochastically occurring discrete electron transfer event changes the electrostatic potential of the electrode. At the nanoscale, the electrode potential is thus also a time-fluctuating quantity VGM1 ; VGM2 and so are the reaction rate coefficients, which depend exponentially on the electrode potential. Moreover, fluctuations of the electrode potential enhance each electrochemical reaction rate compared to the macroscopic limit. This peculiar behavior is closely connected to the potentiostatic control of electrochemical reactions, where the potential difference between the working electrode, i.e., the electrode of interest - in our case the nanoscale electrode - and a reference electrode is held constant. This operation mode imposes a constraint on the evolution of the electrode potential and may lead to non-normal distributions of the variables as well as shifts of the average steady state current densities compared to the macroscopic case VGM2 . The dependence of the rate constant of the electrochemical oscillations on a time-fluctuating quantity complicates theoretical studies on nanoscale electrochemical oscillators considerably. The dynamics of the nanoscale system is governed by an electrochemical master equation (EME) with time-dependent transition rates for the electrochemical reactions and an evolution equation for the electrode potential coupled to the one of the probability distribution of the chemical concentrations. The EME can be solved using an extended Gillespie algorithm VGM1 ; JansenMC .

In this paper, we investigate electrochemical oscillations on nanoscale electrodes, and compare them with chemical oscillations in small systems. We do so with simulations of a prototypical potentiostatic oscillator, the reduction of hydrogen peroxide on Pt electrodes in the presence of halides Mukouyama:2001 . This oscillator belongs to the large class of N-shaped negative differential resistance oscillators. In particular, we determine a critical electrode size below which correlated oscillations do not exist. Furthermore, we discuss the transition of uncorrelated to correlated behavior and test the validity of theoretical results for nano-scale chemical oscillators in the weak noise limit for electrochemical oscillators in regimes where the probability distribution of the period of the oscillations attains a Gaussian shape.

2 Electrochemical master equation for the oscillatory reduction of $\mathrm{H_{2}O_{2}}$ on Pt

2.1 Electrochemical master equation

An electrochemical reaction occurs on the surface of an electrode that is immersed in an electrolyte and involves the transfer of an electron from a species dissolved in the electrolyte to the electrode (oxidation reaction) or from the electrode to a species in the electrolyte (reduction reaction). These electron transfer reactions might be preceded, accompanied or followed by further elementary steps analogous to the ones occurring in a heterogeneous surface reaction, i.e. diffusion to or from the electrode, adsorption on the electrode surface, surface diffusion, reactions with other species and desorption from the electrode surface. Due to the electron transfer steps the rate of an electrochemical reaction depends on the electrostatic potential of the electrode. According to the Butler-Volmer Equation, the potential dependence of a reaction $\rho$ can be expressed through a potential dependent rate constant $k_{\rho}$:

[TABLE]

Here $\phi$ is the electrode potential (with respect to a given reference electrode), $\widetilde{k}_{\rho}^{0}$ the rate constant at equilibrium potential $E^{0}_{\rho}$. $c_{\rho}$ is given by:

[TABLE]

with $\alpha$ the transfer coefficient, $n$ the number of transferred electrons, $F$ the Faraday constant, ${R}$ the ideal gas constant and $T$ the absolute temperature. The value of $\beta_{\rho}$ changes with the form of the reaction $\rho$: $\beta_{\rho}=0$ for reduction reactions and $\beta_{\rho}=1$ for oxidation reactions. Note that for $n_{\rho}=0$, it follows that $c_{\rho}=0$ and the reaction is a purely chemical reaction.

As long as the electrode potential $\phi$ does not change with time, e.g. due to a perfectly conducting connection of the nanoelectrode and the outer electric circuit, $k_{\rho}$ is a constant, and the changes in the number of molecules of the chemical species involved in the reactions are goverend by the chemical master equation:

[TABLE]

Here, $P(\vec{N},t)$ is the probability that the system is at time $t$ in state $\vec{N}=(N_{1},N_{2}...)$ where $N_{i}$ is the number of molecules of species $i$. $W_{\rho}(\vec{N},k_{\rho})$ is the stochastic rate or transition probability of reaction $\rho$, which is expressed in terms of the macroscopic rate constant $k_{\rho}$. Each probabilistic reaction event $\rho$ changes the $N_{i}$s according to their stoichiometric coefficients $\nu_{\rho,i}$ which are equal to the difference of the number of products and reactants of species $i$ in reaction $\rho$.

However, in general $\phi$ changes with time. This is in particular true for all electrochemical oscillators which rely on a direct or hidden negative differential resistance in the current $I-\phi$ curve where $\phi$ is the positive feedback variable Krischer:1999 ; Krischer:2002 . Figure 1 depicts the equivalent circuit of an electrochemical cell. In the deterministic limit the evolution equation for $\phi$ is obtained by applying Kirchhoff’s circuit laws to the equivalent circuit:

[TABLE]

where $c_{dl}$ represents the double layer capacitance per unit area, $i_{F}$ the faradaic current density and $A$ the electrode area. In a potentiostatically controlled experiment the externally fixed potential $U$ and the ohmic resistance $R$ are kept constant.

For meso- or nanoscale electrodes, a faradaic current flows each time an electron transfer reaction takes place. If we assume that after a time $\tau_{\rho}$ $n_{\rho}$ electrons are transferred to the electrode in a reaction $\rho$, the faradaic current density flowing due to this reaction event is given by $n_{\rho}e/A\tau_{\rho}$. Furthermore, since the time difference between two reaction events $\tau_{\rho}$ is very small, the change of the potential due to the potentiostatic control can be expressed as $\frac{U-\phi}{c_{dl}RA}\tau_{\rho}$. The macroscopic evolution equation for $\phi$ (4) thus transforms to the following time-discrete equation for nanoscale electrodes:

[TABLE]

The time dependence of $\phi$ renders $k_{\rho}$ time-dependent, which we will call from now on rate coefficient. Furthermore, due to the potentiostatic control, the potential changes also in between two reaction events, impeding the use of the classical Gillespie algorithm Gillespie:1976 ; Gillespie:1977 for solving the master equation (3). In particular, the formula for the waiting time between two reaction steps of the classical algorithms cannot be employed here. Instead, a generalized approach for the waiting times $\tau_{\rho}$ for each reaction $\rho$ to occur has to be employed VGM1 ; JansenMC .

[TABLE]

where $r_{\rho}$ is a random number between 0 and 1 that is associated with the probability of reaction $\rho$ to occur. The derivation of eq. (6) as well as a summary of the extended Gillespie algorithm employed in our simulations is given in the Appendix. Note that for $c_{\rho}=0$ $\tau_{\rho}$ takes on the value of the expression of the original Gillespie algorithm:

[TABLE]

Furthermore, for convenience we express $W_{\rho}$ through a time-independent, chemical part, $W_{\rho}^{0}$, and a potential dependent term:

[TABLE]

The dynamics of an electrochemical reaction network on a nanoscale electrode is thus determined by eqs. (3) and (5), which are coupled by the potential-dependent transition rates eq. (8) and $\tau_{\rho}$ given by eq. (6).

2.2 Oscillatory hydrogen peroxide reduction on Pt

For our further analysis we will consider a particular oscillatory electrochemical reaction network, the hydrogen peroxide reduction in the presence of halide ions. This system was extensively studied by Nakato’s group both experimentally and through mean field simulations Mukouyama:2001 . The oscillation mechanism involves the following reaction steps:

[TABLE]

The notation $(\ce{Pt-M})$ represents the adsorption of a given molecule $\ce{M}$ on a platinum surface site, Pt. Reactions (10) and (14) are purely chemical, while all others involve an electron transfer step. The adsorption of upd (under potential deposition) hydrogen (forward reaction of eq. (12)) inhibits the adsorption (eq. (9)) and thus also the reduction (eq. (10)) of $\ce{H_{2}O_{2}}$, leading to a negative differential resistance (NDR) in the $I-\phi$ curve. Adsorption of halide ions (forward reaction of (15)) also inhibits $\ce{H_{2}O_{2}}$ adsorption but has an opposite potential dependence; $\ce{X-}$ desorbs when moving to negative potentials. Since the ad- and desorption dynamics of $\ce{X-}$ is much slower than the one of H-upd formation, the NDR is hidden for slow potential changes while it dominates the impedance at higher frequency. The oscillator is thus a classical ’hidden N-shaped negative differential resistance’ (HN-NDR) oscillator Krischer:1999 ; Krischer:2002 . Note that in contrast to the original model in Mukouyama:2001 we neglect transport limitation of $\ce{H_{2}O_{2}}$ as well as the dependence of $k_{1}$ on $\ce{OH}$ and $\ce{X}$ coverage for the sake of simplicity.

The transition rates $W_{\rho}^{0}$ corresponding to reactions (9)-(15) are compiled in table 1, the rate constants and other parameters used in the simulations are given in tables 2(a) and 2(b), respectively, the initial values used in all simulations in table 2(c). The extensivity parameter $\Omega$ is the absolute number of surface sites on the electrode surface:

[TABLE]

Here, $N_{A}$ is the Avogadro constant and $n_{s}$ represents the active site density (given in $\frac{1}{mol\>cm^{2}}$).

The bifurcation diagram of the macroscopic system is shown in fig. 2 in the resistance-applied voltage parameter plane. It was calculated from the deterministic rate equations for the evolution equations for $\phi$ (eq. (4)) and the coverages of upd-H, OH, X, and Hot (’ot’ standing for ’on top’) using AUTO Doedel:1986 . Note that both $\phi$ and $U$ are given versus the saturated calomel electrode (SCE) throughout the paper. The evolution equations for the chemical species can be derived from the transition probabilities in table 1 taking the limit $\Omega\rightarrow\infty$. The solid and dashed lines in fig. 2 show the locations of Hopf and saddle node bifurcations, respectively. Thus, very roughly the system is oscillatory for parameter values between the Hopf line and the right saddle node line. The mesoscopic simulations we discuss in the following were done at $RA=1.4\>\mathrm{\Omega\cdot cm^{2}}$ for different potential values along the red (dark gray) line in fig. 2. At this resistance value the oscillatory region is bounded by a saddle loop bifurcation and a Hopf bifurcation at the low and high voltage borders, respectively. When traversing the oscillatory region from the high to the low voltage border the oscillations change from quasi-harmonic to spike-like relaxation oscillations.

3 Simulations and Discussion

In the following we discuss results obtained from stochastic simulations of the electrochemical master equation for the oscillatory $\ce{H2O2}$ reduction with the extended Gillespie algorithm. Time series of the electrode potential $\phi$ for four different system sizes $\Omega$ but otherwise identical parameters are depicted in Figure 3. It is striking that the dynamics changes qualitatively with system size. From top left to bottom right the system size is increased, the lower right time series, obtained for $\Omega=20000$, being very close to the macroscopic oscillations. For the smallest value of the extensivity parameter shown, $\Omega=1000$ (top left), the system fluctuates around some mean value, however, without exhibiting any variations that would remind of an oscillatory excursion. Rather, it appears that the fluctuations are around a steady state. For intermediate system sizes, the system is still oscillatory, but the period of the oscillations is clearly larger than in the macroscopic limit.

A change of the average period with system size does not occur for purely chemical oscillators Gaspard2002 , which we also verified with the Brusselator model. We therefore attribute it to the peculiar property of electrochemical nanoscale systems that the reaction rate coefficients are fluctuating in time. It was shown that this renders each individual electron transfer event faster on a nanoscopic than on a macroscopic electrode, whereby the slowest steps are enhanced strongest VGM1 . Hence, also the relative rates within the reaction network are changed, and we conjecture that decreasing the system size drives the system effectively through the saddle loop bifurcation which occurs at more negative voltages in the macroscopic system. Figure 4 depicts a projection of three oscillatory trajectories obtained for different system sizes $\Omega$ in the $\phi-\theta_{\ce{X}}$ phase plane, $\theta_{X}$ being the coverage with halide ions. The continuous red line represents the macroscopic oscillation, the blue (dark gray) and green (light gray) lines are the trajectories in smaller systems. It is obvious that for the smaller systems the oscillations are deformed at low $\phi$ values towards higher $\theta_{\ce{X}}$ values, giving further evidence that for small system sizes the trajectories do not just fluctuate around the macroscopic limit cycle but occupy different regions in phase space.

To quantify the trend of the mean oscillation period with $\Omega$, the fluctuating time series had been smoothed and the average return time of the resulting smooth oscillating time series (s.b.) determined. To this end, the return value was set above the noisy slow-passage part of the time series on the rising flank. More precisely, it was set to half the value between the maximum and minimum potential value of the oscillations. In some circumstances, especially for the more harmonic oscillations and small system sizes, this lead to multimodal distributions since here the amplitude fluctuations were so large that not all oscillations reached the set return value. In these cases, a smaller value was chosen by inspection such that the distributions became unimodal. The thus determined average return time, or mean period of the oscillations, $T_{\Omega}$, normalized to the period of the macroscopic limit cycle, $T_{\infty}$, is plotted vs. the system size $\Omega$ in fig. 5 for four different values of $U$. It can be clearly seen that the qualitative trend of an increase in mean period with decreasing system size is independent of the value of $U$. However, quantitatively, the increase is more pronounced close to the saddle loop bifurcation.

Next, let us investigate the robustness of the oscillations with respect to molecular noise. From Figs. 3 and 4 it is obvious that the noise level in the time series, especially around the slow branch of the oscillations, increases with decreasing system size, as one would expect. First and foremost, molecular noise should manifest itself in a spread of the return times stemming from diffusion along the trajectories, i.e. so-called phase diffusion Kurrer:1991 ; Vance:1996 ; Gaspard2002 . Figure 6 displays four distributions of the return times obtained for different systems sizes and $U=-0.30\mathrm{V}$ together with short parts of the respective time series. The distributions of the three smallest systems exhibit asymmetric tails towards higher periods and are thus clearly non-Gaussian. For the two smallest values of $\Omega$ the tails are so pronounced that it does not seem to be reasonable to assign an average or mean period to the time series. This interpretation is also supported by the corresponding autocorrelation functions (not shown) that fall off so quickly that the first maximum lies below $1/\mathrm{e}$. In addition, the time series for $\Omega=1500$ reminds more of a noisy excitable system that fluctuates around a steady state where from time to time the fluctuations exceed a threshold such that an excursion in phase space is triggered. Owing to the nearby saddle loop bifurcation in the deterministic system and the excitable character of the global steady state beyond the saddle loop bifurcation together with the drift of the mean dynamics in phase space with system size, this scenario seems to be indeed conceivable. For increasingly larger systems, the distances between two successive oscillations are less spread, the mean and variance of the distribution being well defined from $\Omega=3000$ on.

For further evaluation, we fitted skew-normal distriubtions $\mathrm{SN}(x)$ to the data sets $x$, which can be represented by twice the product of a standard normal distribution $\varphi(x)$ and its relative cumulative distribution function $\Phi(\alpha x)$, the parameter $\alpha$ describing the asymmetry Azzalini:2013 .

[TABLE]

Considering the so called location $\chi$ and scale $\omega$, which correspond to the mean $\mu$ and the standard deviation $\sigma$, respectively, in the standard normal distribution, and using the fact that the cumulative distribution function depends linearly on the errorfunction $\mathrm{erf}$ of the dataset $x$, we can write the skew-normal distribution as follows:

[TABLE]

Fitting the distribution of the return times with this expression allows us to determine the mean value $\mathrm{E}(x)$, the variance $\mathrm{var}(x)$ and the skewness $\zeta(x)$ from:

[TABLE]

where $\delta=\alpha/\sqrt{1+\alpha^{2}}$.

The blue (dark gray) curves in the right column of fig. 6 represent the best fits of eq. (56) to the distributions. The distribution for $\Omega=3000$ has only a minor contribution of the skewness and can be considered Gaussian, the one for $\Omega=2000$ can well be fitted by a skew-normal distribution, though with a non-negligible skewness of about 0.46, while in the other two cases the distributions are far from being describable by a skew-normal distribution. Thus, in these cases it is also not possible to determine a variance or standard deviation. Note also that whenever we talk about a mean return time or mean period of the nanoscale system that has an assymetric distribution, we refer to the maximum of the distributions, not to the average period.

Figure 6 suggests that the loss of correlation of the oscillations for the spike-like oscillations at $U=-0.30\;\mathrm{V}$ occurs through a lengthening of the interspike distances whereby the Gaussian distribution transforms first to a skew-normal distribution which is then further distorted towards longer distances until the spikes occur essentially uncorrelated. At still smaller system sizes the system fluctuates around a low potential value without exhibiting large amplitude spikes anymore (cf. fig. 3).

This picture changes when the macroscopic base oscillations are more harmonic. Here, the probability distribution of the return times also progressively deviates from a Gaussian distribution when decreasing the extensivity parameter beyond a threshold value, but it stays approximately symmetric with respect to the maximum, exhibiting long-tails towards longer and shorter periods instead (fig. 7). This difference also shows up in the time series. With decreasing system size they exhibit continuously higher noise levels and more strongly fluctuating amplitudes. Also the absolute system sizes at which coherent oscillatory behavior breaks down change considerably with the applied voltage, being much larger for the more harmonic oscillations at less negative potentials than for the spike-like relaxation oscillations at the negative potential border of the oscillations. This can be seen in fig. 8 where the standard deviations normalized to the mean period of those probability distribution functions (pdfs) that could be well fitted by a SN distribution are depicted for $U=-0.26\;\mathrm{V}$ and $-0.30\;\mathrm{V}$. The relative standard deviation of the return times increases much faster with decreasing system size for $U=-0.26\;\mathrm{V}$ than for $U=-0.30\;\mathrm{V}$.

With a given lattice constant $a$, or average distance between two neighboring atoms, the extensivity parameter $\Omega$ can be related to an electrode size. Choosing $a=0.27$ $\mathrm{nm}$ as a typical value, $\Omega_{min}=2000$ as lower limit for correlated electrochemical oscillations at $U=-0.30\mathrm{V}$ corresponds to an electrode area of about 12 x 12 nm2 for flat electrodes or to spherical electrodes with a radius of about 3.5 nm. At $U=-0.26\mathrm{V}$ the lower limit of $\Omega$ supporting regular oscillations was $4000$, which in turn translates to 17 x 17 nm2 flat electrodes or spherical electrodes of 4.9 nm in radius, correspondingly. Compared to sizes reported in literature for heterogeneous surface reactions, these numbers are considerably larger. Also the number of molecules needed for a chemical clock, which was estimated to be between 10 and 100 Gaspard2002 ; Barroo:2015 , is close to 2 orders of magnitude smaller than for our electrochemical oscillator. For an electrochemical oscillator, the latter can be roughly estimated to be 1 to 2 times the value of $\Omega_{min}$, assuming that a typical oscillator has two to three types of adsorbed species with average coverages of 0.3 - 0.8. This implies that the noise level in electrochemical oscillators, where also the electrode potential fluctuates in time, is considerably larger than in chemical oscillators.

Despite the strong dependence of the quality of a nonlinear oscillator on the characteristics of the limit cycle in the macroscopic system, theory predicts that in the low-noise limit where the chemical master equation reduces to a Fokker-Planck equation the product of the relative correlation time $\tau/\left\langle T\right\rangle$ and the relative variance of the first return times $\sigma^{2}/\left\langle T\right\rangle^{2}$ should be equal to $1/2\pi^{2}$ Barroo:2015 . For an electrochemical nanoscale system, the time-dependence of the rate coefficients $k$ hampers the derivation of a Fokker-Plank equation, and in fact makes it impossible to follow the same lines as for the a chemical system VanKampen . Yet, it is possible to test this conjecture numerically. Therefore, we extracted the variance and mean period from the probability distributions of the period and calculated the autocorrelation function $C\left(t-t^{\prime}\right)$ from the fluctuating time series. According to the theory $C\left(t-t^{\prime}\right)$ should exhibit damped harmonic oscillations with an exponentially decreasing envelope, $\mathrm{e}^{-(t-t^{\prime})/\tau}$, where the relative correlation time is the decay rate $\tau$ devided by the mean period $\left\langle T\right\rangle$. In the case of more harmonic oscillations and Gaussian pdfs, i.e. towards the high-potential border of the oscillatory region, the autocorrelation functions exhibited indeed exponentially damped, nearly harmonic oscillations. However, in case of the more spike-like oscillations towards lower potentials, the damped oscillations of the autocorrelation functions were neither harmonic in shape nor was their envelope exponentially decreasing. Rather, the decay of the maxima of the autocorrelation function was stronger than exponential, also in the cases of large $\Omega$ where the pdfs were Gaussian. In this case, we determined the delay time at which the envelope of the autocorrelation function attained the value $1/\mathrm{e}$ by linear interpolation between the corresponding maxima of $C\left(t-t^{\prime}\right)$ and used this delay time in the further evaluation.

A double logarithmic plot of the relative decay time and the relative variance of the period of the oscillations is shown in fig. 9 together with the theoretical predicted curve Barroo:2015 given by:

[TABLE]

which has a slope of $-1$ and an intercept of $\ln(1/2\pi^{2})\approx-2.98$. The data were obtained for two different $U$ values and system sizes that all have Gaussian distributions of the period. Clearly, the data of the two sets deviate significantly from the theoretical line. Furthermore, they separate in two groups corresponding to the $U$ values for which they were obtained. For $U=-0.30\mathrm{V}$ the best linear fit gives a slope of $-1.26\pm 0.05$ and an intercept of $-5.48\pm 0.35$. This means that at a given relative variance of the period the relative correlation time is shorter than predicted, which probably reflects the faster than exponential decay of $C\left(t-t^{\prime}\right)$ and the non-harmonic form of the oscillations. The situation is somewhat different for the more harmonic oscillations at $U=-0.26\mathrm{V}$. Here, the data points obtained for the 5 largest system sizes between $\Omega=20000$ and $12000$ lie on a line with a slope of $-1.03$ and an intercept of $-3.26$, very close to the values anticipated for the weak noise limit of chemical oscillators. However, for small system sizes the increase in relative variance is much stronger than predicted for the respective relative correlation times. Here, clearly the theory does not capture the observed behavior. We also determined the corresponding points for an intermediate value of $U$, where the data points showed a trend somewhat between the two displayed cases. Taken together, the data give evidence that the universal relationship between phase diffusion and correlation time of chemical oscillations is only applicable in very small parameter regions for electrochemical oscillations. At this point we can only speculate that the origin of this different behavior is a consequence of two properties of electrochemical systems, the stochastic time-dependence of the reaction rate coefficients, and the constraint on the dynamics introduced by the potentiostatic operation mode. For further insight more analytical work seems to be necessary. From this perspective, a break through seems to be linked to the success to derive a Fokker-Planck equation for electrochemical systems.

4 Conclusions

In conclusion, we demonstrated that electrochemical nanoscale oscillators feature some peculiar properties which render them less robust to molecular noise than chemical oscillators. The mean period of electrochemical oscillations on mesoscale and nanoscale electrodes depends on the system size, the mean trajectories of the oscillatory time series following different paths in phase space, and close to the macroscopic oscillation border the trajectories are driven out of the oscillatory region, thus confining oscillations to a smaller parameter range. Furthermore, the correlation time decreases strongly with decreasing system size, which might limit the smallness of nanoelectrodes in technological applications. The minimal system size and thus the lower bound to the number of chemical species necessary for time-correlated oscillations is one to two orders of magnitudes larger than for chemical systems, the minimal surface area of the electrode amounting to 100 to 500 nm2. At sufficiently large system sizes the probability distributions of the first return periods are Gaussian, yet also in this region the universal relationship between correlation time and the variance of the period of the oscillations predicted for nanoscale chemical oscillations in the weak noise limit is not found in large parameter regions. Furthermore, for progressively smaller systems we found that the transition from coherent oscillations to noisy fluctuations lacking any defined mean period may occur through different scenarios. For spike-like relaxation oscillations the probability distribution of the next return map becomes more and more asymmetric towards larger periods before the correlation is lost completely, the oscillatory excursions of the smallest systems reminding of those in an excitable system. Thus, here the trajectories in phase space still follow rather defined paths, but at erratic times. For more smooth oscillations in the macroscopic limit, the probability distribution of the return times remains nearly symmetric but obtains long tails for system sizes below a threshold, the trajectories covering an increasingly larger region in phase space with decreasing $\Omega$ and thus the oscillatory cycles become less and less well defined. In the latter case, the system is also more prone to noise, coherent oscillations being lost already at larger system sizes than in the case of more relaxational oscillations.

The decisive difference between the thus far studied chemical systems and the electrochemical system is that in the latter the electrode potential constitutes a further dynamical variable that renders the reaction rate coefficients to depend on time. Any chemical reaction which is not run under isothermal conditions will also involve changes of the temperature. Since all reaction rates depend exponentially on the temperature, a fluctuating temperature causes fluctuations in the reaction rate coefficient, just as the potential does in electrochemical systems. Dust particles in the higher atmosphere at which catalytic reactions occur might be an important realization of such nano- or mesoscale nonisothermal chemical systems. From this perspective, it would be interesting to study the impact of molecular noise in nonisothermal chemical nanoscale oscillators. They might exhibit some of the features we found for electrochemical systems, which would then be of more general nature.

However, also the external control of electrochemcial systems enters into the dynamics in a nontrivial way. The potentiostatic operation mode introduces an external driving which renders all elemenatry reaction steps to be faster on a nanoelectrode than in the deterministic limit VGM2 . A deeper understanding of the peculiar features of the nano-oscillators discussed here requires further studies that allow the impact of the external driving to be separated from that of the time-dependent reaction rate coefficients.

Appendix

Extended Gillespie Algorithm

Derivation of Eq. (6)

The derivation given here follows closely the one given in the supplementary information of [9].

The probability $p_{\rho}(\vec{N},t)$ that a reaction $\rho$ has not yet occurred after a time $t$ is given by the solution of the following differential equation

[TABLE]

According to eq. (5), in a small time interval $\Delta t$, in which no faradaic reaction occurs, the potential changes by

[TABLE]

The probability that in the time interval $\tau_{\rho}$ reaction $\rho$ has not occurred can thus be expressed as

[TABLE]

or, using eq. 59

[TABLE]

With the notation change to discrete potential steps $\phi_{dl}(t_{n})\rightarrow\phi_{dl,j}$, integration of eq. 61 yields:

[TABLE]

Relating the random number $r_{\rho}$ to the probability $P_{\rho}(\tau_{\rho})$, i. e. setting $P_{\rho}(\tau_{\rho})=r_{\rho}$, we can solve equation (62) for the waiting time $\tau_{\rho}$ and obtain:

[TABLE]

which is equal to eq. (6).

Gillespie Algorithm

The electrochemical master equation was solved using Gillespie’s generalized first reaction method. In detail the following steps were carried out:

An input for $\phi$ and the species $N_{i}$ was chosen. 2. 2.

The transition probabilities $W_{\rho}$ were calculated for the initial conditions. 3. 3.

A random number $r_{\rho}$ was generated from the uniform distribution in the unit interval for each reaction $\rho$. 4. 4.

Using eq. (6), a waiting time $\tau_{\rho}$ was calculated for each reaction $\rho$. 5. 5.

The equation with the smallest value of $\tau_{\rho}$ was chosen to advance. All $N_{i}$’s involved in the corresponding reaction are updated according to $N_{i}\rightarrow N_{i}+\nu^{i}_{\rho}$ and $\phi_{j}$ was updated according to eq. (5). 6. 6.

The new propensities were calculated and steps 3 - 6 were repeated until the simulation was stopped, usually after several thousands of oscillation periods.