As the world moves more toward unpredictable renewable energy sources, better energy storage devices are required. Supercapacitors are a promising technology to meet the demand for short-term, high-power energy storage. Clearly, understanding their charging and discharging behaviors is essential to improving the technology. Molecular Dynamics (MD) simulations provide microscopic insights into the complex interplay between the dynamics of the ions in the electrolyte and the evolution of the charge distributions on the electrodes. Traditional MD simulations of (dis)charging supercapacitors impose a pre-determined evolving voltage difference between the electrodes, using the Constant Potential Method (CPM). Here, we present an alternative method that explicitly simulates the charge flow to and from the electrodes. For a disconnected capacitor, i.e., an open circuit, the charges are allowed to redistribute within each electrode while the sum charges on both electrodes remain constant. We demonstrate, for a model capacitor containing an aqueous salt solution, that this method recovers the charge–potential curve of CPM simulations. The equilibrium voltage fluctuations are related to the differential capacitance. We next simulate a closed circuit by introducing equations of motion for the sum charges, by explicitly accounting for the external circuit element(s). Charging and discharging of the model supercapacitor via a resistance proceed by double exponential processes, supplementing the usual time scale set by the electrolyte dynamics with a novel time scale set by the external circuit. Finally, we propose a simple equivalent circuit that reproduces the main characteristics of this supercapacitor.

## I. INTRODUCTION

Supercapacitors have the unique ability to store much more energy than conventional capacitors and to charge and discharge much faster than batteries. Supercapacitors also have other important benefits over batteries, such as a higher efficiency and a longer lifetime.^{1} However, the energy density of supercapacitors is one order of magnitude lower than that of batteries. As renewable, clean energies are replacing fossil fuels, the use of long-term and short-term energy storages is expected to rapidly increase in the near future. Supercapacitors are already deployed in many applications^{2,3} and could potentially play an important role in short-term high-power energy storages. This will only happen if their energy density can be improved, e.g., by combining optimized porous electrodes with ionic liquid electrolytes.^{1} However, the application of ionic liquids comes at the cost of a lower power density due to the slow ion transport in ionic liquids, especially in the electrodes’ pores.^{4}

Molecular dynamics simulations are a powerful tool to investigate the sluggish dynamics of ionic liquids in electrode pores since they provide detailed information at the atomic level.^{5–10} These simulations require dedicated routines to model the redistribution of the charges on the electrodes in response to the dynamics of the ions: the Constant Potential Method (CPM) endows the electrode atoms with narrow Gaussian charge distributions of evolving strengths;^{11–16} the Induced Charge Computation (ICC*) method treats the electrodes as media of infinite dielectric permittivity with evolving charge distributions at their interfaces;^{17} and in the Lorentz–Drude model, a light-weight negative charge representing the fluctuating electron cloud is tethered to a positively charged nucleus.^{18} We refer the reader to recent reviews on simulations of electrodes for more details on these methods and the insights gained by simulations.^{9,10,19} In previous studies, the charging and discharging of supercapacitors were simulated by subjecting the electrodes to a potential difference with a pre-determined evolution in time.^{20–28} In practice, however, one does not control the potential difference between the electrodes during discharging. Rather, in the absence of another power source in a discharging circuit, it is the capacitor itself that determines the potential difference. This difference causes a current to flow in the electric circuit external to the capacitor, which in turn reduces the charges on the electrodes and, thus, co-determines the decay of the potential difference. Similarly, in a simple charging circuit with the capacitor connected to a potential source via a resistance, the capacitor co-determines the voltage difference between its electrodes. As a first step to simulating the (dis)charging process, one would like to simulate a supercapacitor with fixed total charges on its electrodes, i.e., a disconnected supercapacitor or open-circuit. Dufils *et al.*^{29,30} modeled an open circuit by combining CPM with an appropriately tuned electric displacement, which they applied to a single electrode surrounded on both sides by the same periodically continued electrolyte; this setup results in opposing conserved charges on the two surfaces of the electrode. Two recent studies on a “galvanostatic mode” by Zeng *et al.*^{31} and the “ConQ ensemble” by Tee and Searles^{32} simulate two electrodes with conserved total charges, while the charge distribution within each electrode is free to evolve in response to the electrolyte dynamics. These methods result in a fluctuating voltage difference between the electrodes, with a standard deviation that relates to the differential capacitance.^{32} By slowly varying the total charges on the electrodes, one gains access to the supercapacitor's (dis)charging dynamics.^{30,31}

Here, we introduce a more general Constant Sum Charges Method (CSCM) that works with any standard electrostatics solver, whereas the aforementioned methods require dedicated solvers.^{30–33} We next introduce an equation of motion for the sum charges, with the voltage difference between the electrodes (co-)determining the current in a closed circuit. This has a major impact on the (dis)charging dynamics of a supercapacitor: whereas previous simulation studies report relaxation time scales set by the dynamics of the electrolyte^{20,22,24,25,30,33} we find an additional relaxation process due to the external circuit. The ratio between internal and external time scales determines the energy efficiency of the discharging process.

From a macroscopic point of view, complex electronic components—such as supercapacitors—are more conveniently modeled by equivalent circuits, i.e., combinations of elementary ideal electronic components that reproduce the electronic behavior of the complex components.^{34} Equivalent circuit models for supercapacitors typically model the energy storage in the electrode–electrolyte interface as a capacitor and the electrolyte’s viscous resistance to ion transport as a resistance.^{35} For porous electrodes, this model is extended to the transmission line model by approximating the electrode as a stack of slabs, with each slab combining a capacitor for storage in double layers, a resistance for the viscous flow of the electrolyte, and a conductor for the current in the solid fraction.^{20,28} We derive an alternative model in which the volume between two flat electrodes also acts as a capacitor, in parallel to the resistance due to the electrolyte, to account for electrostatic interactions between the two double layers. The analytic solution of the model predicts a double-exponential decay and explains how the relaxation dynamics varies with the value of the external resistance, which is in good agreement with the CSCM simulations.

This paper is organized as follows: The theory underlying the simulation method is presented in Sec. II, along with a phenomenological equivalent circuit. The simulation model, an aqueous electrolyte between flat graphene-like electrodes, is introduced in Sec. III. Simulation results on disconnected capacitors, as well as the evolution of charging and discharging capacitors, are presented in Sec. IV and compared with the equivalent circuit. We summarize our main conclusions in Sec. V.

## II. THEORY

In this section, we derive a microscopic description of a supercapacitor subject to charge conservation on both electrodes separately, derive expressions for the macroscopic properties of this supercapacitor, and propose a phenomenological equivalent circuit model describing the charging and discharging processes.

### A. Constant sum-charge method (CSCM)

^{36,37}atoms are endowed with effective charges located at the atoms’ positions, representing the total charge of the atom due to the positive charge of the nearly point-like nucleus plus the negative charge of the surrounding electron cloud. By the laws of electrostatics, the exact charge distribution within an atom is of minor importance at the interatomic distances relevant to MD, and one conveniently approximates atoms as point charges. We model the electrodes as

*N*atoms with fixed positions

*R*_{i}and variable charges

*Q*

_{i}, while the electrolyte contains

*n*atoms with fixed charges

*q*

_{j}and variable positions

*r*_{j}. The electrostatic energy of this system reads as

**and**

*Q***, respectively. The elements of the symmetric**

*q**N*×

*N*matrix

**A**, the non-symmetric

*N*×

*n*matrix

**B**, and the symmetric

*n*×

*n*matrix

**C**vary with the positions of the atoms. Routine electrostatic solvers determine the energy and forces without evaluating these matrices, using Ewald summation or the Particle–Particle Particle–Mesh (PPPM) method,

^{38–40}but the explicit evaluation of

**A**is crucial to constant potential methods. It is common practice in the CPM field to treat the electrolyte atoms as point charges, while the charges of the electrode atoms are assumed to be Gaussian distributed. The net effect of this narrow point-spread is to endow point-like electrode atoms with chemical hardness,

^{41–43}i.e., an intra-atomic energy associated with non-zero charges,

*K*

_{p}denotes the strength of this energetic “penalty” on non-zero charges. In the above equation, it has been assumed that all electrode atoms are of the same type, which renders

*K*

_{p}a common constant and also eliminates an intra-atomic energy term linear in the charges accounting for atomic eagerness to attract or repel electrons, which is known in a chemical context as electronegativity.

^{41,44,45}The total charge-related energy,

*U*

_{Q}=

*U*

_{e}+

*U*

_{p}, can again be expressed in the above vector-matrix form by introducing the matrix

**A**

_{p}=

**A**+

*K*

_{p}

**1**. Because the electron distribution evolves on a much shorter time scale than the motion of the ions, one may assume that the electron clouds are continuously in equilibrium with the momentaneous ion distribution, i.e., the Born–Oppenheimer approximation. The vector collecting the electrical potentials of the atomic charges on the electrodes—potentials, for short—is obtained by differentiating the energy with respect to the charge distribution,

*V*

_{R}(

*V*

_{L}). The vector collecting the potentials of all electrode atoms, therefore, takes the form

*Y*_{R}(

*Y*_{L}) are unity for atoms on the right (left) electrode and zero for atoms on the opposite electrode. The resulting charges on the electrode atoms are then solved by inverting Eq. (3) and inserting Eq. (4), yielding

*V*

_{R}= (

*V*

_{∗}+ Δ

*V*/2) and

*V*

_{L}= (

*V*

_{∗}− Δ

*V*/2), respectively, where the potential difference Δ

*V*=

*V*

_{R}−

*V*

_{L}is selected by the user and

*V*

