Phase diagrams and dynamical evolution of the triple-pathway electro-oxidation of formic acid on platinum.

Recently, an electro-kinetic model based on a specified reaction scheme for the electro-oxidation of formic acid on platinum was reported. The model evaluated three reaction pathways towards the production of CO2: the dehydrogenation and the dehydration of formic acid, and the third and most active pathway includes fast oxidation of the formate ion. Numerical integrations showed that the model is well-suited to describe the experimental results in voltammetric and oscillatory regimes. In the present paper, we provide detailed stability phase diagrams characterizing the dynamical evolution of this system under galvanostatic and potentiostatic regimes. We find the triple-pathway electro-oxidation of formic acid on platinum to have rather intertwined stability phases and, surprisingly, a total absence of chaotic oscillations. To the best of our knowledge, this is the first study in this direction using a realistic electrochemical model.


Introduction
The electro-oxidation of organic compounds on catalytic surfaces has great interest from fundamental and applied perspectives. 1, 2 The formic acid electro-oxidation reaction (FAEOR) on platinum has been studied several decades ago, 3,4 but the reaction mechanism is not yet completely clear. 5 Recently, a plausible proposal was presented 6 for this electrochemical process, based on different experimental and theoretical observations reported in the literature. 5,[7][8][9][10][11][12][13][14] The reaction scheme proposed for the FAEOR is illustrated schematically in Fig. 1, and considers three reaction pathways towards the CO 2 production.
Each of the reaction steps in Fig. 1 are explicitly described by the following chemical equations: OH ad þ H þ þ e À (r5) CO L þ OH ad ! 6 CO 2 þ H þ þ e À þ 2 Ã (r6) In these equations, asterisks represent a free surface reaction site, the 'B' and 'L' subindices indicate that the chemical species are adsorbed on the electrode surface in a bridge-bonded or linear-bonded configuration, respectively.
The first pathway, the so-called formate pathway, is represented by three chemical eqn (r1)-(r3) for dehydrogenation of formic acid. [10][11][12][14][15][16][17][18][19][20][21][22] The indirect pathway (r4)-(r6) is the CO ad and OH ad formation and the further reaction via Langmuir-Hinshelwood mechanism. 11,12,19,[23][24][25][26][27][28][29][30][31] In this way, we considered a very important aspect in contrast with previous models. 14,21,25,32,33 Our proposal considers that the formation of CO L in platinum can also proceed from a monodentate adsorbed formate, that is, after the reaction steps indicated in eqn (r1) and (r2). This idea was initially suggested for the solid-gas interface 27 and supported by experiments in which a kinetic isotopic effect was attributed to the faster decomposition of HCOO ad than DCOO ad . 7 Recent experiments and calculations using density functional theory also support this idea. 11,12,19,28 In turn, this proposal allows us to consider that the (r4) dehydration process is initiated by an electroreduction process. This can also successfully explain the dependence on the potential of CO L formation. The third pathway, (r7), is the direct one. Unlike a previous model, 14 we consider that the active species in this third way is the anion formate which is quickly oxidized to CO 2 by a dehydrogenation process. 34,35 Some experimental observations corroborate the involvement of the solution formate in the oxidation process. 11,36 The electro-oxidation of organic compounds can be performed under conditions in which unusual and very interesting behaviors are observed. For example, behaviors in which electrochemical variables oscillate in time are rather common, and periodic, mixed mode and chaotic oscillations have been reported. Based on the proposed reaction scheme, we presented a renewed electrochemical model for FAEOR and compared the numerical results with experimental data. 6 The model consists of 4 + 1 differential equations, which will be presented in the next Section, corresponding to the mass balance of the concentration of the reactants on the electrode surface and the charge balance at the electrochemical surface.
The results of the numerical integration showed that the model is well suited to describe the experimental results under voltammetric and oscillatory regimes. Phase diagrams can provide a comprehensive and detailed description of the system's dynamics. Since they display features in relevant twoparameter sections and discriminate simultaneously regions with different types of periodic motions or regions where there is chaotic behavior. 37 These diagrams present a powerful tool to study nonlinear systems because they uncover the fine structure within the oscillatory regions, allows at classifying the dynamics in a broader context, and also permit the identification of emergent universalities found in disparate models. In this work, we report high-resolution stability diagrams in which parametric regions are classified by counting the number of spikes (local maxima) per period of sustained periodic oscillations. For details, see section Computational procedures below. The stability diagrams were computed for both galvanostatic and potentiostatic regimes. Previous studies in electrochemical systems were carried out with minimal and generic models, 37,38 in contrast, herein we use a realistic electrochemical model solidly built in experimental evidences. 6

