Theory of Transport in Highly Concentrated Electrolytes

Ionic liquids are promising candidates for novel electrolytes as they possess large electrochemical and thermodynamic stability and offer a high degree of tunability. As purely-ionic electrolyte without neutral solvent they exhibit characteristic structures near electrified interfaces and in the bulk, both being described theoretically via separate frameworks and methodologies. We present a holistic continuum theory applying to both regions. This transport theory for pure ionic liquids and ionic liquids-mixtures allows the systematic description of the electrolyte evolution. In particular, dynamic bulk-transport effects and interfacial structures can be studied. The theory is thermodynamically consistent and describes multi-component solutions (ionic liquids, highly concentrated electrolytes, water-in-salt electrolytes). Here, we give a detailed derivation of the theory and focus on bulk transport processes of ionic liquids as appearing in electrochemical cells. In addition, we validate our framework for a zinc-ion battery based on a mixture of ionic-liquid and water as electrolyte.


INTRODUCTION
The political and social demand for ecologically friendly, low cost, and rechargeable batteries with high energy densities is an increasing research stimulus for the improvement of materials, and the development of novel batterycomponents. 1 Electrolytes play a key role for the performance of electrochemical systems. 2 Common types of electrolytes include aqueous electrolytes, 3 organic electrolytes, 4 solid electrolytes, 5 and recently, ionic liquids (ILs). 6 ILs have attracted attention in the context of many technologies, including energy management, electrodeposition, bioscience, and biomechanics. 7 A large number of ILs can be mixed from various cations and anions. 8 This allows to tailor-cut them into task specific designer electrolytes. Thus, many ILs exhibit a variety of beneficial properties which makes them promising candidates for novel electrolytes. 9 Their favorable properties include large electrochemical windows, low flammability, or low vapor pressures. Thus, they can be compatible with high voltage electrodes, offer intrinsic safety, or be stable towards air. 10 Furthermore, some ILs suppress dendritic growth during metal deposition. 11 Depending on their composition, many ILs are environmentally friendly. 12 Current theoretical studies concentrate on the bulk, 13 or on interfacial structures. 14 Unifying approaches describing both scales are rare, but deliver remarkable results. [15][16][17][18] However, a complete unified approach for the dynamic description of ILs, and multi-component IL-mixtures, at bulk and interface, is still missing in the literature. In this article, we present such a unified framework based on the concepts of rational thermodynamics (RT) which we previously applied to electrolytes with neutral solvents. 19,20 RT provides a thermodynamically consistent framework to model a great variety of non-equilibrium systems. Here, we use RT to derive a continuum transport theory for strongly correlated electrolytes. In this work, we focus on bulk-transport. In a previous publication, 21 we showed how to supplement the theory with hardcore interactions and describe the interfacial behavior of ILs. We illustrate our holistic framework in fig. 1.
Molecular/atomistic studies of ILs are based on molecular dynamics (MD) simulations and (classical) density functional theory (DFT) simulations. DFT resolves microscopic ion-properties, delivers detailed insights into the molecular arrangement in the electrochemical double layer (EDL) 22 and describes bulk properties like ion-pair formation 23 or smallscale ion-diffusion. 24 MD simulations resolve the complete molecular arrangement and describe the evolution of the nano-structured bulk-landscape of ILs, dependent on external agents like temperature, electric fields and pressure. [25][26][27][28][29] However, DFT/MD simulations are limited due to their computational costs; simulations at length-scales above the nanometer scale are hardly accessible by these atomistic methods.
Continuum theories provide a complementary methodology for dynamic transport simulations of larger systems. Recently, Bazant et al. proposed a phenomenological meanfield-theory for binary ILs at electrified interfaces. 30 Using a generalized Landau-Ginzburg functional, the authors show that higher gradients in the electric potential lead to quasicrystalline structures near electrified interfaces. This finding is confirmed by our theory (see comments below eq. (65)). 21 Yochelis et al. rationalized this approach, and extended it to bulk phenomena and ternary ILs. [15][16][17][18] Their theoretical work allows fast dynamic simulations but cannot resolve detailed transport processes.
Dilute and concentrated electrolytes for lithium-ion batteries drive the development of continuum transport theories. 31 Recently, significant effort was put in the rationalization of consistent theories for neutral-solvent-based electrolytes with high amount of salt. 19,20,32,33 Although neat ILs constitute the extreme limit of concentrated electrolytes, where the neutral solvent vanishes, such transport theories for lithium ion batteries cannot be generalized to solvent-free electrolytes. This is because ILs exhibit some exceptional behaviour. For example, ILs form characteristic quasi-crystalline structures near electrified interfaces. These extend over a couple of nanometers, which is not observed in regular concentrated electrolytes. 34 Our transport theory for ILs, however, also describes stan-dard electrolytes, 35 as the IL-effects decrease under water-dilution, 36 or minor additive salts. 21 Transport theories for highly concentrated electrolytes should contain a condition which prevents "Coulomb collapse" of the ions. 37 The prevalent electrostatic attraction can be counteracted by volumetric constraints ("mean steric effects") or repulsive particle interactions ("atomistic volumetric exclusion"). Here, we impose a mean-volume constraint for multicomponent-incompressibility. This results in a threshold for local ion-concentration due to finite volumes. Near electrified interfaces, the concentration-threshold leads to crowding effects. 21 (Imposing repulsive interactions implies microscopic effects of excluded volume, which lead to overscreening. 21 ) The momentum equation is explicitly used in our derivation. Thereby, we capture the coupling of electric and mechanical stresses in the chemical potentials. This has significant consequences in confined geometries, e.g. electrochemical double layers, where volumetric constraints are in strong competition with Coulomb forces. 21 Electrolyte momentum is given by the center-of-mass convection velocity. Thus, the evolution of convection determines the momentum equation, whereas its variation, expressed by the rate-of-strain, is mandatory for the correct electro-mechanical couplings in the stress tensor. Furthermore, convection is a significant transport-mechanism in solvent-free electrolytes. We take account for these phenomena, derive a convection-equation and consider dissipative, viscous stress in the momentum equation.
We use ILs as validation-template, exemplifying extremely correlated electrolytes. Nevertheless, our theory applies to any liquid multi-component electrolyte, composed of arbitrarily charged (also, neutral) species. Thus it exhibits a very general structure. Our scope of electrolytes includes ILs, IL-salt mixtures, concentrated electrolytes (including aqueous electrolytes) and novel "water-in-salt" electrolytes. 38 In particular, for many applications, neat ILs or solutions with high saltconcentrations are supplemented by additives (e.g., water, 39 organic solvents, 40 or salts 41 ) to improve their performance as electrolytes. The transport theory presented in this work can be tailor-cut to any such heterogeneous mixture.
This article is structured into two main parts. In the first part, we derive our transport theory, which we validate for a zinc-ion battery in the second part. Finally, we discuss the novel aspects of our transport theory in relation to previous works.

