Galilean Bulk-Surface Electrothermodynamics and Applications to Electrochemistry

In this work, the balance equations of non-equilibrium thermodynamics are coupled to Galilean limit systems of the Maxwell equations, i.e., either to (i) the quasi-electrostatic limit or (ii) the quasi-magnetostatic limit. We explicitly consider a volume Ω, which is divided into Ω+ and Ω− by a possibly moving singular surface S, where a charged reacting mixture of a viscous medium can be present on each geometrical entity (Ω+,S,Ω−). By the restriction to the Galilean limits of the Maxwell equations, we achieve that only subsystems of equations for matter and electromagnetic fields are coupled that share identical transformation properties with respect to observer transformations. Moreover, the application of an entropy principle becomes more straightforward and finally helps estimate the limitations of the more general approach based the full set of Maxwell equations. Constitutive relations are provided based on an entropy principle, and particular care is taken in the analysis of the stress tensor and the momentum balance in the general case of non-constant scalar susceptibility. Finally, we summarise the application of the derived model framework to an electrochemical system with surface reactions.


Introduction
The reliable simulation of electrochemical systems on the device level requires thermodynamically consistent continuum models. Standard models such as the Poisson-Nernst-Planck system suffer from well-known limitations, cf., e.g., [1,2], of which the most-obvious deficiency is the missing volume exclusion effects. Several similar remedies for this problem have been proposed; see, e.g., [3][4][5] and the literature cited therein. More recent approaches that also take the dielectric effects of the solvent into account were given, e.g., in [6,7]. The latter one builds on a bulk-surface electrothermodynamic modelling framework summarised in [8], which unifies classical non-equilibrium thermodynamics containing electromagnetic fields as in [9,10] and its extensions to surfaces [11,12]. As an alternative modelling framework, Reference [13] recently developed a multiscale electrothermodynamics approach within the GENERIC formalism. Although not less complex and technical compared to [8], this approach is so far limited to closed bulk systems.
Applying the abstract framework of [8] to electrochemical problems typically first requires strong simplifications of the model, thereby favouring the impression that the framework appears too general and not precisely tailored to the considered class of problems. Moreover, the readers of [8] might be confused by some decisions made there for mainly pragmatic reasons, but leading to a discussion of the stability conditions that are related to the choice of independent variables and that are necessary in addition to the employed entropy principle. Finally, the way that the electromagnetic sub-system and the classical material mixture part are coupled in [8] might not appear fully satisfactory because of the different transformation properties of these sub-systems. The problem of different transformation properties can be remedied in two alternative ways. On the one hand, one can give up mass conservation in favour of a fully relativistic description of matter. On the other hand, one may try to couple a low-velocity Galilean limit system of the Maxwell equations to the classical, i.e., non-relativistic, balance equations of mixtures of matter.
In this work, we proceed with the second of the alternatives, since we aim at engineering applications and, in particular, electrochemical applications. As shown by [14], cf. also [15], there is not one unique Galilean limit of the Maxwell equations, but instead, there are two distinct limits, i.e., the electric-or quasi-electrostatic-limit and the magnetic-or quasi-magnetostatic-limit. While coupling the electric limit to the balance equations of mass, momentum, and energy seems most reasonable for electrochemical applications, the magnetic Galilean limit is more naturally related to magneto-hydrodynamics. We re-examine the derivation of Galilean limit systems of the Maxwell equations here, as we are not aware that this has been carried out before for the case of singular surfaces. This way, we can make sure that we do not miss relevant effects during coupling to the balances of matter later on or during the following derivation of constitutive equations. For each of the two Galilean limit cases considered in this work, there is one corresponding variant of the entropy principle in [8] such that coupling the general full set of Maxwell equations to the classical balance equations of matter and taking the Galilean limit after exploitation of the entropy principle yield the same results as obtained here. However, none of the two alternative variants of the entropy principle in [8] are capable of covering both Galilean limit cases equally well. It turns out that the application of the entropy principle here is much more straightforward compared to [8], as a re-definition of the combined inner energy of field and matter becomes unnecessary and a discussion of the stability of the polarisation relaxation does not arise.
The obtained constitutive equations differ between the two Galilean limit cases. As can be expected, polarisation relaxation is only covered in the quasi-electrostatic case, whereas the relaxation of magnetisation is only present in the quasi-magnetostatic case. Moreover, the generalised driving forces for thermodiffusion differ between the Galilean limits, as does the stress tensor.
We apply the geometrical setting and notation of [8], but neglected the elastic deformations, since these effects do not contribute to the questions discussed here. The results obtained here can be generalised in a safe and straightforward manner to include elastic effects in the same way as in [8].

