A Further Comparison of MPS and TTNS for Nonadiabatic Dynamics of Exciton Dissociation

TL;DR

This study systematically compares MPS and TTNS tensor network methods for simulating exciton dissociation, demonstrating that increasing bond dimensions ensures both methods converge to accurate solutions, thus validating their reliability.

Contribution

It provides a comprehensive comparison of MPS and TTNS, highlighting the importance of bond dimension scaling for accurate quantum dynamics simulations.

Findings

Discrepancies in previous results are mainly due to insufficient bond dimensions.

Increasing bond dimensions reduces differences in electronic state occupancy to below 10%.

Optimized TTNS structures further improve accuracy and efficiency.

Abstract

Tensor networks, such as matrix product states (MPS) and tree tensor network states (TTNS), are powerful ans\"atze for simulating quantum dynamics. While both ans\"atze are theoretically exact in the limit of large bond dimensions, [J. Chem. Theory Comput. 2024, 20, 8767-8781] reported a non-negligible discrepancy in its calculations for exciton dissociation. To resolve this inconsistency, we conduct a systematic comparison using Renormalizer, a unified software framework for MPS and TTNS. By revisiting the benchmark P3HT:PCBM heterojunction model, we show that the observed discrepancies arise primarily from insufficient bond dimensions. By increasing bond dimensions, we reduce the relative difference in occupancy for weakly populated electronic states from up to 60% towards the end of the simulation to less than 10% and the absolute difference from 0.05 to 0.005. We also discuss the impact of tensor network structures on accuracy and efficiency, with the difference further reduced by an optimized TTNS structure. Our results confirm that both methods converge to numerically exact solutions when bond dimensions are adequately scaled. This work not only validates the reliability of both methods but also provides high-accuracy benchmark data for future developments in quantum dynamics simulations.

Tables

Table 1. Table A1: The on-site energy of the CS states.

$n$	${ϵ}_n^{\text{CS}}$ (meV)
1	0.0
2	33.6
3	47.4
4	56.0
5	61.8
6	65.7
7	68.4
8	70.0
9	70.9
10	71.2
11	71.1
12	70.5
13	69.5

Table 2. Table A2: The frequencies and coupling constants of the F F modes and the OT local vibrational modes. The units are meV.

$l$	${\omega }_{F,l}$	${\omega }_{\text{OT},l}$	$g_{F,l}$	$g_{\text{OT},l}$	$g_{\text{OT},l}^{\prime }$
1	200.025	401.283	45.246	7.017	4.035
2	184.269	397.773	65.701	-0.077	2.921
3	177.853	182.714	-40.280	-67.849	-129.712
4	141.11	178.531	-17.511	57.668	46.885
5	93.952	134.550	28.026	-40.145	-32.908
6	79.933	111.848	-13.629	11.68	36.591
7	55.892	42.621	-23.732	-10.784	-20.211
8	33.264	18.316	9.86	-12.309	-7.77

Equations

$$\hat{H}={\hat{H}}^{\text{e}}+{\hat{H}}^{\text{vib}}+{\hat{H}}^{\text{e-vib}}.$$

$${\hat{H}}^{\text{e}}={\hat{H}}_{\text{diag}}^{\text{e}}+{\hat{H}}_{\text{off-diag}}^{\text{e}}.$$

$${\hat{H}}_{\text{diag}}^{\text{e}}={ϵ}^{\text{LE}}\sum _n^{N_o}|{\text{LE}}_n⟩⟨{\text{LE}}_n|+\sum _n^{N_o}{ϵ}_n^{\text{CS}}|{\text{CS}}_n⟩⟨{\text{CS}}_n|.$$

$${\hat{H}}_{\text{off-diag}}^{\text{e}}=J\sum _n^{N_o-1}|{\text{LE}}_n⟩⟨{\text{LE}}_{+1}|+t\sum _n^{N_o-1}|{\text{CS}}_n⟩⟨{\text{CS}}_{n+1}|+\lambda |{\text{LE}}_1⟩⟨{\text{CS}}_1|+\text{H.C.}$$

$${\hat{H}}^{\text{vib}}={\omega }_R{\hat{b}}_R^{†}{\hat{b}}_R+\sum _l^{N_f}{\omega }_{F,l}{\hat{b}}_{F,l}^{†}{\hat{b}}_{F,l}+\sum _n^{N_o}\sum _l^{N_m}{\omega }_{\text{OT},l}{\hat{b}}_{\text{OT},nl}^{†}{\hat{b}}_{\text{OT},nl}.$$

$${\hat{H}}^{\text{e-vib}}={\hat{H}}_R^{\text{e-vib}}+{\hat{H}}_F^{\text{e-vib}}+{\hat{H}}_{\text{OT}}^{\text{e-vib}}.$$

$${\hat{H}}_R^{\text{e-vib}}=g_R\left({\hat{b}}_R^{†}+{\hat{b}}_R\right)|{\text{CS}}_1⟩⟨{\text{CS}}_1|+g_R^{\prime }\left({\hat{b}}_R^{†}+{\hat{b}}_R\right)\left(|{\text{LE}}_1⟩⟨{\text{CS}}_1|+\text{H.C.}\right),$$

$${\hat{H}}_F^{\text{e-vib}}=\sum _n^{N_o}\sum _l^{N_f}{g}_{F,l}\left({\hat{b}}_{F,l}^{†}+{\hat{b}}_{F,l}\right)|{\text{CS}}_n⟩⟨{\text{CS}}_n|,$$

$${\hat{H}}_{\text{OT}}^{\text{e-vib}}=\sum _n^{N_o}\sum _l^{N_f}\left[g_{\text{OT},l}\left({\hat{b}}_{\text{OT},nl}^{†}+{\hat{b}}_{\text{OT},nl}\right)|{\text{CS}}_n⟩⟨{\text{CS}}_n|+g_{\text{OT},l}^{\prime }\left({\hat{b}}_{\text{OT},nl}^{†}+{\hat{b}}_{\text{OT},nl}\right)|{\text{LE}}_n⟩⟨{\text{LE}}_n|\right].$$

$$|\Psi ⟩=\sum _{\{a\},\{\sigma \}}A[1{]}_{{\Lambda }_1,a_1}^{{\sigma }_1}A[2{]}_{{\Lambda }_2,a_2}^{{\sigma }_2}\cdots A[N{]}_{{\Lambda }_N,a_N}^{{\sigma }_N}|{\sigma }_1{\sigma }_2\cdots {\sigma }_N⟩.$$

$$\sum _{{\Lambda }_i,{\sigma }_i}{A}^{†}[i{]}_{{\Lambda }_i,a_i^{\prime }}^{{\sigma }_i}A[i{]}_{{\Lambda }_i,a_i}^{{\sigma }_i}={\delta }_{a_i^{\prime }{a}_i}.$$

$$|\Psi ⟩=\sum _{{\Lambda }_r}A[r{]}_{{\Lambda }_r}\prod _{\{j;a_j\in {\Lambda }_r\}}|a[j{]}_{a_j}⟩,$$

