Interactions of CO2 Anion Radicals with Electrolyte Environments from First-Principles Simulations.

Successful transformation of carbon dioxide (CO2) into value-added products is of great interest, as it contributes in part to the circular carbon economy. Understanding chemical interactions that stabilize crucial reaction intermediates of CO2 is important, and in this contribution, we employ atom centered density matrix propagation (ADMP) molecular dynamics simulations to investigate interactions between CO2 - anion radicals with surrounding solvent molecules and electrolyte cations in both aqueous and nonaqueous environments. We show how different cations and solvents affect the stability of the CO2 - anion radical by examining its angle and distance to a coordinating cation in molecular dynamics simulations. We identify that the strength of CO2 - interactions can be tailored through choosing an appropriate cation and solvent combination. We anticipate that this fundamental understanding of cation/solvent interactions can facilitate the optimization of a chemical pathway that results from selective stabilization of a crucial reaction intermediate.

Introduction

Transforming CO2 into value-added chemicals or fuels is of great interest to create a sustainable carbon neutral cycle to tackle challenges such as climate change.[1−3] However, reducing CO2 is a challenging process due to its thermodynamic stability, poor electron affinity, and large kinetic overpotentials. Additionally, as a nonpolar gas, it is only sparingly soluble in water (0.033 M),[4] a common solvent for electrochemical reduction of CO2.[1] Therefore, it might be advantageous to use nonaqueous (aprotic) solvents due to the higher CO2 solubility, enhanced potential window, as well as lower proton concentrations that can suppress the undesired hydrogen evolution reaction.[5−7]

The first step of the CO2 reduction reaction (CO2RR) involves the formation of CO2–, and although successful studies have been performed to unravel the mechanism of electrochemical CO2RR,[8−17] molecular-level understanding of the dynamics of CO2– interacting with surrounding electrolyte cations and solvents (aqueous vs nonaqueous) remains unclear. Of particular interest is the cation effect,[18−24] which was first introduced by Murata and Hori, and it has been better understood that coupling alkali metal cations (M+) with CO2– cocatalyzes the first step[25,26] in aqueous solvents, i.e., * + CO2 + M+ + e– → M+···*CO2– by a short-range electrostatic interaction. In nonaqueous solvents, however, alkali metals are known to inhibit CO2 reduction due to the formation of a passivation layer containing carbonate species.[27,28] Therefore, electrolyte cations such as NX4+ (X = methyl, ethyl, and butyl groups) are being used, and their catalytic roles were examined.[24,29] In aprotic Li–CO2 batteries, it was argued that a LiCO2 intermediate is crucial to tuning CO2RR toward a solution- versus surface-mediated pathway.[30]

In this contribution, we aim to provide theoretical insights into the fundamental interactions between CO2–, obtained from the electrochemical reduction of CO2 gas, and the supporting electrolyte cations in various solvents, which will help direct the mechanism toward desired products such as CO, formic acid, and/or oxalate, and will lay a foundation to the future study of specific catalyst systems. Here, we focus on the stabilization effect of the CO2– anion radical, as its formation is likely to be the rate-determining step[31] for CO2 reduction in some aqueous[32−34] and nonaqueous[15,35] media, although this is still under debate.[1] We compare water and nonaqueous solvents (mainly dimethoxyethane (DME) in this work), as glyme-ether-based solvents have been successfully used for electrochemistry (e.g., metal-air/Mg battery) and gas separation processes capturing CO2[36,37] but have been rarely studied for CO2 reduction until recently,[27] to the best of our knowledge.

To gain understanding of the specific molecular-level interactions, we used ab initio molecular dynamics (AIMD) simulations.[38,39] These AIMD simulations show the atomic movement within the system, where the extremes and averages of the oscillations can give detailed information on the stability and other thermodynamic properties. We use the atom centered density matrix propagation (ADMP)[40−42] molecular dynamics method, as it gives the accuracy of density functional theory (DFT) and has advantages particularly around handling the dynamics of charged and radical molecular systems.[43−47] This approach allows us to gather accurate information on how the CO2– anion radical interacts with cation and solvent species, based on the computed descriptors such as coordination bond distances (cation–CO2) and bond angles to develop a fundamental understanding of the interactions. The shorter coordination distance and/or more decreased bond angle of CO2– indicate that the intermediate complex can be present for prolonged periods and likely undergo subsequent desired reactions such as chelation, charge transfer interactions, and bond cleavage/formation reactions.[48,49] We further construct a physics-based model to rationalize the relative contributions of electrostatic and covalent/noncovalent interactions to the stabilization effects of CO2–.

Computational Details

