ML-MCTDH: Hierarchical Quantum Dynamics

• Multi-Layer MCTDH is a variational method that recursively expands single-particle functions into a layered, tree-based representation, enabling accurate simulations of many-body dynamics.
• The approach leverages a hierarchical tensor network structure to efficiently capture strong correlations across molecular, transport, and spin systems.
• Enhanced numerical techniques, such as projector splitting and gauge optimization, mitigate instabilities and improve performance in high-dimensional, non-adiabatic, and mixed-statistics applications.

Multi-Layer Multiconfiguration Time-Dependent Hartree (ML-MCTDH) is a variational wave-function method for simulating the time-dependent quantum dynamics of many-body systems with a large number of degrees of freedom. It generalizes standard MCTDH by recursively expanding time-dependent single-particle functions into further layers, yielding a hierarchical tensor representation that is naturally described as a tree tensor network state or hierarchical Tucker format. In the limit of suitably large basis sets, the approach is described as numerically exact, and its scope now spans molecular quantum dynamics, non-adiabatic dynamics, transport, impurity problems, ultracold mixtures, and disordered spin models (Vendrell et al., 2010, Larsson, 2023).

1. Formal structure and variational foundation

The defining feature of ML-MCTDH is the replacement of a single-layer time-dependent Hartree expansion by a recursive, layered expansion. A standard form used across applications is

$$|\Psi (t) \rangle = \sum_{j_1} \sum_{j_2} \ldots \sum_{j_p} A_{j_1j_2\ldots j_p}(t) \prod_{\kappa=1}^{p} |\varphi_{j_\kappa}^{(\kappa)} (t) \rangle,$$

where the single-particle functions are themselves expanded recursively into lower-layer functions. In the multilayer language, the wavefunction is therefore associated with a mode tree: the top node represents the full state, internal nodes represent grouped degrees of freedom, and the leaves are primitive basis functions or occupation-number states (Wang et al., 2017, Vendrell et al., 2010).

This recursive structure is the key extension beyond standard MCTDH. In standard MCTDH, the method is effectively a two-layer scheme: top-layer configuration coefficients and time-dependent single-particle functions expanded directly in a primitive basis. ML-MCTDH allows any number of layers and thus adapts the representation depth to the correlation structure of the problem. The same recursive logic underlies both function-based formulations and the tensor-network description, where each node is a tensor and each edge is a virtual index connecting layers (Larsson, 2023).

The equations of motion are obtained from the Dirac-Frenkel variational principle,

$$\langle \delta \Psi(t) | i\frac{\partial}{\partial t} - \hat{H} | \Psi(t)\rangle = 0,$$

which yields coupled equations for the top-layer coefficients and for all lower-layer expansions. At the root, the equation has the Schrödinger-like form

$$i\,\frac{\partial A_I^1}{\partial t} = \sum_J \langle \Phi_I^1 | \hat{H} | \Phi_J^1 \rangle\,A_J^1,$$

while internal-node equations involve reduced density matrices, mean-field operators, and projectors that enforce orthonormality constraints (Larsson, 2023, Vendrell et al., 2010).

A central conceptual point is that the method is not restricted to a single physical representation. It applies to primitive coordinate bases, occupation-number bases, thermofield-doubled spaces, and mixture-adapted species bases. This flexibility is what permits the same variational architecture to appear in chemically motivated wavepacket propagation, Fock-space fermionic transport, and ultracold-gas many-body dynamics (Sasmal et al., 2020, Cao et al., 2017).

2. Statistical formulations and major variants

A major branch of the theory is ML-MCTDH in second quantization representation, ML-MCTDH-SQR. In this formulation, the wavefunction is represented directly in Fock space, and indistinguishable particles are handled algebraically rather than by explicit construction of Slater determinants or permanents. For fermions, the basic algebra is

$$\{ \hat{a}_P, \hat{a}_Q^\dagger \} = \delta_{PQ}, \quad \{ \hat{a}_P^\dagger, \hat{a}_Q^\dagger \} = 0, \quad \{ \hat{a}_P, \hat{a}_Q \} = 0,$$

