Ion Transport in Organic Electrolyte Solutions for Lithium-ion Batteries and Beyond

The performance of metal-ion batteries at low temperatures and their fast charge/discharge rates are determined mainly by organic liquid electrolyte phase for other alkali metals that are used in batteries (such as sodium, potassium, magnesium, etc.). Atomistic computer simulations can play a primary role and predict ion transport in organic liquids. However, to date, atomistic force fields and models have not been explored and developed exhaustively to simulate such organic liquids in quantitative agreement to experimental measurements.


Introduction
Electrochemical rechargeable batteries have received increasing attention as sustainable and environment-friendly energy storage devices and power sources. The performance of rechargeable metal-ion batteries depends on several parameters, such as the electrode material, energy density, charge/discharge rates [1], storage capacity and cycling stability, electrochemical stability window of the electrolyte [2,3], ion mobility in the electrolyte [4,5], etc. One method of increasing the energy density is by increasing the cation transference number in the electrolyte phase [6]. Fast charging also depends on ion transport and temperature. A key mechanism that influences fast charging is the ion transport: (i) through the electrodes, (ii) across the electrolyte/electrode interface for the anode and the cathode [7], and (iii) through the electrolyte phase, including solvation and desolvation of ions [1]. A low ionic conductivity of the electrolyte phase is a shortcoming that limits its performance and the use of metal-ion batteries in the industry. Such a low conductivity could be due to the low dissociation of ions in the solvents, low cation solvation, [8] or low ion migration (transference number) in the solvent. Although liquid organic (non-aqueous) electrolytes offer high conductivity (1-10 mS/cm) and well-dissociated ions for concentrations > 1 M, they have a cation transference number below 0.5, which implies that the total ionic conductivity results from the anion migration [6]. However, contradictory results have been reported regarding transference number measurements, concentrations, and temperature dependence by different methods [9][10][11][12]. Ion transport is based on ion solvation which is favored by aprotic solvents with a high dielectric constant [13,14], ion dissociation which depends on temperature [15], ion motion through the solvent, long-range electrostatic interactions, and viscosity of the solvent. A direct relation exists between the motion of the ions and the ionic conductivity [16,17]. For instance, if the ions move faster, the ionic conductivity of the solution increases. A low value of the diffusion coefficient of ions can limit the charge/discharge rates and can lead to large concentration gradients across the anode, which will develop side reactions and eventually lead to cell failure.
Although extensive studies and reviews are available in the literature regarding the solvation structure of metal cations (Li + , Na + , K + , etc.) [18,19] and their desolvation [20] energies in several aprotic solvents, a comprehensive understanding (by using both experimental and computer simulation methods) of the ion transport mechanism in electrolyte solution is required to help optimize the ionic conductivity of rechargeable metal-ion batteries. To the best of our knowledge, the number of works on ion transport beyond lithium-ion batteries is rather limited. The higher natural abundance of alkali metals other than lithium (such as sodium [21][22][23][24][25][26] or potassium (K) [27][28][29]) and their lower cost make them potential candidates for future energy storage devices. However, the lack of understanding of the ion transport mechanism can hinder their application.
In this review, we discuss the properties such as ion diffusion (self-diffusion and interdiffusion), dynamics, transference number, and ionic conductivity in the organic liquid (not polymer) electrolyte phases [30][31][32][33][34][35] ( not in the aqueous solutions, ionic liquids, polymer electrolytes, or electrolyte solutions in electrodes) for lithium-ion batteries (LIBs) and beyond, using both experimental and atomistic modeling methods. This review aims to provide fundamental insights into the current understanding of the effect of electrolyte transport on the battery design, its development, and its performance.

Experimental Methods
There are a few experimental methods and techniques that are used for evaluating the ion transport properties, such as the self-diffusion coefficient, the binary diffusion coefficients (interdiffusion), transference number of ions, and ionic conductivity of organic electrolyte solutions of batteries. The key points of these methods have been described very briefly below.

PFG-NMR, DOSY, and E-NMR
Pulsed-field gradient nuclear magnetic resonance (PFG-NMR) spectroscopy is an experimental technique that is used for detecting the self-diffusion coefficients of molecular species [36]. This NMR method overcomes the experimental limitations of the ordinary spin-echo method [37] when detecting small diffusion coefficients. Diffusion ordered spectroscopy (DOSY) is a two-dimensional representation of PFG-NMR [38], which has the advantage of separating multiple species with similar diffusion coefficients through their chemical shifts. Many building blocks can be used for constructing the PFG-NMR experiments [39,40]. For example, the spin-echo based PFG is the pioneer experiment [41] in which a stimulated echo enables a longer diffusion time (required for slow diffusing molecules) compared to the spin-echo [42], the double stimulated echo suppresses the convection artifacts [43], and a longitudinal eddy delay overcomes the echo distortion in the presence of eddy currents [44]. Although the pulse sequences for PFG-NMR experiments are getting large and complex, the principles remain unchanged: the pulsed-field gradients play two roles, namely, encoding and decoding the positional information of the NMR-active nuclei. During the encoding period, the nuclear spins located in the slices of the sample volume experience a gradient phase offset the magnetic field axis and are hence dephased. During the decoding period, the nuclear spins experience an opposite gradient phase offset, which effectively cancels out the dephasing and therefore reconstructs the signal, provided that these spins remain in the same slice. The displacement induced by self-diffusion leads to incomplete decoding of the position information and hence results in the attenuation of the signal. The degree of signal attenuation can be described by the following exponential relation [40]: where D is the diffusion coefficient, is the gyromagnetic ratio of the nucleus, is the duration of the magnetic field gradient pulse, is the strength of the magnetic field gradient, and ∆ ′ is the effective diffusion duration, which depends on the duration, shape, and interval (∆) of the gradient pulses, as well as the chosen pulse sequences [45]. The PFG-NMR spectra are collected with a set of incremental , whereas the remaining parameters are kept fixed. Although the self-diffusion of the ions and solvent molecules can be easily obtained from the PFG-NMR and DOSY experiments, electrophoretic NMR (E-NMR) is an indispensable technique for performing in situ measurements of ionic mobilities [46][47][48][49][50][51] under the influence of an electric field. E-NMR is an advanced version of PFG-NMR experiment, where an additional electric field is applied along the magnetic field. However, in contrast to PFG-NMR, the signal in E-NMR experiences a phase shift rather than attenuation. This contrast is due to the fact that the ensemble of nuclear spins belonging to the same ionic/molecular species experiences a coherent electrophoretic drift along the electric field in E-NMR, whereas the self-diffusion of these species in PFG-NMR is uncorrelated. The extent of phase shift is derived from the following exponential relations (Equations (2)-(4)) [51]: where ∆ is the duration of the electrophoretic drift, the strength of the electric field and is the electrophoretic mobility. The phase shift is calculated as ( ( ) − (0)), where (0) is the phase of the NMR peak in the absence of the electric field. During an E-NMR experiment, the electric field intensity, E, is varied, whereas the PFG related parameters are fixed. As the electric field has a negligible effect on the uncharged species, their NMR signals experience practically no phase shift. Therefore, the signals of these uncharged species can also be used as an internal reference frame to compensate for the flow effect [52,53]. The cation transference number in the dilute limit for a binary salt electrolyte, in which both ions are univalent (1:1 electrolyte) and assuming that the ions are completely dissociated, is given by the following equation [6,54,55]: where D+ and D− are the diffusion coefficients of the cations and anions, respectively. In this dilute limit, the transference number can be considered simply as the fraction of the total ionic conductivity that is carried by the cations [4].
However, in concentrated electrolyte electrolytic solutions [15], where ion association takes place, the diffusion coefficients predicted by the NMR contain contributions from not only the free ions, but also the ion pairs [8,9] since the salt is not dissociated. Thus, the diffusion coefficients take the following form [9]: where α is the degree of dissociation and Dpair is the diffusion coefficient of neutral ion pairs. The cation transference number (Equation 4), calculated using the diffusion coefficients of Equation (5) ( + , − ) when the salts are not completely dissociated, will give erroneous results [9].
However, compared to the electrochemical methods, the NMR method has the advantage of not requiring a binary as a prerequisite.In addition, even multi-component mixtures can be investigated using the NMR method [9,56].

