Reduction of CO2 by a masked two-coordinate cobalt(i) complex and characterization of a proposed oxodicobalt(ii) intermediate

Fixation and chemical reduction of CO2 are important for utilization of this abundant resource, and understanding the detailed mechanism of C–O cleavage is needed for rational development of CO2 reduction methods. Here, we describe a detailed analysis of the mechanism of the reaction of a masked two-coordinate cobalt(i) complex, LtBuCo (where LtBu = 2,2,6,6-tetramethyl-3,5-bis[(2,6-diisopropylphenyl)imino]hept-4-yl), with CO2, which yields two products of C–O cleavage, the cobalt(i) monocarbonyl complex LtBuCo(CO) and the dicobalt(ii) carbonate complex (LtBuCo)2(μ-CO3). Kinetic studies and computations show that the κN,η6-arene isomer of LtBuCo rearranges to the κ2N,N′ binding mode prior to binding of CO2, which contrasts with the mechanism of binding of other substrates to LtBuCo. Density functional theory (DFT) studies show that the only low-energy pathways for cleavage of CO2 proceed through bimetallic mechanisms, and DFT and highly correlated domain-based local pair natural orbital coupled cluster (DLPNO-CCSD(T)) calculations reveal the cooperative effects of the two metal centers during facile C–O bond rupture. A plausible intermediate in the reaction of CO2 with LtBuCo is the oxodicobalt(ii) complex LtBuCoOCoLtBu, which has been independently synthesized through the reaction of LtBuCo with N2O. The rapid reaction of LtBuCoOCoLtBu with CO2 to form the carbonate product indicates that the oxo species is kinetically competent to be an intermediate during CO2 cleavage by LtBuCo. LtBuCoOCoLtBu is a novel example of a thoroughly characterized molecular cobalt–oxo complex where the cobalt ions are clearly in the +2 oxidation state. Its nucleophilic reactivity is a consequence of high charge localization on the μ-oxo ligand between two antiferromagnetically coupled high-spin cobalt(ii) centers, as characterized by DFT and multireference complete active space self-consistent field (CASSCF) calculations.

