Entropy method for generalized Poisson–Nernst–Planck equations

A proper mathematical model given by nonlinear Poisson–Nernst–Planck (PNP) equations which describe electrokinetics of charged species is considered. The model is generalized with entropy variables associating the pressure and quasi-Fermi electro-chemical potentials in order to adhere to the law of conservation of mass. Based on a variational principle for suitable free energy, the generalized PNP system is endowed with the structure of a gradient flow. The well-posedness theorems for the mixed formulation (using the entropy variables) of the gradient-flow problem are provided within the Gibbs simplex and supported by a-priori estimates of the solution.


Introduction
To describe electrokinetic transport occurring in micro-structures in many physical, chemical, and biological applications, a proper mathematical model adhering to the law of conservation of mass is suggested following the approach [5,9]. The reference cross-diffusion system of elliptic-parabolic type is described by nonlinear Poisson-Nernst-Planck (PNP) equations for concentrations of charged species and overall electrostatic potential. For physical consistency, they are generalized with entropy variables associating the pressure and quasi-Fermi electro-chemical potentials.
Based on a suitable free energy (see thermodynamic principles in [24]), in [20] a variational principle was established within the Gibbs simplex, thus preserving the total mass balance and non-negative species concentrations. In [18,19], the generalized PNP problem was stated in two-phase medium composed of pore and particle parts and taking into account for nonlinear interface reactions which are of primary importance in applications. Its rigorous asymptotic analysis was carried our in [7,8]. For a broad class of other relevant transport equations we refer to [6,12,15,22], to [13] for stochastic systems, and to [17,25] for variational principles.
Based on the entropy variables and following the thermodynamic formalism for cross-diffusion systems introduced in [11,16], in the current work we endow the generalized PNP problem with the structure of a gradient flow and analyze it. Within the entropy approach, the question of global solvability of related diffusion problems was investigated in [1,2,4]. For the general theory of linear and quasilinear parabolic equations we refer to [21]. However, the key issue of the entropy approach requires uniformly strongly elliptic property of the governing system. Unfortunately, the ellipticity fails under coupling cross-diffusive phenomena for the PNP problem, thus implying the degenerate case. For a study of degenerate elliptic operators, see [23].
In Sect. 3 we present well-posedness analysis following from the regularization approach by [26]. We set the entropy variables as independent ones. In the fully coupled case, the non-negativity of species concentrations might be lost during the time evolution. Otherwise, when the electro-chemical potentials are well-defined, then the species concentrations are expressed by a normalized canonical ensemble of Fermi-Dirac statistics, thus yielding the non-negativity and the total mass balance. Moreover, in the decoupled case, in Sect. 4 we prove directly well-posedness of the static equilibrium for the underlying problem. A rigorous derivation of energy and entropy estimates is collected in Appendix A.