Initial geometry optimizations were performed using Gaussian 16 software[50] at the ωb97xD/6-31+G(d,p) level of theory[51,52] with ultrafine integration grids and multiple conformers considered.[53] The frequencies were calculated at the same level of theory. Upon geometry optimization, trajectories were calculated using ADMP molecular dynamics at the ωb97xD/jun-cc-pVDZ level of theory[51,54,55] for 5000 time points (500 fs). ADMP simulations were performed with and without the presence of implicit solvents[56,57] (water (H2O), tetrahydrofuran (THF), dimethyl sulfoxide (DMSO), acetonitrile (MeCN), and n,n-dimethylformamide (DMF)). Optimized, maximum, mean, and variation amplitude of the CO2– angle (θ) and the CO2– distance (to the cation, d) were extracted from the AIMD trajectory. As shown in Figure a, the CO2– angle (θ) is measured from one oxygen to the other, with the carbon atom as the angle vertex. The CO2–-to-cation distance (d) is measured from the carbon of the CO2– to the center of the cation (either the cation atom or the nitrogen in the ammonium cations). A typical ADMP trajectory is shown in Figure b, where CO2– oscillates around its equilibrium position over the course of the simulation. Additional details of the simulations are provided in the Supporting Information (SI).

Figure 1

(a) Schematic showing how the change in bond distance (Δd) or change in angle (Δθ) is calculated. Scheme (a) shows a molecular cluster with the dashed black line indicating the bond distance (measured Li+ to CO2– from the center of the cation to the carbon atom) and gray lines showing how the bond angle (θ) is measured (with the carbon as the vertex of the angle). The lithium, carbon, oxygen, and hydrogen atoms are shown in purple, dark gray, red, and light gray, respectively. (b) Computed trajectories of the bond distance (Δd) and CO2– angle (θ) over 500 fs for a Li+–CO2––3H2O molecular cluster. (c) Computed average Δd for various cationic complexes. All complexes in state (i) include three explicit solvent molecules. The dots indicate the difference in bond distances of optimized complexes with and without explicit solvent molecules. Me, Et, and Bu denote methyl, ethyl, and butyl groups, respectively.

Results and Discussion

The interaction strength of two charged particles is inversely related to the separating distance.[58] As shown in yellow bars in Figure c, based on the simulations, we observe that the interaction of dynamic water molecules increases the cation-to-CO2– distance (Δd) for all cations compared to the gas-phase geometry, which manifests less negative complexation energies and binding energies compared to those in the gas phase (ΔHgas, Table ). Note that the equations for the complexation enthalpy calculations are provided in the Supporting Information. In the case of the DME solvent system, the chelation of the DME molecule with the alkali metal cation increases the interaction distance for alkali metal cations with the CO2–, while for NX4+ cations (X = Me, Et, and Bu), the difference is minimal (<0.1 Å). This is largely due to the chemical differences between water and DME molecules. Water molecules are strongly polar and have strong ion-dipole interactions with cations and CO2–, as seen from the computed complexation enthalpies (ΔHcomplex) summarized in Table . DME is weakly polar and has chelating interactions with alkali metal cations, where the unshielded positive charge can interact closely with the lone pairs on the oxygen atoms of DME. On the basis of the computed enthalpies of complexation (Table ), the CO2– species itself weakly interacts with DME as compared to water (by 0.50 eV) because DME molecules do not have any strong localized partial positive charges. The NBu4+ positive charge is highly shielded by the butyl arms and as such does not closely interact with the DME oxygen atoms, which also can be seen from our quantum chemical calculations that the difference between ΔHcomplex and ΔHgas for the NBu4+ complex is much smaller than that for alkali metal cations (Table ). The optimized structures of NBu4+ complexes are shown in Figure S1. The optimized structures of Li+–CO2––solvent with DME and H2O also corroborate the above analysis; i.e., the oxygen atoms of CO2– only bind to Li+ in the DME solvent, as shown in Figure a, whereas in water one of the CO2– oxygen atoms binds to Li+, and the other interacts with water to form an O···H hydrogen bond as shown in Figure b. Overall, this suggests that the identity of the supporting electrolyte cation, as well as the polarity of solvents, both dictate the strength of interactions of the CO2– anion radical, which primarily consist of (electrostatic) ion-ion and ion-dipole interactions and hydrogen bonds. The different types of interactions will be discussed in more detail.

Table 1

Computed Complexation Enthalpies and Binding Enthalpies (ΔHcomplex and BECO2-, in eV) of CO2–, Li+, Na+, K+, and NBu4+ with Three Explicit H2O or DME Molecules