and the fermionic sign structure is enforced through permutation-sign or Jordan-Wigner-type operator constructions. This makes the multilayer expansion compatible with many-electron transport and impurity models (Wang et al., 2011).

For bosons, the corresponding symmetry-adapted extension is ML-MCTDHB. Here the multilayer ansatz combines species states, bosonic permanent states, and time-dependent single-particle functions. The method is described as a variational numerically exact ab-initio method for studying the quantum dynamics and stationary properties of bosonic systems, and the multilayer structure permits accurate treatment of mixed bosonic systems consisting of arbitrary many species, including multi-dimensional and mixed-dimensional settings (Cao et al., 2013).

A further unification is achieved in ML-MCTDHX, which extends the multilayer formalism to mixtures of bosons and fermions. In the binary-mixture form,

$$|\Psi(t)\rangle = \sum_{i=1}^M \sum_{j=1}^M A_{ij}(t)\, |\psi_i^A(t)\rangle \otimes |\psi_j^B(t)\rangle,$$

with each species basis state expanded over number states built from a small set of time-dependent single-particle functions. In this framework, intra-species correlations are encoded within species basis states, whereas inter-species correlations are encoded at the species layer. The only place where bosonic or fermionic statistics enter is in the structure of the number states and the corresponding (anti-)commutation relations (Cao et al., 2017).

Finite-temperature dynamics has also been incorporated by a thermofield-based formulation. In this construction, a mixed thermal state is mapped to a pure state in an enlarged Hilbert space, and the thermal quasi-particle representation of bosonic thermofield dynamics is used to transfer bosonic many-body MCTDH to finite temperature. The thermofield Schrödinger equation,

$$\frac{\partial}{\partial t}|\psi_\beta(t)\rangle = -\frac{i}{\hbar} (\hat{H}-\tilde{H}) |\psi_\beta(t)\rangle,$$

allows finite-temperature wavefunction propagation without stochastic sampling, and agreement with existing studies on the pyrazine model based on the $\rho$MCTDH method is reported (Fischer et al., 2021).

A related hybridization is ML-MCTDH-oSQR, which combines the tensor contraction scheme of ML-MCTDH-SQR with a time-dependent optimized orbital basis. Unlike ordinary MCTDH-X approaches, this formulation does not possess invariance with respect to time-dependent unitary transformations of the orbital basis; explicit gauge operators therefore enter the equations of motion and affect efficiency and accuracy, especially in imaginary-time propagation (Weike et al., 2019).

3. Equations of motion, gauge structure, and singularity-free propagation

The standard ML-MCTDH equations of motion involve inverse reduced density matrices at internal nodes. When the wavefunction is weakly entangled, these matrices can become singular or nearly singular. In the standard formulation, this leads to numerical instability and motivates regularization procedures. The origin of the difficulty is that standard ML-MCTDH expands the state in terms of orthonormal single-particle functions but generally non-orthonormal single-hole functions, so linear dependence can arise during propagation (Lindoy et al., 2021).

A family of alternative equations of motion addresses this by changing the static gauge condition. Instead of enforcing a semi-unitary coefficient tensor while allowing non-orthonormal single-hole functions, the state is re-expanded in terms of orthonormal single-hole functions and arbitrary coefficient tensors. In this representation, the equations of motion never require the evaluation of the inverse of singular matrices. The projector splitting integrator (PSI) and the invariant equations of motion are identified as special cases of this family, obtained from different dynamic gauge choices (Lindoy et al., 2021).

The practical consequence is substantial. In the multi-layer PSI implementation, forward and backward sweeps move the orthogonality center through the tree, and each working equation is linear and free of density-matrix inversion. For spin-boson models with up to $10^6$ bath modes, PSI is reported to require roughly 3–4 orders of magnitude fewer Hamiltonian evaluations and 2–3 orders of magnitude fewer Hamiltonian applications than standard ML-MCTDH, and 2–3/1–2 orders of magnitude fewer evaluations/applications than approaches that use improved regularization schemes. The same framework was applied to multi-spin-boson models with wavefunctions containing up to $\sim 1.3\times 10^9$ parameters (Lindoy et al., 2021).

