Molecular description of electrolyte solution in a carbon aerogel electrode ∗

We develop a molecular theory of aqueous electrolyte solution sorbed in a nanoporous carbon aerogel electrode, based on the replica reference interaction site model (replica RISM) for realistic molecular quenched-annealed systems. We also briey review applications of carbon aerogels for supercapacitor and electrochemical separation devices, as well as theoretical and computer modelling of disordered porous materials. The replica RISM integral equation theory yields the microscopic properties of the electrochemical double layer formed at the surface of carbon aerogel nanopores, with due account of chemical specicities of both sorbed electrolyte and carbon aerogel material. The theory allows for spatial disorder of aerogel pores in the range from micro- to macroscopic size scale. We considered ambient aqueous solution of 1 M sodium chloride sorbed in two model nanoporous carbon aerogels with carbon nanoparticles either arranged into branched chains or randomly distributed. The long-range correlations of the carbon aerogel nanostructure substantially affect the properties of the electrochemical double layer formed by the solution sorbed in nanopores.


Introduction
Nanoporous aerogels and foams have found a striking number of scientific and technical applications.Owing to exceptional physico-chemical properties, these novel nanostructured materials have been proven to be superior as catalysts, sorbers, fuel storage, supercapacitor electrodes, electrosorption cells, sensors, lightweight optics, and many other devices [1].Aerogels are produced in a variety of forms with programmed porosity, catalytic composition, thermal, optical, acoustic, mechanical and electrical features [2,3].In particular, nanoporous carbon aerogels were first developed by Pekala and co-workers [4][5][6] as electrode material of the electrochemical double layer (EDL) aerocapacitor, an advanced high-energy electrochemical storage device [7].It is expected to have specific capacitance approaching 100-200 F/g [8].Carbon aerogels with a pore diameter ranging from 3 to 13 nm have been proven to exhibit the best voltammetry characteristics and the highest capacitance values of 70-150 F/g [9].This provides energy and power densities of respectively 2 Wh/kg and 8 kW/kg for an aqueous electrolyte aerocapacitor with operating voltage of 1 V/cell.Further enhancement of aerocapacitors can be achieved by appropriate selection of the electrolyte, modification of the pores structure and size distribution to increase the accessible surface area, and by functionalization of the carbon surface resulting in pseudocapacitance effects due to quick faradaic charge transfer reactions at the surface of nanopores [9].The latter can be done by special oxidation of the surface [10], attachment of electroactive polymers also providing high charge/discharge rates due to fast access by electrolyte ions and hence high charge/discharge rates for high power density applications [11,12], and insertion of electroactive centers of transition metals oxides in the aerogel [13,14].Attractive also remains further improvement of an aerocapacitor with organic solvent which has an enhanced specific energy for the higher decomposition voltage of the electrolyte (3-5 V) but a poor specific power due to its low conductivity [9].
Another novel process based on carbon aerogel electrodes is electrochemical separation, elaborated by Farmer [15] for removing ionic impurities from aqueous waste streams.As contaminated water flow passes through the electrosorption cell, charged impurities are attracted to the electrode surface and held within the EDL to produce an effluent of purified water.A variety of ions can be scavenged from solution, such as the major species in underground aquifers, seawater, and storage tanks (e.g., Na + , NH + 4 , Cl − , ClO − 4 , NO − 3 ) [8,16,17], and many heavy metals (Cu, Zn, Ni, Cd, Cr, Pb, U) [18,19].The electrosorption on nanoporous carbon aerogel electrodes has been demonstrated to be superior in efficiency and device durability to previous cells with electrodes made of activated carbon powders or packed carbon particles.Unlike such a deionization technology as ion exchange, it requires no acids, bases, or salt solutions for regeneration of the system, which minimizes secondary wastes [8].
The length scale of fractal superstructures of aerogels can range from angströms to hundreds nanometers [20].Therefore, simulating polar molecular fluids and electrolyte solutions sorbed in such a disordered porous material is hampered by the need for a very large simulation box to catch correlations between pores of the material.The nanoporous host material is frequently represented in a simplified manner by a single pore of a spherical, cylindrical, slit-like, or other shape.Rovere et al. [21] performed molecular dynamics simulations for the structure of water confined in a cylindrical nanopore of Vycor glass.Shevade et al. [22] employed grand canonical Monte Carlo (GCMC) simulations to determine the structure and the adsorption isotherms of water-methanol mixture in a slit activated carbon micropore with carboxyl groups periodically grafted to its surface.Cascarini de Torre and Bot-tani obtained the isotherms of adsorption of nitrogen [23] and ethylene fluid [24] on a heterogeneous solid surface exhibiting the adsorption properties similar to real carbonaceous materials.Davies and Seaton [25] proposed a procedure to represent activated carbons as a size distribution of pores of regular shape.They performed a number of GCMC simulations for adsorption of mixtures of methane and ethane in a single slit-like pore of various size, and then determined the pore size distribution to fit the overall adsorption isotherm to experimental data.King and Do [26] described the effect of adsorption centers in activated carbon mesopores on the adsorption isotherms of such polar fluids as water by a semi-phenomenological theory based on a simple picture of an adsorbed phase droplet growing on a single adsorbing spot at the planar pore wall.
A crucial property of aerogel materials is disorder on a microscopic scale, which has a profound effect on the properties of sorbed fluids.MacElroy and Raghavan [27] studied dilute Lennard-Jones (LJ) vapor confined in microporous silica modelled in two ways as a network of randomly interconnected, and randomly distributed silica microparticles.Kaminsky and Monson [28,29] developed a model for adsorption of LJ fluid in microporous silica gel comprising a configuration of microspheres, each considered as a continuum of atomic interaction centers.They derived an analytic interaction potential between such a composite sphere and an atom of sorbed LJ fluid which is an analog of the 9-3 potential model for adsorption on planar solid surfaces.They used GCMC simulations to determine the adsorption isotherms and structure of methane [28] and methane-argon mixtures [28] sorbed in the porous microstructures represented by an equilibrium hard-sphere configuration and a facecentered cubic lattice of adsorbing composite spheres, and demonstrated that the Henry constant and the total amount sorbed are quite sensitive to the fine geometry of pores.The adsorption selectivity of this mixture strongly depends on the molecular size differences and heterogeneity of the adsorbent-fluid potential, resulting in such complex phenomena as the occurrence of azeotropes.Gordon and Glandt [30] employed the Kaminsky-Monson potential in GCMC simulations for adsorption of n-alkanes on a model microporous silica gel, and obtained the adsorption isotherms and a change of adsorption behavior from continuous pore filling by n-octane to layering of n-dodecane molecules at the surface of micropores.Furthermore, Gordon and Glandt [31] modified this model by incorporating either discrete dipoles or a continuum dipole bilayer at the composite sphere surface.Segarra and Glandt [32] modelled microporous activated carbons by an isotropic collection of circular graphitic platelets with dipoles evenly smeared around the edges of platelets to represent surface functional groups attached to the edges of the carbon basal planes.However, they used this advanced model of microporous silica and carbon to simulate the adsorption of methane, ethane, and water only in gas phase.
Very recently, Gavalda et al. [33] have developed a realistic molecular model of nanoporous carbon aerogel and simulated the adsorption isotherms of nitrogen fluid.The connected network of carbon nanoparticles is prepared as a random closepacked structure of slightly overlapping spheres matching a given porosity, and each sphere is represented by a model microporous carbon generated by means of the Reverse Monte Carlo method to fit the X-ray diffraction data for the carbon radial distribution function.Even for such a relatively simple nonpolar molecular fluid as nitrogen, simulation of adsorption in this nanoporous model required a large simulation box of length 36 nm, and hence a very large amount of computational time and the use of a massively parallel GCMC simulation algorithm [33].
An alternative is offered by integral equation theory of liquids which yields a cost-efficient description as well as analytical relations for the fluid properties in the thermodynamic limit [34].A number of works employed the density functional theory for adsorption of LJ fluid in a single slit-like micropore, with further account of the pore size distribution [35] and pore wall heterogeneity [36].Lozada-Cassou and co-workers developed integral equation theories for the structure of the restricted primitive model of electrolyte solutions confined between two plates [37] and in a cylindrical pore [38], as well as between two charged plates [39].Zaini et al. [40] extended the study to electrolyte solutions in a charged cylindrical pore.Yang et al. [41] modelled electrosorption of Na + and F − ions in carbon aerogel from aqueous solution at concentrations up to 0.1 M by treating the carbon aerogel electrodes as planar electrical double-layer capacitors.They used the classical double-layer theory based on the Poisson-Boltzmann equation regarding ions as point charges and solvent as structureless dielectric medium.Besides missing the effects of the ions and solvent molecules sizes, the applicability of this approach to confined spaces is problematic.By comparison with simulations and integral equation theory, it was shown that the results of Poisson-Boltzmann theory are good for dilute monovalent electrolytes in slit-like pores, but become erroneous even at relatively low concentrations of electrolyte in small cylindrical pores with a large overlap of the EDL's [38].
A successful approach to fluids and solutions of complex molecular species in a wide range of density from gas to liquid is provided by the reference interaction site model (RISM) theory pioneered by Chandler and Andersen for nonpolar molecular fluids [42,43] and extended by Hirata et al. to charged species [44][45][46].Kovalenko and Hirata proposed a new closure [47] (below referred to as the KH approximation) which enabled the description of associating molecular fluids and solutions in the whole density range from gas to liquid [48,49].It is a non-trivial optimized coupling of the so-called hypernetted chain (HNC) approximation and the mean spherical approximation (MSA).The former is enabled for repulsive core regions, whereas the latter is automatically applied to long-range enhancement tails of the distribution functions in the critical regime as well as to their high association peaks wherever they arise.The KH closure approximation is adequate for description of charged polyatomic fluids and mixtures over a wide range of density with proper account for their chemical specificities.The RISM/KH theory qualitatively reproduces the vapor-liquid phase diagrams of associating molecular fluids such as water, methanol, and hydrogen fluoride, and their structure in the gas as well as liquid phases, including hydrogen bonding [48,49].It appropriately describes the association structures that sequentially change in a mixture of water and amphiphilic tert-butyl alcohol with their molar fraction, including micromicelles [50].In the context of the present work, the KH approximation is suitable (i) to model carbon aerogel as a nanoporous network of interconnected carbon nanoparticles with both short-and long-range order in their correlations, and (ii) to describe the correlations of aqueous electrolyte solution in charged nanoporous aerogel which are of coupled long-range electrostatic and short-range chemical character for this strongly associating system.Special statistical-mechanical means are required to take proper account of the disorder of the porous host material (usually called "matrix") which comprises immobile particles and thus is not in thermodynamic equilibrium with the sorbed fluid.By using a cluster diagram analysis, Madden and Glandt (MG) [51,52] derived the Ornstein-Zernike (OZ) type [34] integral equations for atomic simple fluid confined in a disordered matrix of spherical obstacles frozen in space.Chandler [53] proposed a RISM generalization of the MG theory to the case of polyatomic fluids sorbed in quenched amorphous materials, in particular a polymer RISM [54][55][56] generalization to fluid of flexible polyatomic molecules sorbed in a quenched matrix.Thompson and Glandt [57,58] employed the polymer RISM-MG theory to describe the structure and thermodynamics of a fluid of freely jointed hard sphere chains in a hard sphere matrix.In the MG approach, the matrix is regarded as a single rigid supermolecule immersed in the fluid [53].However, this approximation does not distinguish between quenched and annealed matrices [59] and thus neglects an additional part of the direct correlation functions between fluid particles arising due to the presence of immobile matrix obstacles.Computer simulations show that these so-called blocking direct correlations become particularly essential for quenched-annealed systems including charged species [60].Employing the so-called replica identity to treat the free energy of the sorbed annealed fluid averaged over the distribution of the quenched matrix, Given and Stell [59,[61][62][63] developed the formalism of the replica OZ equations for fluid and matrix of atomic species.It takes proper account of the blocking correlations in a quenched-annealed system.Hribar et al. [68] used the replica OZ/HNC theory versus simulations to describe 1-1 to 4-1 electrolyte solutions with primitive model solvent confined in hard sphere matrix mimicking solutions of NaCl, Li-Cl, CaCl 2 , LaCl 3 , and ThCl 4 sorbed in microporous silica xerogel.Their calculations showed that electrostatic forces influence the blocking correlation functions which assume high values at small distances, increasing them for equally charged and decreasing for oppositely charged ions.The difference increases with the electrolyte concentration, even though the excluded volume effects contribute significantly at high concentration of the adsorbed electrolyte solution.A crucial merit of the replica integral equation method is that it provides a consistent analytical approach to the thermodynamics of sorbed fluid [64][65][66][67].The matrix distribution functions following from the OZ equation for matrix species are an input to the replica OZ equations for the structure of sorbed fluid.Therefore, the matrix structure can be allowed for in the form of the partial structure factors taken from diffraction experiments or simulation procedures as well as modelled by integral equation approaches.Krakoviack et al. [69] employed the replica integral equation method to study the phase diagram and the structure of LJ fluid sorbed in a highporosity aerogel.The structure factor of the latter was calculated for a connect-ed cluster of microspheres generated by the simulation process of diffusion-limited cluster-cluster aggregation mimicking the synthesis of silica aerogels (see literature in [69]).
Recently, we developed a RISM generalization of the replica method to the theory of molecular quenched-annealed systems comprising both fluid and matrix species described with realistic force fields, such as water sorbed in porous carbon aerogel activated with carboxylic groups [70,71].Along with the distribution functions, we derived closed expressions for the system thermodynamics, in particular the chemical potential and pressure.They acquire new terms in the case of the temperature of the sorbed fluid different from the temperature characterizing the spatial distribution of the quenched matrix species.An important advantage of the replica RISM method is that it can treat molecular fluids and solutions sorbed in disordered porous materials with pore size ranging from micro-to macroscopic scale.Based on the replica RISM equations complemented with the KH closure approximation [48,49] and the modified version [72] of the Verlet closure [73] (MV) applied to the blocking correlations, we considered the fluid phase equilibria of methanol as well as water sorbed in nanoporous carbon aerogel [49].The replica RISM-KH-MV theory predicts that methanol can condense either filling the whole nanopore volume or forming a thin layer physisorbed at the nanopore surface.In the latter case, zigzag chains of hydrogen bonding of methanol bend to adjust to the surface curvature.By contrast, the calculations yield liquid water only continuously filling the nanopores.Apparently, the physisorption on the surface of carbon aerogel nanopores is not strong enough to distort the tetrahedral network of water hydrogen bonds which is much more resistant than zigzag chains of methanol hydrogen bonding.These our findings are in agreement with experimental observations for methanol and water sorbed in microporous carbons [74,75].The latter revealed that methanol molecules can interact with the micropore walls more strongly than water, and form a solid-like localized structure with dimensionality less than three and considerably long-range order not observed in the bulk liquid [74].In contrast, water molecules cannot form uniform adlayers on the hydrophobic pore walls, but are presumed to arrange into small unit clusters up to pentamers uniformly filling the hydrophobic space [75].
Below we employ a dielectrically consistent modification of the replica RISM-KH-MV theory to describe the electrochemistry of a double layer formed by ambient electrolyte solution of sodium chloride at concentration 1.069 mol/l sorbed in nanoporous carbon aerogel with an external electric charge.The aerogel is modelled as a network of partially overlapping compact carbon nanospheres.To illustrate the effect of the aerogel nanostructure, it is represented by two models which are reminiscent of those used by MacElroy and Raghavan [27]: nanospheres that are almost randomly distributed, and nanospheres associating into a network of branched chains.As shown below, the long-range correlations of aerogel nanoparticles strongly affect the structure and thermochemistry of the EDL formed by the sorbed electrolyte solution at the surface of carbon aerogel nanopores.

Replica DRISM-KH-MV integral equations
In the replica RISM approach to the structure and thermodynamics of a quenched-annealed molecular system [49,70,71], the average free energy of molecular fluid (species 1) sorbed in a disordered porous matrix of quenched molecular particles (species 0) is obtained as a statistical average od the fluid free energy at matrix spatial configuration q 0 over the ensemble of matrix realizations, where T 1 is the temperature of fluid.Statistical averaging of a logarithm is difficult, and is obviated by using the so-called replica identity relating the logarithm to the analytic continuation of moments, ln Z 1 = lim s→0 dZ s 1 /ds.For integer s, the moment average takes the form of the equilibrium canonical partition function of a fully annealed (s + 1)-component mixture which comprises N 0 particles of the matrix species and s equivalent replicas of the liquid, each consisting of N 1 particles, with no interaction between the replicas.The average free energy of the mobile molecular fluid is then obtained as the analytic continuation from the free energy A rep (s) of the annealed replicated system, Accordingly, the derivatives of this free energy expression are taken analytically and give consistent access to the thermodynamics of the sorbed fluid [62][63][64]67].
The analytic continuation of the RISM integral equations for the annealed (s + 1)component replicated system to the limit of s → 0 in the assumption of replica symmetry persistence leads to the replica RISM integral equations which yield the radial distribution functions g ij αγ (r) dependent on separation r between interaction sites α and γ of molecular species i, j constituting the quenched-annealed system.Their derivation is given in detail in [70,71].To provide an adequate description of solution consisting of ions at an arbitrary concentration in polar molecular solvent, Perkyns and Pettitt [76,77] developed the dielectrically consistent RISM theory (DRISM) including a so-called bridge function providing consistency of the solution dielectric properties.We adapted the DRISM approach as most appropriate for the system in which charged and polar molecular species are essential, and introduced the replica DRISM integral equations [78].
Written for the total and direct correlation functions, an integral equation of liquid state theory needs an additional, so-called "closure" relation between them.The exact closure is nonlocal and extremely cumbersome, and therefore is replaced with one of the approximations known in liquid state theory [34].Physical phenomena predicted by a theory stem from singularities and asymptotics of a closure approximation.To complement the integral equations for the matrix-matrix, matrix-fluid, and fluid-fluid correlations, we employed the recently developed KH closure approximation [48,49] In the context of the present work, it is suitable (i) to model carbon aerogel as a nanoporous network of interconnected carbon nanoparticles with both short-and long-range order in their correlations, and (ii) to describe the correlations of aqueous electrolyte solution in charged nanoporous aerogel which are of coupled long-range electrostatic and short-range chemical character for this strongly associating system.For simplicity, the temperature of the matrix distribution is assumed to be the same as that of sorbed solution.Different temperatures of the matrix and sorbed fluid lead to new terms in the system thermodynamics [70,71] and can result in non-trivial effects in quenched-annealed systems with charges [60].The use of such linearized closures as the MSA for the blocking correlations neglects the blocking direct correlation function [62][63][64]67].To properly allow for the nonlinearity of the blocking correlations, the integral equation for them is complemented with the closure [78] in the Verlet functional form [73] with the modification to its denominator proposed by Kinoshita et al. [72] (MV closure).

