# Quantum Dynamics of Photoactive Transition Metal Complexes. A Case Study of
Model Reduction

Olga S. Bokareva  Oliver Kühn<EMAIL_ADDRESS>
###### Abstract

Transition metal complexes for photochemical applications often feature a high
density of electron-vibrational states characterized by nonadiabatic and spin-
orbit couplings. Overall, the dynamics after photoexcitation is shaped by
rapid transitions between states of different character and multiplicity. Even
though transient absorption experiments enable characterization in terms of
kinetic rates, the complexity of the systems usually prevents a more detailed
analysis. Quantum dynamics simulations using quantum chemically determined
model Hamiltonians may provide such details. In particular, one is tempted to
pursue a model reduction, such as to identify couplings or vibrational modes
most relevant for the dynamics. Here, we address how such an endeavor is
challenged by the particular nature of transition metal complexes. For that
purpose, we performed quantum dynamics simulations for a recently studied
iron(II) homoleptic complex.

###### keywords:

quantum dynamics , transition metals , linear vibronic coupling , correlation
effects

[1]Leibniz Institute for Catalysis (LIKAT), Albert-Einstein-Str. 29A, 18059
Rostock, Germany [2]Institute of Physics, University of Kassel, Heinrich-
Plett-Str. 40, 34132 Kassel, Germany [3]Institute of Physics, University of
Rostock, Albert-Einstein-Str. 23-24, D-18059 Rostock, Germany

##  1 Introduction

Transition metal (TM) complexes are important for a broad range of
catalytical, medical, biological, and material science applications such as
photocatalysis, molecular memory, cell imaging, etc., see e.g., Refs. [1, 2,
3, 4]. The typical size of many natural and artificial photocatalysts and
their high density of electronic and vibrational states, together with spin-
orbit, relativistic, and environmental effects pose a challenge to quantum
chemical and quantum dynamical simulations [5, 6, 7, 8, 9]. Therefore, most of
the reported simulations of TM excited state dynamics employ the trajectory
surface hopping (TSH) method, see. e.g., Refs. [10, 11, 12, 13, 14, 15]. This
quasi-classical approach enjoys popularity due to its clear concept and the
possibility of studying complex systems with numerous electronic and
vibrational degrees of freedom (DOFs). Of course, one might ask whether all
states and modes are equally important for the dynamics. The answer to this
question often will provide a better mechanistic understanding in addition to
a means for model reduction towards a setup more feasible for quantum
dynamics. Methods for identifying essential nuclear motions are well
established in the field of molecular dynamics [16, 17, 18, 19, 20]. They all
rely solely on the distribution of nuclear geometries, including minimum
energy seams of crossing and minimum energy conical intersections [21].
However, as pointed out by Mai and González [22], displacements alone are not
sufficient to judge on the importance of a particular mode, as the impact of
the amplitude of a given nuclear motion on the wave packet evolution may be
very complex. Instead, information about electronic energy gaps and state
couplings has to be taken into account. A workflow for model reduction
starting from a TSH simulation has been proposed by Gómez et al. [12], which
additionally included feedback on the reduced dimensional quantum simulation,
such as to improve TSH parameters. Trajectory surface hopping was also used by
other authors to identify relevant modes for quantum dynamics simulations of
TM compounds [13].

As excited dynamics takes place on the potential energy surfaces (PESs)
precomputed with the chosen quantum-chemical approach or on-the-fly, the
outcomes of such computations may strongly vary with the applied computational
scheme. An exemplary homoleptic iron(II) complex $\mathbf{FePS}$ with two
push-pull ligands [23] has been recently studied by us using the combination
of TSH excited states dynamics with optimally-tuned long-range separated
functionals [24]. Although the concept of optimal tuning was introduced to
reliably describe charge-transfer properties such as metal-to-ligand charge
transfer (MLCT) excitations, the appearance of multiple low-lying metal-
centered states makes the choice of an appropriate range-separation scheme
more complicated. Being formally local, those MC states may suffer from
uncompensated self-interaction error and, thus, require special attention
beyond the classical optimal tuning. Following this, additional criteria may
be desired to adjust the computational scheme. Assisted by comparison with the
experimental absorption spectra and CASPT2 calculations, two suitable sets of
range-separation parameters were chosen for excited state surface hopping
dynamics simulations. The suggested dynamics mechanisms were notably different
for both sets of parameters. Using the set of optimal parameters from the
classical $\Delta$SCF method (set A) resulted in the population of ^3MLCT
states that remained stable within the 2 ps simulation time, which is at
variance with transient absorption experiments [23]. Applying multi-reference
CASPT2-derived parameter set, TSH calculations led to relaxation processes
into hot singlets states. In the course of the deactivation, the system stayed
in regions of high density of states, reachable for further transfer into
quintet states in accordance with experimental findings [23].

