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.

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.

Let us first recapitulate the main idea of the simplest model used to describe electrolyte transport with the example of a monovalent electrolyte. It is assumed that the two ions move independently of each other. The assumption means that the conductivity, κ, 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 m2 V−1 s−1). For a fully dissociated monovalent electrolyte in water, the conductivity is
κ=Fc(u++u).
(1)
Here 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
τ+=u+u++u,τ=uu++u,
(2)
where τ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=2RTFu+uu++u.
(3)
The dimension of D is m2/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.

In an extended model,2,6,14 the electric conductivity in Eq. (1) changes into
κ=F2(L++L+L++L),
(4)
where Lij 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 Lii = cui/F. Introducing self-diffusion coefficients, Di,self = uiRT/F, while still neglecting cation–anion interactions, the electric conductivity becomes
κ=F2(c+D+,self+cD,self)/RT.
(5)
The different formulas for κ 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.

FIG. 1.

Snapshot of the lithium battery electrolyte mixture consisting of the components shown in Fig. 2. The lithium-ions are colored yellow, the hexafluorophosphate ions are blue/green, the carbon atoms are gray, and the oxygen atoms are red. Hydrogen atoms are not shown. The carbonates shield lithium-ions, but a tendency to cluster positive and negative ions can also be seen.

FIG. 1.

Snapshot of the lithium battery electrolyte mixture consisting of the components shown in Fig. 2. The lithium-ions are colored yellow, the hexafluorophosphate ions are blue/green, the carbon atoms are gray, and the oxygen atoms are red. Hydrogen atoms are not shown. The carbonates shield lithium-ions, but a tendency to cluster positive and negative ions can also be seen.

Close modal

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. 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.

FIG. 2.

Constituents of a typical lithium-ion battery electrolyte. The yellow lithium-ion to the left has the blue–green hexafluorophosphate ion as a counterion. The solvents are ethylene carbonate (EC) and diethyl carbonate (DEC). The carbonate group occupies the center of DEC. It is more exposed in the EC.

FIG. 2.

Constituents of a typical lithium-ion battery electrolyte. The yellow lithium-ion to the left has the blue–green hexafluorophosphate ion as a counterion. The solvents are ethylene carbonate (EC) and diethyl carbonate (DEC). The carbonate group occupies the center of DEC. It is more exposed in the EC.

Close modal

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, 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.

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.

In the mixed component scenario, the entropy production, σ, is6,
σ=JqC1TJLi+B1Tμ̃Li+,TJPF6B1Tμ̃PF6,TJDB1TμD,TJEB1TμE,T.
(6)
Table VI in  Appendix A contains an overview of the quantities used, their symbols and their dimensions. The heat flux in this scenario is JqC (dimension J m−2 s−1), with conjugate force ∇(1/T). Component fluxes, indicated with subscripts Li+ and PF6, have as a driving force the negative electrochemical potential gradient divided by the temperature, i.e., μ̃Li+,T/T and μ̃PF6,T/T,2,3,21 where the electrochemical potential for ion i is defined by30,
μ̃iμi+ziFψ.
(7)
Here, zi 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.
In the neutral component scenario,3,21 the five independent contributions to σ are
σ=JqN1TJLB1TμL,TJDB1TμD,TJEB1TμE,Tj1Tφ.
(8)

The electric current density, j, (in A/m2) 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, JqN, 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 JqN is defined as the energy flux Jq minus the latent heat, Hi (partial molar enthalpy), carried by the three components: JqN=Jqj=13JjHj.12 The heat fluxes in the two scenarios were shown to be related.17 When lithium electrodes are used, we obtain JqC=JqN+TSLi+(j/F). The mass fluxes depend on the frame of reference.

Coefficients are generated with the barycentric frame of reference in simulation programs like the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS).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 m2 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 ui (in m/s) of all species i from the particle velocities vk,i. We have Niui(t)=k=1Nivk,i(t). In the barycentric frame of reference, the isothermal flux of i obeys
xiui=1RTj=1n=4Λijμ̃j,
(9)
where the summation is carried out over 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,
Ji=xicui=ciui=cRTj=1n=4Λijμ̃j=1Tj=1n=4Λijμ̃j.
(10)
The coefficients Λij have dimension m2 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,
Λij=cRΛij.
(11)
In this manner, we find the 4 × 4 = 16 coefficients in the matrix that derive from Eq. (6) for a constant temperature. Six of them are dependent through the Onsager relations, leaving ten unknown coefficients to describe the isothermal system. Fortunately, the number of independent coefficients can be further reduced because the driving forces are dependent through Gibbs–Duhem’s equation. We will reduce the coefficients from ten to six by choosing another frame of reference, and this reduction is our next step.
Gibbs–Duhem’s equation relates the chemical driving forces and offers the possibility of reducing the number of unknown coefficients from ten to six.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
μE,T=xLxEμ̃Li+,TxLxEμ̃PF6,TxDxEμD,T,
(12)
or
μE,T=xLxEμL,TxDxEμD,T,
(13)
where xi is the mole fraction of L or D. We introduce the respective expression in the corresponding expressions for the entropy production and obtain
σ=JqC1TJLi+1Tμ̃Li+,TJPF61Tμ̃PF6,TJD1TμD,T,
(14)
or
σ=JqN1TJL1TμL,TJD1TμD,Tj1Tφ.
(15)
All component fluxes are now in the E frame of reference (not indicated by a superscript). The number of flux-force products is reduced by one, and the flux-force relations have six independent coefficients. The entropy production is still the same, so a relation must exist between the Onsager coefficients in the E frame of reference and those in the B frame of reference.
The FDT is well defined in the mixed component scenario (see Sec. IV). Therefore, our primary target is this set of coefficients. The set is next transformed to the E frame of reference and to the operationally defined set. After some algebra, we find
Lij=cRΛijxjxnΛinxixnΛnj+xixjxn2Λnn.
(16)
The matrix of any set of Onsager coefficients is symmetric because we have used independent variables in each set. The first coefficient set obtained with the help of FDT, Λij, has the B-frame of reference. The set of Lij-coefficients with the E frame of reference is computed from these.
The Onsager coefficients are characteristic of the particular choice of variables. We proceed with the formulation that includes the heat flux to provide a start for thermal transport modeling. The ionic variable set obtains flux equations in the E frame of reference defined by the entropy production in Eq. (14),
JqC=Lqq1TLq+1Tμ̃Li+,TLq1Tμ̃PF6,TLqD1TμD,T,JLi+=L+q1TL++1Tμ̃Li+,TL+1Tμ̃PF6,TL+D1TμD,T,JPF6=Lq1TL+1Tμ̃Li+,TL1Tμ̃PF6,TLD1TμD,T,JD=LDq1TLD+1Tμ̃Li+,TLD1Tμ̃PF6,TLDD1TμD,T.
(17)

