A novel approach to ligand-exchange rates applied to lithium-ion battery and sodium-ion battery electrolytes

A novel approach based on analyzing the forces and velocities of solvents and anions to compute ligand-exchange rates is here presented and applied to lithium-ion battery (LIB) and sodium-ion battery (SIB) electrolytes. By using ab initio molecular dynamics generated data, we find the ligand-exchange rates to increase as functions of electrolyte salt concentration and to be higher in SIB electrolytes as compared to LIB electrolytes. This indicates both that Na + transport will be more non-vehicular in nature and have improved kinetics vs Li + , and that increasing the salt concentration is beneficial. The systems studied were basically the first cation solvation shells of Li/NaPF 6 in propylene carbonate and acetonitrile using three solvent to salt ratios. Overall, the solvation shells are solvent rich at low salt concentrations, and as functions of concentration, the solvents are replaced by anions. As the SIB electrolytes display higher cation coordination and solvation numbers, we also expect an earlier onset of highly concentrated electrolyte behavior for SIB than LIB electrolytes. These observations should all have an impact on the design of electrolytes for optimal bulk properties, but also be useful with respect to interfacial dynamics. Here, we report on a study using ab initio molecular dynamics (AIMD) to study the structure and dynamics of the cation first solvation shell of LIB and SIB electrolytes of lithium/sodium hexaflu-orophosphate (Li/NaPF 6 ) in propylene carbonate (PC) and acetonitrile (ACN). Both PC and ACN are excellent solvents displaying high dielectric constants, wide liquidus ranges, and low viscosities, 13,22–27 and are of interest as HCE solvents. Especially, the dynamics of the solvation shells have been investigated by computing and analyzing the distribution of forces and velocities of solvents/anions to finally allow for a novel method for computing ligand-exchange rates from AIMD data to be presented.


