Mesoscopic theory for a double layer capacitance in concentrated ionic systems

Abstract

The effect of oscillatory decay in charge density in concentrated ionic solutions and ionic liquids on the double-layer capacitance is studied within the framework of a mesoscopic theory. Only Coulomb and steric forces between the ions that are present in all ionic systems are taken into account. We show that the charge oscillations lead to a rescaled distance between the electrode and the virtual monolayer of counterions in the Helmholtz capacitance, and the scaling factor depends on the period of the charge oscillations. Our very simple formula for large density of ions and small voltages can serve as a reference point for the double layer capacitance in concentrated ionic solutions and ionic liquids, and can help to disentangle the universal and specific contributions to the capacitance in particular systems.

---

I. Introduction

Understanding the fundamental properties of a double-layer in concentrated electrolytes and room temperature ionic liquids is of great importance for designing fuel cells, batteries, supercapacitors and energy storage devices. In the early model of the double layer, the counterions form a layer separated from the electrode by a distance L that was first approximated by the ionic diameter a. In this model, the Helmholtz capacitance is simply C_H = ε/(4πa), where ε is the dielectric constant. In a more accurate model of the double layer,¹ thermal motion of the ions is taken into account, and in the screening cloud of ions the charge density decays with the distance from the electrode with the Debye screening length λ_D.^1,2 In the Debye capacitance, the diffuse layer of counterions is taken into account, and L = λ_D. The Debye capacitance C_D = ε/(4πλ_D) is a reference point for dilute electrolytes at very small voltages.

In concentrated electrolytes and ionic liquids (IL), however, the charge density decays with the distance from the electrode in an oscillatory way, and the decay length differs significantly from λ_D. The simple picture of dilute electrolytes where the average distance between the ions is much larger than their diameter is no longer valid, and the assumption of point-like ions is not justified.

A lot of effort has gone into experimental, theoretical and simulation studies of the structure of concentrated electrolytes and IL.^3–22 On the one hand, universal behavior was observed by the Perkin group for the decay length λ_s of the force between charged objects immersed in concentrated electrolytes or IL.^4–6,8λ_s determined in these experiments obeys the scaling λ_s/λ_D ∼ (a/λ_D)ⁿ with n = 3 for all ionic systems confined between crossed mica cylinders.^4–6,8 This result indicates increasing charge–charge correlation length for increasing concentration of ions. In theoretical and simulation studies of concentrated ionic solutions as well as in recent SAXS experiments,^10–23 the increasing λ_s for increasing concentration of ions was confirmed, but different values of n were obtained in different works. The underscreening observed in the above experimental, simulation and theoretical studies was not confirmed by the AFM experiments with ionic systems confined between silica surfaces,²⁴ and the issue remains controversial.

On the other hand, specific effects play an important role in determining the capacitance, especially in the case of polar solvents, where ε may exhibit strong dependence on ρ and on the distance from the electrode.^25–28 In general, the performance of the double layer capacitance is determined by a combination of many factors that include the concentration of ions in the electrolyte solution, size and nature of the ions, solvent polarity, electrode material, electrolyte–surface interactions, pore geometry, temperature, etc.^28–41 In particular, electrode surface morphology modifies the structure of the electric double layer and can either break or enhance overscreening depending on the relationship of the surface roughness to the electrostatic correlation length and the ion size-asymmetry.^42,43

It is not easy to disentangle different effects on the capacitance. In this work, we consider the simplest model of ionic systems with the size of the ions taken into account, namely the restricted primitive model (RPM) of charged hard spheres with equal diameter a and opposite charges in a structureless solvent characterized by the dielectric constant ε. In this model, the Coulomb and steric interactions lead to the oscillatory decay of the charge density with the distance from a charged boundary.^44–46 We can thus determine the effect of the charge layering near the electrode on the capacitance that should be common for many concentrated ionic solutions. In particular cases, specific interactions and the dielectric constant dependence on the local environment can play important roles. Our purpose, however, is to find a reference point for the capacitance in concentrated ionic solutions at very small voltages, in analogy with the Debye capacitance for dilute electrolytes. To achieve this goal, we use the theoretical approach developed in our previous works for the description of systems with spontaneous inhomogeneity in the bulk, near a charged wall and in a slit geometry.^{17,22,45,47,48}

In Section II, we briefly summarize the mesoscopic theory developed in ref. 17,22,45,47 and 48. In Section III.A, we present results for dilute electrolytes to test the mesoscopic theory predictions for the capacitance. In Section III.B, we derive our results for the capacitance of systems with large density of ions. We discuss our results in Section IV, and conclude in Section V.

II. The mesoscopic theory

In this section, we briefly summarize the key concepts, assumptions, approximations and results of the mesoscopic theory that can be applied to dilute as well as to concentrated electrolytes. The theory developed and described in detail in previous studies^17,22,47,48 allows the determination of the differential capacitance of the double layer in terms of the structure of the ionic solution. In principle, different levels of approximation in this theory are possible, but in this work, we limit ourselves to the simplest approximation to highlight the underlying effect of the structure, i.e., the distribution of the ions, on the capacitance.