TRANSPORT THEORY
Our transport theory relies on non-equilibrium thermodynamics, supplemented with elements of electromagnetic theory, and mechanics. 32,[42][43][44] We split this theory-chapter into two sections. In the first section General Transport Theory, we derive the universal framework, which applies to a wide class of systems. We highlight the logical steps of this derivation in four subsections. First, we derive the entropy production rate, eq. (15), from the fundamental assumption of conservation of mass, charge, and energy, by formulating the coupled balance equations for momentum, energy and entropy. We use symmetry arguments to guarantee that our model is independent from the state of the observer. 45,46 Thus, the set of variables is restricted to so-called objective quantities (see section SI 1 A). 47 Our method applies to a large class of materials and leaves a broad tunability for distinct physical systems. Second, we perform our first modeling choices and select a variable-set for strongly correlated electrolytes, eq. (17). Then, we evaluate the entropy production for a general free energy. In the next two steps, we use a linear Onsager-Ansatz to close the system of differential equations. In the third subsection, we couple the thermodynamic fluxes and forces, eq. (31), which ensures a positive entropy production. Fourth, we assume a linear relationship between viscosity and the velocity-gradients, eq. (49). This determines the model-dependent stress-tensor, eq. (50).
In the second section of this theory-chapter, we specify our framework and model the free energy density for ILs, eq. (52). This casts the specific properties of the medium into the formalism. Here, these properties comprise liquid state, polarizability, temperature dependence, viscosity, and multicomponent structure. We further assume incompressibility, and, in the case of bulk transport, electroneutrality. These strong correlations allow to reduce the minimal set of independent variables, eqs. (62) and (63). Finally, we state the equations of motion, eqs. (65), (67) to (70) and (73).

General Transport Theory
Second Law of Thermodynamics: Entropy As introduction to our derivation, we present the general structure of a transport theory based on the concepts of non-equilibrium thermodynamics. Our model is implemented as continuum theory for the fundamental physical quantities mass, charge, momentum, energy, and entropy, described by local volume-specific field densities ψ A (x) at position x = (x 1 , x 2 , x 3 ). We refer to physical quantities with capital letters A, B, C, . . ., to spatial dimensions with lowercase letters i, j, k, . . ., and to dissolved species with Greek letters α, β, γ, . . .. The time evolution of ψ A (x) is governed by continuity equations, 48 The partial time derivative ∂ t ψ A refers to temporal changes in the laboratory system, and the total time derivativeψ A describes temporal changes relative to a co-moving observer at velocity v. Both are related by a convection term,ψ A = ∂ t ψ A + (v · ∇)ψ A . 42 The convection velocity v is defined in the fixed laboratory-frame as center-of-mass average, ρv = N α=1 ρ α v α , describing bulk electrolyte-momentum per unit mass.
We derive expressions for the production rates r A and the non-convective flux densities ξ A for each field variable ψ A in FIG. 1. Scheme of our holistic framework, which captures length-scales from battery-cells to particle interactions. Thus, it describes macroscopic phenomena, like discharging/charging of a battery, mesoscopic effects, like specific electrolyte-dynamics, but also interfacial effects occurring at microscopic scales, like crowding and overscreening. the remainder of this section.
The total mass density ρ = N α=1 ρ α is conserved in the bulk electrolyte, N α=1 M α r α = 0, and it evolves due to convection only, This special role of the total mass density is used to express the continuity equation 2 in a simpler form for the mass-specific field density ψ A = ψ A /ρ. Individual concentrations c α , however, are affected by non-convective fluxes, where N α = c α (v α − v) are the flux densities in the centerof-mass system. We note that these flux densities are related, N α=1 M α N α = 0. As (total) charge = F N α=1 z α c α is conserved, we state the balance equation where J = F N α=1 z α N α is the current density in the centerof-mass system.
Total momentum density ρg comprises kinetic and electromagnetic contributions. The correct form of the electromagnetic contribution is under debate in the literature, and various forms have been suggested. 49 Here, we set ρg = ρv + D ∧ B, using the Minkowski-momentum. 50 We use D and B as a Galilei-invariant description of the electromagnetic field, which is consistent with our choice of variables for the free energy density (see eq. (17)). We formulate momentum conservation using Euler's first law of mechanics, Thus, ρg changes due to body-forces ρb acting upon the system, e.g. gravity, and due to anisotropic surface-forces described with the Cauchy-stress tensor σ. Here, the tensorgradient is defined as (∇σ) i = 3 j=1 ∂ j σ i j . We neglect molecular orientations (internal spin) in the isotropic bulk liquid, which implies a symmetric stress tensor, σ = σ T . 51 In the absence of of gravitation, total momentum is conserved, ρb = 0. Below, we couple energy and momentum balances via bodyforces b. Thus, our derivation proceeds without an explicit model for body-forces. 52 The evolution of the energy density ε is determined by work performed on the system, The first two terms stem from the non-kinematic mechanical work performed on the system, coupling energy and momentum, see also eq. (8). The third term describes heat production h. As h couples energy and reversible entropy, no constitutive equations for h is needed. This is consistent because heat production, like body-forces b, is an external field. The last two heat-fluxes are the material heat-flux due to heat-conduction and the momentum-flux of the electromagnetic field, described by the Poynting-vector E ∧ H. Here, E = E + v ∧ B and H = H − v ∧ D are the Galilei-invariant electric and magnetic fields. 53 We eliminate the mechanical source-term ρvb by substitution of eq. (8) into eq. (9). Using the vector identity ∇(σ T v) = v · ∇σ + σ : κ, we find The Cauchy-stress tensor σ and the strain-rate-tensor κ = (gradv + grad v T )/2 couple via complete contraction σ : κ = i j σ i j κ i j . The trace of the strain-rate-tensor, trκ = ∇v, determines volume-expansion, whereas its anti-symmetric part describes shear of the liquid. 54 According to the second axiom of thermodynamics, entropy is not conserved. We express its evolution with the local Clausius-Duhem inequality, 55 where ξ s is the non-convective entropy flux density. As measure for the deviation from thermodynamic equlibrium, i.e., equality in eq. (11), we define the entropy production rate R, Thus, R/T is the irreversible part of entropy production. We want to find expressions for R and ξ s . To this aim, we replace the energy production h in the entropy production in eq. (12) with the energy evolution in eq. (10), Here, we introduce the internal energy u in the center-of-mass system, and relate it to the total energy ε with the differential relationu =ε − vġ.
Our principal modeling-quantity is the free energy density ϕ H , which is related to the internal energy u by the Legendretransformation ϕ H = u − T s, i.e., The We continue evaluating the entropy poduction R, and the entropy flux density ξ s by modeling the free energy density ϕ H in the next section. The final form of R then describes dissipative processes and determines our transport equations.

