Reproducing capacitive cyclic voltammetric curves by simulation: When are simplified geometries appropriate? Electrochemistry

The usage of multi-physics simulation tools is steadily increasing in the field of electrochemistry. While this is a great opportunity for closing the gap between analytical electrochemists used to simple 1D models and exper-imentalists, there are possible pitfalls that must be avoided. In this work, we raise awareness on numerical ar- tifacts that can mislead the interpretation of cyclic voltammetry experiments through simulations of geometries with different number of spatial dimensions. In particular, we show that one-dimensional simulations can suffer from substantial errors when models go beyond charge neutrality assumption. We exemplify such situations using simple electrolyte/electrode structures with 1D, 2D and 3D geometries. We then show the occurrence of artifacts related to the geometry of the simulation domain on the simulation of cyclic voltammetric curves as those typically performed to characterize conjugated polymer/electrolyte blends. All the models are imple- mented using COMSOL Multiphysics and are accompanied by a detailed description of their implementation. However, geometrical artifacts identified in this work also apply to other simulation approaches.


Introduction
In the field of electrochemistry, several methodologies employing analytical calculus have been successfully employed to predict electric currents at chemically active electrodes, surface reactions kinetics, electrostatic potentials, diffusion mechanisms, etc. [1,2]. These approaches employ idealized geometrical structures that admit analytical closed-form solutions and often cannot represent real systems. The consequence is that one generally resorts to numerical solutions for accurately reproducing complex geometries [3,4]. Among these solutions, finite element simulation (FEM) platforms are gaining popularity, thanks to the increasing computational performance of modern computers, friendly user interfaces, the possibility to handle different physics coupled together (i.e., multiphysics frameworks) and embed them in multidimensional and arbitrary geometries. The applications are manifold. For example, numerical simulations might be useful to validate analytical models that are derived based on approximations whose correctness is not trivial to assess [4][5][6] or to reproduce experimental data while extracting model parameters that are not easy to find directly neither from the experiments nor from analytical models [7]. On the contrary, the main drawback of FEMs is that simulations can be time demanding and may require high computational resources when one has to fit experimental results over a large number of unknown parameters. This is particularly relevant when dealing with electrochemical systems in more than one physical dimension and with simulations that consider multiple physics at once (e.g., chemical reactions, ion transport, fluid flow, heat dissipation, etc.) leading to highly nonlinear models. As a result, time-dependent multiphysics simulations of micrometer-sized structures can easily take up to hours.
One way to greatly shorten the simulation time and model complexity is to reduce the number of spatial dimensions from 3 to 2, or even to 1 by exploiting planar/axial symmetries and neglecting edge effects. Unless symmetries can be exploited mathematically to restore the three-dimensional nature of the problem (e.g., modeling the real geometry as hemispherical or axis-symmetric structure) [8,9], these simplifications may lead to numerical artifacts when diffusion and migration problems are considered. The reason is that reduced geometries do not properly describe electric current dispersion lines, this phenomenon being intrinsically three-dimensional. In electrochemistry, one way to overcome this problem is using the electroneutrality assumption when dealing with concentrated supporting electrolytes [10] where the drift transport of ions is suppressed [11][12][13][14]. Refined models, however, do require that drift and diffusion effects are considered under varying applied voltages, such as in the Poisson-Nernst-Planck modeling framework. In these cases, the use of geometries that oversimplify the real geometry, such as 2D [15][16][17][18] or even 1D [16,19,20], have a considerable risk to be affected by geometrical artifacts.
In this article, we show that the use of simplified geometries in electrochemical systems where diffusion and migration of charges take place under varying applied voltage can lead to wrong results, unless some precautions are considered. To show this, we consider two common situations as case studies. We first analyze an electrolyte contacted by a reference electrode and an ideally polarizable electrode, and show that the same physical model can predict very different results if applied to 3D, 2D or 1D geometries. Simulation results are then explained by means of simplified equivalent circuit analysis. In the second example, we scout geometry-related artifacts in a complex published model [20] to simulate cyclic voltammetry curves of conjugated polymers and polyelectrolyte blends.
The main result of our analysis is that one-dimensional simulations can be critical as the current paths through the bulk of the electrolyte are not correctly taken into account. According to our predictions, these errors are much less critical in 2D geometrical domains thus being a trade-off option to avoid significant errors (such as those in 1D) when (i) it is impossible to identify symmetries and, (ii) when the 3D simulation comes at the expenses of significantly increased simulation time with little additional insights [17,18].
The article is structured as follows. In Section 2 we introduce the models and the simulation platform while in Section 3 we show that the geometrical issues can be studied by means of equivalent circuits, in the case of simple contact/electrolytes systems. Simulation results and discussions are reported in Section 4 which includes a few concluding remarks.