In what follows, we will use the CASPT2-based optimized functional to
illustrate some aspects of quantum dynamics simulations of TM compounds alike
$\mathbf{FePS}$. It should be stressed that the purpose is not to obtain
converged results for the total system, although, in principle, this should be
feasible. Instead, we will highlight some issue which one may encounter when
trying to design a reduced model for a system with a large density of states
and a large number of couplings of comparable magnitude. In what follows, we
first outline the linear vibronic coupling (LVC) model used for the quantum
dynamics simulations. Next, we address the Multi-layer Multiconfiguration
Time-dependent Hartree (ML-MCTDH) method. In the applications section, we
first characterize the LVC model for the $\mathbf{FePS}$ system. Possible
model reduction is investigated using the Time-dependent Hartree (TDH) limit
of MCTDH. For the thus identified reduced model, ML-MCTDH simulations are
performed to address correlation effects. Finally, conclusions are given.

##  2 Theory

###  2.1 Linear Vibronic Coupling

In contrast to trajectory-based methods, quantum dynamics simulations require
a priori knowledge of the PES, including possible electronic state couplings.
For TM complexes, this inevitably means that one has to resort to model PES,
the most prominent one being the vibronic coupling model [25, 26, 27]. It
rests on a harmonic oscillator description of nuclear motion in terms of
normal modes, combined with a low-order expansion for the displacement of the
PESs and the state couplings. Specifically, the linear LVC Hamiltonian reads,
denoting the set of $f$ normal mode coordinates as $\mathbf{Q}$,

$H=\sum_{mn}[\delta_{mn}H_{m}(\mathbf{Q})+V_{mn}(\mathbf{Q})]|m\rangle\langle
n|$ (1)

with the Hamiltonian for electronic state $m$

$H_{m}(\mathbf{Q})=E_{m}+\frac{1}{2}\sum_{\xi}\hbar\omega_{\xi}(-\partial^{2}_{\xi}+Q_{\xi}^{2})+\sum_{\xi}\kappa_{m,\xi}Q_{\xi}$
(2)

and the interstate couplings

$V_{mn}(\mathbf{Q})=V_{mn}^{(0)}+\sum_{\xi}\lambda_{mn,\xi}Q_{\xi}\,.$ (3)

Here, $\kappa_{m,\xi}$ and $\lambda_{mn,\xi}$ are the intra- and interstate
coupling constants, respectively, related to mode $Q_{\xi}$. Further,
$V_{mn}^{(0)}=V_{mn}^{(\rm SOC)}$ denotes the spin-orbital coupling (SOC),
assuming that its coordinate dependence is negligible. The Hamiltonian, Eqs.
(1 \- 3), can be supplemented by

$H_{\rm f}(t)=-{\mathbf{E}}(t)\sum_{mn}{\mathbf{d}}_{mn}|m\rangle\langle n|$
(4)

describing the interaction with an external time-dependent electric field
${\mathbf{E}}(t)$ in dipole approximation (${\mathbf{d}}_{mn}$ is the vector
of the dipole matrix elements).

Note that the LVC Hamiltonian is not only used for quantum dynamics
simulations (see, e.g. Refs. [28, 29]), if applicable, it also provides an
ideal framework for TSH. Fortunately, this often holds true for TM complexes
where on-the-fly generation of PES is computationally very expensive [14, 24].
This model is only valid if no large amplitude motions in the excited state
dynamics are present, for diagnostics of accuracy, see the exemplary study in
Ref. [30]. In cases where bond breaking is simulated, the LVC model for MCTDH
may be complemented by the corresponding dissociative mode, see e.g. studies
of CO-release from TM-complexes [31, 32].

###  2.2 ML-MCTDH Method

The method of choice for high-dimensional quantum dynamics simulations of LVC
Hamiltonians is the MCTDH approach [33, 34, 35, 36, 37], which is briefly
sketched in what follows. Given the expansion of the state vector in the
diabatic electronic basis

$|\Psi({\bf Q};t)\rangle=\sum_{m}\psi_{m}({\bf Q};t)\,|m\rangle\,.$ (5)

The state-specific vibrational wave packets are represented as follows
(skipping the state index $m$ for simplicity)

$\psi({\bf
Q};t)=\sum_{j_{1}=1}^{n_{1}}\ldots\sum_{j_{f}=1}^{n_{f}}A_{j_{1}\ldots
j_{f}}^{1}(t)\cdot\phi_{j_{1}}^{1;1}(Q_{1},t)\ldots\phi_{j_{f}}^{1;f}(Q_{f},t)\,.\,.$
(6)

