Bjerrum Pairs in Ionic Solutions: a Poisson-Boltzmann Approach

TL;DR

This paper models the formation of Bjerrum pairs in ionic solutions using the nonlinear Poisson-Boltzmann equation, revealing their impact on dielectric properties, screening length, and surface interactions at higher concentrations.

Contribution

It introduces a Poisson-Boltzmann framework to analyze Bjerrum pair formation and its effects on electrolyte properties, providing new insights into ion association phenomena.

Findings

Bjerrum pairs reduce dielectric decrement with increasing ionic concentration.

Pairs influence the Debye screening length and surface-force interactions.

They decrease bulk ionic concentration and enhance attraction between charged surfaces.

Abstract

Ionic solutions are often regarded as fully dissociated ions dispersed in a polar solvent. While this picture holds for dilute solutions, at higher ionic concentrations, oppositely charged ions can associate into dimers, referred to as Bjerrum pairs. We consider the formation of such pairs within the nonlinear Poisson-Boltzmann framework, and investigate their effects on bulk and interfacial properties of electrolytes. Our findings show that pairs can reduce the magnitude of the dielectric decrement of ionic solutions as the ionic concentration increases. We describe the effect of pairs on the Debye screening length, and relate our results to recent surface-force experiments. Furthermore, we show that Bjerrum pairs reduce the ionic concentration in bulk electrolyte and at the proximity of charged surfaces, while they enhance the attraction between oppositely charged surfaces.

Equations

$${\lambda }_{\mathrm{D}}=\frac{1}{\sqrt{8\pi l_{\mathrm{B}}{n}_b}},$$

$$Z=1+{\Lambda }_{+}+{\Lambda }_{-}+{\Lambda }_{+}{\Lambda }_{-}{\mathrm{e}}^{\beta J},$$

$$a^3\left(n_w+2n_s+n_p\right)=1.$$

$$\varphi =\left(n_s+n_p\right)a^3=\frac{1}{1+2\Lambda +{\Lambda }^2{\mathrm{e}}^{\beta J}}\left(\Lambda +{\Lambda }^2{\mathrm{e}}^{\beta J}\right),$$

$$\Lambda =\frac{-\left(1-2\varphi \right)+\sqrt{{\left(1-2\varphi \right)}^2+4{\mathrm{e}}^{\beta J}\varphi \left(1-\varphi \right)}}{2{\mathrm{e}}^{\beta J}\left(1-\varphi \right)}.$$

$$n_s{a}^3\simeq \Lambda \simeq \frac{-1+\sqrt{1+4\varphi {\mathrm{e}}^{\beta J}}}{2{\mathrm{e}}^{\beta J}}.$$

$$-{\epsilon }_0{\nabla }^2\psi ={\rho }_f+{\rho }_s+{\rho }_p,$$

$${\epsilon }_0{\psi }^{\prime \prime }(z)$$

$$-p{n}_p\frac{\mathrm{d}}{\mathrm{d}z}\frac{\mathcal{G}\left[\beta p{\psi }^{\prime }(z)\right]}{\mathcal{D}(z)}$$

$$-p_w\left(a^{-3}-2n_s-n_p\right)\frac{\mathrm{d}}{\mathrm{d}z}\frac{\mathcal{G}\left[\beta p_w{\psi }^{\prime }(z)\right]}{\mathcal{D}(z)},$$

$$\mathcal{D}(z)$$

$$+\left(a^{-3}-2n_s-n_p\right)a^{3}g\left(\beta p_w{\psi }^{\prime }\right),$$

$${\epsilon }_{\mathrm{eff}}{\psi }^{\prime \prime }=2\beta e^2{n}_s\psi ,$$

$${\epsilon }_{\mathrm{eff}}(\varphi )$$

$${\color[rgb]{0,0,0}\lambda_{{\rm eff}}(\phi)=\frac{1}{\kappa_{{\rm eff}}(\phi)}=\sqrt{\frac{\varepsilon_{{\rm eff}}(\phi)}{2\beta e^{2}n_{s}(\phi)}}.}$$

$$\frac{1}{{\kappa }_{\mathrm{eff}}^2}\frac{\partial {\kappa }_{\mathrm{eff}}^2}{\partial \varphi }=\frac{1}{n_s}\frac{\partial n_s}{\partial \varphi }-\frac{1}{{\epsilon }_{\mathrm{eff}}}\frac{\partial {\epsilon }_{\mathrm{eff}}}{\partial \varphi }.$$

$$-{\epsilon }_0{\psi }^{\prime }(0)$$

$$+p_w\left(a^{-3}-2n_s-n_p\right)\frac{\mathcal{G}\left[\beta p_w{\psi }^{\prime }(0)\right]}{\mathcal{D}(0)}.$$

$$n_{\pm }(z)=n_s\frac{{\mathrm{e}}^{\mp \beta e\psi (z)}}{\mathcal{D}(z)}.$$

$$\Pi$$

$$+2p_w\left(a^{-3}-2n_s-n_p\right)\frac{\mathcal{G}(\beta p_w{\psi }^{\prime })}{\mathcal{D}}]+k_{\mathrm{B}}T{a}^{-3}\ln \mathcal{D}.$$

