A theoretical study of the time-resolved x-ray absorption spectrum of the photoionized BT-1T cation

The molecule 4-(2-thienyl)-2,1,3-benzothiadiazole (here denoted BT-1T) consists of an electron-accepting benzothiadiazole unit and an electron-donating thiophene unit (see Fig. 1). Several π-conjugated polymers, including either one or both of these units, have been scrutinized as potential organic solar cells^1–4 and generally for organic electronic devices.^5,6 Particularly, applications in organic solar cells are of interest. Such solar cells are becoming a feasible alternative to the more established silicon- or metal-based ones due to their lower production costs^7,8 and with organic solar cells now reaching efficiencies of around 18%.^9–12 An important property of an organic solar cell is its ability to form charge transfer states.^13,14 Thus, this process must be investigated for any potential new solar cell in order to determine its suitability.¹³ 

Although molecules including BT-1T units have already been investigated both experimentally and theoretically,^15,16 the small BT-1T molecule has only recently been studied theoretically on its own by Khalili et al., who considered not only the photoexcitation¹⁷ but also the photoionization of the system.¹⁸ So far, experimental studies of this system have not been published. Since experimental studies do typically rely on theoretical support for a full interpretation of the often complicated experimental spectra,^19–22 further theoretical studies are, nonetheless, still required, especially at higher levels of theory.

Khalili et al.^17,18 studied the BT-1T unit by simulating its time-resolved x-ray absorption spectra (tr-XAS). X-ray spectroscopy and scattering offer time-resolutions on the sub-picosecond timescale characteristic of molecular reactions.^23–28 Furthermore, XAS, especially excited state XAS, has proven sensitive to both electronic and structural changes in small organic molecules.^22,29,30 In Ref. 18, the authors explored, in particular, the possibility to use tr-XAS to follow the charge dynamics of the BT-1T system after photoionization with a vacuum ultraviolet (VUV) pulse. They prepared the BT-1T⁺ cation in a state corresponding to the removal of one of the HOMO-3 electrons in the neutral molecule (see Figs. 2 and 3).

Khalili et al.¹⁸ employed Hartree–Fock (HF) level calculations using a small basis set, which is a rather low level of theory; calculations at a higher level of theory are, therefore, needed to determine the validity of the results. Moreover, the dynamics of the process was studied with trajectory surface hopping (TSH), which treats the non-adiabatic process stochastically. While this is often a good approximation, it is also of interest to investigate the effect of using quantum dynamics. We will, therefore, study the BT-1T system quantum dynamically for the process of tr-XAS (probe step) following the VUV photoionization (pump step) of the neutral species. We will investigate not only whether this type of experiment is suitable for studying the important charge transfer reaction, but also how much the conclusions are affected by resorting to a higher level of theory. Khalili et al.¹⁸ employed ionic potential energy surfaces (PESs) based on Koopmans' theorem (KT). Instead, we will investigate the nuclear dynamics by applying time-dependent density functional theory (TDDFT) directly to the cation. Unlike HF, TDDFT includes some electron correlation,³¹ without increasing the computational cost beyond feasibility. In general, TDDFT has proven a very successful method due to its low computational cost and overall good performance,^32–34 even though the choice of functional can greatly affect the results, and its effect is often difficult to predict.^35,36 Figure 3 pictorially illustrates the differences between the approach used by Khalili et al.¹⁸ and the one used in the present study when considering the molecular orbitals of the neutral molecule.

Following the nuclear dynamics simulations, XAS calculations are then carried out utilizing TDDFT as well as frozen core-valence separated equations-of-motion coupled cluster singles and doubles (fc-CVS-EOM-CCSD).³⁷ This latter method is computationally more demanding, but the quality of the results obtained with CC methods is more predictable and can be systematically improved. Indeed, it is in the CC family of methods that the gold standard for calculations is to be found when the ground state is dominated by a single configuration.^38–40 While not being the gold standard, CCSD still yields results of relatively high quality and has become one of the methods of choice for accurately computing molecular properties at a reasonable cost.⁴⁰ 

The paper is organized as follows. Section II collects all information on the computational approaches utilized for this study. In Sec. III, we present the relevant molecular orbitals (MOs) and valence transitions within the cation, while in Sec. IV, we report the calculated charges on both the sulfur atoms and the two units of the molecule in the electronic states of interest. Next, the dynamics simulation is presented in Sec. V. Here, we perform a TSH simulation⁴¹ for a limited number of trajectories in order to determine the normal modes to be included in a quantum dynamics simulation using the multi-configurational time-dependent Hartree (MCTDH) method.^42,43 The results of these simulations are compared and evaluated, before the XAS calculations are considered in Sec. VI. To simulate the tr-XAS, the predicted average geometries and state populations from the MCTDH simulation are used in Sec. VI C. Here, both the importance of the geometrical changes, which are expected to be small due to the rigidity of the molecule, and the electronic state are evaluated.

Geometries: The molecular geometry of the molecule before ionization was determined by optimizing the molecular structure of the neutral species on the ground state using the ORCA-5.0.0 program package.^44,45 The optimization was carried out at the DFT level of theory employing the BHHLYP functional and the cc-pVDZ basis set.^46,47

Since the photoionization, starting the studied process, is fast, the initial geometry of interest (at a time delay of zero fs) is that of the neutral species. The geometries for non-zero time delays are determined as the average structures (i.e., centroid geometries) found in the quantum dynamics simulation.

