A hierarchical wavepacket propagation framework via ML-MCTDH for molecular reaction dynamics

Abstract

This work presents a computational framework for studying reaction dynamics via wavepacket propagation, employing the multiconfiguration time-dependent Hartree (MCTDH) method and its multilayer extension (ML-MCTDH) as the core methodologies. The core idea centers on the concept of modes that combine several coordinates along with their hierarchical separations because the degrees of freedom are too numerous to be efficiently treated as a single mode. First, the system is partitioned into several fragments within the same layer, and these fragments are further decomposed. Repeating this process, a hierarchical separation of modes emerges, until modes of a manageable size are achieved. Accordingly, the coordinate frame can be designed hierarchically. Second, the kinetic energy operator (KEO) is derived as a sum-of-products (SOP) of single-particle differential operators through a polyspherical approach, while the potential energy surface (PES) is expressed in a similar SOP form of single-particle potentials (SPPs) through (1) reconstruction approaches using an existing PES or (2) direct approaches based on a computed database. Third, the nuclear wave function is expressed in a multi-layer expansion form, where each term is a product of single-particle functions (SPFs) that are further expanded by the SPFs in the deeper layer. This expansion form is also adopted using a variational eigensolver for electronic wave function. Finally, the Dirac–Frenkel variational principle leads to a set of working equations whose solutions reproduce reaction dynamics, say reaction probability and time-dependent expectation. In addition, the hierarchical framework can be rearranged by the mathematical language of tensor networks (TN) or tree tensor networks (TTN). In this work, we compare the methods represented by a function with those in the form of a TN or a TTN. We also discuss the limitations of the present framework and propose solutions, providing further perspectives.

---

Xingyu Zhang is currently pursuing her PhD degree in chemistry at Northwestern Polytechnical University after receiving her BS degree in chemistry at the same university in 2024 (both supervised by Dr Q. Meng). Her research interests focus on chemical dynamics. In the preparation of this manuscript, she was primarily responsible for drafting the initial text, collecting the literature, and organizing figures and tables.

Qingyong Meng received his PhD degree in chemistry from the University of Chinese Academy of Sciences in 2011 (supervised by Prof. Dr M.-B. Huang). After spending years at Ruprecht-Karls-Universität Heidelberg (with Prof. Dr H.-D. Meyer) and Dalian Institute of Chemical Physics of the Chinese Academy of Sciences (with Prof. Dr D. H. Zhang) as a postdoctoral fellow, he became an associate professor (in 2016) at Northwestern Polytechnical University. His research interests focus on chemical dynamics. In the preparation of this manuscript, he contributed through scientific guidance, re-organizing the manuscript and textual refinement.

---

I. Introduction

By solving the Schrödinger equation (SE), i∂Ψ(t)/∂t = ĤΨ(t), we can directly demonstrate chemical dynamics from a microscopic viewpoint. Here, Ψ(t) represents the time-dependent (TD) wave function of a molecular system, and Ĥ = T̂ + V̂ is the Hamiltonian operator, where T̂ and V̂ are the kinetic energy operator (KEO) and the potential energy surface (PES), respectively. For a chemical reaction, we can either examine the transformation between the eigenstates of the reactant and product, known as time-independent (TI) quantum scattering, or investigate time-dependent (TD) wave function of the entire system, known as wavepacket propagation. Usually, the Hamiltonian operator Ĥ of a molecular reaction system is time-independent and hence allows the separation of time t from spatial coordinates making quantum scattering and wavepacket propagation equivalent up to a time-energy Fourier transform. The computational methods of wavepacket propagation, as shown in Fig. 1 have matured into a sophisticated and tightly coupled framework that enables applications to chemical dynamics. However, even elementary reactions exhibit intricate mechanisms, posing significant theoretical and computational challenges. This complexity necessitates the integration and coordination of the quantum dynamics methods to achieve good performance in calculations. A common way to achieve this goal is to express both the wave function and the Hamiltonian operator in a hierarchical expansion form leading to a tree-structure. The tree-structure defined by the connectivity of the isometry single-particle objects plays a crucial role in improving its numerical performance. This is because the tree-structure can effectively enhance flexibility of the wave function and the Hamiltonian.

Fig. 1 Systematic illustration of the quantum dynamics methods. The black words represent computational goals or results, while the red words represent computational techniques. The arrows represent the direction from the input to the computational results. Designing a set of coordinates, the kenetic energy operator (KEO) can be derived through the polyspherical approach. Then, the database is constructed through ab initio energy calculations on a set of sampled geometries. The potential energy surface (PES) is then constructed in an appropriate form. With the Hamiltonian operator, extensive wavepacket propagations can be launched on the basis of the rovibrational eigenstate of the reactant which can be computed by relaxation of a guess function. Then, the initial wave function is computed by assigning a given value of momentum to the initial eigenstate. Having time-dependent wave function, flux analysis and expectation calculations are launched to compute reaction probability and time-dependent expectations. The present work presents a computational framework for wavepacket propagation employing the multiconfiguration time-dependent Hartree (MCTDH) and its multilayer extension (ML-MCTDH) methods, where both the Hamiltonian and wave function are expanded in hierarchical form. The blue arrows belong to construction of the Hamiltonian as given in Sections II and III. The black arrows belong to the wavepacket propagation and analysis of the wave function as give in Section IV.

In this work, we focus on techniques for wavepacket propagation in molecular reaction dynamics, as shown in Fig. 1, including (1) construction of Ĥ (blue arrows), (2) propagation of an initial eigenstate Ψ(t = 0) (black arrows), and (3) analysis of the propagated wave function Ψ(t) (yellow arrows). A hierarchical framework is outlined in this work to review recent developments mainly proposed by this laboratory,^1–4 especially the construction of the Hamiltonian operator Ĥ. The core idea centers on the concept of modes that combine several coordinates along with their hierarchical separations, because the degrees of freedom (DOFs) of high-dimensional systems are too numerous to be efficiently treated as a single mode. Repeatedly breaking the mode, hierarchical separations can be designed resulting in a multilayer representation of the molecular system until modes of a manageable size are achieved. In the simplest case, a d-dimensional system with mass-weighted coordinate set has Ĥ and Ψ(t) in the summation of orthogonal products (also called configurations) of the basis, namely

where and are single-particle operators and functions, respectively, while I = (i₁, …, i_d) means the grid used to represent Ĥ and Ψ(t). In eqn (1) and (2), and are thus the configuration basis defined on grid I. If q^(κ) is the κ-th mode instead of the coordinate, and can be further expanded similarly leading to the multilayer sum-of-products (SOP) form.^1–29 The TD functions in the deepest layer are finally expanded either by a given TI basis set in finite basis representation (FBR) and discrete variable representation (DVR)³⁰ or by a TD basis set of correlation DVR (CDVR)^20,31–34 based on multidimensional time-dependent DVR (TDDVR). For CDVR, eqn (1) and (2) are not strictly required. In this work, we mainly focus methods based on the FBR or DVR that require eqn (1) and (2).

The present framework for wavepacket propagation originates from representation of a molecular fragment by a z-matrix which is well-established. To integrate multiple z-matrices corresponding to individual fragments into a unified representation for the whole reaction system, a z-system was developed.^35,36 For a finite set of points in Euclidean space, this mathematical formalism³⁵ generalizes to a comprehensive theory of polyspherical coordinates. Such coordinates are defined on the orbit spaces of the group of rigid motions.³⁵ This representation enables the construction of hierarchical coordinate frames and facilitates the derivation of KEO expressed as a summation of derivative product terms through a polyspherical approach.^36–41 Following eqn (1), significant developments have been achieved to construct a PES in a sum-of-products (SOP) form either by an existing PES in the general form or by a set of energy values at sampled geometries. Employing eqn (1) and (2), various propagation approaches have been developed over the past few decades. Depending on the number of layers in the expansional ansatz, these approaches are convention (or known as standard) methods, time-dependent Hartree (TDH) products, multiconfiguration TDH (called MCTDH) and multilayer MCTDH (called ML-MCTDH). Due to flexibility of the wave function ansatz, both MCTDH and ML-MCTDH methods have greatly influenced the field of quantum molecular dynamics and continue to do so. In this work, we will mainly discuss the hierarchical wavepacket propagation for molecular reaction dynamics using the MCTDH language which employs the Schrödinger picture with a grid-based coordinate representation. This is a natural choice for chemists who want to directly demonstrate how a chemical reaction occurs. Moreover, grid-based representations of the wave function and PES yield tensor-structured expressions, analogous to a tensor network (TN) or a tree tensor network (TTN) in quantum many-body systems,⁴² such as the matrix product state (MPS) and matrix product operator (MPO). By such tensor-structured representation, the time-dependent variational principle indicates the time-dependent density matrix renormalization group (TD-DMRG) method^43–54 or a variant of MCTDH, called MPS-MCTDH.²⁸ This is not surprising because chemical dynamics is fundamentally a many-body problem and the current ML-MCTDH-based framework employs the single-particle approximation (SPA).

Furthermore, in this work we compare the wavepacket propagation methods with the electronic-structure methods because the antisymmetric wave function in the SOP form is widely adopted by variational approaches for electronic structures.^55,56 These methods are Hartree–Fock (HF), full configuration interaction (CI), multiconfigurational self-consistent field (MCSCF), and so forth. The one-dimensional basis function in the deepest layer is called the molecular orbital (MO) which is further expanded by a given set of basis functions. Since these electronic-structure methods also employ the SPA, substituting the antisymmetry MPS as an ansatz into the variational principle, one can derive working equations for the quantum chemistry (QC) time-independent density matrix renormalization group (TI-DMRG or called QC-DMRG).^56–60 The TI-DMRG method is a special case of time-independent TN (TI-TN) optimizing a one-dimensional (1D) MPS but not a time-independent TTN (TI-TTN). In the following, we collectively refer to TD-DMRG and TI-DMRG as DMRG. Initially, ML-MCTDH and DMRG focused on very different applications. The former is based on quantum molecular dynamics,^19–21 whereas the latter is based on the quantum many-body dynamics.^57,58 Approximately two decades after their inceptions, both ML-MCTDH and DMRG were extended to incorporate the capability of solving each other's problems. The different individual developments of ML-MCTDH and DMRG led to very different mathematical expressions. Recently, Larsson⁵ mediated the communicative barriers between the ML-MCTDH and TD-DMRG formulations by deriving working equations through counterpart's language. Considering the present goal, we have no choice but to limit the present work to omit the details of TI-DMRG and TD-DMRG. We would like to finally emphasize that all of these approaches were developed to solve many-body equations of motion based on the SPA. For solving these equations, the most crucial aspect is handling the correlations among individual particles. These correlations are largely embedded in the interaction term of the Hamiltonian and the symmetry nature of the system.

In this work, to demonstrate the present hierarchical framework for wavepacket propagation, vibrational eigenstates and sticking probability of CO/Cu(100),⁴ dissociative chemisorption probability of H₂O/Cu(111),⁶ and reaction probability of OH + HO₂ → O₂ + H₂O (denoted by HO_x)¹ are taken as examples. Vibrational eigenstates and the probabilities of chemical reactions or physical processes are central problems in chemical dynamics. We refer the reader to the suppplementary information (SI) for implementations and numerical details of these examples. We should mention that the present approaches are not yet universally applicable to a general reaction and require targeted modifications and adaptations. The paper is organized as follows; in Sections II and III, we describe construction approaches of KEO and PES, respectively. Section IV describes function expansions for wavepacket propagation and compares these with those for electronic structures. Section V presents extensions of the present theory. Section VI concludes with a summary.

II. Coordinate set and the kinetic energy operator

To easily derive the KEO, Cartesian coordinate set is the simplest choice to represent the molecular system because all of the coordinates are orthogonal. For a N-atomic system, multiplying the square root of associated mass {m_α}^N_α=1 mass-weighted Cartesian coordinate set [Fraktur X]_N = {x_α,y_α,z_α}^N_α=1 indicates the KEOwhich is automatically given in the SOP form. Unless otherwise specified, in this section mass-weighted coordinates are always used. However, the Cartesian coordinate set is rarely employed in quantum dynamics because chemical processes with breaking and forming of the chemical bonds cannot be intuitively described by such coordinate set. Moreover, Cartesian coordinates always introduce artificial correlations slowing down computational convergence. To overcome this question, generalized coordinates are designed as linear or non-linear functions of the Cartesian coordinates. A set of linear functions of is the rectilinear coordinate set. A typical example is the vibrational normal mode of a bound system. A set of non-linear functions of is curvilinear coordinate set which is able to describe rotations and motions with large amplitude. Both rectilinear and curvilinear coordinates might be involved in describing chemistry dynamics. Therefore, the former is well-suited for quantum dynamics calculations of spectra, while the latter is appropriate for molecular reaction dynamics.

