Several modern technologies for energy storage and conversion are based on the screening of electric charge on the surface of porous electrodes by ions in an adjacent electrolyte. This so-called electric double layer (EDL) exhibits an intricate interplay with the electrolyte’s temperature that was the focus of several recent studies. In one of them, Janssen et al. [Phys. Rev. Lett. 119, 166002 (2017)] experimentally determined the ratio Qrev/Wel of reversible heat flowing into a supercapacitor during an isothermal charging process and the electric work applied therein. To rationalize that data, here, we determine Qrev/Wel within different models of the EDL using theoretical approaches such as density functional theory (DFT) as well as molecular dynamics simulations. Applying mainly the restricted primitive model, we find quantitative support for a speculation of Janssen et al. that steric ion interactions are key to the ratio Qrev/Wel. Here, we identified the entropic contribution of certain DFT functionals, which grants direct access to the reversible heat. We further demonstrate how Qrev/Wel changes when calculated in different thermodynamic ensembles and processes. We show that the experiments of Janssen et al. are explained best by a charging process at fixed bulk density or in a “semi-canonical” system. Finally, we find that Qrev/Wel significantly depends on parameters such as pore and ion size, salt concentration, and valencies of the cations and anions of the electrolyte. Our findings can guide further heat production measurements and can be applied in studies on, for instance, nervous conduction, where reversible heat is a key element.

In a recent experiment,1 the reversible heat flowing into and the electric work applied to a supercapacitor during isothermal charging were measured. This experiment is one of several recent studies, both experimental1–4 and theoretical,5–10 on the intricate interplay between the electrolyte’s temperature and the properties of the electric double layer (EDL). Until now, however, no comparison has been made between the experimental findings of Ref. 1 and theoretical predictions from sophisticated EDL models.

Helmholtz proposed the EDL to be a system where two layers of opposite charges are facing and, thus, screening each other.11 Usually, systems are considered where mobile ions physically screen electric charge, for instance, on a solid electrode’s surface,12 on colloidal particles,13,14 in (biological) ion channels of the plasma membrane,15,16 and near macromolecules such as DNA.17 The resulting diffuse double layer for point-like ions was first described by Gouy and Chapman around 1910 within Poisson–Boltzmann theory,18–20 a framework even nowadays still applied frequently to study electrolyte systems. This simple picture of point charges is refined in more sophisticated models that account for finite ionic volume, where the latter can be important for the microscopic structure of EDLs in narrow geometries and crowded environments. For instance, the finite volume of ions can be crucial for the description of colloidal interactions, capacitances, and understanding certain aspects of screening.21–23 

In modern technologies, EDLs also form the basis for the aforementioned supercapacitors, which can store much more energy than conventional capacitors and can deliver much higher power than batteries.24–26 For practical applications wherein these devices are charged and discharged, it is important to know how the electrolyte temperature can be kept low because increased temperatures are the cause of faster degradation of components. Interestingly, these EDL systems can further be employed to desalinate solutes27 and to harvest energy because concentration28 and temperature8,29,30 can change their capacitance. Accordingly, vast amounts of studies on EDL systems exist and are still performed, but the interplay between ions and the electrolyte’s temperature is still rather unexplored, despite being promising for optimization and new concepts.

Measurements of the temperature of a supercapacitor in operation showed that it heated during charging and cooled during discharging,2 showing an overall trend to warm up during cycling. Such a warm-up is expected as ionic currents in a resistive fluid dissipate Joule heat. The cooling, however, can be understood from an analogy to the adiabatic decompression of an ideal gas: During discharging, ions leave the EDL and their entropy increases. In an isolated supercapacitor, this increase must be balanced by an entropy decrease of the electrolyte, accomplished through a lowering of the electrolyte’s temperature. The opposite happens during charging and causes heating additional to Joule heat. Moving beyond the above ideal-gas analogy, a thermodynamic identity was derived for the temperature rise upon adiabatic EDL formation.29 Predictions of this identity coincided with numerical solutions of the electrokinetic equations for a slow charging process.7 In the latter nonequilibrium framework, reversible heating and irreversible heating are both captured by the heat production term IE in the heat equation5,7,31 with I being the ionic current density and E being the local electric field. While the ionic current density aligns with the local electric field IE in bulk electrolytes, leading to a strictly positive Joule-heating term I2, conversely IE<0 is possible in the EDL when the gradient in electrochemical potential anti-aligns with E, leading to reversible local cooling.5,7 Similar cooling has been observed near an ion-exchange membrane.32,33 As Joule heating is mainly a bulk phenomenon, while reversible heating happens only in the nanometer-wide EDL, a capacitor with a large surface-to-volume ratio is needed to notice appreciable reversible temperature variations. Advanced “microcalorimetry” measurements near flat electrodes, however, can detect much smaller temperature variations.3 

With a setup as sketched in Fig. 1(a), the authors of Ref. 1 studied the temperature of and charge on nanoporous carbon electrodes subject to a suddenly applied potential change. From the difference QtotQirr between the total and irreversible (Joule) heat, for which they had independent measurements, they determined the reversible heat Qrev, i.e., the temporal and spatial integral of the above heat production for slow charging [cf. Eq. (12)]. Moreover, from the system’s capacitance, they determined the electric work Wel during isothermal charging. With Qrev and Wel at hand, the authors of Ref. 1 claimed experimental access to the ratio ΔΩent/ΔΩ (reproduced here in Fig. 4) with ΔΩ being the change of the total grand potential during charging and ΔΩent being its entropic part. Their identification Qrev/Wel=ΔΩent/ΔΩ relied on the identity Qrev=ΔΩent [cf. Eq. (6)], proposed by Overbeek on thermodynamic grounds,34 and on ΔΩ = Wel, which holds for isothermal charging. The linear scaling of both Qrev and Wel with the surface area of the electrodes drops in their ratio, making Qrev/Wel a quantity that can be conveniently compared with theoretical model predictions. In fact, ΔΩent and ΔΩ had been studied before within Poisson–Boltzmann theory34 and extensions thereof accounting for finite-size ions.35,36 In particular, the authors of Ref. 36 found that the Carnahan-Starling bulk chemical potential radically altered ΔΩent/ΔΩ at large voltages: While Poisson–Boltzmann predicts ΔΩent/ΔΩ → 1, their more sophisticated theory suggested ΔΩent/ΔΩ → 0 instead. For this reason, the authors of Ref. 1 speculated about the importance of ionic steric interactions to their measurement of ΔΩent/ΔΩ ≈ 0.25 at ΔΨ = 1 V.

FIG. 1.

Sketch (a) of the experimental setup of Ref. 1. In the experiments, two nanoporous carbon electrodes that carry charges ±Q due to a potential difference ΔΨ are in contact with an aqueous sodium chloride solution and submerged in a thermostatic bath with respect to which its temperature is measured. Within our theoretical description, the pores of one electrode are modeled (b) by two planar walls of area A with equal surface charge density σ. The electrolyte between the walls is modeled by mobile ionic charges.

FIG. 1.

Sketch (a) of the experimental setup of Ref. 1. In the experiments, two nanoporous carbon electrodes that carry charges ±Q due to a potential difference ΔΨ are in contact with an aqueous sodium chloride solution and submerged in a thermostatic bath with respect to which its temperature is measured. Within our theoretical description, the pores of one electrode are modeled (b) by two planar walls of area A with equal surface charge density σ. The electrolyte between the walls is modeled by mobile ionic charges.

Close modal

