Electrolytes structure near electrodes with molecular size roughness
EElectrolytes structure near electrodes with molecular size roughness
Timur Aslyamov ∗ and Iskander Akhatov Center for Design, Manufacturing and Materials, Skolkovo Institute of Science and Technology,Bolshoy Boulevard 30, bld. 1, Moscow, Russia 121205
Konstantin Sinkov † Schlumberger Moscow Research, Leningradskoe shosse 16A/3, Moscow, Russia 125171 (Dated: February 23, 2021)Understanding the electrodes’ surface morphology influence on the ions’ distribution is essen-tial for designing the supercapacitors with enhanced energy density characteristics. We develop amodel for the structure of electrolytes near the rough surface of electrodes. The model describesan effective electrostatic field’s increase and associated intensification of ions’ spatial separation atthe electrode-electrolyte interface. These adsorption-induced local electric and structure propertiesresult in notably increased values and sharpened form of the differential capacitance dependence onthe applied potential. Such capacitance behavior is observed in many published simulations, and itsdescription is beyond the capabilities of the established flat-electrodes theories. The proposed ap-proach could extend the quantitatively verified models providing a new instrument of the electrodessurface-parameters optimization for specific electrolytes.
The supercapacitors are one of the most prospectivemodern energy sources due to the outstanding charg-ing/discharging times, extremely long life-cycle, andecology-friendly process of the charge storage [1]. Thesupercapacitors’ development is mainly focused on en-hancing the capacitance properties, which are notably be-low the batteries characteristics [1]. This research direc-tion started from the vital experimental discovery of thecapacitance increase inside sub-nanoporous electrodes[2]. Further experimental [3] and theoretical [4] stud-ies demonstrated that the capacitance of the nanoporouselectrodes is an oscillating function of the size of the poreswith the largest value corresponding to pore’s size com-parable to ions diameter. These results show the impor-tance of the molecular-scale influence on the electrode-electrolyte interface storing the energy. Therefore, an-other possible option to enhance the capacitance is usingthe electrodes with geometrically heterogeneous surfaces.The systematic molecular dynamic (MD) studies [5–10]demonstrated that the nano-scale patterns on the sur-face of the electrodes result in a significant increase ofthe capacitance.The differential capacitance (DC) C d depending on theapplied potential U is traditionally used as the main char-acteristic to quantify the geometrical influence on theelectrostatic properties in the experiments [11, 12] andsimulations mentioned above. In the case of flat elec-trodes, Kornyshev [13] demonstrated that DC as a func-tion of potential exhibits camel- or bell-shapes in the de-pendence on the electrolyte packing density. The surfaceroughness results in more complicated DCs behavior; forexample, MD simulations [5] showed that the electrodesurface’s geometrical heterogeneity could alter the formof DC from the camel-shape to bell-shape. Also, in work ∗ [email protected] † [email protected] [6] the authors considered significantly rough surfaceswhich induce a sharper form of DC and the existenceof a larger number of peaks. Such properties of DC cannot be described in terms of model [13] and more recentapproaches developed for the flat [14–17] and curved [18]electrodes.In this work, we develop a mean-field model describingthe electrolyte behavior near rough electrode surfaces.We represent the roughness in terms of the correlatedrandom processes accounting for the structures in bothnormal and lateral directions to describe the morphology(see Fig. 1). Such a two-parametric surface model pro-vides height-profiles mimicking the geometry of the realmaterials [19, 20]. The random surface roughness has asubstantial advantage over deterministic well-structuredgeometries used in the MD simulations due to a widerrange of the applications to randomly distributed surface-heterogeneity such as defects, and functional groups thatinfluence ionic-electrode interface [21, 22].We model inner rough surface of porous electrodes asa realization of the correlated random process, where thestandard deviation δ and the correlation length λ definethe roughness in the normal direction and the charac-teristic lateral structure, respectively. More precisely,in our model, the solid surface profile is the correlatedGauss process z s ( x ), where x is the coordinate in thelateral direction. We put the origin of the normal co-ordinate z so that the solid profile average equals zero z s ( x ) = 0. The standard deviation defines the verticaldistribution of the solid media and the correlation prop-erties correspond to the decreasing exponential function z s ( x ) z s ( x + t ) ∼ e −| t | /λ . Therefore, the small λ inducesthe palisade of the solid peaks, the large λ results inthe formation of the sparse structure allowing the fluidmolecules to penetrate into perturbed/rough region ofthe solid medium. Thus, to describe the ions molecules’behavior near rough surfaces, both normal δ and corre-lation λ parameters are crucial. a r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b FIG. 1. The sketch illustrates the characteristic distributionof the asymmetric ions near rough surface at positive appliedpotential.
The electrolyte fills the available space that results inthe inhomogeneous density distributions of the ions mix-ture ρ i = ρ i ( x, z ). Accounting for the applied potential U the electrostatic field ψ = ψ ( x, z ) inside the pore isdefined by the Poisson equation and the boundary con-ditions βe ∆ ψ = − πλ B q in D,ψ = U at ∂D, (1)where β = 1 /k B T , k B is Botlzman constant, T is thetemperature, e is electron charge, ∆ = ∂ xx + ∂ zz is the2D Laplace operator, q = (cid:80) i Z i ρ i , Z i are the valencies, λ B = βe / (4 π(cid:15)(cid:15) ) is the Bjerrum length. The domain D = { x, z s ( x ) < z < H/ } is half-space of the pore ofwidth H above rough surface ∂D = { x, z = z s ( x ) } .Subdividing D into the internal D int = { x, δ < z 2. This behavior contrasts with the molecules distri-butions near rough surfaces that allow ions to penetrateinto the rough region of solid medium. Therefore, thefunction of the ions’ spatial distribution starts from somepoints z i < d i / 2. This starting points z i may be nega-tive for extremely rough surfaces and tends to hard wallsvalue d i / z i ( δ/d i , λ/d i ) [25].The details of the dependence of ions penetration intosolid media on the surface roughness and the diameterof molecules can be found in Appendix B. In comparisonwith the flat surface, the roughness results in the regionwhich is filled by both ions and solid molecules, but atthe same time, the vertical distribution of the solid me-dia decreases the space permitted for the ions. It can beaccounted by excluding the ratio of solid media at eachlevel z from the whole covered surface, then the permit-ted surface area as a function of z has the following form: S ( z ) = S s ( z ) = S (cid:18) − 12 erfc z √ δ (cid:19) (4)where S is the area of the surface projection on the lat-eral plane. Expression (4) depends on the standard devi-ation only and defines the vertical impact of the rough-ness (see derivation in Appendix B).To describe the electrolyte near rough surfaces wetake into account the average properties of the randomrough surface: the averaged electrostatic field; the modi-fied configuration volume; the contact condition with therough solid boundaries. These effects can be written interms of the specific Helmholtz free energy f defined inthe elementary volume V (see details in Appendix C)) f ( ρ i ) = (cid:88) i (cid:2) U ext i ( z ) + ψ ( z ) Z i (cid:3) ρ i ( z ) − k B TV log Z HS (5)where U ext i ( z ) is the hard boundary potential which isinfinity for z < z i and zero for z ≥ z i ; V = S ( z )∆ z is theelementary layer volume accounting for the decreasing ofthe fluids’ permitted space near the rough surface due tothe form of the available area S ( z ); Z HS is the partitionfunction of the hard sphere mixture [26] defined in thelayer volume V as follows: Z HS = ( V − (cid:80) i v i N i ) (cid:80) i N i (cid:81) i N i ! (6)where v i = πd i / i -thcomponent, the number of the molecules is defined fromthe density distribution N i = ρ i ( z ) V . The equilibriumdensity distributions correspond to the constant chemicalpotential µ = µ i across the inhomogeneous profiles µ = ∂F/∂ρ i . As one can see from the detailed calculationsin Appendix C the inhomogeneous density distributionshave the following form: ρ i ( z ) = s ( z ) θ ( z − z i ) ρ i e − Z i βeψ − (cid:80) i γ i + (cid:80) i γ i e − Z i βeψ (7)where ρ i is the bulk density describing the componentfar enough from the surface, γ i = v i ρ i are the modelparameters showing the packing density of the fluids. It isworth noting that in the case of z i = 0, ρ i = ρ and v i = v the result (7) agree with the corresponding expressionfrom work [13] with γ = (cid:80) i v i ρ i .Unlike combinatorial calculations used in work [13] forions with the same size, we implement the partition func-tion (6) which allows us to consider the ions with thedifferent radii. To describe the asymmetric electrolyte,Kornyshev introduced the additional dependence γ ( U )that leads to the result, which is equivalent to the flatlimit of our model (7). Therefore, the proposed approachextends the verified model [13] accounting for the elec-trolyte properties near rough electrodes. Besides model[13] it is possible to consider more recent approaches,for example, the Bazant-Storey-Kornyshev (BSK) the-ory [14] accounting for the electrostatic correlations andvery recent work [15] showing the spatial oscillations ofthe ions density. Both the BSK theory [14] and lineariza-tion of model [15] result in the modified Poisson equationcontaining addition length scale. The proposed averag-ing over the random process could also be applied to themodified Poisson equations.We consider the binary mixture of the hard spheremolecules with the opposite charges Z = − Z = 1 andnon-equal diameters d (cid:54) = d . In the absence of appliedvoltage, the mixture is electrically neutral and compo-sition bulk densities are equal ρ = ρ i . The followingdimensionless variables are introduced U ∗ = eU/k B T , z ∗ = z/d m , H ∗ = H/d m , Q ∗ = Q/ρ d m , λ ∗ B = λ B d m ρ ,where d m = min( d , d ) is the molecular diameter of thesmallest component.First, we investigate the case of a flat pore wall surface δ = 0. The minimal distances between the flat surface FIG. 2. Differential capacitance of ions mixtures with d =2 d (solid lines) and d = 1 . d (dashed lines). MD simula-tions results [27] also shown for reference (disks). and the center of ion are equal to ions’ radii z i = d i / C d = ∂Q ∗ /∂U ∗ as a function of the potential U ∗ forthe case of cations larger than anions. Similarly to themodel [13], the high potential limit of the differential ca-pacitance is C d ∼ / (cid:112) γ | U | and C d ∼ / (cid:112) γ | U | forpositive and negative U , respectively. Since compositionbulk densities ρ = ρ i are equal for Z = − Z = 1,the ratio of right and left wings of C d shown in Fig. 2is defined by the ions sizes as ( γ /γ ) / = ( d /d ) / .As one can see from Fig. 2 the differential capacitanceat negative potential, where a contribution of the largercations to charge dominates, is lower than at a positiveone. Such asymmetric behaviour agrees with publisheddata of MD simulations [27] shown in Fig. 2. Also, Fig. 2demonstrates that the number of capacitance maxima de-pends on the bulk density ρ . The curves correspondingto sufficiently low γ i exhibit two maximum points — thesharp and diffuse peaks at regions of small and large ionsprevailing contribution to the total charge, respectively.The roughness induces more striking changes in thecapacitance properties. To isolate the impact from thesurface geometry, we considered slightly asymmetricalelectrolyte with the following molecular diameters ratio d = 4 / d . The rough surface allows molecules to reachthe external region ( z < δ ), where the electrostatic fieldis defined by (3). Inside the external region, the absolutevalue of electrostatic field | ψ ∗ ( z ) | locally increases (seeFig. 3) due to the sharp decrease of the permitted sur-face S ( z ). Such potential behavior crucially influencesthe ions distributions improving the spatial separationsof the co-/counter-ions. Fig. 3 shows that the counter- FIG. 3. The relative electrostatic potential ψ ∗ /U ∗ (black),cations ρ v (red) and anions ρ v (blue) packing density dis-tributions at two applied potentials U ∗ = 4 (solid lines) and U ∗ = 8 (dashed lines). Molecular diameters ratio d = 4 / d , γ = 0 . γ = 0 . λ ∗ B = 0 . 25. Surface parameters( δ ∗ = 0 . λ ∗ = 1 . z ∗ = − . z ∗ = − . 31) corre-spond to the surface S shown in Fig. 4. Inset shows thedetailed distribution inside inner region using normalized den-sity ρ ∗ i = ρ i /ρ .FIG. 4. The differential capacitance for the binary mixturewith d = 4 / d with γ = 0 . γ = 0 . λ ∗ B = 0 . to signifi-cantly rough S (the surface parameters can be found in Ap-pendix B). Inset shows the larger scale of DC for flat S andslightly rough S surfaces. ions cumulative effect from the external region becomesoverwhelming as the applied potential increases, whilethe electrolyte behavior in the internal domain ( z > δ )remains almost unperturbed. The contribution of Q ext to the total charge is notable and significantly improvesthe capacity properties. This effect of rough electrodesis in agreement with published simulations [6], wherethe capacitance enhancement has been also attributedto the local increase of the electrostatic field near roughsurface and its influence co-/counter-ions distribution in-side electrode-electrolyte interface (cumulative density).Thus, our model combining the roughness effects on boththe electrostatic and ions spatial distribution reflects themajor insights from MD simulations.To describe the total charge dependency on the ap-plied potential, we calculated the differential capacitance C d for the rough electrodes. Fig. 4 shows our results forvarious surface morphology varying from flat to signifi-cantly rough. As one can see from Fig. 4 the roughnessnotably increases the capacitance and, in particular, thevalues of C d ( U ) maxima. Indeed, the quantitative com-parison of the results from Fig. 4 shows that the flat DC isalmost constant while the surfaces with the higher rough-ness exhibit larger C d values and extremely sharp peaks.Such DC behaviour agrees with the observations fromMD simulations [5, 6] for ionic liquids inside rough elec-trodes. The calculated DCs curves from Fig. 4 conservethe number of peaks, while the MD simulations [5, 6] pre-dict the roughness induced appearance/disappearance ofthe DCs curves maxima. The results of work [6] allow usto explain this discrepancy. The authors of [6] demon-strated that the cumulative center-of-mass ion densitiesdescribe DCs curve at a finite range of the potential con-taining only two peaks and the formation of the addi-tional peaks is related to the steric effect of the ionsspatial orientations. Therefore, our model’s applicationis limited by the influences of the electrostatic poten-tial and ion’s distributions of center of mass near roughelectrodes. However, the roughness induced transitionfrom two-peaks to one-peak DC curve simulated in work[5] can be qualitatively described in terms of our modelconsidering the suppressed maximum instead of full ex-tinction. Let us compare the ideal flat geometry and thesurface with the lowest roughness, which are noted as S and S in Fig. 4. The external region of S is filled bythe smallest ions mainly, leaving the other component atthe internal zone. Such asymmetry promotes only oneenhanced peak corresponding to the situation when thesmallest-size component is counter-ions. Thus, as onecan see from inset in Fig. 4 the increase of the surfaceroughness transforms two comparable peaks S of DC toone dominating peak S .Besides the mean-filed models our approach also canbe combined with the classical Density Functional The-ory (c-DFT) which successfully describes static [24] anddynamic [28] properties of the supercapacitors. Despiteseveral versions of c-DFT describes the adsorption of theneutral molecules on the rough uncharged surfaces [29–31] the electrostatic c-DFT has been previously appliedto flat electrodes only. Detailed theoretical investigationof the electrolyte behaviour near a rough surface couldexplain the difference of the roughness impact on theionic liquids and the solutions described in [10]. Thesmall solution molecules may form the thin film on thesurface decreasing the roughness for the consequent ionsadsorption [30, 31]. Then, the solvent molecules com-patible adsorption is the possible origin of the surfaceroughness’s insignificant influence, and the correspondingc-DFT model could provide the quantitative estimations.In conclusion, we developed a theory describing theions distribution structure and accounting for the realis-tic roughness of the porous electrodes. Our model pre-dicts the significant capacitance increase induced by ions-scale roughness. Moreover, we observed that the shapeof the differential capacitance dependency on applied po-tential changes notably with a variation of roughness andbecomes sharper as the roughness increases. Such behav-ior is in agreement with the published MD simulations for the ionic liquids near well-defined heterogeneous surfaces[5, 6]. The proposed method is semi-analytical and in-volves the numerical calculations to solve Poisson equa-tion only. The simplicity of the approach and a wideacceptance of the used roughness model are strong ad-vantages of our theory in comparison with the publishedsimulations. Combining our approach with the quantita-tively approved methods (such as c-DFT) is perspectiveto design the electrodes morphology corresponding to thedesired energy/power characteristics. ACKNOWLEDGMENTS T.A. acknowledges support from the Russian ScienceFoundation (project number: 20-72-00183). K.S. isgrateful to Schlumberger management for their supportand permission to publish this work. T.A. and K.S. con-tributed equally to this work. [1] P. Simon and Y. Gogotsi, Nature Materials , 1151(2020).[2] J. Chmiola, G. Yushin, Y. Gogotsi, C. Portet, P. Simon,and P.-L. Taberna, science , 1760 (2006).[3] C. Largeot, C. Portet, J. Chmiola, P.-L. Taberna,Y. Gogotsi, and P. Simon, Journal of the AmericanChemical Society , 2730 (2008).[4] D.-e. Jiang, Z. Jin, and J. Wu, Nano letters , 5373(2011).[5] J. Vatamanu, L. Cao, O. Borodin, D. Bedrov, and G. D.Smith, The Journal of Physical Chemistry Letters , 2267(2011).[6] L. Xing, J. Vatamanu, G. D. Smith, and D. Bedrov, Thejournal of physical chemistry letters , 1124 (2012).[7] Z. Hu, J. Vatamanu, O. Borodin, and D. Bedrov, Phys-ical Chemistry Chemical Physics , 14234 (2013).[8] J. Vatamanu, L. Xing, W. Li, and D. Bedrov, PhysicalChemistry Chemical Physics , 5174 (2014).[9] D. Bedrov, J. Vatamanu, and Z. Hu, Journal of Non-Crystalline Solids , 339 (2015).[10] J. Vatamanu, O. Borodin, M. Olguin, G. Yushin, andD. Bedrov, Journal of Materials Chemistry A , 21049(2017).[11] Y.-Z. Su, Y.-C. Fu, J.-W. Yan, Z.-B. Chen, and B.-W.Mao, Angewandte Chemie International Edition , 5148(2009).[12] Y. Lauw, M. D. Horne, T. Rodopoulos, V. Lockett,B. Akgun, W. A. Hamilton, and A. R. Nelson, Lang-muir , 7374 (2012).[13] A. A. Kornyshev, “Double-layer in ionic liquids:paradigm change?” (2007).[14] M. Z. Bazant, B. D. Storey, and A. A. Kornyshev, Phys-ical review letters , 046102 (2011).[15] J. P. de Souza, Z. A. Goodwin, M. McEldrew,A. A. Kornyshev, and M. Z. Bazant, arXiv preprintarXiv:2005.04270 (2020). [16] A. Gupta, A. G. Rajan, E. A. Carter, and H. A. Stone,Physical review letters , 188004 (2020).[17] H. Chao and Z.-G. Wang, The journal of physical chem-istry letters , 1767 (2020).[18] M. Janssen, Physical Review E , 042602 (2019).[19] C. Cirac`ı, F. Vidal-Codina, D. Yoo, J. Peraire, S.-H. Oh,and D. R. Smith, ACS Photonics , 908 (2020).[20] E. Gadelmawla, M. Koura, T. Maksoud, I. Elewa, andH. Soliman, Journal of materials processing Technology , 133 (2002).[21] X. Wang, M. Salari, D.-e. Jiang, J. C. Varela, B. Ana-sori, D. J. Wesolowski, S. Dai, M. W. Grinstaff, andY. Gogotsi, Nature Reviews Materials , 787 (2020).[22] S. Evlashin, F. Fedorov, P. Dyakonov, Y. M. Maksi-mov, A. Pilevsky, K. Maslakov, Y. O. Kuzminova, Y. A.Mankelevich, E. Voronina, S. Dagesyan, et al. , The Jour-nal of Physical Chemistry Letters , 4859 (2020).[23] M. Dambrine, I. Greff, H. Harbrecht, and B. Puig, SIAMJournal on Numerical Analysis , 921 (2016).[24] A. H¨artel, Journal of Physics: Condensed Matter ,423002 (2017).[25] A. Khlyupin and T. Aslyamov, Journal of StatisticalPhysics , 1519 (2017).[26] T. Aslyamov and I. Akhatov, Physical Review E ,052118 (2019).[27] M. V. Fedorov and A. A. Kornyshev, The Journal ofPhysical Chemistry B , 11868 (2008).[28] T. Aslyamov, K. Sinkov, and I. Akhatov, arXiv preprintarXiv:2011.04575 (2020).[29] A. V. Neimark, Y. Lin, P. I. Ravikovitch, andM. Thommes, Carbon , 1617 (2009).[30] T. Aslyamov and A. Khlyupin, The Journal of chemicalphysics , 154703 (2017).[31] T. Aslyamov, V. Pletneva, and A. Khlyupin, The Jour-nal of chemical physics , 054703 (2019).[32] G. Caloz, M. Costabel, M. Dauge, and G. Vial, Asymp-totic Analysis , 121 (2006). Appendix A: Electrostatic field near rough surface Derivation of the system (2), (3) relating the ensemble-averaged electrostatic field and charge density distributionsnear rough surface from the system (1) formulated in terms of a single realization of the surface can be divided intwo principle steps. First, one can derive a series of boundary value problems for the deterministic but rough surfaceusing asymptotic expansions with respect to the small parameter characterizing the roughness. At this step we followthe approach presented in [23] and outline derivation below. Second, the boundary value problems can be formallyaveraged over ensemble and used to obtain the problem for the approximation of ensemble-averaged field.In order to employ asymptotic expansion with respect to the surface roughness we recast equations (1) to dimen-sionless form using the pore width H and the applied voltage U as a length and potential scales respectively. Thecharge density scale is then βeU/ πλ B H . The dimensionless Poisson equation is given by∆ ψ = − q in D,ψ = 1 at ∂D, (A1)the domain D = { x, z s ( x ) < z < / } and its boundary ∂D = { x, z = z s ( x ) } .Next, we decompose the domain D with deterministic but corrugated boundary ∂D into two subdomains D = D int ∪ D ext divided by a smooth artificial boundary. Taking into account that realizations of Gaussianrandom process used here to model the pore surface mostly lie in the interval − δ < z < δ , we specify the artificialboundary ∂D int = { x, z = δ } . Thus, the internal domain D int = { x, δ < z < / } , the complementing externaldomain D ext = { x, z s ( x ) < z < δ } and its outer boundary ∂D ext = { x, z = z s ( x ) } .Using the rough surface shape z = z s ( x ) one can introduce the characteristic size of the external domain ε = 2 δ , itsscaled width h ( x ) = ( δ − z s ( x )) /ε and define the coordinate ζ = ( δ − z ) /εh in D ext such that D ext = { x, < ζ < } and the boundaries ∂D int = { x, ζ = 0 } , ∂D ext = { x, ζ = 1 } . Following [23, 32] we adopt double expansions forthe electrostatic potential and charge density with respect to ε using original coordinates x, z in D int and scaledcoordinates x, ζ in D ext u = u int def = ∞ (cid:80) k =0 ε k u int ,k ( x, z ) in D int ,u ext def = ∞ (cid:80) k =0 ε k u ext ,k ( x, ζ ) in D ext , (A2)where u = { ψ, q } . The problem (A1) is reformulated as a transmission problem coupling both values ψ and normalgradients ∂ n ψ of the internal and external solutions for electrostatic field at the boundary ∂D int ∆ ψ int = − q int in D int ,ψ int = ψ ext at ∂D int ,∂ n ψ int = ∂ n ψ ext at ∂D int , ∆ ψ ext = − q ext in D ext ,ψ ext = 1 at ∂D ext . (A3)Rewriting the equations (A3) for the external problem with respect to x, ζ , substituting expansions (A2) andmatching the coefficients of the terms with the same power of ε one can obtain O (1) problem : ∆ ψ int , = − q int , in D int ,ψ int , = ψ ext , at ∂D int , ∂ ζ ψ ext , at ∂D int ,∂ ζζ ψ ext , = 0 in D ext ,ψ ext , = 1 at ∂D ext . Explicit solution for the external domain is given by ψ ext , = 1 (A4)and thus, due to the coupling conditions at ∂D int , ψ int , solves the boundary value problem with the potentialprescribed at the shifted boundary z = δ ∆ ψ int , = − q int , in D int ,ψ int , = 1 at ∂D int . (A5) O ( ε ) problem : ∆ ψ int , = − q int , in D int ,ψ int , = ψ ext , at ∂D int , − h∂ z ψ int , = ∂ ζ ψ ext , at ∂D int ,∂ ζζ ψ ext , = 0 in D ext ,ψ ext , = 0 at ∂D ext . Solution for the external domain in this case is ψ ext , = − ( ζ − h∂ z ψ int , | z = δ (A6)and the coupling conditions at ∂D int imply that ψ int , solves the boundary value problem∆ ψ int , = − q int , in D int ,ψ int , = h∂ z ψ int , at ∂D int . (A7)Although the procedure can be continued for higher orders of ε resulting in more accurate approximation for ψ [23], we restrict consideration to O ( ε ) problem and proceed with averaging over ensemble of random surfaces. Toperform averaging we note that the equations (A5) don’t involve the scaled external domain width h . Accordingly,the ensemble averages ψ int , = ψ int , and h∂ z ψ int , = h∂ z ψ int , . Using the latter equality one can apply ensembleaveraging to both sides of the equations (A5), (A7) and combine the resulting equations to get∆ ψ [1]int = − q [1]int in D int ,ψ [1]int − δ∂ z ψ [1]int = 1 at ∂D int (A8)for the partial sums ψ [1]int = ψ int , + εψ int , and q [1]int = q int , + εq int , . Here we also used that z s = 0 and accordingly h = 1 / ψ [1]ext = ψ ext , + εψ ext , satisfies ψ [1]ext = 1 − ( ζ − εh∂ z ψ [1]int | z = δ = 1 + z∂ z ψ [1]int | z = δ . (A9)The solution of (A8), (A9) approximate the mean electrostatic potential in internal and external domains with thefirst order accuracy ψ = ψ [1] + O ( ε ). Dropping superscripts from (A8), (A9), restoring the dimensional variables andnoting that for homogeneously rough surfaces ψ = ψ ( z ), we readily arrive at the problem (2), (3). Appendix B: Random surface model and hard sphere contact condition Here we briefly discuss the local geometrical properties of the correlated random process which models the realisticsurface roughness. The detailed derivation can be found in work [25] dedicated to the effective molecular potentialbetween a fluid molecule and the solid media with the rough surface. In this approach the properties of the roughsurface geometry corresponds to the Gauss correlated random process defined by the following one- and two-pointsdistribution functions: w ( z ) = 1 √ πδ exp (cid:18) − z δ (cid:19) , (B1) w ( z , x ; z , x + t ) = 12 πσ (cid:112) − K ( t ) exp (cid:18) − z + z − K ( t ) z z σ (1 − K ( t ) ) (cid:19) , (B2)where K ( s ) is the lateral correlation function, δ is the standard deviation. Therefore, the rough solid heigth profile z s can be considered as the realisations of such random process. More precisely we consider the random profiles z s with the zero average and exponentially decreasing lateral correlations: z s ( x ) = 0 , z s ( x ) z s ( x + t ) = K ( s ) = e −| t | /λ (B3)where over-line symbol is the average over the random process, λ is the correlation length. Thus, this approach allowsus to consider two-dimensional rough surface in terms of two parameters δ and λ describing the roughness in thenormal direction and the lateral geometry, respectively.We consider spherical fluid molecule with the diameter d located near rough surfaces at the point (0 , z f ) (the x -axes origin is the location of the fluid molecule). In this case we take into account the process realisation passingbelow the fluid location that approximately equivalent the condition z s (0) ≤ z f − d/ 2. One can account for theproper realisations considering the profiles which pass through the certain starting point (0 , z ) for all z ≤ z f − d/ x -coordinate instead of time. Indeed the random process starts at moment x = 0 and the crosses some level z as thevalue of x -coordinate increases Fig. 5(b). Therefore, the averaging of z -level first hitting defines the average lengthalong the x -coordinate L FPT . In work [25] the authors calculated the analytical expression for L FPT ( z , z ) for theMarkov random correlated process: L FPT ( z , z ) = λδ (cid:90) zz e ξ δ dξ (cid:90) ξ −∞ e − η δ dη = λ (cid:90) dτ e z δ (1 − τ ) − e z δ (1 − τ ) − τ + π (cid:18) erfi z √ δ − erfi z √ δ (cid:19) (B4)As one can see from expression (B4) and Fig. 5(c) the length L FPT ( z , z ) crucially depends on the starting point z . The inset in Fig. 5(c) shows the probability distribution for the profiles at point x = 0 accounting for the fluidmolecules location z f : ρ ( z ) = w ( z ) θ ( z − z f + d/ (cid:82) ∞−∞ w ( z ) θ ( z − z f + d/ dz = w ( z ) θ ( z − z f + d/ (cid:82) z f − d/ −∞ w ( z ) dz (B5)To account for all random process realisations we calculate the average for the length (B4) over the starting point z distribution (B5). The result of the averaging at point z depends on the relative position in respect to the fluidmolecule’s boundary z f − d/ L av ( z f , z ) = (cid:26) L uav ( z f , z ) , z > z f − d/ L bav ( z f , z ) , z ≤ z f − d/ z > z f − d/ 2) the averaging accounts for all starting points in the accordance with the probability FIG. 5. The first step of the random surface averaging for the fluid molecule (black dot) at point (0 , z f ): (a) the characteristicrealisations of the random process with fixed parameters δ and λ , (b) the colored curves correspond to the surface profiles z s ( x )which pass through the point (0, z ), the colored dots define the moments when the profiles exhibit the first crossing of thecertain level z . (c) The colored curves (equation X) are the averaged lengths corresponding to ”the first passage time” for therandom profiles starting from the points z (1)0 and z (2)0 . d ∗ d ∗ δ ∗ λ ∗ z ∗ z ∗ S distribution (B5): L uav ( z f , z ) = (cid:82) z f − d/ −∞ L FPT ( z , z ) w ( z ) dz (cid:82) z f − d/ −∞ w ( z ) dz (B7)The bottom region induces an additional constraint z < z < z f − d/ L bav ( z f , z ) = (cid:82) z −∞ L FPT ( z , z ) w ( z ) dz (cid:82) z f − d/ −∞ w ( z ) dz (B8)The average length L av reflects the local properties of the solid surface geometry and depends on the fluid moleculeslocation Fig. 6(a). Therefore, the function L av ( z f , z ) defines the averaged boundary of the solid media in the vicinityof the fluid molecule at point z f . In work [25] the authors used this geometrical construction for the integration ofthe fluid-solid molecular potential. Here, since we consider hard-sphere solid-fluid potential the average interactioncorresponds to the contact condition shown in Fig. 6(b). As one can see the from Fig. 6(b) the molecule with diameter d touches the upper average length (solid boundary) in some contact point. This contact condition defines the startingpoint z c of the fluid density distribution near the rough surface with known parameters δ and λ .We calculate numerically the contact points z c ( δ/d, λ/d ) as a function of the dimensionless parameters δ/d and λ/d . The contour plot of this calculations is shown in Fig. 6(c). The nonlinear behavior of z c ( δ/d, λ/d ) from Fig. 6(c)demonstrates that the penetration inside solid crucially depends on the molecular size d . We use this rough surfacesproperty to calculate the starting points z i for non-symmetric electrolyte (the hard spheres with the diameters d i ).For example, the parameters z i ( δ/d i , λ/d i ) used in calculations of differential capacity shown in Fig. 4 can be foundin Table I.To calculate the permitted surface area S ( z ) as a function of z we define the ratio of the solid media above level z using the properties of the Gauss random process (B1). Then the surface area permitted for the fluid moleculesis defined as the residue after the excluding the ratio of the area occupied by the solids media: One can obtain thisvalue from the following expression: S ( z ) = S (cid:18) − (cid:90) ∞ z w ( η ) dη (cid:19) = S (cid:18) − 12 erfc z √ δ (cid:19) (B9) FIG. 6. (a) The average length L av (red and blue curves) calculated for various locations of the fluid molecules (red and bluedots). (b) Hard sphere contact condition defining the minimal distance between the fluid and solid media z c . Appendix C: The average distributions of ions near rough surface Here we consider the ions near rough surface accounting for the average electrostatic potential ψ ( z ) and averagedproperties of the Gauss random process. The solid boundaries influence by step-like potential U ext i which prohibitpenetration of the i -th electrolyte’s component below known point z i : U ext i = (cid:26) , z ≥ z i ∞ , z < z i (C1)The thermodynamic contribution of the ions modeled as mixture of the hard spheres with diameters d i contains inthe specific Helmholtz free energy f HS . Therefore the system Helmholtz energy has the following form: f [ ρ i ( z )] = (cid:88) i (cid:2) U ext i ( z ) + ψ ( z ) Z i (cid:3) ρ i ( z ) + f HS (C2)The equilibrium density distributions correspond to the constant chemical potentials across the normal direction: µ = µ i = ∂f∂ρ i (C3)We consider spheres as an ideal gas, but the proper volume accounts the ions size V − (cid:80) i v i N i . In our system thevolume permitted for molecules is a function the distance to the surface V ( z ) = S ( z )∆ z due to the vertical roughness.Inside this layered volume the partition function of the hard sphere ions has the following form: Z HS = ( V − (cid:80) i v i N i ) (cid:80) i N i (cid:81) i N i ! (C4)The correpsponfing Helmholtz free energy is defined as F HS = − k B T log Z HS . After using of the Stieltjes formula forlarge N i the logarithm factor has the following form:log Z HS = (cid:88) i N i log (cid:32) V ( z ) − (cid:88) i v i N i (cid:33) − (cid:88) i N i log N i /e (C5)After differentiation in respect to N i expression (C5) has the following form: ∂ N i log Z HS = log (cid:32) V ( z ) − (cid:88) i v i N i (cid:33) − (cid:88) i N i ( V ( z ) − (cid:80) i v i N i ) − log N i (cid:39) − log ρ i s ( z ) − (cid:80) i v i ρ i (C6)Expression (C6) defines the derivative of the specific energy ∂ ρ i f HS = − k B T ∂ N i log Z HS . Therefore equilibriumcondition (C3) has the following form: ρ i s ( z ) − (cid:80) i v i ρ i = θ ( z − z i ) e βµ − Z i βψ ( z ) (C7)The chemical potential can be derived from the condition that ions are homogeneous far enough from the surface and ψ → e βµ = ρ − (cid:80) i γ i (C8)where a new parameter γ i = v i ρ is introduced. Substituting expression expression (C8) one can solve the system ofequations (C7) in respect to densities: ρ i ( z ) = ρ s ( z ) θ ( z − z i ) e − Z i βeψ − (cid:80) i γ i + (cid:80) i γ i e − Z i βeψβeψ