$$|a[j{]}_{a_j}⟩=\sum _{{\Lambda }_j,{\sigma }_j}A[j{]}_{{\Lambda }_j,a_j}^{{\sigma }_j}|{\sigma }_j⟩\prod _{\{k;a_k\in {\Lambda }_j\}}|a[k{]}_{a_k}⟩$$

$$A[r{]}_{{\Lambda }_r^{\prime },a_j}\stackrel{\text{SVD}}{=}\sum _l^{|a_j|}{U}_{{\Lambda }_{r}l}{s}_l{V}_{l{a}_j}\approx \sum _l^{M^{\prime }}{U}_{{\Lambda }_{r}l}{s}_l{V}_{l{a}_j}.$$

$$1-\text{absolutevalue}{⟨{\Psi }^{\prime }|\Psi ⟩}^2=\sum _{l=M^{\prime }+1}^{|a_j|}{s}_l^2.$$

$$S=-\text{Tr}{\rho }_X\ln {\rho }_X=-\text{Tr}{\rho }_Y\ln {\rho }_Y$$

$${\rho }_X={\text{Tr}}_Y|\Psi ⟩⟨\Psi |=\sum _{{\Lambda }_r^{\prime },a_k,a_k^{\prime }}{A}^{†}[r{]}_{{\Lambda }_r^{\prime },a_k^{\prime }}A[r{]}_{{\Lambda }_r^{\prime },a_k}|a[k{]}_{a_k}⟩⟨a[k{]}_{a_k^{\prime }}|$$

$$S=-\sum _l{s}_l^2\ln s_l^2$$

$$\text{Error}=\frac{{\sum }_{t=150}^{200}\text{absolutevalue}n(t)-n_{\text{ref}}(t)}{200-150}.$$

$${\rho }^{\prime }=\rho +\epsilon \exp -\rho /\epsilon$$

Full text

1 Introduction

Tensor networks have emerged as powerful tools for simulating quantum dynamics. These networks approximate high-order wavefunction tensors as contractions of lower-order tensors, enabling accurate simulations of complex systems orus2019tensor, banuls2023tensor, larsson2024tensor. The accuracy of this approximation is controlled by the size of the low-order tensors, also known as the bond dimension. In the limit of large bond dimensions, tensor network methods converge to the exact solution. Among tensor network methods, matrix product states (MPS) are the foundation of the time-dependent density matrix renormalization group (TD-DMRG) schollwock2011density, paeckel2019time, ren2022time, haegeman2011time, lubich2015tt, lubich2015mctdh, haegeman2016unifying. TD-DMRG has been successfully applied to study systems such as charge transport in organic semiconductors li2021general, non-radiative decay in molecular aggregates wang2023minimizing, and photoinduced ultrafast vibration-coupled electron transfer reactions wang2025numerically, among others borrelli2017simulation, greene2017tensor, ren2018time, ma19, baiardi2019large, borrelli2021finite, peng2023studies. Meanwhile, tree tensor network states (TTNS) form the basis of the multi-layer multiconfiguration time-dependent Hartree method (ML-MCTDH) shi2006classical, wang2015multilayer, which has found wide-spread applications wang2008coherent, wang2010coherent, vendrell2011multilayer, weike2022multi, ke2023tree, reddy2024intramolecular, zheng2025ml such as the study of singlet fission schroder2019tensor and exciton migration binder2018conformational, binder2020first, di2020quantum, popp2021quantum, mondelo2022quantum, brey2024coherent.

While TD-DMRG and ML-MCTDH originated in different research communities meyer1990multi, manthe1992wave, white1992density, white1993density, wang2003multilayer, white2004real, daley2004time, vidal2004efficient, feiguin2005time, it is now well-established that both rely on tensor networks and share fundamental similarities larsson2024tensor. Recent efforts have focused on directly comparing these methods to understand how their different tensor network structures affect simulation accuracy and efficiency gunst2018t3ns, larsson2019computing, ren2022time, li2024optimal, larsson2024tensor, chen2025tree. In principle, both methods should converge to the exact result when the bond dimension is sufficiently large. It is often assumed that TTNS converges faster than MPS with respect to the bond dimension, while MPS offers a better computational scaling in terms of the bond dimension, due to its simpler tensor structure.

However, an intriguing discrepancy emerged in a recent benchmark study dorfner2024comparison. The authors performed an extensive and careful comparison of the two methods with standard implementations, and found that both methods agreed perfectly in most scenarios. Their analysis also connected the observed dynamics with patterns of entanglement entropy and the underlying tensor network structures. For exciton dissociation at poly(3-hexylthiophene):[6,6]-phenyl-\chC61-butyric acid methyl ester (P3HT:PCBM) heterojunctions, TD-DMRG and ML-MCTDH showed small but non-negligible differences at the long time limit. This inconsistency is troubling because both methods are frequently used as reference standards for validating and benchmarking other quantum dynamics approaches peng2022exciton, xu2022taming, li2023efficient, liu2024benchmarking. The comparison is complicated by the use of different software implementations mctdh:MLpackage, fishman2022itensor, making it challenging to isolate the origin of these differences. Very recently, Lindoy and Rungger et al. systematically investigated this system using both multiset and singleset MPS/TTNS ansätze based on the pyTTN package lindoy2025pyttn. They demonstrated that the multiset ansätze yields highly accurate and consistent results, while the singleset ansätze exhibits slow bond dimension convergence. Their analysis suggests that insufficient bond dimension is the most likely source of the discrepancy reported in the benchmark study dorfner2024comparison.

In this work, we revisit the P3HT-PCBM model using the Renormalizer package, which provides a unified framework for both MPS and TTNS li2024optimal, renormalizer. Through extensive benchmark using the projector splitting (PS) algorithm based on the Time Dependent Variational Principle (TDVP) haegeman2011time, lubich2015tt, lubich2015mctdh, haegeman2016unifying, lindoy2021time, lindoy2021time2, ceruti2021time, we conclusively confirm that the observed discrepancies stem from insufficient bond dimensions. In this context, TDVP-PS with MPS corresponds to TD-DMRG, while TDVP-PS with TTNS corresponds to a variant of ML-MCTDH known as Projector Splitting Integrator ML-MCTDH (PSI-ML-MCTDH).

Additionally, guided by the entanglement entropy calculated from TTNS simulations, we propose two new structures for MPS and TTNS, respectively, tailored for this exciton dissociation model. We find that the new TTNS structure is more efficient than the previous MPS and TTNS structures, while the new MPS structure results in higher error.

Our results demonstrate that both methods can achieve numerically exact results, reinforcing their reliability for the simulation of complex quantum dynamics.

2 Methodology

2.1 The Exciton Dissociation Model