A coefficient superscript indicates the interactions among the fluxes in question; cf. Lij. The following dimensions apply to the Onsager coefficients in the mixed component scenario:

  • Lqq has dimension K J m−1 s−1, while the coupling coefficients Lq+, Lq, and LqD have dimension K mol m−1 s−1.

  • L++, L+−, L+D, L−−, L−D, and LDD have dimensions K mol2 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.

In the neutral component scenario, the entropy production of an electrolyte in a cell with lithium electrodes [Eq. (15)] prescribes the following flux-force matrix:
JqN=Lqq1TLqL1TμL,TLqD1TμD,TLqφ1Tφ,JL=LLq1TLLL1TμL,TLLD1TμD,TLLφ1Tφ,JD=LDq1TLDL1TμL,TLDD1TμD,TLDφ1Tφ,j=Lφq1TLφL1TμL,TLφD1TμD,TLφφ1Tφ.
(18)
The dimensions are
  • Lqq has dimension K J m−1 s−1.

  • LqL and LqD have dimensions K mol s−1 m−1.

  • L has dimension K C m−1 s−1.

The coefficient Lqq differs from Lqq because it refers to a different heat flux. The coefficients of the second and third rows have pairwise the same dimensions. The first pair LLq and LDq has dimension K mol m−1 s−1. The second and third pairs, LLL, LLD, and LLD, LDD, have dimensions K mol2 J−1 m−1 s−1. The last coefficients, LLφ and LDφ, 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.

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.

A basis for a definition of diffusion coefficients is obtained when we express ∇φ in terms of j with Eq. (18) and find
JqN=qq1TqL1TμL,TqD1TμD,T+LqφLφφj,JL=Lq1TLL1TμL,TLD1TμD,T+LLφLφφj,JD=Dq1TDL1TμL,TDD1TμD,T+LDφLφφj.
(19)
The uppercase coefficient symbols of Eq. (18) and the lowercase coefficient symbols of Eq. (19) are related by
ij=LijLiφLjφLφφ.
(20)
We use this relation in Sec. V C to obtain the diffusion coefficients ij. The coefficients are measured in the absence of j and ∇φ21 and with the E frame of reference.
The conductivity of the electrolyte is defined for a homogeneous, isothermal electrolyte and is obtained from the last row of Eq. (18) (Ohm’s law),
κ=jφT=0,μL=0,μD=0=LφφT.
(21)
The transference coefficients, ti, are defined by Hittorf experiments from Eq. (18),
ti=Jij/FT=0,μL=0,μD=0=FLiφLφφ.
(22)
The transference coefficient expresses the movement of the component due to electric current, measured in the limit of zero concentration gradients.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 
The Peltier coefficient, πN, is defined as the ratio of the heat flux and the electric current density,
πNJqNj/Fμi,T=0,T=0=LφqLφφ.
(23)
The single Peltier coefficient cannot be measured in isolation. It is accessible from the combination of coefficients, the Peltier 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
πNTT=0=ΔφΔTj=0.
The heat of transfer of components L and DEC, qi*, is best accessible via Onsager relations and the Soret coefficient of the cell; see the equation below. The heat of transfer can be computed from the gradient in chemical potential,
qL*=JqNJLJD=0,j=0,T=0=TμL,TTJL=JD=0,j=0.
(24)
An equivalent expression is obtained for component DEC. The temperature gradient, ∇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 ci and the corresponding thermodynamic factor, cf. Eq. (27).
With knowledge of the above-mentioned quantities, we can model or measure the electric potential gradient across the electrolyte of the cell,
φ=πTTtLμL,TtDμD,Tj/κ.
(25)
The flux-force relations in Eq. (18) can be used to find the profiles of the variables that enter and the single contributions to the cell potential.

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.

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.

The symmetric coefficient matrix can be related to Fick’s diffusion matrix by using the invariance of the entropy production again. We apply the conditions ∇T = 0, j = 0 to Eq. (19) and obtain
JL=LL1TμL,TLD1TμD,T,JD=DL1TμL,TDD1TμD,T.
(26)
Next, we express the gradients in the chemical potential with the gradients in concentration using12,
μi,T=j=12ΓjiRTcjcj,
(27)
where cj is the concentration of the (independent) components j = L and DEC. The non-ideality is expressed here using the thermodynamic factors Γij,33–35,
Γij=ciRTμjciT,V,ci,
(28)
where the subscript ci 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
JL=RcLcL(LLΓLL+LDΓLD)RcDcD(LDΓDD+LLΓDL),JD=RcLcL(DLΓLL+DDΓLD)RcDcD(DLΓDL+DDΓDD).
(29)
The factor 1/T from the force disappears from the coefficients of Eq. (22) when the concentration gradient is introduced. However, the introduction of the gradients in ci and cj 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 Dij,
DLL=R(LLΓLL+LDΓLD)/cL,DLD=R(LLΓDL+LDΓDD)/cD,DDL=R(DLΓLL+DDΓLD)/cL,DDD=R(DLΓDL+DDΓDD)/cD.
(30)
The equations can be used to compute Fick’s law’s diffusion coefficients. The coefficients follow when the concentration gradients are used as driving forces. They need not obey the conditions for the determinant and the main coefficients: that liiljjlijlji ≥ 0, lii > 0. Knowledge of Γij/ci and ij gives the matrix elements Dij. 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.

As pointed out earlier, the descriptions of the mixed and neutral component scenarios are equivalent, and relations exist due to 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:
j/F=JLi+JPF6.
(31)
For lithium reversible electrodes, the observable salt flux is identified by the movement of the anion,
JL=JPF6,
(32)
and the gradient in the measured electric potential is identified by the electrochemical potential gradient of the lithium-ion,
Fφ=μ̃Li+.
(33)
The gradient in the chemical potential of the salt is, furthermore, the sum of the electrochemical potential gradients of the ions,
μL=μ̃Li++μ̃PF6.
(34)
The last identity follows from the definition in Eq. (7).
We introduce the expressions for the ionic fluxes, Eq. (17), into the right-hand side of the equation for the electric current density. Equation (31) gives
j/F=(L+qLq)1T(L++L+)1Tμ̃Li+,T(L+L)1Tμ̃PF6,T(L+DLD)1TμD,T.
(35)
We introduce the expressions for forces in Eqs. (33) and (34) and compare the outcome to the bottom line of Eq. (18). The following relations are then obtained between the coefficients in the mixed (left-hand side) and neutral (right-hand side) component scenarios:
Lφφ=(L++L+L++L)F2,
(36)
LφL=(L+L)F=LLφ
(37)
LφD=(L+DLD)F=LDφ,
(38)
Lφq=(L+qLq)F=Lqφ.
(39)
Equations (36)(39) describe coupling to charge transport, indicated by the subscript φ. For mixtures of multivalent electrolytes, the valency of the respective ion enters. For mass transport, the corresponding coefficient relations are
LLL=L,
(40)
LLD=LD=LDL,
(41)
LDD=LDD.
(42)
Onsager relations have been used. The coefficient relations do not depend on the value of the driving forces. The relations follow from constant entropy production only and appear without assumptions. Approximations can be introduced, but not for one coefficient alone; the whole set must be considered. We have, therefore, given the set of equations the name “Rules for Coupling of Fluxes.” They were derived here for a lithium battery-related case, but similar forms can be found for similar systems. However, they are rules that apply in the solvent E frame of reference.
We can now use the Rules to generalize the formulas for the transport properties that appeared in Eqs. (1)(3). The electric conductivity, κ, prescribed by the rules is the same. We repeat
κ=LφφT=F2(L++2L++L).
(43)
The coupling coefficient L+− can be positive or negative in sign. Coupling may, therefore, enhance or reduce κ.
Component transference coefficients follow. For the salt L,
tL=FLφLLφφ=L+LL++2L++L.
(44)
The transference coefficient is not equal to the transport number, and the transport number of ion 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. 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
tD=FLφDLφφ=L+DLDL++2L++L.
(45)
The coefficient describes the transport of DEC that accompanies charge transport through the electrolyte. Depending on the sign of the coupling coefficients, it will lead to an accumulation of DEC to the right or left of the system. We see that it is zero when the co-solvent interacts equally with the cation and with the anion. In the present electrolyte, these coefficients are unknown. We know from ion-exchange membranes that the number of water molecules transported by cations or anions is high.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.
The first contribution to the main diffusion coefficient for the salt in Eq. (30) can be found in terms of Onsager coefficients,6 
DL=RcΓLLLL=LLLLφL2Lφφ1TμLcL=L(L+L)2L++2L++L1TμLcL.
(46)
This is a common expression, but for binary mixtures. More general expressions can be developed using Eq. (30).

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 JqN is measurable in the neutral component scenario. From the discussion below Eq. (8), we understand that JqC 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.