INTRODUCTION
Lithium-ion batteries (LIBs) are the leading energy storage technology for portable electronics and electric vehicles (EVs). 1,2 The EV market is currently expanding rapidly, and as the world transitions from fossil fuel power to renewable energy sources, such as wind and solar power, the need for large-scale energy storage increases. Therefore, advances in LIB safety, energy and power density, cost, cyclability, and reliable long-term raw material supply are all highly needed.
Sodium-ion batteries (SIBs) show promises of decreased cost, improved power performance, and similar energy densities to LIBs. [3][4][5][6][7] Moreover, SIBs generally make use of abundant raw materials, thus, securing the long-term availability of resources. 5 With these features, SIBs have been suggested as a possible candidate for large-scale energy storage solutions, 8 and there are already some SIBs commercialized aiming toward the market of E-bikes. 9 From a different perspective, highly concentrated electrolytes (HCEs) have garnered great attention in recent years for both LIBs and SIBs. [10][11][12] The local HCE structure results in high rate capability, due to less resistive solid electrolyte interphase (SEI) layers, being derived primarily from anion reduction, an improved cation desolvation rate, and an overall high concentration of charge carriers. [13][14][15] The cation transport mechanism in HCEs is suggested to have the diffusion of the cation decoupled from that of the solvent, and high cation transport numbers have been reported, together indicating that vehicular motion is not the main transport mechanism. Instead, hopping by continuous ligand exchange has been suggested as an alternative transport mechanism. [16][17][18][19][20][21] Here, we report on a study using ab initio molecular dynamics (AIMD) to study the structure and dynamics of the cation first solvation shell of LIB and SIB electrolytes of lithium/sodium hexafluorophosphate (Li/NaPF 6 ) in propylene carbonate (PC) and acetonitrile (ACN). Both PC and ACN are excellent solvents displaying high dielectric constants, wide liquidus ranges, and low viscosities, 13,[22][23][24][25][26][27] and are of interest as HCE solvents. Especially, the dynamics of the solvation shells have been investigated by computing and analyzing the distribution of forces and velocities of solvents/anions to finally allow for a novel method for computing ligand-exchange rates from AIMD data to be presented.
The Perdew-Burke-Ernzenhof (PBE) functional was used throughout, 32 together with the separable dual-space Gaussian pseudo-potential by Goedecker, Teter, and Hutter, [33][34][35] except for Na where the pseudo-potential of Troullier and Martins 36 was used (as the former was not available). A time step of 4 a.u. ≈0.1 fs was used, the plane-wave cutoff was set to 70 Ry, and the temperature was set to 300 K. After equilibrating the systems, a Nosé-Hoover thermostat was turned on and remained on during the production runs. Following the advice in Ref. 29, the thermostat parameters in the PC and ACN systems were set to ωion = 1780 cm −1 and ω electron = 3560 cm −1 , and ωion = 2250 cm −1 and ω electron = 4500 cm −1 , respectively, and the target electron kinetic energy was set to the converged values in the equilibration runs.
The radial distribution function (RDF) is computed as where ni(r) is the average number of atomic species i in a shell of thickness Δr at a distance r from a central cation and ρi is the number density of atomic species i. The first minima in the RDFs are subsequently used to count the number of each atomic species i inside the solvation shell of each cation at each time step, thus, creating a distribution of instantaneous coordination numbers (CNs). The CNs are computed as the average of these distributions, while the variance of the distributions is used as a measure of how many different types of solvation shells are present in the electrolytes. A 95% confidence interval is computed as CNi ± 1.96σ/ √ C, where σ is the standard deviation of the distribution and C is the number of cations in the simulation. This is a conservative estimate of the error, assuming maximal statistical inefficiency. Similarly, partial RDFs are computed with respect to the center of mass of the solvents/anions, and the solvation number (SN) is computed in the same way.
Each solvent and anion are also studied with respect to the position, radial force, and radial velocity on their respective center of masses. In each time step, the radial forces/velocities, in the range [−2.55 nN, 2.55 nN] and [−2.45 km/s, 2.45 km/s], respectively, and distances between the cation and the center of masses, within [2 Å, 6 Å], are computed. A distribution showing the probability of experiencing a radial force/velocity at a given distance from the cation was created by placing each pair of radial forces/velocities and distances into equal size bins of sizes ∼0.07 nN × 0.02 Å and 0.035 km/s × 0.02 Å, respectively. The distributions were then normalized by the total amount of datapoints.
The structure of the cation first solvation shell is studied by computing the angle between the cations and two N/O atoms located within 3.2 Å of the cation using the cosine theorem. Similarly, the angle between the cation and the center of mass of two solvents within the solvation shell was computed.
To analyze the ligand-exchange rate, we propose that for a solvent/anion to be considered inside a solvation shell, their center of mass has to be within a certain critical distance rc, as determined by the minima of the appropriate RDF, from the central cation, while beyond rc, it is no longer a part of the solvation shell and an exchange has occurred, hence, inherently a one-dimensional problem. The radial motion of the solvent/anion center of mass was determined by the radial forces experienced by the solvent/anion and the velocities these forces generate.
The ligand-exchange rate is the inverse of the average time τ a solvent/anion remains in a solvation shell, i.e., the residence time. With ∫ Γ i dr/v(r) being the time it takes to travel along a path Γi from inside the shell to outside the shell, and by averaging over all such paths, we define the residence time of a solvent/anion in the solvation shell as where r is the distance between the cation and the solvent/anion center of mass. This expression can be re-written as (see the supplementary material for details and derivation) where r 0 is the equilibrium distance in the shell and ⟨v⟩ is the average radial velocity of the solvent/anion center of mass within the solvation shell. This is reminiscent of the concept of self-exchange velocity, in which, however, the residence time first is computed and subsequently used to compute an average velocity that subsequently is correlated with internal mobilities. 37,38 A 95% confidence interval for the residence time is computed asτ ± 1.96σ/ √ C × SN solvent , whereτ is an average of the residence time over the last 2 ps of the trajectories, σ is the standard deviation during these 2 ps, and of Chemical Physics ARTICLE scitation.org/journal/jcp SN solvent is the average number of solvents in the solvation shell. Again, this is a conservative estimate of the error, assuming maximal statistical inefficiency.
The residence time has also been calculated using the conventional approach, 39

by defining a function
where i enumerates the cations and s the solvents, and computing the auto-correlation function which is averaged over all cations, solvents, and times t. The autocorrelation function is then fitted to a stretched exponential where τ is the residence time and 0 < β ≤ 1 is a stretch parameter. The conventional method and our new method differ in their definitions; the former assumes the solvent to immediately be part of the first solvation shell once it is within rc of a cation.

RESULTS AND DISCUSSION
We start by describing the local structure of the different LIB and SIB electrolytes before moving on to the force and velocity distributions, and finally, the novel ligand-exchange rate analysis is outlined, applied, and discussed.

