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.

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 applications2,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 Searles32 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 electrolyte20,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.

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.

In this section, we introduce our constant sum-charge method (CSCM) for simulating the charging and discharging dynamics of systems with electrodes, in particular supercapacitors. Following common practice in molecular dynamics simulations,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 Ri and variable charges Qi, while the electrolyte contains n atoms with fixed charges qj and variable positions rj. The electrostatic energy of this system reads as
Ue=12QqTABBTCQq,
(1)
where, for notational convenience, the charges of the electrode atoms and those of the electrolyte atoms have been grouped into column vectors Q and q, respectively. The elements of the symmetric 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,
Up=12Kpi=1NQi2,
(2)
where the positive Kp 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 Kp 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, UQ = Ue + Up, can again be expressed in the above vector-matrix form by introducing the matrix Ap = A + Kp1. 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,
Ψ=UQQ=ApQ+Bq.
(3)
For conducting electrodes, all atoms on the right (left) electrode will have the same potential VR (VL). The vector collecting the potentials of all electrode atoms, therefore, takes the form
Ψ=VRYR+VLYL,
(4)
where the elements of the vector YR (YL) 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
Q=Ap1VRYR+VLYLBq.
(5)
In the constant potential method, the potentials of the two electrodes are set at VR = (V + ΔV/2) and VL = (V − ΔV/2), respectively, where the potential difference ΔV = VRVL is selected by the user and V = (VR + VL)/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, iQi = (YR + YL) · Q = 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.
In an open circuit, the sum charges on the right and left electrodes are conserved independently,
YRQ(t)=QR,
(6a)
YLQ(t)=QL,
(6b)
while the potentials of the two electrodes may vary with time. The challenge now is to find the corresponding electrode potentials, VR(t) and VL(t).31,32 Combining Eqs. (4) through (6) gives the charges on the two electrodes as
QRQL=DVRVLYRAp1BqYLAp1Bq,
(7)
where the four elements of D read as Dαβ=YαAp1Yβ, with α and β taking the values “R” and “L”. Inversion yields the potentials VR and VL for the given electrode charges QR and QL at the momentaneous ion distribution. The charges of the electrode atoms are then solved in CSCM from32,
Q=Ap1GRQR+GLQLHBq,
(8)
where the products of QR and QL with the vectors
GR=D1DLLYRDLRYL,
(9a)
GL=D1DRRYLDRLYR
(9b)
represent the potentials of the electrode atoms in the absence of ions and the projection matrix
H=1GRYRAp1GLYLAp1
(10)
prevents the potential due to the ions, Bq, from altering the sum charges on the electrodes since
YαAp1H=0T.
(11)
For immobile electrode atoms, the products Ap1GR, Ap1GL, and Ap1H are constants and, hence, need to be computed only once. Note that the potential difference between the electrodes now fluctuates with time.29–32 
The charge-related force on the jth ion follows by differentiating the potential energy with respect to the ion’s position,
fjQ=UQrj=UerjQUQQQQrj,
(12)
where the RHS emphasizes that the usual Coulomb force evaluated at the charges Q is extended with a force arising from the position-dependent induced charges of the electrode atoms. The latter term can be rewritten as
UQQQQrj=ΨAp1HBrjq,
(13)
which vanishes identically since ΨAp1H=0T, as follows from combining Eqs. (4) and (11). Hence, after determining the charges on the electrode atoms, the charge-related forces on the ions are readily evaluated using standard routines for electrostatics. Upon coupling this system to a heat bath at temperature T, the resulting model samples the canonical Boltzmann distribution at constant electrode charges,
P(r;QR,QL)dr=1ZeβUQ+UQdr,
(14)
with r denoting the collective coordinates of the electrolyte atoms, Z denoting the normalizing configuration integral, UQ denoting the charge-related energy evaluated at Q, UQ denoting the charge-independent interactions such as the van der Waals terms and covalent contributions, and β = 1/(kBT) with Boltzmann’s constant kB.
The evolving potential difference between the electrodes is readily solved from Eq. (7) as32 
ΔV=gRQR+gLQL+hBq,
(15)
with
gR=DLL+DLR|D|,
(16a)
gL=DRL+DRR|D|,
(16b)
and
h=gRYR+gLYLAp1.
(17)
For immobile electrode atoms, the scalars gR and gL and the vector h are constants and, hence, need to be computed only once.
In the discussion thus far, the sum charges on both electrodes are kept constant, with the distribution of these sum charges over the electrode atoms provided by Eq. (8). The discharging of the supercapacitor is simulated by allowing charge to flow between the electrodes through an external resistance Re, 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
I=dQRdt=dQLdt=ΔVRe.
(18)
This relation is used to update the sum charges QR(t) and QL(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 Ve and a resistance Re. 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 = (ΔVVe)/Re.
We will henceforth limit the discussion to capacitors with an overall neutral electrolyte. For a charge neutral capacitor, the charges on the two electrodes will then be equal in absolute value but opposite in sign,
Q=QR=QL=QRQL2.
(19)
The capacitance C = QV of a supercapacitor will, in general, vary with the potential or charge. Its differential capacitance Cdiff = 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, ⟨ΔVQ.
The charge fluctuations in the CPM ensemble—again assuming, an overall neutral electrolyte—provide a direct route to the capacitance,43,46–48
CdiffCPM=Cdiff,0CPM+βQ2ΔVQΔV2,
(20)
where the bracketed difference measures the variance of the electrode charges due to the electrolyte dynamics and where the capacitance of the empty capacitor reads as
Cdiff,0CPM=Y±SpY±,
(21)
with vector Y± = (YRYL)/2 and the symmetric matrix43,46
Sp=Ap1Ap1YYAp1YAp1Y,
(22)
with vector Y = YR + YL. Note that replacing the vector Y± in Eq. (21) with its projected counterpart,46  Y±Y(Y · Y±)/(Y · Y), does not affect Cdiff,0CPM because the matrix Sp removes homogeneous distributions, SpY = YTSp = 0.
We now derive a similar fluctuation formula in the constant sum charges ensemble, as was also obtained in a recent study by Tee and Searles,32 by focusing on the derivative of the average potential difference,
1CdiffCSCM=dΔVQdQ=gRgL+hBqQQ,
(23)
where the last step follows from Eq. (15). Differentiation of an ensemble average with respect to the electrode charge requires an expression for the corresponding derivative of the average charge-related energy, which is obtained from the Boltzmann distribution of Eq. (14) upon inserting the charge vector of Eq. (8) in the expression for the energy UQ,
UQQQ=GRGLAp1GRGLQHBq+Bq=hG±Q+Bq,
(24)
where in the last step use was made of Eq. (11), of the definitions of Gα, gα, and h, and where G± = GRGL. Differentiating the electrode-induced average potentials on the electrode atoms then gives
BqQQ=βBqQhG±Q+BqQβBqhG±Q+BqQ=βBqQhBqQβBqhBqQ,
(25)
where in the last step we used that h · G±Q is a constant. Insertion in Eq. (23) and noticing that h · Bq represents the only fluctuating contribution to the potential difference ΔV in Eq. (15), it follows that32 
1CdiffCSCM=gRgLβ(ΔV)2QΔVQ2,
(26)
which is the CSCM equivalent to Eq. (20).
An empty capacitor does not show voltage fluctuations in the current approach, nor does it show charge fluctuations in CPM. The capacitance of the empty capacitor does not depend on the voltage in CPM or the charge in CSCM. The above derivations, therefore, offer two expressions for the differential capacitance of an empty capacitor. It can be proven that they agree32 
Cdiff,0CPM=Cdiff,0CSCM=1gRgL=C0,
(27)
by expanding the capacitances in terms of Dαβ=YαAp1Yβ. Both methods are obviously closely related: for a given electrolyte configuration, the charge distribution by CPM for a potential difference ΔV, QCPM = SpVY±Bq), matches that in CSCM for electrode sum charges Qα = Yα · QCPM. 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.
Here, we will derive a simple equivalent circuit for this supercapacitor. We are guided by a simplification of the expected ion distribution; see Fig. 1. In a macroscopic idealization, a supercapacitor consists of two parallel flat electrodes of area 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 − 2d; 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ξ,L=±ξ/(2ε0εrA)êz, 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 êz 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ξ=Q/(ε0εrA)êz, 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ξξ=(ξQ)/(ε0εrA)êz. The total potential energy stored in these electric fields reads as
UE=12ε0εrA2dEQξ2+wEξξ2=12ε0εrA2dQ2+wQξ2.
(28)
There is also free energy stored in the non-uniform distribution of the ions. Using the charge accumulated in the ionic monolayers as a “reaction coordinate” to quantify this non-uniformity, the corresponding free energy is approximated to the lowest order as a quadratic in ξ. We, therefore, model the free energy of the capacitor as
F=12kQQQ2+12kQξ(Qξ)2+12kξξξ2,
(29)
where the three coefficients 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 F/ξ=0 as
ξ=kQξkξξ+kQξQ.
(30)
Back-substitution gives the free energy as
F=kQQkQξ+kξξ+kQξkξξ2kQξ+kξξQ2.
(31)
Hence, the equilibrium capacitance, Ceq, is independent of the charge,
Ceq1=2FQ2=1kQQ1+1kξξ1+kQξ1,
(32)
which corresponds to a capacitor CQQ = 1/kQQ in series with two parallel capacitors Cξξ = 1/kξξ and C = 1/k; 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.
FIG. 1.