Method of Concentration Gradient
This method involves the establishment and relaxation of a concentration gradient in a vertical cell. The concentration profile relaxes exponentially with time. In order to form a concentration gradient, two methods are followed: i) physically, by introducing a high-density solution at the bottom of the cell, having a less dense solution above it, or ii) by placing reversible electrodes on the top and bottom of the cell and applying voltage ( Figure 1) [57]. Ions migrate due to the electric field created and thus form a concentration gradient. This gradient can be analyzed for short times using a galvanostatic polarization method (semi-infinite diffusion) [58] or for long times (restricted diffusion) [59].

Figure 1
Schematic showing the principle of the diffusion coefficient determination. At the start of the experiment, C1 = C2, and after applying a constant current for a predetermined time, C1 > C2. Reprinted with permission from the article by Valøen [57]. Copyright 2005 Institute of Physics and IOP Publishing.
The salt diffusion coefficients can be calculated from the relaxation profiles obtained using the galvanostatic polarization experiments [60]. In particular, the galvanostatic method is comprised of three steps: (i) The cell is allowed to rest in the open circuit condition until equilibration, (ii) a small current is applied for a predetermined time, in which the concentration gradient is built, (iii) the current is switched off and the relaxation time of the open circuit potential is monitored [57]. However, the apparatus for such experiments should be isolated from vibrations [57].
The restricted diffusion was first developed by Harned and French [61]. It was later advanced by Newman and Chapman and was based on the long-time relaxation behavior, following an induced concentration gradient. In particular, Newman and Chapman used the refractive index in order to measure the concentration differences between two points of the vertical cell [62]. The binary diffusion coefficients calculated based on short or long time relaxation can be significantly different [63,64]. Other techniques that are used to monitor the relaxation of the concentration gradient or for measuring the transference number include the following:

Potentiostatic Steady-State Polarization
Such experiments can be conducted using a two-electrode symmetrical lithium polarization cell (two Li metal electrolytes and an electrolyte phase) separated by a porous separator [59,60]. After the potentiostatic steady-state polarization state is achieved, a linear concentration profile is established between the electrodes, and the short galvanic pulse polarization is designed such that the concentrations change only near the electrodes. The steady-state potentiostatic polarization experiment is based on the short-term relaxation behavior of the cell potential. After each individual polarization, an open-circuit voltage (OCV) phase of at least 3 h is applied to ensure complete relaxation of the concentration profile [64]. When a potential is applied, the measured current (I0) is a result of the migration of all charged species, which initially exists at uniform concentration throughout the cell. A constant current is obtained after the steady state (Iss) salt concentration profile is achieved, and the net anion flux is zero. The transference number is calculated as t+ = Iss/I0 [63], assuming that convection and instabilities at the Li electrodes can be neglected, the ions of the electrolyte are completely dissociated, and the applied polarization voltage is small [65]. However, in real cells, charge transfer and conduction through the dynamic passivation layer occur on the electrode. This layer affects the steady-state current by influencing the reactions that occur on the electrode [65].
Bruce and Vincent extended this method to electrolyte solutions with a variable ionic diffusion coefficient [59,61,62], based on the dilute solution theory, which neglects the activity coefficients [66], [67] and is not applicable to concentrated electrolyte solutions [15,62,64,68]. In such cases, the cation transference number,for small polarization potentials (< 10 mV), is given by the following equation [65]: where the subscripts 0 and ss indicate the initial and steady-state values, respectively, R is the sum of the charge transfer resistance and the passivating film resistance, ∆ is the polarization voltage, and I is the current. 0 and can be measured by recording the impedance spectra of the cell before the polarization, after the steady state has been reached, and the direct current (DC) bias potential has been removed [65].
For a fast estimation of the transference number, potentiostatic polarization is the best choice, but it has the disadvantage of being restricted in its applicability to binary electrolytes electrolyte solutions [7] and is not suitable for concentrated electrolyte solutions [63].
More discussion on other experimental methods for measuring the transference number based on polarization cell experiments can be found in references [63,64].

Galvanostatic Pulse Polarization
The polarization time in a galvanostatic pulse experiment has no impact on the long-term relaxation behavior of the concentration gradient [64]. The concentration gradient (concentration changes only in the vicinity of the electrodes [64]) is not measured directly, but the cell potential difference during and after polarization determines the exact potential difference at the current interruption [9]. The polarization currents in galvanostatic experiments are selected such that the current density is below 0.3 mA/cm 2 since this guarantees the absence of lithium dendrite formation [20,67]. Subsequently to these galvanostatic pulse experiments, a steady state potentiostatic polarization is applied. The cation transference number is calculated by (i) measuring the cell potential after galvanostatic polarization (ii) the salt diffusion coefficient (by the restriction diffusion method), (iii) and the determination of the concentration dependence of the potential difference, using the following relation [58]: where ∞ is the bulk concentration of the salt, F is the Faraday constant, D is the salt diffusion coefficient and dΦ/dlnc is the concentration dependence of the potential Φ, and m is the slope of a plot of the cell potential at the time of current interruption, which is obtained from ∆Φ vs jti 1/2 fitting, where j is the current density and ti is the polarization time [58]. The galvanostatic polarization method, for measuring the transference number, is more accurate for concentrated solutions, but it is time consuming since it combines three different measurements [9]. This method is restricted to binary electrolytes with the cation as the active species (anion flux is zero in the steady state), no convection, semi-infinite diffusion, and linear cell geometry. However, such a method can be applied to non-ideal, concentrated solutions [9,68] and can be the most straightforward method to measure binary diffusion coefficients based on long-term relaxation behavior after galvanostatic polarization pulse [64]. It is worth noting that even for low salt concentration when the ion pairing is low, the transference number calculated by the galvanostatic polarization method results a lower value than the one calculated by the NMR method [9,55].

