Editors' Choice—Review—Nanostructured Electrodes as Random Arrays of Active Sites: Modeling and Theoretical Characterization

Alexander Oleinick; Oleksii Sliusarenko; Irina Svir; Christian Amatore

doi:10.1149/2.0302001JES

Nanostructured electrochemical interfaces (NECI) are frequently encountered in practice. For instance, partially blocked electrodes,¹ defects and pinholes in self-assembled monolayers,² 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⁵ 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 site, global diffusion layer propagation toward 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.¹ 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.⁶ 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 and positioned according to an arbitrary probability distribution. Such assemblies 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, i_array, can be written as:

where N_site is the number of active sites in the array and i_cell is a current contribution of a single unit cell.

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** Unit cells for active sites located on a square lattice (a) and top view of the Voronoi unit cells (b).

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 in 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⁸ 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–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, no-flux boundary condition can also be applied at the boundaries of the Voronoi unit cells analogously to the case of regular arrays.¹²

In this situation Eq. 1 rewrites as:

where i⁽ⁿ⁾_cell(t) is the current in the nth unit cell and summation is taken over all of the Voronoi cells. Each i⁽ⁿ⁾_cell(t) 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 to solve 3D problems, 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. 2a 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.¹³ 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. 2a). The resulting bias can be evaluated for regular arrays and was shown to be non detrimental due to automatic compensation between overlapping and avoided solution volumes¹³ (Fig. 2a, bottom).

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** Schematic illustration of circular unit cell approximation for (a) regular array and (b) for random array. In each panel the bottom figure represents the cylindrical unit cells equivalents (a) or their cross-section with the inert substrate plane (b).

This encouraged us^12,14 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.¹² 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^15–17 although they are not relevant to the purpose of our former contributions^1,7,12–14 and those of Compton's group.^10,11

In order to make the circular unit cell approximation fully operative, let us introduce dimensionless area of the Voronoi unit cell ρ = A_cell/A_site = (r_cell/r₀)², A_cell and A_site being geometric areas of the unit cell (being equal for original Voronoi cell and its circular approximation) and active site, r_cell and r₀ are the unit cell and site radii, respectively. Then Eq. 2 rewrites as follows:

where summation is taken over all possible unit cells surface areas, i^(ρ)_cell being the current in each circular unit cell of dimensionless area ρ and N^(ρ)_site the number of such cells. For a large array ρ may be assumed to vary continuously (and not having just discrete values as in Eq. 3). Thus, Eq. 3 can be rewritten in integral form as follows

where f(ρ) is a probability density of unit cell dimensionless area ρ (in other words, f(ρ) is a histogram of normalized N^(ρ)_site 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, 3b 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 (Figs. 3c,3d), f(ρ) 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 f(ρ) distribution is delta-function, viz., f(ρ) = δ(ρ − ρ₀) since all the unit cells have the same dimensionless surface area ρ₀ and, hence, Eq. 4 simplifies into Eq. 1.

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** Statistically generated uniform (a) and Gaussian (b) random arrays (see text for the definitions) along with their Voronoi tesselations and their respective unit cell size distributions f(ρ) (c, d). The average dimensionless surface areas are ρ^uf_avg = 23.2 (uniform) and ρ^g_avg = 32.0 (Gaussian), correspondingly. Solid lines in c) and d) are the fits of the histograms. Red line in c) evaluated via Eq. 16 (see below).

As soon as the ρ-distribution is known the response of an array can be evaluated via Eq. 4. Additionally, since the unit cells have the same uniform circular structure except for their ρ values the contributions i^(ρ)_cell(t) can be evaluated and stored in advance or represented via analytical approximation further simplifying computations. To emphasize this important point, let us examine the chronoamperometric behavior of one of such unit cells.

Chronoamperometric response of a unit cell

Considering diffusion limited conditions, as was indicated above, all unit cells experience two limiting Cottrell regimes irrespective of their sizes, although the time ranges in which this occurs depends on the unit cell dimensions. First, at short times (i.e. when t ≪ r²₀/D, D is the diffusion coefficient of electroactive species) any unit cell gives rise to a Cottrell current governed by active site surface area πr²₀, i.e.:

where n is the number of the transferred electrons per elementary electron transfer (ET) step, F is the Faraday constant, and c_bulk the bulk concentration of electroactive species.

