Bulk-Surface Electrothermodynamics and Applications to Electrochemistry

We propose a modeling framework for magnetizable, polarizable, elastic, viscous, heat conducting, reactive mixtures in contact with interfaces. To this end, we first introduce bulk and surface balance equations that contain several constitutive quantities. For further modeling of the constitutive quantities, we formulate constitutive principles. They are based on an axiomatic introduction of the entropy principle and the postulation of Galilean symmetry. We apply the proposed formalism to derive constitutive relations in a rather abstract setting. For illustration of the developed procedure, we state an explicit isothermal material model for liquid electrolyte|metal electrode interfaces in terms of free energy densities in the bulk and on the surface. Finally, we give a survey of recent advancements in the understanding of electrochemical interfaces that were based on this model.


Introduction
The energy transition from fossil fuels to renewable energy sources gives rise to an increasing demand for more efficient energy storage in stationary and mobile applications, cf. [1]. As current lithium-ion battery technology is expected to reach its theoretical limit soon, several promising future technologies like metal-sulfur and metal-air systems, as well as polymer electrolyte batteries and different other types of solid state batteries, are intensively investigated [2][3][4][5].
Because the demand for both portable water and clean water for industry and agriculture is expected to grow, there is a strong need for better water desalination and water purification technologies. Electrically driven processes like electrodialysis and capacitive deionization are investigated since they promise to be commercially competitive and in particular energy efficient, compared to thermal or pressure driven processes, [6][7][8].
A key ingredient for better understanding of the electrochemical processes in the mentioned technologies is the development of better mathematical models. Standard models like the Poisson-Nernst-Planck system, cf. e.g., [9,10], suffer from deficiencies that are well known. Within a continuum theory, some of these flaws have already been remedied, see e.g., [11,12] and the literature cited therein. However, to develop models for the above mentioned applications, there is the more fundamental problem that additional effects like elastic deformation, stresses or interaction of charge transport with fluid flow have to be included into the models. However, extension of the standard model to cover these further effects and other material properties in a consistent way are far from being obvious.
The aim of this work is the formulation of a general modeling framework that allows the derivation of strictly thermodynamically consistent models for the above applications. With this article, we want to provide a revised and more concise version of [13] and want to make it available to a broader audience. Thereby, we close a gab that remained in our recent works on the improved continuum modeling of electrochemical interfaces [14][15][16][17][18][19]. We here resume the effort of non-equilibrium thermodynamics containing electromagnetic fields and its extensions to surfaces [20][21][22][23], which was combined to a single unified framework in [13] for the first time, together with the derivation of the according constitutive equations. Extending [24] to electromagnetic fields and to surfaces, we build our model on (i) universally valid balance equations in the bulk and on surfaces on the one hand and (ii) the formulation of an entropy principle and the according derivation of constitutive relations on the other hand. In particular, a chemically reacting mixture with neutral and electrically charged components will be described. Other phenomena such as diffusion, heat conduction, viscosity, elasticity, polarization and magnetization are taken into account.
When applying the resulting continuum models to electrochemical systems, typically very different scales in space and time arise and not all of them are necessary to capture the macroscopic relevant features of the considered system. Here, we show by means of dimensional analysis how considerable simplification of the model is possible and still-compared to the standard literature-a much more fundamental approach like [14] is possible to overcome deficiencies of the classical Nernst-Planck model. Thereby, we gain better insight into the internal double layer structure and the electrolyte transport by diffusive fluxes. Moreover, asymptotic analysis provides a mathematical tool to derive more simple reduced models. By this procedure, it is possible to recover on a macroscopic level the double layer capacity of electrode-aqueous electrolyte interfaces and some well established relations of electrochemistry, like Butler-Volmer equations and the Lippmann equation.

Outline
This paper is organized as follows: in the following section, we introduce notation and the geometrical setup. Then, in Section 3, we state the universal balance equations in the bulk and on the surface, containing the balances for the fields of matter and Maxwell's equations and their coupling. The transformation properties of the involved fields and the principle of Galilean symmetry are briefly discussed and stated in Section 4. In Section 5, we formulate the axioms of the entropy principle in bulk and surface and describe the general closing procedure based on the exploitation of the entropy principle. Subsequently, we derive in Section 6 constitutive relations in bulk and surfaces for magnetizable, polarizable, viscous and reactive mixtures and discuss in particular the role of polarization and magnetization. In Section 7, we apply the general model to liquid electrolytes, in particular aqueous electrolytes in contact with metal electodes and give a survey over recent key results. We close with some concluding remarks in Section 8.

Geometrical Setup
We consider a moving orientable surface S(t) dividing a possibly evolving domain Ω(t) ⊂ R 3 into two subdomains Ω ± (t) ⊆ R 3 with S := ∂Ω + ∩ ∂Ω − -see Figure 1. Let θ : [0, t end ) × ω → S be a smooth bijective parametrization of the surface S, where ω ⊂ R 2 is the open parameter domain. The partial derivatives of θ define tangential vectors, the surface normal and the metric As a convention, we chose the orientation of the surface normal such that ν is the inner normal of Ω + . Moreover, we indicate the components of vectors and tensors with respect to Cartesian coordinates by lowercase Latin indices, e.g., i, j, k, whereas we use uppercase Greek indices like e.g., Γ, ∆, Σ for the tangential components. For a vector V defined on the surface, we denote the normal component by V ν and write V ∆ τ , for ∆ = 1, 2, for the tangential components. For the matrix components of the metric tensor g, we use lower indices g ∆Γ and, for the components of the inverse matrix of the metric, we use upper indices g ∆Γ . We apply the convention of implicit summation over coordinate indices appearing twice. The curvature tensor b ∆Γ of the surface S and the Christoffel symbols Γ Σ ∆Γ are defined by a decomposition of the derivatives of the tangential vectors into their tangential and normal components, Then, the mean curvature of S is k M = 1 2 b Γ∆ g Γ∆ .

Covariant Derivatives
Let a : S → R be a scalar and V : S → R 3 a vector field. Then, the covariant derivatives of the tangential components are defined as Let w denote the velocity of the surface S. For a scalar a : [0, t end ) × S → R, we define the time derivative

Traces, Jumps and Mean Values on Surfaces
For a generic function u : [0, t end ) × (Ω + ∪ Ω − ) → R m , we denote for any time t and for each x s ∈ S the trace on either side of the surface by If the function u is not defined in Ω + or Ω − , we set the corresponding trace to zero, i.e., u + = 0 if u is not defined in Ω + and u − = 0 if u is not defined in Ω − . By we define the jump and the mean value of u at the surface S, respectively. The definition of the jump is related to the above convention on the orientation of the surface normal to be the inner normal of Ω + , cf. Figure 1.

Description of Reacting Mixtures
For quantities defined in the domains Ω + or Ω − , there will often be corresponding quantities on the surfaces S. As a convention, the same letters are used for these quantities, but the surface variables are indicated by an underset s.

Constituents and Chemical Reactions
In each of the two domains Ω + and Ω − and on the surface S, we consider a mixture of several constituents, usually referenced as A α . The index set of constituents in Ω ± is denoted by I ± . We assume that the sets I ± are disjoint, i.e., I + ∩ I − = ∅. This assumption takes into account that, even if a certain chemical substance is present in both domains, the functional dependence of the corresponding chemical potentials in general differs between the subdomains. All constituents of the subdomains Ω ± are assumed to be also constituents on the surface S, but, in addition, there may be some constituents that are exclusively present on S. Accordingly, the constituents on the surface S are represented by the set I S with I ± ⊆ I S . For the later construction of constitutive equations, we choose in each subdomain one designated constituent as a reference. We denote these constituents by A 0 − , A 0 + , A 0 in Ω − , Ω + and S, respectively.
For a description on a continuum level (macroscopic level), each of constituent A α is characterized by the (atomic) mass m α and its atomic net charge z α e 0 , where the positive constant e 0 is the elementary charge and z α is the charge number of the constituent. There may be M ≥ 0 chemical reactions among the bulk constituents and M S ≥ 0 chemical reactions on the surface. These reactions may be written in the general form The constants a k α , b k α are positive integers and γ k α := b k α − a k α denote the stoichiometric coefficients of the reactions. The reaction from left to right is called forward reaction with reaction rate R k f > 0. The reaction in the reverse direction with rate R k b > 0 is the backward reaction. The net reaction rate is defined as R k = R k f − R k b . Since charge and mass have to be conserved by each single reaction in the bulk and on the surface, we have Electro-and Thermodynamic State The thermodynamic state of Ω ± at any time t is described by the number densities n α and the velocities υ α for α ∈ I ± and the (specific) internal energy u. Analogously, the thermodynamic state of the surface S is characterized by the number densities of the surface constituents n s α and the surface velocities υ s α for α ∈ I S and the (specific) surface internal energy u s . The electrodynamic state of Ω ± and surface S at any time t is described by the electric field E and the magnetic field B. Alternatively, the electromagnetic behavior of polarizable and magnetizable matter could equally well be characterized by the vectors of polarization P and magnetization M, instead of E and B. Both sets of variables, (E, B) and (P, M), are coupled by constitutive equations and, for the derivation of these relations, it turns out to be favourable to use P and M, for reasons which are discussed in Section 6.4.
Multiplication of the number densities n α by m α gives the partial mass densities The mass density and the barycentric velocity of the mixture are defined by The free charge density is defined as The internal electronic structure of the constituents is not reflected by our macroscopic description, but it has relevance for the overall electromagnetic field. To represent these microscopic effects on the more macroscopic level, we introduce the polarization charge density n P , n s P in Ω ± and S, respectively, and then define the total charge density by Deformation Gradient The barycentric velocity causes deformation of bulk and surface. We formally introduce the deformation gradient F by means of a partial differential equation, viz.
where det(F) > 0 is assumed. For the relation of the so defined F to the deformation of a body under the velocity υ, we refer to the textbook literature, cf. e.g., ([25] p. 328, Equation (67.4)). From the deformation gradient, the unimodular, or volume preseving, deformation gradient F uni is derived as The temporal changes of the surface parametrization θ describes the deformation of the surface. Therefore, the tangential vectors τ 1/2 are the surface equivalents of the deformation gradients. Similar to the bulk, the unimodular tangential vectors are defined as