First Law Of Thermodynamics: Energy
The discussions above are valid for many materials. 44 In this section, we describe the process of modeling material properties and derive constitutive equations for the field variables.
The Helmholtz free energy F is the focal point of our material model. In this work, we study bulk-phenomena and model the corresponding local free energy density ϕ H , Nevertheless, our theory can be extended to non-local phenomena which are important at electrochemical interfaces.
For example, we have recently shown how to supplement F with an interaction-functional, incorporating the particlenature of the constituents. 21 For a multi-component, polarizable and magnetizable liquid, 49 the total differential of the free energy density reads Temperature T and concentrations c α are varied according to standard thermodynamics. 32 The differential energy E · dD + H · dB with the independent, Galilei invariant variables D and B is valid for linear and non-linear materials. 56 Previous works, however, have chosen E 20 or (P, M) 32 as independent variables. Since only unique transformations between the electromagnetic variables exist, this poses no problem, 57 but leads to different expressions for the stress-tensors and the body-foces employed in the balance-laws. 52 Varying the free energy density with respect to the strain-tensor would extend the model to (visco)-elastic materials, which lies beyond the scope of this work. This total differential in eq. (17) consists of products X A · dΥ A , where X A describes a particular reaction of the medium subject to variations of the external quantity Υ A . This motivates using Υ = {Υ A } as variable-set for F and ϕ H . However, the variational expansion of the free energy does not capture the complete variable-set, as it does not comprise viscous, dissipative effects. In particular, viscous liquids cannot sustain shear stress (unlike elastic media) and depend on the rate of strain. Due to symmetry-arguments (see section section SI 1 A), we use the symmetrized, objective strain rate tensor. 58-60 Thus, we model electrolytes with the objective variable set Υ = {T, c α , D, B, κ}.
The variation of free energy density in eq. (17) contains some constitutive equations, These are supplemented by ∂(ρϕ H )/∂κ i j =0. Equation (18) implies a constitutive equation for the internal energy density, ρu = −T 2 · ∂/∂T (ρϕ H /T ). Constitutive equations relate conjugate variable-pairs, (s, T ), (E, D), (µ α , ρ α ), (H, B) and thus complement the universal balance equations with materialspecific properties. We want to determine ρφ H as it enters the entropy production R in eq. (15). To this aim, the standard method of Coleman and Noll, 42 assumesφ H (Υ) = A ∂ϕ H /∂Υ A ·Υ A . However, this expansion would result in ambiguous constitutive equations for the chemical potentials because the variableset Υ is not independent and contains redundancies (see section SI 1 D). Instead, we determine ρφ H from the total differ-ential of the free energy density (see eq. (17)), where we use the mass continuity equation (see eq. (3)). Replacing the concentration changeċ α with the mass balance eq. (5) and using ∇v = Id : κ, we identify the entropy flux density, 44 and the entropy production rate, with the symmetric viscosity tensor To derive eq. (24), we use that E ⊗ D and H ⊗ B is symmetric, see section SI 1 C.
The entropy production rate R comprises correlations between thermodynamic forces and thermodynamic potentials, N α ∇µ α , J ∇Φ and ξ s ∇T . Positivity of R constitutes an important role in our framework, constraining the generalized fluxes. We proceed by determining the fluxes N α , ξ s , and the stress tensor σ in the next two paragraphs, see eqs. (43), (44), (46) and (50).

Flux Densities
In this section, we determine the flux densities N α , ξ s relative to the center-of-mass velocity using the Onsagerformalism. 43 This relates the transport equations to the modeled free energy density ϕ H .
A dimensional analysis of the Maxwell-equations reveals that magnetic effects in the bulk-electrolyte can be neglected under normal operating conditions of electrochemical systems. 32 We thus assume from now on the electrostatic limit (B = 0), and introduce the electric potential Φ such that Using center-of-mass convection as reference velocity in our transport theory implies that the trivial flux-constraint emerges naturally. Thus, only N-1 such fluxes N α are independent, and eq. (25) requires to reduce the number of species. To this aim, we define reduced valences and chemical potentials with respect to the designated species α = 1, In the validation chapter, we demonstrate that the predictions of our theory, if applied correctly, are independent of the choice of the designated species. This reduction can lead to counter-intuitive contributions to the electric current, By construction,z 1 = 0. Thus the designated species does not contribute, even if it is charged, z 1 0. Furthermore, if the designated species is charged, z 1 0, then even neutral species carry effective charge and contribute to J . This is one reason behind recent confusion with respect to solventfree electrolytes, e.g., molten-salts, 61,62 and ILs. 63,64 For electrolytes with neutral solvents, the natural choice for the designated species is the solvent, 20 since thenz α = z α .
The entropy production rate in the reduced formalism becomes Here, we defined the vector of thermodynamic forces X (α) = (Fz α ∇Φ + ∇μ α , ∇T ) and the vector of thermodynamic fluxes The subscripts in brackets denote that only the ionic contributions are speciesrelated.
We make use of this expansion and transfer the products to quadratic terms via a bilinear Onsager-matrix. Thus, this Onsager formalism (i) closes the system of equations, (ii) establishes positivity of R, and (iii) evaluates the flux-force couplings. The phenomenological relations linearly couple thermodynamic fluxes and forces with the Onsager matrix L (α)(β) , The Onsager matrix L (α)(β) is symmetric and positive semidefinite. 65 This guarantees that the corresponding term in the entropy production rate R is always positive (see eq. (30)). Thus, electrolyte transport is a dissipative process, which wants to equilibrate the system. In dilute solutions, the Onsager matrix is diagonal and describes intra-species correlations. In concentrated electrolytes, off-diagonal Onsager coefficients comprise inter-species correlations, e.g., electroosmotic drag through a membrane 66 or asymmetric transference numbers. 67 In total, we define N(N + 1)/2 independent Onsager coefficients, i.e., transport parameters. The thermodynamic fluxes are driven by electric potential gradients, concentration gradients, and temperature gradients, i.e., migration, diffusion, and thermo-electric effects. In the following, we relate the abstract Onsager coefficients to these physico-chemical effects. The electric current density J (see eq. (28)) is expressed as The flux densities N α and the entropy flux density ξ s are often formulated in terms of J instead of Φ. Therefore, we transform eq. (32) by substituting ∇Φ for J with eq. (33) and introducing additional, physically motivated parameters, Transport coefficients introduced in eqs. (33), (34) and (35) are electric conductivity κ, Seebeck coefficient β, and thermal conductivity γ, 31 By construction, κ and γ are positive because the Onsager matrix is positive semi-definite. N-1 transference numbers t α relate particle fluxes and electric current in eq. (34), These N-1 parameters are subject to the normalization condition, N α=2 t α = 1, such that only N-2 transference numbers are independent. Thus, for binary electrolytes all transference numbers are fixed, t N=2 2 = 1, and N 1 relates to J by the specific mass-ratio alone, see section SI 1 H. Apparently, all t α are positive if the Onsager matrix is diagonal (as for dilute solutions).
We proceed by defining the N(N+1)/2 coefficients of the symmetric diffusion matrix D, Since D T T is determined by γ, β and κ, eqs. (36) to (40) yield N(N+3)/2 transport coefficients. Thus, the number of transport parameters defined so far exceeds the number of independent Onsager coefficients N(N + 1)/2. However, further N constraints follow from the relation of J and N α in eq. (28), Thus, only N(N-1)/2 diffusion coefficients are independent, Altogether, N(N-1)/2 independent diffusion coefficients, N-2 independent transference numbers, the electric conductivity, and the Seebeck coefficient constitute the complete set of physically motivated free parameters.
We designate a second species, α = 2, and define the set of N-2 reduced chemical potentials Thus N-1 independent fluxes (N 3 , . . . , N N , ξ s ) exist, N 2 is determined by the electric flux with eq. (28) and N 1 is determined by mass conservation in eq. (25). Our electrolyte potential Φ is the Maxwell potential that appears in the Poisson equation. Typically, 20,31 concentrated solution theory is expressed using the alternative potential If the first designated species is the neutral solvent, ϕ = Φ + µ 2 /Fz 2 corresponds to the electro-chemical potential of the second designated species. With the electrolyte potential ϕ, we can express the current density analogous to eq. (43) and eq. (44), The entropy production rate R in the double reduced formalism is where D red is the diffusion matrix D without the two designated species (see eq. (40)). To summarize, the entropy production comprises three contributions. The first term describes mechanical dissipation. In the next paragraph, we determine τ in such a way that τ : κ ≥ 0 is guaranteed. The second term describes Joule-heating due to migration. Since the electric conductivity is non-negative, J 2 /κ ≥ 0. The last term describes entropy production due to diffusion and heat conduction. Thermodynamic consistency reduces to nonnegativity of the diffusion matrix, D red ≥ 0.