At larger times after beginning of the experiment (t ≫ r²_cell/D) 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.⁷ 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 i^∞_cell(t), 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 quasi-steady state current.⁷ 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, i^ss_cell = 4nFDc_bulkr₀. For smaller unit cells (ρ ≈ 1 − 3) this regime cannot be fully observed resulting in a smooth transition between the two limiting Cottrell 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

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** Amperometric responses of unit cells with ρ = *A_cell*/*A_site* = 9, 25, 49 and 100 plotted in a) linear scale; b) log-log scale; c) normalized current via Eq. 7.

Response of a random array

Normalization 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):¹²

where τ = Dt/r²₀ is the dimensionless time, i⁰_array(τ) and i^∞_array(τ) 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:¹²

where a = 9 × 10^{− 7} and parameters b^ρ₀ and b^ρ₁ are functions of unit cell surface area:

as obtained through fitting of computed φ^ρ_cell(τ) responses for a large set of ρ values. The approximation 9–10 is valid for ρ ⩽ 103, which is amply sufficient for our purpose (compare Figs. 3c, 3d). Evidently, this can be extended to encompass a wider range of ρ if necessary via the same principle.

Knowing the contributions φ^ρ_cell(τ) in analytical form allows a straightforward evaluation of random array responses from Eq. 8 as soon as their distributions of unit cell sizes f(ρ) is available. The average value ρ_avg does not need to be supplied as an independent entry since it can be readily computed from the f(ρ) distribution according to the definition of a mean value:

Note that the lower integration limit in Eqs. 4 and 11 is unity since by construction ρ ≥ 1.

Figure 5 displays the variations of the normalized currents with the dimensionless time evaluated via Eqs. 8–10 for the two types of random arrays shown in Fig. 3. The average normalized surface areas for uniform and Gaussian arrays were ρ^uf_avg = 23.2 and ρ^g_avg = 32.0, correspondingly. In addition, the response of a regular array with ρ^reg_avg = 25.0 is provided for comparison.

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** Normalized current Φ(t) for uniform (red, ρ^uf_avg = 23.2) and Gaussian (blue, ρ^G_avg = 32.0) random arrays shown in Fig. 3 along with response of regular array (black, ρ^reg_avg = 25.0).

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 f(ρ) 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 f(ρ) and represent it as a histogram, H_M(ρ), with M beans:

where B_k(ρ) is a bin function, ρ_k, Δρ and h_k are center, width and height of a k-th bin. Vector of coefficients h = (h₁, ..., h_M) represents thereby a distribution of the unit cell sizes of a random array as soon as the binning parameters (ρ_k, Δρ) are set. Minimizing the difference between the evaluated and experimental currents upon varying the coefficients h_k allows reconstruction of the thought distribution. The latter statement can be written formally as a constrained minimization problem:

where Φ_exp is the normalized experimental current. The constraints 14−15 on coefficients h_k 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 reconstruction procedure one may assume that the currents shown in Fig. 5 are the experimental ones and the distribution in Figs. 3c, 3d 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.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** f(ρ) distributions reconstructed (histograms) from the amperometric responses of (a) uniform and (b) Gaussian random arrays shown in Fig. 5. The corresponding arrays are shown in Figs. 3a, 3b and their original distributions (solid lines) are the same as shown in Figs. 3c, 3d.

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 Fig. 5 evaluated for different kind of arrays confirm this point. Indeed, the range displayed by the Φ(τ) sigmoid along the time axis is largely defined by the average ρ value. The larger average unit cell size the later is transition toward the second Cottrell regime, as evidenced by the shift to larger times of normalized current for Gaussian array (blue curve, Fig. 5) relatively to uniform one (red curve, Fig. 5) in agreement with their relative mean ρ values ρ^uf_avg = 23.2 and ρ^G_avg = 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¹⁴ 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.¹⁴

For the identification of the probability distribution governing locations of the active sites one may rely on evaluation of normalized currents with various f(ρ) 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 probability distribution as shown in Figs. 3a, 3c (i.e. without any specific bias).¹⁸ For our case it writes as:

where Γ(x) 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.¹² 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 f(ρ) 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.

Acknowledgments

This work was supported by CNRS, Ecole Normale Superieure – PSL Research University, Sorbonne University (UMR 8640 PASTEUR). CA acknowledges Xiamen University for his position of Distinguished Professor.

ORCID

Alexander Oleinick 0000-0001-6194-0459

Oleksii Sliusarenko 0000-0003-1667-3552

Irina Svir 0000-0002-6006-3148

Christian Amatore 0000-0002-2561-1137

Editors' Choice—Review—Nanostructured Electrodes as Random Arrays of Active Sites: Modeling and Theoretical Characterization

Author notes

Author notes

Notes

Article metrics

Share this article

Dates

Abstract