Ferrocene Cell Measurements
Another method for measuring the transference number is based on experiments using a concentration cell and determining the thermodynamic factor, (1 + ∂ln ± ( ) ln ), of binary solutions [62], where f± is the mean molar activity coefficient, which is the average of the mean molar activity coefficients of the cation and the anion, by measuring the lithium concentration dependent potential in a lithium electrode versus ferrocene/ferrocenium redox couple used as an internal standard [63,69]. The concentration dependence of the transference number, t+, is usually not known a priori [11]. Alternatively, it can be assumed a constant transference number within a concentration range (δc) about an average concentration, c0 (δc << c0); thus, one can determine an average transference number within the differential concentration range t+(c0+δc).
This method uses only the assumption that the transference method can be assumed to be constant within a differential concentration regime [59]. For a fast and accurate calculation of the transference number, this method is superior to the polarization cell experiments [63].

Electromotive Force Method
Concentration cells having liquids of different molalities (m2, m1: m2 > m1) of the monovalent electrolyte salt (e.g., LiX), have a potential difference, Etrans, and the variation of the potential with concentration, dEtrans, that are given by [9]: where R is the ideal gas constant, T is the temperature, and y is the activity coefficient of the cations. The potentials Etrans and E increase with the increasing concentration difference between the constant high concentration and the varied concentration of the salt.
If the activity coefficients are available, the transference number is determined by the electromotive force measurements.
If the activity coefficients are not available, a measurement of a cell without a liquid junction (with a salt bridge between the half cells with different concentrations) is necessary [9]. The potential is given by The anion transference number is calculated using the following relation: The cation transference number for a binary electrolyte solution is t+ = 1 -t− This method assumes that the transference number is constant in the measured concentration regime and also requires an adequate organic (non-aqueous) salt bridge and reversible electrodes [9].

UV/Vis Absorption
The ultraviolet and visible (UV/Vis) absorption technique can be used for monitoring the relaxation of the concentration gradient [59] that has been introduced physically. In particular, it evaluates the change in concentration at a point that is at the one sixth distance of the liquid column height from the bottom of the column since by doing so, the selection errors on diffusion are reduced [62]. This technique can be useful for evaluating the diffusion coefficients in any liquid solution with a UV active species that have a diffusion coefficient, D, greater than 10 -6 cm 2 /s [70]. However, a limiting factor of this technique is the spatial resolution of the instrument.

Moiré Pattern Method
The Moiré pattern method is used for calculating the binary diffusion coefficient, D±, of salts in organic electrolyte solutions. This technique is based on the optical observation of the timedependent relaxation of the concentration profile after two electrolyte solutions with different concentrations are smoothly and instantly brought into contact in a diffusion cell [70,71]. Thus, a stable interface is formed between the two electrolyte solutions, and mutual diffusion occurs [24].

Rotating Disk Electrode Method
In the rotating disk electrode (RDE) experiments, the idea is to measure the mass transfer of lithium in an organic electrolyte by finding the limiting current densities of lithium electrodeposition (Li + + e -⇒ Li) on the RDE, which is suitable for finding the mass transport of salts in electrolytes [72][73][74]. In the RDE experiments, for systems with limiting current density (presence of a plateau), the diffusion coefficient can be calculated using the following Levich equation [65,66]: where I is the disk electrode current, n is the number of electrons transferred in the elementary reaction step, F is the Faraday constant, A is the disk electrode area (cm 2 ), D is the diffusion coefficient (cm 2 /s), ν is the kinematic viscosity (cm 2 /s), and ω is the angular velocity of the disk electrode (rad/s). The diffusion coefficient is obtained from the slope of the I versus ω 1/2 plot, which according to Equation (11), should be a straight line passing through the origin [43].

Electrical AC Impedance Spectroscopy
In order to measure the ionic conductivity of organic electrolyte solutions, the alternating current (AC) impedance technique can be applied. In particular, from the obtained Cole-Cole plots, the bulk resistance, Rb (Ω), can be determined, and thus the conductivity, σ, can be calculated using the relation σ (c, T) = t/Rb·A [75,76], where t is the thickness of the sample (cm), A is the area of the effective contact with the electrodes (cm 2 ), c is the concentration of the solution, and T is the temperature. The ionic transference number can be measured by using the DC polarization method and is expressed by the relation ti = 1-Ie/It , where Ie, It are the electronic and total currents, respectively [75].

Atomistic Computational Methods
There are certain computational methods that are used to study ion transport (dynamics), diffusion, conductivity, such as the classical molecular dynamics (MD) and ab initio molecular dynamics (AIMD) methods.

