Existence theory for a Poisson-Nernst-Planck model of electrophoresis

A system modeling the electrophoretic motion of a charged rigid macromolecule immersed in a incompressible ionized fluid is considered. The ionic concentration is governing by the Nernst-Planck equation coupled with the Poisson equation for the electrostatic potential, Navier-Stokes and Newtonian equations for the fluid and the macromolecule dynamics, respectively. A local in time existence result for suitable weak solutions is established, following the approach of Desjardins and Esteban [Comm. Partial Diff. Eq., 25 (2000), 1399--1414].


Introduction
Electrophoresis is the motion of charged colloidal particles or molecules through a solution under the action of an applied electrical field. This phenomenon is one of the most powerful analytical tools in colloidal science, being often used in the characterization of colloidal systems [34]. Modern applications include drug delivery and screening, manipulation of particles in micro-/nanofluidic systems, sequencing of genome of the organisms, forensic analysis, micro-chip design and others [8], [32], [20]. Despite the fact that the analytical description of this type of electrokinetic phenomenon is a difficulty task (due to the complicated interplay between hydrodynamic and electric effects [41]), there is a vast literature on the subject. The majority dealing with special problems like electrophoresis of bio-molecules with some geometric symmetry property or small surface potentials (see, for example, [1], [19], [24], [35], [38]). Few of these papers deal with the rigorous mathematical aspects of the models involved.
In this paper we present a model of electrophoresis of a single particle immersed in a viscous and incompressible fluid (a ionic water solution). Unlike the models where the concentrations of charged particles range from colloidal to nano size (see [36]), we deal with the case where the colloidal particle is a very large rigid molecule and obeys the laws of classical mechanics. Also we consider the situation where the particle and fluid are contained in a bounded domain O ⊂ R 3 which represents the enclosure of the system (a tube, in the capillary electrophoresis). These are fairly reasonable assumptions in a number of practical situations such us in the study of electrophoresis of proteins [12]. On the other hand, we assume that the particle is far from the boundary of O. As a consequence we are able to analyze only the local motion of the particle, as we are considering the standard electrokinetics equations. In fact, in the standard models (see [41], [1], [24]), the analysis is restrict to a single particle that is suspended in a fluid which fills all the region in R 3 exterior to the particle. This means that the boundaries are sufficiently far removed from the particle so that they have negligible effects on the electric and fluid velocity fields associated with the particle [23]. According to [37] when the geometry forces the mobile particle into close of proximity with a surface (which may itself be charged) complex phenomena arise and a more detailed analysis should be done.
Although we are interesting mostly in the local motion of the particle, we do not neglect the inertial terms in the Navier-Stokes equations for fluid velocity and pressure. The Reynolds number for fluid flows in a typical electrophoretic motion is small and many authors consider the Stokes equations instead the Navier-Stokes equations (see [41], [32], [3]). However, it is not clear to us under what circumstances the linear theory can be used. With respect to this issue Allison and Stigter [2] suggest that the use of linearization of the constitutive equations with respect to the perturbing electric/flow fields are allowed provided the perturbing fields are weak. Also we believe that the analysis of the non simplified model considered here can be given us a support from the treatment of more general problems where the inertial effects could not be negligible (see [19], [8]).
Another important aspect in the modeling of electrophoretic motion is the electrostatic potential. In more simplified models the total electrostatic potential of the system is given by the applied external electric field (see [3], [38]) and the local potential (due to particle charge distribution and ions) is neglect. On the other hand, the local electrostatic potential is usually modeled by the Poisson-Boltzman equation [5], [4], [7], [42]. This equation is derived from the assumption of thermodynamics equilibrium where the ionic distributions are not affected by fluid flows. According to [30], this is a reasonable assumption for the case of steady-state processes with stationary values of diffusive fluxes, but there are some important cases where convective transport of ions has nontrivial effects [27]. In this context, in the lines of [23] and [41], we consider a more convenient approach: the use of a convective-diffusion type equation, more precisely, the well known Nernst-Planck equation, for the ionic concentration and Poisson equation for the total electric potential. Moreover, we only impose C 2+α regularity on the domain occupied by the particle (without other geometric symmetry restrictions).
Some of the mathematical aspects of a electrophoresis model based on the Poisson-Boltzmann theory have been discussed by the authors in the papers [4] and [5]. However, at least to the knowledge of the authors, there are no rigorous mathematical results on the Poisson-Nernst-Planck model of electrophoresis of the rigid macromolecules. From the mathematical point of view, the main feature relies in the highly coupled equations which must be solved in time dependent domains. We establish the proof of a existence theorem (see Theorem 3.1 in Section 3) for appropriate weak solutions for this system. The proof is based on the approximation technique introduced in the references [14], [15] for the study of the motion of rigid bodies in viscous fluid. A sequence of approximate smooth solutions is construct and the solution is obtained as the limit of this sequence. Our scheme is local in time in the sense that the macromolecule does not touch the boundary of the enclosure in the time interval of existence. Employing a linearization and regularization procedure, this permits us to obtain suitable energy bounds for the electric potential, fluid velocity and ionic concentration (see Lemmas 4.2,4.3,4.5). As a consequence, we establish special results on time compactness (see Lemmas 5.1 and 5.2) that give us strong convergence results for fluid velocity and ionic concentration. Then, the main result follows from the passage to the limit in the approximation equations (see Lemmas 5.3,5.4 and 5.5) and is announced in Theorem 1.
The outline of the paper is as follows. In Section 2 we introduce the coupled system of governing equations; in the Section 3 we obtain formal energy estimates that give us a weak formulation for the system; in the Section 4 we construct a sequence of approximate solutions that are uniformly bounded in energy spaces; in the Section 5 we obtain appropriate compactness results that permit us to pass the limit the approximate solutions.
In a subsequent work we will investigate the effect of proximal boundaries in the electrophoresis of many particles and consider the possible contact between two particles or with the boundary, in the lines of [18].

