Review and Modification of Entropy Modeling for Steric Effects in the Poisson-Boltzmann Equation.

The classical Poisson-Boltzmann model can only work when ion concentrations are very dilute, which often does not match the experimental conditions. Researchers have been working on the modification of the model to include the steric effect of ions, which is non-negligible when the ion concentrations are not dilute. Generally the steric effect was modeled to correct the Helmholtz free energy either through its internal energy or entropy, and an overview is given here. The Bikerman model, based on adding solvent entropy to the free energy through the concept of volume exclusion, is a rather popular steric-effect model nowadays. However, ion sizes are treated as identical in the Bikerman model, making an extension of the Bikerman model to include specific ion sizes desirable. Directly replacing the ions of non-specific size by specific ones in the model seems natural and has been accepted by many researchers in this field. However, this straightforward modification does not have a free energy formula to support it. Here modifications of the Bikerman model to include specific ion sizes have been developed iteratively, and such a model is achieved with a guarantee that: (1) it can approach Boltzmann distribution at diluteness; (2) it can reach saturation limit as the reciprocal of specific ion size under extreme electrostatic conditions; (3) its entropy can be derived by mean-field lattice gas model.

1. Introduction

One of the major limitations of the Poisson-Boltzmann (pan class="Chemical">PB) and Poisson-Nernst-Planck (PNP) models is the assumption of point-like ions without considering their sizes. These models based on mean field theories work well for dilute electrolytes, but break down when the concentration is high and ions are crowded in it. A high concentration would generally cause steric repn>ulsions and additional electrostatic correlations among ions, that cannot be described by classical n>an class="Chemical">PB/PNP models [1]. For example, the concentration of counter-ions, predicted by PB, can be unrealistically high near the electrode surface, when the electrode voltage is large. Another example occurs at the selectivity filter in a potassium channel, where potassium ions are strongly attracted into this extremely narrow filter by the strong negative charges of oxygens on the backbone of the filter. Employing classical PB/PNP models would overestimate the density of potassium inside the filter and give incorrect channel current predictions. Therefore, many researchers have worked on the modification of PB/PNP to include the steric effect of ions.

Steric effect has long been approached in modeling by modifying either the internal energy or entropy in the Helmholtz free energy. Through internal energy, the steric effect has been featured as excess hard-sphere energy either by density functional theory (DFT) [2,3] or Lennard-Jones potential [4]. These energies were all formulated using non-local potentials and cause the resultant modified pan class="Chemical">PB/PNP to produce a series of complicated integro-differential equations, which are hard to compute in higher dimensions. For practical implementations, localization of hard-spn>here potential and simplifying integro-differential equations into pure differential equations has been conducted in [5,6,7,8] for DFT and [9] for Lennard-Jones potential.

Through an entropy approach, Bikerman modified the classical Boltzmann distribution by adjusting bulk and local ion concentrations via the excluded volume concept [10]. Borukhov et al. [11] rigorously derived the same formula independently by adding solvent entropy through excluded volume into the Helmholtz free energy. Although the localized hard-sphere model-based DFT [5,6,7,8] also captures this solvent entropy as one of the terms accounting for excess hard-sphere chemical potential, the Bikerman model [10,11] has been a more popular steric model due to its easiness of application and qualitatively good agreement with experiments [12,13,14,15].

In order to obtain the potential and further derive a neat modified pan class="Chemical">PB equation from free energy, Borukhov et al. [11] treated all ions as having identical size, which has been long criticized for neglecting spn>ecific ion sizes. Researchers have tried to address this shortcoming with spn>ecific ion sizes, and many of them simply extended original Bikerman model by replacing the identical ion size with spn>ecific ones without any rigorous justification. Although the resultant model has a better agreement with expn>eriments than the original Bikerman model [16,17,18], it does not have a Helmholtz free energy to support it. Here modifications of the Bikerman model to include spn>ecific ion sizes have been developed iteratively in Section 4, Section 5 and Section 6, preceded by derivations of classical n>an class="Chemical">PB in Section 2, and the original Bikerman model in Section 3. Finally, in the Discussion and Conclusions section a specific-ion-size Bikerman model is presented with a guarantee that: (1) it can approach Boltzmann distribution at diluteness; (2) it can reach the saturation limit as the reciprocal of specific ion size under extreme electrostatic conditions; (3) its entropy can be derived by a mean-field lattice gas model.

2. Classical Poisson-Boltzmann Model

