Singlet and Triplet Pathways Determine the Thermal Z/E Isomerization of an Arylazopyrazole-Based Photoswitch

Understanding the thermal isomerization mechanism of azobenzene derivatives is essential to designing photoswitches with tunable half-lives. Herein, we employ quantum chemical calculations, nonadiabatic transition state theory, and photosensitized experiments to unravel the thermal Z/E isomerization of a heteroaromatic azoswitch, the phenylazo-1,3,5-trimethylpyrazole. In contrast to the parent azobenzene, we predict two pathways to be operative at room temperature. One is a conventional ground-state reaction occurring via inversion of the aryl group, and the other is a nonadiabatic process involving intersystem crossing to the lowest-lying triplet state and back to the ground state, accompanied by a torsional motion around the azo bond. Our results illustrate that the fastest reaction rate is not controlled by the mechanism involving the lowest activation energy, but the size of the spin–orbit couplings at the crossing between the singlet and the triplet potential energy surfaces is also determinant. It is therefore mandatory to consider all of the multiple reaction pathways in azoswitches in order to predict experimental half-lives.


Contents
Optimizations. In order to investigate the different Z/E thermal isomerization mechanisms, initial guess geometries for Z/E -arylazo-trimethylpyrazole were obtained from conformer searches with crest 2.11 . 1,2 Then, density functional theory (DFT) optimizations were carried out using the ωB97X-D functional. 3,4 This functional was chosen as it was employed in a previous thermal isomerization study of similar heteroaryl azoswitches 5 and a benchmark paper calculating half-lives of Z -azoarenes 6 estimated errors of ca. 2 kcal mol -1 when this functional was used in combination with triple-zeta basis sets. Accordingly, we chose the def2-TZVP basis set. 7,8 In order to replicate the experimental conditions, implicit DMSO (dimethylsulfoxide) solvation using the solvation model based on density (SMD) 9 was employed. The ωB97X-D/def2-TZVP@SMD(DMSO) 3,4,7-9 calculations were carried out using the Gaussian 16 suite. 10 We focused on the conformer 1 (see Figure S1), where the N -methyl group of the pyrazole moiety points away from the aryl moiety in the Z -configuration, and neglected the other conformer (conformer 2) that can be obtained by rotation of the pyrazole ring.
Moreover, these types of photoswitches are known to exhibit chirality, particularly helical chirality, which is usually indicated by the (P) or (M) prefixes. We focus our attention on just one of the enantiomers since they possess degenerate energies and the isomerization rates will likewise be the same. The same level of theory was used for the optimizations of the inversion transition states S6 (TSs: TS iAr , TS iPy ). The T 1 minimum (T min ) was optimized with the Tamm-Dancoff approximation (TDA) 11 based on a restricted Kohn-Sham reference state (S 0 ). According to Ref. 12, it is important to account for multiconfigurational effects to optimize the rotational TS (TS r ). Therefore, the spin-flip time-dependent DFT (SF-TDDFT) approach was used at the ωB97X-D/def2-TZVP@CPCM(DMSO) level as implemented in ORCA5.0. 4,7,8,[13][14][15] The total spin-squared operator,Ŝ 2 , for states 1 to 3 of the TS r was 1.44, 0.78, and 0.12, respectively. States 1 (≈T 1 ) and 2 (≈S 0 ) show spin-contamination between pure triplet (Ŝ 2 =2) and pure singlet (Ŝ 2 =0), which is a common feature in these methods, 16  All optimized geometries (except for MECPs) were verified to be stationary points by frequency analysis (0 imaginary frequencies for minima, 1 imaginary frequency for TSs). The so obtained geometries ( Figure S2) will be referred as DFT geometries in the following and in the main manuscript.
Potential energy scans. Unrelaxed potential energy scans along the Z -TS r -E -(Path r ) and the Z -M 1 -T min -M 2 -E -mechanism (Path rT1 ) were performed. Each scan is an unrelaxed interpolation between the previously DFT-optimized critical-point geometries with 25 single point energy calculations in total for each path at the (TDA-)ωB97X-D/def2-TZVP-@SMD(DMSO) 3,4,7-9 level of theory, as implemented in the Gaussian 16 suite. 10 The results are shown in Figure 3 of the main manuscript and more extensively in Figure S3.   Their results will be denoted as CASPT2 and CASPT2 + , respectively.
The IPEA shift 37 was introduced into CASPT2 38 in order to correct for a systematic error found in CASPT2 when describing dissociation processes. 39,40 This error was described as a general underestimation of the energy of open-shell electronic states. While this general claim seems to not apply to electronically excited states, 41 using the IPEA shift nevertheless corrects the error in the dissociation in the electronic ground state. 37 Since the description of S10 the TS r resembles the electronic structure encountered in ground-state dissociation processes, i.e., a wave function composed mainly of two configurations with similar contributions, it is reasonable to expect the systematic error in non-IPEA-corrected CASPT2 to be present and the energy of TS r to be underestimated as well. Using the IPEA shift with the recommended shift value of 0.25 a.u. 37 should correct for this underestimation. As it will be shown below, the fact that the IPEA-shifted CASTP2 results agree with the MC-PDFT ones, justifies the use of the IPEA in CASPT2 for this case.
MC-PDFT Single Point Calculations. The CASSCF energies were also corrected using hybrid multiconfiguration pair-density functional theory (MC-PDFT) calculations. 42,43 To this aim, the tPBE0 functional (a translated PBE0 18,19 functional), an ultrafine grid, and the PCM(DMSO) implicit solvation model 31,32 was used, as implemented in OpenMOLCAS 21.02. 35 CASSCF and CASPT2 Optimizations. In order to assess the quality of the DFT geometries, we also performed geometry optimization of the critical points using CASSCF and CASPT2 with the same protocol as described above, except that the geometry optimizations were performed in the gas phase. Subsequent single-point calculations were performed using the same solvent models as described previously. Cartesian coordinates of the resulting structures are listed in Appendix A2 and A3. As a starting guess, the DFT-optimized geometries were employed. These calculations were performed using OpenMOLCAS 21.02 35 and the geometry optimization was controlled by the SHARC 2.1 code 17 using the ORCA 5.0 13,14 optimizer.
The CASPT2 optimizations of the Z -isomer and TS iPy did not converge fully using the default convergence criteria in ORCA. Instead, the values listed in Table S1 were used. S11  Figure S4: Orbitals included in the active space of the Z conformer.

