Large-scale quantum-dynamics with matrix product states

TL;DR

This paper introduces a novel TD-DMRG method for simulating large-scale quantum dynamics, enabling accurate real- and imaginary-time evolution of complex molecular Hamiltonians with over 20 degrees of freedom.

Contribution

It presents an exact, sweep-based TD-DMRG algorithm that overcomes previous scaling limitations, allowing simulation of larger and more complex quantum systems.

Findings

Able to simulate systems with more than 20 degrees of freedom

Successfully applied to excitonic and vibronic Hamiltonians

Enables study of real- and imaginary-time quantum dynamics

Abstract

Dynamical electronic- and vibrational-structure theories have received a growing interest in the last years due to their ability to simulate spectra recorded with ultrafast experimental techniques. The exact time evolution of a molecular system can, in principle, be obtained from the time-dependent version of full configuration interaction. Such an approach is, however, limited to few-atom systems due to the exponential increase of its cost with the system dimension. In the present work, we overcome this unfavorable scaling by employing the time-dependent density matrix renormalization group (TD-DMRG) which parametrizes the time-dependent wavefunction as a matrix product state. The time-dependent Schroedinger equation is then integrated with a sweep-based algorithm, as in standard time-independent DMRG. Unlike other TD-DMRG approaches, the one presented here leads to a set of coupled equations that can be integrated exactly. The resulting theory enables us to study real- and imaginary-time evolutions of Hamiltonians comprising more than 20 degrees of freedom that are challenging for current state-of-the-art quantum dynamics algorithms. We apply our algorithm to the simulation of quantum dynamics of models of increasing complexity, ranging from simple excitonic Hamiltonians to more complex ab-initio vibronic ones.

Figures

Equations

$$|\Psi ⟩=\sum _{{\sigma }_1,...,{\sigma }_L}\sum _{a_1,...,a_{\text{L-1}}}{M}_{1,a_1}^{{\sigma }_1}{M}_{a_1,a_2}^{{\sigma }_2}...M_{a_{\text{L-1}},1}^{{\sigma }_L}|{\sigma }_1,...,{\sigma }_L⟩.$$

$$\mathcal{W}=\sum _{\mathbf{\sigma }{\mathbf{\sigma }}^{\prime }}\sum _{b_1,\dots ,b_{L-1}}^{b_1^{\text{max}},\dots ,b_{L-1}^{\text{max}}}{W}_{1b_1}^{{\sigma }_1,{\sigma }_1^{\prime }}\dots W_{b_{i-1}{b}_i}^{{\sigma }_i,{\sigma }_i^{\prime }}\dots W_{b_{L-1}1}^{{\sigma }_L,{\sigma }_L^{\prime }}|\mathbf{\sigma }⟩⟨{\mathbf{\sigma }}^{\prime }|,$$

$$\sum _{{\sigma }_i^{\prime }}\sum _{a_{i-1}^{\prime },a_i^{\prime }}\sum _{b_{i-1},b_i}{W}_{b_{i-1},b_i}^{{\sigma }_i,{\sigma }_i^{\prime }}{L}_{a_{i-1},a_{i-1}^{\prime }}^{b_{i-1}}{R}_{a_i,a_i^{\prime }}^{b_l}{M}_{a_{i-1}^{\prime },a_i^{\prime }}^{{\sigma }_i^{\prime }}=E{M}_{a_{i-1},a_i}^{{\sigma }_i}.$$

$$H_{(a_{i-1}{\sigma }_i{a}_i),(a_{i-1}^{\prime }{\sigma }_i^{\prime }{a}_i^{\prime })}=\sum _{b_{i-1}}\sum _{b_i}{L}_{a_{i-1},a_{i-1}^{\prime }}^{b_{i-1}}{W}_{b_{i-1},b_i}^{{\sigma }_i,{\sigma }_i^{\prime }}{R}_{a_i,a_i^{\prime }}^{b_i}.$$

$$\mathcal{H}|\Psi (t)⟩=i{\partial }_t|\Psi (t)⟩,$$

$$F(t)={\Vert \mathcal{H}|\Psi ⟩-i{\partial }_t|\Psi ⟩\Vert }^2,$$

$$⟨\delta \Psi |\mathcal{H}-i{\partial }_t|\Psi ⟩=0,$$

$$i{\partial }_t|{\Psi }_{\text{MPS}}⟩={\mathcal{P}}_{{\Psi }_{\text{MPS}}}\mathcal{H}|{\Psi }_{\text{MPS}}⟩.$$

$$|{\Psi }_{\text{MPS}}⟩=\sum _{{\sigma }_1,...,{\sigma }_L}\sum _{a_1,...,a_{\text{L-1}}}{A}_{1,a_1}^{{\sigma }_1}\dots M_{a_{l-1},a_l}^{{\sigma }_l}\dots B_{a_{\text{L-1}},1}^{{\sigma }_L}|{\sigma }_1,\dots ,{\sigma }_L⟩.$$

$$\sum _{{\sigma }_l,a_{l-1}}{A}_{a_{l-1},a_l}^{{\sigma }_l}{A}_{a_{l-1},a_l^{\prime }}^{{\sigma }_l}$$

$$\sum _{{\sigma }_l,a_l}{B}_{a_{l-1},a_l}^{{\sigma }_l}{B}_{a_{l-1}^{\prime },a_l}^{{\sigma }_l}$$

$$|{\Psi }_T⟩=\sum _{i=1}^L|{\Psi }_T^{(i)}⟩,$$

$$|{\Psi }_T^{(i)}⟩=\sum _{{\sigma }_1,...,{\sigma }_L}\sum _{a_1,...,a_{\text{L-1}}}{A}_{1,a_1}^{{\sigma }_1}\dots$$

$${\mathcal{P}}_{{\Psi }_{\text{MPS}}}=\sum _{i=1}^N$$

$$\equiv \sum _{i=1}^N$$

$$|a_{i-1}^{(l)}⟩$$

$$|a_i^{(r)}⟩$$

$$|{\Psi }_{\text{MPS}}⟩=\sum _{a_{i-1}^{(l)}{\sigma }_i{a}_i^{(r)}}{M}_{a_{i-1}^{(l)},a_i^{(r)}}^{{\sigma }_i}|a_{i-1}^{(l)}{\sigma }_i{a}_i^{(r)}⟩.$$

$$|a_i^{(l)}⟩=\sum _{a_{i-1}^{(l)},{\sigma }_i}{A}_{a_{i-1}^{(l)},a_i^{(l)}}^{{\sigma }_i}|a_{i-1}^{(l)}{\sigma }_i⟩.$$

$$e^{-i{\mathcal{P}}_{{\Psi }_{\text{MPS}}}\mathcal{H}t}\approx \prod _{i=1}^N{e}^{-i{\mathcal{P}}_i^{(1)}\mathcal{H}t}\prod _{i=1}^{N-1}{e}^{i{\mathcal{P}}_i^{(2)}\mathcal{H}t}.$$

$$\frac{\text{d}}{\text{d}t}{M}_{a_{i-1}{a}_i}^{{\sigma }_i}=-\sum _{a_{i-1}^{\prime },a_i^{\prime },{\sigma }_i^{\prime }}\mathrm{i}{H}_{a_{i-1}{\sigma }_i{a}_i,a_{i-1}^{\prime }{\sigma }_i^{\prime }{a}_i^{\prime }}{M}_{a_{i-1}^{\prime }{a}_i^{\prime }}^{{\sigma }_i^{\prime }},$$

$$|a_i^{(l)}⟩=\sum _{a_{i-1}^{(l)},{\sigma }_i}{A}_{a_{i-1}^{(l)}{\sigma }_i,a_i^{(l)}}|a_{i-1}^{(l)}{\sigma }_i⟩,$$

$$M_{a_{i-1}^{(l)},a_i^{(r)}}^{{\sigma }_i}=\sum _{a_i^{(l{)}^{\prime }}}{A}_{a_{i-1}^{(l)},a_i^{(l{)}^{\prime }}}^{{\sigma }_i}{Q}_{a_i^{(l{)}^{\prime }},a_i^{(r)}}.$$

$$H_{\left(a_i^{(l)}{a}_i^{(r)},a_i^{(l{)}^{\prime }}{a}_i^{(r{)}^{\prime }}\right)}$$

$$=\sum _{{\tilde{a}}_{i-1},{\tilde{a}}_{i-1}^{\prime }}\sum _{{\sigma }_i,{\sigma }_i^{\prime }}{A}_{{\tilde{a}}_{i-1},a_i}^{{\sigma }_i,†}{A}_{{\tilde{a}}_{i-1}^{\prime },a_i^{\prime }}^{{\sigma }_i^{\prime }}⟨{\tilde{a}}_{i-1}^{(l)}{\sigma }_i{a}_i^{(r)}\mid \mathcal{H}\mid {\tilde{a}}_{i-1}^{(l){,}^{\prime }}{\sigma }_i^{\prime }{a}_i^{(r){,}^{\prime }}⟩.$$

