Analysis of a degenerate parabolic cross-diffusion system for ion transport

TL;DR

This paper analyzes a complex cross-diffusion system modeling ion transport through biological membranes, establishing existence, uniqueness, and numerical insights into long-term behavior despite mathematical challenges.

Contribution

It extends the boundedness-by-entropy method to nonhomogeneous boundary conditions and proves global existence and uniqueness for a degenerate, coupled ion transport system.

Findings

Existence of bounded weak solutions is proven.

Uniqueness of solutions is established under certain regularity.

Numerical simulations show slow equilibration rates.

Abstract

A cross-diffusion system describing ion transport through biological membranes or nanopores in a bounded domain with mixed Dirichlet-Neumann boundary conditions is analyzed. The ion concentrations solve strongly coupled diffusion equations with a drift term involving the electric potential which is coupled to the concentrations through a Poisson equation. The global-in-time existence of bounded weak solutions and the uniqueness of weak solutions under moderate regularity assumptions are shown. The main difficulties of the analysis are the cross-diffusion terms and the degeneracy of the diffusion matrix, preventing the use of standard tools. The proofs are based on the boundedness-by-entropy method, extended to nonhomogeneous boundary conditions, and the uniqueness technique of Gajewski. A finite-volume discretization in one space dimension illustrates the large-time behavior of the numerical solutions and shows that the equilibration rates may be very small.

Equations

$$\partial_{t}u_{i}=\operatorname{div}J_{i},\quad J_{i}=D_{i}\big{(}u_{0}\nabla u_{i}-u_{i}\nabla u_{0}+u_{0}u_{i}(\beta z_{i}\nabla\Phi+\nabla W_{i})\big{)}$$

$$-{\lambda }^2\Delta \Phi =\sum _{i=1}^n{z}_i{u}_i+f,$$

$$J_i\cdot \nu =0\text{mbox}on{\Gamma }_N,$$

$$\nabla \Phi \cdot \nu =0\text{mbox}on{\Gamma }_N,$$

$$u_i(\cdot ,0)=u_i^0\text{mbox}in\Omega ,i=1,\dots ,n.$$

$$\partial_{t}u_{i}=\operatorname{div}\bigg{(}\sum_{j=1}^{n}A_{ij}(u)\nabla u_{j}+D_{i}u_{0}u_{i}\nabla F_{i}\bigg{)},$$

$$A_{ii}(u)=D_i{u}_i,A_{ij}(u)=D_i(u_0+u_i),j\ne i.$$

$$H(u)={\int }_{\Omega }h(u)dx,\text{mbox}whereu=(u_1,\dots ,u_n),$$

$$h(u)=\sum _{i=0}^n{\int }_{u_i^D}^{u_i}\log \frac{s}{u_i^D}ds+\frac{\beta {\lambda }^2}{2}|\nabla (\Phi -{\Phi }^D){|}^2+\sum _{i=1}^n{u}_i{W}_i$$

$$\partial_{t}u_{i}=\operatorname{div}\bigg{(}\sum_{j=1}^{n}B_{ij}\nabla w_{j}\bigg{)},\quad i=1,\ldots,n,$$

$$\frac{\partial h}{\partial u_i}=w_i-w_i^D,\text{mbox}where$$

$$w_i=\log \frac{u_i}{u_0}+\beta z_i\Phi +W_i,w_i^D=\log \frac{u_i^D}{u_0^D}+\beta z_i{\Phi }^D,i=1,\dots ,n.$$

$$u_i=u_i(w)=\frac{\exp (w_i-\beta z_i\Phi -W_i)}{1+{\sum }_{j=1}^n\exp (w_j-\beta z_j\Phi -W_j)},i=1,\dots ,n,$$

$$u_{i}(w)\in{\mathcal{D}}:=\bigg{\{}u=(u_{1},\ldots,u_{n})\in(0,1)^{n}:\sum_{i=1}^{n}u_{i}<1\bigg{\}}.$$

$$-{\lambda }^2\Delta {\Phi }^D=f\text{mbox}in\Omega ,\nabla {\Phi }^D\cdot \nu =0\text{mbox}on{\Gamma }_N.$$

$$H_D^1(\Omega )=\{u\in H^1(\Omega ):u=0\text{mbox}on{\Gamma }_D\}.$$

$$u_i{u}_0^{1/2},u_0^{1/2}\in L^2(0,T;H^1(\Omega )),{\partial }_t{u}_i\in L^2(0,T;H_D^1(\Omega {)}^{\prime }),$$

$$\Phi \in L^2(0,T;H^1(\Omega )),i=1,\dots ,n,$$

$$\displaystyle\int_{0}^{T}\langle\partial_{t}u_{i},\phi_{i}\rangle dt+D_{i}\int_{0}^{T}\int_{\Omega}u_{0}^{1/2}\big{(}\nabla(u_{0}^{1/2}u_{i})-3u_{i}\nabla u_{0}^{1/2}\big{)}\cdot\nabla\phi_{i}dxdt$$

$$+D_i{\int }_0^T{\int }_{\Omega }{u}_i{u}_0(\beta z_i\nabla \Phi +\nabla W_i)\cdot \nabla {\varphi }_{i}dxdt=0,$$

$$\displaystyle\lambda^{2}\int_{0}^{T}\int_{\Omega}\nabla\Phi\cdot\nabla\theta dxdt=\int_{0}^{T}\int_{\Omega}\bigg{(}\sum_{i=1}^{n}z_{i}u_{i}+f\bigg{)}\theta dxdt,$$

$$u_0=u_0^D:=1-\sum _{i=1}^n{u}_i^D,u_i{u}_0^{1/2}=u_i^D(u_0^D{)}^{1/2}\text{mbox}on{\Gamma }_D,i=1,\dots ,n,$$

$$\frac{dH}{dt}$$

$$\le -\frac{1}{2}{\int }_{\Omega }\sum _{i=1}^n{D}_i{u}_0{u}_i|\nabla w_i{|}^{2}dx+C(w^D),$$

$$\sum _{i=1}^n{u}_i{u}_0\nabla \log \frac{u_i}{u_0}=4u_0\sum _{i=1}^n|\nabla u_i^{1/2}{|}^2+|\nabla u_0{|}^2+4|\nabla u_0^{1/2}{|}^2,$$

$$\partial_{t}u_{0}=\operatorname{div}\big{(}\nabla u_{0}-u_{0}(1-u_{0})(\beta z\nabla\Phi+\nabla W)\big{)},\quad-\lambda^{2}\Delta\Phi=z(1-u_{0})+f(x),$$

$$\frac{d}{dt}{\int }_{\Omega }(u_0-v_0)(t{)}^{2}dx\le C(\Phi ){\int }_{\Omega }(u_0-v_0{)}^{2}dx,$$

