Transport properties are essential for the understanding and modeling of electrochemical cells, in particular complex systems like lithium-ion batteries. In this study, we demonstrate how a certain degree of freedom in the choice of variables allows us to efficiently determine a complete set of transport properties. We apply the entropy production invariance condition to different sets of electrolyte variables and obtain a general set of formulas. We demonstrate the application of these formulas to an electrolyte typical for lithium-ion batteries, 1M lithium hexafluoro-phosphate in a 1:1 wt. % mixture of ethylene and diethyl carbonates. While simplifications can be introduced, they provide inadequate predictions of conductivity and transport numbers, and we argue that a full matrix of Onsager coefficients is needed for adequate property predictions. Our findings highlight the importance of a complete set of transport coefficients for accurate modeling of complex electrochemical systems and the need for careful consideration of the choice of variables used to determine these properties.

## I. INTRODUCTION

Many central books in electrochemistry, e.g., Ref. 1, contain models for the transport of ions in an electrolyte. Models can also be found in books on irreversible thermodynamics,^{2–4} and they are listed in Handbooks (cf. Ref. 5). In this work, we present a practical procedure, based on existing methods, to find general models for the transport of ions and solvents. Our goal is to provide useful procedures for modeling ion and solvent transport in different electrolyte systems for researchers and practitioners in electrochemistry. To demonstrate the procedures, we shall create a general transport model for the lithium-ion battery electrolyte.

*κ*, is a sum of the single contributions of the positive and negative ions. The mobility

*u*

_{+}(

*u*

_{−}) of the cation (anion) is defined as the absolute velocity of the ion in a unit electric field (dimension m

^{2}V

^{−1}s

^{−1}). For a fully dissociated monovalent electrolyte in water, the conductivity is

*F*is Faraday’s constant, and

*c*is the electrolyte concentration (dimension mol m

^{−3}). The conductivity,

*κ*, has a dimension Ω

^{−1}m

^{−1}. The dimensionless ion transport numbers are accordingly

*τ*

_{i}is the transport number of ion

*i*, defined as the fraction of the electric current carried by

*i*. These expressions result in the diffusion coefficient,

*D*, of the salt in a dilute (ideal) solution,

^{2,6}

*D*is m

^{2}/s. Here,

*R*is the gas constant, and

*T*is the temperature in Kelvin. These equations have been used in numerous applications.

Cation–anion and ion–solvent coupling may, however, be significant, especially when the electrolyte solution is concentrated. The fluctuation-dissipation theorems (FDT), or Green–Kubo relations,^{7–11} are then helpful. These relations take component-self as well as component-component interactions into account. The FDT constitute an essential part of nonequilibrium thermodynamics (NET)—they can be viewed as expressions of Onsager’s regression hypothesis.^{12} The expressions reflect the underlying symmetry implied by the time-reversal invariance of the processes at the particle level. The FDT offers a direct way to systematically extend the simple model of Eqs. (1)–(3), as ion–ion as well as solvent–ion interactions are included. Much work has thus been performed using FDT, particularly on diffusion in isothermal multi-component mixtures (see, e.g., Ref. 13). Work on thermal coupling effects is still scarce, however.

^{2,6,14}the electric conductivity in Eq. (1) changes into

*L*

^{ij}is an Onsager coefficient (to be defined below). Superscript

*ij*refers to the interaction of ions

*i*and

*j*, and Eq. (4) includes coupling between the two ions. The model of independent movement of cations and anions in Eq. (1) is recovered by neglecting cation–anion interactions,

*L*

^{+−}=

*L*

^{−+}= 0, and setting

*L*

^{ii}=

*cu*

_{i}/

*F*. Introducing self-diffusion coefficients,

*D*

_{i,self}=

*u*

_{i}

*RT*/

*F*, while still neglecting cation–anion interactions, the electric conductivity becomes

*κ*have been shown to give very different results when compared to experimental results.

^{15,16}More sophisticated models are, therefore, of interest.

^{11}

In this work, we shall use the fact that the entropy production of the electrolyte is invariant to changes in the variable sets used to describe it and to the choice of frame of reference for the fluxes.^{17} We shall find a set of rules that derive from this, connecting in a systematic way the transport properties of different variable sets. The purpose is to present practical procedures for the determination of two coefficient sets: The first set relates better to the electrolyte structure, while the second set addresses the measurement situation better. It is well known that a change in the frame of reference of fluxes creates a link between coefficient sets (see, e.g., Ref. 6). Here, we shall see that certain links also appear between different sets of electrolyte variables. These links can be combined with FDT to find transport properties in a systematic manner. By exploiting these links, or rules, as we shall call them, we hope to add deeper insight into the processes in a typical electrolyte and set the stage for modeling profiles in intensive and other variables across the electrolyte. As an example to illustrate the theoretical elaborations, we have chosen an important lithium-ion battery electrolyte.

The typical electrolyte consists of 1M of a lithium salt, here lithium hexafluorophosphate, and two solvent molecules, here ethylene carbonate and diethyl carbonate, in a weight ratio of 1:1%^{18} in most cases. We will later in this paper present results from molecular dynamics simulations of this electrolyte, but a simulation snapshot of the complicated mixture is shown already in Fig. 1. As others have also observed, there is a tendency toward solvent-separated ion pairing.^{19,20} Clearly, there are interactions of several sorts in the electrolyte.

Common acronyms for the carbonates are EC (ethylene carbonate) and DEC (diethyl carbonate) respectively, and we shall also use these. When used in subscripts, however, we shall shorten their names further to E and D. The salt as a component is denoted by L only to be distinguished from the single ions Li^{+} and $PF6\u2212$. The components of the electrolyte are illustrated in Fig. 2. We see that DEC is larger than EC. Both molecules have polar carbonate groups with one oxygen pointing out. The carbonates play different roles in the shielding of the lithium-ion. It is shielded by about three DEC molecules in this snapshot, and a relatively large distance between the cation and the anion is noticeable. At the same time, we also see tendencies toward clustering of the two ions.

There is a certain degree of freedom in the choice of variables that can be used to describe the electrolyte, and we consider two different variable sets for our purpose. The first set uses ions and neutral molecules mixed together. The components are Li^{+}, $PF6\u2212$, EC, and DEC. The other set has only electroneutral components, namely components L, EC, and DEC.^{17} The first set is more common. It connects more directly to the structure of the electrolyte and is suitable for use with FDT. The second set contains the independent components according to the phase rule.^{21} This set is better connected to possible measurements, as their operational definitions demonstrate. We refer to the formulas derived for these two sets as the mixed component scenario and the neutral component scenario, respectively.

It is well known that mass transfer limitations and subsequent polarization can occur in the battery electrolyte.^{16,22–25} Therefore, precise knowledge of coefficients for the transport of mass, charge, and heat is central to battery modeling. This work aims to demonstrate how we can find useful transport coefficients in a practical way using specific invariance criteria (rules). The procedure will be illustrated with data for the *isothermal electrolyte*. The principles are general in that the procedures should apply to other electrochemical systems, such as batteries or electrolysis cells. Nonisothermal results are presently unavailable from FDT; however,^{26} we shall indicate how they can be included.

## II. THEORETICAL BASIS

### A. Entropy production in mixed and neutral component scenarios

We assume some familiarity with the theory of NET, the basis of our derivations. For more explanations, see Refs. 6 and 12. We repeat only the basic assumptions as follows:

The theory of NET assumes that the system obeys local equilibrium, meaning that the Gibbs equation (the internal energy written as a total differential of the relevant extensive variables) applies locally to any volume element, even though the total system is not in global equilibrium.

The balance laws for mass, energy, and momentum are written for the independent variables of interest and included in Gibbs’ equation.

The entropy production and the entropy flux are identified by comparing the result to the entropy balance.

The linear laws of transport follow from the entropy production for any volume element in the system.

The fluxes are related to the forces by a set of Onsager coefficients. The linear laws can be used to find the transport coefficients, but they can equally well be found from FDT, which is an integral part of NET.