Universal Balance Equations of Electrothermodynamics
In this section, we introduce a set of balance equations in the bulk and on the surface that we postulate to hold universally, i.e., these balances are assumed to hold independent from the considered material. The surface balances account for the transport in tangential as well as in normal direction. Our approach is oriented on the classical work of Truesdell and Toupin [25] and we apply a notation similar to [21]. In this work, we consider only local balance equations that can be derived from their respective global counterparts, cf. [21].
The set of balance equations is subdivided into the classical balances of matter on the one hand and Maxwell's equations for the electromagnetic field on the other hand. Following [25], Maxwell's equations in the bulk are formulated with the postulation of universally valid Maxwell-Lorentz aether relations. For the surface, there is not such a standard formulation of the equations and, in particular, we here derive surface equations by analogy to the procedure in the bulk. With the later application to electrochemical systems in mind, we use the classical, i.e., non-relativistic form of the balances for the fields of matter.

Balance Equations of Matter
We consider partial mass balances for each of the constituents of the mixture, a single momentum balance of the mixture and an energy balance of matter. Modeling approaches based on different sets of balance equations also exists. For instance, each partial momentum ρ α υ α can be balanced. An introduction to these kinds of models for the bulk systems without surface balances and without electrodynamics can be found in [24].

Balance of Mass
In each of the subdomains Ω ± as well as on the surface S, the partial mass balances are in Ω ± for α ∈ I ± , (16a) Herein, r α and r s α denote the mass production rate of constituent A α . They are defined by the reaction rates, The bulk and surface diffusion flux J α and J s α with respect to the barycentric velocities are defined as The partial mass balances can be combined to derive conservation laws for the total mass density of the mixture in the bulk and on the surface. In later applications, often it turns out that it is advantageous to replace partial mass balances by these conservation laws. From balances (16) together with the constraints (18) and (8), we obtain the conservation of total mass Multiplication of the partial mass balances (16) by z α e 0 /m α and summation over all species together with conditions (8) implies conservation laws for free charge, with the (non-convective) free electric current densities

Balance of Momentum
The total momentum density ρυ of the mixture in bulk changes due to the momentum flux, which is given by the (Cauchy) stress tensor σ, as well as by the gravitational force density ρ f and the Lorentz force k. Accordingly, the surface momentum density ρ The momentum flux σ s is decomposed into its normal and tangential components The tensor with the components S Γ∆ is denoted as surface stress tensor and the vector with the components S ∆ is the normal stress vector. We neglect internal spin. This implies symmetry of the stress tensors and vanishing of the normal stress [21], i.e., for i, j = 1, 2, 3 and Γ, ∆ = 1, 2, we assume

Balance of Energy
The energy density of matter can be split into the internal energy density ρu and the kinetic energy density 1 2 ρ|υ| 2 . The splitting originates from the observation that the kinetic part of the energy can be eliminated by a suitable coordinate transformation whereas the internal energy cannot. This decomposition also implies a decomposition of the fluxes and the external sources into internal and kinetic contributions. The energy balance in the bulk and on the surface read Here, q, q s denote the heat fluxes and π, π s are the Joule heat for bulk and surface, respectively.

Electromagnetic Fields
Maxwell's equations for the electromagnetic field are based on two conservation laws: conservation of electric charge and conservation of magnetic flux, cf. [21,25].

Conservation of Electric Charge
The conservation equations of the total electric charge densities n e and n s e in the bulk and on the surface read ∂ t n e + div(n e υ + J e ) = 0 in Ω ± , (26a) Here, the total electric current flux densities are split into a convective and non-convective part, i.e., j e = n e υ + J e and j s e = n s e υ s + J s e . Then, we introduce the charge potential D and the current potential H as formal solutions of these charge balances, When considering the conservation of the free charge n F , the same argumentation leads to a formulation of Maxwell's equations in terms of the electric displacement field D = ε 0 E + P and the magnetic field density H = 1 µ 0 B − M, instead of D and H. Here, P and M are the polarization and magnetization, respectively, that we introduce later in Section 3.3. In contrast to, e.g., [26], our derivation of Equation (27b) in Appendix A is not based on some averaging technique for the bulk equations.

Conservation of Magnetic Flux
The volume-and surface equations for the electric field E and the magnetic flux density B are derived from the conservation law of magnetic flux, Maxwell-Lorentz Aether Relations The equation system (27) and (28) together constitute the system of Maxwell's equations and boundary conditions. The system is underdetermined such that additional relations between (B, E) and (D, H) are needed. Since all material dependence is incorporated into the total charge density n e , n s e , we can postulate universal valid Maxwell-Lorentz aether relations [21,25,27], Here, ε 0 is the dielectric constant and µ 0 is the magnetic constant. They are related to the speed of light by

Coupling of Equations for Matter and Electrodynamics
The coupling of the equations of matter in Section 3.1 and Maxwell's equations according to Section 3.2 is done in the following two steps. First, we identify the balances for polarization and magnetization. Second, by requiring conservation of total energy and total momentum, we identify the Lorentz force and Joule heat as functions of the electromagnetic fields.

Polarization and Magnetization
The conservation laws (26) for the total charge and (20) for the free charge imply also conservation of the polarization charge n P = n e − n F , n s P = n s e − n s F , i.e., ∂ t n P + div(n P υ + J P ) = 0 in Ω ± , (30a) with the polarization current densities according to Motivated by the introduction of D and H as formal solution of the balance equations (26), we use an analogous approach to introduce polarization P and Lorentz magnetization M by formal solution of the balance equations (30) according to Appendix A. We get This introduction of the fields P and M differs from the classical textbook literature, where often microscopic models are used, e.g., P is introduced by considering microscopic electric dipoles and M is derived from microscopic circular currents. Then, the polarization charge n P and polarization currents J P are formally introduced to couple P and M to the Maxwell's equations, cf. [21,[27][28][29][30][31]. Our approach leads to the same relations between (n P , J P ) and (P, M), but it has the advantages that it is independent from any microscopic model and, most notably, its ability to transfer the concept of polarization and magnetization to the surface. The relation of the fields (P,M) to (E,B) is left to a subsequent material modeling. For the fields P and M in the bulk, there are no according counterparts on the surface in our approach. An extension of the theory seems to be possible. For instance in [32], surface polarization and magnetization are derived from the bulk Maxwell's equation by integration over a singularity and thus for each bulk quantity there exists a corresponding surface quantity. While in our approach P and M are introduced as a particular formal solution of the balance equation for polarization charge, there may be a further formal solution of the surface balance equation that allows the introduction of surface polarization vector and surface magnetization. Although we did not introduce a surface polarization vector or surface magnetization, there exists surface polarization charge n s P and surface flux densities J s P due to polarization and magnetization, cf. Equation (32b).
Due to its importance in electrodynamics and for the upcoming constitutive modelling, we introduce the electromotive intensity and the magnetization

Lorentz Force and Joule Heat
From Maxwell's equations (27) and (28), two more balance equations can be derived: the balance of the electromagnetic momentum and the balance of the energy of the electromagnetic field.

Balance of Electromagnetic Energy
The electromagnetic energy density in the bulk is defined as 1 2 (E · D + B · H). From Maxwell's equations, we get the balance equations, [21,26]: Here, E × H is the electromagnetic energy flux density (Poynting vector) and (n e υ + J e ) · E is the energy gained by the matter due to the electromagnetic field. On the surface, there is no additional electromagnetic surface energy density/flux implied by Maxwell's equations. Thus, the influxes from the bulk are solely balanced by the surface production n s e υ s + J s e ·Ē.

