The modelling of carbon-based supercapacitors: Distributions of time constants and Pascal Equivalent Circuits

Supercapacitors are an emerging technology with applications in pulse power, motive power, and energy storage. However, their carbon electrodes show a variety of non-ideal behaviours that have so far eluded explanation. These include Voltage Decay after charging, Voltage Rebound after discharging, and Dispersed Kinetics at long times. In the present work, we establish that a vertical ladder network of RC components can reproduce all these puzzling phenomena. Both software and hardware realizations of the network are described. 
 
In general, porous carbon electrodes contain random distributions of resistance R and capacitance C, with a wider spread of log R values than log C values. To understand what this implies, a simplified model is developed in which log R is treated as a Gaussian random variable while log C is treated as a constant. From this model, a new family of equivalent circuits is developed in which the continuous distribution of log R values is replaced by a discrete set of log R values drawn from a geometric series. We call these Pascal Equivalent Circuits. Their behaviour is shown to resemble closely that of real supercapacitors. The results confirm that distributions of RC time constants dominate the behaviour of real supercapacitors.


Introduction
Recently, a Universal Equivalent Circuit was proposed for carbonbased supercapacitors [1]. For a single carbon electrode (or two matched carbon electrodes in one cell) the circuit consists of a vertical ladder network of RC components, with each "rung" of the ladder having its own characteristic time constant, see Fig. 1. In the present work, we develop software and hardware realizations of this model, and then compare their responses with real supercapacitors.
Supercapacitors differ from ideal capacitors in a number of important ways. For example, they show voltage decay at open circuit, capacitance loss at high frequency, and voltammetric distortions at high scan rate. Why do these phenomena occur? The short answer is that the complex network of pores inside activated carbon behaves like a vast collection of micro-capacitors wired in parallel, with each wire having its own resistance. As a result, supercapacitors exhibit asynchronous charging. That is, each microcapacitor charges at a different rate. The existence of this phenomenon explains why "mean field" models of carbon electrodes (which assume that all pores charge at the same rate) fail to capture the full transient dynamics of real systems.
Although there are many ways of parameterizing the ladder network shown in Fig. 1A, the most natural choice is to use the admittance y, because admittances in parallel are additive. We therefore adopt that approach here. To begin, we focus on the n th rung of the ladder network considered in isolation. In that case, the parameters R n and C n are constants and the individual rung admittance (siemens) is given by the equation y n ðuÞ ¼ juC n 1 þ juR n C n (1) where the ranges of R n and C n are assumed to be strictly positive.
Evidently, the time constant of each rung is R n C n ¼ t n .
Now we consider all the rungs in parallel. The total admittance of the ladder Y L is obtained by performing the summation Y L ðuÞ ¼ X NL n¼1 y n (2) where N L (dimensionless) is the total number of rungs in the ladder, and R n and C n are discrete random variables. For a few hundred rungs having known values of R n and C n this sum can be computed term-by-term. However, for millions of rungs having broad distributions of R n and C n the evaluation of Y L is more challenging.
To make progress, we assume that the population of rungs in the ladder is so dense that the distributions of R n and C n are effectively continuous. We also assume that the distributions of R n and C n are statistically independent. In that case we can define a joint probability density function qðR; CÞ having a single normalization This then allows us to write Eq. (2) in the compact form where the ranges of R and C are again assumed to be strictly positive.
Ideally, we should like to measure qðR; CÞ directly. However, in the absence of experimental techniques for doing this, we here develop some useful approximations of Eq. (4) valid in cases where the range of either R or C is severely restricted. The first model that we consider (denoted Model I) assumes that the micro-capacitors are all different while the resistors are all the same. The second model (denoted Model II) assumes that the resistors are all different while the capacitors are all the same.

Derivation of Model I
If C is a random variable while R is a constant, then the total admittance Y L of the ladder network as a function of frequency u takes the much simpler form: Here R is the resistance of one rung in the ladder, N L is the total number of rungs in the ladder, and gðtÞ is the normalized probability density function of t, i.e. the spectrum of RC time constants in the system: In the case of a lumped equivalent circuit, the function gðtÞ will have a finite number of discrete "spectral lines" corresponding to the finite number of RC time constants. In the case of a real electrode, however, the function gðtÞ will be essentially continuous.

