Tuning the local chemical environment of ZnSe quantum dots with dithiols towards photocatalytic CO2 reduction

Sunlight-driven CO2 reduction to renewable fuels is a promising strategy towards a closed carbon cycle in a circular economy. For that purpose, colloidal quantum dots (QDs) have emerged as a versatile light absorber platform that offers many possibilities for surface modification strategies. Considerable attention has been focused on tailoring the local chemical environment of the catalytic site for CO2 reduction with chemical functionalities ranging from amino acids to amines, imidazolium, pyridines, and others. Here we show that dithiols, a class of organic compounds previously unexplored in the context of CO2 reduction, can enhance photocatalytic CO2 reduction on ZnSe QDs. A short dithiol (1,2-ethanedithiol) activates the QD surface for CO2 reduction accompanied by a suppression of the competing H2 evolution reaction. In contrast, in the presence of an immobilized Ni(cyclam) co-catalyst, a longer dithiol (1,6-hexanedithiol) accelerates CO2 reduction. 1H-NMR spectroscopy studies of the dithiol-QD surface interactions reveal a strong affinity of the dithiols for the QD surface accompanied by a solvation sphere governed by hydrophobic interactions. Control experiments with a series of dithiol analogues (monothiol, mercaptoalcohol) render the hydrophobic chemical environment unlikely as the sole contribution of the enhancement of CO2 reduction. Density functional theory (DFT) calculations provide a framework to rationalize the observed dithiol length dependent activity through the analysis of the non-covalent interactions between the dangling thiol moiety and the CO2 reduction intermediates at the catalytic site. This work therefore introduces dithiol capping ligands as a straightforward means to enhance CO2 reduction catalysis on both bare and co-catalyst modified QDs by engineering the particle's chemical environment.


