Unraveling the mechanisms of triplet state formation in a heavy-atom free photosensitizer

Triplet excited state generation plays a pivotal role in photosensitizers, however the reliance on transition metals and heavy atoms can limit the utility of these systems. In this study, we demonstrate that an interplay of competing quantum effects controls the high triplet quantum yield in a prototypical boron dipyrromethene-anthracene (BD-An) donor–acceptor dyad photosensitizer, which is only captured by an accurate treatment of both inner and outer sphere reorganization energies. Our ab initio-derived model provides excellent agreement with experimentally measured spectra, triplet yields and excited state kinetic data, including the triplet lifetime. We find that rapid triplet state formation occurs primarily via high-energy triplet states through both spin–orbit coupled charge transfer and El-Sayed's rule breaking intersystem crossing, rather than direct spin–orbit coupled charge transfer to the lowest lying triplet state. Our calculations also reveal that competing effects of nuclear tunneling, electronic state recrossing, and electronic polarizability dictate the rate of non-productive ground state recombination. This study sheds light on the quantum effects driving efficient triplet formation in the BD-An system, and offers a promising simulation methodology for diverse photochemical systems.


I. INTRODUCTION
][10][11] For photosensitizers to function efficiently, the electronic excitation needs to be generated in high yield and persist for a long time.One strategy to achieve this is to engineer the sensitizer to rapidly convert short-lived singlet excited states that are generated through photoexcitation into triplet excited states through intersystem crossing.Relaxation of triplet excited states to the singlet ground state is spin-forbidden, allowing the excitation to persist for orders of magnitude longer than in singlet excited states.In many photosensitizers, efficient intersystem crossing is facilitated by the presence of heavy atoms, such as transition metals, which enhance the spin-orbit coupling between singlet and triplet excited states.3][14][15][16] Understanding how triplet formation happens in these systems is essential for the design of other photocatalysts and photosensitizers.Using explicit molecular simulations of ab initio derived models, we reveal the mechanism by which triplet state formation occurs in a molecule made of only light elements.
8][19] BD -An has recently found applications in synthetic chemistry [20][21][22] and its derivatives have been investigated for phototheraputic applications. 23The competing photophysical processes and the electronic excited states involved are summarized in Fig. 1. a) Electronic mail: tom.patrick.fay@gmail.comb) Electronic mail: dlimmer@berkeley.eduBD -Anuses excited state charge transfer from an anthracenyl (An) group to the photoexcited S BD * forming an S CT state, to enable rapid triplet T BD * formation with a high experimental yield, Φ T = 0.93-0.96 18,24.Naively one might expect excited state charge transfer to reduce the triplet quantum yield, since the charge transfer state provides a charge recombination pathway for relaxation to the singlet ground state.However experiments indicate that charge recombination is suppressed by the large charge recombination free energy change, push- ing this reverse electron transfer deep into the Marcus inverted regime, where increasing the free energy change increases the activation energy. 18This effect is captured qualitatively by Marcus' theory for the reaction rate constant 25,26 where   is the coupling between electronic states  and , Δ → is the free energy change of the reaction and  is the reorganization energy, which encodes how solvent fluctuations and intramolecular vibrations control electronic state transitions, ℏ is Planck's constant and  B  is Boltzmann's constant times the temperature.Spin conserving charge recombination to the ground state is in the Marcus inverted regime, −Δ → ≫ , which requires a significant activation energy to proceed, whilst for the spin-orbit coupled charge transfer to the triplet excited state −Δ → ≈ , the reaction is approximately activation-less and thus this spin-forbidden process is competitive, despite   being much smaller for the spin-forbidden charge recombination.However, Marcus theory is not accurate in the inverted regime due to significant nuclear quantum effects, and alternate triplet formation pathways via high-energy triplet states could contribute, as has been observed in TREPR studies wherein T CT and T An * intermediates were detected at low temperatures. 24 aim to investigate the efficiency of BD -An triplet state generation in solution, going beyond the Marcus picture through first principles computational and theoretical methods, in order to explain how spin-crossover competes with charge recombination and fluorescence in solution.To this end, we interrogate each of the photophysical processes outlined in Fig. 1A by combining electronic structure calculations, molecular dynamics simulations and non-adiabatic rate theories. 26,27Our aim is to develop models that quantitatively predict experimental observables and give physical insight into mechanisms of triplet formation.[33][34] The importance of solvent effects poses a particular challenge in developing a first principle understanding of triplet state formation, because this necessitates the use of explicit solvent models and molecular dynamics. 27However common general force fields for organic molecules are only applicable to describe the ground electronic state of these systems.Previous studies have primarily used gas phase electronic structure calculations to rationalize observed behavior, 18,19 but these have not attempted to quantitatively predict rate constants from first principles.To address these challenges, we have developed a protocol for excited state force field parameterization, enabling us to accurately describe solvent fluctuations that control charge transfer processes in ground and excited states.With these tools, we show that the photophysics of BD -An FIG. 2. A) Absorption and B) emission spectra of BODIPY-Anth comparing calculated and experimental spectra with and without shifts in the excited state energies.The simulated line-shapes are obtained from the spin-boson mapping described in the main-text with bespoke force-fields for the excited states.The energy differences between excited states were obtained from DLPNO-STEOM-CCSD/def2-TZVP(-f) calculations combined with solvation energies from molecular dynamics.Experimental spectra obtained from Ref. 18   can be quantitatively predicted and mechanisms of triplet formation can be understood in detail.We start by providing a brief description of the computational methods used in this study.We then show our results for predicted spectra, free energy changes and rate constants, followed by a discussion of how these can be used to understand efficient triplet formation in BD -An.