$$d(u,v)=\int_{\Omega}\sum_{i=1}^{n}\bigg{(}h(u_{i})+h(v_{i})-2h\bigg{(}\frac{u_{i}+v_{i}}{2}\bigg{)}\bigg{)}dx\geq 0,$$

$$-{\lambda }^2\Delta \Phi =\sum _{i=1}^n{z}_i{u}_i(w,\Phi )+f,u_i(w,\Phi )=\frac{\exp (w_i-\beta z_i\Phi -W_i)}{1+{\sum }_{j=1}^n\exp (w_j-\beta z_j\Phi -W_j)}$$

$$\Vert \Phi {\Vert }_{H^1(\Omega )}\le C(1+\Vert {\Phi }^D{\Vert }_{H^1(\Omega )}),$$

$$\displaystyle\frac{1}{\tau}\int_{\Omega}\big{(}u(v^{k}+w^{D},\Phi^{k})-u(v^{k-1}+w^{D},\Phi^{k-1})\big{)}\cdot\phi dx$$

$$+{\int }_{\Omega }\nabla \varphi :B(v^k+w^D,{\Phi }^k)\nabla (v^k+w^D)dx$$

$$\displaystyle\phantom{xxxx}{}+\varepsilon\int_{\Omega}\bigg{(}\sum_{|\alpha|=m}D^{\alpha}v^{k}\cdot D^{\alpha}\phi+v^{k}\cdot\phi\bigg{)}dx=0,$$

$$\displaystyle\lambda^{2}\int_{\Omega}\nabla\Phi^{k}\cdot\nabla\theta dx=\int_{\Omega}\bigg{(}\sum_{i=1}^{n}z_{i}u_{i}(v^{k}+w^{D},\Phi^{k})+f\bigg{)}\theta dx$$

$${\int }_{\Omega }\nabla \varphi$$

Full text

The authors acknowledge partial support from the Austrian Science Fund (FWF), grants P27352, P30000, and W1245, and from the Austrian-French project of the Austrian Exchange Service (ÖAD), grant FR 04/2016

1. Introduction

The transport of ions through membranes or nanopores can be described on the macroscopic level by the Poisson-Nernst-Planck equations, modeling ionic species and an electro-neutral solvent in the self-consistent field [19]. The equations can be derived in the mean-field limit from microscopic particle models [18] and lead to diffusion equations, satisfying Fick’s law for the fluxes. This ansatz breaks down in narrow ion channels if the finite size of the ions is taken into account. Including size exclusion, the mean-field model, derived from an on-lattice model in the diffusion limit [4, 21] or taking into account the combined effect of the excess chemical potentials [17], leads to parabolic equations with cross-diffusion terms. The aim of this paper is to analyze the cross-diffusion system of [4].

1.1. Model equations

The evolution of the ion concentrations (volume fractions) $u_{i}$ and fluxes $J_{i}$ of the $i$th species is governed by the equations

[TABLE]

for $i=1,\ldots,n$, where $u_{0}=1-\sum_{i=1}^{n}u_{i}$ is the concentration (volume fraction) of the solvent. We have assumed that the molar masses are the same for all species. Varying molar masses are considered in, e.g., [6, 8] in the context of the Maxwell-Stefan theory. The classical Nernst-Planck equations are obtained after setting $u_{0}=1$ [5]. They can be also coupled with fluiddynamical equations; see, e.g., [24]. Modified Nernst-Planck models without volume filling, but including cross-diffusion terms, were suggested and analyzed in [13, 16].

In equations (1), $D_{i}>0$ denotes the diffusion coefficients, $\beta=q/(k_{B}\theta)>0$ is the inverse thermal voltage (or inverse thermal energy) with the elementary charge $q$, the Boltzmann constant $k_{B}$, and the temperature $\theta$, $z_{i}\in{\mathbb{R}}$ is the valence of the $i$th species, and $W_{i}=W_{i}(x)$ is an external potential. Note that Einstein’s relation between the diffusivity $D_{i}$ and the mobility $\mu_{i}=qD_{i}/(k_{B}\theta)=D_{i}\beta$ holds. The electrical potential $\Phi$ is determined by the Poisson equation

[TABLE]

where $\lambda>0$ is the (scaled) permittivity, $\sum_{i=1}^{n}z_{i}u_{i}$ is the total charge density, and $f=f(x)$ is a permanent charge density.

Equations (1)-(2) are solved in the bounded domain $\Omega\subset{\mathbb{R}}^{d}$ ($d\geq 1$). Its boundary is supposed to consist of an insulating part $\Gamma_{N}$, on which no-flux boundary conditions are prescribed, and the union $\Gamma_{D}$ of boundary contacts with external reservoirs, on which the concentrations are fixed. The electric potential is influenced by the voltage at $\Gamma_{E}$ between two electrodes, and we assume for simplicity that $\Gamma_{E}=\Gamma_{D}$. This leads to the mixed Dirichlet-Neumann boundary conditions

[TABLE]

Finally, we prescribe the initial conditions

[TABLE]

Equations (1) can be written as the cross-diffusion system

[TABLE]

where $F_{i}=\beta z_{i}\Phi+W_{i}$ is the effective potential and the diffusion matrix $(A_{ij}(u))$ is defined by

[TABLE]

Mathematically, this system is strongly coupled with a nonsymmetric and generally not positive semidefinite diffusion matrix such that the existence of solutions to (6) is not trivial. A second difficulty is the fact that a maximum principle is generally not available for cross-diffusion systems, and the proof of nonnegativity of $u_{0}=1-\sum_{i=1}^{n}u_{i}$ is unclear. The third problem arises due to the degenerate structure hidden in the equations (see below for details).

For vanishing potentials $F_{i}=0$, the global existence of bounded weak solutions to (6) with no-flux boundary conditions has been shown in [25], based on the boundedness-by-entropy method [14, 15]. The existence of weak solutions to the (easier) stationary problem was proved in [4]. Related models were analyzed recently in [2]. No existence or uniqueness results for solutions to the full transient model (1)-(5) seem to be available in the literature and in this paper, we fill this gap. Compared to the works [14, 25], the novelty here is the inclusion of the electric potential and the mixed Dirichlet-Neumann boundary conditions, which need to be treated in a careful way.

1.2. Key idea of the analysis

We extend the boundedness-by-entropy method [14] to the case of nonconstant potentials and nonhomogeneous boundary conditions. The key observation, already stated in [4], is that (1) possesses an entropy or gradient-flow structure. The entropy or, more precisely, free energy is given by

[TABLE]

and $u_{0}^{D}=1-\sum_{i=1}^{n}u_{i}^{D}$. The free energy is bounded from below if $u_{i}\in L^{\infty}(\Omega)$ and $W_{i}\in L^{1}(\Omega)$. Equations (6) can be written as a formal gradient flow in the sense

[TABLE]