Balance of Electromagnetic Momentum
The electromagnetic momentum density is defined as D × B. Maxwell's equations imply in bulk and surface the balances for electromagnetic momentum [21,26]: The electromagnetic momentum flux density (H ⊗ B + E ⊗ D) − 1 2 (H · B + E · D)I is called the Maxwell stress tensor. On the right-hand sides of (35), the productions n e E + (n e υ + J e ) × B and Maxwell's equations (27) and (28) do not imply surface momentum density/flux, such that the bulk influxes are balanced by the surface production.

Identification of Lorentz Force and Joule Heat
We postulate that, in the absence of gravitation, i.e., for f = f s = 0, the total momentum ρυ + (D × B) and the total energy ρ(u + 1 2 |υ| 2 ) + 1 2 (E · D + B · H) of matter and electromagnetic field are conserved quantities, cf. [21,25,26]. Then, the production in the balance equations of matter have to be canceled by the electromagnetic productions. This implies for the Lorentz force k, k s and the Joule heat While the definition of mixture quantities like the total energy and total momentum in the authors point of view should not be modified, terms in the balance equations can be moved from the fluxes to the productions and vice versa. Discussions in literature on the Lorentz force and Joule heat, cf. e.g., [31], mostly originate from this ambiguity.

Balances of Total Momentum and Internal Energy
As a consequence of (36), we get the total momentum balances where the total stress tensor consisting of Cauchy and Maxwell stress is given by where 1 denotes the identity matrix. In the absence of the gravitational force and assuming stationarity and vanishing velocity, the total momentum balances reduces to In particular, if the surface stress vanishes, then the total stress is continuous at the surface. Therefore, one can only measure the total stress, but not separately the Cauchy stress or the Maxwell stress. This fact can be nicely studied in an experiment with two capacitor plates dipped into a liquid, cf. ( [29], Section 3.6).
For the later exploitation of the entropy principle, we derive an appropriate form of the balance of internal energy. The internal energy balance is derived in three steps: first, the kinetic energy is eliminated in the energy balance (25) by means of the momentum balance (22). Second, the electromagnetic contributions are split into a free and polarization part. Third, the polarization current is replaced by using the identities for the polarization vector and Lorentz magnetization (32b):

Symmetry Principles for Observer Transformations
In this section, we study changes of space and time coordinates of the observer frame of reference, shortly called observer transformations. Their application to the fields of thermodynamics and electrodynamics show diverse transformation properties leading to symmetry principles that particularly restrict the generality of the admissible constitutive functions.
In the previous section, the equations of balance for matter, the conservation law of charge, the magnetic flux conservation and the Maxwell-Lorentz aether relations were formulated with respect to an inertial frame of reference. Time and spatial coordinates refer to an inertial frame if two conditions are met: (1) the mass center of a material body that is not subjected to external forces moves with constant (barycentric) velocity along a straight line, and (2) the Maxwell-Lorentz aether relations hold.

Symmetry Principles
The most fundamental symmetry principle is the principle of relativity. It restricts balance equations as well as constitutive equations by stating invariance with respect to arbitrary observer transformations. However, this principle can only be maintained if time and space and all involved equations are properly combined to four-dimensional objects.
We will not touch the four-dimensional case further and restrict ourselves to a less general symmetry principle suitable for the 1+3-dimensional case. In the following, we introduce and exploit the Galilean symmetry principle stating invariance of the involved equations with respect to Galilean transformations defined below in Equation (41). It is well known that the balance equations of matter (16)- (25) are invariant with respect to Galilean transformations. On the other hand, the Maxwell equations are invariant with respect to arbitrary transformations, including Galilean transformations as well. However, the 1+3-dimensional Maxwell-Lorentz aether relations are only invariant with respect to Lorentz transformations. Nevertheless, in the limit of vanishing barycentric velocity, i.e., v/c → 0, the Galilean transformation is a good approximation of the Lorentz transformation. However, even in this case, one should be aware that by applying the Galilean symmetry principle, important classical effects are ignored. Two examples are: (1) Galilean invariance is intimately linked to the conservation of mass in a chemical reaction, thus we cannot describe the case of atomic fusion in this setting. (2) The unipolar generator cannot be explained with models restricted by the Galilean symmetry principle. A relativistic explanation is needed. For the electrochemical applications that we have in mind, these inconsistencies between the Maxwell-Lorentz aether relations and the Galilean symmetry principle are not relevant and are not further discussed here.
Finally, we mention the Euclidean transformation, which is another classical transformation that generalizes the Galilean transformation. For these transformations, often invariance is proposed with respect to constitutive functions only, and the corresponding symmetry principle is called material frame indifference. In this paper, we do not consider that concept any further.
In summary, we may say that the general essentials of symmetry principles for observer transformation consist of three steps: (i) We choose a class of space-time transformation. (ii) We investigate the transformation properties of the involved quantities under that transformation. (iii) We propose a corresponding symmetry principle stating invariance of certain equations, for example invariance of the constitutive equations.

Galilean Symmetry Principle
For a Galilean transformation, the transformed coordinates are given bȳ where the orthogonal matrix O and the velocity W are time independent. We now introduce the principle of Galilean symmetry that states: The equations of balance for matter, Maxwell's equations and the constitutive functions must be invariant with respect to the Galilean transformation (41).
The fields T i 1 ,i 2 ,...,i N are identified as the components of a Galilean tenor (of rank N) if they transform under Galilean transformations according tō For p = 1, T is called an axial Galilean tensor while p = 0 indicates an absolute Galilean tensor. The special cases N = 0, 1 refer to Galilean scalars and Galilean vectors, respectively.

Classification of the Involved Fields
For the velocity, the deformation and the tangential vectors of the surface, it is possible to derive the transformation properties with respect to Galilean transformations as Next, we classify the involved fields in the balance equations for matter, such that these remain invariant under Galilean transformations, i.e.,

The Entropy Principle
The system of balance equations from Section 3 has to be closed by constitutive equations in the bulk and on the surface. The-thus far undetermined-constitutive quantities are the partial mass fluxes J α , J α s , reaction rates R k , R s k , heat fluxes q, q s stress tensors σ, S, polarization P and magnetization M. The constitutive equations are not uniquely determined, but they are restricted by the second law of thermodynamics, which is here expressed in terms of the entropy principle, and the principle of Galilean symmetry. We here follow the strictly axiomatic approach of [24] where also slightly different forms in the literature are reviewed, cf. e.g., [20,21,33,34]. The formulation of the entropy principle is largely analogous in the bulk and on the surface, and therefore done simultaneously.

Formulation of the Entropy Principle
The entropy principle consists of five axioms, of which the first part consisting of axioms I-III is universal, i.e., independent from the considered material. The second part with the axioms IV-V is material dependent, although in a very general and abstract way, because an appropriate set of independent variables has to be chosen for the considered material. Here, we consider elastic, viscous, magnetizable, polarizable, heat conducting and reactive mixtures.
The factors of each product are either scalars, vectors or tensors with respect to Galilean transformations. Each product N A P A resp. N s A P s A can be associated with exactly one dissipation mechanism. (iii) The entropy productions vanish in equilibrium.
The (absolute) temperature T, T s and the chemical potentials µ α , µ s α are defined as ,

Remarks
Axiom IV contains a specific set of variables. The entropy principle can be formulated with different sets of variables as well, and these different choices will in general lead to different constitutive equations. In particular, we assume that the bulk entropy density depends on the internal energy ρu + M · B instead of ρu. This internal energy allows us to derive an entropy production, which satisfies Axiom III in particular the representation (47). Different other choices for the internal energy that allow a respective formulation of the entropy principle can be found in [31]. In Section 6.4, we show that not every choice of constitutive functions that satisfies Axiom III automatically admits relaxation to the equilibrium.
Axiom IV requires that the set variables for the entropy densities are independent from each other. This is not the case for the deformation gradient F and the tangential vectors τ ∆ , for ∆ = 1, 2, since the mass balances (19) are solved by )]] = 0. Therefore, the deformation gradients and tangential vectors cannot be used in Axiom IV and they are replaced by their unimodular counterparts. There is a difference between the constitutive equations of the entropy density for bulk and surface because in our approach there are no counterparts of the fields P and M on the surface.

Exploitation of the Entropy Principle-General Approach
The general strategy for the design of constitutive relations can be outlined in a material independent way. Using only the structure (47) of the entropy production, we derive thermodynamic consistent constitutive equations from either linear or nonlinear relations between the factors N A and P A , such that the entropy production is non-negative.
For the ease of presentation, we restrict ourselves to the case that the factors of each binary product are scalars. A generalization to vectors and tensors can easily be done.

