Parameter-Fitting-Free Continuum Modeling of Electric Double Layer in Aqueous Electrolyte

Electric double layers (EDLs) play fundamental roles in various electrochemical processes. Despite the extensive history of EDL modeling, there remain challenges in the accurate prediction of its structure without expensive computation. Herein, we propose a predictive multiscale continuum model of EDL that eliminates the need for parameter fitting. This model computes the distribution of the electrostatic potential, electron density, and species’ concentrations by taking the extremum of the total grand potential of the system. The grand potential includes the microscopic interactions that are newly introduced in this work: polarization of solvation shells, electrostatic interaction in parallel plane toward the electrode, and ion-size-dependent entropy. The parameters that identify the electrode and electrolyte materials are obtained from independent experiments in the literature. The model reproduces the trends in the experimental differential capacitance with multiple electrode and nonadsorbing electrolyte materials (Ag(110) in NaF, Ag(110) in NaClO4, and Hg in NaF), which verifies the accuracy and predictiveness of the model and rationalizes the observed values to be due to changes in electron stability. However, our calculation on Pt(111) in KClO4 suggests the need for the incorporation of electrode/ion-specific interactions. Sensitivity analyses confirmed that effective ion radius, ion valence, the electrode’s Wigner–Seitz radius, and the bulk modulus of the electrode are significant material properties that control the EDL structure. Overall, the model framework and findings provide insights into EDL structures and predictive capability at low computational cost.


INTRODUCTION
The electric double layer (EDL) plays a fundamental role in various electrochemical processes, such as colloidal dispersions, surface charging, and charge-transfer reactions.It appears at every electrode/electrolyte interface because of the work function difference between the two materials.The work function difference induces a strong electric field at the interface.The electric field redistributes the solute and solvent molecules in the electrolyte and forms a surface charge and dipole by either attracting or repulsing the molecules to or from the electrode surface.At the same time, electrons in the electrode also redistribute to form a surface charge and dipole by spilling into the electrolyte. 1The surface charge due to the electron spillover is the same value in the opposite sign as that due to the redistribution of the molecules in the electrolyte, which ensures overall electroneutrality.The layered structure formed by the redistribution of electrons in the metal and molecules in the electrolyte is called the EDL.
The increasing focus of EDL studies lies in its impact on reaction kinetics.The electric field strength and electrostatic potential near the interface are determined by the EDL structure, which can impact the concentration and stability of the reactants.Ringe et al. 2 and Shin et al. 26 suggested that the electric field strength at the interface can impact the stability of CO 2 adsorbate and thus the reaction kinetics of CO 2 reduction.Similar effects of EDLs on the reaction kinetics can appear in any catalytic reactions in electrochemical systems.A deeper understanding of these effects will provide guidance for tailoring the electrode/electrolyte interface with improved reaction kinetics.
−5 C d is defined as C d ≔ dQ/ dV, where Q is the surface charge density and V is the electrode potential.This property represents the potential dependence of the ion distribution in the electrolyte.C d can be experimentally obtained by cyclic voltammetry (CV) or electrochemicalimpedance spectroscopy (EIS).In CV measurements, however, separating EDL capacitance from total capacitance can be challenging when a Faradaic reaction occurs in the potential range of interest. 6EIS measurements with appropriate analysis can separate those two currents based on the difference in their time scales. 4−10 This model envisions the EDL as two separate layers: Stern layer and diffuse layer.The Stern layer is a layer between the electrode surface and the closest approach of hydrated ions.In this layer, hydrated ions cannot exist because of the size exclusion.Without specific ion adsorption, the Stern layer does not contain charged species.The dielectric constant in the Stern layer can be much smaller than that of the bulk solvent (∼an order of magnitude) due to impacts of water structure 11−15 and field effects 16−18 near the interface.Unlike in the Stern layer, ions can redistribute, depending on the electrostatic potential in the diffuse layer.The distribution is obtained from the Poisson equation and the conditions for thermodynamic equilibrium.The analytical solution for the ion distribution is called Gouy−Chapman distribution when the chemical potential of ions (μ i ) is described in the form of where k B is the Boltzmann constant, T is the temperature, n i is the number density of solute i, n i 0 is the number density of i in the bulk, z i is the ionic valency, e 0 is the elementary charge, and ϕ is the electrostatic potential.Bikerman 19 improved the accuracy of the diffuse-layer model by modifying the entropy term (first term on the right side) to account for the finite size of the ions.He obtained an analytical expression for the entropy term by assuming that all of the ions in the electrolyte have the same molar volume. 20The continuum model based on GCS theory with Bikerman modification is widely used as a primitive EDL model because of its computational efficiency.
Conventional continuum models, however, are not capable of "predicting" the EDL structure by themselves.They require parameter fitting, especially for the parameters in the Stern layer.The shortcomings of the conventional continuum models stem from the fact that they do not explicitly account for several fundamental microscopic interactions near the interface 21 such as electron spillover and ion-specific adsorption.
Quantum-mechanical calculations like density functional theory (DFT) 22−24 provide the means for predictive analyses of the interface.DFT calculations capture the effect of electron spillover, as well as the specific interactions between the electrode/solute, electrode/solvent, and solute/solvent.The analysis can be further improved by using ab initio moleculardynamics simulations (AIMDs). 25,26AIMDs can include quantum-mechanical interactions, as well as steric restrictions and entropic contributions.However, because of their large computational cost to evaluate the average structure of the statistical equilibrium, these calculations are limited to a small volume that is typically not sufficient to capture the whole picture of the EDL structure accurately.
The advantages of continuum models and quantummechanical calculations can be achieved by coupling these two models.Joint density functional theory (joint DFT) was proposed to combine the quantum DFTs and classical continuum model. 27,28Huang proposed a simple and effective model named density potential functional theory (DPFT). 20his model combines a classical continuum model based on GCS theory with a Bikerman modification and a simplified quantum-mechanical calculation.For the quantum-mechanical model, he employed a jellium model with uniform electron gas approximation. 22A jellium model, which describes the positive charge distribution in the electrode as a uniform background charge, is one of the simplest models that accounts for the spillover of electrons into the electrolyte.Although the jellium model does not provide results with highest accuracy, it has been used for work function analyses 29−33 and EDL modeling 34−36 due to its computational efficiency.−33 Even with Huang's DPFT model, 20 however, parameter fitting is needed to reproduce the experimental differential capacitance measurements.
This study proposes a new modeling framework that gains the advantages of both continuum models and molecular-dynamics simulations, i.e., accurate and predictive analyses with low computational cost.This model computes the distribution of the electrostatic potential, electron density, and species' concentrations by taking the extremum of the total grand potential of the system.The grand potential is calculated from the entropic energy, electrostatic energy, electron energy, solute−solute interactions, and electrode/solute interactions.These energy expressions include not only the terms from previous papers 20 but also the microscopic interactions that are newly introduced in this work: polarization of solvation shells, electrostatic interaction in parallel plane toward the electrode, and ion-sizedependent entropy.All parameters that identify the electrode and electrolyte materials are obtained from independent experiments; in this way, the model is predictive and does not contain additional fitting parameters.

