Poisson-Nernst-Planck equations with steric effects - non-convexity and multiple stationary solutions

TL;DR

This paper investigates the existence and stability of multiple stationary solutions in Poisson-Nernst-Planck equations with steric effects, revealing instability in the classical model and proposing a new stable, multi-solution framework.

Contribution

It demonstrates the instability of multiple solutions in PNP-steric equations and introduces a novel PNP-Cahn-Hilliard model with stable, multi-stable solutions.

Findings

Steric effects lead to multiple stationary solutions in PNP equations.

Classical PNP-steric solutions are unstable under dynamics.

The PNP-Cahn-Hilliard model admits stable, multi-stable solutions.

Abstract

We study the existence and stability of stationary solutions of Poisson-Nernst- Planck equations with steric effects (PNP-steric equations) with two counter-charged species. These equations describe steady current through open ionic channels quite well. The current levels in open ionic channels are known to switch between `open' or `closed' states in a spontaneous stochastic process called gating, suggesting that their governing equations should give rise to multiple stationary solutions that enable such multi-stable behavior. We show that within a range of parameters, steric effects give rise to multiple stationary solutions that are smooth. These solutions, however, are all unstable under PNP-steric dynamics. Following these findings, we introduce a novel PNP-Cahn-Hilliard model, and show that it admits multiple stationary solutions that are smooth and stable. The various branches of stationary solutions and their stability are mapped utilizing bifurcation analysis and numerical continuation methods.

Figures

Equations

$$\begin{array}{rl}& \frac{\partial c_i}{\partial t}+\nabla \cdot J_i=0,i=1,2,\cdots ,N,\\ & {\varphi }_{xx}+\sum _{i=1}^N{z}_i{c}_i=0,\end{array}$$

$$J_i=-c_i\nabla \left(\log c_i+z_i\varphi +\sum _{i=1}^N{g}_{ij}{c}_j\right).$$

$$\log c_i+z_i\varphi +\sum _{i=1}^2{g}_{ij}{c}_j={\lambda }_i,i=1,\cdots ,N,$$

$$z_2<0<z_1.$$

$$\varphi (\pm L)={\varphi }_{\pm L},$$

$$z_1{\bar{c}}_1+z_2{\bar{c}}_2+{\rho }_0=0,$$

$${\bar{c}}_i:=\frac{1}{2L}{\int }_{-L}^L{c}_i(x)dx,$$

$$\mathcal{A}(c_1,c_2,\varphi )={\int }_{-L}^L{k}_{B}T\underset{\mathrm{entropy}}{\underbrace{\sum _{i=1}^2{c}_i\left(\ln \frac{c_i}{{\bar{c}}_i}-1\right)}}+\underset{\mathrm{electrostatic}}{\underbrace{q{\mathbf{z}}^T\mathbf{c}\varphi -\frac{ϵ}{2}|\nabla \varphi {|}^2}}+\underset{\mathrm{steric}}{\underbrace{\frac{1}{2}{\mathbf{c}}^T\mathbf{G}\mathbf{c}}}dx,$$

$$0=\frac{\delta \mathcal{A}}{\delta \varphi }=\epsilon {\varphi }_{xx}+q\left({\mathbf{z}}^T\mathbf{c}+{\rho }_0\right)=0.$$

$$\frac{d{c}_i}{dt}=\breve{D}\nabla \cdot \left(c_i\nabla \frac{\delta \mathcal{A}}{\delta c_i}\right),i=1,2,$$

$$\frac{d{c}_i}{dt}+\nabla \cdot J_i=0,J_i=-\breve{D}{c}_i\nabla \frac{\delta \mathcal{A}}{\delta c_i}=-\breve{D}{c}_i\nabla \left(k_{B}T\log c_i+z_{i}q\varphi +\sum _{j=1}^2{g}_{ij}{c}_j\right).$$

$$\tilde{\varphi }=\frac{q}{k_{B}T}\varphi ,\tilde{x}=\frac{x}{{\lambda }_D},{\tilde{c}}_i=\frac{c_i}{c^o},{\widetilde{\stackrel{ˉ}{c}}}_i=\frac{{\bar{c}}_i}{c^o},\tilde{\mathbf{G}}=\frac{1}{k_{B}T}\mathbf{G},\tilde{t}=\frac{{\lambda }_D^2}{\breve{D}{k}_{B}T}t,$$

$${\lambda }_D:=\sqrt{\frac{k_{B}T\epsilon }{q^2{c}^o}},$$

$$\begin{array}{rl}& \frac{d{c}_i}{dt}+\nabla \cdot J_i=0,J_i=-c_i\nabla \left(\log c_i+z_i\varphi +\sum _{j=1}^2{g}_{ij}{c}_j\right),\\ & {\varphi }_{xx}=-\left({\mathbf{z}}^T\mathbf{c}+{\rho }_0\right),\end{array}$$

$$\varphi (\pm L)={\varphi }_{\pm L},{J_i|}_{\pm L}=0.$$

$$\frac{\delta \mathcal{A}}{\delta c_i}=\log c_i+z_i\varphi +\sum _{j=1}^2{g}_{ij}{c}_j={\lambda }_i(\bar{\mathbf{c}})=\log {\bar{c}}_i+\sum _{j=1}^2{g}_{ij}{\bar{c}}_j,i=1,2,$$

$${\varphi }_{xx}=-\left({\mathbf{z}}^T\mathbf{c}+{\rho }_0\right),$$

$$\varphi (\pm L)={\varphi }_{\pm L}.$$

$$\mathcal{E}(c_1,c_2,\varphi )={\int }_{-L}^L\underset{\mathrm{entropy}+\mathrm{steric}}{\underbrace{h(c_1,c_2)}}+\underset{\mathrm{electrostatic}}{\underbrace{\frac{ϵ}{2}|\nabla \varphi {|}^2}}dx,$$

$$h(\mathbf{c})=\sum _{i=1}^2{c}_i\left(\ln \frac{c_i}{{\bar{c}}_i}-1\right)+\frac{1}{2}{\mathbf{c}}^T\mathbf{G}\mathbf{c}.$$

$$h_{c_1}(\mathbf{c})+z_1\varphi =h_{c_1}(\bar{\mathbf{c}}),h_{c_2}(\mathbf{c})+z_2\varphi =h_{c_2}(\bar{\mathbf{c}}),$$

$$\left[\begin{array}[]{c}\phi\\ E\end{array}\right]_{x}=\left[\begin{array}[]{c}E\\ -{\bf z}^{T}{\bf c}(\phi;\bar{\bf c})-\rho_{0}\end{array}\right],$$

$$\text{mbox}diag(\mathbf{c})\text{mbox}{Hessian}_{\mathbf{c}}(h)=I+\text{mbox}diag(\mathbf{c})\mathbf{G}.$$

$$\mathbf{G}{v}_{-}={\lambda }_{-}{v}_{-},{\lambda }_{-}<0.$$

$$\text{mbox}diag({\mathbf{c}}_{-})\mathbf{G}{v}_{-}=-v_{-}$$

$$E_x=-[z_1(c_1-{\bar{c}}_1)+z_2(c_2-{\bar{c}}_2)],$$

$$\begin{array}{rl}{c}_{1,x}& =-\frac{z_1(1+c_2{g}_{22})-z_2{c}_2{g}_{12}}{c_{2}D}E,\\ c_{2,x}& =-\frac{z_2(1+c_1{g}_{11})-z_1{c}_1{g}_{12}}{c_{1}D}E,\end{array}$$

$$D=|\text{mbox}{Hessian}_{\mathbf{c}}(h)|=(1/c_1+g_{11})(1/c_2+g_{22})-g_{12}^2=\frac{1}{c_1{c}_2}+\frac{g_{11}}{c_1}+\frac{g_{22}}{c_2}+|\mathbf{G}|.$$

$$\frac{d{c}_1}{d{c}_2}=\frac{z_1\left(\frac{1}{c_2}+g_{22}\right)-z_2{g}_{12}}{z_2\left(\frac{1}{c_1}+g_{11}\right)-z_1{g}_{12}},c_1(c_2^0)=c_1^0.$$

$$m{c}_1\le \frac{d{c}_1}{d{c}_2}\le M{c}_1<0,c_2>c_2^0,m=\frac{z_1}{z_2}\left(\frac{1}{c_2^0}+g_{22}\right)-g_{12},M=\frac{z_1{g}_{22}-z_2{g}_{12}}{z_2\left(1+g_{11}{c}_1^0\right)-z_1{g}_{12}{c}_1^0}.$$

$$0<c_1^0\exp [m(c_2-c_2^0)]\le c_1(c_2)\le c_1^0\exp [M(c_2-c_2^0)].$$

$$c_1={\bar{c}}_1-\frac{z_2}{z_1}(c_2-{\bar{c}}_2).$$