Over the past several decades, the polyspherical approach has been developed to design an appropriate set of curvilinear coordinates and to derive the KEO.³⁶ Its core idea involves partitioning the system into multiple subsystems of sufficient simplicity to facilitate coordinate definition, while employing Euler angles to describe the relative orientations between subsystems and the system itself. Moreover, such separation allows the reaction coordinates to be naturally defined as a subset of large-amplitude curvilinear coordinates. Fig. 2 illustrates the hierarchical scheme for designing coordinates through recursive fragmentation of the system into fragments (i.e., subsystems), corresponding to the decomposition of modes into deeper layers. Compared with the previously reported polyspherical approach,^36,39,40 the present hierarchical scheme is given in a multilayer fashion instead of one-layer. As shown in Fig. 2, the translational and rotational motions of the entire system are described in the space-fixed (SF) frame or called the laboratory-fixed (LF) frame, whose axes are fixed at the laboratory. Similarly fixing axes at the system itself, we obtain the E₂ frame in which the internal motions are described. Second, to describe internal motions of individual fragments, body-fixed (BF) frames are established by fixing their axes to specific fragments. The relative orientations between each BF frame and its parent frame in the last layer are specified through three Euler angles {α^(BF),β^(BF),γ^(BF)} that form a relative frame. Repeatedly introducing BF frames with respect to those in the last layer, a hierarchical frame is designed. The deepest layer BF frames describe the smallest inactive fragments during the reaction using minimal coordinates, typically represented by Jacobi or Radau vectors.⁴¹Fig. 3 illustrates the applications of the hierarchical coordinate frame to the present benchmarks, CO/Cu(100),⁴ H₂O/Cu(111),⁶ and HO_x.¹ These coordinate frames are all designed in one layer.

Fig. 2 Hierarchical scheme to design coordinates, where a circle indicates a frame and a box represents the coordinate in the deepest layer. The first layer separates the E₂ frame from relative coordinates with respect to the SF frame, denoted by E₂/SF. Dividing the system into several subsystems, a BF frame is designed for each subsystem. The second layer of the E₂ frame separates all of the BF frames from each other. Moreover, the relative coordinates of BF frame with respect to the E₂ frame (denoted by BF/E₂) are separated from other BF frames in the second layer. Repeating these separation processes, we can obtain all of the coordinate frames. The deepest layer gives (1) either Jacobi or Radau coordinates to describe motions of the fragment and (2) Euler angles or relative positions to describe relative motions.

Fig. 3 Hierarchical coordinate frames of (a) CO/Cu(100), (b) H₂O/Cu(111), and (c) HO_x, which are designed according to Fig. 2. For CO/Cu(100) and H₂O/Cu(111), surface scattering processes are described by the E₂ frame which is separated into (1) BF_CO or BF_wat for molecular motions, (2) relative frame BF_CO/E₂ or BF_wat/E₂ for relative positions, and (3) BF_Cu for surface atoms which is divided into BF₁ and {BF_i} for top and non-top atoms, respectively. The relative frame E₂/SF is further required to describe molecular motions with respect to the surface. For HO_x, if the total angular momentum of the whole system is set to be zero, the Euler angles in the relative frame E₂/SF are all equal to zero which implies that the E₂ and SF frames are identical. The OH and HO₂ fragments are described by BF_OH and BF_HO2, respectively, together with relative frames BF_OH/E₂ and BF_HO2/E₂. The relative frame BF_OH/BF_HO2 is required to define the principal component of the reaction coordinate.

Having mass-weighted coordinates of a d-dimensional system, the KEO can be derived through the metric tensor g^κρ and the Christoffel symbol of the second kind Γ^υ_κρ,^36–40

Both g^κρ and Γ^υ_κρ are determined by the coordinate transformations between {q^(κ)}^d_κ=1 and standard Cartesian coordinates. Several numerical procedures, say the TNUM (means numerical T̂) procedure,^37,38 were developed to obtain g^κρ and Γ^υ_κρ and their derivatives and to finally compute numerical KEOs. Moreover, employing the polyspherical approach, an automatic procedure^39,40 was developed to provide explicit analytical expression of the KEO which can be adopted by MCTDH and ML-MCTDH. In a word, the starting point of these numerical procedures is the relation between the mass-weighted Cartesian coordinates in the BF frame and the curvilinear coordinates. For the (ij⋯pk)-th BF frame (denoted as BF_ij⋯pk in Fig. 2) associated with an N_ij⋯pk-atomic fragment in the deepest layer, the corresponding KEO T̂^(BF_ij⋯pk) can be expressed in terms of 3N_ij⋯pk − 6 momentum operators and three angular momenta {Ĵ_ι}_ι=1³. This operator decomposes into vibrational (i.e.), rotational (i.e.), and Coriolis coupling (i.e.) components,^36–40The parameter tensors include (1) the (3N_ij⋯pk − 6) × (3N_ij⋯pk − 6) matrix Σ, (2) the (3N_ij⋯pk − 6) × 3 matrix σ, and (3) the 3 × 3 matrix Γ. It is often impossible to completely separate these three terms and the splitting itself is artificial.³⁶ This is because the Coriolis coupling and the rotational part depend on the BF_ij⋯pk frame. Following the derivation of T̂^(BF_ij⋯pk), the KEO in the BF_ij⋯p frame can be written as a summation of T̂^(BF_ij⋯pk) and the KEOs from relative frames.³⁶Eqn (5) in a summation of products of coordinate-associated operators naturally leads to a multi-layer expansion form for the final KEO expression. We refer the reader to Section S-I.B of the SI for the KEOs of the present examples.

Main advantages of the hierarchical coordinate frame are (1) its clear descriptions of the chemical and physical insights in a reaction, and (2) its capability to automatically^36–40 derive the KEO in the summation of products of mono-mode operators through the polyspherical approach. First, the reaction dynamics involves large-amplitude motions which are described by curvilinear coordinates. By the hierarchical frame, one can simply define the curvilinear coordinates for large-amplitude motions and hence define the reaction through chemical inspirations. It is generally accepted that a given reaction coordinate maps to a specific reaction channel or mechanism. Consequently, carefully designed reaction coordinates based on chemical inspirations enable both (1) more efficient sampling of geometries in configurational space for PES construction (see Section III) and (2) deeper analysis of the corresponding reaction dynamics. This represents a crucial feature when dealing with large elementary-reactions with complex dynamics mechanisms. For instance, chemical inspirations indicate that the HO_x system has four channels and thus three reactions starting from OH + HO₂,

OH + HO₂ → O₂ + H₂O, OH + HO₂ → H₂ + O₃, OH + HO₂ → O + H₂O₂. (6)

Obviously, coordinates defined in Fig. 3(c) are only suitable to describe dynamics of the OH + HO₂ → O₂ + H₂O reaction. In addition, a description with individual fragments allows us to use direct products of a one-dimensional basis avoiding the singularity problem. This also leads to a reduced coupling between the coordinates. Second, for a one-layer separation scheme, Lauvergnat and co-workers^37,38 as well as Ndong and co-workers^39,40 have developed procedures to derive the KEO in numerical or analytical form. For a hierarchical frame, it is also possible to develop a similar technique. Currently, however, a single-layer scheme is sufficient to describe typical elementary reactions. The hierarchical frame incorporates an interface for future consideration of substituent motions based on the one-layer scheme.

One notable disadvantage of the hierarchical coordinate frame is the redundancy in Euler angles. Specifically, when two relative frames share a common axis (typically denoted by the z-axis), the two Euler angles around this axis depend on each other forming a pair of redundancy azimuth angles, called {α_i}_i=1². To address this issue, one of the authors⁶¹ proposed a solution by applying a Fourier transform to {α_i}_i=1², converting them into their discrete momentum counterparts. Consequently, the corresponding KEO terms {∂²/∂α_i²}_i=1² are replaced by the square of the angle quantum numbers {k_i²}_i=1². Moreover, the action of V̂ on Ψ is accordingly written as an expansion form in the angular-momentum representation,⁶¹

where and satisfies the Fourier transformAnother disadvantage arises from the non-inertial nature of the hierarchical coordinate frame. This issue is not new. The ancient geocentric theory describes planetary motions using a non-inertial frame with the epicycle-on-deferent technique. The epicycle was defined to move in longitude, while the planet moves on the epicycle at a specific velocity.⁶² Similarly, the present hierarchical frame represents arbitrary periodic motions through the superposition of a set of circular motions, which introduces numerical complexity and reduces the intuitiveness of the resulting KEO. A potential solution to mitigate this problem is to adopt an occupancy representation rather than a coordinate representation implemented in a first-quantization (FQ) framework. By the second-quantization (SQ) formalism (shown in Section V later), the atomic motions are represented by excitations of virtual particles (or fields). The kinetic energy term is fundamentally a representation of the single-particle KEO expressed in a chosen basis. In other words, it is constructed from matrix elements of the single-particle KEO in the selected basis. This is helpful in simplifying the KEO term.

III. Decomposition of the potential energy surface

Since the KEO has already been given in the SOP form, as given using eqn (5), now the remaining task to build the Hamiltonian operator Ĥ in the form of eqn (1) is to construct the PES in sum-of-products (SOP) or canonical-polyadic (CP) decomposition (CPD) of single-particle potentials (SPPs),where C_i1⋯id and C_r are the expansion coefficients or called the coefficient tensor, R the CPD rank, and v^(κ)_i(q^(κ)) the i-th SPP for the κ-th coordinate. If q^(κ) means the κ-th mode, v^(κ)_i(q^(κ)) can be further expanded similarly leading to multi-layer expansion form. For the simplest one-layer expansion form of eqn (9), two primary ways can be employed. One way is direct construction from an energy database, such as reproducing kernel Hilbert space interpolation (RKHSI)⁶³ and revised Gaussian process regression (GPR) proposed by this laboratory.^1,3 Another way is reconstruction using an existing PES, say the potential re-fitting (POTFIT) method.^64–66 We collect these methods in Table 1. On the other hand, if the MPS ansatz is applied to represent the wave function, the matrix product operator (MPO)^67,68 as a tensor network operator should be adopted to represent potential tensor on grids. The conventional MPO is designed to handle a SOP of one-body operators such as position operators and second quantization operators. An efficient algorithm for compressing a given PES into a grid-based MPO was recently developed⁶⁹ to work with MPS-MCTDH, where the MPO is constructed from the grid-based numerical integral of the sum of arbitrary many-body coupled operators. It has capability to provide numerical integrals of arbitrary potential operator representations. The potential tensors with the SOP or CPD form as given in eqn (9) represent two of the tensor decomposition approaches^70–81 for the approximation of functions. In Table 2, we briefly summarize tensor decomposition approaches for the PES construction.

Table 1 Comparison of various methods to build decomposed PES with the MCTDH language, including direct construction (upper panel) and reconstruction (lower panel) methods. Abbreviations of these methods are explained in the main text. The third column gives the number of expansional layers of the potential function. The fourth column gives techniques to derive these approaches. The fifth and sixth columns give the rank of the resulting PES and remarks on these methods, respectively. The rightmost column gives references that reported these methods

No.	Method	Layer	Technique	Rank	Remark	Ref.
a The SOP-NN function of eqn (10) is actually in the CPD form. b We refer the reader to the SI for details. c Values of n and nsv are numbers of training data and support vectors satisfying nsv < n. d It is also possible to extend CPD-ws-SVR to its mode-combination version. e It can be extended to the mode-combination or warm-started version. f Using FBR, the wave function is represented by a finite set of time-independent functions. g Parameters {mκ}^dκ=1 are the number of terms in the decomposed function. h Number of expansion terms should be given before interpolation by V^(CP-FBR).
Direct construction
1	SOP-NN^a	One	GLR/NN	F ^(1)	Rank is the number of neurons	82 and 83
2	ML-SOP-NN	Multi	GLR/NN		Rank is the number of neurons	This work^b
3	CPD-GPR	One	KMR/GPR	n	n cannot be optimized	3
4	CPD-mc-GPR	One	KMR/GPR	n	Mode combination is adopted	1
5	CPD-SVR	One	KMR/SVR	n sv	n sv can be optimized	85
6	CPD-ws-SVR^d	One	KMR/SVR	n sv	Warm-started SVR is used	85
7	ML-CPD-GPR^e	Multi	KMR	∼l·n^c	Multi-layer CPD-GPR	This work^b
8	ML-CPD-SVR^e	Multi	KMR	∼l·nsv^c	Multi-layer CPD-SVR	This work^b
9	RKHSI	One	Interpolation	n	Interpolation with kernel function	63
10	SOP-FBR^f	One	Interpolation		See V^(SOP-FBR) in eqn (13)	93
11	CP-FBR^f	One	Interpolation	R	See V^(CP-FBR) in eqn (14)	94
Re-fitting methods
12	POTFIT	One	Regression		Error as target for d ≤ 7 systems	64–66
13	MGPF	One	Interpolation		Interpolation from coarse to fine grids	86
14	MLPF	Multi	Regression		Suitable for d > 7	87
15	MCPF	One	Regression		Monte Carlo for d > 7 SOP	89
16	MCCPD	One	Regression	R	Monte Carlo for d > 7 CPD	4 and 90

---

Table 2 Comparison of various tensor-decomposition methods for approximating function by tensor language, where an existing function is required. For the PES construction, these methods belong to the reconstruction scheme. Abbreviations of these methods are explained in the main text. The third column gives implementation complexity. The fourth column gives the number of expansional layers of potential function. The fifth column gives techniques to derive these approaches. The sixth and seventh columns give the rank of the resulting PES and remarks on these methods, respectively. The rightmost column gives references that reported these methods

No.	Method	Complexity	Layer	Technique	Rank	Remark	Ref.
a Rank of the CP decomposition is hardly determined. b SVD means singular value decomposition. c Rank of the Tucker decomposition is not clearly defined. d The Tucker decomposition requires that dimensionality is smaller than six. e TT + AS is the first application of the TT cross approximation (CA) procedure for the PES construction. f Rank should be determined carefully. g In the active subspaces (AS), the potential function has the most variability. The affine transformation is often used in the construction of the normal modes. h CA means cross approximation, also known as skeleton decomposition.
1	CP	Medium	One	A sum of rank-1 tensors	Low^a	Moderate dimensionality	70 and 75
2	Tucker	High	One	SVD to higher dimensions^b	—^c	Low dimensionality^d	71 and 75
3	FT	High	One	Combines with Chebyshev basis expansions	—^c	Three dimension	72
4	HT	Very high	Multi	Factorization using recursive subspace splitting	—	Suitable for very high-dimensional function	73 and 74
5	TT	Low	One	A sequence of matrix products	Finite	Suitable for very high-dimensional function	76
6	TT + AS^e	Not high	One	An affine transformation of Cartesian coordinates into the active subspaces^g	Finite^f	Coordinates should be carefully determined to get a good low-rank tensor approximation of multidimensional function	78
7	FTT	Not high	One	Combines with basis expansions	Finite^f	Chebyshev basis is a choice	79
8	EFTT	Not high	One	Compressed FTT approximation	Low	Functional low-rank approximations	80
9	CA^h	Medium	One	Interpolation of selected tensor	Finite	Suitable for function with any dimensionality	81

