Double sulfur vacancies by lithium tuning enhance CO2 electroreduction to n-propanol

Electrochemical CO2 reduction can produce valuable products with high energy densities but the process is plagued by poor selectivities and low yields. Propanol represents a challenging product to obtain due to the complicated C3 forming mechanism that requires both stabilization of *C2 intermediates and subsequent C1–C2 coupling. Herein, density function theory calculations revealed that double sulfur vacancies formed on hexagonal copper sulfide can feature as efficient electrocatalytic centers for stabilizing both CO* and OCCO* dimer, and further CO–OCCO coupling to form C3 species, which cannot be realized on CuS with single or no sulfur vacancies. The double sulfur vacancies were then experimentally synthesized by an electrochemical lithium tuning strategy, during which the density of sulfur vacancies was well-tuned by the charge/discharge cycle number. The double sulfur vacancy-rich CuS catalyst exhibited a Faradaic efficiency toward n-propanol of 15.4 ± 1% at −1.05 V versus reversible hydrogen electrode in H-cells, and a high partial current density of 9.9 mA cm−2 at −0.85 V in flow-cells, comparable to the best reported electrochemical CO2 reduction toward n-propanol. Our work suggests an attractive approach to create anion vacancy pairs as catalytic centers for multi-carbon-products.

where n is the number of transferred electrons, F is the Faraday constant, p = 101.3 kPa, Vgas is the volume of gas products, i is the total current detected by the electrochemical workstation, R is the gas constant, c is the molar concentration, Vliquid is the volume of anode electrolyte, and Q is the quantity of applied electric charges during the CO2 reduction. The The partial current densities (j) for products can be calculated as below, where A is the geometric area of the cathode: Computational details.
Based on density functional theory (DFT), structural optimization and energy calculations were performed within GGA-PBE 1 functional as implemented in Vienna ab-initio Simulation Package (VASP) 2 . The projector augmented wave (PAW) method was used to treat exchangecorrelation energy. A plane wave expansion with kinetic-energy cutoff of 450 eV and a residual force tolerance of 0.03 eV/Å were adopted in all calculations to reach sufficiently precise results.
The optimized cell parameters of hexagonal CuS structure are a = b = 3.81 Å and c = 16.51 Å, which is in good agreement with experimental measurement. On the basis of experimentally crystalline characterization, we built (001) and (100) facets of CuS with 3×3×1 and 3×1×1 slabs, respectively, where only the top two layers of atoms were allowed to relax. Accordingly, a 3×3×1 S5 and a 3×2×1 Monkhorst-Pack k-point meshes were used to sample the k-space, respectively. The last number is the supercell number of unit cell in direction of crystal plane. All slab models are at least consisting of four Cu-S-Cu atom layers for CuS (100) facet, while it is three Cu-S-Cu atom layers for CuS (001) facet. A 20 Å of vacuum space in z-direction was selected to remove the effect of slab interactions. To decouple the electrostatic interaction for asymmetric slabs, a dipole correction was included in all cases 3,4 . Additionally, one layer of water molecules and an extra hydrogen atom were placed onto the slab surfaces to explicitly take solvation and electric filed effects into account 4 . Plane-averaged electrostatic potential for a CuS (100) slab with one water layer on one side is shown in Supplementary Fig. 31. The long-range dispersion interactions were described by the D3 correction method 5 . The geometries and barriers of transition state for the CO* coupling process were calculated by climbing-image nudged elastic band (CI-NEB) method 6 .
By employing computing hydrogen electrode (CHE) model, the free energy correction can be written as: where E is the electronic energy, ZPE is the zero-point energy, and S is entropy. The values of these parameters were obtained from vibrational frequency calculations using harmonic approximation for adsorbates and standard thermodynamics tables for gas phase molecules. Since the values of free energy corrections are nearly the same for one adsorbate on different sulfur vacancy sites, energy and free energy change will have little difference for CO* coupling process. Thus, we used energy change to represent the energy barrier of C-C coupling. The charge density difference is defined as Δρ = ρ 100-vac + ρ S -ρ 100 , where ρ 100-vac and ρ 100 refers to the charge density of (100) facet with and without S vacancy, respectively. The iso-surface value for both single and double vacancies is 0.0018 e bohr -3 . ρ S is the charge density of the removed S atoms.

Supplementary Figures
Supplementary Fig. 1 Hexagonal CuS was constructed through the two centrosymmetric sandwich structures, exposing the high-density S atoms on the base (100) (light blue rectangle) or (001) (pink rectangle) crystal planes, and low-density S atoms on (103) Table 3 The distance of Cu-Cu extracted from the final optimized model in Fig.   1b, Fig. 1c, Supplementary Fig. 3 and Supplementary Fig. 5. The reason why SSV-CuSx model cannot dimerize the CO-COCO is not the far distance The attached illustration of basic sites constructed with seven Cu atoms (i→ vii) and two S atoms (or sulfur vacancies), as well as the definition of the measured distance extracted from the final optimized model in Fig. 1b, Fig. 1c, Supplementary Fig. 3 and Supplementary Fig. 5. During the fitting of CuS nanosheets, its coordination numbers was fixed as the same values as that of CuS theory, while the internal atomic distances R, Debye-Waller factor σ 2 , the edge-energy shift ΔE and passive electron reduction factor S0 2 were allowed to run freely. During the fitting of CuS-10 cycles, S0 2 was fixed to be the value of 0.9, and the other four parameters were variable. The error bar in the fitting results for these samples were provided.

S27
Supplementary Table 6 The FEn-PrOH values with the sulfur vacancies concentration calculated based on the double sulfur vacancies/total sulfur atoms in the DFT model in Fig. 2b, and the fitting coordination number obtained from the EXAFS in Supplementary Table 5.  /n-PrOH