$$\underset{c_2\to \infty }{\lim }{c}_1(c_2)=0,$$

$$\underset{c_2\to \infty }{\lim }{c}_1(c_2)-{\Gamma }_E(c_2)=-\infty .$$

Full text

1 Introduction

It is difficult to overstate the importance of electrolyte solutions in biological and electrochemical systems. The pioneering studies of Nernst and Planck in 1889-1890 [1, 2] led to the Poisson-Nernst-Planck (PNP) theory, which is one of the most widely used continuum mean-field model to describe electrolyte behavior. PNP theory, however, is valid only at exceedingly dilute solutions due to various underlying assumptions. In particular, the model assumes point-like ions, while neglecting crowding effects due to the finite-size of the ions. The shortcomings of the PNP model has led to the development of a wide family of generalized Poisson-Nernst-Planck models that take into account ion-ion and/or ion-solvent interactions [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]; for an overview the reader is also referred to [24, 25] and references within. Among these models, we particularly focus here on the PNP-steric [26, 13] model

[TABLE]

The above system is written using non-dimensional variables, where $c_{i}$ are the concentrations of the $i^{\rm th}$ ionic species with corresponding valence $z_{i}$, $\phi$ is the electric potential, and $g_{ij}=g_{ji}$ are non-negative constants depending on ion radii and the energy coupling constant of the $i^{\rm th}$ and $j^{\rm th}$ species ions [13]. When $g_{ij}=0$ for all $1\leq i,j\leq N$, the PNP-steric equations reduces to the classical PNP model that does not account for finite-size effects.

A wide family of generalized PNP models share common features - these models give rise to unique stationary solutions which in the absence of applied external potential reduces to a homogenous state [24, 27]. The uniqueness of a stationary solution, however, is challenged by observations of ionic currents through ionic channels. Indeed, experiments [28] and simulations [29] measuring currents through single protein channels show that single channels switch between ‘open’ or ‘closed’ states in a spontaneous stochastic process called gating. Currents are either (nearly) zero or at a definite level, characteristic of each type of protein, independent of time, once the channel is open. These observations of a bi-stable system strongly suggest that an underlying model for ion channels should have multiple stationary solutions, one corresponding to an ‘open’ state and another to a ‘closed’ state. Following this insight, the recent work [30] due to Lin and Eisenberg considered the existence of multiple stationary solutions of the (non-dimensional) PNP-steric equations (1) in a general (namely, not in an ion-channel) setting.

We note that PNP equations with piecewise constant fixed charge can give rise to multiple solutions, see, e.g., [31, 32, 33]. The work [30], as well as the current study, is concerned whether the PNP-steric equations with spatially constant permanent charges can give rise to multiple steady-states. Namely, whether steric effects themselves are sufficient to give rise to multiple stationary solutions. The study [30] showed that when $N=2$, $z_{1}>0>z_{2}$ and $g_{11}=g_{22}$, the PNP-steric system gives rise to a single stationary solution, but proved the existence of multiple stationary solutions in multi-species systems, $N\geq 3$. This work, however, considered solely the stationary equations satisfying

[TABLE]

where $\lambda_{i}$ is a Lagrange multiplier associated with the conservation of mass, coupled with Poisson’s equation. Particularly, these studies characterized branches of solutions that are functions of $\phi$, i.e. $c_{i}=c_{i}(\phi)$, defined by the algebraic system (2). Dynamics were not considered in these studies. It, therefore, remains an open question whether the PNP-steric model with two species admits smooth multiple stationary solutions, and whether these solutions would be stable under PNP-steric dynamics. These open questions are addressed in this study.

1.1 Main results and paper outline

In this work, we study the PNP-steric system with two counter-charged species. We show that under certain conditions, particularly, for sufficiently concentrated (large $L^{1}$ norm or average ionic concentration) initial conditions, the reduced time-independent system for the stationary states does give rise to smooth multiple solutions. However, at the parameter regime where multiple stationary solutions arise, the PNP-steric equation is found to be ill-posed. The finding that the PNP-steric system is ill-posed for sufficiently concentrated solutions raises a need to develop a model that is valid in this regime. In the second part of this study, we introduce a novel PNP-Cahn-Hilliard (PNP-CH) model which is well-posed also at high bulk concentrations. We show that this model gives rise to multiple stationary solutions that are stable with respect to PNP-CH dynamics. These stationary solutions differ by their bulk structure. Utilizing bifurcation analysis and a numerical continuation study, we map the various branches of stationary solutions.

The paper is organized as follows. In Section 2 we briefly review the PNP-steric model derivation. In Section 3 we show that under certain conditions, the corresponding free energy of the PNP-steric model is convex, and review the theory of generalized PNP models with convex energy. In particular, these models give rise to a unique stationary solution. Section 4 focuses on the case where the corresponding free energy of the PNP-steric system is non-convex. In contrast to the previous works, e.g. [30], we do not focus on the study of the solutions profiles $c_{i}(\phi)$, but rather present a different approach, and consider the solutions trajectories in the $(c_{1},c_{2})$ plane. We show that the $(c_{1},c_{2})$ plane can be divided into two regions, $D^{\rm convex}$ in which the corresponding free energy of the PNP-steric system is locally convex and $D^{\rm concave}$ in which the free energy is locally concave in $(c_{1},c_{2})$. We then present a classification of trajectories based on their path within these two regions. Solution trajectories which are fully contained in $D^{\rm convex}$ (aka type I trajectories) are studied in Section 4.1, and are shown to reflect the case where steric interactions are of perturbative nature, namely, they do not change qualitative properties of the stationary solutions . In contrast, solutions which at the bulk reside in $D^{\rm concave}$ reflect a case where steric effects change the qualitative behavior of the equation and give rise to structure formation in the bulk (aka type III solutions). In Section 4.2 we characterize periodic solutions of the PNP-steric system in this case, and show that they co-exist with the homogenous stationary solution. The above, thus, presents a full characterization of multiple stationary solutions of the PNP-steric model via the study of the reduced time-independent equations. Next, we go beyond the study of the stationary equations and consider stability of the stationary solutions with respect to PNP-steric dynamics. Our results show the PNP-steric system becomes ill-posed for concentrated solutions, and in particular in the regimes where multiple stationary solutions arise. Section 4.3 presents an interim summary of the study of the PNP-steric model, as well as a discussion regarding the reasons for the failure of the PNP-steric system to describe concentrated solutions: Intuitively, in classical PNP theory, it is the entropy that regularizes the solution and drives the system toward a spatially homogenous state. At high bulk concentrations, we show, however, that steric effects dominate entropic effects. Consequently, entropy fails to regularize the solution and the system becomes ill-posed as it is pushed toward higher and higher frequencies. To overcome this problem, in Section 5 we introduce a gradient energy term which regularizes the solution by direct energetic penalty of large gradients, giving rise to the PNP-Cahn-Hilliard (PNP-CH) system. Such gradient energy terms arise naturally in Cahn-Hilliard (Ginzburg-Landau) mixture theory [34], and has been recently used to model steric interactions in ionic liquids [35]. Section 5.1 presents the derivation of the PNP-CH system. Sections 5.2 and 5.3 present the linear and weakly nonlinear stability analysis of the homogenous state in an infinite domain, and characterizes the regimes in which super-critical and sub-critical bifurcations arise. Sections 5.4 and 5.5 consider the case of finite domain with possible applied voltage using numerical continuation techniques. The results present a wide family of stable stationary solutions of the PNP-CH system, showing a that finite domain has a significant effect on the solutions and their stability. Numerical methods are discussed in Section 6, and concluding remarks are presented in Section 7.

2 The PNP-steric Model

We briefly recall the PNP-steric model [26, 13]. For simplicity, we consider only the case of an electrolyte with two ionic species $c_{1},c_{2}$ with valences $z_{1},z_{2}$, respectively, such that

[TABLE]

The solution is bounded between two electrodes which are situated at $x=-L$ and $x=L$, with surface charge density such that the electric potential $\phi$ on the walls equals,

[TABLE]

where the reference potential is taken to be zero voltage at the bulk111Namely, at $L\to\infty$, the average electric potential at the interior of the domain is zero. This implies that the homogenous state $(\bar{c}_{1},\bar{c}_{2})$ is a stationary solution of the PNP-steric system when $\phi_{\pm L}=0$. . The system is globally electroneutral,

[TABLE]

where $\bar{c}_{i}$ is the average density of the $i^{\rm th}$ ionic species

[TABLE]

and where $\rho_{0}$ is a possible permanent ‘background’ charge density. The solvent is described as a dielectric medium with a uniform dielectric response $\epsilon$.

The action potential describing the system is given by

[TABLE]