In order to calculate the capacitance of the double layer, we consider the electrolyte in contact with a planar metallic electrode, and assume that the charge of the electrode is distributed over its flat surface. In the case of concentrated ionic solutions, for example water in the salt electrolyte, IL, or IL mixture with a neutral solvent, the assumption of point charges is not valid. The simplest model that takes the size of the ions into account is the restricted primitive model (RPM) of charged hard spheres with equal diameter a and opposite charges in a structureless solvent characterized by the dielectric constant ε. If the sizes of the positive and negative charge ions, a₊ and a₋, are somewhat different, we assume a = (a₊ + a₋)/2. We adopt this model, and assume in addition monovalent ions with the charge homogeneously distributed over the whole volume of the ion. For the same valency n of the cations and the anions, all our results apply when e is replaced by ne, where e is the elementary charge. For different valences of the anions and the cations, the theory becomes more difficult because of the lack of symmetry.^49,50 We should mention that significantly different sizes of ionic cores together with specific interactions can lead to spontaneous formation of relatively large charged domains,^51,52 as well as an asymmetric shape of the electric double layer capacitance.^53–55 The cases of significant size disparity, different valences of the cations and the anions and strong specific interactions require separate study.

Our theory⁵⁶ is based on the local mesoscopic volume fraction ζ_i of the i^th component of the mixture. The mesoscopic volume fraction is defined by analogy with its macroscopic counterpart, namely, ζ_i(r) is equal to the fraction of a mesoscopic volume with the center at r that is occupied by the particles of the considered species. In general, it depends on the scale of the coarse graining.

We assume that, near a flat electrode, ζ_i depends only on the distance z from the solid–liquid interface. As illustrated in Fig. 1, for the IL in contact with a flat solid surface, we identify the mesoscopic regions with layers of thickness a that are parallel to the solid surface. The mesoscopic volume fraction of the anions or the cations at the center of the layer is equal to the fraction of the volume of that layer that is occupied by the anions or the cations, respectively. With this definition, we obtain continuous functions of the distance from the electrode. Dimensionless mesoscopic densities are defined by ρ_i(z) = 6ζ_i(z)/π, and in the ionic system, i = +, −. It is convenient to introduce the local dimensionless number density of ions, ρ(z) = ρ₊(z) + ρ₋(z) and the local dimensionless charge, c(z) = ρ₊(z) − ρ₋(z). For monovalent ions, the local charge density is ec(z), where e is the elementary charge.

Fig. 1 Cartoon explaining the mesoscopic densities in our theory in the case of IL. The yellow shaded region represents the solid wall with fixed charge or fixed voltage. The red and gray circles represent the negatively and positively charged ions. The mesoscopic density is ρ_i(z) = 6ζ_i(z)/π, where ζ_i(z) is the fraction of the volume of the layer centered at z and of the thickness a that is covered by the i^th type ions. A typical layer with the center at z = z₀ is bounded by the dashed lines. With this construction, ρ_i(z) are continuous functions. The charge at the surface of the solid wall is included in the considered system. We choose z = 0 at the distance a/2 from the surface of the electrode inside the solid, as shown in the cartoon. With this choice, the mesoscopic charge density c(0) = ρ₊(0) − ρ₋(0) contains only the charge from the wall, i.e. the charge included in the liquid does not contribute to c(0). When the charge of the metallic electrode is confined to the surface at z = a/2, a layer of the thickness a with the center at z < 0 contains no charge; therefore, c(z < 0) = 0. With this formulation of the mesoscopic theory, we assume that the charge in the electrode of area A is equal to ec(0)A, and the overall charge neutrality condition is .

There are some ambiguities in defining the mathematical surface representing the solid–liquid interface at the microscopic and mesoscopic levels. We choose the origin of the coordinate frame, z = 0, inside the solid at the distance a/2 from the solid surface at which the electrode charge with the surface charge density eσ₀ is homogeneously distributed (see Fig. 1). From Fig. 1, it can be clearly seen that in our mesoscopic approach, c(0) is equal to the dimensionless surface charge density . For z < 0, we have c(z) = 0 when the charge of the electrode is confined to the surface of the solid, since this surface is now outside the layer with the center at z < 0. For 0 < z < a, both the electrode and the liquid contribute to c(z). Finally, for z > a, the only contribution to the charge density ec(z) comes from the electrolyte. The charge-neutrality condition of the whole system in this theory takes the simple form

The electrostatic potential at z = 0 is given by the Poisson equation in the integral form,

and can be calculated once the shape of c(z) is known. Here and below, the distance is in a units. All the characteristic lengths will be in a units as well, i.e. we will consider dimensionless quantities.

Determination of the equilibrium shape of c(z) is the main difficulty of the theory. From thermodynamics, we know that in a system with fixed volume and fixed temperature T that is in contact with a bulk reservoir of ions, the grand thermodynamic potential Ω takes a minimum. Thus, we should consider the grand potential for different forms of ρ(z) and c(z), and find the functions that minimize the functional

Ω[c,ρ] = U[c] − TS[c,ρ] − μN, (3)

where μ and N are the chemical potential and the number of the ions, respectively. The internal energy U[c] depends only on the charge c if the specific interactions are neglected, and takes in k_BT = 1/β units, the formwhere l_B = βe²/aε is the Bjerrum length in a units, r = |r₁ − r₂|, and g(r) is the pair distribution function. The entropy S[c,ρ] consists of the entropy of mixing of the anions and the cations, and of a contribution associated with packing of the hard spheres representing the ionic cores.

