To capture the long-range part of the electrostatic interaction, we choose a truncated Coulomb potential u r = l B 2 r r r 0 + 1 − r r 0 − 1 , (5) where l B is the Bjerrum length and r 0 is the truncation length. Truncated Coulomb potentials have been employed in the past to cut off the long-range part with the goal to facilitate efficient computer simulations (Linse and Andersen, 1986; Baker et al., 1999). In contrast, we propose to cut off the short-range part of the Coulomb interaction, as shown in Figure 1: the black solid line in both diagrams of Figure 1 shows the truncated potential u ( r ) according to Equation 5. The operator A can be identified conveniently in Fourier space, where the Fourier transformation u ̃ k = 4 π k ∫ 0 ∞ d r r ⁡ sin k r u r = 4 π l B sin k r 0 k 3 r 0 , together with A ̃ u ̃ ( k ) = − 1 , gives rise to A = ∇ 2 4 π l B r 0 ∇ sinh r 0 ∇ = ∇ 2 4 π l B ∑ i = 0 ∞ 2 − 4 i B 2 i 2 i ! r 0 ∇ 2 i = ∇ 2 4 π l B 1 − 1 6 r 0 ∇ 2 + 7 360 r 0 ∇ 4 ∓ ⋯   . (6) In Equation 6, ∇ and B i denote the nabla differential operator and the ith Bernoulli number, respectively. Recall that the set of Bernoulli numbers can be defined via a series expansion of the generating function x / ( e x − 1 ) = ∑ i = 0 ∞ B i x i / i ! . If the truncated Coulomb potential is approximated by the j -th order operator, A j = ∇ 2 4 π l B ∑ i = 0 j 2 − 4 i B 2 i 2 i ! r 0 ∇ 2 i , (7) then the corresponding potential is given by the following equation: u j r = 2 π l B r ∫ 0 ∞ d k sin k r k 1 ∑ i = 0 j 2 − 4 i B 2 i 2 i ! − k 2 r 0 2 i . (8) The function u j ( r ) in Equation 8 is the j -th order approximation of the truncated Coulomb potential. For j = 0 , Equation 8 gives rise to the Coulomb potential u 0 ( r ) = l B / r . In the limit j → ∞ , the truncated Coulomb potential u j → ∞ ( r ) = u ( r ) , as specified in Equation 5, is recovered. As has been pointed out by Lee et al. (2015), the first-order approximation of the truncated Coulomb potential is as follows: u 1 r = l B r 1 − e − 6 r r 0 . (9) The right diagram in Figure 1 shows the truncated Coulomb potential u ( r ) together with the increasingly more accurate approximations u 1 ( r ) , u 2 ( r ) , u 3 ( r ) , and u 4 ( r ) . In the present work, we will focus on the first-order approximation u 1 ( r ) only, leading to a fourth-order Poisson equation A 1 Ψ = − η / ν that utilizes the differential operator A 1 = ∇ 2 [ 1 − ( r 0 ∇ ) 2 / 6 ] / ( 4 π l B ) similar to the BSK model (Bazant et al., 2011). However, Equations 7, 8 provide a method for the incorporation of higher-order approximations in a straightforward manner.

We point out that the BSK model (Bazant et al., 2011) is an effective mean-field approach based on the pair potential u 1 ( r ) . The emergence of u 1 ( r ) from the bare Coulomb potential can be the result of charge discreteness, which we consider in the present work, or of short-range correlations. Our present work considers electrostatic correlations beyond the short-range correlations of the BSK model. Electrostatic correlations between ion pairs are not accounted for by the BSK model.