---

First of all, the direct approaches build the potential function in the SOP or CPD form through a discrete set of energy data. Originating from generalized linear regression (GLR), the neural network (NN) approach has gained widespread adoption in building the PES owing to its universal approximation capability for continuous functions. It is this remarkable property that motivated Manzhos and Carrington⁸² and Koch and Zhang⁸³ to propose a single-layer NN architecture in the SOP/CPD form, subsequently termed SOP-NN,

where b⁽²⁾ and c_k⁽²⁾ are parameters, while b_kκ⁽¹⁾ and w_kκ⁽¹⁾ are biases and weights, respectively. The V^(SOP-NN) formulation in eqn (10) can be extended to a multi-layer architecture by sequentiall0y inserting the output of one neuron into the activation function of the subsequent layer. However, this nested structure inherently loses the CPD characteristic. This limitation can be resolved through a multi-layer SOP-NN approach, called ML-SOP-NN, where q^(κ) in eqn (10) represents the κ-th mode rather than a coordinate. The feasibility of this approach stems from the fact that the utilization of modes and the mode combination¹¹ preserve the essential characteristics required for PES construction. By iteratively incorporating neurons into activation functions following a hierarchical mode structure, one can derive the PES in a multi-layer expansion form, denoted by V^(ML-SOP-NN).

It has been proved that the kernel-model regression (KMR), such as GPR, is equivalent to single-layer NN with infinite neurons.⁸⁴ Inspired by SOP-NN and RKHSI,⁶³ this laboratory proposed a CPD-type GPR function, called CPD-GPR,³ and its revisions with mode-combination (mc) scheme,¹¹ called CPD-mc-GPR,¹

where K^(κ)(·,·) is the kernel function for the κ-th coordinate, ξ_r the r-th optimized coefficient, geometries of n training data. Since the mode combination scheme¹¹ does not affect the construction of PES, we no longer need to distinguish CPD-GPR and CPD-mc-GPR. Eqn (11) indicates that V^(CPD-GPR) has a total of n expansion terms (rank) equal to the number of training data.^1,3 This value is often larger than 10⁴ and leads to low-performance in subsequent calculations.^1,3 To reduce the rank, a method to select decisive data was developed by this laboratory through support vector regression (SVR) in the CPD form, called CPD-SVR.^84,85 The resulting PES is similar to eqn (11) but n is replaced by the number of support vectors n_sv < n,^84,85where intercept b is required by the selection algorithm.^84,85 To easily select support vectors, an existing PES is introduced as the initial decision surface to warm-start (ws) the selection process. This is the CPD-ws-SVR method.⁸⁵ Illustrated in Fig. 4 are comparisons of the numbers of training data n and support vectors n_sv by ratios n_sv/n for various potential energy intervals of H + H₂ and ratios n_sv/n as a function of n for H₂ + H₂ and H₂/Cu(111). These numerical calculations⁸⁵ well demonstrated that SVR reduces the ratio to as low as around 30%–80%, indicating that SVR can effectively lower the rank. We refer the reader to ref. 85 for numerical details of these CPD-ws-SVR calculations. Concerning on the resulting PES expression, we will no longer distinguish CPD-SVR and CPD-ws-SVR, as well as CPD-mc-SVR and CPD-ws-mc-SVR. The above methods based on GPR and SVR in building the decomposed PES in the CPD form are implemented to construct a file named cpd that can be directly used for MCTDH or ML-MCTDH calculations implemented by the Heidelberg MCTDH soft package.¹¹ The code will be open-sourced in the near future. In the meantime, one can contact the authors for access. According to hierarchical coordinate frame, one can repeatedly substitute CPD-GPR/SVR function into V^{(CPD-GPR/SVR)} to obtain multi-layer version of CPD-GPR/SVR, called ML-CPD-GPR/SVR. These methods for the system with l layers predict the ranks of l·n and l·n_sv which are much larger than n and n_sv. Nevertheless, the ML-CPD-GPR/SVR function might be more in line with the multi-layer wave function and easier to simplify the working equations (see Section IV later).

Fig. 4 Comparisons of the numbers of training data n and support vectors n_sv, by (a) ratios n_sv/n for various potential energy intervals of H + H₂ and (b) ratios n_sv/n as a function of n for H₂ + H₂ and H₂/Cu(111). The tolerance error of these SVR calculations is set to be 10⁻⁴ Hartree. (a) The abscissa axis gives seven uniform intervals for the potential energy less than zero that is energy of the asymptotic channel of H + H₂. Each interval represents an energy range of ∼0.62 eV. The left axis represents the numbers of data while the right axis represents the ratio n_sv/n. The red and blue columns give the numbers of support vectors and training data, respectively. The black dots represent the ratios n_sv/n in percentage. (b) The abscissa and ordinate axes represent the number of training data n and ratios n_sv/n, respectively. For simplicity, the abscissa axis is given by a logarithmic scale instead of a linear scale. The ratio n_sv/n as a function of n for the H₂ + H₂ and H₂/Cu(111) systems is given by the blue and red lines, respectively. Similar information of n_sv/n for the H + H₂ system has been illustrated in (a). The support vectors are points that are most difficult to regress such that they have a direct bearing on the optimum location of the target prediction function. Thus, n_sv/n must be smaller than one. We refer the reader for ref. 85 for computational details. Reprinted from ref. 85. Copyright 2024 Zekai Miao, Xingyu Zhang, Qingfei Song, and Qingyong Meng.

Next, we turn to the reconstruction approaches, as collected in Table 1. One of the typical methods is POTFIT^64–66 which transfers the PES with less than six dimensionalities into the SOP form. For a higher-dimension (larger than six), the multi-grid POTFIT (MGPF),⁸⁶ multi-layer POTFIT (MLPF),^87,88 and Monte Carlo POTFIT (MCPF)⁸⁹ methods can be adopted to build the SOP form, while the Monte Carlo CPD re-fitting (MCCPD) method^4,90 is useful to build high-dimensional CPD functions. For instance, Otto and co-workers⁸⁸ implemented the MLPF method through various ways employing MGPF and then compared subsequent quantum dynamics calculations for H₃O₂⁻ with the accuracy of the re-fitting calculations. As shown by Fig. 5, the MLPF re-fitting provides ten times more accurate but only roughly consumes doubles the quantum dynamics runtime. For the second example, the PESs of CO/Cu(100)⁴ and H₂O/Cu(111)³ were re-fitted by MCCPD in this laboratory showing its power in subsequent dynamics calculations. Gatti and co-workers^91,92 performed 75D MCCPD calculations for the PES of a hydrogen atom scattering from graphene. Moreover, Peláez and co-workers^93,94 proposed a novel FBR formulated in the SOP/CPD form, denoted by SOP/CP-FBR. These methods enable the interpolation of a PES from coarse grids to fine grids, provided that the SPPs are represented by analytical functions and expanded through a given set of polynomial series, such as Chebyshev polynomial series {T_μ(·)}^s_μ=1,^93,94

Here, and c^(κ)_μ,r are the expansion coefficients. Due to the sparse character of the coefficient tensors (called the Tucker core tensors), the SOP/CP-FBR methods^93,94 save computational cost of the subsequent propagations. It is not surprising that, if kernel functions of CPD-GPR/SVR are expanded by Chebyshev basis functions with infinite order s_κ → ∞, then CPD-GPR/SVR becomes CP-FBR. Inspired by MLPF,^87,88 Peláez and co-workers⁹⁵ recently proposed and implemented a (binary) tree scheme for hierarchical Tucker (HT) decomposition of a high-dimensional dense target tensor, where a set of auxiliary basis is used to fit the children's factor matrices of eqn (13) and (14). This method is called HT-FBR. Their numerical experiments demonstrated that HT-FBR has faster performance than other approaches resulting in a more compact and optimizable PES and can be interpolated from a small grid sample to a large finer grid.⁹⁵

Fig. 5 Comparison of relaxation results with accuracy of the PES re-fitting calculation for the H₃O₂⁻ system, where (a) illustrates runtime for relaxation versus accuracy of the PES re-fit while (b) illustrates ground state energy versus accuracy of the PES re-fit. For MLPF, the full-grid MLPF (denoted by FG-MLPF) and random-sampling MLPF (denoted by RS-MLPF) methods are performed to compare the MLPF results with those of MGPF. By MGPF + MLPF, it meant that each of the MGPF calculations was then post-processed by MLPF. Comparing MGPF + MLPF with FG-MLPF, the first step of FG-MLPF was replaced by reading in the SOP re-fitting produced by MGPF which computed the SPPs. The runtime is consumed for the MCTDH improved relaxation or the ML-MCTDH regular relaxation. The error of the PES re-fit is estimated by Metropolis-Hastings sampling over 10⁶ points at an energy of 10⁴ cm⁻¹. (b) Data from the RS-MLPF fits are shown by the small grey dots and, where available, their uncertainties are indicated by error bars with one standard deviation. It should be noted that the energy axis of (b) uses linear scaling in the range 6603.29 ± 0.10 cm⁻¹ and logarithmic scaling outside this range. Reprinted with permission from ref. 88. Copyright 2018 Elsevier.

Next, we must turn to the tensor decomposition approaches (see Table 2) for building the potential tensor in the SOP or CPD form. In propagating the wave function, its discretization requires tensor representation of the PES. In this context, the PES is represented by the tensor with a d-way array. For instance, the potential functions shown by eqn (9) can be represented by their values on grids I = (i₁,⋯,i_d),

where means values of the i_κ-th grid for the κ-th DOF. The tensor decompositions V^(SOP)_I and V^(CPD)_I are called Tucker and CP decompositions, respectively. The Tucker decomposition is a higher-order form of principal component analysis while the CP decomposition is a sum of rank-one tensors.⁷⁵ Other tensor decomposition methods include tensor train (TT) decomposition, hierarchical Tucker (HT) decomposition, cross-approximation method, and so forth. In Table 3, we compare the methods given in Table 1 with the tensor decomposition approaches (see Table 2) for fitting the potential tensor in the sum of products of tensors. Baranov and Oseledets⁷⁸ applied the TT cross approximation procedure to the PES construction, in conjunction with an affine transformation of Cartesian coordinates into the active subspaces (AS) where the PES has the most variability. This TT + AS approach⁷⁸ belongs to the reconstruction approach which requires an existing PES or couples to a quantum chemistry code. An alternative way⁷⁸ to perform TT + AS is to sample geometries before, and to minimize the error between the approximation and the true values by coupling the interpolation code and the quantum chemistry software. Avila and Carrington⁷⁷ constructed the PES in the SOP form based on a Smolyak interpolation technique using polynomial-like or spectral basis functions and one-dimensional Lagrange-type functions. Habershon and co-workers⁹⁶ proposed an algorithm for tensor decomposition of existing KMR prediction for on-the-fly MCTDH. Furthermore, like the SOP/CP-FBR methods^93,94 various approaches were proposed combining tensorized Chebyshev interpolation with a function approximation procedure. For instance, Dolgov and co-workers⁷² combined tensorized Chebyshev interpolation with a Tucker decomposition of low multilinear rank to approximate a three-dimensional function defined on a tensor-product domain via function evaluations. Moreover, the functional tensor-train (FTT)⁷⁹ and its extended version (EFTT)⁸⁰ formats were proposed to compress and work with multidimensional functions on tensor product domains combining tensorized Chebyshev interpolation. Both FTT and EFTT are continuous extensions of the TT decomposition using a TT ansatz by replacing the three-dimensional TT cores with univariate matrix-valued functions.

Table 3 Comparison of the PES construction methods with the tensor decomposition approaches for fitting the potential tensor in the sum of direct products of one-dimensional tensors. The former is given in Table 1, while the latter is given in Table 2. Since we focus exclusively on the construction method irrespective of data provenance, no fitting-type distinction of re-construction versus direct construction is required here. Abbreviations of these methods are explained in the main text. The first column gives methods of the tensor decomposition. The second column gives correspondence methods for approximating the decomposed PES. The third column gives expansional expression of the corresponding methods. The rightmost column gives remarks on these methods, in particular the rank of the expanded PES

Tensor decomposition	Function decomposition	Expression	Remark
a MGPF and MCPF were designed for re-fitting high-dimensional (larger than six) PESs to the SOP form. b MLPF was designed for re-fitting high-dimensional (larger than six) PESs to the multilayer SOP form. c It means all variants of the CPD-GPR and CPD-SVR methods given in Table 1.
Tucker	POTFIT	V ^(SOP) in eqn (9)	POTFIT was designed for MCTDH
	MGPF		MGPF is multi-grid variant of POTFIT^a
	MCPF		MCPF is variant of POTFIT^a
HT	MLPF		MLPF is a multilayer variant of POTFIT^b
	ML-SOP-NN		Rank of HT is not clearly defined
	ML-CPD-GPR/SVR^c