Generalized PNP problem
We start with the geometry configuration. Let Ω ⊂ R d (with natural d ∈ N) be a connected domain with the Lipschitz boundary ∂Ω and the normal vector ν = (ν 1 , . . . , ν d ) outward to Ω. Here and in what follows the upper symbol stands for transposition swapping columns and rows. We split ∂Ω into two disjoint parts Γ D and Γ N corresponding to mixed Dirichlet-Neumann boundary conditions. By this consideration we associate Ω to a pore space with a bath boundary Γ D , which is Electric displacement complement to a solid space (bearing in mind possibly disconnected set of microparticles) with the boundary Γ N . For time t ∈ R + and spatial coordinates x = (x 1 , . . . , x d ) ∈ R d , we look for an unknown distribution over the cylinder (0, T ) × Ω (with the final time T > 0) of mass concentrations ρ(t, x) = (ρ 1 , . . . , ρ n ) (natural n ≥ 2) of charged species (ions) with electric charges z = (z 1 , . . . , z n ) , electro-chemical potentials μ(t, x) = (μ 1 , . . . , μ n ) , the overall electrostatic potential φ(t, x), and the pressure p(t, x) according to the generalization that was introduced in [5]. For convenience, all the physical variables and parameters of the model are gathered in Table 1.
Our modeling is based on the general law of cross-diffusion where the vector-valued diffusion fluxes J i (t, x) = ((J i ) 1 , . . . , (J i ) d ) are given by the constitutive law (see [5]) with coupling by means of diffusivity matrices D i j ∈ R d×d , i, j = 1, . . . , n.
Here and in what follows div stands for the divergence, and ∇ for the gradient. Inserting (1b) into (1a) implies a strongly nonlinear equation with respect to ρ and μ. The fluxes have to fulfill the mass conservation law: The electrostatic potential φ is described by the Poisson equation where the electric permittivity A ∈ R d×d . The Navier-Stokes equation (see e.g. [14]) with zero flow velocity results in the force balance The species concentrations should be physically consistent within a Gibbs simplex requiring non-negativity and preserving the total mass C > 0: Introducing the Lagrangian function of a free energy (see [20]) the governing laws (1) are completed with the thermodynamic equilibrium expressed by functional derivatives implying the Gibbs-Duhem equation for the electro-chemical potentials. It is worth noting that substitution of (2b) and (1b) into the diffusion equation (1a) leads to the gradient-flow structure Since p + K is defined by (1e) up to an additive constant K , all the μ i + u C K are determined by (2b) up to the same constant. Taking the gradient of (2b) and using the force balance (1e) leads to formulas [which will be useful later on to calculate the flux in (1b)] for i = 1, . . . , n where the functions Υ 1 , . . . , Υ n are uniformly bounded within the Gibbs simplex Moreover, equating the variation δE δφ of the function E in (2a) to zero leads to the Gauss law in the form of Poisson equation (1d) and the inhomogeneous Neumann boundary condition below in (3c) for φ. From the optimization viewpoint, the pressure p enters (2a) as a Lagrange multiplier to the equality constraint in (1f) implying δE δ p = 0. The elliptic-parabolic system of nonlinear equations in (1)-(2) is endowed with the standard initial condition ρ = ρ 0 as t = 0 ( 3 a ) and mixed Dirichlet-Neumann boundary conditions It is worth remarking that an inhomogeneous condition for the normal diffusion flux in (3c) would be well-posed only when it depends nonlinearly on ρ, this case was investigated in [18][19][20].
In order to guarantee the flux balance identity (1c), it suffices to assume with an elliptic matrix D ∈ R d×d . Indeed, substituting into (1c) the constitutive equations (1b) together with the expression for ρ i ∇μ i from (2c) and using the assumption (4a), after summation of the fluxes J i over i = 1, . . . , n we have since n j=1 (∇ρ j ) = ∇C = 0 and n j=1 Υ j (ρ) = 0 in (2c) due to the total mass balance in (1f). The assumption (4a) is related to quasi-stochastic matrices. In fact, , their sum in every column according to (4a) is equal to the same Such matrices D kl with non-negative entries are called column quasi-stochastic. The standard assumptions for solvability are the ellipticity and boundedness conditions for diffusivities: and for the permittivity: there exist 0 < a ≤ a such that Based on (1)-(4) now we give a weak variational formulation of the generalized PNP problem by excluding the entropy variables μ and p with the help of (2c). Find a pair of functions that satisfy the non-negativity and the total mass balance (1f), the Dirichlet condition (3b), and the following variational equations for i = 1, . . . , n The well-posedness to (5), (3b), and (1f) was investigated earlier in [18][19][20]. When solving problem (5), (3b), the key issue concerns fulfilling explicitly conditions (1f). In the following sections we consider the redundant entropy variable μ as independent one, thus allowing to include conditions (1f) implicitly in the problem formulation.