The entropy production is absolute, but there is a degree of freedom in the choice of variables used to describe it. In addition, the chosen set of fluxes depends on the frame of reference used in their determination.

The NET theory, which originates from Onsager, assumes that particle fluctuations obey the regression hypothesis or, in essence, FDT. Onsager used microscopic reversibility to show that the matrix of coefficients, relating the fluxes to all forces, was symmetric.

^{6,12}

These premises have been shown to hold up exceptionally well. The assumption of local equilibrium was recently proven valid for shock waves.^{27} The assumption can be formulated for reactions using mesoscopic variables.^{28,29}

Electrochemical cells are heterogeneous, so a full system analysis starts by finding the entropy production of three bulk systems and two electrode interfaces. The bulk systems are the bulk of the electrodes and the bulk electrolyte. Kjelstrup and her co-workers showed how to integrate across the system.^{6,17} Flux-force relations are obtained for each part. Electrolyte boundary conditions were given by the processes at the electrode surfaces.

Consider the electrolyte of a lithium battery in a cell with metal lithium electrodes. The electrolyte was pictured in Fig. 1. In the mixed component scenario (properties have a superscript C), there are four component fluxes in addition to the heat flux. The electric current is a dependent variable. In the neutral component scenario (properties have a superscript N), there are three component fluxes and one variable electric current. Its value is measured (and controlled) in the external circuit.^{21} In addition, there is a heat flux. In both cases, there are five contributions to the entropy production of the electrolyte. These descriptions are equivalent in the sense that they give the same entropy production.

Here, we shall take advantage of this invariance of entropy production and examine the rules that bind the particular sets of coefficients together. We apply the condition to the two scenarios mentioned because one is practical for computational reasons (the mixed component scenario) and the other is suited to decompose experimental data (the neutral component scenario). There are rules connecting the sets of coefficients, which can help create models for transport coefficients.

*σ*, is

^{6}

^{,}

^{−2}s

^{−1}), with conjugate force ∇(1/

*T*). Component fluxes, indicated with subscripts Li

^{+}and $PF6\u2212$, have as a driving force the negative electrochemical potential gradient divided by the temperature, i.e., $\u2212\u2207\mu \u0303Li+,T/T$ and $\u2212\u2207\mu \u0303PF6\u2212,T/T$,

^{2,3,21}where the electrochemical potential for ion

*i*is defined by

^{30}

^{,}

*z*

_{i}is the charge number of ion

*i*,

*μ*

_{i}its chemical potential, and

*ψ*is the Maxwell or electrostatic potential. The electrochemical potential gradients in Eq. (6) are all evaluated at a constant temperature. Superscript B indicates from now on that the barycentric frame of reference is used for the mixed component scenario. This choice is natural when we want to use molecular dynamics simulations to provide the set of Onsager coefficients that belong to Eq. (6). Fong

*et al.*

^{11}derived the same entropy production using a different set of premises.

^{3,21}the five independent contributions to

*σ*are

The electric current density, *j*, (in A/m^{2}) is the flux conjugate to −∇*φ*/*T*, where ∇*φ* is the electric potential gradient that we can measure with lithium-reversible electrodes (not to be confused with the Maxwell potential gradient). The heat flux in the neutral component scenario, $Jq\u2032N$, is measurable and has, as the conjugate driving force, the gradient of the inverse temperature 1/*T*. The driving force for components L, DEC, and EC, respectively, is the negative chemical potential gradient divided by the temperature, evaluated at constant temperature.

The measurable heat flux $Jq\u2032N$ is defined as the energy flux *J*_{q} minus the latent heat, *H*_{i} (partial molar enthalpy), carried by the three components: $Jq\u2032N=Jq\u2212\u2211j=13JjHj$.^{12} The heat fluxes in the two scenarios were shown to be related.^{17} When lithium electrodes are used, we obtain $Jq\u2032C=Jq\u2032N+TSLi+(j/F)$. The mass fluxes depend on the frame of reference.

^{31}The FDT is well-defined in the mixed component scenario, where all component positions and velocities are uniquely defined. The FDT technique is well-established, but at present, only for isothermal systems. The system used in this work to give numerical insight is thus isothermal. In Appendix B and below, we explain how to find the Onsager coefficients Λ

_{ij}that belong to the fluxes and forces of Eq. (6) by tracking particle positions and velocities. The coefficients have the dimension m

^{2}s

^{−1}, typical for Fick’s diffusion coefficients. This dimension can be converted to a dimension proper for the flux-force relations (see Appendix A) by multiplying with the factor

*c*/

*R*, as explained by Krishna and van Baten.

^{9}The data reduction procedure starts by computing the average component velocity

**u**

_{i}(in m/s) of all species

*i*from the particle velocities

**v**

_{k,i}. We have $Niui(t)=\u2211k=1Nivk,i(t)$. In the barycentric frame of reference, the isothermal flux of

*i*obeys

*n*= 4 components, and Λ

_{ij}is an element in the matrix of Onsager coefficients. The molar flux of a component is found by multiplying Eq. (9) with

*c*, the total concentration,

_{ij}have dimension m

^{2}s

^{−1}, while the corresponding Λ

^{ij}obtains the dimension needed in the thermodynamic flux-force relations of Appendix A when we multiply with the factor

*c*/

*R*,

### B. Reducing the number of variables: From the barycentric to the solvent frame of reference

^{6}We want to study solvent segregation and choose one of the solvent components as the frame of reference: the EC component. We replace the chemical potential gradient of EC in both scenario formulations. In mixed as well as neutral component scenarios, Gibbs–Duhem’s equation takes the form

*x*

_{i}is the mole fraction of L or D. We introduce the respective expression in the corresponding expressions for the entropy production and obtain

^{ij}, has the B-frame of reference. The set of

*L*

^{ij}-coefficients with the E frame of reference is computed from these.

### C. Flux-force relations for mixed and neutral component scenarios

A coefficient superscript indicates the interactions among the fluxes in question; cf. *L*^{ij}. The following dimensions apply to the Onsager coefficients in the mixed component scenario:

*L*^{qq}has dimension K J m^{−1}s^{−1}, while the coupling coefficients*L*^{q+},*L*^{q−}, and*L*^{qD}have dimension K mol m^{−1}s^{−1}.*L*^{++},*L*^{+−},*L*^{+D},*L*^{−−},*L*^{−D}, and*L*^{DD}have dimensions K mol^{2}J^{−1}m^{−1}s^{−1}.

