Analysis of improved Nernst – Planck – Poisson models of compressible isothermal electrolytes . Part I : Derivation of the model and survey of the results

We consider an improved Nernst–Planck–Poisson model first proposed by Dreyer et al. in 2013 for compressible isothermal electrolytes in non equilibrium. The model takes into account the elastic deformation of the medium that induces an inherent coupling of mass and momentum transport. The model consists of convection–diffusion–reaction equations for the constituents of the mixture, of the Navier-Stokes equation for the barycentric velocity, and of the Poisson equation for the electrical potential. Due to the principle of mass conservation, cross–diffusion phenomena must occur and the mobility matrix (Onsager matrix) has a kernel. In this paper we establish the existence of a global–in–time weak solution for the full model, allowing for a general structure of the mobility tensor and for chemical reactions with highly non linear rates in the bulk and on the active boundary. We characterise the singular states of the system, showing that the chemical species can vanish only globally in space, and that this phenomenon must be concentrated in a compact set of measure zero in time. With respect to our former study [DDGG16], we also essentially improve the a priori estimates, in particular concerning the relative chemical potentials.


Introduction
Increasing the efficiency of actual high-performance energy storage systems requires an exact understanding of their fundamental physical principles. Of particular interest is ion transport in electrolytes, for instance in lithium-ion batteries. Classically this transport is modelled by the Nernst-Planck theory, but it has an important drawback. In the neighbourhood of interfaces it is failing for various reasons (see [DGM13,DGL14]): First of all, the classical Nernst-Planck model neglects the high pressures induced by the Lorentz force that affect the charge transport. Secondly, it does not take into account the interaction between the solvent and the charged constituents. A third drawback of the Nernst-Planck theory is the widely used assumption of local charge neutrality. This assumption completely fails in the vicinity of the boundaries where electric charges accumulate. In the new model, the dissipative diffusion fluxes are not proportional to the concentration gradients like in the Nernst-Planck theory. The reaction rates do not exhibit the structure of products of monomials in the concentrations. The natural estimates available for the model are thus completely different than the ones known for standard diffusion-reaction systems. In this paper, we show that the original modelling yields new PDE structures and progresses in the analysis of models of mathematical electro chemistry. We establish the existence of a global-in-time weak solution for the model, allowing for 1 A general structure of the mobility tensor, with cross-diffusion effects and a nontrivial kernel; 2 Chemical reactions of arbitrary growth rate in the bulk and on an active boundary that represents a one-sided electrochemical interface; 3 A mechanical contribution (pressure) to the diffusion flux; 4 Different reference specific volumes of the constituents (possible solvatisation).
Each of these points represents an absolute breakthrough in the analysis of models in mathematical electro chemistry. In particular, the original thermodynamic modelling for the reaction rates yields a control in Orlicz classes so that global-in-time weak solutions can be defined without the help of the advanced re-normalisation techniques (see [Fis15]) necessary to handle models traditionally called 'of mass action type' (see a. o. [CDF14,HK16a,HK16b,Fis15,HMPW16,BFPR16,MHM15,HHMM16,Mie16]).
In addition to the existence statements, we are able to characterise the singularities of the system associated with the vanishing of species. We show that except for the occurrence of a complete vacuum -which is entirely non physical in the range of validity of the model -the mass density of the species can vanish only globally in the spatial domain and that this phenomenon is concentrated in a compact set of measure zero in time. To the best of our knowledge, this important property is not yet known in the context of Nernst-Planck type models.
Our method relies on the one hand on a priori estimates that result from the thermodynamically consistent modelling, and from the conservation of total mass. The estimates are partly a consequence of known results for the Poisson equation or the Navier-Stokes equations, but we can regard the estimates on the chemical potentials of the mixture constituents, in particular in the presence of chemical reactions, as original. The idea is to make use of the diffusion gradients and the entropy production of the chemical reactions to reveal a relationship between the blow up of differences of chemical potentials and the entire vanishing of certain groups of species. In order to exclude the latter phenonemon, a restriction on the total initial masses of the involved constituents turns out sufficient.
in our study: Particularly the constitutive equations for the pressure, for the diffusion fluxes and for the reaction terms, are different in [MPZ15] and in [DGM13]. The compactness question occurs like in our analysis but is solved assuming a special structure of the viscosity tensor, called Bresch-Desjardins condition. This allows to obtain estimates on the density gradient, a device which is not at our disposal here. A further difference between the two mixture models concerns cross-diffusion, which is described in [MPZ15] and [Zat15] by a special non symmetric choice of the mobility matrix, whereas we allow for general symmetric positive semi definite matrices. Note that the mobility matrix must be symmetric at least in a binary mixture. Among recent less directly related investigations let us mention: In the context of general diffusion, [Bot11], [JS13]; for models with simplified diffusion and pressure laws [FPT08], [BFPR16]; for the analysis of incompressible models of Nernst-Planck-Poisson type [BFS16], [CJ15], [HMPW16].
Due to the length of the investigation, we shall split its presentation into three contributions. In this paper, we derive the model, motivate its reformulation, and propose a survey of the main mathematical results with the fundamental steps of the proofs. The complete proofs of these results follows in [DDGG17a] devoted to the construction of approximate solutions with a priori estimates, and in [DDGG17b] devoted to the investigation of the compactness and the proof of convergence of the approximation scheme. The papers [DDGG17a,DDGG17b] are naturally more technical than the present one, but in our eyes they possess independent interest.
In this paper, in the Section 2, the model will be introduced following [DGM15]. The model is formulated for the normal regime of the system, i.e. it is assumed that the mass densities of the constituents do not vanish. For the mathematical analysis we will derive an equivalent formulation which exhibits more stability against extreme behaviour. The equivalent formulation of the model is the object of Section 6. We then formulate our main results and sketch the proofs.