where $B_{ii}=D_{i}u_{0}u_{i}$, $B_{ij}=0$ if $i\neq j$ provide a diagonal positive semidefinite matrix $(B_{ij})$, and $w_{j}$ are the entropy variables, defined by

[TABLE]

We refer to Lemma 7 below for the computation of $\partial h/\partial u_{i}$. In thermodynamics $\partial h/\partial u_{i}$ is called the chemical potential of the $i$th species. The advantage of formulation (8) is that the drift terms are eliminated and, in this special case, the new diffusion matrix $(B_{ij})$ is diagonal. Note that we have not included the boundary data into the formulation (8). In fact, the free energy is nonincreasing along trajectories to (1)-(5) only if the boundary data are in equilibrium, i.e. if $\nabla w^{D}_{i}=0$. In the general case, the free energy is bounded only; see (12) below.

There is another important benefit of formulation (8). Observing that the relation between $w=(w_{1},\ldots,w_{n})$ and $u=(u_{1},\ldots,u_{n})$ can be inverted explicitly according to

[TABLE]

we see that, if $(w_{1},\ldots,w_{n},\Phi)$ is a solution to (2) and (8),

[TABLE]

This provides positive lower and upper bounds for the concentrations $u_{0},\ldots,u_{n}$ without the use of a maximum principle.

1.3. Main results

We prove (i) the global-in-time existence of bounded weak solutions, (ii) the uniqueness of weak solutions under additional regularity assumptions, and (iii) some numerical results on the large-time behavior of solutions in one space dimension. In the following, we detail these results. First, we specify the technical assumptions.

{labeling}

(A44)

Domain: $\Omega\subset{\mathbb{R}}^{d}$ ($d\geq 1$) is a bounded domain with $\partial\Omega=\Gamma_{D}\cup\Gamma_{N}\in C^{0,1}$, $\Gamma_{D}\cap\Gamma_{N}=\emptyset$, $\Gamma_{N}$ is open in $\partial\Omega$, and $\mbox{meas}(\Gamma_{D})>0$.

Parameters: $T>0$, $D_{i}$, $\beta>0$, and $z_{i}\in{\mathbb{R}}$, $i=1,\ldots,n$.

Given functions: $f\in L^{\infty}(\Omega)$, $W_{i}\in H^{1}(\Omega)\cap L^{\infty}(\Omega)$, and $W_{i}=0$ on $\Gamma_{D}$, $\nabla W_{i}\cdot\nu=0$ on $\Gamma_{N}$, $i=1,\ldots,n$.

Initial and boundary data: $u_{i}^{0}\in L^{\infty}(\Omega)$, $u_{i}^{D}\in H^{1}(\Omega)$, $u_{i}^{0}>0$, $u_{i}^{D}>0$, $1-\sum_{i=1}^{n}u_{i}^{0}>0$, $1-\sum_{i=1}^{n}u_{i}^{D}>0$ in $\Omega$ for $i=1,\ldots,n$, and $\Phi^{D}\in H^{1}(\Omega)\cap L^{\infty}(\Omega)$ satisfies

[TABLE]

Clearly, it is sufficient to define the functions $u_{i}^{D}$, $\Phi^{D}$ on $\Gamma_{D}$. By the extension property, they can be extended to $\Omega$, and we assume in (A4) that the extension of $\Phi^{D}$ is done in a special way. This extension is needed to be consistent with the definition of the free energy (entropy) and the entropy variables; see Lemma 7. We denote these extensions again by $u_{i}^{D}$, $\Phi^{D}$. Furthermore, we introduce the space [23]

[TABLE]

The first result concerns the existence of bounded weak solutions.

Theorem 1** (Global existence of weak solutions).**

Let Assumptions (A1)-(A4) hold. Then there exists a bounded weak solution $u_{1},\ldots,u_{n}:\Omega\times(0,T)\to\overline{{\mathcal{D}}}$ to (1)-(5) satisfying

[TABLE]

and the weak formulation

[TABLE]

for all $\phi_{i}$, $\theta\in L^{2}(0,T;H^{1}_{D}(\Omega))$, $i=1,\ldots,n$. The initial condition is satisfied in the sense of $H^{1}_{D}(\Omega)^{\prime}$, and the Dirichlet boundary conditions are given by

[TABLE]

in the sense of traces in $L^{2}(\Gamma_{D})$.

The proof is based on an approximation procedure, i.e., we prove first the existence of solutions $u_{0}^{(\tau)}$, $u_{i}^{(\tau)}$ to a regularized problem with approximation parameter $\tau>0$ and then pass to the limit $\tau\to 0$. The estimates needed for the compactness argument are coming from a discrete version of the entropy-production inequality (for simplicity, we omit the superindex $\tau$)

[TABLE]

where the constant $C(w^{D})>0$ depends on the $H^{1}(\Omega)$ norm of $w^{D}$. We show in (23) below that

[TABLE]

which yields an $H^{1}(\Omega)$ estimate for $u_{0}^{1/2}$ but not for $u_{i}$ because of the factor $u_{0}\geq 0$. This reflects the degenerate nature of the equations which is more apparent in the component-wise formulation $\partial_{t}u_{i}=\operatorname{div}(D_{i}u_{0}u_{i}\nabla w_{i})$ (see (8)).

To overcome this degeneracy, we employ the technique developed in [3, 25]. We show that $(u_{0}^{(\tau)}u_{i}^{(\tau)})$ is bounded in $H^{1}(\Omega)$ and that the (approximative) time derivative of $u_{i}^{(\tau)}$ is bounded in $H^{1}_{D}(\Omega)^{\prime}$. If $u_{0}^{(\tau)}$ was strictly positive, we could apply the Aubin-Lions lemma to conclude strong convergence of (a subsequence of) $(u_{i}^{(\tau)})$ to some $u_{i}$ which solves (1). However, since $u_{0}^{(\tau)}$ may vanish in the limit, this lemma cannot be used. The idea is to compensate the lack of the gradient estimates for $u_{i}^{(\tau)}$ by exploiting the uniform estimates for $u_{0}^{(\tau)}$. Then, by the “degenerate” Aubin-Lions lemma (see, e.g., [14, Appendix C]), (a subsequence of) $(u_{0}^{(\tau)}u_{i}^{(\tau)})$ converges strongly to $u_{0}u_{i}$, and $u_{0}$, $u_{i}$ solve (1). For details, see Section 2.

Remark 2**.**

1. Theorem 1 also holds when reaction terms $f_{i}(u)$ are introduced on the right-hand side of (1). As in [14], we need that $f_{i}$ is continuous and $\sum_{i=1}^{n}f_{i}(u)(\partial h/\partial u_{i})\leq C(1+h(u))$ holds for some $C>0$ and all $u\in{\mathcal{D}}$.