Double superscripts indicate interacting phenomena. The matrix of coefficients is symmetric, according to Onsager. For isothermal systems, it is common to adsorb the constant 1/*T* into the coefficient.

*L*_{qq}has dimension K J m^{−1}s^{−1}.*L*_{qL}and*L*_{qD}have dimensions K mol s^{−1}m^{−1}.*L*_{qφ}has dimension K C m^{−1}s^{−1}.

*L*

_{qq}differs from

*L*

^{qq}because it refers to a different heat flux. The coefficients of the second and third rows have pairwise the same dimensions. The first pair

*L*

_{Lq}and

*L*

_{Dq}has dimension K mol m

^{−1}s

^{−1}. The second and third pairs,

*L*

_{LL},

*L*

_{LD}, and

*L*

_{LD},

*L*

_{DD}, have dimensions K mol

^{2}J

^{−1}m

^{−1}s

^{−1}. The last coefficients,

*L*

_{L}

*φ*and

*L*

_{D}

*φ*, have dimensions K C mol J

^{−1}m

^{−1}s

^{−1}. The dimension of

*L*

_{φq}is K C m

^{−1}s

^{−1}, the dimension of

*L*

_{φL}and

*L*

_{φD}is K C mol J

^{−1}m

^{−1}s

^{−1}, and the dimension of

*L*

_{φφ}is K C V

^{−1}m

^{−1}s

^{−1}(see Appendix A). We include the common factor 1/

*T*into the coefficients when the temperature is constant.

In order to distinguish the coefficients of the two scenarios, we have used subscripts for the interacting phenomena in the neutral component scenario rather than superscripts as in the mixed component scenario.

### D. Operationally defined properties

The properties relevant for thermodynamic modeling of electrochemical systems are mostly obtained from experiments, i.e., they are operationally defined. The properties include diffusion coefficients, transference coefficients, electric conductivity, heats of transfer, Peltier coefficients, and the chemical potential. The operational definitions are, therefore, briefly recapitulated.

*φ*in terms of

*j*with Eq. (18) and find

*ℓ*

_{ij}. The coefficients are measured in the absence of

*j*and ∇

*φ*

^{21}and with the E frame of reference.

*t*

_{i}, are defined by Hittorf experiments from Eq. (18),

^{21}It must not be confused with the transport number, cf. Eq. (2), the fraction of the electric current carried by one ion. For the co-solvent DEC, the transference coefficient describes the electro-osmosis of DEC and, therefore, solvent segregation. This coefficient is also accessible via Onsager relations and

*emf*-measurements.

^{21}

*coefficient*,

*π*

^{N}, is defined as the ratio of the heat flux and the electric current density,

*heat*of the junction, which constitutes the electrode interface.

^{32}The Peltier heat is accessible via the Onsager relation and the Seebeck coefficient of the cell. The symmetry relation gives

*T*, can be determined in, for instance, nonequilibrium molecular dynamics simulations. The gradient in chemical potential, ∇

*μ*

_{i,T}for

*i*= L or DEC, can be found from knowledge of the composition

*c*

_{i}and the corresponding thermodynamic factor, cf. Eq. (27).

We have seen how the coefficients in the equations earlier are often determined. There are two major routes. In Subsection II B, we referred to the fluctuation-dissipation theorems (FDT) (available for isothermal systems; however, see Ref. 26). In this section, we have outlined how the coefficients can be obtained from their operational definitions. Both methods should give the same results, providing us with a useful check for consistency.

### E. Fick’s law’s coefficients and the chemical potential

The symmetry of the Onsager coefficients on the macroscopic scale reflects the underlying symmetry on the molecular level and is an example of the many symmetry laws in physics.^{8} Therefore, the Onsager coefficients are of primary interest, not only for the present case but for electrochemical systems in general. The diffusion coefficients of Fick are, however, more frequently used in practice,^{13} and we will here comment on their relation to the Onsager coefficients.

*T*= 0,

*j*= 0 to Eq. (19) and obtain

^{12}

^{,}

*c*

_{j}is the concentration of the (independent) components

*j*= L and DEC. The non-ideality is expressed here using the thermodynamic factors Γ

_{ij},

^{33–35}

^{,}

**c**

_{i}indicates that the derivative is taken while holding the concentrations of all components, except

*i*, constant. The assumption of an ideal solution corresponds to Γ

_{ij}=

*δ*

_{ij}, where

*δ*

_{ij}is the Kronecker delta. Introducing the chemical potentials from Eq. (27) into Eq. (26) gives

*T*from the force disappears from the coefficients of Eq. (22) when the concentration gradient is introduced. However, the introduction of the gradients in

*c*

_{i}and

*c*

_{j}and the corresponding thermodynamic factor Γ

_{ij}means that the matrix of coefficients is no longer symmetric.

^{35}We obtain expressions for Fick’s diffusion coefficients

*D*

_{ij},

*l*

_{ii}

*l*

_{jj}−

*l*

_{ij}

*l*

_{ji}≥ 0,

*l*

_{ii}> 0. Knowledge of Γ

_{ij}/

*c*

_{i}and

*ℓ*

_{ij}gives the matrix elements

*D*

_{ij}. In Appendix C, we describe how we calculate the thermodynamic factors Γ

_{ij}from molecular dynamics simulations.

In the present case, it is practical to use FDT [ Appendix B, Eq. (B1)] for the mixed component scenario (because all particle positions are well-definable). The target is the Onsager coefficients of the neutral component scenario. We shall next see how these can be obtained without the knowledge of Fick’s coefficients.

## III. IMPOSING ENTROPY PRODUCTION INVARIANCE

^{6,17,21}The supporting electrodes of the cell in question are both reversible to lithium. The relations follow from the expression for the electric current density:

### A. Rules for coupling of fluxes

*φ*. For mixtures of multivalent electrolytes, the valency of the respective ion enters. For mass transport, the corresponding coefficient relations are

### B. Implications of the rules for coupling of fluxes

*κ*, prescribed by the rules is the same. We repeat

*L*

^{+−}can be positive or negative in sign. Coupling may, therefore, enhance or reduce

*κ*.

*i*in Eq. (2) is the fraction of the electric current carried by the ion in question. In the presence of lithium-reversible electrodes, the transference coefficient of L can be understood as the negative transport number of the anion, $PF6\u2212$. This expresses that salt accumulates on the left-hand side of the cell when a positive electric current is passing from left to right in the cell. Coupling can reduce or increase the transference coefficient. The transference coefficient applies to

*neutral*components. For the DEC component, it becomes

^{21}Therefore, coupling coefficients are likely to be significant in strong electrolyte solutions. The case of the independent movement of ions was referred to in the Introduction.

^{1,2,6}It is a limiting case, obtained by setting

*L*

^{+−}= 0,

*L*

^{−D}= 0, and

*L*

^{+D}= 0.

^{6}

### C. The nonisothermal case. An approximation

Consider again the coefficients that relate to heat transport. While the thermal conductivity of a system can be found with well-established techniques, experimental as well as computational, the procedure to find coupling coefficients for heat and mass and charge, heat, and mass is less well established. A complicating factor is the definition of the heat flux. The heat flux $Jq\u2032N$ is measurable in the neutral component scenario. From the discussion below Eq. (8), we understand that $Jq\u2032C$ in the mixed component scenario is not measurable. Both scenarios are now discussed in view of the Rules for Coupling of Fluxes to gain more insight into the conditions for their determinations.

*J*

_{D}= 0 and $JPF6\u2212=0$. These stationary state conditions can be used to eliminate two driving forces from the flux equations. The heat flux of the mixed component description becomes

*heat of transfer*of the lithium-ion in the mixed component scenario,

*K*and

*M*are lumped coefficients, which are functions of the Onsager coefficients,

*L*

^{+−}= 0,

*L*

^{−D}= 0, and

*L*

^{+D}= 0, and obtain

*K*= 0 and

*M*= 0. The expressions for the heats of transfer simplify to

*L*

_{φq}=

*L*

^{q+}−

*L*

^{q−}, can now be applied. We divide both sides by

*L*

_{φφ}(

*L*

^{++}+

*L*

^{−−}) and obtain

The measurable Peltier heat may obtain a simplified contribution through this.

### D. Summary of procedures for coefficient determinations

We have described the electrolyte in a concentration cell in two equivalent ways: as a set of mixed components and as a set of neutral components only. A procedure for determining practical transport properties has been outlined. A schematic illustration of the procedure is given in Fig. 3, and the procedure consists of four steps, as follows:

Determine the coefficients of the mixed component scenario with FDT in the barycentric frame of reference.

Reduce the number of coefficients by replacing the barycentric frame of reference with the (co-)solvent frame of reference.

Apply the Rules for Coupling of Fluxes to find the coefficients of the neutral component scenario from the coefficients of the mixed component scenario.

Finally, determine the set of practical transport properties by applying the operational definitions to the coefficients obtained in the previous step.

In the following, we shall apply these theoretical steps and demonstrate how they apply to an example. As of yet, the set of Rules can only be used for isothermal systems. FDT for nonisothermal systems is still lacking in the practical literature.

## IV. SIMULATION METHODS

To apply the theoretical framework, we have calculated transport coefficients with equilibrium molecular dynamics simulations, and we provide an overview of our simulation methodology in the following.

### A. Electrolyte model and simulation set-up

The equilibrium MD simulations were performed using the LAMMPS^{31} software. Atomic and molecular interactions were described by the OPLS-AA^{36} potential. The parameters for the solvent molecules were obtained from the LigParGen web service,^{37–39} and the parameters for Li^{+} and $PF6\u2212$ were taken from Jensen and Jorgensen^{40} and Acevedo *et al.*,^{41,42} respectively. The real-space cutoff for the Lennard-Jones and Coulombic forces was set to 13 Å. A Lennard-Jones tail correction was added to the energy and pressure.^{43} Electrostatic forces beyond the cutoff were computed using a particle–particle particle-mesh solver^{44} with a maximum relative error in forces of 1 × 10^{−6}. The ionic charges were scaled by a factor of 0.75 to account for the overestimation of electrostatic interactions between ions in non-polarizable force fields.^{45} Periodic boundary conditions were applied in all directions. Initial configurations of the systems were prepared by randomly placing solvent molecules, Li^{+} and $PF6\u2212$, in a simulation box with the Packmol software.^{46} The composition of the electrolyte corresponded to the weight ratio EC:DEC = 1:1, which is a particle ratio of 5520:4116, and 920 particles of salt were added to reach a concentration of 1M.

The energies of the systems were minimized to avoid particle overlap. We used the routine developed by Molinari *et al.*^{47} for initial equilibration. The systems were further equilibrated at a temperature of 350 K or higher and a pressure of 1 bar in the isobaric-isothermal (*NPT*) ensemble using a time-step of 1.25 fs in order for the potential energy and density of the systems to stabilize. The Nosé–Hoover thermostat and barostat^{48–50} were used to control the temperature and pressure using time constants of 100 and 1000 time-steps, respectively. The final equilibration in the *NPT* ensemble was conducted at a temperature of 300 K while sampling the box volume. The simulation box size was adjusted (to the average volume) after equilibrium to obtain the correct density. The simulation boxes were cubic with sides of 115 Å.

### B. Calculation of transport properties

We sampled the radial distribution functions (RDFs) and transport properties in the canonical ensemble (*NVT*) using a time-step of 1.25 fs during simulations lasting at least 80 ns. The Nosé–Hoover thermostat maintained the temperature, which was set to 300 K. Sampling transport properties in the *NVT* ensemble with the Nosé–Hoover thermostat produces transport properties that are statistically indistinguishable from those obtained using the *NVE* ensemble.^{51} The average pressure was 0.4 bar for these simulations. All RDFs and coefficients were obtained from the OCTP code^{52} by tracking the motion of the central atom of the various components.

We computed the coefficients *L*^{ij} of Eq. (17) in the mixed component scenario and the matrix of *ℓ*_{ij}-coefficients of Eq. (19) in the neutral component scenario. We have assumed monovalent ions in our theoretical framework, where the cation (anion) has a valency of +1 (−1). As described earlier, the ion charges are scaled to *z*_{+} = *z* = 0.75 and *z*_{−} = −*z* = −0.75, respectively, in our simulations. This will modify the equations presented in Sec. III via the inclusion of the valency in Eqs. (31)–(34), as well as in the equations that follow. In particular, each term in the expression for *L*_{φφ} in Eq. (36) must be multiplied by *z*^{2}, while the terms in Eqs. (37)–(39) are multiplied by *z*. This also implies that the transference number in Eq. (44) should be divided by *z*. Finally, we calculate the transport numbers according to the derivation of Gullbrekken *et al.*,^{53} which takes the charge into account.

To determine the thermodynamic factors, we evaluated the Kirkwood–Buff integrals^{33} from the RDFs. The relations between the Kirkwood–Buff integrals and the thermodynamic factors are given in Appendix C. We used the finite-size correction of Ganguly and van der Vegt^{54} for the RDFs and the integration procedure of Krüger *et al.*^{55} in our computations. We refer to Milzetti *et al.* for further details on this procedure.^{56}

## V. RESULTS AND DISCUSSION

### A. Equilibrium properties

The thermodynamic factors give information on attractive and repulsive forces in the electrolyte. Thermodynamic factors (based on mole fractions) of the example electrolyte are presented in Table I (see also Appendix C for a description of the calculation and the concentration-based thermodynamic factors). The thermodynamic factors for the salt and DEC indicate that there are repulsive forces between the components themselves. The conditions are far from ideal, with values for Γ_{LL} and Γ_{DD} of 1.5 and 1.2, respectively (see the two first rows of Table I). The cross terms Γ_{LD} and Γ_{DL} indicate that there are attractive forces between L and DEC embedded in co-solvent EC. It is known that DEC arranges itself around the lithium-ion,^{19,20} as we noted in connection with Fig. 1. A clustering of ion pairs has also been observed. The results in Table I support this.

### B. A comprehensive set of Onsager coefficients

Tables II–V present our comprehensive set of Onsager coefficients for the ternary lithium battery electrolyte that we have considered. This set includes all coefficients for the isothermal electrolyte. Transport coefficients in the nonisothermal electrolyte are yet to be calculated, but we have provided preliminary expressions. Earlier studies have been less comprehensive.^{16,23,57,58}

. | Frame of reference . | ||
---|---|---|---|

Coefficient . | B × 10^{−11} m^{2}/s
. | E × 10^{−11} m^{2}/s
. | E × c/RT 10^{−11} mol^{2}/(J m s)
. |

L^{++} | 0.4 ± 0.04 | 0.8 ± 0.1 | 4.0 ± 0.3 |

L^{−−} | 0.7 ± 0.2 | 1.3 ± 0.2 | 6.2 ± 1.0 |

L^{+−} | 0.1 ± 0.02 | 0.6 ± 0.05 | 2.8 ± 0.2 |

L^{D+} | 0.5 ± 0.1 | 2.4 ± 0.2 | 12.1 ± 1.1 |

L^{D−} | −0.3 ± 0.1 | 1.8 ± 0.2 | 9.2 ± 0.9 |

L^{DD} | 2.8 ± 0.2 | 11.3 ± 1.3 | 55.9 ± 6.2 |

L^{E+} | −0.9 ± 0.1 | ||

L^{E−} | −1.0 ± 0.1 | ||

L^{EE} | 6.1 ± 0.9 | ||

L^{ED} | −3.4 ± 0.4 |

. | Frame of reference . | ||
---|---|---|---|

Coefficient . | B × 10^{−11} m^{2}/s
. | E × 10^{−11} m^{2}/s
. | E × c/RT 10^{−11} mol^{2}/(J m s)
. |

L^{++} | 0.4 ± 0.04 | 0.8 ± 0.1 | 4.0 ± 0.3 |

L^{−−} | 0.7 ± 0.2 | 1.3 ± 0.2 | 6.2 ± 1.0 |

L^{+−} | 0.1 ± 0.02 | 0.6 ± 0.05 | 2.8 ± 0.2 |

L^{D+} | 0.5 ± 0.1 | 2.4 ± 0.2 | 12.1 ± 1.1 |

L^{D−} | −0.3 ± 0.1 | 1.8 ± 0.2 | 9.2 ± 0.9 |

L^{DD} | 2.8 ± 0.2 | 11.3 ± 1.3 | 55.9 ± 6.2 |

L^{E+} | −0.9 ± 0.1 | ||

L^{E−} | −1.0 ± 0.1 | ||

L^{EE} | 6.1 ± 0.9 | ||

L^{ED} | −3.4 ± 0.4 |

Coefficient . | Equation . | Value . | Unit . |
---|---|---|---|

L_{φφ} | (36) | 0.23 ± 0.03 | Ω^{−1} m^{−1} |

L_{φL} | (37) | −2.5 ± 0.8 | 10^{−6} mol C/(J m s) |

L_{φD} | (38) | 2.1 ± 0.8 | 10^{−6} mol C/(J m s) |

L_{LL} | (40) | 6.2 ± 1.0 | 10^{−11} mol^{2}/(J m s) |

L_{DL} | (41) | 9.2 ± 0.9 | 10^{−11} mol^{2}/(J m s) |

L_{DD} | (42) | 55.9 ± 6.2 | 10^{−11} mol^{2}/(J m s) |

ℓ_{LL} | (20) | 3.7 ± 0.2 | 10^{−11} mol^{2}/(J m s) |

ℓ_{DL} | (20) | 11.3 ± 0.6 | 10^{−11} mol^{2}/(J m s) |

ℓ_{DD} | (20) | 53.7 ± 6.4 | 10^{−11} mol^{2}/(J m s) |

Coefficient . | Equation . | Value . | Unit . |
---|---|---|---|

L_{φφ} | (36) | 0.23 ± 0.03 | Ω^{−1} m^{−1} |

L_{φL} | (37) | −2.5 ± 0.8 | 10^{−6} mol C/(J m s) |

L_{φD} | (38) | 2.1 ± 0.8 | 10^{−6} mol C/(J m s) |

L_{LL} | (40) | 6.2 ± 1.0 | 10^{−11} mol^{2}/(J m s) |

L_{DL} | (41) | 9.2 ± 0.9 | 10^{−11} mol^{2}/(J m s) |

L_{DD} | (42) | 55.9 ± 6.2 | 10^{−11} mol^{2}/(J m s) |

ℓ_{LL} | (20) | 3.7 ± 0.2 | 10^{−11} mol^{2}/(J m s) |

ℓ_{DL} | (20) | 11.3 ± 0.6 | 10^{−11} mol^{2}/(J m s) |

ℓ_{DD} | (20) | 53.7 ± 6.4 | 10^{−11} mol^{2}/(J m s) |

Coefficient . | All types of coupling . | No cation–anion coupling . | Self-diffusion coefficients . |
---|---|---|---|

κ/Ω^{−1} m^{−1} | 0.23 ± 0.03 | 0.27 ± 0.03 | 0.41 ± 0.01 |

t_{L} | −0.97 ± 0.12 | −0.81 ± 0.07 | −0.84 ± 0.01 |

t_{D} | 0.90 ± 0.46 | 0.39 ± 0.17 | 0 |

τ_{+} | 0.28 ± 0.09 | 0.39 ± 0.05 | 0.37 ± 0.01 |

τ_{−} | 0.72 ± 0.09 | 0.61 ± 0.05 | 0.63 ± 0.01 |

Coefficient . | All types of coupling . | No cation–anion coupling . | Self-diffusion coefficients . |
---|---|---|---|

κ/Ω^{−1} m^{−1} | 0.23 ± 0.03 | 0.27 ± 0.03 | 0.41 ± 0.01 |

t_{L} | −0.97 ± 0.12 | −0.81 ± 0.07 | −0.84 ± 0.01 |

t_{D} | 0.90 ± 0.46 | 0.39 ± 0.17 | 0 |

τ_{+} | 0.28 ± 0.09 | 0.39 ± 0.05 | 0.37 ± 0.01 |

τ_{−} | 0.72 ± 0.09 | 0.61 ± 0.05 | 0.63 ± 0.01 |

Coefficient . | Equation . | Value × 10^{−11} m^{2} s^{−1}
. | Value × 10^{11} mol^{2} J^{−1} m^{−1} s^{−1}
. |
---|---|---|---|

ℓ_{LL} | (20) | 0.7 ± 0.04 | 3.7 ± 0.2 |

ℓ_{DL} | (20) | 2.3 ± 0.1 | 11.3 ± 0.6 |

ℓ_{DD} | (20) | 10.8 ± 1.3 | 53.7 ± 6.4 |

D_{LL} | (29) | 52.3 ± 2.6 | |

D_{LD} | (29) | 309.4 ± 16.9 | |

D_{DL} | (29) | 222.8 ± 22.4 | |

D_{DD} | (29) | 1453.5 ± 170.2 | |

$DLi+$ | 7.2 ± 0.2 | ||

$DPF6\u2212$ | 12.4 ± 0.2 | ||

D_{D} | 11.6 ± 0.1 | ||

D_{E} | 22.1 ± 0.1 |

Coefficient . | Equation . | Value × 10^{−11} m^{2} s^{−1}
. | Value × 10^{11} mol^{2} J^{−1} m^{−1} s^{−1}
. |
---|---|---|---|

ℓ_{LL} | (20) | 0.7 ± 0.04 | 3.7 ± 0.2 |

ℓ_{DL} | (20) | 2.3 ± 0.1 | 11.3 ± 0.6 |

ℓ_{DD} | (20) | 10.8 ± 1.3 | 53.7 ± 6.4 |

D_{LL} | (29) | 52.3 ± 2.6 | |

D_{LD} | (29) | 309.4 ± 16.9 | |

D_{DL} | (29) | 222.8 ± 22.4 | |

D_{DD} | (29) | 1453.5 ± 170.2 | |

$DLi+$ | 7.2 ± 0.2 | ||

$DPF6\u2212$ | 12.4 ± 0.2 | ||

D_{D} | 11.6 ± 0.1 | ||

D_{E} | 22.1 ± 0.1 |

This paper’s purpose is not to provide coefficients for a particular electrolyte but to present and document a convenient procedure for coefficient determination. As such, the physical-chemical meaning of the results and their application to battery modeling will not be the focus of our discussion here. However, using the lithium battery electrolyte as an example, we will now discuss the single steps of the new procedure and illustrate how they may be used. Most of the tools in this context are familiar and well-established, but a new item has been added: the Rules for Coupling of Fluxes. These Rules allow us to effectively predict measurable transport coefficients from sets of coefficients from FDT, which are relatively easy to compute and better related to electrolyte structure. In this manner, theory, simulation, and experiment go hand-in-hand in one tool.

The set of Onsager coefficients for isothermal conditions, presented in Tables II–V, demonstrate the advantages of such a tool. The first conclusion we can draw is that all coefficients in any of the matrices we have presented are significant. Coupling coefficients should, in general, be accounted for (see Table III).

Name . | Symbol . | Dimension according to Eqs. (17) or (19) . |
---|---|---|

Entropy production | σ | J K^{−1} m^{−3} s^{−1} |

Heat flux | $Jq\u2032$ | J m^{−2} s^{−1} |

Thermal force | ∇(1/T) | K^{−1} m^{−1} |

Mass flux | J_{i} | mol m^{−2} s^{−1} |

Chemical force | −∇μ_{i}/T | J K^{−1} mol^{−1} m^{−1} |

Electric flux | j | C m^{−2} s^{−1} |

Electric force | −∇φ/T | V K^{−1} m^{−1} |

Main heat coefficient | L_{qq} | K J m^{−1} s^{−1} |

Coupling coefficient, heat L | L_{qL} = L_{Lq} | K mol m^{−1} s^{−1} |

Coupling coefficient, heat D | L_{qD} = L_{Dq} | K mol m^{−1} s^{−1} |

Electric conductivity | L_{φφ} | K C V^{−1} m^{−1} s^{−1} |

Coupling coefficient, charge-heat | L_{φq} = L_{qφ} | K C m^{−1} s^{−1} |

Coupling coefficient, charge-mass L | L_{φL} = L_{Lφ} | K C mol J^{−1} m^{−1} s^{−1} |

Coupling coefficient, charge-mass D | L_{φD} = L_{Dφ} | K C mol J^{−1} m^{−1} s^{−1} |

Diffusion main coefficient L | ℓ_{LL} | K mol^{2} J^{−1} m^{−1} s^{−1} |

Diffusion main coefficient D | ℓ_{DD} | K mol^{2} J^{−1} m^{−1} s^{−1} |

Diffusion coupling coefficient LD | ℓ_{LD} = ℓ_{LD} | K mol^{2} J^{−1} m^{−1} s^{−1} |

Diffusion main coefficient ++ | L^{++} | K mol^{2} J^{−1} m^{−1} s^{−1} |

Diffusion main coefficient −− | L^{−−} | K mol^{2} J^{−1} m^{−1} s^{−1} |

Diffusion coupling coefficient +− | L^{+−} = L^{−+} | K mol^{2} J^{−1} m^{−1} s^{−1} |

Diffusion coupling coefficient D^{+} | L^{D+} = L^{+D} | K mol^{2} J^{−1} m^{−1} s^{−1} |

Diffusion coupling coefficient D− | L^{D−} = L^{−D} | K mol^{2} J^{−1} m^{−1} s^{−1} |

Diffusion main coefficient DD | L^{DD} | K mol^{2} J^{−1} m^{−1} s^{−1} |

Name . | Symbol . | Dimension according to Eqs. (17) or (19) . |
---|---|---|

Entropy production | σ | J K^{−1} m^{−3} s^{−1} |

Heat flux | $Jq\u2032$ | J m^{−2} s^{−1} |

Thermal force | ∇(1/T) | K^{−1} m^{−1} |

Mass flux | J_{i} | mol m^{−2} s^{−1} |

Chemical force | −∇μ_{i}/T | J K^{−1} mol^{−1} m^{−1} |

Electric flux | j | C m^{−2} s^{−1} |

Electric force | −∇φ/T | V K^{−1} m^{−1} |

Main heat coefficient | L_{qq} | K J m^{−1} s^{−1} |

Coupling coefficient, heat L | L_{qL} = L_{Lq} | K mol m^{−1} s^{−1} |

Coupling coefficient, heat D | L_{qD} = L_{Dq} | K mol m^{−1} s^{−1} |

Electric conductivity | L_{φφ} | K C V^{−1} m^{−1} s^{−1} |

Coupling coefficient, charge-heat | L_{φq} = L_{qφ} | K C m^{−1} s^{−1} |

Coupling coefficient, charge-mass L | L_{φL} = L_{Lφ} | K C mol J^{−1} m^{−1} s^{−1} |

Coupling coefficient, charge-mass D | L_{φD} = L_{Dφ} | K C mol J^{−1} m^{−1} s^{−1} |

Diffusion main coefficient L | ℓ_{LL} | K mol^{2} J^{−1} m^{−1} s^{−1} |

Diffusion main coefficient D | ℓ_{DD} | K mol^{2} J^{−1} m^{−1} s^{−1} |

Diffusion coupling coefficient LD | ℓ_{LD} = ℓ_{LD} | K mol^{2} J^{−1} m^{−1} s^{−1} |

Diffusion main coefficient ++ | L^{++} | K mol^{2} J^{−1} m^{−1} s^{−1} |

Diffusion main coefficient −− | L^{−−} | K mol^{2} J^{−1} m^{−1} s^{−1} |

Diffusion coupling coefficient +− | L^{+−} = L^{−+} | K mol^{2} J^{−1} m^{−1} s^{−1} |

Diffusion coupling coefficient D^{+} | L^{D+} = L^{+D} | K mol^{2} J^{−1} m^{−1} s^{−1} |

Diffusion coupling coefficient D− | L^{D−} = L^{−D} | K mol^{2} J^{−1} m^{−1} s^{−1} |

Diffusion main coefficient DD | L^{DD} | K mol^{2} J^{−1} m^{−1} s^{−1} |

### C. Reducing the number of variables

The first dataset is the simulation results from FDT, obtained in the barycentric (B) frame of reference. They are listed in Table II, column one. We see a large number of coefficients—ten altogether. A large number of coefficients is clearly a disadvantage, and the first task is to reduce this number by transforming to the solvent (E) frame of reference; see the description in Sec. II B. Here, we choose the more mobile (which has less impact on the electrolyte structure) of the two organic solvents, EC, as the reference. After the change in frame of reference, we obtain the results in the second and third columns of Table II. The number of necessary coefficients changes from ten to six upon the shift to the E frame of reference, and the coefficient values change accordingly.

The frame of reference is always an issue when transport coefficients are concerned. Standard procedures are available on how to change from one frame to another (see, e.g., Ref. 6). Here, the choice of the frame of reference determines the contributions to the *emf* but not the (total) *emf* (which is invariant to the choice^{17,59}).

The results presented in Table II have the same order of magnitude as those published in the literature, 10^{−10} m^{2} s^{−1}.^{58} So far, EC and DEC have been treated as one solvent (however, see Wang *et al.*^{22}) while we distinguish between the co-solvents. A direct comparison of our results with those in the literature is thus not possible. However, all coupling coefficients are large in both frames of reference. This is an interesting point; it means that it will not be correct to neglect coefficients.

Return again to the two last columns of Table II for Onsager coefficients in the E frame of reference. The coefficients in the two columns are equivalent; they differ only by the factor that is used to convert between the dimensions used. The left column has dimensions m^{2}/s as they are directly derived from Eq. (B1). In order to obtain the coefficients with dimension mol^{2}/(J m s), which fit the flux Eq. (17), we have multiplied the left column with the factor *c*/*RT*. This unit gives the coefficient a dimension that fits with the thermodynamic description.