Notation
We consider a geometrical setup where, locally, an orientable surface S divides a domain Ω ⊂ R 3 into two subdomains Ω ± ⊆ R 3 with S := ∂Ω + ∩ ∂Ω − . The domain Ω, as well as the surface S may evolve in time. In addition to quantities defined in the domains Ω + or Ω − , there are, in general, corresponding quantities on the surfaces S, which, in general, do not coincide with the corresponding traces from the subdomains. As a convention, the same letters are used for these quantities, but the surface variables are indicated by an underset s.
We apply the convention of implicit summation over coordinate indices appearing twice. We indicate the Cartesian components of vectors and tensors by lowercase Latin indices, e.g., i, j, whereas at the surface S, we use uppercase Greek indices such as, e.g., Γ, ∆ for the tangential components.
The tangential vectors are defined as the partial derivatives of the smooth bijective parametrisation θ mapping from some open parameter domain to the surface S. In addition, we define the (unit) normal vector and the metric tensor by As a convention, we chose the mapping θ such that ν is the inner normal of Ω + . For the matrix components of the metric tensor g, we use lower indices g ∆Γ , and we use upper indices g ∆Γ for the components of the inverse matrix of the metric. A vector V defined on the surface can, thus, be written as V = V ∆ τ τ ∆ + V ν ν, with the normal component by V ν and the tangential components V ∆ τ = g Γ∆ V τ Γ , for ∆ = 1, 2. 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 Γ∆ . 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 υ s denote the velocity of the surface S. For a scalar a : [0, t end ) × S → R, we define the time derivative: Let u be a generic function defined at least in one of the subdomains Ω ± . We denote the trace of u by whenever u is defined on this side of the surface; otherwise, we set the corresponding trace to zero. We define the jump and the mean value of u at the surface S by

Maxwell Equations and Their Galilean Limits
The Maxwell equations for the electromagnetic field are most generally formulated such that time and space and all involved equations are properly combined into fourdimensional objects. Then, the equations satisfy the most-fundamental symmetry principle, i.e., the principle of relativity, meaning that the balance equations, as well as the constitutive equations remain invariant with respect to arbitrary observer transformations.
Here, 1+3-dimensional Maxwell equations are formulated in a way following the classical work of Truesdell and Toupin [16]. This includes the postulation of Maxwell-Lorentz aether relations that are independent of the considered material. The precise form of these relations depends on the frame of reference. Here and in the following section, we assume an inertial frame of reference such that the Maxwell-Lorentz aether relations take the most-simple form (9) below. In the 1+3-dimensional setting, Maxwell-Lorentz aether relations are only invariant with respect to Lorentz transformations. The Galilean transformation is a good approximation of the Lorentz transformation in the limit of vanishing barycentric velocity, i.e., v/c 0 → 0. However, the derivation of a Galilean limit system of the Maxwell equations is complicated by the fact that there are two different limit systems, cf. [14,15].

General Maxwell Equations
The underlying physical principle for the derivation of the Maxwell equations is the conservation of charge and magnetic flux. The (total) electric charge density n e in each volume domain Ω ± , a well as n s e on the surface satisfy the local conservation equations, viz. ∂ t n e + div(n e υ + J e ) = 0 , (7a) Here, we split the electric current as j e = n e υ + J e . We introduce the charge potential D and the current potential H by means of a formal solution of the charge balances (7); see [8] Therefore, we have The conservation of magnetic flux in the bulk and on the surface reads where E is the electric field. We postulate universal valid Maxwell-Lorentz aether relations, i.e., independent of the considered material, it holds in an inertial frame: where the dielectric constant ε 0 and the magnetic constant µ 0 are related to the speed of light by ε 0 µ 0 = c −2 0 . The remaining constitutive quantity is the electric flux density J e . As an immediate consequence of the universal valid Maxwell-Lorentz aether relations (9) and charge conservation (7a), we conclude that the sum of the electric current and the displacement current is at any time t and location x source-free. This motivates the definition of the total current as Implied Balances From Maxwell's Equation (8), the following two additional balance equations can be derived, cf. [10]. The balance of electromagnetic momentum reads in each subdomain Ω ± , as well as on the surface S where the right-hand sides are the negative Lorentz force in the volume and on the surface, respectively. The Maxwell stress tensor in (11) is where 1 denotes the identity matrix. The balance of electromagnetic energy reads where the right-hand sides are the negative Joule heat in the volume and on the surface, respectively.