Derivation of Model II
If R is a random variable while C is a constant, then the total admittance Y L of the ladder network as a function of frequency u takes the form: Or equivalently Here C is the capacitance of one rung in the ladder, N L is the total number of rungs in the ladder, and f ðtÞ is the normalized probability density function of t, i.e. the spectrum of RC time constants in the system: Model I was developed more than thirty years ago by Brug et al. to quantify inhomogeneous surfaces on uniformly accessible electrodes [2]. Model II, however, is new, and is used here to quantify homogeneous surfaces on non-uniformly accessible electrodes (porous carbon).
Inside porous carbon, the pore radii typically have values over the range rz0:1 nm to rz10:0 mm, which spans five orders of magnitude. This implies that pore resistances may have values that span ten orders of magnitude, since (for cylindrical pores) the pore resistance is given by the formula Here r is the pore resistivity, L is the pore length, and r is the pore radius. By contrast, surface capacitances of carbon typically have values over the range C ¼ 2 mF cm À2 to C ¼ 70 mF cm À2 , depending on whether graphitic basal planes or oxidized edge planes are exposed [3]. Such a range of capacitance values spans less than two orders of magnitude. Overall, therefore, Model II appears much closer to reality. An idealized diagram of the joint distribution of log R and log C values inside an activated carbon electrode is shown in Fig. 1B.

Experimental and results
In order to explore some of the properties of Model II, we now introduce a new family of equivalent circuits that can mimic the wide distributions of resistances found in carbon supercapacitors. We call them "Pascal Equivalent Circuits". Unlike conventional equivalent circuits, Pascal Equivalent Circuits are not "mean field" approximations, and therefore they are not restricted to a small number of RC time constants. Instead, they have 2 nÀ1 RC elements arranged in a vertical ladder network, with a single value of C, but multiple values of R drawn from a geometric series. The number of resistors of a particular value of R is determined by the coefficients of the n th row of Pascal's Triangle ( Fig. 2A).
To illustrate the idea, a "Row 5" Pascal Equivalent Circuit is shown in Fig. 2B. This is a 16-rung vertical ladder network having five sets of RC values arranged in sequence of their time constants. Each set of RC values contains a population of "rungs" determined by the binomial coefficients: 1, 4, 6, 4, 1. Consistent with Model II, all the micro-capacitors are the same, but the sets of resistors differ by factors of ten. As a result, the Pascal Equivalent Circuit displays an approximately lognormal distribution of RC time constants (i.e. a normal distribution of log t values), and thus resembles the very broad distribution of RC time constants found in real supercapacitors.
The theoretical link between Pascal's triangle and the lognormal distribution is provided by the binomial distribution. In the notation of Fig. 2A, the binomial distribution takes the form PrðX ¼ kÞ ¼ Cðn; kÞp k ð1 À pÞ nÀ1Àk (12) where Cðn; kÞ ¼ ðn À 1Þ! k!ðn À 1 À kÞ! (13) and PrðX ¼ kÞ is the probability that the random variable X takes value k. Thus, for the special case p ¼ 1=2, the k coefficients along the n th row of Pascal's triangle are identical to the k coefficients of the binomial distribution. (This may readily be verified by substitution.) Further, by a well-known result in statistics, the binomial distribution approaches the normal distribution for large n; this trend is illustrated in Fig. 2C. [Explicitly, if the distribution of the random variable X is binomial Bðn; pÞ then for sufficiently large n the distribution approaches the normal form Nðm; sÞ with m ¼ np and s 2 ¼ npð1 À pÞ.] Finally, the normal distribution may be transformed into a lognormal distribution by replacing successive values of k with successive values of a geometric series YðkÞ. For example, when constructing trial models in the present work, we found it convenient to use the decadal values Y ¼ 10 k where k ¼ 1; 2; 3…n. The probability density function of a lognormally-distributed parameter Y is given by the equation where m and s are the mean and standard deviation of the underlying normal distribution. A plot of lðYÞ versus Y is shown in Fig. 2D.
Clearly the lognormal density function is strongly skewed to the right, where it is responsible for the very large time constants that appear in practical devices. Interestingly, the skewness is a sensitive function of the underlying standard deviation s,

Voltage decay at open circuit
One feature of real supercapacitors that Pascal Equivalent Circuits reproduce outstandingly well is the phenomenon of Voltage Decay at open circuit. An example is shown in Fig. 3A. This illustrates the effect of charging the "Row 5" Pascal Equivalent Circuit shown in Fig. 2B for different periods of time at 2:0 V. Upon switching to open circuit, the fast-responding micro-capacitors (i.e. those that have completed their charging process at 2:0 V) immediately begin to discharge into the slow-responding micro-capacitors (i.e. those that have not completed their charging process at 2:0 V), thus causing the Voltage Decay.
From a health and safety perspective, it is reassuring to note that the Voltage Decay in real supercapacitors is not due to uncontrolled electrochemical reactions inside the carbon electrodes. Instead, it is due to the wide distribution of RC time constants in the pores, and can therefore be avoided simply by charging devices for longer times. (Even the most resistive pores will become charged eventually, even if it takes 24 h).