### D. Diffusion and charge transfer

#### 1. The mixed component scenario

In the E frame of reference, we find that the co-solvent DEC has the largest mobility of all components. All coefficients are positive, meaning that the movement of one constituent in the mixture always hampers the movement of others. This is intuitively obvious from a visual inspection of the simulation snapshot of the electrolyte in Fig. 1: The clustering tendency of ions and the structuring of DEC around the lithium-ions indicate that interactions obstruct rather than enhance the movements.

All Onsager coefficients are significant. Coupling coefficients have nearly the magnitude of the main coefficients (compare *L*^{+−} to *L*^{++} in Table II). The inequalities *L*_{ii}*L*_{jj} > *L*_{ij}*L*_{ji} and *L*_{ii} > 0 are anyway obeyed. The physical picture that this conveys is clear: One ion does not move without significantly impacting the movement of other ions and other components. All components move significantly relative to one another when one ion moves. This means that the models that assume the independent movement of ions, Eqs. (1)–(3), do not hold. They will, for instance, predict rather different electric conductivities when *L*^{+−} is significant. This is discussed further below.

#### 2. The neutral component scenario

The coefficients *L*_{D}*φ*, *L*_{L}*φ*, and *L*_{φφ} were computed with input from Table II and the stated Rules for Coupling of Fluxes, Eqs. (36)–(42). We find *L*_{L}*φ* = (−2.5 ± 0.8) × 10^{−6} mol C J^{−1} m^{−1} s^{−1}, *L*_{D}*φ* = (2.1 ± 0.8) × 10^{−6} mol C J^{−1} m^{−1} s^{−1}, and *L*_{φφ} = 0.23 ± 0.03 Ω^{−1} m^{−1}. The results are also shown in Table III. Diffusion coefficients, *ℓ*_{ij}, were computed from Eq. (20). The results for *ℓ*_{LL}, *ℓ*_{LD}, and *ℓ*_{DD} are also shown in Table III (bottom rows). The coefficients that describe coupling to charge transport can now be computed from this input. These results are shown in Table IV. The results have numerical accuracy, obtained from three independent simulations.