Polarisation and Magnetisation
When considering a macroscopic non-relativistic continuum description of mixtures of charged matter, there is the free charge density n F defined as a Galilean conserved quantity. In general, this free charge does not coincide with the previously introduced electric charge n e . As an example, n F may solely represent the net charge of molecules in a mixture without taking into account the internal electronic structure. However, also this internal electronic structure in general has relevance for the overall electromagnetic field. To bridge this gap and represent these microscopic effects on the more macroscopic level, we introduce the polarisation charge density n P and the polarisation current density J P as n e = n F + n P , The conservation of the electric charge n e and of the free charge n F then also imply the conservation of the polarisation charge n P = n e − n F . We introduce the vector of polarisation P and Lorentz magnetisation M by the same approach of the formal solution to the conservation equations, viz.
Not all of the introduced quantities so far are Galilean scalars, vectors, or tensors, respectively. Due to its importance in electrodynamics and for the upcoming constitutive modelling, we introduce the Galilean vectors of the electromotive intensity and the magnetisation: The electric current can then be expressed as For later use in the balance of the inner energy, we verify that Upon the definition of the dimensionless constants: we obtain in the volume domains Ω ± the following dimensionless system: For the surface, we introduce an additional reference number density n Then, the surface balances are in dimensionless form: The balance of electromagnetic momentum reads in the volume domains and on the surface Here, the Maxwell stress tensor multiplied by β 2 , viz.
contains terms independent of β 2 , as well as the Lorentz force and the Joule heat do, i.e., The electromotive intensity and magnetisation arȇ The balance of electromagnetic energy reads and the introduction of polarisation and Lorentz magnetisation implies, for the electric current and the surface electric charge flux, We consider in the following two alternative limiting cases for the parameter β, whereas we only assume for the remaining parameters λ, κ, ω, and δ that they remain moderate in size. At this point, we only remark that, later on, it is possible to analyse additional limit processes, such as, e.g., λ → 0, which corresponds to the thin interface limit in electrochemical applications.

Quasi-Electrostatic Limit β 2 1
We assume that β is a small parameter, whereas v ref c 0 and the time derivative of B remain bounded. As a consequence, we obtain the equations of electrostatics for the electric field, which read in the dimensional form A constitutive equation for P, which is independent of B and M, then allows for a given free charge obtaining E and P from (29). The balance of electromagnetic momentum reduces to with the Maxwell stress tensor depending only on E, viz.
In this limit, the electromotive intensity is identical to the electric field. The magnetic field is not vanishing; it does not have any effect on the Lorentz force, viz.
The Joule heat cannot be evaluated from E and P alone, but also requires knowledge of B.
The remaining Maxwell equations, in general, together with a constitutive equation for M, determine the magnetic flux density B, i.e., The electric current and the surface current flux are then

Quasi-Magentostatic Limit
Under the assumption that the derivatives ofȆ +P remain bounded, Maxwell's equations simplify in this limit to A constitutive equation for M, which is independent of E and P, then allows for a given free current to obtain B and M from (35). The balance of electromagnetic momentum reduces to with the Maxwell stress tensor depending only on B, viz.
In the quasi-magnetostatic limit, the magnetisation is identical to the Lorentz magnetisation. The electric field does not vanish; it does not have any effect on the Lorentz force, and the Joule heat is negligible, viz.
The Lorentz force cannot be evaluated from B and M alone, but also requires knowledge of E. The remaining Maxwell equations, together with a constitutive equation for P, determine the electric field, i.e., The electric current and the surface current flux are

Balance Equations of Galilean Electrothermodynamics
The Galilean limit systems of electromagnetics can now be consistently coupled to classical, i.e., non-relativistic balance equations for charged mixtures 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.