Stress Tensor
In contrast to viscous-elastic models, this description does not describe the stress tensor by a constitutive equation in the form of a derivative of the free energy density. 44,68 Instead, we determine σ implicitely via the viscosity tensor. To find τ(κ), we apply the Onsager formalism as in the previous subsection, and assume that τ is linear in the strain-rate tensor. Symmetry demands that τ is an objective tensor, thus the representation theorems for isotropic tensors determine it uniquely up to two scalar transport parameters of viscosity, 69,70 where κ tf is the symmetric trace-free part of κ, λ is the bulkviscosity and η is the shear-viscosity. Since the entropy production rate R is non-negative (see eq. (24)), if the viscosities obey η ≥ 0 and λ ≥ 0. To summarize, we find as constitutive equation for the symmetric stress tensor, 71 The pressure p measured in experiments is derived from the total stress tensor, 32,72 Here, p td = N α=1 c α µ α − ρϕ H + 2ED/3 is the thermodynamic pressure comprising the standard Gibbs-Duhem contribution, and a Maxwell-contribution. The further contribution, λ∇v, describes dissipation due to friction. Since λ ≥ 0, friction can lead to positive and negative pressures, depending on ∇v.
The apparent electromagnetic contribution to the total stress σ in eq. (50) differs from the Maxwell contribution Σ = −ED/2 + E ⊗ D. Note that we recover the correct standard form once we model and specify ϕ H (see eq. (54)).

Model for Correlated Liquid Electrolytes
The transport equations of any electrolyte-model depend on the specific form of the free energy density ϕ H . In our approach, ϕ H is the focal point of modeling. Once ϕ H is specified, all equations follow from plain mathematics. In this section we introduce our model ϕ H , derive the convection velocity from the volumetric electrolyte equation of state, and discuss the dynamic transport equations.

Free Energy Density
As discussed in the previous section, we assume that the free energy density ϕ H depends on temperature T , concentrations c α , and dielectric displacement D (see eq. (17) and the discussion thereafter). 20,32 Our model for liquid electrolytes is generated by The first term comprises the electrostatic energy-density of polarizable media. 56 The eletric field follows from eq. (19), E = D/ε 0 (1 + χ), which describes a linear dielectric medium. We neglect the dependence of susceptibility χ on ion composition such that the chemical potentials do not depend on dielectric displacement (see section section SI 1 I 1).
The second term expresses volumetric contributions with partial molar volumes ν 0 α in a stable reference state, which couple to surface forces acting upon the system. These surface forces are usually expressed by pressure (see eq. (54)). Here, K is the bulk-modulus which acts as a Lagrange-mutliplier in the case of incompressible electrolytes (see eq. (59)). We expand the energy of deformations around a stable reference state in order to easily transfer to incompresssible media, see eq. (59). These volumetric energy penalties encode volume conservation (see eq. (57)), and account for mean steric effects, 14 which are crucial for continuum theories of concentrated electrolytes 73 and ILs. 16 Mean steric effects create a threshold for local ion-concentration and prevent Coulombic collapse of the system. 37 This leads to crowding of ILs near electrified interfaces. 21 We motivate this expression for the elastic energy in the supplementary (section SI 1 F).
The third term is the entropy of mixture for non-interacting systems. We model this term in analogy with ideal gases and neglect contributions from inter-molecular interactions, e.g., solvation effects.
We phenomenologically account for non-ideal interaction contributions in ρϕ int H . The activity-coefficients f α measure the deviation from ideal/dilute electrolytes (f α c=1) via ∂(ρϕ int H )/∂c α = RT ln(f α c). The interaction term ρϕ int H can also describe the heat capacity of the system, C = −T · ∂ 2 (ρϕ int H )/∂T 2 . Since thermal aspects are not our main focus, we refer to the supplementary for a thermal contribution to our free energy model ρϕ H (section SI 1 I 3).
In the following, we calulate the chemical potentials µ α and the stress tensor from the model free energy density in eq. (52). The chemical potentials follow by evaluation of the constitutive eq. (21), The first term comprises entropy of mixture and inter-molecular interaction energies. The stress tensor is determined by eq. (50), Thus, the electrostatic contribution is the expected electrostatic Maxwell-stress tensor Σ. 56 The elastic stress due to compressibility and the stress due to intra-molecular interaction is also diagonal, i.e., isotropic. Shear stress due to viscosity contains the non-isotropic contribution 2ηκ tf .