General Approach.
The EDL is often assumed to be in a quasi-equilibrium 37 because of its rapid formation.The formation time scale can be roughly expressed as 10 l D 2 /D, where l D is the Debye length and D is the diffusion coefficient, 38 and can scale from 10 −9 to 10 −6 s for 1 M to 1 mM of 1:1 electrolytes with a diffusion coefficient of 10 −9 m 2 /s, which is much shorter than time scales of typical electrochemical measurements (e.g., EIS measurement with 100 kHz corresponds to the time scale of 10 −5 s). 1,7Based on the time scale difference, it is reasonable to assume that EDL is in quasi-equilibrium.Also, we focus on potential windows in which the Faradaic reaction rate is small so that the chemical potential gradient of solutes in EDL is negligible.Based on these assumptions, the distribution of electrostatic potential (ϕ), electron density (n e ), and species' number density (n i ) are obtained by taking the extremum of the total grand potential (Ω tot ) of the system. 20,37More specifically, n e and n i are obtained by minimizing Ω tot , while ϕ is obtained by maximizing it. 39Ω tot can be expressed as a functional of these unknown variables (ϕ, n e , and n i ) as  )   where ω tot is the local grand potential.Here, we assume the electrode is flat, and ϕ, n e , and n i distribute uniformly in the plane parallel to the electrode.With this assumption, the grand potential can be simplified in one-dimension (1D) to

Journal of Chemical Theory and Computation
where x describes the position in the direction perpendicular to the electrode.x < 0, x = 0, and x > 0 represent the electrode, electrode/electrolyte interface, and electrolyte, respectively.V is the volume of the system.x = −L metal is the left end of the calculation domain that represents the bulk metal, and x = L sol is the right end of the domain that represents the bulk solution.
The local grand potential can be decomposed into 5 factors.
where ω mix is the mixing entropy, ω els is the electrostatic energy, ω elec is the electron energy, ω ss is the solute−solute interaction, and ω ws is the wall (electrode)-solute nonelectrostatic interaction.The reference potentials of all the grand potential components are based on the electrochemical potential in the bulk electrolyte (denoted by the superscript 0) so that the conditions to take the extremum of Ω̅ tot , which is described in the next section, are satisfied in the bulk.The expressions for ω tot are shown in Section 2.3 and listed in Tables S1−S3 in the Supporting Information (SI).