In this study, we determine Qrev/Wel within the restricted primitive model (RPM), where ions are described as charged hard spheres and the solvent is represented by a homogeneous background. The RPM is easy to simulate and is well described in classical density functional theory (DFT),22,37 a version of the famous quantum DFT adopted to classical systems.38–40 In our study, we apply several theoretical approaches of different sophistication: the (modified) Poisson–Boltzmann theory of Refs. 34 and 35 without and with a Stern layer, a density functional theory with a very accurate description of the hard-sphere interaction,22 and molecular dynamics (MD) simulations.41 For our DFT approaches, we verify Overbeek’s identity Qrev=ΔΩent. Moreover, we discuss the importance of the choice of the thermodynamic process and the corresponding ensemble, which we demonstrate using the analytical Gouy–Chapman solution to the Poisson–Boltzmann equations. Finally, anticipating future experiments with other electrode–electrolyte combinations than used in Ref. 1, we study different pore sizes, ionic radii, valencies, and bulk ion concentrations.

Thermodynamics allows us to draw general conclusions for our model system without using microscopic details. We consider the setup sketched in Fig. 1(a), where two porous electrodes are immersed in an electrolyte held at temperature T. The pores of each electrode have a certain fixed geometry during charging processes. Our system contains an electrolyte with N+ positive ions, N negative ions, and Ns neutral solvent particles. Conjugated to these particle numbers are the intensive chemical potentials μi with i ∈ { +, −, s} that reflect the connection of the pore volume to an (infinitely) large reservoir that the electrodes are immersed in.

Upon connecting the electrodes to a battery that supplies a potential difference ΔΨ, the electrodes acquire electric surface charges Q and −Q. The first law of thermodynamics for this system relates the change of internal energy U, heat Q transferred to the system, and thermodynamic work W done to the system and reads

dU=δQ+δW.
(1)

The electric work performed during charging is given by

Wel=0QΔΨ dQ.
(2)

Using Legendre transforms, we obtain the free energy F(T, Q, Ni) = U(S, Q, Ni) − TS and the grand potential Ω(T, Q, μi) = F(T, Q, Ni) − ∑iμiNi, where the entropy S enters.

According to Eq. (1), a process in, say, a grand canonical system wherein the surface charge on the positive electrode changes from Q1 to Q2 must cause a (reversible) heat flow into the capacitor that reads

Qrev=TΔS.
(3a)

Here, the entropy difference ΔS is given by

ΔS=ΩTμi,Q2+ΩTμi,Q1.
(3b)

Using Eq. (3a) and the Maxwell relation (∂S/∂Q)T = −(ΔΨ/∂T)Q, the heat during the isothermal charging process follows as [see also Eq. (S7) of Ref. 1]

QrevI=TQ1Q2SQT,μidQ
(4a)
=TQ1Q2ΔΨTQ,μidQ.
(4b)

In Eq. (4a), we introduce a superscript I for the reversible heat to distinguish it from the later result in Eq. (12), which uses microscopic information and is valid only in a canonical system. Note that Eqs. (3b), (4a), and (4b) also hold in a canonical ensemble if the grand potential Ω is replaced by the free energy F and the derivatives are taken at constant Ni instead of μi.

For a general process where all the natural variables of a thermodynamic potential except for the surface charge are kept constant, the electric work corresponds to the change in the thermodynamic potential. For instance, in a grand canonical system, the above charging process, where the chemical potentials are kept constant, results in the electric work,

Wel=ΔΩ.
(5)

Likewise, the reversible heat can be expressed as the entropic contribution −TS to the thermodynamic potential. Here, we already mention that both F and Ω contain this term −TS. Thus, knowing expressions of the thermodynamic potentials and being able to identify the contribution −TS would allow us to directly read off the change of entropy and, hence, the reversible heat. Later, we will see that this becomes handy in the framework of DFT. From here on, for convenience, we will use the uncharged electrodes as the reference state, and hence, the change in the grand potential ΔΩ and its entropic contribution ΔΩent vanishes for ΔΨ = 0. Regarding a grand canonical system, for example, this means that the reversible heat during charging satisfies

Qrev=ΔΩent.
(6)

Again, Eqs. (5) and (6) also hold for a canonical system if Ω is replaced by F.

To determine thermodynamic potentials and related state functions for the setup sketched in Fig. 1, we have to model the microscopic details of the capacitor system. To capture the essential physics of EDLs in nanopores, we model the pores of each electrode by two parallel planar walls of surface area A and separation (pore size) L, as sketched in Fig. 1(b); note that, for our MD simulations, we used oppositely charged walls, as explained in Sec. II D and  Appendix A. Both walls combined carry the total charge ±Q of the respective electrodes, leading to a surface charge density = Q/(2A) with the proton charge e.

To benefit from symmetries, we consider each pore wall stretching infinitely in the (x, y) plane of a Cartesian coordinate system such that edge effects are suppressed. The pore walls are positioned at z = 0 and z = L. In this setting, the EDLs at the left and right wall generally overlap. However, if L is sufficiently large, both EDLs can be considered independent, a situation we call free of overlap. In this case, a variation of L does not affect the EDL. If the state is free of overlap and ions are symmetric, the study of only one EDL at one wall is sufficient because the EDLs at all other walls of both electrodes will follow from symmetries.

In this work, we mainly focus on one electrode (two walls) and define the pore volume AL of one of the electrodes to contain N+ positive ions, N negative ions, and Ns neutral solvent particles. The ions of the electrolyte have valencies zi that define the number of positive unit charges e per ion. Unless stated otherwise, we consider z± = ±1. In cases where we do not explicitly account for the volume of solvent particles, we set Ns = 0. While we generally assume solvent particles to have zs = 0, if explicitly present at all, the dielectric nature of the solvent is always accounted for in a dielectric background via a relative permittivity εr.

The experiments of our interest dealt with porous carbon electrodes and aqueous sodium chloride.1 While we usually use states free of overlap by setting L to large values, choosing L = 1.6 nm would result in the same ratio of pore volume to electrode surface area as in the experiments. To account for the steric effects of the ions, we adopt a restricted primitive model (RPM) and describe the sodium and chloride ions as charged hard spheres. As diameter we choose d = d± = 0.68 nm, a value determined in scattering measurements and approximating their effective size in water.22,42 We also study the solvent primitive model (SPM), an extension of the RPM where solvent particles are added as uncharged hard spheres of ds = 0.3 nm such that the total volume fraction is 46.8% (corresponding to pure water in our model).

Apart from our MD simulations, where particles are considered explicitly, our other theoretical treatments handle particle densities ρi(r) of a particle species i ∈ { +, −, s}, i.e., the number of particles at a position r averaged over all states in an ensemble. Due to the infinite extension of the pore walls, the number densities ρi(r) depend only on z. We further denote the respective bulk densities by ρ¯i. They have to satisfy the condition 0=iziρ¯i to ensure electroneutrality in the bulk.

We construct the local unit charge density as

q(z)=σ[δ(z)+δ(zL)]+i{+,}ziρi(z)
(7)

using Dirac δ-distributions δ(z). The Poisson equation now relates q(z) to the electrostatic potential ψ(z) through

ε0εrz2ψ(z)=eq(z)
(8)

with ε0 being the dielectric permittivity of the vacuum. We will frequently use the dimensionless potential ϕ(z) = (z)/kBT with kB being the Boltzmann constant. Using capital letters, we denote the electrode potential by Ψ = ψ(z = 0) and Φ = ϕ(z = 0). As sketched in Fig. 1(b), we set ψ(0) = ψ(L).

To come to a closed set of equations, Eqs. (7) and (8) need to be supplemented with an expression for ρ±(z) in terms of ψ(z). Accordingly, in the following, we present different theoretical approaches, formulated within the framework of classical density functional theory (DFT).

The central quantity in DFT is the grand potential Ω[{ρi}], a functional of the particle densities ρi in the system. While the grand potential functional also depends on T and μi, for readability, we omit these dependencies in our notation. As common, we split up the grand potential functional into

Ω[{ρi}]=Fid[{ρi}]+Fexc[{ρi}]+AiρizVextzμidz
(9)