$$\hat{\rho }(\mathbf{r})$$

$$-\left[s_j^{+}{s}_j^{-}p+\left(1-s_j^{+}\right)\left(1-s_j^{-}\right)p_w\right]{\hat{\mathbf{n}}}_j\cdot \nabla \delta \left(\mathbf{r}-{\mathbf{r}}_j\right)],$$

$$Z$$

$$\times \int \frac{\mathrm{d}{\Omega }_j}{4\pi }\exp \left[-\frac{\beta }{2}\int \mathrm{d}\mathbf{r}\mathrm{d}{\mathbf{r}}^{\prime }\hat{\rho }\left(\mathbf{r}\right)v_c\left(|\mathbf{r}-{\mathbf{r}}^{\prime }|\right)\hat{\rho }\left({\mathbf{r}}^{\prime }\right)\right],$$

$$1$$

$$=\int \mathcal{D}\rho \mathcal{D}\phi \exp \left(i\beta \int \mathrm{d}\mathbf{r}\phi \left(\mathbf{r}\right)\left[\rho \left(\mathbf{r}\right)-\hat{\rho }\left(\mathbf{r}\right)\right]\right).$$

$$Z$$

$$S\left[\phi \right]$$

$$-\frac{k_{B}T}{a^3}\ln (\mathrm{sinc}\left(\beta p_w\left|\nabla \phi (\mathbf{r})\right|\right)$$

$$+{\Lambda }_{+}{\Lambda }_{-}{\mathrm{e}}^{\beta J}\mathrm{sinc}\left(\beta p\left|\nabla \phi (\mathbf{r})\right|\right)$$

$$+{\Lambda }_{+}{\mathrm{e}}^{-i\beta e\phi (\mathbf{r})}+{\Lambda }_{-}{\mathrm{e}}^{i\beta e\phi (\mathbf{r})})].$$

$$F[\psi ]=S[-i\psi ].$$

Full text

I Introduction

Ionic solutions are ubiquitous in electrochemical, colloidal and biological systems. The solution properties are determined by the interplay between the ion mixing entropy and their electrostatic interaction Israelachvily ; VO ; David95 . An important length emerging from this interplay is the Bjerrum length, $l_{{\rm B}}=e^{2}/\left(4\pi\varepsilon k_{\mathrm{B}}T\right)$, where $e$ is the electronic unit charge, $\varepsilon$ the dielectric constant of the solution, and $k_{\mathrm{B}}T$ the thermal energy. At this length, the Coulombic interaction between two unit charges is equal to the thermal energy. For water, $\varepsilon\approx 78\,\varepsilon_{0}$, where $\varepsilon_{0}$ is the vacuum permittivity, and the Bjerrum length is equal to about 0.7 nm, at room temperature.

When the Bjerrum length is comparable with the lattice spacing, it is favorable for a salt crystal (e.g., NaCl) to dissociate, forming an ionic aqueous solution. In such solutions, another important length scale naturally emerges. This is the Debye screening length, $\lambda_{{\rm D}}$, that was introduced in the 1920’s DH by Debye and Hückel (DH). For fully dissociated monovalent salt with bulk concentration $n_{b}$, the Debye screening length is given by

[TABLE]

and at distances larger than $\lambda_{{\rm D}}$, the electric field induced by an ion is exponentially screened.

The apparent dichotomy between salt crystal and fully dissociated free ions in solution is an over-simplification, especially for concentrated solutions, where complexation of ions has been proposed as an alternative structure already in the early works of Bjerrum Bjerrum26 . Bjerrum postulated that oppositely charged ions in solution can associate into ionic pairs (dipoles), which are nowadays referred to as Bjerrum pairs. These pairs were shown, for example, to play an important role on the critical behavior of ionic solutions Fisher93 ; Levin96 . In addition, they arise naturally in two-component hardcore plasma Netz00 .

The formation of ionic pairs has a two-fold effect on the screening length, $\lambda_{{\rm D}}\sim(\varepsilon/n_{b})^{1/2}$. First, pairs reduce the concentration of free ions that participate in the screening to a value lower than $n_{b}$ Zwanikken09 . Second, pairs increase the solvent permittivity, $\varepsilon$, due to their permanent dipolar moment. Furthermore, free ions are known to lower the dielectric constant of ionic solutions Hasted48 ; Barthel95 . This phenomenon was recently attributed to correlation effects Levy12 ; Levy13 and may also be related to the exclusion of water dipoles, as is described below. In this picture, ionic pairs can increase the dielectric constant by excluding less water dipoles.

The simple remarks above are related to recent surface-force experiments Smith16 conducted on ionic liquids and solutions at relatively high ionic concentrations, (up to about 4 M). The screening length fitted from the force profiles between two surfaces was shown to be non-monotonic as function of the ionic concentration. It first decreases for low ionic concentrations, and then anomalously increases for higher values.