Py
Figure S10: Orbitals included in the active space of TS iAr .

S2.1 Energetic Effects on DFT Geometries
Here we discuss the energies obtained from the CASPT2/CASSCF and MC-PDFT/CASSCF calculations as well as the single-reference methods, i.e. the DFT calculations using ωB97X-D, PBE0, PBE0-D3, and B2PLYP as well as MP2. In all cases, the energies were obtained as single point calculations on top of the DFT(ωB97X-D) optimized geometries.
The results are collected in Figure S12, showing the relative energies of the important TS r but also the remaining critical points for comparison. The numerical values are collected in Table S2. At first glance, the relative energy of the TS r (red line) varies substantially   calculations, TS iAr and TS iPy are found at higher energies (6-7 kcal mol -1 ) when using these perturbation methods, even positioning the inversion TSs above the TS r for CASPT2 + . As both TS iAr and TS iPy should be correctly described by a single configuration, the energetic differences stem from the different treatment of dynamical correlation.

S17
The . We are therefore certain that a considerable amount of molecules will isomerize through Path iAr .   Table S3 and the root-mean-square deviations (RMSD) between the atomic coordinates of the differently optimized structures are collected in Table S4. Figure S13 shows the overlap between the three structures. In all cases, it can be seen that the critical points obtained with DFT, CASSCF, and CASPT2 are in good agreement, e.g., with only small differences in the dihedral angle d as well as angles α and α ′ . Most of the RMSD originates from rotations of the methyl groups. The largest deviation (RMSD=0.2Å) is obtained for the T min geometry, due to the relatively flat T 1 energy surface between the MECPs ( Figure S3); however, these larger geometrical differences do not influence the reaction mechanism.
We therefore conclude that, the geometries involved in the potential energy surfaces of the thermal isomerization of PATP do not seem to be particularly affected by the employed level of theory, or in other words, the DFT geometries are robust and in good agreement with the ones obtained at higher levels of theory.    Table S5 and displayed in Figure S14.
Additionally, we compare the DFT(ωB97X-D), MC-PDFT, and CASPT2 + energies against the single-point energies computed at the DFT-optimized geometries (recall Table S2 and Figure S12) in Figure S15. The fact that the DFT and CASPT2 geometries are very similar is translated in very similar energies. Accordingly, when we compare the DFT(ωB97X-D), MC-PDFT, and CASPT2 + energies for geometries optimized with DFT(ωB97X-D) or CASPT2 in Figure S15, we can see that the results within the single point methods are independent of the optimized geometry deviating by only 1.8 kcal mol -1 (DFT), 2 kcal mol -1 (MC-PDFT), and 0.9 kcal mol -1 (CASPT2 + ). The only exception, with deviations of up to 10 kcal mol -1 , is given by the S 0 energy at T min (not shown in Figure S15, see Tables S2   and S5). However, the T 1 energies at the T min only differ by 0.4-2.6 kcal mol -1 . The large difference between the S 0 energies at T min is, thus, a consequence of the flat shape of the potential energy surface of the T 1 state (cf. Figure S3). This feature has no influence on the thermal isomerization mechanism, since it does not involve movement in the ground-state S 0 at the T min . Overall, we can hereby confirm our choice of using DFT geometries to predict the different pathways in the thermal isomerization mechanism of PATP. , and CASPT2 + level of theory. All energies are given relative to their respective Z ground state energy (see Table S5). Dotted lines are only displayed for visual aid.  Tables S2 and S5). Dotted lines are only displayed for visual aid.