Model equations
With the mechanism for the FAEOR described above, 6 we derived the following equations for the reaction rate of each step. During oxidative adsorption of formic acid (r1), formate adsorbs on the electrode surface in the so-called bridge-bonded, HCOO B , configuration, i.e., via two oxygen atoms occupying two platinum sites. Steps (r1) can thus be given as: where y vac 1 À y CO L À 2y CO B À y HCOO L À 2y HCOO B À y OH , C HCOOH = 1.99 Â 10 À4 mol cm À3 , F = 96485.33 C mol À1 , R = 8.314 J (mol K) À1 , T = 298.15 K, and k uni m is unitary but with the appropriate units for v m to be expressed in s À1 . The subscript L in HCOO L and CO L means that these species are occupying a single catalytic site on the electrode surface. Moreover, y CO B is the coverage of CO ad adsorbed in bridge-bonded occupying two catalytic sites. The model considered that the relationship between CO L and CO B is kept near equilibrium, which means 11 y CO B E 0.258y CO L .
After the initial step, the model considers the need for a configurational change of the HCOO B species 11 before the formation reactions of CO L or CO 2 by formate and indirect pathways. This configurational change corresponds to reaction (r2) that forms the HCOO L species in which the C-H bond is oriented with the hydrogen atom pointing to the metal surface. 12 Due to the presence of other neighboring adsorbed species, such as sulfate or other HCOOB molecules 39 a simple rotation around the O-CHO bond is difficult. The need for some neighboring catalytic sites for the interaction of the molecule lying down on the platinum surface is considered. The following equation describes the rate of step (r2): where k 2 = 110 s À1 . The CO 2 formation from the HCOO L species corresponds to a dehydrogenation process in which the hydrogen atom must interact with the electrode surface. This interaction requires a catalytic site that is considered in the reaction (r3): Experimental observations have been reported showing that the CO ad formation requires the presence of at least three contiguous Pt atoms for the dehydration of HCOOH to CO ad . 30 Considering that CO ad is formed from the species HCOO L , (r2) 11,12,19,28 which is already occupying a site, it is clear that two other active sites are involved in this process. Therefore, a reaction order 2 with respect to y vac is considered. Just as it was considered in the case of reaction (2), the HCOO L species must lie down on the electrode surface before or simultaneously to a process in which the adsorption mode should change from a Pt-O to a Pt-C binding mode. 11 The rate of step (r4) is thus Fig. 1 The recently proposed 6 reaction scheme for triple-pathway electro-oxidation of formic acid on platinum.
given as: Reactions steps 5 and À5 consider the oxygenated species formation on the surface of the electrode from water and its respective reverse reaction, and the corresponding rate equations read: where C H 2 O = 55 Â 10 À3 mol cm À3 . The kinetic equation for CO L oxidation to CO 2 , through Langmuir-Hinshelwood mechanism, is expressed as: Finally, the direct pathway considers that the formate in solution is quickly oxidized to CO 2 by a dehydrogenation process with the following rate law: where C HCOO À = 3.55 Â 10 À7 mol cm À3 , and the parameters f m and f Àm were fitted as indicated in Table 1. The exponential terms in the previous equations correspond to the electrochemical rate constants of each reaction step according to the Butler-Volmer theory, i.e. k m ðfÞ ¼ However, to decrease the number of kinetic parameters in the model, the standard rate constants k 0 m and the reaction standard potentials k 0 m have been grouped by the following substitutions: From arbitrarily chosen initial values, the parameters f m and f Àm were varied in intervals of 0.01 units until the best possible congruence between the numerical results and the experimental results (observed oscillations and the voltammetry response) was achieved. The fitting values of f m and f Àm are reported in Table 1.
The five variables that describe the dynamic of the electrooxidation process in our model are y HCOO B , y HCOO L , y CO L , y OH , and f. The core model consist of the equations These equations are supplemented by an additional equation, depending on the experiment under consideration. For galvanostatic experiments, eqn (12)-(15) are supplemented by the equation while for the potentiostatic experiments, instead of the above equation, the supplementary equation used is Here, C d = 24 Â 10 À6 C V À1 cm À2 , R S = 10 O, A = 0.42 cm 2 , and N tot = 2.18 Â 10 À9 mol cm À2 . In this work, j is a parameter that is varied between 0 and 0.002 A cm À2 , and R ext varies between 0 O and 3000 O. Regarding the first model, the current density is a parameter easy controllable in galvanostatic experiment and therefore it is thus an obvious parameter to be explored. Other interesting parameters that may have a determining role in the dynamic behavior of the electrochemical system are the f m terms. 6 Each of these parameters groups together two others: the standard (or formal) potential and the heterogeneous rate constant associated with the respective electrochemical reaction step. 6,14,21,25 Consequently, each f m can be a function that depends on the nature of the electro-catalytic surface, on the chemical environment, and on factors that affect the activation energy of the m-th reaction step. Thus, f m are not factors individually manipulated in traditional electro-chemical experiments. In other words, it is difficult to modify the rate of a single reaction step without altering the rate of the other reactions involved in the FAEOR. In this work, the parameter space j Â f 5 is presented under the premise that it is possible to modulate the intrinsic rate constant of CO L formation through f 5 . The choice of f 5 is justified by the leading role of y CO L in the positive feedback loop associated with HNDR in the pseudo-stationary voltammograms. 40 In potentiostatic experiments, the electrical resistance of the system is a parameter that plays a decisive role, as it gives the Table 1 Kinetic parameters f m and f Àm (in Volts) used in our simulations. In this work, f 5 varies between 0.35 V and 0.55 V in the diagrams below freedom to the interfacial potential differ from the applied one, allowing thus the emergence of oscillations. To modulate the resistance of the system as a bifurcation parameter, an external resistance can be experimentally connected in series with the electrochemical cell. Thus, the change in the value of the electrical resistance allows shifting the dynamics of the process to a region where bistability in the system is favored. 41 Hence, in the potentiostatic case, the isospike diagrams were constructed to represent the parameter space E Â R ext .