In this paper, we systematically incorporate the formation of Bjerrum pairs into the Poisson-Boltzmann (PB) theory. We demonstrate that pairs have a qualitative effect on bulk electrolyte properties. Our model predicts that a non-negligible fraction of salt ions associate into pairs. Considering both free ions and pairs, we are able to account for the nonlinear behavior of the dielectric decrement of ionic solutions on the mean field (MF) level. Furthermore, the screening length is found to be qualitatively modified by Bjerrum pairs, depending on the relative strength of the dipolar moment of the ion pairs, as compared to that of the solvent. Finally, we obtain the effect of Bjerrum pairs on the counterion profiles next to charged surfaces, and on the osmotic pressure between two oppositely charged surfaces.

The outline of our paper is as follows. The model is formulated in Section II, where we calculate the bulk concentration of Bjerrum pairs. In Section III, we describe how ion pairs modify the dielectric constant and screening length. Next, in Section IV, we present results for the ionic profiles at the proximity of a charged surface and elaborate on the local effects of the pair formation, and in Section V, we solve the two-plate problem and calculate the corresponding osmotic pressure. Finally, Section VI offers a discussion that includes a comparison with previous models and relevant experiments.

II Model

Consider an aqueous solution with added electrolyte. For simplicity, we assume a monovalent electrolyte with cations/anions of unit charge $\pm e$. The two ionic species have the same bulk concentration, $n_{b}$, satisfying electro-neutrality. The solvent (water) molecules are modeled as dipoles of concentration $n_{w}$ and permanent dipole moment $p_{w}$. We further assume that a fraction of the cations and anions can associate into Bjerrum pairs, modeled as dipoles with moment $p=be$, and the length $b$ corresponds to a typical separation between the paired ions. As a result of the association, the total number of ions partitions into free ions of concentration $n_{s}$, and ion pairs (dipoles) of concentration $n_{p}$, satisfying $n_{s}+n_{p}=n_{b}$. The two ionic states (free ions and pairs) are in chemical equilibrium that determines the $n_{s}$ and $n_{p}$ values for a given $n_{b}$, as is described below.

II.1 Free-ion and pair concentrations

The concentrations, $n_{s}$ and $n_{p}$, can be calculated, for example, via a lattice-gas model. We model the bulk solution as a cubic lattice with a unit cell of volume $a^{3}$. Each cell is occupied by either solvent, a cation, an anion, or a Bjerrum pair. The ion-pair association energy is $-J$, which accounts for the electrostatic attraction and, possibly, short-ranged ion-specific interactions. As the electrostatic attraction governs in such small length scales, we restrict ourselves hereafter to positive $J$ values, $J>0$. Note that the limit of zero ion pairs corresponds to the limit $J\to-\infty$. On the MF level, no additional interactions are considered in the bulk.

For simplicity, we assign the same lattice size to the solvent, cations, anions, and ion pairs, implying that all species occupy a similar volume. Note that this volume is not necessarily the bare ionic size. Water molecules form hydration shells around ions and dipoles and swell their effective volume Israelachvily . Hence, the length $a$ does not correspond to the ionic diameter but to the diameter of the hydrated ion.

Without any interaction between unit cells, the grand-canonical partition per cell is given by

[TABLE]

where $\beta=1/(k_{\mathrm{B}}T)$ is the inverse thermal energy and the fugacities, $\Lambda_{\pm}=\exp\left(\beta\mu_{\pm}\right),$ depend on the chemical potentials, $\mu_{\pm}$. For convenience, the chemical potential of the solvent (water) is set to zero. This is possible, as the number of solvent molecules, $n_{w}$, is determined from the condition,

[TABLE]

As Eq. (2) is symmetric under the exchange $\Lambda_{-}\leftrightarrow\Lambda_{+}$, the two fugacities are equal, $\Lambda_{+}=\Lambda_{-}\equiv\Lambda$. They are related to the bulk concentration by $2\phi=\Lambda\partial\ln Z/\partial\Lambda$, where $\phi=n_{b}a^{3}$ is the average number of cations/anions in a unit cell, ranging from $\phi=0$ for pure solvent to $\phi=1$, where all lattice sites are occupied by Bjerrum pairs. It then follows that

[TABLE]

where the first term on the right-hand-side is equal to $n_{s}a^{3}$ and the second one is $n_{p}a^{3}$.

Solving Eqs. (4) leads to the following expression for $\Lambda(\phi)$:

[TABLE]

From Eqs. (4) and (5), we obtain the dependence of the free-ion bulk concentration on $\phi$, $n_{s}=n_{s}(\phi)$, as is plotted in Fig. 1.

The limit of $\Lambda\ll 1$ corresponds to dilute ionic solutions ($\phi\ll 1$), where steric effects are negligible and Eq. (5) simplifies to

[TABLE]

The dipolar concentration is then given by $n_{p}=a^{3}n_{s}^{2}\exp(\beta J)$. The result of Eq. (6) can also be obtained straight-forwardly Zwanikken09 by introducing a kinetic constant for the reaction, $K=a^{-3}\exp(-\beta J)$, satisfying $K=\left[n_{s}\right]^{2}/\left[n_{p}\right]$. However, Eq. (5) is more general, and can be applied to higher ionic concentrations, where steric effects must be taken into account.

