Toward Bottom-Up Understanding of Transport in Concentrated Battery Electrolytes.

Bottom-up understanding of transport describes how molecular changes alter species concentrations and electrolyte voltage drops in operating batteries. Such an understanding is essential to predictively design electrolytes for desired transport behavior. We herein advocate building a structure-property-performance relationship as a systematic approach to accurate bottom-up understanding. To ensure generalization across salt concentrations as well as different electrolyte types and cell configurations, the property-performance relation must be described using Newman's concentrated solution theory. It uses Stefan-Maxwell diffusivity, ij , to describe the role of molecular motions at the continuum scale. The key challenge is to connect ij to the structure. We discuss existing methods for making such a connection, their peculiarities, and future directions to advance our understanding of electrolyte transport.

The primary function of a battery electrolyte is to transport[1−5] working ions from one electrode to another to participate in electrochemical reactions at electrode/electrolyte interfaces. Given local charge neutrality and thermodynamic constraints, e.g., the Gibbs–Duhem relation, other electrolyte species also transport along with the working ions (not necessarily in the same direction[6]). Such a coupled transport behavior causes species concentration gradients and voltage drops within the electrolyte. Depending on the battery application, either of these can be rate-limiting. For example, the concentration gradients limit the high current density operation of metal anodes,[7−14] wherein depletion of the electrodepositing ion in the vicinity of the anode can trigger nonuniform growth. Alternatively, in the context of extreme fast charging,[15−19] both concentration gradients and voltage drops influence the behavior of lithium-ion batteries. Here, the local concentration of lithium ions and the local electrolyte voltage regulate the distribution of intercalation reactions in the porous electrodes, cause unwanted side reactions such as lithium plating, and contribute to heat generation[20,21] within the cell.

Fundamentally, the evolution of the concentration and voltage fields is governed by transport properties such as conductivity, transference number, and diffusivity. Hence, the desired transport behavior can be achieved by suitably choosing electrolytes based on transport properties. Since intermolecular forces acting among ions and neutral species govern their motions at the molecular scale, any changes in the molecular constituents influence these motions. When averaged over a statistically significant population of the constituents, these motions manifest as continuum transport properties (Figure a). Thus, in principle, an electrolyte can be designed for the desired transport behavior by tuning its structure.[22]

Figure 1

(a) A general framework for bottom-up understanding of electrolyte transport. (b–d) Comparing Nernst–Einstein conductivities to measured values for three different electrolytes. Self-diffusivity and measured conductivities are obtained from the literature: LiPF6/PC from Hwang et al.,[36] LiPF6/EC/DEC from Feng et al.[37] and LiTFSI/PEO from Pesko et al.[38] Here PC ≡ propylene carbonate, EC ≡ ethylene carbonate, DEC ≡ diethyle carbonate, and PEO ≡ poly(ethylene oxide). (e–g) Changes in electrolyte transport properties by varying the Stefan–Maxwell diffusivities at constant concentrations of ions (c = 1 M) and solvent (c0 = 11.1 M); and fixed temperature, T = 298 K. ○ represents properties corresponding to 1 M LiPF6/PC from Hou and Monroe:[39]+0 = 5.11 × 10–7, –0 = 1.79 × 10–6, and ± = 3.21 × 10–7 cm2/s. Each curve is at a constant value of –0 (in (e), (f)) or ± (in (g)). The dashed curve is 1/10th of the solid curve value, while the dotted curve is 10 times this value. The expressions for t+0, , and κ in terms of are shown alongside. Here c = 2c + c0.