Description of Reacting Mixtures
We use different index sets I ± to refer to the constituents of a mixture in Ω ± and the index set I S for the constituents on the surface S. We apply the non-restrictive assumption that the sets I ± are disjoint, i.e., I + ∩ I − = ∅. All constituents of the subdomains Ω ± are assumed to be also constituents on the surface S, but there may be some additional constituents that are exclusively present on S. Thus, we have I ± ⊆ I S .
There may be several chemical reactions among the bulk constituents, as well as chemical reactions on the surface. Picking some indices k, , then, these reactions in the bulk or surface may be written in the general form: where A α is used as a notation to refer to the different constituents. The constants a k α , b k α are positive integers, and γ k α := b k α − a k α denote the stoichiometric coefficients of the reactions. The net reaction rate is R k , where reactions in the direction from left to right are counted as positive.
Each constituent has the (atomic) mass m α , for α ∈ I ± or α ∈ I S , and the net charge z α e 0 , where z α is the charge number of the constituent. Since charge and mass have to be conserved by each single reaction in the bulk and on the surface, we have To describe the thermodynamic state of the mixture, we use in the volume domains the inner energy density ρu, the particle number densities n α , and partial velocities υ α for α ∈ I ± . On the surface, the surface inner energy density ρ s u s and the number densities and partial velocities of the surface constituents are n s α and υ s α for α ∈ I S . The multiplication of the number densities by m α gives the partial mass densities: The mass density and the barycentric velocity of the mixture are defined by The (non-convective) bulk and surface diffusion flux J α and J s α with respect to the barycentric velocity are defined as The free charge density and free current are are defined as In each of the subdomains Ω ± , as well as on the surface S, the partial mass balances are The partial mass balances can be combined to derive bulk and surface conservation laws for the total mass density of the mixture and for the free charge.

Balance of Momentum
In the absence of electromagnetic fields, the momentum density of the mixture with respect to the barycentric velocity is ρυ. We postulate that, in the absence of gravitation, the total momentum of matter and the electromagnetic field is a conserved quantity. In the Galilean limits of the Maxwell equations, the Maxwell stress σ EM equals the negative Lorentz force, while the electromagnetic momentum density is negligibly small. We introduce the total stress tensor Σ containing σ EM . Then, in either of the Galilean limit cases, the total momentum balances of the matter and electromagnetic field read in Ω ± and on S where f and f s are due to gravitation. The surface momentum flux σ s is decomposed into its normal and tangential components: The tensor with the components S Γ∆ is denoted as the surface stress tensor, and the vector with the components S ∆ is the normal stress vector. We neglect internal spin. Given the symmetry of the Maxwell stress tensor due to the Maxwell-Lorentz aether relations, this implies the symmetry of the total stress tensors and the vanishing of the normal surface stress [10], i.e., the assumption implies, for i, j = 1, 2, 3 and Γ, ∆ = 1, 2, From these balances, several different balances for the momentum ρυ can be derived, depending on which forces should be accounted for on the right-hand side. Subtracting from (48) the momentum balance (30) in the quasi-electrostatic limit, or (36) in the quasimagnetostatic limit, then yields the momentum balance equations of matter with the Lorentz force k, respectively k s , on the right-hand side. In addition to the Lorentz force, one might also want to track the Kelvin polarisation force, electrostrictive force, or Korteweg-Helmholtz force explicitly in the balances (cf. Ref. [7] ), such that the stress tensor then needs to be modified accordingly.

Balance of Inner Energy
The energy of matter consists of the inner energy density ρu and the kinetic energy density ρ|υ| 2 . The balance equations of the mass and momentum of matter imply the balances of the kinetic energy in Ω ± and on S. Moreover, we postulate that, in the absence of gravitation, the total energy of the matter and electromagnetic field is a conserved quantity. This implies the inner energy balances in Ω ± and on S as shown in Appendix A, viz.
We can write Q = q + E × M, such that q coincides with the heat flux in the balance of the inner and kinetic energy when the heat production is given by the Joule heat. Using the particular form of the Lorentz force and the Joule heat in the quasi-electrostatic limit, we obtain for the inner energy balance On the other hand, we obtain in the quasi-magnetostatic limit from the explicit Lorentz force and Joule heat