Illustration of the idealized supercapacitor consisting of two electrodes with charges ±Q, covered by dense counterion monolayers with charges ∓ξ, and a dilute distribution of ions in the center. Also indicated are the electric field strengths resulting from these charges: E within the electric double layers and Eξξ between the electric double layers, with the arrows indicating the directions of these fields for 0 < ξ < Q.

FIG. 1.

Illustration of the idealized supercapacitor consisting of two electrodes with charges ±Q, covered by dense counterion monolayers with charges ∓ξ, and a dilute distribution of ions in the center. Also indicated are the electric field strengths resulting from these charges: E within the electric double layers and Eξξ between the electric double layers, with the arrows indicating the directions of these fields for 0 < ξ < Q.

Close modal
FIG. 2.

Equivalent circuit to the phenomenological discharging model in Eq. (34). The two outermost capacitors of CQQ=2CQQ=2/kQQ represent the energy stored between the electrodes and ionic monolayers, the capacitor C = 1/k accounts for incomplete screening of the electrode charges by the ionic monolayers, and the capacitor Cξξ = 1/kξξ represents the free energy stored in the non-uniform distribution of the ions. The internal resistances Rξ=Rξ/2=1/(2Mξ) oppose the flow of ions within the supercapacitor, and the external resistance Re opposes the external flow of electric current between the electrodes. The arrows indicate the directions of flow of positive charge carriers during discharging, from the positive electrode on the right to the negative electrode on the left, and from the positive ion layer on the left to the negative ion layer on the right. The charging model has an additional external voltage source. The fits of the simulations with this circuit (see Table I) indicate that the capacitor Cξξ is of little consequence; see also Fig. 14.

FIG. 2.

Equivalent circuit to the phenomenological discharging model in Eq. (34). The two outermost capacitors of CQQ=2CQQ=2/kQQ represent the energy stored between the electrodes and ionic monolayers, the capacitor C = 1/k accounts for incomplete screening of the electrode charges by the ionic monolayers, and the capacitor Cξξ = 1/kξξ represents the free energy stored in the non-uniform distribution of the ions. The internal resistances Rξ=Rξ/2=1/(2Mξ) oppose the flow of ions within the supercapacitor, and the external resistance Re opposes the external flow of electric current between the electrodes. The arrows indicate the directions of flow of positive charge carriers during discharging, from the positive electrode on the right to the negative electrode on the left, and from the positive ion layer on the left to the negative ion layer on the right. The charging model has an additional external voltage source. The fits of the simulations with this circuit (see Table I) indicate that the capacitor Cξξ is of little consequence; see also Fig. 14.

Close modal
The free energy model gives the generalized forces on the charges as
FQ=FQ=kQQQkQξ(Qξ),
(33a)
Fξ=Fξ=kQξ(Qξ)kξξξ.
(33b)
Note that the generalized force experienced by the electrode charge reproduces the potential difference between the electrodes, FQ = −ΔV. If the two electrodes are connected by an external resistance Re, as shown in Fig. 2, and assuming we are in the overdamped regime, the generalized velocities of the charges are linear in the forces,
Q̇=MQ(kQQ+kQξ)QkQξξ,
(34a)
ξ̇=Mξ(kξξ+kQξ)ξkQξQ,
(34b)
where Ohm’s law gives the charge mobility as MQ = 1/Re 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 ξ,
Q=AQ,+et/τ++AQ,et/τ,
(35a)
ξ=Aξ,+et/τ++Aξ,et/τ,
(35b)
with expressions for the amplitudes 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,
Q̇=MQkQQ+kQξkξξkQξ+kξξQ=1ReCeqQ.
(36)
In the opposite limit of a short-circuit, Re = 0, the electrode charge quickly drops to a quasi-equilibrium with the ions, solved from F/Q=0 as
Q=kQξkQQ+kQξξ,
(37)
where ξ, and hence Q, relaxes by a single exponential process,
ξ̇=MξkQξkQQkQQ+kQξ+kξξξ=1RξCintξ,
(38)
where the last step, by analogy with Eq. (36), combines an internal resistance Rξ = 1/Mξ and an effective internal capacitance Cint. This equation reflects that for Re = 0, the equivalent circuit reduces to a resistor Rξ in series with both a capacitor Cξξ and the parallel pair of capacitors CQQ and C. 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.
The charging of a capacitor, connected in series with a voltage source and a resistor, is modeled by adding the source term Ve to the force FQ in the equations of motion. The resulting potential difference Ve − ΔV describes the force driving the current through the external resistor, and multiplication by the charge mobility gives this current as
Q̇=MQ(kQQ+kQξ)QkQξξVe.
(39)
The solution in combination with Eq. (34b) again yields a double exponential relaxation, but this time to an equilibrium state with ΔV = Ve; see  Appendix B.
Simulations are performed in the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS),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 nm2 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 2H. 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,
u(r)=4εσr12σr6+ε,r<21/6σ,0,r>21/6σ,
(40)
with strengths ɛ = kBT at temperature T = 298 K and diameters σ = 0.5 nm. As explained elsewhere,43 the potentials on the electrode atoms due to the ions, Bq, 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,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 εr, multiplying potentials with ε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,
mr̈i=γṙi+fi+fiR(t),
(41)
with 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, fi 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.

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 ⟨ΔVQ = 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 C0 = 27.32 e/V, which agrees very well with the idealized value of CH = ɛ0ɛrA/H = 27.10 e/V for this system; upon normalization by the area of the electrode, the former corresponds to C0/A = 17.41 μF/cm2. 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 

FIG. 3.

Charge vs voltage for the supercapacitor containing 156 ion pairs. In CSCM, the average voltage difference between the electrodes varies with their imposed sum charge, while in CPM, the average electrode charge varies with the imposed potential difference. The solid line represents classical density functional theory calculations on a nearly identical system.53 

FIG. 3.

Charge vs voltage for the supercapacitor containing 156 ion pairs. In CSCM, the average voltage difference between the electrodes varies with their imposed sum charge, while in CPM, the average electrode charge varies with the imposed potential difference. The solid line represents classical density functional theory calculations on a nearly identical system.53 

Close modal
FIG. 4.

Natural logarithm of the probability density distribution p of the potential difference between the electrodes ΔV for the system with 156 ion pairs at electrode charge Q = 40.68 e. The straight line represents a Gaussian distribution with an average of 0.401 V and a standard deviation of 0.025 V.

FIG. 4.

Natural logarithm of the probability density distribution p of the potential difference between the electrodes ΔV for the system with 156 ion pairs at electrode charge Q = 40.68 e. The straight line represents a Gaussian distribution with an average of 0.401 V and a standard deviation of 0.025 V.

Close modal
FIG. 5.

Differential capacitance for the system in Fig. 3, as deduced from the potential fluctuations at a constant surface charge in CSCM and the surface charge fluctuations at constant potential in CPM.53 The solid line represents classical density functional theory calculations on a nearly identical system.53 

FIG. 5.

Differential capacitance for the system in Fig. 3, as deduced from the potential fluctuations at a constant surface charge in CSCM and the surface charge fluctuations at constant potential in CPM.53 The solid line represents classical density functional theory calculations on a nearly identical system.53 

Close modal

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.

FIG. 6.

Density profiles for anions and cations, differing only in the sign of their charge, between two identical electrodes. The cation profiles are mirrored to highlight this symmetry. The main plots present profiles for simulations conserving sum charges on both electrodes (CSCM), while the insets show the differences with simulations maintaining a constant potential difference between the electrodes (CPM), ρ±Δ=ρ±CSCMρ±CPM, for (a) 51 ion pairs at ΔV = 0.2 V vs Q = 16.15 e, (b) 156 ion pairs at ΔV = 0.4 V vs Q = 40.68 e, and (c) 156 ion pairs at ΔV = 2.0 V vs Q = 142.70 e.

FIG. 6.

Density profiles for anions and cations, differing only in the sign of their charge, between two identical electrodes. The cation profiles are mirrored to highlight this symmetry. The main plots present profiles for simulations conserving sum charges on both electrodes (CSCM), while the insets show the differences with simulations maintaining a constant potential difference between the electrodes (CPM), ρ±Δ=ρ±CSCMρ±CPM, for (a) 51 ion pairs at ΔV = 0.2 V vs Q = 16.15 e, (b) 156 ion pairs at ΔV = 0.4 V vs Q = 40.68 e, and (c) 156 ion pairs at ΔV = 2.0 V vs Q = 142.70 e.