_{∗}= (

*V*

_{R}+

*V*

_{L})/2 denotes the average potential difference between the electrodes and the environment. In early implementations of CPM,

^{12,13}this potential difference was tacitly set to zero, while in later implementations,

^{15,46}the potential difference is solved by zeroing the total charge of the simulation box, which for a neutral electrolyte is realized by eliminating the total charge on the electrodes,

*∑*

_{i}

*Q*

_{i}= (

*Y*_{R}+

*Y*_{L}) ·

**= 0. Due to the perpetual motions of the ions, the charges on the electrode atoms fluctuate with time. Even the total charge per electrode fluctuates in time, around a value dictated by the imposed potential difference.**

*Q**V*

_{R}(

*t*) and

*V*

_{L}(

*t*).

^{31,32}Combining Eqs. (4) through (6) gives the charges on the two electrodes as

**D**read as $D\alpha \beta =Y\alpha \u22c5Ap\u22121Y\beta $, with

*α*and

*β*taking the values “

*R*” and “

*L*”. Inversion yields the potentials

*V*

_{R}and

*V*

_{L}for the given electrode charges

*Q*

_{R}and

*Q*

_{L}at the momentaneous ion distribution. The charges of the electrode atoms are then solved in CSCM from

^{32}

^{,}

*Q*

_{R}and

*Q*

_{L}with the vectors

**B**, from altering the sum charges on the electrodes since

*q*^{29–32}

*j*th ion follows by differentiating the potential energy with respect to the ion’s position,

**is extended with a force arising from the position-dependent induced charges of the electrode atoms. The latter term can be rewritten as**

*Q**T*, the resulting model samples the canonical Boltzmann distribution at constant electrode charges,

**denoting the collective coordinates of the electrolyte atoms,**

*r**Z*denoting the normalizing configuration integral,

*U*

_{Q}denoting the charge-related energy evaluated at

**,**

*Q**U*

_{≠Q}denoting the charge-independent interactions such as the van der Waals terms and covalent contributions, and

*β*= 1/(

*k*

_{B}

*T*) with Boltzmann’s constant

*k*

_{B}.

^{32}

*g*

_{R}and

*g*

_{L}and the vector

**are constants and, hence, need to be computed only once.**

*h**R*

_{e}, say from a positive right electrode to a negative left electrode. Ohm’s law gives the current, i.e., the time derivative of the sum charges on the electrodes, as

*Q*

_{R}(

*t*) and

*Q*

_{L}(

*t*) after every simulation step based on the momentaneous voltage difference between the electrodes, Δ

*V*(

*t*), as calculated by Eq. (15). Charging is simulated by connecting the supercapacitor in series with an ideal voltage source

*V*

_{e}and a resistance

*R*

_{e}. The current flowing in this circuit and, hence, the time derivatives of the sum charges then follow from the potential over the resistance as

*I*= (Δ

*V*−

*V*

_{e})/

*R*

_{e}.

### B. Capacitance

*C*=

*Q*/Δ

*V*of a supercapacitor will, in general, vary with the potential or charge. Its differential capacitance

*C*

_{diff}=

*dQ*/

*d*Δ

*V*can be deduced from the tangent to a plot of the imposed charge

*Q*vs the average resulting potential difference at that charge, ⟨Δ

*V*⟩

_{Q}.

^{43,46–48}

*Y*_{±}= (

*Y*_{R}−

*Y*_{L})/2 and the symmetric matrix

^{43,46}

**=**

*Y*

*Y*_{R}+

*Y*_{L}. Note that replacing the vector

*Y*_{±}in Eq. (21) with its projected counterpart,

^{46}

*Y*_{±}−

**(**

*Y***·**

*Y*

*Y*_{±})/(

**·**

*Y***), does not affect $Cdiff,0CPM$ because the matrix**

*Y***S**

_{p}removes homogeneous distributions,

**S**

_{p}

**=**

*Y*

*Y*^{T}

**S**

_{p}= 0.

^{32}by focusing on the derivative of the average potential difference,

*U*

_{Q},

*G*_{α},

*g*

_{α}, and

**, and where**

*h*

*G*_{±}=

*G*_{R}−

*G*_{L}. Differentiating the electrode-induced average potentials on the electrode atoms then gives

**·**

*h*

*G*_{±}

*Q*is a constant. Insertion in Eq. (23) and noticing that

**·**

*h***B**represents the only fluctuating contribution to the potential difference Δ

*q**V*in Eq. (15), it follows that

^{32}

^{32}

*V*,

*Q*_{CPM}=

**S**

_{p}(Δ

*V*

*Y*_{±}−

**B**), matches that in CSCM for electrode sum charges

*q**Q*

_{α}=

*Y*_{α}·

*Q*_{CPM}. Similarly, the charge distribution for CSCM for electrode charges ±

*Q*matches that in CPM for a potential difference derived from these sum charges by Eq. (15). This equivalence for a given electrolyte configuration notwithstanding, the ensemble sampled at constant Δ

*V*differs from the ensemble sampled at constant electrode charges ±

*Q*.

### C. Equivalent circuit

*A*with uniform surface charges ±

*Q*/

*A*separated by a height

*H*, two counter-ionic monolayers with uniform surface charges ∓

*ξ*/

*A*at a distance

*d*from the electrodes, and a uniform ion distribution in the central region of width

*w*=

*H*− 2

*d*; see Fig. 1. Note that we do not equate

*ξ*with

*Q*, as these two charges will differ when the system is out of equilibrium and may even deviate for a system in equilibrium. Upon infinitely extending these four monolayers along the two directions parallel to the electrodes, the model becomes one-dimensional. The electric field due to the idealized homogeneous left ionic layer is then readily evaluated as $E\xi ,L=\xb1\xi /(2\epsilon 0\epsilon rA)e\u0302z$, where the plus (minus) sign applies to positions to the right (left) of this ionic layer,

*ɛ*

_{0}is the vacuum permittivity,

*ɛ*

_{r}is the relative permittivity of the background solvent—whose molecules are assumed to be much smaller than the ions—and where the unit vector $e\u0302z$ points from left to right. Similar expressions describe the fields by all four idealized surface charge densities. By linear combinations of these four expressions, the field in the region between an electrode and the adjacent ionic layer is evaluated as $EQ\xi =\u2212Q/(\epsilon 0\epsilon rA)e\u0302z$, where the fields due to the two oppositely charged ionic layers cancel out as both lie to the same side of the region of interest. Similarly, in the region between the two ionic layers, the four fields combine into $E\xi \xi =(\xi \u2212Q)/(\epsilon 0\epsilon rA)e\u0302z$. The total potential energy stored in these electric fields reads as

*ξ*. We, therefore, model the free energy of the capacitor as

*k*will be treated as fit parameters and, in general, will depend on the details of the ions and the electrodes; their extraction from the simulations is discussed below. For a non-connected charged capacitor in thermal equilibrium, the ionic charge accumulated at the electrodes is solved from $\u2202F/\u2202\xi =0$ as

*C*

_{eq}, is independent of the charge,

*C*

_{QQ}= 1/

*k*

_{QQ}in series with two parallel capacitors

*C*

_{ξξ}= 1/

*k*

_{ξξ}and

*C*

_{Qξ}= 1/

*k*

_{Qξ}; a symmetrized version hereof is shown in Fig. 2. The above relation between the charges also applies when the supercapacitor is slowly (dis)charged at a constant temperature. Interesting behavior is expected when the (dis)charging time scale is decreased and becomes comparable to the evolution time scale of the ion distribution.

*F*

_{Q}= −Δ

*V*. If the two electrodes are connected by an external resistance

*R*

_{e}, as shown in Fig. 2, and assuming we are in the overdamped regime, the generalized velocities of the charges are linear in the forces,

*M*

_{Q}= 1/

*R*

_{e}and the ion mobility

*M*

_{ξ}relates to the friction experienced by the ions as they move in the electrolyte. The equivalence of these equations of motion to the circuit in Fig. 2 is proven in Appendix A. The analytical solution of this discharge model gives a double-exponential decay of the charges

*Q*and

*ξ*,

*A*and relaxation times

*τ*

_{±}derived in Appendix B. It then follows from Eq. (33a) that the potential difference Δ

*V*of the supercapacitor also decays double-exponentially, with the same relaxation times. In the limit of a large external resistance, the ions are continuously in quasi-equilibrium with the electrode charges [see Eq. (30)], and the relaxation reduces to a single exponential decay,

*R*

_{e}= 0, the electrode charge quickly drops to a quasi-equilibrium with the ions, solved from $\u2202F/\u2202Q=0$ as

*ξ*, and hence

*Q*, relaxes by a single exponential process,

*R*

_{ξ}= 1/

*M*

_{ξ}and an effective internal capacitance

*C*

_{int}. This equation reflects that for

*R*

_{e}= 0, the equivalent circuit reduces to a resistor

*R*

_{ξ}in series with both a capacitor

*C*

_{ξξ}and the parallel pair of capacitors

*C*

_{QQ}and

*C*

_{Qξ}. This relaxation process is evidently not affected by the external resistance; this decay law also applies to CPM simulations where the potential difference between the electrodes is instantaneously reduced to zero. In general, however, the dynamics of the current external to the capacitor will affect the relaxation time scales within the capacitor.

*V*

_{e}to the force

*F*

_{Q}in the equations of motion. The resulting potential difference

*V*

_{e}− Δ

*V*describes the force driving the current through the external resistor, and multiplication by the charge mobility gives this current as

*V*=

*V*

_{e}; see Appendix B.

## III. SIMULATION SETUP