While such a bottom-up connection is intuitively clear, its usefulness is tied to achievable quantitative accuracy. Since the continuum transport behavior of the same electrolyte, i.e., a given structure, can differ depending on the choice of electrodes (e.g., porous or flat electrodes, intercalation[23] versus precipitation–dissolution electrodes,[24,25] etc.) as well as operating conditions such as current density and temperature, a structure–property–performance mapping as shown in Figure a offers the complete bottom-up understanding. Here performance refers to the spatiotemporal distribution of concentration and voltage fields for a given electrochemical context and how they influence the battery operation. Newman’s concentrated solution theory[23−30] provides an accurate property–performance connection for any electrolyte composition and electrochemical setting.[20,29,31,32] The use of Stefan–Maxwell diffusivity, , to describe species transport is central to the accuracy and universality of this theory.[32] For example, as shown in Figure b–d, unlike , the self-diffusivity-based Nernst–Einstein relation cannot explain the measured electrolyte conductivities in different electrolytes (self-diffusivity measurements are used to compute the Nernst–Einstein conductivity and compare it against the measured conductivities).

For every electrolyte, depending on the identity of the species, multiple ’s jointly represent the macroscopic properties such as conductivity. Since the structural variations influence all , the macroscopic properties cannot be tuned individually. For example, consider a simple electrolyte made up of three species: a cation, +, an anion, −, and a solvent, 0. Three Stefan–Maxwell diffusivities, +0, –0, and ±, jointly describe three macroscopic properties: cation transference number, t+0, ambipolar diffusivity, , and conductivity, κ. Figure e–g illustrate the effect of variations on each of these properties at fixed species concentrations. In addition to each of t+0, , and κ exhibiting different sensitivities to , the corresponding effects can counter each other. For example, an electrolyte with t+0 → 1 and a large is desirable to minimize concentration gradients; while t+0 can be increased by decreasing –0 (Figure e), this simultaneously decreases as shown in Figure f. Alternatively, ± variation can more strongly change κ than variations in (+0 + –0) in Figure g. Thus, instead of universal descriptors,[2,33−35] the property–performance relationship should be used to quantify the usefulness of each variation to the electrolyte transport behavior for the given electrochemical context.

The remaining question is to understand the structure–property relationship. Recently, various studies[37,40−65] have attempted to build such correlations for electrolyte transport. However, most of these have limited applicability given the use of self-diffusivities instead of the Stefan–Maxwell diffusivities. A key bottleneck is to quantify . Existing experiments[26,38,39,57,66−71] are cumbersome, and ’s have to be inferred from macroscopic properties. On the other hand, most theoretical approaches[72−74] do not compute . Either of these[26,38,39,57,66−70,72−74] is not readily extensible to electrolytes with any number of mobile species; e.g., an electrolyte for a Li–air battery[25,28] has four mobile species: cation, anion, solvent, and dissolved oxygen. There are two methods—one proposed by Wheeler and Newman in 2004[75] and another by Fong et al. in 2020[76]—that quantify based on molecular dynamics (MD) simulations and are extensible to any electrolyte system. These methods are similar to calculations from MD in general transport phenomena literature.[77−79] The Stefan–Maxwell transport in electrolytes obeys an additional constraint: charge neutrality, which must be accounted for in estimating from MD trajectories. Given this background, the following discussion explores these two methods to understand their peculiarities in the context of building accurate structure–property relations. Subsequently, we discuss future research directions to advance the bottom-up understanding of electrolyte transport for informing electrolyte design.

Molecular Dynamics (MD) to Quantify Stefan–Maxwell Diffusivities

One of the ideas for connecting ’s to electrolyte structure is to use classical MD to simulate an equilibrium structure for a given electrolyte system. According to Statistical Thermodynamics,[80,81] a state of macroscopic equilibrium has atomic and molecular constituents in a constant state of motion such that the system’s total energy is conserved. In principle, such molecular motions can be interpreted to quantify macroscopic properties defining nonequilibrium interactions,[82−86] such as thermal conductivity for heat transfer. As shown in Figure , given the trajectories computed in an MD calculation, the choice of the Wheeler[75] or the Fong[76] method may give different values since both use different definitions for displacements. To examine these differences, we use LiPF6/PC as an exemplar electrolyte system (PC ≡ propylene carbonate). All transport properties have been measured[39] for this electrolyte, which allows one to compare predictions against measurements.