S3 Transition State Theory
The basis of TST is the supposition that a hypersurface exists in phase space that separates the space into a reactant area and a product area. The trajectories that traverse this 'dividing hypersurface' from the reactant to the product space will not cross the hypersurface again (the no-recrossing/dynamical bottleneck assumption).

S3.1 Conventional Transition State Theory
In conventional transition state theory (cTST), 44,45 the surface of division is situated at the saddle point (TS), making the net rate coefficient equivalent to the one-way flux coefficient.
According to cTST the rate constant of a thermal Z/E isomerization, and therefore the half-life, depends on the Gibbs free energy difference between the stable Z -isomer and the TS of the isomerization pathway en route to the E -isomer (∆G ‡ ). We thus set to investigate the Z/E isomerization mechanism to obtain the associated rate constant k cTST (T ), assuming a first-order reaction using Eyring equation with the Boltzmann constant k B , temperature T , and Planck constant h.

S3.2 Non-Adiabatic Transition State Theory
In case of the S 0 -T 1 -S 0 rotational transition mechanism, we apply the non-adiabatic transition state theory (NA-TST) 46 to obtain the rate constant k NA-TST (T ). NA-TST expands cTST and in our case places two hypersurfaces in the spin-forbidden MECPs instead of the TSs. Moreover, nonrelativistic interactions that cause transitions between states are of importance. For instance, the greater the spin-orbit couplings at the MECPs, the higher the chance of a transition. 46 NA-TST has provided remarkable agreement between experiment in theory in the case of azobenzene. 47 However, similar to cTST, NA-TST is a theory with a S25 plethora of approximations. Most notably, the no-recrossing (also called dynamical bottleneck) assumption, which states that a reactant reaching the energy barrier to product space can not cross back to the reactant space. We also want to note that NA-TST can be reduced to cTST, by setting the transmission coefficient (γ) to unity and only considering one energy barrier (∆G ‡ ). Finally the k NA-TST (T ) can be calculated, as where γ 1 , γ 2 are the transmission coefficients and ∆G ‡ 1 , ∆G ‡ 2 the Gibbs free energy difference between the stable Z -isomer and the S 0 /T 1 MECPs (M 1 , M 2 ) of the rotational isomerization pathway. 12,47,48 The transmission coefficient γ is calculated as where α is defined as H SO is the spin-orbit coupling, ∆F is the singlet-triplet force difference, with the geometric mean of the singlet and triplet forces at the crossing N is the number of atoms, and the reduced mass µ of the reaction coordinate is given as with m n being the mass of the n th atom.
Finally, we calculate the total reaction rate of the thermal isomerization as the sum of the rates for each of the independent parallel paths.