Poisson-Nernst-Planck model
The Poisson-Nernst-Planck (PNP) model is used to simulate ion transport in diluted electrolytes. The electrostatics is governed by the Poisson equation [1]. For monovalent ions only, it reads: where ε el is the dielectric constant of the electrolyte, ψ is the electrostatic potential, F is Faraday's constant and the space charge density is written in terms of sum of the (monovalent) cation and anion concentrations, c + and c − respectively, expressed in molar (M) units. Time enters the continuity equations that provide single ion mass conservation in the space-time domain. Considering ionic fluxes, f, the continuity equation reads: In turn, the flux is expressed using the Nernst-Planck formalism where, neglecting convective effects for simplicity, one has: where D ± is the space-dependent diffusion coefficient of the ionic species, R is the gas constant and T is the temperature of the system, assumed constant. In order to account for steric effects and avoid unrealistic ionic concentration at the charged boundaries of blocking electrodes, we use a Stern layer, namely, an ideal dielectric layer between the electrolyte and the WE with a capacitance C SL = 0.20 F/m 2 [21]. This translate into an additional boundary equation that imposes where σ SL is the surface charge at the electrode and V App is the applied potential.
Eqs. (1)-(4) model an electrolyte in the absence of surface or distributed chemical reactions. The boundary conditions of the Poisson equation allow to fix the potentials at the reference electrode (RE) and at the working electrode (WE). Further details on the implementation in COMSOL Multiphysics [22] are provided in the Supplementary Information S1.