The exact expressions for Ω[c,ρ] (including the precise form of g(r) and S[c,ρ]) are not known, and different approximate theories were developed. In the bulk, the position independent ρ is a function of the chemical potential and temperature. The average charge density is c = 0 because it is equally probable to find an anion or a cation in a given microscopic volume in the absence of external fields. If an ion is kept at a given position, however, then it is more probable to find an oppositely charged ion in its vicinity than an ion of the same charge because the energy in the former case is lower. The charge–charge correlation function is a result of the competition between the energy favoring oppositely charged close neighbors and the entropy favoring the random distribution of the ions in the whole volume.

In our mesoscopic theory,^17,56 we assume that the ions cannot overlap; therefore, g(r) = 0 for r<1 (r is in a units), and g(r) → 1 for r → ∞, since at very large distances the ions are not correlated. In the mean-field (MF) approximation, the correlations are neglected for r > 1, and g(r) ≈ θ(r − 1), where θ(r − 1) = 1 for r > 1 and θ(r − 1) = 0 for r < 1. In this approximation,

whereis the energy per unit amplitude of a charge wave with the wavelength 2π/k excited in the homogeneous system, and ĉ(k) is the amplitude of this wave. V̂_C(k) < 0 and the energy decreases when the charge wave with k > π/2 is excited. It takes a minimum for k₀ ≈ 2.46 (in 1/a units), consistent with energetically favorable oppositely charged close neighbors.

The charge–charge correlation function in Fourier representation takes in this theory the form

where ρ_R = ρ in the MF approximation and for ρ ≪ 1. In concentrated solutions, however, ρ_R < ρ. The renormalized density of ions in the mesoscopic theory follows from the energetically favorable charge waves that play a similar role as the neutral clusters observed in concentrated electrolytes,^16,57,58i.e., they lead to a smaller density of free ions. Equations for ρ_R are developed and discussed in ref.17 and 56.

The charge–charge correlations in real space are obtained by the inverse Fourier transform of 〈ĉ(k)ĉ(−k)〉, and the decay lengths are obtained from simple poles of 〈ĉ(k)ĉ(−k)〉 extended to the complex q plane. In general, the imaginary pole q = ia_i gives a monotonic decay with the decay length 1/a_i, and the complex poles q = iα₀ ± α₁ give an oscillatory decay with the decay length 1/α₀ and the period 2π/α₁.

In the presence of the electrode, the excess density over the bulk value, Δρ(z) = ρ(z) − ρ_b, and c(z) ≠ 0 minimizes the functional ΔΩ[c,Δρ] = Ω[c,ρ_b + Δρ] − Ω_bulk[0,ρ_b]. When the bulk is in thermal equilibrium, the terms linear in Δρ and c in Taylor expanded ΔΩ[c,Δρ] vanish, and ΔΩ[c,Δρ] contains terms of the second and higher orders in the fields,

where the parameters R₀ and R₂ come from the entropic contribution to Ω and from the correlations between fluctuations of the local charge.^22,48 The differential operator is defined through the Taylor expansion, and the last term in eqn (9) is given bywhere βf_h is the free energy density of a mixture of hard spheres in k_BT units. Minimization of the functional (8) leads to the ordinary differential Euler–Lagrange (EL) equationsthat should be solved under boundary conditions (fixed charge or fixed voltage). Note that the above equations are nonlinear. In addition, the second terms on the LHS lead to a coupling between the charge and the number density profiles. This makes the problem technically difficult. However, g_h(ρ_b,c,Δρ) is of a higher order in the fields c and Δρ; therefore, ∂g_h/∂c becomes a negligible correction for a sufficiently small charge, i.e. for small voltages of the electrode. Unfortunately, unless the full set of equations is solved, it is not possible to accurately determine the range of the voltage for which the derivatives of g_h in the above equations can be neglected. We can expect that the voltage is certainly too strong and the derivatives of g_h have to be included when the obtained dimensionless charge is c(z) > 1.

When the derivatives of βg_h can be neglected, the EL equations for Δρ and c are linear and decoupled. The linearized EL equation for c(z) is

A solution of a linear equation is a sum of exponential terms ∝ exp(λz). In our case, λ = −iq, where complex q is a solution of the equation . For our purpose here, it is important that the linearized EL equation leads to c(z) decaying at large distances with the same decay lengths as the charge–charge correlation function (compare eqn (12) and (7)). The solution of eqn (12) must satisfy the charge neutrality condition (1).

Eqn (12), (1), (2) and c(0) = a²σ₀ allow for the calculation of the differential capacitance in the case of small voltages, and with the specific effects of the solvent neglected.

III. Results

III.A Capacitance of dilute electrolytes