Figure 2

A scheme to build structure–property relationships using molecular dynamics (MD). MD trajectories are analyzed differently in Wheeler and Fong’s approaches to compute the Stefan–Maxwell diffusivities, . While the thermodynamic state is specified unambiguously, prescription of the molecular structure and the computational parameters can change the computed for the same electrolyte. Nsalt and Nsolvent are respectively the number of salt and solvent molecules in the MD simulation. c is the continuum concentration, V is the domain volume during NVT calculations, and NA = 6.022 × 1023 #/mol is Avogadro’s constant. One explicitly specifies Nsalt (or Nsolvent) for the MD simulation, and the corresponding volume, V, is identified as a part of the simulation, which in turn prescribes what macroscopic salt concentration the particular simulation represents.

Stefan–Maxwell Diffusivities for LiPF6/PC Electrolyte: an Example

The trajectory data needed to estimate ’s for LiPF6/PC is generated using classical MD simulations with the GROMACS code.[87] Li+ cations, PF6– anions, and PC molecules are modeled using OPLS-AA force fields.[88−90] The ionic charges are scaled to 80%, i.e., Li+0.8 and PF6–0.8, to compensate for the charge transfer, polarization, and binding effects.[91,92] The initial configuration for each simulation (at a given number of ions, Ncation = Nanion = Nsalt, solvent molecules, Nsolvent, and a fixed temperature, T = 298 K) was packed randomly using the Packmol code.[93] The simulation domain is periodic in all three directions. Initially, the energy of each system is minimized using the steepest descent method. The force tolerance is 103 kJ/mol·nm, and the maximum step size is 10–2 nm. Subsequently, each system is relaxed as an NPT ensemble for 50 ns to stabilize the volume and system energy. The pressure is maintained at 1 atm using the Parrinello–Rahman barostat[94,95] with a time constant of 10 ps. Afterward, each system is simulated as an NVT ensemble for 300 ns to generate trajectories for diffusivity calculations. The simulation volume is kept constant and is identical to that identified at the end of NPT relaxation. For both the processes, the temperature is maintained at T = 298 K using the Nosé–Hoover thermostat[96,97] with a time constant of 1 ps. The time step for all simulations is Δt = 2 fs. The electrostatic interactions are computed using the particle mesh Ewald (PME) method,[98] while direct summation is used for the nonelectrostatic interactions, each with a cutoff length of 1.2 nm. All bonds with hydrogen atoms are constrained using the LINCS algorithm.[99]

The trajectories from the NVT run are used to quantify diffusivities. These trajectories are r⃗+,,r⃗-,, and r⃗0 time series where i, j = 1, 2, ..., Nsalt and k = 1, 2, ..., Nsolvent. Each of the coordinates is unwrapped to obtain corresponding position vectors for an infinitely large volume in all directions. Subsequently, the self-diffusivity calculations are straightforward. For every species, Δr⃗2 = (r⃗(t) – r⃗(0))·(r⃗(t) – r⃗ (0)) is computed and then averaged over all its molecules to compute an ensemble averaged value, ⟨Δr⃗2⟩, for example,

Subsequently, the slope of the ⟨Δr⃗2⟩ versus the t curve relates to the self-diffusivity. Here, a long enough time series is needed such that the slope approaches a stable value. Note that r⃗’s are 3D vectors and · represents a dot product. While the thermodynamic state is unambiguously defined by specifying Nsalt/Nsolvent and T for each simulation, one can choose many different Nsalt (and equivalently Nsolvent) values to perform these calculations. Once a large enough Nsalt is used, the slope becomes well-behaved. The time window and the number of molecules are not known beforehand and identified based on such comparisons.[100] Essentially, the large number of molecules simulated over a long enough time history ensures that the computed values represent continuum-scale transport properties.