Imaginary-time propagation introduces additional gauge considerations. In ML-MCTDH-oSQR, the gauge matrix is hermitian in real time but anti-hermitian in imaginary time. A natural gauge that keeps the one-particle density matrix diagonal is possible, but it introduces higher-body terms and an unfavorable scaling of $\sim f^{2N+2}$ compared to $$\langle \delta \Psi(t) | i\frac{\partial}{\partial t} - \hat{H} | \Psi(t)\rangle = 0,$$0 otherwise. A spectral gauge is therefore introduced and recommended for imaginary-time propagation because it minimizes the spectral width of the effective Hamiltonian while retaining the more favorable scaling (Weike et al., 2019).

Particle-number conservation is another central issue in Fock-space formulations. Since ML-MCTDH-(o)SQR operates in Fock space, all particle-number sectors are formally present. This is usually benign in real time but can destabilize imaginary-time propagation by driving the state toward an incorrect lower-energy sector. A constraint operator of the form

$$\langle \delta \Psi(t) | i\frac{\partial}{\partial t} - \hat{H} | \Psi(t)\rangle = 0,$$1

is proposed to enforce the desired particle number efficiently, and the local particle number associated with a given single-particle function is shown to be preserved under the equations of motion. This leads to a sparse block structure and a refined tensor contraction scheme that explicitly utilizes particle number conservation (Weike et al., 2019).

4. Hamiltonian representations, tensor contractions, and basis engineering

ML-MCTDH is most efficient when the Hamiltonian is available in a structured product form. In the conventional molecular-dynamics setting, this usually means a sum-of-products representation,

$$\langle \delta \Psi(t) | i\frac{\partial}{\partial t} - \hat{H} | \Psi(t)\rangle = 0,$$2

which in tensor-network language can be viewed as a matrix product operator or, on general trees, as a tensor network operator. The recursive evaluation of mean fields and reduced matrices then follows the same tree used for the wavefunction (Larsson, 2023).

Potential-energy fitting is therefore a major part of the ML-MCTDH ecosystem. Traditional Potfit supplies a Tucker-type sum-of-products approximation, but its number of terms grows exponentially with dimensionality. Multi-Layer Potfit (MLPF) replaces this by a hierarchical singular value decomposition and yields a potential in hierarchical tensor format with strict upper bounds for the truncation error,

$$\langle \delta \Psi(t) | i\frac{\partial}{\partial t} - \hat{H} | \Psi(t)\rangle = 0,$$3

A recursive scheme for using the hierarchical-tensor potential within ML-MCTDH is derived, and theoretical estimates together with a computational example show that MLPF can reduce the numerical effort for ML-MCTDH by orders of magnitude compared to Potfit; the effect is already beneficial for systems with just four modes (Otto, 2013).

For reaction dynamics, the same structural principle appears in a broader workflow. The kinetic energy operator is derived as a sum-of-products of single-particle differential operators through a polyspherical approach, while the potential energy surface is expressed in a similar sum-of-products form of single-particle potentials through reconstruction approaches using an existing PES or direct approaches based on a computed database. In that framework, the coordinate frame, the Hamiltonian, and the wavefunction are all designed hierarchically, and the whole construction can be rearranged in the mathematical language of tensor networks or tree tensor networks (Zhang et al., 8 Jul 2025).

In transport theory, basis engineering can be as important as operator factorization. For vibrationally coupled electron transport in single-molecule junctions, a scattering-state representation of the electronic Hamiltonian diagonalizes the non-interacting electronic part and greatly reduces the artificial electron correlation present in the original representation. With an appropriate choice of scattering states, one Hartree product is exact in the non-interacting limit, and steady-state currents can be simulated efficiently over a wide physical parameter space. The best performance in the presence of strong electron-vibration coupling is obtained when the bridge energy is shifted by the reorganization energy,

