Nanostructured Electrodes as Random Arrays of Active Sites: Modeling and Theoretical Characterization

This review presents the main principles underlying the theoretical description of the behavior of regular and random arrays of nanometric active sites. It is further shown how they can be applied for establishing a useful semi-analytical approximation of the arrays responses under diffusion limited conditions when they involve the common situation of active sites with identical sizes. This approximation is general and, as exemplified for different type of arrays, can be employed for describing the behavior of any array involving arbitrary distributions of their active sites onto the substrate surface. Furthermore, this efficient approach allows statistical characterization of active sites distributions of any array based on chronoamperometric data.


Introduction Text
Nanostructured electrochemical interfaces (NECI) are frequently encountered in practice. For instance, partially blocked electrodes [1], defects and pinholes in self-assembled monolayers [2], pores in synthetic or natural membranes [3,4], catalytic nanoparticles (NP) or any type of catalytic site distributed at an inert or less reactive substrate [5] etc. are physical representation of NECI. In these systems, the location of the active sites cannot be controlled (excluding some specific situations) thus leading to disordered surface structures referred as random arrays of active sites. In order to model properly the response of such systems one has to deal simultaneously with multiple scales: electron transfer events at the nano-or microsized active site surface, varying distances between the disordered active sites and interaction of the diffusion layers originated at each sites, global diffusion layer propagation towards the bulk solution, eventual homogeneous kinetics etc. In addition, theoretical description of random arrays is complicated by the fact that taking into account exact positions of the active sites with respect to each other is either difficult or impossible, however, the electrochemical response of the system can be modeled employing statistical information about active sites dispersal.
One of the earliest attempts to model disordered systems in the context of partially blocked surfaces was made by one of us [1]. In this work, the behaviors of random arrays with disk-or band-like electrochemically active surface structures were investigated by relying on regular arrays with uniform size of active sites. Such treatment simplified greatly description of the system, but at the same time, it has a good predictive power as confirmed by notable citation of this early work. Nevertheless, due to computer limitations at the time it was difficult to do otherwise. The situation is different nowadays and more sophisticated approaches can be used as detailed below.

Experimental
All the results referred in this work were obtained with home-build programs (Python, Delphi) or using Mathematica 12.0 [6]. The software were run at PC with Intel Xeon E5-2650 v2 processor.

Results
We wish to take advantage of this review to summarize and unify a series of our previous works in which we considered an archetypical situation of arrays of nano-or microdisk active sites embedded in the inert surface of a substrate according to an arbitrary probability distribution. Such systems bear the most general characteristics of the systems mentioned in the Introduction. It should be noted, however, that all considerations herein are restricted to the case of active sites with monodisperse size distribution performing under diffusion limited conditions.
Modelling and simulating such disordered systems is a complex problem, however it will be shown that it can be performed on the basis of principles similar to those used for theoretical description of the regular arrays (e.g. arrays with the active sites located on a squared or hexagonal lattice).

Construction of the unit cells.
One of the main features of the regular arrays is their symmetry. Due to the symmetry the diffusion layers around each active site are identical (except those on the edge of the array, vide infra) both at short times after start of the experiment when diffusion layers of neighboring sites are separate and do not interact, but also at intermediate and long times when the diffusion layers are partially or fully overlap [1,7]. This allows delineating virtual 3D-spaces around each active site, in which all individual diffusion layers behave in an identical way. These spaces are called unit cell and their shape depends on the array lattice (Fig.1a). The limiting planes of the unit cells are the symmetry planes between corresponding neighboring sites, hence the net flux across these boundaries is zero. This significantly simplifies modeling of the regular arrays since it permits to consider each unit cell as performing independently from the rest of array and model mass transport and kinetics inside a cell without influence of its identical neighbors. Hence the total current of the regular array, , can be written as: . .
(1) where is the number of active sites in the array and is a current contribution of a single unit cell.
Two points should be noted with respect to the previous paragraph. First, strictly speaking Eq.(1) is valid only for an infinite regular array or practically for a sufficiently large array with number of active sites at the edge of the array being much smaller than the number of the sites at the interior of the array. Otherwise, in order to take properly into account contribution of the edge active sites to overall array current a correction similar to those used in [8] should be applied. Second, zero net flux boundary conditions are assigned at the unit cells boundaries making them acting independently from each other. However, the virtual boundary of the unit cells are just useful theoretical constructions and, evidently, there exist exchange of the molecules between the neighboring unit cells. Though, these exchanges are fully compensated due to the symmetry of the diffusion layers in adjacent unit cells which results in a net zero total flux across the inter-cell boundary.
Similar construct can be done in the case of a disordered array via Voronoi tessellation [9][10][11][12]. This technique is used in various branches of natural sciences and represents a natural generalization of the approach used for regular arrays. Indeed, in Voronoi tesselation symmetry planes are constructed between the closest active sites as discussed above, but due to the random location of the neighbour sites this results in a non-uniform unit cells, which are then named Voronoi cells. As can be seen in Fig.1b each segment of the oddly shaped unit cell represents the trace on the inert substrate of the vertical symmetry plane between two adjacent sites (see Fig.SM1 in Supplemental Material for construction principle illustration). Hence, noflux boundary condition can also be applied at the boundaries of the Voronoi unit cells analogously to the case of regular arrays [12].

