# Extending Non-Perturbative Simulation Techniques for Open-Quantum Systems to Excited-State Proton Transfer and Ultrafast Non-Adiabatic Dynamics

Brieuc Le Dé, Simon Huppert, Riccardo Spezia, Alex W. Chin

► To cite this version:

Brieuc Le Dé, Simon Huppert, Riccardo Spezia, Alex W. Chin. Extending Non-Perturbative Simulation Techniques for Open-Quantum Systems to Excited-State Proton Transfer and Ultrafast Non-Adiabatic Dynamics. Journal of Chemical Theory and Computation, 2024, 20 (20), pp.8749-8766. 10.1021/acs.jctc.4c00666 . hal-04796864

HAL Id: hal-04796864

https://hal.sorbonne-universite.fr/hal-04796864v1

Submitted on 22 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Extending Non-Perturbative Simulation
Techniques for Open-Quantum Systems to
Excited-State Proton Transfer and Ultrafast
Non-Adiabatic Dynamics

Brieuc Le Dé,† Simon Huppert,† Riccardo Spezia,‡ and Alex W. Chin*,†

†Sorbonne Université, CNRS, Institut des NanoSciences de Paris, 4 place Jussieu, 75005 Paris, France

‡Sorbonne Université, CNRS, Laboratoire de Chimie Théorique, 4 place Jussieu, 75005 Paris, France

E-mail: alex.chin@insp.upmc.fr

Abstract

Excited state proton transfer is an ubiquitous phenomenon in biology and chem-
istry, spanning from the ultrafast reactions of photo-bases and acids to light-driven,
enzymatic catalysis and photosynthesis. However, the simulation of such dynamics
involves multiple challenges, since high-dimensional, out-of-equilibrium vibronic states
play a crucial role, while a fully quantum description of the proton's dissipative, real-
space dynamics is also required. In this work, we extend the powerful Matrix Product
State approach to open quantum systems (TEDOPA) to study these demanding dy-
namics, and also more general non-adiabatic processes that can appear in complex
photochemistry subject to strong laser driving. As an illustration, we initially consider
an open model of a four-level electronic system interacting with hundreds of intramolecular vibrations that drive ultrafast excited state proton transfer, as well as an explicit photonic environment that allows us to directly monitor the resulting dual fluorescence in this system. We then demonstrate how to include a continuous 'reaction coordinate' of the proton transfer that allows numerically exact simulations that can be understood, visualized and interpreted in the familiar language of diabatic and adiabatic dynamics on potential surfaces, while also retaining an exact quantum treatment of dissipation and driving effects that could be used to study diverse problems in ultrafast photochemistry.

# 1 Introduction

Proton transfer is an elementary physical process that underlies a vast number of chemical transformations which play critical roles in biology,¹,² catalysis, and environmental energy harvesting.³,⁴ Excited-State Intramolecular Proton Transfer (ESIPT) is a more specific transfer that is triggered by optical and/or electronic excitation.⁵ This type of proton transfer, which can occur on timescales spanning $10^1 – 10^3$ fs,⁵ has garnered attention in the last 20 years since it holds several interesting characteristics. For example, the extremely large Stokes-shift in optical emission following ESIPT in some molecules (≥ 8000 cm⁻¹) could make it an efficient fluorescent sensor for optical sensing⁶⁻⁸ and biological imaging;⁹,¹⁰ and the related dual-colour emission from some ESIPT systems have been proposed for white-light generation in solid-state (organic) LEDs.⁵ ESIPT is also of interest, as it allows transport and storage of excitations in energy conversion processes,⁵,¹¹ light-driven enzymatic and artificial enzyme biosynthesis,¹² and may provide a way to control ultrafast energy transfer in the nascent field of polaritonic chemistry.¹³ A comprehensive list of other possible applications can be found in a recent review by Joshi *et al*.⁵

Although fundamental, ESIPT is a challenging mechanism to study, both experimentally and theoretically.¹⁴,¹⁵ The hydrogen being a light nucleus, proton transfer can occur on an
ultrafast timescale that is at the limit of many experimental probes, and the identification of
the key 'reaction coordinates' for ESIPT rapidly becomes problematic for increasing molecule
sizes, with some systems suggested to undergo ESIPT under the control of the molecular
'skeletal' vibrations,¹⁶⁻¹⁸ or directly along through-space, hydrogen-bond coordinates that
could themselves be multidimensional, in the presence of protic solvent or protein environ-
ments.¹⁹,²⁰ ESIPT, and more generally proton transfer, can also be highly correlated with
electron transfer, leading to proton-coupled electron transfer (PCET).⁴,²¹ This is a rich and
important field in its own right, but lies beyond the scope of the present article.

Given the wide dynamical range of typical molecular, solvent and host-material vibra-
tions, ultrafast ESIPT will manifestly involve a number of out-of-equilibrium vibronic states
whose dynamics and fluctuations will be temporally entangled with the dissipative (irre-
versible) dynamics of ESIPT. These lead to strong 'memory'-like effects that can cause the
fate of the process to depend on the way the system was initially *prepared*. Indeed, while
describing strongly coupled vibrational 'environments' always entails the consideration of
a large number of degrees of freedom, here, the situation is made more challenging by the
fact that no single dynamical approximation, such as the Markov approximation, can be
invoked for the entire bath: the proton-environment behaves as a quantum many body sys-
tem whose 'reaction' coordinates may be time-dependent and/or collective/hybrid/emergent
modes of the 'environmental' degrees of freedom. However, experimental efforts have led to
approaches based on absorption and fluorescence to probe the environment, characterised
by its spectral density.²²⁻²⁴ A number of time-resolved optical spectroscopy techniques can
also now potentially probe such details of ESIPT dynamics,²⁵⁻³⁰ including X-ray methods
that also report on real-time changes of molecular structure during ESIPT.³¹ Predicting
and understanding the underlying photodynamics of ESIPT at a numerically exact level of
theory – limited only by the quality of model parameterization – is thus a timely problem.

Previous studies of ESIPT have used a variety of theoretical methods to treat particu-
lar aspects of this problem at different levels of approximation. Approaches starting from
first principles, and which highlight the role of the molecular and heterogeneous environment on ESIPT, range from QM/MM free energy simulations invoking semi-classical (adiabatic) approximations, to simulations that predict non-adiabatic reaction rates.<sup>32,33</sup> Recently, a promising *ab initio* approach that blends real-time propagation of TD-DFT with quantum nuclear-electronic orbital dynamics has been presented and applied to ESIPT in condensed-phase environments.<sup>19,34,35</sup> When potential energy surfaces are available, other methods of treating the quantum nature of the proton dynamics and non-adiabaticity include semi-classical trajectory methods such as Surface Hopping (SH),<sup>36</sup> or Exact Factorisation.<sup>37</sup> For model vibronic systems, the Multi-configurational Time-Dependent Hartree (MCTDH) method<sup>38-42</sup> has made a very significant contribution in the last decades, being able to obtain exact quantum dynamics of non-adiabatic phenomena. The studied mechanisms with MCTDH span from exciton transport,<sup>43</sup> ultrafast charge separation<sup>44</sup> to photodissociation.<sup>45</sup> Being based on the treatment of continuous degrees of freedom, the family of MCTDH methods is indeed a suited choice to study proton transfer dynamics onto potential energy surfaces.<sup>46-50</sup> An excellent review of these – and other – theoretical chemistry approaches to ESIPT has been given by Jankowska *et al.*<sup>14</sup>

In this article we attempt another fully quantum mechanical approach, based on the idea that ESIPT can be considered as a complex type of open quantum system.<sup>51,52</sup> Indeed, a recent study of ESIPT of 10-hydroxybenzo[h]quinoline (HBQ) and 2-(2-hydroxyphenyl) benzothiazole (HBT) by Zhang *et al.* examined the non-adiabatic ESIPT dynamics in a diabatic picture of electronic-proton states coupled via a 1D reaction coordinate in the presence of a bath representing the rest of the molecule vibrations.<sup>52</sup> By treating the environment as a simple Markovian bath, they derived a Markovian Lindblad master equation for ESIPT, including the explicit role of the laser excitation, which is problematic for many of the methods mentioned above, such as the commonly used surface hopping methods. Inspired by this viewpoint, here we extend the numerically exact Time-Evolving Density with Orthogonal Polynomial (TEDOPA) method for Non-Markovian open systems to the problem of ESIPT.
TEDOPA has been established over the last decade as a state-of-the-art method for treating the coupling of a quantum system to its environmental excitations in a range of condensed matter, quantum optical and molecular problems, and is particularly efficient in the case of linear-vibronic models.53–59 At the heart of the method is the efficient representation of the evolving system-environment wave function in terms of either Matrix Product States (MPS) or Tensor Network States (TNS),60 with rather interesting connections to the tree-structures found in ML-MCTDH.61 The proximity of these two types of method allows many comparisons. Generally, ML-MCTDH is mainly performed in studying chemical problems while MPS is more established in the physics community. It has to be pointed out that an equivalent of TEDOPA called hierarchical effective-mode chain representation62–67 is also commonly used with ML-MCTDH. The essential properties of these MPS states and why these are particularly powerful for open system problems are presented in Section 2, but a key feature is the access to the complete quantum state of the system *and* its environments. Unlike master equations in which the environment has been removed from the problem (traced out), TEDOPA allows one to track and visualize the environmental dynamics leading to dissipative phenomena in the system.56,57,68

This capability could be particularly useful for photo-chemistry problems, as all modes of the environment are treated on an equal (exact) footing, and interrogation of the environment state can provide ways to identify relevant reaction coordinates on the fly (important in multi-stage dynamics), and/or visualize high dimensional dynamics in a purely diabatic basis.57,59 Indeed, the use of TNS and MPS wave functions, also known as tensor-trains in the applied mathematics literature,69,70 for simulating high-dimensional excited-state dynamics has seen tremendous growth over the last few years.71–76

However, previous work employing TEDOPA has typically explored how environments drive transitions between discrete states of point-like systems, although it is capable of treating spatial-temporal correlations between point-like systems embedded in different regions of a common bath.68 For ESIPT systems possessing well-defined reaction coordinates, describ-
ing dissipative dynamics requires an extension of the current algorithms to treat a system
evolving on one or more *continuous* spatial degrees of freedom. This article presents our first
efforts in this direction, where we use toy models similar to those of Zhang *et al.* to demon-
strate that TEDOPA can be used to visualize dissipative dynamics of ESIPT in real space
with an exact treatment of the vibrational environment, the coupling to the electromagnetic
environment that leads to dual fluorescence, as well as driving by laser pulses of arbitrary
strength.

In this work, we progressively present these capabilities. After introducing the MPS
framework in Section 2, we start from a purely open system model without reaction coor-
dinates (Section 3) to show how access to environmental information allows us to correlate,
for example, the real-time emergence of dual fluorescence with the energy dissipated into
the vibrational environment. This model also allows us to explore the role of strong driving
on the emission properties and dynamics of ESIPT, anticipating future work on how po-
laritonic effects could also alter ESIPT phenomena. We then introduce an explicit reaction
coordinate in Section 4 and present examples where we are able to give a complete spatio-
temporal description of ESIPT, including dissipation and an exact treatment of decoherence
in the non-adiabatic dynamics, which often has to be treated in an *ad hoc* fashion in surface
hopping methods.⁷⁷ Importantly, we also introduce in Section 4.2 a matrix product oper-
ator projector that allows to re-express the dynamics in the chemically familiar language
of adiabatic potential surfaces, allowing for future comparisons with other more established
formalisms, and which also provides a quantitative measure of the 'non-adiabaticity' of the
dynamics. We then summarise our results in Section 5. Overall, this yields a versatile
MPS/TNS framework designed for the study of excited-state chemical problems, as well as
the study of non-adiabatic transitions, which we hope will be useful for researchers working
in these areas.
## 2 A Brief Introduction to Matrix Product States