Consider the first column of Table IV, where no approximations are applied (all types of coupling are considered). Our simulated value is 0.23 ± 0.03 Ω^{−1} m^{−1} (Table IV). Lundgren *et al.*^{23} measured 0.7 Ω^{−1} m^{−1} for a similar electrolyte, while Morita *et al.*^{18} measured 0.8 Ω^{−1} m^{−1}. Newman^{57} measured a conductivity of 0.6 Ω^{−1} m^{−1} for LiPF_{6} in a propylene carbonate electrolyte. Trends in data or ratios are better captured in simulations than absolute values. With the Rules for Coupling of Fluxes available, it is easy to compare the models that we addressed in the Introduction. The neglect of cation–anion coupling has not much to say for the electric conductivity (cf. the second column in Table IV), but there is a large impact on the transference coefficient for co-solvent when this approximation is used. The use of self-diffusion coefficients to estimate ionic transport properties is slightly worse (see the last column of data).

The transport number of Li^{+} in the mixture was obtained from the transference coefficient of the salt.^{53} It is rather small (0.28) for DEC:EC = 1:1. A value less than unity implies that there is always a build-up of salt during battery discharge close to the anode, followed by diffusion. Neglecting the coupling coefficient *L*^{+−} changes the transport number of Li^{+} significantly, from 0.28 to 0.39. Clearly, a smaller value can lead to more concentration polarization at the electrode. The number may be essential for precise interpretations of electrolyte performance.