Variational Analysis. ϕ(x)
, n e (x), and n i (x) are obtained from the condition to let Ω̅ tot [ϕ(x), n e (x),n i (x)] be the extremum.The condition can be written as δΩ̅ tot /δϕ(x) = δΩ̅ tot /δn e (x) = δΩ̅ tot /δn i (x) = 0 for any x.By considering that the local grand potential ω tot depends not only on ϕ, n e , and n i , but also on the first-order gradient of ϕ and n e (∇ϕ and ∇n e ), the conditions are expressed as where ϕ, n e , and n i need to satisfy these equations at any position of x.Equation 5 provides the same number of equations as the number of unknown variables (ϕ(x), n e (x), and n i (x)).proposed a model for the mixing entropy that accounts for the finite size of the species.This expression is frequently employed in continuum models for EDLs because of the availability of analytical solutions.However, the model is associated with the total free energy by assuming that all the species in the solution have the same molar volume, 40−42 which is not always applicable to ions with different molecular size or hydration numbers.
There are several existing models that deal with the asymmetric size of the solute molecules.For example, Gongadze and Iglic 43 proposed a model by modifying the Bikerman model, providing an analytical expression for solute concentrations with molecules' size variation.Another model proposed by Maggs and Podgornik 44 accounted for the entropy due to the solvent.However, a comprehensive formulation of mixing entropy with an asymmetric size effect has yet to be established.The Gongadze's model involves an unphysical segmentation of the solutes to achieve uniform size, while the Maggs' model solely accounts for solvent entropy, neglecting entropy loss of solutes due to the existence of other solutes.
Hence, we derived the mixing entropy (ω mix ) based on the partition function with a lattice model, where ε mix id is the entropy from ideal mixing, ε mix slv is the solvent's entropy considering the finite size of the solute molecules, and ε mix size is the entropy due to the difference in the species' molar volume.μ mix,i is the chemical potential due to the mixing in the bulk solution that is defined to make sure ∂ω mix /∂n i = 0 in the bulk.This mixing entropy is based on the partition function using a lattice model as where N is the number of the solute type, V i is the volume of one molecule of species i, and The derivation of these equations is shown in Section S2.1 in the SI.When all of the solute molecules have the same molar volume, this expression reduces to the conventional Bikerman model. 19.3.2.Electrostatic Interactions, ω els .The grand potential due to electrostatic interactions can be expressed by where ω els ef is the electric field energy, ω els chg is the Coulombic energy of charged species, and ω els pol is the polarization energy of solvent and solute molecules.In this study, ω els ef and ω els chg use conventional expressions based on the Coulombic interactions. 37 where ϵ 0 is the vacuum permittivity, e 0 is the elementary charge, and n e 0 is the valent electron density in the bulk metal.Θ(x) is the Heaviside step function which takes a value of 1 when x > 0 and 0 when x < 0. Θ(x) is used to express the metal properties in x < 0