$$L_{a_i^{(l)},a_i^{(l{)}^{\prime }}}^{b_l}=\sum _{{\sigma }_i,{\sigma }_i^{\prime }}\sum _{a_{i-1},a_{i-1}^{\prime }}\sum _{b_{i-1}}{A}_{{\tilde{a}}_{i-1}^{(l)},a_i^{(l)}}^{{\sigma }_i,†}{A}_{{\tilde{a}}_{i-1}^{(l{)}^{\prime }},a_i^{(l{)}^{\prime }}}^{{\sigma }_i^{\prime }}{W}_{b_{i-1},b_i}^{{\sigma }_i,{\sigma }_i^{\prime }}{L}_{a_{i-1}^{(l)},a_{i-1}^{(l{)}^{\prime }}}^{b_{i-1}},$$

$$H_{\left(a_i^{(l)}{a}_i^{(r)},a_i^{(l{)}^{\prime }}{a}_i^{(r{)}^{\prime }}\right)}=\sum _{b_l}{L}_{a_i^{(l)},a_i^{(l{)}^{\prime }}}^{b_l}{R}_{a_i^{(r)},a_i^{(r{)}^{\prime }}}^{b_l}.$$

$$|{\Psi }_{\text{MPS}}⟩$$

$$\times Q_{a_i^{(l)},a_i^{(r)}}{B}_{a_i^{(r)},a_{i+1}^{(r)}}^{{\sigma }_{i+1}}\dots B_{a_{L-1}^{(r)},1}^{{\sigma }_L}|{\sigma }_1\dots {\sigma }_L⟩$$

$$\equiv \sum _{a_i^{(l)},a_i^{(r)}}{Q}_{a_i^{(l)},a_i^{(r)}}|a_i^{(l)}{a}_i^{(r)}⟩.$$

$$\Delta M_{a_{i-1}{\sigma }_i,b_i{a}_i^{\prime }}=\sum _{b_{i-1}}\sum _{a_{i-1}^{\prime }}\sum _{{\sigma }_i^{\prime }}{L}_{a_{i-1}{a}_{i-1}^{\prime }}^{b_{i-1}}{M}_{a_{i-1}^{\prime }{a}_i^{\prime }}^{{\sigma }_i^{\prime }}{W}_{b_{i-1}{b}_i}^{{\sigma }_i{\sigma }_i^{\prime }},$$

$$\mathcal{W}|\Psi ⟩$$

$$=\sum _{a_1^{\prime }\dots a_{L-1}^{\prime }}\sum _{{\sigma }_1^{\prime }\dots {\sigma }_L^{\prime }}{N}_{1a_1^{\prime }}\cdots N_{a_{L-1}^{\prime }1}|{\sigma }_1^{\prime },\dots ,{\sigma }_L^{\prime }⟩,$$

$$N_{a_{i-1}^{\prime }{a}_i^{\prime }}^{{\sigma }_i}=\sum _{{\sigma }_i^{\prime }}{W}_{b_{i-1}{b}_i}^{{\sigma }_i{\sigma }_i^{\prime }}{M}_{a_{i-1}{a}_i}^{{\sigma }_i^{\prime }}.$$

Full text

\abbreviations

1 Introduction

The vast majority of the algorithms in quantum chemistry aims at solving the time-independent Schrödinger equation. A wide class of phenomena are, however, difficult to target within a time-independent framework. This is the case, for example, for the simulation of spectra obtained with time-resolved spectroscopic techniques that probe a molecular system under non-equilibrium conditions with pulsed light.1, 2, 3, 4, 5 Such spectra are readily obtained from a direct solution of the time-dependent Schrödinger equation with quantum dynamics approaches.6, 7, 8, 9 Similarly, high-order response properties or high-lying excited states can, in principle, be obtained with time-independent approaches, but their calculation is easier within a time-dependent framework.10, 11, 12

Multiconfigurational time-dependent Hartree (MCTDH)13 is currently the reference method for quantum-dynamics simulations with vibrational14, 15 and electronic Hamiltonians.16, 17 As any full-configuration interaction (full-CI) approach, MCTDH suffers from the curse of dimensionality because its computational cost increases exponentially with system size. This unfavorable scaling can be partially overcome with multi-scale formulations of MCTDH, such as the Gaussian MCTDH (G-MCTDH)18, 19 or the multi-layer MCTDH (ML-MCTDH).20, 21, 22, 23, 24, 25 In G-MCTDH, only a subset of the nuclear degrees of freedom are treated with conventional MCTDH, while the rest is described semiclassically with Gaussian wavepackets. Conversely, in ML-MCTDH, strongly correlated vibrational degrees of freedom are first contracted together with MCTDH. The resulting pre-contracted basis is truncated and then employed in a second MCTDH calculation. Systems with several hundreds of vibrational degrees of freedom can be targeted by this hierarchical treatment.24, 26, 27, 28 The efficiency of such multi-layer schemes depends on the choice and definition of the layers, i.e., on the choice of the modes treated with conventional MCTDH for G-MCTDH, and on the precontraction scheme for ML-MCTDH. Well-established schemes to perform this selection are, however, not yet available.

The density matrix renormalization group (DMRG) algorithm has become one of the reference methods for variational calculations on large systems in recent years. Unlike full-CI, the scaling of DMRG is polynomial in the system size and therefore does not suffer from the curse of dimensionality.29, 30, 31, 32 In its original formulation, DMRG is designed for targeting ground states of Hamiltonians of linear systems with nearest-neighbor interactions. It has been then successfully applied to electronic33, 34, 35, 36, 37 and vibrational38, 39, 40, 41 Hamiltonians. DMRG replaces the exact diagonalization of the Hamiltonian of an $N$-body system with a series of smaller eigenvalue problems. Within each problem, only one site (i.e., an orbital for the electronic Hamiltonian) is treated explicitly and the effects of the remaining sites are represented through an effective basis, known as renormalized basis that is iteratively constructed to provide the best representation of the overall wavefunction in a least-squares sense. The dimension of this basis, also known as bond dimension, determines the accuracy of DMRG.

The generalization of DMRG to time-dependent processes faces two problems. The first one concerns the update of the renormalized basis during the propagation. The renormalized basis is known only at the initial time $t$ and must be optimized iteratively to optimize the representation of the wave function over the whole time interval. In the first TD-DMRG formulation, this basis was constructed at the beginning of the simulation and kept constant during the propagation.42 Not surprisingly, such a choice leads to a continuous decrease of accuracy with increasing time, which makes it useful for short-time propagations only.43 More recent TD-DMRG formulations, known as adaptive TD-DMRG,44, 45, 46, 47 build the renormalized basis to provide a balanced least-squares representation of the wavefunction for the whole time-step of the propagation.

For Hamiltonians containing only nearest-neighbor interactions, an efficient solution is to express the time evolution operator through a Suzuki-Trotter factorization48 as a product of terms involving only one and two sites at a time. The time evolution operator can be applied in a sweep-based fashion as in time-independent DMRG. The resulting theory is known as time-evolving block decimation (TEBD).49, 50 Despite its successful applications for model Hamiltonians, TEBD has strong limitations when applied to more complex Hamiltonians, such as the ones of interest in quantum chemistry that contain long-range Coulomb interactions (highlighted by 4-indexed integrals in its second quantization form).

Shortly after its introduction, it was realized that DMRG optimizes wavefunctions encoded as matrix product states (MPSs),51 a parametrization known as tensor train in the mathematical community.52 The MPS form can be generalized to operators to be written in matrix product operator (MPO) format.53, 54, 55, 56 The availability of an analytic parametrization for wavefunctions and operators enables one to simulate the time evolution with the so-called Dirac-Frenkel variational principle.57, 58 This leads to an approximate form of the time-dependent Schrödinger equation that is projected onto the manifold formed by all MPSs with a given bond dimension.59 In practice, this is achieved by scaling the Hamiltonian, expressed as an MPO, by the so-called tangent-space projector.60 For this reason, the resulting theory is known as tangent-space formulation of TD-DMRG.60, 61, 62 We note that other algorithms to propagate wavefunctions encoded as MPSs have been proposed. Greene and Batista introduced the tensor-train spit-operator Fourier Transform (TT-SOFT)63 theory, in which the MPS and the propagator are constructed from their representation in the position/momentum basis. More recently, Frahm and Pfannkuche extended the Runge-Kutta and Lanczos propagation schemes to wave functions encoded as MPS.64 Both approaches share, however, a common limitation because the compact structure of an MPS is not fully exploited and, therefore, the bond dimension increases during the propagation. To avoid an uncontrolled increase of its size, the MPS needs to be truncated after each time-step. This does not hold true for the tangent-space-based scheme that implicitly keep the dimension of the MPS fixed.