In our mesoscopic theory, there are two inverse decay lengths, a₁ > a₂ for large T and small ρ, because there are two solutions of the equation (see eqn (7))and the asymptotic decay of correlations at large distances is monotonic. For ρ → 0, we get , where is the dimensionless Debye screening length (recall that l_B and ρ are dimensionless), and a₁ is of order of unity (in 1/a units). a₁ and a₂ merge at the so-called Kirkwood line on the (ρ,T) diagram,⁵⁹ and at lower T they transform to a pair of complex inverse decay lengths that lead together to an oscillatory decay of the charge–charge correlations. This result agrees with previous theories and simulations,^60–62 and is consistent with dominant roles of the entropy and energy at high and low T, respectively.

Because the decay lengths of the charge density are the same as the correlation lengths, the charge-density profile satisfying the charge-neutrality condition has on the high-T small ρ side of the Kirkwood line the form

where the length unit is the ion diameter a. The dimensionless inverse decay lengths are determined numerically based on ref. 17 and 59. Similar values were obtained in different theories.⁶¹ A representative c(z) is shown in Fig. 2, for , l_B = 2 and ρ = 0.01.

Fig. 2 The dimensionless charge density in the case of dilute electrolytes with ρ = 0.01, l_B = 2 and dimensionless surface charge density as obtained in the mesoscopic theory. The decay lengths satisfy the equation 4πl_Bρ_R cosh a_i = a_i², and ρ_R ≤ ρ satisfies the equations presented in ref. 17. z is in units of the ion diameter a.

The physical meaning of the two decay lengths in the context of the mesoscopic theory can be understood by comparing Fig. 2 for the representative charge density and Fig. 1 illustrating the construction of the theory. For large z, the decay of the charge is monotonic and the sign of c(z) is opposite to the sign of the electrode. The decay length is for ρ → 0, as discussed above. By construction of the mesoscopic theory, however, for z increasing from z = 0 to z = 1, there is a contribution to c(z) from the charge of the electrode and an increasing contribution from the opposite charge of the ions in the electrolyte (see Fig. 1 for 0 < z < 1). Thus, c(z) must change sign for z < 1. The two terms in eqn (14), one with the same sign as the sign of the charge of the electrode and small decay length 1/a₁ and the other one with the opposite sign and large decay length 1/a₂, are consistent with the above physical picture.

The surface-charge density in the mesoscopic theory is

For the potential at z = 0, we obtain from eqn (14) and (2)

The capacitance C = d(eσ₀)/dU is easily obtained from eqn (15) and (16) and is given by

For dilute electrolytes, a₂/a → 1/λ_D that can be easily seen from the equation l_BV̂_C(ia₂) + 1/ρ = 0 for the imaginary simple pole and V̂_C(ia₂) → −4π/a₂² for ρ ≪ 1 (see eqn (7) and (6)), giving a₂² = 4πl_Bρ. Thus, for dilute electrolytes, we obtain

where the dimensionless parameter a₁ in our theory is larger than 2 and increases with increasing (l_Bρ)⁻¹, but remains of order of unity. The precise value of this inverse microscopic length should be determined in a more exact microscopic theory. Thus, in the limit of dilute electrolytes, we obtain the Debye capacitance up to a parameter a₁ = O(1).

In order to understand the origin of the difference between eqn (18) and the Debye capacitance, recall that the Debye capacitance was obtained for point charges. The point charges can be at the distance z = 0 from the electrode surface, and c(z) is monotonic for the whole range of z > 0. If the size a of the ions is taken into account, microscopic details such as the charge distribution over the volume of the ion, the definition of the position of the electrode surface, etc. start to play a role, and some ambiguities in the definition of the capacitance in terms of the charge distribution appear. In our mesoscopic theory, the charge density c(z) is averaged over the layers of the thickness a (Fig. 1), and as discussed above, the charge profile is nonmonotonic and changes sign for z ∼ 1. This leads to a smaller value of Ψ(0) defined in eqn (2), and as a result to a larger value of the capacitance, with the coefficient a₁ associated with the charge distribution in the close neighborhood of the electrode surface. To compare our predictions with the results of simulations or experiments on the quantitative level, we should take into account that formulas (17) and (18) contain the factor a₁ that in the theory with the microscopic structure averaged over the region with the linear size a is 2 ≤ a₁ ≤ 6 rather than a₁ = 1 present in the Debye capacitance for point charges.

III.B Capacitance of concentrated electrolytes

The solution of the EL eqn (12) at the low-T side of the Kirkwood line has the form

c(z) = A_ce^−α₀z sin(α₁z + θ). (19)

The above dimensionless charge density profile agrees with the fixed surface charge density σ₀ and satisfies the charge-neutrality conditions (1) when

A typical dimensionless charge density near a weakly charged electrode is shown in Fig. 3. Predictions of the mesoscopic theory for the charge profile⁴⁵ agree with simulation results of the RPM model in the case of large density of ions.⁴⁴ Namely, the electrostatic potential βΨ(z) was very similar in the simulations and in our mesoscopic theory for two values of the surface charge studied in the simulations and for the same thermodynamic states and parameters of the model. In addition, as shown for example in ref. 46, eqn (19) perfectly fits the charge density in atomistic simulations of IL – alcohol mixture for z > 2π/α₁. The atomistic simulations give different shapes of c(z) for z < 2π/α₁ because the microscopic charge density obtained in the simulations is compared with the charge density averaged over the layers of the thickness a in the mesoscopic theory.

