Effect of Electron Donating/Withdrawing Groups on Molecular Photoswitching of Functionalized Hemithioindigo Derivatives: a Computational Multireference Study

Hemithioindigos are a versatile group of molecules with the ability to undergo double bond isomerisation upon light excitation. Using electronic structure theory, a detailed description of the photoisomerisation of hemithioindigo and five of its functionalized derivatives is reported. We use a combination of Density Functional Theory (DFT), approximate second order Coupled Cluster (CC2), and multireference methods such as the complete active space self-consistent field (CASSCF) method to probe how isomerisation capabilities of functionalized hemithioindigo derivatives are affected by functionalisation with both electron donating and withdrawing groups. We demonstrate that important properties such as the thermal barrier, electronic absorption spectra and, to an extent, the mechanism can all show pronounced differences compared to the parent hemithioindigo molecule. Our results demonstrate how the inclusion of functional groups can have important effects on its photoisomerisation capabilities, which in turn, can influence the fate of these molecules upon light excitation and their utility as molecular switches.


Introduction
The ability to reversibly modify the properties of a material by exposure to external stimuli such as electric fields, light, and temperature changes lies at the heart of many modern materials applications. This includes transistors, sensors, readwrite memory and even drug delivery systems. [1][2][3][4][5] A central goal of molecular nanotechnology is to integrate this ability to the level of single molecules. Molecules that can be reversibly switched between two measurably different states -through isomerism or other intramolecular reactions -are called molecular switches. [6] Whereas nature has produced many examples of responsive molecular switches and machines controlled by temperature, pH, or light; the artificial synthesis and design of molecular switches that can be controlled externally is a challenging field of chemistry. Molecular switches are a class of molecules that has been of significant technological interest over recent decades due to their potential applications in drug delivery, [7] host-guest chemistry, [8] and molecular logic. [3,9] While molecular switches can be activated with a variety of different stimuli, photoswitches -molecules that can reversibly respond to light -have garnered most attention. [10] Examples include diarylethenes, [11] spyropyrans, and azobenzene, [6] the latter being the most studied representative of this class of molecules. [12] In azobenzene, the irradiation of light is used to induce double bond isomerisation between an E and Z isomer. The molecule and its derivatives have been widely applied to create on-surface switches and switchable polymer thinfilms. [13][14][15][16][17] While the mechanism of double bond isomerisation around the diazo double bond in azobenzene is relatively well understood, [18][19][20] it involves certain practical disadvantages. For example, while the Z to E transition can be activated by visible light, the E to Z transition requires ultraviolet (UV) light. [18] This is not ideal as UV light can lead to molecular degradation and is non-penetrative to most media. [21] This is particularly important for medicinal applications and for ensuring that the molecular switch will not degrade over time as a result of frequent switching events. [22] The photophysics and photochemistry of molecular switches can be modified by chemical functionalisation. For example, attempts have been made at constructing azobenzenes with spectrally red-shifted absorption through extending and delocalising the π-system by placing amide groups at the para positions of the phenyl rings. [23] This was effective for enabling switching in response to visible light, but also compromises the thermal stability of the switch, reducing the thermal reisomerisation barrier from the Z to E state. [22] Hemithioindigo (HTI) photoswitches have more recently garnered attention, as they do not suffer from some of the challenges encountered by azobenzene-based photoswitches described above. [24] The structure of HTI is shown in Figure 1, with the thioindigo section to the left of the central carboncarbon and the hemistilbene section to the right. Unlike azobenzene, HTIs can reversibly switch from Z to E and E to Z by absorption of visible light. [25] HTIs and their derivatives also display higher ground state energy barriers providing additional thermal stability when compared with azobenzene photoswitches. This helps to ensure that switching is a result of irradiation by light rather than thermal activation. [26] This may also benefit the application of HTIs in cases where interaction with the environment (e. g. surface adsorption) can reduce the groundstate barrier. [27] The synthesis of HTIs is welldocumented, [28] and shows that functionalized derivatives of the parent molecule can be made with relative ease which is an important prerequisite to streamline functional design of light-responsive nanomaterials. [29] As HTIs can be easily functionalized, they represent a versatile class of molecular switches for the design of compounds with tailored excitations and photodynamics.
Unlike azobenzenes, the central photochromic moiety of HTI is a carbon-carbon double bond, allowing for additional substituents to be present in the molecule due to the additional coordination capabilities of the carbon atoms compared to nitrogen. The change in composition of the switching moiety also is able to directly influence the mechanism of the isomerisation; for example, in azobenzenes the dominant reaction coordinates are dihedral rotation and planar inversion. [30,31] However even for ethylene, the simplest carbon-carbon double bond containing compound, the photoisomerisation cannot be suitably described by a simple angle bending or torsional motion, and requires a more thorough investigation across multiple reaction coordinates (for example including out-of-plane pyramidalization of the sp 2 carbon atoms) and across multiple electronic excited states. [32,33] Theoretical descriptions of the excited state chemistry of the parent HTI have been previously reported by the groups of Dreuw [34] and de Vivie-Riedle. [35] In the study carried out by the latter, multireference quantum chemistry methods are applied across several reaction coordinates, yielding a mechanism that suitably includes conical intersections between the electronic excited and ground states, which was not captured with the single reference description presented by the former. These previous studies have identified the mechanism as being predominantly based on dihedral rotation, yet the molecule also must undergo some early secondary structural distortion (e.g pyramidalization or tilt) in order to reach conical intersections which mediate the photoisomerisation reaction. Whilst extensive experimental studies on functionalized derivatives of HTIs have been reported, illustrating the effect of functionalisation on the rates of switching, there have only been few theoretical studies dedicated to describing the mechanistic origin of these experimental observations. [36][37][38] Cordes et al. observe in their experiment that placement of electron donating groups (EDGs) at the para position of the hemistilbene section leads to an increase in the rate of switching compared to the base compound. [38] Upon replacement of an EDG with an electron withdrawing group (EWG) the opposite effect is observed, namely a decrease in switching rate. [37,39] While it is clear from experiment that the introduction of EDGs and EWGs affects the rate of photoswitching, the underlying mechanisms of this observation are not fully understood.
In this study, we perform a computational study on the effects of para-substitution on the hemistilbene moiety of HTI, varying only a single substituent with our attention fixed on typical EDGs and EWGs in an attempt to build an understanding of how a single functional group contributes to the changes in rate observed in experiment. We start by characterizing the photoisomerisation mechanism of the HTI base molecule in order to derive key design parameters that define the ground state, excited state, and mechanistic details of photoisomerisation for this class of molecule. We find that the isomerisation of the parent molecule is connected to a S 1 ! S 0 electronic transition, which induces a motion that couples dihedral torsion and pyramidalisation, facilitating internal conversion to a dark excited state S 2 . This internal conversion subsequently drives the isomerisation transition before radiationless decay to the ground state S 0 via a conical intersection. Based on this finding, which is in agreement with previous work [35] , we use the same method to investigate how these properties change in response to being functionalized by varying substituents. We introduce dimethylamino and methoxy groups to act as the EDGs (see 2 and 3 in Figure 1) and chloride, cyano, and bromide groups to act as EWGs (see 4-6).
For the EDG-functionalised case, we find energy landscapes that are consistent with a mechanism requiring less pyramidalisation and being more dominated by dihedral rotation around the central double bond. For the EWG-functionalised molecules, the opposite effect is found. We further discuss how these changes in excited-state dynamics connect to the experimentally observed rate changes.