Linear Relations
To account for cross effects, we consider a regular matrix M and set Q = M −1 . Then, the entropy production (47) is rewritten as To have ∑ A N A Q AB proportional to ∑ C M BC P C , we choose positive phenomenological coefficients L B > 0 and set where the coefficient matrix K is symmetric and positive definite by construction. The coefficients of the matrix M and L B in general may be functions of the independent variables as long as they are absolute scalars with respect to Galilean transformations. We want to stress here that, in the generality of the above formulation, the entropy principle provides no information about this dependencies on the independent variables apart from the requirements that L B ≥ 0 and M regular. Instead, theories of different nature or, like e.g., in [35,36], specific assumptions on the considered system are needed to derive such functional dependencies.

Nonlinear Relations
Thermodynamically consistent constitutive equations based on linear relation (51) can also be derived by nonlinear closure relations. Since chemical reactions are often modeled by an Arrhenius type exponential dependency of the reaction rate on the driving force, we thus propose the nonlinear relation where the phenomenological coefficients L B and K B are positive. For a system close to equilibrium, such that 0 ≤ |K B ∑ C M BC P C | 1, the exponential function in (53) may be linearized and the nonlinear closure relations (53) then approach the linear ones in (52).

Galilean Symmetry Principle
The constitutive relations that result from the above closure relations automatically satisfy the Galilean symmetry principle, since the factors of the binary products are scalars, vectors or tensors with respected to Galilean transformations according Axiom III (ii) and for the phenomenological coefficients only absolute scalars according to Section 4 are chosen.