Though the classical pan class="Chemical">PB model is well known, we still derive the model here for review and comparison with its modified versions discussed later. Starting by stating the Helmholtz free energy, internal energy and entropy, we have: where F is Helmholtz free energy; U is internal energy; T is temperature; S is entropy; is electric potential. p, n denote cation/n>an class="Chemical">anion concentrations, and , denote their valence, respectively. e denotes elementary charge. q denotes permanent charge. Permittivity with being the permittivity for vacuum and being the relative permittivity or dielectric constant. is some reference concentration such as bulk concentration of electrolyte. and denote the solvation energies for cations and anions, respectively. Although the traditional PB model generally does not include solvation energy in the expression, it is important when modeling some electrolyte systems involving hydration/dehydration of ions and is therefore it is explicitly included in the energy here. Based on the Born model, the solvation energies for cations and anions are:

Differentiation of F with respect to gives the Poisson equation:

By doing the differentiation of F with respect to p and n, we obtain the chemical potentials for p and n, respectively:

At equilibrium, the chemical potential is uniform everywhere and therefore the local chemical potential must be equal to its bulk value, which is usually known: with the bulk chemical potential for cations and pan class="Chemical">anions: where the subscripn>t b denotes the bulk situation. Equations (8) and (9) can be solved solve for p and n: where , , and:

From (10), as , we obtain . Likewise, as , we obtain . These unrealistic infinite concentrations for p and n are mainly because ions are treated as particles without size in the classical pan class="Chemical">PB model. This pitfall has motivated modifications of the classical n>an class="Chemical">PB/PNP model to account for the finite-size effect, or so-called steric effect, of ions. In reality, the limit of should be at most , where is the particle volume of . This can be derived by considering a volume fully occupied by cation p only, with the number of cation particles being , and then:

Likewise, the limit of is at most , where is the particle volume of .

Substituting (10) into (4), and we obtain the classical pan class="Chemical">PB equation:

For z:z electrolyte without considering solvation energy, the equation above reduces to: with

3. Bikerman Model

As stated earlier, the Bikerman model [10] has been a popular steric-effect model due to its easiness of application and qualitatively good agreement with experimental data. It modifies the free energy of the classical pan class="Chemical">PB (1)–(3) by adding a solvent entropy term. This term also partially represents the excessive energy accounting for overcrowding of ions and solvent molecules in localized hard-spn>here models based on DFT [5,6,7,8]. The free energy in the Bikerman model treats all spn>ecies of ions and solvent molecules with an identical size, and is stated as follows: where is concentration of solvent (such as n>an class="Chemical">water); v is the universal particle volume. Why are all solute and solvent particles treated as having the same size? Why are specific sizes of ions and solvent molecules not used here? This was not explained in the original model [10,11], and the justification of using identical size for all species particles will be addressed later in Section 5.

If we assume that, besides occupation of ions, the rest of space is occupied by solvent molecules (which can be taken as pan class="Chemical">water here). Then: where V is the whole volume of electrolyte; are number of cation, n>an class="Chemical">anion and solvent particles in an electrolyte with volume V, respectively. Equation (18) can then be rewritten as: which simply means the sum of volume fractions of water, cation and anion is one. Note that here we assume that, besides occupation of ions, the rest of space is occupied by water.

Substituting (19) into (17) we can obtain:

Differentiation of F with respect to again gives the Poisson equation:

By doing the derivation of F with respect to p and n, we obtain the chemical potentials for p and n, respectively:

At equilibrium, the chemical potential is uniform everywhere and therefore the local chemical potential must be equal to its bulk value, which is usually known: with:

By substituting (22), (23), (25) and (26) into (24), we can relate the local ion-to-solvent volume fraction ratios (denoted as ) to their counterparts in bulk solution in a Boltzmann manner for p and n, respectively:

Summation of (27) and (28) gives the solute-to-solvent volume fraction ratio as: and we can further obtain the solute volume fraction:

From (29) and (19), we know then volume fraction for p, n and w, respectively. which further gives p and n in terms of their bulk values , , identical particle size v, and local energy , :

Equations (27) and (28) can be re-arranged to obtain: with ionic steric potential expressed as:

The steric potential , first described in [16,17,18,19], characterizes the crowding of ions and their finite-size effect by a bulk-to-local pan class="Chemical">water fraction ratio. Larger local ion concentrations would have a larger steric potential.

Also, by letting , and , Equations (27) and (28) can be simplified as:

Therefore, by (30), we obtain:

Since due to electric neutrality in bulk conditions, therefore: where . Also: then (36) becomes:

Substituting (39) into (21), we obtain the Bikerman-pan class="Chemical">PB equation:

For z:z electrolyte without considering the solvation energy, the equation above becomes: as shown in [11] with .

Two important criteria need to be checked for all modified pan class="Chemical">PB/PNP models accounting for steric effects:

CRITERION I: When ion concentrations p and n are dilute, will they follow the classical Boltzmann distribution?

CRITERION II: As will p and n approach their saturation limits and , respectively?

For CRITERION I when p and n are dilute here, it means their volume fractions are negligible, and therefore and . Steric potential term then vanishes, and by (32) , and , which follows the Boltzmann distribution.

For CRITERION II, let us consider first, and can be derived similarly. As , , and . Therefore, , and by (30), which further means , and . Likewise, as , we can get , and .

4. The Bikerman Model with Specific Ion Sizes

The shortcoming of the Bikerman model is the usage of a universal particle size, denoted by v, for cations, pan class="Chemical">anions and solvents. Using spn>ecific ion and solvent sizes would be closer to reality. Taking n>an class="Chemical">NaCl solution as an example, the spherical diameters for Cl−, Na+ and water are , , and , and then the particle volume ratio is , in which using universal particle volume would be far from reality in the case of high ion concentrations. In appearance, it seems, and many researchers did, we can just simply modify the model to include specific ion sizes by changing and to and ; similarly, and to and for (22) to (41). With this straightforward extension, we obtain p, n as: and the specific-ion-size Bikerman-PB equation:

For z:z electrolyte without considering solvation energy, the equation above becomes:

Let us denote (42) as the specific-ion-size Bikerman model 1 (SISBM1) for convenience of notation. However, we can not find an energy functional like (15)–(17) to support this naive extension, which means chemical potentials (22) and (23) with universal particle size replaced by specific ion sizes cannot be derived. The correct specific-ion-size energy functional and chemical potentials should be derived as follows:

By (47) can be rewritten as:

Differentiation of F with respect to again gives the Poisson equation:

By doing the differentiation of F with respect to p and n, we can obtain the chemical potentials for cations and pan class="Chemical">anions, respn>ectively: where .

At equilibrium, the chemical potential is uniform everywhere and therefore the local chemical potential must be equal to its bulk value, which is usually known: with:

To solve p and n from (52), there is no closed form solution like (31) for p and n due to the nonlinearity, unless some simplified case such as is considered, which is actually reduced to the original Bikerman model with . Like (32), p and n at most can be expressed as: with the steric potential being modified from (33) to include specific ion sizes:

Let us denote (55) as the specific-ion-size Bikerman model 2 (SISBM2) for convenience. Note that a similar model was also obtained in [18,19] without a rigorous derivation.

Again, we need to check criteria I and II for this specific-ion-size model. For CRITERION I, when p and n are dilute, it again means  and . Therefore by (56), and then , and by (55), which follows a classical Boltzmann distribution. For CRITERION II, as in (50), should approach for to be finite. This can only be achieved by and (saturation). Applying the same reasoning for (51), as , should approach for to be finite. Then and (saturation). This specific-ion-size model seems correct and reasonable so far, but actually there is a pitfall. That is its entropy formula (48) cannot be derived by the traditional mean-field lattice gas model. This will be explained in the next section.

5. Mixing Entropy Derivation Based on the Mean-Field Lattice Gas Model

In this section, we would like to derive the entropy in (20) by the traditional mean-field lattice gas model. Consider the entropy for an aqueous electrolyte system: where W is the number of microstates at equilibrium which possess a maximum number of microstates. Mixing entropy in electrolyte studies macrostates through spherical particles’ (solute and solvent) occupation of identical cubic sites is based on the mean-field lattice gas model. The necessity of using identical cubic sites provides a combinatorial basis when computing the maximum number of microstates. The most probable distribution of all solute (ions) and solvent particles, reaching maximum number of microstates for each species, is that each identical cubic site generally would be at most occupied by one solute/solvent particle as depicted in Figure 1a. This is based on the concept that the size of each species’ particle is infinitesimal or finite but dilute. When the actual size for each species’ particle is considered and an aqueous electrolyte is extremely concentrated as depicted in Figure 1b, the most probable distribution above may not be available. The situation in Figure 1b will be addressed in the next section.

Figure 1

(a) Moderately concentrate situation with each solute/solvent particle only occupy one identical site. (b) Extremely concentrate situation with each identical site can be allowed to be occupied by multiple solute/solvent particles of the same species in order to increase packing efficiency in space.

The entropy based on the most probable distribution of K-species solute (ions) and solvent (treated as -th species) particles, under dilute situation, over a total of available identical sites in a system is: where is the particle number of j-species ion. is the particle number of solvent, so the entropy becomes:

Using the Stirling formula with , we can rewrite the entropy as: using the following relations: where is the concentration of j-species particle; V is the volume of system; is the volume of an identical cubic site that composes the volume of system. It is naturally requested that , where is the particle volume of j-species particle. Usually in aqueous electrolyte system, where solute and solvent particles are generally crowded.

Applying (61)–(63) to (60), the entropy density can be expressed as:

For a binary electrolyte, (64) can be expressed as: or: which, without loss of generality, can be augmented as:

Equation (67) is exactly the same as (20) with:

This means the universal particle volume in the Bikerman model is actually the volume of an identical occupation site , which is limited from below by the largest particle size among all solute and solvent particles. The original Bikerman model has long suffered criticism for assuming all ions have the same size instead of using specific ion sizes in the model. The above reasoning explains why specific ion size information is left out mainly due to the need for all cubic sites to be identical in order to support the combinatorial basis demanded by the mean-field lattice gas model. Actually, information of specific ion sizes is still carried but only implicitly as shown in (68). Researchers may prefer to use SISBM2 as illustrated in Section 4, but actually its entropy formula (48) cannot be derived by the mean-field lattice gas model described above. Note that usually solute and solvent particles are treated as spheres in modeling. If are diameters for p, n, and w, respectively and their maximum is for example, then not since the identical occupation is cubic. This is why , instead of , used in [10,11].

In CRITERION II described above, as p and n should approach their saturation limits and , respectively. Here, this would be changed to approach instead of and respectively, although approaching and sounds more physically correct. This paradoxical conclusion is from entropy rigorously derived by the traditional mean-field lattice gas model based on combinatorics requiring identical occupation sites. Can this be fixed to resume the limit approach to and and still holding the ground of combinatorics at the same time? An attempt at this is discussed in the next section.

6. Entropy Fixing for Electrolytes under Extreme Concentration Conditions

Here we hope to construct a steric pan class="Chemical">PB model with entropy able to be derived by the mean-field lattice gas model, and at the same time showing physically correct saturation limits for ions as . The mean-field lattice gas model is fixed here such that each identical cubic site is allowed to be occupied by more than one solute particle of the same spn>ecies as illustrated in Figure 1b. Although this kind of distribution is no more a most probable distribution as stated earlier, it allows more efficient packing when spn>ace is extremely limited and size among spn>ecies varies largely. Again, we consider the entropy for an aqueous electrolyte system:

Let be the particle number of species j and the number of identical sites occupied by j-species particles with . This means that in an extremely concentrated situation an identical site can be occupied by more than one particle of the same species. If an identical cubic site, on average, can allow j-species particles to occupy it, we can then relate and by , or equivalently . Again, is the particle volume of species j. is the volume of an identical cubic site with . The entropy based on the most probable distribution of all ‘grouped’ species particles over a total of available identical sites in a system is:

Note that, after all the ions (in group) are distributed, there are sites that will be filled by solvent molecules in group, so the entropy becomes:

Using the Stirling formula with , we can rewrite the entropy as:

Using the following relations: where is the concentration of species j; V is the volume of the system.

The entropy per unit volume can be expressed as:

Compared with (64), specific ion sizes can now appear explicitly in the entropy formula (76). For a binary electrolyte: or: which, without loss of generality, can be augmented as: by:

Again, at equilibrium, the chemical potential is uniform everywhere and therefore the local chemical potential must be equal to its bulk value, which is usually known:

Usually bulk solutions are dilute and chemical potentials under that condition can be formulated following (80) and (81) with :

By denoting ,

(82) forms and can solve for p and n: or: where , , . Let us denote (87), (88) as the specific-ion-size Bikerman model 3 (SISBM3) for convenience.

Again, we need to check this new model with criteria I and II. For CRITERION II, we can easily deduce from (87) , (saturation). Similarly, from (88), , , (saturation). There is no constraint like as any more, and entropy here can be derived by mean-field lattice gas model.

For CRITERION I, and will not approach a Boltzmann distribution and at diluteness unless . This violation of the Boltzmann distribution at the dilution limit is because we allow multiple ions of the same species to occupy an identical cubic site. This failure and a possible cure will be discussed in next section.

7. Discussion and Conclusions