Dynamics: For all calculations related to the dynamics simulations, the (TD)DFT/BHHLYP/cc-pVDZ level of theory was used. For excited state calculations, the Tamm–Dancoff approximation (TDA) was employed. This approach has been used previously with good results.⁴⁸ 

First, the normal modes were determined analytically for the optimized structure. To determine the relevant normal modes to include in the MCTDH calculation, a TSH calculation was carried out in the SHARC-2.1.1 program package,^41,49,50 interfaced with ORCA-4.2.0.^44,51 The older version of ORCA was utilized here, since the auxiliary python-scripts necessary to perform the SHARC calculation were not adapted to the newer version. In the TSH simulation, 25 trajectories were simulated for 400 fs. This number of trajectories is clearly too small to obtain good statistics for predicting quantum yields and lifetimes. It has, however, previously been shown that a limited number of trajectories can qualitatively determine the normal modes of interest⁴⁸ as long as one is not concerned with rare events. The TSH simulation was carried out from a Wigner distribution of geometries and photoexcitation of the cation to simulate initial conditions. The final states considered are the same with respect to the occupation of neutral MOs, as those obtained by a photoionization of the neutral species (see Fig. 3). The gradient in all states, and hence the PES points, is like the initial conditions, generated by considering photoexcitation of the ion. This way of generating the PES is employed even if initial conditions are chosen based on photoionization of the neutral species (such a simulation can be found in the supplementary material, Sec. S3.2). The (adiabatic) PES in the TSH simulation was determined on the fly. The simulation started in the third excited state (spin state 7 in the diagonal representation). The calculation included spin–orbit couplings and was, therefore, performed with SHARC dynamics,⁵² i.e., based on the diagonal representation of the Hamiltonian. This is an alternative to the molecular Coulomb Hamiltonian (MCH) representation,⁵² which is known as regular surface hopping. Couplings between states were approximated based on wavefunction overlaps, and an energy based decoherence correction was employed.

The MCTDH simulation was carried out along selected normal modes of interest, determined based on the TSH simulation (see Sec. V B and Sec. S3.1 in the supplementary material). The (diabatic) PESs for the MCTDH calculation were fitted to a vibronic-coupling Hamiltonian^53,54 with a harmonic oscillator zero-order term and dimensionless mass-frequency scaled normal coordinates. Harmonic oscillator eigenfunctions were used as a basis for the PES fit. The fit for each normal mode was based on 31 points along the given mode for displacements Q = −15 to Q = 15, as some modes showed displacements of $|Q|>10$⁠. The fitting parameters are given in the supplementary material (Sec. S3.4). All necessary points were calculated with ORCA-5.0.0, and the quantum dynamics simulation was then carried out with the MCTDH-8.4.12 program package.⁵⁵ Input for the latter program was based on a local development version of the VCHAM program,⁵⁶ modified to accommodate output files generated by ORCA-5.0.0.

The MCTDH calculation employed a primitive basis consisting of 45 basis functions for each normal mode chosen. The D₀ state for each mode was described by five single particle functions, while 10 functions were used to describe the remaining states for each mode. The grid populations were investigated and found to be below $10−5$⁠, indicating a converged primitive basis. Likewise, the natural weights were investigated to ensure convergence with respect to the single particle functions. Again, sufficiently small numbers were found (of the orders of magnitude: $10−7$– $10−3$⁠). In addition, multimode single particle functions were used to combine modes in order to reduce computational costs.

Charges: Calculations of CHELPG (CHarges from ELectrostatic Potentials using a Grid based method⁵⁷) point charges were performed in Gaussian 16⁵⁸ at the BHHLYP/cc-pVDZ level of theory, in order to obtain the excited state charges. Such CHELPG charges have previously been successful in applications for classical nuclear dynamics.⁵⁹ 

X-ray absorption spectra: The XAS were calculated at both the TDDFT level of theory with the BHHLYP functional and at the fc-CVS-EOM-CCSD level of theory, adopting the aug-cc-pVTZ basis on sulfur and cc-pVDZ on the remaining atoms. Q-Chem-5.4.2⁶⁰ was used for these calculations.

The XAS spectrum of each electronic state of the cation was determined using TDDFT/fc-CVS-EOM-CCSD. These calculations were based on the Kohn–Sham/CCSD wavefunction of the given state optimized by employing the initial maximum overlap method^61,62 (IMOM) procedure. Here, the excited state wavefunction is generated by solving the self-consistent field equations, imposing as the constraint that the occupied orbitals have the largest possible overlap with an initial set of guess orbitals, rather than utilizing the Aufbau principle. This strategy has previously been employed with success in a number of cases.^30,63–65 In our calculations, the choice of initial guess was based on the orbitals determined for the ground state of the neutral molecule. The occupied orbitals were then specified to be the same with one exception, since the cation has a singly occupied orbital. The choice of this singly occupied orbital for each investigated state is discussed in Sec. III and can also be anticipated based on Fig. 3. The IMOM approach has the disadvantage of generating excited states that are not orthogonal to the ground state. On the other hand, the XAS of the excited states can be straightforwardly computed as one-photon absorption of the IMOM-optimized state. Thus, one can obtain not only the core to singly occupied MO (SOMO) transitions but also transitions from core orbitals to virtual MOs. An alternative approach for simulating excited state XAS is to compute the intensities as the transition strengths between two excited states, namely, the initial valence excited and the final core excited state. The corresponding transition energies are then computed as energy differences between the excitation energies of the two states from the ground state. At the CCSD level, this “quadratic response” approach can, however, only reliably describe core-to-SOMO transitions, as only these correspond to a single photon process from the ground state. An equivalent approach for TDDFT is currently not available in Q-Chem-5.4.2.