developed by Noodleman 8,9 was employed to describe states involving antiferromagnetic coupling. The following thresholds were used for optimizations: energy change tolerance of 5 × 10 −6 a.u., root mean square gradient of 1 × 10 −4 a.u., maximum gradient of 3 × 10 −4 a.u., root mean square displacement of 2 × 10 −3 a.u., maximum displacement of 4 × 10 −3 a.u. The DFT integration grid was always set to Grid6 while the Hartree Fock exchange was evaluated by GridX6 for calculations with the B3LYP functional. Tight convergence criteria (energy tolerance 1 × 10 −8 a.u.) for self-consistent field (SCF) calculations were employed throughout all calculations. The atom-pairwise dispersion correction with Becke-Johnson damping (D3BJ), as implemented in ORCA 10,11 was used throughout to treat non-covalent interactions. The resolution of identity (RI) approximation [12][13][14] was used to accelerate the calculations in conjunction with the auxiliary Coulomb-fitting basis set def2/J to accelerate S-6 the calculations at the BP86 level. To account for relativistic effects, scalar relativistic zeroth order approximation (ZORA) was used. [15][16][17][18] Geometry optimizations were performed without any constraints. Solvent effects during optimization and single point calculations were modeled by employing the polarizable continuum model (CPCM) 19 using benzene (ɛ = 2.3) as the solvent. Harmonic vibrational analyses were undertaken after the geometry optimization of each structure to differentiate local minima from transition states. The local minima were characterized by all positive eigenvalues of the Hessian matrix, whereas the transition states (TSs) have only one negative eigenvalue. The zero-point energies (ZPE), enthalpic corrections and entropic contributions were obtained from the frequency calculations. In order to obtain single-point energies closer to the basis set limit, B3LYP 4,20 and M06L 21 calculations with the larger def2-TZVPP 22 basis set on all elements were carried out for which the auxiliary basis set def2-TZVPP/J was employed for the RIJCOSX or RI approximation as appropriate. 23 Empirical dispersion correction 24,25 with the D3BJ keyword was employed all throughout except for calculations with M06L. The BP86 calculations on TS4 and 11 (at S = 0) leads to an erroneous change in the local spin state of the two cobalt centers, deviating from the desired electronic structure. Since both species are involved in the rate-determining step of the dinuclear dissociative pathway, we re-optimized them (at S = 0 and 3 spin states) at the B3LYP/B1 level of theory. To be consistent, we optimized structures corresponding to rate-determining transition states and the reference intermediates/adducts among various pathways (such as 5, TS2, 6, 10, TS4, 11, 13, 14, and TS7) at B3LYP/B1 and evaluated single point energies at M06L/def2-TZVPP. Reported relative free energies in this work are calculated at the M06L/def2-TZVPP level of theory unless otherwise mentioned.
Finally, highly accurate electronic energies were computed by using open-shell DLPNO-CCSD(T) 26,27 approach in conjunction with the def2-TZVPP basis set and an additional correlation fitting auxiliary basis set, def2-TZVPP/C. For pair natural orbital (PNO) S-7 generation, the default values of TCutPNO, TCutPairs and TCutMKN were used. 28 Localized orbitals according to the Pipek-Mezey scheme were computed before PNO generation. Solvent corrections to coupled-cluster electronic energies were obtained from the M06L calculations. Due to the challenging computational cost, only crucial barriers were computed with coupled-cluster theory and the rest with M06L/def2-TZVPP. The gas phase Gibbs free energy change on addition/binding of CO2 was corrected to the standard state (1 M) by using a correction factor of 1.89 kcal/mol. 29 This value is added to the relative free energy, when one molecule of CO2 interacts with the mononuclear or dinuclear metal complexes to form the corresponding CO2 adducts, which requires transfer of CO2 from the gas phase into benzene solution.
To determine the electronic structure of 3, we employed the multi-reference CASSCF 30,31 approach. 32 The crystal structure of 3 was used as a starting point for geometry optimization, and to reduce computational cost, the two tert-butyl groups in the ligand were replaced with methyl groups. Geometry optimizations were performed using the BP86 functional and the B1 basis set combination. This approximation does not alter the electronic structure of the Co-O-Co moiety. The state-specific CASSCF (in the singlet spin state) calculations employed the triple-z quality def2-TZVPP basis set for all atoms with optimization of the first root. The active space consists of all 3d orbitals of both cobalt centers and the 2p orbitals of the bridging oxo ligand, which results in the active space CASSCF (20,13), 20 electrons distributed into 13 orbitals. Furthermore, to explore the effects of electronic excitations from filled π to vacant π* orbitals of the ligand, larger active spaces were constructed including the bonding and antibonding π-orbitals of the N-C-C-C-N of the β-diketiminate framework. This active space constituted of 24 electrons in 17 orbitals (CASSCF (24,17)). To make the computation of CASSCF (24,17) faster, the configuration interaction here has been treated with the ACCCI approximation. . 1 H NMR spectra of reaction mixture of L tBu Co (1) with 10 equiv CO2 monitored over time at 10 °C in C6D12. Integration with respect to a nickelocene internal standard shows that L tBu Co(CO) (2) is made in 95% spectroscopic and L tBu Co(OCO2)CoL tBu (4) is made in 86% spectroscopic yield after 3 h. Figure S4. 1 H NMR spectrum (25 °C) of crude product mixture from reaction of L tBu Co with 10 equiv of CO2 in C6D12. Figure S5. 1 H NMR spectra of a reaction mixture of 3 before (below) and 2 min after (above) adding 1.5 equiv of CO2 at 10 °C in C6D12. Integration relative to a nickelocene internal standard indicates formation of 4 in 75% spectroscopic yield.

C. NMR kinetics of CO2 reduction by L tBu Co.
A general procedure for one flooding concentration of CO2 is described. Under an Ar atmosphere, a stock solution of 1 in C6D12 (0.200 mL, 13.4 mM, 0.00268 mmol), 0.400 mL of C6D12, and a nickelocene capillary integration standard were added to a J. Young NMR tube. The NMR tube was sealed with a rubber septum and the solution was flash frozen in a liquid N2 bath. A stock solution of CO2 in C6D12 (0.400 mL, 72.9 mM, 0.0292 mmol) was injected anaerobically by syringe. The solution was thawed and injected into an NMR spectrometer pre-cooled to 10 °C. Proton NMR spectra were acquired at 10°C every 10 min for 5 h. This general procedure was followed to measure two different flooding concentrations of CO2 as well as two initial concentrations of 1 as indicated below. The concentration over time profiles of L tBu Co (1), L tBu Co(CO) (2), and L tBu Co(µ-OCO2)CoL tBu (4), as determined by integration of proton NMR data, were fit to first-order kinetics for each run. The value of the rate constant k for each run was then calculated as the average of the three observed rate constants obtained for a given experiment. Using these rate constants in the Eyring equation, a value of DG ‡ =21.0 kcal/mol is obtained.

