Calculation of exchange couplings in the electronically excited state of molecular three-spin systems

Photogenerated molecular three-spin systems, composed of a chromophore and a covalently bound stable radical, are promising candidates for applications in the field of molecular spintronics. Through excitation with light, an excited doublet state and a quartet state are generated, whereby their energy difference depends on the exchange interaction JTR between the chromophore triplet state (T) and the stable radical (R). In order to establish design rules for new materials to be used in molecular spintronics devices, it is of great importance to gain knowledge on the magnitude of JTR as well as the factors influencing JTR on a molecular level. Here, we present a robust and reliable computational method to determine excited state exchange couplings in three-electron-three-centre systems based on a CASSCF/QD-NEVPT2 approach. The methodology is benchmarked and then applied to a series of molecules composed of a perylene chromophore covalently linked to various stable radicals. We calculate the phenomenological exchange interaction JTR between chromophore and radical, which can be compared directly to the experiment, but also illustrate how the individual exchange interactions Jij can be extracted using an effective Hamiltonian that corresponds to the Heisenberg–Dirac–Van-Vleck Hamiltonian. The latter procedure enables a more detailed analysis of the contributions to the exchange interaction JTR and yields additional insight that will be invaluable for future design optimisation.


Introduction
Molecular spintronics is an emerging interdisciplinary research eld that has attracted an increasing amount of attention in recent years since it may allow the development of nanoscale devices with improved performance or new functionalities. [1][2][3] One of the greatest challenges in the eld is to nd new materials that have suitable properties to enable an efficient generation, transport and storage of spin information.
Molecular systems are promising since modern synthesis allows molecules to be tailored with atomic precision. Furthermore, molecular systems have sharply dened electronic states, the manipulation of which could bring functionalities which are not necessarily accessible in the case of solids. [4][5][6] If the desired properties for a possible application are known, suitable molecular systems can thus be developed. However, in order to establish design guidelines for the development of such materials, it is important to know how the material's properties can be inuenced, specically the exchange interaction between spin centres.
Recently it has been shown that photogenerated molecular three-spin systems, composed of an organic chromophore covalently linked to a stable radical, may be ideal candidates to explore the factors governing spin communication on a molecular level. [6][7][8][9] The photophysical processes taking place in these systems aer light excitation are summarised in Fig. 1.
The chromophore is excited to its rst excited singlet state S 1 by absorption of light. This transition corresponds to the HOMO-LUMO transition of the chromophore and exhibits a high transition dipole moment. In the presence of the radical, the triplet ground state T 0 of the chromophore may then be generated by radical-enhanced intersystem crossing (EISC). 10 The requirement for this partially allowed transition is the formation of a correlated doublet state between the chromophore and the stable radical, rather than the coexistence of two isolated systems with different spin multiplicities. The chromophore S 1 / T 0 transition can then be seen as a D 2 / D 1 transition, which is spin-allowed. Finally, the Q 0 state can be generated by intersystem crossing from the D 1 state. [11][12][13][14] If the exchange interaction J TR between the chromophore triplet and the radical doublet surpasses all other magnetic interactions in the system signicantly, the molecular system is said to be in the strong coupling regime. 13,14 Depending on the sign of J TR , the quartet state Q 0 may be lower in energy than the excited doublet state D 1 and vice versa. The energetic order of the doublet and quartet states (sign of J TR ) as well as the magnitude of J TR will determine the magnetic properties of the system and therefore its suitability for different applications. Consequently, to guide the design of these highly modular molecular three-spin systems, detailed knowledge of the factors governing the interaction between the triplet and radical spin centres will be essential.
Experimentally, the excited state exchange interaction in such systems is frequently difficult to determine. Transient electron paramagnetic resonance (EPR) spectroscopy can be used for this purpose, but a reliable value for J TR can only be obtained as long as J TR is smaller or of the same order of magnitude as other magnetic interactions in the system. 15,16 Once the strong coupling regime is reached, where pure quartet and doublet states are formed, the magnitude of J TR has no inuence any more on the shape of the transient EPR spectrum and only a lower bound for J TR is obtained. Experimental data is thus only available for a small subset of the hitherto investigated systems and it seems natural to take recourse to computational methods that allow systematic studies on the properties of these materials to be performed.
In this work, we present a robust computational method for the prediction of the exchange interactions J TR in photoexcited chromophore−radical systems. The approach, based on a CASSCF/NEVPT2 calculation, is generally applicable for any three-spin system and suitable for molecules up to a size of about 250 atoms.
The methodology is applied to a series of perylenes covalently attached to different stable radicals and the trends in J TR are carefully analysed. Further, we illustrate how the phenomenological exchange interaction J TR can be decomposed into three individual exchange interactions J ij (cf. Fig. 1) by using an effective Hamiltonian that corresponds to the Heisenberg-Dirac-Van-Vleck Hamiltonian. 17,18 This approach enables a more detailed analysis of the factors inuencing the magnetic properties of such systems and will be used as a predictive tool to establish guidelines for the optimisation of the design of molecular three-spin systems. Our results further indicate that most of the perylene−radical systems are antiferromagnetically coupled and that the symmetry of the magnetic orbitals is likely to be one of the main factors determining the sign of J TR .