A default convergence criterion of 10⁻⁶ was used for the excitation energy TDDFT calculations, while a convergence threshold of 10⁻⁵ was applied for the computationally more demanding EOM-CCSD calculations. All calculations on the cation employed an unrestricted formalism. TDA was not used for these calculations.

Orbitals: Ground state MOs are plotted based on Gaussian 16 calculations. Excited state MOs were, however, found by utilizing the Q-Chem-5.4.2 program package. All natural transition orbitals (NTOs) were determined based on calculations in Q-Chem-5.4.2. NTOs and MOs were plotted using Molden⁶⁶ and a contour value of 0.01.

The MOs of interest are the HOMO, HOMO-1, HOMO-2, and HOMO-3 displayed in Fig. 4 for both the neutral molecule and the cation. The obtained MOs for the neutral molecule are in good agreement with those presented by Khalili et al.^17,18 The MOs of the cation will be discussed further below. Throughout the text, we will refer to the orbitals for the neutral molecule with a subscript “0” and to those for the cation with a D_n subscript, indicating the considered state. Hence, the HOMO for the neutral molecule is denoted HOMO₀.

To generate the wavefunction of the cation, a β-spin electron was removed from the neutral molecule. Removal of an electron from the HOMO results in a slightly spin contaminated D₀ state with $⟨S2⟩=0.8163$ at the BHHLYP/cc-pVDZ level of theory (a pure doublet spin state has $⟨S2⟩=0.75$⁠).

Taking the D₀ state of the cation as the reference, we also computed the NTOs for its first three valence transitions, which are shown in Fig. 5. These NTOs closely resemble the MOs shown in Fig. 4 (left). The particle NTO for all valence excitations corresponds to HOMO₀, which must, therefore, be singly occupied in the cationic ground state denoted D₀. D₀ corresponds to the hole in HOMO₀, D₁ corresponds to the hole in (HOMO-1)₀, D₂ corresponds to the hole in (HOMO-2)₀, and, finally, D₃ corresponds to the hole in (HOMO-3)₀. For consistency, we shall throughout denote the state associated with a hole in the (HOMO-n)₀MO as D_n.

The MOs and NTOs determined by using the methods and basis sets employed to calculate the XAS are found to agree well with the MOs and NTOs shown here (see the supplementary material, Secs. S2.2 and S2.3).

The four highest-energy occupied MOs of the cation, in each of the four different electronic states considered in the XAS calculations, are shown in columns 2 to 4 of Fig. 4. We only explicitly generate the states D₁-D₃ when calculating the XAS. Therefore, the MOs of the cation are shown for the basis sets also employed in these calculations, i.e., the aug-cc-pVTZ basis on sulfur and cc-pVDZ on the rest. Note that the Unrestricted Hartree–Fock (UHF) MOs used in the CC calculations on the cation (Sec. S2.3) are in good agreement with the ones shown here. Furthermore, the MOs of D₀ clearly resemble those obtained at the BHHLYP/cc-pVDZ level of theory (see Sec. S2.1). We highlight the different energetic ordering of the MO orbitals of equivalent character in the neutral molecule and in each of the four states of the cation. For these cationic states, the orbital with the highest energy (the $HOMODn$⁠) is always the singly occupied one, which corresponds to one of the four highest occupied MOs of the neutral molecule [the (HOMO-n)₀]. In short, the SOMO is in all cases $HOMODn$⁠, which has the same character as (HOMO-n)₀.

Based on the characters of the MOs obtained for D₀ (see Fig. 4), we note that generating, e.g., D₁ corresponds to a hole in $(HOMO‐3)D0$⁠. We inspected the weights of the different transitions contributing to the D₁ excitation and found that the main contribution to this excitation is, indeed, from the transition originating in $(HOMO‐3)D0$ and with final state in $HOMOD0$⁠. We did, however, note from these weights that each state cannot be associated with only one transition (although one is found to dominate in each case). This might partly explain why the D₁ state is associated with a transition originating in the $(HOMO‐3)D0$⁠. Despite the lower energy of the $(HOMO‐3)D0$ MO compared to $(HOMO‐1)D0$⁠, the total energies of the D_n states increase as expected from D₀ to D₃ (see Table III). As previously mentioned, D_n is here generated by moving one electron from the MO with the character of (HOMO-n)₀ to that with the character of HOMO₀. Thus, our chosen naming convention of the states is also in line with the energetic ordering at the ground state optimized structure.

The reordering of MOs indicates that we cannot apply the simple picture of the MOs retaining a predetermined ordering upon ionization, which is generally assumed, e.g., when using Koopmans' theorem. Such reordering of the orbitals upon excitation has also been observed in a recent study by Schmerwitz et al.⁶⁷ The final states obtained in our approach should, however, be noted to correspond well with what would also be expected based on Koopmans' theorem. Our findings do indicate, though, that analyzing the process in terms of one electron moving stepwise down through the MOs might be too simple.