The calculations in this work make use of the Matrix Product State (MPS) Ansatz for many-body quantum states.⁶⁰,⁷⁸,⁷⁹ A general N-body wave function $|ψ⟩$ can be written on the basis of, say, d one-body substates $\{|\phi_i⟩\}$ as $|ψ⟩ = \sum_{\{i\}} C_{i_1 i_2 \dots i_N} |\phi_{i_1} \dots \phi_{i_N}⟩$. The coefficients $C_{i_1 i_2 \dots i_N}$ are the amplitudes of each N-body configuration to the full system state. Clearly, the total number of amplitudes scales as $d^N$, leading to the ‘curse of dimensionality’ that rapidly makes any direct numerical treatment extremely expensive.

However, starting from a well-defined initial configuration, most physical Hamiltonians will only allow a much smaller sub-space of states to be explored in a finite time. Additionally, if the initial state is a product state, then it seems reasonable that the size of this subspace can be *classified* by the amount of the quantum correlation (entanglement) that is generated during the dynamics. Indeed, in the absence of any entanglement, the problem would reduce to a type of mean-field description that would have linear scaling in $N$.

Matrix Products States are essentially low-rank decompositions of the full wave function that are formed by breaking up the multi-linear array (a tensor) $C_{i_1 i_2 \dots i_N}$ into a product of smaller data arrays associated with each particle.⁶⁰,⁸⁰ The MPS decomposition can be represented with the help of diagrams as in Figure 1. In this formalism, the tensor $C_{i_1 i_2 \dots i_N}$ is written as several smaller tensors $T$ (Figure 1 (a)). Each tensor $T$ has three legs (array dimensions), with legs pointing towards the bottom of the diagram ($i_j, j \in [1, N]$) running over physical dimension of the one-particle states on that site (for example the Fock dimension or the number of electronic states). The legs linking tensors represent the entanglement between tensors $T_j$ and $T_{j+1}$, and contain $D_j$ parameters, known as the local ‘bond dimension’ that allow the MPS to describe non-classical correlated states ($D = 1$ would be a product state).

Tensors sharing a common leg implies contraction (Einstein summation, or equivalently a trace) over the shared dimensions of their arrays. Thus for a particular physical state configuration, i.e. if the elements of each downward pointing leg are specified, it can be
seen that the scalar coefficient $C_{i_1 i_2 \dots i_N}$ is obtained from contractions of all the tensors
in the state, which in 1D corresponds to matrix multiplication. The key aspect of the
MPS Ansatz is that the number of parameters required to describe the state is now $N \times$
$d \times D^2$, which has linear scaling with $N$, provided that $D$ is not itself a function of $N$.
Importantly, rigorous results concerning the so-called 'area laws' show that this is indeed the
case for gapped 1D Hamiltonians with short-range interactions,$^{80}$ so that the $\{D_j\}$ can be
used as convergence parameters for numerical simulations of such models.$^{60,80}$ Crucially for
the dynamics of open quantum systems, it is possible to prove that many common system-
environment Hamiltonians are formally equivalent to 1D models with nearest-neighbour
couplings, and that these models can be easily constructed via a unitary transformation
of the environmental degrees of freedom.$^{53,81,82}$ This is the underlying reason for why MPS
are particularly useful for system-environment problems, and an example of the TEDOPA
workflow, i.e. model→chain mapping→MPS simulation is given in Section 3. For the results
presented in this paper, we typically find converged results for chains with $D \approx 10$, $d = 10$
and number of modes going from $N = 150$ to $N = 540$. Only for the less demanding
simulation, it results to a total of $N \times d \times D^2 = 15 \times 10^4$ parameters in comparison with
the full-rank wave function that would have $d^N = 10^{150}$ elements, making the numerical cost
reduction clear. For more details about parameters, see Tables (S1-S7) in the Supporting
Information.