The efficiency of TD-DMRG depends on the possibility of representing the wavefunction as a compact MPS for the whole propagation. The area law65 ensures that such representation is possible for ground states of short-range, gapped Hamiltonians. A time-dependent analog of this theorem is, however, not known, and it has been demonstrated for some spin Hamiltonians66, 67 that the entanglement entropy increases during a time-dependent propagation even with a dynamical update of the renormalized basis. The dimension of the MPS must then also increase during the propagation to represent the wavefunction with constant accuracy. For this reason, particular care is needed when assessing the convergence of TD-DMRG propagations with respect to the renormalized basis set size.

The tangent-space-based TD-DMRG algorithm has been mostly applied to spin Hamiltonians68, 69 and chemical applications have been limited to excitonic models.70, 71, 72 In the present work, we investigate this algorithm on realistic vibrational and vibronic Hamiltonians containing both short- and long-range interactions. Moreover, we improve the flexiblity of the algorithm by allowing for a dynamical selection of the dimension of the MPS, within either a single-site73 and a two-site62 scheme. First, we study the performance of our implementation on an excitonic Hamiltonian74 and compare to alternative TD-DMRG approaches.75, 70, 71 Then, we simulate the full-dimensional excited-state nuclear dynamics of pyrazine by treating all 24 modes at the quantum level. Due to its size, this system represents a challenge for full quantum dynamics approaches and usually requires multi-layer treatments. Finally, the flexibility of the approach will be demonstrated by extending it to the simulation of the imaginary-time propagation of MPSs. This results in an efficient algorithm to optimize ground and excited states encoded as MPSs.

2 Theoretical framework

2.1 MPS/MPO formulation of DMRG

Before presenting the TD-DMRG theory, we describe briefly the time-independent formulation of DMRG to introduce the relevant notation. A more detailed discussion can be found in several recent reviews.31, 32, 33, 37

DMRG is an iterative optimization algorithm to approximate the ground state wavefunction $|\>\Psi_{0}\>\rangle$ of an $L$-sites quantum system. In its MPS/MPO formulation, the wavefunction is expressed as an MPS as51, 53

[TABLE]

where $|\sigma_{1},...,\sigma_{L}\rangle=|\bm{\sigma}\rangle$ are many-particle basis set elements obtained from the direct product of the basis of each site and expressed as occupation number vectors (ONV). The $L$ tensors $\bm{M}^{\sigma_{i}}=\{M_{a_{\text{i-1}},a_{\text{i}}}^{\sigma_{i}}\}$ (usually known as site tensors) have one index ($\sigma_{i}$) associated to the local basis of the $i$-th site, and two additional indices with maximum dimension $m$ (note that $\mathbf{M}^{\sigma_{1}}$ and $\mathbf{M}^{\sigma_{L}}$ are row and column vectors, respectively, of maximum size $m$). A full-CI wavefunction can, in principle, be expressed exactly as an MPS, but its dimension would increase exponentially with system size. The key advantage of DMRG is the fact that ground states of Hamiltonians containing short-range interactions can be encoded by an MPS whose dimension is largely independent of the system size.65 For more general Hamiltonians, with both short- and long-range interactions, MPSs with larger $m$ are usually required to obtain a converged representation of the wavefunction.

Following the site decomposition of an MPS, the matrix product operator (MPO) representation of an operator $\mathcal{W}$ is obtained as53, 32, 76, 55

[TABLE]

where $\bm{W}^{\sigma_{l},\sigma_{l}^{\prime}}=\{W_{b_{l-1},b_{l}}^{\sigma_{l},\sigma_{l}^{\prime}}\}$ are 4-dimensional tensors with two indices ($\sigma_{l}$ and $\sigma_{l}^{\prime}$) associated to the local basis of the $i$-th site.

The energy functional obtained combining Eqs. (1) and (2) is a complex, non-linear function of the entries $\{M_{a_{\text{i-1}},a_{\text{i}}}^{\sigma_{i}}\}$. If the energy is minimized with respect to the $l$-th site tensor by keeping all the other ones fixed, we obtain32

[TABLE]

where $\bm{L}$ and $\bm{R}$ are three-dimensional tensors usually defined as left and right boundary, respectively. The left boundary is obtained by contracting the tensor network representing $\langle\>\Psi\mid\mathcal{H}\mid\Psi\>\rangle$ for all sites ranging from 1 to $i$-1. Similarly, the right boundary is the partial contraction of the MPS with the MPO for the sites ranging from $i$+1 to $L$. Eq. (3) can be recast as a standard eigenvalue equation $\bm{H}\bm{v}=E\bm{v}$, where $\bm{v}$ is a vector collecting all the entries of the $M_{a_{l-1},a_{l}}^{\sigma_{l}}$ tensor. The matrix $\bm{H}$ is defined as

[TABLE]

In DMRG iterations, Eq. (3) is solved sequentially for each site, starting from the first site up to the end of the chain (forward sweep) and then repeating this procedure with the reverse direction (backward sweep). By consistently following the lowest energy eigenpair, an approximation of the ground state is iteratively constructed.

2.2 MPS/MPO formulation of TD-DMRG

To apply the DMRG optimization algorithm of MPS in quantum dynamics simulations, the time-dependent Schrödinger equation,

[TABLE]

must be solved under the condition that the wavefunction $|\Psi(t)\rangle$ is an MPS at any time $t$. It follows from Eq. (5) that the variation of the wavefunction $|\Psi(t)\rangle$ over a time-step $\Delta t$ is proportional to $\mathcal{H}|\Psi\rangle\times\Delta t$. It is known that, if $|\Psi(t)\rangle$ is represented as MPS with bond dimension $m$, an exact representation of $\mathcal{H}|\Psi\rangle$ will require an MPS with bond dimension $b_{i}m$, where bi is the bond dimension of the MPO (see Eq. (2)) and $m$ is the bond dimension of the MPS (see Eq. (1)). A naive integration of Eq. (5) would, therefore, lead to an unbounded increase of the MPS bond dimension.

For approximate wavefunction expressed by some parametrization scheme, Eq. (5) cannot be solved exactly in the most general case. The best approximation to the exact wavefunction $|\Psi(t)\rangle$ is obtained from the minimization of the Dirac-Frenkel functional F(t) defined as57, 58

[TABLE]

where the minimization is performed within the set $S$ of all wavefunctions that can be expressed with the parametrization encoding $|\Psi(t)\rangle$. Eq. (6) is equivalent to

[TABLE]

where $|\delta\Psi\rangle$ is a generic infinitesimal variation of the wavefunction $|\Psi\rangle$ in $S$. The best approximation of the exact, time-dependent wavefunction $|\Psi(t)\rangle$ is therefore orthogonal, at each time step, to the error with respect to the exact solution of the time-dependent Schrödinger equation $(\mathcal{H}-i\partial_{t})|\Psi\rangle$.

Let now $\mathcal{M}_{\text{MPS}}(m)$ be the set of all MPSs with bond dimension $m$. It was proven60 that $\mathcal{M}_{\text{MPS}}(m)$ is a smooth submanifold of the full Hilbert space. This means, in practice, that it is a subset of the full Hilbert space not displaying discontinuities and this facilitates the approximation of $\mathcal{M}_{\text{MPS}}(m)$ in the vicinity of this MPS as a linear subspace of the full configurational space. This linear approximation is known as the tangent space to $\mathcal{M}_{\text{MPS}}(m)$ calculated for the reference MPS $|\>\Psi_{\text{MPS}}\>\rangle$. Eq. (7) can be conveniently expressed in terms of the projector $\mathcal{P}_{\Psi_{\text{MPS}}}$ onto this tangent space,

[TABLE]

$\mathcal{P}_{\Psi_{\text{MPS}}}$ has a non-trivial effect on the right-hand side of Eq. (8) because, as already highlighted, the bond dimension of $\mathcal{H}|\>\Psi_{\text{MPS}}\>\rangle$ is larger than the one of $|\>\Psi_{\text{MPS}}\>\rangle$, $m$. The projection operation ensures that the dimension of the MPS is kept constant during the propagation.

A closed-form expression for $\mathcal{P}_{\Psi_{\text{MPS}}}$ was derived recently77, 61 and applied to the solution of Eq. (8).62 We recall that an MPS is defined as canonically normalized with respect to site $l$ if it is represented as

[TABLE]

where $\bm{A}$ and $\bm{B}$ fulfill the following orthogonality relation,

[TABLE]

An element $|\Psi_{T}\rangle$ of the tangent space to an MPS $|\Psi_{\text{MPS}}\rangle$ of a given rank $m$ can be decomposed as,61

[TABLE]

where,

[TABLE]

$\bm{A}$ and $\bm{B}$ are given in Eq. (9) and ${N}_{a_{i-1},a_{i}}^{\sigma_{i}}$ is a normalized tensor orthogonal to ${M}_{a_{i-1},a_{i}}^{\sigma_{i}}$ introduced in Eq. (9). The orthogonality constraint ensures that the parametrization of Eq. (13) is non-redundant. The projector $\mathcal{P}_{\Psi_{\text{MPS}}}$ may be written as,