Theory
There are different approaches to calculate the exchange interactions. One of the most popular approaches is broken symmetry DFT, which is a single-determinant method. [19][20][21][22] The exchange coupling constants are estimated by only using the energy of the high spin determinant and the energy of the broken symmetry determinant. If the Yamaguchi approach is applied, also the expectation value of the total spin hS 2 i is required. 23 The disadvantage of this method is that an appropriate functional has to be found for a specic problem in order to compute reliable results, which makes this method unsystematic. 24 In addition, the BS-DFT method delivers unphysical spin densities that require projection methods to be used before physical observables can be calculated. 25 Unfortunately, not every system can be benchmarked, since it is frequently not possible to determine the exchange interactions experimentally. As such, BS-DFT can only be applied to electronic ground states, or a limited number of excited states that would not face variational collapse during the selfconsistent eld optimisation. TD-DFT does not alleviate the problem, since it cannot deal with complex spin couplings. However, it is noteworthy that a certain subset of open shell states can be addressed by the spin-ip TD-DFT method. 26,27 Hence, the calculation of exchange interactions is most rigorously addressed theoretical by using genuinely multideterminantal methods to calculate the ground-and possibly also excited states. 18 Commonly used methods of this type include the CASSCF (complete active space self-consistent eld) and CASCI (complete active space conguration interaction) methods. 18,[28][29][30][31][32][33] Both methods have in common that the orbitals are divided into three sub-classes: the internal orbitals, which are occupied exactly twice, the external orbitals, which are not occupied, and the active orbitals, which can have any occupation number between zero and two. In the subspace of the active orbitals, a full conguration interaction is carried out, yielding a qualitatively correct wavefunction if the active space is properly chosen.
If quantitatively accurate results are required, it is necessary to also account for dynamic electron correlation. Very precise values for the exchange interactions can for instance be achieved with multi-reference conguration interaction (MRCI) methods such as DDCI (difference dedicated conguration interaction) or broken symmetry coupled cluster (BS-CC) methods. 24,32,34 However, since the computational effort is very high, MRCI or CC methods can only be applied to very small molecules. Furthermore, MRCI methods are subject to size- consistency errors, which have a strong inuence on the calculated exchange interactions for any larger systems. [35][36][37] An alternative to the mentioned approaches are multireference perturbation theory methods such as NEVPT2 (nelectron valence state perturbation theory) 38,39 or CASPT2 (complete active space perturbation theory), 40 which are less computationally expensive but underestimate the exchange interactions by 60-80% with a minimal active space in comparison to multi-reference CI methods such as DDCI3. 41 Nevertheless, multi-reference perturbation theory methods provide part of the dynamic electron correlation, which leads to better results than a simple CASCI or a CASSCF calculation. 18,24,42 However, based on the analysis of Calzado, Malrieu and co-workers, it is to be expected that a number of physical effects that are relevant for the correct description of exchange couplings will only occur at higher orders of perturbation theory. 43,44