In this section, we present the model for the numerical comparison of MPS and TTNS. The model describes the exciton dissociation at the interface of fullerene molecules and a linear chain of $N_{o}=13$ oligothiophene (OT) molecules, which is originally introduced in the reference paper tamura2013ultrafast, huix2015concurrent, dorfner2024comparison. The model can be considered as a simplified representation of the P3HT:PCBM heterojunction Dang2011P3HT. The reference paper provides a detailed chemical picture of the model. In Fig. 1, we present an interaction-based view of the model.

As shown in Fig. 1, each OT molecule is associated with two electronic states: a local excitation (LE) state and a charge-separated (CS) state. The CS state describes the interaction between the OT molecule and the fullerene cluster. The LE and CS states can hop to their nearest neighbors. Additionally, the LE and CS states at the first OT molecule (LE1 and CS1) can transit between each other. The system starts in the LE1 state due to initial excitation. Besides, there are three types of vibrations in the model, which are shown as filled spheres in Fig. 1. Each OT molecule has 8 local vibrational modes, and they are coupled to both its LE state and the CS state. The fullerene cluster has 8 vibrational modes, denoted as $F$, which couple to all CS states. There is an additional intermolecular vibrational mode, denoted as $R$, which specifically couples to the transition between LE1 and CS1.

The Hamiltonian of the system $\hat{H}$ can be decomposed into three parts

[TABLE]

Here, $\hat{H}^{\textrm{e}}$ represents the electronic energy and interactions, $\hat{H}^{\textrm{vib}}$ represents the vibration energy, and $\hat{H}^{\textrm{e-vib}}$ represents the electron-phonon interaction or the vibronic coupling. We next describe each term in detail.

The electronic Hamiltonian $\hat{H}^{\textrm{e}}$ is further split into diagonal and off-diagonal terms

[TABLE]

The diagonal part consists of on-site energies for the LE and CS states:

[TABLE]

The off-diagonal terms include hopping interactions between neighboring LE and CS states:

[TABLE]

where “H.C.” represents Hermitian conjugation. Following the reference work, the parameters used in this model are $\epsilon^{\textrm{LE}}=100$ meV, $J=100$ meV, $t=-120$ meV and $\lambda=-200$ meV. The values for $\epsilon^{\textrm{CS}}_{n}$ are listed in Table A1.

The vibrational Hamiltonian $\hat{H}^{\textrm{vib}}$ includes three types of vibrations: an intermolecular mode $R$, $N_{f}$ fullerene modes $F$ and $N_{o}\times N_{m}$ local OT vibrational modes. When there is no ambiguity, we refer to the local OT vibrations simply as “OT”. $\hat{H}^{\textrm{vib}}$ is then written as

[TABLE]

Here, $N_{f}=N_{m}=8$, $\omega_{R}=10$ meV, and the frequencies $\omega_{F,l}$ and $\omega_{\textrm{OT},l}$ are listed in Table A2 in Appendix A.

$\hat{H}^{\textrm{e-vib}}$ assumes linear vibronic coupling and can be decomposed to three groups, corresponding to three different types of vibrations:

[TABLE]

The intermolecular $R$ mode has both diagonal and off-diagonal couplings to the LE1 state and CT1 state

[TABLE]

where $g_{R}=-10/\sqrt{2}$ meV and $g^{\prime}_{R}=-30/\sqrt{2}$ meV. The fullerene modes $F$ couple to the on-site energy of the CS states

[TABLE]

where the coupling constants $g_{F,l}$ are included in Table A2. The OT modes couple to the on-site energy of both the LE and CS states

[TABLE]

The coupling constants $g_{\textrm{OT},l}$ and $g^{\prime}_{\textrm{OT},l}$ are included in Table A2. These coupling constants, including $F$ and OT vibrations, are assumed to be independent of $n$, the index of the OT molecule.

Overall, the exciton dissociation model consists of 26 electronic states and $1+8+8\times 13=113$ vibrational modes. While most of the interactions are short ranged, $\hat{H}^{\textrm{e-vib}}_{F}$ introduces long range interactions between fullerene and OT molecules.

2.2 MPS and TTNS

In this section, we briefly introduce the principles of MPS and TTNS, using the tensor network language and with an emphasis on the relationship between their accuracy and the entanglement entropy. For more detailed explanations, readers are encouraged to consult several excellent reviews and introductory literatures beck2000multiconfiguration, schollwock2011density, orus2014practical, chan2016matrix, hauschild2018efficient, ren2022time, larsson2024tensor.

MPS and TTNS are data structures that approximate high-order tensors by the contraction of low-order tensors. Let us consider a quantum system with $N$ degrees of freedom. For each degree of freedom, the corresponding primitive basis is denoted by $\ket{\sigma_{i}}$. The wavefunction of the system can be approximated by MPS and TTNS as

[TABLE]

Here $A[i]$ represents the low-order tensors. $\Lambda_{i}$ is a set of indices that connects to child nodes and $a_{i}$ connects to the parent node. If $\Lambda_{i}$ is an empty set, then $A[i]$ represents a leaf node. The contraction between the tensors occur according to the index $\Lambda_{i}$. While TTNS employs a general tree-like contraction topology, MPS utilizes a linear chain structure, making MPS a special case of TTNS where all tensors are arranged in one dimension. Eq. 10 represents a general TTNS and it reduces to a MPS if $\Lambda_{i}=\{a_{i-1}\}$.

Since tensor contraction operation is undirected, the concept of parent and children nodes in TTNS is arbitrary. In other words, we are free to choose the root of the tree for implementation or formal purposes and it does not affect the contraction result in Eq. 10. If a specific tensor is chosen as the root, we denote it as $A[r]$, where the index $a_{r}$ is not included in any $\Lambda_{i}$ and the dimension of $a_{r}$ is 1. As a result, the index $a_{r}$ is omitted in the following when no ambiguity arises.

In ML-MCTDH or three-legged tree tensor network gunst2018t3ns, there are entirely “virtual” nodes, which are not associated with any physical degree of freedom $\sigma$. For these virtual nodes, an auxiliary physical degree of freedom with a Hilbert space of dimension 1 can be assigned to ensure consistency with Eq. 10. Similarly, there are also nodes where the number of physical degree of freedom is greater than 1 and they can be combined so that formally only one physical degree of freedom presents, which is also known as the mode combination technique. For simplicity, we assume that at the root node the dimension of $\sigma_{r}$ is 1 and thus $\sigma_{r}$ can be ignored. If $\absolutevalue{\sigma_{r}}\neq 1$, we can decompose $A[r]^{\sigma_{r}}_{\Lambda_{r}}$ into two tensors. One that has the indices $\Lambda_{r}$ acts as the new root, and the other whose shape is $\absolutevalue{\sigma_{r}}\times\absolutevalue{\sigma_{r}}$ acts as a child to the new root. An auxiliary physical degree of freedom can then be assigned to the new root and $\absolutevalue{\sigma_{r}}$ becomes 1. In this work $|\sigma|$ or $|a|$ represent the dimension of the index.

Through sequential QR or singular value decomposition, the tensor network in Eq. 10 can be transformed into a “canonical” format. A tensor $A[i]$ is said to be canonical if it satisfies the following condition:

[TABLE]

Eq. 10 is considered canonical if all tensors $A[i]$, except for the root $A[r]$, are canonical. The root node is thus also referred to as the canonical center. Importantly, the canonical center can be moved to any node in the tree, analogous to how any node can serve as the root. In the canonical form, the wavefunction $\ket{\Psi}$ can be expressed as

[TABLE]

where $A[r]$ acts as the coefficient tensor, and $\{\ket{a[j]}\}$ forms an orthogonal basis set. Here each $\ket{a[j]}$ is the orthogonal basis set by one of the children

[TABLE]

The orthogonal relation $\braket{a[j]_{a^{\prime}_{j}}|a[j]_{a_{j}}}=\delta_{a^{\prime}_{j}a_{j}}$ can be derived from the canonical condition Eq. 11 through tree recursion.

The accuracy of the tensor network approximation is controlled by the dimension of the indices $\{a\}$. In MPS and general tensor network literature, this quantity is typically referred to as the bond dimension. In the context of ML-MCTDH, the same quantity is called the number of single-particle functions. In this work, we employ the term “bond dimension” following our previous conventions and use $M$ to denote this quantity.

In principle, the bond dimension can vary for each bond in the tree, and it can be adjusted dynamically during time evolution. In this work, we employ the same fixed bond dimension for every bond, primarily to simplify the setup. Since the initial state of the model introduced in Sec. 2.1 is a product state, Krylov subspace vectors $\hat{H}^{n}\ket{\Psi}$ are added to the product state with small coefficients to expand the bond dimension to the target value $M$.

Specifically, we typically add around 10 such vectors, included together with a small weight of $10^{-10}$, to ensure minimal perturbation while effectively expanding the state space. This enrichment is performed only once during the initial state preparation and is not repeated during the time evolution. We note that here the Krylov vector refers to the global wavefunction, and should not be confused with the Krylov vector method employed for solving the matrix exponential of the time evolution discussed at the end of Sec. 2.3.

The error introduced by the finite bond dimension can be quantified by the singular values of $A[r]$. The square of the singular values corresponds to the natural orbital population in MCTDH language. Suppose $a_{j}\in\Lambda_{r}$ and the other indices of the root node are denoted as $\Lambda^{\prime}_{r}=\Lambda_{r}\setminus\{a_{j}\}$, we can reshape $A[r]_{\Lambda_{r}}$ to a two dimensional tensor, where the first and the second indices correspond to $\Lambda^{\prime}_{r}$ and $a_{j}$ respectively. By performing singular value decomposition on $A[r]_{\Lambda^{\prime}_{r},a_{j}}$, we obtain the singular values $s_{l}$. Let $M^{\prime}$ be a integer smaller than $|a_{j}|$. By keeping only the first (largest) $M^{\prime}$ singular values, the dimension of $a_{j}$ is compressed to $M^{\prime}$

[TABLE]

Let $\ket{\Psi^{\prime}}$ be the wavefunction after compression. Since $U$ and $V$ are unitary matrices and the bases in Eq. 11 are orthogonal, the compression fidelity is $\absolutevalue{\braket{\Psi^{\prime}|\Psi}}^{2}=\sum_{l}^{M^{\prime}}s^{2}_{l}$. Assuming the wavefunction is normalized, $\sum_{l}^{|a_{j}|}s^{2}_{l}=1$. The compression error is the sum of the square of discarded singular values

[TABLE]

Thus, the ideal wavefunction for tensor network approximation are those whose singular values decay rapidly.

The efficiency of the tensor network compression is closely related to the bipartite von Neumann entanglement entropy $S$. In a tree tensor network, cutting an arbitrary bond (edge) $a_{k}$ divides the system degrees of freedom into two parts $X$ and $Y$. The bipartite von Neumann entanglement entropy $S$ is defined as

[TABLE]

where $\rho_{X}$ and $\rho_{Y}$ are the reduced density matrix of $X$ and $Y$, respectively.

In general, $S$ is difficult to calculate for many-body systems, because calculating $S$ involves diagonalizing the reduced density matrices to obtain the eigenvalues. However, the canonical form of tensor networks, as described in Eq. 12, provides an efficient way to calculate $S$. To do so, the canonical center is firstly moved to the bond $a_{k}$ that divides the sub-systems such that $a_{k}\in\Lambda_{r}$. The remaining indices are denoted as $\Lambda^{\prime}_{r}=\Lambda_{r}\setminus\{a_{k}\}$. Supposing $X$ is in the subtree of the bond $a_{k}$, its reduced density matrix $\rho_{X}$ is then

[TABLE]

By performing SVD on $A[r]$ as described earlier in Eq. 14, the entanglement entropy can be calculated by the singular values

[TABLE]

Note that SVD can also be employed to move the canonical center to the neighbouring nodes. Thus, by sweep the canonical center across the tree, we can obtain the entanglement entropy $S$ for each bond.

As indicated by Eq. 18, if the bond $a_{k}$ has bond dimension $M$, the maximumly possible entanglement entropy $S$ is $\ln M$, which occurs when $s_{l}=1/\sqrt{M}$. In other words, for a bipartite system with entanglement entropy $S$, the bond dimension must exceed $e^{S}$ in order to accurately describe the system. As a special case, if $S=0$, then $s_{l}=\delta_{l,0}$, and the required bond dimension is 1. Thus, for systems exhibiting strong entanglement, a large bond dimension is required to achieve an accurate simulation. The choice of tensor network structure, particularly the MPS ordering or more generally the TTNS tree structure, affects the bipartition of the system and consequently the entanglement entropy. Structures that minimize bipartite entanglement entropy are preferred, as they allow for more efficient simulations. Finding the globally optimal structure for a given Hamiltonian is believed to be a challenging problem. In practice, tensor network structures are often designed based on heuristics or the nature of the system’s interactions.

2.3 Tensor Network Structures

In this study, we explore four different tensor network structures, which are illustrated in Fig. 2. The first is a standard MPS, shown in Fig. 2(a), which takes a linear form. The first site of the MPS contains all electronic states, including 13 LE states and 13 CS states. This is followed by the $R$ site, and then 8 $F$ vibration sites. Finally, the local OT vibrations are appended molecule by molecule. The site ordering of the MPS is the same as the site ordering in the reference work dorfner2024comparison.

The second structure is a tree tensor network based on the reference work. As shown in Fig. 2(b), the root node has 3 children. The first child contains all 8 $F$ modes and the $R$ mode. The second child contains the 26 different electronic states. The third child contains all local OT vibrational modes. The OT vibrational modes are divided into two subtrees based on the vibration frequency: low frequency modes and high frequency modes. High frequency modes are defined as the modes whose $\omega_{\textrm{OT},l}>300$ meV. According to Table A2, there are two high frequency modes per OT molecule. The topology of this tree is identical to the one used in the reference work. In the tree structure, a few nodes that are close to the leaves in the OT subtree have three children. Since in this work we use fixed uniform bond dimension across the whole tree, such ternary nodes become too large and lead to an unnecessary computational bottleneck. To address this, we reduce the bond dimension associated with these nodes and their children to $\frac{1}{8}$ of the bond dimension of the entire tree. Since these nodes are close to the leaves, this reduction in bond dimension has minimal impact on accuracy. The bond dimension is the only difference between Tree and the TTNS structure employed in the reference work.