CP	MCCPD	V ^(CPD) in eqn (9)	Rank should be carefully determined
	SOP-NN	Eqn (10)
	CPD-GPR/SVR^c	Eqn (11) and (12)
FTT, EFTT	SOP-FBR	Eqn (13)	Chebyshev basis is used to expand the SPPs in V^(SOP) or V^(CPD)
	CP-FBR	Eqn (14)

---

We now examine the complexity and rank of CPD-GPR/SVR in comparison with SOP-NN.^82,83 As given in Table 1, the ranks of SOP-NN and CPD-GPR are equal to the number of neurons and number of training data, respectively, both of which typically exceed the number of support vectors (rank of CPD-SVR).^84,85 The KMR approach offers distinct advantages through its Bayesian framework, enabling PES construction with reduced training data requirements.^84,97–99 Thus, the rank of CPD-GPR often remains comparable to that of SOP-NN, while CPD-SVR typically demonstrates lower rank due to n_sv < n. Consequently, CPD-GPR and SOP-NN exhibit similar levels of complexity and nonlinearity, with CPD-SVR^84,85 occupying an intermediate position. However, the comparison becomes more nuanced when considering ML-SOP-NN versus CPD-GPR/SVR. Extending this analysis to ML-CPD-GPR/SVR reveals a consistent pattern where SVR maintains its advantage in generating more compact expansion forms compared to both NN and GPR approaches. From the view of tensor decomposition, as indicated in Table 2, these methods exhibit progressively lower implementation complexity in the sequence from HT and Tucker decompositions to CP decomposition and to TT decomposition. This sequence is consistent with the aforementioned scenario if one notes that SOP-NN and CPD-GPR/SVR are based on Tucker decomposition and their multilayer counterparts adopt HT decomposition. The challenge of large-rank essentially arises from incomplete expansion of the PES and bears conceptual similarity to the electron-correlation problem in electronic structure theory. This analogy suggests both the necessity and feasibility of further improvements in building decomposed PES. This laboratory has been actively developing two distinct approaches to address this issue. The first approach leverages the inherent advantage of PES reconstruction in achieving reduced ranks through optimized expansion expression.^4,90 This has led to the development of a hybrid approach that combines direct construction and reconstruction ways to effectively mitigate large-rank challenges. The second innovative approach incorporates unsupervised learning techniques, say the chemistry-informed generative adversarial network (CI-GAN) method² proposed in this laboratory. This method² generates distributions of geometry and energy by n ≲ 10² training data, thereby offering the potential to produce a very small rank of ∼10².

IV. Wave function and propagation

Now, let us consider the expression of nuclear wave function and methods of wavepacket propagation. It has been described elsewhere^5–29 that the nuclear wave function can be expressed in a summation of products of single-particle functions (SPFs) as shown in eqn (2). Repeatedly expanding the SPFs, the wave function is given in the multilayer expansion form. From the point of numerical implementations, the multilayer expansion proves more advantageous than the single-layer expansion for wavepacket propagation exceeding nine dimensions or five atoms, enabling more efficient solutions of the SE equation. This has been widely found by practical calculations.^21–24 The reason might be two-fold, greater flexibility of the basis set and more efficient propagation of the basis functions, enabling a nearly complete expansion of the wavepacket with fewer SPFs (see eqn (2) as the single-layer expansion). First, the multilayer expansion enhances the flexibility of the basis set by hierarchically adapting to local features of the solution, which is critical for resolving the intricate structure of wavepackets. Second, the multilayer expansion enables more efficient propagation of the basis functions through dynamical low-rank approximation or tensor network techniques, effectively reducing the curse of dimensionality. By decomposing the wavepacket into a multilayer representation, such as a nested set of SPFs or a tree tensor network, one can achieve a nearly complete expansion with fewer basis functions than that in single-layer approaches. This not only improves computational performance but also preserves conserved quantities (say, norm of the wavepacket and total energy) more accurately during time evolution. This is because the adaptive basis can dynamically reorganize to capture dominant correlation patterns. In this context, the core advantages of this multi-layer expansion are arisen from mathematics and physics. Mathematically, it provides a more complete basis set for the wavepacket improving accuracy. Physically, the hierarchical structure enables the wavepacket to more easily “capture” the dynamical essence of the time evolution. We refer the reader to the review papers by Meyer and co-workers,^11,18 Wang,²⁵ Manthe,²⁷ and Larsson⁵ for details of this idea together with its historical evolution.

In multi-layer expansion,^{5,18–21,25,27} the SPFs in the (l − 1)-th layer are expanded by products of another series of SPFs in the deeper l-th layer (similar to eqn (2)),

where and , while is a mode. The SPFs of the last layer are finally expanded by a given set of primitive TI basis functions. It is convenient to introduce a diagrammatic notation, called the ML-tree structure (see Fig. 6 for the present three examples) to represent these objects.^19–21 By the number of expansional layers, one can classify propagation methods as collected in Table 4. It is worth noting that wavepacket propagation methods can be classified and named not only by their expansion form but also by the time integrators.^{19–21,100–106} Next, by substituting the wave function and Hamiltonian operator in the time-dependent variational principle (i.e., Dirac–Frenkel variational principle), Meyer, Manthe, Wang, and their co-workers^7,11,19–22 derived the working equations for propagating SPFs and expansion coefficients. Due to the equivalence of all layers, the working equations for each layer are formally identical.^19–21 This feature is particularly advantageous for high-dimensional systems. Thus, by a given time-integrator,^{19–21,100–106} the multi-layer expansion consequently provides a general framework of propagation, enabling the independent propagation of SPFs. Using the ML-MCTDH method, the present benchmarks have been extensively investigated at this laboratory^1–4,6 and the corresponding numerical results are illustrated in Fig. 7. We refer the reader to the SI for numerical details of these ML-MCTDH calculations. Recently, Gatti and co-workers^91,92 reported their 75D ML-MCTDH calculations for a hydrogen atom scattering from graphene, as illustrated in Fig. 8. These numerical results well demonstrate the power and capability of the ML-MCTDH method in computing the quantum dynamic properties of a molecular reaction, such as the spectrum and reaction probability.

Fig. 6 ML-MCTDH wave function structures (ML-tree structures) of the (a) CO/Cu(100),⁴ (b) H₂O/Cu(111),⁶ and (c) HO_x¹ systems. These ML-tree structures are iteratively optimized to indicate fast and well convergence of ML-MCTDH calculations.^22,23 Numbers of SPFs and primitive basis functions are given next to lines. Since the total angular momentum of the HO_x system is set to be zero,¹ its dynamical model unnecessarily incorporates the Euler angles {α^(HO_x),β^(HO_x),γ^(HO_x)} of the relative frame E₂/SF shown in Fig. 3(c).

Table 4 Comparison of wave functions between wavepacket propagation and electronic structure theory. The left and right panels give the properties of the propagation and electronic-structure methods, respectively. The leftmost column gives the number of layers of the expanded wave function. In each panel, the first column gives symbols, while the other columns give the coefficient, configuration, and permutation symmetry requirement of the entire wave function. The TI basis functions for propagation are also indicated here, while those for the electronic wave function is omitted

Layer^a	Wave function propagation				Electronic structure theory
	Symbol^b	Coeff.	Conf.^c	Sym.	Symbol	Coeff.	Conf.^c	Sym.
a Number of layer represents the number of expansions in building the configuration. b The symbol in the form of A/B represents theoretical approach A and the underlying basis set B. c We give the functions whose product consists of the configuration which is the product of a series of one-dimensional functions. d The zero-layer expansion means the configuration itself is the product of the basis functions. e The basis set of discrete variable representation (DVR) is adopted. f The DMRG is typically performed in Fock space. Thus, its antisymmetry is taken care of by the operators and not by the wavefunction.
Zero^d	TDH/FBR	TD function	TI basis	No	HF	Needs optimization	Basis	Antisymmetry
One	Standard/FBR	TD function	TI SPF	No	Full-CI	Needs optimization	MO	Antisymmetry
Two	MCTDH/DVR^e	TD function	TD SPF	No	MCSCF	Needs optimization	Needs optimization	Antisymmetry
Multi	ML-MCTDH/DVR^e	TD function	TD multi-layer SPF	No	TI-DMRG	—	—	Antisymmetry^f

---

Fig. 7 ML-MCTDH results of (a) and (b) the CO/Cu(100) system, (c) the H₂O/Cu(111) system, and (d) the HO_x system. (a) Time-dependent position expectations along out-of-plane coordinates of CO (the solid lines) and the top copper atom (the dashed line), while (b) shows sticking probabilities. The red and blue lines represent ML-MCTDH results from rigid (6D) and nonrigid (21D) surface models, respectively. Reprinted with permission from ref. 4. Copyright 2021 American Chemical Society. (c) Dissociative sticking probabilities computed for the rigid surface model. The colored lines were explained in ref. 6. Reprinted with permission from ref. 6. Copyright 2022 American Chemical Society. (d) Reactive probabilities for the OH + HO₂ → O₂ + H₂O reaction. Reprinted with permission from ref. 1. Copyright 2024 American Chemical Society.

Fig. 8 The 2D reduced density of the wave function for the z directions of the H atom and the q₁ normal mode at 10, 38, 46, 58, 78, 82, 120, and 260 fs during simulation on the PES 2D cut from the DOFs correspondent for the simulation of the H atom with an initial kinetic energy of 0.96 eV. We refer the reader to ref. 91 and 92 for numerical details and explanation. Reprinted with permission from ref. 91. Copyright 2023 American Institute of Physics.

One of the authors^22,23 suggested the necessity of optimizing the ML-tree structure to ensure efficient convergence during propagation. As rules of thumb,^{1,4,6,22–24,26} the following points should be noted when designing an initial ML-tree according to the hierarchical coordinate frames. First, two or three coordinates from a single BF frame in the deepest layer should be combined into one mode. Additionally, distance and angular coordinates should be separated into distinct modes. Second, if more than three coordinates are defined in a BF frame, they should be further divided into separate BF frames in deeper layers. Modes belonging to the same fragment may be grouped into a single mode in the upper layer. Third, the three Euler angles of a relative frame should form a single mode. Distance coordinates defining the reaction coordinate should be separated into individual modes, each containing only one coordinate. This separation is computationally advantageous as it reduces the number of grids required along the reaction coordinate to achieve well-converged propagation. The initial ML-tree is iteratively optimized^22,23 through extensive test calculations where one should note that the larger the mode populations, the larger are the couplings among the modes in the lower layers. For the present benchmarks, the optimized ML-tree structures shown in Fig. 6 differ from the hierarchical coordinate frames shown in Fig. 3. This discrepancy is expected as the coordinate frames are artificially designed. However, a comparison between Fig. 3 and 6 reveals that most of the separations are preserved suggesting physics insights behind the hierarchical coordinate frame. It should be mentioned that the above discrepancy of separating schemes between the wave function and coordinate frame decomposition does not affect the ML-MCTDH results. The tree structures indicate a separating scheme of the modes in either designing the coordinate frame or expanding the wave function. The resulting equations of motion (EOMs) of individual DOFs are essentially independent of the separating scheme. Therefore, the tree-structure difference between the coordinate frame and the wave function is acceptable in solving the working equations for propagation. Later, Vendrell and co-workers¹⁰⁷ proposed an approach to design an optimal ML-tree structure using multivariate statistics, or more specifically, factor analysis and hierarchical clustering. By the TN language, Hikihara and co-workers¹⁰⁸ proposed a tree algorithm to automatically optimize the network structure by local reconnections of isometries to suppress the bipartite entanglement entropy on their legs. This algorithm can be implemented to optimize the tree-structures in DMRG.

We now address the question of whether the multi-layer expansion form is sufficiently flexible for propagation. To this end, we compare the expansion form of the nuclear wave function with that of the electronic wave function,^55,60 as summarized in Table 4. We refer the reader to the SI for such comparison details. Of course, we must mention that the motion of electrons and nuclei during reactions is fundamentally different making electronic-structure theory and wavepacket propagation methods physically distinct. In general, the HF, full-CI, and MCSCF wave functions are expanded in zero-, one-, and two-layer forms, respectively, while the quantum chemistry TI-DMRG state is expanded in the multilayer form. We compare the TI-DMRG or called QC-DMRG^56,60 with ML-MCTDH using various criteria, as given in Table 5. Moreover, numerical comparisons^29,109,110 of ML-MCTDH with TD-DMRG have been reported by representing the state and Hamiltonian in terms of MPSs and MPOs, respectively. We noted that the method spectrum of DMRG (see also Fig. 9 as reference) contains a lot of methods which have been developing rapidly and integrating various numerical techniques (such as deep learning) to obtain numerous approaches tailored for diverse application scenarios. Considering the present goal, we have to omit this portion of the content on TI-DMRG and TD-DMRG. Here, we must caution the reader that the abbreviations TI-DMRG and TD-DMRG represent counterparts of wavepacket-based electronic-structure methods and quantum dynamics methods, respectively. Special attention should be paid to the critical distinction that the projector-splitting integrator developed by Lubich and co-workers^100–102 is exactly what modern TD-DMRG is. These above approaches for either wavepacket propagation or electronic structure employ the wave function in the CI type. This choice is natural because the CI-type wave function, which is expanded by a set of configurations as a basis set, is a general scheme for solving many-body EOMs under the SPA. To do this, the many-body EOMs are firstly reduced to single-particle EOMs by dimensionality reduction. Maximizing the size of the configuration basis set, the correlations among individual particles are incorporated as much as possible. In addition, let us turn to the choice of CC-type nuclear wave function. Loosely, defining reference state Φ₀ of an N-electron system, the CI and CC wave functions are given by,

