Understanding Ion Pairing in High Salt Concentration Electrolytes using Classical Molecular Dynamics Simulations and its Implications for Nonaqueous Li-O$_2$ Batteries
Abhishek Khetan Hamid Arjmandi Vikram Pande Heinz Pitsch Venkatasubramanian Viswanathan
physics.chem-ph


A precise understanding of solvation is essential for rational search and design of electrolytes that can meet performance demands in Li-ion and beyond Li-ion batteries. In the context of Li-O_2 batteries, ion pairing is decisive in determining battery capacity via the solution mediated discharge mechanism without compromising heavily on electrolyte stability. We argue that models based on coordination numbers of the counterion in the first solvation shell are inadequate at describing the extent of ion pairing, especially at higher salt concentrations, and are often not consistent with experimental observations. In this study, we use classical molecular dynamics simulations for several salt anions (NO_3^-, BF_4^-, CF_3SO_3^-, (CF_3SO_2)_2N^-) and nonaqueous solvent (DMSO, DME, ACN, THF, DMA) combinations to improve the understanding of ion paring with the help of a new metric of ion-pairing. We proposed a metric that defines the degree of clustering of a cation by its counterions and solvent molecules on a continuous scale, the limits if which are based on a simple and intuitive condition of charge neutrality. Using these metrics, we identify the extent of ion pairing in good agreement with experimental solvation phase diagrams and further discuss its usefulness in understanding commonly employed measures of salt and solvent donicity such as the Gutmann donor number. 


	
	
	
	
	


	
    < g r a p h i c s >


Keywords: ion pairing, solvation, molecular dynamics






§ INTRODUCTION

Despite having the highest theoretical specific energy among “beyond Li-ion” (BLI) battery chemistries, nonaqueous Li–O_2 batteries are yet to make significant technological strides toward commercialization due to several critical challenges such as low rechargeability and capacity <cit.>. The electronically insulating nature of the primary discharge product Li_2O_2, which passivates the battery cathode as it is formed, puts fundamental limitations on the extent to which discharge reactions can take place at a feasible rate. Recently, strategies to circumvent this issue by enhancing solubility of the reacting intermediates have been reported, which are based upon electrolyte formulations that use solvents <cit.> with high Gutmann Donor Number (DN) and Acceptor Number (AN) for stabilizing the primary reacting ions Li^+ and O_2^- and their solvated intermediates such as LiO_2 and LiO_2^-. However, it has been shown using thermodynamic analysis that aprotic solvents' ability to induce the solubility of these intermediates and thus provide for higher discharge capacity is anti-correlated to their stability against nucleophilic attacks <cit.>, H-abstraction <cit.> and electrochemical oxidation <cit.>. In this regard, the effect of Li-salts <cit.> and additives <cit.> has been also explored to in order to influence the solubility of the intermediates without compromising on the stability of the solvent, with reports of multifold capacity enhancement. For instance, the capacity enhancing effect of NO_3^- anions in comparison to TFSI has been reported widely in literature but uncertainities exist owing to several different explanations for this phenomemon, ranging from higher donicity and ionic association strength of the NO_3^- anion <cit.> to its catalytic effect on the electrodes <cit.> to the deactivation of the usually employed carbon cathode to further degradation <cit.>. 

Using mixtures of two commonly employed Li-salts, LiNO_3^- and LiTFSI, dissolved in dimethoxyethane (DME), we demonstrated in an earlier work that electrolytes containing a high NO_3^-:TFSI ratio resulted in higher discharge capacity <cit.>. Based on the principle that the free energy of the Li^+ ions and thus, the free energy of dissolution of the adsorbed intermediate LiO_2*, is largely dependent on the species that are present in the first solvation shell, we developed a modified Ising model for site occupancy in the shell. The estimated fractional occupation numbers resulting from this work provided a rational basis for selection of the total electrolyte, i.e., solvent and anion, where it was also shown that here is minimal capacity enhancement by changing the electrolyte anion in high-DN solvents. The modified Ising model served a useful tool for phenomenological understanding, however a more refined picture is required to estimate the solvation structure and dynamics to cover the broad spectrum of behavior spawned by varying salt-solvent combinations. 

