Structural interactions in ionic liquids linked to higher order Poisson-Boltzmann equations

TL;DR

This paper derives generalized Poisson-Boltzmann equations from classical binary fluid mixture theories, linking microscopic interactions to macroscopic electrostatic models with potential applications in understanding ionic liquids.

Contribution

It introduces a novel derivation of higher-order Poisson-Boltzmann equations using Legendre transforms, connecting microscopic fluid structure to electrostatic modeling.

Findings

Reduction to a phenomenological equation under symmetry assumptions

Structural interactions influence near-surface ion distributions

Provides a theoretical framework for higher-order electrostatic equations

Abstract

We present a derivation of generalized Poisson-Boltzmann equations starting {from} classical theories of binary fluid mixtures, employing an approach based on the Legendre transform as recently applied to the case of local descriptions of the fluid free energy. Under specific symmetry assumptions, and in the linearized regime, the Poisson-Boltzmann equation reduces to a phenomenological equation introduced by (Bazant et al., 2011), whereby the structuring near the surface is determined by bulk coefficients.

Equations

$${\mathcal{F}}_{electro}[\{c_i\},\mathbf{D}]=\int {\mathrm{d}}^3\mathbf{r}\left[\frac{{\mathbf{D}}^2}{2\epsilon }-\varphi (\text{mbox}div\mathbf{D}-e{z}_i{c}_i)\right]$$

$${\mathcal{F}}_{fluid}[\{c_i\}]=\frac{1}{2}\int {\mathrm{d}}^3\mathbf{r}\left[f(\{c_i\})-{\mu }_i{c}_i+{\kappa }_{ij}\nabla c_i\nabla c_j\right]$$

$${\mathcal{F}}_{fluid}[\{c_i\}]=\int {\mathrm{d}}^3\mathbf{r}f(w(c_i))$$