S3.3 Wigner Tunneling Transition State Theory
With the described version of cTST we neglect the possibility of molecules tunneling through the energy barrier. Because of the relatively large activation energies, a negligible effect on the reactions is expected. However, tunneling might become an important phenomenon for reactions with relatively low energy barriers, since the tunneling probability increases with decreasing barrier height. One possibility to correct cTST 44 reactions rates is to include a one dimensional Wigner tunneling coefficient. 49 We tested the influence of this coefficient using the KiSThelP program 50 (as used by Wang et al. 51 ) on the DFT(ωB97X-D) obtained reactions rates for Path iAr and Path iPy . The differences were negligible with reaction rates of 2.2e-8 s -1 and 1.8e-8 s -1 for Path iAr with and without tunneling, respectively. Similarly, Path iPy has reaction rates of 7.3e-6 s -1 and 6.4e-6 s -1 with and without tunneling, respectively.

S4 Computation of Half-Lives
In  Synthesis of 4-hydroxy-3-(phenyldiazenyl)pent-3-en-2-one (1) Compound 1 (see Figure S16) was prepared according to a modified literature procedure. 52  Synthesis of phenylazo-1,3,5-trimethyl-pyrazole (2) Compound 2 (see Figure S16) was prepared according to a modified literature procedure. 52  the absorption at 339 nm, the E -isomer's absorption maximum, was regularly measured over the course of 15 hours at 23°C. This data was linearized according to equation The linearized absorption data was plotted over time and by applying a linear fit function, rate constant k was gained as slope value ( Figure S18). The half-life of 10.5 days was calculated according to equation whereby κ was assumed to be 1. At a temperature of 23°C the entropic contribution equals -15.6 kJ mol -1 . The enthalpic and entropic parts of activation thereby added up to a total activation energy of 107.2 kJ mol -1 or 25.6 kcal mol -1 , being in the range of the theoretical calculations (see Table 1).
Figure S17: UV/Vis spectra of PATP in the dark (= E -isomer, black line) and after reaching its PSSs for several wavelengths.

S35
Calculation of Z -content of PATP's PSS at 365 nm.
The Z -content of PATP was calculated according to a method by Fischer. 53 With this method, one can calculate PSSs of systems A ⇌ B, in our case E -and Z -isomer, when only A (E -isomer) is known. The experimental data that is necessary for the method by Fischer are absorption spectra of pure species A (equals E -isomer) before irradiation and absorption spectra of PSSs at two different wavelengths. This data was available to us ( Figure S17).
The method by Fischer assumes that the ratio of the quantum yields ϕ E ⁄ϕ Z , ϕ E being the quantum yield of E to Z isomerization and ϕ Z being the quantum yield of Z to E isomerization, does not differ at the two chosen wavelengths, which is generally the case. Therefore, the Z -content of PATP could be calculated for all PSSs that experimental spectra had been measured of. As we were especially interested in the Z -content of the PSS at 365 nm, we used the PSS at 365 nm and 400 nm for the calculations (see Table S6). Thereby we could determine a Z -content of 99% in the PSS at 365 nm and a Z -content of 52% in the PSS at 400 nm. Lasertechnik, set to full power (70 mW) was used for irradiation.