Another interesting case is the limit of small association, $\phi\exp(\beta J)\ll 1$, where the square root term in Eqs. (5) and (6) can be Taylor-expanded, yielding $n_{s}=n_{b}\left[1-\phi\exp(\beta J)\right]\simeq n_{b}$. As expected, almost all ions are free in this case, due to a combination of low ionic strength and small association energy. In the opposite limit of large association, $\phi\exp(\beta J)\gg 1$, Eq. (5) leads to $n_{p}a^{3}=\phi\left(1-\phi\right)$ and $n_{s}a^{3}=\exp\left(-\beta J/2\right)\sqrt{\phi\left(1-\phi\right)}$, i.e., the free-ion concentration is damped by a factor of $\exp\left(-\beta J/2\right)$.

The free-ion concentration, $n_{s}$, with and without steric effects (Eqs. (5) and (6), respectively), is presented in Fig. 1. Both expressions lead to a smaller amount of free ions ($n_{s}<n_{b}$), as compared to the standard theory without pairs ($n_{s}{=}n_{b}$). We also note that pair association is more pronounced for larger $J$ values, leading to a further reduction in free ions. Steric effects become evident for non-negligible concentrations ($\phi\gtrsim 0.1$), and favor pair formation that decreases the ionic volume fraction. As $\phi$ further increases, Eq. (6) breaks down, and only the full expression of Eq. (4) remains valid. As can be seen from Fig. 1a, the free-ion concentration increases with $n_{b}$ and reaches a maximal value always for $\phi=0.5$, and independent of $J$. This can be explained because of the $\phi\leftrightarrow\left(1-\phi\right)$ symmetry of the lattice-gas model.

II.2 Modified dipolar Poisson-Boltzmann (MDPB) equation

So far we discussed the properties of bulk ionic solutions where the electrostatic potential, $\psi$, and the ionic concentrations are fixed throughout the solution. However, a charged surface induces a spatially varying electrostatic potential $\psi({\bf r})$ and ionic concentrations. Such a potential solves Poisson’s equation,

[TABLE]

where $\rho=\rho_{f}+\rho_{s}+\rho_{p}$, is the total charge density. The first term, $\rho_{f}$, is the fixed charge density of any macromolecules and surfaces, $\rho_{s}$ is that of salt ions, and $\rho_{p}$ is related to the polarization field of the solvent and Bjerrum pairs. Note that the vacuum permittivity $\varepsilon_{0}$ is used, because the contribution of the solvent dipolar molecules is taken into account explicitly via $\rho_{p}$. In thermal equilibrium, the densities $\rho_{s}$ and $\rho_{p}$ are further related to the electrostatic potential by their corresponding Boltzmann factor (not shown here for simplicity). Combining these terms in Poisson’s equation, Eq. (7), leads to the modified dipolar PB (MDPB) equation Borukhov00 ; Abrashkin07 ; Iglic11 ; Levy12 ; Levy13 . This equation takes into account steric effects by relating $a$, the lattice constant of Section II.1, with a reference close-packing density, $a^{-3}$.

Consider a planar surface with a fixed charge distribution, $\rho_{f}$, homogeneous in the $xy$-plane. The corresponding MDPB equation depends only on the distance $z$ from the surface and reads (see Appendix):

[TABLE]

where $\psi^{\prime}=\mathrm{d}\psi/\mathrm{d}z$ and the first two terms on the right-hand-side are $\rho_{f}$ and $\rho_{s}$, while the last two are the contributions of dipolar Bjerrum pairs and solvent molecules to $\rho_{p}$, respectively. In Eq. (II.2), $\mathcal{G}(u)=\cosh u/u-\sinh u/u^{2}$ is related to the Langevin function $\mathcal{L}(u)=\coth u-1/u$ by $\mathcal{G}=\left(\sinh u/u\right)\mathcal{L}$. This is the polarization density, written as a product of the dipole density, $\sinh u/u$, and the average dipole moment, given on the MF level by the Langevin function, $\mathcal{L}$.

The denominator in Eq. (II.2), $\mathcal{D}$, is given by

[TABLE]

where the function $g(u)=\sinh u/u$ satisfies $g^{\prime}(u)=\mathcal{G}(u)$. The denominator, $\mathcal{D}$, leads to a saturation of the ionic and dipolar concentrations, bounded from above by the close-packing density, $a^{-3}$. This saturation feature is especially important near highly charged surfaces that attract counterions and dipoles, as is demonstrated in Section IV.

Equation (II.2), together with the relations between $n_{s}$ $n_{p}$, and $n_{b}$, define our model for Bjerrum pairs in ionic solutions. This model can be obtained formally by employing a field-theoretical approach to a lattice-gas model including electrostatic interactions, as is further explained in the Appendix A. Within this framework, Eq. (II.2) is the Euler-Lagrange equation, derived from the variation of the free energy. Note that in the formulation of this section as well as in the Appendix, we account only implicitly for some of the Coulombic interactions in the system via the association energy, $J$, and the hydrated ionic diameter, $a$ (see Section II.1).