Local structure
By the different RDFs, we find that the LIB electrolytes show smaller solvation shells (peak at ∼2.0 Å) than the SIB electrolytes (∼2.4 Å), both with only small shifts as functions of salt concentration (Fig. 1). The peaks in the SIB RDFs are generally broader, supporting that solvents/anions have less strong interactions with Na + than Li + . 40,41 Comparing the solvents, one apparent difference is that while all the ACN RDFs drop to almost zero after the first peak, the PC RDFs never do-due to the bidentate coordination mode of PC made possible by the etheric ring oxygen atoms creating the broader features. 42 The shifts as functions of salt concentration seen in both RDFs originating in the anion, M + -P/F, show how the anions in the cation solvation shell get closer.
From the RDFs, we can derive that Li + has a total coordination number (CN) of ∼4, while for Na + , it is ∼6 ( Fig. 2 and Table I    both cations and solvents, however, the CNF contribution increases drastically with salt concentration, especially for the SIB electrolytes, and, hence, the amount of contact ion-pairs (CIPs). The variance of the CNs is larger for the SIBs, indicating a larger variety of Na + solvation shells, but the variance, indeed, increases for both cations and both solvents upon increasing salt concentration, and hence, more types of solvation shells are present in the HCEs. 42 This trend, however, was only clear once one of the cations was removed from each of the analyses of the LiPF 6 and NaPF 6 in ACN at the 10:1 concentration. The removed cations had solvation shells with up to three anions during the trajectories. These were deemed as highly unlikely structures in reality, and as the systems are small and the trajectories short, they would heavily skew the result. However, the results with those cations included in the analysis are included for completeness (Table S1) and show the same trend in the variance for the PC systems, but the trend is no longer visible in the ACN systems. There is also a distinct difference between the CNs and SNs (Table S2 and Fig. S1); for Li + , both are ∼4, while the SNs of Na + typically are ∼5, which is significantly lower than the CNs, a sign of more bidentate coordination with the solvents/anions for the latter. Finally, the higher SNs of Na + as compared with Li + means that (some of) the HCE properties are likely to emerge at lower salt concentration, i.e., lower solvent to salt ratios for the former.
In addition to the RDFs, CNs, and SNs, the distribution of angles can also be used to reveal the local structure, 42 and here, we find clear differences between the LIB and SIB ACN based electrolytes (Figs. 3 and S2). While the former has a peak corresponding to 105 ○ , consistent with a tetrahedral structure about the cation, the latter has a peak corresponding to 90 ○ and, hence, an octahedral structure. Moreover, the latter distribution is broader, extending all the way up to 180 ○ , while the former has more of a Gaussian distribution-quite consistent with the larger CN variance. Turning to the PC based electrolytes, the angle distributions are extremely similar, except for the lowest salt concentration, where the LIB electrolyte/Li + has a peak corresponding to ∼110 ○ , i.e., the tetrahedral structure, but the dominant feature is at a value corresponding to ∼50 ○ caused by bidentate coordination. 42 If we instead study the angle between the center of mass of two PC molecules and a central cation for PC (Figs. 4 and S3), we can remove the latter peak. Then, the distributions for the two cations are again very similar and very broad, with a main feature at a value corresponding to 80 ○ -90 ○ , showing that the solvation shell structure in the PC based electrolytes is less defined, likely caused by continuous exchange between monodentate and bidentate coordinations.

Force and velocity distributions
Moving to the distribution of forces acting on the solvent center of mass, these can be used to show differences between the first and second solvation shells (Figs. 5 and 6.) In the first solvation shells, at ∼3.3 Å and ∼3.6 Å for the ACN based LIB and SIB of Chemical Physics electrolytes, respectively, the forces are distributed around a very clear linear feature. The linear relationship between the radial force and the distance shows that the center of mass experiences a harmonic oscillator potential. In stark contrast, in the second solvation shell, the forces are evenly distributed around a radial force of zero, indicating that there is no average drive toward or away from the shell. Between the shells, there is a region of very low probability, as this space is only populated by solvents moving between shells.
Comparing the LIB and SIB electrolytes, the slope in the first solvation shell is greater for the former, which is consistent with the interaction strengths, and furthermore, the distribution is wider for the latter showing that the solvents are moving more freely. This can also be seen from the force distribution range being larger for the LIB electrolytes. Moreover, the area between the solvation shells is more populated for the SIB electrolytes, hence, more exchange between solvent shells-which can also be seen in the RDFs for the center of masses (Figs. S4 and S5), and this feature also increases as a function of salt concentration.
The PC based electrolytes show the same features; the forces are distributed around a linear feature in the first solvation shell, while they are centered around a radial force of zero in the second solvation shell (Fig. 6). Furthermore, the same differences are seen between the LIB and SIB electrolytes, but the distributions are much broader, especially for the latter. Moreover, there seems to be some segmentation of the distribution in the first solvation shell; the first peak in the center of mass RDF is far from symmetric, and an additional feature is seen for some of the electrolytes (Fig. S5)indicative of the solvent rotational motion caused by exchanging between monodentate and bidentate coordinations.
In the ACN based electrolytes, the solvent molecules in the first solvation shell experience a harmonic oscillator force and potential (Fig. 7), and similar potentials are obtained by computing the potential of mean force U(r) = −kBT ln g(r) from the center of mass RDFs (Figs. S4-S6), although the potentials are less smooth and a bit shallower. As a solvent enters the region between the first and second solvation shells, it experiences a restoring force, trying to pull it back into the first solvation shell. The restoring force, however, weakens, and very close to the second solvation shell, a weak force drags it into the second solvation shell. In the second solvation shell, the average radial force drops rapidly to zero, showing the solvent to be free from cation influence. Comparing the two cations, the first solvation shell harmonic oscillator potential well is deeper around Li + as compared to that around Na + , consistent with the overall differences in interaction strengths. 40,41 Furthermore, this shows that desolvation occurs more easily in an SIB electrolyte, which partly can explain the high rate capabilities and low cell impedances of SIBs. 3,50 Taking the perspective of a solvent in the second solvation shell, there is a small energy barrier to move inside the first solvation shell, but as it is only about twice the thermal energy, this indicates that a vacancy in the first solvation shell is filled quickly.
In PC, the same general features are seen, although not as clearly due to the high disorder in these electrolytes (Fig. 8). Nevertheless, a linear feature is present in the LiPF 6 in PC electrolyte at ∼4 Å. The same feature is not as clearly seen for NaPF 6 in PC, but there is a hint of a linear feature between 4.2 Å and 4.5 Å. The potentials again show clear, but irregular energy wells. The barriers are lower in the PC based electrolytes as compared to those in the ACN based ones, indicating an increased ligand-exchange rate in these systems, as well as easier desolvation.
In a fashion similar to that for the force distributions, the velocity distributions can also be analyzed (Figs. S7 and S8). They indicate the same general behavior; the solvents are more confined in the Li + solvation shells, and the solvents move faster in the Li + solvation shells.

Ligand-exchange rates
Several ligand exchanges are seen for a majority of the solvation shells. These occur through either of two main processes: an associative, where a solvent/anion enters the first solvation shell and of Chemical Physics another solvent/anion subsequently exits the shell, or a dissociative, where a solvent/anion exits the shell, leaving a vacancy that is subsequently filled by a new solvent/anion. The associative is more prevalent, often preceded by several attempts by the solvent/anion to enter the shell (Fig. 9).
In the further mathematical analyses, we start from the fact that the residence time can be directly computed from the center of mass RDFs and the average radial velocity in the solvation shell. The residence time τ for the solvents has converged for most electrolytes within 10 ps of simulation time, while for a few, a small upward drift by a few % per picosecond is still seen in the trajectories (Figs. 10 and 11). While we find no clear connection between the amount of data and convergence, the HCEs are slower to converge and have more solvation shells. This can be explained by the large CN variance, which indicates that the time needed to adequately sample the phase space grows as a function of salt concentration. Although the anions show a behavior similar to the solvents, they are far more scarce in the first solvation shells, and hence, much longer or larger runs are needed for their residence times to converge.
Using the average of τ during the last 2 ps of each simulation, we find that for the ACN based electrolytes, τ is 60%-80% higher for the LIB electrolytes at each single concentration. As the ligandexchange rate is the inverse of τ, it is 60%-80% greater for the SIB electrolytes (Fig. 10), which is consistent with the lower potential barriers (Fig. 7). As compared to previous MD simulations 35 for the same type of electrolytes, but at 333 K, our averages of τ are somewhat lower, 9-15 ps vs ∼50-100 ps. 26 Another study found a residence time of 25 ps, at 298 K, for ACN in dilute levels of Li + and O − 2 electrolytes, 51 thus, clearly of the same order of magnitude. Finally, compared to the conventional method of computing the residence time, our method is computationally very inexpensive. For the PC based electrolytes (Fig. 11), the τ are extremely similar for all electrolytes, except for the LIB 5:1 ratio electrolyte. Compared to the ACN based electrolytes, however, all τ are a factor of 3-5 lower, and hence, the ligand-exchange rate is much faster. The ligand-exchange rates are three orders of magnitude lower than those previously reported, 21 but the literature study used classical MD without any polarizable force-field. There seems, however, to Due to the small system sizes and the short trajectories, the auto-correlation functions ξ(τ) (Fig. S9), used to compute the residence time in the conventional way, are not well-defined exponentials for all intervals τ, especially for long intervals τ. Hence, the tail ends of ξ(τ) were not included in the fitting of the stretched exponential (Table S3). The residence times (Table II) are of the same order of magnitude as those computed by our method, and the ACN system follows the same trends with decreasing residence time as a function of salt concentration and shorter residence times for Na + than for Li + . There are, however, several discrepancies in the results for the PC systems, especially the outliers, LiPF 6 in PC at 20:1 and 5:1 concentrations and NaPF 6 in PC at 5:1 concentration, which all have stretch parameters far from 1.
Finally, the residence time is analyzed through Kramers' theory, 53 where m is the mass of the solvent, D is the diffusion coefficient of the solvent, ωB = √ U ′′ (xB)/m is the frequency of oscillations at the top xB of the potential energy barrier, ω 0 = √ U ′′ (x 0 )/m is the frequency at the minimum x 0 in the potential energy ( Figs. 7 and 8), and EB is the energy barrier (Table S4). For LiPF 6 in ACN, Kramers' theory yielded the residence times 16.2 ps, 19.3 ps, and 6.0 ps for the 20:1, 10:1, and 5:1 concentrations, respectively, and for NaPF 6 in ACN, 8.0 ps, 4.7 ps, and 5.0 ps for the 20:1, 10:1, and 5:1 concentrations, respectively. In the PC electrolytes, the potential wells are not as well defined, especially around the top of the barriers, making an accurate estimate of ωB difficult. The values in the ACN electrolyte which is the weak friction limit of Kramers' theory. 53 Assuming that the weak friction limit holds in the PC based electrolytes, we find residence times of 4.6 ps, 4.5 ps, and 5.6 ps in LiPF 6 in PC at the 20:1, 10:1, and 5:1 concentrations, respectively, and 10.3 ps, 8.5 ps, and 6.1 ps in NaPF 6 in PC at the 20:1, 10:1, and 5:1 concentrations, respectively. Comparing these results with the previous two methods for the residence time, none of the trends are seen, but the residence times are on the same order of magnitude and closer in absolute values to those arrived at through our new method utilizing the average radial velocity.

CONCLUSIONS
First of all, we find that the solvent center of masses moves in harmonic oscillator potentials inside the first solvation shell and that the potential barriers are lower for Na + vs Li + and for ACN vs PC systems, respectively. Hence, ligand exchange is more prevalent in these systems, facilitating faster desolvation, most likely also at interfaces, and these molecular level phenomena can be excellently connected to the macroscopic observations of lower impedance and better rate capabilities.
From a model point-of-view, the approach of using the computed center of mass RDFs of solvents (and anions) and their average radial velocity in the solvation shell is shown to enable the residence time and subsequently the ligand-exchange rate to be readily calculated. The computed residence times are of the same order of magnitude as in several other studies, and even though the potential barriers show no clear concentration dependence, the residence time decreases with salt concentration-especially in the ACN based electrolyte. Therefore, the transport mechanism can be concluded to be less vehicular in the HCEs, and desolvation is made easier. The latter should, in particular, be important at the interfaces where the ion concentration can be higher than that in the bulk, and this is then consistent with the high rate capabilities of HCEs. Finally, the method itself relies on very few assumptions and is, hence, generally applicable to ligand-exchange phenomena.

SUPPLEMENTARY MATERIAL
See the supplementary material for details on the density calculations, the CNs, including all cations in the analysis, Table S1, the distribution of SNs, Table S2, Fig. S1, and the angle distributions as functions of angle (Figs. S2 and S3). The center of mass RDFs are found in Figs. S4 and S5, as well as the potential of mean force, Fig. S6, derived from these. Figures S7 and S8 contain the solvent radial velocity distributions, and a mathematical proof of the presented method of calculating the ligand-exchange rate and additional computational details are included. The auto-correlation functions as well as the best stretched exponential fits used for computing the ligand-exchange rate using the conventional method are found in Fig. S9 and Table S3. The mean-square displacement and the parameters used for computing the residence time using Kramers' theory are located in Fig. S10 and Table S4.