Modelling of nanoporous carbon aerogel
We represent nanoporous carbon aerogel by the spatial equilibrium distribution of amorphous carbon composite nanospheres of size σ 0 = 32 Å which can overlap up to half the diameter [49].Their interaction potential comprising two Yukawa terms u 00 (r) = 4 00 exp 2 with the depth, size, and width parameters 00 = 17 kcal/mol, d 00 = 15.3Å, and δ 00 = 1 Å, forms a narrow association well at r (min) 00 = d 00 + δ 00 ln 2 ≈ 16 Å.We solved the OZ integral equation with the KH closure and the interaction potential (4) at the temperature T = 298 K and the number density of nanospheres ρ 0 = 7 × 10 −6 Å−3 .The nanospheres associate into a network of branched chains with fractal dimensionality d ≈ 1.5 persisting up to a length of 100 Å and the long-range correlations extending up to 200 Å [49].Previously, we considered a similar fractal structure consisting of interconnected branched chains of associating carbon-like atoms to model microporous carbon material [70,71].Figure 1 depicts the excess ∆n 00 (r) = n 00 − n (unif) 00 of the running coordination number of nanospheres n 00 (r) over the uniform distribution n different chains overlap and interconnect into a network although the long-range correlations in the network persist up to 200 Å.The average number of nanospheres in a cluster in excess over their uniform density thus amounts to ∆n 00 (r clust ) ≈ 16.With distance, the excess coordination number ∆n 00 (r) saturates at a value of approximately 34 extra nanospheres over the uniform distribution.
To reveal the effect of the long-range correlations of associating nanoparticles, we also used another aerogel model with repulsive nanospheres.The attractive part of the interaction potential as defined by the Weeks-Chandler-Andersen prescription [79] is removed from the expression (4) which now contains pure repulsion between nanospheres of size d ( At this density, the radial distribution function of repulsive nanospheres is almost random except for their exclusion volume: ).These almost randomly distributed nanospheres overlap much less than those associating into chains, which considerably reduces the volume of aerogel pores.To keep the porosity equal in the two aerogel models, the number density of non-associating nanospheres is somewhat decreased to ρ 0 = 5.14 × 10 −6 Å−3 .
Porosity and pore size distributions of nanoporous materials are usually tested and characterized by adsorption of such chemically neutral and rigid species as rare gases, nitrogen, or mercury.A simple and reasonable theoretical definition of porosity is a volume fraction of aerogel cavities available for insertion of a test hard sphere into the aerogel ensemble of nanoparticles regarded as hard bodies, p ≡ V cav /V 0 .For the present aerogel models, these are hard nanospheres with effective size σ 0 .According to Widom's test particle method [80], the volume portion V cav available for insertion of a hard sphere into the ensemble of hard particles of the liquid occupying the volume V 0 is related to its chemical potential, µ hs = − ln(V /V 0 ).Relevant in this theorem is that the test hard sphere is being inserted into an ensemble of particles, that is into a frozen spatial distribution representing the equilibrated liquid.This exactly corresponds to the modelling of the aerogel by an ensemble of nanospheres which are first equilibrated and then frozen in space.The aerogel porosity p as tested by hard spheres of size σ hs is thus calculated simply as The pore size distribution of the aerogel (the portion of aerogel cavities fitting hard spheres of size σ hs in the size range dσ hs ) can be determined by differentiation with respect to the size of the test hard sphere, v(σ hs ) = dp/dσ hs .
The Percus-Yevick (PY) closure is known to be a good and quite reliable approximation for hard body molecules.To calculate the chemical potential µ hs (σ hs ), we solve the solute-solvent OZ/PY integral equations for a hard sphere particle at infinite dilution immersed in an equilibrated fluid of hard composite spheres with the given radial distribution function.Figure 2 draws a comparison between the pore  size distributions obtained for the two aerogel models.Both are characterized with the same porosity value of p ≈ 89% at the test particle size σ hs = 3 Å approximately corresponding to the water molecule.Since the separations between composite nanospheres associated into clusters are smaller, the gaps between clusters are larger than those between nanospheres distributed at random.This increases the most probable size of aerogel cavities from 30 Å in the model with repelling composite nanospheres to 42 Å in the aerogel with associating nanospheres, and extends the maximal size of cavities from 100 to 170 Å.
The specific surface area s(σ hs ) of aerogel pores accessible for a test particle of size σ hs can be obtained from a change of the void volume fraction V cav /V 0 with "covering" the surface of aerogel nanoparticles by a thin hard layer.This is equivalent to varying the nanosphere radius R 0 = σ 0 /2 in the aerogel-solution potential (in the general case, simultaneously for all sorts of nanospheres) at the distribution of nanospheres g 00 (r) and the probe size σ hs kept unchanged, which gives Calculation of the surface areas as tested by a 3 Å hard sphere yields a value of s = 66 m 2 per cm 3 of aerogel with nanosphere chains, and s = 88 m 2 /cm 3 for aerogel with random nanospheres.The larger surface area of pores is related to the decrease of overlap of aerogel nanospheres in the latter model.Given the physical density of carbon in the body of an aerogel nanoparticle, the porosity and surface area recalculated per gram of carbon aerogel material are listed below in table 3.
A carbon nanosphere of radius R 0 = σ 0 /2 = 16 Å is treated as an aggregate of uniformly distributed carbon atoms with the 12-6 LJ parameters σ C = 3.4 Å and ε C = 0.0556 kcal/mol, and the number density ρ C corresponding to the physical density of graphite, 2.27 g/cm 3 .This leads to the Kaminsky-Monson potential [28,29].Besides, the interaction potential between interaction site α of solution species and a carbon nanosphere contain the Coulomb interaction of the partial site charge Q 1 α with the carbon nanosphere charge Q 0 induced by the external electric source.The induced charge is considered to be centered on the nanosphere.We neglect electrostatic images of solution species charges at the nanoparticle surface.Simulations reveal a very large degree of cancellation of the surface polarization effects altogether in water near a planar electrode, although the image potential of a single water molecule can be quite large [81].On the other hand, electrostatic image potentials at short range make little sense for a rough surface of a carbon composite nanoparticle.We postpone consideration of such effects which requires a self-consistent explicit description of the electronic structure of amorphous carbon nanoparticles.Thus, the total potential between solution species and a carbon nanosphere is written as The site-site interaction potentials between solution species sorbed in the aerogel are modelled by the Coulomb and LJ terms, The LJ diameter and energy parameters between unlike species or sites in the potentials ( 9) and ( 10) are determined by the standard Lorentz-Berthelot mixing rules respectively as σ ab αγ = (σ a α + σ b γ )/2 and ε ab αγ = (ε a α ε b γ ) 1/2 .Water molecules are represented with the simple point charge (SPC) model [82], and the LJ size and energy of the Na + and Cl − ions are taken from the simple ions parametrization summarized in [83].The interaction site parameters, including carbon atoms of the nanoporous aerogel matrix are shown in table 1.