Close modal

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 Q0 = 40.68 and 142.70 e, corresponding to ΔV ≈ 0.4 and 2.0 V, respectively. The strength Re 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 = ReC0, with C0 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, Re = 0, implies τe = 0. Removing the potential difference from the charge calculation in Eq. (5), or its charge-neutral equivalent obtained by replacing Ap1 with Sp,43,46 does not instantaneously drop the electrode charges to zero. Rather, the pull by the electrolyte ions binds a residual QR ≈ (2/3)Q0 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.

FIG. 7.

[(a) and (b)] The voltage difference and [(c) and (d)] the electrode charge during discharging of the supercapacitor with 156 ion pairs, from (a) and (c) Q = 40.68 e and ΔV ≈ 0.4 V or (b) and (d) Q = 142.70 e and ΔV ≈ 2.0 V, through an external resistor. Charges are expressed in terms of elementary charges, and voltages are in volts. The legend indicates the values of the external relaxation time τe = ReC0, based on the external resistance Re selected in the CSCM simulations and the capacitance of the ion-free capacitor C0. The instantaneous removal of the potential difference in the CPM simulation is equivalent to a short-circuit, Re = 0. The simulation data (solid lines) are fitted with a double-exponential decay, yielding the fit parameters in Table I. These parameters are next used to determine the capacitances and resistances in the equivalent circuit (see Fig. 2) and, finally, to plot the discharge curves of this circuit (dashed lines). The inset shows the same curves on an extended time scale. Voltages are shown until the first time they cross zero. Each solid line is the average of two independent runs.

FIG. 7.

[(a) and (b)] The voltage difference and [(c) and (d)] the electrode charge during discharging of the supercapacitor with 156 ion pairs, from (a) and (c) Q = 40.68 e and ΔV ≈ 0.4 V or (b) and (d) Q = 142.70 e and ΔV ≈ 2.0 V, through an external resistor. Charges are expressed in terms of elementary charges, and voltages are in volts. The legend indicates the values of the external relaxation time τe = ReC0, based on the external resistance Re selected in the CSCM simulations and the capacitance of the ion-free capacitor C0. The instantaneous removal of the potential difference in the CPM simulation is equivalent to a short-circuit, Re = 0. The simulation data (solid lines) are fitted with a double-exponential decay, yielding the fit parameters in Table I. These parameters are next used to determine the capacitances and resistances in the equivalent circuit (see Fig. 2) and, finally, to plot the discharge curves of this circuit (dashed lines). The inset shows the same curves on an extended time scale. Voltages are shown until the first time they cross zero. Each solid line is the average of two independent runs.

Close modal
When the capacitor discharges through a resistor, as in CSCM, the initial instantaneous drop of the potential and charge are replaced by transient decays. Their relaxation times and residual voltage increase with increasing external resistance, while the residual charge is much less affected; see Fig. 7. Following these transient decreases, the various voltage and charge curves display slow declines with approximately parallel slopes (ΔV) or nearly overlapping lines (Q). To quantify the relaxation processes, the data are fitted with a double-exponential decay,
y(t)=Ay,1et/τy,1+Ay,2et/τy,2,
(42)
where y = Q or ΔV, two amplitudes Ay,1 and Ay,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.
TABLE I.

Amplitudes Ay,i and relaxation times τy,i of the double exponential fit to the electrode charges y = Q (in elementary charges e) and potential difference y = ΔV (in volts) while discharging the capacitor containing 156 ions from two initial charges Q0; see Fig. 7. In CSCM, the value of the external resistance is varied to cover a range of external relaxation times, τe = ReC0, with C0 being the capacitance of the ion-free capacitor; in CPM, the two electrodes discharge by a short-circuit, τe = Re = 0. Simulation results, averaged over two independent runs at the same conditions, are fitted over the first 100 ps, i.e., at least ten times the slowest τe. The three energy coefficients k defined by Eq. (29) and the internal friction Mξ of Eq. (34b) are extracted from these fit parameters by fitting against the theoretical prediction; see Eq. (48) and  Appendix B. Their averaged values over five initial charges Q0 (see Fig. 13) are reported in the bottom line of the table.

Q0 (e)yτe (fs)Ay,1 (y)τy,1 (fs)Ay,2 (y)τy,2 (fs)kQQ (V/e)k (V/e)kξξ (V/e)Mξ [e/(V fs)]
40.68 Q 2.98 × 101 2.24 × 105 ⋯ ⋯     
102 2.92 × 101 2.12 × 105 1.13 × 101 1.12 × 102 8.98 × 10−3 2.36 × 10−2 1.06 × 10−4 7.14 × 10−4 
103 2.92 × 101 2.33 × 105 1.07 × 101 1.28 × 103 7.60 × 10−3 2.06 × 10−2 1.06 × 10−4 7.70 × 10−4 
104 3.24 × 101 2.04 ⋅ 105 8.42 ⋅ 100 8.33 ⋅ 103 1.04 ⋅ 10−2 2.78 ⋅ 10−2 1.06 ⋅ 10−4 7.35 ⋅ 10−4 
ΔV 102 6.22 × 10−4 6.87 × 104 6.12 × 10−1 7.95 × 101 2.03 × 10−2 1.36 × 10−2 1.02 × 10−4 1.50 × 10−3 
103 4.97 × 10−3 2.09 × 105 3.73 × 10−1 1.04 × 103 9.27 × 10−3 2.51 × 10−2 9.88 × 10−5 6.88 × 10−4 
104 5.83 × 10−2 1.81 × 105 3.60 × 10−1 8.54 × 103 1.04 × 10−2 2.61 × 10−2 9.91 × 10−5 8.07 × 10−4 
142.70 Q 8.54 × 101 1.40 × 105 ⋯ ⋯     
102 8.50 × 101 1.44 × 105 5.51 × 101 1.24 × 102 1.14 × 10−2 1.80 × 10−2 1.07 × 10−4 9.82 × 10−4 
103 8.63 × 101 1.44 × 105 5.48 × 101 1.17 × 103 1.22 × 10−2 1.87 × 10−2 1.07 × 10−4 9.39 × 10−4 
104 9.52 × 101 1.63 × 105 4.63 × 101 1.06 × 104 1.27 × 10−2 1.81 × 10−2 1.07 × 10−4 9.07 × 10−4 
ΔV 102 2.67 × 10−3 8.64 × 104 1.75 × 100 1.08 × 102 1.34 × 10−2 1.76 × 10−2 1.07 × 10−4 1.33 × 10−3 
103 2.70 × 10−2 9.02 × 104 1.95 × 100 1.02 × 103 1.48 × 10−2 1.83V10−2 1.07 × 10−4 1.24 × 10−3 
104 2.96 × 10−1 9.34 × 104 1.71 × 100 8.77 × 103 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 
Q0 (e)yτe (fs)Ay,1 (y)τy,1 (fs)Ay,2 (y)τy,2 (fs)kQQ (V/e)k (V/e)kξξ (V/e)Mξ [e/(V fs)]
40.68 Q 2.98 × 101 2.24 × 105 ⋯ ⋯     
102 2.92 × 101 2.12 × 105 1.13 × 101 1.12 × 102 8.98 × 10−3 2.36 × 10−2 1.06 × 10−4 7.14 × 10−4 
103 2.92 × 101 2.33 × 105 1.07 × 101 1.28 × 103 7.60 × 10−3 2.06 × 10−2 1.06 × 10−4 7.70 × 10−4 
104 3.24 × 101 2.04 ⋅ 105 8.42 ⋅ 100 8.33 ⋅ 103 1.04 ⋅ 10−2 2.78 ⋅ 10−2 1.06 ⋅ 10−4 7.35 ⋅ 10−4 
ΔV 102 6.22 × 10−4 6.87 × 104 6.12 × 10−1 7.95 × 101 2.03 × 10−2 1.36 × 10−2 1.02 × 10−4 1.50 × 10−3 
103 4.97 × 10−3 2.09 × 105 3.73 × 10−1 1.04 × 103 9.27 × 10−3 2.51 × 10−2 9.88 × 10−5 6.88 × 10−4 
104 5.83 × 10−2 1.81 × 105 3.60 × 10−1 8.54 × 103 1.04 × 10−2 2.61 × 10−2 9.91 × 10−5 8.07 × 10−4 
142.70 Q 8.54 × 101 1.40 × 105 ⋯ ⋯     
102 8.50 × 101 1.44 × 105 5.51 × 101 1.24 × 102 1.14 × 10−2 1.80 × 10−2 1.07 × 10−4 9.82 × 10−4 
103 8.63 × 101 1.44 × 105 5.48 × 101 1.17 × 103 1.22 × 10−2 1.87 × 10−2 1.07 × 10−4 9.39 × 10−4 
104 9.52 × 101 1.63 × 105 4.63 × 101 1.06 × 104 1.27 × 10−2 1.81 × 10−2 1.07 × 10−4 9.07 × 10−4 
ΔV 102 2.67 × 10−3 8.64 × 104 1.75 × 100 1.08 × 102 1.34 × 10−2 1.76 × 10−2 1.07 × 10−4 1.33 × 10−3 
103 2.70 × 10−2 9.02 × 104 1.95 × 100 1.02 × 103 1.48 × 10−2 1.83V10−2 1.07 × 10−4 1.24 × 10−3 
104 2.96 × 10−1 9.34 × 104 1.71 × 100 8.77 × 103 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 
TABLE II.