Kinetics Conditions
• . Reactant and product concentration over time data as determined by NMR spectroscopy in C6D12 at 10 °C; data and first-order fits correspond to Run 1 above. Figure S7. Reactant and product concentration over time as determined by NMR spectroscopy in C6D12 at 10 °C; data and first-order fits correspond to Run 2 above.
S-14 Figure S8. Reactant and product concentration over time as determined by NMR spectroscopy in C6D12 at 10 °C; data and first-order fits correspond to Run 3 above.

D. Attempted synthesis of CO coordinated oxo complex 11
General procedure for CO gas addition to complex 3: Solid 3 (2.7 mg, 0.0023 mmol) was dissolved in 0.6 mL benzene-d6 and transferred to a J. Young NMR tube containing a nickelocene capillary used as an integration standard. On the Schlenk line, the NMR sample solution was frozen using a liquid nitrogen cold bath. The headspace above the frozen solution was evacuated. CO (7 mbar, 12.55 mL, 0.0035 mmol) was condensed on the frozen solution for 20 min and then the NMR sample was sealed and thawed to ambient temperature with vigorous shaking. Proton NMR indicated formation of the known cobalt bis(carbonyl) complex L tBu Co(CO)2 3 in 34% spectroscopic yield in addition to unidentified minor species. Figure S9. Proton NMR spectra of reactions between 3 and CO. S-16

E. X-ray absorbance spectroscopy of L tBu Co(O)CoL tBu (3) and reference compounds
High Energy Resolution Fluorescence Detection X-ray Absorption Spectroscopy (HERFD-XAS) data were collected at the Cornell High Energy Synchrotron Source at C-line end station. The incident Xrays were monochromated using a Si(220) crystal monochromator. Emitted X-rays were analyzed using a set of three spherically bent Si crystals (531 reflection) mounted on a Rowland circle in combination with a Pilatus Area Detector as previously described. 33,34 He-filled plastic bags were used to displace air from the beam flight path to minimize attenuation of the fluorescence. The Co HERFD-XAS spectra were collected by monitoring the intensity of the Ka1 peak as a function of incident energy. Data were collected from 7650 to 7950 eV. A step size and count time of 10 eV (1 sec), 0.3 eV (5 sec), and 2 eV (2 sec) was used from 7650 to 7690 eV, 7690 to 7750 eV and 7750 to 7950 eV respectively. Data were collected at ca. 20 K in a Displex cryostat to minimize photoreduction and to maintain the samples in an inert atmosphere. The data were normalized using PyMCA and six experimental spectra were averaged for each compound. For normalization, the postedge region (E > 7750 eV) was set to an absorbance of 1.0. The first inflection point of Co foil was used as a calibrant, adjusting the energy axis to the reported literature value of 7709.5 eV. All samples were prepared in a nitrogen filled glovebox as 5% Co (wt/wt) solid dilutions in BN. The method of preparation for the last two compounds was recently described. 35 S-17