Governing Equations
We assume that the solvent and macromolecule occupy a bounded open connected set O ⊂ R 3 and the solvent contains J ionic species. At the initial time, the macromolecule occupies a compact region K 0 , where K 0 ⊂⊂ O is a open connected set such that The fluid domain at the initial time is denoted by Ω 0 := O\K 0 . We also assume that K 0 and Ω 0 are C 2+α -domains, 0 < α < 1, and define K t = Q(t)K 0 as the position of the particle at time t, Q(t) is an affine isometry such that Q(0) = I. Also, let us denote Ω t = O\K t , i.e., the fluid domain at time t and assume (formally) that there exists T > 0 such that It is important to remark that, although we have not imposed any symmetry hypothesis on K 0 , the regularity imposed is not satisfactory in a large number of situations. In fact, the determination of the electrostatic interface in biological processes is not a trivial task. According to [11] (see also [12]), any model of molecular surface must follow, in some way, the boundary of the space from which other molecules are excluded. Usually, the molecular surface (the accessibility surface of the solvent) is taken to be the surface described by a point on the surface of an idealized spherical solvent probe as it is rolled around the solvated molecule. This lead us to a C 0 surface with the eventual presence of cusps and with multiple connected components. Accordingly, we are dealing with a idealized situation in order to obtain appropriate regularity results and a priori estimates for the approximate solutions of the system. However, we believe that this condition can be relaxed using a smoothness domain technique. We will discuss this issue in a subsequent paper.
If the dielectric constant of the particle and of the solvent are κ 1 > 0 and κ 2 > 0, respectively, the electrostatic potential ψ, for each t ∈ [0, T ], is governed by the Poisson equation e is the electron charge, Z i and N i are the valence and concentration of the ith ionic species, respectively. Moreover, where ρ 0 is a fixed charge distribution on K 0 such that suppρ 0 ⊂⊂ K 0 . Note that (2.5) implies that the fixed charges of the particle are invariant under rigid motion. As a consequence, we have the conservation of total fixed charges for all t ∈ [0, T ]. These assumptions on ρ corresponds, of course, to an idealized situation. In fact, the determination of the properties of fixed charges in colloidal molecules is a very difficulty task (see [12], [28]).
We assume transmission boundary conditions for ψ where ν(., t) is the exterior normal to ∂K t and Ψ represents a stationary electrostatic potential on ∂O, ψ 1 = ψ| Kt , ψ 2 = ψ| Ωt . The evolution of N i , for each i ∈ {1, . . . , J}, is described by the Nernst-Planck equation where d i is the ionic diffusion coefficient of the ith ionic species, v f is the solvent velocity field, κ B is the Boltzmann constant and θ is the temperature of the system (which we suppose constant). It is important to remark that, for obvious physical reasons, we seek non-negative functions The boundary conditions correspond to no ion diffusion and no ion conduction through the boundaries (see [23]) are given by The system is complemented with the initial condition The solvent velocity field is governed by the Navier-Stokes equations for incompressible flows is the electrical forcing term, p is the pressure, η > 0 and µ f > 0 are the viscosity and the mass density of the fluid, respectively. Let denote by x c (t), v c (t) and w(t) the center of mass, the translational velocity and the rotation vector of the particle in the time t, respectively. If A, v p and µ p > 0 are the symmetric inertial matrix, velocity and the mass density of the particle, we have and Using (2.3) and (2.14), the electrical forcing term F can be written in the di- (see [41]). If M is the mass of the particle, the evolution law for its motion is given by the Newtonian dynamic equations where σ H is the stress tensor of the fluid.
We assume the homogeneous Dirichlet condition for v f on ∂O and require the velocity and normal stress to be continuous at the interface between the solid body and fluid where Σ is the Cauchy stress tensor in the solid. Also, we have the following initial conditions for the velocities In the "zeta" potential approach for electrophoresis modeling, (2.20) is replaced by a "slip" boundary condition on v f : a nonlinear Dirichlet condition that depends on ψ and it is based on a Prandtl boundary layer approximation (see [3], [4]). This approximation is relatively accurate when the Debye screening length of the macromolecule is much smaller comparable with its radius of curvature [38]. Otherwise, in the more realistic cases, this approach is no longer correct [23].