$$\langle \delta \Psi(t) | i\frac{\partial}{\partial t} - \hat{H} | \Psi(t)\rangle = 0,$$4

This basis choice is explicitly motivated by the suppression of representation-induced, rather than physical, correlation (Wang et al., 2013).

A different route to basis design appears in non-adiabatic dynamics without potential energy surfaces. By combining first-quantized nuclei with second-quantized electrons, all degrees of freedom become distinguishable indices, so Tucker and hierarchical Tucker decompositions can be applied directly to the full nuclear-electronic wavefunction. The formalism circumvents explicit potential energy surfaces and non-adiabatic couplings while remaining compatible with MCTDH and ML-MCTDH (Sasmal et al., 2020).

5. Principal application domains

In molecular quantum dynamics, ML-MCTDH has served as a high-dimensional wavepacket propagator for benchmark and realistic Hamiltonians. A fully general implementation was reported for any number of layers and applied to generalized Henon-Heiles models and to pyrazine. For 6D Henon-Heiles the overhead of ML-MCTDH makes the method slower than MCTDH, but for 18D it starts to be competitive; 1458D simulations of the HH Hamiltonian using a seven layer scheme were also reported. In the 24D pyrazine problem, spectra with all the correct features were obtained in calculations needing a fraction of the time and resources of reference MCTDH calculations, and a highly converged ML-MCTDH result used only $$\langle \delta \Psi(t) | i\frac{\partial}{\partial t} - \hat{H} | \Psi(t)\rangle = 0,$$5 coefficients compared with a best available MCTDH benchmark based on $$\langle \delta \Psi(t) | i\frac{\partial}{\partial t} - \hat{H} | \Psi(t)\rangle = 0,$$6 time-dependent coefficients (Vendrell et al., 2010).

Coupled exciton-vibrational dynamics is another mature domain. For the Fenna-Matthews-Olson complex, ML-MCTDH was used for seven- and eight-site models with 518 and 592 harmonic vibrational modes, respectively. The calculations treated all oscillators on an equal footing, analyzed reduced one-exciton density matrices, and identified distinct features due to competing time scales of vibrational and exciton motion and vibronically-assisted transfer. The same broad class of non-adiabatic dynamics was extended to finite temperature in thermofield ML-MCTDH studies of pyrazine, where electronic population decays, vibrational occupations, and absorption spectra were shown to vary systematically with temperature (Schulze et al., 2016, Fischer et al., 2021).

Quantum transport and impurity problems constitute a second major application area. ML-MCTDH-SQR was used for vibrationally coupled electron transport in single-molecule junctions, where the time-dependent formation of a polaron state produced a pronounced suppression of current corresponding to phonon blockade, and later for correlated transport models including both electron-electron and electron-vibrational interaction. In the nonequilibrium Anderson impurity model at zero temperature, two methodological innovations were introduced for the Kondo regime: correlated initial states obtained by imaginary-time propagation of the full Hamiltonian at zero voltage, and logarithmic discretization of the electronic continuum. These improvements strongly suppressed transient oscillations and enabled analysis of linear and nonlinear conductance in the Kondo regime (Wang et al., 2011, Wang et al., 2013, Wang et al., 2017).

The same transport-oriented logic extends to nonequilibrium dynamical mean-field theory. A time-dependent single-impurity Anderson model produced by a Keldysh-contour mapping can be solved by MCTDH and its multilayer extensions, and the self-consistent two-time impurity Green’s function can be computed by propagating wavefunctions in particle-number sectors that differ by $$\langle \delta \Psi(t) | i\frac{\partial}{\partial t} - \hat{H} | \Psi(t)\rangle = 0,$$7. In this setting, MCTDH outperformed exact diagonalization for large baths still within the reach of exact diagonalization and enabled calculations beyond the system size accessible to exact diagonalization (Balzer et al., 2014).