In the Wheeler approach,[75] relative coordinates, R⃗+, = r⃗+, – r⃗0, and R⃗–, = r⃗–, – r⃗0,, are computed. After corresponding t = 0 coordinates are subtracted, the dot product of their ensemble averages, e.g., ⟨ΔR⃗+ ⟩·⟨ΔR⃗– ⟩, is used to compute the slope. Alternatively, the Fong approach[76] uses coordinates relative to the center of mass, e.g., R⃗+, = r⃗+, – r⃗cm, and the dot product precedes the ensemble averaging, for example, ⟨ΔR⃗+ ⟩·⟨ΔR⃗– ⟩. In each method, these averages are multiplied by volume, V, and the slope is computed. Subsequent steps in both approaches only involve information about the thermodynamic state and are independent of explicit V, Nsalt or Nsolvent values. Either of these methods differ from self-diffusivity calculations in the amount of trajectory information required. As shown in Figure for LiPF6/PC electrolyte, V⟨ΔR⃗+ ⟩·⟨ΔR⃗– ⟩ versus t in the Wheeler method or V⟨ΔR⃗+·ΔR⃗– ⟩ versus t in the Fong method must be averaged over multiple 10 ns simulations to obtain stable slopes. In comparison, much fewer simulations are required to reliably compute self-diffusivities. Thus, the Stefan–Maxwell diffusivity calculations require much longer simulations[78,101,102] and add to the complexities of studying them. This difference in Dself and predictions is fundamentally related to how the same trajectory information is analyzed differently. While each of the Δr⃗2 for self-diffusivity calculations is always positive, individual ΔR⃗·ΔR⃗ in Wheeler (or Fong) method can be positive or negative. This cancels out contributions from some of the molecules at a given time step and requires averaging over multiple simulations to obtain stable slopes.

Figure 3

MD calculated (a) self-diffusivities and (b, c) Stefan–Maxwell diffusivities using trajectory information from multiple simulations. Each simulation has a 10 ns trajectory data for all the molecules. ⟨Δr⃗2 ⟩ for self-diffusion, ⟨ΔR⃗⟩·⟨ΔR⃗⟩ for the Wheeler method and ⟨ΔR⃗·ΔR⃗⟩ for the Fong method are computed based on each 10 ns simulation. These ensemble averaged squared displacements are further averaged over # of simulations to estimate their slopes against time. All the simulations are for the identical thermodynamic state (Nsalt/Nsolvent = 1/10 and T = 298 K) of the LiPF6/PC electrolyte. For this electrolyte, Nsalt/Nsolvent = 1/10 is 1 M salt concentration. (d–f) The estimated diffusivity values based on # of simulations are compared against the corresponding values from 30 simulations to show that Dself calculations require much less trajectory information than predictions. Note that the y-scale in (a) is different than (b) and (c), while (d–f) are identically scaled.

As seen in Figure b,c, when applied to the same MD data, Wheeler and Fong methods give identical ’s. While both the methods have different definitions for the displacements and subsequent steps are markedly different, we find identical answers for a few Li+ electrolytes we have analyzed so far. Given the analytical differences, this equivalence may not hold for a different electrolyte class.

Figure plots the computed properties against experimental measurements and shows that ’s are comparable. The measured trend is much better captured for +0 in Figure a than for –0 in (b) and ± in (c). The predicted values change appreciably if the size of MD ensemble is varied at a fixed thermodynamic state (Nsalt/Nsolvent = 1/10 and T = 298 K) as shown by ★ at 1 M salt concentration in Figure a–c. These peculiarities in computed also translate to transport properties as shown in Figure d–f. For example, computed t+0 for Nsalt = 100 in Figure d changes suddenly around 1 M salt concentration, and the magnitude of this change is comparable to the variation with the size of the MD ensemble. Also note that the variation with the MD ensemble size is strongly non-monotonic (unlike the size effect in self-diffusivity estimates[100,103]). Given these variabilities, multiple long trajectory MD simulations are required to discern representative trends in the Stefan–Maxwell diffusivities, adding a considerable computational cost compared to predicting self-diffusivities. Note that, since the same MD data gives well-behaved Dself, the observed trends in Figure a–c are a characteristic of existing methods of analyzing the MD trajectories for . These variabilities are also found while simulating other electrolyte systems.[104−107] It is important to recognize that the computed transport properties in Figure d–f are comparable to the measurements,[39] when accounted for such variabilities.

