1. Introduction
One of the major limitations of the Poisson-Boltzmann (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 repulsions and additional electrostatic correlations among ions, that cannot be described by classical 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 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-sphere 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 PB equation from free energy, Borukhov et al. [
11] treated all ions as having identical size, which has been long criticized for neglecting specific ion sizes. Researchers have tried to address this shortcoming with specific ion sizes, and many of them simply extended original Bikerman model by replacing the identical ion size with specific ones without any rigorous justification. Although the resultant model has a better agreement with experiments 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 specific ion sizes have been developed iteratively in
Section 4,
Section 5 and
Section 6, preceded by derivations of classical 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 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/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 anions:
where the subscript
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 PB model. This pitfall has motivated modifications of the classical 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 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 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-sphere models based on DFT [
5,
6,
7,
8]. The free energy in the Bikerman model treats all species of ions and solvent molecules with an identical size, and is stated as follows:
where
is concentration of solvent (such as 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 water here). Then:
where
V is the whole volume of electrolyte;
are number of cation, 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 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-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 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, anions and solvents. Using specific ion and solvent sizes would be closer to reality. Taking 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 anions, respectively:
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.
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 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 species 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 space is extremely limited and size among species 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 water) would be totally expelled. Taking the K channel selectivity filter as an example, its extreme narrowness and the strong negative 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 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 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.
Above we assume the rest of space after the occupation of ions is exclusively occupied by solvent particles such as water. [
9,
10,
11,
12] have suggested that the rest of space should be occupied by solvent or void, so the
species in (58) and (70) should be interpreted 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 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.