Equation (II.2) can be simplified under several conditions. For small electrostatic potentials and fields, it is possible to linearize the equation, as is further discussed in Section III. This is especially relevant in the bulk electrolyte, far from any fixed charge density.

In addition, as opposed to the picture of Section II.1, where solvent molecules are considered explicitly, it is possible to treat the solvent implicitly as a homogeneous background with a dielectric constant $\varepsilon$ of the solution, by inserting $p_{w}=0$ and substituting $\varepsilon_{0}\to\varepsilon$ in Eq. (II.2). Hereafter, we distinguish between the explicit solvent model of Eq. (II.2) and the implicit solvent model, described here. Throughout this work, we present results for both models and compare between the two.

Finally, for small volume fractions of ions and Bjerrum pairs ($\phi\ll 1$), it is possible to neglect steric effects. Then, the denominator function, Eq. (II.2), simplifies to unity, $\mathcal{D}\simeq 1$ Abrashkin07 ; Levy12 ; Levy13 , and the bulk concentrations of free ions and Bjerrum pairs are given by Eq. (6) Zwanikken09 .

III Dielectric constant and screening length

The association of Bjerrum pairs modifies the permittivity of the solution and decreases the number of free ions responsible for screening. As both the dielectric constant and screening length characterize the bulk electrolyte, the above effects can be captured by linearizing Eq. (II.2), leading to the general form

[TABLE]

where $\varepsilon_{{\rm eff}}$ is the effective dielectric constant, given by

[TABLE]

This expression is typical of effective medium theory, where the contribution of each species is proportional to its volume fraction in the solution.

For pure solvent, Eq. (11) yields $\varepsilon=\varepsilon_{0}+\beta a^{-3}p_{w}^{2}/3$. This is a known MF result for a dilute phase of dipoles and is not expected to yield the correct dielectric constant of pure water. For example, substituting $a=0.5\,$nm, the value $\varepsilon=78$ corresponds to $p_{w}\approx 9.8\,{\rm D}$, which is about five times larger than the physical value $p_{w}=1.85\,{\rm D}$ of water molecules. Therefore, in addition to $a$, the dipole moment $p_{w}$ should be regarded as a parameter of the model. A more physical relation between $p_{w}$ and $\varepsilon_{{\rm eff}}$ can be obtained by incorporating electrostatic correlations and fluctuations beyond MF Levy12 ; Levy13 .

In the absence of pairs ($n_{p}=0$), Eq. (11) describes a dielectric decrement with increased ionic concentration that is known for aqueous solutions Hasted48 ; Barthel95 . In Eq. (11), the decrement originates from the condition of Eq. (3). An increased ionic concentration results in a decreased solvent concentration, reducing the effective permittivity. Moreover, ions in an aqueous solution are surrounded by an hydration shell of water molecules that are not free to rearrange themselves in response to an external field, thus lowering the dielectric constant Dan11 .

Ion-pair association makes the above picture more complex. Bjerrum pairs have a permanent dipole moment that increases the dielectric constant. Each pair adds to the dielectric constant a term proportional to $p^{2}-p_{w}^{2}$, and the behavior of $\varepsilon_{{\rm eff}}(\phi)$, therefore, changes qualitatively according to the sign of the difference $p-p_{w}$. This statement is further discussed below.

The effective screening length, $\lambda_{{\rm eff}}=1/\kappa_{{\rm eff}}$, is obtained from Eq. (10) by converting it to the form $\psi^{\prime\prime}=\kappa_{{\rm eff}}^{2}\psi$. This leads to

[TABLE]

While it is similar to the Debye screening length of Eq. (1), the above Eq. (12) suggests a more intricate dependence on the concentration $\phi$ as both $n_{s}=n_{s}(\phi)$ and $\varepsilon_{{\rm eff}}=\varepsilon_{{\rm eff}}(\phi)$. In particular, the screening length does not necessarily decrease monotonically with $\phi$ Zwanikken09 , as opposed to the classical DH result.

In order to explore this possibility, we examine the derivative of $\kappa_{{\rm eff}}^{2}$ with respect to $\phi$

[TABLE]

The first term on the right-hand-side of Eq. (13) describes whether the free-ion concentration increases or reduces with $\phi$. According to Sec. II.1, it is positive for $\phi<1/2$ and negative for $\phi>1/2$. The second term, on the other hand, depends on whether ions induce a dielectric decrement or increment. As described above, this greatly depends on how large is the ion-pair dipole in comparison to the solvent dipole.

Below, we distinguish between the two cases, $p<p_{w}$ and $p>p_{w}$, and examine the behavior of $\varepsilon_{{\rm eff}}(\phi)$ and $\kappa_{{\rm eff}}(\phi)$ for both cases. For aqueous solutions, $p<p_{w}$ due to the high dielectric constant of water. In order to examine the $p>p_{w}$ scenario, we broaden the scope of our discussion to general solvents rather than focusing exclusively on water.

III.1 The $p_{w}>p$ case