Fig. 3 The dimensionless charge density in the case of the concentrated electrolyte with ρ = 0.4, l_B = 2 and ρ_R ≈ 0.185 (left) and ρ = 0.7, l_B = 4 and ρ_R ≈ 0.154 (right) for dimensionless surface charge density , as obtained in the mesoscopic theory. α₀ and α₁ satisfy the equations (α₁² − α₀²) + 4πl_Bρ_R cosh α₀ cos α₁ = 0 and α₁α₀ − 2πl_Bρ_R sinh α₀ sin α₁ = 0.¹⁷z is in units of the ion diameter a.

The electrostatic potential (2) at z = 0 for the charge density given by eqn (19) and (20) is

The capacitance can be easily calculated using eσ₀ = eA_c sin θ/a² with eqn (20) and (21), and the result is

where ε = 4πε_rε₀, with ε_r and ε₀ ≈ 9 × 10⁻⁶ μF m⁻¹ denoting the relative dielectric constant and the vacuum permittivity, respectively.

The dimensionless wavenumber α₁ ≤ k₀ ≈ 2.46 of the damped charge oscillations decreases slightly for increasing (l_Bρ)⁻¹ in the concentrated electrolyte.⁵⁹ As found recently^4,6,56 in concentrated electrolytes, the dimensionless inverse decay length is α₀ ≪ 1. In the lowest order approximation, we keep only the terms quadratic in Δρ and c, and follow the steps described in ref. 45 and 48. We can assume that in concentrated electrolytes and IL α₀² ≪ α₁², and

Formulas (22) and (23) relating the capacitance with the period 2πa/α₁ of the damped charge oscillations near the electrode are the main result of this work. We should emphasize, however, that the smoothed shape of the mesoscopic c(z) near the electrode leads to a smaller value of the electrostatic potential at z = 0, hence to a larger value of the capacitance, as already discussed in detail in Section III.A.

IV. Discussion

The expressions for the capacitance of dilute and concentrated electrolytes in the RPM near a flat metallic electrode, eqn (17) and (22), are significantly different. They become identical, however, at the Kirkwood line separating the monotonic and oscillatory asymptotic decays of the charge density because at the Kirkwood line a₁ = a₂ = α₀ and α₁ = 0. Thus, we obtained a continuous function for the whole range of the density of ions. Comparison of C with C_H = ε_rε₀/a and C_D = ε/(4πλ_D) for the density of ions 0 < ρ < 0.7 and for fixed Bjerrum lengths l_B = 2 and l_B = 4 is shown in Fig. 4.

Fig. 4 The capacitance C/C_H obtained from (17) and (22) (solid lines) and C_D/C_H (dashed lines) as a function of the dimensionless density of ions for the Bjerrum length l_B = 2 and l_B = 4 in a-units.

Our general analytical formulas (17) and (22) were obtained from the linearized EL eqn (12); therefore, they are valid only for very small voltages. Nevertheless, they highlight the effect of the charge distribution on the capacitance for the whole range of the density of ions on a general qualitative level. For large voltages, the charge density is no longer small, and the coupled nonlinear EL equations for c and Δρ (see eqn (10) and (11) and ref. 56) have to be solved to determine C. This is possible only numerically for particular cases, and will be a subject of our future study.

The formulas for the capacitance become particularly simple for very dilute and very dense electrolytes. In the former case, we obtain the well-known Debye capacitance C_D = ε/(4πλ_D) up to a dimensionless coefficient of order unity (see eqn (18)). For the very concentrated electrolyte or IL, we obtain formula (23) that is strikingly similar to the Helmholtz capacitance C_H(L) = ε/(4πL) in the early model of the double layer, where L is the distance between the electrode and the surface occupied by the counterions. We take into account the whole oscillatory charge profile such as the ones shown in Fig. 3, and find that the alternating oppositely charged layers have the same effect as a single layer of counterions located at the distance L = a/α₁² from the electrode. This shows that the simplest model of the double layer works well even in the case of a rather complex structure, but with the distance of the virtual monolayer of counterions from the electrode, L, depending on the dimensionless wavenumber α₁ ≈ 2.46 of the charge oscillations.

As shown in Fig. 4, C/C_H increases quickly for increasing ρ < 0.15, and slowly for ρ > 0.3, whereas C_D/C_H increases with ρ gradually, with almost constant slope. For the two considered values of l_B, C/C_H > 5 if ρ > 0.3. Moreover, C > C_D, but C − C_D decreases with increasing ρ, and for l_B = 4, C ≈ C_D for ρ ≈ 0.7. Let us compare predictions of the mesoscopic theory for the charge-density profiles shown in Fig. 3 with the classical Debye capacitance in more detail. For ρ = 0.4 and l_B = 2, and the diameter of hydrated ions a ≈ 0.5 nm corresponding to ∼2.65 M NaCl_aq, we get α₀ ≈ 0.93 and α₁ ≈ 2.18, and obtain from eqn (22)

C ≈ 10ε_r μF cm⁻². (24)

From , we get in this case λ_D ≈ 0.16 nm, and the formula valid for dilute electrolytes gives

C_D ≈ 5.6ε_r μF cm⁻². (25)