Within the same formalism, operators can also be expressed as Matrix Product Operators
(MPOs) by writing them as a chain of local tensors $\hat{O} = \sum_{\{i\},\{i'\}} W_{i_1 i'_1 i_2 \dots i_N i'_N}^{i_1 \dots i_N} |\phi_{i_1} \dots \phi_{i_N}\rangle \langle \phi_{i'_1} \dots \phi_{i'_N}|$
(Figure 1 (b)).$^{83}$ Accordingly, an MPO can act on an MPS to lead to a new wave-function
(Figure 1 (c)). Effectively, the MPO form allows us to measure observables of the system –
or the environment – as the action of a local operator acting on the states (tensor) of a given
site, which can be combined with the gauge freedom of the MPS representation to create a
highly efficient procedure for evaluating multi-partite observables.$^{79,80}$

Importantly, this also applies to the Hamiltonian operator in MPO form, allowing for
<figure><img src="image_2.png" /><figcaption>Figure 1: Diagram representing the Matrix Product State Ansatz. (a) The full wave-function represented by the multidimensional coefficient $C_{i_1 i_2 \dots i_N}$ can be approximated by dividing it into smaller tensors. The $i_1 \dots i_N$ coefficients are physical dimensions while $D_1 \dots D_N$ are bond dimensions, i.e. entanglement between tensors. (b) The MPO $\hat{O}$ is composed of local tensors $W_j$ connected by bond dimensions $w_j$. It has twice the number of physical dimension to allow $\hat{O}|\psi\rangle$ and $\langle\psi|\hat{O}$ operations. (c) $W_j$ MPO tensors are acting on MPS tensors $T_j$. Common physical dimensions between the MPS and the MPO are merged to lead to a new wave-function. Legs of dimension one at the extremities of the chains are not represented in the diagrams.</figcaption></figure>
efficient time-evolution of the state by sweeping through a sequence of updates that act on
no more than one or two site tensors at a time.⁶⁰,⁸⁰,⁸⁴ For (bosonic) open system problems,
a particularly efficient and versatile evolution method has been shown to be the one-site
*Time Dependent Variational Principal* (1TDVP) and its recent extensions to allow for dy-
namic updates to the size of the bond dimensions.⁷²,⁸⁴⁻⁸⁸ Further details about TDVP and
its algorithmic implementation can be found in Ref.⁸⁵ In this article, all of the TEDOPA
calculations are performed using our source-source Julia package MPSDynamics.jl⁸⁹,⁹⁰. This
package, along with extensive documentation, tutorial examples and all of the models that
we have been studying for research problems in physics and chemistry can be found at
https://github.com/shareloqs/MPSDynamics. We now return to the particular applica-
tion for this article: ESIPT.

3 A Four-level Open-System Description of ESIPT

Excited-state intramolecular proton transfer (ESIPT) is often described as a four-state cycle,
representing the enol electronic ground state $|e\rangle$, the enol electronic excited state $|e^*\rangle$ and
their keto counterparts ($|k\rangle$, $|k^*\rangle$), with the ESIPT process corresponding to the $|e^*\rangle \rightarrow |k^*\rangle$
transfer. In order to study an excited electronic transfer as well as emission from these
states, the full four electronic-state cycle is simulated. Interactions with environments can
be included along with the main system description. The coupling between a system and
a discrete set of modes can be efficiently implemented within the MPS approach, without
any loss of information. These modes can describe different environments such as molecular
vibrations and photons. In fact, photons are treated in a similar manner as vibrations: both
are represented by a set of discrete harmonic modes described by the associated (bosonic)
creation and annihilation operators. Such environments are characterised by their spectral
density, specifying the interaction between the system and each mode of the environment.
The quantum dynamics of the open system wavefunction (including the environments) allows
to capture the effect of the motion of the molecule on the electronic population dynamics as well as photo-emission and absorption processes. These photonic environments are of particular interest for dynamics of systems and quantized fields have, for instance, been simulated along vibrational environments with MCTDH in order to study polariton dynamics.⁹¹,⁹²

<figure><img src="image_3.png" /><figcaption>Figure 2: Four-state cycle that illustrates ESIPT. Interaction with environments allows transitions between states. Transitions from excited state to ground state correspond to photon coupling ($a_k^\dagger, a_k$ in eq 1.1-1.5) whereas transitions between excited states or ground states are due to vibration coupling ($b_{k'}^\dagger, b_{k'}$ in eq 1.1-1.5). The transfer $|e^*\rangle \rightarrow |k^*\rangle$ represents the excited-state proton transfer. A general ESIPT compound is shown on the bottom with $R_1, R_2$ and $R_3$ as general moieties.</figcaption></figure>

To illustrate the effect of these couplings, the prototype system of Figure 2 is studied. In this model, there is no direct coupling between electronic states. Instead, environments medi-
ate transitions between states: intramolecular vibrations allow to switch from one electronic structure to another and photons allow vertical transitions.

## 3.1 Hamiltonian of Electrons Interacting with Photons and Vibrations

To represent a four-state cycle interacting with environments (Figure 2), a Hamiltonian can be written to account for these couplings. This Hamiltonian is composed of a system part $H_S$, a time-dependent term modelling the effect of an external driving field $H_{drive}(t)$, a bath component $H_B$ and an interaction term $H_{SB}$. The full Hamiltonian reads (here and in the following, we use arbitrary units in which $\hbar = 1$):

$$H(t) = H_S + H_{\text{drive}}(t) + H_B + H_{\text{SB}} \quad (1.1)$$

$$H_S = \omega_e |e\rangle \langle e| + \omega_{e^*} |e^*\rangle \langle e^*| + \omega_k |k\rangle \langle k| + \omega_{k^*} |k^*\rangle \langle k^*| \quad (1.2)$$

$$H_{\text{drive}}(t) = (|e^*\rangle \langle e| + |k^*\rangle \langle k| + \text{h.c.}) \epsilon(t) \sin(\omega_{\text{drive}}t) \quad (1.3)$$

$$H_B = \overbrace{\sum_k \omega_k^{\text{phot}} (a_k^{\dagger} a_k)}^{\text{Photons}} + \overbrace{\sum_{k'} \omega_{k'}^{\text{vib}} (b_{k'}^{\dagger} b_{k'})}^{\text{Vibrations}} \quad (1.4)$$

$$H_{\text{SB}} = (|e^*\rangle \langle e| + |k^*\rangle \langle k| + \text{h.c}) \sum_k g_k (a_k^\dagger + a_k) + (|e^*\rangle \langle k^*| + \text{h.c}) \sum_{k'} y_{k'} (b_{k'}^\dagger + b_{k'}) \quad (1.5)$$

with $a_k^{(\dagger)}$ the annihilation (creation) operator for the $k^{\text{th}}$ mode in the photon bath of frequency $\omega_k$, and $b_{k'}^{(\dagger)}$ the annihilation (creation) operator for the $k'^{\text{th}}$ mode in the vibration bath of frequency $\omega_{k'}$. $\omega_e, \omega_{e^*}, \omega_k$ and $\omega_{k^*}$ are respectively the energy of the enol ground state, enol excited state, keto ground state and keto excited state. We also define $\Delta\omega_e = \omega_{e^*} - \omega_e$ and $\Delta\omega_k = \omega_{k^*} - \omega_k$. The coefficients $g_k$ and $y_{k'}$ are determined for the different $k$ and $k'$ indices defining the weights of the spectral density characterising the coupling between environments and the system for every mode. The spectral densities read $J_{\text{phot}}(\omega) = \sum_k g_k^2 \delta(\omega_k - \omega)$ and $J_{\text{vib}}(\omega) = \sum_k y_k^2 \delta(\omega_k - \omega)$. Although the TEDOPA method (associated with chain
mapping) can handle arbitrary spectral densities that can be tuned to represent a spe-
cific environment,⁹³ here the spectral densities are defined generically, with an Ohmic form.
The time-dependent driving field is characterised by its frequency $\omega_{drive}$ and its (possibly
time-dependent) amplitude $\epsilon(t)$. The amplitude function $\epsilon(t)$ is time-independent for a
continuous-wave driving field, and can be made time-dependent to model pulsed excitation.
Both the driving field and the photon bath can induce optical transitions between electronic
states. The driving field is used to model the (continuous or pulsed) laser excitation of the
system, while the photon bath captures fluorescence.

The vibration bath couples the states $|e^*\rangle$ and $|k^*\rangle$ as well as $|e\rangle$ and $|k\rangle$ whereas photons couple the excited states with their respective ground states. The transfer between $|e^*\rangle$ and $|k^*\rangle$ is the excited state proton transfer and the transfer between $|k\rangle$ and $|e\rangle$ corresponds to ground state proton transfer. In this particular example, the ground state proton transfer is ignored.

<figure><img src="image_4.png" /><figcaption>Figure 3: Star conformation. The electronic system (orange) is interacting with two different environments : *N* photon modes (red) characterised by the spectral density *J*<sub>phot</sub>(*ω*) = ∑<sub>*k*=1</sub><sup>*N*</sup> *g*<sub>*k*</sub><sup>2</sup>δ (*ω*<sub>*k*</sub> − *ω*) on one side and with *M* vibration modes (blue) characterised by the spectral density *J*<sub>vib</sub>(*ω*) = ∑<sub>*k*=1</sub><sup>*M*</sup> *y*<sub>*k*</sub><sup>2</sup>δ (*ω*<sub>*k*</sub> − *ω*) on the other side.</figcaption></figure>

This Hamiltonian corresponds to an electronic system interacting with several vibration
modes of different frequencies. It is represented in Figure 3 in the form known as the 'star-like'
conformation, where each vibration mode is coupled schematically to the system situated at the centre. The star conformation is not directly adapted to the MPS formalism that naturally deals with chain-like systems in which each environment mode interacts with only its nearest neighboring modes as in Figure 1. In this work, we make use of the *Time Evolving Density matrix with Orthogonal Polynomials Algorithm* (TEDOPA) procedure.⁵³,⁸¹,⁸² This approach allows to transform the initial star conformation into a chain, via a change of basis (unitary transformation) based on orthogonal polynomials. In brief, given an environment characterized by a continuous spectral density $J(\omega)$, we express the original environmental modes in a new set of bosonic modes $\{c_n, c_n^\dagger\}$ as $a_k = \sum_{n=0}^{\infty} \sqrt{J(\omega_k)} \pi_n(\omega_k) c_n$, where $\pi_n(\omega)$ is a nth-order orthogonal polynomial w.r.t. to the 'weight function' $J(\omega)$.⁸² Substituting this into the original Hamiltonian, and using the three-term recurrence properties of orthogonal polynomials, transforms $H_B$ into its chain form $H_C$

$$H_C = \sum_{n=0}^{\infty} \epsilon_n c_n^\dagger c_n + t_n c_n^\dagger c_{n+1} + t_n c_n c_{n+1}^\dagger, \quad (2)$$

where all microscopic information about the environment, i.e. $J(\omega)$ is now encoded in the chain parameters $\{\epsilon_n, t_n\}$. In certain cases, such as the widely used Ohmic, sub-Ohmic and super-Ohmic spectral densities, these chain parameters can even be determined analytically,⁸² but efficient numerical procedures implemented in MPSDynamics.jl allow these parameters to be extracted for spectral densities with arbitrary frequency dependencies.

It can further be shown that in the chain representation the system interacts with only the first mode of the chain-like environment and MPS techniques can then be easily applied. The MPS is then finally built as the system surrounded by a first chain of photons and a second chain of molecular vibrations as shown in Figure 4. This practice does not change the simulated wave-function but only the way to calculate dynamics. Transforming it into a chain allows to reach the full power of MPSs for a linear vibronic coupling. Moreover, the inverse mapping can always be performed, allowing us to use the chain results to extract
observables and dynamics of the original bath degrees of freedom.⁶⁸,⁹⁴,⁹⁵ Indeed, we will make extensive use of this back-transformation throughout this article.

<figure><img src="image_5.png" /><figcaption>Figure 4: Chain conformation. The electronic system (orange) is interacting with two linearized environments : a chain of N photons (magenta) characterised by the parameters {$\epsilon_n^{\text{phot}}$, $t_{n-1}^{\text{phot}}$} with $n \in [1, N]$ and a chain of M vibrations (dark blue) characterised by {$\epsilon_m^{\text{vib}}$, $t_{m-1}^{\text{vib}}$} ($m \in [1, M]$). In this conformation, the system is interacting with only two modes and MPS techniques can be used.</figcaption></figure>

## 3.2 Time-dependent Approach

To simulate photo-chemistry processes, it can be necessary to access ultrafast dynamics: initialising the wavepacket onto the excited state at the beginning of the simulation might wipe out rapid features. For this purpose, the time-dependent Hamiltonian $H_{\text{drive}}(t)$ is introduced within the MPS method in this work, representing a laser driving (either continuous or pulsed) that promotes the electronic population to the excited state.

In order to do so, the Hamiltonian is written as a Matrix Product Operator (MPO) and is modified at each time step to let the time-dependent part evolve. Practically, the time-dependent terms of the MPO are added to the time-independent MPO terms at each timestep in order to avoid the calculation of an entire new MPO. In this way, only a matrix sum is computed at each time step.

The wide possibilities of the generic four-level setup described by equations (1.1) to (1.5), treated within the MPS framework and including the time-dependent feature, are illustrated on several model examples in the following paragraphs.
### 3.3 Excited-State Dynamics and Dual-Fluorescence without Laser Driving

First, we consider a simple four-level model where the initial population was prepared in the $|e^*\rangle$ state. The system here interacts with vibrations and photons and is governed by the Hamiltonian of eq. (1.1) to (1.5). In this first example, the Hamiltonian is time independent ($\epsilon(t) = 0$). Examples of this section are provided in arbitrary units and describe very generic models. The corresponding dynamics can be easily rescaled by an arbitrary time factor. all the parameters are reported in Table S1 of the Supporting Information. Beginning the dynamics with the $|e^*\rangle$ state populated, the electronic population is transferred from $|e^*\rangle$ to $|k^*\rangle$ as shown in Figure 5 (a). Ground states $|e\rangle$ and $|k\rangle$ are also populated to a lesser extent via spontaneous photon emission (Figure 5 (c)). Looking at the baths, vibration modes are populated in time mainly at one frequency (Figure 5 (b)) while two frequencies are excited in the photon bath (Figure 5 (d)).

This example illustrates the effect of environments on the electronic system. While there is no direct coupling between electronic states in $H_S$, the coupling to the vibration bath allows transitions from $|e^*\rangle$ to $|k^*\rangle$ as shown in Figure 5 (a). Because of the energy transfer between the system and vibrations, vibration modes are excited at the frequency corresponding to the energy difference between $|e^*\rangle$ and $|k^*\rangle$ (Figure 5 (c)). Coupling to the photon bath also triggers vertical transitions, thereby populating the electronic ground states (Figure 5 (b)). Two emission lines can be observed, increasing in intensity with time (Figure 5 (d)). These peaks are characteristic of dual fluorescence, i.e. simultaneous emission from the enol (at frequency $\Delta\omega_e$) and from the keto forms (at frequency $\Delta\omega_k$). The keto emission is delayed in comparison with the enol one because of the different initial population of the excited states. Although ESIPT happens in a sub-picosecond timescale, this time-resolved photo-emission spectrum can, in principle, be probed experimentally.²⁵,²⁷ This example shows that, by making use of the MPS wave-function method, information can be recovered from environments.
All of these observables together enable the energy transfers to be tracked, both in the system and in the environments.

## 3.4 Excited-State Dynamics under Laser Driving

We now consider a dynamics initiated in the enol ground state $|e\rangle$ and introduce the effect of the driving field via the time-dependent term $H_{\text{drive}}(t)$, in order to populate the excited state (eq. 1.3). The system and bath characteristics are similar to the previous example, while the driving is continuous (i.e. the amplitude $\epsilon(t)$ is time-independent). The field frequency is set up at resonance: $\omega_{\text{drive}} = \Delta\omega_e$ in order to populate the enol excited state and triggers the first step of the ESIPT process. Since the keto energy gap is off resonance, the continuous drive has a negligible effect on the keto states. This type of driving, when at resonance, leads to an almost full transfer. For parameters, see Table S2 of the Supporting Information.

The electronic population dynamics is presented in Figure 6. The $|e\rangle$ and $|e^*\rangle$ populations show damped oscillations in time, with a characteristic period close to 800 arb. units, corresponding to the Rabi period of the drive. In the long term, the total enol population goes to zero, while the keto population $|k^*\rangle$ increases towards 1. Vibrations allow transfers between $|e^*\rangle$ and $|k^*\rangle$ so that $|k^*\rangle$ becomes populated as soon as $|e^*\rangle$ is populated by the drive. In the long term, the population is entirely transferred from enol to keto via the excited states. Environment spectra for this example are not shown but they closely resemble that of Figure 5.

As a last example, we consider the same system under a more intense driving (all parameters are reported in Table S3 of the Supporting Information). It corresponds to strong light-matter coupling, that could be reached by the means of an intense laser excitation or in a resonant optical cavity. Figure 7 shows the excitation spectrum for the photon environment (Figure 7 (a),(b)) and for the vibration environment (Figure 7 (c),(d)). Four main
<figure><img src="image_6.png" /><figcaption>Figure 5: Dual-fluorescence example of a four-level system. (a) Dynamics of the electronic population. The dynamics starts in $|e^*\rangle$ and vibrations allow a transition towards $|k^*\rangle$. (b) Zoom in the dynamics of the electronic population. Photons allow to populate ground states. (c) Vibration average population resolved in time. An horizontal dashed line indicates the frequency difference $\omega_{e^*} - \omega_{k^*}$. (d) Photon average population resolved in time. Two different frequencies are detected, corresponding to the two electronic energy gaps. Time and frequencies are in arbitrary units.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6: Driven four-level system. The electronic dynamics starts in $|e\rangle$ and $|e^*\rangle$ is populated with the continuous driving. Transfer in $|k^*\rangle$ happens due to the coupling between $|e^*\rangle$ and $|k^*\rangle$ with vibrations. Time is in arbitrary units.</figcaption></figure>
photo-emission lines can be identified in the photon bath results, and two different molecular vibration frequencies are excited by the system. The corresponding electronic population dynamics is reported in the Supporting Information (Figure S1).

The formation of these peaks is characteristic of the strong coupling regime and can be explained analytically: the strong coupling between the $|e\rangle \rightarrow |e^*\rangle$ transition and the driving field gives rise to enol dressed states. As a result, the $|e\rangle \rightarrow |e^*\rangle$ emission line splits in three distinct peaks, corresponding to transitions between the different dressed states (Fig. 8). This characteristic signature of strong light-matter coupling is known as Mollow triplet.⁹⁶ The central peak corresponds to the two degenerate transitions that remain at the initial frequency $\Delta\omega_e$, while two other transitions are found at frequencies $\Delta\omega_e + \epsilon$ and $\Delta\omega_e - \epsilon$. Furthermore, the state $|k^*\rangle$ now interacts with the two $|e^*\rangle$ dressed states, which also explains the splitting of the vibration spectrum peak. Finally, though the coupling between the driving field and the keto states is not on resonance, it still has a measurable impact on these states and results in a slight shift of the keto energy gap $\Delta\tilde{\omega}_k$. All the shifts arising from the strong coupling and the emergence of dressed states can be anticipated on the basis of a simple analytical model that can be found in Supplementary Information.

In Figure 7, the lowest frequency photon signal corresponds to the dressed transition $|k^*\rangle \rightarrow |k\rangle$ (Figure 7 (a), (b)) while the three other peaks constitute the Mollow triplets. The highest central peak corresponds to the two degenerate enol transitions at $\Delta\omega_e$ and the two satellite peaks are at $\Delta\omega_e + \epsilon$ and $\Delta\omega_e - \epsilon$. High frequency transitions are favoured by the vibrational spectral density, which explains the slight difference in intensity between the Mollow triplet satellites: the upper enol dressed state being more strongly coupled to $|k^*\rangle$, it depopulates more quickly and contributes less to the photo-emission spectrum. We also note that the keto emission is differed in time from the enol emission, as a result of the delay induced by the $|e^*\rangle \rightarrow |k^*\rangle$ proton transfer step.
The two peaks in the vibration spectrum (Figure 7 (c), (d)) correspond to transitions from the two dressed states of $|e^*\rangle$ to $|k^*\rangle$. Two excited-state transitions are now possible at $\omega_{e^*} - \tilde{\omega}_{k^*} + \epsilon/2$ and at $\omega_{e^*} - \tilde{\omega}_{k^*} - \epsilon/2$ (Figure 8). Finally, strong coupling with the driving field gives rise to Rabi oscillations in time in both the vibration and photon environments. These oscillations are shown in more detail in Supplementary Information (Figure S2).
<figure><img src="image_8.png" /><figcaption>Figure 7: Strongly-driven four-level system. (a) Photon average population resolved in time. A Mollow triplet appears alongside the keto transition. (b) Photon average population at $t_{final} = 400$ arb. units. (c) Vibration modes average population resolved in time. Due to the strong driving, two transitions are now possible between $|e^*\rangle$ and $|k^*\rangle$. (d) Vibration modes average population at $t_{final} = 400$ arb. units. Time and frequencies are in arbitrary units.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8: Four-level system strongly driven on resonance with $\Delta\omega_e$. It creates enol dressed states with a new $\epsilon$ gap and three different optical transitions, known as a Mollow triplet. Two different transitions between the dressed states $|e^*\rangle$ and the state $|k^*\rangle$ are now possible.</figcaption></figure>
# 4 From Discrete States to Non-Adiabatic Dynamics on
Continuous Energy Landscapes

The ESIPT process involves not only an electronic transition between the states $|e^*\rangle$ and $|k^*\rangle$, but also a proton transfer. As the two phenomena are highly correlated, an explicit description of both mechanisms is crucial.⁹⁷ In contrast, in the previous 4-state model, the proton transfer is only described implicitly in the $|e^*\rangle \rightarrow |k^*\rangle$ transition, mediated by the coupling with the vibration environment. The MPS formalism is generally not used to describe continuous reaction coordinates such as the position of the transferring proton. In this section, we develop a model that allows the consistent description, within the MPS Ansatz, of the proton wavefunction along its path from the donor to the acceptor site. The associated continuous reaction coordinate (RC) is described by means of an additional vibration mode, with appropriate couplings to the electronic states in order to generate the double-well energy profile typical of proton transfer processes.

## 4.1 Hamiltonian

To set up the RC feature in the MPS method, additional terms are included in the Hamiltonian. In this section, we focus on the electronically excited states to simplify the description, the ground states $|e\rangle$ and $|k\rangle$ can be treated in exactly the same way. A coupling term $\Delta$ is introduced between the electronic states $|e^*\rangle$ and $|k^*\rangle$, while the RC mode is described by the oscillator annihilation and creation operators $d, d^\dagger$ with a characteristic frequency
$\omega_{RC}$:

$$H = H_S + H_{RC} + H_{int}^{RC-S} \quad (3.1)$$

$$H_S = \omega_{e^*} |e^*\rangle \langle e^*| + \omega_{k^*} |k^*\rangle \langle k^*| + \Delta (|e^*\rangle \langle k^*| + |k^*\rangle \langle e^*|) \quad (3.2)$$

$$H_{RC} = \omega_{RC} \left( d^{\dagger} d + \frac{1}{2} \right) \quad (3.3)$$

$$H_{int}^{RC-S} = g_{e^*} |e^*\rangle \langle e^*| (d + d^{\dagger}) + g_{k^*} |k^*\rangle \langle k^*| (d + d^{\dagger}) \quad (3.4)$$

The coefficients $g_{e^*}$ and $g_{k^*}$ are the interaction weights between the electronic and RC parts of the system. The Hamiltonian above includes two electronic states coupled to an oscillator that can be represented in a truncated Fock states basis. It is in this basis that the Hamiltonian is expressed as an MPO and the dynamics propagated exactly thanks to the TEDOPA method. However, in the analysis of the resulting quantum dynamics, different choices of basis might be used in order to obtain a more expressive representation of the photochemical process. In particular, introducing the RC position operator allows to interpret $g_{e^*}$ and $g_{k^*}$ as spatial shifts of the RC equilibrium position (different for the enol and keto forms), and a transformation of the electronic states towards the adiabatic basis then gives rise to the characteristic double-well potential picture of proton transfer.

**Continuous description in space.** A more intuitive representation of the system is indeed obtained by re-expressing the Hamiltonian on the basis of the RC spatial coordinate. We first recall that the position operator is defined from the annihilation and creation operators as $x = \sqrt{\frac{1}{2m_{RC}\omega_{RC}}}(d + d^{\dagger})$. The Hamiltonian of eq 3.1 can then be rewritten as:

$$H = H^{e^*}(x) |e^*\rangle \langle e^*| + H^{k^*}(x) |k^*\rangle \langle k^*| + \Delta (|e^*\rangle \langle k^*| + |k^*\rangle \langle e^*|) \quad (4)$$

where $H^{e^*}(x) = \frac{1}{2}m_{RC}\omega_{RC}^2 (x - x_{RC}^{eq,e^*})^2 + \omega^{eq,e^*}$ and similarly for $H^{k^*}(x)$. The equilibrium positions $x_{RC}^{eq,e^*}$ and $x_{RC}^{eq,k^*}$ of the excited enol and keto states are directly related to the
interaction coefficients via : $g_i = -\omega_{RC} x_{RC}^{eq,i} \sqrt{\frac{m_{RC}\omega_{RC}}{2}}$ and the energies at the two equilibrium positions follow, $\omega^{eq,i} = \omega_i - \frac{1}{2}m_{RC}\omega_{RC}^2 (x_{RC}^{eq,i})^2$. In this representation, the system can be seen as consisting in two displaced harmonic oscillators for $|e^*\rangle$ and $|k^*\rangle$, with different equilibrium positions and energies and coupled via the tunnelling term $\Delta$. The different parameters enable a rich variety of configurations for the corresponding potential surface, an example is represented by the dashed lines in Figure 9.

## 4.2 The Adiabatic Projector

The inclusion of the continuous RC variable leads to an extended system described in terms of electronic state and spatial coordinate, as shown in Figure 9. The MPS dynamics is performed in the diabatic representation, i.e. with $|e^*\rangle$ and $|k^*\rangle$ as basis states, whereas most photo-chemistry problems are generally described in the adiabatic basis.⁹⁸,⁹⁹ For this reason, we developed an adiabatic projector (expressed in the MPO representation) that enables easy transformations from one basis to the other. The adiabatic projector is here only used as an analysis tool in order to investigate adiabatic observables. It does not affect the potential energy surface or the quantum dynamics propagation. Similar projectors for several degrees of freedom are implemented in the MCTDH¹⁰⁰ and QUANTICS¹⁰¹ packages for MCTDH methods.¹⁰²⁻¹⁰⁵ Although the presented adiabatic projector deals with only one degree of freedom, it demonstrates the new possibility of working in both bases with MPS.

**Diabatic and adiabatic basis.** Two different bases are frequently used to represent this type of systems at a fixed RC position $x$. The initial Hamiltonian of eq 4 is written in the so-called diabatic electronic basis ($|e^*\rangle$, $|k^*\rangle$) in which the MPS dynamics is carried out. The
Hamiltonian is represented by the following matrix:

$$H^{\text{Diab}}(x) = \begin{pmatrix} H^{e*}(x) & \Delta \\ \Delta & H^{k*}(x) \end{pmatrix} \qquad (5)$$

where the anti-diagonal terms are the coupling between the electronic states. The adiabatic basis is then obtained by diagonalizing $H^{\text{Diab}}(x)$ at fixed position $x$:

$$H^{\text{Adiab}}(x) = \begin{pmatrix} H^{\text{LOW}}(x) & 0 \\ 0 & H^{\text{UP}}(x) \end{pmatrix} \qquad (6)$$

The adiabatic energies $H^{\text{LOW}}(x)$ and $H^{\text{UP}}(x)$ are represented in solid lines in Figure. 9. The corresponding eigenvector $|\chi^{\text{LOW}}(x)\rangle$ and $|\chi^{\text{UP}}(x)\rangle$ follow from the diagonalization of the matrix $H^{\text{Diab}}(x)$.

The lower potential surface $H^{\text{LOW}}(x)$ has a characteristic double-well shape. In the following, the RC is treated at the quantum level using the MPS representation, its dynamics therefore describes the evolution of the proton wavepacket, propagating onto the double-well potential energy surface after being photo-excited.

**Construction of the adiabatic projector.** The purpose of the adiabatic projector is to split the system's wavepacket into two parts : the one evolving onto the upper adiabatic surface and the other one onto the lower surface. The projector is built as an MPO, which allows for an efficient numerical treatment.

In order to build the projector, we first consider the normalised eigenstates of the undisplaced RC harmonic oscillator $H_{\text{RC}}$. These are the Fock states $|n\rangle$ over which the MPS dynamics is built, with the associated eigenfunctions $\phi_n(x)$ given by the usual Hermite polynomials. A large number of Fock states should be considered in order to obtain a well-resolved spatial description of the RC. The second ingredient to the adiabatic projector is the
<figure><img src="image_10.png" /><figcaption>Figure 9: Example of potential energy curves parameterised by the Hamiltonian of eq 3.1 in the adiabatic basis. The Hamiltonian is diagonalized leading to two distinct adiabatic surfaces, that are designed as LOW and UP surfaces. Diabatic surfaces are also represented underneath as dashed lines. We define $\Delta\omega = \omega_{e^*} - \omega_{k^*} + \frac{1}{2}m_{RC}\omega_{RC}^2((x_{RC}^{eq,k^*})^2 - (x_{RC}^{eq,e^*})^2)$. The Reaction Coordinate and $\omega$ are in atomic Hartree units.</figcaption></figure>
eigenstates of the diagonalization from the diabatic to the adiabatic electronic basis at fixed $x$ value: $|\chi^{\text{LOW}}(x)\rangle$ and $|\chi^{\text{UP}}(x)\rangle$. A generic system wavefunction can then be decomposed as:

$$ \begin{aligned} |\psi^{\text{Adiab}}(x,t)\rangle &= |\psi^{\text{LOW}}(x,t)\rangle + |\psi^{\text{UP}}(x,t)\rangle \\ &= |\chi^{\text{LOW}}(x)\rangle \phi_{\text{LOW}}(x,t) + |\chi^{\text{UP}}(x)\rangle \phi_{\text{UP}}(x,t) \end{aligned} \quad (7) $$

with $\phi_{\text{LOW}}(x,t)$ and $\phi_{\text{UP}}(x,t)$ the time-dependent adiabatic wavefunctions and $|\chi^{\text{LOW/UP}}(x)\rangle = \chi_{e^*}^{\text{LOW/UP}}(x) |e^*\rangle + \chi_{k^*}^{\text{LOW/UP}}(x) |k^*\rangle$ the local electronic eigenstates. The projector $\hat{P}_{\text{AD}}^{\text{LOW}}$ on the LOW adiabatic surface can then be written as (a similar form holds for the UP projector):

$$ \hat{P}_{\text{AD}}^{\text{LOW}} = \int dx |x\rangle \langle x| \otimes |\chi^{\text{LOW}}(x)\rangle \langle \chi^{\text{LOW}}(x)| \quad (8.1) $$

$$ = \sum_{n,m} \int dx \phi_n(x) \phi_m^*(x) |n\rangle \langle m| \otimes (\chi_{e^*}^{\text{LOW}}(x) |e^*\rangle + \chi_{k^*}^{\text{LOW}}(x) |k^*\rangle) \otimes \quad (8.2) $$

$$ (\chi_{e^*}^{\text{LOW}*}(x) \langle e^*| + \chi_{k^*}^{\text{LOW}*}(x) \langle k^*|) \quad (8.3) $$

In the diabatic basis over which the MPS dynamics is run, the projector can thus be written as:

$$ \hat{P}_{\text{AD}}^{\text{LOW}} = \sum_{n,m,i,j} |n\rangle \langle m| \otimes |i\rangle \langle j| \Delta_{nm}^{ij} \quad (8.4) $$

with

$$ \Delta_{nm}^{ij} = \int dx \phi_n(x) \phi_m^*(x) \chi_i^{\text{LOW}}(x) \chi_j^{\text{LOW}*}(x) \quad (9) $$

where $i$ and $j$ are either $e^*$ or $k^*$ while $n$ and $m$ represent the different Fock states of the RC harmonic oscillator.

**MPO formulation of the adiabatic projector.** From the MPO perspective, the projector can be seen as acting only on two sites: that of the electronic state $\langle e^*|$, $\langle k^*|$ and that of
the RC oscillator. The corresponding tensor therefore combines the electronic tensor with the RC mode tensor. For instance, if we focus on the projection onto the LOW surface with $\hat{P}_{AD}^{LOW}$, and define a discrete spatial grid $x = (x_1, x_2, \dots, x_N)$ and the electronic transition matrix $M^{LOW}$ for each discretized space point $x_i$ ($i \in [1, N]$):

$$ M^{LOW}(x_i) = \begin{pmatrix} \chi_{e^*}^{LOW}(x_i)\chi_{e^*}^{*,LOW}(x_i) & \chi_{e^*}^{LOW}(x_i)\chi_{k^*}^{*,LOW}(x_i) \\ \chi_{k^*}^{LOW}(x_i)\chi_{e^*}^{*,LOW}(x_i) & \chi_{k^*}^{LOW}(x_i)\chi_{k^*}^{*,LOW}(x_i) \end{pmatrix} \quad (10) $$

For a number $M$ of Fock states, we also define the $P(x_i)$ matrix from the products of eigenstates at fixed $x_i$:

$$ P(x_i) = \begin{pmatrix} \phi_1(x_i)\phi_1^*(x_i) & \dots & \phi_M(x_i)\phi_1^*(x_i) \\ \vdots & \ddots & \vdots \\ \phi_1(x_i)\phi_M^*(x_i) & \dots & \phi_M(x_i)\phi_M^*(x_i) \end{pmatrix} \quad (11) $$

Now, the full MPO projector can be written as:

$$ \hat{P}_{AD}^{LOW} = \underbrace{1 \otimes \dots \otimes 1}_{\text{Photons}} \otimes \underbrace{\left( M^{LOW}(x_1) \ M^{LOW}(x_2) \ \dots \ M^{LOW}(x_N) \right)}_{\text{S}} \otimes \underbrace{\left( \begin{array}{c} P(x_1) \\ P(x_2) \\ \vdots \\ P(x_N) \end{array} \right)}_{\text{RC mode}} \otimes \underbrace{1 \otimes \dots \otimes 1}_{\text{Vibrations}} \quad (12) $$

Expressed in this way, the projector combines the adiabatic electronic states obtained for each discrete position $x_i$. The dimension of the associated MPO should therefore scale linearly with the number of grid points $x_i$ as it is clear from eq 12. This should be numerically cumbersome, as a fine grid might be needed for an accurate description of the RC wavepacket. This manoeuvre would dramatically increase the bond dimension of the resulting wavefunction, as the product of a MPS with an MPO is an MPS whose bond dimension is the product of the original MPS' and the MPO's bond dimensions.⁸⁰ One way to circumvent this
unnecessary complexity is to reduce the dimension of the projector, by applying a Singular Value Decomposition (SVD), which reduces the projector operator to its necessary elements, while the result of its application to a wavefunction is not affected. More precisely, the bond dimension between the electronic MPO tensor and the RC mode MPO tensor is the number of discrete grid points $x_i$. This number can be large when Fock states of high index have to be taken into account. In the following examples, corresponding to typical ESIPT four-level systems, the number of $x$ points varies between 1300 and 2800. After SVD, the bond dimension of the projector is reduced to 5 for both the symmetric and non-symmetric double-well cases.

Once expressed in the MPO formalism and reduced using the SVD, the adiabatic projectors $\hat{P}_{AD}^{LOW}$ and $\hat{P}_{AD}^{UP}$ provide a useful and efficient tool to finely analyze the dynamics of the system, as illustrated by the examples of application in the following.

## 4.3 Coupling between Vibrations and the Reaction Coordinate

The MPS approach allows not only to simulate the quantum dynamics of the isolated {Electrons + Reaction Coordinate} system but also to efficiently describe its coupling with the vibrational and photonic environments. In particular, the vibrational environment dissipates energy away from the system, acting as a damping on the continuous RC coordinate that causes the wavepacket to relax to the lower energy states. In this case, vibrations interact only with the RC and do not directly couple $|e^*\rangle$ and $|k^*\rangle$ states, in contrast to the previous model described by equations (1.1) to (1.5).

In practice, the RC oscillator is coupled to a set of harmonic oscillators. The coupling Hamiltonian takes the characteristic form:

$$H_{\text{vib}} + H_{\text{int}}^{\text{vib-RC}} = \left( \sum_k \frac{p_k^2}{2m_k} + \sum_k \frac{1}{2} m_k \omega_k^2 (q_k - \alpha_k x)^2 \right) \quad (13.1)$$

with $q_k = \sqrt{\frac{1}{2m_k\omega_k}} (b_k + b_k^\dagger)$ the position operator of the $k^{th}$ vibration mode of frequency
$\omega_k$, mass $m_k$ and $p_k = i\sqrt{\frac{m_k\omega_k}{2}}(b_k^\dagger - b_k)$ its momentum. The interaction coefficients $\alpha_k$ are potentially different for each mode which enables representing arbitrary spectral densities for the vibration bath. To express this Hamiltonian in the MPS frame, it needs to be converted in second quantization:

$$H_{\text{vib}} = \sum_k \omega_k b_k^\dagger b_k \qquad (13.2)$$

$$H_{\text{int}}^{\text{vib-RC}} = \sum_k \frac{1}{2} m_k \omega_k^2 (-2 \alpha_k q_k + \alpha_k^2 x^2) \qquad (13.3)$$

$$= \sum_k \frac{\alpha_k^2 m_k \omega_k^2}{4 m_{\text{RC}} \omega_{\text{RC}}} (d + d^\dagger)^2 - (d + d^\dagger) \sum_k \frac{\alpha_k \omega_k \sqrt{m_k \omega_k}}{2 \sqrt{m_{\text{RC}} \omega_{\text{RC}}}} (b_k + b_k^\dagger) \qquad (13.4)$$

$$= \lambda_{\text{reorg}} (d + d^\dagger)^2 - (d + d^\dagger) \sum_k z_k (b_k + b_k^\dagger) \qquad (13.5)$$

The final interaction Hamiltonian takes the form a linear coupling between the RC mode and an environment of vibrations, weighted by the coefficients $z_k = \frac{\alpha_k \omega_k \sqrt{m_k \omega_k}}{2 \sqrt{m_{\text{RC}} \omega_{\text{RC}}}}$ corresponding to a spectral density $J(\omega) = \sum_k z_k^2 \delta(\omega_k - \omega)$. A supplementary term arises from the development, designated as the reorganisation energy $\lambda_{\text{reorg}} = \int \frac{J(\omega)}{\omega} d\omega$. This term compensates the modification of the double-well potential shape that would be induced in its absence by the vibrational environment.

This type of coupling leads to the star conformation of Figure 10. The electronic system is coupled to the RC mode which gives rise to a double-well potential energy surface. Environments include photons interacting with the electronic component of the system and molecular vibrations, coupled to the RC mode via equation (13.5). This picture also holds a physical meaning. The RC mode can be interpreted as one particular mode of the molecular vibrations that is strongly coupled to the electronic states. Therefore, this particular mode is treated separately from the other skeleton vibrations of the vibration environment in order to gather more information on its quantum state.
<figure><img src="image_11.png" /><figcaption>Figure 10: {Electrons + Reaction Coordinate} star conformation. The electronic system (orange) is interacting with *N* photon modes (red) characterised by the spectral density $J_{\text{phot}}(\omega) = \sum_{k=1}^{N} g_k^2 \delta(\omega_k - \omega)$. The system is also coupled to the RC mode, forming the double-well. The RC mode is itself interacting with *M* vibration modes (blue) characterised by the spectral density $J_{\text{vib}}(\omega) = \sum_{k=1}^{M} z_k^2 \delta(\omega_k - \omega)$.</figcaption></figure>

## 4.4 Reduced Density Matrix

To take full advantage of continuous description of the reaction coordinate, one can compute the corresponding reduced density matrix. It is obtained by tracing out the two environments in order to obtain a description of the quantum wavepacket corresponding specifically to the RC mode. It is also possible to project the reduced density matrix on one or the other of the adiabatic electronic states, to provide more insight into the adiabatic components of the wavepacket. These measurements allow a rich and visual study of the dynamics, going beyond mean values of observables. The full quantum state of the RC mode can be captured, including the spatial delocalisation of the transferring proton and the associated coherences (via the diagonal and off-diagonal elements of the reduced density matrix, respectively).

## 4.5 Dissipative Tunnelling in a Symmetric Double-Well Potential

To illustrate the different features presented in this section, a first elementary example describes a wavepacket that lies only onto the excited states. The wavepacket is initialised in the
left well of a symmetric excited adiabatic LOW surface (Figure 11), whereas the high energy of the excited UP surface prevents any non-adiabatic transitions in the dynamics. For the sake of the example we call the diabatic parameters $x_{RC}^{\text{eq,e*}} = x_{RC}^{\text{eq,LEFT}}$ and $x_{RC}^{\text{eq,k*}} = x_{RC}^{\text{eq,RIGHT}}$. Parameters are reported in Table S4 and Table S5 of the Supporting Information.

<figure><img src="image_12.png" /><figcaption>Figure 11: Wavepacket (schematised in black) initially localised onto the LOW left well. The chosen Hamiltonian parameters design the adiabatic LOW and UP surface as well as $x_{RC}^{\text{eq, LEFT}}$ and $x_{RC}^{\text{eq, RIGHT}}$. The initial energy of the isolated system is here below the barrier. RC value and $\omega$ are in atomic Hartree units.</figcaption></figure>

Figure 12 shows the dynamics of the isolated system ((a),(b)) and of the system damped by vibrations ((c),(d)). Concerning the isolated system, Figure 12 (a) shows time snapshots of the wavepacket amplitude along the Reaction Coordinate. We can see that the wavepacket is travelling back and forth between the left well and the right well. The mean value of the displacement operator $X$ is also represented in Figure 12 (b) showing oscillations between $\langle X \rangle = x_{RC}^{\text{eq,LEFT}} = -0.25$ a.u. and $\langle X \rangle = x_{RC}^{\text{eq,RIGHT}} = 0.25$ a.u. with smaller oscillations.
Once vibrations are introduced, wavepacket time snapshots (Figure 12 (c)) show a stabilisa-
tion with the wavepacket equally split between the two wells. The ⟨X⟩ dynamics (Figure 12
(d)) is also modified, converging towards ⟨X⟩ = 0 a.u..

Regarding the isolated case, although the initial state has an energy lower than the
potential barrier, the quantum description of the RC mode enables tunnelling to occur.
The adiabatic projector used as a post-treatment tool allows to distinguish between the
wavepacket LOW and a UP components (the latter being essentially zero in this example).
The small oscillations in the ⟨X⟩ dynamics are due to the wavepacket initial position that is
not exactly the equilibrium position of the adiabatic surface. For the case of the system that
is damped by vibrations, the wavepacket can only be represented by the reduced density ma-
trix ρ<sub>RC</sub>(x, x'). The population amplitude in function of the RC value is obtained by taking
the diagonal elements of the reduced density matrix. To obtain the UP and LOW parts, the
reduced density matrix is projected on the fly aside from the dynamics and results in the
reduced density matrix UP ρ<sub>RC</sub><sup>UP</sup>(x, x') and LOW ρ<sub>RC</sub><sup>LOW</sup>(x, x'). Despite the tunnelling process,
the full transition from left to right does not happen here and the system wavepacket loses
energy due to vibrations. Therefore, it stabilises towards the ground state which is charac-
terised by the wavepacket uniformly split between the two wells and a mean displacement
value ⟨X⟩ = 0 a.u..
<figure><img src="image_13.png" /><figcaption>Figure 12: Dynamics of a system initialised in the left well. (a) Time snapshots of the isolated system wavepacket. The LOW and UP wavepacket parts are obtained with the post-analysis adiabatic projector. The tunnelling effect allows passage between the left and the right well. (b) Mean value of the X operator for the isolated system. Oscillations between minima of the left and right wells are shown. (c) Time snapshots of the adiabatic diagonal elements $||\rho_{RC}^{\text{LOW/UP}}(x, x)||$. Vibrations stabilise the system towards the ground state. (d) Mean value of the X operator for the damped system. Vibrations make it converge towards $\langle X \rangle = 0$ a.u.. Time is in atomic Hartree units.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13: Scheme of the energetic landscape of a four-level system described with the RC dimension. Diabatic $|e^*\rangle$ and $|k^*\rangle$ minima are separated by a barrier. The vertical electronic gaps are kept the same as for discretized states. The wavepacket is initialised in $|e\rangle$ and promoted to $|e^*\rangle$ by a pulse $H_{\text{drive}}(t)$.</figcaption></figure>
## 4.6 Excited-State Proton Transfer under Driving: Exact Treatment of Dissipative Non-Adiabatic Dynamics

The second example takes the four-level model of section 3.4 further by incorporating the RC dimension to it (Figure 13). As in Figure 9, the two excited discrete states are directly coupled and become an unique adiabatic surface with its two local minima separated by a barrier. The wavepacket is initialised in the enol ground state and is promoted by a laser pulse of the form:

$$H_{\text{drive}}(t) = (|e^*\rangle\langle e| + |k^*\rangle\langle k| + \text{h.c.}) \epsilon \exp((t-t_0)^2 / (2\tau_{\text{drive}}^2)) \cos(\omega_{\text{drive}}t) \quad (14)$$

with $t_0$ the centre time of the pulse and $\tau_{\text{drive}} = \text{FWHM}/2\sqrt{2\log(2)}$ with FWHM the full width at half maximum. This time-dependent Hamiltonian simulates a laser pulse with a gaussian envelope. The two previous environments (photons and vibrations) are also included in the full wave-function. Parameters are reported in Table S6 of the Supporting Information, their values are chosen to represent a prototypical ESIPT process. In particular, the RC frequency is typical of a stretching mode involving hydrogen.

From these dynamics, diabatic and adiabatic populations are probed. Figure 14 shows the initial population in $|e\rangle$ that oscillates due to the pulse. After the pulse, the electronic population is mainly in $|e^*\rangle$ and $|e\rangle$ with a little bit of $|k^*\rangle$. For the rest of the dynamics, $|e\rangle$ population seems constant while a transfer between $|e^*\rangle$ and $|k^*\rangle$ occurs. After the pulse the excited state population is principally in the LOW excited adiabatic state. Two peaks are detected in the time-resolved photon spectrum at frequency $\Delta\omega_e$ and $\Delta\omega_k$ (Figure 15 (a)), whereas that in the vibration spectrum, the main peak is at $\omega_{\text{RC}}$ but other less energetic excitations can be detected (Figure 15 (b)).

Unlike the previous example of a four-level under strong driving, the gaussian shape of the pulse does not induce dressed states in the dynamics. The pulse populates $|e^*\rangle$ and vibrations induce transition from $|e^*\rangle$ to $|k^*\rangle$, without involving non-adiabatic transitions.
The damping induced by both environments will make the wavepacket to fall down in en-
ergetic minima of the excited manifold while simultaneously slightly drop down towards
ground states by emitting a photon. Although the main vibrational peak is at ωRC, several
frequencies can be detected in the vibration spectrum. In fact, the coupling between the
two wells allows diverse gap energies between vibrational states, broadening the detected
vibrational frequency range around ωRC. A low-frequency feature is also visible due to the
bath displacement. It is important to point out that, although the space dimension has been
added in this example, the dual fluorescence in the photo-emission spectrum is recovered.
This model shows the same main characteristics of the discretized four-level system but with
more information concerning the proton wavepacket.

At this stage, one can really appreciate the potential of the method and the new features.
As a matter of fact, by combining the RC mode with the time-dependent Hamiltonian, the
problematic question in the initialisation of the wavepacket is bypassed. This specific issue
can represent a bottleneck for the use of numerical method such as Surface Hopping when
the Franck-Condon region includes several vibrational modes or in order to include the initial
energy distribution.<sup>106,107</sup>

Dynamics with different parameters were conducted to illustrate non-adiabatic transi-
tions. For the example showed in Figure 16, the system is the same as the previous example
except that the enol ground state is shifted from the enol excited state and the driving has
a higher frequency, leading to a higher energy for the promoted wavepacket. Parameters are
reported in Table S7 of the Supporting Information.

The time evolution of the different states is reported in Figure 17, showing a behaviour
similar to the previous example. The main difference is the UP population that is exhibiting
periodical peaks at short times. With the same period, the $|e^*\rangle$ population undergoes a
<figure><img src="image_15.png" /><figcaption>Figure 14: Diabatic and excited adiabatic population dynamics. The adiabatic populations are calculated with the adiabatic projector. After being populated by the pulse, coupling with vibrations transfer $|e^*\rangle$ population to $|k^*\rangle$. Time is in atomic Hartree units.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 15: (a) Photo-emission spectrum resolved in time. Two peaks are detected at $\Delta\omega_e$ and at $\Delta\omega_k$. (b) Vibration spectrum resolved in time. The main peak is at $\omega_{RC}$ but smaller frequencies are also probed. This frequency range is allowed by the anharmonicity of the double-well. Time and frequencies are in atomic Hartree units.</figcaption></figure>
steep decrease and is transferred to $|k^*\rangle$. These adiabatic population variations highlight non-adiabatic transitions, whereas the steps between diabatic population suggest a tunnelling process. Despite non-adiabatic transitions at short times, the excited-state electronic transfer at longer times is similar to the previous example.

In Figure 18, the reduced density matrix is represented and illustrates a useful tool here. This observable in space gives insights into the position onto the excited states. The diagonal elements of this observable represent the population in space whereas anti-diagonal elements represent coherences. The investigation of the adiabatic reduced density matrix for the excited states allows to observe the oscillation of the wavepacket onto the left well, transferring its weight little by little towards the right well. This space monitoring along the RC helps detecting that non-adiabatic transitions happen when the left wavepacket oscillation gets closer to the barrier at time around $t = 1390$ a.u. and $t = 2110$ a.u. for instance. The frequency of these wavepacket oscillations corresponds to the RC mode frequency $\omega_{RC} = 2\pi/T$ with $T \approx 650$ a.u.. Once vibration dissipation takes effect and stabilises the main wavepacket onto the bottom of the left surface (from $t = 4000$ a.u.), an adiabatic dynamics is shown for the transfer from $|e^*\rangle$ to $|k^*\rangle$, where the excited keto wavepacket is growing at the expense of the excited enol wavepacket. A last interesting observation concerns the anti-diagonal elements of the reduced density matrix. In fact, the space coherences in $||\rho_{RC}^{LOW}(x, x')||$ are probed especially in short times. These features are characteristic of quantum out-of-equilibrium dynamics and they stay significant until the dissipation due to vibrations kills them. These quantum coherences in space are a result of the Hamiltonian constructed with electronic states and a RC quantum harmonic oscillator.
<figure><img src="image_17.png" /><figcaption>Figure 16: Scheme of the energetic landscape of a four-level system described with the RC dimension. The wavepacket is initialised in $|e\rangle$ that is shifted from $|e^*\rangle$ and promoted by a pulse $H_{\text{drive}}(t)$. The promoted wavepacket has a high energy and undergoes non-adiabatic transitions.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 17: Diabatic and excited adiabatic population dynamics. Adiabatic populations are calculated with the adiabatic projector and illustrate non-adiabatic transitions. After being populated by the pulse, $|e^*\rangle$ transfers its population to $|k^*\rangle$ thanks to the coupling with vibrations. Non-adiabatic transitions occur periodically in the short times. Time is in atomic Hartree units.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 18: Time snapshots of the adiabatic reduced density matrix $||\rho_{RC}^{LOW}(x, x')||$. Equilibrium positions of the diabatic excited states are represented by the black crosses. Non-adiabatic transitions occur when the wavepacket becomes close to the barrier. *x* and *t* are in atomic Hartree units.</figcaption></figure>
# 5 Conclusion

In this work, we present two models to tackle the problem of proton transfer with environment interactions using Matrix Product State. A first model represents a four-level electronic system under driving that incorporates intramolecular vibration and photon interactions. Secondly, we extended the TEDOPA algorithm to explicitly include a continuous reaction-coordinate within the system, leading to a real-space double-well description of ESIPT which we studied at a quantum level. This description allows to study the electronic transitions as well as the proton displacement on energy landscapes with a full treatment of non-adiabatic effects. Developing a numerically efficient MPO projector, the resulting dissipative dynamics were analysed in both diabatic and adiabatic bases, demonstrating that the method is an efficient and complete tool to study non-adiabatic processes for open quantum systems. This approach can naturally be extended for further study with thermalised (finite-temperature) environments that could lead to interesting phenomena as shown in Ref. ¹⁰⁸, as well as realistic, structured spectral densities based on *ab initio* calculations, following the ideas recently presented in Refs. ⁵⁹,⁹⁵ Although the energy surfaces presented here were generic toy models, more specific energy landscapes such as a Marcus inverted region, different RC frequencies for $|e^*\rangle$ and $|k^*\rangle$, or multiple surfaces extracted from experimental and/or *ab initio* data could be studied with this methodology, and would provide interesting outlooks. Future work in our team aims to consolidate these developments to model ESIPT in real molecules where ultrafast experimental data is available.

# Acknowledgement

We thank iSiM (Initiative Sciences et ingénierie moléculaires) from the Alliance Sorbonne Université for funding.
## Supporting Information Available

Additional information about the MPO representation of the models. Analytical results and additional observables for the second case of the four-level system under driving. This information is available free of charge via the Internet at http://pubs.acs.org

## References

(1) Reece, S. Y.; Nocera, D. G. Proton-Coupled Electron Transfer in Biology: Results from Synergistic Studies in Natural and Model Systems. *Ann. Rev. Biochem.* **2009**, *78*, 673–699.

(2) Lama, B.; Sarma, M. Unraveling the Mechanistic Pathway for the Dual Fluorescence in Green Fluorescent Protein (GFP) Chromophore Analogue: A Detailed Theoretical Investigation. *J. Phys. Chem. B* **2022**, *126*, 9930–9944.

(3) Ren, G.; Meng, Q.; Zhao, J.; Chu, T. Molecular design for electron-driven double-proton transfer: A new scenario for excited-state proton-coupled electron transfer. *J. Phys. Chem. A* **2018**, *122*, 9191–9198.

(4) Lennox, J. C.; Kurtz, D. A.; Huang, T.; Dempsey, J. L. Excited-state proton-coupled electron transfer: different avenues for promoting proton/electron movement with solar photons. *ACS Energy Lett.* **2017**, *2*, 1246–1256.

(5) Joshi, H. C.; Antonov, L. Excited-State Intramolecular Proton Transfer: A Short Introductory Review. *Molecules* **2021**, *26*, 1475.

(6) Chen, L.; Fu, P.; Wang, H.; Pan, M. Excited-State Intramolecular Proton Transfer (ESIPT) for Optical Sensing in Solid State. *Adv. Opt. Mater.* **2021**, *9*, 2001952.

(7) Sedgwick, A. C.; Wu, L.; Han, H.-H.; Bull, S. D.; He, X.-P.; James, T. D.; Sessler, J. L.;
Tang, B. Z.; Tian, H.; Yoon, J. Excited-state intramolecular proton-transfer (ESIPT) based fluorescence sensors and imaging agents. *Chem. Soc. Rev.* **2018**, *47*, 8842–8880.

(8) Li, Y.; Dahal, D.; Abeywickrama, C. S.; Pang, Y. Progress in Tuning Emission of the Excited-State Intramolecular Proton Transfer (ESIPT)-Based Fluorescent Probes. *ACS Omega* **2021**, *6*, 6547–6553.

(9) Reece, S. Y.; Nocera, D. G. Proton-Coupled Electron Transfer in Biology: Results from Synergistic Studies in Natural and Model Systems. *Ann. Rev. Biochem.* **2009**, *78*, 673–699.

(10) Santos, E. M.; Sheng, W.; Esmatpour Salmani, R.; Tahmasebi Nick, S.; Ghanbar-pour, A.; Gholami, H.; Vasileiou, C.; Geiger, J. H.; Borhan, B. Design of large stokes shift fluorescent proteins based on excited state proton transfer of an engineered photobase. *J. Am. Chem. Soc.* **2021**, *143*, 15091–15102.

(11) Goyal, P.; Hammes-Schiffer, S. Tuning the Ultrafast Dynamics of Photoinduced Proton-Coupled Electron Transfer in Energy Conversion Processes. *ACS Energy Lett.* **2017**, *2*, 512–519.

(12) Liang, M.; Yan, X. Nanozymes: From New Concepts, Mechanisms, and Standards to Applications. *Acc. Chem. Res.* **2019**, *52*, 2190–2200.

(13) Sokolovskii, I.; Groenhof, G. Photochemical initiation of polariton-mediated exciton propagation. *Nanophotonics* **2024**, *13*, 2687–2694.

(14) Jankowska, J.; Sobolewski, A. L. Modern Theoretical Approaches to Modeling the Excited-State Intramolecular Proton Transfer: An Overview. *Molecules* **2021**, *26*, 5140.

(15) Hammes-Schiffer, S. Proton-coupled electron transfer: classification scheme and guide to theoretical methods. *Energy Environ. Sci.* **2012**, *5*, 7696.
(16) Pijeau, S.; Foster, D.; Hohenstein, E. G. Excited-State Dynamics of 2-(2'-Hydroxyphenyl)benzothiazole: Ultrafast Proton Transfer and Internal Conversion. *J. Phys. Chem. A* **2017**, *121*, 4595-4605.

(17) Kim, J.; Kim, C. H.; Burger, C.; Park, M.; Kling, M. F.; Kim, D. E.; Joo, T. Non-Born-Oppenheimer Molecular Dynamics Observed by Coherent Nuclear Wave Packets. *J. Phys. Chem. Lett.* **2020**, *11*, 755-761.

(18) Picconi, D. Nonadiabatic quantum dynamics of the coherent excited state intramolecular proton transfer of 10-hydroxybenzo[h]quinoline. *Photochem. Photobiol. Sci.* **2021**, *20*, 1455-1473.

(19) Xu, J.; Zhou, R.; Blum, V.; Li, T. E.; Hammes-Schiffer, S.; Kanai, Y. First-Principles Approach for Coupled Quantum Dynamics of Electrons and Protons in Heterogeneous Systems. *Phys. Rev. Lett.* **2023**, *131*, 238002.

(20) Schnorr, K.; Belina, M.; Augustin, S.; Lindenblatt, H.; Liu, Y.; Meister, S.; Pfeifer, T.; Schmid, G.; Treusch, R.; Trost, F.; Slavíek, P.; Moshammer, R. Direct tracking of ultrafast proton transfer in water dimers. *Sci. Adv.* **2023**, *9*, eadg7864.

(21) Weinberg, D. R.; Gagliardi, C. J.; Hull, J. F.; Murphy, C. F.; Kent, C. A.; Westlake, B. C.; Paul, A.; Ess, D. H.; McCafferty, D. G.; Meyer, T. J. Proton-coupled electron transfer. *Chem. Rev.* **2012**, *112*, 4016-4093.

(22) Wendling, M.; Pullerits, T.; Przyjalgowski, M. A.; Vulto, S. I. E.; Aartsma, T. J.; van Grondelle, R.; van Amerongen, H. Electron Vibrational Coupling in the Fenna Matthews Olson Complex of Prosthecochloris aestuarii Determined by Temperature-Dependent Absorption and Fluorescence Line Narrowing Measurements. *J. Phys. Chem. B* **2000**, *104*, 5825-5831.

(23) Reimers, J. R.; Cai, Z.-L.; Kobayashi, R.; Rätsep, M.; Freiberg, A.; Krausz, E. As-
segment of the Q-Bands of the Chlorophylls: Coherence Loss via Qx Qy Mixing.
*Scientific Reports* **2013**, *3*.

(24) Pachon, L. A.; Brumer, P. Direct experimental determination of spectral densities of molecular complexes. *J. Chem. Phys.* **2014**, *141*.

(25) Lee, J.; Kim, C. H.; Joo, T. Active Role of Proton in Excited State Intramolecular Proton Transfer Reaction. *J. Phys. Chem. A* **2013**, *117*, 1400–1405.

(26) Schriever, C.; Barbatti, M.; Stock, K.; Aquino, A. J.; Tunega, D.; Lochbrunner, S.; Riedle, E.; De Vivie-Riedle, R.; Lischka, H. The interplay of skeletal deformations and ultrafast excited-state intramolecular proton transfer: Experimental and theoretical investigation of 10-hydroxybenzo[h]quinoline. *Chem. Phys.* **2008**, *347*, 446–461.

(27) Kim, C. H.; Joo, T. Coherent excited state intramolecular proton transfer probed by time-resolved fluorescence. *Phys. Chem. Chem. Phys.* **2009**, *11*, 10266.

(28) Lochbrunner, S.; Wurzer, A. J.; Riedle, E. Microscopic Mechanism of Ultrafast Excited-State Intramolecular Proton Transfer: A 30-fs Study of 2-(2'-Hydroxyphenyl)benzothiazole. *J. Phys. Chem. A* **2003**, *107*, 10580–10590.

(29) Takeuchi, S.; Tahara, T. Coherent Nuclear Wavepacket Motions in Ultrafast Excited-State Intramolecular Proton Transfer: Sub-30-fs Resolved Pump-Probe Absorption Spectroscopy of 10-Hydroxybenzo[h]quinoline in Solution. *J. Phys. Chem. A* **2005**, *109*, 10199–10207.

(30) Wang, C.-H.; Liu, Z.-Y.; Huang, C.-H.; Chen, C.-T.; Meng, F.-Y.; Liao, Y.-C.; Liu, Y.-H.; Chang, C.-C.; Li, E. Y.; Chou, P.-T. Chapter Open for the Excited-State Intramolecular Thiol Proton Transfer in the Room-Temperature Solution. *J. Am. Chem. Soc.* **2021**, *143*, 12715–12724.
(31) Yin, Z.; Chang, Y.-P.; Balčiūnas, T.; Shakya, Y.; Djorović, A.; Gaulier, G.; Fazio, G.; Santra, R.; Inhester, L.; Wolf, J.-P.; Wörner, H. J. Femtosecond proton transfer in urea solutions probed by X-ray spectroscopy. *Nature* **2023**, *619*, 749–754.

(32) Zhao, L.; Wildman, A.; Pavošević, F.; Tully, J. C.; Hammes-Schiffer, S.; Li, X. Excited State Intramolecular Proton Transfer with Nuclear-Electronic Orbital Ehrenfest Dynamics. *J. Phys. Chem. Lett.* **2021**, *12*, 3497–3502.

(33) Zhong, J.; Soudackov, A. V.; Hammes-Schiffer, S. Probing Nonadiabaticity of Proton-Coupled Electron Transfer in Ribonucleotide Reductase. *J. Phys. Chem. Lett.* **2024**, 1686–1693.

(34) Yu, Q.; Roy, S.; Hammes-Schiffer, S. Nonadiabatic Dynamics of Hydrogen Tunneling with Nuclear-Electronic Orbital Multistate Density Functional Theory. *J. Chem. Theory Comput.* **2022**, acs.jctc.2c00938.

(35) Chow, M.; Li, T. E.; Hammes-Schiffer, S. Nuclear–Electronic Orbital Quantum Mechanical/Molecular Mechanical Real-Time Dynamics. *J. Phys. Chem. Lett.* **2023**, *14*, 9556–9562.

(36) Toldo, J. M.; do Casal, M. T.; Ventura, E.; do Monte, S. A.; Barbatti, M. Surface hopping modeling of charge and energy transfer in active environments. *Phys. Chem. Chem. Phys.* **2023**, *25*, 8293–8316.

(37) Agostini, F.; Curchod, B. F. E. Different flavors of nonadiabatic molecular dynamics. *WIREs Comput. Mol. Sci.* **2019**, *9*, e1417.

(38) Meyer, H.-D.; Manthe, U.; Cederbaum, L. The multi-configurational time-dependent Hartree approach. *Chem. Phys. Lett.* **1990**, *165*, 73–78.

(39) Meyer, H.; Gatti, F.; Worth, G. *Multidimensional Quantum Dynamics: MCTDH Theory and Applications*; Wiley, 2009.
(40) Bowman, J. M.; Carrington, T.; Meyer, H.-D. Variational quantum approaches for computing vibrational energies of polyatomic molecules. *Molecular Physics* **2008**, *106*, 2145-2182.

(41) Meyer, H. Studying molecular quantum dynamics with the multiconfiguration time-dependent Hartree method. *WIREs Computational Molecular Science* **2011**, *2*, 351-374.

(42) Manthe, U. The multi-configurational time-dependent Hartree approach revisited. *J. Chem. Phys.* **2015**, *142*.

(43) Popp, W.; Brey, D.; Binder, R.; Burghardt, I. Quantum Dynamics of Exciton Transport and Dissociation in Multichromophoric Systems. *Annu. Rev. Phys. Chem.* **2021**, *72*, 591-616.

(44) Di Maiolo, F.; Worth, G. A.; Burghardt, I. Multi-layer Gaussian-based multi-configuration time-dependent Hartree (ML-GMCTDH) simulations of ultrafast charge separation in a donor-acceptor complex. *J. Chem. Phys.* **2021**, *154*.

(45) Han, S.; Schröder, M.; Gatti, F.; Meyer, H.-D.; Lauvergnat, D.; Yarkony, D. R.; Guo, H. Representation of Diabatic Potential Energy Matrices for Multiconfiguration Time-Dependent Hartree Treatments of High-Dimensional Nonadiabatic Photodissociation Dynamics. *J. Chem. Theory Comput.* **2022**, *18*, 4627-4638.

(46) Perveaux, A.; Lorphelin, M.; Lasorne, B.; Lauvergnat, D. Fast and slow excited-state intramolecular proton transfer in 3-hydroxychromone: a two-state story? *Phys. Chem. Chem. Phys.* **2017**, *19*, 6579-6593.

(47) Anand, N.; Isukapalli, S. V. K.; Vennapusa, S. R. Excited-state intramolecular proton transfer driven by conical intersection in hydroxychromones. *J. Comput. Chem.* **2020**, *41*, 1068-1080.
(48) Craig, I. R.; Thoss, M.; Wang, H. Proton transfer reactions in model condensed-phase environments: Accurate quantum dynamics using the multilayer multiconfiguration time-dependent Hartree approach. *J. Chem. Phys.* **2007**, *127*.

(49) Richings, G.; Polyak, I.; Spinlove, K.; Worth, G.; Burghardt, I.; Lasorne, B. Quantum dynamics simulations using Gaussian wavepackets: the vMCG method. *International Reviews in Physical Chemistry* **2015**, *34*, 269–308.

(50) Craig, I. R.; Thoss, M.; Wang, H. Accurate quantum-mechanical rate constants for a linear response Azzouz-Borgis proton transfer model employing the multilayer multi-configuration time-dependent Hartree approach. *J. Chem. Phys.* **2011**, *135*.

(51) Breuer, H.-P.; Petruccione, F. *The theory of open quantum systems*; OUP Oxford, 2002.

(52) Zhang, L.; Fassioli, F.; Fu, B.; She, Z.-S.; Scholes, G. D. Modeling Excited-State Proton Transfer Using the Lindblad Equation: Quantification of Time-Resolved Spectroscopy with Mechanistic Insights. *ACS Phys. Chem. Au* **2022**, *3*, 107–118.

(53) Prior, J.; Chin, A. W.; Huelga, S. F.; Plenio, M. B. Efficient Simulation of Strong System-Environment Interactions. *Phys. Rev. Lett.* **2010**, *105*, 050404.

(54) Prior, J.; de Vega, I.; Chin, A. W.; Huelga, S. F.; Plenio, M. B. Quantum dynamics in photonic crystals. *Phys. Rev. A* **2013**, *87*, 013428.

(55) Chin, A. W.; Prior, J.; Rosenbach, R.; Caycedo-Soler, F.; Huelga, S. F.; Plenio, M. B. The role of non-equilibrium vibrational structures in electronic coherence and coherence in pigment-protein complexes. *Nature Phys.* **2013**, *9*, 113–118.

(56) Del Pino, J.; Schröder, F. A.; Chin, A. W.; Feist, J.; Garcia-Vidal, F. J. Tensor network simulation of polaron-polaritons in organic microcavities. *Phys. Rev. B* **2018**, *98*, 165416.
(57) Schröder, F. A.; Turban, D. H.; Musser, A. J.; Hine, N. D.; Chin, A. W. Tensor network simulation of multi-environmental open quantum dynamics via machine learning and entanglement renormalisation. *Nature Commun.* **2019**, *10*, 1062.

(58) Dunnett, A. J.; Gowland, D.; Isborn, C. M.; Chin, A. W.; Zuehlsdorff, T. J. Influence of non-adiabatic effects on linear absorption spectra in the condensed phase: Methylene blue. *J. Chem. Phys.* **2021**, *155*.

(59) Hunter, K. E.; Mao, Y.; Chin, A. W.; Zuehlsdorff, T. J. Environmentally Driven Symmetry Breaking Quenches Dual Fluorescence in Proflavine. *J. Phys. Chem. Lett.* **2024**, *15*, 4623-4632.

(60) Bañuls, M. C. Tensor Network Algorithms: a Route Map. *Annu. Rev. Condens. Matter Phys.* **2023**, *14*, 173-191, arXiv:2205.10345 [quant-ph].

(61) Larsson, H. R. A tensor network view of multilayer multiconfiguration time-dependent Hartree methods. *Mol. Phys.* **2024**, e2306881.

(62) Hughes, K. H.; Christ, C. D.; Burghardt, I. Effective-mode representation of non-Markovian dynamics: A hierarchical approximation of the spectral density. II. Application to environment-induced nonadiabatic dynamics. *J. Chem. Phys.* **2009**, *131*.

(63) Zeng, T.; Fedorov, D. G.; Klobukowski, M. Model core potentials for studies of scalar-relativistic effects and spin-orbit coupling at Douglas-Kroll level. I. Theory and applications to Pb and Bi. *J. Chem. Phys.* **2009**, *131*.

(64) Martinazzo, R.; Vacchini, B.; Hughes, K. H.; Burghardt, I. Communication: Universal Markovian reduction of Brownian particle dynamics. *J. Chem. Phys.* **2011**, *134*.

(65) Martinazzo, R.; Hughes, K. H.; Burghardt, I. Unraveling a Brownian particle’s memory with effective mode chains. *Phys. Rev. E* **2011**, *84*.
(66) Popp, W.; Polkehn, M.; Hughes, K. H.; Martinazzo, R.; Burghardt, I. Vibronic cou-
pling models for donor-acceptor aggregates using an effective-mode scheme: Applica-
tion to mixed Frenkel and charge-transfer excitons in oligothiophene aggregates. *J.
Chem. Phys.* **2019**, *150*.

(67) Fischer, E. W.; Werther, M.; Bouakline, F.; Grossmann, F.; Saalfrank, P. Non-Markovian vibrational relaxation dynamics at surfaces. *J. Chem. Phys.* **2022**, *156*.

(68) Lacroix, T.; Lovett, B. W.; Chin, A. W. From non-Markovian dissipation to spa-
tiotemporal control of quantum nanodevices. *Quantum* **2024**, *8*, 1305.

(69) Oseledets, I. V. Tensor-train decomposition. *SIAM J. Sci. Comput.* **2011**, *33*, 2295–2317.

(70) Lubich, C.; Oseledets, I. V.; Vandereycken, B. Time integration of tensor trains. *SIAM J. Numer. Anal.* **2015**, *53*, 917–941.

(71) Borrelli, R.; Gelin, M. F. Finite temperature quantum dynamics of complex systems: Integrating thermo-field theories and tensor-train methods. *Wiley Interdiscip. Rev.: Comput. Mol. Sci.* **2021**, *11*, e1539.

(72) Paeckel, S.; Köhler, T.; Swoboda, A.; Manmana, S. R.; Schollwöck, U.; Hubig, C. Time-evolution methods for matrix-product states. *Ann. Phys.* **2019**, *411*, 167998, arXiv:1901.05824 [cond-mat, physics:quant-ph].

(73) Lyu, N.; Soley, M. B.; Batista, V. S. Tensor-train split-operator KSL (TT-SOKSL) method for quantum dynamics simulations. *J. Chem. Theory Comput.* **2022**, *18*, 3327–3346.

(74) Peng, J.; Hu, D.; Liu, H.; Shi, Q.; Bao, P.; Lan, Z. Studies of nonadiabatic dynamics in the singlet fission processes of pentacene dimer via tensor network method. *J. Chem. Phys.* **2023**, *159*.
(75) Greene, S. M.; Batista, V. S. Tensor-train split-operator fourier transform (TT-SOFT) method: Multidimensional nonadiabatic quantum dynamics. *J. Chem. Theory Comput.* **2017**, *13*, 4034–4042.

(76) Mangaud, E.; Jaouadi, A.; Chin, A.; Desouter-Lecomte, M. Survey of the hierarchical equations of motion in tensor-train format for non-Markovian quantum dynamics. *EPJ ST* **2023**, *232*, 1847–1869.

(77) Shu, Y.; Truhlar, D. G. Decoherence and its role in electronically nonadiabatic dynamics. *J. Chem. Theory Comput.* **2023**, *19*, 380–395.

(78) Cirac, J. I.; Pérez-García, D.; Schuch, N.; Verstraete, F. Matrix product states and projected entangled pair states: Concepts, symmetries, theorems. *Rev. Mod. Phys.* **2021**, *93*, 045003.

(79) Orús, R. A practical introduction to tensor networks: Matrix product states and projected entangled pair states. *Ann. Phys.* **2014**, *349*, 117–158.

(80) Schollwöck, U. The density-matrix renormalization group in the age of matrix product states. *Ann. Phys.* **2011**, *326*, 96–192.

(81) Tamascelli, D. Excitation Dynamics in Chain-Mapped Environments. *Entropy* **2020**, *22*.

(82) Chin, A. W.; Rivas, A.; Huelga, S. F.; Plenio, M. B. Exact mapping between system-reservoir quantum models and semi-infinite discrete chains using orthogonal polynomials. *J. Math. Phys.* **2010**, *51*, 092109.

(83) Strathearn, A.; Kirton, P.; Kilda, D.; Keeling, J.; Lovett, B. W. Efficient non-Markovian quantum dynamics using time-evolving matrix product operators. *Nature Commun.* **2018**, *9*, 3322.
(84) Haegeman, J.; Lubich, C.; Oseledets, I.; Vandereycken, B.; Verstraete, F. Unifying time evolution and optimization with matrix product states. *Phys. Rev. B* **2016**, *94*, 165116.

(85) Dunnett, A. J.; Chin, A. W. Efficient bond-adaptive approach for finite-temperature open quantum dynamics using the one-site time-dependent variational principle for matrix product states. *Phys. Rev. B* **2021**, *104*, 214302.

(86) García-Ripoll, J. J. Time evolution of Matrix Product States. *New J. Phys.* **2006**, *8*, 305-305.

(87) Yang, M.; White, S. R. Time-dependent variational principle with ancillary Krylov subspace. *Phys. Rev. B* **2020**, *102*, 094315.

(88) Rams, M. M.; Zwolak, M. Breaking the entanglement barrier: Tensor network simulation of quantum transport. *Phys. Rev. Lett.* **2020**, *124*, 137701.

(89) Dunnett, A.; Lacroix, T.; Le Dé, B.; Riva, A. shareloqs/MPSDynamics. 2021; https://zenodo.org/record/5106435.

(90) Lacroix, T.; Le Dé, B.; Riva, A.; Dunnett, A. J.; Chin, A. W. MPSDynamics.jl: Tensor network simulations for finite-temperature (non-Markovian) open quantum system dynamics. *J. Chem. Phys.* **2024**, *161*, 084116.

(91) Csehi, A.; Kowalewski, M.; Halász, G. J.; Vibók, A. Ultrafast dynamics in the vicinity of quantum light-induced conical intersections. *New Journal of Physics* **2019**, *21*, 093040.

(92) Triana, J. F.; Herrera, F. Ultrafast modulation of vibrational polaritons for controlling the quantum field statistics at mid-infrared frequencies. *New J. Phys.* **2022**,
(93) Lorenzoni, N.; Cho, N.; Lim, J.; Tamascelli, D.; Huelga, S. F.; Plenio, M. B. Systematic Coarse Graining of Environments for the Nonperturbative Simulation of Open Quantum Systems. *Phys. Rev. Lett.* **2024**, 132.

(94) Schröder, F. A.; Chin, A. W. Simulating open quantum dynamics with time-dependent variational matrix product states: Towards microscopic correlation of environment dynamics and reduced system evolution. *Phys. Rev. B* **2016**, 93, 075105.

(95) Riva, A.; Tamascelli, D.; Dunnett, A. J.; Chin, A. W. Thermal cycle and polaron formation in structured bosonic environments. *Phys. Rev. B* **2023**, 108, 195138.

(96) Mollow, B. R. Power Spectrum of Light Scattered by Two-Level Systems. *Phys. Rev.* **1969**, 188, 1969-1975.

(97) Harshan, A. K.; Yu, T.; Soudackov, A. V.; Hammes-Schiffer, S. Dependence of Vibronic Coupling on Molecular Geometry and Environment: Bridging Hydrogen Atom Transfer and Electron-Proton Transfer. *J. Am. Chem. Soc.* **2015**, 137, 13545-13555.

(98) Förster, T. Diabatic and adiabatic processes in photochemistry. *Pure Appl. Chem.* **1970**, 24, 443-450.

(99) Turro, N. J.; McVey, J.; Ramamurthy, V.; Lechtken, P. Adiabatic Photoreactions of Organic Molecules. *Angew. Chem. Int. Ed.* **1979**, 18, 572-586.

(100) Worth, G. A.; Beck, M. H.; Jäckle, A.; Meyer, H.-D.; Vendrell, O. The MCTDH Package, Version 8.2 (2000), Version 8.3 (2002), Version 8.4 (2007), Version 8.5 (2013), Version 8.6.2 (2022). http://mctdh.uni-hd.de, 2022; Versions 8.5 and 8.6 contain the ML-MCTDH algorithm.

(101) Worth, G. Quantics: A general purpose package for Quantum molecular dynamics simulations. *Comput. Phys. Commun.* **2020**, 248, 107040.
(102) Worth, G. A.; Meyer, H.-D.; Cederbaum, L. S. Relaxation of a system with a conical intersection coupled to a bath: A benchmark 24-dimensional wave packet study treating the environment explicitly. *J. Chem. Phys.* **1998**, *109*, 3518–3529.

(103) Worth, G. A.; Meyer, H.-D.; Köppel, H.; Cederbaum, L. S.; Burghardt, I. Using the MCTDH wavepacket propagation method to describe multimode non-adiabatic dynamics. *International Reviews in Physical Chemistry* **2008**, *27*, 569–606.

(104) Beck, M. The multiconfiguration time-dependent Hartree (MCTDH) method: a highly efficient algorithm for propagating wavepackets. *Physics Reports* **2000**, *324*, 1–105.

(105) Yaghoubi Jouybari, M.; Liu, Y.; Improta, R.; Santoro, F. Ultrafast Dynamics of the Two Lowest Bright Excited States of Cytosine and 1-Methylcytosine: A Quantum Dynamical Study. *J. Chem. Theory Comput.* **2020**, *16*, 5792–5808.

(106) Avagliano, D.; Lorini, E.; González, L. Sampling effects in quantum mechanical/molecular mechanics trajectory surface hopping non-adiabatic dynamics. *Phil. Trans. R. Soc. A* **2022**, *380*, 20200381.

(107) Pieroni, C.; Becuzzi, F.; Creatini, L.; Granucci, G.; Persico, M. Effect of Initial Conditions Sampling on Surface Hopping Simulations in the Ultrashort and Picosecond Time Range. Azomethane Photodissociation as a Case Study. *J. Chem. Theory Comput.* **2023**, *19*, 2430–2445.

(108) Gelin, M. F.; Borrelli, R. Thermo-Field Dynamics Approach to Photo-induced Electronic Transitions Driven by Incoherent Thermal Radiation. *J. Chem. Theory Comput.* **2023**, *19*, 6402–6413.
TOC Graphic

<figure><img src="image_20.png" /></figure>