Influence of anisotropic ion shape, asymmetric valency, and electrolyte concentration on structural and thermodynamic properties of an electric double layer

Grand canonical Monte Carlo simulation results are reported for an electric double layer modelled by a planar charged hard wall, anisotropic shape cations, and spherical anions at different electrolyte concentrations and asymmetric valencies. The cations consist of two tangentially tethered hard spheres of the same diameter, $d$. One sphere is charged while the other is neutral. Spherical anions are charged hard spheres of diameter $d$. The ion valency asymmetry 1:2 and 2:1 is considered, with the ions being immersed in a solvent mimicked by a continuum dielectric medium at standard temperature. The simulations are carried out for the following electrolyte concentrations: 0.1, 1.0 and 2.0 M. Profiles of the electrode-ion, electrode-neutral sphere singlet distributions, the average orientation of dimers, and the mean electrostatic potential are calculated for a given electrode surface charge, $\sigma$, while the contact electrode potential and the differential capacitance are presented for varying electrode charge. With an increasing electrolyte concentration, the shape of differential capacitance curve changes from that with a minimum surrounded by maxima into that of a distorted single maximum. For a 2:1 electrolyte, the maximum is located at a small negative $\sigma$ value while for 1:2, at a small positive value.


Introduction
Development of the electrical double layer (EDL) theory as well as the progress in the numerical technology have resulted in proposition and examination of an increasing number of EDL models. The Gouy-Chapman theory (GC) [1,2] is based on solving the Poisson-Boltzmann equation and applies a mean field approximation to describe the electrostatic interactions. The GC theory can be applied to the simplest model of an electrolyte, which does not take into account the volume of ions. In this model, ions are represented by point electric charges, and the solvent is a dielectric continuum with the relative permittivity r . It is worth noting that the Debye-Hückel theory uses the same model of an electrolyte. The GC theory describes well low density electrolytes at small electrode charges. It breaks down at higher concentrations and charges because the excluded volume and correlation effects are disregarded.
The hard-sphere and electrostatic correlation effects are considered in modern theories of EDL such as the Modified Poisson-Boltzmann Theory (MPB) [3,4], the Mean Spherical Approximation (MSA) [5], the Hypernetted Chain Theory (HNC) [6,7] or more recently the Density Functional Theory (DFT) [8]. They have been used to describe the primitive model of an electrolyte (PM) [9] and the restricted primitive model (RPM) [10]. In the PM model, ions are represented by hard spheres of different diameters with a point electric charge located at the centre. The charged spheres are immersed in a medium with the relative permittivity r characteristic of a solvent. In the RPM model, the ion diameters are the same. For EDL composed of a planar electrode and a PM or RPM electrolyte, the modern theories predict the layering effect at high absolute values of the electrode charges [11,12] and the transition of the capacitance minimum into a maximum at small electrode charges caused by an increasing electrolyte concentration [4]. The GC theory fails to describe these effects.
Due to correlation effects that are included, the MPB and DFT theories have been successfully applied to the solvent primitive model (SPM) electrolyte. SPM is the simplest model of an electrolyte which includes the volume of solvent molecules. In SPM, solvent molecules are modelled by neutral hard spheres whose diameter is the same as [13,14] or different [15] from that of ions. Neutral spheres as well as cations and anions are immersed in the continuous dielectric medium characterised by r . The presence of neutral spheres leads to the generation of density oscillations near the electrode surface and to an increase in the differential capacitance of EDL [16].
Of the formal statistical mechanical theories, the HNC has been applied to the non-primitive model (NP) [17] in which solvent molecules are modelled by hard spheres with a point permanent electric dipole moment located at the centre. However, it leads to lowered values of the relative permittivity.
The DFT theory gives fresh possibilities. Due to the free energy term of intra-molecular interactions, the DFT theory is applicable to ions of more complicated topology than spherical. Here, we must mention the work by Sokołowski [18] who has applied the modified DFT theory to the study of orientation ordering of electrostatic neutral hard dumbbells at the planar surface. Charged dumbbells, called dimers, have aroused great interest in investigation of the EDL properties [19,20]. This advanced model of an electrolyte has been used for modelling systems of high densities like ionic liquids [21].
Torrie and Valleau [10] have introduced the computer simulation technique into the investigation of EDL. At that time computer simulation results were used to confirm the correctness of EDL theories. Now, the tremendous development of numerical technology has opened new area of application inaccessible for theory. Fedorov et al. [22] have studied models of ionic liquids made of one, two or three beads, assuming that one of the hard spheres in the chain had a positive charge. Breitsprecher et al. [23,24] have conducted molecular dynamics simulations for ions with different size and valency, represented by a coarse-grained model. Silvestre-Alcantara et al. [25] have investigated the properties of ELD containing fused dimer electrolyte. Charged hard walls and hard spheres have been replaced by molecular electrodes and soft sphere ions [26][27][28]. Ions of topology characteristic of ionic liquids are investigated. The solvent is no longer a dielectric continuum but is modelled by explicit molecular models [29]. Thus, the present models of EDL have become more realistic.
Recently, we have intensively investigated the EDL containing a dimer electrolyte. Among other effects, we have analysed the influence of concentration of 1:1 electrolyte [21] and of ion valence asymmetries [30] of charged dimers on the properties of EDL. However, we did not consider the properties of EDL for asymmetric ion valencies at different electrolyte concentrations. In particular, we expect new shapes of the differential capacitance. Thus, in this paper, we discuss the influence of anisotropic ion shape on the structural and thermodynamic properties of ELD containing asymmeteric 2:1 and 1:2 dimer electrolytes of different concentrations.