Using the self-diffusion coefficient rather than the coefficient that represents all interactions in the ternary mixture has a smaller impact on the transport number; see the last column of Table IV. The transport number changes to 0.37 in that case. Valøen and Reimers reported 0.38, Lundgren *et al.* 0.18, and Zugmann *et al.* 0.28, all values obtained from concentration cell experiments taking a solvent average frame of reference.^{15,23,60} Neither these nor other authors considered the effect of co-solvent electro-osmosis or an impact on the emf by solvent polarization.

A large positive transference coefficient of D equal to 0.90 is found. This has not been reported before. The positive value means that the co-solvent is moving in the opposite direction of the salt. The solvent transference number is zero only when the coupling coefficients are neglected. On average, 0.90 mol D accumulates on the cathode side when one mol of electric charge passes the electrolyte during discharge. Clearly, the model that assumes independent movement of ions cannot predict the full impact of charge transfer; the fact is that D is carried along when the battery is in operation. The transport of D will cause its own polarization, which has only recently been recognized in the literature as a cause of segregation.^{22}

The transport of EC with respect to the barycenter is opposite to that of DEC (cf. Table II). Large coefficients for solvents mean that the cell emf can have several significant contributions. Such contributions are presently neglected, say, in the formula for concentration polarization, which is used in the determination of transport numbers.