Voltage rebound
Traditionally, it has been difficult to simulate large electrical circuits such as the one illustrated in Fig. 2B. However, equivalent circuit modelling has recently become much easier with the release of Simscape™, a high-level programming language developed by MathWorks, Inc. (Natick, Massachusetts, USA) [4]. Simscape™ not only allows researchers to model electrical equivalent circuits of arbitrary size, but also allows the models to be run on a digital computer. Conveniently, Simscape™ operates in the Simulink ® programming environment, and so is fully compatible with MATLAB ® .
During the present work, Simscape™ was used to simulate the rapid discharge process of supercapacitors, with the goal of equalizing their terminal voltages and thus rendering them safe. (By definition, a fully discharged supercapacitor at open circuit has a voltage of zero volts.) Some model shut-down data are shown in Fig. 3B. To obtain this plot, a Simscape™ model of the "Row 5" Pascal Equivalent Circuit was charged at 1 V for 10 s, and then discharged at 0 V for 1 s. What happened next was very revealing d a Voltage Rebound occurred after the system was switched to open circuit.
Analysis of the transient responses of individual "rungs" provided the explanation. On open circuit, the slowest RC states (i.e. those which had remained full during the 1 s spent at 0 V) were now re-charging the fastest RC states (i.e. those which had emptied during the 1 s spent at 0 V). The spontaneous occurrence of this effect in the Simscape™ model was compelling evidence that it was also occurring in reality. The only remaining question was how broad the distribution of time constants might actually be in a real system.
In theory, a voltage rebound can be generated with as few as two "rungs" in the ladder. However, such a parsimonious allocation of  time constants is inadequate to model the asymptotic long-time behaviour of real electrodes (see below). In the real world, it turns out that a very broad distribution of time constants is needed.

In operando modelling
The successful implementation of Pascal Equivalent Circuits in silico (i.e. as virtual models in a computer) prompted us to explore the possibility of implementing them in operando (i.e. as actual physical models in the real world.) Before embarking on this task, however, it was first necessary to ensure that any commercial resistors and capacitors that we used would respond as closely as possible to their mathematically idealized counterparts. It turned out that near-ideal resistors were readily available from many suppliers, but that near-ideal capacitors were not. Indeed, many offthe-shelf capacitors showed voltage decay on the time scale of seconds, which made them useless for our purposes. Some examples of the voltage decay of commercial capacitors are shown in Fig. 4A. Fortunately, the aluminium/polyethylene-terephthalate (PET) capacitors retained their voltage over the timescale of many seconds, so these were used in all our physical models.
In passing, we remark that the physical origin of voltage decay in ordinary capacitors may be different to that in supercapacitors. For example, in ordinary capacitors, voltage decay may be due to current leakage through the dielectric.
By using aluminium/PET capacitors we were able to implement a fully-functioning hardware model of a "Row 5" Pascal Equivalent Circuit on a printed circuit board (Fig. 4B). For ease of comparison, the component values were selected to be the same as those in Fig. 2B. To test the physical device, we recorded plots of its open circuit voltage versus time, after charging potentiostatically at 2 V for various time intervals (1e10 s). The results are shown in Fig. 5A. The similarity of these results with those of the Simscape™ simulation (Fig. 3A) is noteworthy. In both cases, the voltage decay curves are exponential with time constants equal to the largest time constant in the ladder (2.2 s). No slower states are observed. It is interesting to compare the simulated voltage decay curves with real-world data from a commercial supercapacitor (Fig. 5B). The latter device shows the expected voltage decay, but the voltage-time curves are slightly non-exponential. This is due to the presence of a very long "tail" of ultra-slow states inside the carbon electrodes. The existence of this "tail" means that real supercapacitors must be held at zero volts for 24 h between individual measurements, in order to ensure that their discharge is complete.
In our laboratories, we have found that hardware models of supercapacitors can be very useful for both pedagogical and practical purposes. However, although it is comparatively easy to match the RC products of hardware models with those of real supercapacitors, it is sometimes difficult to match the individual values of R and C. This is because real supercapacitors may have exceptionally high values of C. For example, a real supercapacitor may have a capacitance >10 F, much higher than any normal capacitor.   (10, 20, 30, 40, 50, 60 s). Terminals short-circuited for 24 h between measurements.