The reactivity, selectivity and thermophysical properties of electrolytes depend strongly on the solvation behavior of the reactant ions, which is described by the degree of association of oppositely charged ions in electrolyte solutions to form distinct chemical species called ion pairs <cit.>. It has been shown that solvation and ion pairing behavior has far reaching implications on the electrochemical and chemical stability of both the nonaqueous solvent and the salt anion, which is intricately linked with the coupling between them <cit.>, and which cannot be described by simple atom centered descriptors such as Mulliken charges or electron affinity <cit.>. First order descriptors such as DNs have been successful in explaining trends in capacity and stability but are too coarse to be able to provide reliable trends for use in rational electrolyte design and discovery <cit.>. 

There is much subjectivity as to what forms an ion pair and it is often inferred from spectroscopic experiments or based on a certain distance of separation, r, between two oppositely charged ions in solution, which is smaller than some specified cutoff distance, R, beyond which ions are considered free <cit.>. Electrolytes in which ionic interactions become prominent, also termed often as“solvates” <cit.> remain difficult to characterize, as current theories cannot model their attributes adequately.  In this work we make use of rigorous classical molecular dynamics (MD) simulations to provide a general description of the solvation structure of Li^+ ions and the role of both solvent and anion in ion-pairing for several combinations of nonaqueous solvents and salt anions over a wide range of electron donicities and sizes, as shown in Fig. <ref>. Classical MD simulations help us take into account the effect of size, the averaged charge distribution within a molecule or ion, and their correlation with the intermolecular interactions and the resulting solvation structure.



Extensive studies of the solvation behavior of Li^+ ions using experiments and classical MD simulations for several different salt anions have been performed in solvent like acetonitrile (ACN) <cit.>, dimethyl sulfoxide (DMSO) <cit.>, carbonates <cit.> and glymes of different chain lengths <cit.>. However, computational descriptions of ion pairing typically consider ion-pairing from the solvation structure in the first solvation shell <cit.>, and it can be shown that the conventional choice of the radius defining this shell has limited applicability and cannot account for the extent of aggregation and clustering size.  In this work we propose a novel way to quantify the extent of ion-paring by defining a simple and intuitive metric based on charge neutrality of ionic clusters. Using this metric we offer a new generalized classification scheme for ion pairing and compare the various salt solvent combinations considered in this work. Additionally, we discuss the applicability of the commonly employed donor number in understanding solvation. This proposed method as such will be independent of the accuracy of forcefields or level of theory used for computation.



§ COMPUTATIONAL DETAILS

We make use of the GROMACS code v.5.1.2 <cit.> and OPLS-AA forcefields <cit.>. The forcefields for anions X^- =  NO_3^-, BF^-_4, CF_3SO^-_3 (OTF) anions was taken from the work of Acevedo  <cit.>, for (CF_3SO_2)_2N^- (TFSI) anion was taken from the work of Lopes and Padua  <cit.>, for DME was taken from the work of Watanabe  <cit.>. The forcefields for all other molecules were taken directly from the OPLS-AA forcefield repository <cit.>. The used forcefield parameters, topologies and assigned partial charges can be found in the Supplementary Information

All simulation boxes were chosen to be  5.5 × 5.5 × 5.5 nm^3 in size with periodic boundary conditions such that 100 Li^+ ions and salt anions in the box resulted in a total salt concentration of 1 M. After equilibration for 5 ns in an NVT ensemble and subsequent 10 ns of equilibration in NPT ensemble, the production runs for the simulations were performed in for a further 12 ns. The Parinello-Rahman barostat was used to maintain the system at a pressure of 1 bar and the Nóse-Hoover thermostat was employed to maintain the system at 298.15 K. The Particle-Mesh Ewald (PME) method was used to consider long-range electrostatic interactions, which were truncated at 1.3 nm along with those resulting from non-bonded Lennard-Jones interactions. All the analysis results were averaged over a couple of independent configurations of the same system. In addition to the all the solvent and salt combinations shown in Fig. <ref>, mixtures of  NO_3^- and TFSI of total 1 M concentration were also simulated in order to understand the competition between NO_3^-, which is a small, strongly associating anion to TFSI, which is a large, weakly associating anion.



§ RESULTS AND DISCUSSIONS