^{49}using a slight modification of the CPM routines by Wang

*et al.*

^{13}to calculate the

**A**matrix for point-charge electrode atoms;

^{43}the CPM and CSCM codes are available on GitHub.

^{50,51}We consider an aqueous electrolyte with monovalent ions between two fixed parallel graphene-like layers with a separation of

*H*= 4 nm. Each layer consists of 960 atoms, covering an area of

*A*= 25.15 nm

^{2}at a bond length of 0.142 nm. The simulation box is periodically repeated in the two Cartesian directions parallel to the electrode surfaces; in the third direction, perpendicular to the slab, periodic copies of the box alternate with empty spaces of height 2

*H*. The electrolyte consists of 51 or 156 pairs of monovalent ions. Excluded volume interactions between all atoms are described by the Weeks–Chandler–Andersen (WCA) pair potential,

*ɛ*=

*k*

_{B}

*T*at temperature

*T*= 298 K and diameters

*σ*= 0.5 nm. As explained elsewhere,

^{43}the potentials on the electrode atoms due to the ions,

**B**, are readily calculated using standard electrostatics routines. The charge-related forces on the ions [see Eq. (12)] can also be calculated by standard electrostatics routines. We here use the Particle–Particle Particle–Mesh (PPPM) method,

*q*^{39}with a cutoff of 1.2 nm in real space and a relative tolerance of 10

^{−6}for the long-range forces, including a correction for the slab geometry.

^{52}In the absence of an explicit solvent, the screening of Coulombic interactions by an aqueous electrolyte must be included in the model to prevent the ions from crystallizing. A common expedient is to introduce a water-like relative permittivity of

*ɛ*

_{r}= 78. The resulting scaling of Coulombic energies and forces is appropriate for ions separated by many water molecules, while at shorter distances the interaction will depend on the actual configuration of the water molecules; the variable charges on the electrode atoms very effectively screen Coulombic interactions within the electrodes. This dielectric is implemented here by dividing all charges by $\epsilon r$, multiplying potentials with $\epsilon r$, and multiplying the external resistance by

*ɛ*

_{r}; this scaling is reversed in the simulation values reported below. In the absence of an explicit solvent exposing the ions to friction and perpetual Brownian motion, its impact on the dynamics of the ions is incorporated into the Langevin equation of motion,

*m*= 50 a.m.u. being the mass of the ions, the dots denoting time derivatives,

*γ*= 2.26 a.m.u./fs being the Stokesian drag coefficient for a sphere of diameter

*σ*in water,

*f*_{i}being the conservative force acting on the ion, and $fiR$ being a random force; the latter has zero mean, no time correlations (Markovian), no correlations between ions, and a variance that is related by the fluctuation–dissipation theorem to the drag coefficient and the temperature. The integration time step is set at 0.5 to 5 fs.

## IV. RESULTS AND DISCUSSION

### A. Equilibrium