where ĉ_p and t̂_p introduce the p-excitation term, while N is the expansional order. Eqn (17) clearly demonstrates the approximate equivalence¹¹¹ between Ψ_CI and Ψ_CC for finite N. The CC-type electronic wave function at N → ∞ and the exact solution of the time-dependent SE have equivalence formulation if the nuclear Hamiltonian is time-independent. In this context, the full-CI-type wave function expressed in terms of multilayer expansion is sufficient for propagation. Despite this, time-dependent coupled-cluster (TDCC) theory has been developed for simulating laser-driven electronic dynamics in atoms and molecules and for simulating molecular vibrational dynamics.¹¹² Although TDCC was developed for TD nucleon dynamics^113,114 and TD electronic-structure theory,^112,115,116 it shows applicability in the field of quantum molecular dynamics. Until now, we discussed nothing about perturbation theory, especially Rayleigh–Schrödinger perturbation theory (RSPT), where solutions of variational methods serve as references for further perturbations. When all configurations are considered in perturbation, RSPT evolves into many-body perturbation theory (MBPT). The CC-type wave function generally offers a very convenient resummation of MBPT diagrams offering an infinite-order approximation in selected cluster operators.¹¹¹ Specially, if t̂_p is a connected cluster operator corresponding to p-fold excitations, the CC and MBPT wave functions are approximately equivalent. In this context, the TD perturbation theory (TDPT) is widely employed to understand excitation in the interaction picture with TD perturbation. However, for propagation with a TI Hamiltonian, separating the Hamiltonian into a core term and a perturbation is often impossible, making few, if any, perturbation theory applicable.

Table 5 Comparison and summary of TI-DMRG and ML-MCTDH. The second column gives criteria to compare TI-DMRG with ML-MCTDH, including expansional forms of the Hamiltonian operator and wave function, optimization approaches, time integrators, system features, and application scenarios. Here, the first- and second-generation TI-DMRG methods are given, denoted by o-DMRG and s-DMRG, respectively. Since the DMRG theory was originally developed for field theory of many-body problems (say condensed matter systems), o-DMRG employs the renormalized or complementary Hamiltonian operators. In field theory, renormalized operators focus on regularization and finiteness, while complementary operators focus on symmetry and duality structures. In condensed matter systems, renormalized operators and complementary operators are generally distinct, but their mathematical expressions may overlap in high-symmetry systems

No.	Criteria	TI-DMRG		ML-MCTDH
		o-DMRG^a	s-DMRG^b
a Symbol “o-DMRG” means the original TI-DMRG theory, or first-generation formulation. b Symbol “s-DMRG” means the second-generation TI-DMRG theory defined by MPS and MPO. c Only adopted by MCTDH/DVR. MCTDH/CDVR is not restricted by the SOP or CPD form. d Sometimes called transformation matrices. e The MPSs are named as tensor train (TT) factorization from the view point of numerical analysis. f The sweep-based DMRG optimization uses the alternating least squares (ALS) algorithm in the TT context. g The original ML-MCTDH method is a propagation method. h The original DMRG method is an eigen-value solver. i It means interaction between fragments A and B of a complex system. j Only the original scenarios are given. k DMRG has capability to be applied to solve a wider range of equations but not limited to solve the Schrödinger equation.
1	Hamiltonian	Renormalized or complementary	MPO	SOP or CPD^c
2	Wave function	MPS^d	MPS^e	Multi-layer expansion
3	Optimization	ALS^f	ALS^f	N/A^g
4	Integrator	N/A^h	N/A^h	See ref. 19–21 and 100–106
5	Interaction^i	Short/long-range	Short/long-range	Short/long-range
6	Application^j	Various equations^k	Various equations^k	Nuclear dynamics

---

Fig. 9 Diagrammatic sketch for classification of the methods derived from time-independent (denoted by “TI” in the left panel) or time-dependent (denoted by “TD” in the right panel) variational principle and organized hierarchically with function (upper panel) or tensor (lower panel) representation. The black words indicate methods for molecular dynamics, while the red words indicate electronic-structure techniques. These methods exhibit mathematical tree-structures whether expressed in grid-based function representation or second-quantized representation or tensor representation, demonstrating deep relations among these methods. We would like to emphasize that the classification by “function” and “tensor” is informal serving only as an illustration. For instance, the original way of DMRG views the wave function in terms of renormalized functions and hence should be moved to the upper panel. This comparative framework reveals how the TI or TD variational principle unifies different representations, suggesting potential cross-methodological applications, in particular in addressing systems with correlated DOFs. This implies that insights from the tensor network (TN) or the tree tensor network (TTN) could inform those indicated by wave function, and vice versa, highlighting opportunities for methodological cross-fertilization between the electronic-structure and quantum molecular dynamics domains.

Nevertheless, CI-type, CC-type, and perturbative wave functions necessarily emerge when solving the molecular rovibrational eigenstate problems by the vibrational self-consistent field (VSCF) method and its enhanced variants, as given in Table 6. These methods variationally minimise a trial nuclear wave function composed of a single Hartree product of one-dimensional vibrational wave functions. Therefore, these vibrational-structure methods align more closely with the paradigm of electronic-structure theory than with that of dynamical methods which propagate the wave function in imaginary time to reach the rovibrational eigenstates. Furthermore, due to the wave function in the SOP form, these methods based on VSCF also require the Hamiltonian operator in a similar form, like the MCTDH and ML-MCTDH methods. We refer the reader to ref. 117–122 for VSCF and its variants. Like the relation between MCTDH and TI-DMRG, the vibrational-structure methods also have TI-DMRG counterparts.¹²³ In addition to TI vibrational-structure problems, these methods can be further extended to the TD propagations. For instance, in recent decades, Christiansen and co-workers^124–126 developed coupled cluster theory for time-dependent wave functions for the efficient computation of the quantum dynamics associated with the motion of nuclei, including time-dependent vibrational coupled cluster (TDVCC) and time-dependent modal vibrational coupled cluster (TDMVCC),^124–126 which employ static and adaptive basis sets, respectively. At present, developments of TDVCC and TDMVCC are still in their infancy. The fast configuration-space convergence, flexible choice of basis type and coordinate system, and polynomial-scaling computational cost make TDVCC and TDMVCC promising.^124–126

Table 6 A comparison of the approximations in resolving the vibration-structure problem based on the vibrational self-consistent field (VSCF) method. Resolving the vibration-structure problem, one could obtain ro-vibrational eigen-states of the system. The first column gives various methods for resolving the vibration-structure problem. The second column gives approximations in each method. The third column gives the feature of size consistency. The last two columns give counterparts in electron-structure theory as well as a strategy for resolving the equation of motion. We refer the reader to ref. 117–122 for details of VSCF and its variants

Method^a	Vibration structure		Electron structure
	Approximation^b	Size consistency	Counterpart^c	Strategy
a Abbreviations: “VMP” means vibrational Møller–Plesset perturbation theory, “PT2” and “VPT2” second order perturbation theory, “DC” degeneracy corrected version, “VCI” vibrational configuration-interaction, “ss” state-specific version, “VCC” vibrational coupled-cluster, “VMCSCF” and “VMRCC” multi-configurational version of VSCF and VCC, respectively. b Abbreviation “RSPT” means Rayleigh–Schrödinger perturbation theory. c Abbreviations: “HF SCF” means Hartree–Fock self-consistent field, “MPPT” Møller–Plesset perturbation theory, “MP2” second-order Møller–Plesset perturbation theory, “CI” configuration interaction, “CC” coupled cluster, “MCSCF” multi-configurational self-consistent field, and “MRCC” multi-reference coupled cluster.
VSCF	Hartree product of single-mode	No	HF SCF	Variation
VMP	RSPT based on VSCF wave function	No	MPPT	Perturbation
VSCF-PT2/VPT2	Second-order of VMP	No	MP2	Perturbation
VSCF-DCPT2	Degeneracy corrected VPT2	No	Degenerated MP2	Degenerated RSPT
HDCPT2	Hybrid version VSCF-DCPT2	No	Degenerated MP2	Degenerated RSPT
VCI	Expansion using a set of VSCF wave functions	Yes	Full CI	Variation
ss-VCI	State-specific configuration selected VCI	No	Truncated CI	Variation
VCC	Coupled cluster expansion using VSCF wave functions	Yes	Full CC	Variation
ss-VSCF	State-specific VSCF	No	—	Variation
VMCSCF	Multi-configurational version of VSCF	No	MCSCF	Variation
VMRCC	Multi-configurational version of VCC	No	MRCC	Variation

---

V. Perspectives on quantum molecular dynamics

In the present section, let us discuss the classifications of the methods for propagating wave function (as illustrated in Fig. 9) and then give perspectives on quantum molecular dynamics. In Table 7, we compare the propagation methods for molecular reaction dynamics with those for electronic and nucleon dynamics. All of these methods were designed to propagate wave function of either distinguishable particles (say atoms) or indistinguishable particles (say electrons and nucleons). As indicated in Section IV, substituting a solution ansatz into the time-dependent variational principle, one can derive the working equations. By solving these working equations, one can obtain the resulting dynamics properties, such as reaction probability and the scattering cross section. In these methods, the solution ansatz is expanded by products of one-dimensional (or lower-dimensional) functions, called SPFs or MOs in different scenarios. Correspondingly, the Hamiltonian operator is expanded into a similar summation of products of one-dimensional (or lower-dimensional) operators. This is helpful to reduce dimensionality of the working equations because multi-dimensional integrators are reduced to lower-dimensional ones. Repeatedly expanding, we can finally obtain a hierarchical framework for various dynamics scenarios. In Table 7, we compare the hierarchical frameworks for these scenarios. For instance, the time-dependent Hartree–Fock (TDHF) wave function for either electrons or nucleons has a zero-layer expansion form, similar to the standard propagation for chemical dynamics. Similar to the HF eigensolver, TDHF is an impractical tool due to its computational expense and lack of correlation effects making more layers necessary. Consequently, the time-dependent versions of CI, CC, and density functional theory (DFT) were subsequently developed to account for correlation effects leading to the TDCI, TDCC, and TDDFT methods for electrons or nucleons. Turning to the two-layer cases, the MCTDH-type expansion with permutation asymmetry can be used to propagate wave function of fermions leading to MCTDHF. In addition, there also exists the MCTDHB method for bosons. Therefore, if the hierarchical expansion preserves the permutation symmetry or antisymmetry, it can be used as a solution ansatz to derive a definitive propagation approach. In this context, due to the symmetry, field-theory techniques can be used to represent the wavepacket propagation making the second-quantization representation play an important role in dynamics calculations.

Table 7 Comparison and summary of quantum dynamics methods for various objects in a molecular system, together with those for nucleon dynamics (i.e., dynamics of protons and neutrons). In these methods, a solution ansatz is substituted into the time-dependent variational principle to obtain working equations. The solution ansatz is expanded by products (sometimes called configurations) of one-dimensional (or lower-dimensional) functions, called SPFs or MOs in different scenarios. Correspondingly, the Hamiltonian operator is expanded into a similar summation of products of one-dimensional (or lower-dimensional) operators. Repeating such expansional processes, one can finally obtain a hierarchical framework for various dynamics scenarios. The first column gives the layer of the expansional form of solution ansatz, where the SPFs or MOs in the deepest layer are further expanded by a given set of basis functions. The second and third columns give methods for wave-packet propagation for quantum molecular dynamics (QMD) and TD electronic-structure (TDES) theory, respectively. We refer the reader to Table 4 for comparison of the QMD methods with the TI electronic-structure methods which are essentially eigensolvers. The fourth column gives methods for quantum nucleon dynamics (QND) which are inherently within nuclear physics and hence beyond the scope of this work. The rightmost column gives probable remarks. Here, we would like to emphasize that the propagation methods for electrons or nucleons must preserve the permutation symmetry of TD wave function

Layer	QMD	TDES	QND	Remarks
a TDHF means time-dependent Hartree–Fock. b TDCI means time-dependent configuration-interaction. c TD-EOM-CC means TD equation-of-motion CC. d MCTDHF employs the wave function in the MCTDH form for dynamics of fermion particles.
Zero	Standard	TDHF^a	TDHF^a	An impractical tool due to its computational expense and lack of correlation effects
One	TDH	TDCI^b		Tool for improving time-dependent density functional theory (TDDFT)
		TDCC	TDCC	TDCC with >2 electrons does not converge to the full TDCI limit, and explains the emergence of plasmon behavior, a method for computing the absorption spectrum
		TD-EOM-CC^c
Two	MCTDH	MCTDHF^d	MCTDHF^d	There exists MCTDHB for boson dynamics
Multi	ML-MCTDH	TD-DMRG	TD-DMRG	MPS must preserve asymmetry

---

Following Section II and the above discussions, we now turn to the second-quantization representation (SQR) of quantum molecular dynamics. This is a way to apply quantum field theory to molecular dynamics, where the chemical dynamics is transferred to dynamics of quantum field. The occupation number vector (ONV) specifies which and how many SPFs are occupied by defining a vacuum state |0〉 alongside creation operator and annihilation operator for the κ-th coordinate together with commutators