Ultracold-gas applications use statistics-adapted multilayer formulations. ML-MCTDHB was developed for the non-equilibrium quantum dynamics of multi-species bosonic systems, including tunneling dynamics in one-dimensional double wells, where population imbalances can feature temporal equilibration and bosons of different species can exhibit bunching tendencies. ML-MCTDHX then unified bosonic and fermionic mixtures in a single ab-initio framework, and applications to colliding Bose-Fermi and Fermi-Fermi mixtures showed reflection or transmission processes that are suppressed in the mean-field approximation but dominated by correlations in the many-body treatment (Krönke et al., 2012, Cao et al., 2017).

More recently, ML-MCTDH was applied to quantum spin models with disorder and long-range interactions. Ground states of one- and two-dimensional disordered spin systems were obtained by imaginary-time propagation or improved relaxation, and the hierarchical multilayering was explicitly presented as a way to tackle systems that violate area laws of entanglement entropy. The same formal flexibility that had been exploited in molecular and ultracold physics was thus transferred to disordered spin models (Köhler et al., 2022).

6. Tensor-network interpretation, benchmarking, and current issues

The contemporary understanding of ML-MCTDH is increasingly tensor-network based. The method and DMRG are now routinely described as workhorses built on tensor network states but using different mathematical languages. In this view, ML-MCTDH is a tree tensor network state, standard MCTDH is a Tucker decomposition, and MPS/DMRG correspond to special one-dimensional tensor-network geometries. A major consequence of this reformulation is direct transfer of concepts: optimizing unoccupied single-particle functions in MCTDH corresponds to subspace enrichment in DMRG, and canonicalization or changing the orthogonality center in tensor networks corresponds to changing the root node in ML-MCTDH (Larsson, 2023).

Tree structure is a decisive practical variable. In molecular applications, physically sensible grouping of strongly correlated coordinates is essential; in the 1458D Henon-Heiles calculations, success depended on a seven-layer scheme tailored to the problem structure (Vendrell et al., 2010). In the tensor-network review, tensors of order three in trees were identified as the most optimal scaling for large TTNSs, and a greedy and simulated annealing-based algorithm was proposed for tree optimization (Larsson, 2023). In disordered spin problems, the same issue appears as a choice of logical blockings for one- and two-dimensional geometries, long-range interactions, and disorder (Köhler et al., 2022). This suggests that the mode tree is not a secondary implementation detail but a primary approximation parameter.

Benchmarking against DMRG-style dynamics has clarified a persistent misconception. A recent comparison of TD-DMRG and ML-MCTDH for nonadiabatic exciton dissociation revisited an earlier report of discrepancies up to 60%. Using a unified software framework and systematically increasing bond dimensions, the difference was first reduced to less than 10%, and then to approximately 2% through extrapolation in $$\langle \delta \Psi(t) | i\frac{\partial}{\partial t} - \hat{H} | \Psi(t)\rangle = 0,$$8 and an optimized tensor network structure. The study concluded that both methods converge to numerically exact solutions when bond dimensions are adequately scaled, and that the previously reported discrepancy arose primarily from insufficient bond dimensions rather than from an inherent limitation of either method (Li et al., 30 Aug 2025).

The remaining limitations are structural rather than conceptual. ML-MCTDH can be slower than standard MCTDH in modest dimensionality because of recursive overhead, as seen for 6D Henon-Heiles, and it remains costly when strong entanglement forces large numbers of single-particle functions or large bond dimensions (Vendrell et al., 2010). In spin systems that violate area laws, the method is explicitly presented as promising rather than asymptotically effortless (Köhler et al., 2022). The modern literature therefore treats ML-MCTDH as a systematically improvable hierarchical ansatz whose practical success depends on gauge choice, contraction strategy, conserved-quantity exploitation, Hamiltonian format, and especially tree design (Lindoy et al., 2021, Larsson, 2023).