Before formulating and minimizing a free energy that accounts for charge discreteness and electrostatic correlations between ion pairs, we start with three comments: first, choosing the truncation radius r 0 to be equal or slightly larger than the lattice spacing b enables us to explicitly include the electrostatic interactions between neighboring ions: ω 11 = ω 22 = ω = l B / b between cation–cation and anion–anion pairs and ω 12 = ω 21 = − ω = − l B / b between cation–anion and anion–cation pairs. This introduces the definition ω = l B / b and the notation ω i j that we use in Equation 10. Second, we account for electrostatic correlations between ions located at neighboring lattice sites through QCA (Davis, 1996). Correlations involving more than two ions can, in principle, be incorporated into QCA (Bossa et al., 2015), but they are ignored in the present work. To account for correlated ion pairs, we consider the local fraction ϕ i j = ϕ i j ( r ) of i - j pairs, where the indices i and j adopt the value 1 for a cation and 2 for an anion. Third, we introduce the chemical potentials μ 1 and μ 2 of the cations and anions, respectively, and choose them to ensure coexistence with a bulk phase of the ionic liquid, where ϕ 1 = ϕ 2 = 1 / 2 . We also recall the definition of the lattice coordination number z . With that, we express the free energy of a lattice-based ionic liquid as follows (Downing et al., 2018b): F = 1 ν ∫ d 3 r − ν 2 Ψ A Ψ + 1 − z ϕ 1 ⁡ ln ϕ 1 + ϕ 2 ⁡ ln ϕ 2 + z 2 ∑ i j = 1 2 ϕ i j ⁡ ln ϕ i j + ω i j ϕ i j − μ 1 ϕ 1 − μ 2 ϕ 2 . (10) The four contributions to this free energy account for the electrostatic energy due to the long-range part of the Coulomb potential, the mixing entropy of ion pairs according to QCA, the electrostatic energy associated with neighboring ion pairs, and a Legendre transformation that fixes the chemical potentials of the cations and anions in the ionic liquid. The ϕ i j in Equation 10 must fulfill the three conservation relations ϕ 1 = ϕ 11 + ϕ 12 , ϕ 2 = ϕ 21 + ϕ 22 , and ϕ 12 = ϕ 21 . This leaves one degree of freedom for the local distribution of ion pairs, which we define conveniently as ϕ ̄ = ϕ 12 + ϕ 21 . The difference in mole fractions η = ϕ 1 − ϕ 2 constitutes another degree of freedom. We thus have the following equation: ϕ 1 = 1 + η 2 , ϕ 2 = 1 − η 2 , ϕ 11 = ϕ 1 − ϕ ̄ 2 , ϕ 22 = ϕ 2 − ϕ ̄ 2 , ϕ 12 = ϕ 21 = ϕ ̄ 2 . (11) Note that both ϕ ̄ = ϕ ̄ ( r ) and η = η ( r ) depend on the position r inside the ionic liquid. Calculation of the first variation in the free energy F = F ( η , ϕ ̄ ) leads to δ F = 1 ν ∫ d 3 r δ η Ψ + 1 − z 2 ln 1 + η 1 − η + z 4 ln 1 + η − ϕ ̄ 1 − η − ϕ ̄ − 1 2 μ 1 − μ 2 + δ ϕ ̄ z 4 ln ϕ ̄ 2 1 + η − ϕ ̄ 1 − η − ϕ ̄ − z ω . (12) We have not included boundary terms in Equation 12 because they will be discussed separately below. By symmetry, the chemical potentials μ 1 = μ 2 are equal. Thermal equilibrium demands δ F = 0 , thus leading to the following two equations: e 4 ω = ϕ ̄ 2 1 − ϕ ̄ 2 − η 2 , 0 = 2 Ψ + 1 − z ln 1 + η 1 − η + z 2 ln 1 + η − ϕ ̄ 1 − η − ϕ ̄ , (13) which must be fulfilled at each position r . Equation 13 defines the two relations η = η ( Ψ ) and ϕ ̄ = ϕ ̄ ( Ψ ) , with ω and z as additional parameters. Spatial compositional variations arise from the solutions of Equation 3. However, to simplify calculations, we replace the operator A by its first-order approximation A 1 , as defined in Equation 7. This renders the level of modeling the long-range part of the electrostatic interactions similar to the BSK model (Bazant et al., 2011) yet with different boundary conditions as discussed below. Equation 3 then reads as follows: l 2 ∇ 2 1 − r 0 2 6 ∇ 2 Ψ = − η Ψ , (14) where we have defined the length l = ν / ( 4 π l B ) . For example, l B = 1 nm and ν = 1 nm 3 imply l = 0.3 nm . We will use l as unit length throughout this work. Equation 14 is a self-consistency relationship that takes the form of a fourth-order differential equation with the function η ( Ψ ) defined through a set of two algebraic equations. As pointed out above, electrostatic interactions are represented by the potential u 1 ( r ) in Equation 9. Unlike interpretations in terms of short-range ion correlations (Storey and Bazant, 2012; Moon et al., 2015; Blossey et al., 2017; Shalabi et al., 2019; de Souza and Bazant, 2020), we associate the use of u 1 ( r ) with an approximation of the truncated Coulomb potential u ( r ) in Equation 5, which originates in our goal to model charge discreteness rather than correlations. A complementary interpretation of u 1 ( r ) would be in terms of excluded volume interactions between adjacent ions, and such an interpretation has been proposed recently (Gupta et al., 2020a).

