Selective adsorption of ions in charged slit-systems

We study the selective adsorption of various cations into a layered slit system using grand canonical Monte Carlo simulations. The slit system is formed by a series of negatively charged membranes. The electrolyte contains two kinds of cations with different sizes and valences modelled by charged hard spheres immersed in a continuum dielectric solvent. We present the results for various cases depending on the combinations of the properties of the competing cations. We concentrate on the case when the divalent cations are larger than the monovalent cations. In this case, size and charge have counterbalancing effects, which results in interesting selectivity phenomena.


Introduction
This work was motivated by two previous papers [1,2], and, also, it is a direct continuation of those works. Ion selectivity was a central theme of many of our works for ion channels [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17] which, in turn, motivated our study for the selective adsorption of various ions at highly charged interfaces [1]. In that work, we simulated -using Monte Carlo (MC) -the adsorption of cations at a charged electrode, where the size and charge of cations were varied. We considered the competition of small monovalent (Sm C + ), small divalent (Sm C 2+ ), large monovalent (Lg C + ), and large divalent (Lg C 2+ ) ions. We were interested in the question which species is adsorbed in larger quantity in the double layer (DL) in competition with other species for various electrode charges and mole fractions. This is the question that we pose in this work too, except that the ions now are adsorbed in slits instead of near an isolated wall. As a matter of fact, double layers are commonly simulated between two confining walls, where either of them can be charged and serve as an electrode. If the distance of the walls is large enough, 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, the DLs overlap, and the bulk region disappears. In this case, we talk about a slit. 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 bulk system is a virtual bath where the ions go and where they come from during the insertion/deletion steps of the Grand Canonical Monte Carlo (GCMC) simulations. In this case, however, 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 ground of the potential, therefore, is ill defined. Many studies for electrolytes confined in a slit considered only a single slit in equilibrium with a bulk in the GC ensemble [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32]. The behavior of fluids near surfaces attracted the attention of many workers [33][34][35]. Electrolytes are especially interesting due to the forming DLs. If the wall is charged, well-defined contact value theorems describe the quantitative behavior of ions at the electrode [36][37][38][39][40].
In our previous work [2], therefore, we simulated a slit system, where we allowed the existence of several slits and, more importantly, the existence of two baths on the two sides of the slit system. The slit system was composed of membranes of a fixed width carrying surface charges. The slits between the membranes then adsorbed cations in order to balance the membrane charges. Our main interest in that study [2] was the electrical conditions present inside and outside the slit system, namely, the profiles for the charge density, the electric field, and the mean electrical potential. The electrical potential difference between the inner slits and the bath had two main components. First, DLs were formed outside the slit system. The potential drop associated with these DLs formed the main component of the total drop. Additionally, there was a potential drop across the outer slits. We analyzed the effect of the number of slits, the width of the slits, membrane charge, concentration, and the presence of divalent ions. The electrolyte, however, was a pure electrolyte meaning that only one kind of cation and one kind of anion were present.

Model and method
The electrolyte is modelled by the Primitive Model (PM). 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), R i is the ionic radius, and ǫ 0 is the vacuum's permittivity. We use R A = 1 Å for the radius of the anions, while we use R Sm C = 1 Å and R Lg C = 2.1 Å for the radii of the small and large cations, respectively. The ionic charges are point charges in the centers of hard spheres.
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 (the distance of the two walls forming the membrane: L M = x R α − x L α for every α) in the slit system. The value L M = 10 Å was used throughout this work. In this work, both membrane walls carry σ 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. Ions are not allowed to enter the membranes, so the membrane surfaces behave as hard walls. 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 hand bulk region and on the right hand side of the right hand bulk region. In this study, these walls are uncharged.
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