II. STATE ENERGIES AND SPECTRA
To validate our molecular model, we have computed the BD -An absorption and fluorescence spectra (shown in Fig. 2).We calculated gas phase energies of the excited states using high-level wave-function based the DLPNO-STEOM-CCSD/def2-TZVP(-f) method 35,36 (or DLPNO-CCSD(T)/def2-TZVP(-f) for the T BD * and T An * states 37 ), with geometries for each of the excited states obtained from TDA-TDDFT 38,39 with the ωB97X-D3/def2-SVP functional 40 and basis set. 413][44] We found that wave-function based methods, which account for orbital relaxation in the excited state are required in order to obtain an accurate S 0 -S CT gap.
In the absence of solvation effects, the S BD * state is lower in energy than the S CT state by about 0.5 eV (see SI for list of energies), which is inconsistent with the fluorescence spectrum, which shows a clear peak from the CT state at lower energies than the S BD * peak.Thus in order to predict solvation effects and spectral line-shapes, we constructed bespoke force-fields for the ground and excited states of BD -An, which enabled us to perform molecular dynamics simulations to efficiently compute spectra with the spin-boson mapping. 45Geometries and Hessians from TDA-TDDFT to were used to parameterize intramolecular force-fields 46,47 based on the OPLS-AA forcefield. 48,49Electronic polarizability was accounted for using the Drude oscillator model. 50We used the same procedure to parameterize both polarizable 50 and non-polarizable force-fields for the acetonitrile (ACN) solvent, with further non-bonded parameter refinement targeting the dielectric properties of the solvent.The BD -An molecule was solvated in a box of 512 ACN molecules, and energy gap correlation functions were calculated from NVE trajectories, initial after NPT and NVT equilibration (full details are given in the SI).
From the molecular dynamics (MD) trajectories, the spinboson mapping was constructed, from which spectra were then calculated. 45,51,52In this approach the full anharmonic potential energy surfaces   are mapped onto effective harmonic potential energy surfaces.Observables of this harmonic model are fully determined by the spectral density J  ().We fit the spectral distribution   () = J  ()/ from the energy gap correlation function obtained from molecular dynamics, 29 where Δ =   −   , Δ = Δ − ⟨Δ⟩  and ⟨• • •⟩  denotes the classical phase space average over the equilibrium distribution for state  with dynamics calculated on the same surface.
For the absorption spectrum we use dynamics on  = S BD * and for the fluorescence spectra we use  = S 0 to compute the mapping, and the reorganization energy  is fit from free energy calculations using the same force fields (see below for details).From this mapping the spectra can be calculated from the Fourier transform of correlation function   (), which is given by The absorption,   (), and fluorescence,   (), spectra (with unit area) are then given by Further details of force-field development and the spin-boson mapping are provided in the SI.The unshifted spectra calculated from the spin-boson mapping using DLPNO-STEOM-CCSD/def2-TZVP(-f) gas phase energy gaps are shown in Fig. 2 as dashed lines.Our calculated spectra show good overall agreement in the spectral line shapes, without any additional fitting, capturing the narrow S BD * peak in the absorption and fluorescence spectra, including a small vibronic side band at about 1500 cm −1 from the main peak, as well as the broad S CT fluorescence band.The agreement in the vibronic structure in the S BD * peaks suggests the fitted force fields capture the reorganization energies between excited states relatively well.However we see that the unshifted absorption spectrum calculations underestimates the S BD * energy, which we attribute to the fact that the triple zeta def2-TZVP(-f) basis set is likely still not sufficient for this system.As a result, we shifted all excited states by 805 cm −1 in order to fit the experimental absorption spectrum.This simple shift is justified by the fact that all excited states shift by ∼ 0.15 eV on increasing the basis set size from def2-SVP to def2-TZVP(-f), but differences between excited state energies change by much less (see SI for details).Furthermore it has been found the EOM-CCSD has typical errors of around 0.3 eV≈2400 cm −1 for charge transfer states, so introducing a shift of 805 cm −1 seems justifiable.This shift is also used later in the free energy and rate calculations.
Using the shift from the absorption spectrum, the fluorescence spectrum (Fig. 2B) was calculated as a weighted sum of the S CT and S BD * emission spectra, with weights given by the transition dipole moments from DLPNO-STEOM-CCSD,  2 S BD * ,S 0 = 7.59 a.u. and  2 S CT,S 0 = 0.54 a.u., and equilibrium populations of the two states given by the free energy change of charge separation Δ CS , i.e.
The assumption of equilibrium between the S BD * and S CT states is justified by the fact the time-scale of equilibration of these states is ∼ 10 3 times shorter than the lifetime of these states (as we will discuss shortly).We have also computed the fluorescence spectrum assuming the populations of the S BD * and S CT states are given by the experimental estimate, Δ CS,exp , based on the approximate Weller equation, which is about 0.2 eV larger than our estimate. 18Because Δ CS,exp > 0, the S CT state is significantly less populated relative to the S BD * state and the S CT fluorescence peak is almost completely suppressed, which does not agree with the experimental spectrum.This suggests that the Weller equation cannot be used reliably when free energy changes are close to zero.As an interesting aside, the strongest S CT − S  coupling (see table I) is to the S 0 state, by over a factor of 10, which indicates that the intensity borrowing effect responsible for the S CT emission arises primarily from mixing between S CT and S 0 states, rather than S CT and S BD * states, as has previous been assumed. 24

A. Thermodynamics
Charge separation, the S BD * → S CT process, and charge recombination, the S CT → S 0 process, both play an important role in efficient triplet formation.Efficient charge separation is required to suppress fluorescence from the S BD * state, but slow charge recombination is needed to enable intersystem crossing to occur to generate triplet states.From our excited state force-fields, we have calculated free energy changes associated with these processes from molecular dynamics and the multi-state Bennett acceptance ratio (MBAR). 53,54As discussed above, the calculated free energy change for charge separation is −0.057 eV, thus population of the S BD * state is reduced and fluorescence is suppressed.
We have also calculated the rates of these processes from the same MD simulations, by calculating the probability of two states being at resonance.This probability controls the classical Fermi's Golden rule (FGR) rate for the transition between  and . 55The free energy along the energy gap coordinate, Δ =   −   , is related to the energy gap distribution   () = ⟨(Δ − )⟩  by 27   () = − B  ln(   ()) + (   −   ) for  =  or .In Fig. 3 we show the free energy profiles calculated from MD simulations on each of excited state surfaces with the polarizable ACN model using MBAR.The crossing point of the two curves gives the free energy barrier for the transition, which dictates the classical FGR rate, 27 If the free energy curve is perfectly quadratic, then this reduces exactly to Marcus theory [Eq.(1)]. 25,26,29For the charge recombination the crossing point occurs outside of the sampled region, so we extrapolated to the crossing point using a quadratic polynomial ansatz for the free energy, fitted to the cumulative distribution function.This procedure was found to result in very little loss in accuracy when compared to umbrella sampling/weighted histogram analysis 56 calculations performed using the nonpolarizable ACN model (see SI for details).
The free energy curves for charge separation and charge recombination are shown in Fig. 3A and B, where we see charge separation lies in the Marcus normal regime, whereas charge recombination is deep in the Marcus inverted regime, with a much larger free energy barrier.Using diabatic state couplings calculated from the generalized Mulliken-Hush method 57 with DLPNO-STEOM-CCSD calculations, we can directly calculate the classical FGR rates for these processes (couplings |  | are shown in table I).The classical FGR charge separation rate is 4.8 × 10 10 s −1 , about a factor of 10 smaller than the experimentally observed rate of 5.4 × 10 11 s −1 , however the charge recombination rate is predicted to be 1.1 × 10 −17 s −1 , which is more than 10 24 times too small compared to the experimental estimate of 2.3×10 7 s −1 . 18This enormous discrepancy can be attributed to nuclear quantum effects, in particular the important role of nuclear tunneling in the inverted regime.