2. The approximate solution satisfies a discrete version of the entropy-production inequality; see (17). As explained above, the sequence $(u_{i}^{(\tau)})$ may not converge strongly, such that we are unable to perform the limit $\tau\to 0$ in (17). As a consequence, we cannot prove that the free energy (7) is nonincreasing along trajectories of (1)-(2), and the analysis of the large-time behavior seems to be inaccessible. Therefore, we investigate the decay of $H(u)$ numerically; see Section 4.

3. Since the Neumann boundary condition does not appear explicitly in the weak formulation (10)-(11), we do not need to make expressions like $\nabla\Phi\cdot\nu=0$ on $\Gamma_{N}$ precise. We only mention along the way that terms like $\nabla\Phi\cdot\nu$ on $\Gamma_{N}$ have to be understood in the sense of $H_{00}^{1/2}(\Gamma_{N})^{\prime}$ which is the dual space of $H_{00}^{1/2}(\Gamma_{N})$ consisting of all functions $v$ on $\Gamma_{N}$ such that $v\in H^{1}_{D}(\Omega)$. This space is larger than $H^{-1/2}(\Gamma_{N})$. We refer to [1, Chapter 18] for details. ∎

The second result is the uniqueness of weak solutions.

Theorem 3** (Uniqueness of weak solutions).**

Let Assumptions (A1)-(A4) hold, $\sum_{i=1}^{n}W_{i}\in L^{\infty}(0,T;W^{1,d}(\Omega))$, and let $D_{i}=1$ and $z_{i}=z\in{\mathbb{R}}$ for $i=1,\ldots,n$. Then there exists at most one bounded weak solution to (1)-(5) in the class of functions $u_{i}\in H^{1}(0,T;H_{D}^{1}(\Omega)^{\prime})\cap L^{2}(0,T;H^{1}(\Omega))$, $\Phi\in L^{\infty}(0,T;W^{1,q}(\Omega))$ with $q>d$.

The proof is a combination of standard $L^{2}(\Omega)$-type estimates and the entropy method of Gajewski [9]. In fact, equations (1) partially decouple because of the assumptions $D_{i}=1$ and $z_{i}=z$. Summing (1) over $i=1,\ldots,n$, we find that $(u_{0},\Phi)$ solves

[TABLE]

where $W=\sum_{i=1}^{n}W_{i}$. The uniqueness of solutions is shown by taking two solutions $(u_{0},\Phi)$ and $(v_{0},\Psi)$ and using $u_{0}-v_{0}$ as a test function in the first equation of (13). Then, with the Gagliardo-Nirenberg inequality and the hypothesis $\nabla\Phi\in L^{q}(\Omega)$, we show that

[TABLE]

where $C(\Phi)>0$ depends on the $W^{1,q}(\Omega)$ norm of $\Phi$. Hence, Gronwall’s lemma yields $u_{0}=v_{0}$ and consequently, $\Phi=\Psi$.

The next step is to show, for given $u_{0}$ and $\Phi$, that $u_{i}$ is the unique solution to (1). Since we cannot expect that $\nabla u_{i}\in L^{q}(\Omega)$, $q>d$, for $d\geq 3$, we employ the technique of Gajewski [9] which avoids this regularity. The method seems to work only for linear mobilities $u_{i}$, which is the reason why we cannot apply it to (13). The idea is to introduce the semimetric

[TABLE]

where $h(s)=s(\log s-1)+1$, and to show that $\partial_{t}d(u,v)\leq 0$. Since $d(u(0),v(0))=0$, this implies that $d(u(t),v(t))=0$ for $t>0$ and consequently, $u(t)=v(t)$. Since expressions like $\log u_{i}$ are undefined when $u_{i}=0$, we need to regularize the semimetric. For details, we refer to Section 3.

Remark 4**.**

1.  The regularity $u_{i}\in L^{2}(0,T;H^{1}(\Omega))$ holds if $u_{0}$ is strictly positive. A standard idea for the proof is to employ $\min\{0,u_{0}-me^{-\lambda t}\}^{p}$ as a test function in the first equation of (13), where $\inf_{\Gamma_{D}}u_{0}^{D}\geq m>0$ and $\lambda>0$ is sufficiently large, and to pass after some estimations to the limit $p\to\infty$. We leave the details to the reader; see, e.g., [12] for a proof in a related situation.

2. The regularity condition $\Phi(t)\in W^{1,q}(\Omega)$ with $q>d$ is satisfied if $d\leq 3$, $\partial\Omega\in C^{1,1}$, and the Dirichlet and Neumann boundary do not meet, $\Gamma_{D}\cap\overline{\Gamma}_{N}=\emptyset$ [23, Theorem 3.29]. It is also satisfied in up to three space dimensions if $\partial\Omega\in C^{3}$, $\overline{\Gamma}_{D}\cap\overline{\Gamma}_{N}\in C^{3}$, and $\Phi^{D}\in W^{1-1/q,q}(\Gamma_{D})$, $q>d$ [20]. ∎

The paper is organized as follows. The existence theorem is proved in Section 2, while the uniqueness result is shown in Section 3. The numerical solution in one space dimension and its large-time behavior is illustrated in Section 4. The entropy variables $\partial h/\partial u_{i}$ are computed in the Appendix.

2. Existence of solutions

We consider first the nonlinear Poisson equation

[TABLE]

in $\Omega$ with the boundary conditions (4) for given $w_{i}\in L^{\infty}(\Omega)$. Then $(x,\Phi)\mapsto u_{i}(w(x),\Phi)$ is a bounded function with values in $(0,1)$ and a standard fixed-point argument shows that this problem has a weak solution $\Phi\in H^{1}(\Omega)$. Since $\Phi\mapsto u_{i}(w,\Phi)$ is Lipschitz continuous, this solution is unique. By the maximum principle and $f\in L^{\infty}(\Omega)$, we have $\Phi\in L^{\infty}(\Omega)$. Note that $u(w(x),\Phi(x))\in{\mathcal{D}}$ for $x\in\Omega$. Therefore, the following estimate holds:

[TABLE]

where $C>0$ depends on $\lambda$, $z_{i}$, and $\|f\|_{L^{2}(\Omega)}$.

Step 1: Solution to an approximate problem. Let $T>0$, $N\in{\mathbb{N}}$, $\tau=T/N>0$, and $m\in{\mathbb{N}}$ such that $m>d/2$. Then the embedding $H^{m}(\Omega)\hookrightarrow L^{\infty}(\Omega)$ is compact. Let $v^{k-1}:=w^{k-1}-w^{D}\in H^{1}_{D}(\Omega;{\mathbb{R}}^{n})\cap L^{\infty}(\Omega;{\mathbb{R}}^{n})$, $\Phi^{k-1}-\Phi^{D}\in H^{1}_{D}(\Omega)$ be given. If $k=1$, we set $v^{0}=h^{\prime}(u^{0})-w^{D}$ and let $\Phi^{0}$ be the weak solution to $-\lambda^{2}\Delta\Phi^{0}=\sum_{i=1}^{n}z_{i}u_{i}^{0}+f(x)$ in $\Omega$ with boundary conditions (4). Our aim is to find $v^{k}\in H^{1}_{D}(\Omega;{\mathbb{R}}^{n})\cap H^{m}(\Omega;{\mathbb{R}}^{n})$, $\Phi^{k}-\Phi^{D}\in H^{1}_{D}(\Omega)$ such that