where ${\bf c}=(c_{1},c_{2})$, ${\bf z}=(z_{1},z_{2})$, ${\bf\bar{c}}=(\bar{c}_{1},\bar{c}_{2})$, $q$ is the unit of electrostatic charge, $k_{B}$ is Boltzmann’s constant, and $T$ is the temperature. The action potential consists of the standard entropic and electrostatic free energy together with a correction functional ${\bf c}^{T}{\bf G}{\bf c}$, where ${\bf G}=[g_{ij}]$ is a symmetric $2\times 2$ matrix with non-negative entries, which accounts for steric interactions between ions [26, 13]. Particularly, when ${\bf G}\equiv{\bf 0}$, steric interactions are neglected and the action potential $\mathcal{A}$ reduces to the action potential corresponding to the Poisson-Nernst-Planck model.

Requiring that $\phi$ is a critical point of the action $\mathcal{A}$ yields Poisson’s equation

[TABLE]

Introduction of the non-dimensional variables

[TABLE]

where $\lambda_{D}$ is the Debye length

[TABLE]

and $c^{o}$ is an arbitrary reference density, see remark 1 below, yields the PNP-steric equation (presented after omitting the tildes)

[TABLE]

In particular, the stationary charge density profiles satisfy the Poisson-Boltzmann-steric (PB-steric) equations, composed of generalized Boltzmann equations

[TABLE]

Remark 1**.**

Note the explicit dependence of $\lambda_{i}$ in ${\bf\bar{c}}$, see (11a). Typically, this explicit dependence is eliminated by setting the reference density $c^{o}$, see (9), to be the bulk density, e.g., $\bar{c}=\bar{c}_{1}$. We retain an arbitrary reference density, $c^{o}$, in order to preserve the central role played by the connection between the bulk densities $\bar{c}_{i}$ and the Lagrange multipliers.

3 The case of convex energy and the functions $c_{i}(\phi)$

The free energy is attained by substituting Poisson’s equation (8a) into the action potential (7), and rescaling by (9), yielding

[TABLE]

where

[TABLE]

In [27], we have studied the case of a general function $h$, aka the family of generalized PNP models with corresponding free energy (12). In particular, we focused on the case when the function $h$ is strictly convex. As will be shown, in this case the corresponding generalized PNP model cannot give rise to multiple stationary solutions. Nevertheless, for completeness and to allow comparison, we briefly review this case of a convex free energy.

By Lemma [27, Lemma 1], for any strictly convex function $h$, the generalized Boltzmann system

[TABLE]

uniquely defines the functions $\{c_{1}(\phi,{\bf\bar{c}}),c_{2}(\phi,{\bf\bar{c}})\}$.

[TABLE]

where the dependence of ${\bf c}=(c_{1},c_{2})$ upon $\phi$ depends on $h(c_{1},c_{2})$ via the corresponding generalized Boltzmann system. The system (14), together with (11a), uniquely defines a stationary solution $(c_{1}(x),c_{2}(x),\phi(x))$ of the corresponding generalized PNP equations [27]. Furthermore, in [27], we have shown that, for all strictly convex $h$, the dynamical systems (14) are locally homeomorphic. In this sense, all generalized Poisson-Boltzmann systems with strictly convex $h$, “behave the same”. In particular, at the limit of zero applied voltage $\phi_{\pm L}\to 0$, the solutions of all such generalized Poisson-Boltzmann systems, subject to the constraint (6), tend to the homogenous bulk state ${\bf c}\equiv(\bar{c}_{1},\bar{c}_{2})$ [27, Section 3.1]. Trivially, these stationary solutions are minimizers of the energy $\mathcal{E}$, and therefore are stable under Nernst-Planck dynamics which are (Wasserstein) gradient flows along $\mathcal{E}$.

In the case of the PNP-steric system, see (13),

[TABLE]

Lemma 1**.**

The PNP-steric system (10) has a corresponding strictly convex free energy, $h$, for all ${\bf c}>0$ if and only if ${\bf G}$ is a positive semi-definite matrix.

Proof.

Note that the Lemma also applies to the case of multi-species $N\geq 2$.

If $G$ is a positive semi-definite matrix (p.s.d.), then by (15), $\mbox{Hessian}_{{\bf c}}(h)$ is a positive definite matrix, hence $h$ is strictly convex.

Assume that $G$ is not p.s.d., then ${\bf G}$ has at least one negative eigenvalue $\lambda_{-}<0$ with corresponding eigenvector $v_{-}$, i.e.,

[TABLE]

Let ${\bf c}_{-}=-(1/\lambda_{-},\cdots,1/\lambda_{-})$. Then,

[TABLE]

Hence, ${\bf I}+\mbox{diag}({\bf c}_{-}){\bf G}$ is singular, where $v_{-}\in\ker\{\mbox{Hessian}_{{\bf c}}(h)\}$, implying that $h$ is not strictly convex for all ${\bf c}>0$, see (15). ∎

Figure 2A presents an example of the stationary profiles of the PNP-steric (10) in a case where ${\bf G}$ is positive-definite. It can be seen that in the interior of the domain, and far from the boundaries, the solution reaches a homogenous state. The profiles $c_{i}(\phi)$, corresponding to the solution of Figure 2A, are presented in Figure 2C. As expected, see [27, Lemma 1], these profiles define unique functions of $\phi$ which satisfy (11a). The focus of this work is on multiple stationary solutions. Therefore, in what follows, we will consider the PNP-steric model when $h$ is not convex.

4 Non-convex energetic landscape of the PNP-steric system

Let us consider stationary solutions of the PNP-steric system (10) with counter-charged two species in the non-convex case, i.e., the case $|{\bf G}|<0$, see Lemma 1. In what follows, it becomes useful to consider the differential form of the steric Boltzmann equations (11a). By differentiation of (11a) by $x$ and formal isolation444We consider below the cases where the system is not invertible, i.e., $D=0$, see, e.g., Lemma 5. of $c_{1,x}$ and $c_{2,x}$, the PB-steric system (11) reads as

[TABLE]

In contrast to the convex case, the solutions profiles ${\bf c}(\phi)=\{c_{1}(\phi,{\bf\bar{c}}),c_{2}(\phi,{\bf\bar{c}})\}$ may not be functions of $\phi$. Therefore, rather than focusing on ${\bf c}(\phi)$, we take a different approach, and consider the solutions trajectories in the $(c_{1},c_{2})$ plane. By equations (16b), a solution trajectory passing through a point $(c_{1}^{0},c_{2}^{0})$ satisfies

[TABLE]

The following Lemma characterizes the solution trajectories.

Lemma 2**.**

Let ${\bf c}_{0}=(c_{1}^{0},c_{2}^{0})\in R_{+}^{2}$m and let $z_{1},z_{2}\in R$ that satisfy (3). Then,

There exists a unique solution trajectory $\Gamma_{{\bf c}_{0}}$, defined by (17), that passes through ${\bf c}_{0}$, and which remains in the quarter plane $R_{+}^{2}$. 2. 2.

The solution trajectory $\Gamma_{{\bf c}_{0}}$ intersects with the line $\{{\bf c}|E_{x}=0\}$ at a single point. 3. 3.

Let $D({\bf c}_{0})$ be defined by (16c). If $D({\bf c}_{0})<0$, then the solution trajectory intersects with the curve $\{{\bf c}|D=0\}$ at least at two points.

Proof.

Since $z_{1}>0>z_{2}$, see (3), the right-hand-side of (17) is strictly negative. Therefore, $c_{1}\leq c_{1}^{0}$ for $c_{2}\geq c_{2}^{0}$, implying that

[TABLE]

By (18) and the fundamental theory of ODE, Equation (17) admits a smooth monotonically decreasing solution for any initial condition $c_{1}(c_{2}^{0})=c_{1}^{0}$. Further by Gronwell’s Lemma, for $c_{2}\geq c_{2}^{0}$, the solution is bounded by

[TABLE]

The solution is monotonically decreasing, hence $c_{1}(c_{2})\geq c_{1}^{0}$ for $c_{2}<c_{2}^{0}$. Therefore, $c_{1}(c_{2})>0$ for all $c_{2}>0$. 2. 2.

The curve $\Gamma_{E}(c_{1})$ along which $E_{x}=0$ is given by the strictly monotonically increasing function, see (16a),

[TABLE]

Thus, $c_{1}(c_{2})-\Gamma_{E}(c_{2})$ is a strictly monotonically decreasing function. By (19) and since $m,M<0$,

[TABLE]

implying that

[TABLE]

Applying the same arguments to the inverse function $c_{2}(c_{1})$ yields that $\lim_{c_{2}\to 0^{+}}c_{1}(c_{2})=\infty$.

[TABLE]

Therefore, the solution trajectory $\Gamma_{{\bf c}_{0}}$ intersects with the line $E_{x}=0$ at a single point. 3. 3.

