Size and charge correlations in spherical electric double layers: a case study with fully asymmetric mixed electrolytes within the solvent primitive model

Size and charge correlations in spherical electric double layers are investigated through Monte Carlo simulations and density functional theory, through a solvent primitive model representation. A fully asymmetric mixed electrolyte is used for the small ions, whereas the solvent, apart from being a continuum dielectric, is also treated as an individual component. A partially perturbative density functional theory is adopted here, and for comparison, a standard canonical ensemble Monte Carlo simulation is used. The hard-sphere free energy is treated within a weighted density approach and the residual ionic contribution is estimated through perturbation around the uniform density. The results from both methods corroborate each other quantitatively over a wide range of physical parameters. The importance of structural correlations is envisaged through the size and charge asymmetry of the supporting electrolytes that includes the solvent as a component.


I. Introduction
Modern day scientic and technological developments can be appreciated 1 not only through the improved quality of life for mankind, but also through the miniaturization and productivity of devices, where the focus is mainly on energy, the environment and health and safety. 2 The core research on such developments is based on the innovation of design, process development, delivery, and deployment. 3 The process should be efficient, green or environmentally friendly, and intelligently controlled -requiring minimum human intervention, thereby giving rise to long term stability and exibility in applications. Recent discoveries on a number of health related devices including nano-kits for rst aid, nano-surgical equipment, and multiple user ventilators are a few of these developments that are helping humans continuously. 4 Even personal protective equipment (PPE) and gear that requires specic protection against viruses and bacteria is also dependent on innovations in design and quality that dene their robustness and usability. 5 The organization of present day materials and their application versatility in the fabrication of devices requires thorough knowledge of bulk and interfacial states of matter to obtain the desired properties. 6 Although understanding bulk states of matter is a past problem, current day fabrication routes of more complex materials still require accurate description of the bulk matter, along with the interfacial states for which phase diagrams are still in the edging state. The multiscale modeling of materials and devices is a concrete step in this direction which allows a seamless transition from one scale to the other in terms of length and time. 7 Also, suitable manipulation of the interface with a number of new synthesis strategies, as applied to nanomaterials, are presently in the advanced state of application. 8 The intricate delicacies of the development 9 of nanodevices, as applied to nanotechnological applications, requires new methods of fabrication, including the materials, the medium, and their associated interactions. A delicate modulation of such systems also requires a molecular level description so that the design and behavior can be formulated, 10 even with a minute change of interaction.
The fact that interactions play a major role not only in the design but also in the fabrication of devices, 11 can be directly visualized from a number of experiments that involve ions, radicals, dipoles, and even multipoles. Alongside the direct interactions between the molecular components, the medium 12 or the associated solvent, 13 also plays a major role 14 in the modication of these interactions. 15 For example, crowding effects, 16,17 mostly observed in complex polymeric and biological macromolecular solutions, arise due to the presence of the solvent molecules. Similarly, solvents also act as a stabilizing medium for the ground or excited states for a number of molecules giving rise to different spectrochemical properties known as solvatochromic shis. 18 It also manifests the caging effects, 19 for ions, radicals, and polar molecules, giving rise to a slowing down in solvation time correlation functions. In order to have a direct visualization of the qualitative and quantitative aspects of charge and size correlations arising due to ionic interactions and that of the solvent, 20 in the current work, an attempt is made to understand both the effects in the same system, viz. in colloidal solutions where the supporting electrolyte and solvent is present as individual components.
Although, work including a molecular description of the solvent in an electric double layer (EDL) 21,22 in planar, 23 cylindrical, 24 and spherical geometry, 25 has been reported previously, these are mainly restricted within the restricted primitive model (RPM), 21 mainly because of its computational simplicity. However, these studies could reveal important ndings related to the volume exclusion introduced through the solvent as the competition between energetic and entropic effects gives rise to overcharging (OC) 26 and charge reversal (CR) 26 phenomena in such EDLs. Attempts to include the size asymmetry of the component ions as well as the solvent within the primitive model (PM) have been studied quite extensively in recent times. 25,27,28 However, they are also mostly conned within the system of binary electrolyte that includes the solvent 29 or a mixed electrolyte system without solvent. 30 Although, effects of size asymmetry and charge asymmetry on EDLs are understood in a scattered manner, a concerted effort to include both asymmetries has so far not been attempted.
The structure of colloidal solutions and the surrounding small ions has been improved from time to time from the classical Poisson-Boltzmann (PB) description of point ions, 21 to the modied Gouy-Chapman theory (MGC) for RPM and unequal radius MGC for PM, to include correlations arising due to different contact distances at the Stern layers. 31 A number of studies over the years have been able to conrm that along with charges on the ions and the interface constituting the EDL, size correlations of ions as well as the solvent components also contribute quite substantially in deciding its static structure. These include liquid state analytical theories, viz. modied PB theory, 32 integral equation theory, 33 and density functional theory, 34 and also simulation methods. [35][36][37] In recent times, a number of delicate experiments also point to the detailed structural behavior of ionic clouds of double layers around colloidal macroions. These include atomic force microscopy (AFM), 38 optical tweezer experiments, 39 X-ray photoelectron spectroscopy, 40 electrokinetic measurements, 41 small angle Xray scattering (SAXS), 35 etc. Attempts to include off-center charges of small ions rather than primitive centrally located ones have also been studied in great detail in recent times. 42 In the current work, however, we restrict ourselves to the primitive models but include the solvent as an individual component in the mixed electrolyte system, to understand the effects of charge and size asymmetry within the same spherical double layer (SDL). Apart from an individual component description, the solvent is also included as a continuum with the dielectric constant of water. The theoretical formalism is presented in Section II, the results and discussion part is given in Section III, and nally we offer a few concluding remarks in Section IV.