Influence of the ion friction parameter γ [see Eq. (41)] on the double exponential relaxation [see Eq. (42)] of the charge and potential for the capacitor with 156 ion pairs. Discharge starts at Q0 = 142.70 e, equivalent to ∼2.0 V; the external relaxation time is fixed at τe = 103 fs.

γ (a.m.u./fs)yAy,1 (y)τy,1 (fs)Ay,2 (y)τy,2 (fs)kQQ (V/e)k (V/e)kξξ (V/e)Mξ [e/(V fs)]
0.226 Q 8.72 × 101 1.02 × 106 5.50 × 101 1.03 × 103 1.36 × 10−2 2.18 × 10−2 1.00 × 10−4 1.16 × 10−4 
ΔV 4.03 × 10−3 1.32 × 105 1.83 × 100 1.05 × 103 1.98 × 10−2 2.57 × 10−3 2.32 × 10−4 2.93 × 10−3 
2.26 Q 8.63 × 101 1.44 × 105 5.48 × 101 1.17 × 103 1.22 × 10−2 1.87 × 10−2 1.07 × 10−4 9.39 × 10−4 
ΔV 2.70 × 10−2 9.02 × 104 1.95 × 100 1.02 × 103 1.48 × 10−2 1.83 × 10−2 1.07 × 10−4 1.24 × 10−3 
22.6 Q 9.61 × 101 1.87 × 104 4.62 × 101 9.74 × 102 1.34 × 10−2 2.07 × 10−2 1.08 × 10−4 7.16 × 10−3 
ΔV 1.97 × 10−1 1.72 × 104 1.80 × 100 9.24 × 102 1.40 × 10−2 2.15 × 0−2 1.08 × 10−4 7.31 × 10−3 
γ (a.m.u./fs)yAy,1 (y)τy,1 (fs)Ay,2 (y)τy,2 (fs)kQQ (V/e)k (V/e)kξξ (V/e)Mξ [e/(V fs)]
0.226 Q 8.72 × 101 1.02 × 106 5.50 × 101 1.03 × 103 1.36 × 10−2 2.18 × 10−2 1.00 × 10−4 1.16 × 10−4 
ΔV 4.03 × 10−3 1.32 × 105 1.83 × 100 1.05 × 103 1.98 × 10−2 2.57 × 10−3 2.32 × 10−4 2.93 × 10−3 
2.26 Q 8.63 × 101 1.44 × 105 5.48 × 101 1.17 × 103 1.22 × 10−2 1.87 × 10−2 1.07 × 10−4 9.39 × 10−4 
ΔV 2.70 × 10−2 9.02 × 104 1.95 × 100 1.02 × 103 1.48 × 10−2 1.83 × 10−2 1.07 × 10−4 1.24 × 10−3 
22.6 Q 9.61 × 101 1.87 × 104 4.62 × 101 9.74 × 102 1.34 × 10−2 2.07 × 10−2 1.08 × 10−4 7.16 × 10−3 
ΔV 1.97 × 10−1 1.72 × 104 1.80 × 100 9.24 × 102 1.40 × 10−2 2.15 × 0−2 1.08 × 10−4 7.31 × 10−3 
The slower the capacitor discharges, the closer the evolving ion distributions will adhere to the equilibrium distributions at the prevalent potential difference, or sum charge. Hence, by increasing the external resistance in the simulations to the limit of very slow discharging, τ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
CdiffdΔVdt=dQdt=ΔVRe,
(43)
where the first step follows from the definition of differential capacitance and the second step employs Ohm’s law. This differential equation for Δ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 = 106 fs (green lines). These discharge curves also allow for a calculation of the total energy delivered up to time t,
Eout(t)=0tΔV(t)I(t)dt,
(44)
executed here using matlab’s cumtrapz routine.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 eV/(ReC0). Conventional theory gives the power as P = (ΔV)2/Re, which under the present conditions translates into P(0) ≈ 109 eV/(ReC0), 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.
FIG. 8.

(a) The potential difference and (b) the total electrode charge when discharging the 156 ion pair supercapacitor from Q0 = 142.70 e for several values of the external resistance as indicated by the external relaxation times τe = ReC0 in the legend. The dashed lines for τe = ∞ apply to the ideal situation of continuous equilibrium, as calculated using Eq. (43) and the equilibrium DFT differential capacitance in Fig. 5. Note that the equilibrium curves are not straight lines due to the non-constant differential capacitance of this supercapacitor; the deviations from these curves in the simulations highlight the impact of non-equilibrium conditions on the discharging behavior.

FIG. 8.

(a) The potential difference and (b) the total electrode charge when discharging the 156 ion pair supercapacitor from Q0 = 142.70 e for several values of the external resistance as indicated by the external relaxation times τe = ReC0 in the legend. The dashed lines for τe = ∞ apply to the ideal situation of continuous equilibrium, as calculated using Eq. (43) and the equilibrium DFT differential capacitance in Fig. 5. Note that the equilibrium curves are not straight lines due to the non-constant differential capacitance of this supercapacitor; the deviations from these curves in the simulations highlight the impact of non-equilibrium conditions on the discharging behavior.

Close modal
FIG. 9.

Cumulative energy delivered by the supercapacitor discharging from Q = 142.70 eV ≈ 2 V) as a function of time. The total delivered energy increases with increasing external relaxation time, τe = ReC0, as indicated in the legend. The dashed curve describes the idealized situation for ion distributions in continuous equilibrium with the electrode charges.

FIG. 9.

Cumulative energy delivered by the supercapacitor discharging from Q = 142.70 eV ≈ 2 V) as a function of time. The total delivered energy increases with increasing external relaxation time, τe = ReC0, as indicated in the legend. The dashed curve describes the idealized situation for ion distributions in continuous equilibrium with the electrode charges.

Close modal
The charging process of the supercapacitor is simulated by connecting it in series with a voltage source at a fixed potential difference of Ve = 0.4 or 2.0 V and a resistance of Re. In CPM, the potential differences between the electrodes jump instantaneously to their final values, Ve; 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 Ve = 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(t)=Ay,11et/τy,1+Ay,21et/τy,2,
(45)
yielding the fit parameters of Table III, where the ordinal number 1 again is assigned to the mode with the largest relaxation time. The relaxation times for different external resistances are similar to their counterparts for the discharging process, with the exception of the fit to the rapidly converging potential for the simulations at 0.4 V. The second mode again follows the external relaxation time, τ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.
FIG. 10.

(a) and (b) The voltage difference and [(c) and (d)] the electrode charge for the system with 156 ion pairs during charging. The initially uncharged capacitor, Q = 0, is connected in series with an external resistance Re and a potential source of [(a) and (c)] Ve = 0.4 V or [(b) and (d)] 2.0 V. The legend indicates the values of the external relaxation time τe = ReC0, based on the external resistance Re selected in the CSCM simulations and the capacitance of the ion-free capacitor C0. The instantaneous imposition of the potential difference in the CPM simulation corresponds to Re = 0. The simulation data (solid lines) are fitted with double-exponential curves, yielding the fit parameters of Table III. The charging process of the equivalent circuit (dashed lines) is based on the capacitances and resistances extracted from the discharging simulations. In the potential plots, the thermal fluctuations in ΔV give rise to a non-zero initial voltage for an uncharged capacitor and some irregular behavior that is most pronounced when the electrode charge and the corresponding potential difference are small, i.e., for short times and under the slowest charging conditions (purple lines).

FIG. 10.

(a) and (b) The voltage difference and [(c) and (d)] the electrode charge for the system with 156 ion pairs during charging. The initially uncharged capacitor, Q = 0, is connected in series with an external resistance Re and a potential source of [(a) and (c)] Ve = 0.4 V or [(b) and (d)] 2.0 V. The legend indicates the values of the external relaxation time τe = ReC0, based on the external resistance Re selected in the CSCM simulations and the capacitance of the ion-free capacitor C0. The instantaneous imposition of the potential difference in the CPM simulation corresponds to Re = 0. The simulation data (solid lines) are fitted with double-exponential curves, yielding the fit parameters of Table III. The charging process of the equivalent circuit (dashed lines) is based on the capacitances and resistances extracted from the discharging simulations. In the potential plots, the thermal fluctuations in ΔV give rise to a non-zero initial voltage for an uncharged capacitor and some irregular behavior that is most pronounced when the electrode charge and the corresponding potential difference are small, i.e., for short times and under the slowest charging conditions (purple lines).

Close modal
TABLE III.

Fit parameters for charging [see Eq. (45)] the supercapacitor with 156 ion pairs by coupling to an external voltage source Ve through an external resistance Re.