[TABLE]

with the left ($|a_{i-1}^{(l)}\rangle$) and right ($|a_{i}^{(r)}\rangle$) renormalized bases on site $i$ defined as

[TABLE]

Eq. (14) can be rationalized as follows. The first term is the sum of projectors onto the $|a_{i-1}^{(l)}\rangle\otimes|\sigma_{i}\rangle\otimes|a_{i}^{(r)}\rangle$ basis for each site $i$. The expression of an MPS normalized with respect to site $i$ in this product basis is

[TABLE]

The MPS can be expanded as in Eq. (16) for any site $i$ and this might suggest that the union of all product bases $|a_{i-1}^{(l)}\sigma_{i}a_{i}^{(r)}\rangle$ constitutes a complete basis for a rank-$m$ MPS. However, this basis set would be overcomplete. To show that, consider a left sweep from the $i$-th to the ($i$+1)-th site. The left renormalized basis for the ($i$+1)-th site ($|\>a_{i}^{(l)}\>\rangle$) can be expressed as a function of $|\>a_{i-1}^{(l)}\>\rangle$

[TABLE]

$A_{a_{i-1}^{(l)},a_{i}^{(l)}}^{\sigma_{i}}$ is obtained by singular value decomposition (SVD) or QR factorization of the optimized tensor $M_{a_{i-1}^{(l)},a_{i}^{(l)}}^{\sigma_{i}}$. Eq. (17) highlights that $|a_{i}^{(l)}\rangle$ and $|a_{i-1}^{(l)}\sigma_{i}\rangle$ are linearly dependent basis sets. Therefore, $|a_{i}^{(l)}\sigma_{i}a_{i}^{(r)}\rangle$ are contained in both the spaces spanned by $|a_{i-1}^{(l)}\sigma_{i}a_{i+1}^{(r)}\rangle$ and by $|a_{i}^{(l)}\sigma_{i}a_{i+1}^{(r)}\rangle$. Including only the first sum in Eq. (14) would, therefore, lead to a double counting of the $|a_{i-1}^{(l)}\sigma_{i}a_{i}^{(r)}\rangle$ basis functions. This effect is corrected by the second sum appearing in Eq. (14).

The formal solution of Eq. (8) is $e^{-i\mathcal{P}_{\Psi_{\text{MPS}}}\mathcal{H}t}|\Psi_{\text{MPS}}(0)\rangle$, but a closed form of the exponential operator $e^{-i\mathcal{P}_{\Psi_{\text{MPS}}}\mathcal{H}t}$ cannot be derived. As suggested in Ref. 62, the exponential operator can be approximated by a first-order Lie-Trotter splitting as

[TABLE]

where higher-order terms involving the commutator between $\mathcal{P}_{i}^{(1)}$ and $\mathcal{P}_{i}^{(2)}$ have been neglected. Although a Lie-Trotter decomposition is also employed in TEBD,49 the factorization introduced in Eq. (18) is much more general. On the one hand, TEBD factorizes the time-evolution operator as a product of exponentials, one for each term appearing in the Hamiltonian. Such factorization is, however, efficient only for nearest-neighbor Hamiltonians. On the other hand, the Lie-Trotter splitting in Eq. (18) is based on the expression of the projector $\mathcal{P}_{\Psi_{\text{MPS}}}$ as a sum of site terms, but the factorization of the Hamiltonian operator is avoided.

The order of application of the local operators appearing in Eq. (18) follows naturally from the structure of the MPS. The $|a_{i}^{(l)}\rangle$ basis is easily constructed from its equivalent of the previous site $|a_{i-1}^{(l)}\rangle$, and it is, therefore, convenient to apply the operators $\mathcal{P}_{i}^{(1)}$ and $\mathcal{P}_{i}^{(2)}$ one after the other for increasing values of $i$. The most appealing feature of this MPS/MPO-based formulation of TD-DMRG is that the resulting coupled differential equation can be integrated analytically. The representation of the Hamiltonian in the $|a_{i-1}^{(l)}\sigma_{i}a_{i}^{(l)}\rangle$ basis is given in Eq. (4). The differential equation obtained for the exponential operators containing $\mathcal{P}_{i}^{(1)}$ can then be expressed as

[TABLE]

which is a first-order, constant-coefficient linear differential equation that can be integrated exactly. The integration of the differential equation involving the $\mathcal{P}_{i}^{(2)}$ projector follows a similar strategy, once the representation of the Hamiltonian in the $|a_{i}^{(l)}a_{i}^{(r)}\rangle$ basis is available. We recall that the $|a_{i}^{(l)}\rangle$ basis is obtained from the $|a_{i-1}^{(l)}\sigma_{i}\rangle$ by QR factorization,

[TABLE]

where the $\bm{A}$ tensor, which is left normalized according to Eq. (11), is defined as

[TABLE]

The matrix representation of the Hamiltonian in the $|a_{i}^{(l)}a_{i}^{(r)}\rangle$ basis is then

[TABLE]

By combining Eqs. (4) and (22) and considering that the left boundary for site $i$ can be obtained recursively from the one of site $i$+1 through the following recursive relation,32, 76

[TABLE]

it is easy to show that Eq. (22) can be expressed in a more compact form as

[TABLE]

The representation of an MPS in the $|a_{i}^{(l)}a_{i}^{(r)}\rangle$ basis is obtained by combining Eqs. (9) and (17),

[TABLE]

Hence, the representation is given by the $\bm{Q}$ matrix, which is often referred to as zero-site tensor.62

2.3 Single- and two-site integrators

The theoretical framework introduced above assumes that the bond dimension of the initial wavefunction $|\>\Psi(t=0)\>\rangle$ is kept constant throughout the whole propagation time. This is a major limitation in practical applications of TD-DMRG because the bond dimension $m$ of the wavefunction at the beginning of the propagation is determined by the initial conditions. In Section 3, we will study the nuclear dynamics of a molecule after photoexcitation to an electronically excited state. Potential energy surfaces are modeled with the harmonic approximation and, therefore, the initial state is the vibrational ground state of the lowest electronic state. In the occupation number vector representation, this is a basis vector ($|\>00\cdots 0\>\rangle$) that corresponds to an MPS with $m$=1. A straightforward application of the propagation algorithm described above would, therefore, constrain the bond dimension of the MPS to one. Here, we propose two alternative strategies, adapted from the standard time-independent variant of DMRG, to alleviate this limitation allowing for a dynamical adaptation of $m$.

The first route is based on the subspace expansion algorithm originally devised for tensor trains78, 79 and later applied to single-site time-independent DMRG.73 We will discuss the algorithm for a forward sweep because the extension to a backward sweep is trivial. After its optimization, an MPS tensor $M_{a_{i-1},a_{i}}^{\sigma_{i}}$ is reshaped as a matrix $M_{(a_{i-1}\sigma_{i},a_{i})}$ and $m^{\prime}$ columns are added to the resulting matrix by increasing its bond dimension to $m+m^{\prime}$. The MPS tensor on the next site $M_{a_{i},a_{i+1}}^{\sigma_{i+1}}$ is reshaped as $M_{(a_{i},a_{i+1}\sigma_{i+1})}$ and $m^{\prime}$ null rows are added to match the new dimensions of the preceding tensor. This modification does not change the MPS because the rows added on the $i$-th site tensor are contracted with null columns on the ($i$+1)-th site tensor. Nevertheless, this step enlarges the basis $a_{i}^{(l)}$, and therefore, the optimization of the MPS for the ($i$+1)-th site will be performed in a larger vector space. The efficiency of the subspace expansion step depends strongly on the $m^{\prime}$ vectors that are added to $M_{a_{i-1},a_{i}}^{\sigma_{i}}$ at the end of a microiteration step. In principle,80 these vectors can be chosen randomly without altering the efficiency of the DMRG optimization. A more reliable alternative is to construct the vectors as73

[TABLE]

where the resulting matrix $\Delta\bm{M}$ is to be added to the $\bm{M}^{\sigma_{i}}$ tensor. Eq. (26) can be intuitively motivated as follows. We have already mentioned that applying an operator $\mathcal{W}$ encoded as an MPO to an MPS increases the bond dimension of the MPS. In fact,

[TABLE]

where the index $a_{i}^{\prime}$ combines the indices $b_{i}$ and $a_{i}$, and

[TABLE]

Therefore, for a given site $i$, the left- and right-renormalized bases representing $\mathcal{H}|\>\Psi\>\rangle$ are larger than the ones of $|\>\Psi\>\rangle$. The subspace expansion algorithm enlarges the left-renormalized basis for $|\>\Psi\>\rangle$ with the left-renormalized basis of $\mathcal{H}|\>\Psi\>\rangle$. It is easy to demonstrate that the matrix defined in Eq. (26) expresses this new additional basis in terms of the $|\>a_{i-1}\>\rangle\otimes|\>\sigma_{i}\>\rangle$ product basis.