Here, $\phi_{\lambda}^{1,\kappa}(Q_{\kappa},t)$ are the so-called single-
particle functions (SPFs) which form a time-dependent basis for expanding the
wave packet in terms of different Hartree products of SPFs. Eq. (6) provides a
compact representation such that the number of SPFs per DOF, $n_{\kappa}$,
usually will be well below that required for a static basis. This ansatz can,
in principle, be converged towards the numerically exact limit. On the other
hand side, there is the TDH limit, where only a single SPF per DOF is taken
into account. In practice, the SPFs are expanded into a time-independent
(primitive) basis set, $\chi_{j}^{(\kappa)}(Q_{\kappa})$, according to

$\phi^{1;\kappa}_{\lambda}(Q_{\kappa},t)=\sum^{N_{\kappa}}_{j=1}A^{2;\kappa}_{\lambda;j}(t)\chi_{j}^{(\kappa)}(Q_{\kappa})\,.$
(7)

MCTDH can be seen as a two-layer approach. This view becomes clear upon
interpretation of $A^{2;\kappa}_{\lambda;j}(t)$ as an additional set of
expansion coefficients to represent the time-dependent SPFs of the upper
(first) layer in a time-independent primitive basis in the lower (second)
layer. In $A^{2;\kappa}_{\lambda;j}(t)$ the superscripts 2 and $\kappa$
correspond to the second layer and the $\kappa$th mode, respectively, whereas
the subscript $\lambda$ refers to the indices of the respective SPF of the
first layer. The expansion coefficients in Eqs. (6) and (7) can be obtained
using the time-dependent Dirac-Frenkel variational principle [34].

ML-MCTDH builds on the idea of combining strongly correlated DOFs into $p$
logical coordinates (combined modes or MCTDH particles) yielding the set
$\\{q^{1}_{1},\ldots,q^{1}_{p}\\}$ with
$q^{1}_{\kappa}=\\{Q_{1_{\kappa}},\ldots,Q_{d_{\kappa}}\\}$, where the
superscript refers to the first layer, and $d_{\kappa}$ is the dimension of
the $\kappa$th particle. The wave packet within the combined modes’ picture
for the uppermost (first) layer can be written as

$\displaystyle\Psi(q^{1}_{1},\ldots,q^{1}_{p},t)=\sum_{j_{1}=1}^{n_{1}}\ldots\sum_{j_{p}=1}^{n_{p}}A_{j_{1}\ldots
j_{p}}^{1}(t)\,\phi_{j_{1}}^{1;1}(q^{1}_{1},t)\ldots\phi_{j_{p}}^{1;p}(q^{1}_{p},t)\,.$
(8)

In ML-MCTDH, the multi-dimensional SPFs for the particles,
$\phi_{\lambda}^{1;\kappa}(q^{1}_{\kappa},t)$, are considered as multi-
dimensional wave packets of the second layer and expanded according to

$\displaystyle\phi_{\lambda}^{1;\kappa}(q^{1}_{\kappa},t)$ $\displaystyle=$
$\displaystyle\phi_{\lambda}^{2;\kappa}(Q_{1}^{2;\kappa},\ldots,Q_{d_{\kappa}}^{2;\kappa},t)$
(9) $\displaystyle=$
$\displaystyle\sum\limits_{j_{1}}^{n_{\kappa,1}}\ldots\sum\limits_{j_{d_{\kappa}}}^{n_{\kappa,d_{\kappa}}}A^{2;\kappa}_{\lambda;j_{1}\ldots
j_{d_{\kappa}}}(t)\,\phi_{j_{1}}^{2;\kappa;1}(Q_{1}^{2;\kappa},t)\ldots\phi_{j_{d_{\kappa}}}^{2;\kappa;d_{\kappa}}(Q_{d_{\kappa}}^{2;\kappa},t)\,.$

In the third layer, these wave packets are represented in the primitive basis,
i.e.

$\displaystyle\phi_{\lambda}^{2;\kappa;\sigma}(Q_{\sigma}^{2;\kappa},t)=\sum\limits_{j=1}^{N_{\alpha}}A^{3;\kappa;\sigma}_{\lambda;j}(t)\,\chi_{j}^{(\alpha)}(Q_{\sigma}^{2;\kappa}),\quad\quad\alpha=\sigma+\sum\limits_{i}^{\kappa-1}d_{i}\,.$
(10)

While this example included three layers, the general idea can easily be
applied to a case with more layers containing particles with decreasing
dimensionality when stepping towards the bottom layer. Within this scheme,
electronic states are incorporated by introducing an electronic coordinate,
$Q_{\rm el}$, in the uppermost layer. By construction, ML-MCTDH allows for a
numerical solution of the time-dependent Schrödinger equation with controlled
convergence. Finally we note that the effect of an external laser field, Eq.
(4), can be incorporated into the MCTDH scheme a straightforward manner [38].
For the simulations present below the Heidelberg program package was used
[39].

##  3 Application

###  3.1 LVC Hamiltonian for $\mathbf{FePS}$