with the intrinsic free energy functional Fid[{ρi}] of an ideal gas,40 an excess free energy functional Fexc[{ρi}] that adds contributions due to pair potentials, and a contribution from an external potential Vextr and chemical potentials (the latter enter in the Legendre transform between F and Ω). Importantly, the grand potential functional is minimal for the correct (physical) equilibrium particle densities of the system and its value then equals the value of the actual (thermodynamic) grand potential.43 This property allows us to determine equilibrium density profiles by minimizing a given functional.

While the ideal free energy functional is known exactly, exact excess free energy functionals are only known in a few cases. Nevertheless, many approximations have been tested for specific problems. In this work, we employ three well-established approximations to Fexc: one to describe point-charge particles and two to describe particles that additonally occupy volume in space. For the point charges, we use a mean-field Coulomb functional (cf. Ref. 22) that reads

FC=eA2q(z)ψ(z)dz.
(10)

We refer to this simplest choice by PB because the Euler–Lagrange equations of this functional together with Fid yield the well-known Poisson–Boltzmann equation.20 For the next approach, we extend the above mean-field Coulomb functional by an excess lattice-gas functional Flg that treats the occupied volume of the particles effectively via a maximum local number of allowed particles.44 The respective free energy functional apart from FC reads

Fid+Flg=kBTAi{+,}ρi(z)lnρi(z)ρvac(z)ρMlnρMρvac(z)dz,
(11)

where the density of lattice vacancies is ρvac(z) = ρMρ+(z) − ρ(z) with ρM defining the highest local concentration or, in other words, the number density of accessible lattice sites. The latter is determined from assuming random close packing of hard spheres, resulting in ρMd3π/6 = 0.634.45 We refer to this approach by mPB because the functional in Eq. (11) together with FC yields the modified Poisson–Boltzmann equation35,44 with βμ±=ln(ρ¯±/(ρM2ρ¯±)). Finally, we construct a functional for the RPM and SPM by extending the excess free energy functional Fexc by the non-local “White Bear mark II” functional for hard spheres46 (as in previous work, we additionally apply a correction by Tarazona47). This latter functional allows us to describe hard-sphere interactions between particles and between particles and the walls by employing fundamental measure theory (FMT).48 We refer to this approach by FMT. For explicit expressions and for details on the calculation of the functionals (via Picard iterations and solving the Poisson equation), we refer to previous work.22 Adding a Stern layer for PB and mPB is discussed in  Appendix B. Note that electrostatic interactions beyond the mean-field are still neglected in our Coulomb functional.

To determine the reversible heat, we go back to the previous result in Eq. (4b) now. In a canonical system, the reversible heat produced while charging our model system from surface charge density 0 to σ can also be expressed as (see  Appendix C)

QrevII=A0σ0Lψ(z)ddσiziρi(z) dz dσ.
(12)

We introduce the superscript II to distinguish between our previous result in Eq. (4b) and this result in Eq. (12), where microscopic details enter explicitly through the density profiles of the system. Thus, this method is suitable to determine the reversible heat from DFT data if calculations are performed for a canonical system.

As an additional approach, we study our system of interest through molecular dynamics (MD) simulations. For this purpose, we use the ESPResSo software package49 with the velocity Verlet algorithm for the propagation of the particles in our system. Hence, no real hard-sphere interaction can be used. Instead, we mimic the hard-core interactions by an extremely repulsive Weeks–Chandler–Andersen (WCA) potential,50 essentially a cut and shifted Lennard-Jones potential, that reads

VWCA(r)4ϵ=dLJr12dLJr6+14forr21/6dLJ0otherwise.
(13)

In  Appendix A, we explain our choice of the parameters ϵ and dLJ and verify that this choice yields neutral-sphere density profiles consistent with DFT (FMT) calculations. For the electrostatic interactions, ESPResSo provides the P3M method, a sophisticated Ewald method, as well as an electric layer correction (ELC) method to effectively remove the periodicity in one direction. We use both methods in a three-dimensional simulation box with periodic boundary conditions such that the periodicity in the x- and y-direction accounts for the translational invariance of the system in those directions and the periodicity in the z-direction is suppressed [see Fig. 1(b)].

To model the effects of the charged walls, we first ensured that the EDLs were free of overlap such that we could run simulations with surface charges of opposite sign on both plates (see also the discussion in Sec. II B). Then, we applied an additional constant electric field Ez along the z-direction to all particles, which equals the field induced solely by the surface charges. The corresponding electrostatic potential difference between both walls for a given electric field strength Ez is then obtained by

ΔΨ=mzε0εrAEzL,
(14)

where mz = e∫z(z+ρ+(z) + zρ(z)) dz is the electric dipole moment of the collective distribution of the ions along the z-direction.

Next, we show that Qrev/Wel differs dramatically between charging processes at either fixed μ±, ρ¯±, or N±. For illustrative purposes, we use the Gouy–Chapman solution to the Poisson–Boltzmann equations in this section: This solution allows for (semi) analytical expressions for Qrev/Wel under the three above thermodynamic conditions. Charging processes at constant μ± and N± are most easily treated in the well-known grand canonical and canonical ensembles, respectively. We refer to the charging at fixed bulk densities ρ¯± as semi-canonical, as Qrev generated under this thermodynamic condition turns out to be close to the heat generated in large canonical systems.

1. Recap of the Gouy–Chapman solution

Gouy and Chapman solved the Poisson–Boltzmann equations for a setup of one planar charged hard wall next to an infinite reservoir of 1:1 electrolyte for which ρ¯+=ρ¯. The solution reads18,19

ϕ(z)=4arctanhexpzλDtanhΦ4,
(15a)
   ρ±(z)=ρ¯±expϕ(z),
(15b)

where λD is the Debye length with λD2=4πλBiziρ¯i and λB = e2/(4πε0εrkBT) is the Bjerrum length. Note that Eq. (15) can be easily reformulated for a general |zi|:|zi| electrolyte and so can the results that we derive below with Eq. (15).

From Eq. (15a), the surface charge density with Gauss’s law is as follows:

σ=σ¯sinhΦ2,
(16)

where σ¯=4λDρ¯+. Moreover, inserting Eq. (15) into Eq. (9), we find [cf. Eqs. (24) and (25) of Ref. 34]

ΩGC=ΔΩelGC+ΔΩentGCpV
(17a)

with the bulk pressure p=kBTiρ¯i of an ideal gas and

   ΔΩelGCAkBT=σ¯coshΦ21,
(17b)
ΔΩentGCAkBT=σ¯33coshΦ2+ΦsinhΦ2.
(17c)

Here, ΩGC was partitioned into an entropic contribution ΔΩentGC and an energetic contribution ΔΩelGC.34 The energetic contribution stems from the mean-field Coulomb functional FC, while the contribution ΔΩentGCpV equals FidiμiNi. However, as pointed out by Overbeek, FC only yields purely energetic terms if a constant, in particular, temperature-independent, dielectric constant, is used (see also  Appendix D).

We will demonstrate in Subsection II E 2 that ΔΩentGC as defined in Eq. (17c) does not fulfill the correspondence given in Eq. (6) [cf. Eq. (18)]. However, as it turns out, ΔΩentGC is closely related to the entropic contribution to the free energy in canonical systems (cf. Sec. II E 4).

2. Gouy–Chapman at fixed μ± (grand canonical)

Inserting ΩGC from Eq. (17) into Eqs. (3a) and (3b) to obtain Qrev, we find

QrevGCAkBTσ¯=92μ+kBTcoshΦ21ΦsinhΦ2.
(18)

In  Appendix D, we show that inserting Φ(σ) [as follows from inverting Eq. (16)] into Eq. (4a) yields the same expression for QrevGC.