B. Quantum effects on rates
In order to include nuclear quantum effects in the rate calculations, we employed the same spin-boson mapping approach as was used to compute the spectra.The full FGR rate constant is given by which can be evaluated directly using Eq. ( 3).The reorganization energy  is fitted by requiring that the classical limit of the spin-boson mapping reproduces the exact classical limit rate constant, obtained from the classical free energy barrier calculated from MBAR. 54This approach to calculating the rate can be regarded as a generalization of the commonly used Marcus-Levich-Jortner theory, accounting for the full frequency dependence of the reorganization energy, which is encapsulated in   ().The final rate constant is obtained as a simple average over the rate constants calculated with spectral distributions   () and   ().
The calculated spectral distributions   () can be decomposed into inner sphere, outer sphere and cross-correlated contributions, by decomposing the energy gap into molecular and the remaining environment contributions Δ = Δ mol + Δ env .We find that the cross-correlated contribution is generally negligible for all processes in BD -An, so the reorganization energy is well-described by a simple sum of inner and outer sphere contributions.The inner and outer sphere spectral distributions are calculated with the non-polarizable ACN/solute model, with the outer sphere contribution scaled down to match the polarizable model outer sphere contributions.As can be seen in Fig. 4A, the low frequency proportion of the spectral distribution for the S CT → S 0 transition is dominated by the outer sphere contribution arising from solvent molecule fluctuations, making up ∼ 50% of the reorganization energy, which is well approximated by the Debye model. 32In contrast, the high frequency region of the spectral density is dominated by the inner sphere contribution from changes in equilibrium bond lengths in the BD -An molecule on charge transfer.The inner sphere spectral distribution has contributions over a range of frequencies from around 500 to 1600 cm −1 , all of which contribute to tunneling enhancement of the S CT → S 0 rate, although the dominant mode at ∼ 1400 cm −1 likely corresponds to a C --C stretching motion within the aromatic rings.Qualitatively similar spectral distributions were found for the other charge transfer processes.For processes which do not involve charge transfer the spectral distribution is dominated by the inner sphere contribution, as can be seen for the T An * → T BD * process in Fig. 4B.
When accounting for nuclear quantum effects, the S BD * → S CT rate goes up by a factor of ∼ 3 to 1.46 × 10 11 s −1 , and the Points correspond to free energy curves calculated with MBAR and lines correspond to polynomials fitted to the MBAR cumulative distribution functions (see SI for details).G) A snapshot for molecular dynamics simulations on the S 0 potential energy surface.H) A scheme highlighting the processes in A-F.( = 0) with   ( = 0) (see SI for details). Couplings averaged over gas-phase equilibrium geometries of  and , 18, estimated from spectroscopic measurements. Rate constants from spin boson mapping. Linear response value:  = (⟨Δ⟩  − ⟨Δ⟩  )/2. Radiative rate constant (Eq.( 9)).ℎ with recrossing correction and  without recrossing correction. estimated from spectroscopic measurements. 18 Using reorganization energy from nonpolarizable umbrella sampling calculations (see SI).
S CT → S 0 rate goes up by over 10 24 to 1.0 × 10 8 s −1 , and both calculated rates are now much closer to the experimentally measured values, agreeing much better with the experimental value.Application of Marcus-Levich-Jortner theory with the same inner and outer sphere reorganization energies and a characteristic inner-sphere frequency of 1500 cm −1 also predicts about a 10 24 -fold increase in the rate constant, compared to Marcus theory.This suggests that the large increase is robust to the details of the spectral density.Electronic polarizability is essential to account for in calculating the charge recombination rate.When a non-polarizable model is used instead, the free energy change of the reaction is effectively unchanged but the reorganization energy goes up by nearly 0.1 eV.This lowers the activation energy and accelerates the rate by around a factor of three.
Care should however be taken when using FGR to calculated the charge recombination rate.This is because the diabatic coupling for charge recombination process,   = 1904 cm −1 , is about 20 times larger than  B , and thus higherorder diabatic coupling effects beyond FGR, may be important (although large nuclear quantum effects in the FGR rate have been observed to reduce the importance of higher order effects). 34The large difference in couplings arises from the BD π orbitals involved in the transitions.The S BD * → S CT coupling involves an interaction between π An and π BD (Fig. 5A) orbitals, whereas S CT → S 0 coupling involves the π An and π * BD (Fig. 5B) orbitals.As can be seen in Fig. 5 the π BD has minimal density on the carbon atom bonded to the An, group, whereas the π * BD orbital does.In order to investigate the potential role of higher-order diabatic coupling effects in the S CT→S 0 tran-sition, we have performed Hierarchical Equations of Motion (HEOM) calculations a simple model for this transition.The spectral density for the transition is coarse-grained down to a low-frequency outer-sphere portion described with a Debye spectral density and the inner sphere portion is described with a single under-damped Brownian oscillator spectral density, with a characteristic frequency of 1400 cm −1 .The coarsegrained spectral density is shown in Fig . A. For this simplified model the exact open quantum system dynamics can be obtained using the HEOM method, and from this the rate constant as a function of   can be obtained.These rates are shown in Fig. B. We see that the rate constant is still fortuitously very well described by Fermi's Golden rule for this model, with only a factor of ∼ 0.9 reduction in the rate constant at the calculated value of   .We include this as a correction to the Fermi's Golden rule  S CT→S 0 that we calculated with the full spectral density.
Radiative recombination from the S CT state can also occur in BD -An, either through thermally activated delayed fluorescence via the S BD * state, or directly.The radiative rates can be calculated from the fluorescence spectra obtained from the spin-boson mapping as 26,60 where  () is the fluorescence line-shape computed from the spin-boson mapping.From this we find the fluorescence rate from the S BD * state to be 1.1 × 10 8 s −1 and the fluorescence  rate from the S CT state to be 7.7 × 10 8 s −1 .Assuming a preequilibrium between the S BD * and S CT states, as is justified by the large charge separation rate constant, we find that only 63% of the S 0 re-formation occurs by direct non-radiative recombination, with 15% of recombination events happening by radiative S CT recombination and 22% occurring via S BD * thermally activated delayed fluorescence.

