Monte Carlo simulation of the electrical properties of electrolytes adsorbed in charged slit-systems

We study the adsorption of primitive model electrolytes into a layered slit system using grand canonical Monte Carlo simulations. The slit system contains a series of charged membranes. The ions are forbidden from the membranes, while they are allowed to be adsorbed into the slits between the membranes. We focus on the electrical properties of the slit system. We show concentration, charge, electric field, and electrical potential profiles. We show that the potential difference between the slit system and the bulk phase is mainly due to the double layers formed at the boundaries of the slit system, but polarization of external slits also contributes to the potential drop. We demonstrate that the electrical work necessary to bring an ion into the slit system can be studied only if we simulate the slit together with the bulk phases in one single simulation cell.

If these walls are far enough from each other, the DLs formed at the walls are independent of each other and a charge neutral bulk region is formed in the middle of the simulation cell. The reference point for the electrical potential then can be set in this bulk region. If the walls, however, are close to each other so that the two DLs overlap (slit), the bulk region disappears [24][25][26][27][28][29][30][31][32][33][34]. Slits are generally simulated in the grand canonical (GC) ensemble, where the electrolyte in the slit is in equilibrium with a virtual bulk phase represented by its temperature and the chemical potentials of the ionic species. The ground of the electrical potential cannot be set in this bulk phase, because it is not present in the simulation cell and Poisson's equation cannot be integrated over it.
The studies for electrolytes confined in a slit considered only a single slit (we will call it the "lonely slit") in equilibrium with a bulk in the GC ensemble [24][25][26][27][28][29][30][31][32][33][34]. In this work we are primarily concerned with the electrical properties of the slit. Most importantly, we want to calculate the potential difference between a slit and the bulk phase. Namely, we are concerned with the electrical work needed to bring an ion from bulk into the slit. Therefore, we simulate the slits confined between membranes with bulk electrolytes on the other sides of the membranes.