The function η = η ( Ψ ) that is needed to solve Equation 14 satisfies Equation 13 and is displayed together with ϕ ̄ ( Ψ ) in Figure 2 for different parameters z and ω . The solution of Equation 13 for ω = 0 is η = − tanh ( Ψ ) , irrespective of z , and is shown as a black line in the left diagram of Figure 2. Curves of the same color refer to z = 2 (red), z = 4 (blue), and z = 6 (green) and are calculated for different values of ω , as indicated in the legend. Note that for z = 2 (the red lines in Figure 2), Equation 13 yields the following simple analytic result: η Ψ = − sinh Ψ e 4 ω + sinh 2 Ψ , ϕ ̄ Ψ = e 2 ω − 1 + e 4 ω − 1 η Ψ 2 2 ⁡ sinh 2 ω , (15) which corresponds to the solution of the one-dimensional Ising model (Davis, 1996) in an external field Ψ .

From this point forward, we consider a single charged planar electrode that carries a uniform surface charge density σ and is in contact with the ionic liquid. We locate the origin of a Cartesian coordinate system at the electrode surface, with its x -axis pointing normal into the ionic liquid. The electrostatic potential Ψ = Ψ ( x ) then depends only on the distance x to the electrode, and Equation 14 reads as follows: l 2 Ψ ″ x − l 2 r 0 2 6 Ψ ″ ″ x = − η Ψ x , (16) where a prime denotes the derivative with respect to the position x . We derive solutions of this fourth-order non-linear differential equation for x ≥ 0 subject to the boundary conditions Ψ ( x → ∞ ) = Ψ ′ ( x → ∞ ) = 0 as well as Ψ ′ 0 = 0 ,   Ψ ‴ 0 = 6 s l r 0 2 , (17) where s = ν σ / ( l e ) = 4 π l B l σ / e is a scaled (dimensionless) surface charge density. Different sets of boundary conditions at the electrode surface x = 0 have been proposed in the past for fourth-order differential equations of the type in Equation 16: Ψ ′ ( 0 ) = − s / l and Ψ ‴ ( 0 ) = 0 in the original BSK model (Bazant et al., 2011) and Ψ ′ ( 0 ) = − s / l and Ψ ″ ( 0 ) = 0 based on variational free energy minimization (Gupta et al., 2020b). In addition, Ψ ′ ( 0 ) = − s / l and ( r 0 / 6 ) × Ψ ‴ ( 0 ) = Ψ ″ ( 0 ) have been proposed based on the continuity of the Maxwell stress tensor at the charged electrode (de Souza and Bazant, 2020) and on interpreting u 1 ( r ) in Equation 9 as a specific combination of Coulomb and Yukawa potentials (Bossa and May, 2020). We take the point of view that the truncated Coulomb potential applies to all charges in the system, including those within the ionic liquid and those at the electrode surface. In this case, the energy (per unit area A ) of the system associated with the long-range part of the electrostatic interactions is as follows: U A = l 2 2 ν ∫ 0 ∞ d x Ψ ′ 2 + r 0 2 6 Ψ ″ 2 = − l 2 2 ν { Ψ 0 Ψ ′ 0 − r 0 2 6 Ψ ‴ 0 + r 0 2 6 Ψ ′ 0 Ψ ″ 0 + ∫ 0 ∞ d x Ψ Ψ ″ − r 0 2 6 Ψ ″ ″ } . (18) Its first variation is as follows: δ U A = − l 2 ν { Ψ 0 δ Ψ ′ 0 − r 0 2 6 δ Ψ ‴ 0 + r 0 2 6 Ψ ′ 0 δ Ψ ″ 0 + ∫ 0 ∞ d x Ψ δ Ψ ″ − r 0 2 6 δ Ψ ″ ″ } . (19) The integral in Equation 19 is ( 1 / ν ) ∫ 0 ∞ d x Ψ δ η , which, when combined with the variation in the non-electrostatic contributions to the free energy, will vanish, as demonstrated in Equation 12. Poisson’s equation for the truncated Coulomb potential, l 2 / ν [ Ψ ″ − ( r 0 2 / 6 ) Ψ ″ ″ ] = − η / ν , relates the potential to the local volume charge density e η / ν . Integrating that equation across the surface of the electrode yields l [ Ψ ′ ( 0 ) − ( r 0 2 / 6 ) Ψ ‴ ( 0 ) ] = − s . Because the scaled surface charge density s is fixed, the ensuing relation δ Ψ ′ ( 0 ) − ( r 0 2 / 6 ) × δ Ψ ‴ ( 0 ) = 0 causes the first surface term in Equation 19 to vanish. Vanishing of the second surface term demands Ψ ′ ( 0 ) = 0 . Because of Ψ ′ ( 0 ) = 0 , the term ( r 0 2 / 6 ) × Ψ ′ ( 0 ) Ψ ″ ( 0 ) in the third line of Equation 18 vanishes. The remaining two terms, the electrostatic energy due to the surface charges on the electrode and the electrostatic energy due the ions in the ionic liquid, then render Equation 18 identical to Equation 4. Hence, the boundary conditions at the electrode, Ψ ′ ( 0 ) = 0 and Ψ ‴ ( 0 ) = 6 s / ( l r 0 2 ) , lead to the consistent vanishing of δ F with F defined in Equation 10, including all surface terms.