In order to gauge the ordering of ions and solvent molecules around the Li^+ ions, we first discuss and plot the time averaged normalized radial distribution functions (RDF) of all electrolyte constituents for the cases of 1 M LiNO^-_3 and 1 M LiTFSI in ACN, DME and DMSO in Fig. <ref>. The RDFs for other salt-solvent combinations considered in this work can be found in the Supplementary Information. The RDF or pair correlation function of species j around species i as a function of radius r is defined as:


    g_ij(r) = ⟨ρ_j(r) ⟩/⟨ρ_j⟩_avg

where the term ⟨ρ_j(r) ⟩ is the density of j within an infinitesimal region at r and ⟨ρ_j⟩_avg is the averaged local density of j in the volume contained by r. The variations in g_ij(r) signify the relative probability of finding j at a distance r from the reference i, for e.g., a the closest peak in g_ij(r) indicates a particularly favored separation distance where j is likely to be found around i. The magnitude of g_ij(r) is a function of binning size for r and its absolute value cannot be used to make quantitative comparisons directly. In macroscopic homogeneous systems, g_ij(r) goes to unity at large r. In the present work, g_ij(r) was plotted after normalization with respect to its maximum value for purposes of representation. 



On comparing Fig. <ref>a, Fig. <ref>b and Fig. <ref>c with NO_3^- as anion across the three different solvents, it can be observed that the nearest of all maxima always belonged to the NO_3^- anion at r = 0.198 nm (ACN), 0.196 nm (DME, DMSO) followed very closely by that solvent ∼ 0.01 nm thereafter. This behavior clearly demonstrates that NO_3^-, which is known to be strongly associating, dominates the immediate solvation surroundings of the Li^+ ions irrespective of the solvent. The radial distance marking the first minimum in the RDF, or alternately the radius of the inner solvation shell, is found to be at  r_sh^NO^-_3 = 0.236 nm (ACN), 0.228 nm (DME), 0.226 nm (DMSO). On comparing Fig. <ref>d, Fig. <ref>e and Fig. <ref>f with TFSI as anion across the three different solvents, it can be observed that the first maximum belonged to TFSI only in the case of ACN at r = 0.208 and the first solvent maxima preceded that of the TFSI anion in cases of DME and DMSO. This behavior demonstrates that for the weakly associating TFSI anion the high donor solvents can play an equal or dominant role in the liquid ordering around Li^+ ions. The radius of the inner solvation shell in these cases was found to be  r_sh^TFSI  = 0.298 nm (ACN), 0.292 nm (DME), 0.294 nm (DMSO). 

The RDFs do not directly reveal the occupancy or coordination numbers (CN) of electrolyte constituents around any Li^+ ion. Cumulative values of each coordinating specie can are obtained by the integration of the function g_ij(r), as shown in Fig. <ref> for the cases of 1 M LiNO_3 and 1 M LiTFSI in ACN, DME and DMSO. The CNs are then simply the cumulative values of the respective components at the respective inner shell radii r_sh. In line with the RDFs, the cumulative values show a much larger presence of the NO^-_3 anion closer to the Li^+ ions, which decreases as the DN of the solvent increases in the order ACN < DME < DMSO. This is again in stark contrast to the case of TFSI anions where the solvents dominate in presence around the Li^+ ions and there is not much difference in the calculated cumulative values. In addition to the cumulative numbers of the electrolytes constituents, the total charge q_r contained in the radius r around any Li^+ ion is also shown, which is evaluated simply as the algebraic sum of the forcefield assigned atomic charges belonging to the atoms in the solvating cations and anions. q_r = 1 at r = 0 because of the reference Li^+ ion. Including the atomic charges from solvent molecules served only to smooth the plotted charge curves and was therefore not pursued further. 



The standard method of classifying ion pair formation as solvent separated (SSIP), contact (CIP) or aggregating (AGG) ion pairs is based on the coordination number (CN) of the counterion in the first solvation shell being between 0 to 1, 1 to 2 or more than 2 <cit.> . To examine the applicability of this scheme, we first plot the CNs of the anions and the solvents ACN, DME and DMSO in the inner solvation shell of Li^+ ions as a function of varying NO^-_3:TFSI ratio in the electrolyte composition, as seen in Fig. <ref>. The cumulative values and CNs for other salt-solvent combinations considered in this work can be found in the Supplementary Information.  As can be seen in Fig. <ref>a, the CN of NO_3^- in the inner solvation shell at r_sh increases linearly with increase in NO_3^- concentration and is an order of magnitude higher than the CN of TFSI which falls to 0 with increasing NO_3^- concentration. For both anions, the CNs decrease as the donicity of the solvent increases from ACN < DME < DMSO. The donicity of solvents or anions alone, however, does not decide its CN in the inner solvation shell. As can be seen in Fig. <ref>b, the CN of weakly donating ACN at lower NO_3^- concentrations is higher than those of DME or DMSO, indicating that size effects as well as the anion - solvent compatibility may also play a role. 