[TABLE]

for all $\phi\in H^{1}_{D}(\Omega;{\mathbb{R}}^{n})$ and $\theta\in H^{1}_{D}(\Omega)$. Here, $\alpha=(\alpha_{1},\ldots,\alpha_{n})\in{\mathbb{N}}_{0}^{n}$ is a multi-index, $|\alpha|=\alpha_{1}+\cdots+\alpha_{n}$, $D^{\alpha}=\partial^{|\alpha|}/\partial x_{1}^{\alpha_{1}}\cdots\partial x_{n}^{\alpha_{n}}$ is a partial derivative, and “:” denotes the matrix product with summation over both indices. Since the matrix $B$ is diagonal, we may write the second integral in (15) as

[TABLE]

Lemma 5** (Existence of weak solutions to the time-discrete problem).**

Let the assumptions of Theorem 1 hold and let $w^{D}\in H^{m}(\Omega;{\mathbb{R}}^{n})$. Then there exists a weak solution $v^{k}=w^{k}-w^{D}\in H^{1}_{D}(\Omega;{\mathbb{R}}^{n})\cap H^{m}(\Omega;{\mathbb{R}}^{n})$, $\Phi^{k}-\Phi^{D}\in H^{1}_{D}(\Omega)$ to (15)-(16), and the following discrete entropy production inequality holds:

[TABLE]

where $H$ is defined in (7), $u^{k}=u(w^{k},\Phi^{k})$, $u^{k-1}=u(w^{k-1},\Phi^{k-1})$, and $C_{P}>0$ is the constant of the generalized Poincaré inequality [22, Chap. II.1.4, Formula (1.39)].

Proof.

We employ the Leray-Schauder fixed-point theorem. For this, let $y\in L^{\infty}(\Omega)$ and $\delta\in[0,1]$. Let $\Phi^{k}-\Phi^{D}\in H^{1}_{D}(\Omega)$ be the unique weak solution to the nonlinear problem

[TABLE]

for $\theta\in H^{1}_{D}(\Omega)$. Since $y\in L^{\infty}(\Omega)$, the expression $u_{i}(y+w^{D},\Phi^{k})$ is well-defined. Next, let $X=H^{1}_{D}(\Omega;{\mathbb{R}}^{n})\cap H^{m}(\Omega;{\mathbb{R}}^{n})$ and consider the linear problem

[TABLE]

where

[TABLE]

The bilinear form $a$ and the linear form $F$ are continuous on $X$. Furthermore, using the positive semi-definiteness of the matrix $B$ and the generalized Poincaré inequality with constant $C_{P}>0$ [22, Chap. II.1.4, Formula (1.39)], $a$ is coercive:

[TABLE]

By the lemma of Lax-Milgram, there exists a unique solution $v\in X\subset L^{\infty}(\Omega;{\mathbb{R}}^{n})$ to (18). For later reference, we observe that, since the continuity constant for $F$ does not depend on $y$,

[TABLE]

which gives a bound for $v$ in $H^{m}(\Omega)$ which is independent of $y$ and $\delta$.

This defines the fixed-point operator $S:L^{\infty}(\Omega;{\mathbb{R}}^{n})\times[0,1]\to L^{\infty}(\Omega;{\mathbb{R}}^{n})$, $S(y,\delta)=v$. It clearly holds that $S(y,0)=0$ for all $y\in L^{\infty}(\Omega;{\mathbb{R}}^{n})$. The continuity of $S$ follows from standard arguments; see, e.g., the proof of Lemma 5 in [14]. In view of the compact embedding $H^{m}(\Omega)\hookrightarrow L^{\infty}(\Omega)$, $S$ is also compact. The uniform estimate for all fixed points of $S(\cdot,\delta)$ follows from (19). Thus, by the Leray-Schauder fixed-point theorem, there exists $v^{k}\in X$ such that $S(v^{k},1)=v^{k}$ and $w^{k}:=v^{k}+w^{D}$, $\Phi^{k}$ solve (15)-(16).

It remains to prove inequality (17). To this end, we employ $\tau(w^{k}-w^{D})\in X$ as a test function in the weak formulation of (15). Again, we set $u^{k}=u(w^{k},\Phi^{k})$, $u^{k-1}=u(w^{k-1},\Phi^{k-1})$. Then

[TABLE]

To estimate the first integral, we take $x\in\Omega$ and set

[TABLE]

where we recall that $u_{0}^{D}=1-\sum_{i=1}^{n}u_{i}^{D}$. Then $(\partial g/\partial u_{i})(u)=\log(u_{i}/u_{i}^{D})-\log(u_{0}/u_{0}^{D})$ and $g$ is convex. Hence, $g(u^{k})-g(u^{k-1})\leq g^{\prime}(u^{k})\cdot(u^{k}-u^{k-1})$ or

[TABLE]

Moreover, we infer from the Poisson equation that

[TABLE]

In view of these estimates, the first term in (20) becomes

[TABLE]

We infer from (20) that (17) holds. ∎

Step 2: A priori estimates. Let $(w^{k},\Phi^{k})$ be a weak solution to (15)-(16). Then $u^{k}(x)=u(w^{k}(x),\Phi^{k}(x))\in{\mathcal{D}}$ for $x\in\Omega$, so $(u^{k})$ is bounded uniformly in $(\varepsilon,\tau)$.

Lemma 6** (A priori estimates).**

The following estimates hold:

[TABLE]

where here and in the following, $C>0$ is a generic constant independent of $\varepsilon$ and $\tau$.

Proof.

We need to estimate the second term on the left-hand side of the entropy-production inequality (17). Since $B(w^{k},\Phi^{k})=\mbox{diag}(D_{i}u_{i}^{k}u_{0}^{k})$, we obtain

[TABLE]

where $D_{\rm min}=\min_{i=1,\ldots,n}D_{i}$, $D_{\rm max}=\max_{i=1,\ldots,n}D_{i}$, and we used the fact that $0\leq u_{0}^{k},\,u_{i}^{k}\leq 1$ in $\Omega$. Furthermore, by definition (9) of the entropy variables,

[TABLE]

Inserting these inequalities into (17), it follows that

[TABLE]

We resolve this recursion to find that

[TABLE]