Weak Formulation and Main Result
We can obtain a weak formulation for the problem (2.3)-(2.22) if we take into account the energy framework of the system. First, we need to introduce a global formulation for (2.13), (2.17)-(2.18). Following [14], we define in O × [0, T ] the Eulerian densities the global density µ = µ p + µ f and the global velocity In view of mass conservation, µ must satisfies where we have consider the balance the momentum in the fluid and particle (see [14]). Summing the above equations and using (2.21) we obtain the global formulation for (2.13), (2.17)-(2.18), and that Formally, taking the inner product of (3.24) with u, integrating by parts, using (3.25) and (3.23) we obtain, for all t ∈ (0, T ), (3.27) is the kinetic energy at time t, is the viscous dissipation and is the total work of the electrical force. Then (3.27) corresponds to the energy conservation of the system. Moreover, integrating (3.27) with respect to t we obtain the energy (mechanical) bound (3.28) As a consequence, recalling that there exists a positive constant it is easy to obtain the estimate where µ min = min{2η, µ p , µ f } and µ max = max{µ p , µ f }. Now, supposing that Ψ ∈ H 1 (∂O), we define If we take the product of (2.3) with ψ − Ψ, where Ψ ∈ L Ψ , we obtain, using (2.7) and (2.8), for each t ∈ [0, T ]. Using Cauchy and Young inequalities, we have for each t ∈ [0, T ]. The left side of (3.31) corresponds to the electrostatic energy contribution (the negative of the free energy according to [33]) which we denote by E el ; the first term in E el corresponds to the self-energy of the electric field and the next two terms are the electrostatic energies of the ions and fixed charges (see [12], [28], [33]). As a consequence of (3.31), for known N i and u and each t ∈ [0, T ], for all Ψ ∈ L Ψ . Note that (2.10) (as well as (3.23)) corresponds to a conservation type equation. In fact, from the transport theorem d dt then, integrating (2.10) in Ω t × [0, t], using Gauss theorem, (2.11), (2.19) and (2.12), we obtain the molar conservation for the ith ionic species, (3.33) In particular, using (2.6), we have the total charge conservation for the system. From the above discussion, in particular, from (3.32) and (3.27), we have a natural weak formulation of the above system (see Definition 3.1 below), which is obtained from the following minimization problem with constraints: Find (µ, u, N 1 , . . . , N J , ψ) (in appropriate functional spaces) that minimizes the energy E in O × (0, T ) and such that ψ(., t) minimizes E el [µ, u, N 1 , . . . , N J , .](t) in L Ψ , for each t ∈ [0, T ]; (2.10), (2.11)-(2.12), (3.25)-(3.26) are the constraints; furthermore, the pressure p and the Cauchy stress tensor Σ are the Lagrange multipliers of the problem.
The additional hypotheses on the data are Despite the fact that the ions are concentrated very close to ∂K t , the hypotheses on the support of ρ 0 and N i,0 are reasonable as there is an ion exclusion region close to the boundary of the molecule [12]. In order to obtain the suitable functional spaces of the solutions, we observe that (3.28) and (3.29) can give us H 1 -bounds for u and ψ. However, we need uniform L 2 -estimates for F and for N i . As we will see, for small T (depending on the initial data) uniform L 4 (Ω t )-estimates for ∇ψ and N i can be obtained; as a consequence, we can estimate F and the conduction term in the weak formulation for (2.10).
We need to define the following functional spaces where Definition 3.1. Let us to assume that K 0 and Ω 0 are C 2+α -domains, 0 < α < 1, ρ is invariant under rigid motion (see (2.5)), hypothesis H and (2.1). Then for all w ∈ S and a.e. t ∈ (0, T ), where The main result of the paper is given below As remarked in the introduction of the present work, the strategy of the proof follows the lines of the results in references [14] and [15] for the motion of rigid bodies in incompressible fluids. It is based on the solution of an approximated system, which give us an approximate sequence of solutions (u n , µ n , ψ n , N 1,n , . . . , N J,n ) for (3.34)-(3.38). The solution (u, µ, ψ, N 1 , . . . , N J ) is built up as a limit of these approximations.
The new feature here is the coupling between (3.34)-(3.38) which lead us to consider an appropriate fixed-point schema in order to construct the approximate sequence. The existence results for the Lagrangian version of the parabolic problem (3.38) depend on the suitable estimates on the term containing ∇ψ 2 (see Lemma 4.3). Then, the crucial point of the proof is the obtention of the uniform L 4 -estimates for ∇ψ 2 (see Lemma 4.1 below), in the sense that the bound does not depend on the particle motion. Using results on singular integral operator and elliptic regularity, this is possible if we assume that (2.2) is valid.