The heat of transfer is defined in the mixed component scenario,
qLi+*=JqCJLi+JD=0,JPF6=0,T=0=Tμ̃Li+,TTJD=0,JPF6=0,JLi+=0,
(47)
qPF6*=JqCJPF6JD=0,JLi+=0,T=0=Tμ̃PF6,TTJD=0,JPF6=0,JLi+=0.
(48)
From Eqs. (34) and (24), we obtain
qL*=qLi+*+qPF6*.
(49)
In the mixed component scenario, the conditions for measurement are JD = 0 and JPF6=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
JqC=Lq+LqK+LqDM1Tμ̃Li+.
(50)
The cation flux that corresponds to these conditions is
JLi+=L++L+K+L+DM1Tμ̃Li+.
(51)
It is then possible to define the ratio of these fluxes as the heat of transfer of the lithium-ion in the mixed component scenario,
qLi+*=JqCJLi+JD=0,JPF6=0,T=0=Lq+LqK+LqDML++L+K+L+DM.
(52)
Symbols K and M are lumped coefficients, which are functions of the Onsager coefficients,
K=L+LDLD+/LDDLLDLD/LDD,
(53)
M=LD+LDD+LDLDDK.
(54)
When the assumption of independent movement of ions applies, the complicated expressions reduce. We introduce 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
qLi+*=JqCJLi+JD=0,JPF6=0,T=0=Lq+L++,
(55)
and
qPF6*=JqCJPF6JD=0,JLi+=0,T=0=LqL.
(56)
The transport numbers simplify accordingly. One of the Rules for Coupling of Fluxes, Lφq = Lq+Lq, can now be applied. We divide both sides by Lφφ (L++ + L−−) and obtain
πN=LφqLφφ=Lq+L+++LLqL+++L=τLi+qLi+*τPF6qPF6*.
(57)

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

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:

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

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

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

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

FIG. 3.

Procedure for calculation of diffusion coefficients and electric- and thermal conductivities. Arrows show the sequence of steps. The procedure starts with the calculation of Onsager coefficients for the mixed component scenario using the fluctuation dissipation theorems (FDTs orange block). The coefficients for the mixed component scenario (gray block) are next obtained for the appropriate frame of reference. Using the Rules for Coupling of Fluxes, we next find properties that describe coupled transport of heat, mass, and charge (central block and left bottom corner). These steps will enable us to eventually model electrochemical cells in terms of nonequilibrium thermodynamics (NET).

FIG. 3.

Procedure for calculation of diffusion coefficients and electric- and thermal conductivities. Arrows show the sequence of steps. The procedure starts with the calculation of Onsager coefficients for the mixed component scenario using the fluctuation dissipation theorems (FDTs orange block). The coefficients for the mixed component scenario (gray block) are next obtained for the appropriate frame of reference. Using the Rules for Coupling of Fluxes, we next find properties that describe coupled transport of heat, mass, and charge (central block and left bottom corner). These steps will enable us to eventually model electrochemical cells in terms of nonequilibrium thermodynamics (NET).

Close modal

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.

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.

The equilibrium MD simulations were performed using the LAMMPS31 software. Atomic and molecular interactions were described by the OPLS-AA36 potential. The parameters for the solvent molecules were obtained from the LigParGen web service,37–39 and the parameters for Li+ and PF6 were taken from Jensen and Jorgensen40 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 solver44 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, 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 barostat48–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 Å.

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 code52 by tracking the motion of the central atom of the various components.

We computed the coefficients Lij 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 z2, 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 integrals33 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 Vegt54 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 

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.

TABLE I.

Thermodynamic factors (based on mole fractions) of the lithium battery electrolyte from molecular dynamics simulations (at 300 K and an average pressure of 0.4 bar). The thermodynamic factors are calculated from the radial distribution functions, averaged over three independent simulations. The error estimates were obtained as the standard deviation of the thermodynamic factors calculated directly from the radial distribution functions (without averaging them).

ΓLLΓDDΓLDΓDL
1.50 ± 0.02 1.20 ± 0.04 −0.29 ± 0.02 −0.98 ± 0.03 
ΓLLΓDDΓLDΓDL
1.50 ± 0.02 1.20 ± 0.04 −0.29 ± 0.02 −0.98 ± 0.03 

Tables IIV 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

TABLE II.

Onsager diffusion coefficients for the mixed component scenario of the isothermal electrolyte obtained by molecular dynamics simulations (at 300 K and an average pressure of 0.4 bar). The coefficients are calculated using Eq. (B1) in the barycentric (B) frame of reference and converted to the solvent (E) frame of reference. Coefficients refer to Eq. (17) and are given in two sets of units for the E frame of reference (the total molar concentration, c, is used to convert the units). The values are averaged over three independent simulations, and the standard deviations are taken as the error estimate here.