Remark on Cross Effects and Symmetry of the Phenomenological Coefficient Matrix
The well known Onsager-Casimir reciprocal relations postulate the symmetry of the matrix K that appears in closure relation (52). They are motivated by experimental observations of thermo-electrical effects by Thomson ([37], pp. 237-241) and have been proven for systems of ordinary differential equations describing homogeneous processes in a statistical mechanics by Onsager and Casimir [33,[38][39][40].
By the construction based on the matrix M, the symmetry of K is automatically guaranteed. Moreover, the introduction of M in the entropy production (51) allows in a very simple way the realization of cross effects in the nonlinear case. However, one should note that, as long as the factors N A and P A have not been specified, it is also possible to construct different linear closure relations such that the coefficient matrix K is positive definite, but anti-symmetric(!), cf. [24]. There, the concept of parity is introduced to solve these problems by further characterization of the factors N A and P A . If the quantities N A exclusively consist of objects with physical units containing odd powers of the time unit and all objects of P A have units with even powers of the second, then we have symmetry of K by our construction. Now, it is easy to observe that an exchange of some N A objects with P A objects then yields antisymmetry of K.

Constitutive Relations for the Bulk
To derive the bulk constitutive equations for the mass fluxes J α , reaction rates R k , heat flux q, stress tensor σ, polarization P and magnetization M, we first identify the entropy production. To do so, we insert the entropy density (48a) into the entropy balance (46a) and apply the chain rule of differentiation. Then, the occurring time derivatives are substituted by means of the balance equations, where in particular we apply the balance of the internal energy (40a). At this stage, it is still possible to rearrange the terms in several ways. Once the entropy flux is fixed, the entropy production is uniquely determined. In order to get the structure (47), we choose Other constitutive equations for the entropy flux different from (54) can be admissible as well. In general, such alternative choices will lead to different constitutive equations for the entropy production. For a discussion of alternative choices of the electromagnetic contributions to the entropy density, we refer to Section 6.4 and in [41] a non-trivial example of generalized entropy fluxes in the context of viscous Cahn-Hilliard equations is discussed.
With the choice of the entropy flux (54), the entropy production is given by a sum of six binary products. We can identify six dissipation mechanisms related to their specific entropy production: shear viscosity ξ SV , volume viscosity ξ VV , (bulk-)reactions ξ R , thermodiffusion ξ TD , polarization ξ P and magnetization ξ M . The entropy production is then given by the constitutive equation Here, we defined In the derivation of the entropy production, we used the constraint (18a) to eliminate the flux of the species A 0 , such that only linearly independent mass fluxes appear in the entropy production. Moreover, we used a symmetry condition which originates from the transformation properties of the thermodynamic fields, cf. Appendix B, i.e., we have Furthermore, in (56), the binary products of the entropy production are arranged in a way such that the transformation properties required by Axiom III (ii) are guaranteed. The binary product due to shear viscosity is formulated in such a way that the matrices are symmetric. Thus, the anti-symmetric part of (E ⊗ P)∇υ and (M ⊗ B)∇υ is shifted to the entropy production of polarization and magnetization, respectively.
To derive constitutive equations from (55), we apply both approaches of Section 5.2. In general, different dissipation mechanism can be coupled, e.g., volume viscosity can be coupled with bulk reactions. To reduce the complexity of the constitutive equations, we only consider cross effects related to thermodiffusion. Since the index sets I ± are assumed to be disjoint, the closure procedure can be applied in both subdomains independently.

Thermodiffusion
For the mass fluxes, we choose a linear relation with cross effects, viz.
The coefficient matrix κ L L T M is symmetric and positive definite. In particular, the heat conductivity κ and the mobility matrix M are symmetric and positive definite.

Reactions
For chemical reactions, often an exponential Arrhenius-type dependence on the temperature and some activation energy is expected. Therefore, we choose the nonlinear relation (53) for the chemical reactions, with positive coefficients A k , R k 0 . For simplicity, we neglected cross effects between the different reactions.

Viscosity
We choose linear relations for the volume viscosity and for the shear viscosity, viz.

Polarization and Magnetization
For the vector of polarization P and the magnetization M, we choose linear relations without cross effects, Here, the phenomenological coefficients τ P > 0 and τ M > 0 are the relaxation times of polarization and magnetization.

Constitutive Relations for the Surface
The exploitation of the entropy principle for the surface is analogous to the bulk. On the surface, we have to determine the heat flux q s , the mass fluxes J s α , the stress tensor S and the reaction rates R s k .
Moreover, we have to determine the normal components of the heat flux q ± ν , of the mass fluxes J ± ν,α and of the stress tensor σ ± · ν.
First, we substitute the entropy density (48b) into the entropy balance (46a) and choose the entropy flux in a way that allows for the structure (47) of the entropy production, viz.
Then, we identify six dissipation mechanisms: tangential surface viscosity ξ As in the bulk, the entropy production is formulated such that the factors of the binary products are Galilean scalars, vectors or tensors, respectively. Furthermore, the transformation properties of the thermodynamic fields restrict the entropy density (see Appendix B), and yield the symmetry conditions The surface entropy production (65) shows some similarities to the bulk entropy production (55), but also some differences. While the bulk entropy production contains two dissipation mechanisms due to polarization and magnetization, there are no similar contributions on the surface because we restricted the constitutive functions ρ s η s on a singular surface to be independent from the electromagnetic fields. Moreover, there are on the surface dissipation mechanisms of different types: while the one class contains only fluxes that are tangential to the surface, the second class consists of those related to the fluxes normal to the surface. In general, it is possible to introduce cross effects between these two different classes, cf. [42]. For simplicity, we do not discuss these kinds of cross effects here.

Surface Reactions
The exponential form of the Butler-Volmer equations for the surface reaction rates, cf. e.g., [10], suggests to apply the nonlinear closure (53). Similar to the bulk and neglecting cross effects between the different reactions, we choose the exponential relations Surface Viscosity Analogously to the bulk, we define the shorthand notation The linear closure yields for the trace and for the deviatoric part of the surface stress tensor T s the constitutive equations: The phenomenological coefficients satisfy λ Mass Flux and Stress Normal to the Surface For the entropy contributions coming from the bulk, we here for simplicity choose the linear relations on S. However, we want to mention that the application of the nonlinear closure (53) sometimes can be more appropriate, e.g., in cases where adsorption to the surface becomes the rate limiting step for electron transfer reactions, cf. [17]. We get The coefficients κ

Comparison of Bulk and Surface Equations
Comparison of the constitutive equations for bulk and surface reveals that, in the case of vanishing polarization and magnetization, the bulk and surface equations have an analogous mathematical structure. The differences are due to the fact that we do not consider surface polarization and surface magnetization in constitutive equation (48b) for the surface entropy density. For vanishing polarization and magnetization, the electromagnetic field solely contributes to the mass flux and heat flux in the form of the electromotive intensity E. Although a contribution of the electromagnetic field to the surface mass fluxes (67b), which results from the linear closure we applied here, has to be expected due to the symmetry properties between bulk and surface, it is not present in literature [21,22,43].

Free Energy Density
The temperature T is a quantity of central interest in thermodynamics. It is here defined as the derivative of the entropy density with respect to the internal energy density. For the construction of constitutive equations, it is often beneficial to use the temperature T resp. T s as an independent variable, instead of the internal energy density. To replace the internal energy density as an independent variable by the temperature, the free energy density ρψ is introduced by means of Legendre transformation of the entropy density ρη, viz.
For the free energy density, we thus have representations as functions of the temperature From the definition (73) and using the transformation properties (75), we get relations between the internal energy density and the free energy density: The constitutive equations (62) and (71) then imply that pressure and surface tension depend on deformation, viscosity, polarization and magnetization.
In the case of vanishing electromagnetic fields and vanishing viscosity, the constitutive equations for pressure and surface tension simplify to These equations are the well known Gibbs-Duhem relation and its counterpart on the surface, cf. [20,21].

Adsorption
In chemistry, an adsorption process is often described as a reaction that contributes to the mass production r s α for α ∈ I S . A similar approach is not possible in the context of electrothermodynamics because the constitutive equations for the surface reactions cannot directly couple bulk and surface species. Thus, adsorption of species is here determined by the constitutive equations (72d) and (72c) for the mass fluxes in normal direction onto the surface. The driving forces for the mass fluxes are the differences of chemical potentials between bulk and surface. The adsorption rate of a species A α is thus determined by a kinetic coefficient M s ± α in Equation (72c). If species A α is non-adsorbing, then M s ± α = 0.
In addition, the total mass flux ρ(υ ν − υ s ν ) ± in normal direction onto the surface is determined by the constitutive equation (72d). In particular, it depends on the Cauchy stress tensor. Under the assumption that the viscosity in the bulk is small and the free energy density in the bulk is independent from the deformation gradient, i.e., ρψ = ρψ(T, (ρ α ) α∈I ± , P, M), then relation (72d) simplifies to (due to constitutive equation (62) for σ) If also the electromagnetic contribution (P × B) ν (υ ν − υ s ν ) as well as kinetic contribution 1 2 ρ|υ s − υ| 2 vanishes, then the total mass flux is determined by the adsorption of species A 0 .

Discussion on Polarization and Debye Equation for Dielectric Relaxation
The constitutive equations of this section are consequences of the choice of the variables for the entropy density (48a). In addition, various different sets of variables can be compatible with the axioms I-III of the entropy principle as well. In particular, the choice of P and M is not mandatory, and it might seem more natural to use E and B instead.
In local equilibrium, i.e., ξ P = 0 and ξ M = 0 in the entropy production (55), there is no preference for either choice. The constitutive equations (63) then imply By the Legendre transformation, the variables (P, M) can be replaced by (E, B), cf. e.g., the discussion of the different sets of variables in [31].

Relaxation
The situation changes decisively in non-equilibrium. For simplicity, let us consider a system with only one species and two different constitutive functions for the entropy density In each setting, an according temperature is then defined as The corresponding free energy densities that follow from (82) and (83) are This coincides with the cases (a) and case (c) in ( [31], Section 3.3) if the internal energy density is defined as ρu + M · B, as we did in axiom IV of the entropy principle. The corresponding constitutive functions of the free energy densities are Then, exploitation of the entropy principle yields for the two cases the respective entropy productions due to polarization Both ξ PM P and ξ E B P , are compatible with the structure of the entropy production (47) in the axiom III(ii). Comparing Equations (86a) and (86b), we observe that upon an interchange of the variables P and E the two relations are almost identical. However, while the entropy production ξ PM P begins with a positive sign, there is a minus in ξ E B P . For further analysis, it is necessary to specify the dependence of free energy density on (P, M) respective (E, B). We consider the most simple choice with a constant electric susceptibility χ. When now applying the linear closure to derive the constitutive equations, we get in the case of vanishing velocity, i.e., υ = 0 the constitutive equations where the relaxation constants τ PM , τ E B > 0 are positive. We conclude from the constitutive equation such that in both cases the polarization is proportional to the electromotive intensity, i.e., In non-equilibrium, the situation is different. For simplicity, let us assume that the magnetization is in local equilibrium, i.e., ξ M = 0, and the free energy densities (85) are independent from M and B respectively. Under this simplification M = 0, B = 0 and E = E, and Maxwell's equations simplify to div(ε 0 E + P) = 0 , curl E = 0 .
Then, the divergence of the constitutive equations (88) simplify to Now we see that, in Equation (92a), div(P) vanishes for t → ∞. In contrast, Equation (92b) implies that div(E ) blows up in time with an exponential growth. We conclude that, although both choices of independent variables, (P, M) and (E, B), yield the same equilibrium relation (90), and in both cases the entropy production is non-negative due to the linear closure (88), the latter choice of the relaxation equation (92b) is incompatible with the equilibrium relations. Therefore, the approach based on the independent variables P, M should be strongly preferred. An analogous argumentation also holds for the magnetization. We may recapitulate the described problem by stating: the entropy principle does not necessarily prevent instability of the resulting system of field equations. This phenomenon is already known in a different area. Non-Newtonian fluids can be modeled with differential type constitutive functions for the stress, such that the symmetric part of the velocity gradient and its time derivatives are among the variables. This approach leads likewise to exponential growth in a situation where the fluid is expected to approach equilibrium. The problem is removed by interchanging the role of the velocity gradient and the stress as variable and constitutive quantity. The details are carefully described in [44,45].

Debye Equation
The Debye equation is used to describe for ideal systems the dielectric relaxation, i.e., the response of the electric field inside a dielectric material to the excitation by an oscillatory outer electric field, cf. e.g., [20,46]. In our setting, the Debye equation is given by Equation (88a) and its derivation is an immediate consequence of the entropy principle of Section 5.2 and it is identical to the Debye equation found by deGroot and Mazur in ( [20], p. 400) for matter at rest. In the case of non-vanishing velocities, i.e., υ = 0, the Debye equation found by deGroot and Mazur is not identical to Equation (63a) in this paper.

Application to Electrochemical Systems
The constitutive equations of Section 6 have been derived without making use of any particular specific material properties. They only rely on the universal balance equations and the entropy principle. All material properties of a specific electrochemical system thus have to be incorporated into the constitutive functions of the entropy and the phenomenological coefficients. In this section, we illustrate the application of the developed theoretical framework for the example of liquid electrolytes in contact with metal electrodes.
Our main interest is the accurate description of the charge transport in electrolytes and electrochemical processes at the electrode-electrolyte interface. By interface, we always mean the compound consisting of surface and the adjacent boundary layers from both sides. We summarize recent results [14][15][16][17][18][19], which are related to an improved understanding of the double layer structure that in particular allows quantitative and qualitative prediction of the differential double layer capacity, and enables a better understanding of electrocapillarity effects. Moreover, the above developed framework allows the formulation of extended Nernst-Planck fluxes and constitutive equations for surface electron transfer reactions that allow to recover generalized Butler-Volmer equations.
A particular difficulty for the development of mathematical models for batteries, fuel cells, electrolyzers or other electrochemical systems is due to the complexity that results from very different scales in space and time, the coupling of bulk and surface processes and the interplay between different physical phenomena like electrical and mechanical phenomena. Therefore, we first identify simplifying assumptions that are appropriate to adapt the most general theory developed above to the description of liquid electrolytes and metal-electrolyte interfaces.
we get the following dimension less form of Maxwell's equations in the bulk The non-dimensional bulk definitions of electromotive intensity and magnetization read Depending on the chosen characteristic reference values, the size of the dimensionless quantities may differ by several orders of magnitude, allowing considerable simplifications of the system.
Under the assumption that the derivatives ofȆ +P remain bounded, Maxwell's equations simplify in the asymptotic limit λ → 0. Rescaled to dimensional quantities, we get the equations of magnetostatics for the magnetic flux density B, which, in conjunction with the equation of matter, determine the electric field E.
Electrostatics: β, 1/δ → 0 In this limit, the equations of electrostatics for the electric field are in dimensional form and the electromotive intensity is identical to the electric field, i.e., E = E. The remaining Maxwell's equations determine the magnetic flux density B, i.e., div Equation (98a) implies the existence of an electrostatic potential ϕ, such that Assuming a simple constitutive equation for the polarization, i.e., P = χε 0 E, where χ is the electric susceptibility, we get from Equation (98b) Poisson's equation for the electrostatic potential, viz.
The dimensional analysis can be performed in analogue manner for the surface equations, cf. [13].
In the limit β → 0 and 1/δ → 0, Maxwell's surface equations yield two boundary conditions for the electrostatic potential, where the susceptibility χ can have different values on both sides of the surface. Due to the second condition in Equation (102), the jump of the potential [[ϕ]] is determined up to a constant and, in particular, this constant is independent from the material. Since the equations for the electrostatic potential allow a normalization, we can set this constant to zero. Thus, the electrostatic potential is continuous at the surface, which allows for defining the electrostatic potential of the surface as Application to Electrochemical Systems We consider a system that is varying moderately in time, has a number density in a typical range for solids and liquids, and a magnetic field strength that does not significantly exceed the field of common permanent magnets, viz.
An important feature of electrochemical systems is the formation of narrow double layers at the contact of different materials. The double layer is characterized by a typical width in the range of nanometers and a strong electric field in the range of 1 V nm . We thus choose the reference values implying for the dimensionless quantities The smallness of β 2 and 1/δ 2 compared to λ 2 suggests for the determination of E the use of the electrostatic limit equations (98), or (101) and (102), respectively. In conclusion, for electrochemical systems, the electrostatic approximation is valid, as long as the characteristic values of the system at hand are in the order of magnitude defined in Equations (104) and (105). In particular, in systems where the characteristic time scale is smaller, for instance for impedance measurements, the application of the electrostatic approximation is not appropriate.

Free Energy Models for Liquid Electrolytes and Metal-Electrolyte Interfaces
To develop a complete mathematical model for liquid electrolytes and metal-electrolyte interfaces, we start with some simplifying assumptions. The metal can be described by applying the Sommerfeld free electron model [47]. Then, it turns out that most of the relevant information of the metal-electrolyte interface can already be obtained without solving the model equations on the metal side, cf. [18]. Thus, we here ignore the metal electrode and only consider the electrolyte subdomain that further on is denoted by Ω with the corresponding index set I. For the description of polarizable liquid electrolytes, it is reasonable to neglect the bulk and surface deformation gradients F uni and τ uni 1/2 and, in particular based on the previous dimensional analysis, the dependency of the free energy function on the magnetization M. The constitutive function for the free energy densities simplify to We consider isothermal processes only. Thus, the temperature in the bulk and on surface is given by an appropriate reference temperature T ref that usually be the room temperature: The energy balances (25) then only serve to determine the heat fluxes which are necessary to enable an isothermal process. In the following, free energy functions for bulk and surface are developed. We only sketch the main ideas of the derivation and refer for more details to [13][14][15]18]. In particular, we illustrate here how elasticity contributes to the free energy function and how finite ion sizes can be incorporated into the electrolyte model.