Photocatalytic CO2 reduction
Sample preparation. A ZnSe-BF4 stock solution (64.1 µM in DMF, 23.40 µL), a capping ligand solution (5.0 mM in DMF, typically 30 µL) followed by co-catalyst solution (Ni(cycP, 5 mM in H2O) were added stepwise to a Pyrex glass photoreactor (Chromacol 10-SV, Fisher Scientific) containing a magnetic stirrer bar. The mixture was diluted with ascorbic acid (AA, 0.1 M in water, pH adjusted to 7 with NaOH) to a total solution volume of 3 mL. In the absence of Ni(cycP), NaHCO3 powder (25 mg) was further added to increase the pH to 8.3. The photoreactor was then sealed with a rubber septum and pierced with two needles (inlet and outlet).
Constant flow-setup with automated product quantification. The inlet of the photoreactor was connected to a Mass Flow Controller (Brooks GF040) supplying a stream of CO2 (CP Grade, BOC, humidified with a water bubbler) with a flow rate of 4.0 sccm. The flow rate at the GC outlet was verified prior to the experiment with an Alicat gas flow meter to avoid gas leakage. The outlet of the photoreactor was connected to a flow selection valve controlled by a Shimadzu Tracera GC-2010 Plus gas chromatograph for product quantification of the gaseous reaction products (see below). Six samples (two triplicates of identical conditions) were typically analyzed in parallel. Upon purging with a constant stream of CO2, the solution pH decreased to 6.5 (absence of Ni(cycP)) or 5.5 (presence of Ni(cycP) due to saturation with CO2. The photoreactor was purged for a further 45 min in the dark and sampled via online GC quantification. The first two injections of each sample were used to determine a "background" peak which was subtracted from further injections. The photoreactor was then placed in a water bath maintained at 25 °C, stirred and irradiated by a solar light simulator (Newport Oriel, AM 1.5G, 100 mW cm -2 ). The six samples were evenly distributed within the light simulator to account for possible variations of the light intensity depending on the position in the simulator. UV irradiation was filtered with a 400 nm cut-off filter (UQG).
For 13 C isotopic labelling, photocatalysis experiments were performed as described above, but with accumulating products in the headspace under steady-state conditions and using 13 CO2 as the headspace gas. After 1000 min (16.7 h), the photoreactor headspace was transferred to an evacuated gas IR cell (SpecAc, 10 cm path length, equipped with KBr windows) and a high-resolution gas-phase transmission spectrum was collected.
The Shimadzu Tracera GC-2010 Plus gas chromatograph (GC) used a barrier discharge ionization detector, kept at 300 °C, and was equipped with a Hayesep D (2 m * 1/8'' OD * 2 mm ID, 80/100 mesh, Analytical Columns) pre-column and a RT-Molsieve 5A (30 m * 0.53 mm ID, Restek) main column in order to separate H2, O2, N2, CH4 and CO while hindering CO2 and H2O to reach the Molsieve column. He carrier gas (Grade 5.0, BOC) was purified (HP2-220, VICI) prior to entering the GC. The column temperature was kept at 85°C. The gaseous flow from the flow selection valve was passed through a loop (volume 1.0 mL) and injected approximately every 4.25 min into the GC. Effectively, each individual sample was injected every 25.5 min. The GC calibration was performed with a known standard for H2, CO and CH4 (2040 ppm H2/2050 ppm CO/2050 ppm CH4 in balance gas CO2, BOC) by diluting the mixture with pure CO2. Data analytics. The data was processed and visualized using the statistical programming language R with the tidyverse library. 3,4 First, the flow rates were corrected by subtracting a "background" peak obtained in the dark prior to irradiation (we noticed a marginal CO background peak depending on the residual amount of oxygen present in the sample stream -a feature of the gas chromatograph and not the sample). Second, the momentary product evolution rate corresponding to each injection was calculated using the following formula. ‫﮲‬ = * ̇ * ! where p is the pressure in the photoreactor (ambient pressure, 101325 Pa), is the flow rate (4.0 sccm), R is the universal gas constant, T is the temperature (298 K) and fi is the response factor for each gas determined by the calibration procedure. Third, the total amount of evolved product was calculated using trapezoidal integration of the product evolution rates. The three independent replicates of identical conditions were averaged by calculating the mean and standard deviation over irradiation time and sample. For visual display, the actual values for each sample are plotted as transparent scatter, whereas the mean is represented as a smoothened continuous line. In addition, the standard deviation is visualized by the shaded area surrounding the mean where the transparency is proportional to the standard deviation. Specifically, the calculated standard deviation is used to compute a Gaussian density for that standard deviation, plotting a cloud with the opacity proportional to the density. This appears as a vertical "cloud" of uncertainty 5 . The maximum of the uncertainty cloud is set to 1 standard deviations.

Computational methods
Periodic density functional theory (DFT) calculations reported in this work were carried out using the Perdew-Burke-Ernzerhof (PBE) functional 6 with dispersion corrections through the zero-damping DFT-D3 method by Grimme,7 as implemented in the Vienna Ab Initio Simulation Package (VASP) software, version 5.4.4. 8 The core electrons of the Zn, Ni, Se, S, N, O, C, P and H atoms were described by projector-augmented wave (PAW) pseudopotentials, 9 while their valence electrons were expanded in plane waves with a kinetic energy cut-off of 500 eV. Geometry optimizations were performed with an energy convergence criterion for the electronic steps of 10 -6 eV and a force criterion of 0.015 eV/Å 2 for the ionic steps, using a Gaussian smearing with a width of 0.05 eV.
The ZnSe cubic bulk structure (space group F4 9 3m), displayed in Figure S11A, was retrieved from the Materials Project database (mp-1190). 10 The equilibrium lattice constant was determined by optimizing the structure with lattice parameters ranging ±5% of the initial value and fitting the energies to the Birch-Murnaghan equation of state, while sampling the reciprocal space with Γ-centered k-point grids of 3×3×3, 5×5×5, 7×7×7, 9×9×9 and 11×11×11. Following a convergence criterion of 1 meV atom −1 , the 5×5×5 k-point grid corresponding to a density of 28.7 points Å −1 was selected for surface calculations; this bulk displayed a band gap of 1.17 eV ( Figure S11B), which is in good agreement with the values predicted theoretically at similar levels of theory. 11,12 To describe the QD surface, we employed the 4-layered (200) slab model reported in our previous work 1 with a vacuum space of 15 Å, wherein the two bottom layers were fixed to their bulk positions and the two topmost layers were allowed to relax. Gas molecules were optimized at the Γ-point and with 15 Å of vacuum along the three axes.
Gibbs energy corrections to the electronic energy were computed with the ASE thermochemistry module 12 at the temperature of 298 K and pressure of 1 atm. For gas-phase molecules, these corrections were calculated via the ideal gas limit method, which include translational, rotational, and vibrational contributions. For adsorbates, only vibrational contributions treated harmonically were considered. For the adsorbed Ni(cycP) cocatalyst and dithiol ligands, Gibbs corrections were determined in the systems ZnSe | Ni(cycP) and ZnSe | DT, respectively, and then applied to the combined system ZnSe | Ni(cycP) | DT. This approach was validated by computing the Gibbs corrections for EDT in ZnSe | Ni(cycP) | EDT, where EDT was monodentately adsorbed on site H, or bidentately on sites H-B (see Figure 6C), finding energy differences of only 0.01 eV compared to the energies obtained with the system ZnSe | EDT. Similarly, the Gibbs corrections for the CO2RR intermediates on the Ni(cycP) cocatalyst in ZnSe | Ni(cycP) | DT were adopted from the system without the adsorbed dithiol. The validity of this approach was confirmed by computing the corrections for the *CO2 intermediate in the system Ni(cycP) | ZnSe and comparing them with those calculated in ZnSe | Ni(cycP) | EDT, obtaining energy differences which are deemed negligible (< 1 meV).
Non-covalent interactions (NCIs) were assessed by taking the geometry of the adsorbed species optimized in VASP and computing their electron density via single-point calculations with the dispersion-corrected hybrid exchange-correlation functional ωb97xd 13 and the Gaussian09 software, Revision E.01. 14 In these calculations, Ni atoms were described with the effective core potential Lanl2dz and an additional f-polarization function (exponent = 3.130), 15 while the 6-31g(d,p) basis set was used to describe the C and H atoms. The latter basis set with an additional p-diffuse function was employed to represent the electrons of the more electronegative S, N, P and O atoms. Based on the electron density ρ, the reduced density gradient s(ρ) was calculated according to Eq. 1 using the software Critic2. 16 We note that at regions far from the molecule, both ρ and s(ρ) approach zero exponentially. Similarly, s(ρ) approaches zero at regions where NCIs are present, and becomes identically zero at the covalent bonds, i.e. critical points of ρ. To distinguish between attractive and repulsive interactions, the Laplacian of the density, ∇ 2 ρ, is computed. For repulsive interactions, ρ is expected to be a minimum, and therefore, all components of ∇ 2 ρ along the 3 dimensions are positive. On the other hand, for attractive interactions, ρ is characterized by a maximum along one direction whose eigenvalue in the Hessian matrix of ρ is called λ2. NCI-plots are constructed by plotting s(ρ) against ρ multiplied by the sign of λ2, sign(λ2)ρ. These plots allow to distinguish between bonded (λ2 < 0) and nonbonded (λ2 > 0) NCIs, while the value of ρ indicates the strength of these interactions. 17

Assessment of the ZnSe-QD surface coverage
To model the adsorption of the dithiol ligands and the Ni(cycP) cocatalyst, the ZnSe-QDs coverage was calculated, based on the experimental concentration of species on the QD surface according to the equation: where rQD (2.27 nm) is the experimental radius of the QDs, and a (5.75 Å) and b (4.06 Å) are the lattice constants of the modelled surface unit cell.
For EDT, the experimental coverage for the system ZnSe | EDT corresponded to a single molecule adsorbed in a p(1×1) unit cell, while for HexDT and OctDT this corresponded to a single molecule in a p(1×3) and p(3×1) unit cell, respectively. In the systems ZnSe | Ni(cycP) | DT, the modelled coverage consisted in one cocatalyst and two dithiol molecules adsorbed in a p(3×2) cell.

Modelling of the adsorption modes in ZnSe | Ni(cycP) | DT
In ZnSe | Ni(cycP) | EDT, the monodentate adsorption of EDT was modelled on site H ( Figure 6C) as this is the only position from which this capping ligand can interact with the CO2RR intermediates adsorbed on the Ni(cycP) cocatalyst. To model the bidentate mode of EDT, the sites H-B were considered as the ligand in this position has a minimum distance of 3.31 Å to the Ni(cycP) cocatalyst and its binding energy was expected to be very similar to that on the other available Zn sites further away from the cocatalyst (Table S5). This was verified by adsorbing EDT on the sites D-E, which resulted in a binding energy for EDT that differs only by 0.01 eV from that on H-B.
For HexDT and OctDT, due to the flexibility of the ligand chains, it was possible to identify more than one bidentate adsorption besides all the monodentate ones. Consequently, for all the positions from A to G ( Figure  6C), the dithiol was directed once along the a and b axes, and diagonally, excluding orientations that would bring the dithiols on the Zn sites immediately around the Ni(cycP), as they were deemed unlikely due to steric hindrance.

Neglection of the bidentate ligand in ZnSe | Ni(cycP) | DT
In ZnSe | Ni(cycP) | DT, of the two dithiols in the p(3×2) cell, one was initially modelled as bidentate since this was found to be the most favorable adsorption in the system ZnSe | DT, while the other one was considered to be monodentate to interact with the CO2RR intermediates adsorbed on the cocatalyst. In subsequent reactivity studies, however, the bidentate ligand was removed from the calculation for two main reasons. Firstly, this ligand does not interact with the second coordination sphere of the cocatalyst, and therefore, it may not influence the energetics of the CO2RR intermediates; and secondly, the distance between the neighboring Zn sites (5.75 Å and 4.06 Å along the a and b axes, respectively), is long enough to minimize the interaction between the adsorbed dithiols. Accordingly, the binding energy of two dithiols in the p(3×2) cell is almost equal to the sum of the individual binding energies, as shown in Table S6 for a sample of control combinations computed for HexDT. In short, an HexDT ligand in the most stable bidentate mode was modelled along with another HexDT molecule in the available monodentate configurations. Based on these results, it is reasonable to assume that this approach will also be valid for EDT and OctDT as these ligands differ from HexDT only in the chain length, but the orientations of the adsorption are the same.