IV. TRIPLET STATE FORMATION AND LIFETIME
As with the charge separation and charge recombination processes, we have calculated the free energy changes and free energy profiles for the three triplet formation pathways: from the S CT state to the T CT, T An * and T BD * states (Fig. 3D-F).Free energy calculations reveal that the three pathways are thermally accessible, with all three states lying lower in energy than the S CT state.Furthermore all three pathways are approximately activation-less, which is at first surprising given that each process has a very different free energy change.The T BD * pathways has a larger |Δ → |, than the T An * pathway, but the T An * pathway has smaller inner and outer sphere reorganization energies, so this pathway is also approximately activation-less.The T CT pathway has a very small reorganization energy which is dominated (∼ 90%) by the inner sphere contribution, at only 0.11 eV.This is because the S CT and T CT states have the same orbitals occupied, so the reorganization energy is dictated only by differences in the exchange energy which alters bonds lengths.However the net exchange energy difference between these states is small because the unpaired electrons have low spatial overlap, so overall the reorganization energy is low and this transition is approximately activationless.Much like the spin-conserving charge separation and charge recombination, about 50% of the reorganization energies for the charge transfer processes is outer sphere, with the remaining 50% arising from inner sphere reorganization, although there is a significant range of reorganization energies for the charge transfer processes, from 0.48 eV to 0.58 eV.In contrast, the reorganization energies of processes which do not involve charge transfer are dominated by the inner sphere contribution, 89% for the S CT → T CT spin-crossover and 99% for the T An * → T BD * triplet-triplet energy transfer, as illustrated in Fig. 4B.The triplet-triplet energy transfer still has a reorganization energy comparable to the charge transfer processes, at 0.57 eV, due to a large change in bond order in both the BD and An units in this process.Further analysis of the inner/outer sphere reorganization energies are given in the Supporting Information together with all calcualted spectral densities.
We have also calculated the SOC couplings between the different S CT and triplet states using TDDFT (ωB97X-D3/def2-TZVPP/CPCM(ACN)) and the spin-orbit mean-field (SOMF) treatment of spin-orbit coupling. 61,62The two spin-orbit coupled charge transfer (SOCT) pathways have the largest SOC couplings, at 0.79 cm −1 and 0.63 cm −1 for the T BD * and T An * whilst the formally El-Sayed's rule forbidden pathway has a smaller coupling at 0.21 cm −1 .Using these couplings and the spin-boson mapping, we find that two El-Sayed's rule allowed transitions, via T An * and T BD * , occur at very similar rates, with S CT → T BD * occurring only about 20% faster than the S CT → T An * formation.The triplet-triplet T An * → T BD * energy transfer is also activation-less (Fig. 3C), and has a coupling from fragment energy/charge density (FED/FCD) FIG. 6. A) The coarse-grained model spectral distribution for the S CT→S 0 transition, consisting of a low frequency Debye contribution  D () = (1/2)/(1 + (/  ) 2 ), with  D = 0.1831, and an under-damped Brownian oscillator contribution  BO () = (1/2)Ω 2 /(( 2 − Ω 2 ) 2 +  2  2 ) with  = 4 and Ω = 6.76.The reorganization energy for the Brownian oscillator portion is  = 8.6780 and for the Debye portion is  = 10.1459.B) The rate constant from HEOM calculations for the coarse-grained spectral density as a function of   together with the FGR predictions.The value of   for the S CT→S 0 transition is also indicated.Calculations were performed using the heom-lab code 58 using the HEOM truncation scheme from Ref. 59. calculations 47,63,64 of 2.57 cm −1 , and so occurs about 10 times faster than the triplet formation rate, accelerated by a factor of 1.6 by nuclear quantum effects, so the steady state population of T An * would be difficult to observe spectroscopically at room temperature.The El-Sayed's rule forbidden transition to the T CT state also contributes to triplet formation, although it occurs about 4.5 times slower than T BD * formation.The T CT state very rapidly recombines to the T An * or T BD * states, with these spin allowed transitions occurring at least ∼ 10 4 times faster than the corresponding spin-forbidden transitions, so the T CT state would be very difficult to observe directly at room temperature.Overall the T CT T An * , and T BD * pathways contribute 14%, 47%, and 39% respectively to the overall triplet formation.Surprisingly the most significant pathway is the T An * pathway and not the direct T BD * pathway, which can be rationalized by the lower activation barrier for the T An * spinorbit coupled charge recombination.The observation is consistent with TREPR experiments in which all three triplet states were observed, although at much lower temperatures (80K) in a very different medium (a dichloromethane/isopropanol solid matrix).This work shows that multiple triplet formation pathways, including those forbidden by El-Sayed's rule, can contribute at room temperature in polar solvents.The presence of multiple triplet recombination pathways may also explain the large spread of effective spin-orbit coupled charge transfer rates observed in the family of BD-Aryl molecules studied in Ref. 18.
Using all of the computed rates, we have estimated the observed charge separation and charge recombination rates, as well as the triplet yield.The effective charge separation rate corresponds to the observed equilibration rate between S BD * and S CT states i.e.  CS,eff =  S BD * → S CT +  S CT→ S BD * .Likewise the effective charge recombination rate corresponds to the observed decay rate of the S CT state, which under a pre-equilirbium approximation for the S BD * ⇌ S CT interconversion is given by  CR,eff = S CT ( CR +  F, S CT→S 0 ) + S BD *  F, S BD * →S 0 (10)  where S CT = 1 − S BD * =  CS /(1 +  CS ), with  CS =  −Δ S BD * → S CT / B  ,  F, S BD * is the calculated fluorescence rate from the S BD * state back to the S 0 state and  CR is the total recombination rate from the S CT state, i.e.
CR =  S CT→S 0 +  S CT→ T CT +  S CT→ T An * +  S CT→ T BD * .(11)  The triplet quantum yield Φ T is calculated as Φ T = S CT ( S CT→ T CT +  S CT→ T An * +  S CT→ T BD * ) CR , with  CR = 1/ CR,eff , and the fluorescence yield Φ F is Φ F = S BD *  F, S BD *  CR .We also computed the fraction of nonradiative transitions which produce a triplet state,  CRT = Φ T /(1 − Φ F ), as measured in Ref. 18.
The calculated and experimental values of the rates and yields are summarized in table II.Overall we see excellent agreement between the calculated rates/yields and the experimental measurements from Refs.18 and 24, with less than a factor of 4 error in the charge separation rate and only a factor of ∼ 1.6 error in the charge recombination rate.Similar we only slightly underestimate the triplet yield, with our calculations yielding 0.  between 0.93 and 0.96.If we only included the dominant S CT → T BD * triplet formation pathway, the triplet quantum yield would only be ∼ 0.6, and the error in the rate would be over a factor of 3. We also find that suppression of the charge recombination also plays a large role in efficient triplet formation, which is facilitated by polarizability and recrossing effects.Without including electronic polarizability, the charge recombination rate would be enhanced to ∼ 1.0 × 10 8 s −1 , which would reduce the triplet quantum yield to ∼ 0.63.This corroborates the conclusions drawn in Ref. 18, although we find that multiple triplet pathways also enable the triplet formation to compete with charge recombination, which is suppressed by several effects.The net fluorescence quantum yield from S BD * that we calculate, 0.045, is also in good agreement with the experimental values, between 0.01 and 0.018.These results suggest that the intersystem crossing rates are being slightly underestimated by our models, possibly due to errors in the reorganization energies or the spin-orbit couplings obtained from TD-DFT, which are all less than 1 cm −1 .
The triplet lifetime  T = 1/ T BD * →S 0 plays an important role in determining the utility of a triplet sensitizer or photocatalyst, with longer-lived triplet states allowing more time for diffusive encounters with other molecules enabling more efficient energy transfer.We have also calculated the triplet lifetime for BD -Anusing the methods described above, and we also find good agreement between our calculated value for  T and experimental measurements (table II), with an error of only ∼ 20%.From a simulation perspective, this requires an accurate calculation of the free-energy barrier, which requires enhanced sampling since the transition is very deep in the Marcus inverted regime, since it displays a highly non-quadratic free energy curve.This was achieved using the non-polarizable model with umbrella sampling 65 on the energy gap coordinate Δ sampled with the Fast-Forward Langevin algorithm. 66Use of the non-polarizable model is justified because over 99% of the reorganization energy is inner sphere for both ACN models, and solvent polarizability has less than a 1 meV effect on the free energy of the T BD * state.As with the spin-conserving charge recombination, because the transition is deep in the inverted regime and the spectral distribution is dominated by high frequency inner sphere contributions, there is a very large nuclear quantum effect of over 10 7 in the rate constant.One significant source of uncertainty in this is the validity of the spin-boson mapping, where rates calculated from the spectral distribution obtained from T BD * and S 0 dynamics vary by about 50%.This means that methods that more rigorously account for asymmetry and anharmonicity in the potential energy surfaces, while also accounting for nuclear quantum effects, may be needed to more accurately compute triplet lifetimes for this system and other related systems. 55,67,68However given the simplicity of the spin-boson mapping and its accuracy in this case, it is clearly still useful in prediction of non-adiabatic rates.