species	ΔHcomplex (3H2O)	BECO2- (3H2O)	ΔHcomplex (3DME)	BECO2- (3DME)	ΔHgas (CO2- + M+)
CO2–	–1.72		–1.22
Li+	–3.66	–5.08	–5.73	–3.59	–6.50
Na+	–2.73	–4.96	–4.65	–3.51	–5.66
K+	–2.03	–4.82	–3.55	–3.49	–5.00
NBu4+	–1.17	–3.65	–2.03	–3.28	–3.84

Figure 2

Gas-phase optimized structures of Li+–CO2––solvent with (a) DME and (b) H2O.

To further explore how the NX4+–CO2– interaction responds to the polarity of the solvent, we carried out ADMP simulations with three explicit solvent molecules (in THF, MeCN, DMSO, DMF, and water) and at the same time included the implicit solvent field of corresponding solvent models in Gaussian 16. For these bulky clusters, we repeated each simulation three times with different initial geometries (randomly generated and subsequently optimized) to make sure that the selection of initial structures have no role in dictating the final results. The average and maximum cation-to-CO2– distances are shown in Figure .

Figure 3

(a) Average cation-to-CO2– distance daverage and (b) maximum cation-to-CO2– distance dmax during simulations in various solvents, as plotted against the polarity (ET(30))[62] of solvent. These data points are for clusters of NBu4+, CO2–, three molecules of the solvent in question, and with an implicit solvent field of that solvent. Error bars represent standard deviations of three simulations.

We observe that the strength of interaction between NBu4+ and CO2– decreases with increasing solvent polarity, indicated by the increasing cation-to-CO2– distance, as shown in Figure . Such an observation is consistent with the dominant electrostatic interactions and has experimental implications on salt dissociation or formation of ion pairs.[25,27] The positive charge on the NBu4+ cation is highly shielded by the bulky and electron donating alkyl chains, which means that changing solvents has a relatively minimal effect on the NBu4+. However, the strength of the interaction between CO2– and NBu4+ is affected by how strongly CO2– is coordinated to the solvent and cationic species. Thus, increasing the polarity of the solvent increases CO2– coordination with the solvent, weakening the bond between CO2– and NBu4+. Similar qualitative trends in CO2RR were also reported by Berto et al.[59] and Shi et al.,[60] where voltammetric experiments in the presence of NBu4+ electrolytes have shown larger current density values for CO2RR in highly polar DMF and MeCN, as compared to low polarity THF. Here, DMSO has the largest standard deviation (Figure ), which was also found to exhibit the largest deviation from trends due to coordination via the sulfur atom (as opposed to via the oxygen atom in water and THF) in a study of aprotic solvent (THF, DMF, and DMSO) adsorption on metal surfaces by Nørskov and co-workers.[61] In aprotic solvents and NBu4+ salts, tailoring the lifetime of the CO2– will likely depend on solvent, due to the stabilizing influence of more polar solvents.

The foregoing analysis suggests that both solvents and cations play an important role in the stabilization of the CO2– radical,[59,60] but the precise effects are less clear. On the basis of our simulations, we hypothesize that inner-sphere solvation effects are the result of electrostatic and covalent/noncovalent interactions based around solvent stabilization and cation binding, and to this end, we sought to rationalize these by constructing a simple physics-based model. Compared to gas-phase CO2–, the introduction of a solvent molecule or a cation reduces the angle of CO2–. In Figure a, the largest reductions are observed with cations without solvent molecules, following the trend of Li+ > Na+ > K+ > NX4+ in the gas phase (gray bars), where X = methyl or longer alkyl groups. There is an overall trend of a larger change in the bond angle (Δθ) being associated with more negative complexation enthalpies, as shown in Figure b, which also indicates that increased complex stability is associated with a decreased CO2– angle. This decreased angle is due to anion radical stabilization through ion-ion interactions, ion-dipole interactions, and/or hydrogen bonds. Given the computed trends of solvent polarity and complexation enthalpy simultaneously (Figure and Figure b) and together with known experiment data, we would be able to provide guidance on the electrolyte environments, e.g., with larger stabilizing effects on CO2–[63] toward the desired binding.

Figure 4

(a) Bar chart of the change in the CO2– angle (Δθ) between the uncomplexed CO2– anion radical in the gas phase and the CO2– angle when complexed with cations and/or three solvent molecules. (b) Computed complexation enthalpy plotted against the change in angle from the gas phase to solvated (by DME or H2O).