By (20), $\lim_{c_{2}\to\infty}D(c_{1}(c_{2}),c_{2})=\infty$. Therefore, continuity of $D(c_{1},c_{2})$ implies that the solution trajectory $\Gamma_{{\bf c}_{0}}$ intersects with the curve $\{{\bf c}|D=0\}$ at a point on $\Gamma_{{\bf c}_{0}}$ for which $c_{2}^{0}<c_{2}<\infty$. Similarly, the limit $\lim_{c_{2}\to 0^{+}}c_{\Gamma}(c_{2})=\infty$ implies that the solution trajectory $\Gamma_{{\bf c}_{0}}$ intersects with the curve $D=0$ at a point on $\Gamma_{{\bf c}_{0}}$ for which $0<c_{2}<c_{2}^{0}$.

∎

Figure 2 presents numerous solution trajectories in the $(c_{1},c_{2})$ plane corresponding to different solutions of the same PB-steric system (11). As follows from Lemma 2, the trajectories are monotonically decreasing functions of $c_{2}$ that intersect once with the line $E_{x}=0$. The significance of the intersection with the line $E_{x}=0$ stems from the fact that the solution is locally electro-neutral at this point, and therefore the solution is expected to reach this point in the limit of zero applied voltage or sufficiently far from the charged boundaries.

The shaded region presented in Figure 2 is the region in which $D(c_{1},c_{2})<0$. As will be shown, the sign of $D$ along the solution trajectory dictates key properties of the solutions. Indeed, differentiating (16a) by $x$ and substituting (16b) in the resulting equation yields

[TABLE]

Classical ODE theory implies that the nature of solutions of this second-order linear ODE - oscillatory or not - depends on the sign of the coefficient of $E$. Particularly, the sign of $D(x)$ dictates the sign of the coefficient of $E$ since the coefficient $a(x)$ in equation (21) is strictly positive for all $x$.

Let us distinguish between types of solution trajectories by the sign of $D$ along the trajectory, and at the intersection point with the line $E_{x}=0$:

• •

Type I trajectories correspond to trajectories that remain in the region $\{{\bf c}|D>0\}$, i.e., they do not intersect with the curve $D=0$.

• •

Type III trajectories correspond to trajectories for which $D<0$ in a surrounding of the intersection point with the line $E_{x}=0$, and thus intersect with the curve $D=0$.

• •

Other trajectories are denoted type II trajectories.

Figure 2 presents an example of each of the trajectory type.

4.1 Type I solution trajectories

In this section, we characterize type I solutions, and show that they qualitatively behave as solutions of generalized PNP systems in the case of convex $h$. Particularly, all solutions of the PNP-steric system (10) in the convex case, $|{\bf G}|>0$, are type I solutions, since $D>0$ for all $(c_{1},c_{2})$. Similarly, in the dilute case, $c_{i}(x)\ll 1$, even if $h$ in non-convex, the solution trajectories are type I trajectories. This can be expected since in the dilute case the PB-steric system (11) reduces at leading order to the Poisson-Boltzmann system, i.e., to the system (11) with ${\bf G}\equiv 0$. Lemma 1 implies that the concentrations of solutions corresponding to type I trajectories are bounded by $(c_{1},c_{2})\in[0,-1/\lambda_{-}]\times[0,-1/\lambda_{-}]$ where $\lambda_{-}$ is the negative eigenvalue of ${\bf G}$. The following Lemma quantifies a minimal regime of solution concentrations for which solutions correspond to type I solutions.

Lemma 3**.**

Let ${\bf c}_{0}=(c_{1}^{0},c_{2}^{0})\in(0,c_{1}^{\rm bound})\times(0,c_{2}^{\rm bound})$ where

[TABLE]

and let $g_{11},g_{22},g_{12}>0$ such that $g_{11}g_{22}-g_{12}^{2}<0$. Then, the solution trajectory given by (17) satisfies for all $x$,

[TABLE]

Proof.

It is easy to verify from that if $c_{1}<c_{1}^{\rm bound}$ or $c_{2}<c_{2}^{\rm bound}$ then $D>0$, see (16c). By Lemma 2, $dc_{1}/dc_{2}<0$. Thus, if $c_{2}>c_{2}^{0}$ then $c_{1}<c_{1}^{0}<c_{1}^{\rm bound}$ and $D>0$. Otherwise, $D>0$ since $c_{2}\leq c_{2}^{0}<c_{2}^{\rm bound}$. ∎

Figure 2A presents an example of the concentration profiles of a type I solution. The solution trajectory, corresponding to the solution of Figure 2A is presented in Figure 2B (thick solid curve). In addition, Figure 2B presents the trajectory $\Gamma_{\bf\bar{c}}$ as defined by (17) for any $c_{2}$ (dashed curve). As expected, the trajectory $\Gamma_{\bf\bar{c}}$ remains in the region $D>0$ for all $c_{2}$ implying that $D(c_{1}(x),c_{2}(x))>0$ for any applied potential $V=\phi_{L}-\phi_{-L}$, $\phi_{L}\phi_{-L}<0$.

Intuitively, type I solutions remain in a region in which the free energy is locally convex, are therefore they are expected to behave as solutions in the convex case. For example, as in the convex case, the generalized Boltzmann system (11a) defines unique functions of $\phi$, see, e.g., Figure 2C presenting the profiles $c_{i}(\phi)$ which correspond to the solution of Figure 2A. In addition, as in the convex case, at the limit of zero applied voltage $\phi_{\pm L}\to 0$, type I solutions, subject to the constraint (6), tend to the homogenous bulk state ${\bf c}\equiv(\bar{c}_{1},\bar{c}_{2})$.

Lemma 4**.**

Let ${\bf c}_{0}\equiv(\bar{c}_{1},\bar{c}_{2})$ be a homogenous state satisfying local electroneutrality

[TABLE]

and

[TABLE]

Then, ${\bf c}_{0}$ is linearly stable under steric-Poisson-Nernst-Planck dynamics (10).

Proof.

Consider a perturbation of the homogenous state ${\bf c}_{0}$

[TABLE]

that satisfies the constraint (6). Then, $\mathbf{\lambda}$ satisfies

[TABLE]

By (16c) and (23), $\mbox{Hessian}_{{\bf c}}(h)$ is a positive definite matrix, while ${\bf z}{\bf z}^{T}$ is a positive semi-definite matrix. Therefore, for $k\neq 0$, ${\bf A}$ is positive semi-definite. This implies that $\lambda(k)<0$ for $k\neq 0$. When $k=0$, the above system has one eigenvector ${\bm{\delta}}_{1}=\mbox{diag}({\bf\bar{c}}){\bf z}$ with corresponding eigenvalue

[TABLE]

and one eigenvector ${\bm{\delta}}_{2}=(-z_{2},z_{1})^{T}$ with corresponding zero eigenvalue $\lambda=0$. The latter eigenvalue reflects changes in bulk concentration that maintain local electro-neutrality (22), but violates (6). Therefore, ${\bm{\delta}}_{2}\perp{\bm{\delta}}$. ∎

4.2 Type III solution trajectories

In what follows, we will focus on the case of a concentrated solution such that

[TABLE]

where ${\bar{c}}_{i}$ is the average concentration of each ionic species, see (6).

We first show that in this case the PNP-steric system (10) in free space ($L\to\infty$) admits multiple stationary solutions. To do so, it is useful to consider the related initial value problem (16) for $x>0$ with initial conditions

[TABLE]

rather than the boundary value problem (16,11c). Solutions of the initial value problem (16,26) for $x\geq 0$, extend to the whole domain via the symmetries $E(x)\to-E(-x)$ and $c_{i}(x)\to c_{i}(-x)$. When $c_{i}^{0}=\bar{c}_{i}$, $i=1,2$, the initial value problem (16,26) has the trivial solution $c_{i}(x)\equiv\bar{c}_{i}$, $i=1,2$. We now show that this system also admits periodic solutions:

Theorem 1** (Existence of periodic solutions).**

Let ${\bar{\bf c}}=(\bar{c}_{1},\bar{c}_{2})\in R_{+}^{2}$ such that

[TABLE]

Then, the Poisson-Boltzmann-steric initial value problem (16,26,3) admits smooth periodic solutions.

Proof.

Lemma 2 ensures that there exists a unique solution trajectory $\Gamma_{\bar{\bf c}}=(c_{1}(c_{2}),c_{2})$ that passes through $\bar{\bf c}$. By (27), and due to the continuity of $D$ in $R_{+}^{2}$, there exist a surrounding of $\bar{\bf c}$ in which $\Gamma_{\bar{\bf c}}$ does not cross $D=0$, i.e., there exists $\varepsilon>0$ such that

[TABLE]

Let $\bar{c}_{2}<c_{2}^{0}<\bar{c}_{2}+\varepsilon$. Since $dc_{1}/dc_{2}<0$, see Lemma 2, $c_{1}^{0}:=c_{1}(c_{2}^{0})$ satisfies