II. Theoretical formulation
A. Solvent primitive model The system considered here consists of mixed electrolytes (NaCl/MgCl 2 ) immersed in a solvent around an isolated uniformly charged spherical colloidal macroparticle of radius R with surface charge density (Q) written as where Z M is the valence of the macroion with e as the electronic charge. The solvent primitive model (SPM) is adopted here, where the electrolyte solution consists of the small ions and the solvent, that are taken as the charged and neutral hard spheres, respectively. The schematic representation of the model system is given in Fig. 1 to be noted here that with this single dielectric constant (78.5), the effects due to image charges are completely ignored, which may not be an exact situation considering different Stern layers formed due to the unequal size of the small ions. 43 The interaction potential between any two components (a and b) of the electrolyte solution is given as with the solvent-solvent and solvent-small ion interactions automatically tuned to hard-sphere potentials. Similarly, the interaction potential between any component of the electrolyte (a) and the macroion is written as ; ( with the automatic representation of the solvent for z a ¼ 0.

B. Density functional theory
In DFT, 34 the grand potential, U or the free energy, F of the system is a functional of nonuniform density {r a (r)} and they are related to each other through a Legendre transform given as where {r a } is the singlet density distribution and m a is the chemical potential, for component a. This grand potential, U should be the minimum at equilibrium that gives rise to the nal density prole of any component (a) in the present SDL written as where r 0 a is the uniform counterpart, b 0 ¼ (k B T) À1 , is the inverse temperature with k B being the Boltzmann constant and T the temperature. In eqn (5) is the rst-order correlation function, and j(r) represents the mean electrostatic potential (MEP) arising due to all the ionic charges, and can be written as 45 The expression for the density prole [eqn (5)], although exact, cannot be evaluated directly due to the unavailability of exact quantities of nonuniform c (1) a (r;[{r a ]). Hence a number of approximations are attempted from time to time for the calculation of these quantities. Here, a partially perturbative procedure is adopted and the hard-sphere and electrical components are approximated through different routes. Here, the hard sphere part, c (1) hs a (r;[{r a }]) is evaluated through the Denton and Ashcro (DA) 46 weighted density approach (WDA) given by where, the DA prescription is used to calculate the weighted density, r a (r) representing the uniform uid counterpart.
The direct correlation functions,c (2)hs ab andc (2)el ab required to calculate c (1) hs a and c (1) hs a , [eqn (7) and (8)] are taken from the analytical expressions as given by Blum 48 and Hiroike 49 within the mean spherical approximation (MSA) of the electrolyte solution for size and charge asymmetric mixtures of neutral and charged hard spheres. The calculation of density and the MEP proles are now quite obvious from eqn (5) and (6) using an iterative solution method.

C. Monte Carlo simulations
In the current work, the model system is simulated through canonical Monte Carlo (CMC) simulations (N, V, T). The initial conguration is generated by xing the macroion at the center of a cubic simulation cell and inserting the small ions and the solvent to attain the desired concentration. The cell is sufficiently large so that the interactions with any other cell can be neglected. 50 Considering the spherical symmetry of the problem, in all three perpendicular directions (X, Y, and Z), the periodic boundary conditions are employed. Once the initial conguration is generated, the system is equilibrated through random diffusion of the components in random translational moves. The long range potential between the charges is taken care of through the Ewald summation. 51 The standard Metropolis sampling procedure 52 is adopted for acceptance of the moves. A block averaging procedure is used and the density prole is calculated by dividing the box into spherical bins, and the MEP proles from the respective component densities. The relevance of the simulation is justied with the assumptions based on the model as it will be exact for the given model. For example, the size correlations for individual components were thought to be better manifested in the current method as it resembles hard-sphere mixtures in a spherical geometry. Similarly, the charge correlations will be clearly visualized, once the size correlations are taken out. The beauty of the present simulation also lies in partitioning the contribution of hardsphere solvent in the overall size and charge correlations.

III. Results and discussion
The main strategy behind the current work is to look at the size and charge correlations that are effected in the components that constitute the electrolyte solution on the structure of the SDL, around a colloidal macroion within the SPM representation. To make things more general, an individual component representation of the solvent is considered along with the mixed electrolyte system. The singlet density proles of the small ions and the solvent in the SDL for the mixed salt, i.e., 1 : 2 : 1 (NaCl/MgCl 2 ) with 1 M 1 : 1 (NaCl) electrolyte at varying concentrations of 2 : 1 (MgCl 2 ) salt, are calculated. Unless otherwise stated, the macroion surface charge density is xed at Q ¼ 0.102 C m À2 , with the radius kept constant at R ¼ 1.5 nm. Alongside a wide variation of macroion radii, the radius of coions, the surface charge density on the macroion, and the concentration of the electrolyte and that of the solvent, are also varied to test the effect of all the parameters that inuence the properties of double layers in direct or indirect ways. This also allows direct comparison on DFT and MC predictions to widen their horizons within the SDL.
The density proles of individual components of an SDL having supporting electrolytes, 1 M NaCl with 0.5 M MgCl 2 around a macroion at Q ¼ 0.102 C m À2 of R ¼ 1.5 nm, are depicted in Fig. 2, with the variation of concentration of solvent (C S ) from (a) 15 M to (d) 30 M. As expected, the density proles of all the components show increasing oscillations in moving from lower to higher solvent concentrations. This is the packing entropic effect that dominates with an increase in the total number of molecules. Size correlation is also reected in the multiple layer formation. The separation of charges among the ionic components is manifested with the increasing thickness. Also the double layer shows sufficient increase in its width from the macroion that attends the bulk values of the component densities. At low solvent concentration (C S ), the coion density at contact is lower, however, this also started increasing with increasing solvent concentration, indicating the OC effect in moving from the inner Helmholtz planes (IHPs) that start at (R + s a /2) for a ¼ 1, 2, to the outer Helmholtz plane (OHP) starting at (R + s 3 /2). As is evident, the CR effect starts appearing at the lowest solvent concentration (15 M), and becomes more pronounced with increasing solvent concentration. A clearer picture for this emerges as we plot the MEP proles with variation in solvent concentrations. Thus, Fig. 3, depicts a faster rate of drop of MEP towards the CR phenomena at higher solvent concentrations. A critical examination on the MEP curve for solvent concentration, C S , at 30 M shows substantial oscillations indicating multiple CR effects. This also conrms the larger width of the double layer for higher solvent concentrations.
As we move on to the variation of other important parameters at constant solvent concentration of C S ¼ 27.75 M, for the current SDL, a number of interesting features started appearing in terms of OC and CR effects. Thus, Fig. 4 depicts the density distributions of all the components around the macroion for a mixture of 1 M NaCl and 0.5 M MgCl 2 with varying Q from (a) 0.102 to (d) 0.408 C m À2 . Also plotted are the URMGC ionic proles, that show that the proles are completely monotonic with lower double layer thickness as compared to DFT or MC   predictions. As is clear, the counterions are strongly attracted and the coions get repelled at the interface. With increasing Q, this is seen clearly, as even the solvent density at the interface starts decreasing along with the coions. This is because of the larger size of the counterion, that gets enhanced with the opposite charge. In fact, the OC phenomena starts disappearing when moving towards higher Q. Although the spread of the double layer is reected by an increase in the thickness during passage to higher Q, the width remains more or less constant. The density proles also point to the onset of charge inversion from low Q, which starts on increasing with higher Q. The same has become quite clear in the MEP proles depicted in Fig. 5, which showed an increased diffuse layer potential at contact. The rate of drop of MEP into the CR zone is quite faster at higher Q, although the variation becomes slower in the end, thereby pointing to the same width even at higher Q. The growth of thickness is also reected in the MEP which is due to the continuous decrease of coion density together with increased counterion density at contact.
Since the larger sized counterions always lead to a substantial increase in density at the interface on the macroion, it would be worthwhile seeing what the effect of the larger coion would be on the same system. Thus, Fig. 6 depicts a multiple increase in the contact density of Mg 2+ , when the individual density proles are plotted for a mixture of 1 M NaCl and 0.5 M MgCl 2 with varying Q from (a) À0.102 to (d) À0.408 C m À2 . This is because of its higher valence coupling with the smaller size of the coion. Although the coion density proles show gradual decrease in moving to higher Q, the contact density is still quite large because of the higher concentration of the coion (Cl À ). However, the solvent density shows a gradual decrease because its size is between that of the coion and the counterions. The separation of charges becomes quite large and manifests through the increased thickness, although the dampening of the system within a short distance leads to a contraction of the double layer in its width. The same is also reected in the MEP proles as shown in Fig. 7, where the rate of drop of potential is quite steep in passing to the larger Q, indicating passage to the CR area quite fast. It is also quite clear from the MEP proles that the CR effect started appearing from a lower Q itself.
In order to have a clearer picture of the contribution of charge correlations in the current SDL, we move on to increase the concentration of the electrolyte, whereas the concentration of the solvent still remains constant at C S ¼ 27.75 M. Thus, Fig. 8 depicts the component density proles for a mixture of NaCl and MgCl 2 where the concentration is varied from 0.01 M to 2 M by keeping the ratio of NaCl : MgCl 2 ¼ 2 : 1, around a macroion of R ¼ 1.5 nm, with Q ¼ 0.102 C m À2 . As expected, the counterion density prole shows a lot of enhancement due to the electrostatic attraction and the opposite effect becomes true for the coions. The volume exclusion becomes quite prominent with an increase in the number of layers as we pass on to higher electrolyte concentrations. Charge correlations started contributing leading to an increase in OC phenomena. The interplay between charge and size correlations also leads to an increase in solvent density at the interface. The dampening of densities and the CR effect starts from 1 M electrolyte   concentration. The same is also corroborated through the MEP proles, as depicted in Fig. 9, where it is quite clear that there is a double inversion for 2 M electrolyte concentrations. As far as the spread of SDL goes, it is quite clear that the width and thickness of the double layer continuously decreases in moving towards higher electrolyte concentrations because of effective screening of the charge on the macroion by the individual components of the electrolyte.
The other important parameter through which both the charge and size correlations can be effected is the size of the macroion, hence, it is varied from (a) R ¼ 0.5 nm to (d) R ¼ 6 nm, by keeping the other parameters the same. Thus, Fig. 10 depicts the component density proles for the SDL. In fact, this is also an attempt to study the effect of increasing the overall macroion charge on the SDL. As is clear, there is a slow decrease of coion densities and the same rate of increase of counterion densities at the interface. This is on the expected line as the large sized counterions will be accommodated more compared to small sized coions. In fact, that is the case for solvent density also, which shows a slow increase in its contact density in moving to the case of a large sized macroion. This also leads to substantial increases in layering with multiple layers. The OC effect tends to have a slow decrease and there should be a slow increase in the CR effect. The same effects are also reected in the MEP proles plotted in Fig. 11. There is a sharp drop of the MEP for the larger sized macroion due to the increased CR effect. As far as the spread of the SDL is concerned, the thickness goes on continuously increasing, although the width remains more or less constant.
As steric effects due to the solvent play a crucial role in determining the ionic distributions as well as OC and CR effects, it will be highly desirable to include the bulk water concentration (55.55 M), instead of a lower value. However, this is possible only at a lower concentrations of the electrolyte. Thus, Fig. 12 depicts the density proles of the components at C ¼ 0.1 M and C S ¼ 55.5 M at different Q, as they vary away from    the surface. Since the solvent concentration is quite high, the oscillations seems to be quite large and take a longer distance to diffuse. Both the width as well as the thickness keeps on increasing with increasing Q. The diffusive behavior of the double layers is also envisaged in the MEP proles as can be seen from Fig. 13, where the MEP drops to zero at a large distance from the macroion.
A critical measurement of the width of the double layer is quite important to calculate the spread and the associated inuence of the macroion charge and is dened through the capacitive compactness, given as 53,54 where j 0 is the MEP at the macroion surface. The capacitive compactness (s c ) data is plotted in Fig. 14  The other important quantity that can be accessed through experiment 40 is the integrated charge distribution function P(r), that represents the net charge of the SDL centred around the macroion within the sphere of radius r, and is given as In essence, P(r)Z M < 0 indicates the CR effect, whereas, P(r)Z M > 0 or |P(r)| > |Z M | corresponds to surface charge amplication. The charge amplication as dened through the difference of maximum P(r) and that of the bare macroion charge (Z M ), given as DZ M ¼ P(r) max À Z M is plotted in Fig. 15 against Z M for the SDL with 1 M NaCl and 0.5 M MgCl 2 around a macroion of R ¼ 1.5 nm, with and without the solvent. It is quite clear that the adsorbed charge amplication increased quite drastically in the presence of solvent. The results obtained from DFT can be quantitatively reproduced through simulations. It should be noted here that experiments based on surface potential measurements and the Helmholtz layer composition can be    planned on nanoporous silica or polystyrene particles in the presence of supporting electrolyte solutions that includes solvent. However, specic interpretations for the adsorbed surface charges can be ascertained by coupling molecular modeling with the experimental characterization. 55

IV. Concluding remarks
The major emphasis of the present work is to investigate in detail the charge and size correlations that arise due to the combined effect or charge and size asymmetry on the SDL formed from a mixed electrolyte solution that includes the solvent as an individual component. As a matter of fact, the electrolyte solution consists of four components of different diameters, with a multivalent coion and a neutral solvent. The solvent also makes its presence known as a dielectric medium through the interaction reduced for the bare ions. For simplicity, the solvent is taken as that of water through its dielectric constant, which, however, was kept at the concentration (C S ) at half of that of bulk water at 27.75 M due to constraints in the simulations. However, for the low electrolyte concentration (0.1 M), bulk water concentration can also be included. The DFT adopted here is a partially perturbative scheme, where the DA version 46 of the WDA is used to calculate the hard-sphere part of the free energy contribution and the ionic part through the perturbation with respect to the bulk density. The required DCF to evaluate the hard sphere and the ionic part is taken from the work of Blum 48 and Hiroike 49 on multicomponent electrolytes. For comparison, a canonical ensemble Monte Carlo simulation is also performed on the same system. A wide variation of different parameters are attempted to throw light on the overcharging and charge reversal phenomena.
The size correlations in the SDL is directly reected in the increase of solvent concentration (C S ) that indicates the increase in volume exclusion, formation of multiple numbers of layers, increase of charge separation, and the widening of the charge reversal effects. The increase of thickness of the double layer is also visible in increasing macroion surface charge density (Q), although there is hardly any change in its width. Increasing the coion size larger than the counterion leads to a quick drop of MEP to the CR area within a short distance from the interface indicating a clear increase in the charge separation. This also leads to a decrease in the width of the SDL. The importance of charge correlation is reected through the increase in concentration of the electrolyte, where the OC effect as well as the CR effect shows a gradual increase. This leads to a decrease of the width of the double layer leading to the dampening of the density proles. However, increasing the overall charge on the macroion (Z M ) by increasing the size of the macroion, leads to a slow decrease of the OC effect and a slow increase in the CR effect. In all the parametric variations, it is quite clear that the density and the MEP proles show a quantitative comparison for both the DFT and MC results. The capacitive compactness data clearly indicates that the presence of solvent reduces the compactness of the double layer studied here.
The current study concentrates on the effect of charge and size asymmetry by including solvent in the mixed electrolyte system, that is, having a multivalent cation. Although, a simplistic solvent primitive model with a reduced solvent density than that of the bulk molarity of water is considered here, a rigorous representation 56 of the solvent along with the electrolyte ions is also within the purview of the present study. This description can also be applied extensively to any type of symmetry of the charged interface. In all the results presented here, only one dielectric constant as that of water (3 ¼ 78.5) is considered here. However, due to the presence of three Helmholtz planes, the dielectric constant should not be the same everywhere. 43 Extension of the same model to include image charges will be considered in the next study. The off-centric charge on the small ions will be a simple extension of the current study. 42 The interaction between large colloidal particles can be normalized in terms of solvent, once the SDLs formed from two macroparticles can be considered. The importance of charge and size correlation in such solvent-induced interactions will be quite critical to numerically simulate biological macromolecules 57 and even fuel cells 58 and batteries. 59 The current simulation method can also be applied to predict the Stern layer structure and energetics as accessed through synchrotron X-ray reectivity (XRR) experiments 60 for alkali chloride solutions on mica surfaces. Progress of all this work related to double layers will be reported in future publications.

Conflicts of interest
There are no conicts to declare.