Choice of the molecules and method
The motivation for this study was the determination of the excited state exchange couplings in a series of perylene-based molecules as shown in Fig. 2. In this series, perylene is covalently bound to seven different radicals, which are commonly used in experimental studies on molecular three-spin systems. [45][46][47][48][49][50][51][52][53][54] The goal was to determine the inuence of the nature of the radical on the magnitude of the exchange interaction(s) and to verify whether systematic trends can be identied.
Perylene was chosen as a representative example for a choromophore since its photophysics is well-known, it is highly photostable and can be substituted easily following established protocols. 55 With regard to the computational effort, perylenebased systems are relatively easy to calculate due to the small size and rigidity of the chromophore, which is convenient when wanting to compare a large number of molecules. We also examined a second series of molecules, which contains only perylene-BPNO and perylene-BDPA structures, but with linkers of different lengths. This second series is shown in the ESI † and allows us to comment on the inuence of the linker length on J TR .
In order to nd a suitable computational procedure, we started by calculating the exchange interaction J TR for two previously investigated perylene diimide (PDI) derivatives, which are linked to BPNO radicals, 48 using different approaches, including restricted open-shell CIS (ROCIS) and CASSCF with and without a QD-NEVPT2 correction. The structures are shown in Fig. 3 and only differ with respect to the orientation of the BPNO substituent attached to the imide position of PDI (meta vs. para). Experimentally, a change in sign of J TR between the para-and meta structures is suggested, inferred from an inverted spin polarisation of the central line (m S ¼ +1/2 4 −1/2 transition) in the transient EPR spectra of the formed quartet states. 13,48 Since quartet state formation was observed, also a lower bound for J TR of $0.4 cm −1 can be given, assuming that J TR is about ten times larger than the zero-eldsplitting of a PDI triplet state of $1100 MHz. 56,57 For the comparison of the different computational approaches, the structures were optimised at the B3LYP/def2-TZVP level of theory. [58][59][60] We used the quasi-restricted orbitals from DFT calculations for the active space selection in CASSCF calculations, whereby the active space was dened by the chromophore HOMO/LUMO, and the radical SOMO. For the CASSCF/ QD-NEVPT2 and ROCIS calculations the def2-TZVP basis set was used. Further computational details are given below.
The experimentally observed sign and change in sign of J TR between the para-and meta structures could be reproduced with the CASSCF/QD-NEVPT2 calculation, showing that the results using this method are qualitatively correct. The ROCIS calculations were judged less reliable since they predicted a ferromagnetic coupling for both structures. Using CASSCF/QD-NEVPT2, a J TR value of 0.49 cm −1 was calculated for PDIpara-BPNO, while for PDI-meta-BPNO a value of −1.3 cm −1 was obtained.
Unfortunately, it is more difficult to judge the agreement regarding the calculated magnitude of J TR , since no accurate experimental reference values are available. In the abovementioned experimental study, jJ TR j was estimated to be larger than 3 cm −1 for both the meta and para compounds based on structural comparisons with a number of different biradical compounds for which J could be determined. 48 Assuming that this estimate is correct, this would imply that the CASSCF/QD-NEVPT2 calculation underestimates J TR by roughly a factor of six which is in line with theoretical studies predicting an underestimation of the exchange coupling by up to an order of magnitude. 41 Although J TR is underestimated, the correct prediction of the sign, sign change and trend is highly promising. In addition, any higher-level methods (e.g. MRCI) that could likely yield more accurate results would be computationally unfeasible for molecules of the considered size. Consequently, we settled on the use of CASSCF/(quasi-degenerate)-NEVPT2 for the calculation of the exchange interactions in the perylene series (cf. Fig. 2). 61 The  individual exchange interactions can then be extracted from the ab intio Hamiltonian in the subspace of neutral determinants. The exact procedure is discussed in the following.