[TABLE]

Denote by $(c_{1}^{A},c_{2}^{A},E^{A})$ the solution of (16,26). Since $(c_{1}^{0},c_{2}^{0})\in\Gamma_{\bar{\bf c}}$, the solution trajectory corresponding to the solution $(c_{1}^{A},c_{2}^{A},E^{A})$ is $\Gamma_{\bar{\bf c}}$, see Lemma 2. We first show that there exists a finite interval $[0,x_{A}]$ where $c_{2}^{A}(x)$ monotonically decreases from $c_{2}^{A}(0)=c_{2}^{0}$ to $c_{2}^{A}(x_{A})=\bar{c}_{2}$, see Figure 4 for illustration. Indeed, by (29), $E_{x}^{A}(0)>0$, thus there exists a maximal surrounding $0<x<x_{A}$, possibly $x_{A}=\infty$, such that such tha

[TABLE]

and $E_{x}^{A}(x_{A})=0$. Since the solution trajectory $\Gamma_{\bar{\bf c}}$ intersects with the line $E_{x}=0$ at a single point, see Lemma 2, then by construction $c_{i}^{A}(x_{A})=\bar{c}_{i}$, $i=1,2$, and hence

[TABLE]

Furthermore, since $D(c_{1}^{0},c_{2}^{0})<0$ and due to continuity of $D$ and of the solution $(c_{1}^{A},c_{2}^{A},E^{A})$, there exists a maximal interval $0<x<x_{\delta}\leq x_{A}$ such that

[TABLE]

and either $D(c_{1}^{A}(x_{\delta}),c_{2}^{A}(x_{\delta}))=0$ or $x_{\delta}=x_{A}$. In this surrounding, $E^{A}(x)>0$, thus by (3),

[TABLE]

see (16b), hence

[TABLE]

The inequalities (28) imply that $D<0$ for all $0<x<x_{\delta}$. Hence $x_{\delta}=x_{A}$. Finally, the inequalities (28,30,31) imply that for any $0<\varepsilon<x_{A}$,

[TABLE]

see (16). This implies that $c_{2}^{A}(x)$ monotonically decreases from $c_{2}^{A}(0)=c_{2}^{0}$ to $c_{2}^{A}(x_{A})=\bar{c}_{2}$ at finite time $x=x_{A}<(\bar{c}_{2}-c_{2}^{0})/M$.

At the point $x=x_{A}$, $E^{A}_{x}(x_{A})=0$, hence $E^{A}(x)$ reaches a local maximum point $E^{A}(x_{A})=E^{A}_{\max}$. The solutions of (16,26) with the initial condition $(c_{1}^{0},c_{2}^{0})\in\Gamma_{\bar{\bf}c}$ are parametrized by a single parameter $c_{2}^{0}$. When $c_{2}^{0}=\bar{c}_{2}$, then $c_{2}\equiv\bar{c}_{2}$ and $E^{A}_{\max}=0$. Continuity of the ODE solution with respect to the initial conditions implies that for any $0\leq E_{\max}\leq E^{A}_{\max}$, there exists a solution of (16,26) with (29) satisfying

[TABLE]

Similarly, there exists a solution of $(c_{1}^{B},c_{2}^{B},E^{B})$ of (16,26) with

[TABLE]

such that at a point $x=x_{B}<\infty$

[TABLE]

and for any $0\leq E_{\max}\leq E^{A}_{\max}$, there exists a solution $(\tilde{c}_{1}^{B},\tilde{c}_{2}^{B},\tilde{E}^{B})$ $(c_{1}^{B},c_{2}^{B},E^{B})$ of (16,26) such that

[TABLE]

Let $E_{\max}=\min\{E^{A}_{\max},E^{B}_{\max}\}$, and let $\tilde{E}^{A}$ and $\tilde{E}^{B}$ be solutions of (16,26) satisfying (32) and (33), respectively. Then,

[TABLE]

is a smooth solution of (16,26) in the interval $[0,L]$ where $L=x_{A}+x_{B}$, see illustration in Figure 4. This solution is periodically extended to the whole domain by exploiting the symmetry, $c_{i}(x_{0}-x)=c_{i}(x_{0}+x)$ where $E(x_{0})=0$,

[TABLE]

∎

Figure 2A presents a ‘periodic solution’ of the PB-steric system (16). We note that the solution is not necessarily presented in an interval in which it satisfies periodic boundary conditions. Rather, in what follows, we refer to ‘periodic solutions’ as solution for which there exist an interval in which they are periodic. As follows from Theorem 1, the corresponding solution trajectory lies in the domain $\{{\bf c}|D({\bf c})<0\}$, and intersects with the line $E_{x}=0$, see Figure 2B. The extremum points of the solutions correspond to the end points of the solution trajectory, see $\circ$ and $*$ markers in Figure 2A and 2B. The corresponding solution profiles $c_{i}(\phi)$ have a regions of triple solutions, whereas the branch corresponding to the periodic solutions is bounded between two folds, see Figure 2C. Particularly, the folds in $c_{i}(\phi)$ correspond to points where the solution trajectory $\Gamma_{\bar{\bf c}}$ intersects with the curve $D=0$.

Theorem 1 proves the existence of periodic solutions of the initial value problem (16b,21,26). These solutions maintain $D<0$ for all $x$. We now show, however, that for concentrated solutions satisfying (25), stationary solutions of the PNP-steric equation with large enough applied voltage $|\phi_{L}-\phi_{-L}|$ can only be described by solution which cross the curve $\{{\bf c}|D=0\}$. Indeed, near a sufficiently charged wall, $D>0$: Isolating $\phi$ in the equation for $c_{1}$ in (11a) yields

[TABLE]

Applying the method of dominant balance for $\phi\gg 1$, yields that since $z_{1}>0$ and $c_{i}>0$, the only way to satisfy the above relation is if $z_{1}\phi\approx-\log c_{1}$. The latter implies that $c_{1}\ll 1$, hence $D(c_{1},c_{2})>0$. Similarly, $D(c_{1},c_{2})>0$ when $-\phi\gg 1$ since $c_{2}\ll 1$.

Lemma 5**.**

Under the conditions of Lemma 1, there exists a smooth stationary solution $(c_{1}(x),c_{2}(x),\phi(x))$ of (10) that crosses $D=0$, i.e., $D(c_{1}(x),c_{2}(x))=0$ at a point $x\in(-L,L)$.

Proof.

Let us consider a point ${\bf c}^{*}=(c_{1}^{*},c_{2}^{*})$ on the curve $\{{\bf c}|D=0\}$ which satisfies $c_{1}^{*}>\bar{c}_{1}$ and $c_{2}^{*}<\bar{c}_{2}$. Consider the PB-steric initial value problem (16b,26) with $(c_{1}^{0},c_{2}^{0})=(c_{1}^{*},c_{2}^{*})$. Then,

[TABLE]

where

[TABLE]

Let us choose the intersection point $(c_{1}^{*},c_{2}^{*})$ such that $c_{1}^{*}$ is sufficiently large so that $f(c_{1}^{*},c_{2}^{*})>0$. In this case,

[TABLE]

This implies that $c_{1,x}(c_{1}^{*},c_{2}^{*})$ and $c_{2,x}(c_{1}^{*},c_{2}^{*})$ exists, and thus Equation (16b) has a smooth solution at a surrounding $[-\varepsilon,\varepsilon]$ of $x=0$ which crosses $D=0$. ∎

Lemmas 1 and 5 prove the existence of a family of stationary solutions of the PNP-steric equation (10). Let us now discuss the stability of these solutions. Clearly, solutions that cross $\{{\bf c}|D=0\}$ are unstable. Indeed, let us consider a smooth solution ${\bf u}$ of the initial value problem (16b,21,26) that crosses the curve $\{{\bf c}|D=0\}$ at point $(c_{1}^{*},c_{2}^{*})$, and let us consider a solution ${\bf u}_{\delta}$ of the initial value problem (16b,21,26) with the perturbed initial condition $(c_{1}^{0},c_{2}^{0})=(c_{1}^{*},c_{2}^{*})+\vec{\delta}$ such that $D(c_{1}^{0},c_{2}^{0})\neq 0$. Lemma 2 ensures that the trajectory associates with ${\bf u}_{\delta}$ will intersect with the curve $\{{\bf c}|D({\bf c})=0\}$. However, since $E(x=0)=0$ and $E_{x}(x=0)\neq 0$, then for small enough $\vec{\delta}$, at the intersection point of the solution trajectory with the curve $\{{\bf c}|D=0\}$, the solution will satisfy $E\neq 0$. Thus, as the solution approaches the intersection point where $D=0$, its gradient blows up.