Based on this previosuly employed classification scheme, CIP formation can be inferred only for the cases where [NO^-_3] > 0.75 M in ACN. This assertion, however, is not consistent with room temperature experimental phase diagrams in which NO^-_3, CF_3SO^-_3 are already insoluble in ACN <cit.>, whereas TFSI, BF^-_4  form SSIP/CIP in ACN <cit.>, at ACN:LiX ratio of 1:20 (considered in this work). While the CNs do reveal the fractional ordering in the immediate solvation shell averaged over every Li^+ ion and over time, they do not reveal the extent of ion pairing and formation of aggregates. Besides, the definition of CN and its associated r_sh implies that CN of Li^+ vs. Li^+ is always 0, and hence, they do not contain any information about the self-distribution of Li^+ ions around itself. Therefore, classifications of ion pairs as SSIP, CIP or AGG are rather restricted to strict and limited models which do not include, for instance, clusters of ion pairs with infused solvent molecules, where counterions may not aggregate without necessarily being in direct contact with each other.

The inconsistency pointed out here becomes more evident on analyzing the snapshots from MD simulations in each case, as shown in Fig. <ref>. Several snapshots were taken from every simulation run to select the most representative insights about ion pairing for each of the systems. In every snapshot, we take out the solvent molecules and show the ions only in order to focus on ion-pairing between Li^+ ions and the salt anion. When comparing the snapshots from Fig. <ref>a, Fig. <ref>b and Fig. <ref>c with 1 M NO_3^- anion across the three different solvents, it can be observed that there is a high degree of aggregation in the low DN solvent ACN, which decreases rapidly as the donicity of the solvent increases. On the contrary, not much perceivable difference is observed for when comparing Fig. <ref>d, Fig. <ref>e and Fig. <ref>f, with the weakly associating TFSI as anion across the three different solvents. These observations, however, cannot be deduced using the simple scheme of using CNs for ion pairing classification as has been performed in previous studies.



In order to address this discrepancy, we propose a new metric for understanding ion pairing based on the more general and intuitive condition of charge neutrality. Any distinct cluster of ions and solvents will naturally be charge neutral within a certain tolerance. An analysis of the radial charge q_r contained within r, as seen in Fig. <ref>, reveals the charge distribution of the species dominating in presence in the cluster depending of the sign of the charge curve. Clearly with NO^-_3 as the anion, the fluctuation in the charge is much higher in the low DN ACN. These fluctuations damp down as the DN of the solvent increases from ACN to DME to DMSO. The same is corroborated from the snapshot in Fig. <ref> where the degree of clustering can be seen to drastically decrease. In case of the TFSI anions, the snapshots hardly show any variation in ion pairing across the solvents as do the charge curves from Fig. <ref>.

To quantify these ideas based on the condition of charge neutrality, i.e.  |q_r| ∼ 0, we define the radial cutoff limit as r_cl.  Instead of the actual physical size of the ion pair, r_cl can be understood as the radial distance beyond which all charge fluctuations can be assigned as tail end fluctuations outside of the limits of the cluster under consideration. In Fig. <ref>, a typical RDF and the corresponding fluctuations in the associated charge curve are depicted schematically as a function of the radial distance r with respect to a reference Li^+ ion. As can be observed, such points of charge neutrality usually occur on the charge curve as inflection points between the oppositely charged peaks, indicating the expected ordering where the Li^+ ion is first surrounded by counterions which are followed by compensating Li^+ ions. The same can be corroborated exactly from the RDF plots in Fig. <ref> for each case. Based on these definitions, we propose the cutoff limit of charge neutrality as the outermost inflection or charge neutral point on the charge curve beyond which the amplitude of charge fluctuation or simply the fluctuation width |Δ q_r| is < q_cut = 0.5.