Decomposition of the phenomenological exchange interaction J TR
The magnetic interactions in a quantum mechanical system can be described by the Heisenberg-Dirac-Van-Vleck Hamiltonian: 62-64Ĥ where J ij is the exchange interaction between the electrons i and j, whereby the subscripts 1, 2 and 3 correspond to the chromophore HOMO, radical SOMO, and chromophore LUMO, respectively.Ŝ i andŜ j are the spin operators for the corresponding electrons. In case of a three-electron-three-centre problem, the Hamiltonian reads: Applying the three-electron-three-centre HDVV-Hamiltonian on the neutral determinants, the Hamiltonian can be written in its matrix representation as follows: 17 Diagonalisation gives the eigenvectors, which are the jQ 0 i, jD 1 i and jD 2 i states (with m S ¼ 1/2), and the eigenvalues, which are the energies of these states. 17 with: As can be seen from eqn (6)(7)(8), there are only two linearly independent energy differences, which are dened by three independent constants. As a consequence, it is not possible to calculate the J-couplings only by the energy differences of the eigenstates without any further assumptions. But assumptions could be made on the basis of the structure/symmetry. 17 Assuming that the exchange interaction between the HOMOelectron of the chromophore and the electron of the stable radical J 12 equals the exchange interaction between the LUMOelectron of the chromophore and the electron of the stable radical J 23 , one can dene: If J 13 $ J TR , then the energy difference between the states Q 0 and D 1 can be expressed as: which is the same expression as one would obtain from the Landé pattern using: 41 Unfortunately, these assumptions are not necessarily ful-lled and could lead to a erroneous description of J TR . Instead we will dene J TR as: where the term on the right-hand side corresponds to the energy difference of Q 0 and D 1 , but without the assumption of eqn (10). Here, J TR is directly connected to the energy difference, obtained from the ab initio calculation, whereas according to eqn (11), J TR is equal to J 12 and J 23 or approximately to the average of J 12 and J 23 . In order to extract the individual exchange interactions, an effective ab initio Hamiltonian has to be constructed that corresponds directly to the HDVV-Hamiltonian. Using this approach, one has to project the target states onto a model space, which consists of the neutral determinants jaabi, jabai and jbaai. Then, the projected wavefunctions need to be orthonormalised, such that an effective Hamiltonian can be constructed with the obtained orthonormal wavefunctions. 17,18,65 From the ab initio calculation, the states 4 i and the corresponding energies E i are obtained using the Schrödinger equation within the Born-Oppenheimer approximation: Then, the target states 4 i need to be projected onto the model space S using the projection operatorP S , which is dened as:P where I is an orthonormal basis of the model space S. In this case, it is composed of the neutral determinants jaabi, jabai and jbaai.
Application of the projection operator on the target states 4 i gives the projected states 4 S,i : 4 S;i ¼P S j4 i i: The projected states are not necessarily orthogonal, consequently one needs to orthogonalise them. There are different approaches for this orthogonalisation. In order to obtain a Hermitian effective Hamiltonian, a Löwdin orthogonalisation is applied on the projected basis, using: 66 where S is the overlap matrix with: With the orthonormal basis at hand, the matrix elements of the effective Hamiltonian are calculated using: where I and J are indices for the neutral determinants. Now, the matrix elements of the obtained effective Hamiltonian can be compared directly with the HDVV-Hamiltonian, which allows for the estimation of J 12 , J 13 , and J 23 . The quality of the effective ab initio Hamiltonian can be veried by estimating the shi of the single diagonal elements compared to the diagonal elements of the model Hamiltonian. If the shi remains constant between the single diagonal elements, the Hamiltonian model describes the investigated system well.