The sensitivity of the above solution to changes in the potential $\phi$ is also evident when considering the derivatives of the concentration profiles with respect to $\phi$

[TABLE]

see (16b). Thus,

[TABLE]

Remark 2**.**

Lemma 5 and the discussion below also applies to type II solutions. These correspond to trajectories for which $D>0$ in a surrounding of the intersection point with the line $E_{x}=0$, but do intersect with the curve $D=0$. Therefore, while for low applied voltage type II solutions resemble type I solutions, at higher applied voltages the solution trajectory will cross the curve $D=0$. At which case, the corresponding smooth solution, if exists, will be unstable.

Finally, we consider the stability of homogenous stationary solutions $c_{i}(x)\equiv\bar{c}_{i}$ of the PNP-steric system (10).

Lemma 6**.**

Let ${\bf c}_{0}\equiv(\bar{c}_{1},\bar{c}_{2})$ be a homogenous state satisfying local electroneutrality (22)

[TABLE]

and

[TABLE]

Then, ${\bf c}_{0}$ is linearly unstable under Poisson-Nernst-Planck dynamics (10).

Proof.

Consider a perturbation of the homogenous state ${\bf c}_{0}$

[TABLE]

that satisfies the constraint (6). Following the proof of Lemma 4, $\mathbf{\lambda}$ satisfies (24). In particular, for $k\gg 1$,

[TABLE]

Since $D({\bf c}_{0})=|\mbox{Hessian}_{{\bf c}_{0}}(h)|<0$ and since ${\rm tr}(\mbox{Hessian}_{{\bf c}_{0}}(h))>0$, the Hessian matrix has a negative eigenvalue $\lambda_{-}<0$. Thus,

[TABLE]

see Figure 7, implying that the problem is linearly unstable.

∎

Lemma 6 can also be applied to the periodic solutions described in Theorem 1, in the case they have small enough amplitude $|c_{2}^{0}-\bar{c}_{2}|\ll 1$. We, further verify by direct integration that periodic solutions on finite domain as well as periodic solutions with larger amplitudes are also unstable (data not shown).

4.3 Interim summary - failure of the PNP-steric at high bulk concentrations

The study presented in Section 4 considers solutions of the PNP-steric system (10) in two distinct cases. In the dilute case, particularly when $D(\bar{c}_{1},\bar{c}_{2})>0$, the PNP-steric system admits stable unique solutions, aka, type I solutions. In the limit of zero applied voltage $\phi_{\pm L}\to 0$, these solutions tend to the homogenous bulk state ${\bf c}\equiv(\bar{c}_{1},\bar{c}_{2})$, in qualitative agreement with solutions of the PNP equation and the wide class of generalized PNP equations with a convex free energy. In the concentrated solution case, $D(\bar{c}_{1},\bar{c}_{2})<0$, the PNP-steric system (10) may admit multiple stationary type III solutions, e.g., the homogenous bulk state ${\bf c}\equiv(\bar{c}_{1},\bar{c}_{2})$ and periodic solutions. These solutions, however, are all unstable as higher and higher Fourier-modes become dominant. 555Type II solutions are a borderline case between type I and type III, they are well-defined and stable at low applied voltages, but not at higher ones.

To understand why the PNP-steric system (10) becomes unstable and even ill-posed at high bulk concentrations let us first consider the classical PNP model (1) with $g_{ij}=0$. Intuitively, in classical PNP theory, the entropy regularizes the solution and drives the system toward a spatially homogenous state. Indeed, in the limit of zero charge density ${\bf z}\to 0$, the PNP system (1) with $g_{ij}=0$, reduces to the heat equation and thus tends toward the spatially homogenous solution ${\bf c}\equiv\bar{\bf c}$. Steric effects, however, compete and may even dominate entropic effects. Consider, for example, the PNP-steric model in a domain $[-1,1]$ with $g_{11}=g_{22}=0$, then the free energy takes the form, see (12),

[TABLE]

where

[TABLE]

In the case of the spatially homogenous solution ${\bf c}_{\rm uniform}\equiv\bar{\bf c}$, the entropic and electrostatic contributions of the energy both vanish, while $\mathcal{E}_{\rm steric}=4g_{12}\bar{c}^{2}$. Overall, the PNP-steric free energy of the homogenous solution scales as

[TABLE]

In contrast, the ‘strongly segregated’ pattern ${\bf c}_{\rm segregated}$ with segregation pattern frequency $n$

[TABLE]

yields $\mathcal{E}_{\rm steric}=0$ since $c_{1}(x)c_{2}(x)=0$. Furthermore, direct integration shows that

[TABLE]

so that the overall free energy of the segregated pattern (35) is

[TABLE]

Clearly, the segregated pattern is energetically preferential over the homogenous solution at high concentrations, $\bar{c}\gg\alpha/g_{12}$. Significantly, electrostatic contributions to the energy as well as the overall free energy decrease with increasing pattern frequencies. The above example introduces functions with discontinuities. It is possible to find a smooth solution which is arbitrarily close (e.g., in $L^{2}$ norm) to ${\bf c}_{\rm segregated}$. This solution would clearly have large gradients at the discontinuity points $\frac{k}{2n}$ of ${\bf c}_{\rm segregated}$. Nevertheless, since $\mathcal{E}$ does not impose any energetic penalty for large gradients, one may expect that $\mathcal{E}[{\bf c}_{\rm smooth}]\approx\mathcal{E}[{\bf c}_{\rm segregated}]$. The emerging picture is that at high bulk concentrations, steric effects dominate entropic effects. Consequently, entropy fails to regularize the solution and the system becomes ill-posed as it is pushed toward higher and higher frequencies. This nature of the PNP-steric system becomes apparent in the proof of Lemma 6.

5 The PNP-Cahn-Hilliard (PNP-CH) model

The failure on the PNP-steric (10) model in the case of high bulk concentrations where entropic contributions are insufficient to efficiently regularize the solution suggests that the PNP-steric model is missing terms that become dominant in the presence of large gradients. In this section, we introduce high-order terms that directly regularize the solution giving rise to the PNP-Cahn-Hilliard model.

5.1 Model derivation

Let us consider the action potential

[TABLE]

where $\mathcal{A}_{\rm PNP-steric}$ is the action potential (7) of the PNP-steric model (10), and where $\sigma$ is a gradient energy coefficient. The gradient energy term $\sum_{i=1}^{N}|\nabla c_{i}|^{2}$ introduces a direct energetic penalty to large gradients, and thus regularizes the solution. This term arises naturally in Cahn-Hilliard (Ginzburg-Landau) mixture theory [34], and has been recently used to model steric interactions in ionic liquids [35]. Moreover, while the steric interaction term ${\bf c}^{T}{\bf G}{\bf c}$, see (7), is derived by a leading order local approximation of the Lennard-Jones interaction potential [26, 13], the high-order gradient energy term seems to be tightly related to high-order terms in the expansion of the Lennard-Jones potential [36]. Since the steric and the gradient energy terms appears in the Cahn-Hilliard mixture theory, and since, as will be shown, the resulting system behaves more like a non-local Cahn-Hilliard (Ohta-Kawasaki) system than a PNP system, we will denote the resulting system as the Poisson-Nernst-Planck-Cahn-Hilliard (PNP-CH) model.

Similar to the PNP-steric, requiring that $\phi$ is a critical point of the action yields Poisson’s equation (8a), whereas the Wasserstein gradient flow along the action gives rise to generalized Nernst-Planck equations for the dynamics of charge concentration

[TABLE]

Introduction of the non-dimensional variables (9) and

[TABLE]

yields the PNP-CH equation (presented after omitting the tildes)

[TABLE]

We adopt the no-flux boundary conditions for the species concentration and the Dirichlet boundary conditions for the potential of the PNP-steric equation, see (10b). The high-order terms, however, make it necessary to introduce additional boundary conditions. The newly introduced gradient energy term is expected to become dominant only when the PNP-steric system fails, i.e., when $D<0$, see (16c). Since near a sufficiently charged wall, $D>0$ (see Section 4), it is reasonable to set

[TABLE]

For this choice, the stationary profiles, given by

[TABLE]

where $\lambda_{i}$ are associated with mass conservation, are independent of $\sigma$ at the boundaries.

The arising system (37,10b,38) accounts for the various interactions between the model components in a self-consistent way, while satisfying the second law of Thermodynamics, via the Onsager relation,

[TABLE]

5.2 Linear stability analysis of the homogenous state

Recall that in the absence of applied voltage, $\phi_{\pm L}=0$, the homogenous stationary solution ${\bf c}\equiv\bar{\bf c}$ is stable under dynamics of all generalized PNP models with a convex free energy, see Section 3. In contrast, at large enough concentrations, the homogenous stationary solution is unstable under PNP-steric dynamics, see Lemma 6. Here, we study the stability of the homogenous stationary solution under dynamics of the PNP-CH system (37,10b,38) in free-space, i.e., $L\to\infty$ and $\phi_{\pm L}=0$.