For $p_{w}>p$, ion pairs decrease the dielectric constant and Eq. (11) captures the dielectric decrement described after Eq. (11) for aqueous solutions. This is illustrated in Fig. 2a. At low $\phi$ values, $n_{p}$ is still small, and the dielectric constant decreases linearly with $n_{s}$, similarly to substituting no pairs ($n_{p}=0$ and $n_{s}a^{3}=\phi$) in Eq. (11), which is drawn as the dashed red curve in Fig. 2a. For higher ionic concentrations, Bjerrum-pair dipoles contribute to the permittivity, resulting in a smaller decrease (black curve in Fig. 2a). The diminished decrement for high ionic concentrations is consistent with measurements and was previously accounted for theoretically by electrostatic correlations beyond MF Levy12 ; Levy13 .

The $\kappa_{{\rm eff}}(\phi)$ behavior is illustrated in Fig. 2b. As the dielectric constant decreases in our model, $\kappa_{{\rm eff}}\sim 1/\varepsilon_{{\rm eff}}$ is expected to be larger than in the classical DH theory. This is indeed the case in the absence of pairs (red curves in Fig. 2b). However, the shallow decrement in presence of Bjerrum pairs, alongside the smaller number of free ions, result in a $\kappa_{{\rm eff}}$ (black curve) that is smaller than the DH result. Furthermore, $\kappa_{{\rm eff}}$ increases monotonically, in accordance with Eq. (13). Note that $\kappa_{{\rm eff}}(\phi)$ can decrease for $\phi$ values larger than $1/2$ (not shown in Fig. 2b), where fluctuations are expected to play an important role and the MF approximation breaks down.

III.2 The $p_{w}<p$ case

We consider a solvent such as alcohol with $\varepsilon=20$, smaller than that of water. For such relatively low $\varepsilon$ values, consistent with $p_{w}<p$, we find that the dielectric constant is non-monotonic with $\phi$, as is illustrated in Fig. 3a. While at low $\phi$ values, $\varepsilon_{{\rm eff}}$ follows the linear decrement of the theory without pairs (dashed red curve), the Bjerrum-pair contribution at higher concentrations leads to a dielectric increment (black curve). Such an increment can also be seen in the implicit solvent model, where $\varepsilon_{\rm eff}=\varepsilon+\beta n_{p}p^{2}/3$, and $\varepsilon$ is the dielectric constant of the solution without the ionic pairs.

The $\kappa_{{\rm eff}}(\phi)$ behavior in the $p_{w}<p$ case is smaller than the result of classical DH theory as well as of Eq. (12) without pairs. However, the screening length in Fig. 3b is non-monotonic. We emphasize that this anomalous behavior of $\kappa_{{\rm eff}}(\varepsilon)$ originates from the dielectric increment [Eq. (13)], and is inconsistent with the dielectric decrement that was measured for aqueous solutions.

IV Ionic profiles close to a charged plate

Beside exploring bulk properties of ionic solutions, we also investigate how ion-pair formation influences the ionic concentration near a charged surface. Consider the ionic solution of Section II confined to the $z{>}0$ half-space by a charged surface at $z{=}0$. The surface is homogeneously charged with a surface-charge density, $\sigma$, i.e., $\rho_{f}(z)=\sigma\delta(z),$ where $\delta(z)$ is the Dirac delta function.

The boundary condition at $z{=}0$ is determined by Gauss’ law by integrating Eq. (II.2) over a small interval around $z{=}0$. We assume that the electric field is confined to the upper region $z{>}0$, leading to

[TABLE]

The boundary condition is more complex than in the regular PB model and includes the free surface-charge density, $\sigma$, as well as the polarization-induced bound surface-charge density, which lowers (in absolute value) the total surface charge .

Once Eq. (II.2) is solved and the electrostatic potential is found, the ionic profiles $n_{\pm}(z)$ as function of the distance from the surface, $z$, are evaluated from the corresponding Boltzmann factor,

[TABLE]

As the boundary condition of Eq. (IV) is too complex to be solved analytically, we present in Fig. 4 numerical results for the counterion profile, $n_{+}(z)$. For simplicity, we treat the solvent implicitly as a homogeneous dielectric background with a dielectric constant $\varepsilon$, as opposed to dipoles in vacuum. This is equivalent to substituting $p_{w}\to 0$ and $\varepsilon_{0}\to\varepsilon$ in Eqs. (11) and (IV), as well as in Eqs. (II.2) and (II.2) in Section II.2, and is applicable at low $\phi$ values.

Counterion profiles of positive ions next to a negatively charged surface are illustrated in Fig. 4 for two $J$ values. The effect of ion association is evident from the inset, where the relative difference, as compared to the theory without ion pairs, is presented. We find that the counterion concentration is overall smaller due to the lower bulk value, $n_{s}$, that decreases as function of $J$ (Section II.1).

Furthermore, both free ions and ion pairs accumulate next to the charged surface, where the electrostatic potential and field are the strongest. Although free ions are attracted more strongly to the surface, pairs are also present at the surface proximity due to mixing entropy. As a result of this accumulation, steric effects are pronounced up to a distance $a$ from the surface, as is illustrated in Fig. 4. Far away from the surface, the counterion concentration reaches its bulk value, $n_{s}$.