$$\mathcal{F}[\{c_i\},$$

$$\int {\mathrm{d}}^3\mathbf{r}\left(f(\{c_i\})-({\mu }_i+e{z}_i\varphi +\nabla \cdot {\mathbf{g}}_i)c_i\right)$$

$$+$$

$$\mathcal{F}[\varphi ,\{{\mathbf{g}}_i\}]$$

$$-$$

$$c_i=\frac{\partial P}{\partial {\mu }_i}.$$

$$p=-ϵ\frac{({\partial }_x\varphi {)}^2}{2}+P+\frac{g_i{g}_j}{2{\kappa }_{ij}}-({\partial }_x{g}_i)c_i$$

$$\mathcal{F}=\int {\mathrm{d}}^3\mathbf{r}\left[-{\mu }_i{c}_i+f(w_i\ast c_i)+\frac{{\mathbf{D}}^2}{2\epsilon }-\varphi (\mathrm{div}\mathbf{D}-e{z}_i{c}_i)\right]$$

$$\mathcal{F}=\int {\mathrm{d}}^3\mathbf{r}[-{\mu }_i{c}_i+f(d_i)+$$

$$-$$

$$\mathcal{F}=-\int {\mathrm{d}}^3\mathbf{r}\left[\epsilon \frac{{(\nabla \varphi )}^2}{2}-P(\{w_i^{-1}(\mu -e{z}_i\varphi ))\}\right]$$

$$\mathcal{F}=\int {\mathrm{d}}^3\mathbf{r}\left[-\epsilon \frac{{(\nabla \varphi )}^2}{2}-P\left(\left\{(1-\frac{{\nabla }^2}{k_i^2})({\mu }_i-e{z}_i\varphi )\right\}\right)\right].$$

$$F\approx \int {\mathrm{d}}^3\mathbf{r}\left(-P-c_i\mathrm{div}{\mathbf{g}}_i-ϵ\frac{(\nabla \varphi {)}^2}{2}-\frac{{\mathbf{g}}_i{\mathbf{g}}_j}{2{\kappa }_{ij}}\right)$$

$${\mathbf{g}}_i\approx {\kappa }_{ij}\nabla c_j$$

$$\mathcal{F}=$$

$$=$$

$$\mathcal{F}[$$

$$-$$

$$\displaystyle\overline{g}\cdot\widehat{\kappa}\cdot\overline{g}^{T}=\left(g_{+}\,\,\,g_{-}\right)\left(\begin{array}[]{cc}\widehat{\kappa}_{++}&-\widehat{\kappa}_{+-}\\ -\widehat{\kappa}_{+-}&\widehat{\kappa}_{--}\end{array}\right)\left(\begin{array}[]{c}g_{+}\\ g_{-}\end{array}\right)\,.$$

$$P=\frac{1}{a^3}\ln \left(1+\frac{\gamma }{2(1-\gamma )}\left(e^{\beta m_{+}(x)}+e^{\beta m_{-}(x)}\right)\right)$$

$$m_{\pm }(x)=\pm ez\varphi (x)+{\partial }_x{g}_{\pm }(x),$$

$$-\epsilon {\partial }_x^2\varphi (x)+\frac{\gamma ze}{a^3}\left[\frac{e^{\beta m_{+}(x)}-e^{\beta m_{-}(x)}}{1-\gamma +\gamma (e^{\beta m_{+}(x)}+e^{\beta m_{-}(x)})}\right]=0,$$

$${\hat{\kappa }}_{++}{g}_{+}-{\hat{\kappa }}_{+-}{g}_{-}=\frac{\gamma }{a^3}{\partial }_x\left[\frac{e^{\beta m_{+}(x)}}{1-\gamma +\gamma (e^{\beta m_{+}(x)}+e^{\beta m_{-}(x)})}\right]$$

$${\hat{\kappa }}_{--}{g}_{--}{\hat{\kappa }}_{+-}{g}_{+}-=\frac{\gamma }{a^3}{\partial }_x\left[\frac{e^{\beta m_{-}(x)}}{1-\gamma +\gamma (e^{\beta m_{+}(x)}+e^{\beta m_{-}(x)})}\right].$$

$$-\epsilon {\partial }_x^2\varphi (x)+\frac{\gamma ze\beta }{a^3}(2ez\varphi (x)+{\partial }_x{g}_{+}(x)-{\partial }_x{g}_{-}(x))=0,$$

$${\hat{\kappa }}_{++}{g}_{+}(x)-{\hat{\kappa }}_{+-}{g}_{-}(x)=\frac{\gamma \beta }{a^3}(ez{\partial }_x\varphi (x)+{\partial }_x^2{g}_{+}(x))$$

$${\hat{\kappa }}_{--}{g}_{-}(x)-{\hat{\kappa }}_{+-}{g}_{+}(x)=\frac{\gamma \beta }{a^3}(-ez{\partial }_x\varphi (x)+{\partial }_x^2{g}_{-}(x)).$$

$${\tilde{\kappa }}^2{\partial }_x^4\psi (x)-{\partial }_x^2\psi (x)+\psi (x)=0,$$

$$-\epsilon {\partial }_x^2\varphi (x)+\frac{\gamma ze}{a^3}\frac{\sinh (\beta m(x))}{1+2\gamma {\sinh }^2(\beta m(x)/2)}=0$$

$$\overline{g}(x)=\frac{\gamma }{2a^3\kappa }{\partial }_x\left(\frac{\sinh (\beta m(x))}{1+2\gamma {\sinh }^2(\beta m(x)/2)}\right).$$

Full text

Introduction

Room-temperature ionic liquids display densely packed layers with periodically varying charge next to charged substrates, described in Mezger et al. (2008), and earlier e.g. Atkin and Warr (2007). This phenomenon has been observed in numerous experiments and simulations, see e.g. Perkin et al. (2010); Zhou et al. (2012); Smith et al. (2013); Cheng et al. (2016); Smith et al. (2017); for a review, see Varela et al. (2015).

As argued already in Mezger et al. (2008) the quantitative description of the electrostatic properties of these layers, with their observed unusual packing of cation and anions, cannot be covered by the standard mean-field approach based on the lattice-gas Poisson-Boltzmann equation. What is necessary is its generalization either in terms of structural contributions to the free energy, leading to micro phase separation and charge layering Gavish and Yochelis (2016), or solving the lattice Coulomb gas model beyond the mean-field approximation Démery et al. (2012, 2012). Recently, however, a phenomenological Poisson-Boltzmann equation of fourth-order has been postulated and applied to the structure of double layers in ionic liquids Bazant et al. (2011). The origin of the higher order term is argued to be due to correlations between the ions, an effect known to fundamentally affect the behavior of ionic fluids Naji et al. (2013). Fourth-order Poisson-Boltzmann equations have indeed appeared in other contexts, e.g. by invoking additional ad-hoc parameters in order to separate out the ion correlation effects on different length scales Santangelo (2006). They can also appear in phenomenological nonlocal theories, see e.g. Paillusson and Blossey (2010), where they can arise from specific forms of the wave-vector dependent dielectric function.

The arguments leading to a simple fourth-order Poisson-Boltzmann equation in the context of ionic liquids have been, however, somewhat sketchy. The atomistic simulations of Zhou et al. (2012), compared to this theory, reveal deviations which should not be present if the theory can correctly account for the correlation effects. Despite the substantial body of work that has emerged meanwhile on the topic, a better understanding of the underlying physics is certainly warranted.

Recently, two of us (ACM, RP) have studied (asymmetric) steric interactions in electrostatic double layers with Legendre transform methods Maggs and Podgornik (2016). The starting point of this work is the free energy of a homogeneous binary charged mixture, from which generalized Poisson-Boltzmann equations can be systematically derived in a form, valid for any assumed equation of state of the uncharged fluid. The resulting Poisson-Boltzmann equations, however, remain partial differential equations of the second order, but with different types of local nonlinearities which cover steric effects, see also the earlier work Borukhov et al. (1997, 2000). Rather surprisingly, arbitrary fluid mixtures were shown to give rise to Poisson-Boltzmann equations, in which one requires only knowledge of the non-electrostatic fluid pressure as a function of the chemical potentials of the individual components.

In this work, we present a formal derivation of generalized Poisson-Boltzmann equations, where the appearance of higher-order derivative terms in the local density is directly linked to the spatial variation in the concentration of the ions in the binary mixture. It arises both within a squared-gradient approach or a more general weighted density approximation in the free energy density. Interestingly, the equation governing the electrostatic potential appears, to lowest relevant order, as a fourth-order partial differential equation whose coefficients, however, in general depend on the non-electrostatic equation of state. As a first application of this type of equation we compare its solutions to those of the phenomenological equation employed by Bazant et al. Bazant et al. (2011). We show that even with these generalisations the only input needed by the theory is again the pressure function of the non-electrostatic problem. The phenomenological coefficients of Bazant et al. (2011) are found in terms of this non-electrostatic pressure function.

Theory

The starting point of our discussion is the free energy expression for a charged binary mixture, composed of two contributions, ${\cal F}={\cal F}_{electro}+{\cal F}_{fluid}$, where the electrostatic contribution is given by the expression

[TABLE]

while the fluid term can be either the sum of a local free energy plus a squared gradient

[TABLE]

where the interaction strengths $\kappa_{ij}$ are proportional to the bulk correlation lengths, or can be written in a weighted density form

[TABLE]

where $w$ is supposed to be a linear, nonlocal function Curtin and Ashcroft (1985). Note that, in the expressions, we have employed Einstein’s summation convention and that in what follows, we will be rather careless in the exact symbol (and name) used for the various free energies $\mathcal{F}$. We perform multiple and sometimes non-standard Legendre transforms, while the stationary points of all the objects that we construct have the same value.

The local contribution to eq.xs (2) is the free energy of an isothermal binary mixture with the concentrations $c_{i}$ and was the basis of the discussion that in Maggs and Podgornik (2016) led to a general theory of asymmetric steric interactions in electrostatic double layers; the $\mu_{i}$ are the chemical potentials; ${\bf D}$ is the dielectric displacement field, $\varepsilon=\epsilon\epsilon_{0}$ the relative dielectric permeability, the $z_{i}$ the valencies of the charged species and $\phi$ is the Lagrange multiplier which takes care of Gauss’ law, i.e. the electrostatic potential.

In Maggs and Podgornik (2016), the theory on the mean-field level was transformed into a free energy expression for the electrostatic potential $\phi({\bf r})$. This was achieved by making use of the Legendre transform approach developed in Maggs and Everaers (2006); Pujos and Maggs (2014). The same approach can be applied here, where in addition we further first introduce a vector-valued variable for the concentration gradient given by ${\bf v}_{i}\equiv\nabla c_{i}\,$ with an associated multiplier ${\bf g}_{i}$. An integration by parts and regrouping of terms in eq. (2) leads to the expression

[TABLE]

We now recognise that the stationary point of this functional with respect to the concentration leads to the Legendre transform of the fluid free energy, thus

[TABLE]

where $P$ is the fluid pressure expressed as a function of the set of chemical potentials ${\mu_{i}}$, as discussed in Maggs and Podgornik (2016). It is important to note that the Lagrange multiplier ${\bf g}_{i}$ enters as an argument to the pressure function which is in general a highly nonlinear function. We recall at this point the fundamental Gibbs-Duhem relation linking derivatives of $P$ to the concentration,

[TABLE]

The variational equations found by taking derivatives of the free energy eq. (5) have a simple first integral in one-dimensional geometries. This integral can be found using the standard construction of a Hamiltonian from a Lagrangian, taking the coordinate $x$ as equivalent to the time in a particle system, while the conserved quantity is not an energy, but rather the total pressure in the fluid (where external charges are absent) which can be derived as

[TABLE]

Far from any external sources, where the potential and $g_{i}$ become constant, this reduces simply to the neutral fluid pressure function $P$.

We now take a fluid free energy based on the weighted functional form eq. (3). The weighted density approximation is a more sophisticated approximation than the square gradient approximation. It supposes that one can evaluate the free energy of a non-uniform fluid from the equation of state of a uniform system. However the argument of the uniform state function is a non-local average of the density in some local neighbourhood. Typically such methods give more robust descriptions of fluid behaviour, above all in the presence of strong gradients. In the most sophisticated theories the weighting kernel is also calculated from the correlations of the equilibrium fluid but here we will take a very simplified form of the kernel as a Yukawa form, with range $k_{i}^{-1}$. This specific choice allows us to make considerable simplifications in the theory. The total free energy of the system is then

[TABLE]

Where by $w_{i}\ast c_{i}$ we mean the convolution of the concentration field $c_{i}$ with the Yakawa kernel of range $k_{i}^{-1}$. We pull out the weighting function as an argument of $f$ by introducing a local density $d_{i}=w_{i}\ast c_{i}$:

[TABLE]

where again $\lambda_{i}$ is a Lagrange multiplier. All the concentrations $c_{i}$ occur linearly allowing us to perform the Legendre transform:

[TABLE]

It is here that the utility of the choice of a Yukawa function becomes clear as the inverse operator is then $w^{-1}=1-\nabla^{2}/k_{i}^{2}$. This gives the final Poisson-Boltzmann functional for the weighted density approximation that involves purely the electrostatic potential, without the introduction of any supplementary degrees of freedom:

[TABLE]

Within this local density approximation we have therefore derived a generalisation of the Poisson-Boltzmann equation with a Laplacian within the pressure function, $P$. We note that eq. (10) leads, on expansion in gradients, to a free energy which is a series in $(\nabla^{2}\phi)$. Retaining terms up to second order leads directly to the free energies of the form proposed in Bazant et al. (2011).

Lowest order elimination of the gradient multiplier

In eq. (4) we can expand the pressure function to first order in $({{\rm div}\;}{\bf g}_{i})$ to find:

[TABLE]

We now integrate by parts and solve for ${\bf g}_{i}$

[TABLE]

Substituting back and using Gibbs-Duhem gives an effective action in terms of only a potential variable:

[TABLE]

since $c_{l}$ is clearly itself a derivative of the pressure. The lowest order correction to the Poisson-Boltzmann function has the appearance of an effective shift in the dielectric properties of the fluid or indeed implies that the molecular structure of the system renormalizes its dielectric response.

Application: The symmetric electrolyte and the

square gradient approximation in one dimension

We now develop further the application of eq. (5), which unlike eq. (10) is not directly of the Poisson-Boltzmann form, involving only a functional of the electrostatic potential. We work with a symmetric binary mixture and impose a one-dimensional geometry, having in mind one and two-plate problems. The electroneutrality condition is trivially satisfied since no external charge is considered in the model.

In this case the free energy expression, eq (5) reduces to

[TABLE]

with $i=(+,-)$, where, explicitly, we have for the inverse matrix $\kappa_{ij}^{-1}\equiv\widehat{\kappa}_{ij}$

[TABLE]

The standard expression for the pressure comes from a lattice gas model:

[TABLE]

with

[TABLE]

and, following ref. Maggs and Podgornik (2016), the conditions on the chemical potentials are $e^{\beta\mu_{+}}=e^{\beta\mu_{-}}=\frac{1}{2}[\frac{\gamma}{1-\gamma}]$. We can now vary the functional ${\cal F}[\phi,\{\partial_{x}g_{i}\}]$ with respect to the fields $\phi,g_{+}$ and $g_{-}$. The variation with respect to $\phi$ yields

[TABLE]

while for the two components of $\overline{g}$ one has

[TABLE]

and

[TABLE]

Linearisation. It is instructive to consider the linearized equations — we will see that they contain already the essential physics of this theory. We find

[TABLE]

while for the two components of $\overline{g}$ one has

[TABLE]

and

[TABLE]

Figure 1 shows the resulting solutions to these equations in the case of a one-plate system with constant surface potential, $\phi(0)=V$, for the parameters indicated in the caption. In Figure 1 A, $\kappa_{++}=\kappa_{--}$, while the contrast between the two couplings in increased by a factor of ten in Figure 1 B. The presence of $g_{+}$ and $g_{-}$ leads to oscillations in the profile of the electrostatic potential whose number, amplitude and extension into the bulk depend on the contrast between the coupling parameters. The equation for the electrostatic potential $\phi$ looks similar to that if a mechanical spring (in which the position $x$ plays the role of time), with the damping term given by the functions $g_{+}(x)$ and $g_{-}(x)$, which in their turn essentially obey coupled diffusion equations, with the electric field acting as a source or a sink term. It is therefore not surprising that, as a function of the coupling parameters $\widehat{\kappa}_{ij}$, the solution exhibits damped oscillations.

The variational equations can be reduced to only two fields if one assumes the symmetry condition $\widehat{\kappa}_{++}=\widehat{\kappa}_{--}\equiv\kappa$. We further set $\widehat{\kappa}_{+-}=0$. Subtracting the equations in $g_{\pm}$ and introducing $g\equiv g_{+}-g_{-}$ then yields a coupled system in $\phi(x),g(x)$. Being linear equations, one can eliminate the field $g$ in favor of $\phi$ and ends up with a forth-order equation in the electrostatic field

[TABLE]

where $\psi\equiv(2ez)\phi$, $\tilde{\kappa}^{2}\equiv(2(ze)^{2}(\gamma\beta)^{2}/(\kappa\varepsilon a^{6})$ and the spatial coordinate $x$ has been rescaled by a factor $\sqrt{\epsilon}$ with $\epsilon=\varepsilon a^{3}/(\gamma\beta)$.

Eq. (28) is exactly the linear equation discussed in Bazant et al. (2011), identifying the coupling parameter $\tilde{\kappa}^{2}$ with the parameter $\delta_{c}^{2}$ used in Bazant et al. (2011). This connection also allows us to relate the parameter $\delta_{c}^{2}$ introduced in Bazant et al. (2011) with the interaction strength $\kappa$ in the binary mixture. Furthermore, the linear solution derived in the Supplementary Material of Bazant et al. (2011) can immediately be adopted by just making the replacement $\delta_{c}^{2}$ with $\tilde{\kappa}^{2}$.

The nonlinear case. In the nonlinear case, one has to solve the coupled eqs. (25) - (27). A simplification can again be made if we consider the anti-symmetric subspace $g_{+}(x)=-g_{-}(x)\equiv\overline{g}(x)$, corresponding to antagonistic gradients in concentrations. In this case one has $m_{+}(x)=-m_{-}(x)\equiv m(x)$ and the resulting two equations for $\phi$ and $g$ can be written as

[TABLE]

and

[TABLE]

The full nonlinear equations can easily be solved numerically as a boundary value problem with the conditions $\phi(0)=V$, $\phi(L)=0$, $\overline{g}(0)=\overline{g}(L)=0$. Figure 2 shows the results for the charge density $P^{\prime}(\phi)\equiv(1/\gamma)\tanh(\gamma\phi(x))$, for a value of $\gamma=0.8$ and a value of $\phi(x=0)=V=5$. The results compare favorably with those obtained in Bazant et al. (2011), and therefore indicate that our explicit assumption of antagonistic concentration gradient is implicitly present in the theory developed there.

The two-plate problem. Finally one can also study the two-plate problem for this equation in the similar vein as Gavish and Yochelis (2016). Figure 3 plots an exemplary osmotic pressure, given by eq.(7), between the plates for identical constant potentials at the plates, $\phi(0)=\phi(L)$, based on the linearized two-field theory discussed just before. Also shown is a comparison with the osmotic pressure of the linearized standard Poisson-Boltzmann equation, which yields a parabolic potential $\phi(x)>0$ with an exponentially decaying osmotic pressure for large plate distances. In the structured fluid, upon variation of the two-plate distance $L$, the potential $\phi$ first becomes flat at the center and crosses to negative values (as shown in an insert for $L=10$), resulting in a pronounced minimum in $p(L)$ (not shown). Upon further increase in $L$ the osmotic pressure progresses through a maximum and a second, very shallow minimum appears (see top insert). In this region, the electrostatic potential has a double-well structure next to the plates, as shown in the insert for $L=15$.

Discussion

In this work we have shown that the inclusion of concentration gradients in a charged binary mixture leads in a natural way to the appearance of higher-order terms in the generalized Poisson-Boltzmann equation, which in general turns out to be much more complex than the phenomenological fourth-order PB equation employed in Bazant et al. (2011) due to the coupling of nonlinearities and spatial gradients.

In the model case of a symmetric binary mixture the presence of different coupling constants of the spatial gradient terms lead to pronounced oscillatory behaviour of the electrostatic profiles. Also, on the level of the generalized Poisson-Boltzmann equation, the oscillations near the wall therefore are the consequence of bulk behaviour and not of the surface coupling. Nonetheless, the solution of the corresponding system of coupled nonlinear differential equations for the electrostatic potential and the spatial gradient function yield results that are qualitatively similar to the behaviour found in Bazant et al. (2011), and actually coincide for the linear case, providing a physically more sound justification of the phenomenological approach. At the same time, our results make clear that the present state of theory does not go beyond the mere mean-field level, and a consistent discussion of correlation effects in ionic liquids, in addition to the structure effects, is still missing. A description of the oscillations in electrostatic potential and osmotic pressure in terms of two decay lengths, as discussed in Smith et al. (2017), might therefore be only a special case of a more general situation in which structural, or packing, effects play a crucial role.

Nevertheless, it is worthwhile pointing out that the squared density gradient theory, stemming from the non-electrostatic structural contribution to the free energy, leads at least in the linear regime formally to the same generalized Poisson-Boltzmann equation as one would obtain from the contribution of the non-mean-field ion correlations, approximated via an appropriate length scale separation Santangelo (2006). It would therefore seem difficult to disentangle these two effects without a full theory taking into account the ion correlations as well as the molecular structure of the system.