Ground state optimizations
Initial optimizations of equilibrium geometries for all molecules were performed at the B3LYP/def2-TZVPP level, [40,41] including the D3BJ dispersion correction of Grimme et al. [42] These calculations were performed with ORCA (version 4.1.2), [43] and all utilise the resolution of identity (RI) approximation using the def2-TZVP/C and def2/J auxiliary basis sets in order to improve computational efficiency. [44] For calculations of minimum energy paths in the ground state, the nudged elastic band (NEB) and its climbing image (NEB-CI) variant were used. [45,46] The software package FHIaims (version 200112) was used to calculate energies and forces across the minimum energy path, [47] which was optimized using the NEB implementation in the Atomic Simulation Environment (ASE). [48] For these calculations we use the B3LYP functional, accompanied by the tight numerical basis set.

Calculation of excitation energies
Vertical excitations were calculated with the following methods: time dependent density functional theory (TD-DFT), second order approximated coupled cluster CC (2), and the algebraic diagrammatic construction ADC (2). [49,50] For each of these methods, the def-TZVP basis set was employed throughout. The Tamm-Dancoff approximation was applied in all TD-DFT calculations, [51] for which the functionals used were B3LYP, in addition to the Coulombattenuated variant CAM-B3LYP. [40,41,52] Solvent effects were captured implicitly using the COSMO model. [53] Calculation of 2D potential energy landscapes 2D potential energy landscapes were constructed by performing constrained relaxations across two coordinates of interest -chiefly dihedral rotation and pyramidalization in 10 degree intervals. Using state-averaged complete active space self consistent field (SA-CASSCF) with an active space of 14 electrons in 12 orbitals, geometries for each state were optimized. All calculations also employ the 6-311G* basis set of Pople et al. [54] Further n-electron valence state pertubation theory (NEVPT2) calculations were carried out subsequently to include dynamic correlation effects. [55] In all cases, the strongly contracted NEVPT2 (SC-NEVPT2) formulation was used. For molecule 1 we average across three states; S 0 , S 1 , and S 2 , however for the remainder of the molecules 2-6 an additional state is included, in order to ensure that the electronic character of each state is consistent throughout. Without the inclusion of the additional state, higher lying excitations with low oscillator strength emerge at some geometries across the PES. Input and output files of all computations have been uploaded to the NOMAD repository as a data set. The data can be accessed via following URL: https://dx.doi.org/10.17172/NOMAD/2021.12.16-1.