Classical Molecular Dynamics
Classical molecular dynamics (MD) involves the integration of the classic Newton equation, that governs the motion of particles, according to Equation 12 [10,[77][78][79][80][81][82]: where Fi is the force experienced by the particle i, mi is its mass, and vi is its velocity. This method can be used for predicting the structural as well as dynamical properties of atoms and molecules in the liquid state electrolytes. The dynamics (diffusion) of ions can be calculated from the mean square displacement (MSD) [83][84][85][86]: where 〈| ( ) − (0)| 2 〉 is the mean square displacement of the center of mass. Using Equation (13), classical MD and ab initio MD simulations can be employed for calculating the self-diffusion coefficients and the ion transport mechanism. The cation transference can also be obtained from the MD simulations using Equation (4). It is worth noting that the classical MD simulations allow for the modeling of longer trajectories and larger systems at less computational expense. D can be used for calculating the diffusion activation energy by an Arrhenius type of relation (Ea, similar to the migration barrier obtained from the density functional theory (DFT) nudged elastic band (NEB) simulations) expressed by Equation (14) [87-89]: Here, D0 is a temperature independent pre-exponential term, T is the temperature, and k is the Boltzmann constant. Ea is commonly obtained from the Arrhenius plots, i.e., from the gradient of ln Di versus 1/T plot. Finally, the ionic conduction of a species i ( ) can be obtained through the Nernst-Einstein relationship by assuming that the motion of the ions is uncorrelated, and depends on selfdiffusion [77]: where n = n + + n − is the total number of mobile i ions, and q is the species charge. However, in a real solution, the motion of the anions and the cations is not uncorrelated, and the Nernst-Einstein relation provides an upper boundary prediction for the ionic conductivity with values always higher than the values determined by the electrochemical AC method. Thus, we can calculate the ionic conductivity using the Einstein-Helfand method [77], by extracting it from the slope of the collective translational dipole moment, taking into account the cross correlations of the ions:

Ab Initio Molecular Dynamics
Ab initio molecular dynamics (AIMD) is one of the most precise methods for calculating the ionic diffusivity and conductivity [90,91]. Unlike MD simulations, interatomic interactions in AIMD simulations are implemented by establishing the Schrödinger equation using various approximation methods. These interatomic potentials have a significant advantage over the empirical potentials used in classical MD, in particular, when performing calculations for complex systems such as organic electrolytes. AIMD simulations can realistically explain the active chemical events accurately. For example, AIMD simulations can describe the bond breaking and association process by considering charge transfer and polarization effects during the simulations, which also allow the reproduction of the experimental conductivity of the ions in the electrolyte systems.
In the AIMD simulations, the accuracy of the interatomic potential depends closely on the level of the approximation methods used during the simulation for calculating the electronic structures of the systems. Recently, two major approximations have been considered widely, namely, the Born-Oppenheimer approximation, in which different treatments were used for the ions and electrons, and the Car-Parrinello approximation, in which the degrees of freedom of the ions and electrons are conjugated via fictitious dynamics variables [92]. Notably, the Car-Parrinello approximation requires less time for integrating the equation of motion compared to the Born-Oppenheimer approximation [93][94][95]. The major outcome of the AIMD simulations is the trajectories of the particles, from which the transport properties such as the diffusivity and conductivity can be inferred using the same equations as those detailed for the MD simulations. The primary shortcoming of an AIMD simulation lies in its critical computational cost relative to a classical MD simulation or indeed other ab initio strategies such as NEB calculations. Therefore, organic electrolyte systems have rarely been investigated using AIMD simulations since organic electrolyte systems consist of a large number of atoms.

Experiments
The PFG-NMR technique has been widely used for measuring the diffusivity of ions and solvent molecules in lithium-ion battery electrolytes [16,17,31,[46][47][48][49]96]. An experimental study was performed by Hayamizu [97] (Figure 2) to calculate the self-diffusion coefficient of Li + , PF6 -, and various solvents, such as diethyl carbonate (DEC), propylene carbonate (PC), and ethylene carbonate (EC) for a range of temperatures [87]. From their study, it was concluded that the order of diffusivities was Dsolvent > Danion > Dcation [15,87,98,99]. The self-diffusion and activity coefficient of Li + ion depended greatly on the temperature [87]. The temperature dependence was correlated to the Volger-Fulcher-Tamman (VFT) fitting equation (Equation (14)). The Li + ion was observed to have a smaller diffusion coefficient compared to the PF6anion, although it had a smaller size (van der Waals radii were 0.076 and 0.254 nm, respectively [100]). This was due to the fact that Li + was solvated, and thus the hydrodynamic radius was larger than that of the anion [101]. It was noted that the Nernst-Einstein relation (Equation (15)) was corrected by incorporating a degree of dissociation (α). The value of α was observed to increase from 0.17 to 0.71 as the EC fraction increased from 0 to 100 % [87,88,102]. For a 1 M lithium bis(trifluoromethanesulfonyl)imide (LiTFSI): PC solution, the value of α was reported to be 0.6 [103]. Another experimental effort denoted that Li + ions were free from the strong interaction with anions PF6anion were in a PC solution even at a high salt concentration (2-3 M) [104]. The Li + selfdiffusion coefficient exhibited a monotonically decreasing trend with increasing salt concentration [15,[104][105][106][107][108]. Furthermore, the self-diffusion coefficients of Li + and different spherical anions (PF6 -, BF4 -) in the PC solution were measured for different salt concentrations and were also found to decrease (Figure 3) [105]. The self-diffusion coefficient of the anion showed a slightly stronger concentration dependence than the one of the cation, for LiPF6-PC solutions. It was speculated that this behavior of concentration dependence was attributed to different solvation of cations and anions [9]. It is worth noting that the LiBF4-PC solution [109] exhibited higher diffusivity values than in the LiPF6-PC solution [65]. Thus, the higher conductivity that appeared in the LiPF6-PC solution could not be attributed to the higher ion diffusivities. The Li + ion forms weak ion pairs with PF6 -, but strong ion pairs with BF4anion [99]. In addition, the cation (Li + ) transference number of LiBF4-PC was observed to be higher for most salt concentrations. Thus, the ion association was weaker in LiPF6-PC solution than in the LiBF4-PC solution (Figure 3). The Li + transference number increased from 0.37 to 0.577 when the LiPF6 concentration decreased from 1.5 to 0.25 M [65]. There were variations in the transference number values extracted by the direct electrochemical method and the PFG-NMR measurement [9,35,65]. The differences could be due to the different experimental conditions [65] or due to the fact that PFG-NMR does not take ion association into account [55,110]. In particular, the calculation of the transference number using the NMR method leads to an increasing cation transference number as a function of the salt concentration, whereas the calculation using electrochemical methods (such as galvanostatic polarization), which detect only "free" ions, results in a decreasing cation transference number as a function of the salt concentration, a result that is in agreement with the electrostatic theories [111,112].
Notably, the Li + transference number was found to be almost constant (slight variations) as a function of the salt concentration (and temperature) in a mixture of PC: EC: dimethyl carbonate (DMC) solvents [57]. However, the binary diffusion coefficient was found to greatly depend not only on the temperature, but also on the salt concentration [57,64,69,113]. This was also confirmed by the Moiré pattern method in various PC solutions containing salts such as LiClO4 [70], LiPF6 , LiBF4 , and LiTFSI [71].
The Li + transference number in an EC:ethyl methyl carbonate (EMC) mixture, at 298 K, was observed to decrease with increasing LiPF6 salt concentration, reaching a maximum of 0.37 at a salt concentration of 0.2 M [60]. In such solutions, the dependence of the ionic conductivity with the LiPF6 concentration, c, could be fitted by a polynomial of the form σ = A1c 3 +A2c 1.5 +A3c [60,114]. In another study of EC:EMC solutions with different salts, such as LiPF6, LiBF4, and LiN(SO2C2F5)2, diffusion was investigated using the PFG-NMR method [115]. According to Capiglia et al. [115], the cation transference number ranged between 0.37-0.49, and the maximum value of 0.49 was obtained for 0.5 M solution of LiBF4 in the EC:EMC (2:8) solution. The PF6 − transference number was found to be higher than the one of tLi + , tPF6 − > tLi [115]. The Li + transference number was reported to be 0.34 in γ-Butyrolactone (GBL), independent of the LiPF6 concentration from 0.1-0.8 M [116]. The error in the experimentally measured Li + transference numbers (using the PFG-NMR method) was estimated to be 5 % and 20 % in binary and ternary solutions, respectively [117]. In addition, organic solutions of lithium bis(fluorosulfonyl)imide (LiFSI) in DMC, EC, or mixtures of them [55,118] were also investigated at 298 K using the PFG-NMR method [119]. As in the previous studies, Li + ion diffusion was observed to be the lowest (one order of magnitude lower than bis(fluorosulfonyl)imide(FSI)), despite having the smallest size, due to its strong coordination with the solvent molecules [119]. The diffusivity of ions in DMC was three times higher than those in EC [119]. Furthermore, the transport properties (self-diffusion coefficient and ionic conductivity) of LiPF6 and tris(pentafluoroethane)-trifluorophosphate (LiFAP) were measured in an EC:DMC mixture (50:50 w:w) [120]. The self-diffusion coefficient of Li + ion was very similar for both electrolytes and close to FAPanions. However, it was 1.6 times smaller than that of the PF6 − anions [120]. Following the Nernst-Einstein relation (Equation 15), a higher degree of dissociation was observed for LiFAP [121] compared to that for LiPF6 or lithium 4,5-dicyano-2-(trifluoromethyl)imidazole (LiTDI) [110] in an EC:DMC mixture, due to the higher delocalization of the negative charge in the FAP − anion [120]. The ion dissociation was observed to decrease with decreasing temperature and the neutral ion aggregates formed in the electrolyte solution affected the cation transference number [8,19,121,122].
Furthermore, numerous experiments were performed on the EC:DMC (1:1 w:w), EC:EMC (3:7 w:w), and EC-free electrolyte (EMC:FEC (19:1 w:w)) solutions for measuring their ionic conductivities, binary diffusion coefficients, D±, transference numbers in the temperature range of −10 °C -to 50 °C ), and concentration range of LiPF6 (0.1 M-3.0 M), as can be seen in Figure 4 [69]. In particular, they used a novel analysis for the determination of transference numbers based on the analysis of concentration cell potentials and the short-term potential relaxation after galvanostatic pulses in the symmetric lithium coin cells [63,64,69]. There was a qualitative agreement of the data by Landensfiend [69] with the data by Nyman et al. [60] , however the strong decrease in the Li + transference number observed with increasing LiPF6 concentration and decreasing temperature [69] was observed, in contrast to the data by Valoen and Reimers (measured t+ ≈ 0.3-0.4, though they mentioned the difficulty in achieving a precise result) [57], who assumed temperature-and concentration-independent transference number. Further, the binary diffusion coefficient, D±, of LiPF6-acetonitrile (ACN) and LiPF6-EC:DEC (1:1 w:w) [12] solutions followed the relations D = 3.018×10 -5 exp(-0.357c) and D = 2.582×10 -5 exp(-2.856c), respectively, as they were estimated by UV/Vis absorption [58]. For 0.1 to 1.5 M LiPF6-EC: DEC (1:1 w:w) solutions, using the concentrated solution theory [123,124] and the Darken relation [125] along with the PFG-NMR measurements, the degree of dissociation, α, was measured to be in the range from 0.61 to 0.37 [15]. It was hypothesized that these changes in the ionic conductivity as a function of the salt concentration [126] were caused by the structural changes that switched the conductivity mechanism from vehicular to a Grotthuss-type (ion hopping) or a mixture of both mechanisms [127]. The origin of the concentration exhibiting the highest conductivity, Cmax, was correlated to the eutectic composition in the salt-solvent phase diagram [127].
Furthermore, the self-diffusion coefficient of Li + with different anions (PF6 -, ClO4 -), in PC, EC: DEC, and EC:DMC solutions was estimated from the current density data according to the Levich equation [128]. The diffusion coefficient was observed to increase in the order LiClO4-PC < LiPF6-EC:DEC < LiPF6-EC:DMC with the salt concentration, with a maximum value of 1.39×10 -5 cm 2 /s for 1 M LiPF6-EC:DMC (1:1) solution [128].
In addition, the diffusion coefficients of 14 pure organic solvents, such as EC, tetrahydrofuran (THF), γ-Butyrolactone (GBL) [116] appeared to exhibit a good correlation with viscosity, according to the Stokes-Einstein relation [84,129]: where D is the diffusion coefficient, is the Boltzmann constant, T is the absolute temperature, is the viscosity, and rs is the Stokes radius (Figure 5a). The Stokes-Einstein relation shows that the diffusion coefficient in a liquid is inverse proportional to the viscosity of the liquid. Thus, the diffusion coefficient decreases with the salt concentration since the viscosity of the electrolyte phase increases [101]. The relationship between the diffusion of the electrolyte ions (Li + and TFSI -) in various solvents and pure solvents shows good linearity; the diffusion coefficient of solvents governs to electrolyte ions diffusion (Figure 5b). From the comparison of the diffusion coefficients, it can be seen that the TFSIanion diffuses faster than the Li + cation, although the size of a bare Li + ion is smaller than that of a TFSIion. This is because the Li + ions are solvated by the organic solvents with the average solvation shell containing two solvent molecules, whereas the TFSI − anions are not solvated [31]. The solvents can be classified into three categories depending on the correlation of the Li + and TFSI − diffusion coefficients and the ionic conductivity. One category is the cyclic carbonates and esters, which have a direct correlation to the ionic conductivity and the diffusion coefficients of lithium and TFSI. The second category is triglyme (TG), diglyme (DG), dimethoxyethane (DME), which are the glyme family (CH3O-(CH2CH2O)m-CH3), in which as m increases, the ionic conductivity and the diffusion coefficient increase. The third solvent category includes dimethyl carbonate (DMC), diethoxyethane (DEE), ethylene propionate (EP), and dioxolane (DOx), in which the ionic conductivity is small due to a stable ion-pair formation [8,31].

Simulations
In this section, we will focus on the atomistic modeling studies of electrolytes for LIBs. Several computer simulation studies are available in the literature [77,98,99,106,107,[130][131][132][133][134][135][136][137][138][139] that focus on diffusion in such organic liquid solutions. In particular, the DFT calculations and AIMD simulations have been used for simulating the diffusion of Li + within the liquid EC solvent [140]. Although the work of Pham et al. [133] claimed to be able to simulate the diffusion behavior of Li + in the EC solvent, the short-scale dynamics [141] of Li + ion, up to 10 ps, was actually simulated. In particular in that study [133], the diffusion coefficient of Li + was calculated to be 4.3 cm 2 /µs in a 0.23 M LiPF6 concentration in the EC solvent, which was in a different scale compared to the experimental values of Hayamizu [87]. In another study implementing the ab initio and classical MD simulations, a good agreement between the simulated Li + diffusion was achieved in comparison to the experiments at 300 K [87], a single Li + , PF6pair in 249 PC molecules were simulated in this study [142]. This corresponded to a 0.05 M PC solution, whereas the experimental measurements [87] referred to a 1 M PC solution [87]. Furthermore, in another study of ab initio molecular dynamics [143], the Li + ion dynamics was measured in a liquid EC solvent at a high temperature (450 K). Li + ion solvation and diffusion in PC [100], EC, EMC solvents or mixtures of them [122,144] were calculated from the ab initio [100] and classical reactive MD simulations [141]. In particular, a higher Li + diffusion was observed in a linear carbonate solvent (EMC) compared to that in a cyclic carbonate solvent (PC [100] or EC [141]) .
Only recently, the self-diffusion of Li + for different electrolyte concentrations was calculated in the range of 298-303 K, using the classical MD simulations, in PC solutions [99,107,136,145], EC solutions [130,131,[136][137][138][139]146], DME solutions [132], DMC solutions [137,139,147,148] or even mixtures of them, such as EC:DMC [137,145,146], PC:DME [145] or ternary (EMC:DMC:EC) mixtures [149]. At low LiPF6 concentration (0.1 M), the EMC-LiPF6 electrolyte was observed to exhibit higher ion dynamics than the DEC-LiPF6 and DMC-LiPF6 electrolytes (DMC-LiPF6 electrolyte exhibited the lowest) [147]. A polarizable force field for the LiTFSI-DME electrolyte quantitatively predicted the Li + , TFSI -, and DME diffusivities, for 0.5-1 M concentration, which were compared to the experimental values [132]. The Li + transference number was found to be 0.47-0.49 depending on the salt concentration [132]. The ionic conductivity of that solution increased up to 1 M concentration and then decreased [132]. A similar result for the behavior of the ionic conductivity, as a function of the salt concentration, was observed for the EC-LiPF6 solutions [131], using a nonpolarizable force field [134,150]. The diffusivity of both ions was found to decrease with the salt concentration following an Arrhenius type of relation [131].
In another study, it was found that in PC solutions, there are contact ion pairs and multi ion complexes formed at high concentrations, which impact the transport properties (diffusivity and ionic conductivity) [136]. However, it is well-known from several studies that non-polarizable force fields do not reproduce ion transport accurately. In particular, they underpredict the dynamic (diffusive) properties since they overestimate the ion-ion correlations [151]. Several atomistic computational studies with non-polarizable force fields have failed to predict ion transport (diffusion) in organic liquid electrolyte solutions (PC [99,134,136] and EC [131,134,137,147,150]).
Furthermore, the self-diffusion coefficient and ionic conductivity predictions using a polarizable force field for the LiPF6-EC:DMC electrolyte solution [146] were found to be in good agreement with the experimental measurements [31,152,153]. The simulated degree of ion dissociation was overestimated compared to the experimental values [127]. A monotonic increase in the degree of dissociation and the ion and solvent diffusion coefficients was observed on increasing the EC content of the solvent mixture [127]. An almost perfect agreement of the ionic conductivity dependence on the salt concentration, between a polarizable force field for the LiTFSI-PC and LiTFSI-DME solutions [154] and the corresponding experimental values [155] was observed, as can be seen from Figure 6. Karatrantos et al. [77] modeled PC solutions with a 1 M LiPF6 concentration and calculated the ion and solvent diffusion over a broad temperature range (353-253 K), and compared them with the experimentally measured values [87]. The solvent molecules and atoms were represented by an atomistic model, and the charges of the solvent and the ions were scaled [77]. Not only the PC diffusivities but also the calculated ions diffusivities were found to be in accordance with the experimental measurements performed by Hayamizu [87] over a range of temperatures (353-293 K) (Figure 7). It was shown that DPC > DPF6 -> DLi + following the experimental trend. By appropriate charge scaling, both the generalized amber force field (GAFF) and the optimized potential for liquid simulations all atom (OPLS-AA) force field could predict the ion and PC diffusivities for a range of temperatures. To the best of our knowledge, there have been no studies, even with polarizable force fields [108], that can predict the ion diffusion and PC diffusion that are so in-line with the experimental measurements.  (15)). Open symbols denote the static ionic conductivity measurements (Equation (16)) from the simulations data. Filled diamond symbol is the experimental value of the same 1 M LiPF6-PC solution [105]. Star symbols are the experimental values of 1 M LiPF6-PC solution [76]. Reprinted with permission from the article by Karatrantos et al. Copyright 2020 Elsevier B. V.
In addition, the ionic conductivity of 1 M LiPF6-PC solution at different temperatures was predicted using the two force fields (GAFF: with the PC solvent charges scaled by 90 % and electrolyte ions by 85 % and OPLS-AA : with the PC solvent charges scaled by 80 % and electrolyte ions by 90 %) and compared with the corresponding values obtained using the Nernst-Einstein relation (Figure 7b) (Equation (15)). The results were found to be in accordance with the experimental measurements by Hwang et al. [105] and Ding and Jow [76].
Furthermore, it was found that the DMC solution exhibited a higher diffusivity than EC (or EC:DMC mixtures) solutions and exhibited lower conductivity (contrary to the Nernst Einstein theory) [137]. This was because the clusters observed in the DMC solution had zero net charge, whereas the clusters in EC (or EC: DMC mixture) had a net charge of +1 (or -1) , and thus they were conductive. Such differences in the clustering behavior could be explained by the different dielectric constants of EC and DMC [137]. In ternary (EMC:DMC:EC) solutions, owing to the higher viscosity of a cyclic (EC) co-solvent, the solvent viscosity was increased. Thus, the diffusional motion of ions became slower with increasing EC fraction [137]. Similar behavior was observed for the PC:DME and EC:DMC mixtures [145]. The viscosity of the electrolyte solution and the solvation structure were changed by the EC co-solvent, contributing to the ionic conductivity in the opposite manner, meaning that the ionic conductivity was optimized at a certain amount of EC [137]. In mixtures of cyclic (EC) and linear (DMC) solvents, the exchange dynamics of the solvents depended on the mixture ratio. The rotational dynamics of the linear (but not of cyclic) solvent was found to be sensitive to the mixture ratio [130]. Density variations of the electrolyte solutions could also alter the ion mobility and the Li + solvation structure [106,139]. A comparison of the effect of the type of anions (spherical or linear anions) on the ionic conductivity was investigated for the LiBF4 and LiTFSI salts in adiponitrile (ADN) solvent [156]. The higher affinity of the BF4 − with Li + compared to that of TFSI − linear anion with Li + implied that the ionic conductivity decreased remarkably at various concentrations [156]. However, in all of the above studies, the diffusion coefficients, as predicted by the original non-polarizable GAFF [157,158] and OPLS-AA [134,159] force fields, were underestimated in comparison to the experimentally measured values.
It has been observed from simulations [107,135,136] as well as experiments [99,104] that the increase in the salt concentration beyond 1 M translates into a decrease in diffusion, as can be seen from Figure 8. In particular, the diffusion mechanism of Li + at room temperature was observed to be a mix of vehicular and structural diffusion in the 1-3 M concentration range (Figure 8), exhibiting an increase in the structural diffusive mode in the super-concentrated regime [57,107,127,138]. In addition, a polarizable force field was used for modeling the LiTFSI salts in the EC solutions [138], and the prediction of ions and EC diffusivities by MD were in fair agreement (self-diffusion coefficients predicted by MD were 30 % higher than those obtained from experimental measurements) with the experimental values [31]. Similarly, the ionic conductivity measured from the MD simulations was found to be 30 % higher than the experimental values. The degree of uncorrelated ion motion and the fraction of free ions decreased slightly with the temperature increase, indicating increased ion clustering with temperature increase [138]. Furthermore, ion aggregates [137,146] contributed approximately 50 % to the overall charge transport with another 50 % contribution from the diffusion of free ions [138]. In a similar study, but with the COMPASS II force field, the Li + ionic conductivity and diffusion were measured for various salt concentrations in the PC solvent [135] at room temperature and were found to be in quantitative agreement with the experimental data at 298 K.  [99,104]. Error bars, when smaller than the symbols for computed data are not shown [107]. Reprinted with permission from the article by Self et al. Copyright 2019 American Chemical Society.
Furthermore, classical MD simulations of LiTFSI and LiNO3 electrolytes in mixtures of 1,2dimethoxyethane (DME) and 1,3-dioxolane (DOL) solvents for Lithium-Sulfur batteries [160] showed a strong dependence of the Li + self-diffusion with temperature. In addition, the diffusion mechanism of Li2S4 in DOL was of the mixed vehicular type and hopping character [160,161].

Experiments
Although there are a plethora of works investigating the structure and solvation of Na + and K + ions [20,118,[162][163][164][165], specifically in the aprotic solvents used in Na-ion batteries (NIBs) [3] and Kion batteries (KIBs) [166,167], only a few experimental efforts [168,169] and even fewer computer simulation studies have been done that investigate diffusion and ion transport in the organic liquid electrolyte phase for such batteries. Results from Raman and NMR spectroscopy showed higher mobility of NaPF6 compared to LiPF6 in the EC: EMC mixtures [118]. The diffusion of the TFSI − anions in the above solutions was measured using a pulsed gradient spin-echo NMR (PGSE-NMR) [168]. It was observed that the PC diffusion coefficient was higher in the 1 M NaTFSI-PC solution than that in the 1 M LiTFSI-PC solution. A higher discrepancy in the value of the PC diffusion coefficient was observed for the 2 M solutions. However, the TFSI − diffusion coefficient in 1 M NaTFSI-PC solution was higher than that in the 1 M LiTFSI-PC solution, but only in the 353-313 K range; their diffusion coincided at lower temperatures. The experimentally measured ionic conductivity of the 1 M NaTFSI-PC solutions was higher than that of the 1 M LiTFSI-PC solutions in very dilute (≤ 0.25 M) or highly concentrated solutions (1.5 M, 2 M) [168]. In addition, it was shown that the viscosity (η) of 1 M NaTFSI-PC solutions was lower than that of the 1 M LiTFSI-PC solutions in the temperature range of 353-283 K [168].
Furthermore, the binary diffusion coefficients of NaPF6-EC:DEC solutions, for different concentrations, were calculated and compared to that of LiPF6-EC:DEC solutions by Landesfeind et al. (Figure 9) [169]. It was observed a higher binary diffusion coefficient (D±), Na + transference number, and ionic conductivity for the NaPF6-EC:DEC solution, specifically at higher salt concentrations (1 M), compared to the corresponding values to the LiPF6-EC:DEC solutions, due to the weaker electrostatic interaction of the larger Na + cation [169]. In addition, at high concentrations, higher conductivities were observed for the electrolytes containing larger amounts of Na + cation [73]. This trend was explained by the different viscosities of the electrolyte due to the weaker interaction between the larger Na + ion and the solvent/PF6 − [73]. This was also observed in a previous study for NaClO4 (LiClO4)-PC solutions for concentrations above 1 M [170]. Furthermore, the formation of ion pairs at high salt concentrations could reduce the free charge carriers [15,19,105,171].
Morales et al. studied the viscosity and ionic conductivity of 1 M NaPF6 in the binary EC: PC mixture (1:1), and the ternary 5EC:3PC:2DEC (5:3:2) mixture, including versions with 2 wt% fluoroethylene carbonate (FEC). Both the viscosity and the conductivity were found to follow the VFT relation. The 5EC:3PC:2DEC electrolyte exhibited lower viscosity due to the lower concentration of the more viscous PC solvent and the addition of the DEC solvent with its lower viscosity. The order of the ionic conductivity, in terms of the solvent mixtures, was 5EC:3PC:2DEC + FEC > 5EC:3PC:2DEC ≥ EC: PC ≥ EC: PC + FEC, which was the inverse of the viscosity trend and could depict the role of viscosity in the ion dynamics. However, at high temperatures, the ionic conductivity values were observed to be similar. Further, the Walden products (product of the molar ionic conductivity and solution viscosity), which were used for assessing the ion transport in ionic liquids, were calculated and observed to lie below the Walden ideal line. It was concluded from this interesting set of experiments that by reducing the ion-ion and ion-solvent interactions, the ion dynamics (mobility) would be enhanced [172]. Furthermore, for the KIB solutions, the KN(SO2F)2 salt appeared to have the highest solubility in the carbonate-ester and ether solvents among KN(SO2F)2, KN(SO2CF3)2, KPF6, KBF4, and KClO4 salts. In addition, the KN(SO2F)2-DME solution presented a higher ionic conductivity and lower viscosity in the highly concentrated region compared to the corresponding values for KN(SO2F)2-PC, KN(SO2F)2-GBL, and KN(SO2CF3)2-DME solutions [173].

Simulations
In this section, we will focus on the atomistic modeling studies of electrolytes conducted for NIBs [3] and KIBs [20,173]. It was shown by Pham et al. [133] that the short-scale dynamics, up to 10 ps, of the Na + ion in EC was faster than that of the Li + or K + ion due to its weaker solvation ( Figure 10). Large scale quantum MD simulations of NaFSA-DME electrolyte solutions with salt concentrations of 5 %, 10 %, and 40 % (superconcentrated solution) were investigated by Okoshi et al. [174]. The diffusion coefficients were extracted by simulations of 20 ps, and for the 5 % solution, the DME calculated diffusion agreed with the experimental measurements [31,175]. In the superconcentrated solution, the vehicular type of diffusion was prohibited, and the alternative pathway was the dissociation/association reactions [176] between the carrier ion and the solvent/anion. It was also shown that the ligand exchange reactions were responsible for the carrier diffusion in the superconcentrated solution. Such reactions were between 60-120 ps in the electrolyte solutions.
In a classical MD study of NaPF6-EC solution (for 1 M NaPF6 concentration), it was shown that the self-diffusion of ions increases with temperature and follows an Arrhenius type of relation (Equation (14)) [177]. In an extensive classical MD study of NaClO4 salts in various solvents, the transference number was observed to depend on the relative mobility of Na + to its counterion ClO4and was calculated, by MD, for a range of cyclic, linear solvents, or mixtures of them [178].
The Na + transference number varied from 0.3 to 0.51 depending on the solvent. The EC based solvents exhibited a high Na + transference number of ≈0.5. The ionic conductivity of the solvent depended on its dielectric constant since it dictated the degree association of the salts, and the viscosity of the solvent, which was inversely proportional to the ion mobility [178].
The Na + diffusion coefficient in DME solutions (having TFSIas anions) was predicted quantitatively, corresponding to the experimental values [132], using the classical MD simulations based on a polarized force field at 298 K. The Na + transference number was predicted to be 0.45-0.5 depending on the salt concentration (0.5-2 M) [132]. The Na + transport mechanism was less vehicular than that of Li + in the DME solutions. For a given salt concentration, higher diffusion coefficients were predicted for NaTFSI than for LiTFSI in DME solutions ( Figure 11) [132].

Figure 11
Self-diffusion coefficients of Li + TFSIin DME (circles) and Na + TFSIin DME (squares) as a function of the salt concentration for (a) Li + or Na + and (b) TFSI -. The computational results are represented by the filled symbols. Reprinted with permission from the article by Liyana-Arachchi et al. [132] Copyright 2017 American Chemical Society.
The experimental work conducted by Seki et al. demonstrated that Mg 2+ formed a much stable complex with PC and TFSAanion than Ca 2+ (Figure 12) [168]. Hence, the diffusion coefficient of PC and TFSAwas lower when the cation was Mg 2+ than when it was Ca 2+ . The same study revealed that the interaction energy between the metal cations and PC solvent and TFSAanion increased in the order Na + < Li + < Ca 2+ < Mg 2+ . Thus, for the same salt concentration, the self-diffusion coefficient of PC and TFSAanion exhibited the order Na + > Li + > Ca 2+ > Mg 2+ . However, it was not feasible to experimentally measure the self-diffusion coefficient of divalent and sodium cations [168]. The selfdiffusion coefficients were found to monotonically decrease with increasing salt concentration. A larger salt concentration dependence of the self-diffusion coefficient was observed for solutions with divalent cations than in solutions with monovalent cations [168]. The temperature dependence of the self-diffusion coefficient and ionic conductivity followed an Arrhenius type behavior. The maximum ionic conductivity was observed for 0.5 and 1 M salt concentrations for all electrolyte solutions (Figure 13a). The NaTFSA electrolyte showed the highest ionic conductivity for all salt concentrations since it had the lowest viscosity, , for all the salt systems.  In addition, the relationship between the self-diffusion coefficient of PC (DPC) and ionic conductivity can be seen in Figure 13b. In particular, the order of the ionic conductivity with both solutions have higher conductivity than the monovalent cation solutions for the same viscosity [168].

Simulations
In this section, we focus on the atomistic modeling studies of electrolytes for Mg 2+ , Zn 2+ , Ca 2+ , and Al 3+ ion batteries. A theoretical investigation performed by Okoshi et al. showed that the divalent Mg 2+ ion has a stronger solvation structure with different organic electrolytes compared to the monovalent Li + and Na + [166,174]. Likewise, Mg 2+ and Ca 2+ in the mixture of PC and EC electrolyte solutions consisted of stronger cation-carbonate interactions than Li + and Na + in the same solutions [3]. Classical MD and ab initio simulations of Mg 2+ and Ca 2+ in the EC, vinylene carbonate (VC), PC, butylene carbonate (BC), DMC, EMC, DEC, EC:PC, EC:PC, EC:EMC, and EC:DEC solutions demonstrated that Mg 2+ and Ca 2+ could be more soluble in these electrolytes than Li + and Na + [54,178,189,190]. Therefore, the values of the ionic mobility of the divalent in the same solvents were 1.7 and 1.5 mS/cm, respectively. However, when the concentration was increased to 1 M, the conductivity values changed inversely. At 1.0 M concentration, the ionic conductivity of Li + , Na + , Mg 2+ , and Ca 2+ was 6.3, 6.9, 1.9, and 2.8 mS/cm, respectively. The divalent cation Mg 2+ consisted of a higher polarizing power than Ca 2+ and exhibited lower ionic conductivity. The experimental, classical MD simulations and first principle studies demonstrated that the diffusion coefficient of the Mg 2+ had been changed by changing the anions in the same organic electrolytes, with the order Mg(BH4)2 > Mg(BF4)2 > Mg(TFSI)2 [191]. Furthermore, it was revealed that the solvation of the Mg 2+ ion in the different organic electrolytes was not necessarily related to the dielectric constant of the solvent molecules, although it had an intricate relationship with the molecular size of the solvents. In addition, it has also been reported that increasing the concentration of Mg 2+ in the organic electrolyte solutions could decrease the diffusion coefficient of Mg 2+ by cluster formation [192].
DFT calculation by Li et al. demonstrated that Zn 2+ could adjoin the TFSI − anion more strongly compared to Mg 2+ , and thus, in organic electrolytes, Zn 2+ had more ion-ion association than Mg 2+ [193]. Experimental and MD simulations showed that Mg 2+ and Zn 2+ had almost the same cationoxygen distance with the TFSI anion [191,194]. This was mainly because of the similar cationic radius size of Mg 2+ (0.86 Ǻ) and Zn 2+ (0.88 Ǻ) [195].
The dynamics of Zn 2+ ion in the organic electrolyte systems was directly related to the viscosity, dielectric constant, molecular geometry, and interaction strength between cation-solvent and cation-anion [191,196]. The computational and experimental work by Han et al. showed that the diffusion coefficient of Zn 2+ ions was directly influenced by the viscosity and dielectric constant of the different aprotic solvents, as depicted in Figure 14 [197]. As can be seen from the figure, the cations and the anions of the Zn(TFSI)2 and Zn(CF3SO3)2 salts in the PC solvent showed the lowest diffusion coefficient and conductivity compared to the other aprotic solvents, probably due to their highest dielectric constant and viscosity difference among the various organic electrolytes [198]. Therefore, it was assumed that the reduction in the viscosity difference of the organic solvent could increase the transportation of the Zn 2+ ion. The introduction of DMC in the 0.5 M Zn 2+ in a trimethyl phosphate (TMP) solution increased the conductivity from 4.58 mS/cm to 4.90 mS/cm by reducing the viscosity of the original solvent [199].

Conclusions
In this article, we have reviewed studies on ion and solvent transport in the organic liquid electrolyte phase that are used for lithium batteries and beyond. The diffusion mechanisms and the transference number, which can strongly govern the battery performance, are challenging and nontrivial to access experimentally [6,8]. Although extensive experimental studies focusing on lithium-ion batteries have been conducted, only a limited number of studies focusing on sodium, potassium, and divalent metal ion batteries have been done. This could not only be due to the difficulty of experimentally measuring sodium, potassium, and other cation diffusivities, but also due to the lack of appropriate atomistic computational models since it is well known that nonpolarizable force fields underpredict diffusion by approximately few orders of magnitude. In principle, polarizable force fields could quantitatively predict the diffusion mechanism, the transference number, and the ionic conductivity, but their calculations are computationally demanding and their development and implementation are currently still limited only to lithium electrolytes. An alternative to the polarizable force fields would be the use of non-polarizable force fields with parametrized charge scaling based on the experimental diffusivity and conductivity data. Such models would be successful in shedding light on the cation transport mechanism and the quantitative prediction of diffusion for the simple lithium electrolyte solutions. However, to date, such methodology has not been applied in solutions used in divalent or trivalent (Mg 2+ , Al 3+ , etc.) ion batteries, in which the electrostatic interactions are stronger, and therefore more challenging for the success of this method. However, taking into account the increasing demand for new batteries containing ions that are more abundant in nature than lithium, such as Na + , K + , Mg 2+ , etc. more research efforts should be devoted to studying their ion transport. In addition, the development of appropriate atomistic computational models could prove to help shed light on the diffusion mechanisms and the concentration and temperature dependence of ion transport and thus help in speeding up the design of such next generation batteries. Another issue would be the accurate atomistic modeling of solvent mixtures (binary or high order) that are excessively used as organic solutions. Future efforts should also focus on identifying co-solvents [200] that adequately reduce ion aggregation while simultaneously facilitating rapid metal-ion migration in the electrolyte phase to improve the performance of the battery. Even modest improvements in the transference number, e.g., to t+ ≈ 0.6-0.7, would be beneficial by allowing an increase in the energy density and high charge rates.