V Osmotic pressure between two charged plates

We consider two charged surfaces at $z=\pm d/2$, bounding an ionic solution. The top surface has a surface-charge density $\sigma_{1}$, and the lower $\sigma_{2}$. We focus on the scenario of oppositely charged surfaces, $\sigma_{1}=-\sigma_{2}=\sigma$, where the ion-pair formation is expected to have the most evident effect on the osmotic pressure. The electric field is larger between oppositely charged surfaces, as compared to equally charged ones, leading to higher dipolar concentration and larger deviation from the standard PB theory. As in Section IV, we assume that the electric field is confined to the aqueous solution region, $-d/2<z<d/2$.

The osmotic pressure, $\Pi=p_{{\rm in}}-p_{{\rm out}}$, is defined as the difference between the pressure inside the solution and the one exerted by the bulk electrolyte reservoir. The inner pressure, constant throughout the system, is equal to the first integral Dan09 ; Safinyabook of the differential equation, Eq. (II.2). The outer pressure is obtained by setting the electrostatic potential and electric field to zero. Taking the first integral of the MDPB equation, Eq. (II.2), and subtracting the outer pressure leads to an exact expression,

[TABLE]

In the above equation, the attractive (negative) term originates from the electrostatic energy density $\sim\mathbf{E}\cdot\mathbf{D}$, with $\mathbf{E}=-\psi^{\prime}\hat{z}$ being the electric field, and $\mathbf{D}=D\hat{z}$ the displacement field. The repulsive (positive) terms, on the other hand, originate from the mixing entropy of ions and of dipoles as well as the rotational entropy of dipoles $\sim\ln\mathcal{D}$.

The electrostatic potential can be solved numerically from Eq. (II.2), with the boundary conditions that are obtained in Section IV [Eq. (IV)], and applied to the two surfaces at $z=\pm d/2$. The osmotic pressure is then found via Eq. (V). A characteristic pressure profile between oppositely charged surfaces is illustrated in Fig. 5.

The pressure in Fig. 5 is negative due to the electrostatic attraction between oppositely charged surfaces. It is enhanced by the Bjerrum pairs, as is evident from the relative difference with respect to the theory without pairs, illustrated in the inset. This difference becomes large by an order of magnitude at large separations, where the electrostatic potential and field become small and the linear approximation of Section III is valid. The linearized form leads to an exponentially decaying osmotic pressure, $\Pi\sim-\exp\left(-\kappa_{{\rm eff}}d\right)$, with a screening length that decreases due to ion association (Section III). The relative difference, therefore, is expected to increase exponentially.

We note that for sufficiently small inter-surface separations, the two surfaces effectively neutralize each other, and free ions are expected to be released to the bulk Sam05 ; Dan07 for entropic gain. In such a case, screening effects are negligible, and the above reduction in the pressure should not necessarily hold.

VI Discussion

The present work addresses the effect of Bjerrum pairs on the properties of ionic solutions in the bulk and near charged surfaces. The main result of our model is the MDPB equation in presence of the Bjerrum pairs, Eq. (II.2). This MF equation is written in terms of the free-ion and pair densities, $n_{s}$ and $n_{p}$, respectively, determined by the lattice-gas model of Section II.1. As pairs are predicted to associate for relatively high ionic concentrations, the theory accounts for steric effects and the dielectric decrement of ionic solutions, which are important under such physical conditions.

The ionic mixture considered in our model is symmetric, composed of anions and cations of the same valency and comparable size. However, asymmetric mixtures can also be studied in a similar fashion. For ions and dipoles of different sizes, it is possible to employ an asymmetric lattice-gas model (see, e.g., Ref. maggs16 ). Furthermore, for anions and cations of different valencies, one can consider the association into a charged anion-cation pair dosSantos10 .

The model includes four physical parameters: the association energy $J$, lattice size $a$, the Bjerrum dipolar moment $p$, and the water dipolar moment $p_{w}$. These parameters can be evaluated from measurements of several physical quantities, such as the dielectric constant. While we treat these parameters as independent, they are in fact coupled. Namely, the values of the $J$, $p$, and $p_{w}$ parameters depend on the lattice size, $a$. Furthermore, the quantities are coupled via the effective solvent permittivity that decreases the association energy and increases the hydration shell. Such mutual dependence can be taken into account to some extent via the Bjerrum postulate Bjerrum26 , according to which the kinetic constant $K$ of pair association, $K=\left[n_{s}\right]^{2}/\left[n_{p}\right]$, is proportional to the inverse of the integral $\int_{v}\mathrm{d}{\bf r}\exp\left(l_{{\rm B}}/r\right)$. The integration volume, $v$, is the region within which the ions are considered associated, related to the ionic diameter and Bjerrum length Bjerrum26 ; Fisher93 ; Levin96 ; Zwanikken09 .