∑ (2)
where is the current in the nth unit cell and summation is taken over all of the Voronoi cells. Each value is formally unique since each cell is distinctive in terms of its shape, surface area and active sites position with respect to the boundaries and geometrical center of the cell. It is clear that even if such representation of the random array make use of local symmetries within the array it does not completely simplifies the simulations of the system since either whole array (or a sequence of all the Voronoi cells due to their independence) should be simulated. To overcome this difficulty one may employ approximation conventionally used for simulations of the electrochemical behavior of the regular arrays.

Circular unit cell approximation.
Considering mass transport in a unit cell of the squared or hexagonal regular array requires solving a 3D mathematical model. Nowadays it is easier than before, but still it is a cumbersome, resource and time consuming process especially when a large number of independent cells need to be modelled simultaneously. Therefore each original unit cell (i.e. squared or hexagonal one) is converted into a circular unit cell as shown in Fig.2 as was performed for regular arrays [1,7]. This greatly facilitate simulations since now due to the axial symmetry of the domain the mass transport model can then be formulated in 2D cylindrical coordinates. Such formal conversion of the unit cell into a circular one evidently introduces a bias into the predicted responses of the unit cells. However, we demonstrated earlier that unless the unit cells are extremely small with respect to the active site size, the bias induced in the unit cell current due to transformation of its shape is less than few percent [13]. This validated the widely used cylindrical approximation. It should emphasized that the surface areas of the circular unit cells cross-sections must be equal to the area of the primary unit cell (see below "Chronoamperometric response of a unit cell"), which necessarily involves simultaneously non-negligible overlap and non-negligible avoided areas (Fig.2). The resulting bias can be evaluated for regular arrays (Fig.2a). This was shown to be non detrimental due to automatic compensation between overlapping and avoided solution volumes (Fig.2a, bottom). This encouraged us [12] and Compton group [10,11] to extend the circular unit cell approximation to the realm of randomly located active sites as exemplified in Fig.2b. It should be noted that the unit cells appear intersecting in the figure, nevertheless each of them acts independently due to the no-flux conditions applied at the side circular boundary of the cell as discussed above. Along with conversion of the unit cell shape the active sites were assumed to be at the center of each circular cell. This significantly reduced simulations computational cost while introducing only reasonable bias (1-7%) depending on the relative sizes of the cell and active site [12]. The latter is a good figure since experimental data not often available with accuracy better than 5% for such arrays. Note that several other approaches have been proposed for treating quantitatively electrochemical data gathered at random arrays although they are not relevant to the purpose of our former contributions [1,7,[12][13][14] and those of Compton's group [10,11]. See for example references [15][16][17].
In order to make the circular unit cell approximation fully operative, let us introduce dimensionless area of the Voronoi unit cell / / , and being geometric areas of the unit cell (being equal for original Voronoi cell and its circular approximation) and active site, and are the unit cell and site radii. Then Eq.(2) rewrites as follows: . .
where is a probability density of unit cell dimensionless area (in other words, is a normalization of values).
Equation (4) is general and can be applied to very different random and regular arrays of large dimension. Figure 3 exemplifies two different cases: uniform random array (i.e. the active sites are equiprobably dispersed over the whole surface of the array according to uniform probability distribution) and Gaussian or binormal random array (i.e. sites have tendency to be located near geometric center of the array [left bottom corner in Fig.3b] and distributed according to 2D Gaussian distribution). Figs.3a,b display partial areas of each system for clarity of presentation, while the whole images of the respective arrays can be found in Supplemental Materials. As can be seen, distributions for these two systems have similar structure, but the amplitude and shape are quite different ensuring different electrochemical responses of each of the array. It should be noted that in the case of regular arrays distribution is deltafunction, viz., since all the unit cells have the same dimensionless surface area and, hence, Eq.(4) simplifies into Eq.(1). As soon as the -distribution is known the response of an array can be evaluated via Eq.(4).
where is the number of the transferred electrons per elementary electron transfer (ET) step, is the Faraday constant, and the bulk concentration of electroactive species. At larger times after beginning of the experiment ( ≫ / ) the diffusion layer at each unit cell has fully developed (i.e. when the diffusion layers of adjacent unit cells fully overlap) thus creating a planar diffusion layer expanding into the solution [7]. Hence, a second Cottrell regime is observed with current proportional to the surface area of the unit cell: Since the same surface area of circular unit cell and the original Voronoi cell that it models must yield the same current contribution , Eq.(6) provides the relationship between their surface areas.
The presence (or not) of a third regime depends on the relative size of the unit cell with respect to the size of an active site. Indeed, if the unit cell is sufficiently large its diffusion layer acquires a hemispherical shape before reaching the limit in Eq.(6) which results in a quasisteady state current [7]. Thus for large unit cells ( larger than few units) a perfectly defined transition takes place between the two above Cottrell regimes involving a quasi-steady state current, 4 . For smaller unit cells ( 1 3) this regime cannot be fully observed resulting in a smooth transition between the two limiting currents.
These limiting regimes could be hardly identified at a conventional amperogram (Fig.4a). However, if the same amperogram is plotted in a log-log scale it reveals the two or three modes (Fig.4b) depending on value as just discussed. This representation of the current, which appears as a tilted sigmoid, suggests useful normalization for the problem at hand. Indeed, the current of a unit cell normalized with respect to its (time dependent) Cottrell limits: has a shape of a sigmoid (Fig.4c). It should be emphasized that hereafter this expression is used in the context of a circular unit cell for the sake of simplifying the presentation, however, the above analysis and normalization remain valid for original Voronoi cells with the same outcome [12,14]. (7) has two immediate advantages. First, it can be extended thanks to Eq.(4) to provide the current of the whole random array in the normalized form (see Supplemental Material for the derivation of this expression) [12]:

Response of a random array. Normalization
where / is the dimensionless time, and are the corresponding Cottrell limits of the array. Second, due to the convenient shape of the dimensionless currents (7) they can be easily parametrized via rather simple analytical fitting equation [12]: ( 1 1 ) Note that the lower integration limit in Eqs. (4) and (11) is unity since by construction 1. Eq.(8) proves to be useful also for a more practical and complicated case of reconstruction the distribution of unit cell sizes on the basis of the measured current. Indeed analyzing experimentally is a tedious and difficult exercise, while for some systems such characterization is even impossible. However, theoretically this can be achieved as follows based on the analysis of chronoamperometric currents. First, let us discretize and represent it as a histogram, , with beans: where is a bin function, , Δ and ℎ are center, width and height of a k-th bin. Vector of coefficients ℎ ℎ , … , ℎ represents thereby a distribution of the unit cell sizes of a random array as soon as the binning parameters ( , Δ ) are set. Minimizing the difference between the evaluated and experimental currents upon varying the coefficients ℎ allows reconstruction of the thought distribution. The latter statement can be written formally as a constrained minimization problem: is the normalized experimental current. The constraints (14)-(15) on coefficients ℎ follow from the properties of probability densities, more precisely that it is a non-negative function whose integral over its definition range is equal to unity.
For the sake of demonstration of described reconstructed procedure one may assume that the currents shown in Fig.5 are the experimental ones and the distribution in Fig.3c,d are not known. The corresponding -distributions can then be attempted to be identified by solving optimization problem (13)- (15). The outcome of this is shown in Fig.6 revealing a very good agreement for both cases between the primary (solid lines) and reconstructed (histograms) distributions.

Discussion
As can be seen from Eq.(8) the distribution of unit cell sizes is the only characteristics responsible for unique features of electrochemical response of a given array as soon as the active sites activities and size are monodisperse. The variations of the normalized current in and 32.0, respectively. Conversely, the dispersion of -distribution affects the shape of the Φ . Regular arrays give rise to the fastest transition due to their high symmetry (black curve in Fig.5). Indeed, all their unit cells have the same sizes so that the transition between the different diffusion regimes occurs simultaneously for all of them. On the contrary, for the Gaussian array featuring a long tail in its -distribution the transition is sluggish since unit cells corresponding to a range of larger perform their transitions to the second Cottrell regime at increasingly larger times. It should be noted that transition of the uniform random array starts earlier even if it is slower than the one of the regular array. This, in fact, reflects the unit cells with the smaller surface areas than those of the regular array (see Fig.6a), while the larger unit cells are responsible for sluggish behavior for Φ 0.4 (Fig.5). This sensitivity of Φ to the shape of -distribution make possible the reconstructions shown in Fig.6. Moreover, these results together with those obtained in our previous work [14] demonstrate that with this approach one may characterize, in a statistical sense, the distribution of the active sites and judge if the array under the scrutiny is regular or random as well as deduce information on the type of probability distribution governing the scattering of active sites over the inert substrate. This allows identifying possible biases introduced into the active sites distribution by fabrication procedure or other factors. Additional advantage of this framework that all this information can be extracted from chronoamperograms and does not require any sophisticated equipment, vacuum chambers etc. It should be noted, however, that the reconstruction procedure described above relies on the normalization of the measured current (8), which suppose a good experimental resolution of each Cottrelian limiting regimes (5)-(6). This would be a rare case since the latter implies a good time resolution over several orders of time magnitude as can be seen from the time range in Fig.5. Nevertheless, as we demonstrated the reconstruction performs adequately well when realistic range of times are taken into account with proper approximations performed to renormalize the experimental current [14].
For the identification of the probability distribution governing locations of the active sites one may rely on evaluation of normalized currents with various and then compare with the results of the reconstruction. In this respect, a good reference for this kind of comparison, consists in the empirical analytical distribution established and validated for random arrays with active sites distributed according to a uniform distribution as shown in Fig.3a,c (i.e. without any specific bias) [18]. For our case it writes as: where Γ is the Gamma function. It should be noted, however, that Eq.(16) is valid only in the limit of infinite number of sizeless active sites. In practice this corresponds to sufficiently large numbers of sites whose dimensions are significantly small with respect to the size of the system. This was the case for the array shown in Fig.3a (see Supplementing Material for more details). The solid line shown in Fig.3c was evaluated according to Eq. (16) and resulted in excellent agreement with the histogram of unit cell sizes obtained from this particular uniform random array.
In terms of perspective, our approach can be employed, in principle, for dynamic characterization of nanostructured electrodes with randomly distributed catalytic sites. For example, it is known that during their operation such systems often show some poisoning or clustering of the catalytic sites etc. Theoretical approach reviewed in this paper performed on a series of amperograms acquired sequentially during the system operation may evidence dynamical changes of catalytic sites with respect to their statistical distribution, average cluster size etc.
Finally, it should be stressed that the above described approach is valid under the assumptions indicated at the beginning of 'Results' section, that is all active sites have almost identical dimensions and perform in diffusion limiting regime (i.e. when the concentrations is constant across electrode-solution interface). In case of presence of significant ohmic drop, preceding or following to ET step chemical reactions (of the order higher than one) as well as slow charge transfer the concentrations at the electrode surface will be a function of position as well as the locations of the neighboring electrodes, thus resulting in different time-and space-competition between adjacent Voronoi cells. In simple terms, the pattern and sizes of the Voronoi cells defined at the electrode surface level will not be identical to those prevailing at different distances in the solution. Therefore, the pattern of Voronoi cells used in our works reviewed here and in those of Compton group will not be anymore a geometric characteristics of a given array, but will depend on the experimental time or scan rate for voltammetry as discussed in [12]. We are currently developing theoretical strategies, which seem to offer interesting solution to the problem. These will be reported in a future work after their validation, in particular in the context of voltammetry of fast and slow electron transfer couples.

Conclusions
The approach reviewed in this work allows an efficient modeling and prediction of the electrochemical behavior of nanostructured electrodes of identical sizes dispersed onto an inert substrate. It represents a unified framework to simulate electrochemical responses of regular and random arrays under diffusion limited conditions as soon as the distribution of unit cells dimensionless surface areas is provided. Conversely, the latter can be reconstructed from chronoamperograms giving access to an extremely important information that is excessively difficult to obtain otherwise. In particular, the degree of active sites distribution randomness can be easily evaluated so that different types of random arrays can be differentiated.