Volume Constraint and Convection
In this section, we derive the equation of state for liquids from our model free energy (see eq. (52)), and discuss incompressibility of electrolytes. In addition, we derive a multicomponent incompressibility constraint, and a transport equation for the center-of-mass convection velocity v.
In a homogeneous system we neglect viscosity. In this case p = p td becomes the thermodynamic pressure as discussed for (51). Equation (54) determines, as a function of c α , T , and D. This expression implicitly determines volume V(N α , p, T, D) as function of particle numbers N α , pressure p, temperature T , and dielectric displacement D. We make use of this description, and define partial molar volumes ν α as derivative of V with respect to N α . Implicit differentiation of the function f (c α , T, D) then determines the partial molar volumes (see section section SI 1 E), Thus, the partial molar volumes of the electrolyte follow from the stress tensor. As V(N α , p, T, D) is a homogeneous function of first order in N β , Euler's homogeneous function theorem yields a compact form of the electrolyte equation of state of liquid electrolytes, 16,[73][74][75][76] which takes the form of a volume constraint. From eq. (56), we calculate the partial molar volumes for our model. Liquid electrolytes are hardly compressible, i.e., K p. Repeated application of this approximation yields Here, we used eq. (55) and eq. (54) to replace N β=1 ν 0 β c β . Thus, K is indeed the bulk-modulus, and ν 0 α is the partial molar volume at standard pressure p = 0. For incompressible electrolytes K → ∞, the partial molar volumes ν α = ν 0 α do not depend on pressure, and eq. (57) constitutes an incompressibility constraint. 21 In this case, pressure cannot be determined by eq. (55). Nevertheless, the volumetric contributions in the chemical potentials have to be determined (see eq. (53)). We show in eq. (65), how to solve this challenge with the volumetric contributions in the stress-tensor and momentum conservation.
The volume constraint in eq. (57) remains fulfilled under multi-component transport as the center-of-mass motion balances the volumes. Thus, we can derive an equation for the velocity v from the volume constraint. To this aim, let us first proof a Lemma. Partial molar volumes ν α defined in eq. (56) obey the symmetry-property ∂ν α /∂N β = ∂ν β /∂N α . Applying Euler's homogeneous function theorem to ν α ensures that N β=1 c β ∂ν α /∂N β = 0 holds for all species. These relations proof the following Lemma, where we used the independence of the primary variables c α , T and the volume constraint. Applying the total time derivative to the volume constraint in eq. (57), we conclude from our Lemma that N α=1ċ α ν α = 0. In combination with mass conservation eq. (5) for the concentrations c α , we finally derive the equation for the center-ofmass velocity v, 35,77,78 This equation is generally true for compressible as well as incompressible electrolytes and for constant as well as nonconstant partial molar volumes. The left-hand side of eq. (61) expresses local, isotropic volume-expansion, ∇v = tr(κ). Thus, the first term on the right measures volume-expansion caused by transport. In absence of the second term, and for constant partial molar volumes, the total flux of volume-fractions c α ν α v α measured in the fixed frame is conserved, ∇ N α=1 c α ν α v α = 0. The second term in eq. (61) measures volume-production by chemical reactions, and is an important source for convection in multicomponent systems.
In the case of a single-component liquid, ∇v = 0 according to eq. (61), since by construction the only flux density relative to the center of mass motion vanishes, N 1 = 0. This is a standard-condition for incompressibility of non-reactive media. However, it is often used for complex electrolytemixtures which can be a bad approximation.
From now on we consider electrolytes in the incompressible limit K → ∞, such that ν 0 α = ν α becomes the general partial molar volume (defined for any state).
As we showed above, conservation of mass and charge reduce the number of independent species by two. We assign the designated two physical species to α = 1, α = 2. The volume constraint eq. (57) allows to determine the corresponding concentrations, with the charge density . Furthermore, we transform eq. (61) into the reduced description, whereν α = ν α − M α /M 1 · ν 1 andν α =ν α −z α /z 2 ·ν 2 define the set of independent partial molar volumes. In particular, eq. (63) implies constant convection-profiles for electroneutral, binary electrolytes. Conversely, this also shows that pure ILs cannot sustain concentration gradients. 61

Equations of Motion
The force law follows from momentum-balance in eq. (8) and takes the form of a generalized Navier-Stokes equation. However, we assume highly effective momentum-dissipation in such viscous media. We neglect external body forces, 79 , and assume mechanical equilibrium, 21 Above we assumed constant viscosity-parameters, constant bulk-modulus, and isothermal conditions. Although the force law is not used in most of the battery literature, 20 it plays a fundamental role for highly concentrated electrolytes, 79 and ionic liquids, 21 where Coulomb interactions compete with saturation effects. Stress couples to transport via the volumetric term in the chemical potentials (see eq. (53)). Via this mechanism, we ensure that the complete, coupled set of mechanical and electrostatic stresses is comprised in our theory. By substitution, dissipative electric, viscous, and interaction forces contribute to the chemical potentials, where is the thermodynamic factor. It is determined by the constitutive equation ∇ ln c α /c + N β=1 (δ αβ − ν α c β )∇∂(ρϕ int H )/∂c β , and differs from the standard form, 80 as it accounts for mean steric effects due to nonvanishing molar volumes. Note that the bulk modulus K, which is not a viable material parameter in the incompressible limit, disappears in eq. (65).
Equation (65) constitutes the complete, mechanically coupled form for the chemical potentials in our electrolyte model. These thermodynamic forces obey a viscous Gibbs-Duhem relation, N α=1 c α (∇µ α +Fz a ∇Φ)=(λ+η)∇(∇v)+η∇ 2 v). In equilibrium we can write this with the thermodynamic pressure defined in eq. (51) as ∇p td = − ∇Φ − ∇(ED)/2. As consequence, strong pressure gradients emerge in electrochemical double-layers. 21 We use Poisson's equation for the coupling of electric potential and charge density in the electrolyte, and use the dielectric constant ε R = 1 + χ. This closes the set of complete isothermal equations, Equation (69) comprise N-2 independent equations for the independent concentrations c 3 , . . . , c N (where c 1 and c 2 follow from eq. (62)). Here, we restate the ionic fluxes and the conduction current density, which are subject to the form of the chemical potentials following from the modeled free energy density eq. (65). Note that ϕ = Φ +μ 2 /Fz 2 is the reduced electrochemical potential introduced in eq. (45). Volume-expansion ∇v is determined by eq. (63). Alternatively, the complete set of transport equations can be cast into matrix-form (see section SI 1 J). In section SI 1 I, we complete the set of transport equations by deriving the equation for temperature ("heat equation"),

SIMULATION OF ZINC ION BATTERY CELL
Secondary zinc-ion batteries (ZIB) are an emerging technology for safe and low-cost energy storage. [81][82][83] Here, we model a secondary ZIB with IL-water mixture as electrolyte, based on the experimental work described in 41. Comparison with these experimental results serves as validation for our transport theory. First, we sketch the cell set-up and describe our modeling equations. Second, we present our simulation results.

FIG. 2. Scheme of the simulated zinc-ion battery.
The electrodes consist of a porous zinc-anode (zinc powder) and a Prussian-blue-analogue (PBA) cathode (FeFe(CN) 6 ). The employed electrolyte is composed of a mixture of choline acetate ([Ch]OAc) with 30 wt % water into which 1M of zinc acetate (Zn(OAc) 2 ) is dissolved. Despite the significant amount of water, the [Ch]OAc+water mixture can (still) be denoted "IL with water", and may thus be viewed as highly concentrated. 84,85 This terminology is motivated by the large mass fraction of the salt, (ρ 0 3 + ρ 0 4 )/ρ, see table I. The ZIB is illustrated in fig. 2.

