Solvent-mediated outer-sphere CO2 electro-reduction mechanism over the Ag111 surface

The electrocatalytic CO2 reduction reaction (CO2RR) is one of the key technologies of the clean energy economy. Molecular-level understanding of the CO2RR process is instrumental for the better design of electrodes operable at low overpotentials with high current density. The catalytic mechanism underlying the turnover and selectivity of the CO2RR is modulated by the nature of the electrocatalyst, as well as the electrolyte liquid, and its ionic components that form the electrical double layer (EDL). Herein we demonstrate the critical non-innocent role of the EDL for the activation and conversion of CO2 at a high cathodic bias for electrocatalytic conversion over a silver surface as a representative low-cost model cathode. By using a multiscale modeling approach we demonstrate that under such conditions a dense EDL is formed, which hinders the diffusion of CO2 towards the Ag111 electrocatalyst surface. By combining DFT calculations and ab initio molecular dynamics simulations we identify favorable pathways for CO2 reduction directly over the EDL without the need for adsorption to the catalyst surface. The dense EDL promotes homogeneous phase reduction of CO2via electron transfer from the surface to the electrolyte. Such an outer-sphere mechanism favors the formation of formate as the CO2RR product. The formate can undergo dehydration to CO via a transition state stabilized by solvated alkali cations in the EDL.