In what follows, we will apply MCTDH to study the quantum dynamics of an
exemplary iron(II) homoleptic complex [Fe(cpmp)2]2+, where cpmp =
6,2”-carboxypyridyl-2,2’-methylamine-pyridyl-pyridine, shown as an inset in
Figure 3 and abbreviated as $\mathbf{FePS}$ [23]. The choice is motivated by
our previous investigation of the influence of the parametrization of the
range-separated hybrid functional on the photoinduced dynamics as described by
TSH. However, the present focus will not be on a quantitative comparison with
the TSH dynamics reported in Ref. [24], but on coping with the high
dimensionality of the LVC Hamiltonian for $\mathbf{FePS}$ in general. As an
example, we use the parameter set of the functional derived based on a
comparison with CASPT2 results, which gave a better agreement between
calculated and experimental dynamics. In numbers, taking into account the
lowest 9 singlet and 14 triplet states (energies up to about 2.8 eV), there
are in total 51 electronic states. Further, the molecule has 212 vibrational
DOFs.

Let us start with an analysis of the LVC Hamiltonian reported in Ref. [24].
Assuming an energy threshold of 0.0003 eV, the LVC Hamiltonian contains 4722
$\kappa_{m,\xi}$ terms, 20119 $\lambda_{mn,\xi}$ terms, and 2464 SOC terms.
Figure 1 shows the distribution of LVC parameters, which indicates that the
majority of couplings is rather small. In Figure 2a, we have plotted the
distribution of zeroth-order states, defined as spin-free states without
nonadiabatic state coupling ($\lambda_{mn,\xi}=0$). Apparently, the number of
states increases rather rapidly with energy.

Figure 1: Distribution of LVC and SOC coupling parameters: (a)
$\kappa=|\kappa_{m,\xi}|$, (b) $\lambda=|\lambda_{mn,\xi}|$, and (c) $|{\rm
SOC}|=|V_{mn}^{(\rm SOC)}|$. Figure 2: Analysis of the LVC model for
$\mathbf{FePS}$ [24]. (a) Distribution of diabatic electron-vibrational spin-
free singlet and triplet states (bin width 5 meV). (b) Distribution of state
crossing for all possible transitions between spin-free singlet and triplet
states (bin width 0.5).

Since it is known that TM compounds often feature a rather rapid ISC dynamics,
it is interesting to inspect the crossing points between singlet and triplet
states. According to the LVC Hamiltonian, for a pair of zeroth-order states,
they are given by $q_{\xi}^{*}=(E_{n}-E_{m})/(\kappa_{m,\xi}-\kappa_{n,\xi})$.
The distribution of crossing points for $\mathbf{FePS}$ is shown in Figure 2b.
To put these numbers into perspective, one should note that the mode
elongation at the zero point energy ranges between $\pm$1 (in dimensionless
units). The actual values for $\mathbf{FePS}$ are mostly much larger. The
reason is to be found in the rather small displacements, cf. Figure 1. This
observation is no contradiction to a rapid ISC since the potential curves run
almost parallel until they eventually cross. Since SOC is assumed to be not
coordinate dependent, this facilitates an efficient ISC if energies of spin-
free states are energetically not too far apart.

###  3.2 TDH Dynamics

Obtaining converged MCTDH dynamics for the given model on a few picosecond
time scale should be possible with the multilayer approach but is not the
focus of this contribution. Instead, we will use the TDH limit of the MCTDH
ansatz, which allows for full-dimensional quantum simulations [40]. The TDH
limit facilitates the exploration of overall aspects of the dynamics of
different reduced models on equal footing. Note that for ML-MCTDH
calculations, each model should be converged separately. The analysis in the
previous section showed that the LVC model is characterized by a high density
of states and a large number of couplings terms having a comparable magnitude.
In this section, we will address the question of whether this situation allows
for a model reduction by defining certain cut-offs for the coupling
parameters.

Figure 3: Upper panel: Absorption spectrum of $\mathbf{FePS}$ as calculated by
the TDH approach to the dipole autocorrelation function as well as
experimental results [23]. Also shown is the laser pulse spectrum (gray shaded
area). Lower panels: Decomposition of electronic transitions based on density
matrix analysis [41].

Figure 3 compares the experimental absorption spectrum with the one calculated
with the correlation function approach using TDH dynamics of the full model.
Further, an assignment of spin-free electronic transitions is given. The
overall lineshape of the spectrum is reasonably reproduced; in addition to the
inherent broadening due to the large density of states, we have added a
phenomenological broadening of $\gamma^{-1}=10$ fs to account for solvation
effects. As compared with experiments, there is a shift of about 0.3 eV. Since
this range of the experimental spectrum does not convey many details, model
reduction will be investigated for laser-driven dynamics.

For the following simulations, we will consider the $x$-component of the
dipole operator only (which is of higher magnitude than the other components).
Further we assume a Gaussian pulse shape, i.e.

