Non-equilibrium steady states of electrolyte interfaces

The non-equilibrium steady states of a semi-infinite quasi-one-dimensional univalent binary electrolyte solution, characterised by non-vanishing electric currents, are investigated by means of Poisson-Nernst-Planck (PNP) theory. Exact analytical expressions of the electric field, the charge density and the number density are derived, which depend on the electric current density as a parameter. From a non-equilibrium version of the Grahame equation, which relates the total space charge per cross-sectional area and the corresponding contribution of the electric potential drop, the current-dependent differential capacitance of the diffuse layer is derived. In the limit of vanishing electric current these results reduce to those within Gouy-Chapman theory. It is shown that improperly chosen boundary conditions lead to non-equilibrium steady state solutions of the PNP equations with negative ion number densities. A necessary and sufficient criterion on surface conductivity constitutive relations is formulated which allows one to detect such unphysical solutions.


I. INTRODUCTION
The dynamics of ions in external electric fields determines the properties of numerous important natural and technological processes such as the formation of a membrane potential in biological cells via ion channels [1][2][3][4][5][6][7], the charging and discharging of batteries by charge transfer reactions at electrodes [8,9], the motion of colloids exploiting electrokinetic effects of ionic solvent components [10][11][12] as well as the suppression of electric currents by insulation fluids in high electric fields [13][14][15][16].
Poisson-Nernst-Planck (PNP) theory [10,17,18] is an established and widely used theoretical framework in order to address these questions by considering the distributions of the electric field, the charge density the ion number densities etc.However, in studies of ion channels, high concentrated battery electrolytes and colloid migration steric effects of ions have been noted to be relevant [19][20][21] so that extensions of the PNP theory are currently under investigation [22,23].
In contrast, ideal insulation fluids for high voltage applications would be void of ions and carefully prepared real insulation fluids contain only a few.Hence ionic steric effects should be weak in that type of systems such that PNP theory can be considered a well justified starting point [24].However, in order to sustain the insulation property the ion concentration has to stay low for long times.Hence, understanding the mechanisms of charge generation in insulation fluids is an important topic which has been under investigation for many decades [15,16,[25][26][27][28][29].
Unlike for nano-sized ion channels or colloids, the functioning of insulation fluids requires a truly macroscopic * markus.bier@thws.deE bulk j Q FIG. 1.A semi-infinite electrolyte solution of univalent cations (positive ions, red) and univalent anions (negative ions, green) dispersed in a dielectric continuum (solvent, yellow) is in contact with a planar electrode (blue).Charge transfer processes at the electrode surface and a non-vanishing electric field E bulk in the bulk give rise to a non-equilibrium steady state with a non-vanishing, spatially uniform charge current density jQ (violet arrow).spatial extension of length scales of the order 1 cm and above.Hence a clear separation of length scales occurs: Within distances of molecular size from (typically metallic) surfaces electrochemical processes can occur, which provide the most important sources for ion generation, as electric fields are strongest there.Further away from the surfaces an extended bulk fluid with smooth distributions of the electric field and ions is present.The smoothness of these distributions allows for a description in terms of (partial) differential equations and the smallness of the ion number densities justifies the applicability of the PNP equations (see Sec. II B below).Moreover, the macroscopic character of high-voltage insulation systems quite naturally motivates to study electric currents theoretically in terms of a simple one-dimensional model comprising a planar electrode in contact with a semi-infinite electrolyte solution, similar to the investigations of Gouy and Chapman [10,11,[30][31][32][33] more than a century ago for conditions of thermodynamic equilibrium (see Fig. 1).
It is the purpose of the present work to analytically solve the PNP equations (Sec.II B) for non-equilibrium steady states of a semi-infinite quasi-one-dimensional univalent binary electrolyte solution (Sec.II A).It turns out that the solutions can be represented in terms of elementary functions and the expressions are of similar complexity as the Gouy-Chapman results for thermodynamic equilibrium (Sec.II C).Moreover, a generalisation of the Grahame equation, which for thermodynamic equilibrium expresses the surface charge density in terms of the surface potential [10,11,33], to non-equilibrium steady states is derived (Sec.II D).This allows to discuss the dependence of the distributions of the electric field as well as of charge and ion number density on the electric current (Sec.III A).Moreover, it is shown that the profiles as well as the Grahame relation and the differential capacitance of the space charge region reduce to the well-known Gouy-Chapman results for vanishing electric current (Sec.III B).A novel feature of the considered semi-infinite model, which is shown to occur only out of equilibrium, is the occurrence of solutions of the PNP equations with negative ion number densities (Sec.III C).This leads to the conclusion, that particular care is required in choosing appropriate boundary conditions in order to avoid such unphysical solutions (Sec.IV).

A. Setup
Consider a semi-infinite univalent binary electrolyte solution, which is bounded by a single planar electrode (see Fig. 1).The solvent is described as a dielectric continuum of temperature T and dielectric constant ε.Cations (positive ions, valency Z + = 1, diffusion constant D + , number density ̺ + ) and anions (negative ions, valency Z − = −1, diffusion constant D − , number density ̺ − ) migrate in an electric field E, which is oriented in normal direction of the electrode.Due to the planar symmetry, the electric field E(x) and the number densities ̺ ± (x) are functions of the distance x ≥ 0 from the electrode, but they are independent of the lateral position.Charge transfer processes close to the electrode surface allow for an electric current to occur in the system.
It is assumed that a sufficiently long waiting time has elapsed such that the system has attained a steady state in which the electric field E and the number densities ̺ ± are time-independent.Moreover, at large distances x → ∞ from the electrode these quantities are assumed to approach the constant bulk limits E(x) → E bulk and ̺ ± (x) → ̺ bulk /2.Consequently, for E bulk = 0 the system is in thermodynamic equilibrium, whereas for E bulk = 0 it is in a non-equilibrium steady state.In a non-equilibrium steady state a non-vanishing, spatially uniform electric current density j Q = 0 is present, which requires concomitant charge transfer processes to occur at the electrode surface.

B. Governing equations
In order to quantify the steady states of the system described in the previous Sec.II A the governing equations are derived within PNP theory in the following.
Given ion number densities ̺ ± (x) the total ion number density and the charge density with the elementary charge e are defined.At large distances x → ∞ from the electrode these quantities approach the limits ̺(x) → ̺ bulk and q(x) → 0. By means of Gauss' law [47] the derivative E ′ (x) of the electric field is related to the charge density q(x): where ε 0 is the vacuum permittivity.This shows that approaching a constant bulk electric field E(x) → E bulk in the limit x → ∞ implies local charge neutrality q(x) → 0 in the bulk.The Nernst-Planck equations [10,17,18] where 1/β = k B T with the Boltzmann constant k B denotes the thermal energy, describe the current densities of ions due to diffusion in a density gradient and drift in the electric field.Note that no advection contribution occurs in the present quiescent solvent.
In general the time dependence of the ion number densities ̺ ± (x, t) is given by the continuity equations which describe conservation of the number of ions.However, a steady state is time-independent ( ̺± = 0) such that the current densities j ± (x, t) are spatially uniform and time-independent: The same is true for the total number current density and the charge current density In order to simplify the later calculations the following reduced current densities are introduced: As ̺ ′ (x) → 0 (due to ̺(x) → ̺ bulk ), q(x) → 0 and E(x) → E bulk for x → ∞, Eq. ( 9) implies J N = 0, hence, from Eq. ( 9), Consequently which allows one to express the ionic current densities j ± in terms of the charge current density j Q .This leads to the relation between the charge current density j Q and the reduced charge current density J Q with the average diffusion constant Similarly, as q ′ (x) → 0 (due to q(x) → 0), ̺(x) → ̺ bulk and E(x) → E bulk for x → ∞, Eq. ( 10) leads to By using Eq. ( 13) one infers with the bulk conductivity

C. Analytical solution
In the following the PNP equations ( 3)-( 5) of the considered system are solved analytically.
In a first step inserting Eq. (3) in Eq. ( 9) yields which implies Evaluation of the constant by taking the limit x → ∞ leads to By inserting Eqs. ( 3) and (20) in Eq. ( 10) one obtains the inhomogeneous non-linear ordinary differential equation for the electric field E(x) This equation resembles Eq. (A.9) of Ref. [48] and Eq. ( 48) of Ref. [35].
Introducing the excess electric field ∆E(x) = E(x) − E bulk transforms the differential equation ( 21) to the homogeneous differential equation With the Debye length 1/κ defined by κ 2 = 4πℓ B ̺ bulk , where ℓ B = βe 2 /(4πε 0 ε) is the Bjerrum length [10,11,49,50], one can introduce the dimensionless electric flux parameter which quantifies the deviation of a steady state from thermodynamic equilibrium (η = 0) in terms of the charge current density j Q .The values η = ±1 correspond to a charge current density j Q = ±S bulk κ/(βe) generated by a bulk electric field E bulk = ±κ/(βe).The notion of electric flux η allows one to rewrite Eq. ( 22) in the form Finally the transformation is used to rewrite Eq. ( 24) in the form A first integral of Eq. ( 26) is found by multiplication with E ′ (y) and the integration constant is fixed by using E ′ (y) → 0 and E(y) → 0 for y → ∞: It can be shown, that the expression inside the parentheses is bounded from below by 1/(1 + η 2 ) for any value of E(y), hence Moreover, from Eq. ( 26) one infers that for y → ∞ the asymptotic behaviour of E(y) is monotonic, i.e. constant or exponential, but not oscillatory.Consequently, due to the limit E(y) → 0 for y → ∞, only three cases can occur: E(y) for y ≥ 0 is either (i) constantly zero ( E(y) = 0) or (ii) positive ( E(y) > 0) and strictly monotonically decreasing ( E ′ (y) < 0) or (iii) negative ( E(y) < 0) and strictly monotonically increasing ( E ′ (y) > 0).Otherwise if, say, E(y) approaches the limit E(y) → 0 for y → ∞ from above, but E(y) was not monotonically decreasing for all y ≥ 0, there would be a local maximum at some position y = y * with E ′ (y * ) = 0 but E(y * ) > 0, which contradicts Eq. ( 28).A similar contradiction arises from the assumption of a non-monotonically increasing behaviour when approaching the bulk limit from below.
In the above cases (ii) and (iii) E ′ (y) and E(y) have opposite sign, so that Eq. ( 27) leads to which relates the value E(y) and the derivative E ′ (y).This equation is obviously true also for case (i).
For cases (ii) and (iii) separation of variables leads to the solution [51] with the value E(0) = E 0 at the electrode surface y = 0 playing the role of the integration constant.The solution E(y) ≡ 0 of case (i) is obtained from Eq. ( 30) with E 0 = 0.

D. Grahame equation
An important result in the thermodynamic equilibrium theory of electrolyte interfaces is an expression for the surface charge density in terms of the surface potential, which is commonly referred to as Grahame equation [10,11,33].Here an analogous expression is derived for steady state conditions.
Integration of the excess electric field ∆E(x) leads to the excess voltage which measures deviations from local charge neutrality expressed by the charge density q(x) = ε 0 ε∆E ′ (x).It adds to the voltage required to sustain the bulk electric field E bulk .Using Eqs. ( 25) and ( 29) one obtains [51] βe∆U = 1 Using Eq. ( 3) the total charge per cross-sectional area of the diffuse layer can be calculated by This total charge per cross-sectional area of the space charge region is balanced by an excess surface charge density ∆σ at the electrode surface, which has the same magnitude but the opposite sign: ∆σ = ε 0 ε∆E(0).It expresses the excess of the total surface charge density σ = σ bulk + ∆σ (34) compared to the surface charge density σ bulk = ε 0 εE bulk , which generates the bulk electric field E bulk .By introducing the saturation charge density [52] one obtains the relation by means of which Eq. ( 32) can be rewritten as Solving Eq. ( 37) for the excess surface charge density ∆σ one obtains the Grahame equation The derivative of the excess surface charge density ∆σ w.r.t. the excess voltage ∆U leads to the differential capacitance of the diffuse layer = βeσ sat 4 1 + η 2 cosh βe∆U 2 +η sinh βe∆U 2 .
Hypothetical pure water of pH = 7, i.e. with total ion number density in the bulk of ̺ bulk = 2 • 10 −7 M, at T = 300 K has Bjerrum length ℓ B ≈ 0.7 nm and Debye length 1/κ ≈ 1 µm so that σ sat ≈ 7 • 10 −5 C/m 2 , which corresponds to an electric field E sat = σ sat /(ε 0 ε) ≈ 100 V/mm at the electrode.Furthermore, assuming a bulk conductivity S bulk = 5 • 10 −6 S/m (see Ref. [53]) an electric flux η = 1 corresponds to a charge current density of j Q ≈ 0.1 A/m 2 .However, these values can vary in a wide range for various materials, and the purpose of the  23).These quantities approach their bulk values on a length scale λ(η), which decreases with increasing magnitude |η| of the bulk flux (see Fig. 3).For E bulk = 0 (thermodynamic equilibrium) the electric field E(x) within Gouy-Chapman theory is reproduced, whereas purely exponential solutions are approached for η ≫ 1.
present work is not to discuss a particular system in detail.
The electric field E(x) = E bulk + ∆E(x) is obtained from the analytical solution Eq. ( 30) in conjunction with Eq. (25).The total number density ̺(x) is obtained from Eq. ( 20) and the charge density q(x) is calculated via Eq.(3).
For |η| ≫ 1 Eq.( 30) simplifies to with the length scale i.e. the excess electric field ∆E(x) becomes purely exponential.This is consistent with the fact that Eq. ( 26) reduces to the linear equation It is remarkable that the governing equation ( 26) becomes simple far away from thermodynamic equilibrium.According to Eq. ( 25) the position dependence of the excess electric field ∆E(x) for arbitrary values of the electric flux η is determined by length scale λ(η) in Eq. ( 42): In particular, according to Eq. ( 30), for x ≫ λ(η), ∆E(x) exhibits an exponential asymptotic decay on the length scale λ(η).
In summary, with increasing electric flux |η| ∼ |j Q | the electric field changes from the Gouy-Chapman form in thermodynamic equilibrium to a purely exponential dependence far away from thermodynamic equilibrium.In parallel the corresponding relevant length scale λ(η) decays from the Debye length 1/κ to zero.23)).For η = 0 (thermodynamic equilibrium) it is identical to the Debye length 1/κ.