S1. Computational methods:
CMD model: CMD simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) (version from 22Aug 2018). 1 An aqueous solution of KCl that is a standard electrolyte in electrocatalytic studies, was confined between two charged Ag walls that represented the cathodic and anodic surfaces. The silver slabs with an electrolyte bulk, that consisted of K + , Cland CO 2 and H 2 O molecules were constructed using Atomic Simulation Environment (ASE). The concentrations were set to correspond to the experimentally relevant conditions, namely, 0.86M for K + , 0.86M for Cland 0.06M for CO 2 . The simulations were carried out using a supercell of 33.1x37.2x265.5 Å. The supercell was set to be periodic in the x and y directions and non-periodic in the z-direction. A vacuum spacing of 1.7 Å was set behind the silver slabs. The charging of the electrode surfaces due to externally applied potential was mimicked by introducing negative and positive ghost charges behind the Ag slabs. The calculated magnitudes of the electric field on the cathodic surface due to polarization induced from the ghost charges were 0 V/nm (no ghost charges), 0.05 V/nm, and 0.5 V/nm. The interatomic interactions between the components of the medium were modelled via Lennard-Jones (LJ) potential with the cut-off of 9.0 Å. The cross-terms were obtained using the Lorentz-Berthelot mixing rules. Water was modelled using the parameters of the SPC/E model, 2 that is widely used in simulations of the ions in aqueous solutions. 2b,c,d The LJ parameters for CO 2 , K + , and Clwere taken from the literature. 3,4 The long-range interactions were calculated with the particle-particle-particle-mesh (PPPM) method. The O-H bonds and H-O-H angles were kept rigid within the SHAKE algorithm, whereas the C-O bonds and O-C-O angles were kept rigid via harmonic bond approximation. The temperature of the simulations was maintained using the Nose-Hoover thermostat at the value of 300 K. 3,4 We used a time step of 1 fs. The production runs were 90 ns long, excluding 10 ps of equilibration time in the beginning. The snapshots of trajectory for the further analysis were taken every 10 ps. The post-processing analysis was conducted using the MDAnalysis 20 and maicos_delft (https://gitlab.com/mdopke/maicos_delft) packages. For all the simulations, the system was initially equilibrated in the isothermal regime with different integration time steps (10 -5 , 10 -4 , 10 -3 , 10 -2 , 10 -1 and 1) for 10 3 steps within each time step. However, the snapshots of atom quantities in the final dcd trajectory were included every 10 ps timesteps, giving the first 10 ps as equilibration time. Figure-SI 1 (left) shows that the energy of the production runs starting from the 10 ps is stabilized, with some small fluctuations present. To calculate the self-diffusion coefficient of water, we have calculated the mean-squared displacement (MSD). The slope of the linear regression fit of MSD versus simulation time is proportional to the diffusion coefficient. For calculating the ensemble average, we have considered all the water molecules and multiple time origins. To indicate the difference in self-diffusivity in the bulk and next to the anode/cathode, the model was divided into 8 bins.
In Figure-SI 1 (right), we show the MSD versus time plot of water molecules along the z dimension -the longest in the model. The time frame between 20 -400 ps was used for linear regression of the MSD z versus time. As the first and the last bins were located near the electrodes (1 and 8), the self-diffusion coefficient were calculated for the central bins (2 -7) which represent bulk diffusivity of water. The calculated diffusion coefficient based on the linear regressions resulted in a mean value of 2.78 •10 -9 m 2 /s with the standard error of the mean equal to 0.01 •10 -9 m 2 /s which compares well with the experimental value of 2.3•10 -9 m 2 /s. 21a , and CMD simulations of pure water using the SPC/E model at 298 K (in the range of 2.6 -3.1•10 -9 m 2 /s depedneing on the integrator time-step). 21b,c  The snapshot of the KCl electrolyte with CO 2 confined between two silver slabs simulated at zero potential. The surface on the left represents the anode, with positive ghost charges imposed behind the wall, while the right surface represents the cathode, with negative ghost charges imposed behind the wall. The color code is following: silver is grey, oxygen is red, hydrogen is white, potassium is violet, chlorine is green, carbon of CO 2 is yellow, the ghost atoms are light brown.
The electric field above the silver cathode was imposed by adding the ghost atoms behind the cathode and anode and was calculated according to Gauss's law that relates the distribution of the electric charge to the resulting electric field. The electric flux can be expressed as: where Θ is the electric charge enclosed and is the dielectric constant of water and is the dielectric constant of the vacuum. The values of the total induced charge were calculated according to the following: where is the charge assigned to each ghost atom in LAMMPS simulation (for units style real), is the electron charge equal to 1.6•10 -19 C and 25 is the number of ghost atoms imposed behind each surface. We have assumed that Ag is a perfect conductor and there is no charge density trapped in the bulk volume of Ag. We have also assumed that the dielectric constant of water does not change in the EDL.
The magnitude of the electric field passing through the cathode: where is the surface area of the cathode, which is assumed to be a rectangular prism.

Supporting Information
DFT calculations: we constructed a smaller molecular model of the electrocatalyst system consisting of a 4x4x5 Ag111 slab and performed density functional theory (DFT) calculations for the reactive events. The EDL was composed of 4 Na + cations per supercell which corresponds to ~ -1.0 V vs pzc (vide infra). The supercell was charge neutral. We performed geometry optimizations using periodic DFT calculations using the PBE 5 functional and projector augmented wave (PAW) potentials 19 with the the valence electrons expanded as plane waves (cutoff 450 eV) via VASP 5.4.4 suit of software. 6 In all DFT based geometry optimization was chose a fcc supercell of dimensions 11.53715 Å, 11.53715 Å, and 30.0 Å. The unit cell was constructed in the following manner: first 5 layers of 1x1 fcc cell of Ag111 was created with a lattice constant of 4.079 Å 18 and lattice vectors: (0.70710678, 0.0000000, 0.000000), (-0.35355339, 0.6123724, 0.000000), and (0.000000, 0.000000, 5.1961524). This 1x1 supercell was extended to create the 4x4x5 slab by changing the Bravais Matrix and adding additional atoms. The 4x4x5 Ag111 slab was fully optimized using DFT, and so were "4x4x5Ag111 slab+water molecules", "4x4x5 Ag111 slab+water molecules + Na ions", "4x4x5 Ag111 slab+water molecules + Na ions + CO2", 4x4x5 Ag111 slab+water molecules + CO2". Grimme's D3 method with Becke-Johnson damping (D3(BJ)) was adopted for the van der Waals correction. 13,14 In the DFT based geometry optimizations we considered 24 water molecules for the solvation of the 4Na + cations over the 4x4x5 Ag111 slab. We further placed a CO 2 molecule at the EDLvacuum interface and performed geometry optimizations. Calculations without Na + cations, that is Ag111-(H 2 O) 24 -CO 2 system were also performed where CO 2 was found to adopt a linear configuration upon geometry optimization. For the Ag111-(Na + ) 4 -(H 2 O) 24 -CO 2 system, DFT based geometry optimizations resulted in CO 2 adopting a bent configuration and involved in hydrogen bonding with protons from H 2 O moieties.

AIMD simulations:
We further added 37 additional water molecules to solvate the entire supercell optimized via DFT. This system comprising of Ag111-4Na + -(H 2 O) 61 -CO 2 was used for AIMD simulations in an NVT ensemble at 360 K using the Nose-Hoover thermostat. The BLYP exchange correlation functional was used for AIMD simulations. This combination of the XC functional and simulation temperature provides better description of the structure and dynamics of water, and the lower computational cost. [7][8][9][10][11][12] We reran a part of the constrained AIMD simulations using the PBE functional which resulted in similar barrier as obtained with the BLYP functional (vide infra). Grimme's D3 method with Becke-Johnson damping (D3(BJ)) was adopted for the van der Waals correction. 13,14 The INCAR files used in our DFT and AIMD simulations are provided with the SI and can be referred to for further details on the simulation parameters. To generate a first-guess of the Gibbs free energy profile slow-growth approach (SGA) simulations were performed to traverse the Gibbs free energy surface from the reactant to the product state via the transition state using a reaction coordinate (Q) defined via a collective variable (bond lengths; bond angle etc.). [15][16][17] To further refine the Gibbs free energy profile, several intermediate values of Q were chosen and subjected to long (> 10 ps) constrained (fixed value of Q) AIMD simulations. In all AIMD simulations, constrained and unconstrained, we used an integrator time step of 1 fs. The average force (F) required to maintain the constraint was computed on an equilibrated trajectory via the bluemoon sampling method. The convergence of the force required to maintain the constraint was visually checked by plotting the force profiles for the last 5 ps of the simulation trajectory for each constrained AIMD run. The computed force profiles have been provided with the computational data provided via the 4TU database. The average force calculated on the last 1 ps of the equilibrated trajectory was integrated from the reactant (Q i ) to product (Q f ) state to obtain the corresponding Gibbs free energy change (Δ ): Therefore, C H is estimated to be 18 μF cm −2 assuming that the relative permittivity of water in the EDL is at the minimum value which is a limiting case. We therefore also consider C H = 25 μF cm −2 and 50 μF cm −2 in our estimation of the surface potential. q which corresponds to the total surface charge on the reactive face can be taken as -1.3 (total Bader net atomic charge on the top layer) or -1.46 which is the total Bader net atomic charge on the top 2 layers. The estimated surface potentials are tabulated in Table-SI 2.

S5. Formate formation during homogeneous CO 2 reduction
An outer-sphere ET would reduce CO 2 in the homogeneous phase. To gain insight into this homogeneous CO 2 reduction we carried out AIMD simulations of CO 2 in water (23 H 2 O and 1 CO 2 per supercell). We performed three simulations: neutral, anionic, and dianionic. The neutral simulation consisted of a solvated CO 2 molecule in a neutral supercell. Introduction of 1ein the neural simulation resulted in the anionic system where CO 2 molecule was found to bend (<∠O=C=O> = 134.3°±3.1°). Both O moieties formed hydrogen bonds with water, and the C moiety was found to have a strong H-bonding interaction but did not get protonated during the simulations (Figure-SI 7a). Introduction of a second electron led to a doubly negatively charged dianionic system where formic acid quickly (within 50 fs of simulation time) formed ( Figure-SI 7b,c). From the AIMD simulations of solvated CO 2 we conclude that an outer-sphere 2etransfer to CO 2 would drive the formation of formic acid/formate.   S10. Data availability.

S8. PBE versus BLYP
Complete datasets including sample input files, optimized geometries, and trajectories from DFT, AIMD and CMD simulations including computed forced from constrained AIMD runs will be deposited to 4TU database and will be publically available under the DOI: 10.4121/19142303.