Constitutive Equations for the Quasi-Electrostatic Limit
In the quasi-electrostatic limit, the coupled electrothermodynamic system consists of (29), (47), (48), and (52). To close the coupled system, additional constitutive relations are needed. These can be obtained from the application of an entropy principle as described in [8]. Because it allows a straightforward derivation of the entropy production in the desired form, we assume constitutive functions of the entropy densities of the form: In principle, it would also be possible to include the dependence of ρ η on M, such that, by the application of the entropy principle, one more constitutive equation could be obtained, which then would allow solving (33) and evaluating the electric current j e . In the quasi-electrostatic limit, however, all explicit dependence of the system (29), (47), (48), and (52) on B and M is removed due to the asymptotic smallness of those terms. The application of the entropy principle builds on this system of equations, and therefore, it should not reintroduce asymptotically non-negligible dependencies on B and M. Moreover, the terms in the entropy production depending on B and M would be in the same asymptotic order as the terms that have been neglected in the system (29), (47), (48), and (52). Therefore, the value of the constitutive equations based on such small terms in the entropy production must be questionable. Instead, it seems more advisable to use a constitutive equation relating B and M, which were derived from another source, which might well be the complementary quasi-magnetostatic limit in Section 6. The (absolute) temperature T, T s and the chemical potentials µ α , µ s α are defined as For the free energy density and the entropy density, we, thus, have representations as functions of the temperature as and in an analogous way for the inner energy and entropy. From the construction (56) and the definition of the temperature according to (55), we obtain the well-known thermodynamic relations:

Constitutive Relations for the Bulk
To determine the entropy production in a straightforward manner, the chain rule is applied to the time derivative of ρ η, and the partial mass balances (47) and the inner energy balance (52) are applied. Then, the choice of the entropy flux fixes the entropy production. Nevertheless, the entropy production can still be rewritten in such a way that is most suitable for the derivation of the constitutive relations. We choose the entropy flux as φ = Q T − ∑ α∈I µ α T J α . Moreover, we define the viscous stress tensor and the symmetric velocity gradient as To account in the derivation of the entropy production for the constraint (45) on the mass fluxes, we choose in each subdomain one designated constituent as a reference and denote these constituents by A 0 − , A 0 + , respectively. Moreover, we use a symmetry condition originating from the transformation properties of the thermodynamic fields, cf. Ref. [8]. Then, we obtain the entropy production as the sum of binary products, viz.
We now can identify the dissipation mechanisms related to their specific entropy production.
To reduce the complexity of the constitutive equations, we neglect cross-effects other than thermodiffusion. Thus, we consider the five dissipation mechanisms: shear viscosity ξ SV , volume viscosity ξ VV , (bulk-)reactions ξ R , thermodiffusion ξ TD , and polarisation ξ P . For each dissipation mechanism, we apply closure relations to obtain the constitutive equations.

Thermodiffusion
For the heat flux and the mass fluxes, we choose a linear relation with cross-effects. We chose coefficients κ, L β , M αβ such that the coefficient matrix is symmetric and positive definite. In particular, the heat conductivity κ and the mobility matrix M are symmetric and positive definite. We set for α ∈ I ± \ {A 0± }

Reactions
For simplicity, we neglect cross-effects between the different chemical reactions. We choose positive coefficients A k , R k 0 and and apply a nonlinear closure relation to obtain, for any of the reactions in the volume domains,

Polarisation
We chose a linear closure relation with the relaxation time τ P ≥ 0 of polarisation to obtain

Viscosity
Linear relations for the volume viscosity and for the shear viscosity with phenomenological coefficients η b , η s that satisfy (η b + 2 3 η s ) > 0 and η s > 0 yield This implies the constitutive equation for the symmetric stress tensor:

Constitutive Relations for the Surface
We proceed in analogous manner as in the volume and choose the surface entropy flux as φ We define the tensors T and D with the components: Moreover, we choose one designated constituent A 0 ∈ I S on the surface. Then, we obtain the surface entropy production as a sum of binary products, viz.
We identify six dissipation mechanisms: tangential surface viscosity ξ

Thermodiffusion
We chose phenomenological coefficients κ s , L s β , and M s αβ , such that the coefficient matrix is symmetric and positive definite. We set, for α ∈ I S \ {A 0 },

Surface Reactions
We neglect cross-effects between the different reactions and apply the nonlinear closure relations with positive phenomenological coefficients A s , R s 0 :

Surface Viscosity
Linear closure for the trace and for the deviatoric part of the surface stress tensor This implies, for the symmetric surface stress tensor, the constitutive equation:

Heat Flux Normal to the Surface
We choose coefficients κ s ± and obtain from the linear closure:

Mass Transport Normal to the Surface
The mass transport normal to the surface can be related to the adsorption, which, in the context of electrothermodynamics, is formally different from the mass production due to chemical reactions. However, in experiments, it might not always be clear how to distinguish between adsorption and surface reactions, and we, therefore, here, also applied a nonlinear closure similar to the chemical reactions. We neglect cross-effects and chose for α ∈ I ± \ {A 0± } coefficients β ± α and M s ± α > 0 to obtain, for α ∈ I ± \ {A 0± }, (74b)