Computational procedures
Before presenting stability phase diagrams for the galvanostatic and the potentiostatic models, we briefly describe how such diagrams were obtained. The stability diagrams reported in this work are the so-called isospike diagrams, 37,42-47 namely diagrams discriminating not only periodicity from chaos but, in addition, providing a detailed cartographic representation of the number of spikes per period of the oscillations. Several illustrative examples of isospike diagrams are shown, e.g., in Fig. 3. Such phase diagrams were computed as follows. Given a parameter window of interest, we start by discretizing it into a N Â N grid of equally space points. Then, for each point of the grid, the equations of motion are integrated numerically using the standard fourthorder Runge-Kutta algorithm with fixed time step h = 5 Â 10 À5 . For simplicity, integrations are always started from the same fixed initial condition (y HCOO B ,y HCOO L ,y CO L ,y OH ,f) = (0,0,0,0,0). The first 30 Â 10 6 integration steps are discarded as a transient time needed to approach the attractor. Next, to compute the isospike diagrams, integrations are further continued for 30 Â 10 6 additional time-steps when up to 1600 extrema (maxima and minima) of the five variables are recorded. We always recorded maxima and minima during an interval of time identical to the transient time used. In some cases, to get clean diagrams it was necessary to increase transient and integration times from the above reference value. Below, the figure captions contain the actual values used. From such records of maxima and minima it is possible to determine for each parameter point whether pulses repeated or not. Isospike diagrams are not difficult to construct based on experimental time series. 48 The typical computer time needed to obtain a stability diagram with, say, 400 Â 400 parameter points is of the order of 180 hours when running on 50 processors of a relatively outdated (2011) SGI Altix ICE system with 16 compute blades within the individual rack unit. It is a cluster with 1536 processors and a theoretical peak performance of 16 Tflops. Unless stated otherwise, all phase diagrams of this paper display the analysis of 400 Â 400 parameter points.
All phase diagrams below use 17 colors to represent the number of spikes per period of the oscillations, as indicated by the colorbars in the figures. Oscillations that contain more than 17 spikes per period are plotted recycling modulo 17 the 17 basic colors. Black is used to represent lack of numerically detectable periodicity for a given integration step and transient time used in the Runge-Kutta integrator. To simplify diagrams, white is used to collectively represent divergent, undefined and nonoscillatory solutions, if any. The integration of differential equations for large sets of parameters and initial conditions is a quite demanding numerical task and can only be performed using computer clusters and suitable ad hoc programming.  Fig. 2 one sees that while the number of spikes per period may be easily counted in the panels displaying y HCOO L , y OH and f this is not the case for y HCOO B and y CO L . For these two variables, spikes are rather small and confined to narrow intervals. However, the number of spikes are not difficult to recognize by simply magnifying suitable segments of the periodic signals (not shown here). Signals are never zero but may reach values as small as 10 À4 . These variations in the coverage of the surface by adsorbed chemical species are only accessible through the modeling of the dynamic process. Experimentally, oscillatory patterns of the adsorbed formate coverage and adsorbed carbon monoxide coverage, calculated from the intensity of the IR bands, have been recorded in spectroelectrochemical experiments. 21 However, the current experimental resolution probably would not permit distinguishing those slight variations in the coverage of adsorbate. In contrast, the electrode potential f can be measured. Fig. 2 clearly shows that the number of spikes per period depends sensitively on the variable used to count them. So, it is natural to inquire if stability diagrams derived from such time series will look similar. An answer to this question is provided by the five panels in the top row of Fig. 3, which show parameter regions centered around the experimental values studied in ref. 6. These five panels shows that, indeed, diagrams obtained by counting spikes in y HCOO B and y CO L show greater similarities than the corresponding ones obtained from the spikes in y HCOO L , y OH and f. The second row from the top in Fig. 3 shows magnifications of details contained in the rectangular boxes seen on the panels on the top row. While phase diagrams for y HCOO L , y OH and f contain a rich alternation of stripes of different colors, diagrams for y HCOO B and y CO L display striation of a single color, corresponding to two spikes per period. As j increases, the initially continuous color of such striations starts to be mediated by the blue color of the background, corresponding to periodic oscillation having one spike per period.