The Debye length in this case differs from the physically relevant lengths a/α₀ ≈ 0.54 nm and 2πa/α₁ ≈ 1.44 nm, however. In another example shown in Fig. 3 with ρ = 0.7 and l_B = 4, we obtain assuming a = 0.9 nm

C ≈ 6ε_r μF cm⁻². (26)

The Debye length is λ_D ≈ 0.15 nm, and

C_D ≈ 5.9ε_r μF cm⁻². (27)

The values of C and C_D are very similar in this case, but the interpretation is quite different. For the oscillatory decay of the charge density, the key factor is the period of the charge oscillations near the electrode, and the Debye length is not associated with characteristic lengths of the charge distribution. These examples show that care must be taken in interpreting experiments and simulations because correct numbers can follow from incorrect formulas.

In order to verify the accuracy of C given in eqn (22) and (23) on the quantitative level, we should compare the theoretical and simulation results for the same model. In ref. 44, C was obtained by simulations of the RPM with a = 1 nm, T = 450 K, ε_r = 2 and ρ ≈ 0.6 nm⁻³, and the result of simulations was C/C_D ≈ 0.15. For the above parameters, we obtain λ_D ≈ 0.084 nm and eqn (23) gives C/C_D = α₁²λ_D/a ≈ 0.5, where we used α₁ = 2.45. Recall that the capacitance obtained in the mesoscopic theory for dilute electrolytes was overestimated by the factor 2 ≤ a₁ ≤ 6. In the case of ρ = 0.6 nm⁻³ and large l_B, the theoretical result is about three times larger than the simulation result, i.e. a systematic overestimation of the capacitance by a factor ∼3 is present in the mesoscopic theory for the whole density range.

Our results show the effect on the capacitance of the charge ordering in the region extending to large distances from the electrode. On the quantitative level, however, the capacitance also depends on the details of the microscopic structure in the vicinity of the electrode that should be determined within a more exact microscopic theory. In our mesoscopic theory, the effect of the microscopic structure can be taken into account by additional scaling factor that based on the comparison with simulations is about 1/3. As already discussed at the end of Section III.A, the larger value of the capacitance in our theory follows from the smaller value of the electrostatic potential at z = 0 that in turn is a result of the smoothed shape of the mesoscopic c(z) near the electrode.

Let us discuss consequences of eqn (23) on a general level. The capacitance decreases with increasing size of the ions, in agreement with experimental results for aqueous ionic solutions and IL.^27,63 The dimensionless period of the charge wave in concentrated electrolytes depends rather weakly on density. According to our mesoscopic theory, in IL or highly concentrated electrolytes, the dimensionless α₁ = 2πa/λ_c is 2 ≤ α₁ ≤ 2.46, corresponding to the wavelength of the charge-density wave 2.55a ≤ λ_c ≤ 3.14a. Hence, in IL, we have for eqn (23) the approximation

In the RPM, the dependence of ε_r (and in turn of l_B) on the density of ions and on the distance from the electrode is neglected, whereas in different solvents, especially in water, this dependence can be quite strong. In aqueous solutions, ε_r decreases from about 80 in pure water to about 40 for 5 M solution of NaCl.⁴ In the Stern layer, the orientations of dipoles of water molecules in the hydration shells of ions are almost fixed, and the dielectric constant may decrease to ε_r ∼ 5 or even less.^55,64–66 Thus, quantitative predictions for the capacitance in particular cases are not possible within the RPM, especially for polar solvents such as water. Our results show, however, the general relationship between the capacitance and the period of the damped charge oscillations. The complex charge distribution can be replaced by the simplest model of the double layer, provided that the virtual single layer of counterions is separated from the electrode by the distance equal to the diameter of the ions re-scaled by the coefficient proportional to α₁⁻², where α₁ is the wavenumber of the damped charge oscillations in 1/a units and the proportionality constant is ∼3.

Detailed comparison of our prediction with results of simulations and experiments for particular systems goes beyond the scope of this work because it would be necessary to disentangle the universal properties captured by the RPM and the specific properties such as the size and nature of the ions, polarity of the solvent, the charge distribution and the roughness of the electrode's surface. Such an analysis should be done in future studies.

V. Conclusions

Our goal was to determine on a very general level the effect of charge ordering in concentrated electrolytes and IL on the capacitance of the double layer. We limited ourselves to the restricted primitive model, RPM, where spherical ions with equal diameters and opposite charges are dissolved in a structureless solvent characterized by the dielectric constant ε. This way we can determine the effect of the Coulomb and steric interactions in the absence of specific effects that differ from one system to the other. We obtained very simple expressions for the capacitance in dilute and concentrated electrolytes in the framework of the same mesoscopic theory. The main conclusion is that the simplest early model of the double layer works surprisingly well in the case of large density of ions, provided that the distance between the virtual single layer of counterions and the electrode is equal to the ion diameter re-scaled by a coefficient determined by the period of the damped charge oscillations. This conclusion agrees with recent simulation results.⁶⁶ Our formulas (22) and (23) should play for concentrated ionic systems a similar role as the Debye capacitance plays for dilute electrolytes. These equations can serve as a reference point that allows to disentangle universal and specific features of the capacitance.