Because of the $H^{1}(\Omega)$ estimate (14) for the electric potential and $\tau k\leq T$, the right-hand side is uniformly bounded. Furthermore, using $\sum_{i=1}^{n}u_{i}^{j}=1-u_{0}^{j}$,

[TABLE]

This finishes the proof. ∎

Step 3: Limit $\varepsilon\to 0$. We cannot perform the simultaneous limit $(\varepsilon,\tau)\to 0$ since we need an Aubin-Lions compactness result, which requires a uniform estimate for the discrete time derivative of the concentrations in $H_{D}^{1}(\Omega;{\mathbb{R}}^{n})^{\prime}$ and not in the larger space $X^{\prime}=(H_{D}^{1}(\Omega;{\mathbb{R}}^{n})\cap H^{m}(\Omega;{\mathbb{R}}^{n}))^{\prime}$. Let $k\in\{1,\ldots,N\}$ be fixed and let $u_{i}^{(\varepsilon)}=u_{i}^{k}$ and $\Phi^{(\varepsilon)}=\Phi^{k}$ be a weak solution to (15)-(16). Set $u_{0}^{(\varepsilon)}=1-\sum_{i=1}^{n}u_{i}^{(\varepsilon)}$. By Lemma 6, there exist subsequences of $(u_{i}^{(\varepsilon)})$ and $(\Phi^{(\varepsilon)})$, which are not relabeled, such that, as $\varepsilon\to 0$,

[TABLE]

We have to pass to the limit $\varepsilon\to 0$ in

[TABLE]

We claim that $u_{i}^{(\varepsilon)}(u_{0}^{(\varepsilon)})^{1/2}\rightharpoonup u_{i}u_{0}^{1/2}$ weakly in $H^{1}(\Omega)$. First, we observe that, because of (24) and (26), $u_{i}^{(\varepsilon)}(u_{0}^{(\varepsilon)})^{1/2}\rightharpoonup u_{i}u_{0}^{1/2}$ weakly in $L^{2}(\Omega)$. Then the claim follows from the bound

[TABLE]

using (22). The compact embedding $H^{1}(\Omega)\hookrightarrow L^{2}(\Omega)$ implies that

[TABLE]

and by the $L^{\infty}(\Omega)$ bounds, this convergence also holds in $L^{p}(\Omega)$ for $p<\infty$. This shows that, taking into account (25),

[TABLE]

In fact, since this sequence is bounded in $L^{2}(\Omega)$, the weak convergence also holds in $L^{2}(\Omega)$. Furthermore, by (26), possibly for a subsequence,

[TABLE]

and this convergence holds also in $L^{2}(\Omega)$.

Then, performing the limit $\varepsilon\to 0$ in (15)-(16) leads to

[TABLE]

for all $\phi=(\phi_{1},\ldots,\phi_{n})\in X$ and $\theta\in H_{D}^{1}(\Omega)$, where $u^{k}:=u$ and $\Phi^{k}:=\Phi$. A density argument shows that we may take $\phi\in H^{1}_{D}(\Omega;{\mathbb{R}}^{n})$.

By the trace theorem, $\Phi^{k}-\Phi^{D}\in H^{1}_{D}(\Omega)$. To show that also $u_{i}^{k}-u_{i}(w^{D},\Phi^{D})\in H_{D}^{1}(\Omega;{\mathbb{R}}^{n})$ holds, we observe that $w^{(\varepsilon)}=w^{D}$ on $\Gamma_{D}$ and therefore, $u_{0}^{(\varepsilon)}=u_{0}^{D}$ on $\Gamma_{D}$ in the sense of traces, where $u_{0}^{D}=1-\sum_{i=1}^{n}u_{i}^{D}$ and $u_{i}^{D}:=u_{i}(w^{D},\Phi^{D})$. Since $u_{i}^{(\varepsilon)}(u_{0}^{(\varepsilon)})^{1/2}=u_{i}^{D}(u_{0}^{D})^{1/2}$ on $\Gamma_{D}$ and $\nabla(u_{i}^{(\varepsilon)}(u_{0}^{(\varepsilon)})^{1/2})\rightharpoonup\nabla(u_{i}u_{0}^{1/2})$ weakly in $L^{2}(\Omega)$ (see (28)), the trace theorem implies that $u_{i}u_{0}^{1/2}=u_{i}^{D}(u_{0}^{D})^{1/2}$ on $\Gamma_{D}$.

In Lemma 5, we have assumed that $w^{D}\in H^{m}(\Omega;{\mathbb{R}}^{n})$ since we have taken $w^{k}-w^{D}\in X$ as a test function. We may take a sequence of functions $(w^{D}_{\delta})$ in $H^{m}(\Omega;{\mathbb{R}}^{n})$ approximating $w^{D}$ and then pass to the limit $\delta\to 0$ to achieve the result for $w^{D}\in H^{1}(\Omega;{\mathbb{R}}^{n})$.

Step 4: Limit $\tau\to 0$. Let $u^{(\tau)}(x,t)=u^{k}(x)$ and $\Phi^{(\tau)}(x,t)=\Phi^{k}(x)$ for $x\in\Omega$ and $t\in((k-1)\tau,k\tau]$, $k=1,\ldots,N$, be piecewise in time constant functions. At time $t=0$, we set $u^{(\tau)}(\cdot,0)=u^{0}$. We introduce the shift operator $(\sigma_{\tau}u^{(\tau)})(\cdot,t)=u^{k-1}$ for $t\in((k-1)\tau,k\tau]$. Then, in view of (29)-(30), $(u^{(\tau)},\Phi^{(\tau)})$ solves

[TABLE]

for all piecewise constant functions $\phi_{i}$, $\theta:(0,T)\to H_{D}^{1}(\Omega)$.

Lemma 6 provides the following uniform bounds:

[TABLE]

where $Q_{T}=\Omega\times(0,T)$ and $C>0$ is independent of $\tau$. Moreover,

[TABLE]

We wish to derive a uniform bound for the discrete time derivative of $(u_{i}^{(\tau)})$. To this end, we estimate

[TABLE]

This holds for all piecewise constant functions $\phi_{i}:(0,T)\to H_{D}^{1}(\Omega)$. By a density argument, we obtain

[TABLE]

Summing these estimates for $i=1,\ldots,n$, we also have

[TABLE]

From these estimates, we conclude that, as $\tau\to 0$, up to a subsequence,

[TABLE]

Taking into account (33) and (36), we can apply the Aubin-Lions lemma in the version of [7] to $(u_{0}^{(\tau)})$ to obtain the existence of a subsequence, which is not relabeled, such that $u_{0}^{(\tau)}\to u_{0}$ strongly in $L^{2}(Q_{T})$, and this convergence even holds in $L^{p}(Q_{T})$ for $p<\infty$. As a consequence,

[TABLE]

Thus, by (33), up to a subsequence,