The criterion of q_cut = 0.5 was chosen after careful analysis of all cumulative and charge plots for all the salt-solvent combinations because the fluctuations were found to dampen quite fast with respect to the distance thereafter, as evidenced in the charge plots the Supplementary Information. The value of q_r at r_cl is often close to 0 but not exactly 0, as shown schematically in the ideal case in Fig. <ref>, implying that there will naturally be some errors associated within this scheme. Subject to these requirements and assumptions, we propose that the number of counterions contained within the charge neutral cluster defined by cutoff limit r_cl will decide the degree of ion pairing on a more realistic and continuous scale. It can be intuitively ascertained that the presence of more than 2 counterions, and thus, ion pairs, within the cutoff limit would imply formation of clustered aggregates (Cl. AGG) that also include infused solvent molecules. A value below 2 would imply formation of SSIP or CIP, the distinction between which can be inferred more precisely from the conventional CNs. Besides these, this proposed metric reveals the amount of solvent infused in the clusters. Having established these criteria, which are independent from the choice of salt or solvent in the system, we plot the number of anions and solvents contained within the cutoff limit of the charge-neutral aggregation in Fig. <ref>a and Fig. <ref>b, respectively.



Using these values from the analysis it is easy to ascertain the extent of clustering leading to aggregation for the simulated salt and solvent combinations. In 1 M TFSI, the weakly associating anions form SSIP/CIPS, as evident from the values < 1 as well as the snapshots. Formation of clustered AGGs takes place at starting at [NO^-_3] >= 0.4 M  for the weakly donating. For the strongly donating DMSO, AGGs would probably form only at [NO_3^-] > 1 M. For the moderately donating solvent DME, AGGs begin to form only at [NO_3^-] >= 0.75 M. While it is clear that strongly associating anions will have a higher degree of aggregation, they cannot exert their effect beyond a certain limit in high DN solvents, as evident from the case of DMSO.

The number of solvents in all cases shows a unique feature where a sharp drop in the number of solvents in the cluster is observed  in every case for 0.1 M <= [NO^-_3] <= 0.2 M. This is indicative of a structural transition from a TFSI dominated solvation shell to a NO^-_3 dominated solvation shell, in which the number of solvents are fewer as expected. The number of solvents with increasing NO^-_3 concentration however, varies quite differently. While the number of ACN molecules increases to as much as the number of counterions, the number of DMSO molecules falls off. This is directly correlated to the ascertained cutoff limits for each of the cases, r_cl = 0.910 nm (ACN), r_cl = 0.628 nm (DME) and r_cl = 0.504 nm (DMSO) at [NO^-_3]= 1 M, in which case the immediate solvation shell is dominated by the strongly associating counterion. From these analyses, an estimation of the extent of ion pair formation on a much more continuous scale is obtained, and the same can also be qualitatively correlated to the visualizations in Fig. <ref>. The cutoff limits for all other cases are tabulated in the Supplementary Information. 

Finally, we interpret these estimations with respect to a commonly employed descriptor of electron donicity - the Gutmann donor number (DN). As can be seen in Fig. <ref>, we plot the solvent and anion donor numbers with respect to the two main metrics discussed in the work, CNs and number of anions in the cluster (CL). In Fig. <ref>a and Fig. <ref>b, the CNs of the anion vs. Li^+ in the first solvation shell are described with respect to the anion's and solvent's DNs, respectively. As can be seen, there doesn't exist a well defined correlation between the DNs and the CNs, where despite being the weakest donating anion, BF^-_4 anion has much higher CN in the solvation shell compared to TFSI and OTF. The trends in the CNs with respect to anions is more representative of the often discussed ionic association strength of salt anions <cit.>, which is known to follow the order TFSI < BF^-_4 ∼ OTF < NO^-_3. An anion of higher DN is expected to dominate in presence around the Li^+ ion, however, the coordination behavior in immediate solvation shell is subject to size and steric effects, as well as salt solvent compatibility. The CN of anions are found to have a much more well defined decreasing trend with respect to the solvents.