The electrode charge vs voltage curves of the system with 156 ion pairs, as calculated using CPM at various potentials and CSCM at various electrode sum charges, are shown in Fig. 3. The good agreement between both curves implies that both methods agree on the capacitance and differential capacitance of this capacitor, despite sampling distinct ensembles; we also see excellent agreement with classical Density Functional Theory (DFT) calculations on a system closely resembling the current system.^{53} A CPM simulation at fixed Δ*V* = 0.4 V yields an average electrode charge of ⟨*Q*⟩_{ΔV} = 40.68 *e*; a subsequent CSCM simulation fixing the electrode charge at this value results in an average potential difference of ⟨Δ*V*⟩_{Q} = 0.401 V, in good agreement with the preceding CPM simulation. The distribution of the instantaneous potential relative to this average follows a Gaussian distribution, as illustrated in Fig. 4. Its standard deviation is related to the differential capacitance, as derived in Eq. (26). The ion-free capacitor contributes *C*_{0} = 27.32 *e*/V, which agrees very well with the idealized value of *C*_{H} = *ɛ*_{0}*ɛ*_{r}*A*/*H* = 27.10 *e*/V for this system; upon normalization by the area of the electrode, the former corresponds to *C*_{0}/*A* = 17.41 *μ*F/cm^{2}. The resulting differential capacitance is plotted in Fig. 5 against the average potential, as obtained by simulations at several fixed values of the electrode charge. A good agreement is observed with the differential capacitance deduced from the charge fluctuations at a constant potential difference;^{53} see Eq. (20). The line in Fig. 5 shows the numerical derivative of the charge–voltage characteristic obtained by (DFT) calculations on a system closely resembling the current system.^{53}

In view of the above agreement between CPM and CSCM, it is to be expected that they will also agree on ion distributions between the electrodes. Figure 6 shows three examples illustrating the excellent agreement, even though, as stated before, the two methods sample different ensembles. In each case, a CPM simulation at a fixed potential difference is compared against a CSCM simulation with the surface charge fixed at the average surface charge in the CPM simulation. The resulting average potentials in the CSCM are close to their fixed values in the CPM, at 0.203, 0.401, and 1.999 V vs 0.2, 0.4, and 2.0 V, respectively. The ion distributions of the first and third examples also agree with DFT calculations on nearly identical systems; see Figs. 4 and 5(f) in Ref. 53, respectively. With an increasing potential difference, the ions form layers of increasing compactness against the electrodes. For the 156 ion pair system—which will be the focus of the remainder of this study—at 0.4 V [see Fig. 6(b)], the layer adjacent to the electrode (as defined by the two minima at either side of the density peak) contains about 40 counterions, with the remaining ∼110 counterions fairly homogeneously distributed between the two counterion peaks. Because the counterions next to the electrodes at 2 V form a closely packed hexagonal layer that does not fully compensate for the electrode charge, a non-saturated second layer of counterions forms; see Fig. 6(c). A small dozen of ion pairs are free to roam in the interior of the capacitor, although their non-uniform distribution suggests they are not completely screened from the electrode charges.

### B. Charging and discharging

The capacitor will gradually discharge when the two electrodes are connected by an external resistance. Figure 7 shows the potential difference, Δ*V*, and the charge on the electrodes, *Q*, as functions of time for the 156 ion pair system discharging from equilibrated configurations at initial charges of *Q*_{0} = 40.68 and 142.70 *e*, corresponding to Δ*V* ≈ 0.4 and $\u22482.0$ V, respectively. The strength *R*_{e} of the external resistor has been varied in the CSCM simulations to sample several orders of magnitude in external current relaxation time, estimated as *τ*_{e} = *R*_{e}*C*_{0}, with *C*_{0} the capacitance of the ion-free capacitor. A relaxation time of 100 fs is realized with a resistance of 3.66 V fs/*e* = 22.8 kΩ. The shortest relaxation time is realized using CPM, where the potential difference is instantaneously reduced to zero; this elimination of the driving force of the external current by a short circuit, *R*_{e} = 0, implies *τ*_{e} = 0. Removing the potential difference from the charge calculation in Eq. (5), or its charge-neutral equivalent obtained by replacing $Ap\u22121$ with **S**_{p},^{43,46} does not instantaneously drop the electrode charges to zero. Rather, the pull by the electrolyte ions binds a residual *Q*_{R} ≈ (2/3)*Q*_{0} at the electrodes; see Figs. 7(c) and 7(d). The capacitor is now in a non-equilibrium situation where reduced electrode charges allow ions to diffuse out of the double layer, which in turn allows charge to flow through the external circuit. The resulting near-exponential decay of the electrode charge is represented by the green lines in Figs. 7(c) and 7(d). This process continues until the charges on the electrodes merely fluctuate around zero due to the perpetual thermal motions of the ions.

*V*) or nearly overlapping lines (

*Q*). To quantify the relaxation processes, the data are fitted with a double-exponential decay,

*y*=

*Q*or Δ

*V*, two amplitudes

*A*

_{y,1}and

*A*

_{y,2}, and two relaxation times

*τ*

_{y,1}and

*τ*

_{y,2}. As noted in Sec. II C, this is an analytic solution to our model circuit. The resulting fit parameters are collected in Table I, where the ordinal number 1 is assigned to the mode with the largest relaxation time,

*τ*

_{y,1}>

*τ*

_{y,2}. It can be seen there that the relaxation time of the second mode strongly resembles the external relaxation time,

*τ*

_{y,2}≈

*τ*

_{e}, while the relaxation time of the first mode appears to be independent of the external relaxation time, especially for the fits to the electrode charges (in general, the charge curves yield better fits than the potential curves because the former is much less affected by thermal noise than the latter; see Fig. 4). This suggests that the second mode represents relaxation by the external current, while the first mode relates to the dynamics of the ions. To validate the latter attribution, we altered the mobility of the ions at constant external relaxation time. Upon increasing (decreasing) the friction parameter

*γ*in the Langevin equation of motion by a factor of 10, the relaxation times

*τ*

_{y,1}decrease (increase) by a factor of ∼7.5, while

*τ*

_{y,2}remains essentially constant; see Table II. The modest disagreement between these two scaling factors is due to the conservative force in the second order Langevin equation not being proportional to the friction; varying

*γ*changes the relative importance of conservative and solvent forces rather than uniformly rescaling all forces. The table also shows that the amplitudes are hardly affected by the friction. We furthermore note that the single exponential relaxation process simulated by CPM closely resembles the first mode of the double exponential relaxation process in CSCM, thereby confirming the attribution of internal and external dynamics.

Q_{0} (e)
. | y
. | τ_{e} (fs)
. | A_{y,1} (y)
. | τ_{y,1} (fs)
. | A_{y,2} (y)
. | τ_{y,2} (fs)
. | k_{QQ} (V/e)
. | k_{Qξ} (V/e)
. | k_{ξξ} (V/e)
. | M_{ξ} [e/(V fs)]
. |
---|---|---|---|---|---|---|---|---|---|---|

40.68 | Q | 0 | 2.98 × 10^{1} | 2.24 × 10^{5} | ⋯ | ⋯ | ||||

10^{2} | 2.92 × 10^{1} | 2.12 × 10^{5} | 1.13 × 10^{1} | 1.12 × 10^{2} | 8.98 × 10^{−3} | 2.36 × 10^{−2} | 1.06 × 10^{−4} | 7.14 × 10^{−4} | ||

10^{3} | 2.92 × 10^{1} | 2.33 × 10^{5} | 1.07 × 10^{1} | 1.28 × 10^{3} | 7.60 × 10^{−3} | 2.06 × 10^{−2} | 1.06 × 10^{−4} | 7.70 × 10^{−4} | ||

10^{4} | 3.24 × 10^{1} | 2.04 ⋅ 10^{5} | 8.42 ⋅ 10^{0} | 8.33 ⋅ 10^{3} | 1.04 ⋅ 10^{−2} | 2.78 ⋅ 10^{−2} | 1.06 ⋅ 10^{−4} | 7.35 ⋅ 10^{−4} | ||

ΔV | 10^{2} | 6.22 × 10^{−4} | 6.87 × 10^{4} | 6.12 × 10^{−1} | 7.95 × 10^{1} | 2.03 × 10^{−2} | 1.36 × 10^{−2} | 1.02 × 10^{−4} | 1.50 × 10^{−3} | |

10^{3} | 4.97 × 10^{−3} | 2.09 × 10^{5} | 3.73 × 10^{−1} | 1.04 × 10^{3} | 9.27 × 10^{−3} | 2.51 × 10^{−2} | 9.88 × 10^{−5} | 6.88 × 10^{−4} | ||

10^{4} | 5.83 × 10^{−2} | 1.81 × 10^{5} | 3.60 × 10^{−1} | 8.54 × 10^{3} | 1.04 × 10^{−2} | 2.61 × 10^{−2} | 9.91 × 10^{−5} | 8.07 × 10^{−4} | ||

142.70 | Q | 0 | 8.54 × 10^{1} | 1.40 × 10^{5} | ⋯ | ⋯ | ||||

10^{2} | 8.50 × 10^{1} | 1.44 × 10^{5} | 5.51 × 10^{1} | 1.24 × 10^{2} | 1.14 × 10^{−2} | 1.80 × 10^{−2} | 1.07 × 10^{−4} | 9.82 × 10^{−4} | ||

10^{3} | 8.63 × 10^{1} | 1.44 × 10^{5} | 5.48 × 10^{1} | 1.17 × 10^{3} | 1.22 × 10^{−2} | 1.87 × 10^{−2} | 1.07 × 10^{−4} | 9.39 × 10^{−4} | ||

10^{4} | 9.52 × 10^{1} | 1.63 × 10^{5} | 4.63 × 10^{1} | 1.06 × 10^{4} | 1.27 × 10^{−2} | 1.81 × 10^{−2} | 1.07 × 10^{−4} | 9.07 × 10^{−4} | ||

ΔV | 10^{2} | 2.67 × 10^{−3} | 8.64 × 10^{4} | 1.75 × 10^{0} | 1.08 × 10^{2} | 1.34 × 10^{−2} | 1.76 × 10^{−2} | 1.07 × 10^{−4} | 1.33 × 10^{−3} | |

10^{3} | 2.70 × 10^{−2} | 9.02 × 10^{4} | 1.95 × 10^{0} | 1.02 × 10^{3} | 1.48 × 10^{−2} | 1.83V10^{−2} | 1.07 × 10^{−4} | 1.24 × 10^{−3} | ||

10^{4} | 2.96 × 10^{−1} | 9.34 × 10^{4} | 1.71 × 10^{0} | 8.77 × 10^{3} | 1.49 × 10^{−2} | 1.90 × 10^{−2} | 1.07 × 10^{−4} | 1.30 × 10^{−3} | ||

⟨⋯⟩ | Q | 1.01 × 10^{−2} | 1.98 × 10^{−2} | 1.06 × 10^{−4} | 8.93 × 10^{−4} |

Q_{0} (e)
. | y
. | τ_{e} (fs)
. | A_{y,1} (y)
. | τ_{y,1} (fs)
. | A_{y,2} (y)
. | τ_{y,2} (fs)
. | k_{QQ} (V/e)
. | k_{Qξ} (V/e)
. | k_{ξξ} (V/e)
. | M_{ξ} [e/(V fs)]
. |
---|---|---|---|---|---|---|---|---|---|---|

40.68 | Q | 0 | 2.98 × 10^{1} | 2.24 × 10^{5} | ⋯ | ⋯ | ||||

10^{2} | 2.92 × 10^{1} | 2.12 × 10^{5} | 1.13 × 10^{1} | 1.12 × 10^{2} | 8.98 × 10^{−3} | 2.36 × 10^{−2} | 1.06 × 10^{−4} | 7.14 × 10^{−4} | ||

10^{3} | 2.92 × 10^{1} | 2.33 × 10^{5} | 1.07 × 10^{1} | 1.28 × 10^{3} | 7.60 × 10^{−3} | 2.06 × 10^{−2} | 1.06 × 10^{−4} | 7.70 × 10^{−4} | ||

10^{4} | 3.24 × 10^{1} | 2.04 ⋅ 10^{5} | 8.42 ⋅ 10^{0} | 8.33 ⋅ 10^{3} | 1.04 ⋅ 10^{−2} | 2.78 ⋅ 10^{−2} | 1.06 ⋅ 10^{−4} | 7.35 ⋅ 10^{−4} | ||

ΔV | 10^{2} | 6.22 × 10^{−4} | 6.87 × 10^{4} | 6.12 × 10^{−1} | 7.95 × 10^{1} | 2.03 × 10^{−2} | 1.36 × 10^{−2} | 1.02 × 10^{−4} | 1.50 × 10^{−3} | |

10^{3} | 4.97 × 10^{−3} | 2.09 × 10^{5} | 3.73 × 10^{−1} | 1.04 × 10^{3} | 9.27 × 10^{−3} | 2.51 × 10^{−2} | 9.88 × 10^{−5} | 6.88 × 10^{−4} | ||

10^{4} | 5.83 × 10^{−2} | 1.81 × 10^{5} | 3.60 × 10^{−1} | 8.54 × 10^{3} | 1.04 × 10^{−2} | 2.61 × 10^{−2} | 9.91 × 10^{−5} | 8.07 × 10^{−4} | ||

142.70 | Q | 0 | 8.54 × 10^{1} | 1.40 × 10^{5} | ⋯ | ⋯ | ||||

10^{2} | 8.50 × 10^{1} | 1.44 × 10^{5} | 5.51 × 10^{1} | 1.24 × 10^{2} | 1.14 × 10^{−2} | 1.80 × 10^{−2} | 1.07 × 10^{−4} | 9.82 × 10^{−4} | ||

10^{3} | 8.63 × 10^{1} | 1.44 × 10^{5} | 5.48 × 10^{1} | 1.17 × 10^{3} | 1.22 × 10^{−2} | 1.87 × 10^{−2} | 1.07 × 10^{−4} | 9.39 × 10^{−4} | ||

10^{4} | 9.52 × 10^{1} | 1.63 × 10^{5} | 4.63 × 10^{1} | 1.06 × 10^{4} | 1.27 × 10^{−2} | 1.81 × 10^{−2} | 1.07 × 10^{−4} | 9.07 × 10^{−4} | ||

ΔV | 10^{2} | 2.67 × 10^{−3} | 8.64 × 10^{4} | 1.75 × 10^{0} | 1.08 × 10^{2} | 1.34 × 10^{−2} | 1.76 × 10^{−2} | 1.07 × 10^{−4} | 1.33 × 10^{−3} | |

10^{3} | 2.70 × 10^{−2} | 9.02 × 10^{4} | 1.95 × 10^{0} | 1.02 × 10^{3} | 1.48 × 10^{−2} | 1.83V10^{−2} | 1.07 × 10^{−4} | 1.24 × 10^{−3} | ||

10^{4} | 2.96 × 10^{−1} | 9.34 × 10^{4} | 1.71 × 10^{0} | 8.77 × 10^{3} | 1.49 × 10^{−2} | 1.90 × 10^{−2} | 1.07 × 10^{−4} | 1.30 × 10^{−3} | ||

⟨⋯⟩ | Q | 1.01 × 10^{−2} | 1.98 × 10^{−2} | 1.06 × 10^{−4} | 8.93 × 10^{−4} |

γ (a.m.u./fs)
. | y
. | A_{y,1} (y)
. | τ_{y,1} (fs)
. | A_{y,2} (y)
. | τ_{y,2} (fs)
. | k_{QQ} (V/e)
. | k_{Qξ} (V/e)
. | k_{ξξ} (V/e)
. | M_{ξ} [e/(V fs)]
. |
---|---|---|---|---|---|---|---|---|---|

0.226 | Q | 8.72 × 10^{1} | 1.02 × 10^{6} | 5.50 × 10^{1} | 1.03 × 10^{3} | 1.36 × 10^{−2} | 2.18 × 10^{−2} | 1.00 × 10^{−4} | 1.16 × 10^{−4} |

ΔV | 4.03 × 10^{−3} | 1.32 × 10^{5} | 1.83 × 10^{0} | 1.05 × 10^{3} | 1.98 × 10^{−2} | 2.57 × 10^{−3} | 2.32 × 10^{−4} | 2.93 × 10^{−3} | |

2.26 | Q | 8.63 × 10^{1} | 1.44 × 10^{5} | 5.48 × 10^{1} | 1.17 × 10^{3} | 1.22 × 10^{−2} | 1.87 × 10^{−2} | 1.07 × 10^{−4} | 9.39 × 10^{−4} |

ΔV | 2.70 × 10^{−2} | 9.02 × 10^{4} | 1.95 × 10^{0} | 1.02 × 10^{3} | 1.48 × 10^{−2} | 1.83 × 10^{−2} | 1.07 × 10^{−4} | 1.24 × 10^{−3} | |

22.6 | Q | 9.61 × 10^{1} | 1.87 × 10^{4} | 4.62 × 10^{1} | 9.74 × 10^{2} | 1.34 × 10^{−2} | 2.07 × 10^{−2} | 1.08 × 10^{−4} | 7.16 × 10^{−3} |

ΔV | 1.97 × 10^{−1} | 1.72 × 10^{4} | 1.80 × 10^{0} | 9.24 × 10^{2} | 1.40 × 10^{−2} | 2.15 × 0^{−2} | 1.08 × 10^{−4} | 7.31 × 10^{−3} |

γ (a.m.u./fs)
. | y
. | A_{y,1} (y)
. | τ_{y,1} (fs)
. | A_{y,2} (y)
. | τ_{y,2} (fs)
. | k_{QQ} (V/e)
. | k_{Qξ} (V/e)
. | k_{ξξ} (V/e)
. | M_{ξ} [e/(V fs)]
. |
---|---|---|---|---|---|---|---|---|---|

0.226 | Q | 8.72 × 10^{1} | 1.02 × 10^{6} | 5.50 × 10^{1} | 1.03 × 10^{3} | 1.36 × 10^{−2} | 2.18 × 10^{−2} | 1.00 × 10^{−4} | 1.16 × 10^{−4} |

ΔV | 4.03 × 10^{−3} | 1.32 × 10^{5} | 1.83 × 10^{0} | 1.05 × 10^{3} | 1.98 × 10^{−2} | 2.57 × 10^{−3} | 2.32 × 10^{−4} | 2.93 × 10^{−3} | |

2.26 | Q | 8.63 × 10^{1} | 1.44 × 10^{5} | 5.48 × 10^{1} | 1.17 × 10^{3} | 1.22 × 10^{−2} | 1.87 × 10^{−2} | 1.07 × 10^{−4} | 9.39 × 10^{−4} |

ΔV | 2.70 × 10^{−2} | 9.02 × 10^{4} | 1.95 × 10^{0} | 1.02 × 10^{3} | 1.48 × 10^{−2} | 1.83 × 10^{−2} | 1.07 × 10^{−4} | 1.24 × 10^{−3} | |

22.6 | Q | 9.61 × 10^{1} | 1.87 × 10^{4} | 4.62 × 10^{1} | 9.74 × 10^{2} | 1.34 × 10^{−2} | 2.07 × 10^{−2} | 1.08 × 10^{−4} | 7.16 × 10^{−3} |

ΔV | 1.97 × 10^{−1} | 1.72 × 10^{4} | 1.80 × 10^{0} | 9.24 × 10^{2} | 1.40 × 10^{−2} | 2.15 × 0^{−2} | 1.08 × 10^{−4} | 7.31 × 10^{−3} |

*τ*

_{e}≫

*τ*

_{y,1}, the observed relaxation process ought to converge to the ideal macroscopic discharge curve for a capacitor continuously in equilibrium. The latter is obtained as

*V*(

*t*) is numerically solved using matlab’s Runge–Kutta routine ode45,

^{54}based on a linear interpolation of the smooth differential capacitance curve calculated by DFT; see Fig. 5. Integration of this result, using the second equality in the above equation, then yields the charge evolution,

*Q*(

*t*). Figure 8 shows that the simulation results indeed converge to this theoretical curve with increasing external relaxation time. A close agreement between simulations and the ideal system is realized only when the external relaxation time well exceeds the capacitor’s intrinsic relaxation time

*τ*

_{y,1}, which is reached for

*τ*

_{e}= 10

^{6}fs (green lines). These discharge curves also allow for a calculation of the total energy delivered up to time

*t*,

^{54}This energy increases when the capacitor is discharging through a larger external resistance (see Fig. 9), with the upper limit set by the ideal continuous-equilibrium process. Note that these energy curves nearly coalesce for short rescaled times, with the common initial power delivery

*P*(0) ≈ 70

*e*V/(

*R*

_{e}

*C*

_{0}). Conventional theory gives the power as

*P*= (Δ

*V*)

^{2}/

*R*

_{e}, which under the present conditions translates into

*P*(0) ≈ 109 eV/(

*R*

_{e}

*C*

_{0}), in reasonable agreement with the simulations. The supercapacitor can deliver this power over a period of about 0.5

*τ*

_{e}for small external resistances, or about

*τ*

_{e}for large external resistances, after which time the capacitor is essentially drained.

*V*

_{e}= 0.4 or 2.0 V and a resistance of

*R*

_{e}. In CPM, the potential differences between the electrodes jump instantaneously to their final values,

*V*

_{e}; see Figs. 10(a) and 10(b). The electrode charges likewise jump to plateau values [see Figs. 10(c) and 10(d)], which correspond to about 25% and 40% of their final equilibrium values at

*V*

_{e}= 0.4 and 2.0 V, respectively, followed by a slow rise on a timescale

*τ*

_{y,1}comparable to its counterpart observed under discharging. In CSCM, the transient jumps are replaced by smooth, rapidly rising curves converging to similar plateau values, followed by slow growths with common slopes on the logarithmic scale. The potential difference and electrode charge are reasonably well fitted with the double exponential function

*τ*

_{y,2}≈

*τ*

_{e}, while the relaxation time of the first mode is fairly constant. In the simulation with an external source at 2.0 V, the fitted amplitudes for the charge are remarkably close to their counterparts for the discharging simulations from this potential. This observation suggests that the discharging dynamics, with peaked ion distributions collapsing into a homogeneous distribution, and the charging dynamics, with a homogeneous ion distribution contracting into sharp peaks, can be described on the macroscopic level by the same evolution equations. This is a confirmation of the tacit assumption underlying the phenomenological model in Eq. (34), i.e., the equivalent circuit.

V_{e} (V)
. | y
. | τ_{e} (fs)
. | A_{y,1} (y)
. | τ_{y,1} (fs)
. | A_{y,2} (y)
. | τ_{y,2} (fs)
. |
---|---|---|---|---|---|---|

0.4 | Q | 0 | 2.22 × 10^{1} | 1.03 × 10^{5} | 1.08 × 10^{1} | 0 |

10^{2} | 1.06 × 10^{2} | 7.46 × 10^{5} | 1.09 × 10^{1} | 1.05 × 10^{2} | ||

10^{3} | 1.80 × 10^{1} | 9.33 × 10^{4} | 1.01 × 10^{1} | 9.30 × 10^{2} | ||

10^{4} | 1.61 × 10^{1} | 4.81 × 10^{4} | 5.22 × 10^{0} | 6.96 × 10^{3} | ||

ΔV | 10^{2} | 3.89 × 10^{−1} | 1.01 × 10^{2} | 1.09 × 10^{−2} | 2.19 × 10^{0} | |

10^{3} | 6.93 × 10^{−3} | 1.89 × 10^{4} | 3.90 × 10^{−1} | 8.83 × 10^{2} | ||

10^{4} | 3.54 × 10^{−1} | 9.97 × 10^{3} | 8.61 × 10^{−3} | 3.20 × 10^{−1} | ||

2.0 | Q | 0 | 8.88 × 10^{1} | 8.99 × 10^{4} | 5.53 × 10^{1} | 0 |

10^{2} | 8.90 × 10^{1} | 9.51 × 10^{4} | 5.47 × 10^{1} | 1.01 × 10^{2} | ||

10^{3} | 8.57 × 10^{1} | 8.22 × 10^{4} | 5.27 × 10^{1} | 9.84 × 10^{2} | ||

10^{4} | 1.08 × 10^{2} | 1.02 × 10^{5} | 3.49 × 10^{1} | 7.91 × 10^{3} | ||

ΔV | 10^{2} | 3.55 × 10^{−3} | 8.53 × 10^{4} | 2.00 × 10^{0} | 9.83 × 10^{1} | |

10^{3} | 3.93 × 10^{−2} | 7.45 × 10^{4} | 1.96 × 10^{0} | 9.62 × 10^{2} | ||

10^{4} | 4.01 × 10^{−1} | 1.08 × 10^{5} | 1.61 × 10^{0} | 7.85 × 10^{3} |

V_{e} (V)
. | y
. | τ_{e} (fs)
. | A_{y,1} (y)
. | τ_{y,1} (fs)
. | A_{y,2} (y)
. | τ_{y,2} (fs)
. |
---|---|---|---|---|---|---|

0.4 | Q | 0 | 2.22 × 10^{1} | 1.03 × 10^{5} | 1.08 × 10^{1} | 0 |

10^{2} | 1.06 × 10^{2} | 7.46 × 10^{5} | 1.09 × 10^{1} | 1.05 × 10^{2} | ||

10^{3} | 1.80 × 10^{1} | 9.33 × 10^{4} | 1.01 × 10^{1} | 9.30 × 10^{2} | ||

10^{4} | 1.61 × 10^{1} | 4.81 × 10^{4} | 5.22 × 10^{0} | 6.96 × 10^{3} | ||

ΔV | 10^{2} | 3.89 × 10^{−1} | 1.01 × 10^{2} | 1.09 × 10^{−2} | 2.19 × 10^{0} | |

10^{3} | 6.93 × 10^{−3} | 1.89 × 10^{4} | 3.90 × 10^{−1} | 8.83 × 10^{2} | ||

10^{4} | 3.54 × 10^{−1} | 9.97 × 10^{3} | 8.61 × 10^{−3} | 3.20 × 10^{−1} | ||

2.0 | Q | 0 | 8.88 × 10^{1} | 8.99 × 10^{4} | 5.53 × 10^{1} | 0 |

10^{2} | 8.90 × 10^{1} | 9.51 × 10^{4} | 5.47 × 10^{1} | 1.01 × 10^{2} | ||

10^{3} | 8.57 × 10^{1} | 8.22 × 10^{4} | 5.27 × 10^{1} | 9.84 × 10^{2} | ||

10^{4} | 1.08 × 10^{2} | 1.02 × 10^{5} | 3.49 × 10^{1} | 7.91 × 10^{3} | ||

ΔV | 10^{2} | 3.55 × 10^{−3} | 8.53 × 10^{4} | 2.00 × 10^{0} | 9.83 × 10^{1} | |

10^{3} | 3.93 × 10^{−2} | 7.45 × 10^{4} | 1.96 × 10^{0} | 9.62 × 10^{2} | ||

10^{4} | 4.01 × 10^{−1} | 1.08 × 10^{5} | 1.61 × 10^{0} | 7.85 × 10^{3} |

*t*decreases with increasing external relaxation time, while that transferred to the resistor increases with increasing

*τ*

_{e}. The power provided by the voltage source is given by

*P*

_{e}(

*t*) =

*V*

_{e}

*I*(

*t*), irrespective of the distribution of this power over the resistor and capacitance. The total energy provided to reach equilibrium follows from

*Q*(

*V*

_{e}) is the equilibrium charge at the imposed potential difference; see Fig. 3. Interestingly, this energy is a function of

*V*

_{e}only and does not depend on the external resistor. For a capacitor that has reached equilibrium at

*V*

_{e}, the maximum energy it can deliver is also a function of this

*V*

_{e}. It then follows that the energy lost while charging to equilibrium, i.e., the energy that cannot be retrieved upon discharging, depends on

*V*

_{e}but not on

*R*

_{e}. However, the external resistor does determine the relative distribution of the loss over the two circuit components. The loss in the capacitor is clearly visible in Fig. 11, where the energy absorbed during charging systematically exceeds the maximum energy of 124 eV returned during discharging (see Fig. 9); this difference is lost as heat to the environment, which is modeled in the simulations by the thermostatting property of the Langevin equation of motion. The loss in the capacitor can be made to vanish by maintaining equilibrium, i.e., in the limit of very slow charging. For an ideal conventional capacitor, with

*Q*(

*V*

_{e}) =

*CV*

_{e}, the energy stored in the capacitor reads as $Ec(Ve)=CVe2/2$ and it follows from Eq. (46) that an equal amount is lost as heat in the resistor.

^{55}The energy stored in the supercapacitor can be expressed as

*V*

_{e}

*Q*(

*V*

_{e}) =

*E*

_{e}of a rectangle in the

*Q*–

*V*plane. Equation (47) highlights that the unavoidable energy loss during charging in this setup is given by the integral of

*Q*(

*V*), i.e., the area under the charge–voltage curve in Fig. 3. For an ideal conventional capacitor, this curve is a straight line, equally splitting the supplied energy into stored and lost energy. The concave curve of the present supercapacitor, with a differential capacity that decreases with increasing voltage, makes the total loss of 162 eV exceed the storage of 124 eV; a convex shape, with a differential capacity that increases with voltage, would be advantageous in the present charging procedure. If the capacitor is removed from the charging circuit after a finite time

*t*, the total energy supplied by the external source,

*E*

_{e}(

*t*), will be less than the limit in Eq. (47). In the disconnected state, the sum charges

*Q*on the electrodes remain constant while the ion distributions may still evolve. The final equilibrium of a thermostatted capacitor will be attained for a voltage that can be read off in Fig. 3 at that

*Q*, and said curve also enables a calculation of the stored energy,

*E*

_{c}. Plotting this

*E*

_{c}against

*E*

_{e}(

*t*) (see Fig. 12) highlights that the energy supplied at the start of the simulation, i.e., for Δ

*V*≪

*V*

_{e}, is largely lost; the energy supplied at the end of the simulation, i.e., for Δ

*V*≈

*V*

_{e}, is efficiently stored. Interestingly, the curves for three relaxation times nearly overlap, suggesting that this system is remarkably insensitive to the charging rate. Note that a capacitor can be charged without energy loss to a resistor by directly connecting the capacitor to an external source and smoothly increasing the potential of that source from zero to the desired final value,

^{55}although this will not eliminate losses caused by non-equilibrium conditions in the capacitor.

### C. Equivalent circuit

*A*

_{Q,i}and relaxation times

*τ*

_{Q,i}into the three energy coefficients

*k*and friction parameter

*M*

_{ξ}of the equivalent circuit [see Eqs. (29) and (34)] is involved: the proposed equivalent circuit is readily solved to express the theoretical amplitudes

*A*

_{Q,±}and relaxation times

*τ*

_{±}in terms of the energy coefficients and friction parameter by using the formal solution in Appendix B, but inverting these equations into explicit expressions for the

*k*’s and

*M*

_{ξ}does not look promising. Therefore, we numerically minimize the function

*A*

_{Q,±}and

*τ*

_{±}, where the subscripts in the denominators are to be read as ±(1) = + and ±(2) = −. Minimization in Microsoft Excel results in two solutions for each data set, depending on the initial guesses of the

*k*’s and

*M*

_{ξ}. The two solutions differ in the value of

*k*

_{ξξ}, which converges to either ≈1 × 10

^{−4}V/

*e*or ≈2 × 10

^{−4}V/

*e*, while the values for the other parameters differ only slightly between the two minima. These fitted values of

*k*

_{ξξ}are two orders of magnitude smaller than the fitted values of

*k*

_{QQ}and

*k*

_{Qξ}, indicating that the

*ξ*

^{2}contribution to the free energy of Eq. (29) is negligible. Furthermore, the difference in the resulting amplitudes and relaxation times upon using either value of

*k*

_{ξξ}is negligible, and even inserting

*k*

_{ξξ}= 0 changes their values by a few percent at most. We, therefore, conclude that there is little to choose between the two fits and that the capacitor

*C*

_{ξξ}may be removed from the equivalent circuit in Fig. 2, thereby arriving at the simplified equivalent circuit in Fig. 14. Equivalent circuit parameters extracted from the discharge simulations in Fig. 7 are collected in Table I; the parameters for the aforementioned wider collection of discharge simulations are presented in Fig. 13. They all vary within ∼20% of their respective averages, with the trendlines in the figure suggesting a weak dependence on the initial charge

*Q*

_{0}; there does not appear to be a systematic dependence on the external resistance; see also Table I. The average values extracted from these data are reported in the bottom line of Table I. For the components of the equivalent circuit in Fig. 14, we then obtain $CQQ\u2032=2.0\xd7102e$/V = 1.2 × 10

^{2}

*μ*F/cm

^{2},

*C*

_{Qξ}= 51

*e*/V = 32