Model of the slit system
The electrolyte is modeled with the Restricted Primitive Model (RPM). In this model, the solvent is represented by its dielectric response characterized by the dielectric constant ǫ, while the ions are represented by charged hard spheres interacting through the screened Coulomb + hard sphere pair potential: where r i j is the distance of two ions, q i is the ionic charge (q i = z i e, z i being the valence and e the elementary charge), and R i is the ionic radius. In the RPM, R i is the same for every ionic species (R = 1.5 Å in this work). The ionic charges are point charges in the centers of hard spheres. It was shown that this simple model can reproduce the non-monotonic concentration dependence of the activity coefficient of electrolytes as soon as we assume that the dielectric constant is concentration dependent [50]. Since this work is a model calculation, we do not change the dielectric constant with the concentration, but we fix it at the value ǫ = 78.4. The αth membrane is confined by two hard walls at x L α and x R α , where each hard wall can carry a σ M surface charge. The interaction potential between such a charged hard wall and an ion is where |x| is the distance of the ion from the surface. There are N M membranes of width L M (this is the distance of the two walls forming the membrane: α −x L α for every α) in the slit system . In this work, both membrane walls carry a σ M surface charge. The distance of two membranes, namely, the width of the slit is L S . This distance is kept fixed during the simulation, namely, the slit system is rigid. We should keep in mind, however, that the overlapping DLs between the like-charged macromolecules play a central role in the attractive force acting between them and contribute to the cohesion of such materials [38,[51][52][53].
Ions are not allowed to enter the membranes, so the following Boltzmann factor is used to forbid the ions from the membranes: where x L α and x R α are the coordinates of the left and right walls of the αth membrane. Obviously, the number of slits is N S = N M − 1.
There are two bulk phases of widths L B on the two sides of the slit system. The simulation cell is closed by two hard walls on the left hand side of the left bulk region and on the right hand side of the right bulk region. In this study, these walls are uncharged. The geometry is illustrated in figure 1.

Grand Canonical Monte Carlo simulations
Grand Canonical Monte Carlo (GCMC) simulations have been performed for the system described above. Periodic boundary conditions have been applied in the y − z dimensions. The effect of the ions outside the central simulation cell has been taken into account by the charged sheet method proposed by Torrie and Valleau [2] and developed further by Boda et.al. [54].
In GCMC simulations of the DL, in addition to the usual particle displacement steps, we insert and delete neutral clusters of ions, e. g., ν + cations and ν − anions (ν + and ν − being the stoichiometric coefficients). This way, we make sure that the simulation cell is charge neutral in every instant of the simulation. The acceptance probability of these steps is where N + and N − are the numbers of cations and anions before insertion/deletion, V is the volume into which we insert the centres of the ions, ν = ν + + ν − , χ = 1 for insertion and χ = −1 for deletion, while ∆U is the energy change of the insertion/deletion that contains the pair energies between ions, ions and charged walls, and ions and charged sheets (it becomes infinite in the case of overlap). The quantity

Solution of Poisson's equation
The main output quantities of the simulations are the density profiles of the various ionic species, ρ i (x), from which the ionic charge profile is obtained as

23803-3
Poisson's equation is solved for the mean electrical potential Φ(x) using Neumann boundary conditions. The total charge density q tot (x) contains all charges in the system (including the σ M membrane charges in addition to the ionic charge density q(x)): where σ i surface charges are placed in positions x i . After integration we obtain where the notation i (x i < x) under the sum means that we include only those sheets in the sum whose coordinates are smaller than x. The numerical integration is performed using the rectangle rule. The electric field is defined as (4.5) Since our system is charge neutral: the electric field outside the simulation domain is zero (A = L 2 is the area of the simulation cell in the y − z plane). The first boundary condition, therefore, is from which the value C 1 = 0 is obtained for the first integration constant. After one more integration the electrical potential is obtained where we assumed that C 1 = 0. Numerical integration is performed according to the trapezoidal rule.
Since the system is charge neutral, the Neumann boundary condition on the other side (E (∞) = 0) is automatically fulfilled.
The zero level of the potential is arbitrary. Therefore, we choose the integration constant C 2 so that the potential is zero in the left hand side bulk. This is done by starting with C 2 = 0, averaging the potential over the left hand side bulk (where it is constant), and deducting this value from the potential profile.
Note that this method of solving the Poisson's equation is different from the traditional convolution integral where the boundary condition is set in a bulk region (L ∞ ). This formula is used in theoretical studies of an isolated DL, where the boundary condition is set infinitely far from the electrode (L ∞ → ∞, where the potential and the electric field are zero). This equation can be used in simulation studies if there is a bulk region in the simulation cell and if the statistics is good enough. If there is no bulk region, such as in the case of a lonely slit, the convolution integral is not practical. Furthermore, equation (4.8) provides results for the potential profiles with much better accuracy than the convolution integral. A formula that is equivalent to equation (4.8) was introduced by Kiyohara and Asaka [39] without any reference to the fact that it corresponds to the Neumann boundary conditions.

Results and discussion
In this work, we measure the distances in Å, so particle densities are measured in Å −3 . In the figures, however, we plot concentration profiles that are related to the density profiles through c i (x) = 1660.58 × ρ i (x) (the unit of concentration is mol/dm 3 abbreviated as M). The charge profile is also computed in terms of concentrations and is normalized by the elementary charge: q(x)/e = z + c + (x)+z − c − (x) (the unit is M). The potential is plotted in a dimensionless form as eΦ(x)/kT . The electric field is the derivative of this dimensionless potential, so its unit is Å −1 . We will denote it by E * (x). A typical result is seen in figure 1, where concentration, charge, electric field, and potential profiles are plotted for a 1:1 electrolyte adsorbed in a slit system with parameters indicated in the caption. The concentration profiles of the cations (counter-ions) show peaks in the slits that are necessary to balance the membrane charges. At these parameters, there is enough space in the slits to accept some coions too. The concentration profiles of both ions approach the bulk values in the two bulk regions on the two sides of the slit system. Near the external membranes double layers are formed. These double layers are apparent in the charge profile. They also balance the membrane charge but they extend into the bulk, because space is available to form diffuse layers. These double layers are responsible for the potential drops between the bulk regions and the outmost membranes of the slit system (Φ DL , see figure 1). Let us define the membrane potential, Φ M , as the potential difference between the inner membranes and the bulk. The potential drop associated with the external DLs, Φ DL , dominates the membrane potential, Φ M . The difference is due to the fact that there is less ionic charge in the outer slits than necessary to balance the membrane charge. The missing charge in the outer DLs balances the membrane charges from outside.

23803-5
The fact that the outer slits and the associated membrane charges do not cancel is proved by the electric field profile. The electric field is zero in the inner membranes because inside the slit system the adsorbed ionic charge totally balances the membrane charge. In the outer membranes, however, the electric field is nonzero indicating that there are more counter-ions in the outer DLs and less counter-ions in the outer slits. (The non-zero electric field is reflected by the slope of the potential in the outer membranes.) This polarization adds an additional term to Φ M denoted by Φ POL . Of course, the relation Φ DL +Φ POL = Φ M between these potential values holds. In this case, Φ POL is much smaller than Φ DL . An interesting result is that the behaviour of the lonely slit is quite similar to the behaviour of the innermost slit when bulk phases are present. The shape of the potential profile is similar, but the reference point of the potential profile is not well defined. At least, this cannot be related to the bulk region because the bulk region is not present in the simulation cell when the lonely slit is simulated.

The effect of the number of slits
The shape of concentration profiles also look similar in the lonely slit and in the innermost slits, at least, at the scale of figure 2. Figure 3 magnifies the peaks of the counter-ion profiles. The ionic profiles in the lonely slit show a maximum at contact position: a typical hard sphere effect at the wall. The ionic profiles in the slit system, on the other hand, show a smooth decrease at the wall. This is an electrostatic effect: the dipole field of the polarized charges outside the slit system exerts a repulsive effect on the counter-ions.
This result shows that the lonely slit is a good approximation if we want to study the structure of the density profiles in the slit, but this is insufficient if we want to get information about the electrical work that is needed to bring an ion into the slit. In that case, the slit should be simulated together with the bulk phases outside it.
The effect of the width of slits In the slits, the potential profile declines (in absolute value) to zero, but it does not reach zero because charge neutral bulk regions are not formed in the narrow slits. If the slits were wide enough, independent DLs would form near the membranes and the potential would drop to zero. Figure 4 shows this effect with an increasing slit width.
As L S increases, the average concentration in the slits decreases and the cation and anion profiles in the middle of the slits become more and more similar. The potential profile, therefore, gets closer to zero.
The potential drop in the DL, Φ DL , is the same in all cases: the interior of the slit system does not effect the outer DL. The potential profile inside the slit system, on the other hand, sensitively depends on L S . For large L S , the Φ POL potential drop vanishes, and the potential is the same in all membranes. Polarization of the outer membranes (Φ POL ) appears for narrow slits.    The effect of the membrane charge Increasing the charge of the membranes, σ M , increases the amount of ionic charge that is necessary to balance it. This increases the charge difference between the slits and the membranes, which makes the potential fluctuation in the slit system larger ( figure 5). The amount of charge in the diffuse layers outside the slit systems also increases, which increases the potential drop across these DLs (Φ DL ). The potential increases non-linearly with σ M as seen in figure 6. This phenomenon was already observed in the case of DLs [11]. The DL and POL components of the potential drop are also shown in figure 6. The relative weight of the DL component in Φ M is smaller at small membrane charges, while Φ POL is relatively small at large membrane charges. As a matter of fact, Φ POL shows a non-monotonous behaviour and decreases with increasing σ M at large membrane charges.
The effect of concentration The electrode potential of the DL layer for a given electrode charge (the capacitance of the DL) depends on the concentration of the electrolyte. The electrical properties of the slit system also depend on the concentration. Figure 7 shows The behaviour of the outer slits is also different in the two cases. There are more missing charges in the outermost slit in the 0.05 M case. This is also shown in the non-zero electrical field in the outermost membrane. Therefore, the Φ POL potential term is larger at small concentrations.

23803-8
The structure of the outer DLs is also different. As in the case of a DL at an electrode, the Φ DL potential drop is larger at small concentrations.
The effect of divalent ions The overall situation for 2:1 electrolytes is very similar to that for 1:1 electrolytes. The most notable difference occurs between membrane charges of the opposite sign but the same absolute value. Figure 8 shows results for a 2:1 electrolyte in a slit system with σ M = ±0.05 Cm −2 membrane charges. The divalent ions are more efficient in balancing the membrane charge than monovalent ions because they provide twice the charge occupying approximately the same space. Polarization of the outer slit, therefore, is negligible at a negative membrane charge when the divalent cation is the counter-ion. Since the divalent ions form a more compact diffuse layer, the potential drop Φ DL is smaller in this case.
At positive membrane charge, on the other hand, both Φ DL and Φ POL are larger (in absolute value) than in the case of negative membrane charge.

Conclusions
GCMC simulation results for PM electrolytes adsorbed in a slit system are reported. We conclude that it is necessary to study the slit together with the bulk phase with which it is in equilibrium if we want to get information on the potential difference between the slit and the bulk.
This potential difference has two main components. The double layers outside the slit system produce potential drops analogous to the electrode potential in the case of a separate DL. This component of the potential drop (Φ DL ) depends primarily on the properties of the bulk electrolyte.
If the ions have difficulty in entering the slits, the outermost slits have less charge than necessary to balance the membrane charge. The membrane charge then is neutralized from outside, from the outer DLs. This charge missing from the outer slits results in a polarization of these slits and in a potential difference. This component of the potential drop (Φ POL ) depends on both the electrolyte properties and geometrical parameters of the slit system. The Φ POL component is larger in absolute value (compared to the Φ DL component) (1) if the width of the slits is smaller, (2) if the electrode charge is smaller in absolute value (but not too large), (3) if the concentration is smaller, and (4) if the monovalent ion is the counter-ion (the membrane charge positive for 2:1 electrolytes).
We plan to study (1) competition between counter-ions of different charge and/or size (ionic selectivity -this is inspired by our ion channel studies [57]), (2) ions capable of being adsorbed into the membrane regulated by a Boltzmann-factor as in our previous membrane studies [58], and (3) ion transport using our Local Equilibrium Monte Carlo method [59], with which a chemical potential gradient can be applied across the simulation cell.