Entropy formulation of the PNP problem
When (5), (3b) is solved, multiplying (1e) with the gradient of a smooth test function and integrating the result over the domain, the pressure can be determined as a solution (defined up to a constant K ) to the elliptic equation for all test functions p ∈ H 1 (Ω). The next step is to recover μ from (2b). Taking ∇μ, multiplying it with the gradient of suitable test functions, and integrating over the cylinder, the corresponding electro-chemical potentials can be looked for as solutions to the mutually independent equations holding for all test functions μ i ∈ L 2 (0, T ; H 1 (Ω)), and satisfying the Dirichlet boundary condition In general, this elliptic problem is degenerate, because the operator of (6c) is unbounded due to the presence of factor 1 ρ i . To remedy, we give the following two conditional assertions, which will be justified in Theorem 1 later on.
Proof Indeed, if either condition (i) or condition (ii) is satisfied, then existence and uniqueness of the solution to (6b)-(6d) stated in Proposition 1 follows immediately from a general fact on elliptic systems.
From (2b) we get and summing up these equations over i = 1, . . . , n, due to the total mass balance n i=1 ρ i = C we can express the pressure as Excluding p from (7a) with the help of (7b), we bring the species concentrations in the form of a canonical ensemble of Fermi-Dirac statistics The normalized probabilities obey the inherent behavior Then the non-negativity and the total mass balance (1f) follow straightforwardly from the properties (7d). By this we observe that μ i − z i φ → −∞ would lead to the limit ρ i = 0 in (7c). This is an admissible behavior, as an example, for the function (μ i − z i φ)(x) = − ln | ln |x|| as x → 0, which agrees with the H 1 -spatial regularity of μ and φ stated in (5a) and (6b). In spite of the fact that μ + K is defined up to an additive constant K , the concentrations ρ i = P i (μ + K − zφ) in (7c) are defined uniquely. Based on (7), we reformulate the generalized PNP problem as follows. Find a triple of functions that satisfy the Fermi-Dirac statistics (7c), the Dirichlet condition (3b), and the following variational equations for i = 1, . . . , n for all test functions ρ = (ρ 1 , . . . , ρ n ) ∈ H 1 (0, T ; L 2 (Ω)) n ∩ L 2 (0, T ; H 1 (Ω)) n and φ ∈ H 1 (Ω) such that ρ(T ) = 0; ρ = 0 and φ = 0 at Γ D . The nonlinear parabolic equations (8b) imply a degenerate system because its operator loses the ellipticity property compared to (5b). Firstly, because of the crossdiffusion structure induced by matrices D i j . Second, due to the presence of the factor ρ j , which is not uniformly positive in general. These facts do not allow to apply the boundedness-by-entropy method [16].
Nevertheless, the well-posedness results established in [18,19] guarantee existence of a solution pair (ρ, φ) for the generalized PNP problem in the form of (5), (3b) that satisfies the total mass balance in (1f). From (5b) using (2c) the time derivative ∂ρ ∂t can be defined in L 2 (0, T ; H −1 (Ω)) n as a continuous linear functional   a solution (ρ, φ, μ). By this, (ρ, φ) solve the problem (5), (3b), and μ is from (2b). For such T 0 , the local a-priori estimates hold for all (t, Here we justify the a-priori estimates. In (8e), the bounds of ρ follow from (7c) and (7d), and the bound of φ is proved in Appendix A, Lemma 1. In (8d), the strict bounds of ρ are the consequence of continuity due to the strongly positive initial distribution (3d), and below we prove the estimate of μ. The lower bound of ρ i in (8d) allows us to apply Proposition 1. Based on this fact, we substitute p from (6a) and μ i = μ i into (6c) to obtain where the functions Υ i are defined by (2c). With the help of the Cauchy-Schwarz inequality, from the above equality we get where the upper bound of Υ i (ρ) is from (2d). Applying the estimate of ∇φ from (8e) and the bound of ∇ρ proved in [18][19][20] for arbitrary T > 0 For global in time solvability, we require a stronger than (4a) assumption with the Kronecker δ i j = 1 for i = j, and zero otherwise. The assumption (9a) implies (4a) and imposes decoupling in (8b) as well as (5b). In this case, existence of a solution (ρ, φ) to (5), (3b) satisfying both conditions in (1f) globally in time is proved in [18,19]. Therefore, the existence of the regular entropy variable μ from (2b) is sufficient to state the following theorem (see a relevant work [10]). The global a-priori estimates (8e) hold. If the Dirichlet data in (3b) are constant, i.e.
then the estimate with a constant K μ (T ) > 0 holds (see Appendix A, Lemma 2) In the next section we investigate an equilibrium state when T ∞.

Equilibrium state
Let lim T →∞ g(T , x) =: g ∞ (x) ∈ L 2 (Γ N ) in (3c), and the limit in (3b) be constant independent of x and satisfying according to (3d) the properties We consider a stationary counterpart of the entropy system (8), (7c), and (3b) under the decoupling assumption (9a). Find a triple of functions that satisfy the Fermi-Dirac statistics the Dirichlet condition and the following variational equations for i = 1, . . . , n for all μ ∈ H 1 (Ω) n and φ ∈ H 1 (Ω) such that μ = 0 and φ = 0 at Γ D . (11) exists, which is given by the following relations

Theorem 3 (Existence of equilibrium) A solution of problem
and the quasilinear variational equation for the electrostatic potential for all test functions φ ∈ H 1 (Ω) such that φ = 0 at Γ D .

A Appendix
Here we present auxiliary lemmas proving global a-priori estimates. The proof utilizes the standard boundary trace theorem and Poincare inequality as follows In Theorem 1, (8e) and in Theorem 3 there was used the first lemma.
Proof A solution to the quasilinear elliptic problem (15a) exists by the Schauder-Tikhonov fixed point theorem. Indeed, the principal term in (15a) is linear, coercive and bounded due to (4d). The nonlinear term is uniformly bounded according to (7d) and Lipschitz continuous, because ∂ P k (μ) Testing (15a) with φ = φ − φ 0 provided by zero on Γ D , due to (4d) we evaluate Applying Young's inequality to the following terms where (7d) and Z = n k=1 |z k | were used, after summation of these inequalities, with the help of (14) we arrive at the estimate (15b) and the bound (15c).
In Theorem 2, the estimate (9b) follows from the next lemma.
The time derivative of entropy can be calculated with the help of thermodynamic relations (2b) in the form (16f) Next we remark that the function μ i −z i φ − u C p = k B θ ln(β i ρ i ) is zero on Γ D because of the assumption (16a) and the Dirichlet condition (3b). Thus it can be inserted for a test function into (16b) and we get the following relations thanks to the lower bound in (4c) and (9a). From (2b) and (2c) we easily get Combining this equality with the upper bound in (4c) and (9a), and taking into account the inequality ρ i ≤ C and estimate (2d) for Υ i (ρ) ρ i , we can estimate the former integral in the right-hand side of (16g) from above using Young's inequality as follows: After summation over i = 1, . . . , n, the latter integral in the right-hand side of (16g) agrees with θ T 0 ∂ S ∂t dt due to (16f). Gathering the same terms in (16g)-(16i) and multiplying the resulting inequality with the factor 2 d , we conclude with the a-priori estimate (16c).