Model of conductive polymers for cyclic voltammetric simulations
In cyclic voltammetry (CV), cyclic potential ramps are commonly used to characterize electrochemical cells [1]. In the case of conductive polymers (CPs), subject of many research studies in the last decades [16,[23][24][25][26], CVs reflect the ion-exchange between the conductive polymer/polyelectrolyte blends and the surrounding electrolyte bath during the injection or the subtraction of electronic carriers from the WE [16,20,27]. Experimental evidence reveal that blends of conjugated polymers and polyelectrolytes show an intertwining of polymer strings permeated by water and ions at the microscopic level. However, despite this complicated three-dimensional entangled structure, the coupling between the ionic and electronic phase at the macroscopic level is purely capacitive [16,27] and has been successfully reproduced by PNPinspired models introducing the concept of volumetric capacitance, as proposed in [20] for PEDOT:PSS blends (see Section 2.2).
According to the model in [20], hereafter denoted PNP-CP model, ions in the electrolyte and holes in the CP are governed by two modified sets of PNP equations, where the Poisson equations for the different phases are coupled. For the ionic carriers (i.e., ions), the Poisson equation takes into account the electronic charge (e.g., holes, p, in PEDOT: PSS) distributed in the conductive phase: where c fix is fixed charge distribution, given by the negatively charged PSS − groups. Differently, the electrostatic potential of electronic carriers is affected by the presence of ions and fixed charges in the form of capacitive coupling with the electrostatic potential of ions: In Eq. 6, the coupling factor is the phenomenological volumetric capacitance parameter, C V , that quantifies the 3D electrical double layer between the two immiscible phases in a straightforward manner.
In the PNP-CP model, continuity equations remain the same as those seen in the conventional PNP (see Section 2.1), with the only difference being that the diffusion coefficient of holes depends on their local concentration instead of assuming a fixed value as for ions in each phase [20]. The current at the WE was computed from FEM simulation results in COMSOL by time-differentiation of the volume-integrated hole concentration in the CP: Further details on the boundary conditions and the model implementation in COMSOL Multiphysics can be found in the Supplementary Information S3 and in [20].

A simple equivalent circuit model
To interpret PNP FEM simulations carried out with different spatial L.J. Mele et al.
dimensions, we employ a simple equivalent electrical circuit ( Fig. 1) that describes an electrolyte solution contacted by two electrodes. A more thorough description of this system and related equations can be found in [28]. In the Laplace domain, the voltage at the node V D in Fig. 1, i.e., the potential taken right outside the diffuse layer, at the edge of the quasi-neutral region where charge screening effects are negligible, reads where C DL is the electrical double layer (EDL) capacitance (i.e., the series of the Stern Layer, C SL , and the diffuse layer capacitance, C dl , calculated as is the bulk resistance of the electrolyte, V App is the electrical stimulus applied at the WE and τ = R B C DL is the time constant. C SL is assumed constant in the analytical model used to interpret the FEM simulation, while C dl is calculated numerically by simply integrating the charge in the diffuse layer (as obtained from the PNP simulation) and dividing by the voltage drop between the Stern Layer and the node V D . In this way, we extend the validity of the equivalent circuit to an applied potential much larger than the thermal voltage [28]. Assuming a 1D geometry for the sake of a simpler notation, the diffuse layer capacitance per unit area reads: where T SL is the Stern Layer thickness and T D the screening length. By considering a ramp stimulus in Eq. 8, V App (s) = v/s 2 (where v is the scan rate in V/s), and using the final value theorem for t→∞, one obtains 2 : We observe a direct proportionality between the steady-state potential V D and the bulk electrolyte resistance, which is known to change according to the physical dimensions and its geometry. Table 1 reports the analytical values of R B calculated for idealized structures, such as a line (1D), a circular (2D), and a spherical (3D) crown representative of three different scenarios of approximating more complex geometries. In the first two cases, the value of R B diverges to infinity for increased electrolyte volumes whereas only in the spherical crown it converges to a finite value regardless of the size of the simulated domain. Since R B enters the expression of V D (see Eq. 10), the equivalent circuit analysis suggests that an increased electrolyte volume has different effects in 1D, 2D and 3D geometries, which clearly underlines the possible presence of numerical artifacts when approximating complex structures with simplified 1D or 2D geometries. For example, if one simplifies the structure under study with a 1D geometry, the results would change with T el . In this respect, a 3D domain, even solved in spherical coordinates (that has the same complexity of a 1D simulation) does not suffer from this limitation. These aspects will be analyzed in detail employing both the analytical model and numerical simulations in Section 4.

The impact of geometry on PNP simulations
The analytical theory outlined in the previous section is now The one on the left is the WE and is assumed ideally polarizable whereas the grounded one at the right is the faradaic RE. The potential nodes V SL and V D indicate the electrostatic potential at the Stern Layer thickness, T SL and the distance right beyond the double layer, T D , respectively. T el is the total electrolyte thickness. Note that the horizontal axis in the bottom is not to scale.

Table 1
Sketch of different electrolyte structures in three different spatial dimensions that permit analytical calculus of the liquid resistance. The electrolyte thickness is defined as T el = b − a, A is the cross-section area of the 1D template structure and the resistivity is ρ = RT/(F 2 (D c + + D c − )c 0 ). Note that the three geometries are representations of different template structures.

Spatial dimensions Sketch
2 Notice that when a voltage step is applied instead of a ramp, V D tends to zero following an exponential decay. compared to full PNP numerical simulations taking into account the geometry dependence of R B . The PNP model applied to a simple electrolyte/contact domain is sufficient to highlight the expected geometrical artifacts on the simulation results. For ease of understanding and for better reference with the analytical solutions proposed in Section 3.1, we simulate the structures illustrated in Table 1, where a indicates the size of the WE while b − a is the electrolyte thickness, T el . Unless otherwise specified, the whole external surface acts as RE. In the following numerical simulations (parameters values reported in Table S2), we apply a potential ramp at the WE, using the RE as the ground node. We then record the electrostatic potential 50 nm away from the WE and at t = 1 s, that is much larger than τ. This corresponds to V D in Fig. 1 at the steady-state. The calculated V D are reported in Fig. 2.a vs the thickness of the electrolyte for the different structures shown in Fig. 2.b. Symbols indicate the PNP numerical solution obtained with COMSOL whereas the dashed lines are the calculations using Eq. 10. Despite its simplicity, the equivalent circuit can accurately reproduce the simulation results over a large range of T el values. We observe that V D features a noticeable dependence on the T el in 1D and 2D geometries with respect to the almost flat 3D case. In fact, these variations raise up to almost four orders of magnitude between 1D and 3D for T el ≈ 5 mm, with trends that follow the scaling predictions of R B reported in Table 1. The concerns on geometrical artifacts are therefore worthy of the attention in one-dimensional structures where the electrolyte resistance and so V D (see Eq. 10) increase linearly with the electrolyte volume. Differently, the voltage drop in the electrolyte increases logaritmically in 2D with a dumped increase of the electrolyte resistance with T el . At this point, one could argue that a RE covering the whole external surface is not realistic. To address this concern we simulated the 3D structure with the RE localized at the top of the outer sphere (with 1 μm radius) and restricted the Dirichlet ionic boundary conditions to the same region and so the current from the WE. Results (green triangles in Fig. 2.a) show that the profile is still quite flat as in the case where the RE covers the whole external surface but with increased voltage drop due to higher electrolyte resistance. Interestingly, as shown in Fig. S2 in the Supplementary Information, this voltage drop remains relatively high in the bulk and only goes to zero in close proximity of the RE.
Finally, we highlight that units of millivolts are used in the y-axis of Fig. 2.a. Despite this is tightly related to the choice of the scan rate (see Eq. 10), with typical values of v = 1 V/s, after 1 s the errors due to voltage drop in the bulk of the electrolyte reach up to few millivolts in 1D simulations with a T el of few millimeters. As we show in the following section, these apparently negligible errors are in fact sufficient to cause large variations in the electrochemical cell responses of, e.g., cyclic voltammetric experiments.

Cyclic voltammetry and the "geometrical pitfall"
In this section, we show that CV curves of conductive polymers are remarkably affected by geometrical artifacts in the same way as explained in Section 4.1 for simple electrode/electrolyte systems. We consider CP/electrolyte structures as a vehicle to show that distortions stemming from oversimplified geometries under time-variant applied potential profiles arise in one-dimensional simulations. However, similar considerations also apply to more complex cases.
As an example, we consider the CV measurements of a PEDOT:PSSelectrolyte blend proposed in [20]. In that reference, the authors used a 1 cm 2 gold WE coated with a 600 nm thin CP. Reasonably, given the extremely high aspect ratio of the electrode geometry, the authors adopted a one-dimensional geometry in their simulations, that is equivalent to assuming that the polymer disk and the reference electrode above are the bases of a cylinder filled with the electrolyte. In this work, we implement the same model equations as in [20] and we extend them to different spatial dimensions, namely, we adopt the 1D, 2D, and 3D geometries 3 depicted in Fig. 3.d as case studies, while varying the electrolyte thickness of orders of magnitude higher and lower than the 7.9 mm used in the reference study. The result for the 1D case is shown Fig. 3.a, where the lines with symbols represent our simulations and where T el spans over three orders of magnitude. For the sake of comparison, circles and the dashed line in black report the experimental and simulation results from [20]. We observe that, despite in all cases the same steady-state capacitive current is achieved, the rise time of the I − V App curve increases for larger electrolyte domains. Figure S4 in the Supplementary Information clearly shows that when approximating the CP-Electrolyte structure as an RC series equivalent circuit, a direct proportionality exists between the time constant and the size of the electrolyte. Such an increase in the rise time does not converge for higher values of T el (Fig. 3.a). In reality, one would rather expect that larger bulk electrolyte volumes do not affect the WE nor the interface with the conductive polymer. As shown in Fig. 3.b such expectations are satisfied if the simulation domain is extended to a 2D geometry (which can be seen as a planar electrolyte sheet sandwiched between two insulating layers) where the third physical dimension is not taken into account. Nevertheless, the results are not different from the more accurate 3D case in Fig. 3.c. This is in accordance with Fig. 2.a where the potential drop in the electrolyte for a two-dimensional structure is about three orders of magnitude lower than in 1D and, does not change  Table S2. In Eq. 10, ρ ≈ 8 Ωm.
appreciably for larger electrolyte domains. In fact, in the three simulated structures (Fig. 3.d), the calculated bulk electrolyte resistance R B (that depends on the chosen T el value) lies in the range of 0.05 − 100 Ω for the 1D case, 4 − 14 mΩ for the 2D and ≈ 2 mΩ for the 3D one. We remind the reader that, except for the 3D case in Fig. 3.d, the 1D and 2D geometrical are oversimplifications of the actual 3D system. However, these examples recall choices of approximated geometries found in the literature [15][16][17][18][19][20].
In [20] the authors found the best match with the experiments using the thickness of the electrolyte among the fitting parameters (black dashed line in Fig. 3, where T el = 7.9 mm). In the following, we provide another perspective where this parameter is not so critical. The voltage drop in the electrolyte layer, that in [20] contributed to the best-fit, can instead derive from a (more realistic) voltage drop due to the nonidealities at the contact electrodes, (e.g., charge transfer resistance at the RE or contact resistance at the small WE). This is equivalent to adding a resistance in series to the simulated system as explained in Section S5. Fig. 3.e, compares, i) the results of the best-fit proposed by the 1D model in [20], ii) our results using the same 1D model but setting T el = 10 μm and adding a series resistance of R S = 42 Ω that is consistent with the R B values for millimeters size electrolytes considered previously (blue line with circles) and, iii) the two-dimensional model version of (ii) (red line with diamonds). The perfect match between i) and ii) reveals that 1D structures with a series resistance can reproduce the same fit as in [20] without the need of including a very large electrolyte layer. Then, the good match between ii) and iii) shows that using small electrolyte thickness does not introduce significant voltage drops between 1D and 2D simulations. In fact, using small electrolyte domains mitigates the build up of unrealistic resistance in the bulk of the electrolyte, thus yielding results in 1D which are similar with respect to the (less affected) 2D structure, ultimately almost identical to a threedimensional domain.
When dealing with geometrical artifacts, the example proposed in this section conveys the message that with both awareness of the electrolyte resistance dependence on space and the support of numerical simulations, scientists can determine the extent of the error induced by these artifacts and even choose to adopt 1D or 2D structures over complicated 3D ones when errors are negligible. Knowledgeable decisions such as this one allow one to save computational time and reduce complexity as can be seen in Section S6 of the Supplementary Material, where a benchmark of the simulations run for this work is reported.  Table S4. (a) 1D, (b), 2D and, (c), 3D model results (relative to the geometries in d) are compared with the 1D results in [20] for different choices of the electrolyte length. (d) Simulated structures, that make use of cylindrical coordinates for consistency with the experimental setup used in [20]. In COM-SOL, the 2D and 3D models are solved using twodimensional differential operators for simplicity, since in the 3D structure the axial symmetry cancels out the partial derivatives with respect to the azimuthal angle φ (see Fig. 2.b). (e) The best fits are replicated in 1D and 2D using a relatively limited electrolyte thickness but adding a series resistance at the RE electrode. Note that in (e), the blue line with circles and the black dashed line are almost perfectly overlapped.