V. CONCLUDING REMARKS
Through this study, we have found that triplet formation in the photosensitizer BD -An hinges on a subtle balance of effects.Firstly charge separation occurs efficiently, which suppresses radiative decay from the S BD * state.Secondly multiple triplet recombination pathways can operate, due to the range of reorganization energies and free energy changes associated with the rate-limiting intersystem crossing steps in each pathway, and in fact the high-lying triplet pathways make-up the major contribution to triplet formation, rather than direct SOCT to the ground triplet state.Thirdly, spinconserving charge recombination to the S 0 state is slowed down a high free energy barrier, with the transition being deep in the inverted regime, as well as diabatic recrossing effects, a significant portion of which arises due to electronic polarizability.The S CT state energy plays an important role in triplet formation, since an increase in energy would increase fluorescence from S BD * , but a decrease in its energy would reduce the barrier for spin-conserving charge recombination because this process is in the Marcus inverted regime.Capturing all of these effects depends on a complete description of the photophysics including accurate calculations of electronic state couplings, explicit solvent fluctuations, polarizability, to capture outer sphere reorganization energies, as well as an accurate description of molecular potential energy surfaces and inner sphere contributions to reorganization energies, as well as the nuclear quantum effects arising due to high frequency vibrations, which accelerate some processes by many orders of magnitude.Enhanced sampling techniques are also necessary to obtain accurate free energy barriers for important processes, namely the triplet decay.
The simulation techniques and bespoke force-field parametrization approach developed here paves the way for a quantitative modeling of other triplet photosensitizers and related systems, 69 possibly even enabling straightforward computational screening for properties such as the triplet lifetime.Comparison between simulated and experimental optical spectra indicates that a major source of error is in gas phase energies of excited states.Even the popular wavefunction-based DLPNO-STEOM-CCSD method appears to significantly underestimate transition energies, although the ground-state DLPNO-CCSD(T) method which can be used to calculate the T 1 -S 0 gap seems robust.We also note that whilst the approximate spin-boson mapping seems fairly reliable for these systems, its application to deep inverted regime processes requires scrutiny.1][72][73][74][75][76][77][78][79] The S CT → S 0 transition poses a particular challenge, since it is deep in the inverted regime, nuclear quantum effects are very large and strong diabatic coupling means there may be some effects missed by FGR, which we have estimated using open quantum dynamics simulations.Furthermore in this study we have neglected non-Condon effects 80 and potential spin-vibronic effects, 81 which could also play a role in determine the rates of conversion between excited states in this system.Future investigations into these potential effects could provide further insight into triplet formation in BD -An.
Overall, we believe the mechanistic insights gained from this study, which would be difficult to probe directly with experiment alone, could help light the path towards the development of novel and interesting photochemistry in related systems.The observation that high-energy triplet pathways dominate at room temperature opens the door to the intriguing possibility of engineering triplet anti-Kasha's rule systems, 82 in which higher energy triplet states could be used to drive photochemistry.This could be particularly promising since triplet-triplet energy transfer is strongly distance dependent, 83 so spatial separation of chromophore units could be used to extend the lifetime of high-lying triplet states.In summary, our comprehensive study highlights the intricate balance of factors influencing triplet formation, including the significance of charge separation efficiency, multiple recombination pathways, and nuclear quantum effects.Moving forward, this mechanistic understanding could steer the development of novel photochemical systems, with a wide range of potential applications.For triplet-triplet couplings the mixed fragment-charge density/fragment excitation density method was used to compute diabatic states and spin-orbit couplings. 1Löwdin charges and transition charge were computed using TDA-TDDFT B97X-D3/def2-TZVP at each of the excited state geometries with CPCM treatment of the ACN solvent.From these the fragment charge and excitation density operators were constructed for the adiabatic states corresponding to the states of interest, which were then diagonalized simultaneously using a Jacobi sweep algorithm. 2Spin-orbit couplings were obtained using the TDDFT B97X-D3/def2-TZVPP/CPCM(ACN) using the SOMF method in Orca 5.0.3. 3,4