Active space selection
We conclude this section by justifying the active space used in our calculations. As a first diagnostic for the suitability of the active space, MP2 unrelaxed density calculations for the molecule across the dihedral rotation coordinate were performed. Any unoccupied orbital which shows partial occupancy > 0.05 was included within this active space. Based on previous work in the literature, we choose to include a conservatively large active space, as is used in previous studies, ensuring that the lone pair on the oxygen atom is included as it has been reported that its involvement is an essential contribution to the S 1 ! S 0 excitation. [35] Other active spaces in the literature for the parent molecule utilise an active space of 8 electrons in 8 orbitals; however this active space was only used in order to characterize the ground and first excited state. [56] Attempts to use this smaller active space within our study across three states (ground state, first and second excited states) led to convergence issues across the different reaction coordinates. For this reason we choose to opt for the (14,12) active space for which we consistently achieve convergence across the explored coordinate space. Further details concerning the active space selection are discussed in the supplementary information (Part A.)

Results and discussion
Thermal stability of HTI compounds degrees means that no pyramidalization of carbon atom 3 has occured. Values above and below 180 degrees indicate pyramidalization and partial rehybridisation around carbon 3 from sp 2 to sp 3 . iv. ϕ: Torsional angle of phenyl ring twisting of the hemistilbene section of the molecule. While this coordinate does not typically directly participate in the isomerisation, it can be indicative of hula-twisting motion of the molecule. [57] Hula-twisting refers to the motion achieved by translocating the CÀ H bond which results in the rotation of the two connected single and double bonds. [58] Values for ϕ of 180 or 0 refer to a geometry where hemiindigo and hemistilbene ring systems form parallel planes with respect to each other. A ϕ value approaching 90 or 270 degrees indicates that the two aromatic ring systems within the molecule are perpendicular. The ground state minimum energy path between Z and E corresponds to the barrier associated with thermal reisomerisation. This minimum energy path can be calculated with the NEB method. [46] This, followed by a CI-NEB yields the first order transition state, and the path traversed to reach it. Our results on the energy difference between Z and E and the thermal barrier are reported in Table 1.
For all molecules in the study, the Z isomer is more stable than its respective E counterpart, at least within the electronic ground state. This is in contrast to azobenzene, where E is more stable. [59] This is due to the increased sterical strain of the phenyl ring in HTI as it approaches the carbonyl group which destabilizes the E arrangement, creating an energy difference DE Z=E of about 0.2 eV between the two equilibrium geometries. Table 1 displays the E-Z energy differences of the different molecules in their equilibrium geometries. The inclusion of substituents in any case does not seem to have an immediate effect on the relative energies between the two stable states; however they do appear to affect the thermal isomerisation barrier DG y also reported in Table 1. The introduction of the dimethylamine and methoxy group in (2) and (3), respectively, increases the thermal stability of the E isomer by 0.31 and 0.46 eV compared to the base compound (1) as calculated at the B3LYP/'tight' level of theory. In EWG-functionalized molecules, we also note higher thermal stability; particular in the case of the halogen-containing molecules 4 and 6, which show increases of thermal stability by 0.81 eV and 1.45 eV respectively. A less pronounced effect is shown in the case of 5 which shows a barrier 0.48 eV higher than that of 1.
The trend in barrier height shown in DFT can be reproduced at the CASSCF(14,12)//NEVPT2 level, where the geometries are optimised at the CASSCF level. Seemingly, placing any substituent, whether it be electron donating or withdrawing, will result in an increase of the energetic barrier in the electronic ground state. From this we can judge that thermal reisomerisation after Z!E photoisomerisation will likely be slower in functionalized HTIs with respect to the parent compound, as a result of the increased thermal stability. Key geometric details of each molecule (as denoted in Figure 1 in their equilibrium and transition states have been tabulated in Table SVII within the supplementary information.