[TABLE]

We cannot infer the strong convergence of $(u_{i}^{(\tau)})$ because of the degeneracy occurring in estimate (34). The idea is to employ the Aubin-Lions lemma in the “degenerate” version of [3, 14] (also see the Appendix in [15]). In view of (37), the $L^{2}(0,T;H^{1}(\Omega))$ estimates for $(u_{i}^{(\tau)}(u_{0}^{(\tau)})^{1/2})$ and $((u_{0}^{(\tau)})^{1/2})$ (see (33)-(34)), as well as estimate (35), there exists a subsequence (not relabeled) such that

[TABLE]

Taking into account the uniform bound (34), we also have

[TABLE]

This shows that

[TABLE]

weakly in $L^{1}(Q_{T})$. Furthermore, by (37) and (38),

[TABLE]

These convergences allow us to perform the limit $\tau\to 0$ in (31)-(32) to find that $(u_{i},\Phi)$ solves (10)-(11) for all smooth test functions. By a density argument, we may take test functions from $L^{2}(0,T;H^{1}_{D}(\Omega))$. We can show as in Step 3 that the Dirichlet boundary conditions are satisfied, and the initial condition $u_{i}(\cdot,0)=u^{0}_{i}$ in $\Omega$ follows from arguments similar as at the end of the proof of Theorem 2 in [14].

3. Uniqueness of weak solutions

We prove Theorem 3. For this, we proceed in two steps.

Step 1. Adding (1) from $i=1,\ldots,n$ and taking into account the assumptions $D_{i}=1$ and $z_{i}=z$, we find that $u_{0}=1-\sum_{i=1}^{n}u_{i}$ solves

[TABLE]

in $\Omega$, $t>0$, where $W=\sum_{i=1}^{n}W_{i}$, together with the initial conditions $u_{0}(\cdot,0)=1-\sum_{i=1}^{n}u_{i}^{0}$ and boundary conditions (4) and

[TABLE]

We show that this problem has a unique weak solution $(u_{0},\Phi)$ in the class of functions $\Phi\in L^{\infty}(0,T;W^{1,q}(\Omega))$.

Let $(u_{0},\Phi)$ and $(v_{0},\Psi)$ be two weak solutions to (39) with the corresponding initial and boundary conditions such that $\Phi$, $\Psi\in L^{\infty}(0,T;W^{1,q}(\Omega))$. We take $u_{0}-v_{0}$ as a test function in the weak formulation of the difference of (39) satisfied by $u_{0}$ and $v_{0}$, respectively. Then

[TABLE]

The first integral is estimated using the identity $u_{0}(1-u_{0})-v_{0}(1-v_{0})=(1-u_{0}-v_{0})(u_{0}-v_{0})$ and Hölder’s inequality with $1/p+1/q+1/2=1$, where $q>d$ (and $2<p<\infty$ if $d\leq 2$):

[TABLE]

By the Gagliardo-Nirenberg inequality with $\theta=d/2-d/p\in(0,1)$,

[TABLE]

This shows that

[TABLE]

For the remaining integral, we employ the following elliptic estimate

[TABLE]

such that

[TABLE]

Then, inserting the estimates for $I_{1}$ and $I_{2}$ into (40) leads to

[TABLE]

and we conclude with Gronwall’s lemma that $u_{0}=v_{0}$. Consequently, by the Poisson equation in (39), $\Phi=\Psi$.

Step 2. Next, we show that $u_{1},\ldots,u_{n}$ is the unique weak solution to (1), written in the form

[TABLE]

where $F_{i}=u_{0}+\beta z\Phi+W_{i}$, and $(u_{0},\Phi)$ is the unique solution to (39), together with the corresponding initial and boundary conditions. Since we have assumed that $u_{i}\in L^{2}(0,T;H^{1}(\Omega))$, the formulation (1) can be used instead of (10). The classical uniqueness proof requires that $\nabla F_{i}\in L^{\infty}(0,T;L^{q}(\Omega))$; see the first step of this proof. To avoid this condition, we use the entropy method of Gajewski [9, 10].

Let $u=(u_{1},\ldots,u_{n})$ and $v=(v_{1},\ldots,v_{n})$ be two weak solutions to (41) with initial and boundary conditions (3) and (5). We introduce the semimetric

[TABLE]

where $h_{\varepsilon}(s)=(s+\varepsilon)(\log(s+\varepsilon)-1)+1$ for $s\geq 0$. The regularization with $\varepsilon>0$ is needed to avoid that expressions like $\log(u_{i})$ are undefined if $u_{i}=0$. Since $h_{\varepsilon}$ is convex, we have $h_{\varepsilon}(u_{i})+h_{\varepsilon}(v_{i})-2h_{\varepsilon}((u_{i}+v_{i})/2)\geq 0$ in $\Omega$ and hence, $d_{\varepsilon}(u,v)\geq 0$. Now, using (41), we compute, similarly as in [25],

[TABLE]

Rearranging these terms, we arrive at

[TABLE]

Lemma 10 in [25] shows that the first integral is nonnegative. Therefore, integrating the above identity in time and observing that $d_{\varepsilon}(u(0),v(0))=0$, we obtain

[TABLE]

Arguing as in [25, Section 6], the dominated convergence theorem shows that $d_{\varepsilon}(u(t),v(t))$ $\to 0$ as $\varepsilon\to 0$ (here, we use $\nabla F_{i}\in L^{2}(Q_{T})$). Then, since a Taylor expansion of $h_{\varepsilon}$ gives

[TABLE]

we infer that $u_{i}(t)=v_{i}(t)$ in $\Omega$ for $t>0$, $i=1,\ldots,n$, which finishes the proof.

4. Numerical simulations

We illustrate numerically the behavior of the solutions to (1)-(2) for a specific type of ion channel modeled in [11]. First, our numerical scheme is verified by comparing our stationary solutions to the profiles obtained in [4]. Second, we explore the large-time behavior of the numerical solutions.

4.1. Numerical method

The equations are discretized in time by an implicit Euler method and in space by a finite-volume scheme. We suppose that $\Omega=(0,1)$ and impose Dirichlet boundary conditions.

For the finite volume discretization, the domain is divided into uniform cells of size $h>0$. The concentrations and the potential are piecewise constant in each cell with values $u_{i,m}^{k}$ and $\Phi_{m}^{k}$, respectively, where $i=1,\ldots,n$, $m=1,\ldots,M$, at time $k\triangle t$, $k=1,\ldots,K$. These values are determined by the following system of nonlinear equations:

[TABLE]