B. Force field parameterization
A classical molecular mechanics description of the inter and intramolecular forces was employed for the different electronic states of BODIPY-Anth.For the intermolecular part, a non-polarizable fixed charge model was used for the electrostatic component together with a Lennard-Jones description of repulsion and dispersion effects.Charges for each state were parameterized by taking an average of CHELPG charges computed for each state at the SOS-B2G-GPLY/def2-TZVP(-f) level of theory, with and without CPCM treatment of the acetonitrile solvent as implemented in ORCA. 5,6Lennard-Jones parameters were assumed to be the same for each electronic state and are taken from the OPLS-AA force field with boron parameters taken from Ref. 7.
The intramolecular forcefield is given as a sum of bonds, angles, torsion dihedrals and improper dihedrals, and bond-bond and bond-angle correlation terms, electrostatics and Lennard-Jones terms, The long-range forces are ignored for 1-2 and 1-3 bonded atoms and scaled by 0.5 for 1-4 bonded atoms.Periodic dihedral and improper terms are used, as is standard for the parent OPLS-AA forcefield.Harmonic angle forces are used and bi-linear bond-bond and bond-angle terms were used for all 1-3 bonded triples of atoms, The  ,0 and  ,0 terms were set to the QM equilibrium geometry parameters.A Morse potential truncated at fourth order was used to describe bond stretches, with the Morse  parameter for each bond parameterized from Hessians at displaced geometries in the ground electronic state.These potentials were parameterized by first fitting force-constants and equilibrium bond lengths/angles to the gas phase QM hessian calculated at the B97X-D3/def2-SVP equilibrium geometries, then refined bond lengths and bond angles to reproduce the gas phase equilibrium structures.The first step in parameterization after LJ parameters and charges were assigned was to fit intramolecular force field terms were by minimizing where H  denotes the partial hessian block.For this step the bonds are treated as harmonic.First the equilibrium bond lengths and angles are allowed to vary arbitrarily but in a second step they are fixed, with any lengths more than 0.0001 nm and angles more than 2.5 • away from the QM equilibrium geometry values fixed at these boundaries.The equilibrium bond lengths and angles are then refined to minimize The improper and dihedral parameters were then refined again minimizing L Hess , and the geometry refinement was repeated after this.The forcefield obtained with the lowest value of L geom was taken.The hessian fitting with L Hess with multiple randomly displaced geometries was used to parameterize a re-scaling of the bond force constant   and the Morse   parameter for the ground-state.The same re-scalings and   parameters were then used for the excited states.Symmetry was used to constrain values of the force field parameters for equivalent atoms at each stage of the fitting.
For the polarizable solute model, Drude oscillators are placed on all non-hydrogen atoms.The atom polarizabilities and Thole parameters are fitted by fitting the molecular polarizability for the groundstate for each equilibroum state geometry from SOS-B2GP-PLYP/def2-TZVP(-f) calcualtions.We found that the poalrizability does not change considerably between charged states and singlet/triplet states of the individual components, so we consider this a reasonable approximation.Each atomic polarizability and Thole parameter was regularized vs values fitted to small symmetric molecules, namely C 6 H 6 for aromatic atoms, C 2 H 6 for the methyl groups, CF 4 for F and BF 4 -for B. In the fitting the CHELPG charges are simultaneously readjusted to optimally reproduce gas phase dipole and quadrupole moments, as well as the molecular polarizability tensor.In the fitting the charges are also regularised against the gas phase CHELPG values to avoid overfitting.This procedure cannot account for how charge may flow around the molecule when it is placed in a dielectric environment.To account for this, we also added the difference in density-derived Hirshfeld charges between vacuum and acetontrile CPCM calculations at the SOS-B2GP-PLYP/def2-TZVP(-f) level of theory to the final re-scaled charges obtained from the polarizability fitting.

C. Solvent models
Two acetonitrile force field models were parameterized to describe this system: a non-polarizable and a polarizable model.Both models were constructed by first performing hessian fitting to the RI-MP2/cc-pVQZ geometry and to obtain harmonic bond-stretch and angle-bend force constants.For the non-polarizable model Lennard-Jones parameters were taken from the OPLS-AA force field and charges were taken as a mixture of CPCM and gas phase CHELPG charges,   =  CPCM  + (1 − ) gas  , with the mixing parameter  optimized to reproduce the experimental density at 298 K.
For the polarizable model a single Drude oscillator was added to the nitrile group carbon atom, with an anisotropic polarizability obtained from a RI-MP2/cc-pVQZ calculation.A modified version of the three-step parameterization procedure was then followed to parameterize charges, Lennard-Jones  parameters and  parameters. 8First the charges were re-scaled according to gas  in order to reproduce the experimental dielectric constant at 298 K. Second, the  parameters in the Lennard-Jones potential were re-scaled to reproduce the vaporization enthalpy at 298 K, and third  parameters were re-scaled to reproduce the experimental density at 298 K.The resulting force-fields both accurately reproduce these properties of ACN at 298 K, as well as the dielectric relaxation time.
All fitted free energy curves shown in the main text are fitted by fitting a polynomial approximation to the free energy curves for , Ã, () =  =0     , and , Ã, () = Ã, () +  − Δ → , to to cumulative distribution function, for  =  and .The constraint the the integral of p () = 1 is also added in the fitting.For the fits to the umbrella-sampled free energy curves we only fit to the  curve.

D. Free energy calculations
The free energy change for the transformation from electronic state  to  was calculated for the non-polarizable model using thermodynamic integration runs, where the potential is given by   () = (1 − )  () +   ().Five MD simulations were performed for each transformation with  values equally spaced between 0 and 1.The Fast-Forward Langevin integrator was used to sample configurations with a time-step of 1 fs. 9MBAR was used to compute the Helmholtz free energy change Δ → from these thermodynamic integration runs.WHAM was used with these runs to compute free energy curves and free energy barriers for each electron transfer process.For the S CT → S 0 and T BD * → S 0 processes the crossing region for the two diabatic states was not sampled in the thermodynamic integration runs so umbrella sampling was performed with the the diabatic energy gap Δ =   () −   () used as a biasing coordinate, and again WHAM was used to obtain free energy curves and the free energy barrier for this electron transfer.All fast-forward Langevin dynamics simulations were performed using an in-house code using the OpenMM 7 C++ library for force and energy evaluations. 10he box containing BD -An was set-up with 512 ACN molecules, which was found to provide adequately almost completely decayed dipole-dipole correlations, using Packmol. 11Three independent NPT trajectories were run with the S 0 BD -An force field and the non-polarizable ACN model to find an equilibrated box side length of  = 3.570 nm.In the thermodynamic integration five equally spaced windows were sampled, with each window sampled for 1 ns after 0.1 ns of equilibration with a FF-Langevin friction constant of  = 4 ps −1 and a time-step of 1 fs.For umbrella sampling the same simulation parameters and trajectory lengths were used, with an umbrella sampling force constant of   = 0.01 kJ −1 mol with windows centered at integer multiples of 50 kJ mol −1 .
For the polarizable model free energy changes were obtained from MBAR with configurations sampled on each electronic state.Free energy barriers were obtained by fitting free energy curves obtained from MBAR to a polynomial.For the charge transfer steps a quadratic polynomial was used but for other processes anharmonicity plays a larger role and a higher order order polynomial was used.Error bars were obtained by bootstrapping the six independent runs.All simulations were performed using the python API of OpenMM 7 and were analyzed with in-house scripts.For the NVT runs, trajectories were equilibrated for 5 ps and sampled for 2 ns using the extended Lagrangian approach, with a time-step of 1 fs, a friction coefficient of  = 2 ps −1 for the centers of mass and 20 ps −1 for the Drude oscillators, with the Drude particle temperature set to 1 K.