Photophysics of HTIs
Previous studies have shown that the excitation of 1, which initially causes the photoisomerisation, arises from the lowest electronic transition, the π-π* transition. [39,60] This is verified through steady state irradiation of the molecule at the longest wavelength maxima observed in the absorption spectra, which yields the growth of the opposite isomer. [26] The second lowest excitation is that of n-π* character, which is a dark state with oscillator strength ω osc of virtually zero. Table 2 reports theoretical excitation wavelengths for vertical excitations and oscillator strengths associated with transitions using TD-DFT, CC(2) and CASSCF(14,12)/NEVPT2. The TD-DFT results correspond to calculations carried out with a B3LYP functional in addition to an implict solvation model to capture the conditions of the experiment with which we compare. [39,53] The same implicit solvent correction was also included to simulate absorption spectra at the CC(2) level of theory (shown in Figure 2). We find that CC(2) provides overall better agreement with experiment than TD-DFT at the B3LYP or CAM-B3LYP level. The CC(2) absorption wavelengths are within 15 nm from experiment across all systems for which experimental data exists. Calculated excitations from NEVPT2 have also been included in Table 2, yet this method struggles to accurately predict accurate excitation energies for this class of molecule, which was displayed in a previous study. [35] We were also unable to apply the same solvent correction in our multireference calculations which was used in the TD-DFT and  Table 2. Calculated vertical excitation energies λ (in nm) and corresponding oscillator strengths ω osc for the first two singlet excitations in molecules 1-6.
As can be seen in Figure 2a-c, by placing EDGs such as NMe 2 or OMe on the hemistilbene section (in the para position in both cases -see Figure 1), a spectral red-shift of the excitation wavelength can be observed. This effect is stronger for NME 2 in 2 than for OMe in 3 and is due to the extension of the π-system encountered when adding groups capable of adopting resonance forms. Conversely, upon the introduction of EWGs (4-6) a spectral blue-shift is instead observed. In 4 the inductive effect felt from the electron withdrawing Cl will pull electrons out of the π-system, which further stabilizes the HOMO, thus decreasing the excitation wavelength of the optical transition. In 5 we observe a competing effect from the strong electron withdrawing inductive effect against the extension of the π-system through resonance forms, resulting in a negligible shift of S 1 ! S 0 excitation wavelength overall. In order to achieve efficient photoconversion rates between Z and E isomers, the absorption maxima of the two should ideally be well separated. In the case of the parent compound, the maxima as predicted by CC(2) are separated by Dl max = 19 nm. The wavelength separation of the absorption maxima is increased in all functionalised molecules to between a Dl max value of 27 nm and 30 nm, respectively. Another interesting observation is that the wavelength of the S 2 ! S 0 excitation across the three EWG-containg molecules 4-6 is significantly blue shifted compared to the parent compound.
We can also see in Figure 2 that the introduction of the substituents has a remarkable effect on oscillator strength exhibited by the functionalized molecules. Through inclusion of the EDGs, one can observe approximately a five fold increase in oscillator strength as shown when comparing 1 and 2. In order to make photoswitching as responsive as possible, the oscillator strength should be maximised if and whenever possible in order to ensure the light drives the reaction reliably. The S 1 ! S 0 excitation is chiefly composed of a HOMO!LUMO transition. The HOMOs of EDG functionalized molecules show a higher spatial overlap with the corresponding LUMO, as shown in Figure S4 in the supplementary information where the frontier molecular orbitals for each molecule are visualized. The extension of the π-system caused by substitution in molecules 2-3 cause both of these frontier orbitals to also become extended, resulting in a higher overlap between HOMO and LUMO for these molecules when compared with the parent compound. On the other hand, EWG-containing molecules 4-6 are fairly similar in their S 1 ! S 0 transition strength when compared with the parent molecule 1.