A direct application of the subspace expansion algorithm would increase the bond dimension by a factor of $b_{i}$ per microiteration, where $b_{i}$ is the bond dimension of the MPO introduced in Eq. (27). To avoid an uncontrolled increase of $m$, it was suggested in Ref. 73 to scale the matrix $\Delta\bm{M}^{\sigma_{i}}$ by a parameter $\mu<1$ and then truncate $\bm{M}^{\sigma_{i}}$ with SVD. Even if adaptive schemes to select $\mu$ have been proposed,73 in the present work $\mu$ is kept constant for $n_{\mu}$ sweeps and then set to 0 because already such simple choice provides reliable and converged results.

A second, alternative algorithm to dynamically control the bond dimension $m$ is the two-site optimizer,32, 76 in which the energy is minimized with respect to the simultaneous variation of two consecutive sites. In practice, a two-site tensor is defined as,

[TABLE]

and the optimization is carried out for $\bm{T}^{\sigma_{i}\sigma_{i+1}}$ instead of $M_{a_{i-1}a_{i}}^{\sigma_{i}}$. To return to the MPS factorization, $\bm{T}$ is first reshaped as a matrix as $T_{(a_{i-1}\sigma_{i},a_{i+1}\sigma_{i+1})}$ and then factorized by SVD,

[TABLE]

where $\bm{U}$ and $\bm{V}$ are orthogonal matrices and $\bm{S}$ is a diagonal matrix containing the singular values. The upper bound $m^{\prime}$ to the sum over the index $a_{i}^{\prime}$ is equal to the rank of $T_{(a_{i-1}\sigma_{i-1},a_{i+1}\sigma_{i})}$ and can, in general, be different from the rank of the matrix before the optimization. To avoid an uncontrolled increase of the updated bond dimension $m^{\prime}$, the SVD can be truncated. For $S_{\text{keep}}$ being the set of diagonal elements that are kept after the truncation, it is possible to demonstrate81 that the error $\eta_{\text{err}}$ introduced in the approximation is given by

[TABLE]

It is, therefore, possible to fix a threshold on the error $\eta_{\text{thresh}}$ and keep the minimal number of vectors in the SVD truncation so that the error is smaller than $\eta_{\text{thresh}}$. This idea is at the heart of the dynamical block state selection scheme81 and can be extended also to time-dependent propagations.

2.4 MPO representation of vibrational and vibronic Hamiltonians

The efficiency of DMRG relies on the ability of representing the Hamiltonian as a compact MPO. Different algorithms have been proposed to build MPO representations of general operators.82, 54, 83, 76, 56 Among them, we chose the one introduced in Ref. 76 for electronic structure problems which we recently generalized to vibrational Hamiltonians.39, 40 This algorithm requires as input parameter a representation of the Hamiltonian in terms of local operators. The MPO is then represented as a graph, where each local operator is associated with a different node. A single term in the Hamiltonian is represented by a connection between the vertices. The major advantage of this construction scheme over other ones available in the literature is its generality. In fact, it supports strings of second-quantized creation and annihilation operators with arbitrary length and is not limited to Hamiltonians with short-range interactions.

Even if TD-DMRG is not limited to any specific form of the Hamiltonian, for the present work we chose to focus on vibrational and vibronic Hamiltonians. We will turn to electronic Hamiltonians in future works. Regarding the former class, we already discussed in Ref. 39 how to generate an MPO representation of the Watson Hamiltonian,

[TABLE]

where Cartesian normal modes $\bm{q}$ have been used as a reference coordinate set and the potential operator $\mathcal{V}(\bm{q})$ has been expanded as a Taylor series in powers of $\bm{q}$. A second-quantization representation if Eq. (31) is obtained by applying the canonical quantization rules to the position operator $\hat{q}_{i}$ and its conjugate momentum $\hat{p}_{i}$,

[TABLE]

where $\hat{b}_{i}^{+}$ and $\hat{b}_{i}$ are bosonic creation and annihilation operators, respectively. Within this canonical quantization representation, the basis set is represented in terms of occupation number vectors (ONVs) $|\>n_{1},n_{2}\ldots n_{L}\>\rangle$ where $n_{i}$ labels the occupation number of mode $i$.

Vibronic Hamiltonians describe the coupled dynamics of nuclei and electrons and, therefore, their encoding as MPOs is less straightforward. A wavefunction including nuclear and electronic degrees of freedom can be parametrized as

[TABLE]

where we label the electronic wavefunction obtained within the Born-Oppenheimer approximation at a given reference geometry as $|\>\alpha\>\rangle$. $|\>\Psi_{\text{vib}}^{\alpha}\>\rangle$ is the vibrational wavefunction associated to the nuclear motion on the $\alpha$-th potential energy surface. For each value of $\alpha$, $|\>\Psi_{\text{vib}}^{\alpha}\>\rangle$ can be expanded as linear combination of harmonic oscillator eigenfunctions for each degree of freedom $i$. In the ONV representation, Eq. (33) can be rewritten as

[TABLE]

In the first line of Eq. (34), the electronic degrees of freedom are labeled through an ONV with a single index, $|\>\alpha\>\rangle$, as in Eq. (33), where different electronic states are associated to different values for $\alpha$. In the second line, a separate list of indices $|\>\alpha_{1},\ldots,\alpha_{N_{\text{ex}}}\>\rangle$ is employed for each electronic state in the ONV representation. The indices $\alpha_{i}$ can be either 0, if the molecule is not in electronic state $i$, or 1, if electrons are solely described by the $i$-th electronic wavefunction. Therefore, the following constraint must be fulfilled

[TABLE]

Eq. (35) reflects the conservation of the number of excited electronic states (usually known as excitons) and induces a block structure of the MPS. This structure can be exploited to reduce the computational cost of DMRG.84

The MPS form of the expansion given in the second line of Eq. (34) is obtained as,

[TABLE]

In principle, it is possible to derive an MPS representation of the first line of Eq. (34). However, the resulting parametrization would be less flexible. As already discussed in the literature, the efficiency of DMRG depends on the sorting of the sites on the lattice and different strategies to derive an optimal sorting have been proposed85, 86, 87, 88 of which the ordering according to orbital entropies has been demonstrated to be very efficient. In Eq. (36), each electronic state is associated to a site and it is, therefore, possible to order the sites so that the electronic states are not close to each other. This is not possible in the first representation given in the first line of Eq. (34), in which only one site is associated to the electronic degrees of freedom of the molecule.

This second, more flexible occupation number vector representation of the electronic degrees of freedom can be associated to second-quantizated creation and annihilation operators $\hat{a}_{\alpha}^{+}$ and $\hat{a}_{\alpha}$, respectively, such that

[TABLE]

The operators defined in Eq. (37) satisfy fermionic anticommutation rules but, due to the constraint of Eq. (35), only one electronic site is occupied. Hence, no antisymmetry constraint must be imposed. For cases in which multiple excitations are present, the fermionic symmetry can be included as discussed, for example, in Ref. 76.

3 Results

We apply the TD-DMRG theory developed here to three examples. The first one is an aggregate of perylene bisimide-based dye molecule, whose excited-state dynamics can be modeled with an excitonic Hamiltonian. This molecular system has been already studied with both MCTDH89, 90, 74 and DMRG,75 and therefore, it is a suitable starting point for analyzing our TD-DMRG implementation against alternative formulations.

The second molecule studied in this work is pyrazine. The excited state dynamics for the S1 $\leftarrow$ S2 non-radiative decay simulated on the potential energy surface of Ref. 91 is the most common test case for new quantum dynamical approaches. Among the various TD-DMRG formulations proposed in the literature, only the TT-SOFT algorithm63 presented in the introduction has been applied, up to now, to pyrazine. A full-dimensional quantum dynamics study of pyrazine, performed by treating all 24 modes at the same level of accuracy, is challenging for most state-of-the-art approaches. We show here how TD-DMRG allows us to obtain converged spectra for this system by treating all vibrational degrees of freedom at the same footing without the need of introducing any hierarchical treatment of the modes, as is done in ML-MCTDH. The pyrazine example will also demonstrate the reliability of DMRG in targeting general vibronic Hamiltonians lacking a regular interaction pattern.

Finally, we extend TD-DMRG to imaginary time propagations to optimize ground and excited vibrational levels of ethylene. Imaginary-time propagation is the idea underlying a diffusion Monte Carlo approach and has been employed to optimize ground and excited states of electronic92 and vibrational93, 94, 95 Hamiltonians.

3.1 Computational details

We implemented the tangent-space-based TD-DMRG theory in our QCMaquis-V program designed to perform DMRG calculations on vibrational Hamiltonians39, 40 and derived from the QCMaquis program for electronic structure calculations.34, 76