The following three comments complete our discussion of the boundary conditions. First, for x < 0 , the potential is constant, Ψ ( x ) = Ψ ( 0 ) . Hence, because the region x < 0 does not contribute to U , the integration in Equation 18 may start from x = 0 instead of also including negative values of x . Second, short-range electrostatic interactions of the ionic liquid with the electrode surface contribute another surface term − η ( 0 ) s l / b to U / A in Equation 18. Because the variation of this term vanishes, it does not add to the boundary conditions. Third, our approach recovers the boundary condition Ψ ′ ( 0 ) = − s / l , used previously to solve the second-order self-consistency differential equation l 2 Ψ ″ ( x ) = − η ( Ψ ( x ) ) that results from our model for r 0 = 0 (Downing et al., 2018b).

The solution Ψ ( x ) of Equation 16 reveals how the surface potential Ψ 0 = Ψ ( x = 0 ) depends on the scaled surface charge density s , allowing us to compute the differential capacitance as follows: C diff = ϵ ϵ 0 l d s d Ψ 0 = ϵ ϵ 0 l C ̄ diff . Calculation of the scaled differential capacitance C ̄ diff = d s / d Ψ 0 is a major focus of the present work.

For small potentials, Ψ ≪ 1 , only the linear part of the relationship η ( Ψ ) is significant. We use Equation 13 to perform a series expansion of η ( Ψ ) up to the linear order. The result, η = − Ψ / 1 + ( z / 2 ) × e 2 ω − 1 , can conveniently be expressed as η = − ( 3 / 2 ) × ( l / r c ) 2 × Ψ , where we have defined the length: r c = l 3 2 1 + z 2 e 2 ω − 1 . (20) This length will be shown to play a role in separating different regimes of the solution Ψ ( x ) for the linear problem. We thus solve the equation Ψ ″ x − r 0 2 6 Ψ ″ ″ x = 3 2 r c 2 Ψ x , subject to Ψ ′ ( 0 ) = 0 , Ψ ‴ ( 0 ) = ( 6 / r 0 2 ) × ( s / l ) , Ψ ( x → ∞ ) = Ψ ′ ( x → ∞ ) = 0 . The solution of the linear problem, denoted in the following by Ψ lin ( x ) , can generally be expressed as the sum of two contributions: Ψ lin x = s a 1 e − λ 1 x r 0 + a 2 e − λ 2 x r 0 , (21) with λ 1 3 = 1 + 1 − r 0 r c 2 , λ 2 3 = 1 − 1 − r 0 r c 2 . (22) The two constants a 1 and a 2 follow from the boundary conditions in Equation 17, a 1 = − r 0 3 l 1 1 − r 0 r c 2 1 + 1 − r 0 r c 2 , a 2 = r 0 3 l 1 1 − r 0 r c 2 1 − 1 − r 0 r c 2 . (23) In the following, we discuss how the structure of the solution changes as r 0 is varied. For r 0 = 0 , the corresponding potential, Ψ lin x = 2 3 r c l s e − 3 2 x r c , (24) is characterized by a single decay length. For 0 < r 0 < r c , there is a double-exponential decay with the two decay lengths r 0 / λ 1 and r 0 / λ 2 . At r 0 = r c , the decay length r 0 / λ 2 diverges, and the potential becomes as follows: Ψ lin x = s r c l 1 3 + x r c e − 3 x r c . (25) Finally, for r 0 > r c , the quantities λ 1 = λ 3 + i λ 4 and λ 2 = λ 3 − i λ 4 , as well as a 1 = a 3 − i a 4 and a 2 = a 3 + i a 4 , adopt complex conjugate values with λ 3 = 3 2 r 0 r c + 1 ,   λ 4 = 3 2 r 0 r c − 1 and a 3 = r c 6 l 1 r 0 r c + 1 ,   a 4 = r c 6 l 1 r 0 r c − 1 . The potential Ψ lin x = 2 s e − λ 3 x r 0 a 3 ⁡ cos λ 4 x r 0 + a 4 ⁡ sin λ 4 x r 0 (26) then exhibits exponentially decaying oscillations. Figure 3 shows Ψ lin ( x ) / s for ω = 1 (upper three curves) and ω = 0 (lower three curves) and three different ratios of r 0 / r c , namely r 0 = 0 (black), r 0 = r c (red), and r 0 = 2 r c (blue). Figure 3 suggests the tendency of ω to increase the magnitude and decay length of the potential. The truncation length r 0 only moderately modifies this tendency. Irrespective of the ratio r 0 / r c , the same linear relationship, Ψ 0 lin = s 2 3 r c l 1 1 + r 0 r c , (27) between surface potential Ψ 0 lin = Ψ lin ( x = 0 ) and scaled surface charge density s results. Hence, Equation 27 is valid within the linearized model for any choice of z , ω , and r 0 . This leads to a prediction for the differential capacitance C ̄ diff ( s = 0 ) = C ̄ diff lin in the limit of an uncharged electrode, where both s = 0 and Ψ 0 = 0 , C ̄ diff lin = d s d Ψ 0 lin = 1 + r 0 r c 1 + z 2 e 2 ω − 1 . (28) Recall that r c depends on z and ω through Equation 20. Clearly, increasing ω lowers C ̄ diff lin , whereas larger r 0 increases C ̄ diff lin .