Conclusions
In this paper, we have shown that the use of an approximated geometry can result in numerical artifacts when electrolyte models include the Poisson-Nernst-Planck formalism without the assumption of charge neutrality at the bulk of the electrolyte. For the case of 1D simulations, an appreciable potential drop can arise in large domains since the geometry imposes constraints on the current density lines. The consequence is that a non-zero electric field can develop also in regions of the device that are supposed to be equipotential. Despite these voltage drops appear relatively small (even when considering one-dimensional systems), their impact on ionic concentration profiles is remarkable as it is shown for cyclic voltammetric simulations of conductive polymers.
In the authors' opinion, despite the greatly simplified simulation experience that multiphysics simulation platforms such as COMSOL provide to users, these tools should be handled with care, as other works also suggest [5,6]. Obviously, the best modeling recommendation to avoid geometrical artifacts is to replicate actual 3D system geometries. When symmetries can be identified, the computational burden can otherwise be reduced by using lower dimensional structures and bearing that equations should reflect the existence of a rotational axis. When this is not possible, an approximated 2D version can also be satisfactorily used as the distortions may still be minimal. Finally, when 1D is the only option, we recommend modelers using electrolyte domains that do not go far beyond the charge screening length, as layers larger than tens of micrometers result in incorrect and non-negligible series resistances.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.