In order to investigate the charge transfer reaction, we computed the CHELPG point charges of the cation in the ground state as well as in the three lowest lying excited states (see Table I). We observe only small changes in the charge on the sulfur atom located in the benzothiadiazole unit (labeled S_BT) across the considered states. A large change is, however, observed on the sulfur atom in the thiophene unit (labeled S_T) in D₁. We note that electron density is, thus, moved away from S_T in D₁, as was also concluded by Khalili et al.¹⁸ One could, therefore, in principle, observe the migration of electron density away from S_T by monitoring the de-excitation process from D₃ to D₁. Another observable indication of the electron density migration might be the strength of the core transitions originating on the different sulfur atoms, since a high intensity is generally associated with high electron density on the atom on which the core transition originates. This will be further considered in Sec. VI C. It should be noted, however, that the changes in charge density on S_BT do not correspond to the changes of the charge density on the entire benzothiadiazole unit, although a good qualitative agreement is observed for the sulfur atom in the thiophene unit. While the charge density on S_BT, as mentioned, does not display any significant changes, it is clear that negative charge is moved from the thiophene unit to the benzothiadiazole unit from D₃ to D₁. Hence, while S_T might be the most important atom to consider for the charge transfer in the thiophene unit, one cannot simply inspect the S_BT atom of the benzothiadiazole unit. Instead, the entire unit should be scrutinized.

When performing the MCTDH calculation (as well as the TSH simulation), the PESs are calculated based on excited state energies of the cation. The starting point of the simulation is the D₃ excited state. The zero point energy of the calculation is set to the ground state energy of the neutral molecule at its equilibrium geometry, as this is the starting point of the corresponding experiment.

The normal modes that appear to be of interest in the TSH simulation are evaluated further in order to choose the modes to be included in the MCTDH calculation (see the supplementary material, Secs. S3.1–S3.3).

The evolution of populations in the TSH simulation is studied by considering the states based on energy ordering (i.e., adiabatic states), which is the natural choice for TSH simulations. It is also possible to consider the states based on oscillator strengths. The latter can be useful for excited states that are close in energy but have oscillator strengths of different magnitudes. This is the case here (see Table II), however, the oscillator strengths are seen to change significantly during the reaction (cf. Tables S3, S20, and S21). Hence, it is difficult to define the states in terms of non-overlapping oscillator strength intervals. This approach has, thus, been deemed inappropriate for the system considered here. The populations based on energies are determined by the fraction of trajectories in each state at a given time. Moreover, to better compare the evolution of the populations determined here to the one determined using MCTDH, the populations are also transformed to a diabatic basis using auxiliary SHARC scripts. These rely on overlap matrices as well as wavefunction coefficients in the diabatic basis.

When simulating XAS later on, one must be careful to verify whether the first and second excited states do, indeed, correspond to a hole in (HOMO-1)₀ and (HOMO-2)₀, respectively. At the geometries associated with our investigated time delays, the NTOs for the three valence excitations appear unchanged compared to the NTOs at the initial geometry (see Figs. S3, S46, and S47). This indicates that we can associate the D_n state with a hole in (HOMO-n)₀ throughout.

The TSH simulation gave rise to three trajectories jumping to the D₀ state. Plots of the evolution of the populations can be seen in Fig. 6. The diabatized populations were also determined using an auxiliary SHARC script, and the results are shown in the plots in Fig. 7. The agreement between Figs. 6 and 7 is excellent, although Fig. 7 shows a smoother evolution of the populations compared to Fig. 6. This result indicates that, in this case, the diabatic and adiabatic pictures give rise to the same states, further justifying our choice of denoting the different states as D_n for an empty (HOMO-n)₀MO throughout.

The TSH simulation shown here exhibits the same rapid decrease in the D₃ population and simultaneous increase in the populations of D₂ and D₁ as also shown by Khalili et al.¹⁸ We, however, obtain a significantly lower yield of D₀ compared to the previous study. This discrepancy could be a result of the different levels of theory used. Furthermore, the wavefunction of the ionized species was determined by Khalili et al.¹⁸ based on Koopmans' theorem,^68,69 without including additional relaxation. The otherwise good agreement leads us to conclude that the normal modes of interest for the reaction at this level of theory can, indeed, be determined based on the presented TSH simulation. We also note that the evolution of the populations appears very similar in the two additional TSH simulations presented in the supplementary material, Secs. S3.2 and S3.3. In Sec. S3.2, photoionization is used to generate initial conditions, while in Sec. S3.3, the nuclear dynamics was carried out in the MCH representation for initial conditions generated based on photoexcitation of the cation. This indicates that the results presented are relatively robust with respect to small changes in the initial condition generation and the state representation, in which the nuclear dynamics is considered.

When determining which normal modes to include in an MCTDH calculation based on a TSH simulation, many types of analyses can be made;⁷⁰ one can consider the activity of the modes as well as coherence and correlation. The linear vibronic-coupling (LVC) parameters can also be investigated. Based on all these types of analyses (see the supplementary material, Sec. S3.1), ten normal modes of interest were identified: modes 1, 2, 8, 16, 24, 33, 44, 46, 47, and 48 (labeled based on frequencies of increasing magnitude). These ten modes were also found to be of importance in the two additional TSH simulations shown in Secs. S3.2 and S3.3.

The PESs of the ten chosen normal modes were determined using the VCHAM procedure. Here, the PESs were fitted as harmonic oscillators to best reproduce the adiabatic PESs obtained by quantum chemistry calculations. The fitting parameters for the vibronic-coupling Hamiltonian considered^53,54 can be found in the supplementary material (Sec. S3.4).