Figure 2(a) shows the ratio QrevGC/Wel>0 of reversible heat and electric work for grand canonical charging. Here, we used the atomic mass 23u corresponding to sodium to determine QrevGC and we used Eqs. (5) and (17a) to determine Wel from ΔΩGC. We observe QrevGC/Wel>0 up to ΔΨ ≈ 0.9 V, suggesting that heat flows into the system during charging, contradicting the experimental findings of Ref. 1. This positive ratio is caused by the net ion adsorption in both electrodes ΔN=4Aσ¯coshΦ/210 within the system during charging when more counterions are attracted than coions are expelled [see also Eq. (19) and  Appendix E]. Yet, an entropy contribution from increasing particle numbers is unlikely to have occurred in the experiments of Ref. 1: The system of porous electrodes and electrolyte reservoir used there, though certainly large, was closed and, hence, canonical. We conclude that the reversible heat QrevGC with its uncommon explicit dependence on the ionic chemical potentials (hence on Planck’s constant and ionic mass alike) is not relevant for the experimental setup of Ref. 1.

FIG. 2.

Ratio QrevGC/Wel of negative reversible heat and electric work during isothermal charging as a function of applied potential ΔΨ calculated with the Gouy–Chapman solution (a) in a grand canonical setting at fixed μ± and (b) in a canonical setting at fixed N± for several pore widths L. Panel (b) also shows the ratio ΔΩentGC/ΔΩGC in a semi-canonical setting, as follows from Eq. (17), which corresponds to fixed ρ¯±.

FIG. 2.

Ratio QrevGC/Wel of negative reversible heat and electric work during isothermal charging as a function of applied potential ΔΨ calculated with the Gouy–Chapman solution (a) in a grand canonical setting at fixed μ± and (b) in a canonical setting at fixed N± for several pore widths L. Panel (b) also shows the ratio ΔΩentGC/ΔΩGC in a semi-canonical setting, as follows from Eq. (17), which corresponds to fixed ρ¯±.

Close modal

3. Gouy–Chapman at fixed N± (canonical)

Going from a grand canonical to a canonical description, the total numbers of particles per species N±tot are kept fixed during charging rather than the chemical potentials μ±. Hence, tracing the system states in a two-dimensional (N±tot,μ±) diagram during charging, grand canonical systems move along lines of constant μi, whereas canonical systems move along lines of constant N±tot. As our model of the supercapacitor consists of two charged hard walls for each electrode, the total numbers of both cations N+tot and anions Ntot read

N±tot=A0Lρ±(z) dzat pos.electrode+A0Lρ±(z) dzat neg.electrode.
(19)

We consider systems whose EDLs are free of overlap: The smallest L = 10 nm used in this subsection is still much larger than λD ≈ 0.31 nm. Then, we determined Wel during canonical charging with the Gouy–Chapman solution as follows: For a given Φ, we inserted ρ±(z) from Eq. (15) into Eq. (19) and varied ρ¯± until the prescribed N±tot was attained. Clearly, the bulk densities ρ¯± decrease, while Φ increases at fixed N±tot.51 For each combination of Φ and ρ¯±, we find σ with Eq. (16), after which Wel follows from Eq. (2) straightforwardly (see also  Appendix E). Next, to determine the heat Qrev produced during canonical charging, we are confronted with the problem that DFT is formulated in the grand canonical ensemble. However, since our system is assumed to be infinite along the in-plane directions, the equivalence between the thermodynamic potentials (here Ω and F) holds true. We may thus perform a Legendre transform to obtain the free energy of our system as F = Ω + ∑iμiNi, where we use ΩGC [Eq. (17)] in place of Ω. Numerically calculating the derivative of the free energy with respect to temperature then yields the reversible heat.

Figure 2 shows ratios Qrev/Wel of reversible heat and electric work for canonical systems of different lengths L. Here, the expected sign, corresponding to heat flowing out of the system during charging, is obtained. We also note that, though all systems considered are free of EDL overlap, Qrev/Wel depends markedly on L. This is because the smaller the system, the faster ρ¯± decreases during canonical charging. The connected reservoir in the experiment of Ref. 1 being large in comparison to the volume filled by EDLs and desalination of the bulk being negligible during charging brings up the question as to how the ratio Qrev/Wel behaves in the limit of large systems where L/λD → ∞. Based on the arguments of Subsection II E 2, we do not expect Qrev/Wel to be the same in grand canonical and canonical processes in this limit.

4. Gouy–Chapman at fixed ρ¯± (semi-canonical)

The authors of Ref. 1 found an expression for Qrev using Eq. (12), which holds for canonical systems only, inserting, however, Gouy–Chapman density profiles pertaining to a grand canonical system. Interestingly, their expression for QrevII also follows from the combination of Eqs. (6) and (17c). Now, an identical expression for QrevI can be obtained from Eqs. (3) and (4)if the partial derivatives therein are carried out not at fixed μi (imperative in grand canonical settings), but for constant ρ¯±, that is, ρ¯± independent of T, σ, and hence Φ.

Importantly, at fixed ρ¯±, the chemical potentials μi vary with T and the particle numbers N± vary with Φ. Hence, fixed-ρ¯± charging is neither grand canonical nor canonical, and we call it “semi-canonical” instead. Meanwhile, as the ionic density profiles for given ρ¯± are the same in grand canonical and semi-canonical systems, they have the same Φ(σ)-relation and, through Eq. (2), the same Wel as well.

In Fig. 2, we plot the ratio of reversible heat and electric work obtained from the Gouy–Chapman solution via Eq. (17). As discussed above, the ratio ΔΩentGC/ΔΩGC in this semi-canonical system indeed represents the limiting ratio in the canonical system for increasing amounts of connected bulk and, thus, increasing L. It is astonishing that calculations carried out for this semi-canonical process at fixed ρ¯± do not only reflect the conditions of the experiment much better but also simplify calculations [e.g., Eqs. (D3) and (D4)].

5. Conclusion

As demonstrated using the Gouy–Chapman solution, the ratio Qrev/Wel depends strongly on the used thermodynamic conditions. Next to conventional grand canonical and canonical charging processes, we introduced a third process, namely, a semi-canonical charging process at constant ρ¯±. This process mimics a charging process in a system connected to an infinite bulk such that the system in combination with the bulk is canonical. Even though the density profiles for different σ or Φ but constant T are the same as in a grand canonical system, the reversible heat produced during semi-canonical charging resembles the heat generated in a canonical system instead. Importantly, this semi-canonical process reflects the conditions of the experiment, as described in Ref. 1 best and, accordingly, it is used in the following.

We consider a parameter set corresponding to the experiment of Ref. 1: T = 300 K, ρ¯±=1 M, d = 0.68 nm, and εr = 80 (in all approaches, including SPM), resulting in a Bjerrum length of around λB = 0.696 nm. Note that d only enters in Fhs and Flg. Moreover, we now set L = 10 nm and discuss narrower pores later in Sec. III C.

In Fig. 3, we show results from FMT RPM for the reversible heat QrevI and QrevII. Both calculations via Eq. (3a) (dots) and Eq. (12) (solid line) clearly yield the same result numerically. We further check numerically whether the free-energy contribution of the hard-sphere interaction within FMT adds to the entropic contribution to the grand potential. For this purpose, we calculate the contributions Fhs and Fid with the density profiles from FMT; we calculate the particle numbers via Eq. (19). The two resulting terms ΔFhs and Δ(FidiμiNi), shown with colored areas in Fig. 3, add up to precisely ΔΩent=Qrev [Eq. (6)]. Thus, as expected, the contribution of the hard-sphere excess term goes completely into the entropic part. We performed the same checks in the mPB model. Again, grouping the volume exclusion term (Flg) into the ΔΩent term, we verified ΔΩent=Qrev also for mPB.

FIG. 3.

Negative reversible heat Qrev per area A obtained via Eq. (3a) (dots) and Eq. (12) (solid line) from FMT RPM during an isothermal charging process. In addition, we show the change of the contributions Fhs (upper orange area) and FidiμiNi (lower blue area) to the grand potential.

FIG. 3.

