Electrolyte solutions at curved electrodes. II. Microscopic approach

TL;DR

This paper employs density functional theory to analyze the microscopic structure of electrolyte solutions near curved electrodes, revealing how particle size and surface charge influence the electric double layer and capacitance.

Contribution

It introduces a microscopic approach using the civilized model to study electrolyte-electrode interactions, extending beyond mesoscopic models and highlighting the impact of microscopic details.

Findings

Small electrode radii lead to charge-independent results.

Large radii show significant effects of particle size and surface charge.

Electrode shape determines the necessity of detailed microscopic modeling.

Abstract

Density functional theory is used to describe electrolyte solutions in contact with electrodes of planar or spherical shape. For the electrolyte solutions we consider the so-called civilized model, in which all species present are treated on equal footing. This allows us to discuss the features of the electric double layer in terms of the differential capacitance. The model provides insight into the microscopic structure of the electric double layer, which goes beyond the mesoscopic approach studied in the accompanying paper. This enables us to judge the relevance of microscopic details, such as the radii of the particles forming the electrolyte solutions or the dipolar character of the solvent particles, and to compare the predictions of various models. Similar to the preceding paper, a general behavior is observed for small radii of the electrode in that in this limit the results become independent of the surface charge density and of the particle radii. However, for large electrode radii non-trivial behaviors are observed. Especially the particle radii and the surface charge density strongly influence the capacitance. From the comparison with the Poisson-Boltzmann approach it becomes apparent that the shape of the electrode determines whether the microscopic details of the full civilized model have to be taken into account or whether already simpler models yield acceptable predictions.

Tables

Table 1. Table 1: Two sets of values for those independent, dimensionless parameters which we keep fixed for each case studied numerically (see Sec.  II.4 ). The choice of values for the remaining dimensionless parameters [ κ ​ r 1 𝜅 subscript 𝑟 1 \kappa r_{1} , κ ​ r 2 𝜅 subscript 𝑟 2 \kappa r_{2} , κ ​ R 𝜅 𝑅 \kappa R , and σ / ( e ​ κ 2 ) 𝜎 𝑒 superscript 𝜅 2 \sigma/(e\kappa^{2}) ] will be quoted for the corresponding numerical results. The choices of the values of the dimensionless parameters given in the table are guided by adopting realistic values for the corresponding dimensional quantities (see the main text). The values in case A assume the Debye length as 1 / κ ≈ 9.600 × 10 − 10 ​ m 1 𝜅 9.600 superscript 10 10 m 1/\kappa\approx 9.600\times 10^{-10}\,\text{m} and 1 / κ ≈ 3.036 × 10 − 9 ​ m 1 𝜅 3.036 superscript 10 9 m 1/\kappa\approx 3.036\times 10^{-9}\,\text{m} in case B. Both cases correspond to the energy scale β ≈ 2.414 × 10 20 ​ J − 1 𝛽 2.414 superscript 10 20 superscript J 1 \beta\approx 2.414\times 10^{20}\,\text{J}^{-1} , i.e., T = 300 ​ K 𝑇 300 K T=300\,\text{K} . The parameter values for case B emerge from those of case A by multiplying the latter by ( 1 / κ A ) / ( 1 / κ B ) 1 subscript 𝜅 𝐴 1 subscript 𝜅 𝐵 (1/\kappa_{A})/(1/\kappa_{B}) and [ ( 1 / κ B ) / ( 1 / κ A ) ] 3 superscript delimited-[] 1 subscript 𝜅 𝐵 1 subscript 𝜅 𝐴 3 [(1/\kappa_{B})/(1/\kappa_{A})]^{3} , respectively, and by replacing I A subscript 𝐼 𝐴 I_{A} by I B = I A / 10 subscript 𝐼 𝐵 subscript 𝐼 𝐴 10 I_{B}=I_{A}/10 .

parameter	case A	case B
$\kappa r_0$	$0.1617$	$0.05114$
${\varrho }_0^{\text{b}}/{\kappa }^3$	$29.36$	$931.8$
$I/{\kappa }^3$	$0.05329$	$0.1685$
$\kappa m/e$	$0.04012$	$0.01269$
$\beta u$	$-1.500$	$-1.500$
$\kappa r_c$	$0.5821$	$0.1841$
${\varrho }_w/{\kappa }^3$	$9.843$	$311.3$

Equations

$$\begin{array}{r}\mathbf{\omega }=\left(\begin{array}{c}\sin (\vartheta )\cos (\phi )\\ \sin (\vartheta )\sin (\phi )\\ \cos (\vartheta )\end{array}\right)\end{array}$$

$$\begin{array}{rl}{\bar{\varrho }}_0(\mathbf{r})& :=\int {\text{d}}^2\omega {\varrho }_0(\mathbf{r},\mathbf{\omega })\\ & :=\underset{0}{\overset{\pi }{\int }}\text{d}\vartheta \sin (\vartheta )\underset{0}{\overset{2\pi }{\int }}\text{d}\phi {\varrho }_0(\mathbf{r},\vartheta ,\phi ).\end{array}$$

$$\begin{array}{rl}\beta \Omega [\varrho ]=& \int {\text{d}}^{3}r\int {\text{d}}^2\omega {\varrho }_0(\mathbf{r},\mathbf{\omega })\left[\beta V_0(\mathbf{r})-\beta {\mu }_0\right]\\ & +\sum _{i=1}^2\int {\text{d}}^{3}r{\varrho }_i(\mathbf{r})\left[\beta V_i(\mathbf{r})-\beta {\mu }_i\right]\\ & +\beta {\mathcal{F}}^{\text{id}}[\varrho ]+\beta {\mathcal{F}}^{\text{hs}}[\varrho ]+\beta {\mathcal{F}}^{\text{el}}[\varrho ]+\beta {\mathcal{F}}^{\text{att}}[\varrho ].\end{array}$$

$$\begin{array}{rl}\beta {\mathcal{F}}^{\text{id}}[\varrho ]=& \int {\text{d}}^{3}r\int {\text{d}}^2\omega {\varrho }_0(\mathbf{r},\mathbf{\omega })\left\{\ln \left[{\Lambda }_0^3{\varrho }_0(\mathbf{r},\mathbf{\omega })\right]-1\right\}\\ & +\sum _{i=1}^2\int {\text{d}}^{3}r{\varrho }_i(\mathbf{r})\left\{\ln \left[{\Lambda }_i^3{\varrho }_i(\mathbf{r})\right]-1\right\},\end{array}$$

$$\begin{array}{rl}\beta {\mathcal{F}}^{\text{el}}[\varrho ]=& \frac{1}{2}\sum _{i,j=1}^2\int {\text{d}}^{3}r\int {\text{d}}^3{r}^{\prime }\frac{\beta }{4\pi {ϵ}_0{ϵ}_{\text{ex}}}\frac{q_i{q}_j}{|\mathbf{d}|}{\varrho }_i(\mathbf{r}){\varrho }_j({\mathbf{r}}^{\prime })\Theta [|\mathbf{d}|-(r_i+r_j)]\\ & +\sum _{i=1}^2\int {\text{d}}^{3}r\int {\text{d}}^3{r}^{\prime }\int {\text{d}}^2{\omega }^{\prime }\frac{\beta }{4\pi {ϵ}_0{ϵ}_{\text{ex}}}{q}_i\frac{m{\mathbf{\omega }}^{\prime }\cdot \mathbf{d}}{|\mathbf{d}{|}^3}{\varrho }_i(\mathbf{r}){\varrho }_0({\mathbf{r}}^{\prime },{\mathbf{\omega }}^{\prime })\Theta [|\mathbf{d}|-(r_0+r_i)]\\ & +\frac{1}{2}\int {\text{d}}^{3}r\int {\text{d}}^2\omega \int {\text{d}}^3{r}^{\prime }\int {\text{d}}^2{\omega }^{\prime }\frac{\beta }{4\pi {ϵ}_0{ϵ}_{\text{ex}}}{m}^2\left[\frac{\mathbf{\omega }\cdot {\mathbf{\omega }}^{\prime }}{|\mathbf{d}{|}^3}-3\frac{(\mathbf{\omega }\cdot \mathbf{d})({\mathbf{\omega }}^{\prime }\cdot \mathbf{d})}{|\mathbf{d}{|}^5}\right]{\varrho }_0(\mathbf{r},\mathbf{\omega }){\varrho }_0({\mathbf{r}}^{\prime },{\mathbf{\omega }}^{\prime })\Theta (|\mathbf{d}|-2r_0)\\ & +\sum _{i=1}^2\underset{\mathbf{r}\in \mathcal{A}}{\int }{\text{d}}^{2}r\int {\text{d}}^3{r}^{\prime }\frac{\beta }{4\pi {ϵ}_0{ϵ}_{\text{ex}}}\frac{\sigma q_i}{|\mathbf{d}|}{\varrho }_i({\mathbf{r}}^{\prime })\Theta (|\mathbf{d}|-r_i)\\ & +\underset{\mathbf{r}\in \mathcal{A}}{\int }{\text{d}}^{2}r\int {\text{d}}^3{r}^{\prime }\int {\text{d}}^2{\omega }^{\prime }\frac{\beta }{4\pi {ϵ}_0{ϵ}_{\text{ex}}}\sigma \frac{m{\mathbf{\omega }}^{\prime }\cdot \mathbf{d}}{|\mathbf{d}{|}^3}{\varrho }_0({\mathbf{r}}^{\prime },{\mathbf{\omega }}^{\prime })\Theta (|\mathbf{d}|-r_0),\end{array}$$