The fitted PESs are shown in Fig. 8. It is observed that, while the fit for normal modes 8, 16, 46, 47, and 48 appears to be in excellent agreement with the calculated points, the fits for the remaining normal modes, however reasonable, show some discrepancy. For normal mode 1, the fit shows a large discrepancy beyond $|Q|=5$⁠, since the shape is not a parabola as assumed in the harmonic oscillator approximation. From the TSH calculation, however, the displacements along normal mode 1 are found to be small (⁠ $|Q|<0.8$⁠), and, thus, the fit can, nonetheless, be utilized. The fit of the PES for normal mode 2 is also off between $|Q|=5$ and $|Q|=10$⁠, since the predicted PES is narrower than the calculated points would indicate. Again, however, only small displacements were observed in the TSH simulation (⁠ $|Q|<1.0$⁠). If a better agreement at larger Q had been required, one might employ a different model for the fit, instead of the harmonic oscillator. For the last modes, 24, 33, and 44, the general shapes of the PESs are in good agreement with the calculated points, except for small shifts of the minima compared to the quantum chemistry calculations.

When determining which normal modes to include in the MCTDH calculation, we look for two types of modes. First, we are interested in modes that appear to displace the molecule away from the equilibrium (tuning modes). Second, modes that potentially give rise to crossing of states (coupling modes) must be considered. Modes 1 and 2 both display a symmetric set of PESs with the minimum energy centered around a zero displacement. While mode 1 does show close-lying PESs for the four states, neither mode shows indication of PESs crossing for the different states until displacements larger than those found in the TSH simulations (⁠ $|Q|≥7$⁠), and, thus, these modes are of less interest. Neither of the remaining modes is completely symmetric about their energy minima, although modes 46 and 48 appear almost symmetric. In addition, the minimum of the PESs occurs slightly displaced from equilibrium. Thus, these might all be of interest. They all display close-lying states; however, only modes 33, 44, 46, and 47 appear to display degeneracies between states at low Q. Investigating the selected modes in terms of their effect on the S_BT–N, S_T–C, and C₂–C₃ bonds as well as the C₁–C₂–C₃ and C₂–C₃–C₄ angles, and the C₁–C₂–C₃–S_T dihedral angle (see Fig. 1) shows that all modes affect at least one of these geometrical parameters (see Sec. S3.6). These parameters, also considered by Khalili et al.,¹⁸ were all found to be affected by the chosen set of modes. The characters of the modes were further investigated by visualization. It was found that they account for in plane stretches of the molecule as well as in plane movement of the thiophene- and benzothiadiazole-units toward each other (see the supplementary material, Sec. S3.7).

An MCTDH simulation was carried out including all eight modes of interest (8, 16, 24, 33, 44, 46, 47, and 48), where multimode single particle functions were used to combine modes 44 and 48 and modes 46 and 47. These modes can be combined, since they have similar vibrational frequencies. Moreover, their coupling parameters (see Sec. S3.4), determined by the VCHAM program, were observed to be larger than for other mode combinations. In Fig. 9, the evolution of the populations over time can be seen for the MCTDH calculation.

The MCTDH simulation shows the same general evolution of the populations as observed in the TSH simulation (Fig. 7). We do, however, note that the population in D₃ decreases much more slowly, and that the evolution overall appears smoother. The latter is likely a result of the few simulated trajectories in the TSH calculation. The predicted half-life of the D₃ state in our quantum dynamics simulation is found to be approximately 41.6 fs based on an exponential fit. This is significantly longer than the approximately 5 fs obtained from the TSH simulation. However, the overall trends in the evolution of the populations with the initial drop in D₃ population and a rise in D₂ and D₁ are in good agreement with the TSH simulation. The fact that the D₁ population then appears to become dominant is also observed in both the previous and current TSH simulations. While the population in D₀ was also observed to be low in the TSH simulation carried out in this study, the D₀ population was found to increase to approximately 90% in 400 fs in the study by Khalili et al.¹⁸ This discrepancy must be attributed to differences in the simulation methods. We emphasize here the differences regarding the on-the-fly generation of the PES, which, in their study, was based on Koopmans' theorem. Furthermore, the choice of quantum chemistry method and, particularly, the lack of electron correlation in their calculations might have affected the result. The evolution of the D₁–D₃ states found by Khalili et al.¹⁸ appears, however, to be well reproduced by the MCTDH calculation except for the half-life of D₃. This was found to be 8 fs by Khalili et al.¹⁸ (in agreement with our TSH simulation). The longer lifetime of the D₃ state in the MCTDH calculation is attributed to differences in methodology: the PES generation and approximations of the non-adiabatic couplings as well as the dynamic treatment of the transitions. Another reason is the fact that the MCTDH simulation is not full dimensional. Including additional degrees of freedom (DOFs), and hence restricting the movement of the molecule less, would, to some extent, decrease the lifetime. To determine the effect of restricting the degrees of freedom, a TSH simulation was performed using an LVC Hamiltonian based on the linear fitting parameters used for the MCTDH simulation. Thus, this simulation was carried out in a reduced-dimension along 8 degrees of freedom only. This simulation showed the same trends with respect to lifetimes as the MCTDH simulation (see Sec. S3.5). Thus, despite the relatively rigid molecular structure, the difference in lifetime of the D₃ state is attributed mainly to the reduced DOFs in the MCTDH simulation. We do, however, note that while the overall trends of all the TSH simulations are the same, the evolution of the populations differs, particularly in the adiabatic representation. This could be related to the few trajectories simulated. Additionally, differences between the TSH simulation based on the LVC Hamiltonian and the MCTDH simulation could partly be related to differences in methodology. Since the evolution of the populations of D₃–D₁ in the MCTDH simulation is in better agreement with the simulation by Khalili et al.,¹⁸ where 100 trajectories were considered, we have more faith in the MCTDH based populations although we emphasize that the timescale is here too long. We, therefore, choose to base our tr-XAS calculations on the MCTDH simulation, to better reflect the population of states. Since the overall trends in the evolution of states are the same in both TSH and MCTDH based simulations, the 8 normal modes chosen do reflect the most important nuclear movement of the system although additional DOFs were shown to speed up the process. We, therefore, conclude that the MCTDH based geometries are sufficient to investigate the effect of small geometrical changes of the system on the XAS compared to the changes induced by a change in the state populations.