43601-2
was taken into account by the charged sheet method proposed by Torrie and Valleau [71] and developed further by Boda et al. [72].
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 assure that the simulation cell is charge neutral in every instant of the simulation. The acceptance probability of these steps is found in our previous paper [2]. The chemical potentials needed for the GCMC simulations' input were determined by the Adaptive-GCMC method of Malasics et al. [73,74].
The dimensions of the simulation cell in the y − z dimensions is in the range of 120-150 Å. System-size checks indicated little sensitivity of the concentration profiles on the y − z dimensions of the cell. In a typical simulation, the sample contained several hundreds of millions (10 8 ) configurations.
The main output quantities of the simulations are the density profiles of various ionic species, ρ i (x), from which the ionic charge profile is obtained as follows: Poisson's equation is solved for the mean electrical potential Φ(x) and the electric field E (x) = −dΦ(x)/dx using Neumann boundary conditions. The total charge density that contains all charges in the system [including the σ M membrane charges in addition to the ionic charge density q(x)] is integrated once to get E (x) and once more to get Φ(x) with the boundary condition E (−∞) = 0. The system is always charge neutral due to the design of the GCMC insertion/deletion steps. The other condition E (∞) = 0, therefore, is automatically satisfied. The zero level of the potential can be chosen freely. The equations are found in [2].

Results and discussion
In this work, we measure 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 charge profile is also computed in terms of concentrations and is normalized by the elementary charge: q(x)/e = i z i c i (x) (the unit is M, which is mol/dm 3 ). The potential is plotted in a dimensionless form as eΦ(x)/k B T . The electric field is the derivative of this dimensionless potential, so its unit is Å −1 . We will denote it by E * (x). The temperature is T = 298.15 K. We show the results for N S = 3 slits; other numbers gave similar results. We have several variables that can be changed: valences, radii, and concentrations of cations, width of slits (L S ), and membrane charge (σ M ). We used slit widths L S = 5, 10, and 20 Å, but most results will be shown for 5 Å. We changed the membrane charge in the interval −0 We change ion concentrations in a way that the ionic strength is constant in a series of calculations. The composition of the mixture can be expressed in terms of either mole fraction or in terms of a fraction of ionic strength where it was assumed that the anion is monovalent. If species i is monovalent, If the two cations have the same valence, the two fractions are the same. If one of the species is monovalent and the other is divalent, they can be transformed to each other through We will plot our results as functions of x I i , but we must be aware of the difference from the classical mole fraction. To characterize the selectivity behavior of the slits, we show selectivity (or binding) curves, in which we plot occupancies against a quantity characterizing composition, x I 1 , in this work. We define occupancy as the average number of an ionic species in the central slit divided by the cross section of the simulation cell. Occupancy, therefore, is a surface density; its unit is Å −2 . There are other ways to characterize selectivity, however. We can plot ion exchange isotherms, where the mole fraction of one of the competing species inside the slit is plotted as a function of the mole fraction of the same species in the bulk. Jamnik

43601-4
Selective adsorption in slit-systems and Vlachy [61] showed the results for the separation factor  In agreement with our ion channel studies, the slit is selective for the smaller ion as figures 1 (c)-(d) show. This is hardly a surprise. In the other two cases, monovalent ions compete with divalent ions. The case shown in figure 1 (a) was called the "selective" case, where the small divalent has the advantages of both the larger charge and the smaller size. The other case of figure 1 (b) was called the "competitive" case, where large divalent ions compete with small monovalent ions. In this case, the smaller size is an advantage for the monovalents, while the larger charge is an advantage for the divalents. We will focus on this case hereinafter.  For σ M = −0.15 C m −2 , the situation is the opposite, we "overcharge" the outer slits. In this case, |Φ DL | > |Φ M |, which results in a peculiar behavior of the potential; the potential is deeper in the outer DL than in the inner membranes.
An interpretation of the potential profile can be given as follows. Gillespie suggested breaking the excess chemical potential into various components (hard sphere, electrical, mean field, screening) in his Density Functional theoretical study of the Ryanodine Receptor calcium channel [75,76]. This made an energetic analyzis of ion selectivity in the selectivity filter (acting as a binding site) of calcium channels possible. In a GCMC study [15], we extended this approach to the presence of inhomogeneous dielectric. It was also used to analyze selective adsorption at highly charged interfaces [1]. The mean-field term appearing in this analysis is the interaction of the ion with the mean potential, z i eΦ(x). This term provides information on how the ion interacts with the average electrical potential. The negative well of the mean electrical potential means that the divalent cation has an electrical advantage by interacting with the negative mean potential. There are, however, other terms that describe hard sphere exclusion or ionic correlations. These terms counterbalance the mean field term and make the chemical potential constant. All these terms together make it possible for the ions to be in equilibrium with the bath. The mean electrical potential is only one term of the many. Figure 1 shows the selectivity curves for a specific membrane charge, σ M = −0.15 C m −2 . Selectivity, however, depends on the charge of the slit system. Figure 3 shows the selectivity curves for three smaller membrane charges. The first panel shows the data for an uncharged slit system. The charge of the ions has no effect in this case. Therefore, the small ion will enter the slits in a larger quantity (Sm C + ) along with the counterions (A − ).