$E(t)=E_{0}\cos(\Omega t)\exp\\{-(t-t_{0})^{2}/2\sigma^{2}\\}\,.$ (11)

Here, $\Omega$, $t_{0}$, and $\sigma$ are the carrier frequency, the pulse
center, and the pulse width, respectively. The transition energy at the peak
of absorption due to the $x$-component of the dipole is at 2.27 eV. For the
width parameter, we use 24 fs yielding a full-width a half maximum of the
pulse intensity of 40 fs. The field amplitude has been chosen as 0.001
Hartree/Bohr, which gave a population transfer of about 25-30% to the excited
singlet states for the TDH dynamics.

Figure 4: TDH dynamics of summed excited singlet state populations after
excitation with a $x$-polarized Gaussian laser pulse centered at $t_{0}=100$
fs ($\hbar\Omega=2.27$ eV, $\sigma=24$ fs, $E_{0}=0.001$ Hartree/Bohr). The
different curves correspond to the full model (denoted ’ref’) and various
reduced modes as explained in the text.

Figure 4 compares the population of all singlet states (except the ground
state) for different models. The reference case (full model) is characterized
by an initial rise due to the laser pulse and an immediate rapid decay. (Note
that this decay is not due to stimulated emission.) Subsequently, one observes
a slower decay with a partial revival around 600-700 fs as well as
superimposed rapid oscillations. The dynamics of the other three models shown
in the figure essentially agree during the first 200 fs. Putting the threshold
for the nonadiabatic couplings to $|\lambda_{mn,\xi}|=1$ meV (I, blue line),
the decay after 200 fs is a bit slower, and the partial revival doesn’t occur.
Using, in addition, the threshold $|\kappa_{m,\xi}|=5$ meV (II, orange line)
brings back the revival. Putting in addition a threshold for SOC of
$|V_{mn}^{(\rm SOC)}|=0.5$ meV (III, red line), there is again no pronounced
revival. It should be stressed that these parameters have been taken after
consideration of certain larger cut-off values such as $|\lambda_{mn,\xi}|=5$
meV, which gave a pronounced difference in the dynamics. (Note that for the
purpose of illustration, it was not attempted to pinpoint these parameters in
more detail.) At first glance, the small magnitude of these cut-off parameters
is a bit surprising. In fact, inspecting Figure 1 one could have been tempted
to select just a few modes having with the largest couplings. Nevertheless, in
the present case of $\mathbf{FePS}$, the dynamics after laser excitation
appears to be shaped by the combined effect of many small-magnitude couplings.
We note that the absorption spectrum is essentially indistinguishable for
these four models. The resulting LVC Hamiltonian with all three coupling
parameters reduced contains 2941 $\kappa_{m,\xi}$ terms, 13243
$\lambda_{mn,\xi}$ terms, and 2166 SOC terms.

Figure 5: Comparison of the population dynamics for model III of Fig. 4 (212
modes) with various models having a reduced number of vibrational modes as
indicated; for selection details see text. The inset shows the absorption
spectrum for the four modes according to the $x$-component of the dipole.

Next, we explore a possible reduction of the number of vibrational DOFs.
Intuitively one would be tempted to compare the coordinate expectation values
with the values of the crossing points $q_{\xi}^{*}$ between singlet and
triplet states. However, according to Figure 2b, the majority of $q_{\xi}^{*}$
takes rather large values which are not reached during the dynamics. As a
technical note, we would like to mention that for the simulation, the single-
set method has been used where a given SPF is optimized for all electronic
states and coordinate expectation values are calculated accordingly. Figure 5
compares the reduced model identified above with models having different cut-
offs for the average normal mode elongation. Namely, the cut-off has been set
to 0.5 (24 modes), 0.2 (70 modes), and 0.1 (131 modes). First, we note that in
contrast to Figure 4, the initial dynamics during the first 200 fs is
influenced by the cut-off. In fact, inspecting the absorption spectrum shown
in the inset, one notices marked changes upon reduction of the number of
modes. The change in cut-off modifies the overlap with the pulse spectrum and,
thus, the wave packet created on the excited state manifold. To have
reasonably accurate initial excitation and, at the same time, a considerable
reduction of the number of modes, we identified the model with 70 modes as a
compromise. Again, we emphasize that the thresholds are rather small and not
that much different, i.e. it is the net effect of a large number of small-
amplitude vibrations, which is responsible for the observed dynamics.

###  3.3 ML-MCTDH Dynamics