These points together underline once more that the assumption of independent movement of ions fails to describe component transport inside batteries. The neglected coefficients will lead to polarization, which eventually reduces cell efficiency. The assumption that the two solvents behave like one is also incorrect, and solvent segregation will reduce the cell voltage. These conclusions have been reached straightforwardly, using FDT and the Rules for Coupling of Fluxes.

#### 3. Relation to Fick’s law’s coefficients

Fick’s diffusion coefficients are used more frequently than the Onsager coefficients. They are often regarded as being more accessible than the cumbersome Onsager coefficients, so they deserve a special comment. Fick’s law’s coefficients for the present electrolyte, given by Eq. (29), are compared to the Onsager *ℓ*_{ij}-coefficients in Table V. We first note that these Onsager coefficients vary by one order of magnitude. Clearly, DEC is more mobile than salt. The self-diffusion coefficients (the four bottom rows of Table V) do not address inter-diffusion.

The set of Fick’s coefficients describes the same reality as the Onsager coefficients. However, their interpretation is not as straightforward as for the Onsager coefficients. They contain part of the driving force, i.e., the thermodynamic factor. As described in Appendix C, the values of these factors depend on the chosen ensemble conditions and on the concentration metric used (e.g., mole fractions or molar concentrations). This is a source of some ambiguity, especially when comparing different works. In Appendix C, we provide the thermodynamic factors for *NVT* and *NPT* conditions and for using either molar concentration or mole fraction. The different factor sets lead to significantly different diffusion coefficients. Such ambiguities can be avoided using Onsager coefficients. Our aim has been to make the Onsager coefficients more accessible and avoid the ambiguities that arise with Fick’s law.

#### 4. Coupling to heat transport

One of the Rules for Coupling of Fluxes, presented in this work, connected the coupling coefficients for charge and heat transport for the mixed and neutral component scenarios. We have not been able to exploit this rule fully yet. Fluctuation-dissipation theorems have been formulated for this case but have been used less. Clearly, this opens the road for future work. The need for better methods and models to deal with thermal transport phenomena in electrochemical cells is evident and has been pointed out.^{17} With the results of Table II, we can compute the value of *K* [see Eq. (53)] and find that the term *L*^{+−} × *K* is 20% of the value of *L*^{++}. In other words, the transport coefficients that describe thermal phenomena must also be taken into account to get a complete picture.

### E. Conclusions and perspectives

A set of rules that connect the Onsager transport coefficients of a mixed component scenario to a neutral component description were obtained here for a ternary electrolyte mixture using entropy production invariance. No assumption was involved, and the rules are general in this sense. All types of coupling between ion and solvent transport processes were systematically accounted for. The relations could be reduced to well-known expressions for the independent movement of ions, Eqs. (1)–(3), when the interaction between particles was set to zero. The rules allow for the systematic introduction of assumptions when information is lacking. A coupling coefficient cannot be neglected without also neglecting its reciprocal Onsager coefficient. The main Onsager coefficients are always positive. All Onsager coefficients obey symmetry criteria, which reflect the nature of molecular fluctuations. Fick’s law has its place, mostly for historical reasons.

We have obtained these results using a theoretical basis that takes advantage of NET for heterogeneous systems. Coupled transport phenomena in the example electrolyte were dealt with systematically. The procedure is applicable to electrochemical cells in general, and the Onsager coefficients play a central role. These coefficients have a solid foundation in the FDT and the entropy production of the system. Their symmetry property reflects the underlying molecular events, which should be brought out.

We have taken as an example the electrolyte of a typical lithium battery, containing a concentrated salt solution in two organic carbonate solvents, and presented in detail the rules that connect transport coefficients for a set of mixed variables, including ions, and a set of neutral variables that express coupled transport of heat, charge, and mass on a thermodynamic level. Solvent segregation, as well as concentration polarization, can now be described. The name “Rules for Coupling of Fluxes” has been used because the relations follow from entropy production invariance for the transformations in the (co-)solvent frame of reference.

## ACKNOWLEDGMENTS

The Research Council of Norway is acknowledged for Project No. 262644, PoreLab, by S.K., A.F.G., S.K.S., and A.L. S.K.S. acknowledges financial support from the NRC through Project No. 275754. Ø.G. acknowledges the Research Council of Norway for its support of the Norwegian Micro- and Nano-Fabrication Facility, NorFab, Project No. 295864. Dick Bedeaux is acknowledged for a critical review of the manuscript in its last stages and for explicit formulations of the fluctuation dissipation theorems.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Signe Kjelstrup**: Conceptualization (lead); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Supervision (equal); Validation (equal); Writing – original draft (lead); Writing – review & editing (equal). **Astrid Fagertun Gunnarshaug**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). **Øystein** **Gullbrekken**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). **Sondre K. Schnell**: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Supervision (equal); Validation (equal); Writing – review & editing (equal). **Anders Lervik**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (supporting); Validation (equal); Writing – original draft (supporting); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: FLUXES, FORCES, TRANSPORT COEFFICIENTS, AND THEIR DIMENSIONS

Table VI summarizes the thermodynamic transport coefficients for the ternary electrolyte, consisting of components L (the lithium salt), EC (ethylene carbonate), and DEC (diethyl carbonate), along with their dimensions.

### APPENDIX B: COEFFICIENTS FROM THE FLUCTUATION DISSIPATION THEOREMS

^{9}gave expressions for Onsager coefficients that are convenient for data reduction in molecular dynamics simulations at isothermal conditions. The Onsager coefficients Λ

_{ij}were expressed in terms of particle position vectors

**r**

_{k,i}, suitable for the simulation box,

*N*

_{i}and

*N*

_{j}are the particle numbers of components

*i*and

*j*, respectively, and

*N*is the total number of particles. The particles can be charged or neutral. The ⟨⋯⟩ brackets indicate an ensemble average, and the period ⋅ indicates a contraction of two vectors. The dimension of Λ

_{ij}is m

^{2}/s, and Onsager symmetry is obeyed Λ

_{ij}= Λ

_{ji}. The expression is particularly useful for the mixed component scenario, as particle positions are then well defined. Figure 4 displays the mean square displacements in Eq. (B1) from our simulations.

*fluctuating contributions*to the molar fluxes that constitute entropy production. These give

*δ*is the Dirac delta function.

**r**,

**J**

_{i,R}(

**r**,

*t*), is given by

*N*

_{A}is Avogadro’s number. By introducing Eq. (B4) into Eq. (B3), we obtain

**u**

_{i}in the E frame of reference, we obtain them for the

