Understanding the formation of the solid–electrolyte interphase (SEI) in lithium-ion batteries is an ongoing area of research due to its high degree of complexity and the difficulties encountered by experimental studies. Herein, we investigate the initial stage of SEI growth, the reduction reaction of ethylene carbonate (EC), from both a thermodynamic and a kinetic approach with theory and molecular simulations. We employed both the potential distribution theorem and the Solvation Method based on Density (SMD) to EC solvation for the estimation of reduction potentials of Li+, EC, and Li+-solvating EC (s-EC) as well as reduction rate constants of EC and s-EC. We find that solvation effects greatly influence these quantities of interest, particularly the Li+/Li reference electrode potential in EC solvent. Furthermore, we also compute the inner- and outer-sphere reorganization energies for both EC and s-EC at the interface of liquid EC and a hydroxyl-terminated graphite surface, where total reorganization energies are predicted to be 76.6 and 88.9 kcal/mol, respectively. With the computed reorganization energies, we estimate reduction rate constants across a range of overpotentials and show that EC has a larger electron transfer rate constant than s-EC at equilibrium, despite s-EC being more thermodynamically favorable. Overall, this manuscript demonstrates how ion solvation effects largely govern the prediction of reduction potentials and electron transfer rate constants at the electrode–electrolyte interface.
I. INTRODUCTION
Lithium-ion batteries (LIBs) play an integral role in today’s technologies due to their high energy and power densities, making them an optimal choice for a large variety of applications, ranging from small handheld devices up to electric vehicles and even grid-scale storage. However, despite the wide array of use-cases for LIBs, they are not without their drawbacks. Among other downsides, LIBs suffer from capacity fade or gradual loss of total storage capacity over time due to many factors: pulverization of electrodes, separation from current collectors, dissolution of active material, formation of the solid–electrolyte interphase (SEI), etc. Among these many modes of capacity fade, it is seen that the formation of this solid–electrolyte interphase (SEI) layer is a leading contributor.1 Thus, it is important to understand the mechanisms by which this degradation process occurs such that it can be controlled.
While the formation of this SEI layer can result in capacity fade, it is also important to understand that the SEI protects the anode from solvent intercalation. Without it, the graphitic anode would delaminate, giving way to rapid failure of the battery.2 Furthermore, the SEI can form a passivating layer, thereby preventing its continual growth. However, controlling the chemistry behind this growth process is difficult because the governing chemistry is not well understood. Thus far, it is known that the SEI growth process is driven by the reductive decomposition of the battery’s electrolyte. There are many studies in the literature that explore the possible reaction mechanisms that follow from electrolyte reduction,3–10 but there are limited studies of formulating this process from the vantage point of electron transfer (ET) rate theory in the condensed phase.11–13 Notably, Kim et al. approached a similar problem but focused on Li+ reduction at a platinum electrode interface with a solid polymer electrolyte.13 The method they employed allowed for the evaluation of reorganization energies as a function of applied electrode potential and found that there was a negligible change in reorganization energy as potential was increased.
The complexity of electron transfer in the vicinity of interfaces (charged or uncharged) in which the interface can act as either a donor or an acceptor is an active area of research.7,11,14–17 Viewing interfacial processes through the lens of Marcus theory provides a potentially powerful theoretical framework to isolate the important role of local solvent fluctuations as a descriptor for a quantitative understanding of the solvent reduction. One aspect of this work is focused on generating a methodology for computing the rate of electron transfer (ET) for electrolyte reduction at a graphitic anode interface [shown in Fig. 1(a)] within the formalism of Marcus.18–22 We focus our attention on two electron acceptors: ethylene carbonate (EC) (in pure EC solution) and s-EC [EC in the solvation shell of a Li+, Fig. 1(b)].
(a) Snapshot of the hydroxyl-terminated graphite surface in contact with liquid EC. (b) Molecular structure of EC and Li+. Side views of the relaxed geometries of (c) EC and (d) EC−, highlighting inner-sphere reorganization. Carbon, oxygen, hydrogen, and lithium are denoted by gray, red, white, and green atoms, respectively.
(a) Snapshot of the hydroxyl-terminated graphite surface in contact with liquid EC. (b) Molecular structure of EC and Li+. Side views of the relaxed geometries of (c) EC and (d) EC−, highlighting inner-sphere reorganization. Carbon, oxygen, hydrogen, and lithium are denoted by gray, red, white, and green atoms, respectively.
Briefly, Marcus theory describes the rate of ET through interactions with solvent degrees of freedom in Marcus’s equation,
where VAB is the electronic coupling strength, λ is the reorganization energy, β = 1/kBT, and ΔA‡ is defined as
where ΔArxn is the free energy of reaction. Given that most interfacial reduction reactions at an electrode surface likely occur within the adiabatic regime of ET,23 we do not compute electronic coupling values.
The free energy of reaction is directly related to the reduction potential, which is a key quantity that describes the electrochemical stability of electrolytes,
where ΔA°rxn is the standard free energy of reaction, n is the number of electrons being transferred, F is Faraday’s constant, and Φ°abs is the standard reduction potential in the absolute scale (note that in this manuscript, the subscript of reduction potentials denotes the reference and subsequent parentheses contain the redox couple of interest). The inverse relationship in Eq. (3) implies that higher Φ° values correspond to a greater propensity for reduction. It is often valuable to report reduction potentials relative to a standard reference (Li+/Li in our case) rather than the absolute scale, given the difficulty of measuring absolute free energies experimentally. This is often achieved by using the standard hydrogen electrode (SHE) potential of 4.42 V in the absolute scale along with the Li+/Li reference potential of −3.05 V vs SHE,6,24–26
This conversion by subtracting 1.4 V is commonly used when computing reduction potentials for LIB electrolyte species, such as EC, dimethyl carbonate (DMC), vinylene carbonate (VC), and fluoroethylene carbonate (FEC).24–26 However, we discuss later in this paper the limitations of the conversions used in Eq. (4) for the estimation of reduction potentials.
In both thermodynamic and kinetic descriptions of EC reduction, the solvation free energy plays an important role that must be carefully addressed for the accurate prediction of reduction potentials and ET rate constants. The chemical potential of a species X in a solvent is represented by
where is the single-ion solvation free energy at a fixed point, ΛX is the deBroglie wavelength, γX is the activity coefficient, and ρX is the concentration of X. The first term in Eq. (5) describes the local interactions of solute X with a solvent, whereas the second term provides a statistical definition of the standard states used in this study.27
The single-ion solvation free energy is computed in a classical potential by
where UXS is the solute–solvent interaction energy and is defined as UXS = UXS,Ns − UNs, which is the energy difference for a given solvent structure with and without the ion. is the energy of X in vacuum, and ⟨⋯⟩0 denotes that averaging is done in the ensemble with no solute–solvent interactions. Equation (6) can be simplified into the following form using the potential distribution theorem (PDT):28–30
where is the free energy required to form a cavity with a radius matching species X and is the charging free energy that describes the solvent’s response to the creation of a solute’s charge. Equation (7) would also include contributions from the surface potential, but in this work, those contributions were computed to be negligible. A more in-depth derivation of can be found in Refs. 29 and 30. In the solvent model based on density (SMD)31, the solvation free energy can be approximated directly by taking the difference in energy between each state with and without a solvent model,
The charging free energy can be computed by incrementally adding small amounts of charge to the ion,29,30
where M is the total number of windows for charging; Uiδq and U(i+1)δq are the average potential energies of the system in windows i and i + 1, respectively; and δq is the incremental charge step size (δq = q/M). To ensure that the step size is small enough, can be computed in reverse via
and should match the profile generated by Eq. (9).
The focus of our work is on identifying the impact of solvation on the thermodynamics and kinetics of ethylene carbonate reduction in the condensed phase by computing reduction potentials for Li+, EC and s-EC, and ET rate constants for EC and s-EC. Multiple early works have studied ET at an interface23,32–34 and provide a good foundation for our study. We closely followed the work of Rose and Benjamin32 but with an increased complexity of the systems studied. In this paper, we use the potential distribution theorem (PDT) to evaluate solvation for condensed phase calculations of reduction potentials and compared to both an implicit solvent model (SMD) and literature values. We also compute reorganization energies for EC and s-EC reduction at a simple electrode interface for the prediction of ET rate constants across a range of overpotentials.
II. METHODS
Two main quantities are computed in this work for Li+, EC, and s-EC: reduction potentials and electron transfer (ET) rate constants. For the calculation of reduction potentials, several steps are involved (shown in Fig. 2), including the calculation of solvation free energies and free energies of reaction. The free energies of reaction are extracted from thermodynamic cycles (Fig. 3), which are converted to reduction potentials via Eq. (3). The key quantity of focus in this work is the solvation free energy, which largely affects the free energy of reaction and consequently the reduction potential. The ET rate constant calculations are governed by reorganization energies, which are discussed in Sec. II D.
Workflow for computing EC reduction potentials from solvation free energies. The same general workflow is also applied when computing Li+ and s-EC reduction potentials. This highlights the connection between solvation free energies and predicted reduction potentials.
Workflow for computing EC reduction potentials from solvation free energies. The same general workflow is also applied when computing Li+ and s-EC reduction potentials. This highlights the connection between solvation free energies and predicted reduction potentials.
Thermodynamic cycles for (a) Li+ + e− → Li, (b) EC + e− → EC−, and (c) Li+⋯EC + e− → Li+⋯EC− redox reactions. μ*(X) denotes the single ion solvation free energy in EC, and ΔAsub° denotes the standard free energy of sublimation. EIP and EEA denote the ionization potential and electron affinity in vacuum, respectively. X(solv), X(g), X(s) denote solvated, gaseous, and solid phases, respectively, for species X. ΔAb°(X, Y) denotes the standard binding energy between X and Y. The signs of thermodynamic quantities in the cycles arise from the transitions depicted by the arrows.
Thermodynamic cycles for (a) Li+ + e− → Li, (b) EC + e− → EC−, and (c) Li+⋯EC + e− → Li+⋯EC− redox reactions. μ*(X) denotes the single ion solvation free energy in EC, and ΔAsub° denotes the standard free energy of sublimation. EIP and EEA denote the ionization potential and electron affinity in vacuum, respectively. X(solv), X(g), X(s) denote solvated, gaseous, and solid phases, respectively, for species X. ΔAb°(X, Y) denotes the standard binding energy between X and Y. The signs of thermodynamic quantities in the cycles arise from the transitions depicted by the arrows.
A. Classical molecular dynamics (MD)
All classical simulations were run using GROMACS35 in the NVT ensemble at 300 K and a step size of 0.5 fs. The Nosé–Hoover thermostat36,37 was used to maintain the system’s temperature and the all-atom optimized potentials for liquid simulations (OPLS-AA) forcefield38–40 was used to describe all interatomic interactions. The OPLS-AA forcefield has been used previously to accurately compute transfer free energies for similar systems.41 All restraints (unless otherwise specified) were applied via the PLUMED library.42
B. Potential distribution theorem (PDT)
Solvation free energies are computed from condensed phase classical MD simulations using the potential distribution theorem (PDT),28–30 which breaks the solvation free energy into multiple contributions that can be computed separately, shown in Eq. (7). Charging free energy simulations were done in a 3 nm cubic box for the Li and EC systems solvated in bulk EC with no interface following Eq. (9) in increments of 0.1 |e|. Simulations of cavity formation were performed in the same conditions as charging, but instead of adjusting charge, the σ Lennard-Jones parameter of a neutral Li atom was incremented from 0 to 2.2 Å in 50 steps of 0.044 Å.
C. Reduction potentials
Quantum density functional theory (DFT) calculations for reduction potentials were performed with Gaussian 16 Rev. A.03.43 The B3PW91 hybrid functional44 was used with the 6-311++G(d,p) basis set45,46 for the calculation of reduction potentials, electron affinities, and ionization potentials. When an implicit solvent was utilized, the SMD method31 was employed with ɛ = 89.78 to match that of liquid EC.47 Free energies computed with DFT were obtained via frequency calculations, which compute entropic contributions and zero-point energy corrections to enthalpy at a reference state (1 atm and 298.15 K) using the rigid rotor harmonic oscillator (RRHO) approximation.
Binding energies were computed via umbrella sampling48 in the same simulation environment as charging using classical MD. The binding energies in Fig. 3(c) are derived from the following equation:
where Kb is the equilibrium binding constant. Kb is computed from a potential of mean force (PMF), w(r), via r2-weighted integration,
where r* denotes a radial distance within the bulk region and C° is the standard concentration (1 mol/l). We have defined the binding site region to begin at r = 0 and end at the first peak in the PMF.49
D. Reorganization energies
Reorganization energies for EC and s-EC are computed with respect to their reduced, cyclic structure (i.e., not the ring-opened structure). Although this ring-opened structure of EC− is lower in energy than the cyclic structure, there still exists an energy barrier separating these two states.50 Therefore, we have chosen to compute reorganization energies with respect to the cyclic structure of EC−.
Inner-sphere reorganization energies were calculated by first optimizing EC and EC− structures with q = 0 and −1 |e| net charges, respectively, with static DFT in CP2K.51 The Perdew–Burke–Ernzerhof exchange–correlation functional52 was employed with double zeta basis sets and GTH pseudopotentials53–56 along with a 400 Ry cutoff. Singlet and doublet multiplicities were used for the neutral and reduced states, respectively. Optimized structures for EC and EC− were then used to evaluate their energies with swapped charges (q = −1 and 0 |e|, respectively). The energy difference between the bent and flat configurations for each charge state was calculated and averaged to yield the inner-sphere reorganization energy in Table III. Alternatively, we evaluated the inner-sphere reorganization energy using ab initio molecular dynamics (AIMD) simulations of liquid EC in CP2K. AIMD simulations were run with a 0.5 fs time step and at the same level of theory as the inner-sphere reorganization energy calculations. EC− structures were sampled in liquid EC and then substituted into a separate AIMD simulation of pure liquid EC (no net charge) by choosing a single EC molecule and aligning their structures by minimizing the root-mean-square deviation (RMSD). The trajectory with the substituted structure was resampled (still with 0 net charge) and the resulting energy difference was averaged to give the AIMD value in Table III.
Marcus parabolas were constructed with a hydroxyl-terminated graphite surface contacting liquid EC and box size of 7.10 × 2.93 × 2.60 nm3 in GROMACS. EC and s-EC were held at the electrode–electrolyte interface (EEI) using a harmonic restraint. Regions beyond the free energy minima were sampled by directly adjusting the oxidant’s charge, which was then reweighted back to the two charge states (oxidized and reduced) in a manner similar to umbrella sampling (also discussed in Ref. 32), allowing us to reconstruct the full parabolas for high energy regions. Partial charges for the reduced state of EC were estimated via a restrained electrostatic potential (RESP) calculation in CP2K with the same level of theory used for inner-sphere reorganization energy calculations.
The reaction coordinate used in the Marcus parabolas was the electrostatic potential on the geometric center of the ion,32,57,58
where denotes the configuration of all atoms in the system, ke is Coulomb’s constant, qi is the charge on solvent atom i, and is the distance between atom i and the geometric center of the EC.
The order parameter in Eq. (13) is not easily restrained for umbrella sampling windows. Instead, we used the molecular net charge as an effective restraint along the order parameter. In this way, we can simply take the static bias potential as
where δq is the change in the net charge of EC, is the potential energy for the m-th window with net charge δq, and Uq=0,−1 is the potential energy for the same molecular configurations but with q = 0 or −1 |e|. This allows us to reconstruct both Marcus parabolas for the reactant and product diabatic states from the same series of simulations, in which the net charge of an EC molecule is changed from 0 to −1 in 0.2 |e| charge increments. The free energy, A, along the solvent coordinate for window m is then computed by
where δ[⋯] is the Dirac delta function, u is the reaction coordinate along which the free energy profile is established, and ⟨⋯⟩δq denotes that averaging is done in the ensemble with δq partial charges on EC. This process produces a PMF for each window along the order parameter following the procedure of umbrella sampling.48 Each PMF segment is manually aligned, so they create a single, smooth PMF. This manual alignment was allowable because the shape of the full PMF was largely expected to follow a parabolic trend, which was ultimately observed. The PMFs are then fitted to symmetrical parabolas for the reactant and product states.
III. RESULTS AND DISCUSSION
A. Solvation free energy
The free energy of reaction describes the driving force behind ET and can be approximated from the thermodynamic cycles depicted in Fig. 3. The free energy of ET reactions can also be represented as the reduction potential, which describes the likelihood of a species to be reduced. In Figs. 3(a) and 3(b), the oxidant is desolvated into the vacuum phase where the electron is added. Following reduction, the species is then deposited to the solid phase [Fig. 3(a)] or solvated again [Fig. 3(b)]. However, in Fig. 3(c), the Li+ and EC are coordinated and must be separated in the solvent phase first, after which the EC can follow the same cycle as Fig. 3(b) before finally re-coordinating with a Li+.
The key quantities in each of the thermodynamic cycles are the solvation free energies, which are represented as single-ion solvation free energies, . We posit that isolating this local27 contribution to the chemical potential will provide a useful descriptor for important collective effects.
In the case of EC, the full single-ion solvation free energy was not required, given that the difference in solvation energies between its reduced and neutral states was sufficient. The term in Eq. (7) cannot be ignored for EC because it is a polyatomic ion with a non-zero vacuum energy, whereas it can be neglected in the case of lithium. The difference in solvation free energy is then
Since the thermodynamic pathway for EC− solvation contains the EC solvation pathway, the cavity formation free energies are the same and cancel. denotes the charging from EC with partial charges in its neutral state to partial charges of EC−.
We have computed the charging free energy profiles, , for lithium and EC using Eq. (9), which are shown in Fig. 4(a). Lithium yields a large charging free energy at −106.7 kcal/mol, whereas EC yields a much smaller charging free energy of −21.7 kcal/mol. In Fig. 4(b), we compare the μ* and Δμ* values for Li+ and EC, respectively, when computed with SMD (lighter shade) and PDT (darker shade). For Li+, the single-ion solvation free energy was computed to be −89.3 and −104.8 kcal/mol for SMD and PDT, respectively. The PDT value was computed using Eq. (7) where was taken from Fig. 5 to be 1.9 kcal/mol. This free energy profile of cavity formation in liquid EC increases at a rate approximately twice as fast as that in water,59 which is supported by the difference in their densities (1.32 g/cm3 vs 0.997 g/cm3). The SMD value for was computed using Eq. (8) and reports a significantly lower value than PDT, which can be attributed to the limited applicability of the SMD model for ionic species with such small cavities.25
(a) Charging free energy profiles of EC and Li. From left to right, the charge of EC transitions from 0 to −1, whereas Li transitions from 0 to +1. (b) Comparison of charging free energies computed with SMD (lighter shades) and PDT (darker shades).
(a) Charging free energy profiles of EC and Li. From left to right, the charge of EC transitions from 0 to −1, whereas Li transitions from 0 to +1. (b) Comparison of charging free energies computed with SMD (lighter shades) and PDT (darker shades).
Free energy of cavity formation for a hard sphere in ethylene carbonate up to the radius of Li+.
Free energy of cavity formation for a hard sphere in ethylene carbonate up to the radius of Li+.
There is a 15.5 kcal/mol difference between the two solvation methods, but our computed value for the single-ion solvation free energy of lithium in EC is supported by quasi-chemical theory (QCT) calculations that also predict a +15 kcal/mol transfer free energy for Li+ between EC and water, with water being the favored solvent.41 Comparing our value to the free energy of solvation for Li+ in water (−120 kcal/mol),30 we see very close agreement in the transfer free energy from Ref. 41. Again, we note that these comparisons only apply to the single-ion solvation free energy and do not account for chemical activities in Eq. (5).
For EC, SMD predicts that kcal/mol from Eq. (8), whereas PDT predicts a value of −26.1 kcal/mol. The term of Eq. (16) was computed to be −4.4 kcal/mol. SMD predicts a more favorable solvation free energy for EC over PDT by ∼10 kcal/mol.
It is important to note that these reduction reactions occur near the electrode interface, but the single-ion solvation free energies computed in this work are for the bulk phase. An attempt to adjust the free energy from the bulk phase to the interfacial region did not yield noticeable changes, which was due largely to non-polarizability of the surface used. An improved description of the interface with a polarizable surface would yield a better estimate of the interfacial solvation region. The charging of species directly at the interface would not be required since bulk phase solvation energies can be adjusted by the difference in interfacial binding energies. Recently, Angarita-Gomez and Balbuena have demonstrated the destabilization of the solvation shell of a Li+ as it approaches an anode interface, with changes in the solvation free energy of ∼1 eV for a Li+ in EC.60 Their work further highlights the impact of interfacial effects on solvation energies.
B. Reduction potentials
In Tables I and II, we have used the single-ion solvation free energies in Fig. 4(b) to compute reduction potentials for Li+, EC, and s-EC. The values for −EIP(Li) and −EEA(EC) in Fig. 3 were computed to be −128.6 and 7.1 kcal/mol, respectively. The two binding energies in Fig. 3(c), −ΔAb°(Li+, EC) and , were computed to be −0.5 and −5.7 kcal/mol, respectively. The handbook61 value for Li sublimation, ΔA°sub(Li), is 38.1 kcal/mol.
Reduction potentials for Li+, EC, and s-EC with solvation free energies computed via PDT. Values in the first column are reported in volts in the absolute scale, whereas all other columns are reported in volts vs Li+/Li. The column corresponds to the Li+/Li scale as computed in this work, whereas the column represents reduction potentials in the Li+/Li scale using the commonly applied relationship . Columns with an additional 1.3 V represent the inclusion of the ring-opening reaction as discussed in the text. Although the last column closely matches the experimental value (0.8 V vs Li+/Li),62 the origin of this lacks the self-consistency of using the Li+/Li scale as computed in this work as discussed in the main text.
Species . | Φabs° . | . | . | . | . |
---|---|---|---|---|---|
Li+ | 2.60 | 0.00 | ⋯ | ⋯ | ⋯ |
EC | 0.82 | −1.78 | −0.48 | −0.58 | 0.72 |
s-EC | 1.10 | −1.50 | −0.20 | −0.30 | 1.00 |
Species . | Φabs° . | . | . | . | . |
---|---|---|---|---|---|
Li+ | 2.60 | 0.00 | ⋯ | ⋯ | ⋯ |
EC | 0.82 | −1.78 | −0.48 | −0.58 | 0.72 |
s-EC | 1.10 | −1.50 | −0.20 | −0.30 | 1.00 |
Reduction potentials for Li+, EC, and s-EC with solvation free energies computed via the SMD implicit solvent model. Values in the first column are reported in volts in the absolute scale, whereas all other columns are reported in volts vs Li+/Li. The column corresponds to the Li+/Li scale as computed in this work, whereas the column represents reduction potentials in the Li+/Li scale using the commonly applied relationship . Columns with an additional 1.3 V represent the inclusion of the ring-opening reaction as discussed in the text.
Species . | Φabs° . | . | . | . | . |
---|---|---|---|---|---|
Li+ | 3.35 | 0.00 | ⋯ | ⋯ | ⋯ |
EC | 1.28 | −2.07 | −0.77 | −0.12 | 1.18 |
s-EC | 1.56 | −1.79 | −0.49 | 0.16 | 1.46 |
Species . | Φabs° . | . | . | . | . |
---|---|---|---|---|---|
Li+ | 3.35 | 0.00 | ⋯ | ⋯ | ⋯ |
EC | 1.28 | −2.07 | −0.77 | −0.12 | 1.18 |
s-EC | 1.56 | −1.79 | −0.49 | 0.16 | 1.46 |
For comparison, we used both implicit (SMD) and explicit (PDT) solvent models to compute the charging free energy. For each respective calculation method, we first computed the Li+/Li reduction potential in the absolute scale, Φabs°(Li+/Li) (reported in Tables I and II), which were then used to convert the reduction potentials of EC and s-EC from the absolute scale to the Li+/Li scale via
where ΔA°rxn is the standard free energy of reaction, n is the number of electrons transferred, and F is Faraday’s constant.
The reported value for Φabs°(Li+/Li) ranges from 1.3 to 1.79 V, with 1.4 V being commonly used as it aligns with the conversion described in Eq. (4).26 However, this range of reported values do not correspond to Li+ solvation in EC but rather with water, DMSO, or acetonitrile. While the solvent is expected to only shift the reduction potential by a couple tenths of a volt,26 we computed Φabs°(Li+/Li) to be 3.35 and 2.60 V with SMD and PDT, respectively. In Tables I and II, we have included additional columns that correspond to the commonly used 1.4 V value and/or the addition of 1.3 V, which corresponds to the inclusion of the ring-opening reaction of EC (discussed later in this section). The additional columns do not include values for Li+ because the new conversions/modifications do not apply to Li+ reduction. These columns are included to demonstrate how we can align with literature values when using similar conversion methods from the absolute scale. However, a major focus of this study is to compute reduction potentials self-consistently [i.e., using Φabs°(Li+/Li) computed in this work, along with Eq. (17)].
Importantly, the solvation free energy of Li+ is seemingly the only quantity that affects Φabs°(Li+/Li) in Fig. 3(a) when changing solvents. The value of 1.4 V that is often used6,24–26 corresponds to a solvation free energy of ∼−130 kcal/mol for Li+, which is even greater than the solvation free energy for water by 10 kcal/mol. However, Li+ is expected to have a less favorable solvation free energy in EC vs water by ∼15 kcal/mol (−105 vs −120 kcal/mol).41 The resulting 25 kcal/mol difference in between our computed value (−105 kcal/mol) and the purported −130 kcal/mol leads to a 1.08 V difference in Φabs°(Li+/Li) values.
Additionally, the common 1.4 V corresponds to aqueous systems and thus is expected to differ for non-aqueous systems such as EC. Busch et al. presents an approach for computing the standard hydrogen electrode potential for aqueous and non-aqueous and reports a Φabs°(SHE) value of ∼4.1 V in water and 4.9 V in polar, non-aqueous solvents (e.g., DMSO, acetonitrile, and DMF).63 This apparent 0.8 V difference demonstrates the disparity between aqueous and non-aqueous electrode potentials.
Comparing the reduction potentials for EC and s-EC when converted via 1.4 V (“” columns in Tables I and II) to those in the literature,6,24 we see a difference of 0.5–0.8 V for the PDT-derived values and 0.1–0.2 V difference for the SMD-derived values. The deviation in PDT-derived values is largely due to the difference in methods for computing solvation free energies. The close agreement with the SMD-derived values is unsurprising given that approaches in the literature often employ implicit solvent models.6,24 These reduction potentials are included to allow for the direct comparison with the literature, but we maintain the importance of computing the reduction potentials self-consistently.
Compared to PDT, SMD reports a lower reduction potential for the Li+/Li reference due to its less favorable solvation free energy but predicts larger Φ°abs values for EC and s-EC. The primary difference between EC and s-EC is the difference in binding energies for the s-EC cycle in Fig. 3(c), and thus, ion concentration would likely impact Φ°(s-EC) values through the binding energies computed with Eqs. (11) and (12).
An important aspect to consider is the overall reaction that is measured experimentally. Zhang et al. previously measured the reduction potential of EC in THF to be 1.36 V vs Li+/Li, which matched closely with their computed value of 1.46 V vs Li+/Li.64 However, the putative reaction in their thermodynamic cycle also included the ring-opening degradation reaction of EC−, whereas we are only considering the ET reaction. An inclusion of the ring-opening reaction would increase the s-EC and EC reduction potentials by ∼1.3 V (shown in Tables I and II) if the ring-opening free energy of reaction in solvent is −30 kcal/mol.50 The inclusion of the ring-opening reaction does bring our computed reduction potentials (“” columns in Tables I and II) for EC and s-EC closer to literature values, but they still differ by ∼1 V from the experimental value of 0.8 V vs Li+/Li.62 However, compared to theoretical values computed in the literature,6,24 our PDT-derived values only differ by ∼0.5 V. This discrepancy is further lessened by using the common 1.4 V value along with the ring-opening reaction (“” columns). However, the goal of this study is to present a workflow where reduction potentials can be computed self-consistently, i.e., using a computed value for Φabs°(Li+/Li) to convert to the Li+/Li scale. Although the rightmost column in Table I closely matches the experimental value for EC reduction (0.8 V vs Li+/Li), we argue that this result is confounded by its lack of self-consistency in its approach.
Another possible contribution to this discrepancy is related to the limited knowledge of the (electro)chemical environment at the electrode–electrolyte interface (EEI). Electrolyte reduction is expected to occur at the EEI within the electric double layer,6 wherein the structuring of dipoles and ions is highly ordered and electric potentials are not fully screened, which can greatly impact the solvation free energy. The thermodynamic cycles in Fig. 3 do not discriminate between interfacial and bulk solvation and thus introduce a potential source of error that is not easily resolved. Baskin and Prendergast have developed continuum models to describe the speciation at such interfaces and would prove useful in future studies.65 Another useful tool would be the constant potential method (CPM),66 which allows for a constant potential, polarizable electrode that can provide an accurate description of the chemical environment at the EEI.
C. Reorganization energy
With the free energies of reaction computed, we shift our focus to another major factor that governs ET rates, the reorganization energy, which describes the entire system’s response to a transfer of charge. This reorganization energy, λ, can be decomposed into two components,
where λi and λo are the inner- and outer-sphere reorganization energies, respectively. Inner-sphere refers to the internal reconfiguration associated with a change in charge, e.g., extension of a bond or loss of planarity. Outer-sphere reorganization refers to the solvent’s response to the change in the charge state of the oxidant.
EC exhibits both an inner- and outer-sphere response to charge transfer. Once reduced, the sp2 carbon in the carbonate group bends out of plane to account for the newly acquired electron. The inner-sphere reorganization energy is the energy difference between the bent and flat conformations for a given charge state, which we have computed to be 46.6 and 35.8 kcal/mol when EC’s charge is 0 and −1, respectively. The average of the two quantities gives us an inner-sphere reorganization value of 41.2 kcal/mol, which closely matches what has been computed in the literature previously.11 The inner-sphere reorganization energy was also computed with AIMD calculations in the condensed phase to produce a similar value of 43 kcal/mol. Given that there is little difference between the two computed λi values, we conclude that there is a limited impact on inner-sphere response by the solvent degrees of freedom. All inner-sphere reorganization energy values are shown in Table III. The inner-sphere reorganization energy is applied to both EC and s-EC in this work.
PBE/DZVP energies of bent and flat configurations of EC with net charges of 0 and −1 [relative to flat, neutral EC, i.e., Fig. 1(c)]. Inner-sphere reorganization energies for each charge state are shown on the right. All energies are in kcal/mol.
Charge (|e|) . | Flat . | Bent . | λi . |
---|---|---|---|
0 | 0.0 | 46.6 | 46.6 |
−1 | 39.9 | 4.1 | 35.8 |
Mean λi,DFT | 41.2 | ||
Mean λi,AIMD | 43 |
Charge (|e|) . | Flat . | Bent . | λi . |
---|---|---|---|
0 | 0.0 | 46.6 | 46.6 |
−1 | 39.9 | 4.1 | 35.8 |
Mean λi,DFT | 41.2 | ||
Mean λi,AIMD | 43 |
Outer-sphere reorganization energies, λo, were computed by constructing the full Marcus parabolas along the solvent coordinate, shown in Fig. 6. EC and s-EC were held directly at the interface during the construction of the Marcus parabolas. In the classical potential, no new force-field parameters were created to account for the EC− bent structure [Fig. 1(d)] and instead neutral EC parameters were used. Thus, the reorganization energies from Fig. 6 represent only outer-sphere response. The hydroxyl-terminated graphite surface [Fig. 1(a)] that was used in this study is not polarizable; however, we are interested in computing reorganization energies at the interface, which have been shown to remain stable across a range of voltages.13
Marcus parabolas along the solvent coordinate for (a) EC and (b) s-EC. The black and red curves in (a) and (b) denote the reactant and product states, respectively, with outer-sphere reorganization energies labeled. Parabolas are shown for the case of , which also corresponds to η = 0. (c) Relationship between the predicted rate constant for electron transfer and overpotential, η, for EC and s-EC.
Marcus parabolas along the solvent coordinate for (a) EC and (b) s-EC. The black and red curves in (a) and (b) denote the reactant and product states, respectively, with outer-sphere reorganization energies labeled. Parabolas are shown for the case of , which also corresponds to η = 0. (c) Relationship between the predicted rate constant for electron transfer and overpotential, η, for EC and s-EC.
For both EC and s-EC, we see large contributions from the solvent (35.4 and 47.7 kcal/mol, respectively) that go into the full reorganization energy. Leung previously employed an implicit solvent model to estimate the full s-EC reorganization energy at 47.3 kcal/mol, which is nearly the same amount as the value computed in this work for purely the outer-sphere contribution. Overall, the full reorganization energies for EC and s-EC are estimated at 76.6 and 88.9 kcal/mol, respectively. This result demonstrates the strength of the solvent’s response to EC reduction and highlights the importance of explicit solvation.
D. Prediction of electron transfer rate constants
With both free energies of reaction and full reorganization energies computed for EC and s-EC, the last remaining unknown quantity in Eq. (1) is the electronic coupling, VAB. However, it is expected that the reduction of EC during SEI growth occurs in the adiabatic regime of ET, implying that the electronic coupling is high. To approximate this, we assume that |VAB| = 100 meV. We also tested a range of values and found that there was a negligible impact on the trends reported herein. Although we strived for self-consistency to obtain values for the ET rate constant within the molecular framework presented in our manuscript, computing the diabatic couplings remains a challenge. Currently, we subscribe to the theory of Ref. 67 for the condensed phase. However, further modifications of this theory have been proposed in Ref. 68 and will be the focus of a future study.
With each quantity computed or defined, we can introduce a dependence on an overpotential that adjusts the system’s free energy of reaction,32
where η is the overpotential and Δq is the total change in charge (i.e., −1).32 We estimated the ET rate constant as a function of overpotential, which is shown in Fig. 6(c). An overpotential of 0 is assigned for the case of .32
For both EC and s-EC, there is an apparent peak in the rate constants of ET between 3 and 4 V in overpotential, wherein the rate constants decrease as driving force increases due to the inverted region. However, systems do not typically operate at such large overpotentials and the inverted region is not expected to be reached. EC seemingly exhibits a larger ET rate constant than s-EC for overpotentials up to the inverted region, but both systems have been scaled to η = 0 when . Since the overpotential describes the deviation from theoretical reduction potential values, Fig. 6(c) shows that EC has a slightly larger ET rate constant than s-EC when at equilibrium, but s-EC is still more thermodynamically favorable due to its higher reduction potentials in Tables I and II. The primary reason for EC having a higher ET rate constant at equilibrium is due to its lower reorganization energy over s-EC because in the limit of ΔArxn = 0, ΔA‡ = λ/4.
If we assume that η represents the deviation from the theoretical reduction potential,69 then setting the η = 0 locations to the respective Φ° for each curve would demonstrate that ket for EC and s-EC are nearly overlapping in the region of their reduction potentials, demonstrating the competition between kinetics and thermodynamics for this system. Kim et al. also demonstrated that changes in electrode potential had a minimal effect on the reorganization energy using a constant potential method for a model platinum metal electrode.13 Although this approach would likely impact the absolute value of outer-sphere reorganization energies, the trend between EC and s-EC would not be expected to change.
IV. CONCLUSIONS
In conclusion, we have utilized the potential distribution theorem to describe solvation and compared it to the SMD implicit solvent model in DFT in the prediction of reduction potentials for Li+, EC, and s-EC. Our results highlight the impact of solvation free energy regarding the disparity between the commonly used 1.4 V for Φabs°(Li+/Li) and our computed values, which were lower by 1–2 V. Moreover, we have established a workflow for the prediction of rate constants of electron transfer in heterogeneous reduction reactions that require treatment of both inner- and outer-sphere reorganization. We find that s-EC is more likely to be reduced over pure EC due to their reduction potentials. However, at equilibrium, EC has a faster ET rate constant than s-EC due to its lower reorganization energy. This implies that the kinetics of EC reduction may compete with the favorable thermodynamics of s-EC reduction.
Although we argue the importance of explicit solvation, we acknowledge that the approach in this work is not without limitations. Equation (7) is missing a quantum mechanical contribution term, whereas real solvation free energies should account for quantum mechanical interactions according to Ref. 30. Therefore, the difference in real solvation free energies () would also include the difference in quantum mechanical contributions to solvation.
Exciting future works involve an improved description of the chemical environment at the EEI to better mimic an electric double layer with the use of a constant potential, polarizable electrode.66 In conjunction with more expensive DFT calculations of ion solvation free energies, the constant potential approach would facilitate improved estimations of reduction potentials beyond what are provided in the commonly used thermodynamic cycles in Fig. 3, in addition to better estimates of reorganization energies in the electric double layer. Furthermore, reduced models such as a generalized Langevin equation can also be constructed to estimate electron transfer rates with reactive flux calculations. Experimentally, this work can be supported by measuring Li+ ion solvation free energy in EC for direct comparison. The measurement of the Li+/Li reference electrode in EC solvent would also greatly benefit this study.
ACKNOWLEDGMENTS
This material is based upon work supported by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences, CPIMS Program (Award No. DE-SC0019483) and the U.S. Department of Energy, Office of Science, Office of Workforce Development for Teachers and Scientists, Office of Science Graduate Student Research (SCGSR) program. The SCGSR program is administered by the Oak Ridge Institute for Science and Education for the DOE under Contract No. DE-SC0014664. C.J.M. acknowledges funding support from the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences. This research used the resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231. PNNL is a multiprogram national laboratory operated for the DOE by Battelle under Contract No. DE-AC05-76RL01830. Computing resources were generously allocated by PNNL’s Institutional Computing program. L.D.G. would like to thank Tim Duignan and Eric Stuve for helpful discussions.The authors would also like to thank David Prendergast and Artem Baskin for helpful discussions.
AUTHOR DECLARATIONS
Conflict of Interest
The authors declare no competing financial interest.
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.