F. X-ray crystallography
Crystals were placed onto the tip of a 0.1 mm diameter glass capillary tube or fiber and mounted on a Bruker SMART APEX II CCD Platform diffractometer for a data collection at 100.0(5) K. 1 The full data collection was carried out using MoKα radiation (graphite monochromator). A randomly oriented region of reciprocal space was surveyed: six major sections of frames were collected with 0.50° steps in ω at six different ϕ settings and a detector position of -38° in 2θ. The intensity data were corrected for absorption. 2 Final cell constants were calculated from the xyz centroids of >3800 strong reflections from the actual data collection after integration. 3 Structure (3) was solved using SIR97 4 and refined using SHELXL-2012. 5 The space group P-1 was determined based on intensity statistics. A direct-methods solution was calculated which provided most non-hydrogen atoms from the E-map. Full-matrix least squares / difference Fourier cycles were performed which located the remaining non-hydrogen atoms. All non-hydrogen atoms were refined with anisotropic displacement parameters. All hydrogen atoms were placed in ideal positions and refined as riding atoms with relative isotropic displacement parameters. The structure suffered both from non-merohedral twinning and highly disordered solvent. After the non-merohedral twin law, [ -0.510 0.497 -0.020 / 1.470 0.490 -0.061 / -0.490 -0.497 -0.980 ], a 180 degree rotation about reciprocal lattice [ -1 -3 1 ], was determined, 6 the data were re-integrated, 3 and a new absorption correction was applied. There were 21623 unique reflections solely in the first component, 21604 unique reflections solely in the second component, and 4867 unique overlapping reflections. The mass ratio of the two components refined to 60:40. However, since the structure also contained highly disordered solvent (toluene and perhaps hexane), the overlapping reflection contributions from the minor mass component of the twin were removed. 7 Once the data were modified for twinning, they were further modified to remove reflection contributions from the disordered solvent. 8 There were 4 sites of ~200 Å 3 each in which around 75 electrons were removed, for a total of 297 electrons in 829 Å 3 removed per unit cell. Since the exact amount and identity of solvent were unknown, no solvent was included in the atom list or molecular formula. Structure (4) was solved using SIR97 4 and refined using SHELXL-97. 9 The space group Cc was determined based on systematic absences and intensity statistics. A direct-methods solution was calculated which provided most non-hydrogen atoms from the E-map. Full-matrix least squares / difference Fourier cycles were performed which located the remaining non-hydrogen atoms. All nonhydrogen atoms were refined with anisotropic displacement parameters. All hydrogen atoms were placed in ideal positions and refined as riding atoms with relative isotropic displacement parameters. The structure was refined as an inversion twin, with a final twin occupancy ratio of 87:13.  Figure S11. Overlay of the optimized geometry of L tBu Co and L'Co with the crystal structure calculated at BP86/B1 with D3BJ correction (see below for basis set information). Table S2. Crystallographically obtained and DFT calculated interatomic distances (Å) and angles ( o ) for L tBu Co and L'Co. Deviations from crystal structure are shown by Δ and written in parenthesis.
Comparing the geometrical parameters of L tBu Co and L'Co (Table S2) Table S2). This also validates the reliability of the method chosen for optimization. Moreover, the κ-N,η 6 -arene ligation is well reproduced in the initial geometry of the complex, i.e. L'Co. The t Bu substituents on the crystal structure and the methyl replacers in L'Co are almost collinear. This justifies the validity of the truncated model L'Co to be used in place of the actual L tBu Co framework. However, whenever necessary, the original L tBu Co framework has been used for calculations; for example in comparing the calculated and experimental structure for dinuclear carbonate (shown in Table 1, main text) and also for µ-oxo dicobalt species (shown in Table  2, main text).