Let us consider a perturbation of the homogenous state $c_{i}\equiv\bar{c}_{i},\,i=1,2$,

[TABLE]

Substitution of (39) into the PNP-CH system (37) yields the linearized system ${\bf M}v=\lambda v$, where

[TABLE]

and $v=[u_{1},v_{1}]^{T}$. The trace of the above matrix ${\bf M}(k)$ is negative for all $k\in R$. Therefore, at least one of its eigenvalues in strictly negative. Let us denote by $\lambda(k)$ the larger eigenvalue of ${\bf M}(k)$.

When $k=0$, the matrix ${\bf M}$ has one zero eigenvalue $\lambda(k=0)=0$ with corresponding eigenvector $\vec{v}_{0}:=(-z_{2},z_{1})^{T}$. This eigenvector corresponds to changes in bulk concentration that do not effect global charge neutrality, see (5). For $k\neq 0$, $\lambda(k)$ satisfies

[TABLE]

Varying $\sigma$, the instability onset is obtained by seeking $\sigma=\sigma_{c}$ and $k_{c}$ such that $\lambda(k=k_{c},\sigma=\sigma_{c})=0$, $\lambda(k\neq k_{c})<0$ and $\frac{d\lambda}{dk}(k_{c})=0$. At these conditions, $\sigma>\sigma_{c}$, implies that $\lambda(k)<0$ for all $k\neq 0$, hence the homogenous solution is stable. Substituting $\lambda=0$ in (41) yields that $\sigma_{c}$ satisfies

[TABLE]

Since all terms in the above equality are strictly positive, except possibly the fourth term, a necessary condition for instability about $(\bar{c}_{1},\bar{c}_{2})$ is that

[TABLE]

or

[TABLE]

This condition is nothing but the condition $D(\bar{c}_{1},\bar{c}_{2})<0$, see (16c). When (42) has a solution $k_{c},\sigma_{c}>0$, the solution of the linearized system ${\bf M}v=\lambda v$ at the bifurcation point $\sigma=\sigma_{c}$ takes the form

[TABLE]

where $\varepsilon$ is the magnitude of the $k=0$ and $k=k_{c}$ modes of the initial condition, $\alpha,\,\beta$ are constants, and $\vec{v}_{0}$, $\vec{v}_{k_{c}}$ are the normalized eigenvectors corresponding the zero eigenvalue of the matrices ${\bf M}(0)$ and ${\bf M}(k_{c})$, respectively,

[TABLE]

The expressions for $(k_{c},\sigma_{c})$ are analytically available but, in the general case, they are cumbersome. In the symmetric case, $g:=g_{11}=g_{22}$, $\bar{c}:=\bar{c}_{1}=\bar{c}_{2}$ and $z_{1}=-z_{2}=1$, the expressions for the instability onset reduce to

[TABLE]

5.3 Weakly nonlinear analysis

To understand the nature of instability at $\sigma=\sigma_{c}$, we conduct a weakly nonlinear analysis. For simplicity, we restrict to the case $\bar{c}:=\bar{c}_{1}=\bar{c}_{2}$. Let us consider the slow dynamics of a solution slightly below the instability onset $\sigma=\sigma_{c}-\varepsilon$.

[TABLE]

Motivated by (43), we consider the ansatz

[TABLE]

Note that we limit this study to spatially uniform amplitudes $\alpha_{0}$ and $\beta_{0}$. In this case, conservation of mass implies that

[TABLE]

Hence, $\alpha_{0}(T)\equiv 0$.666In the general case where $\alpha_{0}$ is spatially varying, the nonlinear interaction between the two modes, $k=0$ and $k=k_{c}$, specifically via the spatial derivatives of $\alpha_{0}$, must be taken into account, see, e.g., [37, 38, 39].

At $O(\varepsilon)$, the solution satisfies

[TABLE]

whose solution is of the form

[TABLE]

where $\vec{\gamma}$ satisfies

[TABLE]

At $O(\varepsilon^{3/2})$, the equations take the form

[TABLE]

Since ${\bf M}(k_{c})$ is a singular matrix with $\vec{v}_{0}\in\ker{\bf M}(k_{c})$, Equation (44) has a solution only when

[TABLE]

In particular, the stationary solution $\beta_{0}(T)\equiv\beta_{0}$ satisfies

[TABLE]

If $-\frac{<\vec{v}_{0},\vec{a}>}{<\vec{v}_{0},\vec{b}>}>0$, the bifurcation is super-critical leading to the solution

[TABLE]

otherwise the bifurcation is sub-critical. In the symmetric case, $g_{11}=g_{22}$ and $z_{1}=-z_{2}=1$, the analytic expression is available

[TABLE]

implying that bifurcations in the symmetric case are always super-critical.

Figure 8A presents the regions in parameter space in which the system undergoes supercritical and subcritical phase transitions. As expected, in the symmetric case $g_{11}-g_{22}=0$, phase transitions are always super-critical, while subcritical bifurcations emerge in the asymmetric case $g_{11}\neq g_{22}$, near the curve $g_{12}=g_{12}^{\rm crit}$. Below the curve $g_{12}<g_{12}^{\rm crit}$, the only stable solution is the homogenous state ${\bf c}\equiv\bar{\bf c}$, see Figure 8(1). The bifurcation diagram of the system at point ‘a’ residing in the supercritical regime is presented in Figure 8B. At $\sigma<\sigma_{c}$, the homogenous state becomes unstable, and rather a branch of stable periodic solutions whose amplitude is proportional to $\sqrt{\sigma_{c}-\sigma}$, see Figure 8(2), emerges. Finally, Figure 8C presents the bifurcation diagram of the system at point ‘b’ which resides in the subcritical regime. In this case, there exists a region $\sigma_{c}<\sigma<\sigma_{\rm fold}$, where $\sigma_{\rm fold}$ is the fold point, of multiple stationary which are all stable. For example, as verified by direct simulation, at point $\sigma=\sigma^{*}$, both the homogenous state (Figure 8(1)) and the periodic state presented in Figure 8(3) are stable under PNP-CH dynamics with periodic boundary conditions for $\phi$.

5.4 Effect of finite domain size

To study the effect of finite domain size, we conduct a numerical continuation study of the PNP-CH system (37,10b,38) in a finite domain $[-L,L]$ using pde2path 2.0 [40, 41]. Unlike a typical continuation study which considers stability of the solution with respect to the (elliptic) stationary equation, we test here stability with respect to PNP-CH dynamics by linking the continuation method to a Comsol 5.2a finite element solver. When a solution is found to be unstable, the solver returns the computed (stable) stationary solution it had converged to. Therefore, in addition to stability information, this method allows to efficiently map distinct solution branches by studying the stationary solutions computed by the solver. See Section 6 for numerical details.

The bifurcation diagram of the PNP-CH (37,10b,38) in the domain $[-L,L]$ where $L=3\pi/2k_{c}$, and otherwise with the same parameters as in Figure 8, is presented in Figure 9A. As expected, we recover the branch of periodic solutions with period $\pi/k_{c}$, see blue branch in Figure 9A, and a plot of a solution in point ‘a’ along this branch in Figure 9a. As noted, in free-space, solutions along this branch become stable under PNP-CH dynamics after the fold point $\sigma=\sigma^{*}$, see Figure 8C. In contrast, we observe that, in a finite domain, all solutions along this branch are unstable under PNP dynamics. Rather, additional branches with corresponding stable solutions emerge. The black branch in Figure 9A corresponds to periodic solutions with period $3\pi/k_{c}$ (point ‘b’, see Figure 9b) or with period $3\pi/2k_{c}$ (point ‘e’, see Figure 9e). Two additional branches (green and red with points ‘c’ and ‘d’, respectively) correspond to asymmetric solutions, see, e.g., Figure 9c. Note that since $\phi=V=0$ on both boundaries, there is no preferred direction for the solutions. Hence, every point along the green and red branches corresponds to two solutions which are reflections of each other. Hence, for example, at $\sigma=0.003$, we find that the system has six different stable stationary solutions: the homogenous state, the symmetric solution corresponding to point ‘b’, two solutions corresponding to a point on the red branch and two additional ones corresponding to a point on the green branch. Particularly, we observe multiple stable stationary solutions over the a wide range of $\sigma$.

5.5 Effect of applied voltage

Following the study of finite-domain effect in Section 5.4, we now apply numerical continuation methods [40, 41] to study the effect of applied voltage.