S36
Samples of 50 µM PATP in dry DMSO were used for the measurements and photosensitizer methylene blue (MB) was added in either 1, 10 or 100mol% prior to first irradiation.
After sample preparation and initial recording of a dark-state adapted E -isomer UV/Vis spectrum, the samples were photoswitched by irradiation with 365 nm for 15 seconds. A UV/Vis spectrum was recorded of that state, equaling maximum Z -content. Afterwards, the cuvette with the sample was placed underneath the 660 nm laser (3 mm distance between laser and cuvette top) and irradiated (see Figure S23 for experimental setup). After certain time intervals of 660 nm irradiation, intermediary UV/Vis spectra were recorded. Those time intervals were: 5 min, 5 min (total 10 min), 5 min (total 15 min), 15 min (total 30 min), 15 min (total 45 min), 15 min (total 60 min), 60 min (total 120 min). Those intervals were, if necessary, adjusted to the specific sample.
As can be seen in Figure S25, the Z/E isomerization can be significantly enhanced by providing MB as photosensitizer together with constant irradiation with 660 nm. Moreover, we could observe a strong correlation between the content of MB in the sample and the speed of the Z/E isomerization (see Figure S24 to Figure S26). In addition to those experiments, negative controls were measured. Therefore a 50 µM sample of PATP in dry DMSO was prepared and no MB was added, while keeping the irradiation protocol with 660 nm as mentioned above. We observed an accelerated Z/E isomerization similar to the one observed in the sample with 1mol% MB (see Figure S24 and Figure S27). This meant that irradiation with 660 nm led to some increase in the speed of the Z/E isomerization, however not comparable to the effect observed when we applied 10mol% or 100mol% of MB. This also meant, that addition of only 1mol% MB had a neglectable effect on the Z/E isomerization, as the observed faster relaxation for that sample was mostly induced by the 660 nm laser irradiation itself. Additionally, one sample was prepared with 1mol% MB that was not irradiated with 660 nm after initial photoswitching towards the Z -isomer with 365 nm. The half-life of that sample was determined according to the procedure mentioned in Section S5.2 (see Figure   S28). The observed half-life was 10.8 days and thereby very close to the determined half-life S37 without a photosensitizer. As our earlier results had already foretold that 1mol% MB hardly influences the Z/E isomerization even with 660 nm irradiaton, we additionally prepared a sample with 100mol% MB that was kept in the dark for 30 minutes after the initial irradiation with 365 nm. Without irradiation with 660 nm, no significant Z/E isomerization could be observed in that sample during that time interval (see Figure S28).
Finally, two key experiments were repeated with 50 µM PATP samples in degassed (=oxygen-free), dry DMSO keeping the sample solutions as inert as possible by establishing argon layers. The repeated experiments were the ones with addition of either 1 or 10mol% of MB. Thereby, we sought to eliminate potential quenching of the triplet state by oxygen in the solution. Under oxygen-depleted conditions and addition of 10mol% MB, we observed an even more accelerated Z/E isomerization than under ambient conditions (see Figure S30).
This showed that the ISC-enhanced Z/E isomerization benefits from oxygen-depleted conditions due to less triplet state quenching but is not limited by the presence of ambient oxygen. The oxygen-depleted experiment with 1mol% MB on the other hand showed very similar behavior to the sample with 1mol% MB under ambient conditions (see Figure S31). This might be due to some residual oxygen in the solution or disruptions of the argon layer during the experiments. As the overall concentration of MB in the 1mol% MB experiment equals 0.5 µM, such traces might be enough to still cause photobleaching. The MB concentration dependency of the Z/E isomerization acceleration was again highlighted as well as that 660 nm irradiation alone has a significantly lower effect on the Z/E isomerization than irradiation together with 10mol% MB.
S38 Figure S23: Experimental setup for continuous 660 nm irradiation of the sample solution in the cuvette (here empty).