In comparison, the number of anions in the clusters as estimated by the metric proposed in this work show much better proportionality with respect to DNs of both the anion and solvent, as can be seen in Fig. <ref>c and Fig. <ref>d, respectively. This behavior can be rationalized on the grounds that first order descriptors are more representative of larger scale behavior and are not apt for direct comparison with solvation phenomena in the immediate solvation shell. For the considered solvent:Li ratios, as mentioned in the Supplementary Information, our results show good agreement with respect to solubility and aggregation behavior observed in all the anions for ACN <cit.> and DME <cit.>. Although clear trends can be inferred from the analysis of coordination numbers as well, the estimation of extent of ion-pairing and clustered aggregation is not directly possible from them.

The application of this scheme is independent of the level of theory or the accuracy of force fields used for simulating the molecular dynamics, although, it must be noted carefully that the classical MD simulations presented here have been performed with non-polarizable forcefields which generally tend to overestimate CNs <cit.>. In comparison to first principles methods, classical MD simulations are very inexpensive and can reveal trends on the basis of which electrolyte components can be analyzed and compared. Besides, this scheme can be applied to both dilute and concentrated electrolytes.



§ SUMMARY AND CONCLUSIONS


In the present work, classical MD simulations were employed to simulate highly concentrated 1 M Li salts in various nonaqueous solvents, both spanning over a wide range of electron doncity, polarizability and size. As a test case, 1 M mixtures of LiNO^-_3 + LiTFSI with varying NO^-_3:TFSI ratios were also simulated to understand change in solvation structure transitioning from a large weakly donating TFSI anion to a small strongly associating NO^-_3 anion. RDFs of each of the systems were analyzed to understand the ordering of anions and solvents around Li^+ ions and CNs were estimated by analyzing the cumulative distributions obtained from the integration of the RDFs.  The evaluated CNs were further employed to differentiate and predict the type of ion pairing for the test cases in ACN, DME and DMSO. It was argued that based on the estimated CNs, only 1 M NO^-_3, which is the most strongly associating salt considered in the work, would form CIPs in ACN, which is the most weakly donating solvent considered in this work.

In order to resolve this discrepancy, a new metric based on charge neutrality of ionic clusters was detailed. Considering the fact that within some nominal deviation all ionic clusters will essentially have net null charge, a new cutoff limit was devised on the basis of which new estimates for the number of anions, and hence ion pairs, in these charge neutral clusters were obtained. The numbers of ion pairs in such clusters were further used to classify and understand ion pairing categories and thus explore aggregation in a more generalized way which allowed inclusion of solvent molecules in those clustered aggregates. These estimates were then pitted against the donor numbers of both the solvent and the salt anions to show that CNs of anions are rather directly related to commonly understood association strength of anions rather than their donor numbers, which were found to be better correlated to the longer length scale estimated number of ion pairs in charge neutral clusters.

The present work offers a new strategy to understand solubility and ion pairing which can potentially be combined with current efforts to identify novel, stable electrolytes, including those in which organic molecules are entirely removed from the electrolyte <cit.>, to develop an electrolyte that could enable high-energy, long-life Li-O_2 batteries. Using these descriptions of ion pairing to model and calculate solvation thermodynamics using cluster models will be the focus of out next research efforts.






  * A. Figure S1: Plots for normalized radial distribution function, cumulative values and radially contained ionic charge for all solvent-anion combination considered in the present study.

  * B. Table S1: Cumulative numbers of Li^+ ions, solvent molecules and anions in the ionic cluster (CL) defined by the cutoff limit r_cl for 1 M mixtures of NO_3^-+TFSI in solvents ACN, DME and DMSO.

  * C. Table S2: Coordination numbers (CN), cumulative numbers of Li^+ ions, solvent molecules and anions in the ionic cluster (CL) defined by the cutoff limit r_cl for all salt anion-solvent combinations considered in this work.

  * D. GROMACS forcefield parameters for all considered solvents.

  * E. GROMACS forcefield parameters for all considered salt-anions.




A.K. and H.P. acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG) within the graduate school mobilEM. A.K. also acknowledges the grant of computational time for the project by the Jülich Aachen Research Alliance (JARA) via Project No. 99000538(JHPC4401). A.K. and V.V. acknowledge support from Convergent Aeronautics Solutions (CAS) project under the NASA Aeronautics Research Mission Directorate.





[pages=1-26]supp_info_f.pdf