Figure 4

Comparing MD predictions and measurements of (a–c) Stefan–Maxwell diffusivities and (d–f) corresponding transport properties for LiPF6/PC electrolyte. In each of the plots, ○ represent measurements,[39] and • are MD calculations at Nsalt = 100. For these calculations, Nsolvent is varied to obtain different salt concentrations. Additionally, at Nsalt/Nsolvent = 1/10 (i.e., 1 M salt concentration), Nsalt and Nsolvent are simultaneously varied to examine the effect of ensemble size on the computed properties. These values are shown as ★.

Variabilities Associated with the Attributes of the Electrolyte Structure

For a given electrolyte, the transport properties are ideally only a function of the thermodynamic state. Hence, the variations with Nsalt in Figure are numerical artifacts. Note that the time step, Δt, is chosen to be small enough to conserve the system energy and not introduce additional variabilities related to the numerical time integration. Apart from these, the attributes of the electrolyte structure in Figure represent an additional source of subjectivity.

In these classical MD simulations, force fields describe how the neighboring atoms interact with each other (a cutoff distance is typically prescribed to identify neighbors). Although simulations with general nonpolarizable force fields can often capture the basic physics of these electrolytes, unfortunately, most of the force fields are only benchmarked for selected structural and thermodynamic quantities.[88] Consequently, simulated and measured transport behavior in the concentrated electrolyte can greatly differ when electronic polarization or charge transfer becomes significant.[108,109] If such polarizability is ignored, the electrostatic interactions among the species intensify and slow down molecular motions. These limitations can be overcome by scaling up the partial charges or using polarizable force fields.[103,108,110] The polarizable force fields[103,108,110] are computationally expensive for the long time history MD calculations, as is the case for calculations in Figure . Alternatively, use of scaled ionic charges in nonpolarizable force fields, for example, Li+0.80 instead of Li+1.00, is meant to provide a mean-field representation of charge screening. Typical scaling factors[111−113] range from 0.5 to 0.8. While charge scaling is an adjustable parameter, it is not always sufficient to accurately capture species transport and can require additional adjustments, e.g., tuning solvent dipole moment, for simulating different electrolytes.[61,114−117]

Another structural aspect is coarse-graining, i.e., using one pseudo atom to represent a group of atoms.[117,118] An example is to represent the PF6– anion as one entity rather than explicitly treating individual P and F atoms. Such coarse-graining is popular for simulating large ensembles and long time scales. It is also valuable in simulating large polymer chains with tractable molecular details.[117] Besides simplifying the force-field description, larger Δt (often 20 times larger) compared to all-atom MD calculations can be used. While such a simplified treatment allows longer simulations,[119] one should verify if the relevant dynamics is captured.

When switching from polarizable to nonpolarizable to coarse-grained force fields, the number of force calculations decreases considerably at the cost of transferability of the corresponding force fields.[103,108] Ideally, one would expect an identical interatomic force field, for example, between Li+ cations and H atoms in solvent molecules, across multiple Li+ electrolytes with different solvents. Such force fields are referred to as transferable since they do not require any adjustments for describing interactions between the same pair of atoms when simulating across different molecules and electrolytes. The nonpolarizable force fields often require different charge scaling for Li+ cations if the anion or solvent structures vary widely across the electrolyte space of interest. Alternatively, the coarse-grained force fields would have to be explicitly tuned for every electrolyte. Note that the polarizable force fields are believed to be transferable but are not always available.[103,108]