The linearized model is based on the relationship η ∼ Ψ . In this section, we employ a perturbation approach to analytically model the nonlinear region, where the potential deviates only slightly from Ψ lin . To this end, we expand the relationship η ( Ψ ) up to the third order, − η = 3 2 l r c 2 Ψ + b 3 Ψ 3 . (29) From Equation 13, we find the third-order coefficient as follows: b 3 = z 2 e 2 ω + 2 e 2 ω − 1 2 − 2 6 1 + z 2 e 2 ω − 1 4 . (30) The left diagram in Figure 2 displays the third-order approximation, as specified in Equations 29, 30 (color-matching thick transparent lines). Hence, we aim to solve the nonlinear differential equation Ψ ″ x − r 0 2 6 Ψ ″ ″ x = 3 2 r c 2 Ψ x + b 3 l 2 Ψ 3 x , (31) subject to Ψ ′ ( 0 ) = 0 , Ψ ‴ ( 0 ) = ( 6 / r 0 2 ) × ( s / l ) , and Ψ ( x → ∞ ) = Ψ ′ ( x → ∞ ) = 0 . We express the potential Ψ ( x ) = Ψ lin ( x ) + b 3 Ψ per ( x ) as the sum of the solution for the linear problem ( Ψ lin ) and a small perturbation ( b 3 Ψ per ) . Substituting Ψ ( x ) into Equation 31 leads to the linear, inhomogeneous differential equation Ψ p e r ″ x − r 0 2 6 Ψ ″ ″ p e r x = 3 2 r c 2 Ψ per x + 1 l 2 Ψ lin 3 x (32) for the perturbation contribution of the potential Ψ per ( x ) , where Ψ lin ( x ) is given by Equation 21. The solution Ψ per ( x ) = Ψ inh ( x ) + Ψ hom ( x ) of Equation 32 can be expressed as the sum of a specific solution for the inhomogeneous equation (index “inh”): Ψ inh x = s 3 r 0 l 2 a 1 3 e − 3 λ 1 x r 0 3 λ 1 2 − 1 6 3 λ 1 4 − 3 2 r 0 r c 2 + 3 a 1 2 a 2 e − 2 λ 1 + λ 2 x r 0 2 λ 1 + λ 2 2 − 1 6 2 λ 1 + λ 2 4 − 3 2 r 0 r c 2 + 3 a 1 a 2 2 e − λ 1 + 2 λ 2 x r 0 λ 1 + 2 λ 2 2 − 1 6 λ 1 + 2 λ 2 4 − 3 2 r 0 r c 2 + a 2 3 e − 3 λ 2 x r 0 3 λ 2 2 − 1 6 3 λ 2 4 − 3 2 r 0 r c 2 , where λ 1 , λ 2 , a 1 , and a 2 are defined in Equations 22, 23, and a solution for the homogeneous equation (index “hom”): Ψ hom x = s 3 c 1 e − λ 1 x r 0 + c 2 e − λ 2 x r 0 . (33) The two constants c 1 and c 2 in Equation 33 can be determined such that besides Ψ per ( x → ∞ ) = Ψ per ′ ( x → ∞ ) = 0 , Equation 32 also satisfies the boundary conditions Ψ per ′ ( 0 ) = 0 and Ψ per ‴ ( 0 ) = 0 . This results in an expression for the surface contribution of the perturbation potential Ψ per ( 0 ) = s 3 ⁡ B with B = − 1 4 1 + z 2 e 2 ω − 1 5 2 × 1 + 25 6 r 0 r c + 41 12 r 0 r c 2 1 + r 0 r c 5 2 1 + 5 3 r 0 r c . (34) The total surface potential Ψ 0 = Ψ lin ( 0 ) + b 3 Ψ per ( 0 ) = s / C ̄ diff lin + s 3 B b 3 can be used to calculate the scaled differential capacitance, C ̄ diff = 1 d Ψ 0 d s = C ̄ diff lin − 3 B b 3 C ̄ diff lin 2 s 2 = C ̄ diff lin + Ω 2 s 2 , (35) analytically up to quadratic order in s . Using the expressions for C ̄ diff lin , b 3 , and B in Equations 28, 30, 34, respectively, to calculate the quadratic-order coefficient Ω = − 6 B b 3 ( C ̄ diff lin ) 2 yields the following: Ω = z 2 e 2 ω + 2 e 2 ω − 1 2 − 2 4 1 + z 2 e 2 ω − 1 5 2 × 1 + 25 6 r 0 r c + 41 12 r 0 r c 2 1 + r 0 r c 3 2 1 + 5 3 r 0 r c . (36) The quadratic expression for C ̄ diff ( s ) in Equation 35 with the analytic results for C ̄ diff lin and Ω given in Equations 28, 36 is the central outcome of this work, accounting (albeit on an approximate level) for both charge discreteness and electrostatic correlations.