Negative reversible heat Qrev per area A obtained via Eq. (3a) (dots) and Eq. (12) (solid line) from FMT RPM during an isothermal charging process. In addition, we show the change of the contributions Fhs (upper orange area) and FidiμiNi (lower blue area) to the grand potential.

Close modal

In Conclusion, in two new cases, we have numerically verified that ΔΩent=Qrev only holds if the excess functional accounting for steric interactions is grouped into the entropic contribution ΔΩent of the grand potential. From here on, we prefer to speak about Qrev/Wel instead of ΔΩent/ΔΩ, although the latter has been used in previous work.1,36 This is because both Qrev and Wel are unambiguously defined in Eqs. (2), (3a), and (12) and can be measured experimentally.

In Fig. 4, we show the ratio Qrev/Wel obtained from the different methods introduced before as well as the experimental data of Ref. 1. Explicitly, we compare the ratios obtained via PB without [red dotted line, same data as in Fig. 2(b)] and with a Stern layer (thick blue dotted line), mPB without (purple dashed-dotted line) and with a Stern layer (thick pink dashed-dotted line), FMT RPM (cyan solid line), FMT SPM (teal dashed line), and MD simulations (orange crosses).

FIG. 4.

Ratio Qrev/Wel of negative reversible heat and electric work during isothermal charging as a function of applied potential ΔΨ for PB without (red dotted line) and with a Stern layer (thick blue dotted line), mPB without (purple dashed-dotted line) and with a Stern layer (thick pink dashed-dotted line), FMT RPM (cyan solid line), FMT SPM (teal dashed line), and MD simulations (orange crosses), all determined for the semi-canonical process at fixed ρ¯±=1 M with L = 10 nm (EDLs are free of overlap for L ≥ 4 nm), d = 0.68 nm, and |z±| = 1. Also shown with black dots are the experimental results from Ref. 1.

FIG. 4.

Ratio Qrev/Wel of negative reversible heat and electric work during isothermal charging as a function of applied potential ΔΨ for PB without (red dotted line) and with a Stern layer (thick blue dotted line), mPB without (purple dashed-dotted line) and with a Stern layer (thick pink dashed-dotted line), FMT RPM (cyan solid line), FMT SPM (teal dashed line), and MD simulations (orange crosses), all determined for the semi-canonical process at fixed ρ¯±=1 M with L = 10 nm (EDLs are free of overlap for L ≥ 4 nm), d = 0.68 nm, and |z±| = 1. Also shown with black dots are the experimental results from Ref. 1.

Close modal

First, we notice that PB, the only description where all steric interactions among the particles are neglected, deviates from all other curves in the way that QrevWel at large potentials. At small applied potentials, we have Qrev=Wel/2. This is all in perfect agreement with earlier descriptions by Overbeek.34 Second, as speculated in Ref. 1, mPB theory describes the experimental data much better than PB theory. PB and mPB coincide at small applied potentials, which is understood from their equal leading-order expansion in ΔΨ, the Debye–Hückel equation.35 Conversely, for ΔΨ much beyond the thermal voltage, Qrev/Wel radically changes: PB predicts the ratio to rise to 1 for large potentials, while mPB predicts this ratio to decrease with increasing potential instead. A similar qualitative change upon accounting for steric repulsions was found in Ref. 36. Similar conclusions for PB and mPB also hold when we add a Stern layer, as explained in  Appendix B. Accounting for a Stern layer, however, dramatically alters the ratio Qrev/Wel at low applied potentials. Interestingly, the value of Qrev/Wel around 0.24 agrees well with the more sophisticated FMT approaches that we discuss now. We see that the predictions of FMT both for RPM and SPM are almost equal and agree with all experimental data within two standard deviations. The similarity of the RPM and SPM results, however, does not mean that solvent properties do not affect Qrev/Wel. For example, the SPM does not account for dipolar interactions within water, which might influence Qrev/Wel. At large potentials, RPM and SPM predictions for the ratio Qrev/Wel are similar to those from mPB, but RPM and SPM predict this ratio to be roughly constant, while mPB predicts this ratio to still decrease with increasing potential. Interestingly, PB predictions with a Stern layer are also similar to those from FMT, even for large potentials. At small potentials, RPM and SPM deviate significantly from Qrev=Wel/2, as predicted by the approaches without a Stern layer (PB and mPB). In both RPM and SPM, particles cannot get closer to the wall than half a particle diameter, introducing a Stern-like layer, whereas in PB and mPB particles can get arbitrarily close to the wall.

Furthermore, we performed MD simulations for ρ¯±=1 M with 600 particles per species in a box of 10 nm3. The MD simulations give access to the (equilibrium) internal energy U and the electric work Wel done to the system [see Eq. (2)], from which the heat flowing into the system follows as Q=ΔUWel. In Fig. 4, we see that MD (orange crosses) predicts slightly larger values for Qrev/Wel than FMT and mostly describes the experiment worse. We performed a convergence analysis for the parameter L by regarding this ratio for different system lengths. We found that the statistical error for potential differences ΔΨ ≥ 0.2 V is smaller than the used marker size and hence negligible compared to the experimental uncertainty.

The deviation between MD and FMT predictions can have different reasons: Most probably, the approximate mean-field functional used for the electrostatic interactions in the FMT approach simply does not capture crucial contributions. For instance, it is known that inaccurate approaches, like the FMT approach, predict qualitatively wrong adsorption to weakly charged walls in the RPM.52 Another similar reason could be the treatment of image charges in both MD and FMT. While image charges are not captured in our MD, it is not yet understood whether they are captured in the ensemble averaged DFT approach. Nevertheless, from the agreement between FMT and PB with a Stern layer, we conclude that the complex structure of the ion density profiles near the surface of the charged hard wall, as predicted only by the FMT approaches and MD, is less important than the steric interaction between the charged wall and the ions represented by the Stern layer.

Using the FMT RPM approach, we discuss how different parameters affect Qrev/Wel.

1. Pore size L

Figure 5 shows Qrev/Wel for several L and all other parameters as in Fig. 4, for which, in particular, λD ≈ 0.31 nm. For interacting EDLs obtained for small pore sizes around LλD, one can see a rapid increase in Qrev/Wel with decreasing L. This finding is relevant to many supercapacitor experiments, where pores in electrodes are nanometer sized, and hence, strong EDL overlap can be expected. For systems larger than L = 4 nm, no further effects from L on Qrev/Wel can be seen because the EDLs decay almost completely within half a system length. As stated earlier, to get the same ratio of pore volume to surface area as in the experiment, a pore size of L = 1.6 nm must be used. For such a pore size, the agreement between the experiment and the FMT RPM curve would be much worse. However, the effective diameter of the ions used here includes some contributions due to hydration shells. These hydration shells might be partly shed when ions get adsorbed at the electrode, causing some effective increase in pore size due to shrinking effective ion sizes. This issue cannot be completely resolved here with the excess functional Fexc that we use and needs further investigation.

FIG. 5.

Ratio Qrev/Wel for FMT RPM as in Fig. 4, but for different system lengths L.

FIG. 5.

Ratio Qrev/Wel for FMT RPM as in Fig. 4, but for different system lengths L.

Close modal

2. Ionic diameter d

Figure 6 shows Qrev/Wel for several d and all other parameters as in Fig. 4. We also show the analytical Gouy–Chapman solution to the Poisson–Boltzmann equations for point charges that follows from Eq. (17). We see that, with d → 0, the RPM results move progressively toward the PB predictions. Notably, however, even for (unphysically) small d, the RPM qualitatively differs from PB as it does not approach 1 but rather decreases at large ΔΨ. We interpret these results as follows: The change in entropy upon charging is associated with the increasing order in the system when ions separate. As PB predicts larger Qrev/Wel than RPM in Fig. 4, the change in entropy upon charging is strongest for point charges. Steric interactions of the hard spheres in the RPM model counteract this trend to order and thus decrease the ratio Qrev/Wel.