where + and − mean anticommutators for fermions and commutators for bosons, respectively. If and are defined for the κ-th mode, the SPFs in the same layer can be expanded by products of states represented by acting a creation operator on |0〉. Repeating this process, the quantum state is finally represented in the multi-layer formalism. Since the occupation and coordinate representations are equivalent, the SQ Hamiltonian has the same decomposed structure in the matrix elements . Such hierarchical SQR might be helpful in simulating many-body open quantum systems¹²⁷ which have strong quantum correlations in both space and time. For instance, Tremblay and co-workers¹²⁸ implemented a Monte Carlo wavepacket (MCWP) algorithm in the MCTDH framework, called stochastic MCTDH (sMCTDH), and then simulated multi-dimensional dissipative quantum dynamics of O₂/Pt(111) in the Markovian regime through a Hamiltonian model of the Lindblad operator. In this MCWP implementation,¹²⁸ a bath-temperature (denoted by T) dependent dissipation operator is used for the κ-th DOF. Their sMCTDH simulations¹²⁸ found that thermalization can be approximately achieved if the intramolecular coupling is weak. In addition, Burghardt and co-workers¹²⁹ also reported numerous contributions on diffusive dynamics. We refer the reader to ref. 129 for further details.

Turning to SQR, over the past few years, Wang and Thoss,¹³⁰ Christiansen and co-workers,¹³¹ Manthe and Weike,^132,133 and Sasmal and Vendrell¹³⁴ reported their SQR for MCTDH/ML-MCTDH (called MCTDH-SQR/ML-MCTDH-SQR) or other methods for computing (i) a strongly correlated electronic model,¹³⁰ (ii) vibrational eigen-energy,¹³¹ (iii) imaginary time propagation,¹³³ and (iv) non-adiabatic spin dynamics.¹³⁴ Furthermore, many developments in MCTDH-SQR/ML-MCTDH-SQR were reported mainly concerned with dynamics and properties of the condensed matter rather than quantum molecular dynamics. This is because the occupation number representation in second quantization naturally describes collective behaviors of identical particles (see also eqn (18)). Employing the SQR, for instance, Christiansen and co-workers¹²⁶ developed the time-dependent full vibrational configuration interaction (TDFVCI), TDVCI[n], and TDVCC[n] methods where n means the number of the Hartree products in the solution ansatz. Usually, convergence of TDVCI[n] is slower compared to TDVCC[n], highlighting the advantage of the non-linear CC parameterization. We refer the reader to ref. 126 for details of TDVCI and TDVCC under the SQR and to Fig. 9 for comparison of TDVCI and TDVCC with other methods. In addition to the above field-theory techniques for particles, as given in Sections III and IV, another way to implement quantum molecular dynamics for large system,^29,109,110 such as material simulations, focuses on the methods with MPS and MPO, which is a fundamental form of TN or TTN and the most widely used representation for quantum many-body theory of condensed matter systems. In general, employing the grid-based representation the wave function and Hamiltonian operator in the decomposition form are represented by tensor decomposition. In this context, the methods with MPS and MPO should be equivalent to those with functional forms, as derived by Larsson.⁵ A popular method with MPS and MPO is the TD-DMRG method for propagating the nuclear wavepacket, while another approach is MPS-MCTDH.^28,69 Over the past decade, various comparisons^29,109,110 of such methods with ML-MCTDH were reported by modeling vibrational or photoelectron spectra because the DMRG method was originally developed to address challenges in simulating strongly correlated many-body systems. The TD-DMRG method was found to be highly accurate for dynamics of one-dimensional chain systems¹⁰⁹ and has been extended to higher-dimensional models.^29,110

The concepts of ONV, creation, and annihilation operators introduce virtual particles as spatial points along paths in a path-integral framework. This demonstrates an equivalence between quantum dynamics and path-integral molecular dynamics (PIMD).^135–140 It is a classical dynamics formulation in imaginary-time mathematically equivalent to a Wick rotation of path-integral quantum dynamics. It should be emphasized that the current PIMD and RPMD methods focus on temperature-dependent kinetics simulations, which corresponds to relaxation rather than propagation. By relaxation, the PIMD method enables the simulation of temperature-dependent chemical kinetics, extending beyond classical molecular dynamics to capture quantum statistical effects. To understand this point, for a one-dimensional system with action S = mẋ²/2 − V(x) and Hamiltonian Ĥ_ori = −∂_x²/(2m) + V(x), the propagator is given by

where the path from x to x′ is divided into P + 1 parts by points {x_i}^P_i=1, and the time from 0 to t is divided into ε + 1 parts such that Pε = t. By introducing imaginary time β = it, interpreted as inverse temperature, the partition function becomesyielding the effective Hamiltonian, known as the ring-polymer Hamiltonian,The ring-polymer Hamiltonian Ĥ^(P)_eff describes an ensemble of P replicas of the original system where two adjacent replicas interact through a harmonic-oscillator potential with frequency ω_P, as shown in Fig. 10. Over the past few decades, the power of ring-polymer molecular dynamics (RPMD) in simulating kinetics has been demonstrated by this laboratory and other groups.^135–140 When P → ∞, eqn (19) and (20) indicate equivalence between the d-dimensional quantum system and the (d + 1)-dimensional classical system, referred to as an equivalent classical system. However, we must mention that, this terminology is appropriate only in cases where all the path weights are positive. It will be true for systems for which quantum Monte Carlo (QMC) calculations can be performed (see also eqn (20)) in a practical implementation. In some cases, the (d + 1)-dimension effective model is an anisotropic analog of the classical model in d + 1 dimensions. In most cases, the equivalent classical system does not simply correspond to the original model in one higher dimension (comparing Ĥ_ori with Ĥ^(P)_eff). Rather, the path integral typically maps to a fundamentally different statistical mechanics problem in d + 1 dimensions, bearing no obvious resemblance to the original d-dimensional quantum system. In the last case, the d-dimensional quantum system is equivalent, on large length and time scales, to the same classical system in d + 1 dimensions.

Fig. 10 Diagrammatic sketch of the ring-polymer potential given by eqn (21). The blue symbols represent P virtual systems which are all identical and share the same Hamiltonian operator. The gray line represents the potential function V(x), where the abscissa and ordinate axes represent coordinate x and potential V, respectively. The red lines represent harmonic interaction between adjacent virtual systems while P is the number of virtual systems. The virtual systems can be seen as those assumed in the path integral. These interacted virtual systems form an ensemble if P → ∞.

Finally, we would like to mention the new developments in the application of generative artificial intelligence (AI) models to dynamics simulations. A typical example is the CI-GAN method.² In Section III, we mentioned that CI-GAN has the capability to directly build the PES in the CPD form with a smaller rank. In this case, the training database only contains geometries and associated energies without any time-dependent parameters. However, if the training data are computed in time-series order, the generated geometries are, in principle, in the time-series order. In other words, if the training data contain time or time-dependent quantity (say trajectory or wave function) it is possible to generate its time for each generated quantity (say geometry or wave function). With time and associated geometry one can actually obtain a trajectory. We would like to emphasize that this is not a simple and direct extension of the CI-GAN method² because the time index is certainly different from geometries and associated energies. Setting the beginning time to be zero, t₀ = 0, based on geometry X(t₀) as well as its small neighbourhood, one can generate geometries X(t) by CI-GAN. As shown above, these time-dependent geometries are expected to replace a classical trajectory. But, the predicted distributions of both time and geometry still require further developments and considerations. Very recently, we derived a generator for generation of the time index. However, if we turn to the time-dependent wave function, it is far more complex than the above discussions and the CI-GAN method.² One can expect that it is very difficult to implement for generating quantum dynamics.

VI. Conclusions

In this work, approaches to propagate nuclear wave functions are reviewed through a hierarchical framework under the coordinate representation. Its core concept is how to combine several coupled coordinates into a single mode making the coordinates defined in a hierarchical frame. By the hierarchical coordinate frame, the Hamiltonian operator and wave function are expanded by products of single-particle operators and functions, respectively. By substituting the multilayer solution ansatz to the Dirac–Frenkel variational principle, a set of coupled working equations for the expansion coefficients and the single-particle functions can be derived. First, the coordinate frame is hierarchically designed enabling the KEO to be expressed in a SOP form. Defining the coordinate frame and interested reaction coordinate, one can sample geometries based on which electronic-structure energy calculations are performed to compute the training dataset. Second, the PES is decomposed in two ways, direct construction by the training dataset and reconstruction through the existing PES. Thus, the Hamiltonian operator has been given in summation of products of single-particle operators. Third, the wave function is expanded repeatedly by a series of configurations which are products of SPFs. Since the Hamiltonian operator and wave function are both expressed in multilayer summation form, the working equations can be proved to be identical for all layers. This implies that the present hierarchical framework is useful for propagating a high-dimensional wave function. To show this point, the present benchmarks were extensively investigated at this laboratory^1–4,6 by the ML-MCTDH method. The numerical results of CO/Cu(100), H₂O/Cu(111), and HO_x are illustrated in Fig. 7. Moreover, 75D ML-MCTDH calculations for a hydrogen atom scattering from graphene^91,92 are also illustrated in Fig. 8 well demonstrating the power and capability of the ML-MCTDH method. In addition, the advantages and disadvantages of the present hierarchical framework are discussed to give perspectives on molecular reaction dynamics.

Conflicts of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

All data have been reported in this work. The data supporting this article have been included as part of the SI. Supplementary information: (1) Implementation details of the present multi-layer unified theory on the benchmarks of CO/Cu(100), H₂O/Cu(111), and HO_x; (2) details of the approaches to build the PES in general form, (3) details of the SOP-NN, ML-SOP-NN, CPD-GPR/SVR and ML-CPD-GPR/SVR methods, and (4) detailed perspectives on (i) comparisons to electronic-structure approaches, (ii) Fourier transformation, (iii) multi-layer SQR and PIMD, and (iv) the RPMD method. See DOI: https://doi.org/10.1039/d5cp02270c.

Acknowledgements

The authors gratefully acknowledge financial support from the National Natural Science Foundation of China (Grant No. 22273074), the Centre national de la recherche scientifique (CNRS) and International Research Network (IRN) “MCTDH”. The authors are grateful to Prof. Dr H.-D. Meyer and Dr M. Schröder (Universität Heidelberg) and Prof. Dr F. Gatti and Prof. Dr D. Peláez (Université Paris-Saclay) for helpful discussions and also grateful to anonymous reviewers for critical discussions.

References