H. Benchmarking the various density functional methods with DLPNO-CCSD(T).
It is observed from Figure S12 that single point energy calculations on optimized geometries (optimized at BP86(D3BJ)/B1) along the reaction coordinate leading to transformation of 1 to 1' followed by CO2 addition shows variations in the relative barriers depending on the choice of wavefunction based method and density functional. In particular, it is found that B3LYP with and without empirical dispersion correction deviates from the DLPNO-CCSD(T) predicted relative free energies. B3LYP functional predicted that the C-O bond cleavage in the dinuclear dissociative pathway via TS4 is the rate determining step. However, it would indicate a buildup of the dinuclear carboxylate bridging intermediate 10, in contrast to the observed reaction. Repeated experimental attempts could not observe the κ 2 N,N' isomer (1′) of the starting material 1 which indicates that the RDS should be the ligand isomerization. Also kinetics studies showed that the rate does not depend on the concentration of CO2 (as has been discussed in section C of this supporting information), which further emphasized on the experimental identification of only the κN,η 6 -arene isomer. B3LYP on the other hand underestimated the metal-aryl bonding and predicted the κ 2 N,N' isomer to be lower in energy by -10 kcal/mol than the κN,η 6 -arene isomer, in disagreement to the experiments (see Figure  S12). Interestingly, M06L/def-TZVPP delivers a similar barrier for TS1 as DLPNO-CCSD(T)/def2-TZVPP. This observation suggests that the B3LYP functional even with the non-covalent corrections cannot properly describe the metal aryl π-bonding in κ-N,η 6 -arene ligated L'Co, whereas M06L is designed to cover the weak interactions. Hence based on this observation, we predict all the important barriers and reaction free energies with DLPNO-CCSD(T)/def2-TZVPP and M06L/def2-TZVPP in the main text. Figure S12. Gibbs free energy (ΔG) profile for isomerization of 1 (L'Co) to 1' followed by CO2 activation at different levels of theory.

S-21
I. Different bridging binding modes of CO2 with two metal complexes and their geometrical parameters. Figure S13. The two different binding modes of bridging CO2 possible with two metal complexes. Table S3. Binding parameters and geometries of CO2 bridged adducts. Distances in Å and angles in degrees.
The κ 1 C, κ 2 O,nding mode for the S=3 and S=0 (BS33) has been discussed in the main text.  The overlap between nonorthogonal magnetic orbitals of broken symmetry UHF solution has been marked with the notation "S". The subscript represents the number of unpaired electrons locally on a center. From Figure SI4 and Table S4 we find that the high spin septet intermediate with S=3/2 on both the cobalt (II) centers ( 7 10) possesses the same energy as the broken symmetry singlet, 1 10BS33. Also, the ferromagnetically coupled quintet 5 1013 and its corresponding antiferromagnetic triplet counterpart 3 10BS13 have slightly higher energies than the high spin states. Hence, to reduce computational effort, we have chosen the lowest energy septet and the broken symmetry singlet to explore the reaction mechanism. It is to be noted that geometries of the chosen intermediate on these two states are also similar.

O. Comparison of diketiminate bond metrics between 3 and L tBu CoCl
For additional evidence that oxo complex 3 contains cobalt(II) ions, we compared the diketiminate bond metrics of the oxo (right ORTEP diagram) to the known three-coordinate cobalt(II) complex L tBu CoCl (left ORTEP diagram; from ref. 1 of ESI). All ligand backbone bonds are within error between the two compounds, arguing against ligand redox activity within the diketiminate supporting ligand.
S-31 P. Electronic structure of 3 from multi-configurational calculations.  Table S5 shows that the first root of the CAS(20,13) calculation on intermediate 3 has several competing electronic configurations and has significant multiconfigurational character. However, all the CSFs reveal that the oxygen-based ligand orbitals, which are lowest in energy, are never singly occupied; they have predominantly oxygen nπ orbital character. This analysis confirms that the oxygen center has all p-orbitals filled and is best described as O 2-, i.e, intermediate 3 has relatively polar bonds and is best described as Co 2+ -O 2--Co 2+ .
We constructed a larger active space considering the doubly-filled π and vacant π* orbitals of the βdiketiminate ligand in order to allow relaxation of the wave-function and explore excitation between these orbitals. The active space now consists of all metal based d-orbitals, the O centered p orbitals, and two high-lying π orbitals: one on each of the β-diketiminate ligands and their corresponding antibonding π* orbitals. This yielded an active space consisting of 24 electrons in 17 orbitals. In Figure  S18, we have compared how the active space changes before and after calculation of CASSCF (24,17). Importantly, the π* orbitals of the β-diketiminate ligand were essentially vacant with occupation numbers in the range 0.04-0.06. Moreover, the π-orbitals of β-diketiminate ligand are not correlated to the Co-O-Co core ( Figure S18) and the bonding between Co and the ligand N-C-N is primarily ionic in nature as evident from strong localization on the N-C-N motif (80-85%). Thus, inclusion of bonding and antibonding π orbitals of the β-diketiminate ligand does not change the overall description of the electronic structure of the cobalt centers as high-spin Co(II). This justifies the suitability of CASSCF(20,13) calculation for complex 3.

S-35
Q. Frontier orbital analysis during dinuclear associative pathway Figure S19. Important orbitals during C-O cleavage mediated by dinuclear associative pathway.
The dinuclear associative pathway involves heterolytic C-O bond cleavage via TS7. Here charge neutralization around the target oxygen center occurs due to interaction with the nearest metal center. The Lewis acidity of the metal center therefore plays a key role in facilitating the C-O bond scission. The Co-O interaction polarizes the target C-O bond, which may be responsible for the lower barrier pathway in the bimetallic species.

R. Thermodynamics of additional pathways with mononuclear species
The formation of 3 from 6 or 7 in the presence of 1, involving release of a CO molecule in 6 in the former reaction, was computed to have favorable driving forces; however, in the mononuclear dissociative pathway the generation of 6 and 7 has to transvers a high barrier for the C-O bond cleavage. Hence, these two pathways were not considered in the free energy profile in the main text. S-36

S. Why the mononuclear associative pathway is not predicted
CO2 is a linear non-polar molecule featuring very short C-O bond lengths (1.16 Å). Due to the electronegativity difference between C and O, it has partially positive carbon and partially negative oxygen. Activation of CO2 is manifested by bending of the O-C-O angle (~ 130 °) and a decrease in the C-O bond order. To be noted that CO2 has a carbon-localized LUMO and the large HOMO-LUMO gap in linear molecule is not conducive for nucleophilic/electrophilic attack. Bending of the molecule lowers the energy of its LUMO, primarily of the in-plane orbital, and increases the carbon weight in its LUMO and oxygen content in its HOMO. These factors enhance the electron accepting ability of CO2. Hence reduction of CO2 requires a synchronous bi-functional activation. 36 To the best of our knowledge, the mononuclear associative pathway has not been discussed in the literature in the context of a (CO2)2 dimer formation. Agarwal and others have proposed a mononuclear associative pathway consisting of association of formic acid with CO2. 37 In that paper, the presence of a proton on HCOOH ensures strong hydrogen bonding with oxygen of incoming CO2 which should have promoted overall electron transfer process in the bent CO2. However, the dissociative pathway is a more common prediction and there are several examples of both mononuclear and dinuclear dissociative pathways studied theoretically by several groups with the βdiketiminate ligand based metal complexes of Fe, Co and Mn. [38][39][40][41] However, only a single study by Agarwal and others discusses the dinuclear associative pathway for insertion of CO2 into a rhenium carboxylate dimer, [Re(dmb)(CO)3]2, where dmb = 4,4′-dimethyl-2,2′-bipyridine, to release CO and form a carbonate-bridged rhenium dimer. 42 Hence to say associative pathway is a unique mechanistic channel studied in the present report which has very few precedence both for mononuclear and dinuclear cases. In the various CO2 activation pathways that we have summarized above in Figure S20, we found that in each of the feasible pathways, CO2 activation is manifested by cooperative effects of a Lewis acid (like a second metal center) and a Lewis base, electron-rich negatively charged oxo group or a two electron-reduced [CO2] 2-. In intermediate 5, there is only partial electron transfer from Co I to bound CO2 which makes it [CO2] -1 in charge. Therefore, the nucleophilicity or basicity of this bound CO2 is not conducive to promote electron transfer to another incoming CO2. Moreover, the lack of a neighbouring Lewis acid group which might help in bending of the incoming CO2 prior to activation is also a reason that inhibits facile association of a second CO2 unit to the mononuclear adduct 5 and hence a mononuclear associative pathway is not possible. However, during dinuclear associative pathway, the bridging CO2 which is actually [CO2] 2has conducive basicity for nucleophilic attack to The choice of functional can have important influences on the result of a computational study, and therefore we used several different density functionals: BP86 (which is a local gradient corrected functional consisting of Becke's exchange and Perdew's correlation), M06L (the local meta-GGA functional developed by Truhlar and co-workers) and the most popular hybrid functional, B3LYP (consisting of 20% Hartree Fock exchange) to evaluate the electronic structure of intermediate 3. High-spin models using each of these functionals gave quasi-restricted DFT orbitals consistent with a d 7 configuration at cobalt, with 4 electrons in the doubly-occupied region and 3 electrons in the singly occupied region. The doubly occupied orbitals have > 90% Co d-parentage while the singly occupied orbitals have ~ 70-90% contribution from metal-based d-orbitals. This is also true for the broken-symmetry models, where we also found three anti-ferromagnetic coupling (for e.g. small overlap in the range of 0.5 to 0.2 at BP86/def-TZVPP level) between the non-orthogonal orbitals of either metal site, in addition to four uncoupled doubly occupied orbitals. Therefore, analysis of electron occupancy with different HF exchange indicate unanimously the dicobalt(II) oxo electronic structure. The schematic representation of the electronic structure of species 3 (both with S=0 and S=3) along with occupancy and atomic contributions to each orbitals with these three different functionals containing different % HF exchange is provided above in Figure S21.
The MEP plot in Figure S22 reveals distribution of negative electrostatic potential (blue colored isosurface) centrally on oxygen and positive electrostatic potential (red colored isosurface) on some parts of the ligand. The cobalt centers are in the borderline region (green colored isosurface). The MEP mapping clearly shows that the oxygen in species 3 is enveloped by a negative electronic potential and explains the reason for its high reactivity.