FIG. 6.

Ratio Qrev/Wel for FMT RPM as in Fig. 4, but for different particle diameters d. For small d, this ratio approaches PB predictions (red dotted line).

FIG. 6.

Ratio Qrev/Wel for FMT RPM as in Fig. 4, but for different particle diameters d. For small d, this ratio approaches PB predictions (red dotted line).

Close modal

3. Bulk density ρ¯±

Figure 7 shows Qrev/Wel for several ρ¯± and all other parameters as in Fig. 4. Clearly, Qrev/Wel decreases with increasing ρ¯±. We have verified that |Qrev| decreases with increasing ρ¯± as well, as was found in experiments.53 Note that, in the RPM, a phase transition would be expected at a packing fraction roughly above 0.45 that corresponds to a bulk concentration of ρ¯±=5.4666 nm−3 = 9.1M,54 which is far from our system at 1M. Furthermore, the reduced temperature T* = d/λB = 0.977 in our system is much larger than the temperatures where gas–liquid phase separation occurs.54 

FIG. 7.

Ratio Qrev/Wel for FMT RPM as in Fig. 4, but for different bulk concentrations ρ¯±.

FIG. 7.

Ratio Qrev/Wel for FMT RPM as in Fig. 4, but for different bulk concentrations ρ¯±.

Close modal

4. Valencies zi

Figure 8 shows Qrev/Wel for several binary multivalent electrolytes. In order to get an electrically neutral bulk, the bulk densities ρ¯± must be changed accordingly. We choose bulk densities such that ρ¯++ρ¯=2 M to get the same total bulk particle density as before. All other parameters are as in Fig. 4. We see that Qrev/Wel decreases with the amount of total charge defined as ei|zi|ρ¯i. Next to integer charges, we also studied one case of fractional charges, relevant for effective charges of larger molecules as they occur, for instance, in the description of ionic liquids.55 

FIG. 8.

Ratio Qrev/Wel for FMT RPM as in Fig. 4, but for different valencies |z+|:|z|.

FIG. 8.

Ratio Qrev/Wel for FMT RPM as in Fig. 4, but for different valencies |z+|:|z|.

Close modal

Note that, due to the simple form of the Coulomb interaction, a multiplicative factor for the valencies can be theoretically mapped onto different εr or T. For example, the data with valencies zi = 0.5 correspond to one with zi = 1 at a four times higher temperature. Note that also the surface charges, and hence the electric potential, would have to be rescaled.

In this study, we calculated the ratio Qrev/Wel between the negative of the reversible heat emerging during isothermal charging processes and applied electric work for the RPM, SPM, and PB-type models of electrolytes. Our findings agree with the experimental results from Ref. 1 and demonstrate the importance of ionic steric interactions with the charged wall to explain the experiments theoretically.

First, we found that the ratio of reversible heat and electric work depends sensitively on the thermodynamic conditions. A semi-canonical process explains the experiments from Ref. 1 best. This process describes charging in a system connected to a very large bulk reservoir such that the reservoir has constant bulk density during particle exchange as in a grand canonical system. Simultaneously, the bulk density does not change with temperature such that the chemical potential does change as in a canonical system. We demonstrated and discussed this finding for different ensembles and processes; for simplicity, we used the Gouy–Chapman solution to PB theory.

For the calculation of the reversible heat and electric work in the RPM, we used different theoretical approaches, namely, classical DFT and MD simulations. To describe the ionic volume in DFT, we used a modified PB approach as well as a sophisticated approach using FMT. Furthermore, we performed calculations for point-like ions within Poisson–Boltzmann theory with and without a Stern layer to emphasize the importance of steric wall–ion interactions for a description of the experiments. While mPB and FMT predictions for Qrev/Wel are similar for large ΔΨ, they deviate significantly around ΔΨ = 0 V. Such differences would in principle be experimentally testable though the data of Ref. 1 for ΔΨ < 0.2 V have large uncertainty. Conversely, the addition of a Stern layer to PB and mPB yields predictions in excellent agreement with FMT even at low potentials. We point out that deviations between MD and FMT most probably arise from the inaccurate treatment of the electrostatic interactions in our used functional, but other sources, such as the different treatment of image charges, are possible as well. Meanwhile, discrepancies between our theoretical calculations and the experimental data could have different causes: In our model, we neglected adsorption and faradaic reactions and, further, the geometry in our model is an oversimplification of real porous electrodes. Moreover, the treatment of the solvent as a constant dielectric background means that we cannot describe the shedding of an ion’s hydration shell when it enters an ultranarrow pore. Notwithstanding these reservations, our results point toward the important role of the finite size of particles to heat production experiments of capacitive systems.

For this reason, we have further tested the importance of volume effects of the electrolyte solvent that is not contained in the RPM: We performed additional FMT calculations in the SPM, where steric interactions of solvent particles are treated explicitly; solvent particles are described as neutral hard spheres, while we retained the dielectric background of the RPM. We could not find significant deviations between our calculations for RPM and SPM.

Anticipating more experimental data and having at hand a predictive theoretical approach, we further studied how the ratio Qrev/Wel changes with pore size L, ion size d, salt concentration ρ¯±, and ionic valencies z±. We found that the ratio depends sensitively on all these parameters. Experimental heat production measurements should be able to pick these trends up.

In the future, experiments for the ratio Qrev/Wel could become a valuable tool to test aspects of EDL theories. However, the large experimental uncertainty of the available data at small applied potentials1 hinders this at the moment. We thus hope that our study inspires more experimental work on this topic. Furthermore, our findings are of interest for applications, where EDLs change cyclically or where they are exploited in combination with temperature changes, for instance, in heat-conversion processes.

Finally, DFT has the advantage to provide thermodynamic potentials innately. While DFT is limited to equilibrated systems, dynamical DFT allows us to study non-equilibrium processes, for instance, fast (dis)charging of supercapacitors where Joule heating comes into play as well.56–58 Such a framework could be of interest for several applications. For instance, recently, the heat produced in nervous conduction has been found to contain a large fraction of reversible heat, while irreversible contributions are small, if existing at all.9 Now, our work could guide studies, where the theoretical description of nervous conduction goes beyond ideal solutions and homogeneous bulk concentrations, as typically applied, and, thus, could shed new light on this fundamental neuroscience process.

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

F.G. and A.H. acknowledge support by the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through Grant No. INST 39/963-1 FUGG (bwForCluster NEMO). The research leading to these results has received funding from the European Union’s Horizon 2020 research and innovation program under Marie Skłodowska-Curie Grant Agreement No. 801133. Finally, F.G. acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Project No. 430195928—and A.H. acknowledges funding by the DFG—Project No. 406121234.

For the MD simulations, we had to set the WCA parameters entering Eq. (13). For the strength of the repulsion, we used ϵ = 100kBT, as in previous work.37 Note that, in the literature, many different methods to obtain an effective radius can be found.50,59 We set dLJ such that density profiles from DFT and MD agree well for a system of neutral hard spheres which happens when 21/6dLJ is about 1.4% larger than the hard-sphere diameter. To obtain this value, we demand that the first Mayer f-bond contribution for a hard-sphere system equals the one for the WCA potential.40 

As the counterion density near the walls shoots up with the applied potential, one needs to check if the simulation box is sufficiently large to capture the spatial correlations along the lateral directions. Accordingly, we calculated the radial distribution function projected on the lateral plane for particles close to the wall. We checked that the projected radial distribution function decays to zero within half a lateral box length. Furthermore, we checked that the length of the system along the normal direction is long enough such that, as for DFT, the different EDLs do not interact nor desalinate the bulk. We found that a box of 10 nm3 with 600 particles per species meets the desired conditions. The particle numbers in the simulations are chosen such that we find the bulk densities from DFT calculations in the center of the simulation box. This allows us to use oppositely instead of equally charged walls in our simulations. Hence, a homogeneous electric field can be used to mimic the effects of surface charges.