EDL potential
An external electric source induces opposite specific charges ±q ex on the two carbon aerogel electrodes of the supercapacitor, q ex = Q 0 ρ 0 .They attract ions of the corresponding sign to the EDL inside the nanoporous electrodes, and cause separation of cations and anions between the two electrodes until the bias of their concentrations inside each electrode satisfies the local electroneutrality condition The statistically averaged distribution of charge density around a carbon nanosphere due to both aerogel matrix and sorbed solution species, creates the average electrostatic potential ψ 0 (r) satisfying the Poisson equation Neglecting the effects of the surface corrugation and internal microstructure of a carbon nanosphere and regarding it as a classical uniform conducting sphere of radius R 0 with charge Q 0 , the total average electrostatic potential around the nanosphere is written as The value of the electrostatic potential ( 14) inside the carbon nanosphere, φ C (q ex ) = φ 0 (r R 0 ; q ex ), corresponds to the potential level of the carbon conducting framework of the nanoporous electrode.Its value at a distance from the nanosphere, φ av (q ex ) = φ 0 (r = ∞; q ex ), by definition gives the average potential level in the whole electrode including both carbon nanospheres and sorbed solution.It should be emphasized that the change from the nanosphere potential φ C (q ex ) to the electrode bulk level φ av (q ex ) comprises not only its own EDL voltage but also the electric field of all other surrounding carbon nanoparticles and their EDL's distributed with the spatial correlations of porous material.The average electrostatic potentials inside each of the two oppositely charged nanoporous electrodes with sorbed solution, φ av (q ex ) and φ av (−q ex ), are related by the macroscopic electric field in the permeable membrane separating the electrodes.The chemical potential of sorbed solution species s comprises the ideal gas contribution µ id s = k B T 1 ln (ρ s Λ s ), the excess term ∆µ s due to the microscopic intermolecular interactions inside the electrode, and the energy in the macroscopic electrostatic field between the supercapacitor electrodes, where ρ s and Q s are the number density of species s and their total charge, and Λ s = (2π 2 /(m s k B T 1 )) 1/2 is the de Broglie thermal wavelength of ideal monatomic particles with molecular weight m s .The rotational, vibrational, and internal electronic terms omitted in the ideal term µ id s (ρ s ) are implicitly included into the reference potential level.They should be taken into account to obtain the values of the full Gibbs free energy to be compared with experimental data.Notice that in the higher-order approximation, these terms can be coupled with the spatial structure of solution.In general, the excess chemical potential ∆µ s (q ex ) of each solution species depends on the nanoporous electrode charge.It is strongly different for ions sorbed in the two oppositely charged electrodes.This causes diffusion of ions creating a gradient of the ionic concentrations in the permeable membrane between the electrodes until it is balanced by the ideal terms difference as well as the electric field of the macroscopic diffuse double layer in the membrane to equalize the chemical potentials in the two electrodes, The potential difference (16) does not depend on the effects at the boundaries of the nanoporous electrodes, in particular on the structure of electrode pores openings towards the permeable membrane and on the membrane shapes.The chemical potentials of sorbed solution are bulk characteristics of the nanoporous electrode, since the relative contribution of the boundary to the free energy vanishes for the electrode of macroscopic size.The relation ( 16) must be satisfied with the same average electrostatic potential bias φ I av −φ II av for each ionic species, subject to the electroneutrality condition (11) for the solution densities.A discrepancy between this bias for cations and anions causes disbalance of the chemical potentials and simultaneous diffusion of both cations and anions, subject to (11), to the electrode with the lower chemical potential.Similar transport can occur for neutral solvent molecules with Q s = 0 if their excess chemical potentials in the two electrodes do not satisfy the condition (16).This process of electroosmosis continues until the excess chemical potentials ∆µ s changing with the solution being "pumped" between the electrodes meet the balance ( 16) for all solution species.
The supercapacitor voltage is found by adding up the electrostatic potential differences inside each of the nanoporous electrodes ∆φ I C (q ex ) = φ I C − φ I av and ∆φ II C (q ex ) = φ II C − φ II av calculated from the distribution functions by equations ( 12)-( 14), and that between the average levels of the electrodes φ I av − φ II av given by equation (16), The integral specific capacitance calculated per volume of one nanoporous electrode of the supercapacitor reads c = 2q ex /U (q ex ).The expression (17) reveals the effect driving the voltage and hence specific capacitance of a nanoporous EDL supercapacitor.It appears that a significant contribution to U (q ex ) comes from the change of the excess chemical potentials of ions ∆µ s (q ex ) with charging the electrodes.The latter are determined jointly by the specific geometry of the EDL inside a nanoporous electrode and the short-range structure of solution species.The excess chemical potentials of the fluid sorbed in a quenched matrix split up into the matrix-fluid, fluid-fluid, and blocking terms [49,70,71].The first two are obtained in the KH approximation which is applied to these correlations.An important feature of the KH closure to integral equation theory is that it has an exact differential of the free energy functions [47][48][49].The considerably smaller contribution of the blocking correlations described in the MV approximation can be obtained in an analytical form too, using the unique functionality assumption [84].This yields the excess chemical potentials of solution species ∆µ s in a closed analytical form in terms of the total and direct correlation functions of nanoporous carbon and sorbed solution [49,70,71,78].Their values to be inserted into equation (17) are conveniently calculated by integrating the radial correlation functions of the system on solving the replica DRISM-KH-MV integral equations.