Improved Nernst-Planck-Poisson model
Throughout the paper, the bounded domain Ω ⊂ R 3 is representing an electrolytic mixture. The boundary of Ω possesses a disjoint decomposition ∂Ω = Γ ∪ Σ: The surface Γ represents an active interface between an electrode and the electrolyte, where chemical reactions and adsorption may occur. The other surface Σ represents an inert outer wall with no reactions and no adsorption.
The compressible mixture consists of N ∈ N species denoted by A 1 , . . . , A N . A species A i is a carrier of electric charge, z i ∈ Z, and of molecular mass, m i ∈ R + .
We assume that the system is isothermal, so that the absolute temperature denoted by θ is a positive constant. Under the isothermal assumption the thermodynamic state of the mixture at time t ∈ [0, T ] is described by the mass densities ρ 1 , . . . , ρ N of the species, the barycentric velocity v of the mixture and the electric field E. As usual in electro chemistry, a quasi-static approximation of the electric field is considered, i.e. the magnetic field is constant and the electric field satisfies E = −∇φ. The scalar function φ is called electrical potential.
The active boundary Γ can be viewed as a mixture of N Γ = N + N ext constituents denoted by A 1 , . . . , A N Γ , where the additional N ext constituents take into account the species of the adjacent exterior matter, i.e. electrode species. Thus we only consider surface chemical reactions with participating species that also exist in the adjacent bulk domains. The surface constituents have the surface mass densities ρ Γ 1 , . . . , ρ Γ Γ, respectively. The kth chemical reaction in the bulk (k ∈ {1, · · · , s}) and on the boundary (k ∈ {1, · · · , s Γ }) possesses the general structure The constants a k i , b k i and a k Γ,i , b k Γ,i are positive integers called stoichiometric coefficients. For k = 1, . . . , s, we define a vectorial coefficient associated with the kth bulk reaction via (1) Due to the inclusion of the molecular mass, γ 1 , . . . , γ s are not as commonly the stoichiometric vectors, but this will simplify the notation. The forward reaction rate of the kth reaction is R f k > 0, and the backward reaction rate is rate R b k > 0. The net reaction rate of the kth reaction is defined as The same definitions hold for the surface reactions on Γ. Here the vectorial coefficients are defined and the surface reaction rates are R Γ k = R Γ,f k − R Γ,b k for k = 1, . . . , s Γ . Since charge and mass are conserved in every single reaction (3)