Ve (V)yτe (fs)Ay,1 (y)τy,1 (fs)Ay,2 (y)τy,2 (fs)
0.4 Q 2.22 × 101 1.03 × 105 1.08 × 101 
102 1.06 × 102 7.46 × 105 1.09 × 101 1.05 × 102 
103 1.80 × 101 9.33 × 104 1.01 × 101 9.30 × 102 
104 1.61 × 101 4.81 × 104 5.22 × 100 6.96 × 103 
ΔV 102 3.89 × 10−1 1.01 × 102 1.09 × 10−2 2.19 × 100 
103 6.93 × 10−3 1.89 × 104 3.90 × 10−1 8.83 × 102 
104 3.54 × 10−1 9.97 × 103 8.61 × 10−3 3.20 × 10−1 
2.0 Q 8.88 × 101 8.99 × 104 5.53 × 101 
102 8.90 × 101 9.51 × 104 5.47 × 101 1.01 × 102 
103 8.57 × 101 8.22 × 104 5.27 × 101 9.84 × 102 
104 1.08 × 102 1.02 × 105 3.49 × 101 7.91 × 103 
ΔV 102 3.55 × 10−3 8.53 × 104 2.00 × 100 9.83 × 101 
103 3.93 × 10−2 7.45 × 104 1.96 × 100 9.62 × 102 
104 4.01 × 10−1 1.08 × 105 1.61 × 100 7.85 × 103 
Ve (V)yτe (fs)Ay,1 (y)τy,1 (fs)Ay,2 (y)τy,2 (fs)
0.4 Q 2.22 × 101 1.03 × 105 1.08 × 101 
102 1.06 × 102 7.46 × 105 1.09 × 101 1.05 × 102 
103 1.80 × 101 9.33 × 104 1.01 × 101 9.30 × 102 
104 1.61 × 101 4.81 × 104 5.22 × 100 6.96 × 103 
ΔV 102 3.89 × 10−1 1.01 × 102 1.09 × 10−2 2.19 × 100 
103 6.93 × 10−3 1.89 × 104 3.90 × 10−1 8.83 × 102 
104 3.54 × 10−1 9.97 × 103 8.61 × 10−3 3.20 × 10−1 
2.0 Q 8.88 × 101 8.99 × 104 5.53 × 101 
102 8.90 × 101 9.51 × 104 5.47 × 101 1.01 × 102 
103 8.57 × 101 8.22 × 104 5.27 × 101 9.84 × 102 
104 1.08 × 102 1.02 × 105 3.49 × 101 7.91 × 103 
ΔV 102 3.55 × 10−3 8.53 × 104 2.00 × 100 9.83 × 101 
103 3.93 × 10−2 7.45 × 104 1.96 × 100 9.62 × 102 
104 4.01 × 10−1 1.08 × 105 1.61 × 100 7.85 × 103 
The energy provided by the external potential source during charging is absorbed by the capacitor and the external resistor. Both contributions are shown in Fig. 11 for three values of the external resistor, where the cumulative energy transferred to the capacitor up to time 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 Pe(t) = VeI(t), irrespective of the distribution of this power over the resistor and capacitance. The total energy provided to reach equilibrium follows from
Ee=0Pe(t)dt=Ve0I(t)dt=VeQ(Ve),
(46)
where Q(Ve) is the equilibrium charge at the imposed potential difference; see Fig. 3. Interestingly, this energy is a function of Ve only and does not depend on the external resistor. For a capacitor that has reached equilibrium at Ve, the maximum energy it can deliver is also a function of this Ve. It then follows that the energy lost while charging to equilibrium, i.e., the energy that cannot be retrieved upon discharging, depends on Ve but not on Re. 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(Ve) = CVe, 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
Ec(Ve)=0Q(Ve)V(Q)dQ=Ee0VeQ(V)dV,
(47)
where the last step follows from the graphical interpretation of the two integrals as together describing the area VeQ(Ve) = Ee of a rectangle in the QV 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, Ee(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, Ec. Plotting this Ec against Ee(t) (see Fig. 12) highlights that the energy supplied at the start of the simulation, i.e., for ΔVVe, is largely lost; the energy supplied at the end of the simulation, i.e., for ΔVVe, 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.
FIG. 11.

Cumulative energies Ein delivered to the supercapacitor (solid lines) and the external resistor (dashed lines) up to time t during charging of the supercapacitor by an external potential source at Ve = 2.0 V. The legend denotes the external relaxation times.

FIG. 11.

Cumulative energies Ein delivered to the supercapacitor (solid lines) and the external resistor (dashed lines) up to time t during charging of the supercapacitor by an external potential source at Ve = 2.0 V. The legend denotes the external relaxation times.

Close modal
FIG. 12.

Energy stored in the supercapacitor, Ec, as a function of the total energy supplied by a Ve = 2.0 V external source, Ee(t), if the charging circuit is broken at time t and the supercapacitor is next allowed to equilibrate under thermostatted conditions. The legend denotes the external relaxation times during charging. The nine markers for τe = 104 fs are equidistant in time, while for the two smaller τe only the end points after 250 ps are marked. While Ec grows to its limiting value, the slope of the curve steadily increases from zero to unity.

FIG. 12.

Energy stored in the supercapacitor, Ec, as a function of the total energy supplied by a Ve = 2.0 V external source, Ee(t), if the charging circuit is broken at time t and the supercapacitor is next allowed to equilibrate under thermostatted conditions. The legend denotes the external relaxation times during charging. The nine markers for τe = 104 fs are equidistant in time, while for the two smaller τe only the end points after 250 ps are marked. While Ec grows to its limiting value, the slope of the curve steadily increases from zero to unity.

Close modal
In Sec. IV B, we established that the best fits are obtained from the charge during the discharging process. We, therefore, simulate this process for three external relaxation times and five values of the initial electrode charge, corresponding to increasing the initial potential difference from 0.4 to 2.0 V in steps of ≈0.4 V. Converting the extracted amplitudes AQ,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 AQ 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
f=i=i2AQ,iAQ,±(i)12+i=i2τQ,iτ±(i)12
(48)
with respect to the four model parameters entering the analytical expressions for AQ 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 kQQ and k, 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 Q0; 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=2.0×102e/V = 1.2 × 102 μF/cm2, C = 51 e/V = 32 μF/cm2, Cξξ = 0, and Rξ=5.6×102 V fs/e = 3.5 × 105 Ω. The first of these values agrees remarkably well with the capacitance of Cd = ɛ0ɛrA/d = 2.2 × 102 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/2=62μF/cm2, similar to the differential capacitance in the low-voltage range in Fig. 5. While the simple model of Eq. (28) predicts C to be equal to
Cw=ε0εrAwCH,
(49)
where the last holds under the condition Hwd, 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 Dion = 6.56 × 10−6 cm2/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.
FIG. 13.

Fit parameters for the equivalent circuit model extracted by discharging the capacitor from five initial charges Q0 for the three external relaxation times indicated in the legend. The corresponding initial potential differences increase with Q0 from 0.4 to 2.0 V in steps of 0.4 V. Cubic functions are fitted as guides to the eye.

FIG. 13.

Fit parameters for the equivalent circuit model extracted by discharging the capacitor from five initial charges Q0 for the three external relaxation times indicated in the legend. The corresponding initial potential differences increase with Q0 from 0.4 to 2.0 V in steps of 0.4 V. Cubic functions are fitted as guides to the eye.

Close modal
FIG. 14.

Simplification of the equivalent circuit in Fig. 2, based on fits to the simulations. The two outermost capacitors of CQQ=2CQQ=2/kQQ represent the energy stored between the electrodes and ionic monolayers; the capacitor C = 1/k accounts for incomplete screening of the electrode charges by the ionic monolayers. The internal resistance Rξ = 1/Mξ opposes the flow of ions within the supercapacitor, and the external resistance Re opposes the external flow of electric current between the electrodes. In this reduced model, Ceq=CQQ=CQQ/2 and Cint = CQQ + C.

FIG. 14.

Simplification of the equivalent circuit in Fig. 2, based on fits to the simulations. The two outermost capacitors of CQQ=2CQQ=2/kQQ represent the energy stored between the electrodes and ionic monolayers; the capacitor C = 1/k accounts for incomplete screening of the electrode charges by the ionic monolayers. The internal resistance Rξ = 1/Mξ opposes the flow of ions within the supercapacitor, and the external resistance Re opposes the external flow of electric current between the electrodes. In this reduced model, Ceq=CQQ=CQQ/2 and Cint = CQQ + C.

Close modal

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 × 102, 4.62 × 104, and 2.20 × 105 fs, with corresponding amplitudes of 55.0, 23.0, and 64.4 e. The first time again matches the fast external relaxation time of 102 fs, while the last time lies remarkably close to the slow internal time τQ,1 = 2.12 × 105 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 × 105 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 AQ,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.

FIG. 15.

Six discharge simulations at τe = 102 fs, starting from distinct microscopic configurations corresponding to the same macroscopic initial charge Q0 = 142.70 e. The average over the six simulations (black solid line) is fitted with a triple exponential function (black dashed line). The inset illustrates the fast initial decay related to the external resistor; the main figure shows a slow second exponential decay until ∼100 ps and a slower third exponential decay for longer times. Although there is a substantial electrode charge left after 100 ps, the change in the slope is of little practical relevance because this supercapacitor has completely released its energy content by t = 5τe = 0.5 ps; see Fig. 9.

FIG. 15.

Six discharge simulations at τe = 102 fs, starting from distinct microscopic configurations corresponding to the same macroscopic initial charge Q0 = 142.70 e. The average over the six simulations (black solid line) is fitted with a triple exponential function (black dashed line). The inset illustrates the fast initial decay related to the external resistor; the main figure shows a slow second exponential decay until ∼100 ps and a slower third exponential decay for longer times. Although there is a substantial electrode charge left after 100 ps, the change in the slope is of little practical relevance because this supercapacitor has completely released its energy content by t = 5τe = 0.5 ps; see Fig. 9.

Close modal

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 Ve = 0.4 V. The fit parameters extracted from the double exponential discharge curves, for relaxation times from 102 to 104 fs, yield a slow mode that is more pronounced and slowed down relative to its counterpart in Table I, with AQ,1 ≈ 38 e and τQ,1 ≈ 5 × 105 fs, while the fast mode has been reduced to AQ,2 ≈ 3 e with relaxation times about three times shorter than those in Table I. The extracted CQQ for the equivalent circuit remains unchanged, indicating that the electric double layers are not affected by the distance between the electrodes, while the C 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 ≤ 104 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, AQ, are also independent of Re in this limit, as are the measured relaxation amplitudes, AQ,i; see Fig. 16(b). We note that the normalized amplitudes AQ,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 = 102 fs is an outlier.

FIG. 16.

(a) The relaxation times τy,i of the charge (blue markers) and voltage (green markers) in the simulations, along with the theoretical relaxation times τ± of the fitted equivalent circuit (solid lines), plotted against the external relaxation time τe. Triangles pointing up and left denote the slowest i = 1 and fastest i = 2 modes, respectively, when discharging from Q(0) = 40.7 e or ΔV(0) = 0.4 V; triangles pointing down and right represent their counterparts when discharging from Q(0) = 142.7 e or ΔV(0) = 2.0 V. For large τe, a third relaxation mode grows in importance, i = 3 (asterisks), as discussed in the main text. The solid lines, black for the slowest τ+ mode and red for the fastest τ mode, are calculated using Eq. (B3). The dashed lines are extrapolations of the low τe limit, and the dashed–dotted lines those of the high τe limit. The normalized amplitudes of (b) the charge and (c) the voltage plotted against the external relaxation time using the same marking conventions; the crosses denote the sum of the two slowest modes, Ay,p = Ay,1 + Ay,3.

FIG. 16.

(a) The relaxation times τy,i of the charge (blue markers) and voltage (green markers) in the simulations, along with the theoretical relaxation times τ± of the fitted equivalent circuit (solid lines), plotted against the external relaxation time τe. Triangles pointing up and left denote the slowest i = 1 and fastest i = 2 modes, respectively, when discharging from Q(0) = 40.7 e or ΔV(0) = 0.4 V; triangles pointing down and right represent their counterparts when discharging from Q(0) = 142.7 e or ΔV(0) = 2.0 V. For large τe, a third relaxation mode grows in importance, i = 3 (asterisks), as discussed in the main text. The solid lines, black for the slowest τ+ mode and red for the fastest τ mode, are calculated using Eq. (B3). The dashed lines are extrapolations of the low τe limit, and the dashed–dotted lines those of the high τe limit. The normalized amplitudes of (b) the charge and (c) the voltage plotted against the external relaxation time using the same marking conventions; the crosses denote the sum of the two slowest modes, Ay,p = Ay,1 + Ay,3.

Close modal

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 Re to being independent of Re, 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 AQ,− and AΔV,− in Figs. 16(b) and 16(c), respectively, highlight that in the equivalent circuit only the slowest mode exists for large Re, 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 AQ,+ for the charge, but for the voltage we again observe a modest overshoot.

In mapping the MD simulations onto an equivalent circuit, we have equated the sum charges on the electrodes 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,
qR=j=1nqjθzjξ,
(50)
where θ denotes the Heaviside step function and zj is the position of the jth 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 qR differs slightly between the simulations at τe = 102 f and τe = 104 f, in agreement with τ+ being nearly independent of Re for these τe and matching the corresponding small difference between the τQ,1 values in Table I.
FIG. 17.

Evolution of the parameter ξ of the equivalent circuit model (dashed lines) and (minus) the ionic charge qR in the right half of the capacitor in the MD simulations [solid lines, see Eq. (50)] while discharging from Q0 = 142.70 e. The qR curves display the results of two independent simulations at the external relaxation times denoted in the legend.

FIG. 17.

Evolution of the parameter ξ of the equivalent circuit model (dashed lines) and (minus) the ionic charge qR in the right half of the capacitor in the MD simulations [solid lines, see Eq. (50)] while discharging from Q0 = 142.70 e. The qR curves display the results of two independent simulations at the external relaxation times denoted in the legend.

Close modal

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 routines30–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 Δ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ΔVQ. Under the latter condition, the fluctuating potential difference in CSCM agrees well with the constant potential difference in the constant potential method (CPM), ΔVQΔ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 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.

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).