Structure of the sorbed solution
In this work, we considered the aerogel filled with aqueous solution of sodium chloride at temperature T = 298 K and concentration 1.069 M in the bulk.At zero charge, the number densities of sorbed solution species ρ External charge ±q ex imposed on the two aerogel electrodes of the supercapacitor attracts the corresponding ions of opposite charge to ensure the electroneutrality condition (11), and causes separation of sodium and chlorine ions of the sorbed solution by diffusion through the semipermeable membrane between the electrodes.For the external charge per unit volume q ex = ±0.9613C/cm 3 , we specify the ionic concentrations as biased respectively to ρ 1 Na+ (q ex ) = (1 ∓ 0.005236)ρ 1 Na+ (0) and ρ 1  Cl− (q ex ) = (1 ± 0.005236)ρ 1 Cl− (0).In general, due to the different sizes and solvation structures of the cations and anions this should leads to a change in the water solvent density too, subject to the condition of constant pressure of the sorbed solution.We imply this effect to be of a higher order of magnitude and postpone its consideration to further development.
We solved the replica DRISM-KH-MV integral equations with the interaction potentials ( 9), (10), and (4) or (5) for either of the aerogel models on a radial grid of 8192 points with a linear resolution of 0.1 Å.The convolutions are calculated by multiplication in reciprocal space using the fast Fourier transform (FFT) technique, with the Coulomb singularities of the correlation functions separated out and treated analytically [86].We converged the equations on a Pentium PC to a root mean square accuracy of 1 × 10 −10 by using the modified direct inversion in the iterative subspace (MDIIS) method [86,87].As compared to the procedure of Picard iterations (or damped iterations), the MDIIS provides drastic acceleration of solving these highly nonlinear integral equations.
Figure 3 shows the site-site radial distribution functions (RDF's) g 11 αγ (r) between the species of the NaCl aqueous solution sorbed in the nanoporous aerogel with branched chains of carbon composite nanospheres.The peaks of the RDF's reflect the hydration structure of ions and the formation of ionic pairs and ionic clusters consisting of several ions of the opposite sign.They are similar to those observed in the bulk ambient aqueous solution of NaCl at 1.069 M concentration which are discussed in detail in [88].Highest is the first peak of the Na + and Cl − ions associating in a contact ionic pair.The second peak of the Na + -Cl − RDF can be attributed both to an ionic pair separated by water molecules of their hydration shell, and to the farthest opposite ions in a cluster comprising three consecutively flipped ionic pairs: ⊕ ⊕ ⊕ .Unlike the bulk solution, the Na + -Cl − RDF has a noticeable third peak which might be related to the formation of even larger clusters of opposite ions layered at the surface of a carbon aerogel nanoparticle.The first peak of the Cl − -Cl − RDF is substantially higher at a finite concentration of NaCl than at infinite dilution.It corresponds to the Cl − ions in a cluster comprising several pairs of opposite ions in  the former case, and to a single hydrated ionic pair Cl − -Cl − in the latter.In both cases water molecules make hydrogen bonds to the Cl − ions in contact and form stabilizing bridges between them, but the stability of the Cl − ions in the cluster is enhanced by attractive bridges through Na + ions.Because of the smaller size, Na + ions included in the cluster are preferably separated by Cl − ions in contact.Similarly, a single hydrated ionic pair Na + -Na + is preferred in the water-separated rather than contact arrangement.Therefore, the first peak of the Na + -Na + RDF is somewhat lower than the second one.The peaks of the ion-water oxygen and hydrogen RDF's recapitulate the structures of the ion hydration shells.As is well known, water molecules associated to a Na + ion orient with the oxygen towards and their dipole moments straight outwards the ion, whereas water molecules around a Cl − ion are tilted forming hydrogen bonds to the ion with one of their hydrogens.The peaks of the water-water RDF's are similar to those of the bulk NaCl solution.
We do not depict the distributions g 11 αγ (r) of the solution sorbed in the second model of carbon aerogel with randomly distributed carbon composite nanospheres, since they are visually close to those in figure 3. The essential difference consists in the shorter range of their asymptotic tails.The connected part of the site-site correlation functions of the sorbed solution solely determines its compressibility [64,70,71] which is in fact a local physical property of a solution portion confined in an aerogel cavity, and primarily shapes its short-range chemical structure, with such chemical specificities as hydrogen bonding.For both the aerogel models, the connected part of the total correlation functions of the solution, h (c) αγ (r), falls off by two orders of magnitude already at r ≈ 13 Å and quickly decays with distance.On the other hand, their blocking part h (b) αγ (r) bears long-range correlations "imprinted" in the structure of the sorbed solution by those of the porous matrix material.Figure 4 draws a comparison between the blocking total correlation functions of the sorbed solution for the aerogel models with the long-range order of branched chains of carbon nanospheres, and with randomly distributed nanospheres.The sitesite correlations run rather close to each other, although showing some splitting into the ion-ion, water-ion, and water-water groups for the former aerogel model.Important is a large difference in the range of the blocking correlations which last to a characteristic correlation length of the corresponding matrix structure.It is the clustering size 100 Å for the former aerogel model (figure 1), and the carbon composite nanosphere size σ 0 = 32 Å for the latter.Thus, the blocking parts h (b) αγ (r) constitute the long-range tails of the distributions g 11 αγ (r) and can be interpreted as interpore correlations of sorbed fluid.They manifest in the structure of sorbed solution measured in diffraction experiments which show strong enhancement of the intensity in a small angle scattering region, typically below 0.5 Å−1 (see discussion in [70,71]).These interpore correlations also affect the structure of the EDL formed by the sorbed solution at the surface of carbon aerogel nanoparticles.
Figure 5 shows the RDF's of sorbed solution species around a carbon composite nanosphere of the aerogel electrode.For the aerogel with branched chains of carbon nanospheres, the nanosphere-solution RDF's acquire the long-range depletion tails extending up to 100 Å (displayed in part to better resolve the short-  range structure).This can be ascribed to the exclusion effect of clustering of aerogel nanospheres resulting in the enhancement of their local density.In the aerogel of random nanospheres, the nanosphere-solution RDF's quickly decay at the distance larger than the nanosphere size, r > 32 Å.The short-range peaks of the RDF's for the two aerogel models differ in height but are qualitatively similar.Their positions reveal the peculiarity of the adsorption character of the cation and anion and their hydration structure at the carbon nanosphere surface.Near the surface, the carbon nanosphere-ion RDF's develop oscillations with alternating peaks of sodium and chloride ions.Figure 6 sketches the ion adsorption arrangements corresponding to the first maxima of the nanosphere-ion potentials of mean force.We obtained that Cl − ions are preferably adsorbed in the inner Helmholtz layer, in contact with nanospheres, whereas Na + ions are found in the outer Helmholtz layer, separated from the nanosphere surface by water molecules of their hydration shell.Some fraction of Na + (smaller than Cl − ), however, comes in contact with the surface, forming a shoulder at the separation smaller than the first Cl − peak.Water molecules compete with chlorine ions for adsorption in the inner Helmholtz layer.These arrangements are in agreement with the commonly accepted picture [9] observed also in simulations of NaCl aqueous solution in a charged slit-like pore [89].The first peaks of the carbon nanosphere-water oxygen and hydrogen RDF's are close in position, with the latter half as high but wider.This corresponds to adsorbed water molecules with the hydrogens oriented preferably along the surface but also tilted both towards and outwards it, which is typical of a hydrophobic hydration shell.The short-range solvation structure is similar for the two aerogel models.However, in the aerogel with the long-range correlation of carbon nanospheres the peaks of the RDF's are reduced approximately by a factor of two along with their long-range depletion.Notice that in our theory, the preferential adsorption of Cl − does not involve any specific adsorption force, but follows from the interplay of the anion hydration structure and the steric effect of the hydrophobic carbon surface.

Potential of the electrochemical double layer
With imposing external charge q ext on the carbon aerogel electrode, the RDF's of solution ions around a carbon nanosphere respond by depletion or enrichment according to the signs of the ion charges, whereas the aerogel-water RDF's scarcely change.As is seen in figure 5, the deviation is most noticeable in the short-range peaks of the aerogel-ion RDF's, and is substantially larger for Na + than for Cl − .The hydration shell of adsorbed Na + ions thus appears to be softer than that of Cl − ions in the Helmholtz layer.The aerogel charge also affects the long-range tails of the aerogel-ion RDF's for the former aerogel model.The changes in the RDF's of solution species modify the average electrostatic potential created by the EDL around a carbon nanosphere.Figure 7 shows the run of the electron energy in the average electrostatic potential with distance from the nanosphere center, −eφ (es) α (r).The electrostatic potential oscillations stem from spatial and orientational ordering of ions and water molecules at the surface of carbon nanoparticles.At nonzero q ext , the first bending of the curves for r R 0 corresponds to the almost bare Coulomb potential between the charged carbon nanosphere surface and the charges of the centers of solution species in the inner Helmholtz layer.Then follows the strong potential rise in the surface dipole layer formed by water oxygens in the inner Helmholtz layer and their hydrogens oriented towards the surface, and on the other hand, the surface dipole of ions in the inner Helmholtz layer, the Cl − ions and a portion of the smaller Na + ions in contact with the surface.The next potential drop almost compensating the increase originates from the opposite surface dipole created by oxygens of adsorbed water molecules and their hydrogens directed outwards the surface, as well as by Cl − ions in the inner and Na + ions in outer Helmholtz layers (see figure 5).The two consecutive diminishing oscillations originate from the ordering of the second and the third solvation shells around carbon nanospheres.
The external charge induced on carbon nanoparticles of the aerogel electrode causes a bias of the ionic concentrations in aerogel pores, and causes long-range tails of the EDL potential.The long-range ionic distributions in the aerogel with branched chains of carbon nanospheres lead to the long-range tail of the EDL potential around a carbon nanosphere, even at the potential of zero charge.The long-range tail is absent in the non-charged aerogel with random nanospheres.Notice that this region of the ionic concentration bias and the electrostatic potential around a carbon nanosphere slowly decaying up to a distance of hundred nanometers should not be confused with a diffuse layer.In an EDL with the slit-like geometry, the diffuse layer thickness is characterized with the Debye length which amounts in ambient aqueous solution of NaCl at 1.069 M to about 3 Å.As distinct, the long-range tails of the RDF's and the electrostatic potential of solution species around an aerogel nanosphere result from the statistical and orientational averaging of its EDL as well as the EDL's around all other nanoparticles forming nanopores of various shapes.It is clear from figure 7 that the response of the EDL potential on the external charge, φ C (q ex ) − φ C (0), comprises the contributions from the inner and outer Helmholtz layers, and from the aerogel nanosphere-solution ion correlations at long range.The latter essentially dominate the EDL potential response in the carbon aerogel with branched chains of nanoparticles.In contrast, with the random distribution of nanoparticles, the EDL response is essentially determined by the electric field in the Stern layer.
To satisfy the electroneutrality (11) for the nanoporous carbon aerogel electrode with charge q ex = 0, the densities of Na + and Cl − ions were biased symmetrically, neglecting at the moment the electroosmosis effects which are of a higher order of magnitude.As is seen from table 2, the changes in the chemical potentials of Na + and Cl − with charging the electrodes are rather similar, and that of water is comparably small.Worthwhile noting also is that the change in the excess chemical potential ∆µ s is substantially larger than in the ideal term µ id s .That is, osmosis effects are of a higher order of magnitude as compared to the electrostatic and short-range structure forces.An additional EDL emerges in the solution at the macroscopic boundary of each electrode to counterbalance the difference between the ion chemical potential µ id s + ∆µ s inside each of the two electrodes.These macroscopic boundary EDL's do not have effect on the supercapacitor voltage, but can influence its charge/discharge characteristics.A salient feature of a nanoporous carbon aerogel supercapacitor is that its voltage U (q ex ) comprises the response of not only the EDL potential step φ C −φ av at a particular carbon nanoparticle, but also of the solution chemical potentials µ id s +∆µ s inside the nanoporous electrode which are a collective effect of all surrounding nanopores filled with the electrolyte solution.These two factors add up in the supercapacitor voltage given by equation (17).A comparison of the EDL shifts in figure 7 with the differences of the chemical potentials in table 2 shows that the latter constitute half the supercapacitor voltage in the carbon aerogel electrodes with randomly distributed nanoparticles.The integral specific capacitance calculated per gram of one nanoporous electrode of the supercapacitor is given in table 3. It follows from the above decomposition of the supercapacitor voltage that the long-range correlations between carbon aerogel nanoparticles decrease the specific capacitance of the supercapacitor electrode due to the dominating long-range tails introduced into the EDL  5).As is seen in figure 1, the carbon aerogel with the long-range correlation of carbon nanoparticles is characterized with their enhanced clustering.The capacitance decrease might be attributed to the overlap of EDL portions between densely clustered carbon nanoparticles, where nanopores shrink to micropores.The latter constitute a significant portion of the pores surface in this aerogel model, because the pore surface area per its volume increases for micropores as it scales inversely to the pore size.Such an EDL overlap strongly affects the EDL properties [38,40], and thus diminishes the effective surface area of the EDL in micropores.

Conclusion
We have developed the replica RISM integral equation theory for molecular quenched-annealed systems.With our KH closure approximation, the replica RISM theory is suitable for realistic modelling of electrolyte solutions sorbed in disordered porous materials with pore size scale ranging from a few angströms to hundreds of nanometers.This microscopic theory allows for chemical functionalities of both sorbed solution species and porous host material, including activating chemical groups or catalytic centers.It properly and consistently accounts for the effect of disorder of porous host material on the structural, thermodynamic, and electrochemical properties of sorbed electrolyte solution.To illustrate the theory, we used two carbon aerogel models that comprise composite nanospheres of amorphous carbon, either randomly distributed or clustered into a network of branched chains with fractal dimensionality 1.5 up to a size of 100 Å.We calculated the structure and the specific capacitance of the electrochemical double layer in a nanoporous carbon aerogel electrode filled with ambient aqueous solution of sodium chloride.We also showed the effect of the long-range order of nanoporous aerogel on the properties of the electrochemical double layer.
We aimed in this work to demonstrate that this theory is capable of microscopically describing an electrochemical double layer in realistic nanoporous electrodes.The specific capacitance obtained exhibits the magnitude comparable to that achieved in carbon aerogel supercapacitor devices [9], which is millions times greater than that of the conventional dielectric capacitor.Further development will involve elaborated carbon aerogel models with such characteristics as porosity, pores surface area, and the pore size distribution maximum fitted to experimental data.It will also include activating chemical groups and catalytic centers grafted on to the surface of aerogel pores, which significantly contribute to pseudocapacitance.Moreover, it is straightforward to perform calculations for a number of molecular electrolytes of various size and chemical properties.
ed ch ain s o f carb on n a n o sp h eres ran d o m ly d istrib u ted carb on n a n o sp h eres

Figure 2 .
Figure 2. Distribution v(σ hs ) of aerogel pores size (volume fraction of pores per unit range of pore size) as measured by insertion of a hard sphere of size σ hs .Aerogel with interconnected branched chains of carbon composite nanospheres (solid line) and with almost randomly distributed nanospheres (dashed line).

1 α
are set to those in the bulk ambient solution [85] reduced by the factor of matrix porosity p taken at test particle size of 3 Å.For p = 0.890 in the first model of aerogel with branched chains of nanospheres, this gives the values of water density ρ 1 water = 2.912 × 10 −4 Å−3 and equal ionic densities ρ 1 Na+ = ρ 1 Cl− = 5.730 × 10 −4 Å−3 .The second model of aerogel with almost randomly distributed nanospheres has slightly smaller porosity p = 0.8875, resulting in the corresponding adjustment of the solution densities.

Figure 3 .
Figure 3. Radial distribution functions g 11 αγ (r) between species of ambient aqueous NaCl solution at concentration 1.069 Mol sorbed in nanoporous carbon aerogel of porosity 89% comprising interconnected branched chains of carbon nanospheres.

Figure 4 .
Figure 4. Interpore (blocking) total correlation functions g (b) αγ (r) between species of the NaCl aqueous solution sorbed in nanoporous carbon aerogel that comprises carbon composite nanospheres forming interconnected branched chains (upper part) and nanospheres distributed almost randomly (lower part).

Figure 5 .Figure 6 .
Figure5.Radial distribution functions of solution species around a carbon nanosphere of radius R 0 = 16 Å for NaCl aqueous solution sorbed in a nanoporous carbon aerogel electrode with the specific charge q ext = 0.96; 0; -0.96 C/cm 3 (dash-dotted, solid, dashed lines, respectively).Carbon nanospheres clustered into branched chains (upper part) and almost randomly distributed (lower part).

Figure 7 .
Figure 7. Statistically averaged EDL potential around a carbon nanosphere for the NaCl aqueous solution sorbed in the charged aerogel electrode.The models of nanoporous carbon aerogel and nomenclature of lines are the same as in figure 5.

Table 1 .
Interaction site charges and Lennard-Jones potential parameters of NaCl aqueous solution sorbed in carbon aerogel.

Table 2 .
Chemical potentials µ id s + ∆µ s of the solution species sorbed in a nanoporous carbon aerogel electrode with specific electric charge q ex .

Table 3 .
Properties of the nanoporous aerogel electrode models.For the carbon density in a nanoparticle equal to the graphite density 2.27 g/cm3, and the aerogel porosity measured by a zero size probe.b As measured by a hard sphere probe of size σ hs = 3 Å.potential step (upper part in figure a