*μ*F/cm

^{2},

*C*

_{ξξ}= 0, and $R\xi \u2032=5.6\xd7102$ V fs/

*e*= 3.5 × 10

^{5}Ω. The first of these values agrees remarkably well with the capacitance of

*C*

_{d}=

*ɛ*

_{0}

*ɛ*

_{r}

*A*/

*d*= 2.2 × 10

^{2}

*e*/V expected for two homogeneous surface charges separated by a distance

*d*[see Fig. 1 and Eq. (28)] when this distance is approximated as the radius

*σ*of the WCA potential. The equilibrium capacitance of the simplified equivalent circuit is then evaluated as $Ceq=CQQ\u2032/2=62\mu $F/cm

^{2}, similar to the differential capacitance in the low-voltage range in Fig. 5. While the simple model of Eq. (28) predicts

*C*

_{Qξ}to be equal to

*H*≈

*w*≫

*d*, the actual value is almost twice the capacitance of the ion-free capacitor. The macroscopic ion mobility

*M*

_{ξ}is related to the microscopic mobility of individual ions, which in turn is related to the ionic self-diffusion coefficient of

*D*

_{ion}= 6.56 × 10

^{−6}cm

^{2}/s for the bulk electrolyte, as discussed in Appendix C. The discharge curves of the equivalent circuit using these average parameters match reasonably well with the simulations at three external resistances; see the dashed lines in Fig. 7. The same parameters also yield a reasonable agreement with the charging simulations (see the dashed lines in Fig. 10), although these simulations are not included in the extraction of the circuit components.