Finally, it is worth mentioning that since our theory is suitable for describing the structural properties of concentrated electrolytes and their effect on capacitance, it is worthwhile to extend it to nonequilibrium properties, because in energy storage devices, the charging/discharging dynamics plays an important role.^67,68 The extension can be performed by analogy with the DFT extension to the DDFT.⁶⁹

Author contributions

Conceptualization: AC; formal analysis: OP and AC; investigation: OP and AC; methodology: OP and AC; visualization: OP and AC; writing – original draft: AC and OP; writing – review and editing: AC and OP.

Data availability

All formulas derived in this work are included in the article. The figures were obtained from these formulas. No more data were used.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We gratefully acknowledge the financial support from the European Union Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 734276 (CONIN).

References

1. P. Debye and E. Hückel, Phys. Z., 1923, 24, 185–206.
2. J. L. Barrat and J.-P. Hansen, Basic Concepts for Simple and Complex Liquids, Cambridge University Press, 2003.
3. M. A. Gebbie, M. Valtiner, X. Banquy, E. T. Fox, W. A. Henderson and J. N. Israelachvili, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 9674–9679.
4. A. M. Smith, A. A. Lee and S. Perkin, J. Phys. Chem. Lett., 2016, 7, 2157.
5. A. Lee, C. S. Perez-Martinez, A. M. Smith and S. Perkin, Phys. Rev. Lett., 2017, 119, 026002.
6. A. A. Lee, C. S. Perez-Martinez, A. M. Smith and S. Perkin, Faraday Discuss., 2017, 199, 239–259.
7. R. Kjellander, J. Chem. Phys., 2018, 148, 193701.
8. T. Groves and S. Perkin, Faraday Discuss., 2024, 253, 193–211.
9. F. Coupette, A. A. Lee and A. Härtel, Phys. Rev. Lett., 2018, 121, 075501.
10. B. Rotenberg, O. Bernard and J.-P. Hansen, J. Phys.: Condens. Matter, 2018, 30, 054005.
11. R. M. Adar, S. A. Safran, H. Diamant and D. Andelman, Phys. Rev. E, 2019, 100, 042615.
12. S. W. Coles, C. Park, R. Nikam, M. Kanduč, J. Dzubiella and B. Rotenberg, J. Phys. Chem. B, 2020, 124, 1778–1786.
13. J. Zeman, S. Kondrat and C. Holm, Chem. Commun., 2020, 56, 15635–15638.
14. J. Zeman, S. Kondrat and C. Holm, J. Chem. Phys., 2021, 155, 204501.
15. P. Cats, R. Evans, A. Härtel and R. van Roij, J. Chem. Phys., 2021, 154, 124504.
16. S. Safran and P. Pincus, Soft Matter, 2023, 19, 7907–7911.
17. A. Ciach and O. Patsahan, J. Mol. Liq., 2023, 377, 121453.
18. J. Yang, S. Kondrat, C. Lian, H. Liu, A. Schlaich and C. Holm, Phys. Rev. Lett., 2023, 131, 118201.
19. A. Härtel, M. Bültmann and F. Coupette, Phys. Rev. Lett., 2023, 130, 108202.
20. S. Wang, H. Tao, J. Yang, J. Cheng, H. Liu and C. Lian, J. Phys. Chem. Lett., 2024, 33, 7147–7153.
21. G. R. Elliott, K. P. Gregory, H. Robertson, V. S. Craig, G. B. Webber, E. J. Wanless and A. J. Page, Chem. Phys. Lett., 2024, 843, 141190.
22. A. Ciach and O. Patsahan, J. Phys.: Condens. Matter, 2021, 33, 37LT01.
23. M. Dinpajooh, E. Biasin, E. T. Nienhuis, S. T. Mergelsberg, C. J. Benmore, G. K. Schenter, J. L. Fulton, S. M. Kathmann and C. J. Mundy, J. Chem. Phys., 2024, 161, 151102.
24. S. Kumar, P. Cats, M. B. Alotaibi, S. C. Ayirala, A. A. Yousef, R. van Roij, I. Siretanu and F. Mugele, J. Colloid Interface Sci., 2022, 622, 819–827.
25. L. Fumagalli, A. Esfandiar, R. Fabregas, S. Hu, P. Ares, A. Janardanan, Q. Yang, B. Radha, T. Taniguchi, K. Watanabe, G. Gomila, K. S. Novoselov and A. K. Geim, Science, 2018, 360, 1339–1342.
26. M. H. Motevaselian and N. R. Aluru, ACS Nano, 2020, 14, 12761–12770.
27. A. D’Alessandro, S. Bellani, A. Gamberini, V. Mastronardi, M. I. Zappia, M. Abruzzese, S. Thorat, E. Calcagno and F. Bonaccorso, Batteries Supercaps, 2024, 7, e202300458.
28. E. Gongadze and A. Iglič, Electrochim. Acta, 2015, 178, 541–545.
29. O. K. Coskun, M. Muñoz, S. Dongare, W. Dean and B. E. Gurkan, Langmuir, 2024, 40, 3283–3300.
30. A. Kasina, E. Cocklin, S. Horswell, Y. Grunder and C. A. Lucas, J. Phys. Chem. C, 2024, 128, 13318–13332.
31. K. Wippermann, Y. Suo, C. Rodenbücher, C. Korte and A. A. Kornyshev, Electrochim. Acta, 2023, 469, 143207.
32. K. Zhang, C. Wei, M. Zheng, J. Huang and G. Zhou, Molecules, 2024, 29, 1246.
33. R. Costa, I. V. Voroshylova, M. N. D. Cordeiro, C. M. Pereira and A. F. Silva, Electrochim. Acta, 2018, 261, 214–220.
34. A. A. Kornyshev, J. Phys. Chem. B, 2007, 111, 5545–5557.
35. M. V. Fedorov and A. A. Kornyshev, Chem. Rev., 2014, 114, 2978–3036.
36. J. P. de Souza, Z. A. Goodwin, M. McEldrew, A. A. Kornyshev and M. Z. Bazant, Phys. Rev. Lett., 2020, 125, 116001.
37. J. Wu, Chem. Rev., 2022, 122, 10821–10859.
38. P. Cats, R. S. Sitlapersad, W. K. den Otter, A. R. Thornton and R. van Roij, J. Solution Chem., 2022, 51, 296–319.
39. G. Jeanmairet, B. Rotenberg and M. Salanne, Chem. Rev., 2022, 122, 10860–10898.
40. E. Gongadze and A. Iglič, Bioelectrochemistry, 2012, 87, 199–203.
41. C. Cruz, A. Ciach, E. Lomba and S. Kondrat, J. Phys. Chem. C, 2019, 123, 1596–1601.
42. A. Khlyupin, I. Nesterova and K. Gerke, Electrochim. Acta, 2023, 450, 142261.
43. I. Nesterova, N. M. Evstigneev, O. I. Ryabkov, K. M. Gerke and A. Khlyupin, J. Colloid Interface Sci., 2025, 677, 396–405.
44. M. V. Fedorov and A. A. Kornyshev, Electrochim. Acta, 2008, 53, 6835–6840.
45. A. Ciach, J. Mol. Liq., 2018, 270, 138–144.
46. J. M. Otero-Mato, H. Montes-Campos, O. Cabeza, D. Diddens, A. Ciach, L. J. Gallego and L. M. Varela, Phys. Chem. Chem. Phys., 2018, 20, 30412–30427.
47. A. Ciach and O. Patsahan, Condens. Matter Phys., 2012, 15, 23604.
48. O. Patsahan, A. Meyra and A. Ciach, J. Mol. Liq., 2022, 363, 119844.
49. A. Ciach, W. T. Góźdź and G. Stell, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 051505.
50. O. Patsahan and A. Ciach, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 86, 031504.
51. O. Borodin, L. Suo, M. Gobet, X. Ren, F. Wang, A. Faraone, J. Peng, M. Olguin, M. Schroeder, M. S. Ding, E. Gobrogge, A. von Wald Cresce, S. Munoz, J. A. Dura, S. Greenbaum, C. Wang and K. Xu, ACS Nano, 2017, 11, 10462–10471.
52. O. Patsahan and A. Ciach, ACS Omega, 2022, 7, 6655–6664.
53. M. V. Fedorov and A. A. Kornyshev, J. Phys. Chem. B, 2008, 112, 11868–11872.
54. A. Velikonja, V. Kralj-Iglič and A. Iglič, Int. J. Electrochem. Sci., 2015, 10, 1–7.
55. A. Iglič, E. Gongadze and V. Kralj-Iglič, Acta Chim. Slov., 2019, 66, 534–541.
56. A. Ciach and O. Patsahan, J. Mol. Liq., 2023, 382, 121949.
57. Z. A. Goodwin and A. A. Kornyshev, Electrochem. Commun., 2017, 82, 129–133.
58. Z. A. M. Chen, Z. A. H. Goodwin, G. Feng and A. A. Kornyshev, J. Electroanal. Chem., 2018, 819, 347–358.
59. A. Ciach, W. T. Góźdź and R. Evans, J. Chem. Phys., 2003, 118, 3702–3710.
60. P. Attard, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 1993, 48, 3604.
61. R. L. de Carvalho and R. Evans, Mol. Phys., 1994, 83, 619–654.
62. G. M. Torrie and J. P. Valleau, J. Chem. Phys., 1980, 73, 5807–5816.
63. V. Lockett, M. Horne, R. Sedev, T. Rodopoulos and J. Ralston, Phys. Chem. Chem. Phys., 2010, 12, 12499–12512.
64. A. Velikonja, E. Gongadze, V. Kralj-Iglič and A. Iglič, Int. J. Electrochem. Sci., 2014, 9, 5885–5894.
65. A. V. Dubtsov, S. V. Pasechnik, D. V. Shmeliova, A. S. Saidgaziev, E. Gongadze, A. Iglič and S. Kralj, Soft Matter, 2018, 14, 9619–9630.
66. S. Park and J. G. McDaniel, J. Chem. Phys., 2024, 160, 164709.
67. S. Kondrat, G. Feng, F. Bresme, M. Urbakh and A. A. Kornyshev, Chem. Rev., 2023, 123, 6668–6715.
68. C. Noh and Y. Jung, Phys. Chem. Chem. Phys., 2019, 21, 6790–6800.
69. A. J. Archer, J. Chem. Phys., 2009, 130, 014509.

---