We also note that although the root node in Fig. 2(b) and Fig 2(d) has three children, it is not the bottleneck of the computation. This is because as the root of the tree, the total number of virtual indices it carries is three, which is comparable to other nodes that have two children and one parent. Therefore, a root node with three children ensures balanced computational cost over all nodes in the tree.

Next, we describe the tensor network structures proposed in this work, which are designed based on the calculation result from the Tree structure. The third structure, shown in Fig. 2(c), is another MPS with different configuration compared to Fig. 2(a). For clarity, in this paper, we use “Chain” and “ChainX” to denote the MPS configurations shown in Fig. 2(a) and Fig. 2(c), respectively. In ChainX, the LE state, CS state and the vibrations of each OT molecule are firstly grouped together. Then, the OT1 unit is placed at the center of the the chain, with OT units of even indices on the left side and OT units of odd indices on the right side. The $R$ mode is placed next to the OT1 unit and the F modes are positioned at the beginning of the chain. The structure is designed so that the degrees of freedom with the strongest entanglement entropy are placed in the center of the chain.

The last structure is a tree tensor network, shown in Fig. 2(d). We use “Tree” and “TreeX” to denote the tree structures shown in Fig. 2(b) and Fig. 2(d), respectively. In “TreeX”, the electronic states and vibration states for each OT molecule are first grouped to form a subtree. The subtree is exemplified by the OT2 subtree shown in Fig 2(c). The corresponding electronic node for OT2 contains three bases: the vacuum state, the LE2 state, and the CS2 state. Quantum number constraints are applied to restrict the wavefunction to the single-exciton manifold. For discussions on grouping versus distributing electronic states in tensor networks within the MPS context, we refer readers to existing literature baiardi2019large, sheng2024td. The overall tree is then constructed based on the 13 subtrees. The root node of the TreeX structure has 3 children. The first child contains OT1 and OT2 subtree, as well as the $R$ vibration. The second child contains the OT subtrees 3 to 6. The third child contains the rest of the subtrees, as well as the $F$ vibrations. TreeX is designed to minimize the entanglement entropy at the top layer of the tree.

For all tensor structures, the harmonic oscillator eigenbasis is employed for the primitive basis of the vibrations, unless otherwise stated. Following the reference work dorfner2024comparison, the number of the oscillator states for a given vibrational mode is set to $(g/\omega)^{2}+3g/\omega+N_{b}$, where $g$ is the maximum coupling constant across all types of couplings, and $N_{b}$ is an adjustable offset. In this work we employ $N_{b}=18$, which according to the reference paper should be enough for converged result dorfner2024comparison. In Appendix C, we show that using the discrete variable representation (DVR) basis does not significantly affect the calculated dynamics.

The multiset formulation employs multiple, independent tensor network states, which are then linearly combined to represent the total wavefunction fang1994multiconfiguration, kloss2019multiset. This approach has demonstrated good performance for the exciton dissociation model studied in this work lindoy2025pyttn. The multiset wavefunction allows different electronic states to couple with entirely different nuclear wavefunctions, which reduces the bond dimension for each individual tensor network state. In this work, we focus on singleset tensor network states where all electronic and vibrational states are encoded in one tensor network. However, we are actively developing a multiset implementation, and a systematic benchmark comparing multiset and singleset tensor network states will be the focus of our future research.

Another important aspect of MPS and TTNS algorithm is the corresponding operators, namely matrix product operators (MPO) and tree tensor network operators (TTNO) schollwock2011density, chan2016matrix, szalay2015tensor. Similar to MPS and TTNS, MPO and TTNO are tensor network representations of quantum operators.

Within ML-MCTDH, the Multi-Layer operators introduced by ML-Potfit are closedly related to TTNOs otto2014multi. ML-Potfit itself is a generalization of the Potfit algorithm jackle1996product, adapted for the multi-layer framework. Its primary function is to numerically construct the potential energy operator by fitting a continuous analytical potential energy surface into the required ML operator form.

In this work, MPO and TTNO are constructed exactly given the sum-of-product form of the operator. This distinguishes our approach from standard ML-MCTDH, where the direct sum-of-product form is generally used for model systems. In our previous work, we demonstrated that employing MPO provides a computational scaling advantage over the sum-of-product approach ren2022time. We utilize the automatic MPO/TTNO construction algorithm based on bipartite graph theory ren2020general, li2024optimal, which efficiently generates the most compact MPO/TTNO with negligible computational cost. Here, “most compact” means the exact MPO and TTNO with minimal bond dimension, which will directly affect the computational cost. The maximum MPO/TTNO bond dimension for Chain, Tree, ChainX and TreeX are 29, 29, 29, 14, respectively.

For comparison, the corresponding Hamiltonian operator contains 503 terms in the sum-of-product form.

For both Chain and Tree structures, the maximum MPO/TTNO bond dimension is associated with the electronic node, because all 26 electronic states are represented in a single leaf node.

The 29 operaters include 26 diagonal operators for each electronic state, one identity operator to be coupled with $\hat{H}^{\textrm{vib}}$, one electronic Hamiltonian operator $\hat{H}^{\textrm{e}}$ to be coupled with the identity operator of the vibrations, and one off-diagonal operator $\ket{\textrm{LE}_{1}}\bra{\textrm{CS}_{1}}+\textrm{H.C.}$ to be coupled with the $R$ mode. In TreeX, the maximum bond dimension is reduced to 14 because the diagonal electronic operators are coupled with the vibrations and merged with $\hat{H}^{\textrm{e}}$ in the subtrees. This is not possible in ChainX, where all diagonal electronic operators must individually couple with the $F$ mode placed at the boundary of the chain. For more details of MPO/TTNO construction, we refer readers to our prior works, where we provide the specific matrix elements of MPO and TTNO for several typical models ren2018time, li2024optimal.

The time evolution of MPS/TTNS is carried out using the TDVP-PS algorithm, which represents another key difference between our approach and standard ML-MCTDH. The TDVP-PS algorithm for TTNS (PSI-ML-MCTDH) and its variants have been known for over a decade lubich2015mctdh, lindoy2021time, lindoy2021time2, weike2021symmetries, however, the variable mean-field (VMF) method remains predominant for time integration in the ML-MCTDH community. The VMF method is also employed in the reference work. The PS algorithm features a sequential sweep over the tensor network to perform the time evolution, while VMF updates all nodes in the tensor network synchronously. We shall discuss the difference between the TDVP-PS and the VMF in Appendix C. In this work, unless otherwise stated, the 1-site TDVP-PS algorithm is employed. The product between a matrix exponential and a vector is solved by the Krylov solver, or more specifically the short iterative Lanczos algorithm for the Hermitian system studied here. For all results reported, the time evolution step is set to 1 fs. During initial benchmarks, we found that TDVP-PS supports a much longer time step than 1 fs, yet for the purposes of this paper, we choose 1 fs as the time step for denser data points. Benchmark data for the time step is provided in Appendix B.