Several arguments may explain the modest discrepancies between the equivalent circuit and the molecular dynamics simulations. The model predicts a constant differential capacitance [see Eq. (32)] in contrast to the simulations; see Fig. 5. Evidently, the phenomenological equations underlying the equivalent circuit are merely crude approximations, ignoring any details on the actual ion distributions. The derivation in Appendix C suggests that the ion mobility *M*_{ξ} will vary with the ion distribution. The argumentation behind Eqs. (28) and (29) may constitute a justifiable model for ions forming a monolayer adjacent to each electrode, but in the simulated capacitor, an incomplete second layer develops at the higher electrode charges; see Fig. 6(c). A closer look at the supercapacitor discharging from the highest initial charge (see Fig. 15) reveals signs of a triple-exponential decay. Upon fitting the charge relaxation curve over a longer time span and a larger set of simulations, we obtain relaxation times of 1.03 × 10^{2}, 4.62 × 10^{4}, and 2.20 × 10^{5} fs, with corresponding amplitudes of 55.0, 23.0, and 64.4 *e*. The first time again matches the fast external relaxation time of 10^{2} fs, while the last time lies remarkably close to the slow internal time *τ*_{Q,1} = 2.12 × 10^{5} fs when discharging from a low initial charge of 40.68 *e* (see Table I) with only a partial monolayer adjacent to the electrodes; see Fig. 6(b). It is, therefore, tempting to speculate that this largest relaxation time is related to the ionic monolayer nearest the electrode, while the intermediate relaxation time applies to the weaker bound next-nearest ionic layer. The *τ*_{Q,1} = 1.44 × 10^{5} fs reported in Table I for discharging from 142.70 *e*, which lies between the two largest times of the triple fit, would then appear as an effective relaxation time accounting for these two layers collectively, and the sum of their amplitudes, 87.4 *e*, matches with the *A*_{Q,1} in Table I. However, whereas Fig. 15 locates the transition between the second and third regimes at around 100 ps, visual inspection of instantaneous density profiles suggests that the second layer has disappeared by about 30 ps.