The finite size of ions affects both the ion–ion interaction and the ion–wall interaction. A simple way to account for finite ion size (in the ion–wall interaction) is through a Stern layer, which is a charge-free region reaching from the electrode surface into the electrolyte over the ionic radius d/2. From Eq. (8) and Gauss’s law, zϕ = −4πλBσ follows the potential difference over the Stern layer as ΦStern = 2πλBσd. Note that the same potential drop ΦStern applies to PB and mPB theories.

In Eqs. (17b) and (17c), we expressed ΩelGC as a function of the potential Φ. Using, instead of Eq. (16),

Φ=2 arcsinhσσ¯+2πλBσd,
(B1)

we can express ΩelGC as a function of σ as

ΩelGC+Stern(σ)=ΩelGC(σ)+e2σ2Ad4ε0εr.
(B2)

As ΩentGC(σ) is unaffected by the Stern layer, we obtain the ratio Qrev/Wel as a function of σ as

QrevWel(σ)=ΩentGC(σ)ΩentGC(σ)+ΩelGC+Stern(σ).
(B3)

Using this result together with Eq. (B1) to obtain the potential Φ, we obtain the thick dotted blue line shown in Fig. 4.

For a thermodynamical (dis)charging process, the difference in heat δQ can be calculated from the difference in internal energy dU and the work δW done during the process as

δQ=dUδW.
(C1)

If the only non-zero contribution to the internal energy comes from the electrostatic interaction, we have

dU=dε0εrA20L(E(z))2dz
(C2)
=ε0εrA0LdE(z)E(z)dz
(C3)
=eA0Ldq(z)ψ(z)dz.
(C4)

In the last step, an integration by parts is performed. Furthermore, one should keep in mind that the charge distribution q(z) contains contributions from both the ions ∑iziρi(z) and the surface charges on the electrodes σ.

For a process of duration T, the change in internal energy follows as

ΔU=0TdUdtdt
(C5)
=eA0T0Ldq(z)dtψ(z)dz dt
(C6)
=A0T0LI(z)E(z)dz dt,
(C7)

where the continuity equation is used and another integration by parts is performed.

Similarly, the electric work during (dis)charging can be written as

ΔWel=0TdQdtΔΨdt.
(C8)

Subtracting Eq. (C8) from Eq. (C7) [or, respectively, Eq. (C6)], one is left with the ionic currents. However, one should keep in mind that this difference yields the corresponding heat only if there are no other work terms. In a grand canonical charging process, for example, one would also get a work term due to a particle flux into/out of the system,

ΔWch=iμiΔNi.
(C9)

If one wants to calculate the reversible heat, an infinitely slow (dis)charging process must be regarded where the system is in equilibrium at every time. Variable substitutions in Eqs. (C6) and (C8), replacing the time integrals by integrals over surface charge density, yield

QrevII=A0σ0Lψ(z)ddσiziρi(z)dz dσ.
(C10)

Note that we introduced index II solely for conformity with the main text. This equation is very useful to calculate the reversible heat from DFT data.

For the Gouy–Chapman solution, the heat flow near a single charged wall follows from Eqs. (4a) and (16) and the definition of σ¯ as

QrevI,GCAkBT=20σTsinh1σ/σ¯Tα,σ dσ
(D1)
=20σsinh1σσ¯σσ2+σ¯2lnσ¯lnTα dσ
(D2)
=σ¯ΦsinhΦ2+21coshΦ2×1+12ln[εr(T)T]lnT+12lnρ¯+lnTα.
(D3)

Here, α stands for the variable(s) kept fixed during the partial T-derivative. If we consider the ensemble of fixed ionic concentration (α=ρ¯+=ρ¯), the last term in Eq. (D3) drops out. If, in addition, (εr(T)T)/∂T = 0, which is accurate for water (cf. p. 69 in Ref. 34), Eq. (D3) yields

QrevI,GCAkBT=σ¯22coshΦ2+ΦsinhΦ2
(D4)
=ΔΩentGC+ΔΩelGC,
(D5)

the same as Eq. (S11) in Ref. 1 (up to a typo in their Ω subscript). This, however, would mean that at every point of charging, the amount of electric work put into the system would flow out of the system in the form of reversible heat. Thus, the internal energy of the system would remain constant during charging, which would be surprising. This may be resolved by including some explicit model for the solvent responsible for the T dependence of εr(T). The explicit model would yield further entropic contributions and may also resolve the problem that, for T-dependent εr(T), one gets QrevIQrevII, as was found in Ref. 1. If (εr(T))/∂T = 0 is considered instead, we find

QrevI,GCAkBT=σ¯33coshΦ2+ΦsinhΦ2
(D6)
=ΔΩentGC.
(D7)

In the grand canonical ensemble (α = μ+ = μ), we make use of μ±=kBTlnΛ±3ρ¯± to find

lnρ¯±lnTμ±=32μ±kBT.
(D8)

For the case that [εr(T)]/∂T = 0, one now finds

QrevI,GCAkBT=σ¯ΦsinhΦ2+1coshΦ292μ+kBT,
(D9)

which is Eq. (18) of the main text.

In a canonical charging process, the total number of ions must be conserved. To derive the corresponding equation for ρ¯±, a system free of overlap is assumed. For simplicity, every electrode is assumed to be one charged hard wall in order to get rid of an additional factor of 2 that cancels anyway for the ratio Qrev/Wel.

Following Ref. 51 [Eqs. (7)(8)(9)(10)], we see that the total number of ions per species in our system can be written as

N±A=ΔNA+ρ¯±L,
(E1)

where ΔN is

ΔNA=σ¯coshΦ21.
(E2)

This equation can be solved for ρ¯± where one should keep in mind that σ¯ is ρ¯± dependent.

In combination with Eq. (16), we have thus two equations to numerically search for combinations of σ, Φ, and ρ¯± that solve these equations. From ρ¯±, one can simply calculate μi, and thus, we can calculate Wel and F then.