3 Results and Discussion

We first compare the dynamics of $\braket{\hat{n}_{\textrm{LE}_{1}}}$ by Chain and Tree, calculated by Renormalizer, in Fig. 3. The employed bond dimensions $M$ for Chain and Tree are 512 and 256, respectively. The full dynamics in Fig. 3(a) agrees with the reference paper dorfner2024comparison, which describes the depopulation of the LE state to other LE and CS states. Moreover, the difference between the Chain and Tree results is almost negligible. In Fig. 3(b) we zoom in on region marked by the rectangle in Fig. 3(a) for a clear visualization of the difference at the long time limit. From Fig. 3(b), we estimate that the difference between the Chain and Tree results is at the order of 0.005. In Fig. 3(c), we illustrate the relative difference of $\braket{\hat{n}_{\textrm{LE}_{1}}}$, using the Tree results as the reference. In general, the relative difference increases over time, and the maximum difference is approximately 8%. We note that the same MPS ordering and TTNS structure as in the reference work are employed to produce Fig. 3. Overall, Fig. 3 shows that Chain and Tree produce consistent dynamics, with the differences being negligible for most practical purposes.

We then provide further evidence that in our calculation Chain and Tree are converging to the exact limit. In Fig. 4 we show how $\braket{\hat{n}_{\textrm{LE}_{1}}}$ evolves when $M\rightarrow\infty$. More specifically, we choose $t=60,90,160$ and 190 fs as several representative time frames, and plot $\braket{\hat{n}_{\textrm{LE}_{1}}}$ against $1/M$. For reference, we include the full dynamics calculated with different $M$ in Fig. B1 in Appendix B. In Fig. 4, we observe that as $M$ increases and $1/M$ decreases, the results from Chain and Tree converge. We then estimate the value of $\braket{\hat{n}_{\textrm{LE}_{1}}}$ in the $1/M\rightarrow 0$ limit by linear extrapolation. The technique is inspired by the extrapolation scheme employed in large scale static DMRG calculations white2005density, li2019electronic, larsson2025benchmarking. To exclude the noise induced by the calculations with small $M$, we employ only the two data points with the smallest $1/M$ for the extrapolation. Fig. 4 shows that this simple extrapolation scheme reduces the difference by Chain and Tree when $t=160$ and 190 fs. We note that such extrapolation is only valid when the bond dimension is sufficiently large. Fig. 4 shows that the dependence of $\braket{\hat{n}_{\textrm{LE}_{1}}}$ on $1/M$ is highly nonlinear. As a result, extrapolation using the datapoint with $M<100$ (for Tree) or even $M<200$ (for Chain) will likely result in even larger errors compared to the raw data. This behavior is fundamentally different from static DMRG calculations. There, the energy is directly minimized, resulting in a monotonic relationship with $1/M$ that enables reliable extrapolation.

Empirically, we found that the extrapolation improves the consistency of results across different tensor network structures particularly at the long time limit. Therefore, we consider the extrapolated results to be more reliable at the long time limit. However, in other computational tasks where such extensive benchmarks are unavailable, the extrapolation should be applied with caution.

Having confirmed that the dynamics in Fig. 3 is numerically exact, we investigate the origin of the discrepancy reported in the reference work. Following the prescription described in the reference paper, we reproduced these results using Renormalizer, which are shown in Fig. C1 in Appendix C. We find that the bond dimensions of both Chain and Tree are the key to the difference. On the one hand, the bond dimension of Chain is expanded using 2-site TDVP-PS with singular value truncation, and the bond dimension becomes fixed after the maximum bond dimension reaches a target value. While effective for short-time dynamics, this method creates a bond dimension distribution biased toward early-time entanglement. This becomes the source of the error for the studied model in which the entanglement spreads from the OT1 molecule to other molecules over time. In our approach for the results in Fig. 3, a fixed uniform bond dimension across the whole MPS chain is employed, avoiding the problem. On the other hand, in the reference work, the maximum bond dimension of the Tree calculation is $M=40$. However, our convergence analysis in Fig. 4 shows that further increasing $M$ can enhance the long-time accuracy. As a result, both the Chain and Tree results deviate from the exact solution by a small yet non-negligible margin. Further details are provided in Appendix C.

We note that despite the relatively small bond dimension employed in the reference work, the short-time dynamics agree nearly perfectly. Additionally, the maximum absolute difference in electronic state occupancy at the long-time limit is only 0.05, which is often sufficient for many practical applications.

The residual difference between Chain and Tree structures with large bond dimension motivated the development of optimized tensor network structures. For MPS and TTNS, the new structures are termed ChainX and TreeX, respectively, and are shown in Fig. 2 (c) and (d). ChainX and TreeX are designed based on the entanglement analysis of Tree calculations, shown in Fig. 5. The bonds with the largest entanglement entropy are shown as they are the bottleneck of the computation. Fig. 5 reveals that the low frequency OT modes exhibit significantly stronger entanglement compared to the $F/R$ vibrations, the electronic states, and the high frequency modes. High-frequency modes show negligible entanglement, allowing flexible placement in the network without affecting the calculated dynamics. Furthermore, the low frequency modes associated with the smallest OT index have the strongest entanglement. Guided by these observations, TreeX groups OT 1-2, OT 3-6 and OT 7-13 as the top subtrees, to distribute the entanglement evenly. ChainX is inspired by TreeX, where the OT 1 unit with the strongest entanglement is placed in the middle of the chain. We note that although our numerical experiments confirm that TreeX reduces errors compared to other tensor network structures, it is unlikely the optimal tree structure for this model and further structural optimization could lead to additional enhancements larsson2019computing, li2022fly.

Our ChainX/TreeX calculations, performed with bond dimensions up to $M=512$ and $M=256$ respectively and combined with the extrapolation scheme from Fig. 4, are shown in Fig. 6(a).

Chain, Tree, ChainX and TreeX results demonstrate exceptional agreement. The corresponding maximum absolute difference is estimated to be approximately 0.002, with an average error of approximately 0.001. Such precision is remarkable for the complex quantum system and highlights the reliability of numerically exact tensor network methods. The complete set of extrapolated results for all 26 electronic states is provided in Appendix D, and all source data are available in our repository data. We believe the data will serve as accurate benchmark data for the further development of quantum dynamics algorithms.

In Fig. 6(b) we show that TreeX shows the fastest convergence compared to Chain, Tree and ChainX structures. Using the extrapolated values from Fig. 6(a) as reference, we quantify the error as the time-averaged absolute deviation between $t=150$ and 200 fs:

[TABLE]