Earlier, we noted that a ten-fold increase (decrease) of the friction parameter *γ* in the equation of motion resulted in a ∼7.5-fold decrease (increase) of the relaxation times *τ*_{y,1}, but hardly affected the other fit parameters. From the analysis in Table II, especially from the charge relaxation, it is clear that these changes of *γ* generate a ≈8-fold increase (decrease) in *M*_{ξ}, with no systematic impact on the other three circuit parameters. This confirms the identification of *M*_{ξ} as related to the mobility of ions in the electrolyte.

An additional set of simulations is undertaken using a distance between the electrodes tripled to *H* = 12 nm. For this system, the number of ion pairs is increased to 520 to recover the previous ion densities in the center and at the electrodes of the reference system of 156 ion pairs at *V*_{e} = 0.4 V. The fit parameters extracted from the double exponential discharge curves, for relaxation times from 10^{2} to 10^{4} fs, yield a slow mode that is more pronounced and slowed down relative to its counterpart in Table I, with *A*_{Q,1} ≈ 38 *e* and *τ*_{Q,1} ≈ 5 × 10^{5} fs, while the fast mode has been reduced to *A*_{Q,2} ≈ 3 *e* with relaxation times about three times shorter than those in Table I. The extracted *C*_{QQ} for the equivalent circuit remains unchanged, indicating that the electric double layers are not affected by the distance between the electrodes, while the *C*_{Qξ} due to interactions between the two double layers is reduced by a factor of about four, in line with the earlier expectation that this capacitance scales inversely with the distance. The increased distance raises the internal resistance by a factor of 3.5, thereby suggesting that *R*_{ξ} ∝ *H* in agreement with Eq. (C8).