*n*− 1 components,

*n*− 1. By multiplication with

*c*/

*R*, we obtain the flux of component

*i*,

^{−2}s

^{−1}) relative to component E is then

*L*

^{ij}and Λ

_{ij}in the B- and E-frames of reference.

### APPENDIX C: THE THERMODYNAMIC FACTORS

In the main text, we used the thermodynamic factors (Γ_{ij}) to account for non-ideality when expressing the chemical potential gradients as concentration gradients for the Fickian diffusion picture. The present section will explain how we calculated the thermodynamic factors from molecular dynamics simulations.

*NVT*ensemble, where

^{34,35}

*c*

_{i}=

*N*

_{i}/

*V*is the concentration of species

*i*. This derivative is carried out at constant

*T*,

*V*, and concentrations of all components except

*i*(indicated by the subscript “

**c**

_{i}”), or equivalently, at constant particle numbers

*N*

_{i}for all components except

*i*(indicated by the subscript “

**N**

_{i}”).

^{33,34}related the thermodynamic factors to integrals of the radial distribution functions of the species. These Kirkwood–Buff integrals,

*G*

_{ij}, are

*g*

_{ij}is the radial distribution function for species

*i*and

*j*and they are related to the thermodynamic factors via

**B**is the matrix with elements

*B*

_{ij}=

*c*

_{i}

*δ*

_{ij}+

*c*

_{i}

*c*

_{j}

*G*

_{ij}and

*C*

_{ij}= (−1)

^{i+j}det(

**B**

_{ij}) is the

*ij*-cofactor in the cofactor expansion of the determinant of

**B**(

**B**

_{ij}is obtained from

**B**by deleting the

*i*th row and the

*j*th column). Since

*G*

_{ij}=

*G*

_{ji}, we have

**B**=

**B**

^{⊤}and

*C*

_{ij}=

*C*

_{ji}so that Γ

_{ij}

*c*

_{j}= Γ

_{ji}

*c*

_{i}. In practical applications, it might be more straightforward to calculate the thermodynamic factors directly from Γ

_{ij}=

*c*

_{i}

*A*

_{ij}, where

*A*

_{ij}are the elements of the matrix

**=**

*A*

*B*^{−1}. For a ternary system, we get

*NPT*ensemble and how we can introduce mole fractions. The transformation to

*NPT*conditions is

^{33,34}

*i*,

*κ*

_{T}denotes the isothermal compressibility,

*NPT*conditions (which are typically encountered in experimental settings). Analogous to Eq. (C1), we define a thermodynamic factor for

*NPT*conditions as

*a*

_{ij}=

*a*

_{ji}, which means that $\Gamma ijNPTNj=\Gamma jiNPTNi$, and from the Gibbs–Duhem relation at d

*T*= 0 and d

*P*= 0, we have

*j*. For a ternary system, with the following definitions:

*c*

_{1}

*a*

_{11}+

*c*

_{2}

*a*

_{21}+

*c*

_{3}

*a*

_{31}= 0.

^{61}

^{,}

*γ*

_{i}is the activity coefficient, and we use the superscript “

*x*” to indicate that $\Gamma ijx$ is based on the mole fraction. The subscript Σ indicates that the derivative is taken with the constraint that the mole fractions sum to one, and it is carried out at a constant mole fraction for all components except the last (the

*n*th). The required derivatives are

^{61–63}

^{,}

^{13,61,63}

*η*is the same as defined in Eq. (C12). We refer to Ruckenstein and Shulgin

^{63}for the expressions for $\Gamma 31x$ and $\Gamma 32x$.

_{ij}

*N*

_{j}= Γ

_{ji}

*N*

_{i}, and use the numerical values from Table VII to confirm the relations in Eq. (C11),

Thermodynamic factor . | NVT Eq. (C4)
. | NPT Eqs. (C10) and (C13)
. | NPT (mole fraction) Eq. (C16)
. |
---|---|---|---|

Γ_{LL} | 1.678 ± 0.004 | 1.43 ± 0.02 | 1.50 ± 0.02 |

Γ_{DD} | 47 ± 2 | 0.81 ± 0.03 | 1.20 ± 0.04 |

Γ_{EE} | 18.7 ± 0.5 | 0.40 ± 0.01 | ⋯ |

Γ_{LD} | 1.30 ± 0.02 | −0.30 ± 0.01 | −0.29 ± 0.02 |

Γ_{DL} | 5.8 ± 0.1 | −1.32 ± 0.04 | −0.98 ± 0.03 |

Γ_{LE} | 0.850 ± 0.008 | −0.018 ± 0.005 | ⋯ |

Γ_{EL} | 5.10 ± 0.05 | −0.11 ± 0.03 | ⋯ |

Γ_{DE} | 24.7 ± 0.7 | −0.38 ± 0.01 | ⋯ |

Γ_{ED} | 33 ± 1 | −0.51 ± 0.02 | ⋯ |

Thermodynamic factor . | NVT Eq. (C4)
. | NPT Eqs. (C10) and (C13)
. | NPT (mole fraction) Eq. (C16)
. |
---|---|---|---|

Γ_{LL} | 1.678 ± 0.004 | 1.43 ± 0.02 | 1.50 ± 0.02 |

Γ_{DD} | 47 ± 2 | 0.81 ± 0.03 | 1.20 ± 0.04 |

Γ_{EE} | 18.7 ± 0.5 | 0.40 ± 0.01 | ⋯ |

Γ_{LD} | 1.30 ± 0.02 | −0.30 ± 0.01 | −0.29 ± 0.02 |

Γ_{DL} | 5.8 ± 0.1 | −1.32 ± 0.04 | −0.98 ± 0.03 |

Γ_{LE} | 0.850 ± 0.008 | −0.018 ± 0.005 | ⋯ |

Γ_{EL} | 5.10 ± 0.05 | −0.11 ± 0.03 | ⋯ |

Γ_{DE} | 24.7 ± 0.7 | −0.38 ± 0.01 | ⋯ |

Γ_{ED} | 33 ± 1 | −0.51 ± 0.02 | ⋯ |

Γ_{ij}/Γ_{ji}
. | NVT
. | NPT
. | N_{i}/N_{j}
. |
---|---|---|---|

Γ_{LD}/Γ_{DL} | $1.305.80=0.224$ | $\u22120.30\u22121.32=0.227$ | $NL/ND=9204116=0.223$ |

Γ_{LE}/Γ_{EL} | $0.855.10=0.167$ | $\u22120.018\u22120.11=0.164$ | $NL/NE=9205520=0.167$ |

Γ_{DE}/Γ_{ED} | $24.733=0.748$ | $\u22120.38\u22120.51=0.745$ | $ND/NE=41165520=0.746$ |

Γ_{ij}/Γ_{ji}
. | NVT
. | NPT
. | N_{i}/N_{j}
. |
---|---|---|---|

Γ_{LD}/Γ_{DL} | $1.305.80=0.224$ | $\u22120.30\u22121.32=0.227$ | $NL/ND=9204116=0.223$ |

Γ_{LE}/Γ_{EL} | $0.855.10=0.167$ | $\u22120.018\u22120.11=0.164$ | $NL/NE=9205520=0.167$ |

Γ_{DE}/Γ_{ED} | $24.733=0.748$ | $\u22120.38\u22120.51=0.745$ | $ND/NE=41165520=0.746$ |

## REFERENCES

*Nonequilibrium Thermodynamics in Biophysics*

*Non-Equilibrium Thermodynamics for Engineers*

*CRC Handbook of Chemistry and Physics*

*Non-Equilibrium Thermodynamics of Heterogeneous Systems*

*Fluctuation, Relaxation and Resonance in Magnetic Systems*

^{7}Li NMR imaging

*Irreversible Thermodynamics: Theory and Applications*

_{6}in EC:DEC

*Nanothermodynamics: Theory and Applications*

*Understanding Molecular Simulation*

*Computer Simulation Using Particles*

*n*algorithm in LAMMPS

_{6}-based Li-ion battery electrolytes