TDH is the mean-field approximation to MCTDH and, as such, neglects
correlations in the dynamics of the different DOFs. One could argue that in
cases of a large number of small couplings combined with a large density of
states, the mean-field approximation should work well. However, since
$\mathbf{FePS}$ already gave some surprises, it will be better to compare TDH
with MCTDH dynamics, in particular, since ML-MCTDH provides a means to obtain
numerically converged results. For the purpose of illustration, we will use
the final 70 modes model obtained from the TDH simulations in the previous
section. ML-MCTDH requires to identify a ML-tree to structure the wave packet.
This procedure is a nontrivial task, and different trees might lead to rather
different CPU usage for obtaining the same accuracy [42]. Presently, we did
not attempt to optimize the tree structure. Instead, we started with two
branches, one including all modes with an average elongation beyond 0.5 (21
modes) and one including the remaining 49 modes. The first branch was further
subdivided into three branches, each containing two two-dimensional and one
three-dimensional particle. The second branch was first subdivided into four
subbranches, each of them containing three further branches. On the lowest
level, there were 23 two-dimensional and one three-dimensional particle. The
number of SPFs in each layer has been chosen such that the largest population
of the least occupied natural orbital was smaller than 2%.

Figure 6: Comparison between TDH and ML-MCTDH dynamics after laser pulse
excitation for the reduced model introduced in Section 3.2. (a) Total
population of excited singlet states. (b) Contributions of individual types of
states to the ML-MCTDH dynamics as indicated.

Figure 6a compares TDH and ML-MCTDH dynamics for the total excited singlet
state population. Even though the linear absorption spectrum is essentially
identical, the difference in population dynamics is dramatic. First, ML-MCTDH
does not show the pronounced initial peak during and immediately after the
laser pulse as is observed for TDH. More notably, however, the decay of the
total singlet state population is much slower in the ML-MCTDH case. In fact,
this brings the present simulation in closer agreement with the TSH results
reported in Ref. [24]. Of course, one might argue that at this point, one
should go back to the model reduction performed in Section 3.2 and to repeat
it using ML-MCTDH. However, we don’t expect a qualitative difference, i.e.,
the illustration of the challenges to model reduction performed with TDH will
still be valid.

In Figure 6b, ML-MCTDH populations of individual sets of states are shown. The
rapid initial population of triple states is typical for TM compounds and
essentially occurs during the laser pulse action. Here, we notice that the
most/least population is in the 1MLCT/1MC states. The population of 3MLCT/3MC
states is about equal and in between the 1MLCT and 1MC cases. This general
behavior is also in accord with the TSH results of Ref. [24]. Note that one
would not expect a quantitative agreement, not at least because the
description of the laser pulse excitation differs in the two simulations (cf.
systematic study of this aspect in Ref. [43]).

##  4 Conclusions

We have studied the laser-driven dynamics of the exemplary FePS complex whose
photophysical and quantum chemical properties were previously reported along
with TSH simulations. FePS shows the high density of electron-vibrational
states typical for such TM complexes. Using an LVC Hamiltonian, the reduction
of the number of electronic coupling parameters was explored, followed by a
reduction of the number of vibrational modes. The structureless linear
absorption spectrum was not useful in this respect, and, therefore, the TDH
dynamics was analyzed. It was found that the dynamics is influenced by a large
number of small couplings, which prevents the identification of a
representative model with just a few DOFs. For the thus established model,
still containing 70 vibrational modes and thousands of coupling terms, the
effect of correlations on the dynamics was investigated. Here one would have
expected that given the nature of the model, effects beyond the mean-field
description play a minor role only. Surprisingly, only the inclusion of
correlations brings the dynamics in good agreement with previously reported
TSH simulations and experiments.

We expect that the challenges for model reduction discussed here for an
exemplary case are typical for similar TM complexes.

## Acknowledgement

The authors are grateful to Dr. J.-P. Zobel for providing the LVC Hamiltonian
and L. Kleindienst for his support with preparing MCTDH input files from LVC
data. This work has been financially supported by the Deutsche
Forschungsgemeinschaft [Priority Program SPP 2102 “Light controlled reactivity
of metal complexes” (grant KU952/12-1)].