Non-exponential voltage decay curves
The Voltage Decay curves of real supercapacitors are similar to those of the hardware and software models, except the curve shapes tend to be non-exponential rather than exponential. This subtlety arises from the fact that the RC time constants in real supercapacitors are arbitrarily close together, whereas those in the models are a finite distance apart. As a result, real supercapacitors tend to display multi-exponential decays on all time scales. Some typical experimental data are shown in Fig. 5B.
As noted in Section 2.3, powerful evidence for the existence of a distribution of RC time constants in real supercapacitors is the appearance of a Voltage Rebound after a period at short circuit, coupled with non-exponential voltage-time responses. An example is shown in Fig. 6A. In this case, an NEC/Tokin activated carbon/ sulfuric acid supercapacitor (nominally 10 mF) was charged at 1 V for 2 s then short-circuited for 2 s. A rising non-exponential voltage-time curve then appears. This indicates that a redistribution of electrical charge is taking place from the slower RC states to the faster RC states inside the carbon electrodes.
In order to curve-fit the voltage-time data we tested two wellknown models of dispersed kinetics, namely the stretched exponential model and the lognormal distribution model. The stretched exponential model is obtained by inserting a constant b into the exponential function: Typically, the choice of the constant b is restricted to the range ð0; 1Þ so that a graph of log f versus t has a "stretched" appearance. Among physicists, the stretched exponential function is also known as the Kohlrausch function, after Rudolf Kohlrausch, who used it to model the residual voltage on a Leyden jar [5].
On the other hand, the lognormal distribution model is obtained by inserting the exponential function into a Fredholm integral equation of the first kind: Here m and s are the mean and standard deviation of the random variable k, and 4ðkÞ is its probability density function [6,7]. In other words, the distribution of the argument k is assumed to be normal. Each of the above models added just one degree of freedom to the classic case of exponential decay. In the case of the stretched exponential function the extra degree of freedom was the exponent b, while in the case of the lognormal distribution the extra degree of freedom was the standard deviation s. However, in both cases the improvement in the goodness-of-fit was dramatic. Indeed, both distributions could be fitted to the experimental data within experimental error at asymptotically long times, although the lognormal distribution fared slightly better at short times (Fig. 6B). The obvious conclusion is that several theoretical models can be fitted to the experimental data within practical accuracy. That is, the fitting problem is "ill-posed" in the sense of Tikhonov [8].
The experimental observation of smooth, non-exponential, voltage-time plots in real supercapacitors over more than three orders of magnitude in time is powerful evidence in favour of a broad distribution of time constants. If a distribution of time constants is arbitrarily truncated, the models show that the asymptotic behaviour of the voltage at long times is dominated by the largest time constant remaining. Some examples are shown in Fig. 7A and B. These figures provide plots of open circuit voltage versus time for the Simscape™ simulation and for the hardware model, after charging potentiostatically for fixed periods of time. The curves are all exponential functions with time constants of 2.2 s (the slowest in the ladder). No such exponential responses are observed in the real world. Instead, the non-exponential responses indicate the existence of a highly skewed distribution like that in Fig. 2D.

Discussion
The vertical ladder network shown in Fig. 1A reproduces many of the characteristic features of carbon-based supercapacitors. However, the exact distribution of RC time constants remains an unsolved problem. Clearly, it would be very desirable to measure this distribution with accuracy and precision. At present, however, all we can say is that the distribution is broad and highly skewed to the right.
One method of estimating the distribution of RC time constants is by means of molecular dynamics simulations. As recently shown by Pean et al. [9] these techniques generate molecular level properties that can then be inserted into equivalent circuit models. Similarly, Kondrat and Kornyshev [10] have used Monte Carlo simulations to probe the remarkably complex phenomena that occur inside nanopores. In particular, they have shown how ions and solvent molecules compete for space [11]. In future, it would be instructive to combine these molecular approaches with equivalent circuit theory.

Conclusions
In the present work, two probability models of supercapacitor electrodes have been developed based on the assumption that the internal resistance and internal capacitance are random variables. The first model that we considered (denoted Model I) assumed that the capacitance values were widely different, while the resistance values were all the same. The second model that we considered (denoted Model II) assumed that the resistance values were widely different, while the capacitance values were all the same. It was then noted that, in activated carbon electrodes, there was a much wider spread of log R values than log C values, suggesting that Model II was more appropriate.
Model II was then transformed into a new kind of equivalent circuit, which we called a "Pascal Equivalent Circuit". Unlike conventional equivalent circuits, the Pascal Equivalent Circuit was not restricted to a small number of RC time constants, but instead had 2 nÀ1 RC elements arranged in a vertical ladder network. The number n could be arbitrarily large. By implementing Pascal Equivalent Circuits in both software and hardware several counterintuitive responses of supercapacitors were explained. These included Voltage Decay after charging, Voltage Rebound after discharging, and Dispersed Kinetics at long times. All of these were previously mysterious.
Finally, it is interesting to compare the ladder model of Fig. 1A with the literature standard, which is the de Levie transmission line model for a single pore [12e17]. Evidently, the models are mutually indistinguishable (degenerate) if the de Levie model is generalized to many different pores in parallel. The fact that neither model is definitive therefore serves as a warning against over-interpretation of individual resistance and capacitance values derived from one model alone.