1. Q. Song, X. Zhang, Z. Miao and Q. Meng, Construction of Mode-Combination Hamiltonian under the GridBased Representation for the Quantum Dynamics of OH + HO₂ → O₂ + H₂O, J. Chem. Theory Comput., 2024, 20, 597–613.
2. Z. Miao, X. Zhang, Y. Zhang, L. Wang and Q. Meng, Chemistry-Informed Generative Model for Classical Dynamics Simulations, J. Phys. Chem. Lett., 2024, 15, 532–539.
3. Q. Song, X. Zhang, D. Peláez and Q. Meng, Direct Canonical-Polyadic-Decomposition of the Potential Energy Surface from Discrete Data by Decoupled Gaussian Process Regression, J. Phys. Chem. Lett., 2022, 13, 11128–11135.
4. Q. Meng, M. Schröder and H.-D. Meyer, High-Dimensional Quantum Dynamics Study on Excitation-Specific Surface Scattering Including Lattice Effects of a Five-Atom Surface Cell, J. Chem. Theory Comput., 2021, 17, 2702–2713.
5. H. R. Larsson, A Tensor Network View of Multilayer Multiconfiguration Time-Dependent Hartree Methods, Mol. Phys., 2024, 122, e2306881.
6. Q. Song, X. Zhang, F. Gatti, Z. Miao, Q. Zhang and Q. Meng, Multilayer Multiconfiguration Time-Dependent Hartree Study on the Mode-/Bond-Specific Quantum Dynamics of Water Dissociation on Cu(111), J. Phys. Chem. A, 2022, 126, 6047–6058.
7. H.-D. Meyer, U. Manthe and L. S. Cederbaum, The Multi-Configurational Time-Dependent Hartree Approach, Chem. Phys. Lett., 1990, 165, 73–78.
8. U. Manthe, H.-D. Meyer and L. S. Cederbaum, The Multi-Configurational Time-Dependent Hartree Approach, J. Chem. Phys., 1992, 97, 3199–3213.
9. H.-D. Meyer, U. Manthe and L. S. Cederbaum, The multi-configuration Hartree approach, in Numerical Grid Methods and their Application to Schrödinger's Equation, ed. C. Cerjan, Kluwer Academic Publishers, Dordrecht, 1993, pp. 141–152.
10. H.-D. Meyer, Multiconfiguration time-dependent Hartree method, in The Encyclopedia of Computational Chemistry, ed. P. V. R. Schleyer, N. L. Allinger, T. Clark, J. Gasteiger, P. A. Kollman, H. F. Schaefer III and P. R. Schreiner, John Wiley and Sons, Chichester, 1998, vol. 5, pp. 3011–3018.
11. M. H. Beck, A. Jäckle, G. A. Worth and H.-D. Meyer, The Multi-Configuration Time-Dedependent Hartree (MCTDH), Phys. Rep., 2000, 324, 1–105.
12. H.-D. Meyer and G. A. Worth, Quantum Molecular Dynamics: Propagating Wavepackets and Density Operators Using the Multiconfiguration Time-Dependent Hartree (MCTDH) Method, Theor. Chem. Acc., 2003, 109, 251–267.
13. A. Viel, R. P. Krawczyk, U. Manthe and W. Domcke, The Sudden-Polarization Effect and its Role in the Ultrafast Photochemistry of Ethene, Angew. Chem., Int. Ed., 2003, 42, 3434–3436.
14. T. Wu and U. Manthe, A Potential Energy Surface Construction Scheme for Accurate Reaction Rate Calculations: General Approach and a Test for the H + CH₄ → H₂ + CH₃ Reaction, J. Chem. Phys., 2003, 119, 14–23.
15. T. Wu, H.-J. Werner and U. Manthe, First-Principles Theory for the H + CH₄ → H₂ + CH₃ Reaction, Science, 2004, 306, 2227–2229.
16. H. Wang and M. Thoss, Numerically Exact Quantum Dynamics for Indistinguishable Particles: The Multilayer Multiconfiguration Time-Dependent Hartree Theory in Second Quantization Representation, J. Chem. Phys., 2009, 131(2), 024114.
17. R. Wodraszka and U. Manthe, A Multi-Configurational Time-Dependent Hartree Approach to the Eigenstates of a Multi-Well System, J. Chem. Phys., 2012, 136, 124119.
18. H.-D. Meyer, Studying Molecular Quantum Dynamics with the Multiconfiguration Time-Dependent Hartree Method, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2012, 2, 351–374.
19. H. Wang and M. Thoss, Multilayer Formulation of the Multiconfiguration Time-Dependent Hartree Theory, J. Chem. Phys., 2003, 119, 1289–1299.
20. U. Manthe, A Multilayer Multiconfigurational Time-Dependent Hartree Approach for Quantum Dynamics on General Potential Energy Surfaces, J. Chem. Phys., 2008, 128, 164116.
21. O. Vendrell and H.-D. Meyer, Multilayer Multiconfiguration Time-Dependent Hartree Method: Implementation and Applications to a Henon-Heiles Hamiltonian and to Pyrazine, J. Chem. Phys., 2011, 134, 044135.
22. Q. Meng, S. Faraji, O. Vendrell and H.-D. Meyer, Full Dimensional Quantum-Mechanical Simulations for the Vibronic Dynamics of Diflurorbenzene Radical Cation Isomers using the Multilayer Multiconfiguration Time-Dependent Hartree Method, J. Chem. Phys., 2012, 137, 134302.
23. Q. Meng and H.-D. Meyer, A Multilayer MCTDH Study on the Full Dimensional Vibronic Dynamics of Naphthalene and Anthracene Cations, J. Chem. Phys., 2013, 138, 014313.
24. Q. Meng and H.-D. Meyer, A Full-Dimensional Multilayer Multiconfiguration Time-Dependent Hartree Study on the Ultraviolet Absorption Spectrum of Formaldehyde Oxide, J. Chem. Phys., 2014, 141, 124309.
25. H. Wang, Multilayer Multiconfiguration Time-Dependent Hartree Theory, J. Phys. Chem. A, 2015, 119, 7951.
26. Q. Meng and H.-D. Meyer, Lattice Effects of Surface Cell: Multilayer Multiconfiguration Time-Dependent Hartree Study on Surface Scattering of CO/Cu(100), J. Chem. Phys., 2017, 146, 184305.
27. U. Manthe, Wavepacket Dynamics and the Multi-Configurational Time-Dependent Hartree Approach, J. Phys.:Condens. Matter, 2017, 29, 253001.
28. Y. Kurashige, Matrix Product State Formulation of the Multiconfiguration Time-Dependent Hartree Theory, J. Chem. Phys., 2018, 149, 194114.
29. M. F. X. Dorfner, D. Brey, I. Burghardt and F. Ortmann, Comparison of Matrix Product State and Multiconfiguration Time-Dependent Hartree Methods for Nonadiabatic Dynamics of Exciton Dissociation, J. Chem. Theory Comput., 2024, 20, 8767–8781.
30. Z. Bačić and J. C. Light, Theoretical Methods for Rovibrational States of Floppy Molecules, Ann. Rev. Phys. Chem., 1989, 40, 469.
31. R. van Harrevelt and U. Manthe, Multidimensional Time-Dependent Discrete Variable Representations in Multiconfiguration Hartree Calculations, J. Chem. Phys., 2005, 123, 064106.
32. U. Manthe, Layered Discrete Variable Representations and their Application within the Multiconfigurational Time-Dependent Hartree Approach, J. Chem. Phys., 2009, 130, 054109.
33. R. Ellerbrock and U. Manthe, A Non-Hierarchical Correlation Discrete Variable Representation, J. Chem. Phys., 2022, 156, 134107.
34. R. Ellerbrock, H. Hoppe and U. Manthe, A Non-Hierarchical Multi-Layer Multi-Configurational Time-Dependent Hartree Approach for Quantum Dynamics on General Potential Energy Surfaces, J. Chem. Phys., 2024, 160, 224108.
35. D. B. Dix, Polyspherical Coordinate Systems on Orbit Spaces with Applications to Biomolecular Shape, Acta Appl. Math., 2006, 90, 247–306.
36. F. Gatti and C. Iung, Exact and Constrained Kinetic Energy Operators for Polyatomic Molecules: The Polyspherical Approach, Phys. Rep., 2009, 484, 1–69.
37. D. Lauvergnat and A. Nauts, Exact Numerical Computation of a Kinetic Energy Operator in Curvilinear Coordinates, J. Chem. Phys., 2002, 116, 8560.
38. D. Lauvergnat, E. Baloïtcha, G. Dive and M. Desouter-Lecomte, Dynamics of Complex Molecular Systems with Numerical Kinetic Energy Operators in Generalized Coordinates, Chem. Phys., 2006, 326, 500–508.
39. M. Ndong, L. Joubert Doriol, H.-D. Meyer, A. Nauts, F. Gatti and D. Lauvergnat, Automatic Computer Procedure for Generating Exact and Analytical Kinetic Energy Operators Based on the Polyspherical Approach, J. Chem. Phys., 2012, 136, 034107.
40. M. Ndong, A. Nauts, L. Joubert-Doriol, H.-D. Meyer, F. Gatti and D. Lauvergnat, Automatic Computer Procedure for Generating Exact and Analytical Kinetic Energy Operators Based on the Polyspherical Approach: General Formulation and Removal of Singularities, J. Chem. Phys., 2013, 139, 204107.
41. Z. Zhao, J. Chen, Z. Zhang, D. H. Zhang, D. Lauvergnat and F. Gatti, Full-Dimensional Vibrational Calculations of Five-Atom Molecules Using A Combination of Radau and Jacobi Coordinates: Applications to Methane and Fluoromethane, J. Chem. Phys., 2016, 144, 204302.
42. Y.-Y. Shi, L.-M. Duan and G. Vidal, Classical Simulation of Quantum Many-Body Systems with A Tree Tensor Network, Phys. Rev. A, 2006, 74, 022320.
43. S. R. White and A. E. Feiguin, Real-Time Evolution Using the Density Matrix Renormalization Group, Phys. Rev. Lett., 2004, 93, 076401.
44. A. J. Daley, C. Kollath, U. Schollwöck and G. Vidal, Time-Dependent Density-Matrix Renormalization-Group Using Adaptive Effective Hilbert Spaces, J. Stat. Mech., 2004, 2004, P04005.
45. G. Vidal, Efficient Simulation of One-Dimensional Quantum Many-Body Systems, Phys. Rev. Lett., 2004, 93, 040502.
46. A. E. Feiguin and S. R. White, Time-Step Targeting Methods for Real-Time Dynamics Using the Density Matrix Renormalization Group, Phys. Rev. B:Condens. Matter Mater. Phys., 2005, 72, 020404.
47. J. José García-Ripoll, Time evolution of Matrix Product States, New J. Phys., 2006, 8, 305.
48. K. A. Al-Hassanieh, A. E. Feiguin, J. A. Riera, C. A. Büsser and E. Dagotto, Adaptive Time-Dependent Density-Matrix Renormalization-Group Technique for Calculating the Conductance of Strongly Correlated Nanostructures, Phys. Rev. B:Condens. Matter Mater. Phys., 2006, 73, 195304.
49. E. Ronca, Z. Li, C. A. Jimenez-Hoyos and G. K.-L. Chan, Time-Step Targeting Time-Dependent and Dynamical Density Matrix Renormalization Group Algorithms with ab initio Hamiltonians, J. Chem. Theory Comput., 2017, 13, 5560–5571.
50. D. Iouchtchenko and P.-N. Roy, Ground States of Linear Rotor Chains via the Density Matrix Renormalization Group, J. Chem. Phys., 2018, 148, 134115.
51. Q. Shi, Y. Xu, Y. Yan and M. Xu, Gaussian-Based Multiconfiguration Time-Dependent Hartree: A Two-Layer Approach. III. Application to Nonadiabatic Dynamics in a Charge Transfer Complex, J. Chem. Phys., 2018, 148, 174102.
52. A. Baiardi and M. Reiher, Large-Scale Quantum Dynamics with Matrix Product States, J. Chem. Theory Comput., 2019, 15, 3481–3498.
53. X. Dan and Q. Shi, Theoretical Study of Nonadiabatic Hydrogen Atom Scattering Dynamics on Metal Surfaces Using the Hierarchical Equations of Motion Method, J. Chem. Phys., 2023, 159, 044101.
54. S. Bai, S. Zhang, C. Huang and Q. Shi, Hierarchical Equations of Motion for Quantum Chemical Dynamics: Recent Methodology Developments and Applications, Acc. Chem. Res., 2024, 57, 3151–3160.
55. A. Szabo and N. S. Ostlund, Modern Quantum Chemistry, Dover, New York, 1996.
56. S. Szalay, M. Pfeffer, V. Murg, G. Barcza, F. Verstraete, R. Schneider and Ö. Legeza, Tensor Product Methods and Entanglement Optimization for ab initio Quantum Chemistry, Int. J. Quantum Chem., 2015, 115, 1342–1391.
57. S. R. White, Density Matrix Formulation for Quantum Renormalization Groups, Phys. Rev. Lett., 1992, 69, 2863.
58. S. White, Density-Matrix Algorithms for Quantum Renormalization Groups, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 48, 10345.
59. G. K.-L. Chan, Low Entanglement Wavefunctions, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2012, 2, 907–920.
60. S. Wouters and D. Van Neck, The Density Matrix Renormalization Group for ab initio Quantum Chemistry, Eur. Phys. J. D, 2014, 68, 272.
61. Q. Meng and H.-D. Meyer, MCTDH Study on Vibrational States of the CO/Cu(100) System, J. Chem. Phys., 2013, 139, 164709.
62. N. R. Hanson, The Mathematical Power of Epicyclical Astronomy, Isis, 1960, 51, 150–158.
63. T. Hollebeek, T.-S. Ho and H. Rabitz, A Fast Algorithm for Evaluating Multidimensional Potential Energy Surfaces, J. Chem. Phys., 1997, 106, 7223–7227.
64. A. Jäckle and H.-D. Meyer, Reactive Scattering using the Multiconfiguration Time-Dependent Hartree Approximation: General Aspects and Application to the Collinear H + H₂ → H₂ + H Reaction, J. Chem. Phys., 1995, 102, 5605.
65. A. Jäckle and H.-D. Meyer, Product Representation of Potential Energy Surfaces, J. Chem. Phys., 1996, 104, 7974.
66. A. Jäckle and H.-D. Meyer, Product Representation of Potential Energy Surfaces II, J. Chem. Phys., 1998, 109, 3772.
67. F. Verstraete, J. J. García-Ripoll and J. I. Cirac, Matrix Product Density Operators: Simulation of Finite-Temperature and Dissipative Systems, Phys. Rev. Lett., 2004, 93, 207204.
68. B. Pirvu, V. Murg, J. I. Cirac and F. Verstraete, Matrix Product Operator Representations, New J. Phys., 2010, 12, 025012.
69. K. Hino and Y. Kurashige, Encoding a Many-Body Potential Energy Surface into a Grid-Based Matrix Product Operator, J. Chem. Theory Comput., 2024, 20, 3839–3849.
70. J. D. Carroll and J.-J. Chang, Psychometrika, 1970, 35, 283–319.
71. L. Tucker, Analysis of Individual Differences in Multidimensional Scaling Via an N-way Generalization of “Eckart-Young” Decomposition, Psychometrika, 1966, 31, 279.
72. S. Dolgov, D. Kressner and C. Strössner, Functional Tucker Approximation Using Chebyshev Interpolation, SIAM J. Sci. Comput., 2021, 43(3), A2190–A2210.
73. W. Hackbusch and S. Kühn, A New Scheme for the Tensor Representation, J. Fourier Anal. Appl., 2009, 15, 706–722.
74. C. Lubich, T. Rohwedder, R. Schneider and B. Vandereycken, Dynamical Approximation by Hierarchical Tucker and Tensor-Train Tensors, SIAM J. Matrix Anal. Appl., 2013, 34, 470–494.
75. T. G. Kolda and B. W. Bader, Tensor Decompositions and Applications, SIAM Rev., 2009, 51, 455–500.
76. I. Oseledets and E. Tyrtyshnikov, TT-Cross Approximation for Multidimensional Arrays, Linear Algebra Appl., 2010, 432, 70–88.
77. G. Avila and T. J. Carrington, Using Multi-Dimensional Smolyak Interpolation to Make a Sum-of-Products Potential, J. Chem. Phys., 2015, 143, 044106.
78. V. Baranov and I. Oseledets, Fitting High-Dimensional Potential Energy Surface Using Active Subspace and Tensor Train (AS + TT) Method, J. Chem. Phys., 2015, 143, 174107.
79. A. Gorodetsky, S. Karaman and Y. Marzouk, A Continuous Analogue of the Tensor-Train Decomposition, Comput. Methods Appl. Mech. Eng., 2019, 347, 59–84.
80. C. Strössner, B. Sun and D. Kressner, Approximation in the Extended Functional Tensor Train Format, Adv. Comput. Math., 2024, 50, 54.
81. S. A. Goreinov, On Cross Approximation of Multi-Index Arrays, Dokl. Math., 2008, 77, 404–406.
82. S. Manzhos and T. Carrington, Jr., Using Neural Networks to Represent Potential Surfaces as Sum of Products, J. Chem. Phys., 2006, 125, 194105.
83. W. Koch and D. H. Zhang, Separable Potential Energy Surfaces from Multiplicative Artificial Neural Networks, J. Chem. Phys., 2014, 141, 021101.
84. Q. Song, X. Zhang, Z. Miao, Q. Zhang and Q. Meng, Unified Regression Models in Fitting Potential Energy Surfaces for Quantum Dynamics, J. Math. Chem., 2022, 60, 1983–2012.
85. Z. Miao, X. Zhang, Q. Song and Q. Meng, Canonical-Polyadic-Decomposition of the Potential Energy Surface Fitted by Warm-Started Support Vector Regression, arXiv, 2024, preprint, arXiv:2410.23529.
86. D. Peláez and H.-D. Meyer, The Multigrid POTFIT (MGPF) Method: Grid Representations of Potentials for Quantum Dynamics of Large Systems, J. Chem. Phys., 2013, 138, 014108.
87. F. Otto, Multi-Layer POTFIT: An Accurate Potential Representation for Efficient High-Dimensional Quantum Dynamics, J. Chem. Phys., 2014, 140, 014106.
88. F. Otto, Y.-C. Chiang and D. Peláez, Accuracy of Potfit-based potential representations and its impact on the performance of (ML-)MCTDH, Chem. Phys., 2018, 509, 116–130.
89. M. Schröder and H.-D. Meyer, Transforming High-Dimensional Potential Energy Surfaces into Sum-of-Products form using Monte Carlo Methods, J. Chem. Phys., 2017, 147, 064105.
90. M. Schröder, Transforming High-Dimensional Potential Energy Surfaces into a Canonical Polyadic Decomposition using Monte Carlo Methods, J. Chem. Phys., 2020, 152, 024108.
91. L. Shi, M. Schröder, H.-D. Meyer, D. Peláez, A. M. Wodtke, K. Golibrzuch, A.-M. Schönemann, A. Kandratsenka and F. Gatti, Quantum and Classical Molecular Dynamics for H Atom Scattering from Graphene, J. Chem. Phys., 2023, 159, 194102.
92. L. Shi, M. Schröder, H.-D. Meyer, D. Peláez, A. M. Wodtke, K. Golibrzuch, A.-M. Schönemann, A. Kandratsenka and F. Gatti, Full Quantum Dynamics Study for H Atom Scattering from Graphene, J. Phys. Chem. A, 2025, 129, 1896–1907.
93. R. L. Panadés-Barrueta and D. Peláez, Low-Rank Sum-of-Products Finite-Basis-Representation (SOP-FBR) of Potential Energy Surfaces, J. Chem. Phys., 2020, 153, 234110.
94. N. Nadoveza, R. L. Panadés-Barrueta, L. Shi, F. Gatti and D. Peláez, Analytical High-Dimensional Operators in Canonical Polyadic Finite Basis Representation (CP-FBR), J. Chem. Phys., 2023, 158, 114109.
95. Z. Qiu, F. Magoulés and D. Peláez, Single-Entry Computation of Analytical Hierarchical (Binary) Tree Structures, Adv. Eng. Software, 2025, 203, 103873.
96. G. W. Richings and S. Habershon, Predicting Molecular Photochemistry Using Machine-Learning-Enhanced Quantum Dynamics Simulations, Acc. Chem. Res., 2022, 55, 209–220.
97. A. Kamath, R. A. Vargas-Hernández, R. V. Krems, T. Carrington and S. Manzhos, Neural Networks vs. Gaussian Process Regression for Representing Potential Energy Surfaces: A Comparative Study of Fit Quality and Vibrational Spectrum Accuracy, J. Chem. Phys., 2018, 148, 241702.
98. Q. Song, Q. Zhang and Q. Meng, Neural-Network Potential Energy Surface with Small Database and High Precision: A Benchmark of the H + H₂ System, J. Chem. Phys., 2019, 151, 114302.
99. Q. Song, Q. Zhang and Q. Meng, Revisiting the Gaussian Process Regression for Fitting High-Dimensional Potential Energy Surface and Its Application to the OH + HO₂ → O₂ + H₂O Reaction, J. Chem. Phys., 2020, 152, 134309.
100. C. Lubich, A Variational Splitting Integrator for Quantum Molecular Dynamics, Appl. Numer. Math., 2004, 48, 355–368.
101. C. Lubich, I. V. Oseledets and B. Vandereycken, Time Integration of Tensor Trains, SIAM J. Numer. Anal., 2015, 53, 917–941.
102. B. Kloss, I. Burghardt and C. Lubich, Implementation of a Novel Projector-Splitting Integrator for the Multi-Configurational Time-Dependent Hartree Approach, J. Chem. Phys., 2017, 146, 174107.
103. M. Bonfanti and I. Burghardt, Tangent Space Formulation of the Multi-Configuration Time-Dependent Hartree Equations of Motion: The Projector–Splitting Algorithm Revisited, Chem. Phys., 2018, 515, 252–261.
104. L. P. Lindoy, B. Kloss and D. R. Reichman, Time Evolution of ML-MCTDH Wavefunctions. I. Gauge Conditions, Basis Functions, and Singularities, J. Chem. Phys., 2021, 155, 174108.
105. L. P. Lindoy, B. Kloss and D. R. Reichman, Time Evolution of ML-MCTDH Wavefunctions. II. Application of the Projector Splitting Integrator, J. Chem. Phys., 2021, 155, 174109.
106. H. Wang and H.-D. Meyer, On Regularizing the ML-MCTDH Equations of Motion, J. Chem. Phys., 2018, 149, 044119.
107. D. Mendive-Tapia, H.-D. Meyer and O. Vendrell, Optimal Mode Combination in the Multiconfiguration Time-Dependent Hartree Method through Multivariate Statistics: Factor Analysis and Hierarchical Clustering, J. Chem. Theory Comput., 2023, 19, 1144–1156.
108. T. Hikihara, H. Ueda, K. Okunishi, K. Harada and T. Nishino, Automatic Structural Optimization of Tree Tensor Networks, Phys. Rev. Res., 2023, 5, 013031.
109. S. Mainali, F. Gatti, D. Iouchtchenko, P.-N. Roy and H.-D. Meyer, Comparison of the Multi-Layer Multi-Configuration Time-Dependent Hartree (ML-MCTDH) Method and the Density Matrix Renormalization Group (DMRG) for Ground State Properties of Linear Rotor Chains, J. Chem. Phys., 2021, 154, 174106.
110. H. R. Larsson, M. Schröder, R. Beckmann, F. Brieuc, C. Schran, D. Marx and O. Vendrell, State-Resolved Infrared Spectrum of the Protonated Water Dimer: Revisiting the Characteristic Proton Transfer Doublet Peak, Chem. Sci., 2022, 13, 11119–11125.
111. R. J. Bartlett and M. Musiał, Coupled-Cluster Theory in Quantum Chemistry, Rev. Mod. Phys., 2007, 79, 291–352.
112. B. Sverdrup Ofstad, E. Aurbakken, Ø. Sigmundson Schøyen, H. E. Kristiansen, S. Kvaal and T. B. Pedersen, Time-Dependent Coupled-Cluster Theory, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2023, 13, e1666.
113. P. Hoodbhoy and J. W. Negele, Time-Dependent Coupled-Cluster Approximation to Nuclear Dynamics. I. Application to A Solvable Model, Phys. Rev. C, 1978, 18, 2380–2394.
114. P. Hoodbhoy and J. W. Negele, Time-Dependent Coupled-Cluster Approximation to Nuclear Dynamics. II. General Formulation, Phys. Rev. C, 1979, 19, 1971–1982.
115. J. J. Goings, P. J. Lestrange and X. Li, Real-Time Time-Dependent Electronic Structure Theory, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2018, 8, e1341.
116. X. Li, N. Govind, C. Isborn, A. E. DePrince III and K. Lopata, Real-Time Time-Dependent Electronic Structure Theory, Chem. Rev., 2020, 120, 9951–9993.
117. J. M. Bowman, S. Carter and X. Huang, MULTIMODE: A Code to Calculate Rovibrational Energies of Polyatomic Molecules, Int. Rev. Phys. Chem., 2003, 22, 533–549.
118. O. Christiansen, A Second Quantization Formulation of Multimode Dynamics, J. Chem. Phys., 2004, 120, 2140–2148.
119. L. Pele and R. B. Gerber, On the Mean Accuracy of the Separable VSCF Approximation for Large Molecules, J. Phys. Chem. C, 2010, 114, 20603–20608.
120. P. Seidler, M. Sparta and O. Christiansen, Vibrational Coupled Cluster Response Theory: A General Implementation, J. Chem. Phys., 2011, 134, 054119.
121. J. A. Faucheaux and S. Hirata, Higher-order Diagrammatic Vibrational Coupled-Cluster Theory, J. Chem. Phys., 2015, 143, 134105.
122. J. A. Faucheaux, M. Nooijen and S. Hirata, Similarity-Transformed Equation-of-Motion Vibrational Coupled-Cluster Theory, J. Chem. Phys., 2018, 148, 054104.
123. N. Glaser, A. Baiardi and M. Reiher, Flexible DMRG-Based Framework for Anharmonic Vibrational Calculations, J. Chem. Theory Comput., 2023, 19, 9329–9343.
124. M. G. Højlund, A. Zoccante and O. Christiansen, Time-Dependent Coupled Cluster with Orthogonal Adaptive Basis Functions: General Formalism and Application to the Vibrational Problem, J. Chem. Phys., 2024, 160, 024105.
125. R. B. Jensen and O. Christiansen, Unitary Vibrational Coupled Cluster: General Theory and Implementation, J. Chem. Phys., 2025, 162, 084112.
126. M. G. Højlund, A. Zoccante, A. B. Jensen and O. Christiansen, Time-Dependent Vibrational Coupled Cluster Theory With Static and Dynamic Basis Functions, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2025, 15, e70001.
127. L. Ye, Y. Wang and X. Zheng, Simulating Many-Body Open Quantum Systems by Harnessing the Power of Artificial Intelligence and Quantum Computing, J. Chem. Phys., 2025, 162, 120901.
128. S. Mandal, F. Gatti, O. Bindech, R. Marquardt and J.-C. Tremblay, Multidimensional Stochastic Dissipative Quantum Dynamics Using A Lindblad Operator, J. Chem. Phys., 2022, 156, 094109.
129. W. Popp, D. Brey, R. Binder and I. Burghardt, Quantum Dynamics of Exciton Transport and Dissociation in Multichromophoric Systems, Annu. Rev. Phys. Chem., 2021, 72, 591–616.
130. H. Wang and M. Thoss, A Multilayer Multiconfiguration Time-Dependent Hartree Study of the Nonequilibrium Anderson Impurity Model at Zero Temperature, Chem. Phys., 2018, 509, 13–19.
131. N. K. Madsen, M. B. Hansen, A. Zoccante, K. Monrad, M. B. Hansen and O. Christiansen, Exponential Parameterization of Wave Functions for Quantum Dynamics: Time-Dependent Hartree in Second Quantization, J. Chem. Phys., 2018, 149, 134110.
132. U. Manthe and T. Weike, On the Multi-Layer Multi-Configuration Time-Dependent Hartree Approach for Bosons and Fermions, J. Chem. Phys., 2017, 146, 064117.
133. T. Weike and U. Manthe, The Multi-Configurational Time-Dependent Hartree Approach in Optimized Second Quantization: Imaginary Time Propagation and Particle Number Conservation, J. Chem. Phys., 2020, 152, 034101.
134. S. Sasmal and O. Vendrell, Non-Adiabatic Quantum Dynamics without Potential Energy Surfaces based on Second-Quantized Electrons: Application within the Framework of the MCTDH Method, J. Chem. Phys., 2020, 153, 154110.
135. Y. V. Suleimanov, R. Collepardo-Guevara and D. E. Manolopoulos, Bimolecular Reaction Rates from Ring Polymer Molecular Dynamics: Application to H + CH₄ → H₂ + CH₃, J. Chem. Phys., 2011, 134, 044131.
136. J. W. Allen, W. H. Green, Y. Li, H. Guo and Y. V. Suleimanov, Full Dimensional Quantum Rate Coefficients and Kinetic Isotope Effects from Ring Polymer Molecular Dynamics for a Seven-Atom Reaction OH + CH₄ → CH₃ + H₂O, J. Chem. Phys., 2013, 138, 221103.
137. Q. Meng, J. Chen and D. H. Zhang, Rate Coefficients of the H + CH₄ → H₂ + CH₃ Reaction from Ring Polymer Molecular Dynamics on a Highly Accurate Potential Energy Surface, J. Chem. Phys., 2015, 143, 101102.
138. Q. Meng, J. Chen and D. H. Zhang, Ring Polymer Molecular Dynamics Fast Computation of Rate Coefficients on Accurate Potential Energy Surfaces in Local Configuration Space: Application to the Abstraction of Hydrogen from Methane, J. Chem. Phys., 2016, 144, 154312.
139. Q. Meng, K. M. Hickson, K. Shao, J.-C. Loison and D. H. Zhang, Theoretical and Experimental Investigations of Rate Coefficients of O(¹D) + CH₄ at Low Temperature, Phys. Chem. Chem. Phys., 2016, 18, 29286–29292.
140. C. Li, Y. Li and B. Jiang, First-Principles Surface Reaction Rates by Ring Polymer Molecular Dynamics and Neural Network Potential: Role of Anharmonicity and Lattice Motion, Chem. Sci., 2023, 14, 5087–5098.

---