Model Equations
We neglect hydration of the salt-ions by water and assume a complete dissociation of the salt [Ch]OAc − −− −− − Ch + + OAc − . In principle, various different ionic zinccomplexes may form under bulk-reactions, depending on the pH-value of the solution. 86,87 However, such a detailed investigation lies beyond the scope of this paper. Thus, we assume complete complexation of the zinc-salt according to 88 We model the electrolyte as quaternary mixture of the three charged constituents Ch + , OAc − , [Zn(OAc) 3 ] − , and water. [89][90][91] As reference, we designate water, thusz α = z α . This choice corresponds to transport theories developed for aqueous solutions, based on water as solvent. 19 Since length scales above some nanometers suffice for the description of batteries, 19 we can safely assume our system to be electroneutral, viz.
= 0. 31 As a consequence, the charge density is not a free variable and we solve ∂ t = 0 in eq. (68) for ∇J instead of Poisson's eq. (67). As thermal aspects play a minor role in elementary battery cells, 20 we assume thermal equilibrium at room temperature, i.e., ∇T = 0. Hence, the temperature appears only as constant a parameter. Here, we set the activity coefficients to f α = 1 and neglect non-ideal interactions ϕ int H = 0. Furthermore, we assume constant viscosity-coefficients (η = 25.3 mPa s and λ = 0) 92 in the chemical potentials.
The only independent specific concentrations are c 3 and c 4 , which determine c 1 and c 2 via eq. (62), subject to the condition = 0.
Due to the porous cell structure, the set of transport equations must be modified by means of volume-averaging, as described by porous electrode theory. 93 Volume-averaging over liquid and solid phases implies volumetric (ε) and dynamic (ε β ) transport-corrections, where ε = V l /V is the liquid-phase volume-fraction or porosity. β is the Bruggemann-coefficient, which depends on the microstructure, Reactions occur only at the electrode surfaces and appear as Since water is the designated species,z α = z α . The columns show initial concentrations c 0 α , mass fractions ρ 0 α / ρ, and volume fractions c 0 α ν α . source-terms in our model equations, defined as This sum includes all reactions k at all electrode-interfaces Γ, involving species α. Here, ν Γ k;α are the stoichiometries of the reactions, and the specific surfaces a Γ measure the surface-tovolume ratio of the electrode Γ. The specific surface-reactionrate i Γ constitutes interface-conditions between the solid and liquid phase, and depends upon the reaction overpotential. We model this quantity using a Butler-Vollmer-Ansatz, 94 see section SI 1 L.
The half cell ractions occuring at the Zn-anode and the PBA-cathode (FeFe(CN) 6 ) are (see section SI 2 A) [Zn(OAc) 3 ] − + 2 e − + FeFe(CN) 6 Thus, zinc dissolves from the anode and intercalates into the cathode (and vice versa). We neglect solid-state diffusion and model the cathodic reactions as deposition-processes (see section SI 2 B 1).

Simulation Results
In this section, we present the simulation results of the asmodeled ZIB and compare them with the experimental results (we specify our numerical methods in section SI 2 C). 41 In this way, we study the spatial profile and temporal evolution of concentrations and convection velocity. We highlight the relevance of convection for cell performance and compare possible definitions of transference numbers.

Competition of Diffusion, Migration, and Convection
First, we validate our model. For this purpose, we simulate the galvanostatic discharge and charge of the ZIB by applying a moderate current density (I = 0.1 mA cm −2 ). We compare our simulation results with experimental observations in fig. 3. Apparently, the specific capacities and cell voltages obtained from simulation are in good agreement with the experimental values. However, the discharge profile obtained from experiment exhibits two discharge phases with a transition at approximately 20 mA h g −1 . Endres et al. 41 attribute this effect to two different electro-reactivities, stemming from low-/ and high-spin states of Fe(III) in the PBA. Because such atomistic processes lie beyond the scope of our model in this paper, this transition is not found in simulations.
Our model provides a complete description of the electrolyte dynamics. In the following, we discuss the ion dynamics and their influence on the overall cell performance in order to understand the interplay between the different transport mechanisms (migration, diffusion, and convection), and the reactions at the electrode surfaces.
For a proper discussion of the electrolyte evolution during discharge of the battery, we designate some characteristic moments. Figure SI 1 shows the discharge-profile of the cell-voltage and the designated moments. We set one focus on processes occuring during the initial discharge-phase (t = 250 s, 420 s, 520 s). Furthermore, we designate two intermediate moments, and, finally, the moment of complete discharge of the cell. The significance of the designated times becomes apparent below when evaluating the electrolyte-quantities Φ, v, and c α at these moments. Figure 4 shows the temporal evolution of the convection velocity profile. Apparently, the direction of the motion of the center-of-mass switches from towards the anode (negative convection), to towards the cathode (positive convection) right after the initial phase, and, then, the profile quickly becomes quasi-stationary. As discussed below, v is not dominant in our system. In the supplementary, we observe a similar behavior for the ionic concentrations (see fig. SI 1c)-f)). Initially (first three designated moments), concentration-differentials grow throughout the cell, followed by quasi-stationary profiles during the rest of the discharge. This suggests that the dynamics of the concentration and the convection profile arise from a common initiation-mechanism. Due to the interplay of multiple transport mechanisms, i.e., diffusion, migration, and convection, the dynamics of the concentrations and of the convection velocity are coupled.
We explain these observations with the push of the electrolyte out of its initial equilibrium-state by the application of the discharge current through the reactions at the electrodes. As the system relaxes towards the new stationary state under galvanostatic discharge, the concentrations approach stationary profiles. Spontaneous electrode reactions directly influence the profile of the ion concentrations near the electrodes. In addition, concentration gradients are built up, due to convective transport of zinc from anode to cathode. As zinc is dissolved as complex [Zn(OAc) 3 ] − , a net loss of available electrolyte volume occurs at the cathode, and a net gain of electrolyte volume at the anode. This pushes electrolyte from cathode to anode according to eq. (76) and results in negative convection velocities. We give a detailed analysis of the electrolyte-dynamics in the SI (see section SI 3 A).
The evolution of the electric potential in the electrolyte is shown in Figure SI 1b). As expected, the cathode is more electronegative than the anode (∆Φ cell ≈ −0.1 mV) at all times during discharge. This implies a migrational pull of the anionic species towards the anode, and of the Ch + -ions towards the cathode. Thus, migration hinders cell operation, which depends on negative [Zn(OAc) 3 ] --ions. As a consequence, diffusion and convection must overcompensate the electric forces in the electrolyte to sustain cell operation. Figure 5 shows normalized, quasi-stationary concentrationprofiles of all electrolyte species at end-of-discharge. Inter- estingly, the gradients of the two anionic species (OAcand [Zn(OAc) 3 ] -) have opposite orientations. The profile for the positively charged Ch + -ions moderately increases at the more electronegative cathode. Likewise, water shows a small gradient, despite being neutral (z H 2 O = 0 in this reference-frame) and thus not being susceptible to migration. Most importantly, the large concentration gradient of the [Zn(OAc) 3 ] --ions is necessary to overcome the migration pull. This is mandatory to maintain cell operation, as zinc is intercalated into the cathodic PBA-structure. OAc --ions have to move back to the anode, which is realized by a small concentration gradient and the migration pull. This motion of bulky OAc --ions induces the positive convection velocity during the stationary phase of galvanostatic discharge seen in fig. 4 (see eq. (76)). Figure 6 illustrates the relevance of convective transport for the electrolyte-species. The ratio ε β N α /εc α v shows the relation between flux densities within the center-of-mass system and the center-of-mass velocities. This ratio is larger than one (roughly one), if convection is negligible (dominant) in transport of the respective species. Apparently, convective flux contributions are negligible for the two negative species, whereas convection plays a significant role for the dynamics of water and Ch + -ions. Thus, convection is important for those species that do not contribute to the half-cell reactions. We infer from the sign of the ratio ε β N α /εc α v that OAc --ions move towards the anode, whereas [Zn(OAc) 3 ] --ions moves towards the cathode. This confirms our previous finding for the [Zn(OAc) 3 ] --ions that diffusion overcompensates migration.