To simulate the sulfur K-edge XAS before the photoionization process, we calculated the spectrum for the neutral molecule at its equilibrium structure. Experimental XAS spectra have not been reported. The results are shown in Fig. 10. The numerical values of the excitation energies and oscillator strengths are tabulated in Tables S12 (BHHLYP) and S13 (CCSD). The transitions are characterized based on their dominant NTOs in Figs. S30 and S31/S32.

Relatively modest differences are observed between the XAS spectral profiles obtained at the fc-CVS-EOM-CCSD and BHHLYP levels of theory. There is, however, a large overall energy shift. Such large energy shifts are commonly observed for TDDFT,⁷¹ whereas fc-CVS-EOM-CCSD generally requires relatively small shifts to align with experiment.^37,72 Thus, throughout, the TDDFT results have been shifted to align with the fc-CVS-EOM-CCSD results.

According to both electronic structure methods, four electronic transitions are responsible for the main band around 2479 eV. The first three transitions almost overlap at the BHHLYP level, whereas the fourth is separated by roughly 0.5 eV. This yields the double peaked shape of the first spectral band in the simulated spectrum. The four transitions are more regularly separated from each other at the fc-CVS-EOM-CCSD level, thus giving rise to one single peak when broadened with the chosen half width at half maximum (HWHM).

The characters of the four transitions are similar for the two methods, according to the NTOs in Figs. S30 and S31/S32. The first transition is assigned as a core excitation from the S_T atom to an antibonding orbital. The electron density of this orbital is predominantly around the thiophene moiety. Furthermore, a noticeable contribution from an in-plane 3p orbital on the S_T atom is observed. We label this transition $ST→σ1,T*$⁠. The second transition is 1s(⁠ $SBT)→π1,BT*$⁠, where the $π1,BT*$ orbital is clearly localized on the benzothiadiazole unit. An equivalent 1s(⁠ $ST)→π2,T*$ transition, with $π2,T*$ mainly on the thiophene unit, is found as the third transition. The fourth transition is the S_BT equivalent of the first one, i.e., 1s(⁠ $SBT)→σ2,BT*$⁠, with noticeable contribution of an in-plane 3p orbital of S_BT. The fc-CVS-EOM-CCSD calculation, furthermore, predicts an additional low intensity signal at higher energy, which is not captured by BHHLYP within the given number of computed transitions.

The strongest transition of the main spectral band is the one originating on S_BT, but a transition of comparable intensity originating on S_T also contributes greatly to the signal. Hence, the signals from the two S atoms might be very difficult to distinguish, when only considering the ground state XAS.

The S K-edge XAS of the cation has been determined for the four states of interest, D₀–D₃, at the two considered levels of theory, see Fig. 11. As the computations are based on UHF and Kohn–Sham calculations, the results show spin contamination. The $⟨S2⟩$ values can be found in Table III. Here, we note that spin contamination in both the BHHLYP and CCSD level calculations is modest. The NTOs of the main peaks can be seen in Sec. S4.2. It is noted that while the D₀ and D₃ states might be difficult to distinguish, due to the low intensity of the core-to-SOMO transitions, the D₁ and D₂ states are clearly distinguishable. While their core-to-SOMO transitions do fall in the same energy region, there is a difference not only in the intensity of the signals but also in the number of signals visible.

In Fig. 11, the effect of the MOs reordering, observed in our approach for generating excited states, can be seen. We note that the separation between the core-to-SOMO signals varies significantly for the different states. If relaxation of the orbitals is not considered following photoionization, one would expect this separation to remain constant, as the ordering (and relative energies) of the MOs would here be unchanged. Likewise, the separation between the main signal (transitions from core to virtual orbitals) and the core-to-SOMO transitions would, in this case, be expected to increase from D₁ to D₃. We, however, observe the opposite trend, which underlines the changes observed in the MOs when they relax to the excited state. Hence, the core-to-valence energetics in the spectra are reversed in our approach, compared to an approach relying on Koopmans' theorem due to the hole orbital always appearing as the HOMO in all ionic states. Our observations are, thus, a consequence of the calculation strategy utilizing IMOM to generate the excited states from the neutral ground state by choosing an MO associated with the desired character to be singly occupied. It should be noted that the computations of transitions from core orbitals to virtual orbitals are spin incomplete. This might affect the separation between the core-to-SOMO and the main signal. An experiment might provide information on whether the shown trends provide a physical description or if an alternative approach for computing excited states is required. For now however, we content ourselves with the present calculations, noting that they are consistent with our observations of the MOs exchanging positions in different states (see Fig. 4). While the energy shifts between the core-to-SOMO and core to virtual MOs might not be reliable, we expect the calculation to reflect reliably, which transitions are visible, since the calculation reflects the character of the SOMOs and fully occupied MOs.