B. Space charge
At distances x ≫ λ(η) the electric field and the densities become spatially uniform (see Fig. 2).Depending on the processes taking place at the electrode surface deviations from spatial uniformity can occur there, which are commonly referred to as the formation of space charges.In the theory of electrolyte solutions in thermodynamic equilibrium this space charge region is traditionally called the diffuse layer [10,11,33,54].
In the notation of Sec.II D the amount of space charge per cross-sectional area is given by −∆σ, where ∆σ denotes the excess surface charge density, which is the charge per cross-sectional area on the electrode surface induced by the space charge.The relation of ∆σ to the excess voltage ∆U of the space charge region (see Eq. ( 31)) is called the Grahame equation (38), in analogy to the equation of the same name within the theory of electrolyte solutions in thermodynamic equilibrium [10,11,33].
Figure 4(a) displays this relation for electric flux η ∈ {0, 2, 4, 6, 8, 10}.For η = 0 (thermodynamic equilibrium, j Q = 0) the traditional Grahame equation is found, whereas for a given excess voltage ∆U > 0 the amount of space charge increases upon increasing the electric flux η.According to Eq. ( 38), for η → ∞ this increase is asymptotically proportional to η. Hence the mean capacitance ∆σ/∆U increases upon increasing the electric flux, which can be attributed to the decrease of the diffuse layer thickness λ(η) upon increasing |η|.
one infers an asymptotically proportional dependence of C on η for η → ∞.This dependence on η can again be attributed to the decrease of the diffuse layer thickness λ(η) upon increasing |η|.
Note that within PNP theory no packing effects due to finite molecular sizes of the ions are taken into account.This precludes the decrease of the differential capacitance C for large values of the excess voltage ∆U , which otherwise would set in once the inner Helmholtz plane is fully occupied such that additionally adsorbed ions have to be accommodated at larger distances from the electrode surface [19,21,55,56].
For η = 0 (thermodynamic equilibrium, j Q = 0) condition Eq. ( 49) is fulfilled for all excess surface charge densities ∆σ, i.e. negative number density solutions within PNP theory can occur for non-equilibrium steady states, but not for states in thermodynamic equilibrium.
Relations between excess surface charge density ∆σ and electric flux η (see Eq. ( 23)) for three cases of surface conductivity models: (1) linear relation with surface-to-bulk conductivity ratio S surf /S bulk = 2 (blue line), (2) diffusion limited process with small saturation current (green curve) and ( 3) diffusion limited process with large saturation current (violet curve).The grey regions, bounded by red curves, correspond to unphysical conditions where solutions of the PNP equations ( 3)-( 5) occur which exhibit negative values of the number densities ̺±(0) close to the electrode.
In order to demonstrate the applicability of the previous result consider the particular case of surface processes which give rise to a strictly linear constitutive relation between charge current j Q and total electric field E(0) at the electrode surface: The proportionality factor S surf is called the surface conhere.Using Eqs.(23) and (36) one can rewrite Eq. ( 52) of the linear surface conductivity model as The case of a surface-to-bulk conductivity ratio S surf /S bulk = 2 is represented by the blue line labelled with "(1)" in Fig. 5.The fact that this line crosses over into the grey regions, where the conditions Eqs. ( 50) and ( 51) are violated, shows that the purely linear surface conductivity model cannot be applied under these conditions for too large electric fluxes.Obviously the same argument applies to any system with S surf > S bulk (high surface conductivity), because then the slope of the corresponding line is negative so that intersections with the unphysical grey regions occur for sufficiently large electric flux |η|.It should be noted that in the grey regions of Fig. 5 no mathematical problems arise: Equations ( 30), ( 45) and ( 46) are the solutions of the PNP equations for steady states (see Sec. II B).But these steady state solutions of the PNP equations may be physically meaningless due to the occurrence of negative number densities.
In systems with S surf ≤ S bulk (low surface conductivity) the linear model Eq. ( 52) leads to a line Eq.( 53) in Fig. 5 with non-negative slope, which does not intersect the grey regions, i.e. no negative number densities occur under these conditions.However, it is possible that other quantities exist, for which the PNP solutions exhibit unphysical properties.Moreover, whether the linear surface conductivity model Eq. ( 52), even if no unphysical values occur, is able to quantitatively describe real systems is an unrelated question.
In fact, the linear surface conductivity model Eq. ( 52) can be expected to be an acceptable description for sufficiently small surface fields only, because for large surface fields saturation of the charge density current j Q sets in due to an exhaustion of ions.Such diffusion limited surface processes can be described by the constitutive relation [8,9] where j Qsat > 0 is the saturation charge current density and S surf is the differential conductivity in the limit of infinitesimal surface electric fields.In Fig. 5 the case of S surf /S bulk = 2 for two different values of the saturation charge current density j Qsat are shown by the curves labelled "(2)" and "(3)".The green curve "(2)" demonstrates that for sufficiently small |j Qsat | no negative ion number densities occur, although a highly conductive surface is present at weak surface fields.However, the violet curve "(3)" for too large |j Qsat | does exhibit unphysical solutions inside the grey regions.Hence, great care is required to choose appropriate surface conductivity models, which, in conjunction with the PNP equations, lead to physical solutions.