Here $n$ represents the occupation of the LE1 state. The reason for choosing the 150-200 fs window is two-fold. The first is that the error in this time window dominates the total error. The second is that throughout the manuscript we often focus on the dynamics in this time range for detailed analysis.

Our analysis reveals that TreeX achieves comparable accuracy at $M=32$ to what Chain attains at $M=128$. And when $M=128$ the average error of TreeX is already as small as 0.002, approaching the uncertainty limit of our reference values. This enhanced performance is consistent with the convergence trends shown in Fig. B1, and demonstrates that the optimized tree structure more effectively captures the entanglement pattern of the system. Meanwhile, ChainX exhibits the largest error compared with the other three tensor network structures. This is because, although the units with the strongest entanglement are placed in the middle, the ordering of the sites in ChainX does not follow the original site ordering of the model. This misalignment typically leads to large entanglement entropy and consequently a large error in MPS structures. In comparison, TreeX is able to be mapped onto a one dimensional topology while preserving site ordering and more critically to distribute entanglement entropy evenly across subtrees. ChainX cannot achieve this balanced distribution simply by placing the most entangled units in the middle. The contrast of the different convergence rates highlights the importance of careful structural design for achieving optimal computational efficiency.

In order to gain deeper insight into the different convergence rate of Chain, Tree, ChainX and TreeX, in Fig. 7 we analyze the time evolution of the von Neumann entanglement entropy $S$ for the four methods. The maximum $S$ across all virtual bonds in the correponding tensor network structure is reported. Since in this work we have employed a uniform bond dimension distribution, the maximum $S$ determines the error of the calculation. As shown in Eq. 15, the tensor network compression error at a particular bond is the sum of the square of the discarded singular values. Thus, the bond with the largest $S$ will have the largest absolute compression error when a fixed bond dimension is applied uniformly. In other words, this single largest error becomes the dominant and limiting source of inaccuracy for the entire calculation and determines its overall error.

The spatial distribution of maximum $S$ in the tensor network varies between methods. For Chain, the maximum $S$ appears between the local vibrations of OT1 and OT3 molecules, which is shown in Fig. C2. For Tree, the maximum $S$ is associated with the bond that connects the local vibrations modes of OT1 to OT7 to the rest of the system, which is show in Fig. 5. For ChainX, the maximum $S$ appears between the $R$ mode and the OT1 unit, which aligns with the design principle. For TreeX, the maximum $S$ is designed to appear in the top layer.

In Fig. 7, we first find that the convergence of maximum $S$ is much slower than the convergence of the occupation for all methods. Thus, the maximum $S$ may serve as a strict criteria for numerical convergence. Additionally, while Tree initially shows slower entanglement growth than Chain, its entanglement increases rapidly thereafter, ultimately reaching comparable values by $t=200$ fs. ChainX shows the highest entanglement entropy, which leads to the larger error observed in Fig. 6(b). In contrast, TreeX shows the slowest entanglement growth. As a result, when $t=200$ fs, TreeX has the smallest $\max_{\rm{bonds}}S$ across the four methods. This suppressed entanglement accumulation directly correlates with TreeX’s superior accuracy shown in Fig. 6(b).

Lastly, we discuss the efficiency of the algorithms. For the same bond dimension $M$, the MPS structures require less memory, which allows us to perform $M=512$ calculation with Chain/ChainX. Both Tree and TreeX have similar memory requirement. However, the TTNO of Tree has a larger bond dimension, leading to larger intermediate tensors and a higher computational cost compared to TreeX. MPS calculation with $M=256$ and $M=512$ and TreeX calculations with $M=256$ are carried out on a V100 (32GB) GPU card in combination with 4 cores of an AMD EPYC 74F3 CPU li2020numerical. Tree calculation with $M=256$ is carried out on a A100 (80GB) GPU card in combination with Intel Xeon Gold 6226R CPU, due to its higher memory consumption.

The thermal design power (TDP) for the AMD EPYC 74F3 CPU and the NVIDIA GPUs is rated at around 300 W. The TDP value represents their maximum thermal output under full load. Our calculations account for a more typical usage scenario where only one-sixth of the CPU cores are active, and GPU utilization averages 50%. Based on this reduced load, we estimate the total power consumption for this workload to be approximately 200 W.

Fig. 8 presents the tradeoff between the average error and computational efficiency for Chain, Tree, ChainX and TreeX. The wall times reported here include the time evolution and the calculation of all physical observables, such as the bond singular values and the RDMs of all degrees of freedom in the system. As shown, the TreeX structure demonstrates the best balance between accuracy and computational cost, achieving the smallest average error while requiring less wall time. Chain is more efficient than Tree, because of its lower computational scaling. This figure highlights that tree structure design is crucial to the success of the tree tensor network algorithms.

4 Conclusion and Outlook

In this work, we develop from a previous work and perform systematic benchmark of MPS and TTNS for a challenging exciton dissociation model. Using the TDVP-PS algorithm, where MPS time evolution corresponds to TD-DMRG and TTNS to PSI-ML-MCTDH, we demonstrate exceptional agreement between MPS and TTNS, with an average difference in electronic state occupancy of only 0.001. This is achieved through: (1) the use of sufficiently large bond dimensions ($M=512$ for MPS and $M=256$ for TTNS), which alone reduces the difference to $\sim 0.005$; (2) the development of an optimized tree structure termed TreeX that shows superior convergence properties. The exceptional efficiency of TreeX is demonstrated by its achievement of an average occupancy error $<0.005$ with only $M=32$ bases. We note that TreeX is not designed a priori, but rather based on the calculated entanglement entropy from TTNS simulations. While we also designed a similar new MPS structure termed ChainX, the error was higher in comparison. Our results show that the structural flexibility of the tree, compared to a linear topology, offers more flexibility to reduce entanglement entropy and enhance computational efficiency.

Our calculation confirms that both MPS and TTNS are numerically-exact ansätze that are capable of achieving high accuracy for the simulation of complex chemical systems with moderate cost. We hope the high-accuracy benchmark data presented here will facilitate future developments of not only tensor network algorithms but also other types of quantum dynamics methods. At the same time, our study highlights practical challenges in tensor network simulations. High-quality results require extensive convergence tests for parameters such as bond dimension and tree structures. To overcome this, we need new descriptors to quantify errors in time-evolved observables, as well as algorithms that can efficiently and reliably generate near-optimal tensor network structures without extensive trial calculations. We expect continued cross-fertilization between TD-DMRG and ML-MCTDH will enable the simulation of even more challenging quantum systems with unprecedented accuracy and efficiency.

Appendix Appendix A Model Parameters

In this section, we list the specific model parameters for the exciton dissociation model studied in this work. All parameters are directly adopted from the reference work dorfner2024comparison. In Table A1 we list the on-site energy of the 13 CS states, used in Eq. 3. In Table A2 we list the frequencies and coupling constants of the $F$ modes and the OT local vibrational modes.

Appendix Appendix B Convergence Benchmark