for $i=1,\ldots,n$, $m=1,\ldots,M$, and $k=1,\ldots,K$. The Dirichlet boundary conditions are accounted for by setting $\Phi^{k}_{0}=\Phi^{D}(0)$ and $\Phi^{k}_{M+1}=\Phi^{D}(1)$, and similarly for the concentrations. Furthermore, we set $f_{m}=\frac{1}{h}\int_{(m-1)h}^{mh}f\,dx$, and the fluxes $J^{k}_{i,m\pm 1/2}$ from cell $m$ to cell $m\pm 1$ are given by

[TABLE]

The concentrations at the cell borders are determined by the logarithmic mean of the cell values:

[TABLE]

for $i=0,\ldots,n$. An advantage of this choice is that the fluxes can be reformulated in terms of the entropy variables

[TABLE]

at least if the concentrations are strictly positive. (We do not use this formulation in the numerical approximation.) The above scheme is implemented using MATLAB, version R2015a. The nonlinear discrete system (42)-(43) is solved by a full Newton method in the variables $u_{i}^{k}$ and $\Phi^{k}$.

4.2. Simulation of a calcium-selective ion channel

We consider a model for an L-type calcium channel described in [11] and used for numerical simulations also in [4]. We choose a simple geometry, where the channel is made of an impermeable cylinder opening up symmetrically into two baths, where Dirichlet boundary conditions are prescribed. For the simulations, three different types of ions are taken into account: calcium (Ca2+, $u_{1}$), sodium (Na+, $u_{2}$), and chloride (Cl-, $u_{3}$). The selectivity filter of the channel consists in eight confined oxygen ions (O*-1/2*), which contribute to the permanent charge density $f=-u_{\text{O}}/2$ as well as to the sum of concentrations in the channel, so that $u_{0}=1-\sum_{i=1}^{3}u_{i}-u_{\text{O}}$. Since these ions are confined, their concentration is assumed to be constant in time. The concentration profile used in our simulations is a simple piecewise constant function, $u_{\text{O}}(x)=0.89$ for $0.45<x<0.55$ and zero else.

In order to obtain results comparable to [4], we use the same one-dimensional approximation of the three-dimensional model that is based on the assumption that the longitudinal extension of the considered domain is much larger than the cross section of the channel. This leads to the reduced system of equations

[TABLE]

where $a(x)$ is the cross-sectional area of the domain at $x\in(0,1)$. It is given by $a(x)=\pi r(x)^{2}$, where the radius $r(x)$ is determined by the piecewise linear function

[TABLE]

For our simulations, we use the parameters given in [4, Section 5.1, Table 1]. The initial concentrations are linear functions connecting the Dirichlet boundary conditions. The initial potential is then computed from the corresponding Poisson equation. The simulations are carried out until the stationary state is reached approximately, which we determine by computing the $L^{2}$ error between the solution at two consecutive time steps:

[TABLE]

The simulation is terminated as soon as $\text{err}_{k}<10^{-13}$. We use the time step size $\triangle t=0.001$ and the mesh size $h=0.01$.

Figure 1 shows the three ion concentrations and the electric potential at various time instances. The scaled concentration values are multiplied by 61.5 mol/liter to obtain physical values. For small times, there is more sodium than calcium present inside the channel region, due to the higher bath and initial concentration of sodium. After some time, the sodium inside the channel is replaced by the stronger positively charged calcium. For higher initial calcium concentrations, the calcium selectivity of the channel acts immediately. The steady-state solution from our simulation coincides with the stationary profile computed in [4, Figure 5], which confirms our numerical scheme. The steady state is reached after 749 time steps, which corresponds to about 23.7 nanoseconds.

4.3. Numerical study of the large-time behavior of the solutions

We investigate numerically the large-time behavior of the solutions and their decay rates to the equilibrium state. First, we consider the setup of the previous subsection. Figure 2 (left) shows the evolution of the relative entropy (7), where the boundary data is replaced by the steady-state solution $(u^{\infty},\Phi^{\infty})$ (see the previous subsection). The right figure displays the $L^{1}$ errors $\|u_{i}^{k}-u_{i}^{\infty}\|_{L^{1}}$ and $\|\Phi^{k}-\Phi^{\infty}\|_{L^{1}}$ versus the number of time steps $k$. We observe that the relative entropy converges exponentially fast to the equilibrium state. By the Csiszár-Kullback inequality (see, e.g., [15] and references therein), the convergence rate in the $L^{1}$ norm is expected to half of that one for the relative entropy, and this is confirmed by Figure 2 (right).

Because of the degeneracy at $u_{0}=0$ in the entropy-production inequality (12), a general proof of exponential convergence rates seems to be not feasible when the solvent concentration $u_{0}$ vanishes locally. Our second numerical example confirms this statement. For this, we choose the oxygen concentration

[TABLE]

All other parameters are kept unchanged. This choice leads to a solvent concentration $u_{0}$ that nearly vanishes in a large part of the computational domain. Consequently, the entropy production in (12) becomes “small” and we may expect a rather slow convergence to equilibrium. Figure 3 illustrates this behavior. After a short initial phase and for the first 20 000 time steps, the convergence rate is very small. This comes from the fact that the values of $u_{0}$ are of the order $10^{-6}$ in the channel region $x\in[0.4,0.6]$, causing the solution to remain nearly unchanged. After about 20 000 time steps, the values of $u_{0}$ increase up to approximately $10^{-3}$ inside the channel region, which initiates the strong exponential decay to equilibrium. These results indicate that exponential decay rates cannot be expected when the solvent concentration vanishes.

Appendix A Entropy variables

The appendix is devoted to a (formal) computation of the entropy variables.

Lemma 7**.**

Let

[TABLE]

Then

[TABLE]

Proof.

It is clear that

[TABLE]

Set $H_{\rm el}(u)=(\beta\lambda^{2}/2)\int_{\Omega}|\nabla\Psi[u]|^{2}dx$, where $\Psi[u]=\Phi-\Phi^{D}$. Recall that $\Phi^{D}$ solves

$-\lambda^{2}\Delta\Phi^{D}=f$ in $\Omega$, $\nabla\Phi^{D}\cdot\nu=0$ on $\Gamma_{N}$. Then $\Psi[u]$ satisfies $-\lambda^{2}\Delta\Psi[u]=\sum_{i=1}^{n}z_{i}u_{i}$ in $\Omega$ together with homogeneous mixed boundary conditions and, by the Poisson equation (2),

[TABLE]

Set $h_{\rm el}(u)=(\beta/2)\sum_{i=1}^{n}z_{i}u_{i}\Psi[u]$. It remains to show that $\partial h_{\rm el}/\partial u_{i}=\beta z_{i}\Psi[u]$. For this, we observe that for any (smooth) functions $u=(u_{i})$, $v=(v_{i})$,

[TABLE]

Let $e_{i}$ be the $i$th unit vector in ${\mathbb{R}}^{n}$ and $w$ be a smooth scalar function. Then, using the linearity of $u\mapsto\Psi[u]$ and (47),

[TABLE]

which shows the claim. ∎