The equivalent circuit enables the extension of the analysis of relaxation times beyond the range manageable by MD simulations. In Fig. 16(a), the relaxation times measured in the simulations, *τ*_{y,i}, are compared with their theoretical counterparts in the equivalent circuit, *τ*_{±}. It becomes clear that the simulation data used to fit the equivalent circuit, with external relaxation times *τ*_{e} ≤ 10^{4} fs (see Table I and Fig. 13) are in the limit of low external resistance. Here, the short relaxation time *τ*_{y,1} is linear with the external resistance, while the long relaxation time *τ*_{y,2} is independent of this resistance; the detailed formulas are collected in Appendix D. The predicted amplitudes of the charge decay modes, *A*_{Q,±}, are also independent of *R*_{e} in this limit, as are the measured relaxation amplitudes, *A*_{Q,i}; see Fig. 16(b). We note that the normalized amplitudes *A*_{Q,i}/*Q*(0) for *Q*(0) = 40.68 and 142.70 *e* differ systematically—compare the triangles pointing up (left) with those pointing down (right)—unlike the prediction by the linear model; it, therefore, appears that the simulations also contain an additional or non-linear contribution not captured by the present equivalent circuit model. For low external resistance, the voltage amplitudes *A*_{ΔV,±} predict the existence of a single voltage decay mode, namely the fast mode; see Fig. 16(c). The simulations confirm the near-absence of the slow mode, thereby making it difficult to extract reliable model parameters by fitting the voltage relaxation, as reported earlier. The amplitude of the slow mode in the simulations, *A*_{ΔV,1}, agrees with or modestly exceeds the theoretical prediction. Here too, the normalized amplitudes hint at a dependence on the initial electrode charge; the consistency observed across the fitted amplitudes suggests that the relatively high value of *A*_{ΔV,2} at *τ*_{e} = 10^{2} fs is an outlier.

The analytic solution of the equivalent circuit shows an interesting transition for high external resistances. Figure 16(a) shows that with increasing external resistance, the fast relaxation time *τ*_{−} switches from being linear in *R*_{e} to being independent of *R*_{e}, while the slow relaxation time *τ*_{+} switches in the opposite direction; the relevant expressions are provided in Appendix D. Extracting these two relaxation times from the simulations requires longer simulations and enlarging the range of the fit, with the latter making the fit more sensitive to slow processes; at the same time, Figs. 16(b) and 16(c) indicate that the fastest mode is predicted to have a low amplitude. Fitting with a double exponential decay did not pick up the fastest mode; hence, we here resort to triple exponential fits. We continue to denote the slowest and fastest modes with ordinal numbers 1 and 2, respectively, and use the number 3 for the intermediate mode. As Fig. 16(a) shows, the relaxation time of the slowest mode agrees well with theory, while that of the fastest mode shows a more modest agreement. The relaxation times of the intermediate mode lie near the extrapolation of *τ*_{−} from the low *τ*_{e} limit; interestingly, the intermediate mode detected in Fig. 15 lies close to the extrapolation of *τ*_{−} from the high *τ*_{e} limit—exploring these observations in more depth exceeds the objectives of this study. The small theoretical amplitudes of *A*_{Q,−} and *A*_{ΔV,−} in Figs. 16(b) and 16(c), respectively, highlight that in the equivalent circuit only the slowest mode exists for large *R*_{e}, as the resulting slow decay of the electrode charges causes the electrolyte to be in continuous quasi-equilibrium with the electrodes. The sum of the amplitudes of the two slowest modes, $Ay,1+=Ay,1+Ay,3$, recovers the predicted amplitude *A*_{Q,+} for the charge, but for the voltage we again observe a modest overshoot.

*Q*in both models. The mapping does not employ a corresponding one-to-one relation for the coordinate

*ξ*, initially introduced in Sec. II C as the charge in the monolayer adjacent to the electrode and subsequently generalized to a reaction coordinate quantifying the distribution of the ions. In the macroscopic model,

*ξ*has effectively been reduced to an auxiliary coordinate that co-determines the evolution of the charge

*Q*, while its exact meaning remains unspecified. By fitting the instantaneous counterion distribution with Gaussian functions, we extracted the evolving sum charges of the ion layers at both electrodes. Comparing the sum charge in the first layer with the value of

*ξ*in the equivalent circuit does not produce a satisfactory agreement (data not shown), nor does the sum of the two nearest layers yield a match. An agreement is observed, however, between

*ξ*and (minus) the total ionic charge in the right halve of the capacitor,

*θ*denotes the Heaviside step function and

*z*

_{j}is the position of the

*j*th ion relative to the supercapacitor's midplane. Both are shown in Fig. 17 to follow a nearly mono-exponential decay, even though the corresponding charges on the electrodes show a double-exponential decay in the onset to Fig. 7(d). The circuit model explains this observation as due to the small amplitude of the fast relaxation of

*ξ*at the employed external relaxation times. The slow decay time of

*q*

_{R}differs slightly between the simulations at

*τ*

_{e}= 10

^{2}f and

*τ*

_{e}= 10

^{4}f, in agreement with

*τ*

_{+}being nearly independent of

*R*

_{e}for these

*τ*

_{e}and matching the corresponding small difference between the

*τ*

_{Q,1}values in Table I.

## V. CONCLUSIONS

We presented a general method, the constant sum-charge method (CSCM), for molecular dynamics simulations of (dis)connected supercapacitors, accounting for the fluctuating charges on the individual electrode atoms induced by the perpetually changing configuration of the ions in the electrolyte, while the evolution of the sum charge on the electrodes is co-determined by the electric circuit connected to this supercapacitor. CSCM is relatively easy to implement in comparison with alternative constant sum-charge routines^{30–32} because all calculations of Coulombic forces and potentials are outsourced to standard electrostatic solvers.

First, in the simulations of a disconnected supercapacitor with constant sum charges *Q* on the electrodes, the potential difference between the electrodes Δ*V* samples a Gaussian distribution around an average $\Delta VQ$, with a standard deviation that is related to the differential capacitance. The ensemble averages extracted from these simulations, such as the ion distribution and differential capacitance, are a close match to their counterparts obtained by simulations at the constant potential difference Δ*V* for which $Q\Delta V\u2248Q$. Under the latter condition, the fluctuating potential difference in CSCM agrees well with the constant potential difference in the constant potential method (CPM), $\Delta VQ\u2248\Delta V$.

Second, we consider the charging and discharging of a supercapacitor in a closed circuit using CSCM. Here, an equation of motion for the sum charge on the electrode is introduced, allowing current to flow to or from the electrodes. The charge and potential during (dis)charging show double-exponential relaxation, with one relaxation time related to the dynamics of the ions within the capacitor and the other to the current dynamics in the external closed circuit. The slower the latter process, the more energy the capacitor can deliver while discharging—although slowing down obviously reduces the delivered power. A higher external relaxation time also reduces the amount of heat produced in the capacitor while charging. For a capacitor charged in series with a resistor, however, this merely shifts the heat loss to the resistor. The modeled system sees its differential capacitance decrease with potential, causing the energy loss in this charging setup to exceed the energy stored. In conventional simulations of (dis)charging by CPM, the external potential source—rather than the capacitor—is in control over the potential difference between the electrodes; they lack the external relaxation process and its impact on the internal dynamics of the capacitor. Our simulations with CSCM show that it can handle both CPM’s limit of vanishingly short external relaxation times and the continuous equilibrium limit of very long external relaxation times, as well as the entire intermediate range.

A simple phenomenological model with two coupled linear first order differential equations of motion, one for the charge on the electrodes and one for a single coordinate description of the ion distribution, provides a satisfactory agreement with all simulation results reported above. Endowing the volume between the electric double layers with a capacitance, i.e., the *C*_{Qξ} in Figs. 2 and 14, is crucial to recovering the double exponential relaxation observed in the simulations. Physically, this capacitance can be understood as representing interactions between the left and right electrode–electrolyte double layers—including the shoulders of the ion distribution extending into the capacitor’s interior—when these layers are not neutral. This capacitance is, therefore, important in simulations, where the slit width is typically just one order of magnitude larger than the thickness of the electric double layer, while it will be less relevant under experimental conditions with much wider slits. Note that the distances between electric double layers within a porous electrode can be as small as the simulated slit widths; hence, the effect observed here will probably also occur within porous electrodes. The equivalent circuit agrees with the molecular dynamics simulations on the influence of the external resistance on the relaxation times and amplitudes of the double exponential decay. This equivalent circuit approaches its limits when the high charge on the electrode induces the formation of a second layer of counterions; this happens in the current MD model for potentials exceeding ∼1.0 V. While the contribution of a third relaxation process is still small in our simulations, it would be interesting to see what extensions to the model are required when the volume between the flat electrodes contains higher ion concentrations or ionic liquids.

## ACKNOWLEDGMENTS

This work forms part of the Data-driven Science for Smart and Sustainable Energy Research program, with Project No. 16DDS014, by the Netherlands Organization for Scientific Research (NWO), which is funded by the Dutch Ministry of Education, Culture and Science (OCW).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Ranisha S. Sitlapersad**: Conceptualization (supporting); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (supporting); Project administration (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). **Anthony R. Thornton**: Conceptualization (lead); Funding acquisition (lead); Methodology (supporting); Project administration (equal); Resources (equal); Software (supporting); Supervision (equal); Writing – review & editing (equal). **Wouter K. den Otter**: Conceptualization (lead); Formal analysis (supporting); Methodology (lead); Project administration (equal); Resources (equal); Software (supporting); Supervision (equal); Validation (supporting); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: EQUIVALENT CIRCUIT

*Q*

_{Qξ}, and

*C*

_{ξξ}will be denoted as $IQQ\u2032$,

*I*

_{Qξ}, and

*I*

_{ξξ}, respectively. Kirchhoff’s current law then says $IQQ\u2032=IQ\xi +I\xi \xi $. The charges on the right electrodes of the four capacitors will be denoted as $QQQ\u2032$ (twice),

*Q*

_{Qξ}, and

*Q*

_{ξξ}, with the opposite charges residing on their left electrodes. Recall that the charges on the electrodes “inside” the supercapacitor represent the ions of the electrolyte, where the numbers of ions are conserved for a closed capacitor. Equivalently, the total charge on either three-legged network combining one “interior” electrode of $CQQ\u2032$ with one electrode of

*C*

_{Qξ}and one electrode of

*C*

_{ξξ}and one resistor $R\xi \u2032$ is conserved, adding up to zero,

*C*

_{Qξ}and the other combining the two internal capacitors, to yield

*C*

_{ξξ}, but it does not provide additional information. Solving the currents from these equations, using the conservation laws $IQQ\u2032=\u2212Q\u0307QQ\u2032$ and $I\xi \xi =\u2212Q\u0307\xi \xi $ and employing Eq. (A1) to eliminate

*Q*

_{Qξ},

*ξ*=

*Q*

_{ξξ}. The remaining identifications are read off as

*M*

_{Q}= 1/

*R*

_{e}, $kQQ=2/CQQ\u2032$,

*k*

_{Qξ}= 1/

*C*

_{Qξ}, $M\xi =1/(2R\xi \u2032)$, and

*k*

_{ξξ}= 1/

*C*

_{ξξ}.

*Q*

_{Qξ}by the definition of

*ξ*yields

### APPENDIX B: ANALYTICAL RELAXATION

**x**, with matrix

**x**= (

*Q*,

*ξ*)

^{T}. The two relaxation times

*τ*

_{±}of this equation are related to the eigenvalues of the matrix,

*λ*

_{±}, by

*k*’s and

*M*’s, the eigenvalues are negative, with

*λ*

_{+}<

*λ*

_{−}< 0, and the relaxation times are positive, with 0 <

*τ*

_{−}<

*τ*

_{+}. The double-exponential decay of

**x**is then solved by employing a similarity transformation or by inserting the expected double-exponential decay function (with unknown amplitudes) in the differential equation, yielding for the two components of

**x**the relaxation functions

Upon using that for an equilibrated charged capacitor the initial values of *Q*(0) and *ξ*(0) are related by Eq. (30), one readily evaluates the amplitudes of the double exponential decays of *Q*(*t*) and *ξ*(*t*). Insertion of this evolution in the linear transformation of Eq. (33a) produces the corresponding driving force on the electrode charge, i.e., the voltage difference between the electrodes, Δ*V*(*t*), whose decay is, therefore, characterized by the same pair of relaxation times.

The charging process described by Eqs. (39) and (34b) reduces to an ordinary differential equation of the form $y\u0307=My+e$, with vector **y** = (*Q*,*ξ*)^{T}, matrix **M** as given above, and vector $e=(MQVe,0)T$. This equation is readily solved by noting that the transformation **x** = **y** + **M**^{−1}**e** recovers an expression of the above form. The resulting vector **y** starts from zero at *t* = 0 and evolves by a double-exponential converging to **y** = −**M**^{−1}**e** for large times, as in Eq. (45).

### APPENDIX C: ION MOBILITY

*M*

_{ξ}of the ion distribution to the microscopic motions of the constituent ions. Consider the density distributions

*ρ*

_{±}(

*z*) of monovalent ions between two electrodes, −

*H*/2 ≤

*z*≤

*H*/2, as shown in Fig. 6. Following Eq. (50) and assuming symmetry around

*z*= 0, the reaction coordinate follows as

*ρ*

_{±}(0), and the mean velocities of the ions crossing that midplane,

*v*

_{±}(0),

*F*

_{±}, as

*M*

_{ion}denotes the mobility of an ion and

*E*(0) is the electric field at the midplane. Assuming that the charge distribution is uniform parallel to the electrodes, the electric field in the midplane of the capacitor is readily obtained by integrating over the charges to the left and right of the midplane,

*Q*

_{L}= −

*Q*and

*Q*

_{R}=

*Q*being the charges on the left and right electrodes, respectively, and where Eq. (C1) was used in the last step. Combining the above expressions yields

*k*

_{ξξ}≈ 0 for the present system, gives the electrolyte mobility as

*D*

_{ion}= 6.6 × 10

^{−6}cm

^{2}/s, the ion mobility is obtained as

*M*

_{ion}=

*D*

_{ion}/(

*k*

_{B}

*T*). The ion density in the midplane is estimated in Fig. 6 as

*ρ*

_{±}(0) = 1 nm

^{−3}. In combination with the

*k*

_{Qξ}of the equivalent circuit, we find

*M*

_{ξ}= 1.9 × 10

^{−7}C/(V s) or 1.2 × 10

^{−3}

*e*/(V fs), in remarkably good agreement with the value of 8.93 × 10

^{−4}

*e*/(V fs) extracted from the MD simulation. At the same time, the explicit dependence of the mobility

*M*

_{ξ}on

*ρ*

_{±}(0) and its implicit dependence on these densities through

*M*

_{ion}highlight that the actual relaxation process is more complicated than the simplified model of Eq. (34) and Fig. 14. Inserting the approximation for

*C*

_{Qξ}in Eq. (49) gives

*R*

_{ξ}= 1/

*M*

_{ξ}is proportional to this distance.

### APPENDIX D: ANALYTICAL LIMITS

*R*

_{e}→ 0 and

*M*

_{Q}→ ∞, the two relaxation times in Eq. (B3) converge to

*t*= 0, the values of

*Q*(0) and

*ξ*(0) are related by Eq. (30), yielding for the electrode charges,

*Q*(

*t*) follows double exponential decay, while

*ξ*(

*t*) reduces to single exponential decay. For the potential difference then follows, using Eq. (33a),

*R*

_{e}→ ∞ and

*M*

_{Q}→ 0, the two relaxation times converge to

*ξ*is enslaved to

*Q*by Eq. (30).

## REFERENCES

*Nanoscience and Technology: A Collection of Reviews from Nature Journals*

*Handbook of Nanocomposite Supercapacitor Materials I: Characteristics*

**2022**(9),

*Understanding Molecular Simulations. From Algorithms to Applications*

*Computer Simulation of Liquids*

*Introduction to Solid State Physics*

*Computer Simulation Using Particles*