If not otherwise specified, TD-DMRG propagations were performed with a time step of 1 fs and an overall propagation time of 4000 fs. Single-site simulations employ the subspace expansion variant presented in Ref.73. The scaling factor $\mu$ for the MPS $\Delta\bm{M}$ (recall Eq.(26 for the definition of $\Delta\bm{M}$) was set to 10*-8* for the first 100 sweeps (note that each sweep corresponds to one time step), and to 10*-15* up to the 200-th sweep, and then set to 0. Other combinations of values for $\mu$ and numbers of sweeps have been tested and provide results equivalent to the ones reported in the next section.

For excited-state nuclear dynamics, vibronic spectra were obtained from the Fourier transformation of the autocorrelation function $C(t)=\langle\>\Psi(t)\mid\Psi(0)\>\rangle$. Broadening effects were reproduced by scaling the autocorrelation function by an exponential damping function $f(t)=e^{-t/\tau}$, corresponding to a Lorentzian broadening function in the frequency domain. The parameter $\tau$ determines the decay rate of the autocorrelation function, and therefore, the broadening of the spectrum in the frequency domain. Low values of $\tau$ produce large broadenings, whereas large values of $\tau$ show sharp spectra.

The parameters of the excitonic Hamiltonian describing the energy transfer in perylene bisimide aggregates were taken from Ref. 74. The parameters for the S1/S2 vibronic Hamiltonian of pyrazine were taken from Ref. 91, employing the data calculated with multi-reference CI. All spectra were shifted by 35609 cm*-1* to match the experimental vertical excitation energy. As in our previous works,39, 40 the potential energy surface for ethylene was taken from Ref. 96 by expressing it as a quartic force-field in terms of the Cartesian-based normal modes employing the PyPES library.97

3.2 Excitonic dynamics of perylene-bisimide aggregates

We first apply our TD-DMRG theory to the simulation of the photochemical properties of a perylene bisimide-based dye, namely PBI-I. Perylene-bisimide derivatives can form molecular aggregates through $\pi$-stacking interactions. Upon absorption of UV light, the excitation is transferred efficiently between the constituents of the aggregate. As models for light-harvesting compounds, PBI-I aggregates have been extensively studied both computationally98, 89, 90, 99 and experimentally.100, 101

PBI-I features 400 vibrational modes and a full quantum treatment is prohibitive already for the isolated monomer in its electronic ground state. A simplified form of the vibronic Hamiltonian of molecular aggregates is obtained within the excitonic approximation. First, the electronic states of the aggregate are described within a diabatic representation, in which the diabatic states of the aggregate are assumed to coincide with the electronic states of the isolated monomers. Then, couplings between vibrations localized on different monomers are neglected. In the following, we label the ground state of the $i$-th monomer as $\alpha_{0}^{(i)}$ and its $n$-th excited state as $\alpha_{n}^{(i)}$. The vibronic Hamiltonian can be modeled as the so-called Holstein Hamiltonian

[TABLE]

where the electronic ONVs

[TABLE]

are given in the representation reported of Eq. (34). The normal modes of the $i$-th monomer are labeled as $\bm{q}_{i}$, $\mathcal{V}_{i}^{\text{GS}}$ is the ground-state PES of the $i$-th monomer and $\mathcal{V}_{i}^{\alpha}$ is its $\alpha$-th excited state PESs. Eq. (38) can be further simplified by approximating $\mathcal{V}_{i}^{\text{GS}}$ as a harmonic surface

[TABLE]

We further assume that the harmonic frequencies $\omega_{j}$ and the normal modes are the same for all the monomers in the aggregate ($q_{i,j}=q_{j}$). Finally, we approximate $\mathcal{V}_{i}^{\alpha}\left(\bm{q}_{i}\right)$ with the linear coupling model,

[TABLE]

$k_{\alpha,j}$ being the $j$-th component of the gradient of the PES associated to the $\alpha$-th electronically excited state. Following the protocol proposed in Ref. 74, we included in Eqs. (40) and (41) only 10 modes associated to the largest value of the gradient. To include excitation transfer effects, a coupling term $\mathcal{V}_{\text{c}}$ is introduced

[TABLE]

where the sum $\langle\alpha,\alpha^{\prime}\rangle$ includes only neighboring terms.

A second quantization form of Eqs. (40) and (41) is trivially obtained from the relations given in Eq. (32) for the nuclear part and Eq. (37) for the electronic one. We sorted the sites on the DMRG lattice as in Ref. 75 so that, by employing the notation introduced above, the sorted ONV is $|\>\alpha_{0}^{(1)}\alpha_{1}^{(1)}n_{1}^{(1)}\cdots\alpha_{0}^{(2)}\alpha_{1}^{(2)}n_{1}^{(2)}\ldots\>\rangle$. Hence, the sites are sorted based on the monomer on which they are centered. Within each monomer, the sites associated to an electronically excited state appear before the nuclear ones.

We first assess the stability of our MPS/MPO formulation of TD-DMRG with a single-site optimizer with respect to the integration time step $\Delta t$ and the bond dimension $m$. Figure 1 reports the absorption spectrum of a single PBI-I molecule obtained with $m$=10, 20, and 60, a time step of 1 fs and an overall propagation time of 4000 fs. The absorption spectrum is fully converged already for $m$=10 and increasing $m$ to 20 or 60 does not modify the overall bandshape. The same trend is observed also for the dimer, as reported in the upper panel of Figure 2, where a fully converged spectrum is obtained with $m$=20. We find that, for the dimer, the initial state of the propagation is a coherent superposition of excitations localized on each monomer. Also for the hexamer, which produces a DMRG lattice with 66 sites, a fully converged spectrum is obtained with $m$=20. This indicates that the convergence of TD-DMRG can be largely independent on the system size. This is not surprising, since the Holstein Hamiltonian given in Eq. (38) contains long-range interactions only between electronic degrees of freedom, while no long-range interactions are present for the vibrational part. These results suggest that the same holds true also for the time-dependent wavepacket during the whole propagation. We note that the $J$-coupling is classified as long-range because the sites are sorted in the lattice such that electronically excited states are not placed close one to the other. As will be discussed in more detail below, other sortings can be devised in which the $J$-coupling is short-ranged and that lead to a different convergence rate of DMRG.

For a given value of $m$, the accuracy of TD-DMRG depends on the step size $\Delta t$ chosen to discretize the propagation. The Lie-Trotter factorization of Eq. (18) is valid in the limit of short time steps. Therefore the accuracy is expected to deteriorate with increasing $\Delta t$. The error introduced in the Lie-Trotter factorization is quadratic in the time step $\Delta t$, proportional to the commutator of the local operators and can, in principle, be estimated.102, 103 Here, we only compare the spectra obtained with different time steps. Figure 3 reports the real (upper panel) and imaginary (lower panel) part of the time-dependent autocorrelation function of the PBI-I dimer for values of the time-step between 1 and 16 fs. For the first 600 fs of the propagation, the autocorrelation function is the same independent of the time step, suggesting that the Lie-Trotter approximation is reliable and time steps that are significantly larger than the ones employed in other MCTDH and TD-DMRG studies on this system can be safely employed.89, 75 We note that the time step affects the numerical Fourier transform of the autocorrelation function. The maximum frequency that can be detected with a numerical Fourier transform is $1/\Delta t$,104 and therefore, a large value for the time step impedes the calculation of the absorption spectrum in a wide energy range.

The graphical representation of the autocorrelation function for the hexamer obtained with different time steps is reported in Figure 4. We note that larger differences between the results obtained with different time steps are observed in this case. The autocorrelation function obtained with $\Delta t$=8 fs agrees with the one obtained with $\Delta t$=1 fs only up to 300 fs, and then they diverge. The error introduced in the Lie-Trotter factorization is proportional to the sum of all the possible commutators between the terms into which the Hamiltonian is partitioned. The number of these terms increases with system size, and therefore, a smaller time step is expected to be required for larger DMRG lattices.

We now investigate the reliability of the two-site propagation scheme for the dimer of PBI-I. As already mentioned in the theoretical section, the main advantage of a two-site TD-DMRG algorithm over a single-site one is the possibility of adapting dynamically the bond dimension based on the accuracy that is required in the representation of the wavefunction. We calculate the absorption spectrum of the dimer using different values for the threshold $\eta$ on the truncation error. We report in the upper panel of Figure 5 a graphical representation of the bond dimension $m$ obtained with two-site TD-DMRG for different values of $\eta$ ranging from 10*-5* to 10*-9*. For each value of $\eta$, the bond dimension $m$ increases continuously with time. As highlighted in Ref. 67, the entanglement of a quantum system grows during a non-equilibrium time evolution, and therefore, the bond dimension $m$ required to represent the wavefunction with a given accuracy (that is proportional to $\eta$) increases as well. Our calculations suggests that this is also the case for the PBI-I dimer. We note that, even if a steep increase for $m$ is observed by reducing $\eta$ from 10*-5* to 10*-9*, the choice of $\eta$ will have a negligible effect on the autocorrelation function (lower panel of Figure 5). This suggests that a rather high value of $\eta$ can be sufficient to obtain a fully-converged autocorrelation function, even if obtained for non-converged wavefunctions.

We recall that the autocorrelation function was obtained from the overlap of the time-dependent wavefunction $|\>\Psi(t)\>\rangle$ with the initial state $|\>\Psi(0)\>\rangle$. The latter state is a product state that can be represented with an MPS with $m$=1. It is, therefore, possible that the overlap with such a compact MPS converges faster than the representation of the full-dimensional wavefunction. This hypothesis agrees with the analysis reported in Ref. 105, where it was shown that even if, for fixed $m$, the truncation error grows with increasing time in single-site TD-DMRG, the calculation of averaged quantities is stable.

Besides the autocorrelation function, the population of each excited state as a function of time is easily extracted from TD-DMRG propagations through the average of the number operator $\mathcal{N}_{\alpha}$, defined as

[TABLE]

The number operator is local and can be easily encoded as an MPO. The dynamics of the population of excited states is particularly interesting if, as an initial state, we consider an excitation localized on one terminal molecule of the aggregate. In this case, tracking the average of $\mathcal{N}_{\alpha}$ enables one to calculate the rate of the transfer of excitation along the aggregate. We show in Figure 6 the population of the two terminal monomers of the PBI-I hexamer, $\langle\mathcal{N}_{1}\rangle$ and $\langle\mathcal{N}_{6}\rangle$, as a function of time $t$. The first discrepancies between the results obtained for different values of $m$ are observed after 25 fs, the population of the first site being overestimated with $m$=10. The populations obtained with with $m$=20 and 60 agree qualitatively and the results obtained with $m$=60 are identical to the $m$=100 ones.

The ordering of the DMRG lattice given in Ref. 75, which we have employed so far, is not unique and can be expected to affect the convergence of TD-DMRG. To assess the relevance of the ordering, we analyze the population dynamics with an alternative ordering obtained as $|\>\alpha_{0}^{(1)}\alpha_{1}^{(1)}\alpha_{0}^{(2)}\alpha_{1}^{(2)}n_{1}^{(1)}\cdots n_{1}^{(2)}\cdots\>\rangle$. In practice, all the electronic degrees of freedom are placed before the vibrational ones. In the Holstein Hamiltonian, two interaction terms are present, the hopping term, given in Eq. (42), and the linear vibronic term of Eq. (41). The relative strengths of these coupling terms affects the efficiency of the two orderings. We show in Figure 7 the population of the excited monomer, at one end of the aggregate, as a function of time for the two orderings described above. If the electronic sites were put close to each other, converged results had been obtained with $m$=20 up to 200 fs whereas, with the other ordering, deviations had been observed already at 100 fs. This suggests that the $J$ coupling is, in this specific case, predominant and that the electronic sites are more strongly entangled than the vibrational ones. Clearly, a quantitative analysis of the entanglement would require the calculation of the von Neumann entropy.86, 88, 106 However, such an analysis is beyond the scope of the present work.

3.3 Excited-state dynamics of pyrazine

We now calculate with TD-DMRG the low-energy vibrationally-resolved electronic spectrum of pyrazine for the S2 $\leftarrow$ S0 excitation. The S2 excited state of pyrazine is strongly coupled to S1 by non-adiabatic couplings, and therefore, an accurate absorption spectrum is obtained only with a vibronic Hamiltonian treating simultaneously these two states. For this work, we chose the vibronic Hamiltonian derived in Ref. 107 with the updated parameters reported in Ref. 91,

[TABLE]

where $\mathcal{V}_{S_{1}}(\bm{q})$ and $\mathcal{V}_{S_{2}}(\bm{q})$ are the potential energy surfaces of the S1 and S2 electronic states, respectively, and are obtained as a second-order Taylor expansion around the equilibrium geometry of the ground state. The off-diagonal coupling $\mathcal{V}_{12}(\bm{q})$ is also expressed as a second-order Taylor expansion. The Hamiltonian of Eq. (44) has become a standard benchmark system to probe the efficiency of new semiclassical108, 109, 110 and quantum91, 18, 23, 111, 112 dynamics algorithms. A full-dimensional quantum simulation of the dynamics of pyrazine, including all 24 vibrational degrees of freedom treated as independent and at the same level of accuracy is prohibitive for most quantum dynamics algorithms available in the literature, and reduced-dimensionality schemes must be adopted.91, 18, 112 Alternatively, the computational costs can be reduced by adopting ML-MCTDH.23 In the present work, we demonstrate how TD-DMRG can be applied to the simulation of the full-dimensional quantum dynamics of pyrazine without the need of introducing any hierarchical treatment of the degrees of freedom.

We begin our analysis with the simulation of the dynamics of pyrazine upon photoexcitation for a reduced-dimensionality model including only four modes, namely $\nu_{10a}$, $\nu_{6a}$, $\nu_{1}$, and $\nu_{9a}$ (following the notation of Ref. 91), which are the most relevant modes governing the dynamics. In the occupation number vector representation, the initial state of the propagation is $|\>0_{S_{1}}1_{S_{2}},0000\>\rangle$ since, under the Franck-Condon approximation, the molecule is excited to the S2 electronic state without modifying its vibrational wavefunction.

We report in Figure 8 the real and imaginary part of the autocorrelation function $C(t)=\langle\>\Psi(t)\mid\Psi(0)\>\rangle$ obtained by TD-DMRG propagation with different values of $m$. In all calculations we chose $N_{\text{max}}=10$ and a time step $\Delta t$ of 0.5 fs. For the first 50 fs of the propagation, the autocorrelation function converged already with $m$=10 with negligible variations observed by increasing $m$ up to 100. Starting from $t=50$ fs, small deviations are observed and only with $m$=20 converged results are obtained. We report in Figure 9 the vibronic spectrum obtained from these propagations. We include broadening effects with an exponential damping function with $\tau$=50 fs*-1* and 200 fs*-1*. For the lowest-resolution spectrum ($\tau$=50 fs*-1*) significant discrepancies are visible between the spectra obtained with $m$=10 and higher values of $m$. By increasing the resolution ($\tau$=200 fs*-1*), slight discrepancies between the spectra are observed also between the $m$=20 and $m$=40 spectra. However, no additional changes emerge by further increasing $m$.

In Figure 10, we compare the real part of the time-dependent autocorrelation function obtained for $m$=20 and different values for the integration time step. For the first 200 fs of the propagation, we note that the integration time can be increased up to 2 fs without a significant loss of accuracy, whereas strong deviations are observed by setting $\Delta t$ to 4 fs. Hence, a smaller time step is sufficient compared to the one needed to reach convergence for the previous example. This could be due to the fact that the Hamiltonian contains both short- and long-range interactions, and therefore, its MPO representation requires a higher bond dimension to maintain accuracy. Applying this MPO on an MPS increases the bond dimension faster than in the previous example, and therefore, a smaller time step is needed to ensure that the wavefunction remains on the manifold of the MPSs with a given value of $m$. Nevertheless, the time step is still higher than the one employed in other TD-DMRG formulations.75

We conclude our analysis of this 4-mode reduced-dimensionality Hamiltonian of pyrazine by investigating the two-site integration scheme. Figure 11 shows the graphical representation of the real part of the time-dependent autocorrelation function obtained with $N_{\text{max}}$=10 and different values of the truncation error $\eta$. With values of $\eta$ larger than 10*-10*, the decay of the autocorrelation function is slower than with $\eta$ smaller than 10*-10*, with constructive interference approximately every 50 fs. Comparing the results reported in Figure 11 with Figure 8 we note that a truncation error of $10^{-11}$ is sufficient to reach converged results.

Considering the absorption spectrum of pyrazine based on the 24-modes full-dimensional vibronic Hamiltonian, we report in the upper panel of Figure 12 the vibronic spectrum obtained with different values for the bond dimension $m$ employing a damping factor $\tau$ of 50 fs*-1*. We note that, despite the strong increase in system size (the DMRG lattice is now composed by 26 sites), the overall bandshape has already been qualitatively reproduced for 20 renormalized block states. By increasing $m$ to 40 and 60, only slight deviations are observed in the low-energy tail of the higher-energy band. By increasing the resolution of the spectrum through a damping factor of 200 fs*-1*, larger differences are observed as can be seen in the lower panel of Figure 12. For the lower-energy band, an overall blue-shift is observed for $m$=20. Nevertheless, the spectrum obtained for $m$=60 is superimposable to that for $m$=100. A similar effect is observed for the more intense S2 $\leftarrow$ S0 band. Increasing $m$ to 60 corrects the band position, but does not provide a fully converged spectrum, since a significant redistribution of the different bands of the vibronic spectrum is observed by further increasing $m$ to 100.

3.4 Imaginary time propagation: vibrational states of C2H4

The range of application of TD-DMRG can be extended by generalizing it to imaginary-time propagations. For example, the time evolution of thermal ensembles can be obtained within the so-called ancillla state approach by finite-time imaginary time evolution of a matrix product state.113, 32, 114 This approach, which was been generalized to an MPS/MPO formulation,115 was applied to the calculation of finite-temperature spectra of the aggregate studied at the beginning of this section.75 Imaginary-time propagations are also an alternative route to the sweep-based scheme to optimize MPSs. In fact, the following property holds:

[TABLE]

where $|\>\Psi_{\text{MPS}}\>\rangle$ is a random MPS and $|\>\Psi_{0}\>\rangle$ is the MPS representation of the ground state of the Hamiltonian $\mathcal{H}$. Eq. (45) is the basis of Diffusion Monte Carlo methods116 that solve the differential equation stochastically to avoid an uncontrolled increase in the basis size. In vibrational calculations, imaginary-time propagation of the MCTDH equations was applied to the calculation of the vibrational ground and excited states of several strongly anharmonic systems.93, 94, 95, 117 As for real-time propagations, the computational costs of MCTDH can be limited with TD-DMRG.49 Even if the sweep-based optimization is in most cases efficient enough, imaginary-time propagation can be easily generalized to excited states through,

[TABLE]

where $\omega$ is an energy-shift parameter and $t$ is squared to match the units of the squared Hamiltonian. It can be easily shown92 that the limiting state $|\>\Psi_{\omega}\>\rangle$ is the excited state of $\mathcal{H}$ with energy closest to $\omega$. Time-evolution governed by the propagator $e^{-t^{2}\left(\mathcal{H}-\omega\right)^{2}}$ can be simulated based on the algorithm described in Section 2 by employing the MPO representation of $\left(\mathcal{H}-\omega\right)^{2}$ instead of the one of $\mathcal{H}$. As discussed in our previous work,40 the calculation of matrix elements of powers of the Hamiltonian is more efficient and accurate within an MPS/MPO framework.

Here, we apply imaginary-time TD-DMRG to the optimization of the vibrational levels of ethylene, which we considered in previous work in the context of vDMRG.39, 40 We employ the same potential energy surface as in Ref. 40, taken from Ref. 96 and converted to a fourth-order Taylor expansion in Cartesian normal modes by neglecting first-order $\pi$-$\pi$ coupling terms in the kinetic energy operator. The upper panel of Figure 13 reports the results obtained with N${}_{\text{max}}$=6 and $m$=10 for the non-shifted imaginary-time propagator obtained from Eq. (45). The energy converges within three sweeps, and therefore, the algorithm is as efficient as the standard time-independent optimization scheme.39 In Ref. 62, it was noted that the back-propagation step of the zero-site tensor (the terms involving the $\mathcal{P}_{i}^{(2)}$ projector operators in Eq. (14)) is not needed in imaginary-time propagations, in which the wavefunction for finite values of $t$ is not of interest and only the asymptotic limit of the propagation is needed. As reported in the upper panel of Figure 13, skipping the back-propagation step improves the rate of convergence of the energy, therefore confirming the previous observation. Based on these results, if not otherwise specified, imaginary-time evolutions are performed without back-propagating the zero-site tensor in the following. As reported in the lower panel of Figure 13, the convergence is fast also for higher values of $m$. Moreover, the asymptotic value of the energy that corresponds to the zero-point vibrational energy of ethylene agrees with the data reported in our previous work.39, 40

We extend our previous analysis to low-lying vibrational excited states of ethylene through the imaginary-time evolution with the square-and-shift operator defined in Eq. (46). We set the energy shift to 11800 cm*-1* for the first state, to 11900 cm*-1* for the second one, and to 11930 cm*-1* for the third one. As shown in Figure 14, in all cases the optimization converges smoothly within several sweeps to the energy of the excited state, which is in good agreement with our reference data reported in Ref. 40.

We have already noted in the context of our recently introduced time-independent energy-specific formulation of vDMRG40 that an estimate of the excitation energy is required to choose the energy shift to ensure that the optimization converges to the desired excited state. The choice of $\omega$ is particularly critical in regions with a high density of states, in which root-flipping events can occur. In time-independent simulations, we employed a root-homing criterion (also known as maximum overlap118, 119, 120, 121) to consistently follow, within an iteration of the optimization, the solution with the largest overlap with a predetermined MPS. Root-homing techniques cannot be adopted in imaginary-time propagation because only one state is available in each sweep. Nevertheless, it is possible to exploit a partially converged time-independent DMRG calculation to build a reliable guess for the energy and the wavefunction of the target state to be then employed in an imaginary-time propagation scheme.

4 Conclusions

In the present work, we extended the tangent space-based formulation of TD-DMRG62 and applied it to the simulation of the vibrational and vibronic dynamics of molecular systems. Due to the intrinsic complexity of the Hamiltonians describing vibrational and vibronic motion in molecules, we chose an MPS/MPO approach62 considered so far only for studying the dynamics of one-dimensional Hamiltonians.75, 70 Our approach has several advantages over alternative TD-DMRG schemes. First, it supports Hamiltonians with general structures containing both short- and long-range interactions. Moreover, it leads to a set of coupled differential equations that can be integrated exactly. This is achieved by constraining the time evolution of the wavefunction on the manifold composed by MPSs with a given rank. We extend the range of applicability of our TD-DMRG formulation by two alternative strategies to dynamically adapt the dimension of the MPS during the propagation. The first one employs a subspace-expansion scheme73 and the other one relies on a two-site integration.81 We first investigated the reliability of our implementation on the excitonic Hamiltonian describing the excitation energy transfer in a perylene bisimide aggregate. The application of TD-DMRG for a 24-mode vibronic Hamiltonian, describing the S1 and S2 electronic states of pyrazine, demonstrated the reliability of the algorithm when applied to large systems, which are difficult to access with current quantum dynamics approaches. We further generalized our TD-DMRG formulation to imaginary-time propagation, which represents an alternative to the standard, sweep-based scheme for optimizing ground and excited states as MPSs.

On the one hand, the results of the simulations suggest that the bond dimension of an MPS should be increased constantly during the propagation to keep the truncation error, i.e. accuracy of the MPS representation, fixed for the whole propagation. On the other hand, our tests show that employing a fixed value for the bond dimension $m$, and therefore letting the truncation error grow with time, enables one to obtain accurate results. This is probably due to the fact that absorption spectra are mainly determined by short-time dynamics and, therefore a small value of the bond dimension is sufficient to represent them accurately. However, additional investigations are required to check if this hypothesis is correct.

Due to its flexibility, our framework can be easily extended to the electronic Hamiltonian, the only additional hurdle being the inclusion of the spin symmetry. This hurdle can be overcome based on the strategies developed for standard time-independent DMRG.76, 122 The resulting theory represents a cost-efficient alternative to time-dependent CI123, 124, 12, 125 and allows one to simulate the electron dynamics for systems that are currently targeted by approximate approaches such as real-time time-dependent density functional theory.126, 127, 128 In the present formulation of the TD-DMRG algorithm, the local basis set is fixed during the whole propagation. A simultaneous optimization of both the CI coefficients, expressed as MPS, and of the local basis in a self-consistent-field fashion could enhance the efficiency of TD-DMRG. If applied to the electronic Hamiltonian, the resulting theory will constitute a cost-effective alternative to the time-dependent complete active space selft-consistent-field algorithm recently proposed in the literature,129, 130, 131 applied only to atoms and few-electron systems so far. For vibrational Hamiltonians, the coupled optimization of CI coefficients and local basis functions will lead to a DMRG approach analog of MCTDH.132, 70

Besides being a cost-effective alternative to MCTDH, TD-DMRG can also drive the definition of strongly-interacting degrees of freedom for ML-MCTDH.24, 25 For electronic structure theory, entanglement-based measures86, 133, 106, 134, 135, 136 obtained from DMRG calculations have been demonstrated to be reliable metrics to detect strong-correlation between orbitals. When generalized to vibrational wavefunctions,137, 138, 139 the same metrics could support the identification of strongly interacting vibrational degrees of freedom.

All results reported in the present work are obtained from dynamics performed on fitted potential-energy surfaces obtained as a Taylor expansion around a reference geometry. This is a major limitation, especially for long-time dynamics in which the wavefunction explores regions of the potential energy surface that cannot be efficiently represented by a power series in terms of some coordinate set. Some attempts are available in the literature to perform on-the-fly dynamics with MCTDH.140, 141 These algorithms find a natural extension in DMRG through the framework discussed in this work.

{acknowledgement}

This work was supported by ETH Zürich through the ETH Fellowship No. FEL-49 18-1.