Bulk Free Energy Density
The free energy density is not directly measurable. Therefore, a derivation of a free energy density must be based on constitutive relations that are in simple equilibrium situations backed by sufficient empirical evidence or that are derived from a microscopic theory. For liquid electrolytes, we consider free energy functions which consist of four contributions that are related to polarization, mechanical stresses, mixing entropy and a temperature dependent reference contribution, i.e., For the derivation of the individual contributions of the free energy, it is convenient to introduce the total mole density of particles n and the mole fractions y α by n = ∑ α∈I n α and y α = n α n with ∑ α∈I y α = 1 .
The entropy of mixing accounts for the number of possible arrangements of ions and solvent molecules for a give macroscopic state. It is determined by statistical thermodynamics by means of the Boltzmann formula. Thus, the corresponding free energy contribution reads The mechanical part of the free energy density function is derived by integration of the identity , which in turn is derived from a generalization of the Gibbs-Duhem equation (78), cf. [13,48]. Equation (78) is well known for non-polarizable materials. Since the polarization (116) is independent from the number densities, the definition of the Gibbs-Duhem equation (78) can be easily extended to the case handled here (see [14,15]). We assume a simple linear constitutive relation for the pressure p, Herein, p ref is the pressure of the reference state and K denotes the bulk modulus of the electrolyte.
By v ref α , we denote the partial specific volume of the constituent A α under the pressure p ref and temperature T ref . The function H is the mean specific volume of the mixture and accounts for volume changes due to a local variation of the mixtures composition. The mechanical contribution to the free energy density is derived from the constitutive equation (112) as The reference contribution to the free energy is assumed to be where ψ ref α denotes the reference free energy of each individual constituent. In the reference values ψ ref α , the specific heat is encoded, but this is not further outlined here because of the assumption of isothermal processes.
In equilibrium, we want to maintain the simple constitutive relation between electric field and polariziation An integration of relation (63a) in equilibrium and in the electrostatic approximation, i.e., 0 = − ∂ρψ ∂P + E then yields the corresponding contribution to the free energy density due to polarization Here, the susceptibility χ is assumed to be independent from the number densities, but may be dependent on the temperature. The assumption that the susceptibility is independent from number densities is not necessary, but this leads to simplified constitutive relations, in particular the chemical potentials are independent from polarization. From the free energy contributions (109)-(114) and the definition (49a), we get the chemical potentials Here, we used the constitutive relation (112) to express the pressure dependency of the chemical potentials.

Surface Free Energy Density
A simple surface free energy model of a metal-electrolyte interface was proposed in [18] that we will summarize here. On the surface, there is no contribution from the electric field, but, as in the volume, we assume additive mechanical, entropic and reference contributions, i.e., where the reference contributions are of the form The reference values may depend on the temperature and in general also on the crystallographic orientation of the surface. The set of constituents I S on the surfaces contains all electrolytic bulk species, the metal bulk constituents, which are the electrons e and metal ions M, and in addition a certain number of reaction products of the bulk species. We assume that electrons interact neither in entropic nor in elastic manner with the other surface species and thus only contribute to the reference energy. This assumption allows to consider the specific volume of the electrons to be vanishing. Let a ref M and a ref α denote the specific area of a metal ion and adsorbates of species α ∈ I S , respectively. We assume that the surface is built from an one atomic layer of metal ions that offers adsorption sites to the electrolyte species and the reaction products. Since some of the sites may be empty, we introduce the surface number density of vacancies by The entropic contribution to the surface free energy is determined by statistical thermodynamics as in the bulk and takes into account the mixing entropy of the adsorbates and reaction products, The mechanical contribution of the free surface energy is based on a simple constitutive model for the surface tension. Let γ Here, the surface tension is assumed to be a function of the metal ion density n s M only. However, in general, adsorption and surface reaction may change the metal density and therefore the surface tension. Since in general K s is large for metals, small changes in n s M lead to significant changes in the surface tension. Thus, the surface tension is also indirectly related to the surface coverage with adsorbates. Insertion of the definition (123) into constitutive equation (78) right and integration after an appropriate variable transformation, see [18], Appendix A, yields for the mechanical contribution to the free energy, The resulting chemical potentials are In contrast to the bulk, the surface chemical potentials of the adsorbates are independent from the stress tensor, i.e., surface tension, but we have contributions from the vacancies on the underlying lattice.

The Incompressible Limit
For liquid electrolytes, in particular for aqueous electrolytes, we expect that the volume does not change significantly if the pressure is varied. We can incorporate this behavior by means of the asymptotic limit K/p ref → ∞. In an analogous way, we consider on the surface the limit K s /γ s ref → ∞.
In the incompressible limit, the relations (112) and (123) cannot be used anymore to determine the pressure p and the surface tension γ s , respectively. Therefore, p and γ s become new independent unknowns of the system and, instead of the constitutive equations (112) and (123), there are the so-called incompressibility constraints The chemical potentials become linear in the pressure respective surface tension, i.e.,

Solvation
In polar solvents, in particular in water, the microscopic dipoles of the solvent molecules give rise to a microscopic electrostatic interaction between the solvent and the ionic species. This interaction leads to clustering of a finite number of solvent molecules around an ion, which is known as solvation. In Figure 2, the solvation of ions is illustrated. The solvation concept can be easily transferred from the volume to the surface as also illustrated in Figure 2. Solvation has a profound impact on the mixing entropy and the specific volume of the ions within the electrolyte model: solvent molecules that are bounded by the ions do not participate in the entropic interaction with the other constituents of the electrolytic mixture. Therefore, we handle an ion and its solvation shell as a single entity called solvated ion. Binding of solvent molecules into solvated ions decreases the amount of remaining free solvent molecules in the electrolyte. Let

The Electrochemical Double Layer in Equilibrium
At the electrode-electrolyte interface, ionic as well as electronic species can accumulate, forming boundary layers of only few nanometer thickness on both sides of the surface. This structure is commonly known as electric double layer [10,49]. It is a key feature of many electrochemical applications. For the theory presented here, the double layer is of particular interest because it allows valuable insight into the free energy which is not directly accessible to measurements. Since the width of the double layer is often much smaller than the macroscopic length scales of electrochemical components or cells, it can be asymptotically approximated by planar or radially symmetric (semi-)infinite systems.