Limiting Discharge Currents
Now, we investigate the power limiting mechanisms of this ZIB with IL electrolyte. For this purpose, we simulate the cell discharge at high current densities. Figure 7 illustrates the impact of increased discharge cur-rents on the overall cell performance. The shape of the discharge profile is preserved for moderate discharge current densities (up to I = 2 mA cm −2 ), despite being shifted to smaller discharged capacities. At higher current densities, a strong capacity fade is observed and the discharge profiles exhibit steep voltage drops. These voltage drops suggest that diffusion becomes too slow to supply the cathodic interfacial reaction mechanism with sufficient amount of [Zn(OAc) 3 ] --salt. Such a behavior is typical for ILs, which are often highly viscous. 95 This explanation is confirmed by fig. 8, in which the concentration profiles of the [Zn(OAc) 3 ] --ions for the different discharge currents at the end of discharge is shown. We observe steep concentration gradients throughout the cell. The salt concentration at the cathode decreases with increasing discharge currents, and, finally, [Zn(OAc) 3 ] --ions get depleted.  3 ]concentration at the anode increases with discharge current density. Interestingly, it starts to decrease for very high current densities. In this case, the cell fails so quickly due to diffusion limitation that a quasi-stationary concentration-profile cannot establish.

Discussion of Transference Numbers
Finally, we demonstrate the consistency of our transport theory. In our model, we choose two reference species because of momentum and charge conservation. The model predictions should be independent of this choice. Here, we compare simulations with different reference species.
In the theory-chapter, we showed that all N(N+1)/2 transport-parameters, appearing in a N-component electrolytemixture, follow from the Onsager matrix L, see eqs. (36) to (40). Although the designated species does not appear explicitly, L comprises the complete set of inter-speciescorrelations, including correlations involving the designated species. This becomes apparent in the closure-relation for the independent fluxes N α | α≥2 , eq. (31), which implicitly determines N 1 . Thus, the Onsager-matrices with respect to different designated species are mutually coupled, and cannot be stated independently from each other. In section SI 1 K we show that simple conversion relations exist, which allow to transfer between the different choices of reference species.
In section SI 3 B we prove consistency of our framework. Figure SI 3 shows a comparison of the discharge curves and the final concentration profiles of the [Zn(OAc) 3 ] --ions, obtained by simulations using the three different sets of reference species comprised in table II. Apparently, the deviations lie within numerical accuracy (relative error ≈ O(10 −6 )). This is consistent, since both cell voltage and concentrations must not depend upon the choice of reference. Thus, they give, up to numerical accuracy, identical results.
In contrast to the electric conductivity, some transport parameters, e.g., transference numbers, depend on the choice of reference species. Figure 9 shows the transference numbers at the end of discharge for two different reference species. In addition, table II summarizes spatially averaged values of the transference numbers at end of discharge, for three different sets of reference species. In quaternary electrolytes, only two transference numbers are independent and only three transference numbers are well-defined.
We observe in fig. 9 that sign and magnitude of the transference numbers depends on the reference species chosen due to momentum conservation. Interestingly, t OAc − is rather similar for both reference species H 2 O, and Ch + , respectively. A significant discrepancy is observed for t [Zn(OAc) 3 ] − . When Ch + as charged species is designated, water acquires an effective charge, and thus contributes to the electric current. In particular, t H 2 O is negative in the Ch + -frame.
Transference numbers are particularly intuitive for the neutral reference species H 2 O. In this setting, the signs of transference numbers endorse our interpretation of the overall electrolyte dynamics, discussed above.   (43)).
The relations between the numerical results for the transference numbers obtained in the three different referenceframes agree with the analytical predictions, see eqs. (SI 122a) to (SI 122f).
Since, there is an ongoing debate regarding the sign and magnitude of transference numbers in concentrated electrolytes and ionic liquids, 63,64 we shall briefly comment on this topic here.
Various, differing definitions for transport parameters exist in the literature, which bears the potential for confusion when comparing results from different authors. Thus, a complete characterization should state all transference numbers and make their definitions clear. In our theory, we define transference numbers and diffusion coefficients relative to the bulk-motion of the center-of-mass. As it was shown earlier for molten salts, 62 such an internal description predicts that transference numbers can take up any real value if the Onsager matrix contains off-diagonal terms. Some works in the literature perform a less general approach and choose one specific ionic species as internal reference, 96 or choose an external description using fixed coordinates as reference. 97 Since each approach bears its own experimental significance, 98 we resolve this ambiguity in section SI 1 G. There, we derive simple relations, which transform between the different referencesystems.
Of course, the theoretical concepts must be probed by experiments. However, the experimental determination of transport parameters in concentrated electrolytes is challenging. In principle, suitable methods must take into account all interspecies-correlations and the formation of ionic clusters and ion-aggregates. 99 In particular, applicability of NMR/PFG-NMR experiments for determination of transport parameters is limited for concentrated electrolytes, as they (i) neglect formation of ionic complexes, (ii) derive transference numbers from diffusion coefficients via Nernst-Einstein relations strictly valid only for ideal electrolytes and (iii) provide only averaged values. 64 Recently, the measurement of ionic electrophoretic mobilities was proposed to overcome these obstacles, either using non-blocking electrodes, 67 or experiments based on electrophoretic NMR (eNMR). 100,101 Unfortunately, the first technique does not apply to neat IL, as it is not suitable for non-metal ions. 102 In contrast, eNMR applies to both, neat ILs and Li-IL-mixtures. 63 In principle, this method is evaluated with the external fixed laboratory-frame as reference. However, this is correct only if convection is negligible (see eq. (SI 27)).

