The linearized Poisson-Nernst-Planck system as heat flow on the interval under non-local boundary conditions

The linearized of the Poisson-Nernst-Planck (PNP) equation under closed ends around a neutral state is studied. It is reduced to a damped heat equation under non-local boundary conditions, which leads to a stochastic interpretation of the linearized equation as a Brownian particle which jump and is reflected, at Poisson distributed time, to one of the end points of the channel, with a probability which is proportional to its distance from this end point. An explicit expansion of the heat kernel reveals the eigenvalues and eigenstates of both the PNP equation and its adjoint. For this, we take advantage of the representation of the resulvent operator and recover the heat kernel by applying the inverse Laplace transform.


Introduction
The Poisson -Nernst-Planck (PNP) system [7] is a fundamental model for electrodiffusion and is one of the main tools in modeling ion channels in cell membrane(see, e.g.[12]).In one of its simplest forms, it contains a pair of drift-diffusion equations for positive and negatively charged ions, coupled with the equation for the electric field induced by the charges.
We concentrate on the case of two types of ions (positively (C + ) and negatively (C − ) charged) and a closed channel, where the flux of C ± is zero at the ends of the channel, hence the number of ions of each type (and, in particular, the total charge C + − C − ) is preserved in time.
The physical model behind the PNP is a drift diffusion for charged particles, where the diffusion is due to independent Brownian motions of the ions, and the drift is due to the external field induces by the potential difference between the ends of the channel, and the mean electric field generated by the moving ions.
Thus, the PNP can be considered as a system of Kolmogorov forward equation, [9,10], whose solutions represent the probability distribution of a test particle for each type of ions in the system.
In the case of zero external field, the neutral case (C + = C − ) induces a steady, uniform distribution for both charges.A linearization of this equation around this constant neutral case is reduced, up to a re-scaling of the time, into the naive looking damped diffusion equation [11,12] 1 for the local charge u ≈ C + − C − : where the interval [0, 1] is the channel, u(x, t) is the local charge at x ∈ [0, 1], t ≥ 0 and κ 2 is the inverse Debye screening length.This looks like a fairly naive equation.However, the boundary conditions are where the electric field E is given by the Poisson equation driven by the voltage difference V across the end points x = 0, 1.These are non-local boundary conditions.Indeed, we show that (2, 3) can be reduced to the following The steady state for the linearized problem can easily be obtained: so we can subtract it from the solution u of equation to get a homogeneous boundary conditions 1 I wish to thank Dr. Doron Elad for turning my attention to this formulation Some versions of parabolic equations under non-local boundary conditions were studied by several authors (see, e.g.[1,6]).A stochastic interpretation of linear diffusion equations under non-local boundary conditions goes back to Feller [5].In that paper he extended his seminal paper [4] to non-local boundary conditions, and interpreted the diffusion equation in terms of a Brownian particle which may undergo a jump from a point on the boundary of the interval to a distributed position at the interior.This extension was later studied by several authors, see e.g [6,8,13].However, in all these cases the process is allowed to jump from a boundary point to the interior, and not the other way around.This will be the case if, e.g., κ2 is replaced by −κ 2 in (5).The boundary conditions ( 5) associated with the operator d 2 /dx 2 − κ 2 suggests a diffusion process which jump at a random Poisson time of mean κ 2 from an inner point x ∈ (0, 1) and reflected at the endpoint x = 0 with probability 1 − x, and at the endpoint x = 1 with probability x.In this sense, it is a forward Kolmogorov equation representing the evolution of a probability distribution of the charge. 2 The heat kernel K(x, y, t) of such an equation generates the solutions This kernel represents the probability of the particle to be at position x at time t + s, conditioned that it was at point y at time s.In particular K(x, y, t) ≥ 0, for (x, y) ∈ (0, 1) × (0, 1), t > 0, and lim t↓0 K(x, y, t) = δ x−y .The adjoint equation, then, represents the backward Kolmogorov equation of the process.In general, it is also a diffusion equation of the same form and adapted boundary condition, whose kernel K * is given by the interchange of x and y: K * (x, y, t) = K(y, x, t).However, in the case of b.c (5), an explicit form of the adjoint operator is not clear.
In this paper we attempt to calculate the spectral expansion of the heat kernel.The information encoded in this expansion contains the eigenvalues, as well as the eigenfunctions of both the operator and its adjoint.
To obtain this, we take advantage on the explicit solutions of the resulvent R = R(λ, x, y) where λ ∈ C, (x, y) ∈ (0, 1) × (0, 1): where δ x−y is the Dirac delta function and R satisfies the boundary conditions (5) in the x variable.These solutions can be expressed locally as a combination of the trigonometric functions sin(λ 1/2 x) and cos(λ 1/2 x).It turns out that the solution of the resulvent equation exists whenever Re(λ) < −κ 2 .This resulvent R can also be written in terms of the heat kernel K (see [3], and also the review in the Appendix): It turns out that R is a meromorphic function of λ in the complex plane, analytic if Re(λ) < −κ 2 , and admits a countable number of simple poles in the half plane Re(λ Under some conditions which we can verify (c.f.Appendix) we can recover the heat kernel form (6) using the inverse Laplace transform via where Γ is a contour enclosing all the poles of R.Then, we use the Residue theorem [2] to evaluate the contour integral.
The main results are summarized below: Theorem 1.1.The heat kernel of (1, 5) is given by where λ n are the roots of Det is given by (15) and In particular, is the heat kernel for the operator L 0 = d 2 /dx 2 on the domain (5 ), where • ψ (2) are the eigenstates of L 0 , and are the eigenstates of its adjoint L † 0 .This poses a challenging question regarding the formulation of this problem, since (except of the constant), φ n are not solutions of φ xx + λφ = 0 for any λ ∈ C. It seems that the adjoint operator may not be given by a differential one, and the non locality of the boundary conditions leaks into the operator itself.(c.f Section 5).

The linearized PNP system
The one dimensional PNP equation takes the form [11] , where C ± is the concentration of positive/negative ions, and E is the electric field given in terms of the concentrations C ± and the potential difference V : The special case of non penetrating charges corresponds to zero flux on the boundary In the neutral case C + = C − = η and E = 0. We linearize this system : and ignore all terms of second order in Ec ± to obtain Here we concentrate in the case B = 0 which reduces to a single equation on u.Without loss of generality we also assume D = 1: Lemma 2.1.The three b.c ( 9), together with the constraint −ǫE x = u can be introduced as a pair of non-local conditions: Proof.By the field equation and the boundary condition (9, 5) admits a classical C 2 solution and u(x, 0) ≥ 0 and so 3 Properties of the linrarized PNP We start from the following We start from the following c.For any t > 0, u(•, t) ∞ ≤ cosh(κ/2) u(•, 0) ∞ .

Eigenvalues and eigenfunctions
The eigenfunctions of the operator d 2 /dx 2 are given by a sin(λ 1/2 x)+b cos(λ 1/2 x).Substitute this in (12) we get The system (13, 14) is a linear system for the coefficients a, b.The determinant of this system is Lemma 3.1.λ 1/2 Det(λ) is a meromorphic function on the complex plane.The roots of Det(λ) = 0 are given by (2kπ) 2 where k ∈ Z.In addition λ m , m ∈ N where {λ m } are the roots of In addition, λ = 0 is a root of Det(λ) only if κ 2 = 12, and it coincides with a root λ 1 (κ) of ( 16) as κ 2 → 12.
In addition, λ = −κ 2 is the only negative root of Det, and it is a simple one.
Proof.The case of non-zero roots follows directly from (15).
At this stage we can only show that there exists R(κ) > 0 such that all roots of (16) outside the disc {|z| < R(κ)} are real and simple.Numerical test verifies, for all selected values of κ, that all roots inside the disc are real and simple as well. in the interior of the square whose boundary is Γ n .Since g is a polynomial of order 3, the number of zeroes of g inside the square is 3 for any n large enough.Since all the poles of h are identical to the poles of f , which are given by (±k + 1/2)π on the real line, k ∈ N ∪ {0}, there are exactly two poles of h in inter(Γ k+1 ) − inter(Γ k ) for all k large enough, namely at x = (±k + 1/2)π.Thus, there are exactly two zeroes (or one zero of order 2) in inter(Γ k+1 ) − inter(Γ k ).Evidently, there are two real roots in this domain, one in each interval (kπ, (k + 1)π) and (−(k + 1)π, −kπ).Thus these are the only roots in inter(Γ k+1 ) − inter(Γ k ).It follows, in particular, that h has only real roots outside a large enough square.
• If κ 2 = 12 then λ 1 = 0 and ψ (2) Proof.The proof follows by Lemma 3.1 and (13, 14).If κ 2 = 12 then 0 is not an eigenvalue, even though it is a root of Det.The reason is that the coefficients of (13,14) are degenerate in that case.However, if κ 2 = 12 then the first root λ 1 of ( 16) is zero, and the eigenfunction follows by substitution.

The resolvent
The resolvent operator for Neumann problem on [0, 1] is expressed in terms of the eigenvalues and eigenfunctions of the operator: where The resolvent R corresponding to the boundary condition ( 12) can be written as From the boundary conditions of ( 12) and ( 17) we obtain We can now solve (19,20) for any λ = 0 which is not a root of Det, a(λ, y) = κ 2 Det −1 (λ) After some trigonometric manipulations on (21, 15) we obtain 4 The heat kernel To obtain the heat kernel corresponding to the equation (12) we use (33) to obtain K(x, y, t) = (23) We now recall that the second term above is just the heat kernel of the Neumann problem.This can be expanded in eigenfunctions: Then we calculate the residues of (R − R N )e −λt using (22).The residue at the pole λ = −κ 2 due to the second term in (22) is Let us now evaluate the other poles of ( 22).The poles of the first term (the coefficients of sin(λ 1/2 (x − 1/2)) are originated by two sources: Since sin(λ 1/2 /2)/Det(λ) has no singularity at λ = (2mπ) 2 , the only singularity due to Det(λ) are the roots of (16), i.e at λ = λ m .The residue Theorem at this singularities yield where A(λ, y) as given in (18).However, the first term of (22) contains also the poles at λ = (2k + 1)π due to the singularity of A(•, y) at these points.

k
(y) = A(λ k , y) are not trigonometric functions.In particular we cannot identify L † 0 with d 2 /dx 2 on a certain domain D † , as we did for L 0 .Open question: Find an explicit expression for the generator to the adjoint operator L † 0 and its domain D † .t = e −γt δ t=0 = δ t=0 as a distribution.The heat kernel can, then, be written as K(x, y, t) e −λt R(λ, x, y)dλ (34)