Approximate solutions
Let us fix 0 < γ < γ 0 . We begin with a technical result Lemma 4.1. Suppose t > 0 and Q(t) an affine isometry such that Q(0) = I and K t = Q(t)K 0 , Ω t = O\K t satisfy (2.2). Consider the situation in which is replaced by f (., t) ∈ L 2 (O) and Ψ is replaced by g ∈ H 1 (∂O). Then the problem (3.37) has a unique solution ψ(., t) ∈ H 1 (O) that satisfies Proof. From a theorem due to Stampacchia (see [6]) it is a routine to check that there exists an unique solution in H 1 (O) for the problem (3.37). Let us to extend f to be zero outside of O.
and there exists C = C(O) such that For the above results see Theorem 9.9 in [21].
and consider the following problem (4.43) From well known results (see Theorem 2.2 in Chapter 4 in [29]), this problem has a unique solution Φ(., where C = C(O, κ 2 , γ) and we have extended Φ to be zero in O\B γ . Now, let us to consider 0 < σ < γ/4 and define the set Note that, from (2.2), we have dist(A σ (t), B γ ) > γ/4. We also define ω(ξ, t) = ω(Q(t)ξ, t), ∀ξ ∈ K 0 and consider the following auxiliary problem This problem has a unique solution where we have extended φ to be zero outside A σ (0) ∪ K 0 . In fact this follows from suitable modifications in the arguments of the Theorem 7.2 in the reference [40] (see also [22]). Furthermore, from the results on singular integral operators of [10] and [40] (for related results see [4]), there exists a constant C > 0, depending on K 0 , κ 1 , κ 2 and γ such that This follows from an elementary calculation using the properties of Q(t). We claim that (4.50) In fact, the transmission boundary conditions (2.7)-(2.9) are satisfied in the sense of traces, in particular, (2.7) implies that ψ(., t) ∈ H 1 (O); also it it easy to check that (3.37) is valid. The result then follows from the uniqueness of the problem (3.37). Now, note that, it is easy to check that Henceforth, from (4.51), standard Sobolev embedding, (4.48), the trace theorem, (4.52) and (4.42), we have where C = C(K 0 , O). Similarly, Consequently, (4.39) follows from (4.50), (4.45), the trace theorem, Sobolev embedding, (4.42), (4.53) and (4.54).
In what follows we describe how to construct the sequence of approximations. We adapt the proof in [15]. The idea is to introduce an approximation scheme which consists in solving a system of regularized equations: an "almost" linear problem related to the velocity field as well as appropriate linear problems for N i and ψ. The existence of these regularized solutions is obtained from the use of a fixed point type theorem. As described below, this is done in small time intervals, chosen in a such way that the advecting vector field is close to the identity, (2.2) is valid in each time interval, and using a Lagrangian Galerkin type method for the linear problems related to u and N i (a similar approach may be seen in the study of the free surface problems [39]).
Let us consider where d min = min{1/2, d 1 /8, . . . , d J /8} and From the Sobolev embedding theorem, there exists a positive constant Let us define L 1 > 0 such that and Following [15] we observe that it is possible to show the existence of a homeomorphism Θ from the space of incompressible vector fields and which corresponds to a rigid motions in K t into the representation space Let us also define W T = L ∞ (0, T ; L 2 (Ω 0 )). Here we will consider the natural corresponding norms in Y T and W T . Now, using standard regularization techniques, we consider Let us also consider, for f ∈ L 2 (0, T ), where {g ǫ } ǫ>0 is a family of regularizing kernels with respect to time such that where R ǫ is a regularization operator (see [14]), such that R ǫ (v) is analytical in time and smooth in space (in particular, Lipschitz in space), 3 , as ǫ → 0 + . Furthermore, v ǫ is divergence free in O and corresponds to a rigid motion in the particle domains Also, we denote, for each t ∈ [0, T ], Ω ǫ,t = O\K ǫ,t . From the construction of Θ (see [15] for ǫ > 0 sufficiently small. Then, let consider ǫ ′ > 0 such that (4.60), (4.64), (4.62) and are valid for all 0 < ǫ < ǫ ′ . As a consequence of (4.64) and (4.59), it is easy to check that The inequality (4.67) implies that, in each time interval (t k−1 , t k ), the advecting vector field is close to the identity. This is an important point in the proof and it is based on [15]. As we will see, from (4.67), the Lagrangian forms of certain operators that appear in the Lagrangian formulation for (3.37) and (3.38) are uniformly elliptic. In fact, from (4.67), it is easy to check that, for each m ∈ {0, 1, . . . , N } and t ∈ [mt 0 , (m + 1)t 0 ], As a consequence, we have existence results for the related linearized problems, in each time interval (t k−1 , t k ).
Let us consider L 2 , L 1 , T , ǫ ′ > 0, N and t 0 as above. For each 0 < ǫ < ǫ ′ and in each time interval (t k−1 , t k ), we want to solve a set of linearized problems and to apply the Schauder's fixed point theorem. Below we give the details for the first time interval (0, t 0 ).
A direct calculation using (4.68) gives us that the problem (4.76) is uniformly parabolic. Using the regularity of the coefficients we can obtain the existence of a unique nonnegative solution N i,ǫ ∈ C 2+α,1 (Ω 0 × [0, t 0 ]) of (4.76)-(4.77) (see Theorem 5.3, Chapter IV in [25]) that corresponds to the unique solution of (4.78). For sufficiently small t 0 , we can obtain uniform estimates on N i,ǫ . In fact, we have the following Lemma.
where L 2 is given in (4.55).

Fixed Point Procedure
Let us recall that we have chosen (v, ϑ 1 , . . . , ϑ J ) ∈ B × X × . . . × X . Now we observe that for (x, t) ∈ K ǫ,t × (0, t 0 ), x ǫ (t) is the center of mass of K ǫ,t . Let us define the map t 0 w ǫ (τ )dτ and Θ −1 2 (y ǫ ) denotes the incompressible component in Ω 0 of Θ −1 (y ǫ ); as described in [15], y ǫ is a divergence-free vector field constructed from u ǫ , v ǫ and from suitable rigid current functions. In order to show that this map has a fixed point in B × X × . . . × X we need to prove some technical results.
Proposition 4.1. The operator G ǫ has a fixed point in B × X × . . . × X .
Proof. First we observe that G ǫ is a compact operator. In fact, this follows from Lemma 4.6, the rigidity properties in K 0 and from the Sobolev embeddings in Ω 0 for N i,ǫ (as a consequence of the estimates (4.79) and (4.85)) and u ǫ (from the estimates (4.96) and (4.97)). Then, the map G ǫ has a fixed point in B, using Lemma 4.6 and Schauder's theorem.
The first time step on (0, t 0 ) is completed. From interpolation results (see Theorem 1.3.8 in [31]), u ǫ (., t 0 ) ∈ H 1 (O) 3 . Therefore, we can proceed similarly on the interval (t 0 , t 1 ), considering (u ǫ (., t 0 ), N ǫ (., t 0 )) as the initial conditions for the problems (4.95) and (4.78), respectively. We give below some details related to the obtention of the bound (4.79) in this second step of the procedure.
Let us consider ψ ǫ and N i,ǫ the solutions of the problems P1 and P2 in (0, t 0 ), respectively. We extend these functions in (t 0 , t 1 ) to be the solutions of P1 and P2 in this interval. First, we observe that (4.73) is valid in (t 0 , t 1 ) using the Lemma 4.1. From a similar argument as in (4.80), we have, for all t ∈ [t 0 , t 1 ], and, for all t ∈ [0, t 0 ], (4.108) Considering t = t 0 in (4.108), we obtain the validity of inequality (4.108) for all t ∈ [0, t 1 ]. As a consequence, we obtain a version of (4.79) in the interval (0, t 1 ), following the lines of the proof of Lemma 4.3. In particular, (4.79) is valid if we consider the interval (t 0 , t 1 ). A similar argument can be used in order to obtain the estimate (4.96) in (0, t 1 ). Repeating all the process in each step (t i−1 , t i ) (recall that the number of steps N = T /t 0 only depends on ǫ and L 1 ) we obtain, for each ǫ > 0, ( u ǫ , N ǫ , ψ ǫ ) as the solutions, respectively, of the problems (4.95), (4.78) and (3.37) in (0, T ) and µ ǫ that satisfies (4.71). If we consider u ǫ and N i,ǫ in terms of the Eulerian coordinates, we see that, for each 0 < ǫ < ǫ ′ , (u ǫ , µ ǫ , N 1,ǫ , . . . , N J,ǫ , ψ ǫ ) is a solution of the problems P1, P2 and P3 and we can obtain an approximate sequence (u n , µ n , N 1,n , . . . , N J,n , ψ n ) n∈N of solutions for (3.34)-(3.38). Using  In order to establish convergence properties for N n it is convenient to consider its extension to fixed domains. So we consider the standard extension operator E : H 1 (Ω 0 ) → H 1 (R 3 ) such that Ef = f a. e. in Ω 0 , ∇(Ef ) 2 0,2,R 3 ≤ C ∇f 2 0,2,Ω0 , Ef 2 0,2,R 3 ≤ C f 2 0,2,Ω0 and Ef ≥ 0 a. e. in R 3 as f ≥ 0 a. e. in Ω 0 ; the constant C depends only on ∂K 0 and ∂O (see Chapter 5 in [17]). Then, from (4.109) and (4.68), we have ess sup where N * i,n (x, t) = N * i,n (X n (x, t, 0), t) and N * i,n = E N i,n .
It is a routine to check that (µ, u) satisfies (3.35) and (3.36) and u satisfies (4.109). In particular, from (4.111), this implies that (2.2) is valid for the domains K t and Ω t that correspond to u. Below we establish some compactness results.
Proof. Following [13], we have a scalar version of (5.116). More precisely, (5.116) is valid for functions in H 1 (R 3 ). In this case, Let us consider i ∈ {1, . . . , J} fixed and denote Then, for each n ∈ N, we bound by where, Using (4.112) and Lemma 1 of the reference [13], we have where C = C(L 2 , C, T ). Now, extending R n (u n ) to be zero in R 3 \O, we see that I n (x, t) = I Ωn,t (x) is the solution of the transport problem ∂ t I n + ∇ · (I n R n (u n )) = 0, I n (., 0) = I Ω0 .
in D ′ ((0, T ) × R 3 ). Then, using (5.118), we have where C = C(L 2 , L 1 , T, C). As a consequence, Denoting the first an the second terms in the right side above by P ′ 3 and P ′′ 3 , respectively, we have  (3.38). Then, we have, where we have extended [ψ n ] n (., τ ) to be zero if τ > T . Subtracting the above expressions, we obtain, after some calculation, Hence,   In particular this implies that N * i ≥ 0 a. e. in R 3 , as it is easy to check. Now, let us consider the problem (3.37) with the data (K t , Ω t , N * 1 , . . . , N * J , ρ, Ψ). From Lemma 4.2 this problem has a unique solution ψ(., t) ∈ H 1 (O) that satisfies the bound (4.73). As a consequence we can prove the following convergence result for ψ n (., t).
We have obtained T > 0 and a weak solution for the problem in (0, T ). From the bound (4.109), Q(t) is continuous in [0, T ] (see [14]), so that lim t→T − dist(K t , ∂O) ≥ γ > 0 and we can iterate the existence result in order to obtain Theorem 3.1.