Frame of reference
CoefficientB × 10−11 m2/sE × 10−11 m2/sE × c/RT 10−11 mol2/(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 
LD+ 0.5 ± 0.1 2.4 ± 0.2 12.1 ± 1.1 
LD− −0.3 ± 0.1 1.8 ± 0.2 9.2 ± 0.9 
LDD 2.8 ± 0.2 11.3 ± 1.3 55.9 ± 6.2 
LE+ −0.9 ± 0.1   
LE− −1.0 ± 0.1   
LEE 6.1 ± 0.9   
LED −3.4 ± 0.4   
Frame of reference
CoefficientB × 10−11 m2/sE × 10−11 m2/sE × c/RT 10−11 mol2/(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 
LD+ 0.5 ± 0.1 2.4 ± 0.2 12.1 ± 1.1 
LD− −0.3 ± 0.1 1.8 ± 0.2 9.2 ± 0.9 
LDD 2.8 ± 0.2 11.3 ± 1.3 55.9 ± 6.2 
LE+ −0.9 ± 0.1   
LE− −1.0 ± 0.1   
LEE 6.1 ± 0.9   
LED −3.4 ± 0.4   
TABLE III.

Onsager diffusion coefficients for the neutral component scenario of the isothermal electrolyte. The top six coefficients are computed from the Rules for Coupling of Fluxes using the indicated equation and values in Table II. The bottom three values are obtained from Eq. (20). Conditions are otherwise the same as for Table II.

CoefficientEquationValueUnit
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) 
LLL (40) 6.2 ± 1.0 10−11 mol2/(J m s) 
LDL (41) 9.2 ± 0.9 10−11 mol2/(J m s) 
LDD (42) 55.9 ± 6.2 10−11 mol2/(J m s) 
LL (20) 3.7 ± 0.2 10−11 mol2/(J m s) 
DL (20) 11.3 ± 0.6 10−11 mol2/(J m s) 
DD (20) 53.7 ± 6.4 10−11 mol2/(J m s) 
CoefficientEquationValueUnit
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) 
LLL (40) 6.2 ± 1.0 10−11 mol2/(J m s) 
LDL (41) 9.2 ± 0.9 10−11 mol2/(J m s) 
LDD (42) 55.9 ± 6.2 10−11 mol2/(J m s) 
LL (20) 3.7 ± 0.2 10−11 mol2/(J m s) 
DL (20) 11.3 ± 0.6 10−11 mol2/(J m s) 
DD (20) 53.7 ± 6.4 10−11 mol2/(J m s) 
TABLE IV.

Transport coefficients computed from the Rules for Coupling of Fluxes applied to the isothermal ternary battery electrolyte. The effect of neglecting cation–anion coupling, L+− = L−+ = 0, is shown in the third column. The effect of using self-diffusion coefficients alone is shown in the fourth column. The results are computed from the values in Table III. Symbols are defined in the text in connection with the equations mentioned. Conditions are otherwise the same as for Table II.

CoefficientAll types of couplingNo cation–anion couplingSelf-diffusion coefficients
κ−1 m−1 0.23 ± 0.03 0.27 ± 0.03 0.41 ± 0.01 
tL −0.97 ± 0.12 −0.81 ± 0.07 −0.84 ± 0.01 
tD 0.90 ± 0.46 0.39 ± 0.17 
τ+ 0.28 ± 0.09 0.39 ± 0.05 0.37 ± 0.01 
τ 0.72 ± 0.09 0.61 ± 0.05 0.63 ± 0.01 
CoefficientAll types of couplingNo cation–anion couplingSelf-diffusion coefficients
κ−1 m−1 0.23 ± 0.03 0.27 ± 0.03 0.41 ± 0.01 
tL −0.97 ± 0.12 −0.81 ± 0.07 −0.84 ± 0.01 
tD 0.90 ± 0.46 0.39 ± 0.17 
τ+ 0.28 ± 0.09 0.39 ± 0.05 0.37 ± 0.01 
τ 0.72 ± 0.09 0.61 ± 0.05 0.63 ± 0.01 
TABLE V.

Onsager’s ij’s and Fick’s interdiffusion coefficients, Dij, and self-diffusion coefficients, Di. The electrolyte conditions were described in Table II. The lithium salt L and its ions, and the co-solvent DEC, diffuse in a co-solvent of E.

CoefficientEquationValue × 10−11 m2 s−1Value × 1011 mol2 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 
DLL (29) 52.3 ± 2.6  
DLD (29) 309.4 ± 16.9  
DDL (29) 222.8 ± 22.4  
DDD (29) 1453.5 ± 170.2  
DLi+  7.2 ± 0.2  
DPF6  12.4 ± 0.2  
DD  11.6 ± 0.1  
DE  22.1 ± 0.1  
CoefficientEquationValue × 10−11 m2 s−1Value × 1011 mol2 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 
DLL (29) 52.3 ± 2.6  
DLD (29) 309.4 ± 16.9  
DDL (29) 222.8 ± 22.4  
DDD (29) 1453.5 ± 170.2  
DLi+  7.2 ± 0.2  
DPF6  12.4 ± 0.2  
DD  11.6 ± 0.1  
DE  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 IIV, 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).

TABLE VI.

Thermodynamic transport coefficients of the ternary electrolyte of components L, EC, and DEC and their dimensions according to the entropy production. Component EC serves as the frame of reference. Common for isothermal systems is to divide by the temperature T. Coefficient units are then multiplied by 1/K.