Mechanism of HTI photoisomerisation
It is possible to describe geometric transformation from one isomer to another by sampling over a single internal degree of freedom of the molecule; for example, dihedral rotation around the ω torsion angle or via an in-plane inversion around the α angle (see Figure 3a and b). Monitoring the potential energy landscape across these reaction coordinates gives a good indication into how the molecule behaves during the isomerisation. In reality, however, sampling a simple one dimensional Figure 2. Simulated absorption spectra for HTI molecules 1-6 calculated at the CC(2)/def2-TZVP level with implicit solvation for CH 2 Cl 2 . Calculated excitations have been artificially broadened with a Lorentzian fitting using a FWHM of 50 nm. Experimentally reported wavelengths have also been indicated for molecules with vertical bars where available. [39] reaction coordinate is likely to be quite different from the actual mechanism observed in experimental studies. Particularly, for carbon carbon double bond rotation, it is known that multiple internal degrees of freedom need to be considered.
The most commonly modelled reaction coordinate in HTI is the dihedral rotation across the central carbon-carbon double bond. It is generally accepted in the literature that this is the dominant coordinate in operation. In one of their earlier experimental studies, Cordes et al. find through transient electronic absorption spectroscopy (TEAS) that upon excitation, the structure relaxes extremely quickly within 0.3 ps during which a spectral red-shift is observed. This is followed by a significant decrease in oscillator strength, which is indicative of a slight structural distortion caused by reaching a state with significant charge transfer character. [61] The exact structural distortion was not actually determined at the time. Numerous mechanisms have been proposed in literature; including what seems likely to be a combination of dihedral rotation and pyramidalization (see Figure 3(c)). [35,62] It is worth noting that the pyramidalization coordinate itself is not sufficient to complete the transformation between states, but rather facilitates an internal conversion to a different electronic state, which is followed by dihedral rotation. [63,64] Pyramidalization is a common coordinate studied for carbon-carbon double bond rotations as it seems that even in ethylene, a pure onedimensional dihedral rotation does not exhibit a conical intersection between the first excited state and the ground state to enable a radiationless transition. [65] The pyramidalization coordinate has also been found to be important in the photoisomerisation mechanism of stilbene. [66] Another reaction coordinate often discussed for some HTI compounds is that of the "hula-twist" (d). This mechanism commonly occurs in the isomerisation of sterically hindered molecules (for example by functionalisation or in solvent), in an attempt to conserve volume and can be thought of as a cooperating rotation and inversion only in the plane of the C=C bond. [57,67] Such a pathway has been previously reported for azobenzene in solution, where the isomerisation can occur through a rotation of the N 2 moeity with both phenyl rings fixed in place. [68] HTI has shown tendencies to form hula-twist products, however this occurs in extremely low yields and can be notably influenced by temperature. [69] A potential energy surface for the hula-twist in the excited state has been calculated previously, but it was reported that the energy increased too steeply in this coordinate to be feasible. [34] Another prevalent coordinate used to monitor C=C isomerisation is that of the stretching of the isomerisable bond itself. This is most notably shown by the isomerism of retinal studied by Olivucci and co-workers, where they demonstrate that an elongation of this C=C bond (in addition to the overall skeletal expansion of the molecule) is needed as a complementary motion to the torsion coordinate to access the conical intersection. [70,71] This coordinate, important in molecules with a conjugated carbon backbone such as retinal, does not appear to be relevant to enable non-radiative decay in the more compact, ring-containing structure of HTIs, which has been illustrated in Figure S6. Here we show that a relaxed scan in both the C=C stretch and rotation coordinates does not seem to create any conical intersections or avoided crossings between the S 1 and S 2 states.
Previous multi-reference quantum chemical calculations have been carried out by different groups which take into account the reaction coordinates displayed above in Figure 3. [35] Nevertheless, a conclusive picture has not yet been reached. Zhu et. al. report a detailed joint experimental and computational study for the isomerisation of 1. [56] In their study, the authors have only probed the first singlet excited state, proposing to additionally calculate triplet states in order to determine whether intersystem crossing is likely to occur. However, carbon-carbon double bond photoisomerisation requires the consideration of multiple crossing excited states across a multi-dimensional energy landscape; [32] in this study we revisit the potential energy surfaces of the ground and the first two excited states of HTI. Figure 4 depicts the ground state and first and second singlet excited state energy landscapes of 1 along the pyramidalisation and dihedral rotation coordinates θ and ω computed with CASSCF(14,12)//NEVPT2. Roman numerals denote key events across the isomerisation path, these points of interest have been indicated on the contour plots of each Figure which illustrate the approximate location of the following: i. Franck Condon region of molecule in Z state immediately after S 1 ! S 0 excitation ii. Local basin in S 1 associated with Z!E isomerisation iii. Basin in S 2 which forms conical intersection with S 0 iv. Franck Condon region of molecule in E state immediately after S 1 ! S 0 excitation v. Local basin in S 1 associated with E!Z isomerisation The first of these (i) is shown on the S 1 surface and corresponds to the Franck-Condon or vertical excitation (S 1 S 0 ) of the molecule in the Z state. The curvature of the S 1 state induces a pyramidalisation motion to a metastable state (ii) at a pyramidalisation angle of about q �20°and a dihedral angle of w �50°. This is consistent with the location of the transition state calculated by de Vieve Riedle and colleagues. [35] At (ii), the S 1 to S 2 transition can occur. As shown in Figure 4, the S 2 state becomes lower in energy than S 1 towards increasing ω, leading to a dihedral rotation towards the w ¼90°p oint (iii), where the radiationless transition to the ground state occurs with subsequent relaxation to the E state (iv). Note that minor pyramidalisation is only required to access point (ii) on S 1 , whereas the dynamics on the dark S 2 surface lead to a rehybridisation back to q �0.
Considering the reverse E!Z reaction, the photoisomerisation mechanism will likely involve a similar process although the S 1 does not feature a similarly prounounced local minimum as point (ii). Photoisomerisation from the E state will therefore require little to no pyramidalisation to reach the dark S 2 state, which intersects with the S 1 state at a dihedral angle of about 140°. From this point, relaxation on the S 2 surface into the region of the conical intersection at (iii) can occur to perform the full switch in reverse.
The photoisomerisation of the parent HTI compound 1 therefore involves a minor pyramidalisation of carbon atom 3 (cf. Figure 1), followed by dihedral rotation around ω. The CASSCF(14,12)//NEVPT2 S 1 energy landscape suggests that the dynamics to reach the S 1 /S 2 conical intersection may have to overcome a shallow barrier. The pyramidalization step is likely to be extremely short lived and difficult to observe experimentally, as has been shown for other compounds e.g nucleobases [72] and in other HTI derivatives where the geometric changes in the excited state are the rate limiting step. However, the dynamics associated with excited-state motion display lifetimes at the sub-picosecond timescale so exact structural distortions remain difficult to detect and characterize. [36] Experimental measurements for molecules of this class state that the forward Z!E isomerisation reaction proceeds slower than the reverse E!Z reaction. [35,38,39] In the parent molecule 1 this has been reported to be t Z=E = 38 ps and τ E/Z = 21 ps, however the decay of signal for the E!Z reaction takes 2.3 ps. [39] Furthermore, some functionalized derivatives, perform the Z!E reaction significantly slower. With only minor functionalization -adding one methyl group -the observed time constants for switching are measured at τ Z/E = 216 ps and τ E/Z = 10 ps. [37] For this reason, in the following sections, we will study the effect of functionalisation on the excited-state pathways.