We consider the PNP-CH system (37,10b,38) in the domain $[-5,5]$ with applied voltage $\phi_{L}=-\phi_{-L}=1$, and otherwise with the same parameters as in Figure 8. The corresponding bifurcation diagram is presented in Figure 10 where solid curves correspond to solutions which are stable with respect to PNP-CH dynamics, and dashed curve correspond to unstable solutions. Figure 11 presents concentration profiles of solutions corresponding to points ‘a’-‘o’ along various branches presented in Figure 10. We observe that the bifurcation diagram is significantly more complex than the one presented in Figure 9 for solutions of a PNP-CH (37,10b,38) in the domain $[-L,L]$ where $L\approx 1.49$ and with no applied voltage $V=0$. As we shall see, this is mainly due to the larger domain size in Figure 10. Indeed, as expected, the solutions profiles are essentially identical near the charged boundaries, see Figure 11. The reason is that near a sufficiently charged wall, the concentration profiles near the boundaries are independent of the high-order gradient-energy terms. Namely, boundary layers of $O(1/\sqrt{\sigma})$ width do not develop in the concentration profiles, and rather dominant balance shows that

[TABLE]

see discussion before Lemma 5. The different branches, therefore, describe distinct behavior in the bulk, namely, different numbers, amplitudes and locations of peaks. Particularly, while the lowest branch in the bifurcation diagram, see Figure 10, corresponds to solutions with a uniform bulk, see Figure 11a, solutions ‘f’ and ‘o’ have one peak in the bulk, solutions ‘c’, ‘l’ and ‘m’ have two peaks, and solution ‘h’ with three peaks. Finally, solution ‘j’ has four peaks, see Figure 11. The maximal number of peaks is dictated by the domain size. Thus, the bifurcation diagram on a larger domain will have more branches. To distinguish between the number of peaks, and to distinguish between two solutions that have the same number of peaks but at different locations, e.g., solutions ‘l’ and ‘m’, we choose a weighted arc-length solution norm

[TABLE]

to plot bifurcations in the $(\sigma,\|c_{1}\|)$ plane.

6 Numerical methods

The stationary solutions of the PNP-steric equation, presented in Sections 3 and 4, were computed by solving the initial value problem (16,26) using Matlab’s ode45. As noted, this is the reason many of the boundary-value problem parameters in Figures 2-2 are not round numbers.

The numerical continuation study of the stationary solutions of the PNP-CH system (37,10b,38) in a finite domain $[-L,L]$ was conducted using pde2path 2.0 [40, 41]. Resolving the bifurcation diagram for these system, however, is a demanding problem for standard continuation procedure. Therefore, we utilize a combined method which integrates the numerical continuation procedure with simulations of the full PNP-CH system. To test the stability of the stationary solutions computed by pde2path, we use LiveLinkTM for Matlab to solve the PNP-CH system, with initial condition attained from pde2path output, via a Comsol 5.2a finite element solver. Whenever the stationary solution computed by pde2path is found to be unstable, the simulation of the PNP-CH system convergences to a new stable stationary solution. If the new stationary solution does not already reside on a known branch, it is used to find a new branch using numerical continuation. Thus, there is a mutual flow of data between the continuation procedure and the PDE solver.

The pseudo-code is as follows.

Find stationary solution using numerical continuation, starting with a given (possibly approximate) stationary solution ${\bf u}_{0}^{\rm p2p}$.

1. (a)

At step i

1. i.

Compute stationary solution ${\bf u}_{i}^{\rm p2p}$. 2. ii.

Solve the PNP-CH system with initial condition ${\bf u}(0,{\bf x})={\bf u}^{\rm p2p}_{i}$ until it reaches steady state at time $T$. Denote ${\bf u}^{\rm stationary}_{i}:={\bf u}(T,{\bf x})$. 3. iii.

The stationary solution ${\bf u}_{i}^{\rm p2p}$ is stable if

[TABLE]

and otherwise unstable. 4. iv.

If ${\bf u}_{i}^{\rm p2p}$ is unstable, and if ${\bf u}^{\rm stationary}_{i}$ does not reside on the same branch as the previously computed solution ${\bf u}^{\rm stationary}_{i-1}$, i.e.,

[TABLE]

then add ${\bf u}^{\rm stationary}_{i}$ to list $\mathcal{L}$ of solutions on potentially new branches. 2. 2.

Save branch data $\mathcal{B}_{k}=\{\sigma_{i},\|{\bf u}_{i}^{\rm p2p}\|\}$ 3. 3.

For each point in list ${\bf u}\in\mathcal{L}$ of solutions on potentially new branches,

1. (a)

verify that ${\bf u}$ does not reside on known branches $\{\mathcal{B}_{1},\mathcal{B}_{2},\cdots,\mathcal{B}_{k}\}$ (i.e., $\|{\bf u}\|$ is at distance larger than TOL from branch curve). 2. (b)

Repeat stage 1, starting the numerical continuation from ${\bf u}$, i.e., setting ${\bf u}_{0}^{\rm p2p}={\bf u}$

The combined method allows to map distinct solution branches more efficiently than considering solely the equations for the stationary states. For example, we found that numerical continuation from point ‘a’ retrieved the blue branch (associated with point ‘c’) but not the red branch (associated with point ‘’e’). Therefore, mapping the three unstable branches emanating from point ‘a’, see Figure 10, using numerical continuation from this point is demanding. Using the combined method, such issues are avoided. For example, solving the PNP-CH system with initial conditions corresponding to solutions along the blue branch connecting points ‘a’ and ‘b’, gives rise to a stable stationary solution of the PNP-CH system on the red branch connecting points ‘e’ and ‘f’. The red branch is then traced back from this point to point ‘a’ using continuation.

7 Conclusions

In this study, we have considered the existence and stability of stationary solutions of the PNP-steric system, and have shown that the corresponding stationary system gives rise to smooth multiple stationary solutions. The PNP-steric equation, however, turns out to be ill-posed at the parameter regimes where multiple solution arise. Following this finding, we introduced a novel PNP-Cahn-Hilliard model which is valid at high concentrations. We have shown that this model gives rise to multiple stationary solutions that are stable with respect to PNP-CH dynamics, and utilized bifurcation analysis and numerical continuation to map the various branches of stationary solutions. In particular, we observe that the stationary solution profiles are essentially the same near the interfaces, and differ only by their behavior in the bulk, namely in regions where the electric field applied from the boundaries is screened.

The PNP-steric model stems from the PNP-Lennard Jones (PNP-LJ) model which accounts for ion-ion steric interactions via the singular Lennard-Jones potential [23]. The complicated structure of the PNP-LJ equations, however, gives rise to an highly demanding computational problem and to an intractable analytic problem. Therefore, the PNP-steric model was derived using a leading order approximation of the Lennard-Jones potential [13]. The fact that PNP-steric model fails at high concentrations, as shown in this work, implies that higher order approximations of the Lennard-Jones potential are required. It is plausible that the second order approximation of the Lennard-Jones potential will give rise to gradient coefficient terms as in the PNP-CH model. Derivation of a high-order PNP-steric model by utilizing more accurate approximations of the LJ potential is left for future work.

The PNP-CH model differs from the broad family of generalized PNP models as it gives rise to bulk pattern formation. From a physical point of view, emergence of periodic patterns or pattern formation, in general, does not seem natural in electrolyte systems - salt dissolved in a cup of water is not typically observed to spontaneously self-assemble into spatial nano-structures. Self-assembly of bulk nano-structure, however, is well documented in room-temperature ionic liquids, see [42] and references within. Room-temperature ionic liquids are solvent-less mixtures of ions which remain liquid at room temperature. In a sense, they provide a limiting example of the most concentrated electrolyte solutions. Hence, sufficiently concentrated electrolyte solutions may also self-assemble into spatial nano-structures in the bulk [43].

Electrolyte behavior near charged boundaries, aka the electric double layer structure, has been at the spotlight of electrolyte research due to its importance within a wide range of areas including electrochemistry, biochemistry, physiology, and colloidal science. Recent works emphasize the significance of the coupling between bulk and interfacial nano-structures in such systems, [44, 45, 46, 47]. A study addressing the effect of bulk nano-structure on electrolyte behavior, and its coupling to the electric double layer structure will be presented elsewhere.

Finally, one of the questions that motivated the study of bi-stability in modified PNP systems is whether bi-stability can explain gating phenomena in ionic channels. This study does not address this question directly as the PNP-CH model with the no-flux boundary conditions considered in this study does not describe ion channels. Nevertheless, the multiple solutions presented in this study are all associated with multiple bulk states, and are not driven by boundary conditions, see, e.g., Figure 11. Thus, it is plausible that an appropriate PNP-CH model of ionic channels will give rise to similar multiple bulk states. A systematic study of the stationary solutions of the PNP-CH model in an ion channel setting, and their corresponding current levels, will be presented elsewhere.

Acknowledgements

This research was supported by the Technion VPR fund and from EU Marie–Curie CIG Grant 2018620. The author thanks Hannes Uecker for valuable advice in the application of pde2path.