DISCUSSION
In this section, we discuss the relevance of the derived transport model. We give a detailed comparison with continuum models previously presented in the literature, before we quantitatively asses the role of center-of-mass convection and mechanical forces for a consistent model of highly concentrated electrolytes.
The so-called Newman-model is standard for the mathematical modeling of batteries on the cell-level. 31,103,104 This approach relies on Stefan-Maxwell theory to relate flux densities and chemical potentials. The canonical Newman model describes a ternary system composed of a cation, an anion, and a neutral solvent with three parameters, diffusivity, conductivity, transference number. It is referred to as concentrated solution theory, as it takes correlations between anions and cations into account. However, the dynamics of the solvent is usually neglected. 103 Monroe et al. extend the Newman model to locally non-neutral multi-component concentrated electrolytes. 105 They model flux densities relative to a designated species-velocity, which serves as convection velocity. In accordance with this description, the transport parameters relate to the dynamics of the designated species ("Hittorftransference" numbers). However, this approach does not take into account momentum conservation and mechanical couplings. We show below that these become relevant in highly concentrated electrolytes.
These Newman-type models can be generalized to account for the electrolyte equation of state, eq. (57) . 74,106 We use it to derive the convection equation in our framework, see eq. (61). Monroe et al. conclude that the relative magnitude of the molar volumes is a crucial factor (see also eqs. (80) to (82)). 106 Nevertheless, their theory is evaluated only for ternary systems with neutral solvent.
A systematic approach for the description of liquid electrolytes has recently been introduced by Latz and coworkers. 19,20 Their approach uses a thermodynamically consistent framework, 43,69 which is based on thermodynamics principles and balancing laws to derive the general structure of the transport equations with an Onsager Ansatz. 107,108 They discuss a ternary system in comparison to the stan-dard Newman model. Mutual couplings between the ionspecies are evaluated, but the dynamics of the neutral solvent is neglected in the description of transport. Furthermore, the electro-mechanical forces following from momentum conservation are not further evaluated.
In a series of publications, 32,75,109 Guhlke et al. presented a similar framework. Their theory is developed using the same underlying rigorous assumptions and accounts for multicomponent mixtures. Furthermore, similar to our approach, all mechanical couplings are incorporated into the transport equations. One highlight of this framework is the thermodynamic description of singular surfaces, which is used for the consistent description of electrochemical double layers. The role of convection for such thermodynamically consistent frameworks is discussed by the authors in great generality in Ref. 65. However, similar to the description of Latz and coworkers, the role of convection is assumed negligible in the discussion of a ternary electrolyte with neutral solvent.
We extend the approach of Latz and coworkers, and derive a more general, and more consistent framework, as we do not restrict the number of species, their valences, nor neglect transport of neutral solvent. At the same time, we make sure that the theory intuitively connects to the standard Newman model by using comparable parameter definitions. We highlight two advances of our model: First, center-of-mass convection plays a key role in our theory. Thus, molar masses appear in the definition of transport parameters, see for example eq. (26). Second, we take into account the coupling of mechanics and transport. Thus, molar volumes appear in the chemical potentials and affect the transport dynamics of the electrolyte, see for example eq. (66). Due to the consistent reduction of independent transport parameters, the abstract form of the general transport equations simplify when specified to particular systems. For example, the electric conductivity is the only transport parameter needed to describe a binary IL. In the following, we quantitatively discuss the importance of these two fundamental extensions for highly concentrated electrolytes.
Equation (63) determines the spatial inhomogeneity of the center-of-mass velocity. We calculate the ratio of the variation of convective flux density and the variation of non-convective flux density. In a ternary electrolyte with neutral solvent (z 1 = 0) and binary salt (c 2 = c 3 , z 2 = z 3 ), it holds Thus, the relative mass density of the solvent ρ 1 /ρ and the relative volume density of the solvent c 1 ν 1 determine the relative variation of convection velocity. If the solvent dominates electrolyte mass and volume, i.e., ρ 1 ∼ ρ and c 1 ν 1 ∼ 1, the convection velocity is constant. If solvent mass is negligible, i.e., ρ 1 ρ, the variation of the center-of-mass convective flux outweighs the relative variation of the salt flux. Thus, we define the term highly concentrated salts based on mass fractions, e.g., ρ 1 /ρ.
The absolute magnitude of the convection velocity is determined by the electrochemical reactions at the electrodeelectrolyte interfaces. The non-convective species-fluxes N α + c α v = c α v α are subject to flux boundary conditions, c α v α | Γ = R Γ ν Γ α , at each electrode-interface Γ. Here, we model the reaction source-terms via interfacial currents as defined by eq. (77) such that R Γ = dx r α /ν Γ α , where ν Γ α denotes the stoichiometry of the reactions. Thus, the relevance of convective fluxes at the interface is determined by the masses / mass-densities, and the stoichiometries, We conclude that the mass fractions ρ α /ρ determine the relevance of the convective flux density. We apply this analytic result to the zinc ion electrolyte specified in table I. At the anode, we find for the reacting species These values show that in this electrolyte center-of-mass convection is relevant, but not dominant. This is in agreement with the mass fraction of water, see table I, which is comparable to the mass fraction of the salt, but not negligible. Thus, in the bulk the electrolyte behaves as concentrated electrolyte. Our simulation results in fig. 6 validate our estimates for inhomogeneity and magnitude of the center-of mass velocity in eq. (80) and eq. (81), respectively. We note that the apparent complexity of the electrolyte depends on the choice of reference species as shown in fig. 9. Furthermore, electrolyte composition is significantly altered at electrified interfaces, where ions can accumulate. 21 It is the unique strength of our consistent transport theory that it describes both, the interfacial and the bulk behavior of concentrated electrolytes independent from the choice of reference species.
We predict an additional contribution to the thermodynamic factor in eqs. (65) and (66), based on the consistent coupling of mechanical forces. For ideal electrolytes (f α = 1 for all species α), it holds ∇µ mixing Thus, we find cross-couplings between species with different molar volumes. We discuss the relevance of speciesasymmetry with two limiting cases, where we assume that the concentrations lie in the same order of magnitude, i.e., c α ∼ c β . First, if all species have the same molar volume, eq. (83) reduces to the standard ideal form, where all species are decoupled, ∇µ mixing α f γ =0 = RT ∇c α /c α . Second, if the first species is much larger than the others, ν 1 ν α | α≥2 , then the forces on the small species effectively decouple, ∇µ mixing α | α≥2 = RT · ∇c α /c α , whereas the designated species experiences only inter-species correlations ∇µ mixing 1 = −RT ν 1 N β 1 ∇c β . This analyis exemplifies the relevance of the consistent coupling of transport and mechanics.

CONCLUSIONS
To the best of our knowledge, we have developed the first thermodynamically consistent transport theory applicable for pure ILs. Our theory describes the complete set of coupled transport equations for composition (concentration), temperature, charge density, electric potential, and convection. We make explicit use of the force law to include all mechanical couplings. Our detailed equations for the electrolyte dynamics describe the individual contributions of all species to transport of mass, charge, heat, and to convection.
The present work completes our two-fold validation procedure. In a first publication, the theory was validated in a joint experimental/theoretical investigation of the doublelayer-structures formed by a binary IL near an electrified interface, and of how a ternary salt influences these characteristic IL-structures. 21 Here, we apply the electroneutral transport theory to a complete secondary battery cell.
Our simulations determine the concentration dependent electrolyte dynamics during a complete cycle of dischargingcharging of a secondary zinc-ion battery based on an ILelectrolyte. The resulting discharge profile is in good agreement with the experimental results. However, the cell performance is limited by the negativity of the zinc ions. These must overcome a Coulombic potential barrier to maintain cell operation by the formation of concentration gradients. In con-trast to binary ILs, convection does not play a significant role compared to the competing mechanisms of migration and diffusion.
Our theory provides a rigorous framework for the determination of the transport parameters, which all follow from the fundamental Onsager matrix. Furthermore, we enrich the present discussion regarding the interpretation of transport parameters in solvent-free electrolytes, as they depend on different reference-frames. Our model clarifies the exact number of independent parameters and contains simple equations for their conversions.
The generality of the presented framework covers a huge variety of physico-chemical systems and it can be customized by appropriate free energy models. In particular, it applies to electroneutral cell-systems as well as to confined geometries near electrified interfaces. Our transport theory thus provides a valid framework for the description of complex, multicomponent battery-electrolytes in general (concentrated electrolytes, water-in-salt-electrolytes, ILs) as well as for the description of surface-effects of ILs and IL-mixtures.