The authors have no conflicts to disclose.

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).

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

We here demonstrate that the equations of motion in Eq. (34) are reproduced by the equivalent circuit in Fig. 2. The currents running from left to right through the capacitors CQQ, Q, and Cξξ will be denoted as IQQ, I, and Iξξ, respectively. Kirchhoff’s current law then says IQQ=IQξ+Iξξ. The charges on the right electrodes of the four capacitors will be denoted as QQQ (twice), 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 with one electrode of C and one electrode of Cξξ and one resistor Rξ is conserved, adding up to zero,
QQQQQξQξξ=0.
(A1)
Kirchhoff’s voltage law can be applied to two loops, one combining the external resistor with the internal capacitor C and the other combining the two internal capacitors, to yield
0=2QQQCQQ+QQξCQξIQQRe,
(A2a)
0=QξξCξξ2IξξRξQQξCQξ,
(A2b)
respectively; there is a third loop, combining the external resistor with the internal capacitor Cξξ, but it does not provide additional information. Solving the currents from these equations, using the conservation laws IQQ=Q̇QQ and Iξξ=Q̇ξξ and employing Eq. (A1) to eliminate Q,
Q̇QQ=1Re2QQQCQQ+QQQQξξCQξ,
(A3a)
Q̇ξξ=12RξQξξCξξQQQQξξCQξ.
(A3b)
One readily recognizes in this pair of equations the structure of Eq. (34), with Q=QQQ and ξ = Qξξ. The remaining identifications are read off as MQ = 1/Re, kQQ=2/CQQ, k = 1/C, Mξ=1/(2Rξ), and kξξ = 1/Cξξ.
For the reduced equivalent circuit in Fig. 14, the above two Kirchhoff’s voltage laws become
0=2QQQCQQ+QQξCQξIQQRe,
(A4a)
0=IξξRξQQξCQξ.
(A4b)
Upon defining the difference between capacitor charges, ξ=QQQQQξ, differentiation gives
ξ̇=Q̇QQQ̇Qξ=IQQ+IQξ=Iξξ,
(A5)
where Kirchhoff’s current law was used in the last step. Solving currents from the above two voltage laws and eliminating Q by the definition of ξ yields
Q̇QQ=1Re2QQQCQQ+QQQξCQξ,
(A6a)
ξ̇=1RξQQQξCQξ.
(A6b)
The identification with the elements of Eq. (34) proceeds as before.
To solve the dynamics of the equivalent circuit of Eq. (34) in a more compact notation, we consider the ordinary differential equation ẋ=Mx for the two-component vector x, with matrix
M=abcd,
(B1)
with diagonal elements (left) and off-diagonal elements (right)
a=MQkQQ+kQξ,d=Mξkξξ+kQξ,b=MQkQξ,c=MξkQξ,
(B2)
and x = (Q,ξ)T. The two relaxation times τ± of this equation are related to the eigenvalues of the matrix, λ±, by
λ±=12a+d±ad2+4bc=1τ±.
(B3)
For positive 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
x1(t)=1D(D+bc)x1(0)b(λ+d)x2(0)eλ+t+1Dbcx1(0)+b(λ+d)x2(0)eλt,
(B4a)
x2(t)=1Dc(λa)x1(0)bcx2(0)eλ+t+1Dc(λa)x1(0)+(D+bc)x2(0)eλt,
(B4b)
where
D=λ+dλabc,
(B5a)
=λ+λλa,
(B5b)
=λ+λdλ+.
(B5c)

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 ẏ=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−1e 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−1e for large times, as in Eq. (45).