E. Spin-boson mapping
For each non-adiabatic process a spin-boson mapping was constructed where an effective potential for each state is constructed as The reorganization energy  was fitted to reproduce the free energy barrier to reach the crossing point for  → , as obtained from the free energy calculations.The rate constants and spectra for the spin-boson mapping depend only on the spectral density This can be parameterized directly from the Fourier transform of the energy gap correlation function obtained from classical MD.In what follows, we parameterize this as where the spectral distribution () is calculated from MD as where Δ =   −   , Δ = Δ − ⟨Δ⟩  and ⟨• • •⟩  denotes the classical phase space average over the equilibrium distribution for state  with dynamics calculated on the same surface.By fitting the reorganization energy based off of the free energy barrier we ensure the spin-boson mapping we construct becomes exact in the high-temperature limit.There are two possible () distributions for each transformation, obtained from dynamics on either  or , and in general these distributions will differ.The similarity of the distributions, and in particular the quantities of interest computed from the two different distributions, gives an indication of how accurate the spin-boson mapping can be expected to be.For all computed rate constants, we use an average of rate constants computed with   () and   (), although we have found that first averaging the spectral distributions and computing rates from this averaged distribution yields almost identical final rate constants.The correlation function appearing in the full FGR rate constant can be evaluated as where Δ  = −Δ → .This is evaluated by discretizing the spectral density using the method in Ref. 12.The spectra can also be directly obtained from the Fourier transform of this function.The correlation function is multiplied by an exponential decay  −/ when computing the spectra to account for instrument broadening.A time-constant  = 350 fs is used, as has been used in other studies.Spectral densities for the non-polarizable model were obtained from 20 NVE trajectories each 30 ps long with a time step of 0.5 fs.The cross correlated portion of the spectral distribution was found to have a negligible effect of the calculated rate constants (as long as the total reorganization energy was held fixed at the fitted value).Furthermore, given that the Debye relaxation time for the non-polarizable and polarizable models are almost identical (see table S1) and that this model very accurately describes the outer sphere contribution, we constructed the polarizable model spectral distribution as (S14) Fig. S1.Full polarizable ACN spectral distribution and the approximate form given by Eq. (S14) for the S CT → S 0 process, calculated with dynamics on the S CT surface.
This avoids the need to run expensive self-consistent field Drude oscillator calculations to obtain the spectral density for the non-polarizable model.The proportion of inner sphere reorganization energy was set as which was found to agree well with values obtained from the co-variance based estimates.In Fig. S1 we show the approximated spectral distribution based off of the re-scaled non-polarizable ACN calculations to the full spectral distribution calculated from three NVE trajectories with the full polarizable ACN potential with SCF Drude oscillator integration.We see excellent agreement (within the uncertainty of the polarizable model simulation) between the approximate and full spectral distribution.The coupling   for the rate constant calculation was taken as the root-mean-square of the couplings between the two gas-phase equilibrium geometries of  and , i.e.   = √︁ (| ,  | 2 + | , | 2 )/2.This ensures a symmetric definition of the rate constants, but does not account for non-Condon effects.

F. Hierarchical equations of motion calculations
The hierarchical equations of motion (HEOM) method, as implemented in the heom-lab Matlab code 13,14 was used to obtain exact rate constants for the model spectral density for the S CT→S 0 transition.The system was initialised in the state ρ(0) = | ⟩⟨| ρb and the long time dynamics of the  state population were fitted to a function in order to obtain the rate constant.The calculations used the HEOM truncation scheme from Ref. 15 with  = 80, Γ = 6 and a Matsurbara expansion was used for the bath correlation functions with the  = 1 Matsubara mode included explicitly in the calculations.The remaining Matsubara terms were treated with the low temperature correction scheme in Ref. 14. Dynamics were run out to  = 10ℏ.

II. INTENSITY BORROWING
The effective transition dipole moments used in the spectral calculation are taken directly from the adiabatic DLPNO-STEOM-CCSD calculations.From the diabatization we can analyze the origin of the relatively large transition dipole moment for the S CT state.The adiabatic states are given by first order perturbation theory in the diabatic couplings   as From this we find the mixing coefficient for the S CT and S 0 states is 0.12 compared to 0.053 for the S CT and S BD * states.The transition dipole moment is given to lowest order in perturbation theory by and for the S CT → S 0 transition this is given approximately by The second term in this expression, arising from mixing of the S CT and S 0 states dominates.Thus it is primarily S CT-S 0 mixing that enables intensity borrowing in the S CT emission.

A. Solvent model properties
In table S1 we show a summary of solvent properties calculated from molecular dynamics for the solvent models.All simulations used a box of 512 ACN molecules.All simulations were performed at 298 K.These properties were considered as good target properties for accurate simulations of electron transfer in ACN at 298 K.

B. Gas phase energies
A comparison of gas phase minimum energy geometry energies, calculated with different basis sets, relative to the S 0 state energy is shown in table S2.

C. Couplings
We have computed both the spin-conserving diabatic couplings as well as spin-orbit couplings between singlet and triplet states.Table S3 shows the computed values of these   at the minimum energy geometries of state  and .The spin-orbit couplings are of course significantly smaller than the spin-conserving diabatic couplings.We note that couplings between the S CT and S 0 state are significantly larger than the S CT-S BD * couplings, and the SOC couplings are all small, with the coupling for the El-Sayed's rule forbidden S CT→ T CTtransition not being significantly smaller than the other SOC mediated charge transfer processes, despite there not being a change in orbital occupancy between these states.This can likely be attributed to the fact that although the SOC-mediated CT processes involve a change in orbital occupancy, the orbitals involved are localised on different fragments of the molecule, so the poor orbital overlap reduces the SOC between these states.Overall we see that all three S CT→T n pathways, which are all feasible according to the free energy changes, are weakly allowed by the SOC interaction.