Journal of Chemical Theory and Computation
and solution properties in x > 0. ϵ op m and ϵ op s are the optical dielectric constant of the metal and solvent, respectively.ω els pol accounts for the dielectric saturation 16,17 and the polarizability of hydrated ions. 45,46In the calculation, we assume that dielectric saturation of ion solvation shell takes place in a similar way to water molecules 17,46 where a pol is a constant that controls the significance of dielectric saturation: large a pol increase the effect of electric field on the effective dielectric constant.a pol is obtained by fitting to an ab initio molecular dynamic analysis. 17ϵ eff ∇ϕ = 0 is the dielectric constant without electric field and is defined as , where β i is the polarizability of solute i, ϵ s 0 is the solvent's dielectric constant, and N avo is the Avogadro's constant.The electrostatic potential in the bulk solution is set to zero.With these expressions of ω els , δΩ̅ tot / δϕ(x) = 0 gives the Poisson equation with the effective dielectric constant ϵ eff as described in eqs S31 and S32 in the SI.In the bulk solution (|∇ϕ| → 0, x > 0, Hence, β i can be obtained by measuring concentrationdependent dielectric constant as β i = −N avo dϵ eff /dn i 045 for the bulk solution. 2.3.3.Electron Energy, ω elec .This study employs electron energy to account for the effect of a potential-dependent electron spillover.We employ a jellium model with uniform electron gas approximation 22 for its simple and low-cost computation.In the jellium model, positive charges in core atoms of the metal (n m ) are expressed as the uniform background charge, and the density distribution of the valence electrons (n e (x)) is evaluated based on the kinetic and exchangecorrelation energy of electrons with the mean-field approximation.Based on Smith's expression, 29 our model evaluates the kinetic energy with the first-order density-gradient expansion, while it uses the local-density approximation for exchangecorrelation energy.Because of the simplified interaction between the core electrons and the valence electrons, jellium models need modifications to improve their accuracy.Perdew et al. 32 introduced a structureless pseudopotential that represents the Madelung energy and the repulsive interaction from the core electrons.The accuracy of Perdew's model in predicting the work function and the bulk modulus of metals is, however, limited mainly because of the assumption of uniform background positive charge in the jellium model.The discrepancy becomes larger when the model is applied to transition and noble metals.Russier and Badiali attributed the discrepancy on transition and noble metals to the contribution of d-electron. 47o correct the errors, we employ the structureless pseudopotential in a semiempirical manner.We employ two parameters in the model that are determined based on basic metal properties.One parameter, structureless pseudopotential 30,32 (μ ps m ), was determined from the bulk modulus of the metal (see eq S29 in the SI).The other parameter, Δϕ WF , was determined from the work function in vacuum and assumed to be independent of the electrode potential (see Section 2.5.1).We also assumed that the electrons that spill into the electrolyte phase feel a non-negligible potential due to the interaction with solvent molecules.This potential was expressed as the pseudopotential in the electrolyte phase.It was determined from the potential of zero charge of Ag(110) 48 and used as a constant value for other metals.
The grand potential for the electron energy is expressed as where ε elec txc is the summation of kinetic, exchange, and correlation energy of the electron based on Smith's expression 29 where e au is the Hartree energy and a 0 is the Bohr radius.ε elec ps is the structureless pseudopotential that takes different values in the metal (x < 0) and in the solution (x > 0), where μ ps m and μ ps s are the pseudopotentials in the metal and solution, respectively, and μ elec is the chemical potential of the electron and is defined as where E WE is the absolute potential of the working electrode and Δϕ WF is a constant used to correct the potential to match the vacuum work function.E WE is the control parameter in the calculation.

Solute−Solute Interaction
, ω ss .We account for the electrostatic solute−solute interactions by assuming that the solvation shells prevent the solutes from getting close to other solute molecules to feel chemical interactions.Although the electrostatic interaction (ω els ) accounts for the electrostatic interaction along the x direction; it does not include the interaction in the y−z plane because it is a one-dimensional model.Ions in the solution interact with each other and redistribute due to electrostatic interaction, which affects the total energy.This y−z plane electrostatic interaction is accounted for as solute−solute interaction.The solute−solute interaction is obtained by analytically solving the Poisson− Boltzmann equation.The governing equations are the same as that of Debye−Huckel theory; however, we used a cylindrical coordinate in the y−z plane, instead of spherical coordinates, to avoid double counting electrostatic interactions in the x direction.The interaction is thus expressed in the form of where ε ss ion is the ion−ion interaction in the y−z plane and μ ss ref is the reference chemical potential expressed as

Journal of Chemical Theory and Computation
where r ave is the average radius of the solute molecules, j is the imaginary unit, Y 0 and Y 1 are the Bessel functions of the second kind of order 0 and order 1, respectively, and x ss is an i n t e r m e d i a t e v a r i a b l e d e s c r i b e d b y In the derivation of this expression, we employed the first-order approximation |z i e 0 ϕ/ k B T| ≪ 1.Although this approximation becomes less accurate when the ion effective radius is small, it enables analytic expression, which is needed to include the effect of y-z direction distribution in the one-dimensional grand potential, ω tot .This interaction is newly introduced in this model, whose derivation is shown in Section S2.2 in the SI.
2.3.5.Wall-Solute/Solvent Interaction, ω ws .To simplify the analysis, this work focuses on a system without a significant specific interaction between the electrode and solute or solvent.Hence, in this study, ω ws only accounts for solute ions' steric restriction for closest approach due to their finite size where μ ws,i is the chemical potential for wall-solute/solvent interactions and is defined as μ ws,i := μ cut Θ(r i − x), μ cut is the cutoff chemical potential (set to 1000 eV), and r i is the species' effective radius.

Boundary Conditions.
Since the problem to be solved is one-dimensional with physics described by second-order differential equations, we need two boundary conditions for each variable.For electrostatic potential, the boundary conditions are set as ϕ(L sol ) = 0 and ∇ϕ(−L metal ) = 0.For electron density, n e (L sol ) = 0 and n e (−L metal ) = n e 0 .From the definition of the reference energy, ∇n e (L sol ) = ∇n e (−L metal ) = 0 is satisfied when L slv and L metal are large enough.For species' density n i , the boundary conditions are ∇n i (L sol ) = ∇n i (0) = 0, which automatically satisfies n i (L sol ) = n i 0 because of the definition of the reference energy.With these boundary conditions, the working electrode (E WE ) was controlled to estimate how the system reacts with the electrode potential change.
2.5.Parameter Setting.Here we describe the parameter setting protocols for electrode-or electrolyte-specific parameters.The values of the parameters for the electrode and electrolyte are listed in Tables 1 and 2, respectively.

Electrode-Specific Parameters.
−53 The Wigner−Seitz radius is employed to evaluate the electron density in the bulk metal, n e 0 := 3/(4πr ws 3 ).The bulk modulus is used to estimate the pseudopotential in the metal (μ ps m ) using the relation B = V(∂ 2 E)/(∂V 2 ) N , where V, E, and N are the volume, energy, and number of electrons in the metal, respectively (see Section S2.3 in the SI).Δϕ WF was assumed to be a constant and was determined from the work function in vacuum as , where Φ vac exp is the experimental work function and Φ vac calc is the calculated work function when Δϕ WF = 0.In the calculation of work function in vacuum, we set the effective dielectric constant to be 1, species' density n i to be 0, and pseudopotential outside of the metal to be 0.
2.5.2.Electrolyte-Specific Parameters.The model employs six electrolyte-specific properties: the ionic valence of the anion and cation (z a and z c ), effective hydrated radius of anion and cation (r a and r c ), and the polarizability of the hydrated shell of anion and cation (β a and β c ). z i is obtained from the ionic formula.β i (:= dϵ eff /dc i ) is from the experimental results for concentration-dependent dielectric potential. 45r i is calculated from the distance between the ion's center and center of the nearest water molecule, 54 including up to the first water layer for anions and second water layer for cations. 55The electrons' pseudochemical potential in the solvent (μ ps s ) also needs to be determined using an experimental result.We employed the potential of zero charge (PZC) on Ag(110) in 5 mM NaF (−0.731V vs SHE 48 ) to get μ ps i = −0.043eV.We fixed this value for all of the calculations using aqueous electrolyte because we assume that it stems from the interaction between the electron from metal and the water molecules. 36In the calculation, the absolute potential for standard hydrogen electrode (SHE) is set as 4.44 V. 56

Evaluation of the Calculated Results.
The model calculates the distribution of ϕ, n e , and n i (cations: i = c, anions i = a) at the given electrode potential (E WE ) as shown in Figure 1.The material-specific parameters are listed in Tables 1 and 2. The electrostatic potential (ϕ) increases under a higher electrode potential, while the electron density (n e ) in the electrolyte phase decreases for higher E WE .The anion density (n a ) is higher when a higher potential is applied, while the opposite trend is observed for cations (n c ).The anion and cation density decrease to almost zero, respectively, at x ≤ 0.51 nm and x ≤ 0.65 nm because of the steric wall-solute interaction (ω ws ) defined in eq 22, corresponding to the Stern layer.These trends are physically representative and reasonably consistent with the previous analyses for the EDLs. 8,9,20  a For anions, only the first hydration shell was included, while second hydration shell was included for cations assuming the stronger hydration affinity of cations. 55For water radius, 0.138 nm was used. 54

Journal of Chemical Theory and Computation
The first term in the parentheses represents the positive background charge in the metal.At a potential of zero charge (PZC), σ ion = σ metal = 0.For comparison with the experiment, the differential capacitance, C d , was calculated by taking the derivative of σ ion with respect to the electrode potential The differential capacitance is experimentally measurable with EIS or CV by using the relation of dσ ion /dt = −I, where I is the current (positive current means an oxidative current).The surface charge density and differential capacitance are shown in Figure 2. The surface charge density is positive below the PZC and negative above the PZC.In the varied potential range, the surface charge density is monotonic.The double-layer capacitance exhibits a local minimum at the PZC because this is where the interface is least charged, and then there are two humps below and above PZC that are not symmetric.This asymmetry arises from anions being smaller in effective hydrated radius compared to the cations, which is discussed in Section 3.3.

Comparison with Experimental Results
. At first, we analyzed Ag(110) in NaF and NaClO 4 electrolytes and compared the differential capacitance predicted by the model with the experimental results by Valette, 48 whose analysis suggests the effect of specific ion adsorption in this system is not significant.Figure 3a,b shows that our model captures the trends in the experiments: the differential capacitance has two local maxima and one local minimum.The local minimum increases when the ion concentration increases, and the height of the first local maximum (−0.9 to −0.8 V vs SHE), in the potential range below the PZC, is almost the same for NaClO 4 and NaF, while the second local maximum (−0.65 to −0.55 V vs SHE), above the PZC, depends strongly on the anion species identity.Because the first local maximum can be attributed to cation adsorption (when the metal is negatively charged) and the second one is due to anion adsorption (when the metal is positively charged; see Section 3.5.), the change in the anion species affects only the second local maximum.The calculated height of the maxima and minima also agrees well with the experiment.Overall, the calculated results demonstrate good agreement with the experimental ones for all the electrolyte concentration and species without parameter fitting, suggesting the predictivity of the model.
Next, we checked the applicability of the model to different electrodes.As a system without significant specific adsorption, we calculated the differential capacitance on Hg in NaF solutions and compared it with the experimental results by Grahame, 57 whose analysis suggests the specific ion adsorption is not significant in this system.The comparison between calculation and experimental results is shown in Figure 3c.Although the agreement is not perfect, the calculation successfully reproduces the different features of Hg compared to Ag: lower differential capacitance with significantly suppressed capacitance maxima.Considering that the model does not use any adjustable parameters, the agreement is satisfactory.The quantitative discrepancy between the experiment and calculation for Hg can probably be attributed to water chemisorption on Hg, 58 which is not accounted for in the present model.The difference in C d between Ag(110) and Hg can be attributed to the difference in the stability of the electrons in the metals.In our model, the metal bulk modulus is associated with the pseudopotential of electrons in the metal.Here, the bulk modulus of Hg is 0.267 Mbar, 50 while that of Ag is 1.01 Mbar. 33s described in eq S29 in the SI, a larger bulk modulus results in higher pseudopotential in the metal because of stronger repulsive interactions from the core electrons.With higher pseudopotential, the electrons become less stable in the metal, and the electron density profile (n e (x)) becomes more sensitive to the applied electrode potential (compare Figures 1b and 4a).The larger sensitivity of n e (x) on E WE results in a larger dσ metal / dE WE , which equals C d (eq 25).When we focus on the Helmholtz capacitance (C H ), the difference between the two electrodes becomes clearer as shown in Figure 4b.Here, C H is calculated by where L zeta = min(r a ,r c ) is the closest approach of ions.Also, we compared the calculated differential capacitance on Pt(111) in 100 mM KClO 4 with the experimental result by Pajkossy and Kolb 4 to check the model applicability to a system with specific adsorption. 59,60Figure 3d shows the results for the Pt(111) electrode.The figure shows a significant discrepancy between the experiment (the blue solid line) and the calculation (the blue dotted line) in terms of peak height, location, and number of peaks.The discrepancy is thought to be due to the specific interaction between Pt(111) surface and ClO 4 − ion, 59 which can be corrected by using an appropriate function for chemical potential for wall-anion interaction (μ ws,a ) (see the red dotted line in Figure 3d).Here we used the form of Morse potential as where G is the adsorption free energy (0.65 eV), α 1 is the relaxation parameter (4.72 nm −1 ), and r a ads is the effective radius of adsorbed anions (0.4 nm).Also, we accounted for the anione effective radius shift as r a (x) = r a 0 + (r a ads − r a 0 )Θ(x − 0.5(r a 0 +r a ads )), where r a 0 is the effective anion radius in the bulk solution.Because of the lack of quantum-mechanical analysis (e.g., DFT calculation) on this interaction, we cannot conclude that the assumed interaction is physically reasonable.However, this result confirmed that the ion-electrode-specific interaction can significantly alter the differential capacitance and needs to be accounted for when applying the model to a system with specific adsorption.Also, it should be noted that quantum-mechanical analysis can be used to assess the significance of the specific interaction term in a certain system in case sufficient experimental data is unavailable.By running a DFT simulation with various ion-metal surface distances, one can obtain a free energy profile.This profile can serve as the wall-solute interaction term (μ ws,a (x)) in a similar manner to eq 27.If the adsorption free energy is sufficiently large, it will affect the calculated differential capacitance, which suggests the importance of including a specific interaction term in the model.
−13 Although our model includes the polarization energy that accounts for the stabilization of solvent dipole for the electric field, it does not include the explicit expressions for the two-dimensional hydrogen bonding in the vicinity of the interface.The model will be capable of including these effects by setting the parameters based on further analyses.The former effect, stabilization of water, prevents solute molecules from approaching the surface.This blocking effect can be included by obtaining the wall-solute interaction energy (μ ws ) from DFT calculations with explicit water molecules, which evaluates the energy to put the solute molecule near the interface, including the reorganization energy of the water hydrogen bond.Also, the latter effect, a change in the dielectric properties, can be implemented by making the dielectric parameters (e.g., ϵ s 0 and a pol ) position-dependent.However, this analysis is beyond the scope of the current study and requires parameter fitting that is not readily available and/or is computationally cost-prohibitive.

Effect of Grand Potential Components.
The effect of each grand potential component is analyzed by enabling them one by one in the calculation for Ag(110) in 100 mM NaF by substituting the total local grand potential (ω tot ) in eq 5 with ω tot (k) in Table 3.
The initial model (ω tot (1) ) corresponds to the standard GCS theory.It accounts for the ideal mixing entropy (ε mix id ), electrostatic energy (ω els ), and wall-solute interaction (ω ws ).Since the electron energy (ω elec ) is missing in this model, the partial derivative of ω tot with n e and ∇n e in eq 5 does not provide any meaningful information.Hence, instead of solving for the differential equation for n e , we manually set n e (x) = n e 0 Θ(−x), which results in neglecting the effect of the surface dipole due to electron spillover.This model also neglects the dielectric saturation due to the polarization of the solvent and solutes by setting a pol → 0. In addition, the solute molecules are assumed to have the same size to exclude the effect of ion size differences.The second model (ω tot (2) ) corresponds to the Bikerman model, 19 which accounts for the effects of finite size of solute ions on mixing entropy (ε mix slv ) by assuming all ions have the same size.The third model (ω tot (3) ) adds the electron energy term (ω elec ) to the second model that enables the evaluation of n e (x), which corresponds to the density potential functional theory  For k = 1 and 2, the electron density is set as n e (x) = n e 0 Θ(−x).
developed by Huang. 20The fourth model (ω tot (4) ) activates the dielectric saturation by setting a pol = 6 nm/V, a fitted value to an ab initio simulation result. 17The fifth model (ω tot (5) ) enables the solute−solute interaction (ω ss ) to visualize its effect on the calculation results.Finally, the sixth model (ω tot (6) ) includes the effect of ion size difference between anions and cations on the mixing entropy by activating ε mix size and setting r a ≠ r c . Figure 5 shows that ω tot (1) gives the differential capacitance with one local minimum, which saturates when the potential is far from the PZC as expected from the GCS theory.The line for ω tot (2) shows that the addition of ε mix slv results in the two local maxima, suggesting that they are due to the ions' finite size.ε mix slv reduces the capacitance far from the PZC because the limited availability of the space for solvent near the interface mitigates further ion accumulation.The line for ω tot (3) shows that the electron energy (ω elec ) increases the capacitance around PZC, and the effect is more significant in negative potential versus PZC than in positive potential.This is because the electron spillover becomes more significant in negative potential, where the electron is relatively unstable in the electrode (thus the Helmholtz capacitance C H becomes larger in negative potential as shown in Figure 4b); this result is consistent with the discussion by Huang. 20The polarization of solvent and solute ions (a pol ≠ 0) reduces the differential capacitance in the potential far from the PZC (the line for ω tot (4) ).This is because the polarization of the ions and solvents induces dielectric saturation, which reduces the effective dielectric constant near the interface and thus the Stern layer capacitance when the electric field is strong. 16The addition of solute−solute interaction (ω ss ) slightly increases the height of the peak (the line for ω tot (5) ) because of the favorable interactions between the ions at higher ion strength.Although the effect of this interaction is not significant in case of monovalent electrolyte, it becomes larger as the ionic valence increases (see Figure S1 in the SI).The line for ω tot (6) shows increased differential capacitance especially above the PZC.This is because the smaller size of the anion enhances the anion accumulation near the interface.Another analysis of the effects of the interaction terms, which uses a different order, also demonstrate significant effects of polarization of solvent and solute ions, as well as the size-dependent entropy (see Figure S2 in the SI).This analysis suggests that all the interactions introduced in this study, polarization of solvent and solute ions, size-dependent entropy, and solute−solute interactions, are responsible for predicting differential capacitance.
3.4.Parameter Sensitivity.Parameter sensitivity analysis is conducted with the model to identify the critical material properties that impact the EDL structure.The parameters to be analyzed are the material properties introduced in Section 2.5: the effective hydrated radius of anion and cation (r a and r c ), polarizability of the hydrated shell of anion and cation (β a and β c ), ionic valence of anion and cation (z a and z c ), Wigner−Seitz radius of the electrode (r ws ), bulk modulus of the electrode (B), and vacuum work function of the electrode (Φ vac ).By changing these material properties by ±5% independently from the base parameters (values on Ag(110) in 100 mM NaF), we evaluated 4 properties of EDL structures: PZC (denoted asE PZC ), C d at E PZC , C d at E PZC −0.1 V, and C d at E PZC + 0.1 V.The difference of ±5% is chosen as a value small enough to obtain a numerical derivative and large enough to neglect the effect of the numerical error due to solver tolerance.Then, the sensitivity of EDL property k on the material property j (S k,j ) was calculated by where m j is the base value of the material property j, Δe k,j is the shift of the EDL property k due to the shift in the material property j, and Δm j is the shift in the material property j.The calculated sensitivities of the EDL properties on the material properties are shown in Figure 6.
First, the effects of the electrolyte properties, r a , r c , β a , β c , z a , and z c are discussed.Figure 6 demonstrates that none of these parameters have a recognizable effect on PZC, which confirms that PZC is independent of electrolyte properties in this model with no specific interactions.The properties of the anion (r a , β a , and z a ) mainly affect C d (E PZC + 0.1 V), while the properties of the cation (r c , β c , and z c ) show similar effects but in C d (E PZC − 0.1 V).Above the PZC, the anion accumulates near the interface (Figure 1c).Since the anion properties affect the anion's accumulation affinity, they can impact the surface charge profile and thus the differential capacitance in high potential.The  cation, on the other hand, accumulates under the PZC (Figure 1d) and thus the cation properties have stronger effects on the differential capacitance in the negative potential.The calculation results show that the electrolyte properties: r a , r c , z a , and z c have large effects on the differential capacitance.Hence, one must take care of the effective ion radius and ionic valence to control the EDL structures.
In terms of the sensitivity of the electrode properties: r ws , B, and Φ vac , Figure 6 demonstrates that r ws has a strong effect on all EDL properties.This strong effect can be attributed to the thirdorder effect of r ws on n e 0 (n e 0 := 3/(4πr ws 3 )) and the impact of n e 0 on the electron spillover (see the shift from ω tot (2) to ω tot (3) in Figure 5).Φ vac significantly changes the PZC but shows no recognizable effects on differential capacitance.In this model, Φ vac is used to evaluate Δϕ WF , which is a constant that simply shifts the absolute potential of the electrode.The shift in the absolute potential directly changes the PZC but does not change the shape of the potential-dependent differential capacitance.The bulk modulus (B) affects all of the EDL properties, but the impact is smaller than r ws .Even with the small sensitivity, it was the main reason for the difference between Ag(110) and Hg, as discussed in Section 3.2, because of the large difference in the bulk modulus (1.01 Mbar for silver 33 and 0.267 Mbar for mercury 50 ).The positive sensitivity of C d on B is consistent with the smaller C d on Hg than that on Ag(110).Although the Wigner−Seitz radius is the most sensitive property of the electrode, one needs to also consider the bulk modulus as an important property to predict the differential capacitance.

Potential Applications of the Model.
Since this model provides a means to predict the structure of EDLs without parameter fitting, it can be applied in various applications, including different materials and reaction microenvironments, whereas conventional continuum EDL models require parameter fitting based on the experimental results.For example, Huang's model 20 calibrates the optical dielectric constant based on the experimental differential capacitance.Since the model presented herein does not require parameter fitting, it can be used to predict the EDL structure without experimental data, which gives us insights into material selection to control the surface properties.
Also, the continuum descriptions developed in this study can be extended into higher dimensions to address the inhomogeneous nature of the electrode surface, which plays a significant role in electrocatalytic processes. 15−64 The present model would improve the accuracy of DFT simulations by providing a more accurate description of the interactions in the electrolyte.Although our model might increase the computational cost due to the need for iterative calculations for solute molecule distributions (n i ), it would improve the calculation accuracy of the implicit water model because of the improved descriptions in the size-dependent entropy term and the polarization model.One should note that in this application, the calculation becomes three-dimensional so that we no longer need to include the solute−solute interaction (ω ss ), which is derived by integrating the Coulomb interaction in the y−z direction as a special treatment for a one-dimensional model.

CONCLUSIONS
In this study, a predictive multiscale continuum model was developed for the electric double layer that eliminates the need for parameter fitting due to the incorporation of microscopic interactions and independent material properties.The model reproduced the trends in the experimental differential capacitance with multiple noninteracting electrode and electrolyte materials (Ag(110) in NaF, Ag(110) in NaClO 4 , and Hg in NaF), which verifies the accuracy and predictiveness of the model.However, the poorer predictions for Pt(111) in KClO 4 suggest the necessity of further study to capture the effect of electrode-ion-specific interactions in a parameter-fitting-free manner.The difference in the differential capacitance between Hg and Ag(110) was attributed to the difference in the electron stability in the metal.Sensitivity analyses confirmed that all the newly incorporated interactions added in this study play a role in predicting the differential capacitance.It also demonstrated the effective ion radius, the ionic valence, the electrode's Wigner− Seitz radius, and the bulk modulus of the electrode are significant material properties that control the EDL structure.The model framework and findings provide insights into the EDL structures and enable predictive evalumation of EDLs with low computational cost.
■ ASSOCIATED CONTENT research process, including literature survey, theory development, summarizing results, and manuscript preparation.Additionally, I.V.Z.provided an initial model for the study.

D
diffusion coefficient (∼10 −9 m 2 /s) Ω tot total grand potential ( eq 2 ) V volume of the system N number of the solute species Ω̅ tot volumetric total grand potential (Ω tot /V) G ads adsorption free energy for ClO − adsorption on Pt(111) (0.40 nm, fitted) S k,j sensitivity of EDL property k on the material property j (S k,j ≔ m j (Δe k,j /Δm j )) m j base value of the material property j Δe k,j shift of the EDL property k due to the shift in the material property j Δm j shift in the material property j Functions Y 0 Bessel function of the second kind of order 0 Y 1 Bessel function of the second kind of order 1 Θ Heaviside step function (Θ(x > 0) = 1 and Θ(x < 0) = 0)

Figure 1 .
Figure 1.Calculated distributions of (a) electrostatic potential (ϕ), (b) electron density (n e ), (c) anion density (n a ), and (d) cation density (n c ) for Ag(110) in 100 mM NaClO 4 .The horizontal axis shows the position in the x direction (x < 0 is in the electrode phase and x > 0 is in the electrolyte phase), and the different line styles represent different electrode potentials as shown in the legend in (a).The inset in (b) is the enlarged view of the domain inside the gray dashed box.N avo is the Avogadro number.The black dashed line indicates the interface between the electrode (x < 0) and the electrolyte (x > 0).

Figure 2 .
Figure 2. Calculated surface charge density of ions (σ ion , blue line, left axis) and differential capacitance, (C d , red line, right axis) with respect to the electrode potential (E WE ).The parameter set for Ag(110) in 100 mM NaClO 4 is used.The black dashed line corresponds to the potential of zero charge (PZC).σ ion , E WE , and E PZC represent the surface charge in the electrolyte (eq 23), the working electrode potential, and the PZC, respectively.

Figure 3 .
Figure 3.Comparison between the calculated (dotted lines) and experimental (solid lines) differential capacitance; (a) Ag(110) in NaF, (b) Ag(110) in NaClO 4 , (c) Hg in NaF, and (d) Pt(111) in 100 mM KClO 4 (solid blue line: experimental, dotted blue line: calculated without specific interaction, dotted red line: calculated with specific interaction (eq 27)).Experimental results for (a, b) are extracted from ref 48, data for (c) are from ref 57, and data for (d) are from ref 4. The black dashed lines indicate the potential of zero charge (PZC).

Figure 4 .
Figure 4. Calculated (a) electron density distributions on Hg in 100 mM NaClO 4 and (b) Helmholtz capacitance (eq 26) of Hg (blue) and Ag(110) (red) electrodes in 100 mM NaClO 4 .For (a), the axis and line styles are the same as those in Figure 1b.The inset in (a) is the enlarged view of the domain inside the gray dashed box.The black dashed lines in (a, b) indicate the electrode−electrolyte interface and potential of zero charge (PZC), respectively.

Figure 5 .
Figure 5.Comparison of calculated differential capacitance with ω tot (k) : The different line styles represent different expressions for ω tot as shown in the legend.The parameter set for Ag(110) in 100 mM NaF was used as the base parameter.

Figure 6 .
Figure 6.Sensitivity (eq 28) of EDL properties on the material properties.The different colors represent different material properties, as listed in the legend.The parameter set for Ag(110) in 100 mM NaF is used as the base parameters.E PZ is in the unit of mV.

Table 1 .
By integrating the charge density of ions in the electrolyte phase, one can obtain the surface charge density σ ion Input Properties of the Electrodes a sol +(23)when L sol and L metal are large, the total charge density σ tot := σ ion + σ metal becomes zero, where

Table 2 .
Input Properties of the Electrolyte r i /nm (ref 54) a β i /M −1 (ref 45) z i

Table 3 .
Expressions for ω tot(k)to Analyze the Effects of the Major Potential Components a

■
ACKNOWLEDGMENTS This work was partially supported by the Hydrogen and Fuel Cell Technologies Office (HFTO), Office of Energy Efficiency and Renewable Energy, U.S. Department of Energy (DOE) through the Million Mile Fuel Cell Truck (M2FCT) consortia, technology managers G. Kleen and D. Papageorgopoulos, and a CRADA with Toyota Central R&D Laboratories, Inc. Additionally, this work was supported by the DOE Office of Science Energy Earthshot Initiatives as part of the Center for Ionomerbased Water Electrolysis at Lawrence Berkeley National Laboratory under contract number DE-AC02-05CH11231.Also, the authors are thankful to Dr. Ryosuke Jinnouchi in Toyota Central R&D Laboratories., Inc. for fruitful discussions on microscopic interactions in the double-layer structure.