To decompose the quantum chemistry-calculated complexation enthalpies into stabilization effects by electrostatic ion-ion interactions and/or ion-dipole interactions explicitly, we have constructed a simple physics-based model (see Supporting Information). This model takes in the geometry of a cluster and computes the different types of interactions based on fundamental physics equations. The model-predicted complexation energies against the known DFT-calculated complexation energies are shown in Figure , where red triangles (“ion-CO2–”, e.g., K+–CO2–) represent systems that have only electrostatic ion-ion interactions (also see ΔHgas in Table ), green squares (“ion-solvent”, e.g. NBu4+–2H2O) represent systems that have ion-dipole and dipole-dipole interactions, and purple circles (“ion-solvent-CO2–”, e.g., Li+–CO2––3DME) represent systems that have ion-ion, ion-dipole, and dipole-dipole interactions. The mean absolute error of the model-predicted complexation energy relative to the DFT-calculated value is about 0.12 eV per ion and solvent molecule in a cluster. For all the “ion-solvent-CO2–” clusters (purple circles), we calculate that on average 73% of the predicted complexation energy results from the electrostatic ion-ion interactions; likewise, the ion-dipole interaction is much stronger than the dipole-dipole interaction. This indicates that CO2 activation could be achieved through electrolyte stabilization of the CO2– anion product due to a stronger interaction with CO2– than the neutral and nonpolar CO2 gas.[64,65]

Figure 5

Physics-based model-predicted complexation energy against the known DFT-calculated energy. The dashed line indicates a perfect agreement between the model-predicted complexation energy and the known DFT-calculated complexation energy.

Finally, we show the average CO2– bond angles in DME and water as a function of the number of solvent molecules (Figure ) from our simulations. The average cation-to-CO2– distances are shown in Figure S2. Here, a small and computationally tractable number of solvent molecules could be used to simulate a reasonable solvation behavior,[66,67] where key interactions between cation/anion and solvent (ion-dipole) are captured. We note that there is a different behavior for Li+–CO2––nDME where the CO2– angle keeps increasing as we increase the number of solvent molecules (Figure a). Specifically, for Li+–CO2––3DME, only one CO2– oxygen binds to Li+ (Figure a), resulting in a weak interaction and thus increased angle and increased Li+-to-CO2– distance (by 0.9 Å, see Table and Figure c). Such a weakening interaction from bidentate to monodentate geometry is also captured by the physics model. For Li+–CO2––nH2O the increased angle and Li+-to-CO2– distance is rather minimal (Table and Figure b). These observations are important to keep in mind, as fully solvated systems indicate that the Li+–CO2– interaction is the weakest among alkali metal cations in DME but the strongest in water (Figure S4). However, we know that hard shell Li+ cation coordinates poorly to an adsorbed CO2 on the surface,[26] and thus we also compare the relative stability of two configurations of Li+–CO2––4H2O (Figure S6) as a function of the partial charge of CO2, assuming the dominant interaction comes from the electrostatic interaction. We find that for a partial charge of −1e, Li+ preferably coordinates with CO2–, but only within about 0.16 eV compared to the other configuration according to DFT. As the charge of CO2 is decreased to −0.6e, which mimics its charge state on the surface, we find that Li+ now preferably coordinates with 4H2O according to our physics-based model (Figure S6).

Figure 6

Average CO2– angle of different cation–CO2––solvent (DME or H2O) systems from AIMD simulations: (a) cation–CO2––nDME and (b) cation–CO2––nH2O (n = 0, 1, 2, 3).

Table 2

Computed Li+ Coordination Distances for Complexes with O Atoms (of DME/H2O) and C Atoms (of CO2–)

	d(Li–Osolvent)/Å	d(Li–CCO2–)/Å
Li+–CO2–		1.92
Li+–CO2––1DME	2.00	2.19
Li+–CO2––2DME	2.09	2.39
Li+–CO2––3DME	2.20	3.07
Li+–CO2––1H2O	1.86	2.63
Li+–CO2––2H2O	1.93	2.76
Li+–CO2––3H2O	2.01	2.89

Conclusions

We have used ADMP simulations to investigate the interactions of the CO2– anion radical in various chemical environments, including both aqueous and nonaqueous electrolytes. Simulations show the effect of cations and solvent molecules on the fundamental properties of bond distance and angle of CO2, e.g., DME solvation increases the interaction distance for small cations with CO2–, while for NX4+ cations the difference is minimal. This correlates well with the complexation enthalpy of CO2–, offering a guide to which combinations of solvent and supporting electrolyte may provide the desired experimental conditions. Moreover, we have shown that the identity of the supporting electrolyte cation will likely matter more in polar solvents than in solvents of less polarity. Therefore, we can tailor the strength of interactions between CO2– and supporting electrolytes through the judicious choice of electrolyte cations and solvents. Our simulation and physics-based model provides a general and suitable approach to studying solvation effects particularly in an aprotic electrolyte environment, and it was shown that the electrochemical CO2RR activity could be correlated to the bulk solvation properties in DME and DMSO.[27] We note the limitation of our model not yet including the electrode surface as well as the pH effect at the interface,[68] and future experimental and computational efforts are necessary to elucidate the interactions of the CO2 anion radical with the electrode surface.