Results and discussion
The mixed-mode oscillations depicted in Fig. 2 consist of at least for two types of oscillators with different frequencies.
In one case, the positive feedback loop can be assigned to the process in which an increase in potential from 0.56 V accelerates the production of y OH , which blocks active sites and causes further increase of the electrode potential. When the potential reaches the high limiting value (around to 0.85 V), the amount of y OH on the surface decreases, releasing active sites, and the potential decreases in order to maintain the current level set in the experiment, this process is defined as the negative feedback loop. The y OH participation in this oscillator was evidenced through the time series analysis, where large amplitudes are observed at higher frequency oscillations. This process has already been reported in the literature and associated with the origin of the NDR during the voltammetric experiments of the FAEOR. 31 The second oscillator mainly involves the variables y HCOO B and y CO L . In this case, the positive feedback loop can be rationalized as follows: an increase of the electrode potential accelerates the production of y HCOO B , and consequently the y HCOO L and y CO L also increase, reactions (r1), (r2) and (r4), respectively. Those species block the active sites on the surface and, in order to keep the current level, a further rise of the potential results. Then, when the oxidation of CO L at high potentials releases surface active sites, the negative feedback loop takes part. Once the sites are available, the oxidation reaction rates rise and the potential decrease. Hence, the amplitudes of the oscillations of y HCOO B and y CO L seem determined by the lower frequency oscillator.
The lower frequency oscillator determines the oscillation period while the higher frequency oscillator provides the number of peaks computed in the isospike diagrams, Fig. 3. Therefore, it is to be expected that the isospike diagrams of y HCOO B and y CO L are similar to each other but clearly different from that for y OH . Regarding y HCOO L and f, they are affected by the changes in the other variables, responding to the behavior of both oscillators. Thus, the number of spikes counted from those variables should correspond to the number of spikes produced by the higher frequency oscillator, i.e. the number of peaks for y HCOO B and f is close to the number of peaks for the variable y OH , which explains the similarity between those isospike diagrams. Now, before describing the mechanism underlying all these phase diagrams, we mention an interesting feature in Fig. 3 that, in fact, is a significant finding of our work. As defined in Section 3, our phase diagrams use black to represent solutions, which, although not divergent, lack any numerically detectable periodicity. This type of behavior is characteristic of 'chaotic' oscillations. Accordingly, the upper left corner in most panels  Table 2. of Fig. 3 seem to indicate periodic oscillations to end up abruptly into a phase of chaotic oscillations. That this is not the case is demonstrated by comparing the black phases seen in the three bottom rows of Fig. 3. The phase diagrams in the second row from the top were obtained after removing a transient of 30 Â 10 6 time steps, the third after 60 Â 10 6 steps, and the forth after 120 Â 10 6 . As it is not difficult to realize, the size of the black phases decrease continuously when the transient time removed increases. This tendency indicates that, rather than chaos, the black phase merely records parameter regions that require removing considerably longer transients before finally revealing that they correspond to periodic oscillations. Several additional ad hoc tests corroborated this to be the case. Chaotic oscillations were never observed in our simulations. Incidentally, we observe that similar chaos-free oscillations were recently reported in a rather different context, namely, in the low-frequency limit of externally driven systems: in a driven Duffing oscillator, a popular model of micro-and nanomechanical devices, and for the so-called Brusselator, 49 which were shown to display wide mosaics of stability phases arising from periodic oscillations of ever increasing periodicity, but with no chaos. Here, in contrast, we find chaos-free phases to arise in non-driven systems, and to emerge only in the form of stripes, not as mosaics. Fig. 4(a) shows a magnification of essentially the same box shown on the top row of Fig. 3, obtained by counting spikes of the y HCOO L variable. Basically, Fig. 4(a) shows a division of the control parameter window into three phases displaying characteristic behaviors. First, for low f 5 and high j the galvanostatic model displays divergent solutions. Second, for  Table 2. Black arrows mark location of small spikes. high f 5 and low j one sees a black phase, indicating that although numerical solutions do not diverge in this region, no periodic oscillations are observed. Third, between the two phases above, one finds a wide intermediary region, an 'oscillatory phase', characterized by a large number of colored stripes, indicating the presence of stable periodic oscillations containing an apparently unbounded number of spikes (local maxima). Fig. 4(b) shows a magnification of the box in Fig. 4(a), where one sees the continuation of the sequence of points c i . The stability diagram in Fig. 4(c) shows the same parameter window depicted in Fig. 4(a), but displaying the continuous variation of the period of the y HCOO L oscillations. Fig. 4(c) complements the information provided in Fig. 4(a) and conspicuously shows the division into the three aforementioned phases. Noteworthy in this figure is the rather fast increase of the period length towards the boundary with the black phase. As already mentioned, the transition to the black phase is simply pointing out an increase in the number of spikes per period and a decrease of their amplitudes. By adjusting the integration step and neglecting longer transient times we found that it was possible to decrease continuously and systematically the size of the black phase. As before, no chaos was ever found. All these tests consume large amounts of computer time. Therefore, we have not attempted to eliminate completely the black phase. At any rate, periodic oscillations in the black phase contain large amounts of fast-oscillating small-amplitude peaks that are likely to be out of reach in experiments.
In Fig. 4(a) there are two vertical sequences of points labeled i i and c i . Under the stability diagrams, Fig. 4 shows typical time series observed when increasing the value of f 5 while maintaining j constant. Such time series make plain the meaning of the stripes in Fig. 4(a): as f 5 increases, there is a smooth and regular complexification of the time series through the acquisition of more and more spikes. The mod 17 recycling of colors seems to indicate that their number of spikes may continue to increase continuously without bound. It is worth mentioning that the rate of CO ad formation on the electro-catalytic surface decreases when f 5 grows. This relationship increases the period of the oscillations and consequently a larger number of spikes is obtained.
Next, we investigate what happens when computing phase diagrams using a variable different from y HCOO L , used in Fig. 4. Experimentally, the variable f is directly accessible in conventional electrochemical experiments, while y HCOO B and y CO L can be measured through the coupling of in situ infrared spectrophotometer with the electrochemical setup. 21 The variable y OH , however, while an important component of the model, is rather difficult to study experimentally. Fig. 5 reproduces the same numerical experiment used to obtain Fig. 4, but now based on f (instead of y HCOO L ). A cursory  Table 2 Coordinates ( j,f 5 ), period T, and number of peaks p HCOO L and p f for fixed values of the current j, for points labeled c i , d i and i i in Fig. 4 and 5. Transients for i 8 and i 9 are 60 Â 10 6 and 90 Â 10 6 , respectively, and 30  comparison between the phase diagrams in Fig. 4 and 5 may give the impression that all three diagrams coincide.
However, a closer look shows that there are some differences among them. For instance, the divergent phase in Fig. 5 connect directly to the one-spike blue phase, not displaying the two-spikes phase which mediates these phases in Fig. 4.
Furthermore, although Fig. 4(b) and (b) display the same parameter window, their inner content is slightly different. As a consequence, in Fig. 5 we consider a series d i of points, distinct from the c i in Fig. 4. Fig. 5 contains arrows indicating the location of very small spikes. As illustrated by this sequence of time series, the arrows also indicate were  new spikes are added to the oscillatory patterns when f 5 increases.
Taking into account the rather distinct patterns of the time series y HCOO L and f, it seems remarkable that the stability phases agree so well over such wide parameter window. Noteworthy is the fact Fig. 4(c) and 5(c) display practically the same distribution of the periods, independently of the local number of spikes per period of the variables used to detect them. Fig. 6 compares time series of y HCOO L and f for the points i c represented in Fig. 4 and 5. This figure shows clearly that, as f 5 grows, y HCOO L develops more and more high-frequency and low-amplitude oscillations within a period. Similarly, the difference between adjacent peaks in one oscillation period of y HCOO L becomes smaller and smaller for increasing values of f 5 . As time grows, peaks grow further and further apart. At the end of this growth, a new peak emerges at the location indicated by the arrow in each panel. Thus, the sustained oscillations recorded for y HCOO L and f display the same regular dynamical evolution as a function of both f 5 and j. Of course, according to two groups of behaviors discussed in connection to Fig. 3, we expect the dynamics described in detail for y HCOO L and f to be always the same, independently of the variable used to record local peaks. Table 2 collects the coordinates for the six points i c considered in Fig. 6 as well as for the points i 7 , i 8 , and i 9 located in the black phases in Fig. 4 and 5. Also recorded in Table 2 are the period and the number of spikes per period for each individual point. Periodic points for i 8 and i 9 were obtained after discarding transients of 60 Â 10 6 and 90 Â 10 6 , respectively. As before, after suitable transients, both points located in the black phase display periodic oscillations. Fig. 7 complements our results for the galvanostatic model, describing the organization of the sustained oscillations as recorded for y HCOO B . To understand the origin of the 'finger-like' intertwined wedges corresponding to one and two, two and three, etc spikes, we computed bifurcation diagrams along the three lines in Fig. 7(a). Such diagrams display the evolution of the local maxima, spikes, of y HCOO B . Fig. 7(b) shows the bifurcation diagram obtained along the vertical line at j = 0.001, highlighting the situation for the eight arbitrary points belonging to the 1-2-1-2 spikes alternation. This bifurcation diagram shows the curious repetition of a simple pattern. As f 5 increases, the initial one-spiked oscillation eventually acquires a second spike. Upon further increase of f 5 , the former one spike disappears with what was the second spike then becoming the single spike in existence for an interval of f 5 . As this happens, the amplitude of the y HCOO B oscillations grow steadily. This cycle repeats again, with the intervals that contain one and two spikes getting narrower and narrower.
Additional bifurcation diagrams were computed along other two lines, namely for the left line running between ( j,f 5 ) = (0.0005,0.542) to (0.00045,0.43), and the right line going between (0.0005,0.542) and (0.00165,0.462). Note that such diagrams involve the simultaneous variation of two control parameters, not the common single parameter variation. From the bifurcations along the L line one sees that the above 1-2-1-2 now turns into a 2-3-2-3 repetitive pattern, occurring also for narrower and narrower parameter intervals and for increasing amplitudes of y HCOO B . Moving to higher values of j, the bifurcations along the R line show that the intervals where two spikes coexist become quite narrow, in agreement to the organization displayed in Fig. 7(a). The last two rows in Fig. 7 display the time series proving explicit details of the regular complexification of the time evolutions for the eight selected points in Fig. 7(a). Note the different vertical scales of these latter plots, showing the steady amplitude increase of y HCOO B . Fig. 8 presents more detailed diagrams illustrating the complex but regular unfolding of the 2-3-2-3 pattern. As seen, from this figure, it is possible to discern an apparently unbounded sequence of analogous m-n-m-n spike patterns, where m and n are integers differing by one unit.
So far, we have discussed only the galvanostatic model. Now we report stability diagrams obtained for the potentiostatic model. Similarly to Fig. 3, Fig. 9 compares stability diagrams illustrating the E Â R ext control parameter planes obtained by counting peaks for the five variables of the model. For the potentiostatic model, the integration step used is h = 10 À5 , five times smaller than the one used for the galvanostatic model. As it is not difficult to see, both models produce phase diagrams, which largely agree. Furthermore, as for the galvanostatic model, Fig. 9 shows that black phases of the potentiostatic model decrease systematically in size as the transient time increases. This indicates that, here also, no chaotic oscillations are present in the system. Fig. 10 illustrates the stripes observed for the potentiostatic model when counting spikes from y HCOO L while moving towards smaller values of E. The organization parallels the one found in Fig. 4 for the galvanostatic model. Remarkably, the time series obtained for the potentiostatic model look quite similar to the corresponding ones in Fig. 4. This seems to indicate that, at least in this region of control parameters, the dynamics of both models is the same. This corroborates the great similarity already transparent from the comparison of the time series presented in Fig. 2.