Note, however, that, in the D₁-state of the ion, i.e., the hole in the orbital resembling (HOMO-1)₀ (and, thus, the $HOMOD1$—see Fig. 4), the predicted NTOs differ between the methods. At the BHHLYP level of theory, the first bright transition is found to be from 1s(S_T) to an orbital resembling (HOMO-1)₀, as expected. The remaining transitions fall in the same region as the D₀ signal and involve many of the orbitals also observed for that signal. At the fc-CVS-EOM-CCSD level of theory, the dominating NTO pair for the first transition appears to be from 1s(S_T) to the (HOMO-2)₀. Naturally, this orbital should be occupied, and, indeed, upon inspection of the weights of the transitions contributing to the excitation, the dominating transition is found to be from a core orbital into the MO resembling (HOMO-1)₀. Many transitions do, however, contribute significantly to the core excitation, hence to the NTO pair. The appearance of the NTO pair is, therefore, attributed to this fact. The computed spectrum is in good agreement with the TDDFT spectrum, which, as mentioned, shows the types of NTO pairs we expect. Moreover, a similar spectrum is obtained at the fc-CVS-EOM-CCSD level of theory by utilizing the orbitals of the cation (D₀) as the initial guess for the IMOM procedure. Here, the expected core to (HOMO-1)₀ transition is predicted by the NTOs for the first core transition (see Sec. S4.2.2). We, therefore, conclude that the spectral shape can be predicted by employing IMOM, but that the assignment of NTOs must be done with care. When the calculation methods yield the same results, however, it is considered possible to make character assignments. Furthermore, we find that one can use either the MOs of S₀ as the initial guess for the IMOM procedure or the MOs of D₀ and obtain similar results.

In the ground state of the ion, D₀ as well as in the D₂ and D₃ states, the NTOs determined for the main transitions are similar in both the BHHLYP and fc-CVS-EOM-CCSD calculations (see Sec. S4.2). In D₀, the two main transitions are found to originate on different sulfur atoms. As also observed for the neutral species, their separation is small, and, thus, only, with a very good experimental resolution, could these possibly be distinguished. For the D₂ state, we observe that the first two bright transitions are found to be from the 1s orbitals of S_T and S_BT into an orbital resembling the (HOMO-2)₀. Likewise, for D₃, the first bright transition is found to be from 1s(S_T) to an orbital resembling (HOMO-3)₀. This signal, however, has very low intensity and is, therefore, unlikely to be observable in the experiment. The remaining parts of the spectra resemble the D₀ spectrum as also expected.

The XAS calculations are simulated at three different time delays after photoionization determined in terms of D₃ population: 100%, 50%, and 20%. In the MCTDH simulation, this corresponds to time delays of 0.0 fs (100% D₃), 50 fs (52% D₃, 19% D₂, 28% D₁, and 1% D₀), and 150 fs (23% D₃, 30% D₂, 44% D₁, and 3% D₀), while in the TSH simulation, where timescales are expected to better reflect the actual lifetime of states, the corresponding time delays are 0.0 fs (100% D₃), 4.5 fs (52% D₃, 20% D₂, 28% D₁, and 0% D₀), and 18.0 fs (24% D₃, 28% D₂, 48% D₁, and 0% D₀). We recall that the NTOs for the three lowest lying valence transitions of the cation were investigated at the geometries associated with these time delays (see Secs. S2.2 and S4.3.1). They showed that, at these time delays, D_n does, indeed, correspond to a hole in (HOMO-n)₀. Furthermore, to ensure the generation of the correct states in the IMOM procedure, the MOs of the neutral molecule at the geometries associated with the non-zero time delays were determined. For both time delays, they were found to retain the ordering seen for the ground state optimized geometry (see Sec. S4.3 and Fig. S44).

The spectra are computed by weighting the XASs associated with the four doublet states by the state populations of the MCTDH simulation. Furthermore, they are based on the (overly optimistic) assumption of 100% yield of the initial photoionization. Therefore, the total spectra [Fig. 12(A)] are not expected to completely match experiment, as the experimental signal will be a mix of the GS spectrum in Fig. 10 and of the shown spectra. On the other hand, the difference signal [Fig. 12(B)] will not be changed except for an overall scaling. This signal is determined as the S₀ spectrum at the Franck–Condon geometry subtracted from the total spectrum. Since it is the difference spectrum that is of main interest, this procedure is employed here, as the expected yield of the photoionization is not known. The fc-CVS-EOM-CCSD calculations showed poor convergence in the IMOM procedure at displaced geometries. Thus, since we previously found that the XAS are equally well described at the BHHLYP level of theory at a lower computational cost, the tr-XAS is only calculated at this level of theory.