Effect of electron donating groups on the mechanism
We have also calculated the ground-and excited-state energy landscapes for 2 and 3, which are derivatives that feature an EDG. In both cases, a dark state with low oscillator strength influences the potential energy surface. This was remedied by including an additional state which we believe to be of little importance to the isomerisation mechanism. We include this state to create complete separation from the S 3 state which does seem to cut through some high energy sections of the potential energy surface of the S 1 and S 2 states.
As can be seen in Figure 5, the internal conversion between S 1 and S 2 states for 2 occurs at lower values of ω compared to compound 1. In this case, no pyramidalization is required to reach a crossing of the two surfaces. This occurs at a dihedral angle of approximately 20°. The subsequent dynamics on the S 2 state lead via dihedral rotation to the S 0 /S 2 crossing at w �90°. As observed for 1, after proceeding back to the ground state via the conical intersection, relaxation across the S 0 surface to yield the E isomer at point (iv) can occur. Experimental studies find the Z!E isomerisation of this molecule to complete in a timescale of τ Z/E = 10 ps, which is significantly faster than what is measured for 1. [37] When considering the backwards E!Z reaction, we again observe there to be very little to no pyramidalization required in order for S 1 /S 2 internal conversion to occur. Instead this appears to happen very close to the Franck-Condon region, with a ω value of � 160°. Experimental studies report a τ E/Z = 2.3 ps, which again is extremely fast, but shows no further increase in rate when compared to the base compound 1. Figure 6 displays the potential energy surfaces for 3, which appears to behave in a similar manner to that of 1, with a clear presence of a basin which facilitates the conversion from the S 1 surface towards the S 2 which features a conical intersection with the ground state S 0 . This basin is observed at a dihedral angle of 20°and pyramidalization angle of 20°, whereas crossing between the two states occurs at ω = 40°. Experimental studies report a time constant of τ Z/E = 3.6 ps which is an almost ten-fold speed-up compared to 1. [39] In contrast to 1, both 2 and 3 feature a pathway for relaxation on S 1 that is barrierless, and the vibrational relaxation is likely to occur quickly, facilitating swift internal conversion into the S 2 surfaces. In addition to this, the gradient of the S 2 surface in 3 towards the conical intersection with the S 0 is significantly steeper than that in 1, which is likely to cause the pronounced increase in rate.
The reverse reaction in 3 seems to proceed without the need for any pyramidalization, and an S 1 /S 2 crossing can be observed at around ω = 120°. This crossing point is easily accessed by the excited molecule in the S 1 surface, from which the reaction proceeds via the same S 2 /S 0 conical intersection, followed by relaxation to the Z equilibrium geometry. The experimental measurement determines that this process occurs with τ E/Z = 0.9 ps.
As aforementioned, barriers in the S 1 state grow significantly in the EDG functionalized HTIs 2-3 compared to the parent compound 1. This is likely due to additional conjugation and resonance that results from inclusion of the EDGs, yielding a higher energy penalty being incurred upon distortion of the molecule during pyramidalization. It is also likely that this large