Ab initio molecular dynamics[120−124] (AIMD) may seem like a promising avenue to resolve the subjectivity related to the structural representation. In AIMD, the force fields are computed “on the fly” from electronic structure calculations. However, it is computationally more expensive and limits the simulations to smaller ensembles and shorter time histories (tens of picoseconds), thus incurring different challenges in predicting since a representative set of ensemble configurations is not sampled. In turn, capturing the correlated transport[125] becomes challenging, especially in slow or dynamically trapped systems.[126,127]

Recently, a hybrid approach[128] is gaining popularity wherein AIMD is used to tune the force fields for the subsequent MD calculations. However, this would have to be performed for every electrolyte combination and can become computationally prohibitive. Separately, several accelerated MD approaches, e.g., temperature-accelerated dynamics, parallel replica dynamics, parallel trajectory splicing, etc., have been developed to sample more representative configurations.[129−131] However, their general usefulness for electrolyte transport predictions is yet to be verified.

Open Questions for Predicting using MD Simulations

The question of predicting for any electrolyte system is only partially resolved since

as compared to self-diffusivity calculations, multiple and much longer MD trajectories are required to obtain meaningful and account for associated variabilities;

accurately prescribing force fields for any electrolyte is ambiguous when exploring an electrolyte space.

While the later aspect is common to classical MD predictions of various properties, it is exacerbated in the context of predicting for electrolytes. For instance, polarizable force fields for ionic liquids[103,108] were tuned to ion self-diffusivities. However, we cannot directly tune force fields to ij since we only have measurements for a limited set of electrolytes,[26,38,39,49,57,66−69,71] and the MD calculations for are considerably more expensive. Most importantly, MD loses much of its usefulness if we tune the force fields against the very property we wish to predict. We require a consistent approach to tune the force fields against a property that is both inexpensive to predict and measure.[88,132] Such a property must be linked to transport since is a transport property. Note that as long as the force fields accurately predicting transport for every electrolyte system can be identified, they need not be transferable.

Thus, at present, the need for large MD data and having to identify a meaningful force field from multiple representations make computations expensive and can benefit from leveraging Machine Learning techniques[60,133−135] that do not sacrifice the physics, for example, tuning the force field without having to explore all possible combinations of the coefficients appearing in force expressions. Or, if the space of electrolyte exploration is kept systematic, e.g., glymes as the solvent (monoglyme, diglyme, triglyme, etc.), one may explicitly tune the force fields to a few representative candidates and then exploit the systematic variation to predict force fields for the rest of the combinations.

Outlook

Historically, electrolyte development has been a heuristic process that focused on increasing conductivity to improve transport. Our understanding of ion transport has advanced to show that the electrolyte transport is governed by multiple properties, not all of which improve with conductivity. Hence, we require an alternate approach that systematically connects the electrolyte structure to its performance in the cells, and in turn, makes it possible to accurately design electrolytes by tuning its structure. The critical piece connecting the structure to the performance is the Stefan–Maxwell diffusivity, .

The use of leads to accurate predictions of polarized electrolyte behavior via Newman’s concentrated solution theory. While this connection is clear and precise, our understanding of how relates to the structure is imperfect. At present, it is infeasible to measure for every electrolyte. Alternatively, the most promising route to predict uses molecular dynamics calculations. In principle, this route is universal across all possible electrolytes; however, the predictability is limited due to expensive calculations and ambiguities related to force fields describing the structure.

We believe that building an accurate structure–Stefan–Maxwell diffusivity relationship is the key question limiting bottom-up understanding of electrolyte transport in various battery chemistries. In turn, future efforts, experimental and theoretical, should focus on elucidating this connection. This question also has broader usefulness beyond batteries since multicomponent Stefan–Maxwell diffusion[136,137] is critical to many applications ranging from biological to engineering systems.

We acknowledge that the present discussion only focuses on the bottom-up understanding of electrolyte transport. The electrolyte also serves additional functions in batteries, especially defining reactions at the electrode/electrolyte interfaces.[5,138−142] The bottom-up understanding for each of these functions will require corresponding structure–property–performance mappings with yet to be agreed upon properties to connect molecular structure to continuum performance.