NameSymbolDimension according to Eqs. (17) or (19)
Entropy production σ J K−1 m−3 s−1 
Heat flux Jq J m−2 s−1 
Thermal force ∇(1/TK−1 m−1 
Mass flux Ji 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 Lqq K J m−1 s−1 
Coupling coefficient, heat L LqL = LLq K mol m−1 s−1 
Coupling coefficient, heat D LqD = LDq K mol m−1 s−1 
Electric conductivity Lφφ K C V−1 m−1 s−1 
Coupling coefficient, charge-heat Lφq = L K C m−1 s−1 
Coupling coefficient, charge-mass L LφL = LLφ K C mol J−1 m−1 s−1 
Coupling coefficient, charge-mass D LφD = LDφ K C mol J−1 m−1 s−1 
Diffusion main coefficient L LL K mol2 J−1 m−1 s−1 
Diffusion main coefficient D DD K mol2 J−1 m−1 s−1 
Diffusion coupling coefficient LD LD = LD K mol2 J−1 m−1 s−1 
Diffusion main coefficient ++ L++ K mol2 J−1 m−1 s−1 
Diffusion main coefficient −− L−− K mol2 J−1 m−1 s−1 
Diffusion coupling coefficient +− L+− = L−+ K mol2 J−1 m−1 s−1 
Diffusion coupling coefficient D+ LD+ = L+D K mol2 J−1 m−1 s−1 
Diffusion coupling coefficient D− LD− = L−D K mol2 J−1 m−1 s−1 
Diffusion main coefficient DD LDD K mol2 J−1 m−1 s−1 
NameSymbolDimension according to Eqs. (17) or (19)
Entropy production σ J K−1 m−3 s−1 
Heat flux Jq J m−2 s−1 
Thermal force ∇(1/TK−1 m−1 
Mass flux Ji 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 Lqq K J m−1 s−1 
Coupling coefficient, heat L LqL = LLq K mol m−1 s−1 
Coupling coefficient, heat D LqD = LDq K mol m−1 s−1 
Electric conductivity Lφφ K C V−1 m−1 s−1 
Coupling coefficient, charge-heat Lφq = L K C m−1 s−1 
Coupling coefficient, charge-mass L LφL = LLφ K C mol J−1 m−1 s−1 
Coupling coefficient, charge-mass D LφD = LDφ K C mol J−1 m−1 s−1 
Diffusion main coefficient L LL K mol2 J−1 m−1 s−1 
Diffusion main coefficient D DD K mol2 J−1 m−1 s−1 
Diffusion coupling coefficient LD LD = LD K mol2 J−1 m−1 s−1 
Diffusion main coefficient ++ L++ K mol2 J−1 m−1 s−1 
Diffusion main coefficient −− L−− K mol2 J−1 m−1 s−1 
Diffusion coupling coefficient +− L+− = L−+ K mol2 J−1 m−1 s−1 
Diffusion coupling coefficient D+ LD+ = L+D K mol2 J−1 m−1 s−1 
Diffusion coupling coefficient D− LD− = L−D K mol2 J−1 m−1 s−1 
Diffusion main coefficient DD LDD K mol2 J−1 m−1 s−1 

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 choice17,59).

The results presented in Table II have the same order of magnitude as those published in the literature, 10−10 m2 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 m2/s as they are directly derived from Eq. (B1). In order to obtain the coefficients with dimension mol2/(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.

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 LiiLjj > LijLji and Lii > 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 LDφ, LLφ, and Lφφ were computed with input from Table II and the stated Rules for Coupling of Fluxes, Eqs. (36)(42). We find LLφ = (−2.5 ± 0.8) × 10−6 mol C J−1 m−1 s−1, LDφ = (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. Newman57 measured a conductivity of 0.6 Ω−1 m−1 for LiPF6 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.

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.

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.

The authors have no conflicts to disclose.

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).

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

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.

Krishna and van Baten9 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 rk,i, suitable for the simulation box,
Λij=16NlimΔt1Δtk=1Nirk,i(t+Δt)rk,i(t)l=1Njrl,j(t+Δt)rl,j(t).
(B1)
The expression is called the Einstein approximation. Here, Ni and Nj 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.
FIG. 4.

Mean squared displacement (MSD) of the various molecular and ionic correlations as a function of time difference. The slopes of the curves are the Onsager diffusion coefficients in Table II. The region used for linear regression is chosen as the data points that are closest to a slope of 1 in the log–log plot. At least four data points are used in the linear regression. The time intervals used for linear regression will vary due to the different correlations between the species.

FIG. 4.

Mean squared displacement (MSD) of the various molecular and ionic correlations as a function of time difference. The slopes of the curves are the Onsager diffusion coefficients in Table II. The region used for linear regression is chosen as the data points that are closest to a slope of 1 in the log–log plot. At least four data points are used in the linear regression. The time intervals used for linear regression will vary due to the different correlations between the species.

Close modal
Equation (B1) can be derived from fluctuation-dissipation theorems. Consider isothermal conditions. The start is the entropy production of the isothermal system, obtained from Eq. (17) as σ=1Tj=1nJiμ̃j. The fluctuation-dissipation theorem applies to the fluctuating contributions to the molar fluxes that constitute entropy production. These give
Ji,R(r,t)Jj,R(r,t)=2RLijδ(rr)δ(tt).
(B2)
A bold symbol means a vector. The symbol δ is the Dirac delta function.
The displacements in Eq. (B1) are first written as integrals over the velocities. We find
Λij=16NlimΔt1Δttt+Δtdttt+Δtdt×k=1Nivk,i(t)l=1Njvl,j(t).
(B3)
The fluctuating contributions to the velocities are the same in and away from equilibrium. The fluctuating contribution to the molar flux at position r, Ji,R(r, t), is given by
NAJi,R(r,t)=k=1Nivk,i(t)δ(rk,i(t)r),
(B4)
where NA is Avogadro’s number. By introducing Eq. (B4) into Eq. (B3), we obtain
Λij=NA26NlimΔt1Δttt+Δtdttt+ΔtdtVdrVdr×Ji,R(r,t)Jj,R(r,t).
(B5)
We next introduce Eq. (B2) in this formula and find that
Λij=RcΛij.
(B6)
This result is exactly the same that we find from the linear laws, see Eq. (11). We, therefore, confirm agreement with Eq. (B2), the fluctuation-dissipation theorem.
For the component velocities ui in the E frame of reference, we obtain them for the n − 1 components,
xiui=1RTj=1n1ΛijxjxnΛinμ̃j.
(B7)
The summation is over n − 1. By multiplication with c/R, we obtain the flux of component i,
Ji=xicui=ciui=cRTj=1n1ΛijxjxnΛinμ̃j.
(B8)
This flux (in mol m−2 s−1) relative to component E is then
JiE=ci(uiun)=cRTj=1n1ΛijxjxnΛinxixnΛnj+xixjxn2Λnnμ̃j=1Tj=1n1LijEμ̃j   for   i=1,..,n1.
(B9)
These equations give the relation between Lij and Λij in the B- and E-frames of reference.

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.

The thermodynamic factor is particular to the ensemble in question. In the present work, equilibrium molecular dynamics simulations are performed in the NVT ensemble, where34,35
ΓijcikBTμjciT,V,ci=NikBTμjNiT,V,Ni,
(C1)
and ci = Ni/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 “ci”), or equivalently, at constant particle numbers Ni for all components except i (indicated by the subscript “Ni”).
Kirkwood and Buff33,34 related the thermodynamic factors to integrals of the radial distribution functions of the species. These Kirkwood–Buff integrals, Gij, are
Gij=4π0gij(r)1r2dr,
(C2)
where gij is the radial distribution function for species i and j and they are related to the thermodynamic factors via
1kBTμicjT,V,cj=Cijdet(B)=Γjicj,
(C3)
where B is the matrix with elements Bij = ciδij + cicjGij and Cij = (−1)i+j det(Bij) is the ij-cofactor in the cofactor expansion of the determinant of B (Bij is obtained from B by deleting the ith row and the jth column). Since Gij = Gji, we have B = B and Cij = Cji so that Γijcj = Γjici. In practical applications, it might be more straightforward to calculate the thermodynamic factors directly from Γij = ciAij, where Aij are the elements of the matrix A = B−1. For a ternary system, we get
Γ11ζ=1+G22c2+G33c3+c2c3G22G33G23G32,Γ12ζ=c1G12+c3G12G33+G13G32,Γ13ζ=c1G13+c2G12G23G13G22,Γ22ζ=1+G11c1+G33c3+c1c3G11G33G13G31,Γ23ζ=c2G23+c1G11G23+G13G21,Γ33ζ=1+G11c1+G22c2+c1c2G11G22G12G21,
(C4)
where
ζ=det(B)c1c2c3.
(C5)
For completeness, we will detail how we can adapt the thermodynamic factors to the isothermal-isobaric NPT ensemble and how we can introduce mole fractions. The transformation to NPT conditions is33,34
aji=μiNjT,P,Nj=μiNjT,V,NjV̄iV̄jVκT,
(C6)
where V̄i represents the partial molar volume of species i,
V̄i=kckCkii,jcicjCij,
(C7)
and κT denotes the isothermal compressibility,
κT=det(B)kBTi,jcicjCij.
(C8)
Together, these equations give
aji=μiNjT,P,Nj=kBTVdet(B)CijkckCkikckCkji,jcicjCij.
(C9)
This equation can be employed to compute thermodynamic factors under NPT conditions (which are typically encountered in experimental settings). Analogous to Eq. (C1), we define a thermodynamic factor for NPT conditions as
ΓijNPT=NikBTμjNiT,P,Ni=NikBTaij.
(C10)
From Eq. (C6), we have that aij = aji, which means that ΓijNPTNj=ΓjiNPTNi, and from the Gibbs–Duhem relation at dT = 0 and dP = 0, we have
i=1nciaij=0,
(C11)
for all j. For a ternary system, with the following definitions:
η=1c1c2c3ijcicjCij,Δij=Gii+Gjj2Gij,
(C12)
one may show that
a11=kBTVηc2+c3+c2c3Δ23c1,a22=kBTVηc1+c3+c1c3Δ13c2,a33=kBTVηc1+c2+c1c2Δ12c3,a12=a21=kBTVη1+c3G33+G12G13G23,a13=a31=kBTVη1+c2G22+G13G12G23,a23=a32=kBTVη1+c1G11+G23G12G13.
(C13)
The relations in Eq. (C11) can be verified from these equations; for instance, c1a11 + c2a21 + c3a31 = 0.
The thermodynamic factors are also commonly expressed as derivatives with respect to the mole fractions,61,
Γijx=xikBTμixjT,P,Σ=δij+xilnγixjT,P,Σ,
(C14)
where γi is the activity coefficient, and we use the superscript “x” to indicate that Γ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 nth). The required derivatives are61–63,
μixjT,P,Σ=NμiNjT,P,NjNμiNnT,P,Nn,
(C15)
and with Eqs. (C15) and (C14), explicit expressions for Γijx are13,61,63
Γ11x=1ηc1+c2+c3+c1c2G12+G13+G22G23+c2c3Δ23,Γ12x=c1ηc2G12+G13+G22G23+c3G12+G13+G23G33,Γ21x=c2ηc1G11G12G13+G23+c3G13+G23G12G33,Γ22x=1ηc1+c2+c3+c1c2G11G12G13+G23+c1c3Δ13,
(C16)
where η is the same as defined in Eq. (C12). We refer to Ruckenstein and Shulgin63 for the expressions for Γ31x and Γ32x.
The calculated thermodynamic factors from our simulations are given in Table VII. The salt is treated as one component, “L.” In Table VIII, we examine the symmetry relations, ΓijNj = ΓjiNi, and use the numerical values from Table VII to confirm the relations in Eq. (C11),
ΓLLNPT+ΓDLNPT+ΓELNPT=1.431.320.110,ΓLDNPT+ΓDDNPT+ΓEDNPT=0.30+0.810.510,ΓLENPT+ΓDENPT+ΓEENPT=0.0180.38+0.400.002.
(C17)
Regarding the last check, we note that the value of ΓLENPT is on the same order of magnitude as the uncertainties in ΓDENPT and ΓEENPT.
TABLE VII.