IV. CONCLUSIONS
In the present work the analytical solution Eq. (30) of the PNP equations for steady states (see Sec. II B) of a semi-infinite univalent binary electrolyte solution in contact with a planar electrode (Fig. 1) has been derived.It can be expected to play a similar role as the Gouy-Chapman solution [10,11,[30][31][32][33] for thermodynamic equilibrium, to which the derived solutions reduce for the case of a vanishing charge current density (Fig. 2).The characteristic length scale of the electric field as well as the number and charge density profiles, which in thermodynamic equilibrium is given by the Debye length, decreases for non-equilibrium steady states upon increasing the magnitude of the charge current density (Fig. 3).The Grahame equation, which expresses the surface charge density at the electrode in terms of the voltage [10,11,33], is generalised to nonequilibrium steady states (Eq.( 38)).The excess surface charge density at the electrode and the differential capacitance of the space charge layer for given excess voltage are found to vary with the current charge density of nonequilibrium steady states (Fig. 4).Finally it is found that, in contrast to the case of thermodynamic equilibrium within Gouy-Chapman theory [10,11,[30][31][32][33], the excess surface charge density may not take an arbitrary value for a given non-vanishing charge current density of a non-equilibrium steady state: Steady state solutions of the PNP equations exist which give rise to physically meaningless negative ion number densities (Fig. 5).A concise criterion is formulated which can serve to identify such unphysical solutions (Eq.( 49)).
The most remarkable observation of the present work, i.e. the existence of steady state solutions of the PNP equations which are physically meaningless, calls for further investigation.Two approaches are conceivable: First, as Gauss' law Eq.( 3) and the continuity equation (5) are unexceptionable for fundamental physical reasons, one could suggest to modify the Nernst-Planck equation ( 4) in order to avoid unphysical solutions.Second, one could try to keep the Nernst-Planck equation (4) unchanged, but require boundary conditions to fulfil Eq. (49).Whereas the second approach is the more pragmatic one, it will be the first one which provides more fundamental insight into the non-equilibrium properties of electrolyte solutions.

FIG. 2 .
FIG.2.Steady state distributions of (a) the electric field E, (b) the cation number density ̺+ and (c) the anion number density ̺− as functions of the normal distance x from the electrode surface for fixed total surface charge density σ = 2.5 σsat and various values of the electric flux η defined in Eq. (23).These quantities approach their bulk values on a length scale λ(η), which decreases with increasing magnitude |η| of the bulk flux (see Fig.3).For E bulk = 0 (thermodynamic equilibrium) the electric field E(x) within Gouy-Chapman theory is reproduced, whereas purely exponential solutions are approached for η ≫ 1.

FIG. 4 .
FIG. 4. (a) Excess surface charge density ∆σ and (b) differential capacitance of the diffuse space charge region close to the electrode as functions of the corresponding voltage drop ∆U for various values of the electric flux η in Eq.(23).For η = 0 (thermodynamic equilibrium) these quantities are identical to those within Gouy-Chapman theory.