D. Umbrella sampling
For the S CT→S 0 and T BD * →S 0 processes the transitions lie very deep in the Marcus inverted regime, which means that enhanced sampling is required to obtain accurate free energy profiles.We have performed Umbrella sampling directly on the energy gap coordinate Δ for these two transitions, in order to sample rare fluctuations of the solvent and molecule in which Δ = 0 and subsequently compute the free energy barrier.Efficient molecular dynamics sampling on the biased potentials with the polarizable ACN model is not possible, because the extended Lagrangian approach cannot be applied, so we only use the non-polarizable ACN model in these simulations.In Fig. S2 we show these free energy curves, together with a quadratic fit to the  state curve data (obtained by fitting the corresponding cumulative distribution function, see SI for details), excluding the portion of the curve where the energy gap, , is less than 1 eV.We see that in the S CT→S 0 (Fig. S2A) the quadratic fit provides an excellent fit, which extrapolates accurately to the crossing point.The effective fitted Marcus reorganization energy,  fit , obtained by equating the Marcus-theory (Gaussian) probability density at the crossing point and the Umbrella sampling probability density, agree very well, with  fit = 0.590 eV compared to the quadratic estimate  quadratic = 0.585 eV.This suggests that obtaining  fit for the polarizable model by fitting the free energy curves (which do not sample the  = 0 region) to a quadratic and extrapolating should introduce an error of only ≃ 0.005 eV.In stark contrast, for the T BD * and S 0 states there is a large deviation between a quadratic fit and the computed free energy curves for the T BD * →S 0 process (Fig. S2B).This can be attributed to the fact that ∼ 98% of the reorganization energy for the T BD * →S 0 arises from intra-molecular inner sphere contributions (see below for how this is calculated), for which asymmetry in the vibrational frequencies and anharmonicity in the molecular potential energy surfaces plays a significant role.However for the S CT→S 0 process the inner sphere contribution is only about 40% of the total reorganization energy, so the outer sphere contribution arising from fluctuations in the solvent polarization dominates, which is very well described by a Gaussian field theory.

E. Reorganization energies
It is interesting to quantify how different the reorganization energies are between the different state transitions, but as observed above, not all of the free energy curves are perfectly quadratic.However, as stated above, we can fit   ( = 0) to the Gaussian functional form in Marcus theory  G  () = (4 B ) −1/2  − (Δ → −  +) 2 /4 B  , to obtain a effective reorganization energy  fit .Predictably it is smallest for the S CT→ T CT transition, which involves no change in orbital occupancy or significant charge  In order to further investigate the role of non-Gaussian fluctuations in transitions between different excited states, we have computed several other estimates of the the Marcus theory reorganization energy for each of the processes.In addition to the fitted value, we can also estimate the reorganization energy from linear-response theory as 12,20  LR = 1 2 (⟨Δ⟩  − ⟨Δ⟩  ) (S20) or from the energy gap fluctuations on a given surface  =  or , 12,20  var,  = (Δ − ⟨Δ⟩  ) 2 If the linear response approximation in Marcus theory is valid, then all of these estimates will agree, and therefore deviations between these estimates are indicative of non-Gaussian effects.We see that for the S BD * → S CT and S CT→S 0 the four estimates are in very close agreement, but for the other processes there are more significant deviations between the estimates, although these deviations are still relatively small.The largest deviations are for the processes dominated by the inner-sphere contribution to the reorganization energy, namely the S CT→ T CT, T An * → T BD * and T BD * →S 0 processes, which is unsurprising given the greater importance of asymmetric vibrational frequencies and anharmonicity in the inner sphere contribution.We see however that the linear-response estimate for the reorganization energy generally agrees very well with the fitted value, which suggests that asymmetry in the energy gap fluctuations between surfaces  and  contribute primarily to the free energy differences between these states, and they have only a small effect on the barrier height.The variance based estimators for the reorganization energy also enable us to decompose the reorganization energy in inner sphere, outer sphere, and inner-outer-sphere cross-correlation contributions.We decompose the energy gap Δ into an intramolecular contribution and a remaining environment contribution, Δ = Δ mol + Δ env .We take the molecular contribution

S
FIG. 1. A) Scheme showing the excited state interconversion processes we consider in this work.B) The chemical structure of BD -An.C) Difference densities for each of the excited states calculated at with TDDFT using the SOS-ωB2GP-PLYP functional and def2-TZVP(-f) basis set.

FIG. 3
FIG.3.A-F) Free energy curves for the six A − −− → B processes considered with the reaction A − −− → B labeled on each figure.Points correspond to free energy curves calculated with MBAR and lines correspond to polynomials fitted to the MBAR cumulative distribution functions (see SI for details).G) A snapshot for molecular dynamics simulations on the S 0 potential energy surface.H) A scheme highlighting the processes in A-F.

FIG. 4 .
FIG.4.Spectral distribution () for A) the S CT→S 0 process computed from dynamics on the S CT potential energy surface and B) the T An * → T BD * process computed from dynamics on the T An * potential energy surface.The decomposition into inner and outer sphere contributions and the Debye approximation for the outersphere component is also shown,  D () = (2/)/(1 +  2 D  2 ), where  D = ( ∞ /  ) rel , and  rel is the solvent dipole-dipole autocorrelation relaxation time, and  ∞ /  are the optical and static dielectric constants of the ACN model.

Fig. S2 .
Fig.S2.Free energy curves for the CT→S 0 (A) and T BD * →S 0 (B) processes obtained from umbrella sampling on the S CT and T BD * surfaces.The dashed line shows a quadratic fit to the curves with data c at energy gaps > 1 eV, i.e. close to the  minimum.

TABLE II .
Calculated and experimental rates, quantum yields and triplet lifetime for the photophysics of BD -An.

Table S1 .
Physical properties of ACN calculated for the polarizable and non-polarizable models and experimental values for comparison.

Table S2 .
Gas phase energies of the different excited states of BD -An at their respective minimum energy geometries (from gas phase TDA-TDDFT ωB97X-D3/def2-SVP geometry optimizations) computed at the DLPNO-STEOM-CCSD level of theory with different basis sets. Geometry obtained from calculation with CPCM treatment of ACN solvent. Energies computed from  = 1 ground state DLPNO-CCSD(T) calculations.