Conclusions and outlook
This work reported sets of high-resolution stability diagrams for parameter windows around recently studied experimental values. 6 By classifying the spikes per period 37,42-48 of the sustained periodic oscillations, stability diagrams were computed both for the galvanostatic and for the potentiostatic control modes.
The major findings distilled from our phase diagrams are briefly as follows. From Fig. 3, it is possible to recognize two types of phase diagrams, obtained either from y HCOO L , y OH , f or from y HCOO B and y CO L . The difference between these two types of phase diagrams is related with the involvement of those variables in two oscillators present in the electrooxidation of the formic acid. One of these oscillators involves the y OH species that blocks the active sites from 0.56 V up to 0.85 V approximately when y OH is consumed. The other oscillator involve the y HCOO B and y CO L species in blocking and unblocking cycles of the active sites of the platinum. This oscillator has a lower frequency than the first one.
The first type of diagrams, illustrated in Fig. 4 and 5, displays a regular series of striations, which arise from periodic oscillations acquiring more and more spikes when f 5 increases. The second type of diagrams, illustrated in Fig. 7 and 8, contains an intertwining of finger-like structures which arise from the peculiar way in which local maxima arise when f 5 is increased. For the potentiostatic model, Fig. 9 shows that the distribution of periodic oscillations closely resembles the ones obtained for the galvanostatic model (Fig. 3). Phase diagrams displaying the period of the oscillations were also computed. They reveal the existence of certain accumulation limits of the period characterized large period following a regular acquisition of a large number of small-amplitude and fast oscillations local maxima. Such large-period oscillations exist in the black regions of our phase diagrams.
A startling finding is that careful numerical experimentation seems to indicate that, in the parameter windows investigated, chaotic oscillations are nowhere to be found in the control parameter space of both models. Similar wide chaos-free phase were recently observed, 49 but only in externally driven systems. In the present case, we find chaos-free oscillations to arise in distinct setup, namely in the selforganization of the oscillations in an autonomous (non driven) system. The striations observed in the galvanostatic and potentiostatic models resemble similar one previously seen in nested arithmetic progressions of oscillatory phases in Olsen's enzyme reaction model. 50 However, in the present cases they are never mediated by chaotic oscillations. The systematic complexification of the oscillations of both models was found to occur in a rather systematic way in control parameter space. It would be now interesting to characterize the metric properties of such regular unfolding of the time series, to investigate their eventual dependence with other parameters of the models, and to compare all these predictions with experimental results.

Conflicts of interest
There are no conflicts to declare.