Discussion of a General Free Energy Model
The obtained constitutive equations are general, in the sense that they were 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. We now restrict the constitutive function for the free energy to the structure: where the first part is independent of the electromagnetic field and the second part depends only on the absolute value of P. Different free energy models can be used for the fieldindependent free energy contribution ρψ, depending on the considered material. For liquid electrolytes, the model of simple mixtures of solvated ions (cf. [5,17]) appears to be appropriate. The general construction of suitable free energies and their incompressible limit is analysed in [18]. The analysis emphasizes the importance of linear approaches with respect to the partial molar volume and illustrates the possibility of recovering the nonlinear behaviour of the mixtures based on the nonlinear reaction behaviour. As an example of a different material model, we refer to [19], which models yttria-stabilised zirconia.

Dielectric Susceptibility and Debye Relaxation
Since it might appear more familiar to formulate constitutional equations in terms of E instead of P, we introduce a conjugate variable to P which coincides with E in the case of polarisation equilibrium. We define Then, we consider a change of variables by means of the Legendre transformation: such that we have the following representation and relations: We introduce the scalar dielectric susceptibility as such that, by (78), this implies, for the vector of polarisation, that Considering the case of vanishing velocity, i.e., υ = 0, the constitutive Equation (63) reduce to where the relaxation constant τ P > 0 is positive. From div(ε 0 E + P) = n F , we conclude that, for the divergence of the constitutive Equation (81), it holds We conclude that, if n F = 0 in Equation (82), then div(P) vanishes for t → ∞, which is compatible to stable polarisation relaxation. We remark that the stability of relaxation does not depend on the applied change of variables.

Pressure, Surface Tension, and Stress Tensor
In the context of electrothermodynamics, the proper definition of pressure is not obvious. We define here pressure based on only these contributions, which are also present in the absence of electromagnetic fields. By the Gibbs-Duhem relation known from classical thermodynamics and its counterpart on the surface (cf. [9,10]), we define pressure and surface tension aŝ where the splitting of the free energy in (75) implies a similar splitting of the chemical potentials. The Legendre transformation (56) implies With the above definitions of pressure and surface tension, we can rewrite the stress tensor in the volume in (65) and on the surface in (71) as Using the scalar susceptibility χ according to (79) and assuming polarisation relaxation equilibrium, i.e., E Eq = E, such that P = χε 0 E, according to (80), we infer with the Maxwell stress tensor according to (31) that Due to the structural similarity to (83), one might refer to the last bracket in (86) as polarisation pressure. Nevertheless, we prefer to keep this term separate from the pressure contributionp to the stress tensor and relate the last terms in (83) to forces in the momentum balance. In the particular case of constant χ, we check ρψ pol = − 1 2 χε 0 |E| 2 such that the stress tensor further simplifies to Σ = −p 1 + T + (1 + χ)σ EM .