Recent surface-force experiments Smith16 suggest a non-monotonic dependence of the screening length on the bulk ionic concentration for ionic concentrations corresponding to a range of $\phi$ values as low as $0.1$. This behavior was recently attributed to solvent molecules acting as defects in an ionic crystal Lee17 . Our prediction for the effective screening length [Eq. (12)] deviates from the usual PB result. We obtain different results, depending on the relative dipolar moment $p-p_{w}$. While for $p<p_{w}$, the screening length is monotonically decreasing, it can increase with ionic concentration for $p>p_{w}$. For aqueous solutions and simple salt, the $p>p_{w}$ scenario is not appropriate , and expected for lower dielectrics (e.g., alcohol) and/or polarizable ions. We demonstrate in Eq. (13) that in light of the dielectric decrement of ionic solutions, the screening length can increase with the ionic concentration only due to a reduction in the amount of free ions.

The reduction of free-ion concentration at low $\phi$ values can have two origins. Steric effects can be more important than what the lattice-gas framework conveys. For example, it is possible to take steric effects into account by using a virial expansion of hard spheres, where the maximal value of the free ionic concentration, $n_{s}$, is calculated at $\phi\approx 0.1$ Sam17 . In addition, correlations beyond MF between ions and dipoles have a crucial role in highly concentrated ionic solutions and ionic liquids, and may promote the reduction of $n_{s}$. This possibility will be explored in the future by employing a loop expansion of the system free-energy, leading to predictions beyond the MF theory.

*Acknowledgments. * This work was motivated by discussions with S. Safran. We thank him and R. Podgornik for numerous helpful suggestions and comments. This work was supported by the Israel Science Foundation (ISF) under grant No. 438/12, the U.S.- Israel Binational Science Foundation (BSF) under grant No. 2012/060, and the ISF-NSFC joint research program under grant No. 885/15. D.A. thanks Alexander von Humboldt Foundation for a Humboldt research award, and T.M. acknowledges support from the Blavatnik postdoctoral fellowship programme.

Appendix A Field-theoretical approach

Consider the aqueous solution of Section II modeled via a cubic lattice with cells of size $a$. Each cell is occupied either by a solvent molecule with dipole moment $p_{w}$, a cation of charge $e$, an anion of charge $-e$, or a Bjerrum pair with dipole moment $p$. The association energy of pairs is $-J,$ where $J$ is a positive energy parameter.

The partition function can be written in terms of spin-like variables assigned to each cell, accounting for its occupation. Each cell $j,$ located at $\boldsymbol{r}_{j},$ is described by a pair of variables, $\left(s_{j}^{+},s_{j}^{-}\right),$ counting the cations and anions in the cell, respectively, where $s_{j}^{\pm}=0\,{\rm or}\,1$. With these variables, the charge density operator is given by

[TABLE]

where $\rho_{f}({\bf r})$ (the first term) describes any possible fixed charge density, the 2nd term corresponds to the cations and anions, and the 3rd term accounts either for the presence of a water dipole or a Bjerrum-pair dipole in the $j$-cell, pointing in the direction of the unit vector, $\boldsymbol{\hat{n}}_{j}$.

Incorporating the Coulombic energy and the ion-association energy yields the following grand-partition function:

[TABLE]

where $\mu_{\pm}$ is the chemical potential of the positive and negative ions, respectively, $\beta=\left(k_{B}T\right)^{-1}$ is the inverse thermal energy, $\Omega_{j}$ is the solid angle of $\boldsymbol{\hat{n}}_{j}$, and $v_{c}\left({r}\right)=1/\left(4\pi\varepsilon_{0}\left|{r}\right|\right)$ is the Coulomb interaction. For sake of convenience, we have set the chemical potential of solvent molecules to be zero. The Coulombic self-energy can be formally absorbed in the chemical potentials and $J$, overcoming the problematic divergence of $v_{c}\left(0\right).$

We replace the occupation degrees of freedom, $\{s_{j}\}$, with a spatially fluctuating field by introducing a density field $\rho\left(\boldsymbol{r}\right)$ and its conjugate field, $\varphi\left(\boldsymbol{r}\right),$ via the functional identity:

[TABLE]

Substituting Eq. (A) in Eq. (A) and after further manipulations (for more details, see Ref. Borukhov00 ), the partition function can be written as

[TABLE]

with the field action

[TABLE]

In the mean field (MF) approximation, the partition function is approximated by its value at the saddle point, $\varphi=\varphi_{0}.$ We denote $i\varphi_{0}=\psi$, and by using the relation $F=-k_{B}T\ln Z$, we obtain the free energy

[TABLE]

By examining Eqs. (A)-(22) with $\Lambda_{\pm}=0$ and $p_{w}=0$, it can be seen that $\psi$ is the electrostatic potential. Furthermore, the fugacities are constant throughout the system and are related to the total number densities of positive and negative ions, as is described in Section II.1.

The potential $\psi$ solves the Euler-Lagrange equation obtained from the variation $\delta F/\delta\psi=0.$ Assuming that $\rho_{f}$ depends only on the $z$-coordinate, the Euler-Lagrange equation for the free energy, Eq. (22), coincides with Eq. (II.2), and is a more formal way to obtain the MDPB equation.