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.

## I. INTRODUCTION

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 experimental^{1–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 solutes^{27} and to harvest energy because concentration^{28} and temperature^{8,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 $I\u2192\u22c5E\u2192$ in the heat equation^{5,7,31} with $I\u2192$ being the ionic current density and $E\u2192$ being the local electric field. While the ionic current density aligns with the local electric field $I\u2192\u221dE\u2192$ in bulk electrolytes, leading to a strictly positive Joule-heating term $\u221dI\u21922$, conversely $I\u2192\u22c5E\u2192<0$ is possible in the EDL when the gradient in electrochemical potential anti-aligns with $E\u2192$, 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 $Qtot\u2212Qirr$ 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 *W*_{el} during isothermal charging. With $Qrev$ and *W*_{el} 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=\u2212\Delta \Omega ent/\Delta \Omega $ relied on the identity $Qrev=\u2212\Delta \Omega ent$ [cf. Eq. (6)], proposed by Overbeek on thermodynamic grounds,^{34} and on ΔΩ = *W*_{el}, which holds for isothermal charging. The linear scaling of both $Qrev$ and *W*_{el} 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 theory^{34} 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.

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=\u2212\Delta \Omega 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.

## II. THEORY

### A. Thermodynamics

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 *N*_{s} 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

The electric work performed during charging is given by

Using Legendre transforms, we obtain the free energy *F*(*T*, *Q*, *N*_{i}) = *U*(*S*, *Q*, *N*_{i}) − *TS* and the grand potential Ω(*T*, *Q*, *μ*_{i}) = *F*(*T*, *Q*, *N*_{i}) − ∑_{i}*μ*_{i}*N*_{i}, 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 *Q*_{1} to *Q*_{2} must cause a (reversible) heat flow into the capacitor that reads

Here, the entropy difference Δ*S* is given by

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]

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 *N*_{i} 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,

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

### B. Microscopic model setup

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 *eσ* = *Q*/(2*A*) 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 *N*_{s} neutral solvent particles. The ions of the electrolyte have valencies *z*_{i} 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 *N*_{s} = 0. While we generally assume solvent particles to have *z*_{s} = 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 *d*_{s} = 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 $\rho i(r\u2192)$ of a particle species *i* ∈ { +, −, s}, i.e., the number of particles at a position $r\u2192$ averaged over all states in an ensemble. Due to the infinite extension of the pore walls, the number densities $\rho i(r\u2192)$ depend only on *z*. We further denote the respective bulk densities by $\rho \xafi$. They have to satisfy the condition $0=\u2211izi\rho \xafi$ to ensure electroneutrality in the bulk.

We construct the local unit charge density as

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

with *ε*_{0} being the dielectric permittivity of the vacuum. We will frequently use the dimensionless potential *ϕ*(*z*) = *eψ*(*z*)/*k*_{B}*T* with *k*_{B} 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*).

### C. EDL modeling within classical density functional theory

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

with the intrinsic free energy functional $Fid[{\rho i}]$ of an ideal gas,^{40} an excess free energy functional $Fexc[{\rho i}]$ that adds contributions due to pair potentials, and a contribution from an external potential $Vextr\u2192\u2009$ 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

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

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 *ρ*_{M}*d*^{3}*π*/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 equation^{35,44} with $\beta \mu \xb1=ln(\rho \xaf\xb1/(\rho M\u22122\rho \xaf\xb1))$. 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 spheres^{46} (as in previous work, we additionally apply a correction by Tarazona^{47}). 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)

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.

### D. MD simulations

As an additional approach, we study our system of interest through molecular dynamics (MD) simulations. For this purpose, we use the ESPResSo software package^{49} 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

In Appendix A, we explain our choice of the parameters *ϵ* and *d*_{LJ} 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 *E*_{z} 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 *E*_{z} is then obtained by

where *m*_{z} = *e∫z*(*z*_{+}*ρ*_{+}(*z*) + *z*_{−}*ρ*_{−}(*z*)) d*z* is the electric dipole moment of the collective distribution of the ions along the *z*-direction.

### E. $Qrev$ depends sensitively on the boundary conditions of the charging process

Next, we show that $Qrev/Wel$ differs dramatically between charging processes at either fixed *μ*_{±}, $\rho \xaf\xb1$, 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 $\rho \xaf\xb1$ 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 $\rho \xaf+=\rho \xaf\u2212$. The solution reads^{18,19}

where *λ*_{D} is the Debye length with $\lambda D\u22122=4\pi \lambda B\u2211izi\rho \xafi$ and *λ*_{B} = *e*^{2}/(4*πε*_{0}*ε*_{r}*k*_{B}*T*) is the Bjerrum length. Note that Eq. (15) can be easily reformulated for a general |*z*_{i}|:|*z*_{i}| 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:

where $\sigma \xaf=4\lambda D\rho \xaf+$. Moreover, inserting Eq. (15) into Eq. (9), we find [cf. Eqs. (24) and (25) of Ref. 34]

with the bulk pressure $p=kBT\u2211i\rho \xafi$ of an ideal gas and

Here, Ω^{GC} was partitioned into an entropic contribution $\Delta \Omega entGC$ and an energetic contribution $\Delta \Omega elGC$.^{34} The energetic contribution stems from the mean-field Coulomb functional $FC$, while the contribution $\Delta \Omega entGC\u2212pV$ equals $Fid\u2212\u2211i\mu 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 $\Delta \Omega entGC$ as defined in Eq. (17c) does not fulfill the correspondence given in Eq. (6) [cf. Eq. (18)]. However, as it turns out, $\Delta \Omega 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 $\mu \xb1$ (grand canonical)

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 *W*_{el} 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 $\Delta N=4\u2009A\sigma \xafcosh\Phi /2\u22121\u22650$ 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.

#### 3. Gouy–Chapman at fixed *N*_{±} (canonical)

Going from a grand canonical to a canonical description, the total numbers of particles per species $N\xb1tot$ are kept fixed during charging rather than the chemical potentials *μ*_{±}. Hence, tracing the system states in a two-dimensional $(N\xb1tot,\mu \xb1)$ diagram during charging, grand canonical systems move along lines of constant *μ*_{i}, whereas canonical systems move along lines of constant $N\xb1tot$. 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 $N\u2212tot$ read

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 *W*_{el} during canonical charging with the Gouy–Chapman solution as follows: For a given Φ, we inserted *ρ*_{±}(*z*) from Eq. (15) into Eq. (19) and varied $\rho \xaf\xb1$ until the prescribed $N\xb1tot$ was attained. Clearly, the bulk densities $\rho \xaf\xb1$ decrease, while Φ increases at fixed $N\xb1tot$.^{51} For each combination of Φ and $\rho \xaf\xb1$, we find *σ* with Eq. (16), after which *W*_{el} 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}*μ*_{i}*N*_{i}, 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 $\u2212Qrev/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, $\u2212Qrev/Wel$ depends markedly on *L*. This is because the smaller the system, the faster $\rho \xaf\xb1$ 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 $\u2212Qrev/Wel$ behaves in the limit of large systems where *L*/*λ*_{D} → ∞. Based on the arguments of Subsection II E 2, we do not expect $\u2212Qrev/Wel$ to be the same in grand canonical and canonical processes in this limit.

#### 4. Gouy–Chapman at fixed $\rho \xaf\xb1$ (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 $\rho \xaf\xb1$, that is, $\rho \xaf\xb1$ independent of *T*, *σ*, and hence Φ.

Importantly, at fixed $\rho \xaf\xb1$, the chemical potentials *μ*_{i} vary with *T* and the particle numbers *N*_{±} vary with Φ. Hence, fixed-$\rho \xaf\xb1$ charging is neither grand canonical nor canonical, and we call it “semi-canonical” instead. Meanwhile, as the ionic density profiles for given $\rho \xaf\xb1$ are the same in grand canonical and semi-canonical systems, they have the same Φ(*σ*)-relation and, through Eq. (2), the same *W*_{el} 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 $\Delta \Omega entGC/\Delta \Omega 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 $\rho \xaf\xb1$ 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 $\rho \xaf\xb1$. 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.

## III. RESULTS

We consider a parameter set corresponding to the experiment of Ref. 1: *T* = 300 K, $\rho \xaf\xb1=1\u2009M$, *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.

### A. Check of $Qrev=\u2212\Delta \Omega ent$ [Eq. (6)]

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 $\Delta Fhs$ and $\Delta (Fid\u2212\u2211i\mu iNi)$, shown with colored areas in Fig. 3, add up to precisely $\Delta \Omega ent=\u2212Qrev$ [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 $\Delta \Omega ent=\u2212Qrev$ also for mPB.

In Conclusion, in two new cases, we have numerically verified that $\Delta \Omega ent=\u2212Qrev$ 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 *W*_{el} are unambiguously defined in Eqs. (2), (3a), and (12) and can be measured experimentally.

### B. $Qrev/Wel$ within theoretical approaches

In Fig. 4, we show the ratio $\u2212Qrev/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).

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 $Qrev\u2192\u2212Wel$ at large potentials. At small applied potentials, we have $Qrev=\u2212Wel/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, $\u2212Qrev/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 $\u2212Qrev/Wel$ at low applied potentials. Interestingly, the value of $\u2212Qrev/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 $\u2212Qrev/Wel$. For example, the SPM does not account for dipolar interactions within water, which might influence $\u2212Qrev/Wel$. At large potentials, RPM and SPM predictions for the ratio $\u2212Qrev/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=\u2212Wel/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 $\rho \xaf\xb1=1$ M with 600 particles per species in a box of $10\u2009nm3$. The MD simulations give access to the (equilibrium) internal energy *U* and the electric work *W*_{el} done to the system [see Eq. (2)], from which the heat flowing into the system follows as $Q=\Delta U\u2212Wel$. In Fig. 4, we see that MD (orange crosses) predicts slightly larger values for $\u2212Qrev/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.

### C. Influence of pore size, ionic diameter, bulk density, and valencies on $Qrev/Wel$

Using the FMT RPM approach, we discuss how different parameters affect $\u2212Qrev/Wel$.

#### 1. Pore size *L*

Figure 5 shows $\u2212Qrev/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 $\u2212Qrev/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 $\u2212Qrev/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.

#### 2. Ionic diameter *d*

Figure 6 shows $\u2212Qrev/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 $\u2212Qrev/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 $\u2212Qrev/Wel$.

#### 3. Bulk density $\rho \xaf\xb1$

Figure 7 shows $\u2212Qrev/Wel$ for several $\rho \xaf\xb1$ and all other parameters as in Fig. 4. Clearly, $\u2212Qrev/Wel$ decreases with increasing $\rho \xaf\xb1$. We have verified that $|Qrev|$ decreases with increasing $\rho \xaf\xb1$ 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 $\rho \xaf\xb1=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}

#### 4. Valencies *z*_{i}

Figure 8 shows $\u2212Qrev/Wel$ for several binary multivalent electrolytes. In order to get an electrically neutral bulk, the bulk densities $\rho \xaf\xb1$ must be changed accordingly. We choose bulk densities such that $\rho \xaf++\rho \xaf\u2212=2\u2009M$ to get the same total bulk particle density as before. All other parameters are as in Fig. 4. We see that $\u2212Qrev/Wel$ decreases with the amount of total charge defined as $e\u2211i|zi|\rho \xafi$. 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}

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 *z*_{i} = 0.5 correspond to one with *z*_{i} = 1 at a four times higher temperature. Note that also the surface charges, and hence the electric potential, would have to be rescaled.

## IV. DISCUSSION AND CONCLUSIONS

In this study, we calculated the ratio $\u2212Qrev/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 $\rho \xaf\xb1$, 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 potentials^{1} 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.

## DATA AVAILABILITY

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

## ACKNOWLEDGMENTS

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.

### APPENDIX A: SIMULATION DETAILS

For the MD simulations, we had to set the WCA parameters entering Eq. (13). For the strength of the repulsion, we used *ϵ* = 100*k*_{B}*T*, 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 *d*_{LJ} such that density profiles from DFT and MD agree well for a system of neutral hard spheres which happens when 2^{1/6}*d*_{LJ} 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\u2009nm3$ 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.

### APPENDIX B: ADDING A STERN LAYER TO THE PB AND mPB APPROACH

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 $\Omega elGC$ as a function of the potential Φ. Using, instead of Eq. (16),

we can express $\Omega elGC$ as a function of *σ* as

As $\Omega entGC(\sigma )$ is unaffected by the Stern layer, we obtain the ratio $\u2212Qrev/Wel$ as a function of *σ* as

### APPENDIX C: REVERSIBLE HEAT FROM INTEGRATED HEAT PRODUCTION

For a thermodynamical (dis)charging process, the difference in heat $\delta Q$ can be calculated from the difference in internal energy d*U* and the work $\delta $*W* done during the process as

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

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 ∑_{i}*z*_{i}*ρ*_{i}(*z*) and the surface charges on the electrodes *σ*.

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

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

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

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,

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

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.

### APPENDIX D: DERIVATION OF $QrevI$

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

Here, *α* stands for the variable(s) kept fixed during the partial *T*-derivative. If we consider the ensemble of fixed ionic concentration ($\alpha =\rho \xaf+=\rho \xaf\u2212$), 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

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 $QrevI\u2260QrevII$, as was found in Ref. 1. If *∂*(*ε*_{r}(*T*))/*∂T* = 0 is considered instead, we find

In the grand canonical ensemble (*α* = *μ*_{+} = *μ*_{−}), we make use of $\mu \xb1=kBT\u2061ln\Lambda \xb13\rho \xaf\xb1$ to find

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

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

### APPENDIX E: ADSORPTION IN THE CANONICAL GOUY–CHAPMAN SOLUTION

In a canonical charging process, the total number of ions must be conserved. To derive the corresponding equation for $\rho \xaf\xb1$, 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 $\u2212Qrev/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

where Δ*N* is

This equation can be solved for $\rho \xaf\xb1$ where one should keep in mind that $\sigma \xaf$ is $\rho \xaf\xb1$ dependent.

In combination with Eq. (16), we have thus two equations to numerically search for combinations of *σ*, Φ, and $\rho \xaf\xb1$ that solve these equations. From $\rho \xaf\xb1$, one can simply calculate *μ*_{i}, and thus, we can calculate *W*_{el} and *F* then.