Computational details
All structures were optimised at the B3LYP/def2-SVP level of theory. 58,59,67 The structures shown in Fig. 2, except for the tetrathiaryl trityl compound and the BDPA compound, were also optimised at the B3LYP/def2-TZVP level of theory, 60 in order to estimate the inuence of the quality of the structure on the exchange interactions.
Aer every geometry optimisation, a frequency calculation was carried out to verify that the optimisation converged to the ground state structure. All optimisations were carried out using the Gaussian 16 program. 68 For the selection of the active space, a TD-DFT calculation with the RIJCOSX approximation was performed for every structure at the CAM-B3LYP/def2-TZVP level of theory using the ORCA 5.0.3 program. 69,70 Using the results from the TD-DFT calculation, the orbitals dening the three-electron-threecentre problem were determined. We considered the energetically lowest possible transition, which shows the strongest transition dipole moment. Typically, this excited state is composed of transitions within the HOMO/LUMO orbitals of the chromophore and the SOMO of the stable radical and is related to the HOMO-LUMO transition of the chromophore. The photophysical mechanism on which the selection of the active orbitals is based is also illustrated in Fig. 1.
For the calculation of the excited states, which were used to construct the effective Hamiltonian, we carried out a state averaged CASSCF(3,3) calculation with a QD-NEVPT2 calculation on top, in order to account for dynamic electron correlation. The calculations were sped up by the RI-JK approximation for the coulomb and exchange integrals. 71 As starting orbitals we used the orbitals obtained from the TD-DFT calculation. The optimised active orbitals were localised by a Foster-Boys localisation, 72 which allows for an easier interpretation of the excited states. All excited state calculations were also performed using the ORCA 5.0.3 program.
Regarding the choice of the active space, we would like to point out here, that, in studying such systems, it is important to rst understand the nature of the low-lying excited states by carefully assessing their orbital contributions. Initially, this may involve a limited amount of trial and error. However, once the relevant orbitals for the state of interest have been established, we are of the school of thought that the smallest active space that leads to a qualitatively correct description of these states is the preferred one.

Results and discussion
In this section, we will rst focus on the phenomenological exchange coupling constant J TR of eqn (13), the value of which can be obtained almost directly from the ab initio calculation.
We will then consider the individual exchange interactions J 12 , J 13 , and J 23 extracted from the effective Hamiltonian, which corresponds to the HDVV-Hamiltonian. The methodology is presented taking the perylene-BPNO radical as an example. The individual exchange interactions of the remaining compounds of our series can be found in the ESI †. Fig. 4 shows the calculated exchange couplings J TR for all compounds of the series. The values are shown for the optimisation performed using the def2-SVP basis set. The calculations on the structures optimised using a larger basis set (def2-TZVP) provide very comparable results for all structures. With the exception of the perylene-TEMPO compound, a slightly smaller absolute value is predicted for J TR using def2-TZVP. The corresponding values can be found in the ESI †.

Calculation of J TR
As can be seen from Fig. 4, a relatively strong dependence of the exchange interaction on the radical type is predicted. The absolute values of the computed exchange interaction jJ TR j range from close to zero (e.g. perylene-TEMPO) to >10 cm −1 for perylene-BPNO.
Most perylene−radical systems are antiferromagnetically coupled. Only the perylene-proxyl system and the perylene-eTEMPO system are ferromagnetically coupled. Interestingly, both of these radicals are asymmetric towards the chromophore−radical bonding axis, which is not the case for the remaining antiferromagnetically coupled compounds. This suggests that the symmetry of the magnetic orbitals might be crucial for the resulting sign (and magnitude) of the exchange interaction J TR . Compared to the inuence of symmetry, the extent of electron delocalisation appears to play a minor role when visually comparing the tetrathiaryl compound and the BDPA compound, or the proxyl compound and the TEMPO compound. As shown in the ESI, † a dependence on the linker length can also be observed, whereby a larger linker length causes a reduction in the exchange interaction. This trend is to be expected from the exponential distance dependence of exchange interactions since the main effect of a longer linker will be the increase in distance between the interacting spin centres. An exception would only occur here if the radical electron were also delocalised over the linker, which is conceivable, since the energy difference between the linker HOMO and the radical SOMO decreases with increasing linker length, making the mixing of these orbitals more likely.