The objective here will be to relate the macroscopic 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 ≤ zH/2, as shown in Fig. 6. Following Eq. (50) and assuming symmetry around z = 0, the reaction coordinate follows as
ξ=eA0H/2ρ+ρdz=eAH/20ρ+ρdz,
(C1)
where the last step holds true for an overall neutral electrolyte. The time derivative ξ̇ requires the ion densities in the middle of the capacitor, ρ±(0), and the mean velocities of the ions crossing that midplane, v±(0),
ξ̇=eAρ+(0)v+ρ(0)v.
(C2)
Assuming that the ions are in the overdamped regime, their mean velocities follow from the mean forces acting on these ions, F±, as
v±(0)=MionF±(0)=±eMionE(0),
(C3)
where Mion 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,
E(0)=1ε0εrQLA+e0H/2ρ+ρdzeH/20ρ+ρdzQRA
(C4)
=2ε0εrAξQ,
(C5)
with QL = −Q and QR = 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
ξ̇=2e2ε0εrρ+(0)+ρ(0)MionξQ.
(C6)
Comparison to Eq. (34b), with kξξ ≈ 0 for the present system, gives the electrolyte mobility as
Mξ=2e2ε0εrρ+(0)+ρ(0)MionkQξ.
(C7)
Using the ionic self-diffusion coefficient in the bulk electrolyte of Dion = 6.6 × 10−6 cm2/s, the ion mobility is obtained as Mion = Dion/(kBT). The ion density in the midplane is estimated in Fig. 6 as ρ±(0) = 1 nm−3. In combination with the k 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 Mion highlight that the actual relaxation process is more complicated than the simplified model of Eq. (34) and Fig. 14. Inserting the approximation for C in Eq. (49) gives
Mξ=2e2ρ+(0)+ρ(0)AHMion.
(C8)
If the ion density in the center of the supercapacitance is kept constant while increasing the distance between the electrodes, as is the case in Sec. IV C, the internal resistance Rξ = 1/Mξ is proportional to this distance.
The limiting behaviors of Eq. (34) for discharging through small and large external resistances are obtained by Taylor expansion of the solutions in  Appendix B. For a small external resistance, Re → 0 and MQ → ∞, the two relaxation times in Eq. (B3) converge to
τ+=RξCint,
(D1a)
τ=ReCQQ1+CQξ11.
(D1b)
This limit applies when
ReCintCQQ1+CQξ1Rξ.
(D2)
The amplitudes are solved from Eq. (B4), where for a capacitor in equilibrium at time t = 0, the values of Q(0) and ξ(0) are related by Eq. (30), yielding for the electrode charges,
AQ,+=kQξ2(kQQ+kQξ)(kξξ+kQξ)Q(0),
(D3a)
AQ,=Q(0)AQ,+,
(D3b)
and for the internal charge distribution,
Aξ,+=kQξkQξ+kξξQ(0),Aξ,=0.
(D4)
Hence, Q(t) follows double exponential decay, while ξ(t) reduces to single exponential decay. For the potential difference then follows, using Eq. (33a),
AΔV,+=0,AΔV,=kQQkξξkξξ+kQξQ(0),
(D5)
and, therefore, a single exponential decay.
In the limit of large external resistance, Re → ∞ and MQ → 0, the two relaxation times converge to
τ+=ReCeq,
(D6a)
τ=RξCQξ1+Cξξ11.
(D6b)
This limit applies when
ReCeq1CQξ1+Cξξ11Rξ.
(D7)
Solving the corresponding amplitudes yields single-exponential decay for the electrode charges, internal charge distribution, and potential difference,
AQ,+=Q(0),AQ,=0,
(D8a)
Aξ,+=kQξkQξ+kξξQ(0),Aξ,=0,
(D8b)
AΔV,+=kQQkξξkQξ+kξξQ(0),AΔV,=0,
(D8c)
as expected when ξ is enslaved to Q by Eq. (30).
1.
P.
Simon
and
Y.
Gogotsi
, “
Materials for electrochemical capacitors
,” in
Nanoscience and Technology: A Collection of Reviews from Nature Journals
(
World Scientific
,
2010
), pp.
320
329
.
2.
S.
Banerjee
,
B.
De
,
P.
Sinha
,
J.
Cherusseri
, and
K. K.
Kar
, “
Applications of supercapacitors
,” in
Handbook of Nanocomposite Supercapacitor Materials I: Characteristics
, edited by
K. K.
Kar
(
Springer International Publishing
,
Cham, Switzerland
,
2020
), pp.
341
350
.
3.
T.
Younas
,
K. M. U.
Rehman
,
M. T.
Khan
, and
T.
Inayat
, “
Supercapacitors and its enactment for renewable energy resources
,”
3C Tecnología
2022(9),
65
95
.
4.
A.
Brandt
,
S.
Pohlmann
,
A.
Varzi
,
A.
Balducci
, and
S.
Passerini
, “
Ionic liquids in supercapacitors
,”
MRS Bull.
38
,
554
559
(
2013
).
5.
R.
Burt
,
G.
Birkett
, and
X. S.
Zhao
, “
A review of molecular modelling of electric double layer capacitors
,”
Phys. Chem. Chem. Phys.
16
,
6519
6538
(
2014
).
6.
A. C.
Forse
,
C.
Merlet
,
J. M.
Griffin
, and
C. P.
Grey
, “
New perspectives on the charging mechanisms of supercapacitors
,”
J. Am. Chem. Soc.
138
,
5731
5744
(
2016
).
7.
Z.
Bo
,
C.
Li
,
H.
Yang
,
K.
Ostrikov
,
J.
Yan
, and
K.
Cen
, “
Design of supercapacitor electrodes using molecular dynamics simulations
,”
Nano-Micro Lett.
10
,
33
(
2018
).
8.
K.
Xu
,
H.
Shao
,
Z.
Lin
,
C.
Merlet
,
G.
Feng
,
J.
Zhu
, and
P.
Simon
, “
Computational insights into charge storage mechanisms of supercapacitors
,”
Energy Environ. Mater.
3
,
235
246
(
2020
).
9.
L.
Scalfi
,
M.
Salanne
, and
B.
Rotenberg
,
Annu. Rev. Phys. Chem.
72
,
189
212
(
2021
).
10.
G.
Jeanmairet
,
B.
Rotenberg
, and
M.
Salanne
, “
Microscopic simulations of electrochemical double-layer capacitors
,”
Chem. Rev.
122
,
10860
10898
(
2022
).
11.
J. I.
Siepmann
and
M.
Sprik
, “
Influence of surface topology and electrostatic potential on water/electrode systems
,”
J. Chem. Phys.
102
,
511
524
(
1995
).
12.
S. K.
Reed
,
O. J.
Lanning
, and
P. A.
Madden
, “
Electrochemical interface between an ionic liquid and a model metallic electrode
,”
J. Chem. Phys.
126
,
084704
(
2007
).
13.
Z.
Wang
,
Y.
Yang
,
D. L.
Olmsted
,
M.
Asta
, and
B. B.
Laird
, “
Evaluation of the constant potential method in simulating electric double-layer capacitors
,”
J. Chem. Phys.
141
,
184102
(
2014
).
14.
A.
Coretti
,
C.
Bacon
,
R.
Berthin
,
A.
Serva
,
L.
Scalfi
,
I.
Chubak
,
K.
Goloviznina
,
M.
Haefele
,
A.
Marin-Laflèche
,
B.
Rotenberg
,
S.
Bonella
, and
M.
Salanne
, “
MetalWalls: Simulating electrochemical interfaces between polarizable electrolytes and metallic electrodes
,”
J. Chem. Phys.
157
,
184801
(
2022
).
15.
L. J. V.
Ahrens-Iwers
and
R. H.
Meißner
, “
Constant potential simulations on a mesh
,”
J. Chem. Phys.
155
,
104104
(
2021
).
16.
S. R.
Tee
and
D. J.
Searles
, “
Fully periodic, computationally efficient constant potential molecular dynamics simulations of ionic liquid supercapacitors
,”
J. Chem. Phys.
156
,
184101
(
2022
).
17.
S.
Tyagi
,
M.
Süzen
,
M.
Sega
,
M.
Barbosa
,
S. S.
Kantorovich
, and
C.
Holm
, “
An iterative, fast, linear-scaling method for computing induced charges on arbitrary dielectric boundaries
,”
J. Chem. Phys.
132
,
154112
(
2010
).
18.
I. L.
Geada
,
H.
Ramezani-Dakhel
,
T.
Jamil
,
M.
Sulpizi
, and
H.
Heinz
,
Nat. Commun.
9
,
716
(
2018
).
19.
S.
Kondrat
,
G.
Feng
,
F.
Bresme
,
M.
Urbakh
, and
A. A.
Kornyshev
, “
Theory and simulations of ionic liquids in nanoconfinement
,”
Chem. Rev.
123
,
6668
6715
(
2023
).
20.
C.
Péan
,
C.
Merlet
,
B.
Rotenberg
,
P. A.
Madden
,
P.-L.
Taberna
,
B.
Daffos
,
M.
Salanne
, and
P.
Simon
, “
On the dynamics of charging in nanoporous carbon-based supercapacitors
,”
ACS Nano
8
,
1576
1583
(
2014
).
21.
A.
A Lee
,
S.
Kondrat
,
G.
Oshanin
, and
A.
A Kornyshev
, “
Charging dynamics of supercapacitors with narrow cylindrical nanopores
,”
Nanotechnology
25
,
315401
(
2014
).
22.
Y.
He
,
J.
Huang
,
B. G.
Sumpter
,
A. A.
Kornyshev
, and
R.
Qiao
, “
Dynamic charge storage in ionic liquids-filled nanopores: Insight from a computational cyclic voltammetry study
,”
J. Phys. Chem. Lett.
6
,
22
30
(
2015
).
23.
K.
Breitsprecher
,
M.
Abele
,
S.
Kondrat
, and
C.
Holm
, “
The effect of finite pore length on ion structure and charging
,”
J. Chem. Phys.
147
,
104708
(
2017
).
24.
K.
Breitsprecher
,
C.
Holm
, and
S.
Kondrat
, “
Charge me slowly, I am in a hurry: Optimizing charge–discharge cycles in nanoporous supercapacitors
,”
ACS Nano
12
,
9733
9741
(
2018
).
25.
C.
Péan
,
B.
Rotenberg
,
P.
Simon
, and
M.
Salanne
, “
Multi-scale modelling of supercapacitors: From molecular simulations to a transmission line model
,”
J. Power Sources
326
,
680
685
(
2016
).
26.
K.
Breitsprecher
,
M.
Janssen
,
P.
Srimuk
,
B. L.
Mehdi
,
V.
Presser
,
C.
Holm
, and
S.
Kondrat
, “
How to speed up ion transport in nanopores
,”
Nat. Commun.
11
,
6085
(
2020
).
27.
A. M.
Sampaio
,
G. F. L.
Pereira
,
M.
Salanne
, and
L. J. A.
Siqueira
, “
Comparing the performance of sulfonium and phosphonium ionic liquids as electrolytes for supercapacitors by molecular dynamics simulations
,”
Electrochim. Acta
364
,
137181
(
2020
).
28.
Z.
Gan
,
Y.
Wang
,
M.
Wang
,
E.
Gao
,
F.
Huo
,
W.
Ding
,
H.
He
, and
S.
Zhang
, “
Ionophobic nanopores enhancing the capacitance and charging dynamics in supercapacitors with ionic liquids
,”
J. Mater. Chem. A
9
,
15985
15992
(
2021
).
29.
T.
Dufils
,
G.
Jeanmairet
,
B.
Rotenberg
,
M.
Sprik
, and
M.
Salanne
, “
Simulating electrochemical systems by combining the finite field method with a constant potential electrode
,”
Phys. Rev. Lett.
123
,
195501
(
2019
).
30.
T.
Dufils
,
M.
Sprik
, and
M.
Salanne
, “
Computational amperometry of nanoscale capacitors in molecular simulations
,”
J. Phys. Chem. Lett.
12
,
4357
4361
(
2021
).
31.
L.
Zeng
,
T.
Wu
,
T.
Ye
,
T.
Mo
,
R.
Qiao
, and
G.
Feng
, “
Modeling galvanostatic charge–discharge of nanoporous supercapacitors
,”
Nat. Comput. Sci.
1
,
725
731
(
2021
).
32.
S. R.
Tee
and
D. J.
Searles
, “
Constant potential and constrained charge ensembles for simulations of conductive electrodes
,”
J. Chem. Theory Comput.
19
,
2758
2768
(
2023
).
33.
L. J. V.
Ahrens-Iwers
,
M.
Janssen
,
S. R.
Tee
, and
R. H.
Meißner
, “
ELECTRODE: An electrochemistry package for atomistic simulations
,”
J. Chem. Phys.
157
,
084801
(
2022
).
34.
J. I.
Hidalgo-Reyes
,
J. F.
Gómez-Aguilar
,
R. F.
Escobar-Jiménez
,
V.
Alvarado-Martínez
, and
M.
López-López
, “
Classical and fractional-order modeling of equivalent electrical circuits for supercapacitors and batteries, energy management strategies for hybrid systems and methods for the state of charge estimation: A state of the art review
,”
Microelectron. J.
85
,
109
128
(
2019
).
35.
M. Z.
Bazant
,
K.
Thornton
, and
A.
Ajdari
, “
Diffuse-charge dynamics in electrochemical systems
,”
Phys. Rev. E
70
,
021506
(
2004
).
36.
D.
Frenkel
and
B.
Smit
,
Understanding Molecular Simulations. From Algorithms to Applications
, 2nd ed. (
Academic Press
,
San Diego, CA
,
2002
).
37.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
, 2nd ed. (
Oxford University Press
,
Oxford
,
2017
).
38.
C.
Kittel
,
Introduction to Solid State Physics
(
John Wiley & Sons, Inc.
,
New York
,
1986
).
39.
R. W.
Hockney
and
J. W.
Eastwood
,
Computer Simulation Using Particles
(
Adam Hilger
,
Bristol
,
1989
).
40.
A. Y.
Toukmaji
and
J. A.
Board
, Jr.
, “
Ewald summation techniques in perspective: A survey
,”
Comput. Phys. Commun.
95
,
73
92
(
1996
).
41.
H.
Nakano
and
H.
Sato
, “
A chemical potential equalization approach to constant potential polarizable electrodes for electrochemical-cell simulations
,”
J. Chem. Phys.
151
,
164123
(
2019
).
42.
L.
Scalfi
,
T.
Dufils
,
K. G.
Reeves
,
B.
Rotenberg
, and
M.
Salanne
,
J. Chem. Phys.
153
,
1174704
(
2020
).
43.
R.
Sitlapersad
,
A. R.
Thornton
, and
W. K.
den Otter
, “A simple efficient algorithm for molecular simulations of constant potential electrodes,”
J. Chem. Phys.
160
, 034107 (2024).
44.
A. K.
Rappe
and
W. A.
Goddard
III
, “
Charge equilibration for molecular dynamics simulations
,”
J. Phys. Chem.
95
,
3358
3363
(
1991
).
45.
M.
Buraschi
,
S.
Sansotta
, and
D.
Zahn
, “
Polarization effects in dynamic interfaces of platinum electrodes and ionic liquid phases: A molecular dynamics study
,”
J. Phys. Chem. C
124
,
2002
2007
(
2020
).
46.
L.
Scalfi
,
D. T.
Limmer
,
A.
Coretti
,
S.
Bonella
,
P. A.
Madden
,
M.
Salanne
, and
B.
Rotenberg
, “
Charge fluctuations from molecular simulations in the constant-potential ensemble
,”
Phys. Chem. Chem. Phys.
22
,
10480
10489
(
2020
).
47.
D. T.
Limmer
,
C.
Merlet
,
M.
Salanne
,
D.
Chandler
,
P. A.
Madden
,
R.
van Roij
, and
B.
Rotenberg
, “
Charge fluctuations in nanoscale capacitors
,”
Phys. Rev. Lett.
111
,
106102
(
2013
).
48.
J. B.
Haskins
and
J. W.
Lawson
, “
Evaluation of molecular dynamics simulation methods for ionic liquid electric double layers
,”
J. Chem. Phys.
144
,
184707
(
2016
).
49.
A. P.
Thompson
,
H. M.
Aktulga
,
R.
Berger
,
D. S.
Bolintineanu
,
W. M.
Brown
,
P. S.
Crozier
,
P. J.
in't Veld
,
A.
Kohlmeyer
,
S. G.
Moore
,
T. D.
Nguyen
,
R.
Shan
,
M. J.
Stevens
,
J.
Tranchida
,
C.
Trott
,
S. J.
Plimpton
, and
S. J.
Plimpton
, “
LAMMPS—A flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales
,”
Comput. Phys. Commun.
271
,
108171
(
2022
).
50.
Z.
Wang
,
Y.
Yang
,
D. L.
Olmsted
,
M.
Asta
, and
B. B.
Laird
, https://github.com/zhenxingwang/lammps-conp,
2014
.
51.
R. S.
Sitlapersad
,
A. R.
Thornton
, and
W. K.
den Otter
, https://github.com/r-s-sitlapersad/lammps,
2023
.
52.
I.-C.
Yeh
and
M. L.
Berkowitz
, “
Ewald summation for systems with slab geometry
,”
J. Chem. Phys.
111
,
3155
3162
(
1999
).
53.
P.
Cats
,
R. S.
Sitlapersad
,
W. K.
den Otter
,
A. R.
Thornton
, and
R.
van Roij
, “
Capacitance and structure of electric double layers: Comparing Brownian dynamics and classical density functional theory
,”
J. Solution Chem.
51
,
296
319
(
2022
).
54.
MATLAB version: 9.13.0 (R2022b),
The MathWorks, Inc.
,
Natick, MA, USA
.
55.
K.
Mita
and
M.
Boufaida
, “
Ideal capacitor circuits and energy conservation
,”
Am. J. Phys.
67
,
737
739
(
1999
);
K.
Mita
and
M.
Boufaida
,
Am. J. Phys.
68
,
578
(
2000
).