Balance equations in the bulk
In the isothermal case the evolution of the thermodynamic state is described by the equations of partial mass balances, of momentum balance, and by the Poisson equation.
In ]0, T [×Ω the mixture obeys for i = 1, . . . , N the partial mass balances In (4), v denotes the barycentric velocity of the mixture. The quantities J 1 , . . . , J N are called the diffusion fluxes. We use upper indices in their case because they are vector fields of R 3 and not scalars. The mass production r i of constituent A i is related to the reaction rates via The total mass is defined as = N i=1 ρ i . Together with v it has to satisfy the continuity equation Thus the conservation of total mass yields additional constraints on the diffusion fluxes and the mass productions: While the constraint (7) 2 is a consequence of (5) and the conservation of mass in every chemical reaction (3), the side-condition on the diffusion fluxes has to be guaranteed by an appropriate constitutive modelling.
The principle momentum balance possesses the expression Herein σ denotes the Cauchy stress tensor, b is the force density due to gravitation, and the symbol n F E stands for the Lorentz force due to the electric field. The quantity n F represents the free charge density that is defined via Throughout the paper, we are going to neglect the gravitational force that plays no role in the analysis.
In the electrostatic setting the balance equation for the electric field reduces to the Poisson equation for the electrical potential, Here χ > 0 is the constant susceptibility of the electrolyte.

Constitutive equations
The constitutive equations for the mass fluxes, the reaction rates and the stress tensor can be derived from a single free energy density ψ of a general form ψ = h(θ, ρ 1 , . . . , ρ N ) .
The derivatives of the free energy function with respect to the mass densities are called chemical potentials, In the isothermal setting, the balance equations and the free energy density yield a local entropy production ξ = ξ D +ξ R +ξ V ≥ 0 with three contributions due to diffusion, ξ D , reaction, ξ R , and viscosity, ξ V (see [MR59,BD15,DGM15]). A constitutive model that relies on the free energy function (11) implies explicit expressions for the three entropy productions as binary products. From these expressions we may derive constitutive equations that yield three separate non-negative entropy productions. For more details regarding the derivation of the entropy production we refer to [MR59,dM63,BD15]. In [BD15] its is shown how cross-effects revealing the Onsager symmetry can be introduced.

Diffusion fluxes
The entropy production due to diffusion reads where D 1 , . . . , D N are the thermodynamic driving forces for diffusion, The simplest constitutive ansatz for the diffusion fluxes J 1 , . . . , J N that implies ξ D ≥ 0 is given by The proportionality factor M ∈ R N ×N sym is called the mobility matrix. It is positive semi definite and may depend on ρ. Moreover, the side- For instance, following the paper [DGM13], one can construct M from an empirical mobility matrix M emp (ρ) and a linear operator P : where d 1 , . . . , d N −1 > 0 are diffusion constants, and the lines of the matrix P are given by the differences e i − e N of standard basis vectors for i = 1, . . . , N . In fact, any operator P that satisfies for i = 1, . . . , N the condition N j=1 P i,j = 0 can be chosen in (15) in order that (14) is valid. However our analytical results do not rely on the particular structure (15) of the matrix M .

Reaction rates
The entropy production due to chemical reactions assumes the form The driving forces D R 1 , . . . , D R s are defined via To achieve ξ R ≥ 0, we assume that the vector of production rates are derived from a convex, nonnegative potential This choice is in fact more general than in [DGM15], where the following potential is employed, with positive constants α 1 , . . . , α s and constants β 1 , . . . β s ∈]0, 1[, C ∈ R arbitrary. By this choice Dreyer et al. achieve for the rates an ansatz of Arrhenius type, which is widely used in chemistry, Stress tensor The entropy production due to viscosity is represented via ,...,3 , and Id denotes the identity matrix. We split the Cauchy stress tensor into a viscous part S visc and the pressure p via σ = −p Id + S visc . Then the material pressure can be calculated from the free energy function (11). The resulting representation is called Gibbs-Duhem equation and reads The simplest constitutive choice for the viscous stress tensor S visc satisfying ξ V ≥ 0 describes a Newtonian fluid. It reads where η > 0 is the coefficient of shear viscosity, and the coefficient λ of bulk viscosity satisfies λ + 2 3 η ≥ 0.

Choice of the free energy function
The constitutive model is derived from a free energy density of the general form (11). However, for the analysis of the model, we need in some extent to specify the choice of the free energy function. To this end the free energy density ψ is additively split into three contributions, The constants µ ref i (i = 1, . . . , N ) are related to the reference states of the pure constituents. The contribution h mech is the mechanical part of the free energy that is neglected in the classical Nernst-Planck theory. The function h mix represents the mixing entropy.
In the presentation of [DGM13,DGL14], the contributions h mech and h mix are naturally given as functions of the number densities (molar concentrations) n 1 , . . . , n N of the constituents. These are defined via n i := ρ i /m i (i = 1, . . . , N ). The number fractions (concentration fractions) are defined via y i := n i N j=1 n j for i = 1, . . . , N .
The mechanical free energy is associated with the isotropic elastic deformation of the mixture. It takes into account different reference specific volumes V 1 , . . . , V N ∈ R + of the constituents. Assuming a constant bulk compression modulus K > 0 the mechanical free energy according to [DGL14] is a function of the mixture volume For the sake of generality, we express h mech in the form Dreyer et al. use F (x) := x ln x + C 1 for an simple mixture, whereas the Tait equation corresponds The contribution h mix results from the entropy of mixing and is given by where k B is the Boltzmann constant.

The model for the boundary Γ
The active boundary Γ represents an interface between the electrolyte and an external material. In the most important application the external material is an electrode which is likewise a mixture of N ext ∈ N constituents. Here we have analogous quantities to those that occur in the electrolyte, namely the barycentric velocity, and diffusion fluxes and so on. To distinguish between the electrolyte and the external material we provide the external quantities the suffix ext .
In this paper we assume for simplicity that the constituents occurring on Γ also exist in the bulk, either in the electrolyte or in the external material. Thus the interface Γ is a mixture of N Γ = N + N ext constituents.
The equations of an interface representing a surface mixture are derivable in the context of surface thermodynamics, and we refer the interested reader to [ABM75, DGM15,Guh14]. As in the bulk there are universal surface balance equation and material depending surface constitutive equations.
To simplify the surface equations we assume on ]0, T [×Γ: (a) Time variations of the surface mass densities and tangential transport are negligible in comparison to mass transfer across the surface and to chemical surface reactions. Then the surface balance equations become stationary.
(b) The interface is fixed in space, i.e. the interfacial normal speed is zero.
(c) There is no velocity slip and the normal barycentric velocity is equal to the interfacial normal speed,

Surface mass balances and surface reaction rates.
In what follows the interfacial unit normal ν points into the external material. Under the assumptions above, the surface mass balance equations Here we make use of the convention that the N first species on Γ are the electrolyte constituents, while the constituents with indices N + 1, . . . , N + N ext are the external ones.
It remains to specify the surface mass production r Γ due to surface reactions. As in the bulk, the production r Γ is related to the surface reaction rates R Γ by The interfacial entropy production ξ Γ R due to chemical reaction is (see [DGM15]) The entropy production of the surface has the same structure as the corresponding entropy production in the bulk (16). Thus in order to satisfy the entropy inequality an ansatz similar to (17), (18) may be used. We assume the existence of a potential Ψ Γ so that Here, µ Γ 1 , . . . , µ Γ N are the surface chemical potentials of the electrolytic species, whereas the quantities µ Γ N +1 , . . . , µ Γ N +N ext are associated with the external species. These equations describe the adsorption of a constituent from the bulk to the surface. The kinetics of this process is controlled by positive semi definite matrices Simpler form of the transmission conditions. In the general thermodynamic setting, the surface chemical potentials are derivatives of a surface free energy. Due to the assumption of stationary surface equations, and that the boundary is fixed, we are able to formulate all surface equations in terms of the bulk chemical potentials. From a mathematical viewpoint the equation system (24), (27), (28) only serves to eliminate the surface chemical potentials µ Γ in order to calculate the external fluxes of the electrolytic species.
Following the Appendix of the paper [DDGG16], we can introduce linear spaces via (Image and it turns out that the interactions on the boundary Γ take place tangentially to the linear subspace We defineŝ Γ := dim U and we choose basis vectorsγ 1 , . . . ,γŝ Γ for U. In the paper [DDGG16], it is proved that the conditions (24), (27), (28) allow to represent the flux of the electrolytic species via Moreover, the external response J 0 is a mapping from [0, T ] × Γ into U. The modified reaction term r is a mapping from [0, T ] × Γ × Rŝ Γ into U and possesses the structurê Moreover, there is a convex nonegative potentialΨ Γ : The data µ ext , M Γ and M Γ,ext are absorbed in the position dependence of J 0 andΨ Γ . In the analysis we can for simplicity assume that J 0 andΨ Γ are the boundary data, even if in applications the solution of a few additional nonlinear algebraic equations are necessary to compute them.
Electrical potential. The boundary condition for the electrical potential can be derived from Maxwell's equations for surfaces, which are satisfied in the quasi-static stetting by a continuous electrical poten- where φ 0 is the given electric potential at ]0, T [×Γ.

Summary model equations
Domain Ω. Summarising, the evolution of the state (ρ, v, ϕ) in ]0, T [×Ω is described by the PDE- Here n F is given by the formula (9) while the fluxes J 1 , . . . , J N obey (13). The reaction rates R 1 , . . . , R s obey (17) where the external chemical potentials µ ext , the external potential φ 0 and the kinetic matrices M Γ and M Γ,ext are given. The reaction rates r Γ obey (25) together with (38) and (39), withr of the structure (31) and J 0 a given mapping into U.
Initial conditions. Initial conditions are prescribed for the variables ρ 1 , . . . , ρ N . We denote them ρ 0 i , i = 1, . . . , N . Moreover, an initial state v 0 is also given for the velocity vector.

Notation
To get rid of overstressed indexing, we simplify the notation by making use of vectors. For instance we denote ρ the vector of mass densities, n the vector of number densities i.e.
Moreover we define the vector 1 := 1 N := (1, 1, . . . , 1) ∈ R N , and the vectors of quotients of charge and mass, and of volume and mass z m Using these conventions, we have a. o. the identities The diffusion fluxes J 1 , . . . , J N span a rectangular matrix J = {J i j } ∈ R N × R 3 . The upper index corresponds to the lines of this matrix. Vectors of R N are multiplicated from the left, as for instance in 1 · J = N i=1 J i which is an identity in R 3 . The vectors γ 1 , . . . , γ s span a rectangular matrix γ = {γ k i } ∈ R s × R N . The upper index corresponds to the line of the matrix. Vectors of R s are multiplicated from the left, as for instance in the identity r = R · γ = s k=1 R k γ k in R N . Analogously the vectors γ 1 Γ , . . . , γ s Γ Γ span a rectangular In order to describe the reactions, we shall further make use of the abbreviationsR : R s → R s ,R := −∇Ψ andR Γ : R s Γ → R s Γ ,R Γ := −∇Ψ Γ .
For a convex, nonnegative potential Ψ ∈ C(R s ), s ≥ 1 (and for its conjugate Ψ * ), the vectorial Orlicz classes L Ψ (Q; R s ) and L Ψ * (Q; R s ) are well known. We make use of the notation ForΨ Γ ∈ L ∞ (S; C 2 (Rŝ Γ )), we define a vectorial Orlicz class LΨ Γ (S; Rŝ Γ ) as the set of all mea-

Mathematical assumptions on the data
From the thermodynamic viewpoint, the free energy function h, the mobility matrix M and the potentials Ψ, Ψ Γ are the essential objects determining the properties of the constitutive equations. Mathematical results can be obtained under suitable restrictions to these objects. In addition, restrictions are as usual necessary concerning the geometry and the quality of boundary and initial states. In this investigation we are not concerned with pointing at intrisic classes for the second type of data.
Assumptions on the free energy function. Our estimates on the (relative) chemical potentials moreover require the special form (21), where the mixing entropy obeys the precise representation (23). We allow for a certain generality only at the level of the function h mech which we assume of the form (22). Recall that K is a positive constant, of which we will keep track due to its importance.
We assume that F belongs to C 2 (R + ) ∩ C(R 0,+ ) and is a convex function. We assume that there are 3 2 < α < +∞ and constants 0 < c 0 , c 1 such that In the neighbourhood of zero, we assume that F (s) behaves like s ln s: There are constants positive constants k 0 < k 1 and s 0 > 0 such that In fact, in order to obtain an unconstrained PDE system, we crucially need that F : R + → R is surjective. This not satisfied for instance by the pure polynomial ansatz according to Tait (proof in the paper [DDGG17a]), but it always follows from (42).
Assumptions on the mobility matrix. The mobility matrix M is symmetric and positive semi definite. Throughout the paper, we assume that M is mass conservative, that is Moreover we assume that the entries of M are linear-growth, continuous functions of the vector ρ of the partial mass densities.
Except for these few points, the exact structure of the mobility matrix is a delicate topic (in particular there are connections to the Maxwell-Stefan theory, see [BD15]). In this paper we restrict ourselves to the assumption that M has rank N − 1 independently on ρ. In other words, denoting 0 = λ 1 (M ) ≤ λ 2 (M ) ≤ . . . ≤ λ N (M ) the eigenvalues of the matrix M , we assume that there are positive constants 0 < λ ≤ λ such that λ ≤ λ i (M (ρ)) ≤ λ (1 + |ρ|) for all i = 2, 3, . . . , N, ρ ∈ R N + .
Let us remark that due to this assumption, only regularisations of the original ansatz of the paper [DGM13] are included in the analysis: In the formula (15) we must, for example, apply a cut-off from below to the entries of the empirical matrix M emp .
Assumptions on the reaction rates. We assume that the reaction rates are derived from a strictly convex, non-negative potential 1 Ψ ∈ C 2 (R s ). Moreover, Ψ satisfies Similarly, we require that the boundary reaction rates are derived from a strictly convex, non-negative For simplicity we explicitly require at least linear growth of the reaction rates via It is well known that r(Ω, Γ) > 2 in general (see [Grö89] a. o.), but there are numerous situations where, depending on the boundary of the domain and the structure of the surface Γ, the optimal exponent satisfies r(Ω, Γ) > 3 (see [DKR15] for results and discussions on this topic). We require with α from (41). This of course might be a restriction only if α < 2.
1 It is always possible to achieve the non negativity because the modelling only requires that Ψ has a global minimum at zero Assumptions on the remaining boundary data. We consider only non degenerate initial and boundary data. This means that Moreover we assume as a compatibility condition the validity in the weak sense of Note that the last assumption in (51) garanties that J 0 ∈ L ∞ (]0, T [×Γ; U) (see [DDGG16], Appendix for the elementary proof).
The selection S will be called uncritical if dim(W S ) = |S| and critical otherwise.
is the finite union of sub manifolds of dimension at most N − 1. We say that the initial compatibility condition is satisfied if the initial vector of the total massesρ 0 := Ω ρ 0 dx ∈ R N + satisfiesρ 0 ∈ M crit .

Identification of natural variables in the equations of mass transfer
In this section we formulate three essential remarks that affect the solution concept and the mathematical analysis of the equations (32) of mass transfer taken independently: There are state-constraints (ρ ≥ 0) on the solution vector; The mobility matrix has a kernel (M 1 = 0) so that the structure of the system is not entirely parabolic; In the context of weak solutions, the continuity equation ∂ t + div( v) = 0 might generate a local vacuum that affects the regularity of the solution.

State-constraints
The pair of vector fields (ρ, µ) : [0, T ] × Ω → R N + × R N is subject to the algebraic relation µ = ∇ ρ h(ρ) (cf. (12)). Obviously: The vector of mass densities ρ must belong to the domain of the gradient of the free energy function, while the vector of chemical potentials belongs to its image.
Meaningful choices of the function h must in general guaranty that the domain of ∇h is R N + . Indeed, the algebraic constraint on ρ must be equivalent with the physical non negativity restriction on the mass densities. The vector of chemical potentials µ must belong to the image of ∇h. There are models, for instance the special constitutive assumption (21) in the case of the Tait equation, in which this image is a true subset of R N .
These algebraic state-constraints are a fundamental obstacle to the application of a functional analytic method to prove the resolvability of the model.
In order to overcome this difficulty we are going to exploit a particular observation: For the special constitutive assumption (21), we can show that ∇h : R N + → R N is a bijection if the first derivative of the function F is surjective onto R. At least for a relevant particular choice of h, the PDE system is unconstrained in µ, and the chemical potentials are a favourable set of variables to perform existence theory or numerical approximation. This was already noted in the context of less complex models, for instance in [JS13].

An 'hyperbolic' component
Diffusion and chemical reactions are the dissipative structure that provide a control on the vector µ. But the fluxes J 1 , . . . , J N and the functions r 1 , . . . , r N occurring in the system (32) in fact only depend on the projection of the vector µ on the subspace 1 ⊥ := {ξ ∈ R N : ξ · 1 = 0} (see the side conditions (14) for the diffusion flux, and to the restriction (3) on the vectors γ 1 , . . . , γ s ).
Thus, natural estimate can be obtained only for a (N − 1)−dimensional projection of the vector µ.
Due to this observation, it was noted in [DGM13] that a change of variables is necessary in order to define the solution. We keep as main variables: (a) On the one hand, one coordinate of the vector field ρ, namely the total mass density ρ · 1 that we shall denote throughout the paper. This is the 'hyperbolic' component subject to the continuity equation (6); (b) On the other hand, N −1 coordinates of the vector of chemical potentials µ defined via a projection onto the linear space 1 ⊥ ⊂ R N .
The possibility of these choices relies on an algebraic result that we want to afore mention here.

q) .
A proof of the elementary result of Proposition 6.1 is given in the second part of the article [DDGG17a].
We here note two additional remarks concerning the choice of the projector Π.

Remark 6.2.
It might be of importance to allow for a general Π, as its choice can be suited to the structure of the mobility matrix (e. g. (15)) in order to simplify the structure of the diffusion.

Vacuum oscillations
Although the presuppositions of the free energy model (21) in fact completely fail if the total mass density is below a lower critical value, the mathematical analysis cannot exclude the occurrence of a complete vacuum.
From the viewpoint of our analytical treatment, a vacuum is characterised by the fact that the variables ρ and q are 'decoupled', in the sense that the mapping q → R( = 0, q) is trivial on the entire R N −1 .
A concrete technical difficulty is raised concerning the compactness. Time-derivative estimates are available only for the ρ−variables and do not transfer one to one to the q variables, since a sequence of mass densities ρ n = R( n , q n ) (n ∈ N) such that n → 0 would converge strongly even if the corresponding q n exhibit oscillatory behaviour.
Since the reaction densities are non-linear expressions of q 1 , . . . , q N −1 , the vacuum-oscillations affect the concept of the solution at this level: The representation or the production rates R = −∇Ψ(γ 1 · Eq, . . . , γ s · Eq) is restricted to the set where is strictly positive. An analogous situation occurs at the boundary ]0, T [×Γ whenever it is in contact with a vacuum.
In order to include the possibility of this extreme behaviour, we relax the concept of a solution to (32), For the representation of the bulk reactions, we require r = s k=1 γ k R k with the following weaker condition: R =R(γ 1 · Eq, . . . , γ s · Eq) in Q + ( ) . For the concept of the solution, we ask for the representationr = ŝ Γ k=1γ k R Γ k together with R Γ =R Γ (t, x,γ 1 · Eq, . . . ,γŝ Γ · Eq) in S + ( ) . (57)

The mathematical results
For t > 0, we denote Q t :=]0, t[×Ω the space-time cylinder, and if T > 0 is the final time of the process, we abbreviate Q := Q T . We denote S t :=]0, t[×Γ and S = S T .

The solution class
Exploiting the preliminary considerations of the Section 6, a solution vector to the entire system (32), An essential property of solutions is the mass and energy conservation.
Definition 7.1. We say that ( , q, v, φ, R, R Γ ) satisfies the (global) energy (in)equality with free energy function h and mobility matrix M if and only if the associated fields and variables (58) satisfy for almost all t ∈]0, T [ We say that ( , q, v, φ, R, R Γ ) satisfies the balance of total partial masses if the vector field is subject toρ The conservation of energy provides the natural bounds that will allow to define the weak solution. We introduce what one could call a natural class B, because this class naturally arises from the global energy and mass conservation identities associated with the model. The class B reflects the regularity of the weak solution and essentially depends on several parameters The final time T > 0, the domain Ω and the partition Γ ∪ Σ of its boundary (see the condition (49)); The choice of the free energy function h and in particular the growth exponent of (41); The mobility matrix M , in particular the number rk M ; The choice of the potentials Ψ and Ψ Γ for the reaction densities.
The variables , φ and v will satisfy the conditions with the exponents α > 3/2 and r(Ω, Γ) > 2 of the conditions (41) and (49), and with β := min r(Ω, Γ), 3α For the variables R and R Γ we consider the conditions For the variable q, a control is achieved on the spatial gradient thanks to the assumption (44). However in the context of flux boundary conditions, the bound on the L 1 norm is a non trivial problem. We shall nevertheless obtain the integrability in time via complex estimates involving the diffusion gradient, the reactions and the conservation of total mass. Under the assumption (47) of at least quadratic potentials, a natural class for the variable q is then q ∈ L 2 (Q; R N −1 ) , ∇q ∈ L 2 (Q; R (N −1)×3 ) .
The natural class B also encodes an information concerning the conservation of global mass (integration of (4), (32) over Ω). We additionally introduce a non-negative function Φ * ∈ C([0, T ] 2 ), Φ * (t, t) = 0 constructed from the functions Ψ, Ψ Γ (and thus from R and R Γ ) via for all 0 ≤ t 1 ≤ t 2 ≤ T . Here C 0 is an appropriate constant that we will choose later. For a function u ∈ C 1 ([0, T ]), we define a weighted modulus of uniform continuity via We are finally in the position to introduce the solution class.
The concept of weak solution is well defined owing to standard estimates.
Assume moreover that the vectorρ 0 of the total initial masses has positive distance to the manifold M crit of (53).
We can characterise the singular states of the system associated with species vanishing. (2) V t = {x ∈ Ω : (t, x) = 0} is a set occupied by a complete vacuum; (3) λ 3 (N t ) = 0.
If one starts with total initial masses on the critical manifold, then it is possible that certain groups of species are completely consumed after finite time, and the solution then exists only up to this time. Afterwards, it might be necessary to restart the system with a smaller number of species. Then, there are a time 0 < T 0 depending only on the data and a time T 0 ≤ T * ≤ +∞ such that there is a weak solution ( , q, v, φ, R, R Γ ) ∈ B(t, Ω, α, N −1, Ψ, Ψ Γ ) in the sense of Definition 7.3 to (P t ) for all t < T * . Moreover the following alternative concerning T * is valid: (1) Either T * = +∞; (2) Or there is ( , q, v, φ, R, R Γ ) ∈ B(T * , Ω, α, N − 1, Ψ, Ψ Γ ) that weakly solves (P t ) for all t < T * , such that min i=1,...,Nρi (t) > 0 for all t ∈ [0, T * [, and such that

Main steps of the proof
For the existence proof, we first regularise the problem in order to obtain equations that are easier to solve. It is important that the regularisation is thermodynamically consistent: it must preserve the fundamental dissipation mechanisms and the natural estimates. For the larger class of (regularised) problems, we are then able to derive the energy identity and the global balance of partial masses in a rigorous way. This entails a priori estimates in the natural class. The a priori estimates for the variable q is one of the most original part of the analysis. In order to pass to the limit with the numerous nonlinearities of the system, it is at last necessary to derive compactness statements: This is the second supporting pillar of the analysis.

The approximation scheme
The regularisation strategy, though not mass conservative, will be chosen thermodynamically consistent, since it consists in two essential steps: (1) A positive definite regularisation of the mobility matrix M ; (2) A convex regularisation of the free energy function h.
For δ, σ, τ > 0, we call weak solution to the problem (P τ, σ, δ ) a vector (µ, v, φ) subject to the energy inequality and such that the quantities satisfy the identities (70), (72), and instead of (71) Since the second definition in (74) in general implies that N i=1 J i = 0, it is necessary to add this term in the momentum equation (75) in order to preserve the energy identity.
The existence of approximate solutions for the regularised scheme is not a trivial problem, because there is no monotone or pseudo-monotone structure inherent to the diffusion. We carry over this technical step in the paper [DDGG17b] by means of a time-continuous Galerkin method. It turns out that for parameters δ, σ, τ > 0, weak solutions to (P τ,σ,δ ) exist and develop no vacuum.

The a priori bounds
Via standard techniques (mainly the Gronwall Lemma), we obtain from the energy identity natural uniform bounds, in particular From the equation (72), we obtain via standard elliptic theory (3−α) + }. Moreover, if β ≥ α a control on the Lorentz force is achieved: From the Navier-Stokes theory (for (71), (75) we obtain improved bound on the pressure: If 3/2 < α ≤ 3, r(Ω, Γ) > α and 1 · J ≡ 0, then p A major challenge is the proof of the following estimate (to find in [DDGG17a]). We give here the formal argument.