Extraction of the individual exchange interactions
In order to better understand the rather large differences in magnitude and sign of the calculated exchange interactions J TR within the perylene series and to be able to interpret them in a meaningful way, it is advisable to consider the individual contributions to J TR . We will illustrate the procedure using the perylene-BPNO system as an example.
First, the target states of the ab initio calculation will be projected onto the subspace of the neutral determinants and normalised (see also section 5 in the ESI † for additional details and an excerpt of the corresponding output le). We obtain: The projected target states are not necessarily orthogonal, consequently a Löwdin orthogonalisation is applied on these states using eqn (17). The states are then renormalised, which gives the orthonormalised states: The quality of the orthonormalisation can be veried by simply calculating the overlap matrix S. The off-diagonal elements of S should ideally be zero, or at least close to zero.
For the calculation of the effective ab initio Hamiltonian, we need to consider the calculated energies of the target states, whereby we will set the energy of the Q 0 state to zero.
Using the coefficients of the target states and their energies, the effective Hamiltonian can be constructed according to eqn (19), which can also be represented as: where c are the coefficients for the determinants I of the states i.
Using this Equation, we obtain the following effective Hamiltonian: All values are given in cm −1 . Now the individual exchange couplings can be extracted from this effective Hamiltonian, since it has a one-to-one correspondence to the Heisenberg-Dirac-Van-Vleck Hamiltonian. The following exchange couplings are obtained: It should again be mentioned, that the subscripts refer to the predened spin centres. Fig. 5 shows a visualisation of the corresponding localised spin centres. In the case of the perylene-BPNO system, the size of the individual exchange interactions can be explained by such a visualisation of the spin centres in a relatively straightforward manner. Electrons from spin centres that are closer together show a larger exchange interaction compared to those that are farther away. By visual analysis of the orbitals, without numerical considerations, it can be stated that J 12 will likely be larger in magnitude than J 23 . The interpretation of the sign of the individual exchange interactions, on the other hand, is more difficult. However, as already noted above, the sign is likely to depend on the symmetry of the individual spin centres.
We would like to emphasise here, that the orbitals shown in Fig. 5 are their localised representations. In order to ensure a consistent assignment of the spin centres for all molecules in our series, the spin centres were not assigned according to their occupancy number but according to their orbital compositions, which allows for a better comparison.
The advantage of the decomposition method is that we are now in a position that allows us to better analyse and understand the different contributions to the exchange interaction J TR , their individual importance and how we can modify them, if necessary. If J 13 $ J TR , then we can approximately express J TR as the average of J 12 and J 23 : As a consequence, if we want to modify the exchange interaction J TR , we have to focus mainly on J 12 and J 23 . In order for the exchange interaction J TR to become minimal, J 12 and J 23 must either trivially approach zero or they must cancel each other out. However, in order for the exchange interaction J TR to be maximised, J 12 and J 23 must either have the same sign or one of the two exchange interactions must be signicantly larger. The latter can presumably be achieved by lowering the molecular symmetry, e.g. by appropriate asymmetric substitution.

Conclusions
To our knowledge, this is the rst study dealing with the problem of calculating the exchange interaction J TR in the excited state of molecular three-spin systems, which turned out not to be an easy task. For example, the frequently used broken symmetry DFT method would fail on such systems, since they cannot be sufficiently described with a single ground state determinant. Furthermore, high-level methods as for example the DDCI3 method are only feasible for very small structures and therefore not applicable for the molecules investigated here. We have found a robust method that is relatively fast and easy to apply while providing reliable trends for the exchange interaction in excited triplet−doublet systems.
We could show that most of the investigated excited perylene−radical systems are predicted to be antiferromagnetically coupled, whereby asymmetric chromophore−radical systems (regarding the chromophore−radical bonding axis) exhibit a ferromagnetic coupling. By extracting the individual exchange couplings using an effective Hamiltonian that corresponds to the Heisenberg-Dirac-Van-Vleck Hamiltonian, we could analyse the phenomenological exchange interaction J TR in more detail. We showed that the exchange interaction J TR can be expressed as the average of the two chromophore-radical exchange interactions (J 12 and J 23 ). As a consequence, any future optimisation of J TR should mainly focus on controlling these two contributions.
Future investigations in our group will focus on the decomposition of the individual exchange couplings, which are still effective parameters. From this decomposition we will get valuable information on the direct exchange and the kinetic exchange contributions, 18,43,73 which is needed in order to understand the signs of the individual exchange interactions. Further, we will examine the symmetry of the magnetic orbitals and its inuence on the exchange interactions.

Data availability
The data supporting the ndings of this study are available within the article and in the ESI †.

Author contributions
Quantum chemical calculations, investigation, methodology, formal analysis, data validation and visualisation, original dra writing, review and editing M. F.; Supervision, review and editing of the manuscript F. N.; Conceptualisation, project administration, funding acquisition, supervision, original dra writing, review and editing S. R.

Conflicts of interest
There are no conicts to declare.