Application to Momentum Balance
We again assumed E Eq = E and P = χε 0 E. The divergence of the Maxwell stress tensor σ EM equals in the quasi-electrostatic limit the Lorenz force n e E, and it seems appropriate to treat also the divergence of the last two terms in (86) as forces in the momentum balance. Whenever χ is independent of the concentrations, i.e., independent of all ρ α for α ∈ I, we apply ∇ρψ pol = −(∇E Eq ) P = −χ(∇E) ε 0 E, to conclude, for the divergence of the stress tensor in the momentum balance, if χ is independent of all ρ α : div(Σ) = −∇p + div(T ) + div ((1 + χ) Due to its similarity, the term n F E is also referred to as the Lorentz force. On the other hand, if ρψ pol is a homogeneous function of degree one with respect to all ρ α for α ∈ I, then χ is a homogeneous function as well, and the terms in the brackets in (86) cancel. We conclude if ρψ pol is a homogeneous function of ρ α : div In addition to the Lorenz force, the last term in (88) is referred to as the Kelvin polarisation force. In the momentum balance on the surface, we have

Constitutive Equations for the Quasi-Magnetostatic Limit
The construction of constitutive relations for the quasi-magnetostatic limit follows largely the same lines as described in the previous section. The entropy principle is now based on constitutive functions of the entropy density in the form: The (absolute) temperature T, T s and the chemical potentials µ α , µ s α are defined in an analogous way as in (55), and we introduce the free energy in an analogous way as in (56). More specifically, we restrict the constitutive function for the free energy to the structure: such that the splitting of the free energy implies a similar splitting of the chemical potentials. As in (83), we definê

Constitutive Relations for the Bulk
Taking the same entropy flux as above, the entropy production mostly looks the same as in (60), and we identify, again, five dissipative mechanism, where only polarisation is now replaced by magnetisation. We proceed in the same way as in Section 5 and apply the closure relations to each dissipative mechanism. For the reactions, there is no change to (62).
With the definition of pressure according to (92), the free energy as in (91) with the scalar magnetic susceptibility according to (96), and assuming magnetisation relaxation equilibrium, i.e., M Eq = −M, we obtain where, in the case of constant χ M , we have ρψ mag = − χ M 1+χ M 1 2µ 0 |B| 2 and the stress tensor further simplifies to Σ = −p 1 + T + 1 1+χ M σ EM . The divergence of the Maxwell stress tensor σ EM equals in the quasi-magnetostatic limit the Lorenz force (n e υ + J e ) × B, and it seems appropriate to treat also the divergence of the last two terms in (101) as forces in the momentum balance.

Constitutive Relations for the Surface
We use the definitions of T s and D s as in (66) to derive the surface entropy production in an analogously way as before. For the surface reactions, there is no change to (69), and for the viscous stress tensor on the surface, there is no change to (70).

Thermodiffusion
With the same choice of the symmetric and positive-definite coefficient matrix, we set, as in (68) for α ∈ I S \ {A 0 }, (102b)

Summary and Discussion
By asymptotic considerations, we motivated the formulation of two distinct Galilean limit systems of the Maxwell equations. Both Galilean limits are well known for volume domains and, here, were transferred to surfaces. Due to their transformation properties, each one of these limit systems can easily be coupled in a consistent manner to standard non-relativistic balance equations of matter that are well suited for most engineering applications. Then, constitutive equations to close the coupled system can be obtained from an entropy principle.

Comparison to the Results of [8]
Considering first the quasi-electrostatic limit, the same results as in Section 5, in principle, can also be been obtained by [8], where the full Maxwell equations are employed, and the quasi-electrostatic limit, then, can be taken after the application of the entropy principle. However, this procedure is less straightforward, because in [8], the inner energy density of the coupled system is redefined in order to obtain the entropy production related to magnetisation as a binary product satisfying a Galilean symmetry principle. Such a modification of the inner energy as ρu + M · B seems admissible because, in coupled electrothermodynamics, it is not a priori evident what is the correct definition of the inner energy. The non-dimensionalisation of Section 3.3 reveals that this modification made to the inner energy is asymptotically small in the quasi-electrostatic limit and can, therefore, be neglected.
An alternative application of the entropy principle based on E as the independent variable was discussed in [8]. It relies on a different definition of the inner energy as ρu − E · P, and therefore, it is fundamentally different from the change of variables employed here in Section 5.3. The term E · P is not asymptotically small compared to ρu in the quasi-electrostatic limit and, thus, implies a non-negligible difference in the definition of the temperature. Moreover, the application of a linear closure relation then yields an unstable relaxation of polarisation. In contrast, the change of variables applied here in Section 5.3 and in the same way also in [7], changing the variables from P to E Eq , does not alter the stability of polarisation relaxation, and in equilibrium, we have E Eq = E.
In addition, the entropy principle in [8] also covers magnetisation. We omitted magnetisation here in Section 5, because, in the quasi-electrostatic limit, the corresponding contributions to the entropy production are in the asymptotic order of the terms that we neglected in the balance equations. Without a constitutive equation relating B and M, we are in general not able to evaluate the electric current j e , although, in many experiments, this is the observable quantity of pivotal importance. Therefore, either an additional constitutive law has to be obtained from a different source, or we have to apply an additional asymptotic limit such as the thin double-layer limit, where the contribution to j e due to M vanishes. The relaxation of magnetisation was not discussed in [8]. We note that an analogous approach as performed here in Section 6, but applied to the main approach of [8] based on the variables (P, M), would imply the blow up of curl(M) for a diamagnetic material.
In the quasi-magnetostatic limit, we assumed in Section 6 the entropy density to depend on the independent variable E. This is different than the main approach of [8], but compatible with the alternative approach discussed there and with [10]. Within the quasi-magnetostatic limit, the definition of the inner energy as ρu − E · P is a negligible modification, whereas M · B is not asymptotically small. Accordingly, we observed in Section 6 the stable relaxation of magnetisation, whereas the approach based on the variable M would lead to unstable relaxation for a diamagnetic material.
Each of the modifications of the inner energy in [8] implies according changes to the stress tensor when it is expressed in terms of free energy. Non-constant susceptibility requires careful treatment of the momentum balance; cf. [7].

Conclusions for Coupling with General Maxwell Equations
Coupling the full Maxwell equations to the non-relativistic balance equations of matter appears to be a more general approach than the one presented here. For this general case, two alternative formulations of the entropy principle were analysed in [8], one based on the variables (P, M) and the other one based on the variables (E, B). While the (P, M)variant allows reaching the quasi-electrostatic Galilean limit where only P remains as a variable, the (E, B)-variant allows reaching the quasi-magnetostatic Galilean limit with only B as the remaining variable. Neither of the two variants is capable of covering both limit cases equally well, as for the (P, M)-variant, magnetisation relaxation is unstable for a diamagnetic material, whereas, for the (E, B)-variant, polarisation relaxation is unstable. However, this should not be considered a source of inconsistency for the limit cases, as the two Galilean limits do not have any relevant overlap. We remark that building the entropy principle on (P, B) as independent variables does not provide a remedy for the general case, since P and B do not form an antisymmetric four-tensor in the four-dimensional formulation of the Maxwell equations. Therefore, we conclude that, if the electric and magnetic effects both are relevant in the application, then one is well advised to consider, in addition to the general Maxwell equations, also a relativistic description of the mixture of matter. This leaves as an open question how to find appropriate models for, e.g., liquid metal batteries or magneto-hydrodynamic forces in the electrolyte transport of the Hall-Heroult process.

Electrochemical Model in Polarisation Equilibrium
We summarise here the complete system of the model equations for electrochemical applications under the simplifying assumption of fast polarisation relaxation. An important feature of electrochemical systems is the formation of double-layers at the contact of different materials. The double-layer is characterised by a typical width in the range of nanometres, wherein the electric potential may vary in the order of one Volt. The magnetic field strength is assumed to be below a guideline value for electromagnetic fields in electrical household appliances. These reference values: imply, for the dimensionless quantity, The smallness of β 2 suggests the use of the quasi-electrostatic limit, where (29b) implies the existence of an electrostatic potential ϕ, such that where the constant is frequently chosen as zero.
For each volume domain Ω ± , we use index sets I ± for the constituents of the mixture. The constitutive equations in the volume are built on the free energy density of the structure: The chemical potentials for α ∈ I ± , pressure, and the inner energy are given in terms of the free energy as The dielectric susceptibility is defined as Since polarisation relaxation in liquid electrolytes is typically fast with a time constant in the range of 10 −8 s (cf. [20]), we can assume that polarisation relaxation is in quasi-equilibrium, i.e., P = −χε 0 ∇ϕ. In the following, we distinguish two cases: χ is independent of ρ α , indicated as ( * ) , homogeneous of degree one in ρ α , indicated as ( * * ) . (110) On the surface with the index set I S for the constituents, the free energy density does not depend on the electric field. The chemical potentials for α ∈ I S , surface tension, and the inner energy are given in terms of the free energy as We choose a constituent A 0 ∈ I S and let α ∈ I S \ {A 0 }. Then, the surface balance equations are  Moreover, there are the boundary conditions: 0 , ( * * ) .
Finally, we remark that the evaluation of the electric current, in general, also requires knowledge of the magnetic field, although the magnetic field does not have any influence on the solution of the above system. In order to compute B for a given solution of the above model equations, we assume that magnetisation relaxation is fast and, hence, apply the equilibrium relation obtained in the quasi-magnetostatic case, viz. M = χ M 1+χ M · 1 µ 0 B. Then, we have to solve div(B) = 0 , and can evaluate the electric current as For a paramagnetic and a diamagnetic material, typically, |χ M | 1, and thus, the Lorentz magnetisation is small, such that it often might be appropriate to neglect the magnetic term, making the evaluation of the electric current j e available already without the solution of (117). Subtracting (A1) from (A2) and setting q EM = E × M then yield the inner energy balance (51). The choice of q EM is motivated by the fact that it makes the production terms in the inner energy balance directly computable within either of the Galilean limits of the Maxwell equations, as demonstrated below, without having to rely on the additional remaining Maxwell equations.