In this section, we present a convergence analysis of both the bond dimension and the time step. The dynamics of $\braket{\hat{n}_{\textrm{LE}_{1}}}$ with different bond dimension $M$ is shown in Fig. B1. While all four tensor network structures appear converged at the full scale in the first column, closer inspection of specific time windows reveals differences. The $50\sim 100$ fs region in the second column shows satisfactory convergence for large $M$ values across all methods. However, in the long time limit ($t=150\sim 200$, third column), Chain, Tree and ChainX exhibit significant convergence challenges, while TreeX maintains stable convergence behavior. From Fig. B1(c) and Fig. B1(f), we see that if $M$ is not sufficiently large, Chain consistenly overestimates the occupation, while Tree tends to underestimate the occupation. These opposing trends explain the discrepancies observed between the two methods when small bond dimensions are used. Nonetheless, Chain and Tree are converging to the same result as $M$ increases. The convergence patterns support our extrapolation approach shown in Fig. 4.

In Fig. B2 we show the effect of the time step on the dynamics obtained by Tree with $M=128$. We find that using a time step of 4 fs will leads to an error smaller than the estimated uncertainty caused by insufficient bond dimension. Thus for the benchmark data reported in this work in principle we can use 4 fs as the time step. However, as shown in Fig. B2(b), using such large time step will cause sparse data points and discontinuous dynamics. Thus in the main text we choose to employ 1 fs as the time step for denser data points.

We note that the large time step does not cause large error despite the high frequency modes in the model. This is because these high frequency modes shows neglegible entanglement entropy with the rest of the system, as shown in Fig. 5. Therefore, one may take the mean-field picture to deal with these high frequency modes. In other words, their fast motion constitudes a mean-field to the dynamics of the electron. For the purpose of calculating the electron occupations, it’s not necessary to exactly track the high-frequency nuclear motion. On the hand hand, the projector-splitting integrator empoloyed in this work features sequential update of each node in the tensor network. The integrator is then able to adaptively perform the time evolution of each node with different number of steps (the number of Krylov vectors in matrix exponential). The symplectic feature of the integrator also helps to ensure the numerical accuracy at long time scale. In our previous work, we have found that TDVP-PS tolerant exceptionally long time step in MPS simulations li2020numerical.

Appendix Appendix C Reproducing the Reference

Using our Renormalizer implementation, we successfully reproduce the difference between Chain and Tree structures reported in the reference work dorfner2024comparison. Fig. C1 demonstrates that our reproduced results and the original data points extracted from the reference work are in excellent agreement. This validation confirms both the reliability of our implementation and the reproducibility of the earlier findings. Besides, from Fig. C1(b), the Tree simulation in the reference paper agrees better with the results in this work compared with the Chain simulation.

The key factor enabling this reproduction is the bond dimension setup. While maintaining nearly identical computational parameters to those used in our main text, where we achieved 10% relative difference and 0.005 absolute difference, we specifically matched the bond dimension of the reference work. More specifically, for MPS calculation, starting from a Hartree product state where $M=1$, we perform 2-site TDVP-PS for a few steps, with a very small singular value truncation threshold of $10^{-6.5}$. While this procedure allows rapid growth of the bond dimension, it creates an uneven distribution where regions of initially strong entanglement develop large bond dimensions while other areas remain constrained. When the maximum bond dimension is above a target value, which in this case is 250, the time evolution algorithm is switched to 1-site TDVP-PS, and the time evolution proceeds with fixed bond dimension. The transition happens at very early stages of the time evolution ($t<20$ fs). A comparison of the bond dimension obtained using this approach and the fixed uniform bond dimension employed in this work is included in Fig. C2(a). The direct result of this setup is that the region where the initial entanglement is small has constrained bond dimension throughout the time evolution. In Fig. C2(b), we plot the accurate bipartite entanglement entropy at each bond using our uniform bond dimension. We can see that the region where initially has small entanglement developes moderate amount of entanglement entropy at later times. The small bond dimension in this region creates persistent bottlenecks in the calculation and limits the overall accuracy regardless of maximum bond dimension elsewhere. Increasing the maximum bond dimension has very slow convergence because the maximum bond dimension is increasing much faster over time than the bond dimension in the low-entanglement region.

On the other hand, for the Tree structure calculation, we replicate the tree structure and bond dimension reported in the reference work, including its maximum bond dimension of 40. As shown in Fig. 4, $M=40$ will underestimate the occupation in the long times, which is consistent with the underestimation in Fig. C1(a). The combined inaccuracies from both MPS and TTNS implementations account for the observed discrepancy in the reference work.

To investigate the potential effect of methodological differences beyond bond dimension, we examine two additional factors from the reference study. The first is the use of different time evolution schemes. The reference work employs the PS method for MPS and the VMF beck2000multiconfiguration scheme for TTNS, as VMF is the traditional integrator for ML-MCTDH. To investigate their effect, we perform further benchmarks based on the Tree structure and bond dimension in Fig. C3. We find that the VMF and PS schemes show excellent agreement, with $\sim 2\%$ relative difference (0.001 absolute difference), which is consistent with our previous MPS benchmarks  li2020numerical. Besides, in the reference work, the harmonic oscillator eigenbasis and DVR basis are employed in MPS and TTNS calculation respectively. We perform VMF time evolution with sine DVR basis beck2000multiconfiguration and the results are also included in Fig. C3. We find that the two factors collectively contribute to 3% of relative difference (0.002 of absolute difference). In Fig. C4 we show the effect of different regularization parameter $\varepsilon$. In our VMF integration the density matrix $\rho$ is regularized using the following equation

[TABLE]

where by default $\varepsilon=10^{-10}$. As shown in in Fig. C4, using a small $\varepsilon$ will lead to larger error in the initial dynamics. In the long time limit, the difference caused by $\varepsilon$ is smaller than 6% (0.003 of absolute difference). We conclude that neither factor significantly contributes to the discrepancy observed in the reference work, and that bond dimension remains the primary determinant of accuracy in these calculations.

Appendix Appendix D Extrapolated Dynamics

In this section, we show the extrapolated dynamics of all 26 electronic states in Fig. D1. The results by Chain, Tree, ChainX and TreeX are all shown in each panel. The excellent agreement across all states validates the consistency of MPS and TTNS, with most cases showing nearly indistinguishable curves between methods. While minor deviations appear for certain states, such as in Fig. D1(m), these occur primarily where the occupation is small in magnitude. For example, in Fig. D1(m), LE7 shows visible difference between Tree and TreeX, yet the absolute difference is approximately 0.004.

Acknowledgement

Weitang Li and Jun Yan are supported by the Guangdong Basic Research Center of Excellence for Aggregate Science and the Shenzhen Science and Technology Program (No. KQTD20240729102028011). Jiajun Ren is supported by the National Natural Science Foundation of China (Grant No. 22273005 and No. 22422301).

Code Availability

The complete code that supports reproducing all data in this study is deposited on Zenodo data.

Data Availability

The source data for all figures and the Python scripts to produce all figures are deposited in Zenodo data.

Competing interests

The authors declare no competing interests.