The spectra are evaluated based not only on the structures determined in the MCTDH simulation but also on the ground state optimized structure (i.e., the Franck–Condon geometry). This is done in order to evaluate the effect of small structural changes on the XAS. The structures based on the MCTDH simulation are average structures, which is only a valid approximation as long as the wavepacket does not show significant density in two different areas. This was not the case for the three time delays considered here (see Sec. S3.9) but might yet be the case for other time delays. From the difference signal in Fig. 12(B), it is observed that a decrease in intensity over time occurs in the energy region associated with transitions originating on S_T. Here, a negative difference signal is observed to become more negative over time. An increase in intensity, on the other hand, is observed in the region associated with transitions originating on S_BT. Here, a positive difference signal that decreases from a time delay associated with 100% D₃ to 50% D₃ but increases slightly from 50% D₃ to 20% D₃ is seen. This suggests a decreased electron density around S_T and increased density around S_BT, as also noted in Sec. IV. In addition, increased intensity is observed in the region where the core-to-SOMO transitions occur, as expected based on an increased population in D₁ and D₂, which both display bright core-to-SOMO transitions. Note also that, when computing the D₂ spectrum at geometries of the time delays associated with D₃ populations below 100%, the intensity of the core-to-SOMO transition originating on S_T is observed to decrease. This observation corroborates the conclusion that a charge transfer occurs (see Sec. S4.3). Thus, this study indicates that an electron density transfer from S_T to S_BT might, in principle, be tracked by XAS as also found by Khalili et al.¹⁸ Due to the short timescale of the process, however, this is not expected to be observable in a tr-XAS experiment with the currently available temporal resolution. This is in contrast to the conclusions of Khalili et al.,¹⁸ who reported that the $D1−D0$ decay might be observable in 400 fs. It should, however, be possible to obtain a difference signal corresponding to the time delay associated with 20% D₃.

The changes in geometry at the considered time delays are not seen to affect the total or difference spectra significantly despite the above-mentioned differences of the D₂ spectrum, since the core-to-SOMO transitions have an overall low intensity. In fact, as can be seen in Figs. 11 and 12, as well as Fig. S43, it appears that while structural changes might affect the spectra to some extent, the main changes, particularly when considering the difference spectra, are caused by changes in the electronic structure. Thus, for systems where only small structural changes can occur, one might save computational resources by simply computing excited state XAS based on the ground state optimized geometry, when the evolution of the excited states is the main focus. It should be emphasized that a dynamics simulation is still required to determine the state populations, and hence the tr-XAS. One might, however, compute the excited state spectra only once per state, rather than at each investigated geometry. This approach will not be valid for spectroscopies that are strongly affected by structural changes. Likewise, for molecules more prone to nuclear motion, this approach needs further investigations. Finally, of course, these theoretical findings must be verified experimentally for different systems and processes to ensure the validity of the approach also in practice. The photoexcitation of BT-1T might, for instance, be studied as longer lifetimes are indicated for this process in previous studies.¹⁷ 

A computational study of the tr-XAS of the BT-1T⁺ cation has been presented. It was found that the relevant normal modes for simulating the nuclear dynamics can be identified based on a small set of surface hopping trajectories, as the investigated process to the best of our knowledge does not rely on rare events. A full quantum simulation of the dynamics in reduced space could subsequently be performed based on the determined normal modes. This simulation provided information not only on the evolution of the molecular structure but also on the evolution of the populations of relevant doublet states. We did, however, find that while the evolution of populations could be well described using a reduced space calculation, the lifetime of the process was significantly overestimated by this approach.

Our calculation strategy for simulating XAS resulted in the excited state spectra displaying a trend in the separation between core-to-SOMO peaks and the main band that is counter-intuitive, if the orbitals are considered to be unchanged upon photoionization. The observations were, however, consistent with our observations of orbital changes after photoionization.

When simulating the XAS at different time delays, it was found that the relative populations of the involved states affected the difference signals more than the small structural changes did. Thus, the total tr-XAS might be constructed from a single geometry, when studying processes that rely more on electronic transitions than on nuclear movement. This approximation allows for significant computational savings, which, in turn, makes the approach accessible even when larger systems are considered. It must be emphasized, however, that this approximation should only be considered when investigating simple non-radiative decay processes close to the Franck-Condon point. Despite this simplification, a non-adiabatic nuclear dynamics simulation cannot be avoided as the information it provides on the evolution of the state populations is paramount for the computation of the spectra. As the BT-1T⁺ cation investigated is rather rigid, the approach of considering only one geometry might not be applicable for all molecules, and further investigations are, therefore, needed.

Finally, the computed XAS at different time delays revealed increased intensity in energy regions associated with transitions originating on the S_BT moiety. At the same time, decreased intensities in regions associated with transitions originating on S_T were found. Despite this indication of the electron transfer process from S_T to S_BT, the determined lifetime of the process indicates that it cannot be studied experimentally with the currently available temporal resolutions.

See the supplementary material for additional calculated MOs and NTOs, the initial geometry of the system, and the detailed normal mode analysis as well as all data tables for the plotted spectra. Movies of the eight normal modes considered in the MCTDH simulation are also available.

The research leading to the presented results has received funding from the Independent Research Fund Denmark, Grant Nos. 8021-00347B (KBM) and 7014-00258B (SC), and from the Hungarian National Research, Development and Innovation Fund, Grant No. NKFIH PD 134976 (MP). A.K.S.P. acknowledges the DTU Partnership Ph.D. programme for funding. M.P. acknowledges support from the János Bolyai Scholarship of the Hungarian Academy of Sciences. The authors acknowledge the DTU computational resources,⁷³ on which all calculations were performed.

The authors have no conflicts to disclose.

Anna Kristina Schnack-Petersen: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Funding acquisition (equal); Investigation (lead); Methodology (equal); Project administration (lead); Software (lead); Supervision (lead); Visualization (lead); Writing – original draft (lead). Mátyás Pápai: Conceptualization (supporting); Data curation (lead); Formal analysis (supporting); Investigation (supporting); Methodology (equal); Software (equal); Supervision (supporting); Writing – original draft (equal). Sonia Coriani: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (lead); Project administration (equal); Software (lead); Supervision (equal); Writing – original draft (equal). Klaus Braaagaard Møller: Conceptualization (lead); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (lead); Software (lead); Supervision (lead); Writing – original draft (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.