Thermodynamic factors of the lithium battery electrolyte from molecular dynamics simulations (at 300 K and an average pressure of 0.4 bar). The thermodynamic factors are calculated from the radial distribution functions, averaged over three independent simulations. The error estimates were obtained as the standard deviation of the thermodynamic factors calculated directly from the radial distribution functions (without averaging them). The superscript “NPT” from Eq. (C10) and the superscript “x” from Eq. (C14) are excluded here.

Thermodynamic factorNVT 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 factorNVT 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 ⋯ 
TABLE VIII.

Consistency check for the thermodynamic factors, ΓijNj = ΓjiNi, numerical values from Table VII.

ΓijjiNVTNPTNi/Nj
ΓLDDL 1.305.80=0.224 0.301.32=0.227 NL/ND=9204116=0.223 
ΓLEEL 0.855.10=0.167 0.0180.11=0.164 NL/NE=9205520=0.167 
ΓDEED 24.733=0.748 0.380.51=0.745 ND/NE=41165520=0.746 
ΓijjiNVTNPTNi/Nj
ΓLDDL 1.305.80=0.224 0.301.32=0.227 NL/ND=9204116=0.223 
ΓLEEL 0.855.10=0.167 0.0180.11=0.164 NL/NE=9205520=0.167 
ΓDEED 24.733=0.748 0.380.51=0.745 ND/NE=41165520=0.746 
1.
J. S.
Newman
,
Electrochemical Systems
,
2nd ed.
(
Prentice Hall
,
1991
).
2.
R.
Haase
,
Thermodynamics of Irreversible Processes
(
Addison-Wesley
,
1969
).
3.
A.
Katchalsky
and
P.
Curran
,
Nonequilibrium Thermodynamics in Biophysics
,
2nd ed.
(
Harvard University Press
,
1975
).
4.
S.
Kjelstrup
,
D.
Bedeaux
,
E.
Johannessen
, and
J.
Gross
,
Non-Equilibrium Thermodynamics for Engineers
,
2nd ed.
(
World Scientific
,
2018
).
5.
CRC Handbook of Chemistry and Physics
,
96th ed.
, edited by
W.
Haynes
(
CRC Press
,
Boca Raton
,
2015
).
6.
S.
Kjelstrup
and
D.
Bedeaux
, in
Non-Equilibrium Thermodynamics of Heterogeneous Systems
,
1st ed.,
Series on Advances in Statistical Mechanics (
World Scientific
,
2008
), Vol.
16
.
7.
R.
Kubo
, “
The fluctuation-dissipation theorem
,”
Rep. Prog. Phys.
29
,
255
(
1966
).
8.
H.
Callen
, “
The fluctuation-dissipation theorems and irreversible thermodynamics
,” in
Fluctuation, Relaxation and Resonance in Magnetic Systems
, edited by
D.
ter Haar
(
Scottish Universities’ Summer School, Oliver and Boyd
,
Edinburgh
,
1962
), pp.
15
23
.
9.
R.
Krishna
and
J. M.
van Baten
, “
The Darken relation for multicomponent diffusion in liquid mixtures of linear alkanes: An investigation using molecular dynamics (MD) simulations
,”
Ind. Eng. Chem. Res.
44
,
6939
6947
(
2005
).
10.
X.
Liu
,
T. J. H.
Vlugt
, and
A.
Bardow
, “
Predictive Darken equation for Maxwell-Stefan diffusivities in multicomponent mixtures
,”
Ind. Eng. Chem. Res.
50
,
10350
10358
(
2011
).
11.
K. D.
Fong
,
H. K.
Bergström
,
B. D.
McClosky
, and
K. K.
Mandadapu
, “
Transport phenomena in electrolyte solutions: Nonequilibrium thermodynamics and statistical mechanics
,”
AIChE J.
66
,
e17091
(
2020
).
12.
S.
de Groot
and
P.
Mazur
,
Non-Equilibrium Thermodynamics
(
Dover
,
London
,
1984
).
13.
X.
Liu
,
A.
Martín-Calvo
,
E.
McGarrity
,
S. K.
Schnell
,
S.
Calero
,
J.-M.
Simon
,
D.
Bedeaux
,
S.
Kjelstrup
,
A.
Bardow
, and
T. J. H.
Vlugt
, “
Fick diffusion coefficients in ternary liquid systems from equilibrium molecular dynamics simulations
,”
Ind. Eng. Chem. Res.
51
,
10247
10258
(
2012
).
14.
J.
Newman
and
T. W.
Chapman
, “
Restricted diffusion in binary solutions
,”
AIChE J.
19
,
343
348
(
1973
).
15.
S.
Zugmann
,
M.
Fleischmann
,
M.
Amereller
,
R. M.
Gschwind
,
H. D.
Wiemhöfer
, and
H. J.
Gores
, “
Measurement of transference numbers for lithium ion electrolytes via four different methods, a comparative study
,”
Electrochim. Acta
56
,
3926
3933
(
2011
).
16.
M.
Klett
,
M.
Giesecke
,
A.
Nyman
,
F.
Hallberg
,
R. W.
Lindström
,
G.
Lindbergh
, and
I.
Furó
, “
Quantifying mass transport during polarization in a Li ion battery electrolyte by in situ 7Li NMR imaging
,”
J. Am. Chem. Soc.
134
,
14654
14657
(
2012
).
17.
S.
Kjelstrup
,
K. R.
Kristiansen
,
A. F.
Gunnarshaug
, and
D.
Bedeaux
, “
Seebeck, Peltier and Soret effects: On different formalisms for transport equations in thermogalvanic cells
,”
J. Chem. Phys.
158
,
020901
(
2023
).
18.
M.
Morita
,
Y.
Asai
,
N.
Yoshimoto
, and
M.
Ishikawa
, “
A Raman spectroscopic study of organic electrolyte solutions based on binary solvent systems of ethylene carbonate with low viscosity solvents which dissolve different lithium salts
,”
J. Chem. Soc., Faraday Trans.
94
,
3451
3456
(
1998
).
19.
D. M.
Seo
,
S.
Reininger
,
M.
Kutcher
,
K.
Redmond
,
W. B.
Euler
, and
B. L.
Lucht
, “
Role of mixed solvation and ion pairing in the solution structure of lithium ion battery electrolytes
,”
J. Phys. Chem. C
119
,
14038
14046
(
2015
).
20.
S.-K.
Jeong
,
M.
Inaba
,
Y.
Iriyama
,
T.
Abe
, and
Z.
Ogumi
, “
Surface film formation on a graphite negative electrode in lithium-ion batteries: AFM study on the effects of co-solvents in ethylene carbonate-based solutions
,”
Electrochim. Acta
47
,
1975
1982
(
2002
).
21.
K.
Førland
,
T.
Førland
, and
S.
Kjelstrup Ratkje
,
Irreversible Thermodynamics: Theory and Applications
(
Wiley
,
Chichester
,
1988
).
22.
A. A.
Wang
,
S.
Greenbank
,
G.
Li
,
D. A.
Howey
, and
C. W.
Monroe
, “
Current-driven solvent segregation in lithium-ion electrolytes
,”
Cell Rep. Phys. Sci.
3
,
101047
(
2022
).
23.
H.
Lundgren
,
M.
Behm
, and
G.
Lindbergh
, “
Electrochemical characterization and temperature dependency of mass-transport properties of LiPF6 in EC:DEC
,”
J. Electrochem. Soc.
162
,
A413
(
2015
).
24.
R.
Chandrasekaran
, “
Quantification of contributions to the cell overpotential during galvanostatic discharge of a lithium-ion cell
,”
J. Power Sources
262
,
501
513
(
2014
).
25.
K.
Shigenobu
,
F.
Philippi
,
S.
Tsuzuki
,
H.
Kokubo
,
K.
Dokko
,
M.
Watanabe
, and
K.
Ueno
, “
On the concentration polarisation in molten Li salts and borate-based Li ionic liquids
,”
Phys. Chem. Chem. Phys.
25
,
6970
6978
(
2023
).
26.
S.
Viscardy
,
J.
Servantie
, and
P.
Gaspard
, “
Transport and Helfand moments in the Lennard-Jones fluid. II. Thermal conductivity
,”
J. Chem. Phys.
126
,
184513
(
2007
).
27.
T. W.
Maltby
,
B.
Hafskjold
,
D.
Bedeaux
,
S.
Kjelstrup
, and
Ø.
Wilhelmsen
, “
Local equilibrium in liquid phase shock waves
,”
Phys. Rev. E
107
,
035108
(
2023
).
28.
J. M.
Rubi
and
S.
Kjelstrup
, “
Mesoscopic nonequilibrium thermodynamics gives the same thermodynamic basis to Butler-Volmer and Nernst equations
,”
J. Phys. Chem. B
107
,
13471
13477
(
2003
).
29.
S.
Kjelstrup
and
A.
Lervik
, “
The energy conversion in active transport of ions
,”
Proc. Natl. Acad. Sci. U. S. A.
118
,
e2116586118
(
2021
).
30.
E. A.
Guggenheim
, “
The conceptions of electrical potential difference between two phases and the individual activities of ions
,”
J. Phys. Chem.
33
,
842
849
(
1929
).
31.
A. P.
Thompson
,
H. M.
Aktulga
,
R.
Berger
,
D. S.
Bolintineanu
,
W. M.
Brown
,
P. S.
Crozier
,
P. J.
in ’t Veld
,
A.
Kohlmeyer
,
S. G.
Moore
,
T. D.
Nguyen
,
R.
Shan
,
M. J.
Stevens
,
J.
Tranchida
,
C.
Trott
, and
S. J.
Plimpton
, “
LAMMPS—A flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales
,”
Comput. Phys. Commun.
271
,
108171
(
2022
).
32.
L.
Spitthoff
,
A. F.
Gunnarshaug
,
D.
Bedeaux
,
O.
Burheim
, and
S.
Kjelstrup
, “
Peltier effects in lithium-ion battery modeling
,”
J. Chem. Phys.
154
,
114705
(
2021
).
33.
J. G.
Kirkwood
and
F. P.
Buff
, “
The statistical mechanical theory of solutions. I
,”
J. Chem. Phys.
19
,
774
777
(
1951
).
34.
A.
Ben-Naim
,
Molecular Theory of Solution
(
OUP Oxford
,
2006
).
35.
D.
Bedeaux
,
S.
Kjelstrup
, and
S. K.
Schnell
,
Nanothermodynamics: Theory and Applications
(
World Scientific
,
2023
).
36.
W. L.
Jorgensen
,
J. D.
Madura
, and
C. J.
Swenson
, “
Optimized intermolecular potential functions for liquid hydrocarbons
,”
J. Am. Chem. Soc.
106
,
6638
6646
(
1984
).
37.
W. L.
Jorgensen
and
J.
Tirado-Rives
, “
Potential energy functions for atomic-level simulations of water and organic and biomolecular systems
,”
Proc. Natl. Acad. Sci. U. S. A.
102
,
6665
6670
(
2005
).
38.
L. S.
Dodda
,
J. Z.
Vilseck
,
J.
Tirado-Rives
, and
W. L.
Jorgensen
, “
1.14*CM1A-LBCC: Localized bond-charge corrected CM1A charges for condensed-phase simulations
,”
J. Phys. Chem. B
121
,
3864
3870
(
2017
).
39.
L. S.
Dodda
,
I.
Cabeza de Vaca
,
J.
Tirado-Rives
, and
W. L.
Jorgensen
, “
LigParGen web server: An automatic OPLS-AA parameter generator for organic ligands
,”
Nucleic Acids Res.
45
,
W331
W336
(
2017
).
40.
K. P.
Jensen
and
W. L.
Jorgensen
, “
Halide, ammonium, and alkali metal ion parameters for modeling aqueous solutions
,”
J. Chem. Theory Comput.
2
,
1499
1509
(
2006
).
41.
S. V.
Sambasivarao
and
O.
Acevedo
, “
Development of OPLS-AA force field parameters for 68 unique ionic liquids
,”
J. Chem. Theory Comput.
5
,
1038
1050
(
2009
).
42.
B.
Doherty
,
X.
Zhong
,
S.
Gathiaka
,
B.
Li
, and
O.
Acevedo
, “
Revisiting OPLS force field parameters for ionic liquid simulations
,”
J. Chem. Theory Comput.
13
,
6131
6145
(
2017
).
43.
D.
Frenkel
and
B.
Smit
,
Understanding Molecular Simulation
,
2nd ed.
(
Academic Press
,
San Diego
,
2002
).
44.
R. W.
Hockney
and
J. W.
Eastwood
,
Computer Simulation Using Particles
(
Hilger
,
Bristol
,
1988
).
45.
I.
Leontyev
and
A.
Stuchebrukhov
, “
Accounting for electronic polarization in non-polarizable force fields
,”
Phys. Chem. Chem. Phys.
13
,
2613
2626
(
2011
).
46.
L.
Martínez
,
R.
Andrade
,
E. G.
Birgin
, and
J. M.
Martínez
, “
PACKMOL: A package for building initial configurations for molecular dynamics simulations
,”
J. Comput. Chem.
30
,
2157
2164
(
2009
).
47.
N.
Molinari
,
J. P.
Mailoa
, and
B.
Kozinsky
, “
Effect of salt concentration on ion clustering and transport in polymer solid electrolytes: A molecular dynamics study of PEO–LiTFSI
,”
Chem. Mater.
30
,
6298
6306
(
2018
).
48.
W.
Shinoda
,
M.
Shiga
, and
M.
Mikami
, “
Rapid estimation of elastic constants by molecular dynamics simulation under constant stress
,”
Phys. Rev. B
69
,
134103
(
2004
).
49.
W. G.
Hoover
, “
Canonical dynamics: Equilibrium phase-space distributions
,”
Phys. Rev. A
31
,
1695
1697
(
1985
).
50.
S.
Nosé
, “
A molecular dynamics method for simulations in the canonical ensemble
,”
Mol. Phys.
52
,
255
268
(
1984
).
51.
J. E.
Basconi
and
M. R.
Shirts
, “
Effects of temperature control algorithms on transport properties and kinetics in molecular dynamics simulations
,”
J. Chem. Theory Comput.
9
,
2887
2899
(
2013
).
52.
S. H.
Jamali
,
L.
Wolff
,
T. M.
Becker
,
M.
de Groen
,
M.
Ramdin
,
R.
Hartkamp
,
A.
Bardow
,
T. J. H.
Vlugt
, and
O. A.
Moultos
, “
OCTP: A tool for on-the-fly calculation of transport properties of fluids with the order-n algorithm in LAMMPS
,”
J. Chem. Inf. Model.
59
,
1290
1294
(
2019
).
53.
Ø.
Gullbrekken
,
I. T.
Røe
,
S. M.
Selbach
, and
S. K.
Schnell
, “
Charge transport in water–NaCl electrolytes with molecular dynamics simulations
,”
J. Phys. Chem. B
127
,
2729
2738
(
2023
).
54.
P.
Ganguly
and
N. F. A.
van der Vegt
, “
Convergence of sampling Kirkwood–Buff integrals of aqueous solutions with molecular dynamics simulations
,”
J. Chem. Theory Comput.
9
,
1347
1355
(
2013
).
55.
P.
Krüger
,
S. K.
Schnell
,
D.
Bedeaux
,
S.
Kjelstrup
,
T. J. H.
Vlugt
, and
J.-M.
Simon
, “
Kirkwood–Buff integrals for finite volumes
,”
J. Phys. Chem. Lett.
4
,
235
238
(
2013
).
56.
J.
Milzetti
,
D.
Nayar
, and
N. F. A.
van der Vegt
, “
Convergence of Kirkwood–Buff integrals of ideal and nonideal aqueous solutions using molecular dynamics simulations
,”
J. Phys. Chem. B
122
,
5515
5526
(
2018
).
57.
J.
Newman
,
K. E.
Thomas
,
H.
Hafezi
, and
D. R.
Wheeler
, “
Modeling of lithium-ion battery
,”
J. Power Sources
119–121
,
838
843
(
2003
).
58.
J.
Landesfeind
and
H. A.
Gasteiger
, “
Temperature and concentration dependence of the ionic transport properties of lithium-ion battery electrolytes
,”
J. Electrochem. Soc.
166
,
A3079
A3097
(
2019
).
59.
S. K.
Ratkje
,
H.
Rajabu
, and
T.
Førland
, “
Transference coefficients and transference numbers in molten salt mixtures relevant for the aluminium electrolysis
,”
Electrochim. Acta
38
,
415
423
(
1993
).
60.
L. O.
Valøen
and
J. N.
Reimers
, “
Transport properties of LiPF6-based Li-ion battery electrolytes
,”
J. Electrochem. Soc.
152
,
A882
(
2005
).
61.
R.
Fingerhut
,
G.
Herres
, and
J.
Vrabec
, “
Thermodynamic factor of quaternary mixtures from Kirkwood–Buff integration
,”
Mol. Phys.
118
,
e1643046
(
2020
).
62.
D. A.
Jonah
and
H. D.
Cochran
, “
Chemical potentials in dilute, multicomponent solutions
,”
Fluid Phase Equilib.
92
,
107
137
(
1994
).
63.
E.
Ruckenstein
and
I.
Shulgin
, “
Entrainer effect in supercritical mixtures
,”
Fluid Phase Equilib.
180
,
345
359
(
2001
).