If we wish to obtain a model for electrolytes such that: (1) it can approach a Boltzmann distribution at diluteness; (2) it can reach the saturation limit as the reciprocal of specific ion size under extreme electrostatic conditions; (3) its entropy can be derived by a mean-field lattice gas model. The only options here is SISBM3 with , since SISBM2 satisfies (1) and (2) but not (3). How can we justify for SISBM3 here? Interpreting all ion sizes as being about the same is certainly not acceptable. Remember SISBM3 is designed for extremely high ion concentrations motivated by the more efficient packing shown in Figure 1b. Actually for situations that would give rise to extremely high ion concentrations and make the steric effect not negligible, such as the Stern layer in the electric double layer (EDL) of a charged wall (discussed next) and the selectivity filter of a K channel [20], there would be ‘locally’ one species only, which is the counter-ion of the local electrostatic environment, since co-ions (and even pan class="Chemical">water) would be totally expn>elled. Taking the K channel selectivity filter as an example, its extreme narrowness and the strong negative n>an class="Chemical">oxygen charges inside it would definitely justify only one species being inside the selectivity filter, which is definitely potassium. This implies to be 1 locally in the filter, and we can justify inside the filter as well since anions would be extremely dilute there due to strong electrostatic repulsion. For the rest of the K channel where ions are at most moderately concentrated, the steric effect is much less significant, and basically the original Bikerman model would be appropriate for it. Since SISBM3 with would be a very good approximation of the original Bikerman model under mild ion concentrations, we can use SISBM3 with globally for the whole K channel then. Notice, under , SISBM3 is actually same as SISBM1, but with a rigorous derivation now. This model has been useful and proven to fit the experimental data quite well [16,17,18]. Although here we just discussed steric-effect modifications for PB, modifications for PNP can be likewise derived.

Here we compare SISBM1 and pan class="Chemical">PB by computing ion distributions in a 1D charged wall problem. Many researchers have used this physical model to investigate the surface differential capacitance of electrodes adjacent to electrolyte solutions [12,13,14,15]. Here (14) for n>an class="Chemical">PB and (44) for SISBM1 were used to calculate the ion distributions of a binary KCl electrolyte solution without considering the solvation energy and permanent charges. The associated boundary conditions are and . The bulk concentration of KCl as is set to be mM, and dielectric constant is set to 80 for the whole domain . The Debye length, featuring the order of thickness of EDL, is . The simulation result is shown in Figure 2 with Figure 2a being the distributions of for SISBM1 and PB under V and 2 V. Figure 2b is the counterpart plot of Figure 2a for . In Figure 2a, the distributions for SISBM1 and PB are very close to each other and almost indistinguishable in the graph at a weak wall voltage V. When the wall voltage increases to V, calculated by SISBM1 reaches its saturation limit mM right adjacent to the wall, but unrealistically increases beyond the saturation limit when computed by the PB model. The main effect of SISBM1 is to offer a saturation limit for counter-ions (K here) when electrostatic attraction from electrode is strong enough, while it is very close to the result of PB when the electrostatic attraction is weak. In Figure 2b, distributions calculated by SISBM1 and PB are very close to each other for both strong and weak wall voltages due to the diluteness caused by electrostatic repulsion to the co-ion (Cl here) of the electrode. Note that, corresponding to a saturation layer of adjacent to wall at V (see Figure 2a), almost vanishes at that layer as well (see Figure 2b). This implies a total exclusion of Cl over there due to the saturation of K, and justifies the locally one-species argument above. If we use the original Bikerman model (41), in which ion sizes are universal, a similar saturating phenomenon for counter-ion concentration can still be obtained. However, specific ion sizes are particularly desired when electrolyte solutions are ternary, like a mixture of KCl and NaCl solutions since K and Na have different sizes, which would saturate at different limiting concentrations. These would be otherwise indistinguishable if using the original Bikerman model.

Figure 2

(a) K distributions in charged wall problem computed by SISBM1 and PB models under different wall voltages. (b) Cl distributions with conditions same as (a).

Above we assume the rest of space after the occupation of ions is exclusively occupied by solvent particles such as pan class="Chemical">water. [9,10,11,12] have suggested that the rest of spn>ace should be occupied by solvent or void, so the spn>ecies in (58) and (70) should be interpn>reted as solvent or void. This may make more sense. Taking the selectivity filter of a K channel as an example, more and more evidences have shown the selectivity filter of a K channel is exclusively occupied by n>an class="Chemical">potassium and voids, and water is not allowed there due to the strong solvation energy barrier. [9,10,11,12] even explicitly separate water and voids as two species in their modeling. However, that means the species transport equation (Nernst-Planck equation) of water needs to be modeled explicitly when constructing a PNP type model. This water equation is generally hard to model due to its physical complexity, especially for its electrostatic behavior since water is a dipole.