ChemPhotoChem
Research Articles doi.org/10.1002/cptc.202100290 energetic penalty observed in S 1 is the reason that dark states appear for 2 and 3, as the higher energies observed in the S 1 state approach energies of other higher-lying states, even in some cases showing higher potential energies than the S 3 state.
Pyramidalization acts as a coordinate to facilitate the S 1 !S 2 conversion in both 1 and 3, which in the case of 2 is not necessary due to the already lower energy surface of the S 2 compared to that of the S 1 . NMe 2 is a more potent EDG than OMe; this seems to promote traversal across the dihedral rotation pathway and allows it to isomerise without any need for preliminary pyramidalization, which is reflected in the reduction in time constant. No pyramidalization is required to perform the E!Z isomerisation in all three molecules 1 to 3, therefore only subtle changes to the rate are observed upon introducing an EDG.

Effect of electron withdrawing groups on the mechanism
Experimental studies have observed an acceleration of switching -particularly in the Z!E isomerisation -upon inclusion of EDGs. The opposite effect is displayed in studies for derivatives including EWGs. [36,37,39] Figure 7 displays the potential energy surface for the Clcontaining derivative 4. Upon vertical excitation into the S 1 of the molecule in the Z state, relaxation of the molecule into the nearby basin can occur. However, in order to perform internal conversion to the S 2 state, there is a modest barrier that must be overcome before a crossing between the two surfaces occurs. The crossing seam includes the point marked (ii) which lies at w �30°, and q �25°. Once on the S 2 surface, relaxation into the minimum on this surface (iii) can occur. Similar to all previously described molecules in this study, this minimum is located at w �90°. Contrarily to molecules 1-3, this minimum is located at higher values of pyramidalization, q �25°. From here, transition into the S 0 can begin and relaxation along this surface brings about the switch to the E isomer.
Concerning the reverse E!Z reaction, there are two distinct features that can be observed that are not present in the molecules 1-3. The first concerns the energies of each excited state at the Franck-Condon region (iv): the S 2 surface remains lower in energy than the S 1 . Secondly, there does not seem to be a crossing between the two excited states in the space of the two coordinates at any value of ω greater than 90°. To access the conical intersection that facilitates the isomerisation, other internal degrees of freedom might be required. A relatively slow internal conversion may have to compete with fluorescence or possibly a route involving intersystem crossing to a triplet state, resulting in a slower transition and overall loss of quantum yield. The minimum in S 1 has been denoted by (v) on the potential energy surface.
The potential energy surface for 5 is displayed in Figure 8. Excitation into the S 1 at the Franck-Condon region enables vibrational relaxation into a crossing point (ii) which is accessed by relaxation into the pyramidalization coordinate, which displays a comparatively steeper gradient to the minima over the rotation coordinate. This crossing point between S 1 and S 2 exists at a geometry with w �20°and q �30°. In the S 2 , traversal across the dihedral coordinate into the large minima displayed at w �90°can occur. In 5, similarly to 4, the minimum in the S 2 is located at higher values of θ, which appear to be stabilised by the presence of EWGs. A conical intersection between S 2 and S 0 enables the return to the ground state S 0 , from which relaxation to the equilibrium geometry for the E-state at (iv) can occur.
The E!Z reaction occurs through excitation into the S 1 which remains lower than the S 2 around the Franck-Condon region. Relaxation of molecule 5 along the potential surfaces requires distortion in both coordinates until a crossing point is reached at w �110°, q �20°. Following internal conversion to S 2 , relaxation into the minimum here allows traversal into S 0 through the conical intersection, which yields the Z isomer following relaxation in the ground state.
The final molecule that we study, 6, displays features akin to the other halogenated derivative 4, although the effects of the electron withdrawing capabilities of the Br are modest compared to the more electronegative Cl. Figure 9 shows the potential energy surfaces for 6, which again begins with excitation into the S 1 , followed by relaxation into the pyramidalization and rotation coordinates where a crossing point between S 1 and S 2 is encountered; with w �20°and q �20°. Like both the previously described EWG-containing derivatives, the minima in the S 2 surface occurs at higher values of θ compared to that of molecules 1-3. Relaxation into this minimum enables the conical intersection to be accessed at point (iii), with w �90°, q �20°. From here, the full Z!E isomerisation can be completed by relaxation across the rotation coordinate into point (iv).
In the reverse reaction, molecule 6 shows a S 1 /S 2 crossing point at (v), which was not observed for molecule 4. Light excitation of the molecule in the E-state leads to excitation into the S 1 state. This time, there is a radiationless pathway into the S 2 surface which occurs approximately at w �120°, with q �25°. In this molecule in particular, the stabilizing effect of the EWGs on the S 2 state at higher values of θ is clearly visible. Only at low values of θ does the S 2 surface appear at higher energies, however upon introducing the pyramidalization to the molecule, this state becomes much more favourable. The conical intersection between the S 2 and S 0 is displayed at w �90°, q �20°.
The inclusion of EWGs seems to cause the molecule to access higher values of pyramidalization angles, which become necessary to facilitate the conversion in the larger barrier experienced in the rotation coordinate, particularly in the S 1 . Dihedral rotation occurs with significantly smaller barriers in the case of EDGs, allowing them to isomerise at extremely fast rates compared to 1. The reverse effect appears to be observed for EWG-containing molecules, which cause the molecule to energetically prefer the pyramidalization, which results in a delayed switching from Z!E due to the more convoluted path