Boundary Layer Structure
We consider a half space problem, where the domain x > 0 is occupied by an incompressible liquid electrolyte. The region close to the boundary surface at x = 0 represents the boundary layer. A potential difference is prescribed between x = 0 and x → ∞. From the constitutive relations derived in Section 6.1, we deduce that the equilibrium equations in the bulk resulting from the mass and momentum balance are, cf. [14], ∇(m α µ α + z α e 0 ϕ) = 0 for α ∈ I , (128a) Integration of Equation (128a) and using the constitutive functions (127a) for the chemical potentials provides an implicit representation of the mole fractions, viz.
where the mole fractions y ∞ α and pressure p ∞ at infinity are specified by bulk values of the electrolyte. Using the incompressibility constraint (126a), n F can be expressed as a function n F (ϕ, p). For given boundary values of ϕ and p, the equilibrium state of the double layer is determined by the coupled system The solution of the system (130) provides the spatial profiles of the ionic number densities in the double layer, cf. Figure 3. We observe that, in front of an electrode with a potential larger than in the electrolyte bulk, anions accumulate and for sufficiently large voltages the anion concentration saturates at a certain level that is related to the specific volume of the anions. Assuming that the anions have the same volume as the solvent, the saturation level is given by the number density of the pure solvent. In comparison, solvated ions are much larger then the pure solvent molecules and thus the saturation level is much lower and the saturated zone is much wider. Accordingly, the region where the cations are depleted is much wider for the solvated ions. In contrast, classical Nernst-Planck/Poisson-Boltzmann theory is based on the assumption of a dilute solution and representations of the number densities as Poisson-Boltzmann: where 0 denotes the solvent. Compared to representation (129), here the dependence on the local pressure is missing and instead of the coupled system (130), only the single Poisson-Boltzmann equation −(1 + χ)ε 0 ∆ϕ = n F (ϕ) needs to be solved. However, Figure 3 demonstrates that already for moderate or even small applied voltages, classical Poisson-Boltzmann theory leads to unphysical, almost unlimited accumulation of ions, far beyond the number density of the pure solvent and thereby it violates the underlying strong dilution assumption. These limitations and inconsistencies of classical Poisson-Boltzmann theory are well known and we refer to [12] for a review of the extensions made in the literature. Our approach based on non-equilibrium thermodynamics provides an additional major advantage over the classical theory and its extensions in the literature. From system (128) and the Gibbs-Duhem equation (78), we recover the momentum balance in equilibrium The Lorentz force k = −n F ∇ϕ is balanced by the pressure gradient. Thus, pronounced pressure gradients have to be expected in charged boundary layers that screen the electric field. Figure 3 depicts the spacial profiles of potential and pressure in the boundary layer. We observe drastic increase of the pressure towards the electrode surface, where p grows up to several GPa. Although the solvated ion case shows considerably smaller pressure than the case of all particles having the same size, the pressure at the electrode is still orders of magnitude larger than atmospheric pressure of 0.1 MPa. However, one has to keep in mind that the mechanical stress, which is controlled in experiments, and which can be measured, is the total stress (38). In the computations above, this total stress equals the outer pressure of 0.1 MPa everywhere in space. Also for the classical Poisson-Boltzmann model, the local pressure can be calculated from (132). For the so computed pressure, we observe in Figure 3 a growth near the interface that is orders of magnitude stronger than in our model. In the one-dimensional setting we use here, a combination of Equations (128b) and (132) yields the relation p − p ∞ = 1+χ 2 ε 0 |∂ x ϕ| 2 . This reveals that in the Poisson-Boltzmann case, the almost unlimited accumulation of ions that shield the electric field at the surface leads to a steeper slope of the potential and thus to higher pressure, compared to our model with the solvated ions.
The electric charge of the boundary layer can be defined in the one-dimensional setting as A remarkable relation from ( [18], Equation (150)) states This directly shows that classical Poisson-Boltzmann theory predicts the storage of an unrealistic huge charge in the boundary layer while our improved models considerably reduce the stored charge and the solvated ions approach can provide quantitatively meaningful results.

Double Layer Capacity
In addition to the charge accumulation in the double layer, also adsorption from the electrolyte to the surface and electron transfer reactions may take place, depending on the electrode metal and the ionic species. In a similar way like above, a surface charge Q S can be defined as function of the surface particle densities n s α . The double layer charge is then Q = Q BL + Q S . In equilibrium, the solution of the boundary value problem for the double layer, in particular the double layer charge, becomes a function of the applied voltage U between electrode and electrolyte. This allows for the introduction of the differential capacity or double layer-capacity [9,10,49] as As the charge Q, C also depends on the salt concentration, the temperature and the parameters of the free energy like specific volumes and adsorption energies. The double layer capacity can be measured directly and provides specific characteristics of the electrolyte and electrode-electrolyte interface at hand [50][51][52]. In Figure 4, measured and calculated double layer capacities are depicted for different salt concentrations. All simulations are done with the same set of parameters, independent from the salt concentration. For further details, we refer to [18]. The simulation shows the typical camel shape of the double layer capacity of aqueous electrolytes and a transition from a two-maximum curve for low salt concentrations to a one-maximum curve at high salt concentrations. Altogether, remarkable agreement between simulated and measured double layer capacities is reached.  [50] with permission from Elsevier). Right: comparison of the computed capacity (reprinted from [18] with permission from Elsevier).

Electrocapillarity-Lippmann Equation
Electrocapillarity describes the relationship between the interfacial tension γ and the applied voltage U. In a thin capillary tube, it can be observed that when a potential difference U is applied, a mercury-aqueous electrolyte interface moves according to the pressure difference while the curvature k M of the interface remains almost constant. The interfacial tension γ can be determined by the Young-Laplace equation, Moreover, the Lippmann equation [9,10,49] relates the slope of the surface tension with respect to the applied voltage to the double layer charge. When discussing electrocapillarity in the context of our model, some subtle differentiation is necessary. The interfacial tension γ in Equations (136) and (137) must not be identified with the surface tension γ s of Equations (78) or (123). In the context of non-equilibrium electrothermodynamics, the surface balance of momentum (22b) does not simplify to the Young-Laplace equation (136) due to the non-zero electromagnetic field.
In [19], we showed that in the asymptotic limit of thin double layers the Young-Laplace equation as well as the Lippmann equation can be derived. However, instead of the surface tension γ s , these equations contain the newly defined interfacial tension Here, the two boundary layer contributions γ ± BL of the corresponding electrode and electrolyte phase are non-negative functions of the electric field in the space charge layers. Thus, charging of the double layer always causes non-negative contributions γ ± BL which lower the interfacial tension. This directly explains the U-or parabola-shaped electrocapillary curves which are observed in experiments-see Figure 5 for experimental measurements and comparison with simulation.  Figure 5. Left: measured electrocapillarity curves for various salts according to [53], Figure 1 (reprinted with permission from the American Chemical Society). Right: computed interfacial tension as described in [19], Figure 1 (reprinted with permission from Cambridge University Press).

Electrochemical Systems in Non-Equilibrium
Generalized Nernst-Planck Flux To derive an explicit representation of the partial mass fluxes, we substitute the chemical potentials (117) into the constitutive equations (59b). With the simplifying choice of a diagonal mobility matrix, where the diagonal elements are of the form M αα = M α T m α ρ α , we get generalized Nernst-Planck fluxes. Denoting the solvent in Ω by the index 0, they read in the incompressible case for the ionic species α ∈ I \ {0}, Compared to the standard Nernst-Planck model, cf. [9,10], there are three additional terms highlighted here in blue. The first of these three terms is due to the solvent-ion interaction. It originates from the construction of the entropy production (55) where we incorporated the constraint ∑ α∈I J α = 0. The second term results from the incompressibility constraint and the possibly different specific volume of the constituents. The third one contains the pressure contribution that is required for thermodynamically consistent fluxes. The pressure contribution only vanishes if the atomic masses and specific volumes of all species are equal. The classical Nernst-Planck equations are derived under the dilute solution assumption, cf. [54] and sometimes the additional assumption of locally electroneutral solution, cf. [55]. The validity of both assumptions cannot be guaranteed inside the double layer, as we have seen above. For a dilute solution, i.e., n α n 0 , the first two additional terms in the flux (139) can be neglected. Moreover, we will discuss next that, outside of the double layers, the electrolyte bulk can be considered locally electroneutral and in particular isobaric, causing also the third additional term in flux (139) to vanish. Thus, the generalized Nernst-Planck flux in a dilute electrolyte bulk reduces to the classical one. However, in the double layer, the difference in general is significant and even might dominate the overall behavior of the considered system.

Asymptotic Analysis and Reduced Models
The partial mass balance equations with the fluxes (139) have to be coupled to the total mass and momentum balance, the Poisson equation and the surface (jump) conditions, leading to a rather complex system that contains several strongly different scales. In order to derive simpler models, formal asymptotic analysis can be applied, cf. [16,56] for a survey of the literature. An important characteristic parameter is λ according to definition (93). It relates the Debye length λ 2 x ref = k B T ε 0 e 2 0 n ref , that describes the width of the double layer, to characteristic macroscopic length x ref of the system. For the asymptotic analysis, all dimensionless parameters are related to powers of λ. Then, any function is approximated by two different formal expansions with respect to λ, the outer expansion in the bulk and an inner expansion near the surface. For the inner expansion, the space coordinate is rescaled by λ. Figure 6 illustrates the asymptotic method. For each polynomial power in the parameter λ, the corresponding terms in the model equations are collected. Finally, both expansions are connected by matching conditions relating the boundary values of the outer expansion to the far field of the inner expansion. In [14,16,19], reduced models for the generalized Poisson-Nernst-Planck system are derived. They are characterized by following properties: • In leading order, the bulk domain is locally electroneutral and pressure is constant in the bulk.