$$\begin{array}{r}\Theta (x)=\{\begin{array}{ll}0,& x<0\\ 1,& x>0\end{array}\end{array}$$

$$\begin{array}{r}\beta {\mathcal{F}}^{\text{att}}[\varrho ]=\frac{1}{2}\sum _{i,j=0}^2\int {\text{d}}^{3}r\int {\text{d}}^3{r}^{\prime }\beta u\Theta (r_c-|\mathbf{d}|)\\ \times {\varrho }_i(\mathbf{r}){\varrho }_j({\mathbf{r}}^{\prime })\Theta [|\mathbf{d}|-(r_i+r_j)].\end{array}$$

$$\begin{array}{r}\beta V^{\text{att}}(\mathbf{r})=\underset{\mathcal{V}}{\int }{\text{d}}^3{r}^{\prime }{\varrho }_w\beta u\Theta (r_c-|\mathbf{d}|).\end{array}$$

$$\begin{array}{r}\beta V_i(\mathbf{r})=\{\begin{array}{ll}\infty ,& \exists {\mathbf{r}}^{\prime }\in \mathcal{V}:|\mathbf{d}|<r_i,\\ \beta V^{\text{att}}(\mathbf{r}),& \text{otherwise}.\end{array}\end{array}$$

$$\displaystyle\frac{\delta(\beta\Omega)}{\delta\varrho_{0}(\bm{r},\bm{\omega})}\bigg{|}_{\varrho^{\text{eq}}_{0},\varrho^{\text{eq}}_{1},\varrho^{\text{eq}}_{2}}=0,\quad\frac{\delta(\beta\Omega)}{\delta\varrho_{1,2}(\bm{r})}\bigg{|}_{\varrho^{\text{eq}}_{0},\varrho^{\text{eq}}_{1},\varrho^{\text{eq}}_{2}}=0.$$

$$\displaystyle\begin{aligned} \varrho_{0}&(\bm{r},\bm{\omega})=\frac{\varrho_{0}^{\text{b}}}{4\pi}\exp\Big{\{}-\beta V_{0}(\bm{r})+c_{0}^{\text{hs}}(\bm{r})-c_{0}^{\text{hs,b}}\\ &+c_{0}^{\text{att}}(\bm{r})-c_{0}^{\text{att,b}}+\beta m\bm{\omega}\cdot\left[\bm{E}(\bm{r})-\bm{E}^{\text{aux}}(\bm{r})\right]\Big{\}}\end{aligned}$$

$$\displaystyle\begin{aligned} \varrho_{i}&(\bm{r})=\varrho_{i}^{\text{b}}\exp\Big{\{}-\beta V_{i}(\bm{r})+c_{i}^{\text{hs}}(\bm{r})-c_{i}^{\text{hs,b}}\\ &+c_{i}^{\text{att}}(\bm{r})-c_{i}^{\text{att,b}}-\beta q_{i}\left[\Phi(\bm{r})-\Phi_{i}^{\text{aux}}(\bm{r})\right]\Big{\}}.\end{aligned}$$

$$\displaystyle c_{i}^{\text{x}}(\bm{r}):=-\frac{\delta(\beta\mathcal{F}^{\text{x}})}{\delta\varrho_{i}(\bm{r})},\quad c_{i}^{\text{x,b}}:=-\frac{\delta(\beta\mathcal{F}^{\text{x}})}{\delta\varrho_{i}(\bm{r})}\bigg{|}_{\text{bulk}},$$

$$\begin{array}{r}\mathbf{P}(\mathbf{r}):=\int {\text{d}}^2\omega m\mathbf{\omega }{\varrho }_0(\mathbf{r},\mathbf{\omega }),\end{array}$$

$$\begin{array}{rl}& \mathbf{E}(\mathbf{r}):=\sum _{i=1}^2\int {\text{d}}^3{r}^{\prime }\frac{q_i}{4\pi {ϵ}_0{ϵ}_{\text{ex}}}\frac{\mathbf{d}}{|\mathbf{d}{|}^3}\left[{\varrho }_i({\mathbf{r}}^{\prime })-{\varrho }_i^{\text{b}}\right]\\ & -\int {\text{d}}^3{r}^{\prime }\frac{1}{4\pi {ϵ}_0{ϵ}_{\text{ex}}}\left\{\frac{\mathbf{P}({\mathbf{r}}^{\prime })}{|\mathbf{d}{|}^3}-\frac{3\mathbf{d}[\mathbf{P}({\mathbf{r}}^{\prime })\cdot \mathbf{d}]}{|\mathbf{d}{|}^5}\right\}\\ & +\underset{{\mathbf{r}}^{\prime }\in \mathcal{A}}{\int }{\text{d}}^2{r}^{\prime }\frac{\sigma }{4\pi {ϵ}_0{ϵ}_{\text{ex}}}\frac{\mathbf{d}}{|\mathbf{d}{|}^3},\end{array}$$

$$\begin{array}{rl}& {\mathbf{E}}^{\text{aux}}(\mathbf{r}):=\\ & \sum _{i=1}^2\int {\text{d}}^3{r}^{\prime }\frac{q_i}{4\pi {ϵ}_0{ϵ}_{\text{ex}}}\frac{\mathbf{d}}{|\mathbf{d}{|}^3}\left[{\varrho }_i({\mathbf{r}}^{\prime })-{\varrho }_i^{\text{b}}\right]\Theta (r_0+r_i-|\mathbf{d}|)\\ & -\int {\text{d}}^3{r}^{\prime }\frac{1}{4\pi {ϵ}_0{ϵ}_{\text{ex}}}\left\{\frac{\mathbf{P}({\mathbf{r}}^{\prime })}{|\mathbf{d}{|}^3}-\frac{3\mathbf{d}[\mathbf{P}({\mathbf{r}}^{\prime })\cdot \mathbf{d}]}{|\mathbf{d}{|}^5}\right\}\\ & \times \Theta (2r_0-|\mathbf{d}|)\\ & +\underset{{\mathbf{r}}^{\prime }\in \mathcal{A}}{\int }{\text{d}}^2{r}^{\prime }\frac{\sigma }{4\pi {ϵ}_0{ϵ}_{\text{ex}}}\frac{\mathbf{d}}{|\mathbf{d}{|}^3}\Theta (r_0-|\mathbf{d}|),\end{array}$$

$$\begin{array}{rl}& \Phi (\mathbf{r}):=\sum _{j=1}^2\int {\text{d}}^3{r}^{\prime }\frac{1}{4\pi {ϵ}_0{ϵ}_{\text{ex}}}\frac{q_j}{|\mathbf{d}|}\left[{\varrho }_j({\mathbf{r}}^{\prime })-{\varrho }_j^{\text{b}}\right]\\ & +\int {\text{d}}^3{r}^{\prime }\frac{1}{4\pi {ϵ}_0{ϵ}_{\text{ex}}}\frac{\mathbf{P}({\mathbf{r}}^{\prime })\cdot \mathbf{d}}{|\mathbf{d}{|}^3}+\underset{{\mathbf{r}}^{\prime }\in \mathcal{A}}{\int }{\text{d}}^2{r}^{\prime }\frac{1}{4\pi {ϵ}_0{ϵ}_{\text{ex}}}\frac{\sigma }{|\mathbf{d}|}\end{array}$$

$$\begin{array}{rl}& {\Phi }_i^{\text{aux}}(\mathbf{r}):=\\ & \sum _{j=1}^2\int {\text{d}}^3{r}^{\prime }\frac{1}{4\pi {ϵ}_0{ϵ}_{\text{ex}}}\frac{q_j}{|\mathbf{d}|}\left[{\varrho }_j({\mathbf{r}}^{\prime })-{\varrho }_j^{\text{b}}\right]\Theta (r_i+r_j-|\mathbf{d}|)\\ & +\int {\text{d}}^3{r}^{\prime }\frac{1}{4\pi {ϵ}_0{ϵ}_{\text{ex}}}\frac{\mathbf{P}({\mathbf{r}}^{\prime })\cdot \mathbf{d}}{|\mathbf{d}{|}^3}\Theta (r_0+r_i-|\mathbf{d}|)\\ & +\underset{{\mathbf{r}}^{\prime }\in \mathcal{A}}{\int }{\text{d}}^2{r}^{\prime }\frac{1}{4\pi {ϵ}_0{ϵ}_{\text{ex}}}\frac{\sigma }{|\mathbf{d}|}\Theta (r_i-|\mathbf{d}|).\end{array}$$

$$\begin{array}{rl}\Delta \Phi (\mathbf{r})& =-\frac{1}{{ϵ}_0{ϵ}_{\text{ex}}}\sum _{j=1}^2{q}_j{\varrho }_j(\mathbf{r})+\frac{1}{{ϵ}_0{ϵ}_{\text{ex}}}\nabla \cdot \mathbf{P}(\mathbf{r}),\\ & \underset{\mathbf{r}\in \mathcal{A}}{\int }{\text{d}}^{2}r\mathbf{n}(\mathbf{r})\cdot \mathbf{E}(\mathbf{r})=\frac{\sigma }{{ϵ}_0{ϵ}_{\text{ex}}}|\mathcal{A}|,\\ & \underset{\lambda \to \infty }{\lim }\underset{\mathbf{r}\in \mathcal{B}(\lambda )}{\int }{\text{d}}^{2}r\mathbf{n}(\mathbf{r})\cdot \mathbf{E}(\mathbf{r})=0,\\ & \mathcal{B}(\lambda ):=\left\{\mathbf{r}+\lambda \mathbf{n}(\mathbf{r})\in {\mathbb{R}}^3|\mathbf{r}\in \mathcal{A}\right\}.\end{array}$$

$$f(\mathbf{r},\mathbf{\omega }):=\frac{{\varrho }_0(\mathbf{r},\mathbf{\omega })}{{\bar{\varrho }}_0(\mathbf{r})}=\sum _{l=0}^{\infty }\sum _{m=-l}^l{f}_{lm}(\mathbf{r})Y_{lm}(\vartheta ,\phi ).$$

$$\int {\text{d}}^2\omega f(\mathbf{r},\mathbf{\omega })=1,$$

$$P(\mathbf{r})=\sqrt{\frac{4\pi }{3}}m{\bar{\varrho }}_0(\mathbf{r})f_{1,0}(\mathbf{r}),$$

$$\displaystyle\begin{aligned} \bar{\varrho}_{0}&(\bm{r})=\varrho_{0}^{\text{b}}\Big{\{}\exp\Big{[}-\beta V_{0}(\bm{r})+c_{0}^{\text{hs}}(\bm{r})-c_{0}^{\text{hs,b}}\\ &+c_{0}^{\text{att}}(\bm{r})-c_{0}^{\text{att,b}}\Big{]}\Big{\}}\frac{\sinh\{\beta m[E(\bm{r})-E^{\text{aux}}(\bm{r})]\}}{\beta m[E(\bm{r})-E^{\text{aux}}(\bm{r})]},\end{aligned}$$

$$\begin{array}{r}{f}_{1,0}(\mathbf{r})=\sqrt{\frac{3}{4\pi }}\mathcal{L}\{\beta m[E(\mathbf{r})-E^{\text{aux}}(\mathbf{r})]\}\end{array}$$

$$k_i^{\Phi }:=-\underset{\mathbf{r}\to \infty }{\lim }\frac{{\Phi }_i^{\text{aux}}(\mathbf{r})}{\Phi (\mathbf{r})},k^E:=-\underset{\mathbf{r}\to \infty }{\lim }\frac{E^{\text{aux}}(\mathbf{r})}{E(\mathbf{r})},$$

$${\varrho }_i(\mathbf{r})$$

$${\bar{\varrho }}_0(\mathbf{r})$$

$$f_{1,0}(\mathbf{r})$$

$$\Delta \Phi (\mathbf{r})\simeq \frac{\frac{e^{2}I\beta }{{ϵ}_0{ϵ}_{\text{ex}}}\left(2+k_1^{\Phi }+k_2^{\Phi }\right)}{1+\frac{{\varrho }_0^{\text{b}}\beta m^2}{3{ϵ}_0{ϵ}_{\text{ex}}}\left(1+k^E\right)}\Phi (\mathbf{r})={\kappa }^2\Phi (\mathbf{r}).$$

$$\kappa :=\sqrt{\frac{2e^{2}I\beta }{{ϵ}_0ϵ}},$$

$${ϵ}_{\text{ex}}=ϵ\left[1+\frac{1}{2}\left(k_1^{\Phi }+k_2^{\Phi }\right)\right]-\frac{{\varrho }_0^{\text{b}}\beta m^2}{3{ϵ}_0}\left(1+k^E\right).$$

Full text

I Introduction

In the first part Reindl2016 of our study electric double layers have been discussed in detail on mesoscopic scales within the Poisson-Boltzmann (PB) theory, which was pioneered by Gouy Gouy1910 and Chapman Chapman1913 . The focus of our analysis in Ref. Reindl2016 has been on electrodes of spherical or cylindrical shape surrounded by an electrolyte solution. Such kind of setups serve as models for certain parts of so-called double layer capacitors, i.e., electrical energy storage devices which are promising candidates for supporting sustainable energy systems. In part I Reindl2016 the differential capacitance $C$, which is the change of the surface charge density upon varying the electrostatic potential at the wall and thus is experimentally accessible, has been analyzed in particular concerning its dependence on the surface charge density and on the radius of the curved electrode. The focus has been, within the PB model, on a thorough discussion of the dependence of the capacitance on the geometry of the system.

Especially for large and small curvatures of the electrodes the behavior has been discussed systematically. On one hand, the simplifying assumptions underlying the PB approach facilitate such a detailed discussion, and, on the other hand, they are also the reason for the approach to be only reliable for weak ionic strengths (below $0.2\,\text{M}=0.2\,\text{mol}/\ell$) and low electrode potentials (below $80\,\text{mV}$) in the case of monovalent salts in aqueous solutions Butt2003 . But already in these ranges deviations from experimental data are observable: the predicted differential capacitance is larger than the measured one Schmickler2010 .

Therefore Stern introduced a model which accounts for that shortcoming by combining the Gouy-Chapman description of the diffuse layer, i.e., the charge in the fluid is distributed continuously following the PB equation, with the model of a Helmholtz layer of counterions, i.e., the charge in the fluid is directly attached to the electrode surface within a molecular layer Stern1924 . Accordingly, the pure Gouy-Chapman description, which does not consider the granular character of the fluid, has been improved by introducing a Stern layer in between the electrode and the diffuse layer, i.e., the actual nonzero particle volumes are taken into account only close to the wall. As a consequence the total capacitance of the system can be considered as a circuit of two capacitors (i.e., the Stern layer and the diffuse layer) in series Schmickler2010 ; Butt2003 and, following the rules for electric circuits, the total capacitance is smaller than the Gouy-Chapman capacitance. Practically, one is left to fit the capacitance of the Stern layer to experimental results which, on one hand, might lead to good agreements with measurements but which, on the other hand, provides only a coarse microscopic picture of the electrical double layer. One possible interpretation of the results of the Stern theory is that close to the wall the permittivity is reduced Butt2003 ; Schmickler2010 .

More sophisticated models have been developed in order to describe the structure of an electrical double layer more precisely than within the mesoscopic PB approach according to which the ions are treated as pointlike charges dissolved in a homogeneous background. Three classes of corresponding models are common in the literature. Within the so-called “primitive model” (PM) the ions are considered as charged hard spherical particles which are dissolved in a solvent which is taken into account only via the permittivity. In the so-called “solvent primitive model” (SPM) the solvent particles exhibit also a nonzero volume and are often described as hard spheres. Sometimes these models are labelled in conjunction with the attribute “restricted” which means that all particle radii are equal. Yet more elaborate theories incorporate the electrostatics between ions and solvent particles by providing the latter with a dipole.

In 1980 Carnie and Chan used the so-called “civilized model” in order to model an electrolyte solution Carnie1980 . Therein, as opposed to a primitive model, both the ions and the solvent are treated on an equal basis. The ions are represented by hard spheres with charges whereas the solvent is represented by hard spheres with an embedded dipole. In Ref. Carnie1980 exact and approximative results within the mean spherical approximation are presented for the structure of the electrolyte solution at a charged planar surface. For low ionic concentrations the surface potential has the Stern layer form, i.e., the expression of the surface potential is the sum of the diffuse part of the electrical double layer and an additional part which can be interpreted as the contribution of a Stern layer. Thus the result is regarded as a derivation of the Stern layer behavior. Analytic expressions reveal that both the nonzero ion size and the dipolar solvent alter the capacitance, which otherwise reduces to the double-layer capacitance of the linearized Gouy-Chapman theory. However, this change is not very large Carnie1980 .

In Ref. Blum1981 a similar approach is used in order to complement the results of Ref. Carnie1980 . Hard core and solvent effects, both of which are taken into account, lead to a reduction of the differential capacitance, i.e., the trend indicated by the Stern model is confirmed. However, the alignment of the dipoles close to the planar wall persists for several layers into the fluid. This ordering is not confined to a single Stern-like layer next to the electrode.

The so-called reference hypernetted-chain theory is used in Refs. Torrie1989 ; Torrie1989_2 in order to examine the double layer at the surface of large spherical macroions within a multipolar hard sphere model of electrolyte solutions, i.e., the ions are charged hard spheres and the solvent molecules are hard spheres carrying point multipoles. The structure of the double layer is discussed in terms of, e.g., the potential of mean force between the macroion and a counterion, ion or solvent number density profiles, and mean electrostatic potentials. The dependences of the structure on the surface charge density, the ionic concentration, and the (macro-)ionic radii are examined.

In Ref. Biben1998 a possible realization of a microscopic model for electric double layers within density functional theory (DFT) is proposed. The functional is formulated for a mixture of spheres with embedded point charges or dipoles. Consequently, the excess part of the DFT functional splits into a Coulombic and into hard sphere parts where for the latter the fundamental measure theory (see Refs. [12, 13] in Ref. Biben1998 ) is adopted. Within this approach various density profiles near a charged planar wall are presented, which, due to the granular nature of the fluid, exhibit pronounced oscillations. The authors conclude by comparison that the primitive model is able to reproduce well ionic density profiles in low concentration regimes and beyond a certain distance from the wall. However, within the microscopic description and for higher concentrations the number density oscillations become increasingly pronounced such that a primitive model description turns out to be inadequate.

Within the microscopic model used in Ref. Lamperski1999 the Coulombic interactions are taken into account in a mean-field-like kind while the excluded-volume effects are accounted for by the Percus Yevick approach. The radii of the various particle species, i.e., anions, cations, and solvent molecules, are chosen differently. This qualitatively affects the dependence of the differential capacitance $C$ on the surface charge density $\sigma$ of the planar electrode: within the Gouy-Chapman theory $C$ is an even function of $\sigma$, whereas for unequal values of the particle radii this symmetry disappears.

Oleksy and Hansen have used DFT with an explicit solvent description in order to investigate the wetting and drying behavior of ionic solutions in contact with a charged solid substrate Oleksy2010 . All particles are treated as hard spheres and the corresponding excluded-volume correlations are taken into account by fundamental measure theory. The other interactions, i.e., the electrostatic interaction and the cohesive Yukawa attraction, are treated within mean field theory in order to ensure full thermodynamic self-consistency. The key finding is the remarkable agreement between this version of a civilized model and a previous model in which the dipoles are not explicitly taken into account.

Henderson and coworkers proposed a nonprimitive model which differs from the aforementioned models with respect to the shape of the solvent particles Henderson2012 . Within DFT hard spherical ions are considered to be dissolved in a solvent which is composed of neutral dimers, i.e., touching positively and negatively charged hard spheres. The model predicts a larger electrode potential than in the case of an implicit model like the restricted primitive model. However, for technical reasons, results could only be obtained for comparatively low potential values.

Recently, electrolyte aqueous solutions near a planar wall have been discussed within a so-called polar-solvation DFT Warshavsky2016 . Therein all particles are considered as hard spheres via the fundamental measure theory and the mean spherical approximation is used in order to calculate the remaining part of the direct correlation functions. The comparison between the results of this polar-solvation DFT, within which the solvent particles are assumed to carry an embedded dipole, and the results of the unpolar-solvation DFT of Refs. Medasani2014 ; Ovanesyan2014 , within which the solvent particles are taken into account as hard spheres only, yields, e.g., a discrepancy in the density profiles which increases with increasing electrode potential. Due to technical reasons the model in Ref. Warshavsky2016 could not be used to study ion concentrations above $10\,\text{mM}$.

For many years a lot of efforts have been spent in order to understand the structure of electric double layers either at planar walls or around macroions. Therefore various realizations of microscopic approaches, which in general include an explicit description of dipolar solvents, have been proposed Carnie1980 ; Blum1981 ; Torrie1989 ; Torrie1989_2 ; Biben1998 ; Lamperski1999 ; Oleksy2010 ; Henderson2012 ; Warshavsky2016 . Here we study the case of a double layer at a spherical electrode of arbitrary radius $R$. In the spirit of part I Reindl2016 of this study, the structure of the electrolyte solution is captured in terms of the differential capacitance $C$, which facilitates to judge the relevance of various system parameters. As compared to the mesoscopic PB approach used in part I Reindl2016 , in the present microscopic description the size of the spherical electrode affects the layering behavior of the particles due to steric effects. In addition, the charges or dipoles, embedded in the particles, interact with each other and with the charge on the electrode. Due to the interplay of various interactions, structural features, such as the spatially varying dipole orientation of the solvent particles, are expected to exhibit a comparatively complex dependence on the electrode radius. The model used within the present study is inspired by the work of Oleksy and Hansen Oleksy2010 . It incorporates the aforementioned features (non-vanishing and distinct particle volumes as well as spatially varying dipole orientations), which contribute to the differential capacitance in ways that are, due to the influence of the electrode size, not yet well understood. Apart from dealing with the hard spheres the interactions are taken into account within a random phase description Evans1979 , the relatively simple structure of which allows one to conveniently generalize the planar geometry to the spherical one. Moreover, the concept of introducing an additional cohesive attraction amongst the particles corresponds to our intention to discuss an electrolyte solution in the liquid state under realistic conditions. To that end the respective model parameters are chosen to mimic liquid water as the solvent at room temperature and ambient pressure. The differential capacitance is discussed as function of the remaining inherent system parameters, in particular, the wall curvature. Here, especially the dependence on the latter contains mechanisms which are not revealed by more primitive models such that quantitative or even qualitative differences between our model and the more primitive ones can be expected to occur.

In Sec. II.1 the present version of the density functional is discussed in detail. The corresponding Euler-Lagrange equations (ELE) are presented in Sec. II.2. The method to account for the boundary conditions in the process of obtaining the numerical solution is described in Sec. II.3. Section II.4 summarizes the chosen parameters and the notation used here. Technical details are discussed in Appendices A, B, and C. The results of the calculations are presented in Sec. III where the structures of the electric double layers are illustrated via spatially varying number density profiles. Subsequently the capacitance data obtained for the planar wall are shown for various choices of system parameters and models. Finally, the capacitance data of spherical electrodes are discussed as function of the electrode curvature and of the surface charge density. The influence of various choices of system parameters is discussed and various models are compared with each other.

II Model

II.1 Density functional

Our microscopic approach follows the one of Oleksy and Hansen in Ref. Oleksy2010 which is sometimes referred to as the so-called civilized model Carnie1980 according to which both the ions and the solvent are treated on equal footing, in contrast to the primitive model. We consider an electrolyte solution composed of three species: solvent particles (hard spheres with radius $r_{0}$ and an embedded dipole of strength $m$), monovalent cations (hard spheres with radius $r_{1}$ carrying a charge $q_{1}=e>0$, $e$ denoting the absolute value of the elementary charge), and monovalent anions (hard spheres with radius $r_{2}$ carrying a charge $q_{2}=-e<0$). The spatial extent of the electrolyte solution is restricted due to the presence of an electrode occupying the volume $\mathcal{V}\subseteq\mathbb{R}^{3}$. In the present study the focus is on two types of electrodes or walls: planar electrodes correspond to the half-space $\mathcal{V}=\{(x,y,z)\in\mathbb{R}^{3}\;|\;z<0\}$ with surface $\mathcal{A}=\{(x,y,z)\in\mathbb{R}^{3}\;|\;z=0\}$ whereas spherical electrodes of radius $R$ occupy the domain $\mathcal{V}=\{(x,y,z)\in\mathbb{R}^{3}\;|\;x^{2}+y^{2}+z^{2}<R^{2}\}$ with surface $\mathcal{A}=\{(x,y,z)\in\mathbb{R}^{3}\;|\;x^{2}+y^{2}+z^{2}=R^{2}\}$. The surfaces $\mathcal{A}$ of the walls may be homogeneously charged with a surface charge density $\sigma$. The, in general inhomogeneous, distribution of solvent particles is given by the number density $\varrho_{0}(\bm{r},\bm{\omega})$, i.e., the number of particles per volume with dipole orientation $\bm{\omega}$, $|\bm{\omega}|=1$, at position $\bm{r}\in\mathbb{R}^{3}\backslash\mathcal{V}$. In a fixed, but otherwise arbitrary, coordinate system the dipole orientation can be represented by

[TABLE]

with polar angle $\vartheta$ and azimuthal angle $\varphi$. The number density of all solvent particles at a point $\bm{r}$ irrespective of the orientation is given by

[TABLE]

In the following the integration over all orientations, i.e., over all possible values of the two angles $\vartheta$ and $\varphi$ [see Eq. (1)], is denoted as $\int\text{d}^{2}\omega$ [see, e.g., Eq. (2)]. The number density of ion species $i\in\{1,2\}$ at position $\bm{r}$ is $\varrho_{i}(\bm{r})$.

Density functional theory (DFT) is a particularly useful approach to determine these density profiles and with them the structure of the electrolyte solution in contact with the wall Evans1979 . To this end we consider the following approximation for the grand potential functional $\Omega[\varrho_{0},\varrho_{1},\varrho_{2}]=:\Omega[\varrho]$ of the number densities $\varrho_{0}(\bm{r},\bm{\omega})$ and $\varrho_{1,2}(\bm{r})$:

[TABLE]

In Eq. (3) one has $\beta=(k_{B}T)^{-1}$ with the Boltzmann constant $k_{B}$ and the absolute temperature $T$. $V_{i}$ and $\mu_{i}$ denote the external and the chemical potential of species $i$, respectively. The external potential $V_{0}(\bm{r})$ acting on the solvent particles is taken to be independent of their dipolar orientations. Unless indicated differently, volume integrals $\int\text{d}^{3}r$ run over the space $\mathbb{R}^{3}\backslash\mathcal{V}$. $\mathcal{F}^{\text{id}}$ is the ideal gas contribution,

[TABLE]

with the thermal wave lengths $\Lambda_{i}$.

The hard sphere interaction between the particles is taken into account by means of the functional $\mathcal{F}^{\text{hs}}$ which, in the present case, is the White Bear version of the fundamental measure theory (see Ref. Roth2002 ). For the contribution $\mathcal{F}^{\text{hs}}$ {and likewise $\mathcal{F}^{\text{att}}$ [see Eq. (7) below]} the orientations $\bm{\omega}$ of the dipoles do not matter which is why it is a functional of $\bar{\varrho}_{0}(\bm{r})$ [see Eq. (2)].

The electrostatic interactions, both amongst the particles and between the particles and the wall, are captured by the functional

[TABLE]

where $\epsilon_{0}$ is the vacuum permittivity,

[TABLE]

denotes the Heaviside step function, and $\bm{d}:=\bm{r}-\bm{r^{\prime}}$ is the spatial offset between positions $\bm{r}$ and $\bm{r^{\prime}}$. Within the present approach, the dielectric properties of the solvent are described in terms of an explicit mean-field-like contribution due to the solvent particles with embedded dipole moments $\bm{p}:=m\bm{\omega}$ and an implicit excess contribution given by the excess relative permittivity $\epsilon_{\text{ex}}$ which captures dielectric properties beyond the mean-field description [see also Sec. II.3 and in particular Eq. (31)]. Alternatively, one could use descriptions based on the mean spherical approximation Warshavsky2016 or model the solvent molecules as dimers in the style of Ref. Henderson2012 . However, the present approach has technical advantages and, moreover, it has turned out that the explicit dipolar contribution affects the results only weakly.

The influence of Coulomb correlation contributions has been examined in Ref. Oleksy2006 , where semi-primitive model electrolytes have been described both in terms of a mean-field density functional, which neglects Coulomb correlations, and in terms of a more complex density functional including such correlations. The outcome of these different approaches has been compared with each other and with Monte Carlo simulations. The authors have found that Coulomb correlation corrections alter the mean-field results significantly only for high surface charges in the presence of divalent cations. Based on that finding, the model in the present study is expected to accurately describe electrolytes consisting of monovalent ions. In addition, we have compared the results of the present model for an ionic strength of $0.1\,\text{M}$ with Figs. 5 and 6 in Ref. Patra2009 , which provide simulation results for a charged spherical macroparticle surrounded by an electrolyte solution within the molecular solvent model. All profiles show good agreement with the simulation data which supports the aforementioned argument that the present model is a reliable description for monovalent ions. Another consistency check has been carried out with respect to the bulk limit. To that end number density profiles around a spherical electrode, with the same size and charge as an ion, have been calculated. These density profiles, which correspond to pair distribution functions, have been compared with Fig. 2 in Ref. Wu1999 where simulation data for an electrolyte within a solvent primitive model are presented. Our results almost lie on top of the curves denoted by “HNC” and are in good agreement with the simulation data.

An attractive interaction between the particles is taken into account by the contribution $\mathcal{F}^{\text{att}}$. This interaction is rationalized by a van der Waals attraction which enables the formation of a liquid state under realistic conditions, i.e., room temperature and ambient pressure. Amplitude and range of this potential, chosen to be square-well-like, are given by $u$ and $r_{c}$, respectively:

[TABLE]

[Note that in Eq. (7) $\varrho_{0}(\bm{r},\bm{\omega})$ reduces to $\bar{\varrho}_{0}(\bm{r})$]. We make the simplifying assumption that this kind of attractive potential is the same for the interaction among the fluid particles and for their interaction with the wall particles. Thus for a homogeneous number density $\varrho_{w}$ of the wall particles, the attractive van der Waals interaction gives rise to the substrate potential

[TABLE]

The attractive van der Waals potential $V^{\text{att}}$ together with the hard repulsive interaction between wall and fluid particles comprise the external potentials $V_{i}$, $i\in\{0,1,2\}$:

[TABLE]

Further contributions, e.g., due to electrostatic image forces in the case of a dielectric contrast between wall and fluid, would have to be added to the second line of Eq. (9).

II.2 Euler-Lagrange equations

In accordance with the variational principle underlying density functional theory Evans1979 the equilibrium densities $\varrho^{\text{eq}}_{i}$ minimize the functional in Eq. (3) and thus fulfill the Euler-Lagrange equations (ELE)

[TABLE]

Their forms will be discussed in more detail below (see also Ref. Oleksy2010 ). For clarity the superscript ${}^{\text{eq}}$ is omitted and subsequently the focus is on equilibrium densities. In the bulk and thus in the absence of inhomogeneities, the density profiles entering the ELE (10) are uniform and isotropic: $\varrho_{0}(\bm{r},\bm{\omega})=\varrho_{0}^{\text{b}}/(4\pi)$, $\varrho_{1}(\bm{r})=\varrho_{1}^{\text{b}}\equiv I$, and $\varrho_{2}(\bm{r})=\varrho_{2}^{\text{b}}\equiv I$ with the ionic strength $I$. Note that $\varrho_{1}^{\text{b}}=\varrho_{2}^{\text{b}}$ implies local charge neutrality. This also holds at points sufficiently far away from the wall. By subtracting from the ELE (10) the respective expressions in the bulk, the chemical potentials $\mu_{i}$ and the lengths $\Lambda_{i}$ drop out of the equations:

[TABLE]

and, $i\in\{1,2\}$,

[TABLE]

Here the following quantities have been introduced: one-point direct correlation functions, $\text{x}\in\{\text{hs},\text{att}\}$,

[TABLE]

the polarization

[TABLE]

electric fields

[TABLE]

and

[TABLE]

as well as electric potentials

[TABLE]

and

[TABLE]

The original Heaviside functions inherited from the functional $\mathcal{F}^{\text{el}}$ in Eq. (5) are split according to $\Theta(x)=1-\Theta(-x)$. This is the reason for the appearance of the auxiliary quantities, denoted with the superscript ${}^{\text{aux}}$ [see Eqs. (16) and (18)]. Since these quantities act like electric fields or potentials, respectively, we use the same corresponding notions for them. The advantage of this separation is of technical nature: on one hand, the integration domains in $\bm{E}^{\text{aux}}$ [Eq. (16)] and $\Phi_{i}^{\text{aux}}$ [Eq. (18)] are bounded and therefore numerically manageable. On the other hand, the total electric field $\bm{E}=-\bm{\nabla}\Phi$ [Eq. (15)] and the total electric potential $\Phi$ [Eq. (17)] are determined by electrostatics and fulfill Poisson’s equation with Neumann boundary conditions:

[TABLE]

Hence, $\bm{E}$ and $\Phi$ are the solution of the boundary value problem posed by Eq. (19). Since the differential equation (19) has to be evaluated only locally, this route is technically more convenient than to perform the integrals over the whole space in Eqs. (15) and (17). The unit vectors $\bm{n}$ in the second and in the last line of Eq. (19) point into the radial direction away from the wall and are locally perpendicular to the respective surface. $|\mathcal{A}|$ denotes the area of the wall surface $\mathcal{A}$. Poisson’s equation (19) determines the potential $\Phi$ up to an additive constant which is chosen such that $\lim\limits_{\bm{r}\rightarrow\infty}\Phi(\bm{r})=0$, i.e., the potential $\Phi(\bm{r})$ at any position $\bm{r}$ corresponds to the voltage with respect to the bulk at large distances from the wall $(\bm{r}\rightarrow\infty)$. We note that the equations for the density profiles [see Eqs. (11) and (12)] depend only on the differences $\bm{E}-\bm{E}^{\text{aux}}$ and $\Phi-\Phi_{i}^{\text{aux}}$.

Due to the dependence of the ELE (11) on both the position $\bm{r}$ and the orientation $\bm{\omega}$, in general the problem has to be solved in a high-dimensional space which is difficult to handle. However, for certain geometries of the electrode, this dimension can be reduced to a large extent Oleksy2010 . To this end the distribution function $f$ of the dipole orientations is introduced and expanded in terms of spherical harmonics:

[TABLE]

Due to the definition of the orientation independent number density $\bar{\varrho}_{0}(\bm{r})$ in Eq. (2), $f$ is normalized, i.e.,

[TABLE]

which determines the value of the coefficient $f_{00}(\bm{r})=(4\pi)^{-1/2}$. In principle, the expansion in Eq. (20) leads to a dependence of the polarization $\bm{P}$ [Eq. (14)] on the coefficients $f_{lm}$ of order $l=1$, i.e., on $f_{1,0},f_{1,1}$, and $f_{1,-1}$. However, for planar and spherical electrodes the orientation of $\bm{P}$ is perpendicular to the electrode surface everywhere and to all corresponding parallel surfaces. The respective normal component $P$ is

[TABLE]

i.e., the polarization depends only on the coefficient $f_{1,0}$, provided that the polar axis of the spherical harmonics is chosen perpendicular to the electrode surface pointing away from the wall. Likewise the electric fields $\bm{E}$ and $\bm{E}^{\text{aux}}$ are perpendicular to the surface $\mathcal{A}$ with components $E$ and $E^{\text{aux}}$, respectively. This facilitates integration over the orientations $\bm{\omega}$ such that the ELE (11) for $\varrho_{0}(\bm{r},\bm{\omega})$ can be split into two equations (see Appendix A): one for the orientation independent density

[TABLE]

and another one for the coefficient

[TABLE]

with the Langevin function $\mathcal{L}(x)=\coth(x)-1/x$. Altogether the model is described by four equations: one for each of the number densities $\varrho_{i}(\bm{r})$ of the ion species $i\in\{1,2\}$ [see Eq. (12)], one for the solvent number density $\bar{\varrho}_{0}(\bm{r})$ independent of the dipole orientation [see Eq. (23)], and one for the orientation coefficient $f_{1,0}(\bm{r})$ of the dipoles [see Eq. (24)]. The entire dependence on the orientation $\bm{\omega}$ is covered by the latter quantity. Therefore and because in the case of planar and spherical electrodes the four profiles $\bar{\varrho}_{0}$, $\varrho_{1}$, $\varrho_{2}$, and $f_{1,0}$ vary only along the direction perpendicular to the surface, the dimensionality of the original problem has been reduced considerably.

II.3 Behavior at large distances from the wall

Equations (12), (16), (18), (19), (22), (23), and (24) form a complicated system of coupled nonlinear integro-differential equations, which can be solved only numerically. This requires discretization of the various profiles on a large but finite grid along the radial direction. This approach requires assumptions concerning the profiles outside the numerical grid at large distances from the wall. Here it is assumed that the one-point direct correlation functions $c_{i}^{\text{hs}}(\bm{r})$ and $c_{i}^{\text{att}}(\bm{r})$ decay rapidly towards their bulk values $c_{i}^{\text{hs,b}}$ and $c_{i}^{\text{att,b}}$, respectively, such that in Eqs. (11) and (12) the differences $c_{i}^{\text{hs}}(\bm{r})-c_{i}^{\text{hs,b}}$ and $c_{i}^{\text{att}}(\bm{r})-c_{i}^{\text{att,b}}$ can be neglected outside the numerical grid. (See Appendix C for a detailed discussion of the decay behavior of these one-point direct correlation functions.) Global charge neutrality requires the vanishing of the electric field $\bm{E}$ infinitely far away from the wall [see the third line of Eq. (19)]. Since a numerical grid can span only a finite distance from the wall, the solution of Poisson’s equation (19) in the asymptotic range outside the grid has to be determined, e.g., in terms of a linearized theory and matched with the numerical solution inside the grid. Therefore the electric fields $\bm{E}(\bm{r})$ and $\bm{E}^{\text{aux}}(\bm{r})$ as well as the potentials $\Phi(\bm{r})$ and $\Phi_{i}^{\text{aux}}(\bm{r})$ are not required to vanish outside the numerical grid. Instead it is assumed that at large distances from the wall the electric fields and potentials are sufficiently small to allow for a linearization of the exponential function in Eqs. (11) and (12) such that the ELEs are given by Eqs. (44) and (45) in Appendix B. It can be shown by means of Eqs. (15)–(18) that for equally-sized ions and weakly charged walls the quantities $E_{x}-E_{x}^{\text{aux}}$, $E_{y}-E_{y}^{\text{aux}}$, $E_{z}-E_{z}^{\text{aux}}$, $\Phi-\Phi_{1}^{\text{aux}}$, and $\Phi-\Phi_{2}^{\text{aux}}$ (i) exhibit the same asymptotic decay behavior at large distances from the wall $(\bm{r}\rightarrow\infty)$ and (ii) are proportional to the surface charge density (see Appendix B for details). We make use of this property by introducing constants $k^{E}$, $k_{1}^{\Phi}$, and $k_{2}^{\Phi}$,

[TABLE]

for systems with weakly charged planar walls, where $E$ and $E^{\text{aux}}$ denote the radial components of the electric fields. Since both numerators and denominators in Eq. (25) exhibit the same decay behavior, as well as the asymptotic proportionality to $\sigma$ in the limit $\sigma\rightarrow 0$, to leading order the constants do not vary spatially and do not depend on the surface charge density. It has turned out numerically that the constants $k^{E}$, $k_{1}^{\Phi}$, and $k_{2}^{\Phi}$ as determined for a weakly charged planar wall are valid for all curvatures and surface charges used in the present study. Moreover, we have found that this procedure works also for small differences between the particle radii, which are at most as large as the ones considered in the following. By using Eq. (25) the asymptotically leading contribution to the auxiliary fields in Eqs. (44) and (45) in Appendix B can be expressed in terms of the constants $k^{E}$, $k_{1}^{\Phi}$, and $k_{2}^{\Phi}$ as well as the total electric potential $\Phi$ and the radial component $E$ of the total electric field. [Note that the radial component of the electric field is the relevant one (see Appendix A).] A treatment analogous to the one in Sec. II.2 and in Appendix A leads to the ELEs for the four relevant profiles

[TABLE]

which correspond to simplified versions of Eqs. (12), (23), and (24). Equations (26)–(28) together with Poisson’s equation (19) lead to a linearized, modified Poisson-Boltzmann equation

[TABLE]

The requirement to recover the Debye length $1/\kappa$ in Eq. (29) with

[TABLE]

and with the total relative permittivity $\epsilon$ defines the excess relative permittivity

[TABLE]

Note that in the case of vanishing particle volumes, i.e., $k_{1}^{\Phi}=k_{2}^{\Phi}=k^{E}=0$, and vanishing dipole moment $m=0$ the excess relative permittivity equals the total relative permittivity $\epsilon_{\text{ex}}=\epsilon$. The linearized modified Poisson-Boltzmann equation (29) can be solved analytically within our geometries.

II.4 Choice of parameters

If lengths, charges, and energies are measured in units of the Debye length $1/\kappa$ [Eq. (30)], the elementary charge $e$, and the thermal energy $1/\beta=k_{B}T$, respectively, the present model of a monovalent salt solution is specified by the following eleven independent, dimensionless parameters: $\kappa r_{0}$, $\kappa r_{1}$, $\kappa r_{2}$, $\varrho_{0}^{\text{b}}/\kappa^{3}$, $I/\kappa^{3}$, $\kappa m/e$, $\beta u$, $\kappa r_{c}$, $\kappa R$, $\sigma/(e\kappa^{2})$, and $\varrho_{w}/\kappa^{3}$. This implies that those systems are equivalent, which exhibit the same values for these dimensionless parameters. We note that for this choice of forming dimensionless ratios the relative permittivity $\epsilon$ is not an independent parameter but is absorbed in the expression of the Debye length [Eq. (30)]. The present study is focused on examining the influence of the electrode geometry. Therefore and for illustration purposes, in the following some of these parameters are fixed to certain, realistic values, i.e., they are chosen such that, at best, they describe a realistic system.

Table 1 provides an overview of the corresponding dimensionless parameters for two cases A and B. The choices of the parameter values within each case and for the two cases relative to each other are guided by adopting realistic values for the corresponding dimensional quantities. These would, for example, refer to an aqueous electrolyte solution at room temperature $T=300\,\text{K}$ and ambient pressure $p\approx 1013\,\text{hPa}$. The dipole moment $m=1.85\,\text{D}\approx 6.171\times 10^{-30}\,\text{Cm}$ of the model solvent particles is chosen corresponding to the literature value of $m$ for the water molecule Moore1976 ; CRC2017 . The relative permittivity of water in static fields takes the value of $\epsilon=77.7003$ for our chosen temperature CRC2017 . The equation of state for the pure solvent is derived from the functional Eq. (3) in the bulk and is matched to the saturation properties of liquid water at $T=300\,\text{K}$, i.e., its saturation number density and its pressure NIST2016 . This fixes the particle radius of the solvent to $r_{0}\approx 1.552\times 10^{-10}\,\text{m}$. In addition the amplitude $u$ and range $r_{c}$ of the attractive interaction [Eq. (7)] are adjusted in order to obtain the best possible accordance between the first peak of the structure factor of the present model and the corresponding data for water determined by X-ray scattering (see Ref. Franks1972 and Ref. [783] therein). This way one obtains the values $u\approx-1.500\times k_{B}T$ and $r_{c}=3.6\times r_{0}$; we recall that within our approach the interaction potentials between the substrate particles and the fluid particles are chosen to be the same square well ones as the ones among the fluid particles. Finally, in our model the homogeneous number density $\varrho_{w}$ of the particles forming the electrode enters into the strength of the attractive interaction [Eq. (8)] between the wall and the fluid particles. Its value is estimated from the number density profile of the pure solvent in contact with an uncharged planar wall. The choice $\varrho_{w}\approx 1.112\times 10^{28}\,\text{m}^{-3}$ ensures that the number density peak closest to the wall matches that of water in contact with a single graphene layer Zhou2012 . Although the value of $\varrho_{w}$ is expected to depend on the chosen electrode material, the aforementioned value leads to a surprisingly good agreement also with the data corresponding to an aqueous electrolyte solution at a charged Ag-surface (see Fig. 4a in Ref. Toney1994 ). Case A corresponds to an ionic strength $I=0.1\,\text{M}$, i.e., to a Debye length $1/\kappa\approx 9.600\times 10^{-10}\,\text{m}$, whereas case B corresponds to $I=0.01\,\text{M}$, i.e., $1/\kappa\approx 3.036\times 10^{-9}\,\text{m}$. The pressure $p$ follows from the equation of state derived from the functional Eq. (3) in the bulk with the equilibrium number densities $\varrho_{0}^{\text{b}}$ and $I$. The ELEs in the bulk relate the chemical potentials $\mu_{i}$ [Eq. (3)] and the thermal wave lengths $\Lambda_{i}$ [Eq. (4)], $i\in\{0,1,2\}$, with bulk quantities which have already been quoted at the beginning of the current Subsec. II.4. Therefore, in this sense $\mu_{i}$ and $\Lambda_{i}$ are not independent variables. The solvent number density $\varrho_{0}^{\text{b}}$ has to be adjusted in order to render the required value of the pressure $p$ for all examined ionic configurations. However, these variations are marginal such that for given $\kappa$ the numerical value of $\varrho_{0}^{\text{b}}/\kappa^{3}$ in Tab. 1 is valid with the precision of four significant digits. We do not claim that the present model is able to accurately describe liquid water because it lacks crucial properties such as hydrogen bonds and the tetrahedral shape of the water molecules. Nevertheless, this procedure precludes one from choosing arbitrary parameter values which correspond to “exotic” or even unrealistic systems. In the following this type of system, corresponding to the civilized model introduced in Secs. II.1–II.3, is abbreviated by “CIV”, possibly in conjunction with additional parameter specifications or modifications (e.g., a vanishing dipole moment, $m=0$). In the following the remaining dimensionless parameters $r_{1}/r_{0}$, $r_{2}/r_{0}$, $\sigma/(e\kappa^{2})$, and $\kappa R$ are varied and their influence on the structure of the electrolyte solution is studied. The radii $r_{1}$ and $r_{2}$ of the ions are given in units of the radius $r_{0}$ of the solvent particles, which is equivalent to providing them in units of $1/\kappa$, and we choose either $r_{1}\leq r_{0}\leq r_{2}$ or $r_{2}\leq r_{0}\leq r_{1}$. By choosing special values for some of the parameters, other well known models can be obtained within the described framework. For the restricted primitive model (RPM) one has $r_{1}=r_{2}$, $m=0$, $\beta u=0$, and $\varrho_{0}^{\text{b}}=0$, and for the Poisson-Boltzmann (PB) description one has $r_{1}=r_{2}=0$, $m=0$, $\beta u=0$, and $\varrho_{0}^{\text{b}}=0$.

III Discussion

In Fig. 1 various profiles relevant for the electrostatics are displayed as functions of the distance $z$ from a charged planar wall. The CIV model, within which all particles have the same radius, and PB (see Sec. II.4) are compared with each other.

The electrode is positively charged and consequently the electrostatic potential $\Phi$ [Fig. 1(a)] has a positive value at the wall. Qualitatively there are no significant differences in $\Phi$ between PB and CIV and even quantitatively both models lead to similar results. This is in contrast to the charge density [Fig. 1(b)]. Within the microscopic CIV the centers of the fluid particles cannot get closer to the wall than their own radius. Hence, there is a discontinuity at the distance of contact. Furthermore, again due to the non-vanishing particle volumes, the charge density exhibits a layered structure close to the wall.

For clarity of the presentation the solvent profile $\bar{\varrho}_{0}(z)$ is not shown here. However, it is noteworthy that the high density of the solvent particles contributes considerably to the pronounced layering of the charge density. It is remarkable that the oscillating behavior of the charge density corresponds to a rather smooth potential $\Phi$. In contrast, in the case of PB, within which the particle volumes are neglected, the charge density exhibits a monotonic behavior. The polarization, the component $P$ of which in the direction normal to the wall [Eq. (22)] is shown in Fig. 1(c), is identically zero in PB. Within CIV also this profile has a layered structure close to the wall. Its positive value is in accordance with expectation because it corresponds to dipoles, which, on average, point away from the positively charged wall. The discontinuity at contact with the wall causes the slight kink in $\Phi$ at the same distance [see Eq. (19)]. At large distances from the wall both models exhibit a monotonic exponential decay on the scale of the Debye length $1/\kappa$.

Figure 2 shows results for the same models as in Fig. 1 but for a larger value of the reduced surface charge density. As expected, the larger surface charge density leads to an increase of the absolute values of the shown profiles. In addition nonlinear effects are more pronounced than in Fig. 1. This is clearly visible in the insets of panel (b): in Fig. 1 the PB result for the reduced charge density is almost a straight line which is in accordance with the linearized PB equation; in Fig. 2 deviations from an exponential behavior occur. The layering of the reduced charge density [Fig. 2(b)] and of the component of the reduced polarization in the direction normal to the wall [Fig. 2(c)] within the CIV model is less pronounced in the case of high surface charges. That is, in comparison with Fig. 1, the peak closest to the wall is large relative to the subsequent peaks. For large distances from the wall, PB predicts a larger value for both the potential [Fig. 2(a)] and the absolute value of the charge density [inset of Fig. 2(b)] than the CIV model does. However, at contact with the electrode the order is reversed.

Figure 3 shows results for the CIV model for various particle radii. The profiles are shown as functions of the distance $z$ from an uncharged planar wall. If all radii are equal, i.e., $r_{1}/r_{0}=r_{2}/r_{0}=1$, there is no electric field present and the profiles of the charge density and the polarization in Fig. 3 vanish identically due to symmetry reasons. This changes in the case of different values of the radii. Figure 3(c) provides an enlarged view of the charge density close to the wall. (Due to the large zoom factor there the points of the numerical grid become visible.) The arrows pointing downwards indicate the reduced positions of closest approach of the ions ($\kappa r_{1}$ and $\kappa r_{2}$). The colors of the arrows and their labels correspond to the colors of the keys. The arrow pointing upwards indicates the reduced position of closest approach of the solvent particles ($\kappa r_{0}$) which is the same for all systems shown there. The space between the electrode surface $z=0$ and the point of closest approach ($\kappa r_{2}$) of the smaller ions (here negative) cannot be penetrated by any particle; in this region the charge density is identically zero. Subsequently, in the direction away from the wall the charge density is negative because only negative ions can approach that space. This holds up to the point beyond which the positive ions are able to penetrate that space ($\kappa r_{1}$); there the charge density becomes positive. The further behavior is visible in Fig. 3(b) which shows a high but narrow positive peak. The inset of Fig. 3(b) reveals that this positive charge is compensated by the subsequent wide region of negative charge density such that global charge neutrality is fulfilled. The polarization [Fig. 3(a)] has a discontinuity at the position $z=r_{0}$, i.e., at the point of closest approach of the solvent particles. In Fig. 3(c), the arrow pointing upwards indicates this position $\kappa r_{0}$. Because the value of the radius of the solvent particles is chosen to be in between the values of the radii of the ions, the discontinuity of the polarization $(\kappa r_{0})$ is located to the right of the point of closest approach of the negative ions $(\kappa r_{2})$ and to the left of the point of closest approach of the positive ions $(\kappa r_{1})$.

In the following we discuss the properties of an electrolyte solution in contact with an electrode in terms of the differential capacitance Schmickler2010

[TABLE]

which is the change of the surface charge density $\sigma$ upon varying the (constant) potential at the wall $\Phi(\bm{r})|_{\mathcal{A}}$ taken relative to its bulk value. Within the present study the ELEs in Sec. II.2 are solved for various values of the surface charge density $\sigma$. Together with the solutions of the ELEs the electric potential $\Phi$ [Eq. (17)] is known. The relation between the potential at the wall and $\sigma$ is used in order to determine $C$ [Eq. (32)] numerically. The differential capacitance is experimentally accessible, e.g., via cyclic voltammetry, chronoamperometry, and impedance spectroscopy Butt2003 . It contains integrated properties of the structure of the electrolyte solution (see, e.g., Figs. 1–3). This facilitates the comparison with other models and the analysis of the influence of the parameters $r_{1}/r_{0}$, $r_{2}/r_{0}$, $\sigma/(e\kappa^{2})$, $I/\kappa^{3}$, and $\kappa R$.

Figure 4 depicts results of various models in which the ionic radii are chosen to be equal ($r_{1}=r_{2}$, see Sec. II.4). The electrolyte solution is in contact with a planar electrode the rescaled surface charge density $\sigma$ of which is the horizontal axis. Here and in the following, the differential capacitance is plotted in units of the double-layer capacitance $\epsilon_{0}\epsilon\kappa$ which facilitates comparison with Gouy-Chapman results (see also Ref. Reindl2016 ). Within the range shown, all displayed curves exhibit the same characteristics as the PB result: for small $\sigma$ the differential capacitance attains a constant with zero slope. Upon increasing $\sigma$ also the capacitance increases; the main differences between the models are borne out within this range. The two indicated CIV results differ with respect to the strength of the dipole moment $m$: in one case (black solid line) the latter is chosen according to Sec. II.4 and in the other case (orange dashed line) it is set to zero. The two curves almost coincide and only for large values of $\sigma$ a small deviation is visible. Hence, for these two systems the influence of the dipole moment is relatively weak. This finding is in accordance with previous studies of the differential capacitance Carnie1980 ; Warshavsky2016 as well as of wetting phenomena Oleksy2010 for which the explicit dipole description turned out to have a relatively small effect. The agreement between the simple PB and the comparatively complex CIV models throughout the studied range is remarkable, in particular when taking into account that RPM, endowed with an intermediate degree of complexity, clearly shows deviations from the otherwise common trend. An explanation for this observation could be that within CIV and PB all particles of the electrolyte solution are described consistently on the same footing whereas within RPM they are not: within PB all particles are pointlike and within the displayed cases of CIV the particles are treated as hard spheres of equal radii. In contrast, within RPM the solvent is a structureless continuum and the ions are described as hard spheres of finite size.

In Fig. 5 the influence of various choices for the ionic radii is examined within CIV: the solid curve corresponds to the case in which all radii are equal whereas they are unequal for the other curves. In the limit of an uncharged ($\sigma\rightarrow 0$) electrode the capacitance increases with increasing difference between the ionic radii. Curves with the same shade of gray correspond to the same set of ion size ratios $\{r_{1}/r_{0},r_{2}/r_{0}\}$ where the cations are the larger ions $(r_{1}/r_{0}>r_{2}/r_{0})$ for dotted curves whereas the anions are the larger ions $(r_{1}/r_{0}<r_{2}/r_{0})$ for dashed curves. As expected, curves of the same shade concur at the vertical axis: For $\sigma\rightarrow 0$ swapping the ion radii is equivalent to flipping the sign of all charges and the latter does not change the differential capacitance $C$ [see Eq. (32)]. However, for a charged electrode ($\sigma>0$) this equivalence does not hold and for the differential capacitance two branches occur. If the positive ions are the smaller (larger) particles, the capacitance curve exhibits (non-)monotonic behavior within the investigated interval of $\sigma$. Compared to the models in Fig. 4, where the ionic radii are chosen to be equal, the decrease of $C$ for small values of $\sigma$ is a new feature in Fig. 5 for $r_{1}/r_{0}>1$ and $r_{2}/r_{0}<1$. On the other hand, for $r_{1}/r_{0}<1$ and $r_{2}/r_{0}>1$ the capacitance increases rapidly thus leading to values of $C$ larger than those for the cases shown in Fig. 4. Due to the symmetries of the present model the differential capacitance fulfills the relation $C(r_{1},r_{2},\sigma)=C(r_{2},r_{1},-\sigma)$, where the first, second, and third argument are the cation radius, the anion radius, and the surface charge density, respectively. That is, if the capacitance is known for a particular system, the same capacitance is obtained for a system in which the values of the ionic radii are swapped and the surface charge density is taken to be opposite. This explains why curves with the same shade of gray meet at $\sigma=0$ in Fig. 5. Moreover the above relation enables one to extend the curves in Fig. 5 to negative values of $\sigma$. For unequal values of the ionic radii the resulting curve has a minimum at a certain nonzero value of $\sigma$ and the curve is not symmetric with respect to this minimum. The shape of such a curve is in better qualitative agreement with experimental findings than the curve for $r_{1}=r_{2}$ which is symmetric around the minimum at $\sigma=0$ (see, e.g., Ref. Valette1981 ).

So far the focus has been on planar electrodes (see Figs. 1–5). In Fig. 6 the same models as in Fig. 4 are used in order to investigate electrolyte solutions in contact with spherical electrodes of various radii $R$. The surface charge density of $\sigma/(e\kappa^{2})\approx 1.527$ is chosen such that the planar limit, i.e., for zero curvature $1/(\kappa R)=0$, corresponds to that planar electrode in Fig. 4 with the largest surface charge density. In this planar limit PB and CIV (dotted and solid black line in Fig. 4, respectively) yield almost the same value for the differential capacitance. However, in Fig. 6 differences between the two models appear for nonzero curvatures, i.e., finite electrode radii: PB predicts larger values for the capacitance than CIV. Hence, compared with the situation at planar electrodes, where PB is a surprisingly accurate approximation for CIV with equal particle radii, at curved electrodes larger deviations occur. It is likely that these differences originate from the hard-sphere character of the particles within CIV. As within PB, the charged electrode interacts with the charges of the ions and denies them access to a certain $R$-dependent volume. However, in the case of CIV in addition the layering of the particles is influenced by varying the radius of the electrode. A mechanism of such kind is not present within PB. This might explain, why between CIV and PB there are differences in the curvature dependences and why the importance of microscopic details hinges on the geometry of the electrode. Again (as in Fig. 4) the two CIV results shown in Fig. 6 are close to each other indicating that within CIV the dipole moment has no significant effect on the capacitance. Since already at a planar wall RPM exhibits clear deviations from the other models (see Fig. 4), it does not come as a surprise, that the curve predicted by it deviates considerably also at spherical electrodes (see Fig. 6). For small wall radii the various models seem to attain a linear dependence on curvature and, compared with intermediate values of the curvature, these lines are relatively close to each other.

Figures 7–13 display the differential capacitance $C$ of spherical electrodes as a function of their curvature $1/(\kappa R)$. The corresponding data are obtained within CIV (see Sec. II.4) and each curve corresponds to a fixed value of the surface charge density $\sigma$. For $1/(\kappa R)=0$ the capacitance values reduce to the corresponding ones for a planar wall. For the largest curvatures considered, the radius of the electrode approximately equals the radii of the fluid particles. It is remarkable that all systems studied exhibit a common behavior for large curvatures: irrespective of the value of $\sigma$ all curves converge to the graph corresponding to $\sigma\rightarrow 0$. However, the dependence on $\sigma$ becomes non-trivial for small curvatures. A similar general behavior for large curvatures has been observed in part I Reindl2016 of our study, where the PB model is discussed in detail. There it is possible to show analytically, that for sufficiently small radii (i.e., large curvatures) of the electrode the linearized version of the PB equation is a reliable description. Within the linearized PB theory the electrode potential is proportional to the surface charge density and hence the differential capacitance is independent of $\sigma$. This explains within PB, why at large curvatures the dependence of the capacitance $C$ on $\sigma$ disappears. It is not possible to analytically analyze the CIV model as detailed as PB. However, the data from the CIV model in the present study reveal the same behavior as in the PB model, i.e., the dependence of $C$ on $\sigma$ weakens for large curvatures. Furthermore, Fig. 14 demonstrates that the capacitance $C$ at large curvatures is also independent of the radii of the particles. Hence the capacitance exhibits a general behavior for large curvatures which is independent of $\sigma$ and of the particle radii.

Figures 7–9 are related because in the cases studied there the radii of all particles are chosen to be equal. Compared with Fig. 7, in Fig. 8 the ionic strength is reduced, and in Fig. 9 the dipole moment is set to zero. Qualitatively these three systems exhibit similar curves and their shapes resemble the results obtained within the pure PB approach in part I Reindl2016 of this study. Quantitatively, however, the distinct models render deviations which are visible most clearly in Fig. 6, where various approaches are compared for one fixed value of the surface charge density. Again, the results are only weakly affected by the strength of the dipole moment: the graphs in Fig. 7 (with the dipole moment chosen as in Sec. II.4) and Fig. 9 (no dipole moment) differ only slightly. Also in the case of planar walls (see Fig. 4 and Refs. Carnie1980 ; Oleksy2010 ; Warshavsky2016 ) the explicit dipole description turned out to have only a small effect. It would be interesting to counter-check this finding with alternative approaches such as computer simulations. Thereby it might be possible to clarify, whether the aforementioned small differences in the capacitance originate from an insufficient description of the dipoles or whether already simpler models are capable to capture sufficiently accurately the relevant structure of an electrolyte solution. In this case it might be justified to skip the comparatively sophisticated description of the dipoles. Simulations for models, which take dipoles explicitly into account, have already been carried out.

In Ref. Boda1998 results for mixtures of hard spherical ions and dipoles in contact with charged walls are presented in terms of spatially varying profiles. Reference Pegado2016 summarizes several simulation studies concerning the effective interaction between two charged surfaces separated by a solution described by ions dissolved in a Stockmayer fluid which is a Lennard-Jones fluid with an embedded point-dipole. However, computer simulations of ion-dipole mixtures are regarded to be technically difficult Boda1998 . Possibly, this is the reason why, to our knowledge, numerical capacitance data derived from models with and without explicit dipole description had not yet been compared with each other, as it is done, e.g., in Fig. 4.

Unequal particle radii can give rise to qualitatively different behaviors which may be discussed according to the sizes of the ionic species. In the case that the positive ions are the smaller particles (Figs. 10 and 11) the capacitance in the planar limit $1/(\kappa R)\rightarrow 0$ increases with increasing difference in the particle radii. Moreover the graphs become more concave (from below) for small curvatures and in particular for large surface charge densities. In the case that the positive ions are the largest particles (Figs. 12 and 13) the capacitance in the planar limit shows a more complex behavior (see also Fig. 5): for $\sigma\rightarrow 0$ the capacitance increases with increasing difference in the particle radii. However, for intermediate and large values of $\sigma$, the capacitance decreases upon increasing the particle size difference. As a consequence, some curves approach the graph for $\sigma\rightarrow 0$ from below (see, e.g., the (green) graph for $\sigma/(e\kappa^{2})\approx 0.3824$ in Fig. 13), whereas in the most other cases the convergence is from above (see Figs. 7–11).

Already in the case of planar electrodes it has become apparent that a variation of the particle radii has a relatively strong effect on the shape of capacitance data (see Figs. 4 and 5). This finding is confirmed for the case of spherical electrodes when comparing the data corresponding to unequal particle radii (Figs. 10–13) with the data corresponding to equal particle radii (Figs. 7–9), upon varying the ionic strength $I$ or the dipole moment $m$. The difference in particle radii has a strong influence on the charge distribution because the smallest species can approach the wall closest. For a planar wall this behavior is captured in Fig. 3. For spherical electrodes the interplay of these steric effects with electrostatic interactions is influenced additionally by the radius $R$ of the electrode which increases the complexity and gives rise to the various shapes of the presented capacitance data. Figure 14 facilitates the comparison of distinct data sets. The solid curve corresponds to the case in which all radii are equal whereas the radii are unequal for the cases corresponding to the other curves. For large curvatures the curves approach each other and exhibit a common behavior independent of the chosen particle radii. In view of the limiting behavior at large curvatures, as shown in Figs. 7–13, the common behavior at large curvatures is also independent of the surface charge density. At this stage it is already known that the simple PB model is a rather good approximation for CIV with equal particle radii in the limit of small electrode radii (see Fig. 6). Furthermore, in part I Reindl2016 of this study it is shown that in the limit of large curvatures these results are in accordance with the linearized PB description. Combined with the insight obtained from Fig. 14 it seems that for $1/(\kappa R)\gg 1$ the linearized PB model might be an adequate approximation for all systems displayed in Fig. 14. This finding might be interesting for describing small electrodes or highly curved parts of electrodes.

IV Summary

We have analyzed the electric double layer (EDL) of an electrolyte solution in contact with electrodes of planar or spherical shape. Inspired by the study of Oleksy and Hansen Oleksy2010 the electrolyte solution is described in terms of density functional theory (DFT) based on the functional given in Eq. (3). This approach, which is a certain version of the so-called civilized model (CIV, see Sec. II.4), takes into account all particle species on equal footing. All particles are modelled as hard spheres with non-vanishing volumes, embedded charges (in the cases of the monovalent anions or cations) or point-dipoles (in the case of the solvent molecules), and with an attractive interaction amongst all particles which enables one to discuss an electrolyte solution in the liquid state under realistic ambient conditions. This microscopic model is a possible extension of the mesoscopic Poisson-Boltzmann (PB) approach, which was used in part I of our study Reindl2016 in order to discuss EDLs at curved electrodes. Close to the wall the microscopic description gives rise to a layering behavior of the charge density and of the polarization (see Figs. 1–3) whereas the PB approach renders monotonic profiles only. As in part I Reindl2016 the structural features of the EDL enter into the differential capacitance $C$ [Eq. (32)] which facilitates the comparison of various models with each other or to evaluate the influence of various system parameters such as particle radii, dipole moment of the solvent molecules, ionic strength, surface charge density, and electrode radius. At the planar wall and for equal radii of all particles, PB and CIV lead to similar values for the capacitance (see Fig. 4). Since compared with CIV (see Sec. II.4) PB neglects many microscopic details, this finding is not obvious. Against this background, in its turn it is remarkable, that in the case of spherical electrodes of finite radii $R$ the agreement between the predictions of the two models deteriorates (see Fig. 6), i.e., the relevance of microscopic details, captured by the CIV model, depends on the geometry of the electrode. The restricted primitive model (RPM), in which the particles are not treated on equal footing, clearly exhibits a different trend in comparison with the other models (see Figs. 4 and 6). In the case of spherical electrodes the capacitance data obtained within CIV for equal particle radii are qualitatively similar to the PB results of part I Reindl2016 of this study (see Figs. 7–9). Nevertheless there are quantitative differences (see Fig. 6). Considering the dipoles explicitly has no large effect (compare the two curves labelled with CIV in Figs. 4 and 6 or compare Fig. 7 with Fig. 9). Qualitative and relatively large quantitative differences occur if the particle radii are unequal. This is the case both for planar electrodes [see Fig. 5 where the PB result turns out to be close to the solid black curve (see Fig. 4)] and for spherical electrodes (compare Fig. 7 with Figs. 10–13). However, the differences are borne out only for small and intermediate curvatures $1/(\kappa R)$. For large curvatures the capacitance curves of all considered cases exhibit a common behavior and converge to the limiting graph valid for small surface charge densities $\sigma\rightarrow 0$, i.e., in this limit the behavior becomes independent of $\sigma$ (see Figs. 7–13). Moreover, this behavior becomes also independent of the choice of the particle radii (see Fig. 14). For $1/(\kappa R)\gg 1$ the simple linearized PB model appears to be an adequate approximation of the relatively complex CIV. We conclude that the geometry of the electrode determines the relevance of microscopic details. Apart from the limit of small electrode radii, for which a general behavior is observed, PB provides acceptable estimates in the case of equal particle sizes and large electrode radii.

Appendix A Derivation of the ELEs for the solvent in the form of Eqs. (23) and (24)

Equations (23) and (24) follow from the ELE (11) which contains all information needed about the solvent. It can be written as

[TABLE]

Due to the dependence of Eq. (33) on both the position $\bm{r}$ and the orientation $\bm{\omega}$, in general the equation has to be solved in a high-dimensional space. In order to reduce the dimensionality of the problem we focus only on the relevant information contained in Eq. (33). With Eq. (33) and the definition of the orientationally integrated number density $\bar{\varrho}_{0}(\bm{r})$ of the solvent [Eq. (2)] one has

[TABLE]

In order to carry out the angular integration in Eq. (34) the orientation vector $\bm{\omega}$ [Eq. (1)] is represented in a coordinate system the polar axis of which points into the radial direction away from the wall, i.e., the polar axis is parallel to the electric fields $\bm{E}(\bm{r})$ and $\bm{E}^{\text{aux}}(\bm{r})$. Therefore the scalar product in Eq. (34) reduces to

[TABLE]

so that

[TABLE]

which is equivalent to Eq. (23).

For this orientation of the coordinate system, the normal component $P$ [Eq. (22)] of the polarization only depends on the coefficient $f_{1,0}(\bm{r})$ of the expansion in Eq. (20):

[TABLE]

where the asterisk ∗ denotes the complex conjugate. In order to carry out the angular integration in Eq. (37) the scalar product therein is treated like in Eq. (35) such that one obtains

[TABLE]

which is equivalent to Eq. (24).

Appendix B Asymptotic behavior at large distances from the wall

In Sec. II.2 the full ELEs are presented in Eqs. (11) and (12) which provide the most general description of the model. From them the relevant reduced Eqs. (12), (23), and (24) are derived in Sec. II.2. They have to be solved numerically on a large but finite grid along the radial direction. However, beyond the finite grid radial cutoff, the position of which is characterized by the length $r_{g}$ in this Appendix, assumptions concerning the profiles have to be made. Here we focus on these assumptions and on the resulting behavior at large distances from the wall, where the external potentials [Eq. (9)] are identically zero and where it can be assumed that the quantities

[TABLE]

are sufficiently small to allow for linearization of the exponential functions in Eqs. (11) and (12). The following considerations focus on the general case of a spherical wall; the corresponding results for a planar wall can be obtained from taking the limit of infinite wall radius. With this, at large distances from the wall the ELEs (11) and (12) take the form

[TABLE]

Within this Appendix the focus is on the spatial decay of the electrostatic quantities $\Delta\Phi_{i}(\bm{r})$ and $\Delta\bm{E}(\bm{r})$ which follow from Eqs. (15)–(18). These expressions are independent of the contribution $\Delta c_{0}^{\text{hs,att}}(\bm{r})$. The latter can enter $\Delta\Phi_{i}(\bm{r})$ and $\Delta\bm{E}(\bm{r})$ only via the polarization $\bm{P}(\bm{r})$. However, due to the angular integration in Eq. (14), terms in $\varrho_{0}(\bm{r},\bm{\omega})$ which are independent of the orientation $\bm{\omega}$, such as $\Delta c_{0}^{\text{hs,att}}(\bm{r})$ in Eq. (42), do not contribute to $\bm{P}(\bm{r})$. In the limit of equal radii of the ions, $r_{1}=r_{2}$, the contributions $\Delta c_{i}^{\text{hs,att}}(\bm{r})$, $i\in\{1,2\}$, are equal and drop out of the sums in the first terms of Eqs. (15)–(18) due to the same absolute value of the ionic charges $q_{1}=-q_{2}$. In this limit, as it will turn out in the following, the electrostatic contributions $\Delta\bm{E}(\bm{r})$ and $\Delta\Phi_{i}(\bm{r})$ exhibit the same decay behavior. The corresponding length scale is set by the Debye length $1/\kappa$ [Eqs. (29) and (30)]. In Appendix C the asymptotic decay behavior of $\Delta c_{i}^{\text{hs,att}}(\bm{r})$ is discussed. There it is shown that in the limit of equal particle radii the one-point direct correlation functions decay on the length scale of the bulk correlation length $\xi$ emerging due to the presence of the hard spherical and attractive interactions. Except very close to critical points, at which $\xi$ diverges, $\xi$ is typically much smaller than $1/\kappa$ such that the contributions $\Delta c_{i}^{\text{hs,att}}(\bm{r})$ decay rapidly towards zero and can be neglected in Eqs. (42) and (43) for large distances from the wall. This can be expected to hold also for small differences between the particle radii which is the limit we focus on in the present study. Accordingly, our approximation for the ELEs at large distances from the wall reads

[TABLE]

Equations (15)–(18) are discussed in detail now. For this purpose the original integration domain $\mathbb{R}^{3}\backslash\mathcal{V}$ is extended to $\mathbb{R}^{3}$. Note that within $\mathcal{V}$ the integrands are identically zero for $r_{1}=r_{2}$. The extended domain $\mathbb{R}^{3}$ is split into a region $|\bm{r}|\leq r_{g}$ close to the wall, where the full solutions are known, and into a region $|\bm{r}|>r_{g}$ further away from the wall, where the number densities are approximated either as in Eqs. (42) and (43) in the limit of equal ionic radii, or as in Eqs. (44) and (45) under the assumption of negligible one-point direct correlation function differences. This leads to

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

In the limit $r_{1}=r_{2}$ the expressions of $S_{i}^{(0)}$ and $S_{\alpha}^{(1)}$ are odd functions of $\sigma$. That is, for small surface charge densities $\sigma\rightarrow 0$ both contributions are of the order $O(\sigma)$. We will refer to this property at the end of this Appendix. Note that $\delta_{\alpha\gamma}$ is the Kronecker symbol and the abbreviation $\bm{d}=\bm{r}-\bm{r^{\prime}}$ for the spatial offset is still valid. The Greek indices $\alpha,\gamma$ denote the vector components $x,y,z$ whereas the Latin indices $i,j$ denote the ion species $1,2$. For example, the first term of Eq. (46) can be derived as follows: from Eqs. (15) and (16) one obtains

[TABLE]

In the next step, the asymptotic expression for the density profiles at large distances $|\bm{r}|>r_{g}$ from the wall [Eq. (43) or (45)] is used. We note that $q_{1}\varrho_{1}^{\text{b}}+q_{2}\varrho_{2}^{\text{b}}=0$ due to local charge neutrality. The equation is written in terms of each component $\alpha\in\{x,y,z\}$ of the electric field:

[TABLE]

In the limit of equal radii for the ions the contribution $\Delta c_{i}^{\text{hs,att}}(\bm{r^{\prime}})$ drops out of the sum. The integrations with respect to positions $|\bm{r^{\prime}}|\leq r_{g}$ close to the wall are collected in $S_{\alpha}^{(1)}$ [Eq. (52)] and will not be considered further in this example. Finally by using Eqs. (49) and (53) one has

[TABLE]

In terms of the Fourier transform

[TABLE]

Eqs. (46) and (47) lead to a system of five linear equations for the components $\Delta\hat{E}_{\alpha}(\bm{q})$, $\alpha\in\{x,y,z\}$, of the vector $\Delta\hat{\bm{E}}(\bm{q})$ and for $\Delta\hat{\Phi}_{i}(\bm{q})$ with $i\in\{1,2\}$:

[TABLE]

and

[TABLE]

With the 5-component vectors $\bm{v}$ and $\bm{s}$ as well as the $5\times 5$ matrix $\bm{M}$,

[TABLE]

defined by

[TABLE]

the system of linear equations (58) and (59) has the form

[TABLE]

with the identity matrix $\bm{1}$. With the abbreviations

[TABLE]

one obtains

[TABLE]

For various combinations of the indices $\alpha,i$ and $\alpha,\gamma$ the Fourier transforms $\hat{K}_{\alpha i}^{(1)}(\bm{q})$ and $\hat{K}_{\alpha\gamma}^{(2)}(\bm{q})$ are identically zero. Therefore the matrix $\bm{1}+\bm{M}$ exhibits this kind of block structure. Equation (68) relates the response of the fields $\bm{v}$ to an external perturbation $\bm{s}$, i.e., it is a multidimensional analogue of Yvon’s equation in Fourier space Hansen1976 . Therefore, up to irrelevant factors, the matrix $\hat{\bm{G}}(q):=(\bm{1}+\bm{M})^{-1}$ is the Fourier transform of the bulk two-point direct correlation functions of the field components $\Delta E_{x,y,z}$ and of the potential differences $\Delta\Phi_{1,2}$. With the inverse Fourier transform

[TABLE]

one obtains $\bm{G}(r)$, which determines the asymptotic spatial dependence of $\Delta E_{x,y,z}(\bm{r})$ and $\Delta\Phi_{1,2}(\bm{r})$. The integral in Eq. (78) can be studied by using the residue theorem. As a consequence the exponential decay of $\bm{G}(r)$ is determined by the poles of $\hat{\bm{G}}(q)$ Evans1993 ; Evans1994 . The pole $q^{\prime}+\text{i}q^{\prime\prime}\in\mathbb{C}$, with $q^{\prime},q^{\prime\prime}\in\mathbb{R}$, of $\hat{\bm{G}}(q)$ with the smallest imaginary part $|q^{\prime\prime}|$ determines the asymptotic decays of $\Delta E_{x,y,z}(\bm{r})$ and $\Delta\Phi_{1,2}(\bm{r})$ on the length scale $1/|q^{\prime\prime}|$ away from the charged wall. Since, according to Cramer’s rule, the inverse matrix $(\bm{1}+\bm{M})^{-1}=\hat{\bm{G}}(q)\propto 1/\text{det}(\bm{1}+\bm{M})$ is proportional to the reciprocal of $\text{det}(\bm{1}+\bm{M})$, the poles of $\hat{\bm{G}}(q)$ are given by the roots of the determinant $\text{det}(\bm{1}+\bm{M})$:

[TABLE]

Equation (79) can be solved numerically. We find for our parameter choices (see Sec. II.4) purely imaginary roots $\text{i}q^{\prime\prime}$, i.e., the asymptotic decay of $\Delta E_{x,y,z}(\bm{r})$ and $\Delta\Phi_{1,2}(\bm{r})$ is monotonic. However, the important finding is that all electric field components and electric potentials decay on the same length scale $1/|q^{\prime\prime}|$ and are proportional to $\sigma$ in the limit of equal radii $r_{1}=r_{2}$ of the ions and for $\sigma\rightarrow 0$. This finding is relevant for Sec. II.3 and in particular for Eq. (25).

Appendix C Asymptotic decay of the one-point direct correlation functions

Within this Appendix the asymptotic decay behavior of the one-point direct correlation functions $c_{0,1,2}^{\text{hs}}(\bm{r})$ and $c_{0,1,2}^{\text{att}}(\bm{r})$ is examined at positions far away from the wall. This behavior is related to the decay behavior of the number densities $\bar{\varrho}_{0}(\bm{r})$ and $\varrho_{1,2}(\bm{r})$ which fulfill the ELEs (23) and (12). We introduce the deviations of the number densities from their respective bulk values, $\Delta\varrho_{0}(\bm{r}):=\bar{\varrho}_{0}(\bm{r})-\varrho_{0}^{\text{b}}$ and $\Delta\varrho_{1,2}(\bm{r}):=\varrho_{1,2}(\bm{r})-\varrho_{1,2}^{\text{b}}$, and use the notation in Eqs. (39)–(41) in order to rewrite Eqs. (23) and (12) as

[TABLE]

Within this Appendix, a spherical wall is discussed, as in Appendix B, and we consider the case that all particle species have the same radius: $r_{0}=r_{1}=r_{2}$. As a consequence the one-point direct correlation functions and the differences of the electrostatic potentials are the same for all species, i.e., $\Delta c_{0,1,2}^{\text{hs,att}}(\bm{r})=:\Delta c^{\text{hs,att}}(\bm{r})$ and $\Delta\Phi_{1,2}(\bm{r})=:\Delta\Phi(\bm{r})$. In the following the equation for the solvent [Eq. (80)] is discussed in detail so that the presented procedure can be applied analogously to the two equations for the ionic species in Eq. (81). The Fourier transform [Eq. (57)] of Eq. (80) is given by

[TABLE]

We introduce the length $r_{g}$ large enough such that

[TABLE]

The integration in Eq. (82) over the whole space can be split into two domains, one with $|\bm{r}|\leq r_{g}$, where the full ELE (80) is integrated, and another one with $|\bm{r}|>r_{g}$, where the ELE (80) is approximated according to Eqs. (83)–(85). We make use of the fact that the equilibrium density $\bar{\varrho}_{0}(\bm{r})$ fulfills the ELE (80) locally which is why the integral in the domain $|\bm{r}|\leq r_{g}$ vanishes. Therefore, with the approximations in Eqs. (83)–(85), the Fourier transform in Eq. (82) reads

[TABLE]

The one-point direct correlation function $\Delta c^{\text{hs,att}}(\bm{r})$, which is a functional of the number densities $\varrho_{0,1,2}$ [see Eq. (13)], can be expressed in terms of a Taylor series expansion with respect to the bulk value:

[TABLE]

Since the integral in Eq. (86) is restricted to positions far away from the wall, i.e., $|\bm{r}|>r_{g}$, where the number densities are close to their respective bulk values and hence the deviations $|\Delta\varrho_{0,1,2}|$ are small, in Eq. (87) only terms up to and including linear order in $\Delta\varrho_{0,1,2}$ are taken into account. Note that $\Delta c^{\text{hs,att}}$ [Eq. (41)] measures the difference of the one-point direct correlation functions from their respective bulk values. Therefore evaluation in the bulk leads to $\Delta c^{\text{hs,att}}=0$ which is the zeroth order of the expansion in Eq. (87). The quantity $c^{(2)\text{hs,att,b}}$ denotes the bulk two-point direct correlation function governed by the hard spherical (hs) and the attractive (att) interaction, respectively. In Eq. (86) a conveniently chosen zero is added such that the original integration over $|\bm{r}|>r_{g}$ is written in terms of a Fourier integral and an integration over the domain $|\bm{r}|\leq r_{g}$ in order to obtain

[TABLE]

The same procedure can be applied to the sum of the two equations in Eq. (81). The contribution of the electric potential drops out of the sum because the ion species carry a charge of the same absolute value $q_{1}=-q_{2}$, thus leading to

[TABLE]

The functions $F_{0}(\bm{q})$ and $F_{1,2}(\bm{q})$ are entire, i.e., they do not possess any poles in $\bm{q}\in\mathbb{C}^{3}$, because the outer-most integrations in Eqs. (88) and (89) are those of continuous integrands which are entire in $\bm{q}$ over a compact domain. The sum of Eqs. (88) and (89),

[TABLE]

is an analogue of Yvon’s equation in Fourier space Hansen1976 because it relates the number densities with an external perturbation given by the numerator on the right hand side of Eq. (90). The expression $S(\bm{q}):=[1-\hat{c}^{(2)\text{hs,att,b}}(\bm{q})(\varrho_{0}^{\text{b}}+2I)]^{-1}$ is the bulk structure factor Hansen1976 . Hence, following the line of argument in Appendix B and recognizing that the numerator in Eq. (90) does not have poles, the asymptotic decay of the sum $\Delta\varrho_{0}(\bm{r})+\Delta\varrho_{1}(\bm{r})+\Delta\varrho_{2}(\bm{r})$ is given by the pole $q^{\prime}+\text{i}q^{\prime\prime}$ of $S(\bm{q})$ with the smallest imaginary part $|q^{\prime\prime}|$, which sets the length scale of the decay. Here the length scale can be identified as the bulk correlation length $\xi=1/|q^{\prime\prime}|$ emerging from the hard spherical and attractive interactions. Far away from the wall, i.e., at positions $\bm{r}$ at which the approximations in Eqs. (83)–(85) can be applied, the ELE for the solvent Eq. (80) and the sum of the ELEs for the ions in Eq. (81) are given by

[TABLE]

From Eqs. (91) and (92) it follows that the decay of the one-point direct correlation function difference $\Delta c^{\text{hs,att}}(\bm{r})$ is given by the decay of the number densities $\Delta\varrho_{0,1,2}(\bm{r})$. That is, $\Delta c^{\text{hs,att}}(\bm{r})$ decays on the length scale of the bulk correlation length $\xi$. This result is used in Sec. II.3 and in Appendix B in order to justify the neglect of the one-point direct correlation function differences in Eqs. (44) and (45).