ChemPhotoChem
Research Articles doi.org/10.1002/cptc.202100290 required to be taken by the molecule across the ω coordinate of the potential energy surface.

Conclusion
Using multireference electronic structure calculations we have studied the photoisomerisation of hemithioindigo derivatives. Our characterization of the parent compound 1 is in good agreement with previous theoretical descriptions, and builds upon them through the inclusion of additional excited states. Application of our model to the functionalized forms of the parent HTI molecule 1 reveals that the isomerisation mechanism is indeed altered in response to the substitution at the para-position. We find changes in the energy landscape, e. g. the elimination of a small barrier in the S 1 , that are consistent with experimentally measured increases in rate for EDGcontaining derivatives of HTI compared to the base compound. This concurs with various earlier experimental studies. [36,39] We also find that for EWG-containing derivatives, pyramidalization along the S 1 surface becomes more and more favourable, and becomes more prominent, particularly for halogenated derivatives 4 and 6. The increased involvement of the pyramidalization coordinate and the larger barriers in the S 1 that need to be overcome to reach the S 1 /S 2 crossing points are consistent with experimental findings of diminished reaction rates. [39,73] Furthermore, we show that subtle modifications to the molecule can result in large differences in the propensity of the molecules to act as molecular photoswitches; for example the boosted thermal stability shown by 4 or the increased oscillator strength shown for 2. We also demonstrate the sensitivity that absorption spectra display towards functionalization. In order to achieve red-shifted absorption peaks, EDGs capable of resonance such as in 2 are required. In 2, we also found a significant boost in oscillator strength. EWGs such as Cl and Br in 4 and 6 lead to a blue shift of the absorption maxima.
In this study we have shown that photoswitches based on hemithioindigo are sensitive to chemical substitution and properties such as thermal stability, absorption spectra and even the mechanism can all be tuned depending on the chosen functional group. We have interrogated the potential energy landscapes at high levels of theory, across the key reaction coordinates and found that pyramidalization in the excited state becomes more prevalent in derivatives including EWGs, whereas EDGs appear to promote a simpler mechanism dominated by dihedral rotation. Whilst it has been previously observed experimentally that substitution on the para-position introduces significant effects on the rate, the geometric distortions that occur have not been explicitly determined. [36] In this study we provide a mechanistic theoretical explanation for these experimental observations. As a caveat, we note that molecule substitution will likely also affect the interaction with the solvent during the photodynamics, which we have not considered in this study. A more detailed understanding of photoisomerization rates will require nonadiabatic molecular dynamics simulations, which have provided valuable findings for reactions of other photoisomerization mechanisms. [74,75] However at present, due to the size of HTI molecules in addition to the size of the active space required in the calculations, application of these methods is not straightforward. The emergence of machine learning methods in the context of excited-state dynamics simulations will support future works in this direction. [76] All molecules in this study featured substitution on the para-position, yet a whole library of HTI derivatives -both mono-and multi-substituted -has been reported and is likely to show different behavior to the molecules explored in our study. [29] HTI is an extremely promising photoswitching system and a complete characterization of how the mechanism may change upon the inclusion of substituents is imperative in order to design new and improved switching molecules. Future studies should assess the contributions of individual functional groups across multiple substitution sites in order to effectively capture the extent to which switching capabilities change in response to not only the nature of a functional group, but also the substitution position. Such a study would be extremely helpful for the design of HTI photoswitches and could assist in a targeted tuning of properties.