•
The double layer is globally electroneutral and is in quasi-equilibrium. Thus, the results of the preceding section can be applied.

•
Boundary layer charge and surface charge are both quantities of first higher order in λ.

•
Analysis of the inner equations allows the formulation of new boundary conditions in terms of the outer variables. Figure 6. Generic field variable u for an electrochemical system of two substances, here metal and electrolyte, in contact with the surface S (left). Boundary layers due to a small parameter λ in the model equations and decomposition of u in bulk part u λ and boundary layer part u λ (middle). In a simplified model for λ → 0, modified jump conditions contain the relevant information of the double layer that is not resolved any more (right). Figure reprinted from [16].

Surface reactions-Butler-Volmer Equation
Surface reactions, e.g., electron transfer reactions, have been intensely studied by experiments and there is a strong empirical basis for a macroscopic relation in which the surface reaction rate R s is driven by a potential difference at the interface that is called the surface overpotential η S . This relation is known as the Butler-Volmer equation and is considered to be "the central equation in phenomenological electrode kinetics" ( [49], p. 1053). It can be written as, cf. [9,10,49], where R f /b are the forward and backward exchange rates, which may depend on the temperature T and bulk number densities n α . The transfer coefficients α f and α b are considered as constant phenomenological coefficients.
In non-equilibrium electrothermodynamics, we are confronted with the observation that the surface Maxwell's equations (102) require the electric potential ϕ to be continuous at an interface, at least in the electrostatic setting. Therefore, no natural potential difference exists, which could be used to define an overpotential η S . By means of the asymptotic analysis above, it is possible to derive a general Butler-Volmer equation in the context of electrothermodynamics [17]. The derivation relies on two necessary conditions: (i) the boundary layers behave quasi-static such that the electrochemical potentials µ α + z α e 0 m α ϕ are constant and (ii) all adsorption processes are fast compared to the surface reactions, i.e., M s α → ∞ in the constitutive equations (72c).
Rather simple representations of the quantities defined in (140) in terms of the electrothermodynamic quantities can be found in the case of a single surface reaction where set of relevant surface constituents can be restricted to the constituents of the bulk phases, i.e., I S \ I = ∅. In the thin double layer limit λ → 0, we have In the asymptotic limit λ → 0, we denote the thin double layer interface with its internal structure by I, in order to distinguish it from the physical surface S in the complete model (see Figure 6). The constant β ∈ (0, 1) is known as the symmetry factor, which fosters either the forward or backward reaction.
The transfer coefficients α f /b depend on β and the stoichiometric coefficients γ s α . The quantities supplied with an overbar represent equilibrium quantities that are defined by the Nernst equation In accordance with usual definitions in electrochemistry [10,49], the overpotential η S describes the deviation of the actual potential difference from the equilibrium voltage between the bulk phases. The exchange rates R f /b depend on temperature and in an indirect way via the bulk chemical potentials µ α on the bulk number densities.
To summarize: there exists a thermodynamic consistent basis for the Butler-Volmer equation and the thermodynamic origin is the constitutive relation (68) for arbitrary surface reactions. In contrast to the Butler-Volmer equation of electrokinetics, the general constitutive relation for surface reaction rates (68) can be applied to electrochemical systems where it is necessary to spatially resolve the electrical double layer and where rate limiting adsorption processes are involved. For a detailed derivation and analysis of general Butler-Volmer equations, we refer to [16,17].

Conclusions
Better theoretical understanding of many modern electrochemical applications demands extensions of classical continuum models. However, deriving such extensions in a thermodynamically consistent way is a non-trivial task, due to the coupled physical phenomena and the multi-scale nature of the considered systems. The derivation of the mathematical models greatly profits from the application of an entropy principle that restricts the modeling freedom. In the literature, several different variants or flavours of entropy principles can be found that often only differ in details. We chose an entropy principle oriented on [24] with the postulation of a specific structure of the entropy production as the sum of binary products. In this work, we focused on the coupling of electrodynamics and thermodynamics and the coupling of bulk and surface equations. We introduced classical balance equations for matter and Maxwell's equations for the electromagnetic field with emphasis on the analogies between bulk and surface.
We restrict the constitutive modeling by a symmetry principle. Due to the different transformation properties of the (non-relativistic) balance equations of matter and Maxwell's equations with the (1+3-dimensional) Maxwell-Lorentz aether relations for the electromagnetic field, the Galilean symmetry principle is chosen. By this choice, some relativistic effects are excluded, but typically these effects are negligible in electrochemical applications.
Another point one should draw attention to is the choice of independent variables for a specific considered material. This set of independent variables is not uniquely determined and different choices are possible. However, different choices imply different definitions of temperature, and they furthermore may lead to different stability properties, even though the entropy principle is satisfied, as we illustrated in Section 6.4. Thus, how can the appropriate set of independent variables be chosen? Typically, the choice is guided by experience and justified a posteriori by the derivation of the entropy production and the additional conditions for the relaxation to equilibrium.
Exploitation of the entropy principle reveals numerous cross relations between the derived constitutive equations. This is a consequence of having one pivotal ingredient to encode material behavior, i.e., the entropy density or the free energy density, respectively. In particular, bulk chemical potentials are defined as derivatives of the bulk entropy density and they appear e.g., in the bulk mass fluxes as well as in the adsorption conditions. Thus, these phenomena of a very different nature cannot be modeled independently.
How does the existing theory for charge transport in electrolytes and electrochemical processes at the electrode-electrolyte interfaces relate to the described general framework? As a first step to establish such relations, a dimensional analysis of general model is used to identify the relevant time and space scales for the application to electrolytes. This justifies substantial simplifications of the general model, in particular with respect to Maxwell's equations. Next, an explicit model for the free energy density for liquid electrolytes is derived. The resulting liquid electrolyte model generalizes many of the already existing improved Nernst-Planck models with finite ion size effects such as e.g., [11,12]. Finally, by applying formal asymptotic analysis to the charged boundary layer at the electrode, it is possible to recover fundamental equations of electrochemistry like Butler-Volmer equations or the Lippmann equation.
While in principle it is possible to extend classical models by constitutive equations different from the ones derived here, it is in general not trivial to prove or disprove thermodynamic consistency of the obtained models. On the contrary, by applying the framework developed in this paper, the resulting models are guaranteed to be thermodynamically consistent. Moreover, starting from the more general framework, we gain deeper insights into the structure of electrochemical interfaces and the coupling of different phenomena.

Outlook
Now that this very general modeling framework is developed, we are ready to apply it to many different electrochemical systems other than liquid electrolytes. To model e.g., solids and polymer electrolytes, appropriate material models have to be formulated in terms of a free energy density. In solids, one has to abstain from the simplifying assumption of the free energy being independent from the unimodular deformation gradient F uni and one has to incorporate lattice velocities which generate further side conditions on the mass fluxes. For polymers, one has to incorporate the chain length of the polymers into the entropy of mixing and their impact on the elasticity.
With the introduction of the Galilean symmetry principle, the 1+3-dimensional Maxwell-Lorentz aether relations cannot be satisfied, but Lorentz invariance would be needed instead. The remedy is the consequent application of a relativistic setting for the equations of matter.
Author Contributions: All authors contributed equally to this work.
Funding: This research received funding by Einstein Foundation Berlin within the MATHEON Project CH11 "Sensing with Nanopores". The publication of this article was funded by the Open Access fund of the Weierstrass Institute and the Open Access fund of the Leibniz Association.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Formal Solution of Surface Balances
Let I be some time interval and Ω an evolving domain that is intersected by the moving (singular) surface S into Ω + and Ω − . By υ and υ s , we denote the barycentric velocity in the bulk and on the surface, respectively. Given the vector fields F, G : I × Ω ± → R 3 , define the density n : I × Ω ± → R and the flux J : I × Ω ± → R 3 by n = div(F) and J + nυ = −∂ t F + curl(G) in Ω ± . To show the assertion, three identities are derived from (A3): First, (A3) 1 is differentiated with respect to time t. With the definition (4), we infer For the second identity, we define Then, we differentiate (A3) 2 with respect to the surface parameter u Σ Finally, (A3) 1 is differentiated with respect to the surface parameter u ∆ ,

Bulk
Let ρu + M · B, (ρ α ) α∈I ± be scalars, P, M vectors and F uni a tensor with respect to Galilean transformation. Let h be an (absolute) scalar, i.e., h is invariant with respect to Galilean transformations according to Section 4. Moreover, let h be a constitutive function such that h = h(ρu + M · B, (ρ α ) α∈I ± , F uni , P, M) . (A8) Then, there holds the following symmetry relation for i, j = 1, 2, 3 Here, δ, β and γ are the respective angles of rotation. Substituting O = R x into (A10) followed by differentiation with respect to δ yields Then, we conclude for δ = 0 that Applying the same arguments with the matrices R y and R z yields the assertion.