43601-6
tivity appears to be present in this case. This is a result, however, of the cancellation of two competing effects. Its small size is advantageous for the monovalent cation, while its double charge is advantageous for the divalent cation. At this membrane charge, the two effects balance each other.
Increasing (in magnitude) the membrane charge to σ M = −0.1 C m −2 increases the importance of electrostatic effects and favors the divalent ion. The selectivity curves are not linear and cross each other at larger [SmC + ] concentrations. That is, more Sm C + ions are needed to replace the Lg C 2+ ions in the slits. Figure 1 (b) shows even stronger divalent selectivity for σ M = −0.15 C m −2 .
The selectivity behavior of this electrolyte mixture, therefore, depends on the membrane charge. Membrane charge determines which effect will dominate: the small size (entropic advantage) or the large charge (electrostatic advantage) of the ion. That is why we called this case the "competitive" case.
We define the crosspoint as the [SmC + ]/I value at which the two curves cross each other, namely, at which there are equal numbers of Sm C + and Lg C 2+ ions in the slits. We can characterize selectivity with this value. If it is close to 0, the system is selective for the Sm C + ions. If it is close to 1, the system is selective for the Lg C 2+ ions. Figure 4 shows the dependence of this crosspoint on membrane charge for two different slit widths. The figure clearly indicates that the divalent ions are favored at large membrane charges. The effect is more pronounced for narrow slits. This can be plotted in a different way by fixing the bath composition and showing the results as functions of the membrane charge. This is shown in the top panel of figure 5 (b). Sm C + ions dominate at small membrane charges, while Lg C 2+ dominate at large membrane charges.
Interestingly, this effect is exactly the opposite of what we found at a highly charged electrode [1].
There, we found (see figure 6 of that paper) that monovalents were favored at very large (|σ| > 1 C m −2 ) surface charges. We did not use such large membrane charges in this study. The space in a slit is limited; at this charge there would not be enough space for the cations to balance the membrane charges inside the slits. This being not possible, they would balance the charge of the slit system from outside, from the outer DLs. That practically would be the case considered in our previous work for the isolated interface [1]. The explanation of our earlier result that small monovalent ions dominate near the highly charged surface is that the density is so large close to the surface that it favors the small ions. There is, however, enough space in the diffuse layer for the divalent ions farther from the electrode so they can contribute to screening from afar. Figure 5 plots the selectivity curves as functions of the electrode charge. Figure 5   The relative amount of the two competing cations depends on slit width; the narrower the slit is, the more selective it is for the divalent ion.
Selectivity for the divalent ion against the monovalent ion was also observed by Jamnik and Vlachy [61], who simulated the competitive ion partitioning between a charged micropore and the bulk. The radius of the pore in that study, however, was quite large (40 Å), so the DLs did not overlap. Divalent vs. monovalent selectivity is weaker in this case as confirmed by figure 5 (b).
Finally, we return to the bottom panel of figure 2 (b), where the mean potential curves are plotted for different membrane charges. This figure also defines the potential differences Φ M and Φ DL . These

43601-8
potential values are shown as functions of the membrane charge in figure 6. The different shapes of the potential curves in figure 2 have the effect on the curves in figure 6. It is seen that |Φ M | > |Φ DL | at small membrane charges (in absolute value), while the reverse is true for large membrane charges. The explanation of this behavior has been given at the discussion of figure 2 (b). The basis is that we "overcharge" the outer slits at large membrane charges. This is because the large membrane charge favors the divalent ions in the slits. The DLs outside the slit system, on the other hand, is less selective for the divalent ion because strong competition is not enforced by a limited space. The excess of divalents in the outer slits, then, "overcharges" these slits. According to charge neutrality, there is a shortage of ionic charge in the outer DLs.
Our simulations for selective adsorption of various competing cations in charged slits revealed that the basic behavior is the same as in the case of calcium channels: the narrow slits favor small ions and divalent ions. The interesting case of the "competitive" case (Sm C + vs. Lg C 2+ ) was discussed in more detail.