Model and methods
The electric double layer is composed of a planar electrode and an electrolyte, which is a mixture of spherical anions and anisotropic cations in the shape of dimers. The dimer consists of two tangentially tethered hard spheres, one of which has a point electric charge immersed at the centre and the other is neutral. The diameters of the two spheres of a dimer and the sphere of a monomer anion are the same. The ion valencies are asymmetric with the ions being immersed in a homogeneous medium of the relative permittivity, r . The electrode is modelled by a hard planar wall with a uniformly distributed charge of surface density, σ. The image effect is not considered which means that r of the electrode material and solvent are the same. The ion-ion and electrode-ion interactions are given by and respectively. Here, e is the magnitude of the elementary charge, and Z s is the valency of the particle of species s. Also, 0 is the vacuum permittivity, r is the separation between the centres of the two hard spheres, and x is the perpendicular distance from the electrode surface to the centre of a hard sphere. The local number density ρ s (x) of the species s at a distance x is the first average quantity obtained from our simulations. The reduced local density or the singlet distribution function g s (x) = ρ s (x)/ρ 0 s , (ρ 0 s is the corresponding bulk number density) is used to describe the structure of ELD. Also, the mean The differential capacitance C d is the property which can be compared with the experimental results.
The differential capacity is defined as In practice, C d was calculated from the interpolation polynomials method introduced by Lamperski and Zydor [31].
and the local net charge per unit area, σ Σ (x) are the properties that can be used for indicating the charge inversion (CI) and charge reversal (CR) phenomena. The CI effect takes place when the electrode charge and the charge density of the second layer and, less commonly, subsequent layers of ions next to the electrode are of the same sign. When the charge of the first layer of counter-ions overcomes the charge of an electrode, the electric field reverses its direction. The function σ Σ (x) has the sign opposite to that of σ. This effect is called the charge reversal [32,33]. Essentially at some x, the integrated charge overcompensates or overscreens the electrode charge.
Hence, the often used term charge overscreening occurs in the literature.
The second average quantity obtained from our simulations is the mean orientation function 〈cos θ〉. It depends on the distance x from the electrode surface. The function 〈cos θ〉 is the average value of cos θ, where θ is the angle between the normal to the electrode surface and the straight line joining the centres of hard spheres constituting a dimer. The origin is located at the centre of the charged sphere.
The GCMC technique was applied to calculate the local densities, ρ s (x), and orientation 〈cos θ〉 profiles. This technique is recommended for inhomogeneous systems [34] as it eliminates difficulties associated with the determination of the bulk concentration which appears when MC simulations are carried out in the canonical ensemble. The details of the GCMC technique and its implementation to our investigation have been described in the previous papers [4,35]. The procedure FLIP3 [36] was used to rotate the dimer while the long-range electrostatic interactions were estimated with the method proposed by Torrie and Valleau [10]. The ionic activity coefficients required by the GCMC technique as input were obtained from the inverse GCMC method [37].