1.
M.
Janssen
,
E.
Griffioen
,
P. M.
Biesheuvel
,
R.
van Roij
, and
B.
Erné
,
Phys. Rev. Lett.
119
,
166002
(
2017
).
2.
J.
Schiffer
,
D.
Linzen
, and
D. U.
Sauer
,
J. Power Sources
160
,
765
(
2006
).
3.
J.
Lindner
,
F.
Weick
,
F.
Endres
, and
R.
Schuster
,
J. Phys. Chem. C
124
,
693
(
2020
).
4.
H. P.
Rodenburg
, “
Measuring the ion-specific heat of electrical double layer formation in porous carbon
,” Master’s thesis,
Faculty of Science
,
Utrecht University
,
2020
.
5.
A.
d’Entremont
and
L.
Pilon
,
J. Power Sources
246
,
887
(
2014
);
A. L.
d’Entremont
and
L.
Pilon
,
J. Power Sources
273
,
196
(
2015
).
6.
R.
Kumar
,
J. P.
Mahalik
,
E.
Strelcov
,
A.
Tselev
,
B. S.
Lokitz
,
S.
Kalinin
, and
B. G.
Sumpter
, “
Microscopic theory for electrocaloric effects in planar double layer systems
,” arXiv:1503.09141 (
2015
).
7.
M.
Janssen
and
R.
van Roij
,
Phys. Rev. Lett.
118
,
096001
(
2017
).
8.
C.
Cruz
,
A.
Ciach
,
E.
Lomba
, and
S.
Kondrat
,
J. Phys. Chem. C
123
,
1596
(
2019
).
9.
A. C. L.
de Lichtervelde
,
J. P.
de Souza
, and
M. Z.
Bazant
,
Phys. Rev. E
101
,
022406
(
2020
).
10.
A.
Alizadeh
and
M.
Wang
,
Electrophoresis
41
,
1067
(
2020
).
11.
H.
Helmholtz
,
Ann. Phys.
165
,
211
(
1853
).
12.
B. E.
Conway
,
Electrochemical Supercapacitors. Scientific Fundamentals and Technological Applications
(
Springer
,
1999
).
13.
E. J. W.
Verwey
,
J. T. G.
Overbeek
, and
K.
Van Nes
,
Theory of the Stability of Lyophobic Colloids: The Interaction of Sol Particles Having an Electric Double Layer
(
Elsevier Publishing Company
,
1948
).
14.
B.
Derjaguin
and
L.
Landau
,
Prog. Surf. Sci.
43
,
30
(
1993
).
15.
R.
Roth
,
D.
Gillespie
,
W.
Nonner
, and
R. E.
Eisenberg
,
Biophys. J.
94
,
4282
(
2008
).
16.
A.
Peyser
,
D.
Gillespie
,
R.
Roth
, and
W.
Nonner
,
Biophys. J.
107
,
1841
(
2014
).
17.
A. A.
Kornyshev
,
D. J.
Lee
,
S.
Leikin
, and
A.
Wynveen
,
Rev. Mod. Phys.
79
,
943
(
2007
).
18.
M.
Gouy
,
J. Phys. Theor. Appl.
9
,
457
(
1910
).
19.
D. L.
Chapman
,
London, Edinburgh Dublin Philos. Mag. J. Sci.
25
,
475
(
1913
).
20.
J.-L.
Barrat
and
J.-P.
Hansen
,
Basic Concepts for Simple and Complex Liquids
(
Cambridge University Press
,
Cambridge
,
2003
).
21.
J.-P.
Hansen
and
H.
Löwen
,
Annu. Rev. Phys. Chem.
51
,
209
(
2000
).
22.
A.
Härtel
,
J. Phys.: Condens. Matter
29
,
423002
(
2017
).
23.
F.
Coupette
,
A. A.
Lee
, and
A.
Härtel
,
Phys. Rev. Lett.
121
,
075501
(
2018
).
24.
P.
Simon
and
Y.
Gogotsi
,
Nat. Mater.
7
,
845
(
2008
).
25.
K. V. G.
Raghavendra
,
R.
Vinoth
,
K.
Zeb
,
C. V. V.
Muralee Gopi
,
S.
Sambasivam
,
M. R.
Kummara
,
I. M.
Obaidat
, and
H. J.
Kim
,
J. Energy Storage
31
,
101652
(
2020
).
26.
P.
Simon
and
Y.
Gogotsi
,
Nat. Mater.
19
,
1151
(
2020
).
27.
M. E.
Suss
,
S.
Porada
,
X.
Sun
,
P. M.
Biesheuvel
,
J.
Yoon
, and
V.
Presser
,
Energy Environ. Sci.
8
,
2296
(
2015
).
28.
D.
Brogioli
,
Phys. Rev. Lett.
103
,
058501
(
2009
).
29.
M.
Janssen
,
A.
Härtel
, and
R.
van Roij
,
Phys. Rev. Lett.
113
,
268501
(
2014
).
30.
A.
Härtel
,
M.
Janssen
,
D.
Weingarth
,
V.
Presser
, and
R.
van Roij
,
Energy Environ. Sci.
8
,
2396
(
2015
).
31.
L.
Landau
,
J.
Bell
,
M.
Kearsley
,
L.
Pitaevskii
,
E.
Lifshitz
, and
J.
Sykes
,
Electrodynamics of Continuous Media
, Course of Theoretical Physics (
Elsevier Science
,
2013
).
32.
S.
Porada
,
H. V. M.
Hamelers
, and
P. M.
Biesheuvel
,
Phys. Rev. Res.
1
,
033195
(
2019
).
33.
P. M.
Biesheuvel
and
J. E.
Dykstra
,
Physics of Electrochemical Processes
(
published online by authors
,
2020
), Chap. 9, ISBN: 978-90-9033258-1.
34.
J. T. G.
Overbeek
,
Colloids Surf.
51
,
61
(
1990
).
35.
V.
Kralj-Iglič
and
A.
Iglič
,
J. Phys. II France
6
,
477
(
1996
).
36.
P. M.
Biesheuvel
and
M.
van Soestbergen
,
J. Colloid Interface Sci.
316
,
490
(
2007
).
37.
A.
Härtel
,
M.
Janssen
,
S.
Samin
, and
R.
van Roij
,
J. Phys.: Condens. Matter
27
,
194129
(
2015
).
38.
C.
Ebner
,
W. F.
Saam
, and
D.
Stroud
,
Phys. Rev. A
14
,
2264
(
1976
).
39.
R.
Evans
,
Adv. Phys.
28
,
143
(
1979
).
40.
J. P.
Hansen
and
I. R.
McDonald
,
Theory of Simple Liquids: With Applications to Soft Matter
(
Elsevier Science
,
New York
,
2013
).
41.
D.
Frenkel
and
B.
Smit
,
Understanding Molecular Simulation: From Algorithms to Applications
, Computational Science (
Elsevier Science
,
2001
).
42.
E. R.
Nightingale
,
J. Phys. Chem.
63
,
1381
(
1959
).
43.
N. D.
Mermin
,
Phys. Rev.
137
,
A1441
(
1965
).
44.
I.
Borukhov
,
D.
Andelman
, and
H.
Orland
,
Phys. Rev. Lett.
79
,
435
(
1997
).
45.
C.
Song
,
P.
Wang
, and
H. A.
Makse
,
Nature
453
,
629
(
2008
).
46.
H.
Hansen-Goos
and
R.
Roth
,
J. Phys.: Condens. Matter
18
,
8413
(
2006
).
47.
P.
Tarazona
,
Phys. Rev. Lett.
84
,
694
(
2000
).
48.
R.
Roth
,
J. Phys.: Condens. Matter
22
,
063102
(
2010
).
49.
F.
Weik
,
R.
Weeber
,
K.
Szuttor
,
K.
Breitsprecher
,
J.
de Graaf
,
M.
Kuron
,
J.
Landsgesell
,
H.
Menke
,
D.
Sean
, and
C.
Holm
,
Eur. Phys. J.: Spec. Top.
227
,
1789
(
2019
).
50.
J. D.
Weeks
D.
Chandler
H. C.
Andersen
J. Chem. Phys.
54
5237
1971
51.
N.
Boon
and
R.
van Roij
,
Mol. Phys.
109
,
1229
(
2011
).
52.
D.
Gillespie
,
M.
Valiskó
, and
D.
Boda
,
J. Phys.: Condens. Matter
17
,
6609
(
2005
).
53.
X.
Zhang
,
W.
Wang
,
J.
Lu
,
L.
Hua
, and
J.
Heng
,
Thermochim. Acta
636
,
1
(
2016
).
54.
A.-P.
Hynninen
,
M. E.
Leunissen
,
A.
van Blaaderen
, and
M.
Dijkstra
,
Phys. Rev. Lett.
96
,
018303
(
2006
).
55.
D.
Roy
and
M.
Maroncelli
,
J. Phys. Chem. B
114
,
12629
(
2010
).
56.
M.
Schmidt
and
J. M.
Brader
,
J. Chem. Phys.
138
,
214101
(
2013
).
57.
M.
Schmidt
,
Phys. Rev. E
84
,
051203
(
2011
);
J. G.
Anero
,
P.
Español
, and
P.
Tarazona
,
J. Chem. Phys.
139
,
034106
(
2013
).
[PubMed]
58.
A. A.
Lee
,
S.
Kondrat
,
D.
Vella
, and
A.
Goriely
,
Phys. Rev. Lett.
115
,
106101
(
2015
).
59.
J. A.
Barker
and
D.
Henderson
,
J. Chem. Phys.
47
,
4714
(
1967
).