## References

  * [1] A. J. Esswein, D. G. Nocera, Hydrogen Production by Molecular Photocatalysis, Chem. Rev. 107 (2007) 4022–4047.
  * [2] R. E. Blankenship, D. M. Tiede, J. Barber, G. W. Brudvig, G. Fleming, M. Ghirardi, M. R. Gunner, W. Junge, D. M. Kramer, A. Melis, T. A. Moore, C. C. Moser, D. G. Nocera, A. J. Nozik, D. R. Ort, W. W. Parson, R. C. Prince, R. T. Sayre, Comparing Photosynthetic and Photovoltaic Efficiencies and Recognizing the Potential for Improvement, Science 332 (2011) 805–809.
  * [3] J. Jiao, J. He, M. Li, J. Yang, H. Yang, X. Wang, S. Yang, A porphyrin-based metallacage for enhanced photodynamic therapy, Nanoscale 14 (2022).
  * [4] M. A. Halcrow (Ed.), Spin-Crossover Materials: Properties and Applications, Wiley, Chichester, West Sussex, United Kingdom, 2013.
  * [5] L. González, D. Escudero, L. Serrano-Andrés, Progress and Challenges in the Calculation of Electronic Excited States, ChemPhysChem 13 (2012) 28–51.
  * [6] T. J. Penfold, E. Gindensperger, C. Daniel, C. M. Marian, Spin-Vibronic Mechanism for Intersystem Crossing, Chemical Reviews 118 (2018) 6975–7025.
  * [7] C. Daniel, Ultrafast processes: Coordination chemistry and quantum theory, Phys. Chem. Chem. Phys. 23 (2020) 43–58.
  * [8] P. S. Wagenknecht, P. C. Ford, Metal centered ligand field excited states: Their roles in the design and performance of transition metal based photochemical molecular devices, Coord. Chem. Rev. 255 (2011) 591–616.
  * [9] J. P. Zobel, L. González, The Quest to Simulate Excited-State Dynamics of Transition Metal Complexes, JACS Au 1 (2021) 1116–1140.
  * [10] J. C. Tully, Molecular dynamics with electronic transitions, The Journal of Chemical Physics 93 (1990) 1061–1071.
  * [11] S. Mai, P. Marquetand, L. González, Nonadiabatic dynamics: The SHARC approach, WIREs Comput. Mol. Sci. 8 (2018) e1370.
  * [12] S. Gómez, M. Heindl, A. Szabadi, L. González, From Surface Hopping to Quantum Dynamics and Back. Finding Essential Electronic and Nuclear Degrees of Freedom and Optimal Surface Hopping Parameters, J. Phys. Chem. A 123 (2019) 8321–8332.
  * [13] M. Pápai, Toward Simulation of Fe(II) Low-Spin → High-Spin Photoswitching by Synergistic Spin-Vibronic Dynamics, J. Chem. Theory Comput. 18 (2022) 1329–1339.
  * [14] J. P. Zobel, M. Heindl, F. Plasser, S. Mai, L. González, Surface Hopping Dynamics on Vibronic Coupling Models, Acc. Chem. Res. 54 (2021) 3760–3771.
  * [15] F. Plasser, S. Mai, M. Fumanal, E. Gindensperger, C. Daniel, L. González, Strong Influence of Decoherence Corrections and Momentum Rescaling in Surface Hopping Dynamics of Transition Metal Complexes, J. Chem. Theor. Comp. 15 (2019) 5031–5045.
  * [16] A. M. Virshup, J. Chen, T. J. Martínez, Nonlinear dimensionality reduction for nonadiabatic dynamics: The influence of conical intersection topography on population transfer rates, J. Chem. Phys. 137 (2012) 22A519.
  * [17] J. P. P. Zauleck, S. Thallmair, M. Loipersberger, R. de Vivie-Riedle, Two New Methods To Generate Internal Coordinates for Molecular Wave Packet Dynamics in Reduced Dimensions, J. Chem. Theor. Comp. 12 (2016) 5698–5708.
  * [18] G. Capano, T. J. Penfold, M. Chergui, I. Tavernelli, Photophysics of a copper phenanthroline elucidated by trajectory and wavepacket-based quantum dynamics: A synergetic approach, Phys. Chem. Chem. Phys. 19 (2017) 19590–19600.
  * [19] X. Li, Y. Xie, D. Hu, Z. Lan, Analysis of the Geometrical Evolution in On-the-Fly Surface-Hopping Nonadiabatic Dynamics with Machine Learning Dimensionality Reduction Approaches: Classical Multidimensional Scaling and Isometric Feature Mapping, J. Chem. Theor. Comp. 13 (2017) 4611–4623.
  * [20] S. R. Hare, L. A. Bratholm, D. R. Glowacki, B. K. Carpenter, Low dimensional representations along intrinsic reaction coordinates and molecular dynamics trajectories using interatomic distance matrices, Chem. Sci. 10 (2019) 9954–9968.
  * [21] Y. Harabuchi, J. Eng, E. Gindensperger, T. Taketsugu, S. Maeda, C. Daniel, Exploring the Mechanism of Ultrafast Intersystem Crossing in Rhenium(I) Carbonyl Bipyridine Halide Complexes: Key Vibrational Modes and Spin–Vibronic Quantum Dynamics, J. Chem. Theor. Comp. 12 (2016) 2335–2345.
  * [22] S. Mai, L. González, Identification of important normal modes in nonadiabatic dynamics simulations by coherence, correlation, and frequency analyses, J. Chem. Phys. 151 (2019) 244115.
  * [23] J. Moll, R. Naumann, L. Sorge, C. Förster, N. Gessner, L. Burkhardt, N. Ugur, W. Nürnberger, P. Seidel, C. Ramanan, M. Bauer, H. K., Pseudo-ocathedral iron(II) complexes with near-degenerate charge transfer and ligand field states at the Franck-Condon geometry, Chem. Eur. J. 28 (2022) e2022ß1858.
  * [24] J. P. Zobel, A. Kruse, O. Baig, S. Lochbrunner, S. I. Bokarev, O. Kühn, L. González, O. S. Bokareva, Can range-separated functionals be optimally tuned to predict spectra and excited state dynamics in photoactive iron complexes?, Chem. Sci. 14 (2023) 1491–1502.
  * [25] H. Köppel, W. Domcke, L. S. Cederbaum, Multimode molecular dynamics beyond the Born-Oppenheimer approximation, Adv Chem Phys 57 (1984) 59.
  * [26] T. J. Penfold, E. Gindensperger, C. Daniel, C. M. Marian, Spin-Vibronic Mechanism for Intersystem Crossing, Chem. Rev. 118 (2018) 6975–7025.
  * [27] V. May, O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems, Fourth Edition, Wiley-VCH, Weinheim, 2023.
  * [28] M. Fumanal, E. Gindensperger, C. Daniel, Ultrafast Intersystem Crossing vs Internal Conversion in $\alpha$-Diimine Transition Metal Complexes: Quantum Evidence, J. Phys. Chem. Lett. 9 (2018) 5189–5195.
  * [29] M. Pápai, T. Rozgonyi, T. J. Penfold, M. M. Nielsen, K. B. Møller, Simulation of ultrafast excited-state dynamics and elastic x-ray scattering by quantum wavepacket dynamics, J. Chem. Phys. 151 (2019) 104307.
  * [30] T. J. Penfold, J. Eng, Mind the GAP: Quantifying the breakdown of the linear vibronic coupling Hamiltonian, Physical Chemistry Chemical Physics 25 (2023) 7195–7204.
  * [31] K. Falahati, H. Tamura, I. Burghardt, M. Huix-Rotllant, Ultrafast carbon monoxide photolysis and heme spin-crossover in myoglobin via nonadiabatic quantum dynamics, Nat. Comm. 9 (2018) 4502.
  * [32] M. Fumanal, C. Daniel, E. Gindensperger, Excited-state dynamics of [Mn(im)(CO) 3 (phen)] + : PhotoCORM, catalyst, luminescent probe?, J. Chem. Phys. 154 (2021) 154102.
  * [33] H. D. Meyer, U. Manthe, L. S. Cederbaum, The multi-configurational time-dependent Hartree approach, Chem. Phys. Lett. 165 (1990) 73–78.
  * [34] M. H. Beck, A. Jäckle, G. A. Worth, H.-D. Meyer, The multiconfiguration time-dependent Hartree method: A highly efficient algorithm for propagating wavepackets, Phys. Rep. 324 (2000) 1–105.
  * [35] H. Wang, M. Thoss, Multilayer formulation of the multiconfiguration time-dependent Hartree theory, J. Chem. Phys. 119 (2003) 1289–1299.
  * [36] H.-D. Meyer, Studying molecular quantum dynamics with the multiconfiguration time-dependent Hartree method, WIREs Comput. Mol. Sci. 2 (2011) 351–374.
  * [37] O. Vendrell, H.-D. Meyer, Multilayer multiconfiguration time-dependent Hartree method: Implementation and applications to a Henon-Heiles Hamiltonian and to pyrazine, J. Chem. Phys. 134 (2011) 044135.
  * [38] H. Naundorf, G. A. Worth, H.-D. Meyer, O. Kühn, Multiconfiguration Time-Dependent Hartree Dynamics on an ab Initio Reaction Surface: Ultrafast Laser-Driven Proton Motion in Phthalic Acid Monomethylester, J. Phys. Chem. A 106 (2002) 719–724.
  * [39] G. A. Worth, M. H. Beck, A. Jäckle, H.-D. Meyer, The MCTDH Package, Version 8.2, (2000). H.-D. Meyer, Version 8.3 (2002), Version 8.4 (2007). O. Vendrell and H.-D. Meyer, Version 8.5 (2013), used Version 8.6.3 (2023), See http://mctdh.uni-hd.de/.
  * [40] GK. Paramonov, H. Naundorf, O. Kühn, Ultrafast multidimensional dynamics of strong hydrogen bonds - A Cartesian reaction surface approach, Eur. Phys. J. D 14 (2001) 205–215.
  * [41] F. Plasser, TheoDORE: A toolbox for a detailed and automated analysis of electronic excited state computations, J. Chem. Phys. 152 (2020) 084108.
  * [42] M. Schröter, S. D. Ivanov, J. Schulze, S. P. Polyutov, Y. Yan, T. Pullerits, O. Kühn, Exciton-vibrational coupling in the dynamics and spectroscopy of Frenkel excitons in molecular aggregates, Phys. Rep. 567 (2015) 1–78.
  * [43] M. Heindl, L. González, Validating fewest-switches surface hopping in the presence of laser fields, The Journal of Chemical Physics 154 (2021) 144102.