Results
Simulations were carried out for the asymmetric ion valencies 1:2 and 2:1 at three electrolyte concentrations 0.1, 1.0, and 2.0 M in the range of the electrode surface charges varying from −1.0 to +1.0 C m −2 .
Diameters of positively charged and neutral spheres of dimers and of spherical anions are equal to d = 425 pm. The other physical parameters were temperature T = 298.15 K, and the relative permittivity of water, r = 78.5. Because of the anisotropic shape of cations, the simulations were expected to require significantly longer configurational sampling. We sampled 2 billion configurations to obtain high precision averages. Figure 1 shows the singlet distribution profiles of charged dimers and spherical anions for three elec- As a result, the internal structure of the distant cation has little influence on the double layer structure with the dimer behaving like a spherical ion (see for example, reference [30]). The surface properties of spherical ions are well known, so we do not discuss them here. The upper panel of figure 1 shows the results for the 1:2 electrolyte, while the lower one for the 2:1 case. The density profiles of charged spheres are similar to those observed for spherical counter-ions. A sharp peak corresponds to the contact distance, d /2. As expected, its height decreases with an increasing electrolyte concentration. The peak is higher and thinner for di-valent dimers as we have observed earlier [30]. The influence of a neutral sphere is hardly visible. The neutral sphere density profiles have two maxima. The first corresponds to the contact distance while the second is at x/d ≈ 1.45. The height of the second maximum relative to the height of the first one increases with a decreasing electrolyte concentration. The first maximum indicates that the large fraction of dimers take parallel orientation to the electrode surface. Assuming the perpendicular orientation of a dimer with charged sphere adjacent the electrode surface, the centre of the neutral sphere is located at x/d = 1.5. A good agreement of this prediction with the simulation results  indicates that the perpendicular orientation is also very probable and its probability increases with a decreasing electrolyte concentration. This conclusion confirms the mean orientation 〈cos θ〉 results.
The mean orientation 〈cos θ〉 profiles are shown in figure 2. The positive values of 〈cos θ〉 mean that dimers prefer the orientation with the charged sphere towards the electrode. The function 〈cos θ〉 has a nearly linear course from the contact distance to the position of the second maximum on the neutral sphere density profile. At this position, the 〈cos θ〉 function has a minimum. The minimum is negative for higher electrolyte concentrations and larger dimer valencies. At longer distances, the minimum transforms into a positive maximum. After leaving the maximum, the 〈cos θ〉 function tends to zero and the orientation randomisation takes place. Thus, the behaviour of 〈cos θ〉 suggests a generally perpendicular orientation of the dimers near the electrode with the charged head nearer to it than the neutral tail, a pattern that is consistent with our earlier studies. The contact values of 〈cos θ〉 are nearly independent of the dimer valency, but at some distance from the electrode surface they are lower for the divalent dimers. However, the course of curves remains similar. In the vicinity of the electrode surface, the rotation of the neutral sphere around the charged one is hindered mainly by the hard wall of the electrode. The wall effect is independent of electrolyte concentration. The value of 〈cos θ〉 varies from 0.5 at a contact distance to 0 at x/d = 1.5. The simulation contact results oscillate at around 0.5. They are greater than 0.5 for c = 0.1 M and lower than 0.5 for c = 2.0 M. It means that, as we have stated earlier, the low concentration electrolytes support the perpendicular orientation, while the electrolytes of high concentration enhance the parallel one. In this range of distances, the rotation of the neutral sphere around the charged one is hindered by steric and electrostatic interactions with neighbouring ions, only. The second maximum of 〈cos θ〉 shows that dimers form the second layer have their charged spheres oriented towards the electrode. It is worth noting that this layer is not visible at the charged sphere density distributions for c = 1.  of the CI and CR phenomena. The increase in electrolyte concentration or the presence of a higher valency counterion intensifies these features. Indeed, at the same electrolyte concentration, the effect is more pronounced for the 2:1 system than for the 1:2 system. The CI effect is observed closer to the electrode surface than to the CR one. The mechanism of CI and CR is not completely clear [33].
The mean electrostatic potential profiles calculated from equation (3) are presented in figure 4. For  The remaining curves have a positive maximum characteristic of divalent electrolytes [38]. This maximum is related to the the CI effect. With an increasing electrolyte concentration and dimer valency, the absolute value of the mean potential and the width of EDL decrease.

13804-6
The dependence of the mean electrostatic potential of the electrode, ψ(0), on the electrode charge density σ calculated for the investigated electrolyte concentrations and on ion valencies is shown in figure 5. The curves intersect in the vicinity of σ = 0. The slope of the curves in the range of adsorption of divalent ions is smaller than in the range of adsorption of monovalent ones. The ψ(0) results are used for the calculation of the differential capacitance of EDL. Figure 6 shows the dependence of the differential capacitance, C d , of EDL on the electrode charge density calculated for different asymmetric valencies and electrolyte concentrations. For c = 0.1 M, the C d curves have a minimum corresponding to σ ≈ 0 and surrounded by two asymmetric maxima. An increase in the electrolyte concentration to c = 2.0 M transforms the minimum into a distorted maximum. For the 1:2 electrolyte, the maximum is located at a positive σ, while for the 2:1 electrolyte, at a negative σ value.
The Gouy-Chapman theory [1,2] predicts that the differential capacitance curve of EDL has the parabolalike shape with a minimum at σ = 0. Transition of the capacitance curve from that having a minimum to that having a maximum with an increasing electrolyte concentration has been explained for the RPM electrolyte by an increase in the thickness of EDL due to the formation of bi-or multi-layer structures of counterions [13,39]. In the case of charged dimers, the increase in thickness of EDL is additionally caused by the neutral spheres of dimers which separate the two layers of charged spheres: adsorbed on the electrode surface and form the second layer. Finally, the increase in valency of counterions makes the EDL thinner. Explanation of the behaviour of C d − σ curve shown in figure 6 requires consideration of all these effects.

Conclusions
A charged dimer is a simple and useful model of molecules of charged surfactants and ionic liquids. Charged surfactants are composed of a charged head and a neutral, hydrophobic tail. In a charged dimer, the head is represented by a charged sphere while the tail by a neutral one. Surfactants are diluted in a solvent (water). Therefore, more advanced investigation of surfactants modelled by a charged dimer should include solvent molecules. The interaction of the charged head and the neutral tail with molecules of solvent depends strongly on the structural and polar properties of solvent molecules. Presumably they can influence the density and orientation profiles. Ionic liquids are solvent free molten crystals composed of large organic ions with the electric charge located off the centre of the ion. Charged dimers can model ionic liquids when their size is large and the concentration high. Current investigation of EDL composed of ionic liquids as an electrolyte is extremely interesting. However, the investigation of EDL formed by charged dimers by the method of molecular simulations breaks down because of high density of a system. The problem can be solved by replacing hard potentials with soft ones.
In this paper we have investigated the structural and thermodynamic properties of an electric double layer composed of a planar charged hard wall, charged dimers and spherical anions at different electrolyte concentrations and asymmetric ion valencies, using the grand canonical Monte Carlo simulation. The density profiles of neutral spheres have two maxima corresponding to parallel and perpendicular orientations of the dimer against the electrode surface. The height of the second maximum relative to the height of the first one increases with a decreasing electrolyte concentration. Thus, we have concluded that the probability of the perpendicular orientation increases with a decreasing electrolyte concentration.
The variation in the electrolyte concentration and ion valency leads to a diversity of shapes of differential capacitance curves. At low concentrations, the capacitance curves have a minimum surrounded by two asymmetric maxima. With increasing concentrations, the minimum transforms into a distorted maximum. For the 1:2 electrolyte, the maximum is located at positive σ, while for the 2:1 electrolyte, at a negative σ value.