In a recent paper [H.-D. Meyer and H. Wang, J. Chem. Phys. 148, 124105 (2018)], we have examined the regularization of the equations of motion (EOMs) of the multiconfiguration time-dependent Hartree (MCTDH) approach. We could show that the standard regularization scheme used by almost all researchers in the field is not optimal. The improved regularization allows for larger values of the regularization parameter ϵ, is less sensitive to the actual choice of ϵ, and performs the rotation of initially unoccupied single-particle functions into the “correct” direction in Hilbert space much faster than the old scheme. The latter point increases both the accuracy and efficiency of time propagation for challenging problems. For simple problems, the new scheme requires some additional numerical work as compared with the old scheme, ranging from negligible to almost doubling the total numerical labor. For demanding problems, on the other hand, the additional numerical work of the new scheme is often overcompensated by less steps taken by the integrator. In the present paper, we generalize the new regularization scheme to the multi-layer (ML) extension of MCTDH. Although the principle idea of the new regularization scheme remains unaltered, it was not obvious how the new scheme should be implemented into ML-MCTDH. The ML-MCTDH EOMs are much more complicated than the MCTDH ones, and for optimal numerical performance it was necessary to derive a recursive algorithm for implementing the new regularization scheme.

Over the last two decades, the multiconfiguration time-dependent Hartree (MCTDH) approach1–6 has been very successful. Several problems of molecular quantum dynamics can be solved in full dimensionality only with this method, for example, the computation of the photoexcitation spectrum of pyrazine,7 or other photoexcitation or photoionisation spectra of polyatomic molecules with vibronically coupled electronic states,8–10 the computation of the IR spectrum of the Zundel cation (H5O2+),11–14 and studies on the tunneling splitting in malonaldehyde.15–18 Furthermore, difficult reactive scattering problems have been successfully tackled with MCTDH; see, e.g., the computation of reaction rates and cross sections for the H + CH4 → H2 + CH3 reaction.19–22 Moreover, studies on electron transfer reactions in condensed phases23–26 have been performed by employing high-dimensional models for the solid.

In 2003, Wang and Thoss introduced a multi-layer (ML) extension to the MCTDH algorithm,27 and in 2008 Manthe provided a recursive formulation28 of the rather intricate equations; the latter formulation simplifies the implementation of the cumbersome ML-equations. The ML-MCTDH method27–30 has been shown to be very successful when treating systems with many degrees of freedom (DOFs)—between 10 and 1000 DOFs, say—and when many of these DOFs are only weakly coupled. Examples of successful applications of ML-MCTDH are, e.g., the study of photo-dissociation of a molecule embedded in a complex,31 vibronic dynamics,32 reactive scattering,33,34 exciton vibrational dynamics,35–37 singlet fission,38 surface scattering,39 quantum transport,40–45 and charge transfer processes.46–59 

Despite the great success of MCTDH and ML-MCTDH, there remains a problem with their equations of motion (EOMs). If there are unoccupied single particle functions (SPFs)—and starting from a Hartree product all SPFs, expect the first one in each degree of freedom, are unoccupied—then there are two problems. First, the initially unoccupied SPFs are undefined and they can be chosen at random. Second, the equations of motion (EOMs) are singular and need to be regularized to make them numerically solvable. The first problem was addressed by Manthe60,61 (see also the discussion of Lee and Fischer62 on bosonic systems), who showed how to define optimal initially unoccupied SPFs, where the optimal choice increases the accuracy of the MCTDH propagation.61 The second problem is usually solved by regularizing the density matrix, as the singularities of the density matrix cause the singularities of the EOM. However, we recently have shown63 that much better results are obtained when the A-coefficients (the expansion coefficients) are regularized, which in turn regularizes the density. The new regularization scheme quickly rotates the unoccupied SPFs toward their optimal form and thus lessens the need for starting with optimal SPFs. This could be demonstrated for MCTDH calculations on a series of examples. In the present paper, we extend our approach to ML-MCTDH. For ML-MCTDH calculations, the new regularization scheme is in general more important than for MCTDH ones, as in ML-MCTDH singularities are more frequent, because they may appear at every layer. As ML-MCTDH is a rather intricate formalism, this extension of the regularization is not simply straightforward, although the general idea of the new regularization scheme is identical to the one of our previous work.63 

We make use of the nomenclature introduced by Manthe,28 with slight modifications introduced in Ref. 29, and write an MCTDH wavefunction as

Ψ(q11,,qp11,t)=j1n11jp1np11A1;j1,,jp11(t)×φj1(1;1)(q11,t)φjp1(1;p1)(qp11,t),
(1)

where qκ1 are multi-dimensional mode-combined logical coordinates. There are p1 logical coordinates and nκ1 denotes the number of SPFs of the κth mode. The upper index 1 indicates that the quantities are for the layer number 1, also called the top layer. In the MCTDH method, which is equivalent to a two-layer ML-MCTDH approach,64 the multi-dimensional SPFs are given by the second layer equation

φm1;κ1(qκ1,t)=j1N1κ1jp2;κ1Np2;κ1κ1Am;j1,,jp2;κ12;κ1(t)×χj1(κ1,1)(q12;κ1)χjpκ(κ1,p2;κ1)(qp2;κ12;κ1),
(2)

where χ denotes a primitive basis function and qκ22;κ1 is here the physical coordinate of the κ2th DOF of the node (2; κ1), which has p2;κ1 branches.

Adding additional layers, one has to specify which node of the ML-tree is under consideration. For this, it is convenient to introduce the symbol

z=l;κ1,,κl1 and z1=l1;κ1,,κl2
(3)

and similar for z + 1, etc. Here l denotes the layer and κ1, …, κl−1 denote the modes of the layers above which leads to the node under consideration. Hence, z precisely defines one particular node in the ML-tree.28,29

The ML expansion can now be written in full generality,

φmz1,κl1(qκl1z1)=j1n1zjpznκlzAm;j1,,jpzzκl=1pzφjκlz,κl(qκlz)
(4a)
=JAm;JzΦJz(qκl1z1),
(4b)

which implicitly defines the configurations ΦJz and the composite index J, and where the coordinates are combined as

qκl1z1={q1z,,qpzz}.
(5)

Equation (4) holds for the top layer if one defines q0 as the combination of all DOFs and sets Ψ(q0,t)=φ10(q0,t), and it also holds for the bottom layer if one identifies φz,κl with the primitive function χκl, when z points to a bottom node. Hence Eq. (4) holds for all layers.

The EOMs for the propagation of the SPFs are formally the same for all layers

iφnz,κlt=(1P^z,κl)j,m(ρz,κl)nj1Ĥjmz,κlφmz,κl,
(6)

where P^z,κl=j|φjz,κlφjz,κl| is the projector onto the space spanned by the φjz,κl SPFs, ρz,κl is a density matrix, and Ĥz,κl is a matrix of mean-field operators acting on the φjz,κl functions. Note that this equation is identical to the MCTDH-EOM, except that all quantities carry in addition the node-index z.

The EOMs for the A-vectors follow from Eq. (6) by projecting the SPFs onto their corresponding basis,

An;Jzt=ΦJz|φnz1,κl1t.
(7)

There is no contribution from ΦJz/t because, due to the projector, the SPFs are orthogonal to their time derivatives. However, the top A-vector obeys a different equation

iA1;I1t=JΦI1|Ĥ|ΦJ1A1;J1.
(8)

The density and the mean fields are given by

ρjmz,κl=Ψjz,κlΨmz,κl
(9)

and

Ĥjmz,κl=Ψjz,κlĤΨmz,κl
(10)

and a single-hole function (SHF), Ψjz,κl, is defined by erasing a SPF from the full wavefunction,

Ψ=m=1nκlzΨmz,κlφmz,κl.
(11)

The top-layer SHF has a simple appearance

Ψm1,κ1=Jκ1A1;Jκ1,m1ΦJκ11;κ1,
(12)

where Jκ1 is a composite index similar to J, but it misses the κ1th entry. Similarly, ΦJκ11;κ1 denotes a product of SPFs except the function φj1,κ1, i.e., in general, ΦJκlz,κl=ΦJz/φjz,κl. And A1;Jκ1,m1 is similar to A1;J1 but with the κ1th index replaced by m. (Compare with Ref. 29.) The SHFs of the deeper layers are more complicated and they obey the recursion28 

Ψmz,κl=j=1nκl1z1Ψjz1,κl1Ψ̃m,jz,κl,
(13)

where

Ψ̃m,jz,κl=JκlAj;Jκl,mzΦJκlz,κl.
(14)

The second layer SHF thus reads

Ψm2;κ1,κ2=Jκ1,Jκ2jA1;Jκ1,j1Aj;Jκ2,m2,κ1ΦJκ11;κ1ΦJκ22;κ1,κ2.
(15)

It is convenient to abbreviate the product of A-vectors by B, starting with the top layer

BJ1;κ1,m1;κ1=A1;Jκ1,m1      
(16)

and continuing recursively

BJz,κl,mz,κl=j=1nκl1z1BJz1,κl1,jz1,κl1Aj;Jκl,mz,
(17)

where Jz,κl=(Jκ1,,Jκl) is a composite index, which can be extended recursively as Jz,κl=(Jz1,κl1,Jκl).

We consider BJz,κl,mz,κl as a matrix and, due to the orthonormality of the SPFs, write the densities as

ρz,κl=Bz,κlBz,κl.
(18)

From Eqs. (17) and (18) follows the top-down recursion for the densities,

ρijz,κl=a,b=1nκl1z1ρabz1,κl1JκlAa;Jκl,iz*Ab;Jκl,jz.
(19)

To derive a recursive generation of the mean-fields, we first introduce a sum-of-products (SOP) form of the Hamiltonian

Ĥ=r=1scrĤr=r=1scrκ=1fĥr(prim;κ),
(20)

where ĥr(prim;κ) operates on the κth DOF of the bottom layer. The operator summands in Eq. (20) can be written as

Ĥr=Ĥrz,κlĥrz,κl
(21a)
=Ĥrz,κlκl+1=1pκlzĥrz+1,κl+1,
(21b)

where ĥrz,κl acts on the logical coordinate qκlz, while Ĥrz,κl acts on the corresponding remaining space. As seen from Eq. (21b), ĥrz,κl can as well be written as a direct product of operators ĥrz+1,κl+1, which act on logical coordinates qκl+1z+1. In fact, both operators Ĥrz,κl and ĥrz,κl can be broken down to a simple product of operators which act on one primitive coordinate only [see Eq. (20)].

Using the SOP form, we write the mean field operator as

Ĥjmz,κl=r=1sHrz,κlĥrz,κl,
(22)

where the mean-field matrix H is given by

Hrz,κl=Bz,κlMrz,κlBz,κl
(23)

with

Mr;Jz,κl,Jz,κlz,κl=crΦJκ11;κ1ΦJκlz,κlĤrz,κlΦJκ11;κ1ΦJκlz,κl
(24)

For an efficient implementation of the EOMs, the matrix elements of the operators ĥrz,κl with respect to the SPFs are needed because the SPFs of one layer (z) are expanded in the SPFs of the layer below (z + 1). We follow the discussion of Ref. 29 and define

ĥrz,κlaκl,bκl=φaκlz,κl|ĥrz,κl|φbκlz,κl=φaκlz,κl|κl+1=1pzĥrz+1,κl+1|φbκlz,κl=κl+1=1pz+1Jκl+1i,j=1nκl+1z+1Aakl;Jκl+1,iz+1*ĥrz+1,κl+1i,jAbkl;Jκl+1,jz+1.
(25)

Here the last line defines the recursion of the h-matrix elements, starting from the bottom layer

ĥrlmax;κ1,,κlmaxiκ,jκ=χiκ(κ)|ĥrprim;κ|χjκ(κ),
(26)

where lmax denotes the maximal layer (i.e., lmax = 2 for normal MCTDH) and proceeds bottom-up through the tree structure. In a practical implementation of the EOMs (6), one replaces ĥrz1,κl1 with

ĥrz1,κl1κl=1pκl1z1iκl,jκlφiκlz,κlĥrz,κliκl,jκlφjκlz,κl
(27)

because an SPF φmz1,κl1 is expanded into the SPFs φjκlz,κl.

From Eq. (21b) follows the top-down recursion for Ĥrz,κl,

Ĥrz,κl=Ĥrz1,κl1νκlpz1ĥrz,ν
(28)

and inserting this into Eq. (24) gives the recursion for M,

Mr;Jz,κl,Jz,κlz,κl=Mr;Jz1,κl1,Jz1,κl1z1,κl1ΦJκlz,κl|νκlpz1ĥrz,ν|ΦJκlz,κl=Mr;Jz1,κl1,Jz1,κl1z1,κl1νκlφjνz,κl|ĥrz,ν|φjνz,κl=Mr;Jz1,κl1,Jz1,κl1z1,κl1νκlĥrz,νjν,jν,
(29)

where jν and jν are the indices which appear in Jκl and Jκl, respectively. With this result, we finally arrive at the desired recursion for the mean-field matrix [compare with Eqs. (17) and (23)],

Hr;mnz,κl=a,b=1nκl1z1Hr;abz1,κl1Jκl,JκlAa;Jκl,mz*ΦJκlz,κl|νκlpz1ĥrz,ν|ΦJκlz,κlAb;Jκl,nz=a,b=1nκl1z1Hr;abz1,κl1νκlpz1Jκl,νjν,jνnνzAa;Jκl,ν,m,jνz*×ĥrz,νjν,jνAb;Jκl,ν,n,jνz,
(30)

where the symbol Jκl,ν denotes a composite index similar to J, but with the κlth and νth index missing. And Ab;Jκl,ν,n,jνz is to be interpreted as Ab;Jz, but with the κlth and νth entry of J replaced with n and jν.

In an actual implementation, only the A-tensors of the SPFs are propagated. Using Eqs. (4), (6), (7), (22), and (27), the EOMs are given by

iAn;Jzt=KδJKaAa;JzAa;K*zj,mρz1,κl1n,j1×rHr;j,mz1,κl1Iκl=1pκl1z1ĥrz,κlkκl,iκlAm;Iz,
(31)

where the indices kκl and iκl are the ones which appear in the composite indices K and I, respectively.

1. Formal expression

Using Eqs. (18), (22), and (23), one can write the EOM (6) as

iφz,κlt=(1P^z,κl)Bz,κlBz,κl1Bz,κl×r=1sMrz,κlBz,κlĥrz,κlφz,κl.
(32)

To regularize this EOM, one may proceed similarly as done in Ref. 63 and perform a singular value decomposition (SVD) of the matrix B,

Bz,κl=Uz,κlλz,κlVz,κl,
(33)

where V is a unitary (nκlz×nκlz) matrix, λ is a (nκlz×nκlz) diagonal matrix with real non-negative singular values, and U is a (n¯κlz×nκlz) matrix with orthonormal column vectors. Here

n¯κlz=j=0l1νκljpzjnνzj
(34)

is the size of the index space {Jz,κl}.

Inserting Eq. (33) into Eq. (18), we notice

ρz,κl=Vz,κlλz,κl2Vz,κl,
(35)

which shows that the singular values are the square roots of the eigenvalues of ρ, i.e., of the so-called natural populations.

Proceeding as discussed in Ref. 63, we insert the SVD (33) into Eq. (32) (except for the very last matrix B), let several terms cancel each other, and obtain

iφz,κlt=(1P^z,κl)Vz,κlλz,κl1Uz,κl×r=1sMrz,κlBz,κlĥrz,κlφz,κl.
(36)

To arrive at a singularity-free EOM, one has to regularize the singular values λ. For this, we may replace λjz,κl with λj,regz,κl=max(λjz,κl,ϵ1/2) [or use a variant of replacing λjz,κl with λj,regz,κl=λjz,κl+ϵ1/2exp(λjz,κl/ϵ1/2)], where ϵ denotes a small positive regularization parameter.63 The conventional approach, where only the density matrix is regularized, is equivalent63 to replacing (λz,κl)1 in Eq. (36) with (λregz,κl)2λz,κl. The jth diagonal element of the latter product vanishes if the jth singular value vanishes. Adopting the new regularization scheme, however, the jth element of (λregz,κl)1 assumes in this case the large value ϵ−1/2. If there is an unoccupied SPF (more precisely an unoccupied natural orbital), then the old scheme will not move this SPF until it has obtained an occupation of the order of ϵ. The new scheme, on the other hand, will immediately start to rotate the unoccupied SPF with a high speed into its “optimal” direction in Hilbert space.63 

However, the procedure (33) and (36) is not practical for two reasons. First, n¯κlz, the size of the index space of Jz,κl, may become very large, in particular, for the lower layers. An SVD of the matrix B in general will become very expensive, if not impossible. Second, as discussed in Sec. II B, the matrices B and M are not explicitly built in practice, but the densities and the mean-field matrices are built recursively starting from their simple top-layer forms.

2. Practical implementation

To derive a practical recursive evaluation of the SVD, we insert the diagonal form of the density, Eq. (35), into the recursion formula for the density, Eq. (19),

k=1nκlzVikz,κlλkz,κl2Vjkz,κl*=Jκla,b,c=1nκl1z1Aa;Jκl,iz*      Vacz1,κl1λcz1,κl12×Vbcz1,κl1*Ab;Jκl,jz=JκlaÃa;Jκl,iz,κl*Ãa;Jκl,jz,κl=m=1nκlzvimz,κlμmz,κl2vjmz,κl*,
(37)

where in the first step we have introduced the definition

Ãa;Jκl,jz,κl=c=1nκl1z1λaz1,κl1Vcaz1,κl1*Ac;Jκl,jz
(38)

and in the second step made use of a singular value decomposition of Ã, where (a;Jκl) is considered as a composite index,

Ãa;Jκl,jz,κl=k=1nκlzua;Jκl,kz,κlμkz,κlvjkz,κl*.
(39)

Comparing the right- and left-hand side of Eq. (37) shows that

vjkz,κl=Vjkz,κl,μkz,κl=λkz,κl
(40)

holds.

This is a quite remarkable result, as it shows how to compute the right-hand side SVD matrix V and the singular values λ of the matrix B without performing an SVD of B, but by performing an SVD of the much smaller matrix à instead. The matrix Ãz,κl has the dimension (ñκlz×nκlz) with ñκlz=nκl1z1νκlpznνz. This number is usually not prohibitively large and is much smaller than the corresponding B index n¯κlz, because in an ML calculation the numbers of logical modes, pz, are usually small. (For a binary tree, pz = 2 holds throughout.) Finally we note that when building Ã, Eq. (38) in practice, we employ the regularized singular values rather than the original ones. This improves the stability of the SVD of Ã, Eq. (39).

What is still missing is a recursion for U. For this, we combine the SVD of B, Eq. (33), with the recursion of B, Eq. (17),

BJz,κl,mz,κl=k=1nκlzUJz,κl,kz,κlλkz,κlVmkz,κl*=a,b=1nκl1z1UJz1,κl1,az1,κl1λaz1,κl1Vbaz1,κl1*Ab;Jκl,mz=aUJz1,κl1,az1,κl1Ãa;Jκl,mz,κl=a=1nκl1z1k=1nκlzUJz1,κl1,az1,κl1ua;Jκl,kz,κlλkz,κlVmkz,κl*.
(41)

Comparing the first and the last line shows the recursion

UJz,κl,kz,κl=a=1nκl1z1UJz1,κl1,az1,κl1ua;Jκl,kz,κl,
(42)

where, as before, Jz,κl=(Jz1,κl1,Jκl).

Of course, one does not want to generate the huge matrix Uz,κl explicitly; one rather computes the modified mean-field [compare with Eq. (23)]

H̃rz,κl=Uz,κlMrz,κlBz,κl
(43)

recursively. Proceeding similarly as when deriving Eq. (30), and making use of Eqs. (17), (29), (42), and (43), one arrives at the desired recursion for the modified mean-field matrix

H̃r;mnz,κl=a,b=1nκl1z1H̃r;abz1,κl1νκlpz1Jκl,ν×jν,jνnνzua;Jκl,ν,m,jνz,κl*ĥrz,νjν,jνAb;Jκl,ν,n,jνz,
(44)

where, as before, the symbol Jκl,ν denotes a composite index similar to J, but with the κlth and νth index missing. And tensor Ab;Jκl,ν,n,jνz is to be interpreted as Ab;Jz, but with the κlth and νth entry of J replaced with n and jν. And a similar notational convention is adopted for tensor u. Finally, the EOM for the coefficients [compare with Eq. (31)] reads

iAn;Jzt=KδJKaAa;JzAa;K*zj,mVn,jz1,κl1λj,regz1,κl11×rH̃r;j,mz1,κl1Iκl=1pκl1z1ĥrz,κlkκl,iκlAm;Iz.
(45)

The working equations for the recursions are (25), (38)–(40), and (44). First, the h-matrix elements are built bottom-up (25) and then the top-level mean-field matrices are built. These are given by the usual MCTDH formulas and their generation requires a SVD of the top-layer coefficient tensor, A1;Jκ1,m1. Then the top-down recursions for the mean-fields (44) are performed, which requires the building of à and its singular-value decomposition (38)–(40). After that, the EOM for the coefficients, Eqs. (8) and (45), is to be solved.

Similar as in the previous work,63 we test the new regularization scheme on the spin-boson model.65,66 This is a system-bath type Hamiltonian that has been widely used in the context of electron transfer theories. The model has two electronic states (|ϕ1⟩ and |ϕ2⟩) linearly coupled to a bath of harmonic oscillators. Using mass-weighted coordinates, the Hamiltonian reads

H=Dσz+Δσx+12j=1N(pj2+ωj2qj2)+σzj=1Ncjqj,
(46)

where σx and σz are Pauli matrices

σx=ϕ1ϕ2|+|ϕ2ϕ1,
(47a)
σz=ϕ1ϕ1||ϕ2ϕ2.
(47b)

The properties of the bath that influence the dynamics of the two-state subsystem are specified by the spectral density function65,66

J(ω)=π2jcj2ωjδ(ωωj).
(48)

In the test examples below, we use an Ohmic (linear) spectral density with an exponential cutoff

J(ω)=π2αωeω/ωc,
(49)

where α is the dimensionless Kondo parameter that characterizes the system-bath coupling strength and ωc is the cutoff frequency of the bath. The continuous bath spectral density of Eq. (49) can be discretized to the form of Eq. (48) via the relation24,27

cj2=2πωjJ(ωj)ρ(ωj),
(50)

in which the density of frequencies ρ(ω) is defined from the integral relation

0ωjdωρ(ω)=j,j=1,,N.
(51a)

In this work, ρ(ω) is chosen as

ρ(ω)=N+1ωceω/ωc,
(51b)

where ωN+1 = is removed from the simulation. The number of bath modes N in the discretization scheme is a convergence parameter and needs to be sufficiently large to represent the continuum within the time scale of interest. Unlike in our previous tests using MCTDH,63 ML-MCTDH is a much more powerful method that can easily converge N. The examples below employ a few hundred to a few thousand bath modes and treat physical regimes that have intermediate to very strong coupling strength. The number of layers ranges from 3 to 6 in the simulation. These examples represent challenging cases of realistic ML-MCTDH simulations.

In the examples below, we choose a simple initial state with the bath starting at its ground state. The dynamical property examined here is the time-dependent population (difference) of the two-level subsystem

P(t)σz(t)=Ψ(t)σzΨ(t).
(52)

We first consider an example with parameters D/Δ = 0, ωc/Δ = 25, and an intermediate coupling strength, α = 0.5. Figure 1 illustrates convergence of the result with respect to the regularization parameter ϵ for a bath of 250 modes. A three-layer ML-MCTDH was used in the simulation. The top layer has five single particle groups (logical coordinates) with six SPFs per group and a mode with two SPFs for the electronic DOF, resulting in a total of 65 × 2 = 15 552 configurations. For the second layer, there are four single particle groups per each top layer group and five SPFs are used per each second layer group. For the bottom layer, up to 13 modes are combined through use mode combination3 with the primitive basis functions (three per DOF) adiabatically contracted.23 Among the six SPFs for each top layer group and five SPFs for each second layer group, the first one represents the initial state. The remaining SPFs are randomly selected and orthonormalized against the first SPF. A variable stepsize predictor-corrector ordinary differential equation (ODE) solver was used in the time integration. Both the matrix diagonalization (of the reduced density matrix) in the standard scheme and the SVD in the new scheme are handled by the corresponding LAPACK (Linear Algebra Package) subroutines.

FIG. 1.

Time dependence of P(t) = ⟨σz(t)⟩ for D/Δ = 0, α = 0.5, and ωc/Δ = 25. A bath of 250 modes is used in a three-layer ML-MCTDH simulation. Convergence of P(t) with respect to the regularization parameter ϵ using the standard regularization scheme is shown by the various dashed, dotted, and dashed-dotted lines (see the inset). The ϵ = 10−8 curve lies on top of the ϵ = 10−10 one.

FIG. 1.

Time dependence of P(t) = ⟨σz(t)⟩ for D/Δ = 0, α = 0.5, and ωc/Δ = 25. A bath of 250 modes is used in a three-layer ML-MCTDH simulation. Convergence of P(t) with respect to the regularization parameter ϵ using the standard regularization scheme is shown by the various dashed, dotted, and dashed-dotted lines (see the inset). The ϵ = 10−8 curve lies on top of the ϵ = 10−10 one.

Close modal

From Fig. 1, it looks like that for the standard regularization of the reduced density matrices, ϵ ≤ 10−8 is required to converge the result. A more quantitative measure of the convergence can be obtained by introducing a relative cumulative deviation (hereafter simply referred to as “deviation”), defined as

deviation=1PmaxPmin1τ0τ|P(t)Pref(t)|dt,
(53)

where τ is the time scale of the simulation, Pref(t) is the reference result obtained by using the smallest ϵ in the new regularization scheme, and Pmax and Pmin are the maximum and minimum values of Pref(t). The error tolerance of the predictor-corrector ODE solver was set such that there is ∼10−5 roundoff error in the deviation. Any deviation below, say 5 × 10−6, can be considered as allowable numerical noise. For the numerical convergence of P(t), we set the criterion that the deviation be less than ∼10−4, which is stringent enough as compared with other convergence parameters and is still sufficiently larger than the roundoff error.

In Table I, the convergence behavior using the new regularization is compared with the one using the standard regularization scheme. As discussed in our previous work,63 reducing ϵ may thus have two opposite effects: (1) a smaller ϵ results in a stiffer set of EOMs and thus will increase the computational cost; (2) by reducing ϵ, the correct SPF subspace is more quickly reached and the integrator will work at a true physical time scale thereafter and thus the overall effect is that the computational cost is reduced. Table I examines these two opposite effects by displaying the number of calls to the ML-MCTDH derivative subroutine in which time derivatives of the A-vectors of all layers are evaluated. The new scheme quickly rotates the unoccupied SPFs to the “correct” direction and is thus more effective, under which convergence is reached when ϵ ≤ 10−8. The standard scheme requires ϵ ≤ 10−9 to achieve the same convergence. Results with a larger ϵ in the new regularization scheme do not deviate as much as those obtained in the standard scheme. This shows that the new scheme is more robust.

TABLE I.

Number of calls to the ML-MCTDH derivative subroutine and deviation versus the regularization parameter ϵ. The physical parameters are given in Fig. 1. A three-layer ML-MCTDH is used to treat a bath of 250 modes.

Standard regularizationNew regularization
ϵNumber of callsDeviationNumber of callsDeviation
10−4 18 238 2.4 × 10−2 17 744 6.1 × 10−3 
10−6 17 895 1.2 × 10−2 20 840 2.8 × 10−3 
10−7 21 411 4.1 × 10−3 21 314 1.8 × 10−4 
10−8 17 265 4.3 × 10−4 16 541 1.4 × 10−4 
10−9 15 730 1.6 × 10−4 16 198 7.5 × 10−5 
10−10 15 798 9.0 × 10−5 15 758 2.9 × 10−5 
10−11 12 588 3.7 × 10−5 12 555 1.2 × 10−5 
10−12 10 483 1.5 × 10−5 10 282 4.9 × 10−6 
10−13 10 885 8.4 × 10−6 10 259 2.7 × 10−6 
10−14 19 356 4.2 × 10−6 10 502 1.8 × 10−6 
10−15 44 565 2.1 × 10−6 10 782 1.2 × 10−6 
10−16 191 629 2.7 × 10−6 12 583 9.3 × 10−7 
10−18 … … 13 305 2.4 × 10−6 
10−20 … … 17 841 Reference 
Standard regularizationNew regularization
ϵNumber of callsDeviationNumber of callsDeviation
10−4 18 238 2.4 × 10−2 17 744 6.1 × 10−3 
10−6 17 895 1.2 × 10−2 20 840 2.8 × 10−3 
10−7 21 411 4.1 × 10−3 21 314 1.8 × 10−4 
10−8 17 265 4.3 × 10−4 16 541 1.4 × 10−4 
10−9 15 730 1.6 × 10−4 16 198 7.5 × 10−5 
10−10 15 798 9.0 × 10−5 15 758 2.9 × 10−5 
10−11 12 588 3.7 × 10−5 12 555 1.2 × 10−5 
10−12 10 483 1.5 × 10−5 10 282 4.9 × 10−6 
10−13 10 885 8.4 × 10−6 10 259 2.7 × 10−6 
10−14 19 356 4.2 × 10−6 10 502 1.8 × 10−6 
10−15 44 565 2.1 × 10−6 10 782 1.2 × 10−6 
10−16 191 629 2.7 × 10−6 12 583 9.3 × 10−7 
10−18 … … 13 305 2.4 × 10−6 
10−20 … … 17 841 Reference 

Table I illustrates the similarity between the ML-MCTDH and MCTDH regularizations: a quick rotation (small ϵ) of the SPF subspaces is necessary in the beginning of the integration to increase both the accuracy and efficiency. Once ϵ is small enough, the number of calls to the derivative subroutine remains approximately the same (and in fact slightly goes down as for the MCTDH case) within a certain range. As ϵ becomes very small (ϵ = 10−16), the ML-MCTDH EOMs become very stiff under the standard regularization scheme and the computational cost increases by an order of magnitude. By contrast, under the new regularization scheme, the cost remains approximately the same over a large range of ϵ. This suggests robustness of the algorithm in terms of both accuracy and efficiency. Again, the difference lies on the fact that though ϵ formally represents the same parameter in the two schemes, it is used differently in the numerical implementation: in the standard scheme, ϵ−1 is the effective regularization parameter whereas in the new scheme it is ϵ−1/2. An even smaller ϵ = 10−24 may still be handled in the new regularization scheme.

We now consider a longer time dynamics than Fig. 1, which also requires a larger number of modes. Figure 2 illustrates convergence of P(t) for twice as long as in Fig. 1 and for a bath of 1000 modes. A four-layer ML-MCTDH was used in the simulation. The top layer still has five single particle groups, but eight SPFs are used per each group, resulting in a total of 85 × 2 = 65 536 configurations. For the other layers, there are up to four lower single particle groups per each upper layer group and five SPFs are used per each group. Mode combination and adiabatic basis contraction are again used for each bottom layer group, which holds up to 13 modes.

FIG. 2.

Time dependence of P(t) = ⟨σz(t)⟩ for D/Δ = 0, α = 0.5, and ωc/Δ = 25. A bath of 1000 modes is used in a four-layer ML-MCTDH simulation. Convergence of P(t) with respect to the regularization parameter ϵ using the standard regularization scheme is shown by the various dashed, dotted, and dashed-dotted lines (see the inset). Note that the ϵ = 10−6 curve is completely off and that the ϵ = 10−10 and ϵ = 10−12 curves are on top of each other.

FIG. 2.

Time dependence of P(t) = ⟨σz(t)⟩ for D/Δ = 0, α = 0.5, and ωc/Δ = 25. A bath of 1000 modes is used in a four-layer ML-MCTDH simulation. Convergence of P(t) with respect to the regularization parameter ϵ using the standard regularization scheme is shown by the various dashed, dotted, and dashed-dotted lines (see the inset). Note that the ϵ = 10−6 curve is completely off and that the ϵ = 10−10 and ϵ = 10−12 curves are on top of each other.

Close modal

Since numerical singularities occur at all layers in ML-MCTDH, convergence requires smaller ϵ as the number of layers and/or SPFs increases. Figure 2 shows that under the standard regularization scheme, the simulation with ϵ = 10−6 predicts a completely erroneous dynamics—a coherent motion of P(t) near its initial value—that contradicts the converged result even qualitatively. The simulation with ϵ = 10−7 also predicts an incorrect initial delay. A more detailed analysis in Table II shows that a large ϵ also requires more integration steps due to its wrong dynamics.

TABLE II.

Number of calls to the ML-MCTDH derivative subroutine and deviation versus the regularization parameter ϵ. The physical parameters are given in Fig. 2. A four-layer ML-MCTDH is used to treat a bath of 1000 modes.

Standard regularizationNew regularization
ϵNumber of callsDeviationNumber of callsDeviation
10−6 132 677 6.6 × 10−1 63 404 5.6 × 10−2 
10−7 50 285 7.2 × 10−2 49 341 1.2 × 10−2 
10−8 45 403 1.7 × 10−2 49 635 1.8 × 10−3 
10−9 41 964 2.7 × 10−3 47 184 2.2 × 10−4 
10−10 31 844 3.1 × 10−4 31 747 7.5 × 10−5 
10−11 28 998 8.1 × 10−5 30 976 3.7 × 10−5 
10−12 29 100 4.0 × 10−5 29 056 1.5 × 10−5 
10−13 29 422 1.4 × 10−5 28 965 5.2 × 10−6 
10−14 33 765 9.2 × 10−6 29 205 4.4 × 10−6 
10−15 69 370 7.3 × 10−6 30 641 3.8 × 10−6 
10−16 277 915 7.1 × 10−6 29 996 4.4 × 10−6 
10−18 … … 31 292 3.6 × 10−6 
10−20 … … 35 347 4.9 × 10−6 
10−22 … … 40 092 Reference 
Standard regularizationNew regularization
ϵNumber of callsDeviationNumber of callsDeviation
10−6 132 677 6.6 × 10−1 63 404 5.6 × 10−2 
10−7 50 285 7.2 × 10−2 49 341 1.2 × 10−2 
10−8 45 403 1.7 × 10−2 49 635 1.8 × 10−3 
10−9 41 964 2.7 × 10−3 47 184 2.2 × 10−4 
10−10 31 844 3.1 × 10−4 31 747 7.5 × 10−5 
10−11 28 998 8.1 × 10−5 30 976 3.7 × 10−5 
10−12 29 100 4.0 × 10−5 29 056 1.5 × 10−5 
10−13 29 422 1.4 × 10−5 28 965 5.2 × 10−6 
10−14 33 765 9.2 × 10−6 29 205 4.4 × 10−6 
10−15 69 370 7.3 × 10−6 30 641 3.8 × 10−6 
10−16 277 915 7.1 × 10−6 29 996 4.4 × 10−6 
10−18 … … 31 292 3.6 × 10−6 
10−20 … … 35 347 4.9 × 10−6 
10−22 … … 40 092 Reference 

Table II shows that the standard regularization scheme requires ϵ ≤ 10−11 to achieve convergence, whereas the new scheme requires ϵ ≤ 10−10. Results with a larger ϵ in the new regularization scheme do not deviate as much as those obtained in the standard scheme. Furthermore, the computational cost under the new scheme remains almost the same spanning a large range of the regularization parameter, ϵ = 10−20–10−10. All this suggests again the stability of the new regularization scheme.

Examples above are typical applications of ML-MCTDH. Due to the fact that the new regularization scheme performs SVD on the matrix unfoldings of the A-vectors of all layers, instead of diagonalizing smaller density matrices in the standard scheme, it is expected that the new scheme is slightly more expensive if the numbers of calls to the derivative subroutines are the same under the two schemes. Indeed we find that the new scheme takes 25% more time than the standard scheme, e.g., for ϵ = 10−10 in the example above. However, for ϵ ≤ 10−15, the new scheme is clearly more efficient than the standard scheme. As shown below, only the new scheme is practical for investigating difficult regimes.

We now consider challenging applications. Figure 3 shows such an example for D/Δ = 0, α = 2, and ωc/Δ = 25. Dynamics of P(t) exhibit localization after a short transient relaxation. Convergence is shown with respect to the regularization parameter ϵ for a four-layer ML-MCTDH simulation on a model bath of 500 modes. The top layer has four single particle groups, with eight SPFs per each group. For the other layers, there are up to four lower single particle groups per each upper layer group and five SPFs are used per each group. Mode combination and adiabatic basis contraction are used for each bottom layer group. Due to the stronger coupling strength, four primitive basis functions are needed for each mode. Hence each bottom layer group contains only up to 8 modes.

FIG. 3.

Time dependence of P(t) = ⟨σz(t)⟩ for D/Δ = 0, α = 2, and ωc/Δ = 25. A bath of 500 modes is used in a four-layer ML-MCTDH simulation. Convergence of P(t) with respect to the regularization parameter ϵ using the standard regularization scheme is shown by the various dashed, dotted, and dashed-dotted lines (see the inset). Note that the ϵ = 10−8 curve is qualitatively wrong. The ϵ = 10−16 and ϵ = 10−17 curves are on top of each other.

FIG. 3.

Time dependence of P(t) = ⟨σz(t)⟩ for D/Δ = 0, α = 2, and ωc/Δ = 25. A bath of 500 modes is used in a four-layer ML-MCTDH simulation. Convergence of P(t) with respect to the regularization parameter ϵ using the standard regularization scheme is shown by the various dashed, dotted, and dashed-dotted lines (see the inset). Note that the ϵ = 10−8 curve is qualitatively wrong. The ϵ = 10−16 and ϵ = 10−17 curves are on top of each other.

Close modal

As shown in Fig. 3 and Table III, convergence of P(t) is slow, and sometimes even oscillatory, versus ϵ. To achieve deviation less than 10−4 under the standard scheme, ϵ ≤ 10−16 is required. Such a small ϵ results in more than an order of magnitude more computational cost in the standard scheme. On the contrary, not only does the new scheme give almost the same convergence with ϵ ≤ 10−14, but it also requires much less computational effort. Even pushing to a very small ϵ = 10−22, the cost is only three times as expensive.

TABLE III.

Number of calls to the ML-MCTDH derivative subroutine and deviation versus the regularization parameter ϵ. The physical parameters are given in Fig. 3. A four-layer ML-MCTDH is used to treat a bath of 500 modes.

Standard regularizationNew regularization
ϵNumber of callsDeviationNumber of callsDeviation
10−8 7 166 3.4 × 10−1 7 168 3.7 × 10−2 
10−9 5 713 1.3 × 10−2 6 181 3.1 × 10−2 
10−10 5 658 4.0 × 10−2 5 943 1.3 × 10−2 
10−11 5 673 1.9 × 10−2 6 205 4.6 × 10−3 
10−12 5 644 7.0 × 10−3 6 274 1.5 × 10−3 
10−13 6 370 2.4 × 10−3 6 065 4.9 × 10−4 
10−14 12 252 7.8 × 10−4 5 541 1.6 × 10−4 
10−15 43 570 2.7 × 10−4 7 560 4.5 × 10−5 
10−16 217 340 1.0 × 10−4 6 865 1.8 × 10−5 
10−17 1 153 089 5.0 × 10−5 8 200 2.4 × 10−6 
10−18 … … 9 867 3.6 × 10−6 
10−19 … … 10 943 8.3 × 10−6 
10−20 … … 13 284 4.8 × 10−6 
10−21 … … 16 019 3.2 × 10−6 
10−22 … … 19 501 Reference 
Standard regularizationNew regularization
ϵNumber of callsDeviationNumber of callsDeviation
10−8 7 166 3.4 × 10−1 7 168 3.7 × 10−2 
10−9 5 713 1.3 × 10−2 6 181 3.1 × 10−2 
10−10 5 658 4.0 × 10−2 5 943 1.3 × 10−2 
10−11 5 673 1.9 × 10−2 6 205 4.6 × 10−3 
10−12 5 644 7.0 × 10−3 6 274 1.5 × 10−3 
10−13 6 370 2.4 × 10−3 6 065 4.9 × 10−4 
10−14 12 252 7.8 × 10−4 5 541 1.6 × 10−4 
10−15 43 570 2.7 × 10−4 7 560 4.5 × 10−5 
10−16 217 340 1.0 × 10−4 6 865 1.8 × 10−5 
10−17 1 153 089 5.0 × 10−5 8 200 2.4 × 10−6 
10−18 … … 9 867 3.6 × 10−6 
10−19 … … 10 943 8.3 × 10−6 
10−20 … … 13 284 4.8 × 10−6 
10−21 … … 16 019 3.2 × 10−6 
10−22 … … 19 501 Reference 

ML-MCTDH simulations become more challenging for larger systems. Table IV analyzes the result with the same physical parameters as in Table III, but for a bath of 2000 modes. A five-layer ML-MCTDH is employed here. Under the standard scheme, convergence can be barely reached with an ϵ as small as 10−18. More severely, the computational effort increases here by more than two orders of magnitude. Things were quite different when the new regularization scheme was used. Convergence can be achieved within a broad range of ϵ = 10−24–10−16 and with very reasonable computational effort.

TABLE IV.

Number of calls to the ML-MCTDH derivative subroutine and deviation versus the regularization parameter ϵ. The physical parameters are given in Fig. 3. A five-layer ML-MCTDH is used to treat a bath of 2000 modes.

Standard regularizationNew regularization
ϵNumber of callsDeviationNumber of callsDeviation
10−8 55 007 5.1 × 10−1 16 298 4.9 × 10−1 
10−9 42 682 5.0 × 10−1 14 523 1.0 × 10−1 
10−10 15 014 3.0 × 10−1 11 534 4.3 × 10−2 
10−11 13 456 1.7 × 10−2 9 669 3.0 × 10−2 
10−12 12 112 3.8 × 10−2 8 484 1.2 × 10−2 
10−13 12 124 1.9 × 10−2 8 234 4.4 × 10−3 
10−14 15 801 7.1 × 10−3 8 939 1.5 × 10−3 
10−15 30 988 2.5 × 10−3 9 985 4.8 × 10−4 
10−16 149 421 8.4 × 10−4 10 711 1.5 × 10−4 
10−17 635 887 2.7 × 10−4 9 616 3.5 × 10−5 
10−18 2 279 470 1.9 × 10−4 9 635 1.1 × 10−5 
10−20 … … 14 895 2.2 × 10−5 
10−22 … … 20 842 5.0 × 10−6 
10−24 … … 32 983 Reference 
Standard regularizationNew regularization
ϵNumber of callsDeviationNumber of callsDeviation
10−8 55 007 5.1 × 10−1 16 298 4.9 × 10−1 
10−9 42 682 5.0 × 10−1 14 523 1.0 × 10−1 
10−10 15 014 3.0 × 10−1 11 534 4.3 × 10−2 
10−11 13 456 1.7 × 10−2 9 669 3.0 × 10−2 
10−12 12 112 3.8 × 10−2 8 484 1.2 × 10−2 
10−13 12 124 1.9 × 10−2 8 234 4.4 × 10−3 
10−14 15 801 7.1 × 10−3 8 939 1.5 × 10−3 
10−15 30 988 2.5 × 10−3 9 985 4.8 × 10−4 
10−16 149 421 8.4 × 10−4 10 711 1.5 × 10−4 
10−17 635 887 2.7 × 10−4 9 616 3.5 × 10−5 
10−18 2 279 470 1.9 × 10−4 9 635 1.1 × 10−5 
10−20 … … 14 895 2.2 × 10−5 
10−22 … … 20 842 5.0 × 10−6 
10−24 … … 32 983 Reference 

For an even larger bath of 5000 modes, Table V shows that the standard scheme is unable to achieve convergence in a practical sense (because of the two orders of magnitude more computational effort): for ϵ = 10−17, the simulation under the new regularization scheme takes less than 4 h (14 438 calls to the ML-MCTDH derivatives) and almost reaches convergence (deviation ∼10−4), whereas under the standard scheme, ϵ = 10−18 is required to achieve similar convergence and the simulation takes a month (3 046 997 calls to the ML-MCTDH derivatives). In such a situation, the new regularization scheme is the only viable choice.

TABLE V.

Number of calls to the ML-MCTDH derivative subroutine and deviation versus the regularization parameter ϵ. The physical parameters are given in Fig. 3. A six-layer ML-MCTDH is used to treat a bath of 5000 modes.

Standard regularizationNew regularization
ϵNumber of callsDeviationNumber of callsDeviation
10−8 77 628 5.1 × 10−1 80 769 5.0 × 10−1 
10−9 88 948 5.1 × 10−1 28 753 4.8 × 10−1 
10−10 89 720 4.8 × 10−1 25 155 7.7 × 10−2 
10−11 23 541 2.0 × 10−1 20 498 4.4 × 10−2 
10−12 23 581 3.7 × 10−2 14 441 2.8 × 10−2 
10−13 16 476 3.5 × 10−2 11 120 1.1 × 10−2 
10−14 17 973 1.5 × 10−2 9 946 3.8 × 10−3 
10−15 38 166 5.4 × 10−3 8 907 1.3 × 10−3 
10−16 176 864 1.8 × 10−3 10 738 4.2 × 10−4 
10−17 664 745 6.2 × 10−4 14 438 1.2 × 10−4 
10−18 3 046 997 1.3 × 10−4 13 250 3.0 × 10−5 
10−20 … … 13 698 1.5 × 10−5 
10−22 … … 19 282 7.8 × 10−6 
10−24 … … 41 213 Reference 
Standard regularizationNew regularization
ϵNumber of callsDeviationNumber of callsDeviation
10−8 77 628 5.1 × 10−1 80 769 5.0 × 10−1 
10−9 88 948 5.1 × 10−1 28 753 4.8 × 10−1 
10−10 89 720 4.8 × 10−1 25 155 7.7 × 10−2 
10−11 23 541 2.0 × 10−1 20 498 4.4 × 10−2 
10−12 23 581 3.7 × 10−2 14 441 2.8 × 10−2 
10−13 16 476 3.5 × 10−2 11 120 1.1 × 10−2 
10−14 17 973 1.5 × 10−2 9 946 3.8 × 10−3 
10−15 38 166 5.4 × 10−3 8 907 1.3 × 10−3 
10−16 176 864 1.8 × 10−3 10 738 4.2 × 10−4 
10−17 664 745 6.2 × 10−4 14 438 1.2 × 10−4 
10−18 3 046 997 1.3 × 10−4 13 250 3.0 × 10−5 
10−20 … … 13 698 1.5 × 10−5 
10−22 … … 19 282 7.8 × 10−6 
10−24 … … 41 213 Reference 

Finally we consider a variant of the spin-boson model, the polaron-transformed Hamiltonian of (46),

ĤPol=T^+ĤT^=Dσz+12j=1N(p^j2+ωj2q^j2)+Δϕ1ϕ2j=1Ne2ip^jcjωj2+H.c.,
(54)

where T^ is a translational operator for different electronic states

T^=j=1Neip^jcjωj2ϕ1ϕ1+j=1Neip^jcjωj2ϕ2ϕ2.
(55)

This is a pathological example and has been shown earlier by us63 that the standard regularization scheme completely fails in MCTDH for a bath of 40 modes and the new regularization scheme needs be used to obtain the correct result. The phase in the translational operator (55) requires a quick rotation of the single particle space to capture the correct dynamics. As shown in (55), this effect is more pronounced for a stronger coupling strength (larger cj’s) and/or a larger bath (more terms in the product and thus the exponential summation).

Here we consider a similar parameter set as in our previous MCTDH work,63 but with a larger bath of 512 modes. A four-layer ML-MCTDH was used in the simulation. The top layer contains five single particle groups, with eight SPFs used per each group, resulting in a total of 85 × 2 = 65 536 configurations. For other layers, there are 3–4 lower single particle groups per each upper layer group and eight SPFs are used per each group. Mode combination and adiabatic basis contract are again used for each bottom layer group, which holds up to 9 modes.

Under the standard scheme, P(t) remains at its initial value of 1 as time evolves even with the smallest possible ϵ that does not crash the ODE solver, indicating the inability of rotating the single particle space to the “correct” direction. On the contrary, the new scheme is capable of achieving convergence when ϵ ≤ 10−18. As shown in Table VI, the computational cost of the new regularization scheme is quite reasonable within a large range of ϵ. As the coupling strength and number of modes increase, this model will become so ill-conditioned that the new scheme will eventually fail. For practical simulations using ML-MCTDH, one needs to seek a better representation or reformulation of the problem when possible (e.g., in this case, go back to the standard spin-boson model). The fact that the new scheme can handle a much broader regime than the standard scheme provides more flexibility for realistic applications and reformulations of the problem.

TABLE VI.

Number of calls to the ML-MCTDH derivative subroutine and deviation versus the regularization parameter ϵ for the polaron-transformed spin-boson model. A four-layer ML-MCTDH is used to treat a bath of 512 modes. The physical parameter are D/Δ = 0, α = 2, and ωc/Δ = 2.5. Only the results under the new regularization scheme are shown since the standard scheme completely fails.

ϵNumber of callsDeviation
10−10 2747 6.7 × 10−1 
10−12 2077 1.6 × 10−1 
10−14 2498 1.6 × 10−2 
10−16 3130 1.6 × 10−3 
10−18 4386 1.6 × 10−4 
10−20 6376 4.5 × 10−5 
10−22 6144 2.9 × 10−5 
10−24 6221 Reference 
ϵNumber of callsDeviation
10−10 2747 6.7 × 10−1 
10−12 2077 1.6 × 10−1 
10−14 2498 1.6 × 10−2 
10−16 3130 1.6 × 10−3 
10−18 4386 1.6 × 10−4 
10−20 6376 4.5 × 10−5 
10−22 6144 2.9 × 10−5 
10−24 6221 Reference 

In passing, we note that Manthe has recently also discussed61 the polaron-transformed spin-boson model for a somewhat weaker coupling strength. He correctly solved the dynamical problem by using optimized60 initially unoccupied SPFs. Hence this approach may serve as an alternative to the use of the new regularization scheme. It is also possible, and presumably more effective, to combine Manthe’s choice of initial SPFs with our new regularization scheme.

In this work, our previously proposed new regularization of the MCTDH equations of motion is generalized to ML-MCTDH. In this scheme, the expansion coefficients of all layers, instead of the density matrices, in ML-MCTDH are regularized through their matrix unfoldings. This automatically regularizes the density matrices. Despite the principle idea remaining the same, the ML-MCTDH equations are much more complicated and this generalization is not obvious. For practical purposes, we have derived a recursive implementation of the new regularization scheme.

We have tested the new regularization scheme on several examples of the spin-boson model whose parameter set ranges from intermediate to very strong couplings and the size of the system ranges from a few hundred to 5000 degrees of freedom. These are very realistic applications of ML-MCTDH. Since singularities occur more frequently in ML-MCTDH, it is more important to adopt the new regularization scheme.

In general, it is found that the new scheme is much more robust with respect to changing the regularization parameter ϵ. For weak to intermediate coupling strength, it is still possible to achieve convergence in the standard scheme by reducing ϵ. However, this becomes prohibitively expensive in challenging applications where the coupling strength becomes strong. In those situations, the new regularization scheme becomes the only practical choice. Furthermore, it is often preferred not to check so many parameters in real applications of ML-MCTDH. The fact that the new regularization scheme is stable over a large range of the regularization parameter ϵ makes it possible to use a conservative default value and not worry about artifacts induced otherwise.

Because the coefficient tensors in ML-MCTDH are usually smaller than in MCTDH, the cost of SVDs in the new scheme is often only moderately more expensive than the density matrix diagonalization in the standard scheme for routine applications (for challenging applications, the latter is neither reliable nor practical.) Still, owing to the fact that orbital rotations are the most important at the beginning, it is possible to combine the new regularization scheme with the standard scheme such that the latter can be switched on when the natural populations become large enough.

In some previous work, the singularity for the short time (ML)-MCTDH propagation raises concerns and effort has been made to remove it. Our work shows that this singularity is needed to rapidly rotate the unoccupied SPFs into their optimal direction in Hilbert space. If this is not done properly, not only the accuracy but also the efficiency of the ML-MCTDH propagation will be affected. Many times, a more rapid orbital rotation at the beginning will actually reduce the overall computational cost (see, e.g., Tables II, IV, and V) because of the exploration of the correct regions of the Hilbert space. The new regularization offers an efficient way of achieving this goal.

H.W. acknowledges the support from the National Science Foundation Grant No. CHE-1500285 and resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. H.-D.M. acknowledges financial support from the Deutsche Forschungsgemeinschaft (DFG).

1.
H.-D.
Meyer
,
U.
Manthe
, and
L. S.
Cederbaum
, “
The multi-configurational time-dependent Hartree approach
,”
Chem. Phys. Lett.
165
,
73
78
(
1990
).
2.
U.
Manthe
,
H.-D.
Meyer
, and
L. S.
Cederbaum
, “
Wave-packet dynamics within the multiconfiguration Hartree framework: General aspects and application to NOCl
,”
J. Chem. Phys.
97
,
3199
3213
(
1992
).
3.
M. H.
Beck
,
A.
Jäckle
,
G. A.
Worth
, and
H.-D.
Meyer
, “
The multi-configuration time-dependent Hartree (MCTDH) method: A highly efficient algorithm for propagating wave packets
,”
Phys. Rep.
324
,
1
105
(
2000
).
4.
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.
109
,
251
267
(
2003
).
5.
Multidimensional Quantum Dynamics: MCTDH Theory and Applications
, edited by
H.-D.
Meyer
,
F.
Gatti
, and
G. A.
Worth
(
Wiley-VCH
,
Weinheim
,
2009
).
6.
H.-D.
Meyer
, “
Studying molecular quantum dynamics with the multiconfiguration time-dependent Hartree method
,”
WIREs: Comput. Mol. Sci.
2
,
351
374
(
2012
).
7.
A.
Raab
,
G.
Worth
,
H.-D.
Meyer
, and
L. S.
Cederbaum
, “
Molecular dynamics of pyrazine after excitation to the S2 electronic state using a realistic 24-mode model Hamiltonian
,”
J. Chem. Phys.
110
,
936
946
(
1999
).
8.
C.
Cattarius
,
G. A.
Worth
,
H.-D.
Meyer
, and
L. S.
Cederbaum
, “
All mode dynamics at the conical intersection of an octa-atomic molecule: Multi-configuration time-dependent Hartree (MCTDH) investigation on the butatriene cation
,”
J. Chem. Phys.
115
,
2088
2100
(
2001
).
9.
S.
Mahapatra
,
G. A.
Worth
,
H. D.
Meyer
,
L. S.
Cederbaum
, and
H.
Köppel
, “
The Ã2E B̃2B2 photoelectron bands of allene beyond the linear coupling scheme: An ab initio dynamical study including all fifteen vibrational modes
,”
J. Phys. Chem. A
105
,
5567
5576
(
2001
).
10.
A.
Markmann
,
G.
Worth
,
S.
Mahapatra
,
H.-D.
Meyer
,
H.
Köppel
, and
L.
Cederbaum
, “
Simulation of a complex spectrum: Interplay of five electronic states and 21 vibrational degrees of freedom in C5H4+.
,”
J. Chem. Phys.
123
,
204310
(
2005
).
11.
O.
Vendrell
,
F.
Gatti
,
D.
Lauvergnat
, and
H.-D.
Meyer
, “
Full dimensional (15D) quantum-dynamical simulation of the protonated water dimer. I. Hamiltonian setup and analysis of the ground vibrational state
,”
J. Chem. Phys.
127
,
184302
(
2007
).
12.
O.
Vendrell
,
F.
Gatti
, and
H.-D.
Meyer
, “
Full dimensional (15D) quantum-dynamical simulation of the protonated water dimer. II. Infrared spectrum and vibrational dynamics
,”
J. Chem. Phys.
127
,
184303
(
2007
).
13.
O.
Vendrell
,
M.
Brill
,
F.
Gatti
,
D.
Lauvergnat
, and
H.-D.
Meyer
, “
Full dimensional (15D) quantum-dynamical simulation of the protonated water dimer. III. Mixed Jacobi-valence parametrization and benchmark results for the zero-point energy, vibrationally excited states and infrared spectrum
,”
J. Chem. Phys.
130
,
234305
(
2009
).
14.
O.
Vendrell
,
F.
Gatti
, and
H.-D.
Meyer
, “
Full dimensional (15D) quantum-dynamical simulation of the protonated water dimer. IV. Isotope effects in the infrared spectra of D(D2O)2+, H(D2O)2+ and D(H2O)2+ isotopologues
,”
J. Chem. Phys.
131
,
034308
(
2009
).
15.
M.
Schröder
,
F.
Gatti
, and
H.-D.
Meyer
, “
Theoretical studies of the tunneling splitting of malonaldehyde using the multiconfiguration time-dependent Hartree approach
,”
J. Chem. Phys.
134
,
234307
(
2011
).
16.
M.
Schröder
and
H.-D.
Meyer
, “
Calculation of the vibrational excited states of malonaldehyde and their tunneling splittings with the multi-configuration time-dependent Hartree method
,”
J. Chem. Phys.
141
,
034116
(
2014
).
17.
T.
Hammer
and
U.
Manthe
, “
Intramolecular proton transfer in malonaldehyde: Accurate multilayer multi-configurational time-dependent Hartree calculations
,”
J. Chem. Phys.
134
,
224305
(
2011
).
18.
T.
Hammer
and
U.
Manthe
, “
Iterative diagonalization in the state-averaged multi-configurational time-dependent Hartree approach: Excited state tunneling splitting in malonaldehyde
,”
J. Chem. Phys.
136
,
054105
(
2012
).
19.
G.
Schiffel
and
U.
Manthe
, “
Quantum dynamics of the H + CH4 → H2 + CH3 reaction in curvilinear coordinates: Full-dimensional and reduced dimensional calculations of reaction rates
,”
J. Chem. Phys.
132
,
084103
(
2010
).
20.
G.
Schiffel
and
U.
Manthe
, “
Communications: A rigorous transition state based approach to state-specific reaction dynamics: Full-dimensional calculations for H + CH4 → H2 + CH3
,”
J. Chem. Phys.
132
,
191101
(
2010
).
21.
G.
Schiffel
and
U.
Manthe
, “
A transition state view on reactive scattering: Initial state-selected reaction probabilities for the H + CH4 → H2 + CH3 reaction studied in full dimensionality
,”
J. Chem. Phys.
133
,
174124
(
2010
).
22.
R.
Ellerbrock
and
U.
Manthe
, “
Communication: Reactivity borrowing in the mode selective chemistry of H + CHD3 → H2 + CD3
,”
J. Chem. Phys.
147
,
241104
(
2017
).
23.
H.
Wang
, “
Basis set approach to the quantum dissipative dynamics: Application of the multiconfiguration time-dependent Hartree method to the spin-boson problem
,”
J. Chem. Phys.
113
,
9948
(
2000
).
24.
H.
Wang
,
M.
Thoss
, and
W.
Miller
, “
Systematic convergence in the dynamical hybrid approach for complex systems: A numerical exact methodology
,”
J. Chem. Phys.
115
,
2979
(
2001
).
25.
M.
Thoss
and
H.
Wang
, “
Quantum dynamical simulation of ultrafast photoinduced electron transfer processes in a mixed-valence compound
,”
Chem. Phys. Lett.
358
,
298
(
2002
).
26.
H.
Wang
and
M.
Thoss
, “
Theoretical study of ultrafast photoinduced electron transfer processes in mixed-valence systems
,”
J. Phys. Chem. A
107
,
2126
2136
(
2003
).
27.
H.
Wang
and
M.
Thoss
, “
Multilayer formulation of the multiconfiguration time-dependent Hartree theory
,”
J. Chem. Phys.
119
,
1289
1299
(
2003
).
28.
U.
Manthe
, “
A multilayer multiconfigurational time-dependent Hartree approach for quantum dynamics on general potential energy surfaces
,”
J. Chem. Phys.
128
,
164116
(
2008
).
29.
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.
134
,
044135
(
2011
).
30.
H.
Wang
, “
Multilayer multiconfiguration time-dependent Hartree theory
,”
J. Phys. Chem. A
119
,
7951
(
2015
).
31.
T.
Westermann
,
R.
Brodbeck
,
A. B.
Rozhenko
,
W.
Schoeller
, and
U.
Manthe
, “
Photodissociation of methyl iodide embedded in a host-guest complex: A full dimensional (189D) quantum dynamics study of CH3I@resorc[4]arene
,”
J. Chem. Phys.
135
,
184102
(
2011
).
32.
Q.
Meng
and
H.-D.
Meyer
, “
A multilayer MCTDH study on the full dimensional vibronic dynamics of naphthalene and anthracene cations
,”
J. Chem. Phys.
138
,
014313
(
2013
).
33.
R.
Welsch
and
U.
Manthe
, “
Reaction dynamics with the multi-layer multi-configurational time-dependent Hartree approach: H + CH4 → H2 + CH3 rate constants for different potentials
,”
J. Chem. Phys.
137
,
244106
(
2012
).
34.
R.
Welsch
and
U.
Manthe
, “
The role of the transition state in polyatomic reactions: Initial state-selected reaction probabilities of the H + CH4 → H2 + CH3 reaction
,”
J. Chem. Phys.
141
,
174313
(
2014
).
35.
J.
Schulze
and
O.
Kühn
, “
Explicit correlated exciton-vibrational dynamics of the FMO complex
,”
J. Phys. Chem. B
119
,
6211
(
2015
).
36.
J.
Schulze
,
M. F.
Shibl
,
M. J.
Al-Marri
, and
O.
Kühn
, “
Multi-layer multi-configuration time-dependent Hartree (ML-MCTDH) approach to the correlated exciton-vibrational dynamics in the FMO complex
,”
J. Chem. Phys.
144
,
185101
(
2016
).
37.
M. F.
Shibl
,
J.
Schulze
,
M. J.
Al-Marri
, and
O.
Kühn
, “
Multilayer-MCTDH approach to the energy transfer dynamics in the LH2 antenna complex
,”
J. Phys. B
50
,
184001
(
2017
).
38.
J.
Zheng
,
Y.
Xie
,
S.
Jiang
, and
Z.
Lan
, “
Ultrafast nonadiabatic dynamics of singlet fission: Quantum dynamics with the multilayer multiconfigurational time-dependent Hartree (ML-MCTDH) method
,”
J. Phys. Chem. C
120
,
1375
(
2016
).
39.
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.
146
,
184305
(
2017
).
40.
K. A.
Velizhanin
,
H.
Wang
, and
M.
Thoss
, “
Heat transport through model molecular junctions: A multilayer multiconfiguration time-dependent Hartree approach
,”
Chem. Phys. Lett.
460
,
325
(
2008
).
41.
H.
Wang
,
I.
Pshenichnyuk
,
R.
Härtle
, and
M.
Thoss
, “
Numerically exact, time-dependent treatment of vibrationally coupled electron transport in single-molecule junctions
,”
J. Chem. Phys.
135
,
244506
(
2011
).
42.
K. F.
Albrecht
,
H.
Wang
,
L.
Mühlbacher
,
M.
Thoss
, and
A.
Komnik
, “
Bistability signatures in nonequilibrium charge transport through molecular quantum dots
,”
Phys. Rev. B
86
,
081412
(
2012
).
43.
H.
Wang
and
M.
Thoss
, “
Numerically exact, time-dependent study of correlated electron transport in model molecular junctions
,”
J. Chem. Phys.
138
,
134704
(
2013
).
44.
H.
Wang
and
M.
Thoss
, “
Multilayer multiconfiguration time-dependent Hartree study of vibrationally coupled electron transport using the scattering-state representation
,”
J. Phys. Chem. A
117
,
7431
(
2013
).
45.
E. Y.
Wilner
,
H.
Wang
,
G.
Cohen
,
M.
Thoss
, and
E.
Rabani
, “
Bistability in a nonequilibrium quantum system with electron-phonon interactions
,”
Phys. Rev. B
88
,
045137
(
2013
).
46.
D.
Mendive-Tapia
,
E.
Mangaud
,
T.
Firmino
,
A.
de la Lande
,
M.
Desouter-Lecomte
,
H.-D.
Meyer
, and
F.
Gatti
, “
Multidimensional quantum mechanical modeling of electron transfer and electronic coherence in plant cryptochromes: The role of initial bath conditions
,”
J. Phys. Chem. B
122
,
126
136
(
2018
).
47.
M.
Thoss
,
I.
Kondov
, and
H.
Wang
, “
Theoretical study of ultrafast heterogeneous electron transfer reactions at dye-semiconductor interfaces
,”
Chem. Phys.
304
,
169
(
2004
).
48.
I.
Kondov
,
M.
Thoss
, and
H.
Wang
, “
Theoretical study of ultrafast heterogeneous electron transfer reactions at dye-semiconductor interfaces: Coumarin 343 at titanium oxide
,”
J. Phys. Chem. A
110
,
1364
1374
(
2006
).
49.
H.
Wang
and
M.
Thoss
, “
Quantum dynamical simulation of electron-transfer reactions in an anharmonic environment
,”
J. Phys. Chem. A
111
,
10369
(
2007
).
50.
I.
Kondov
,
M.
Thoss
, and
H.
Wang
, “
Quantum dynamics of photoinduced electron transfer reactions in dye-semiconductor systems: Description and application to coumarin 343-TiO2
,”
J. Phys. Chem. C
111
,
11970
11981
(
2007
).
51.
M.
Thoss
,
I.
Kondov
, and
H.
Wang
, “
Correlated electron-nuclear dynamics in ultrafast photoinduced electron-transfer reactions at dye-semiconductor interfaces
,”
Phys. Rev. B
76
,
153313
(
2007
).
52.
I. R.
Craig
,
H.
Wang
, and
M.
Thoss
, “
Proton transfer reactions in model condensed-phase environments: Accurate quantum dynamics using the multilayer multiconfiguration time-dependent Hartree approach
,”
J. Chem. Phys.
127
,
144503
(
2007
).
53.
H.
Wang
and
M.
Thoss
, “
Nonperturbative quantum simulation of time-resolved nonlinear spectra: Methodology and application to electron transfer reactions in the condensed phase
,”
Chem. Phys.
347
,
139
(
2008
).
54.
H.
Wang
and
M.
Thoss
, “
From coherent motion to localization: Dynamics of the spin-boson model at zero temperature
,”
New J. Phys.
10
,
115005
(
2008
).
55.
D.
Egorova
,
M. F.
Gelin
,
M.
Thoss
,
H.
Wang
, and
W.
Domcke
, “
Effects of intense femtosecond pumping on ultrafast electronic-vibrational dynamics in molecular systems with relaxation
,”
J. Chem. Phys.
129
,
214303
(
2008
).
56.
K. A.
Velizhanin
and
H.
Wang
, “
Dynamics of electron transfer reactions in the presence of mode-mixing: Comparison of a generalized master equation approach with the numerically exact simulation
,”
J. Chem. Phys.
131
,
094109
(
2009
).
57.
H.
Wang
and
M.
Thoss
, “
From coherent motion to localization. II. Dynamics of the spin-boson model with sub-ohmic spectral density at zero temperature
,”
Chem. Phys.
370
,
78
86
(
2010
).
58.
Y.
Zhou
,
J.
Shao
, and
H.
Wang
, “
Dynamics of electron transfer in complex glassy environment modeled by the Cole-Davidson spectral density
,”
Mol. Phys.
110
,
581
(
2012
).
59.
Y.
Xie
,
J.
Zheng
, and
Z.
Lan
, “
Full-dimensional multilayer multiconfigurational time-dependent Hartree study of electron transfer dynamics in the anthracene/C60 complex
,”
J. Chem. Phys.
142
,
084706
(
2015
).
60.
U.
Manthe
, “
The multi-configurational time-dependent Hartree approach revisited
,”
J. Chem. Phys.
142
,
244109
(
2015
).
61.
U.
Manthe
, “
Optimized unoccupied single-particle functions in the (multi-layer) multi-configurational time-dependent Hartree approach
,”
Chem. Phys.
(in press).
62.
K.-S.
Lee
and
U. R.
Fischer
, “
Truncated many-body dynamics of interacting bosons: A variational principle with error monitoring
,”
Int. J. Mod. Phys. B
28
,
1550021
(
2014
).
63.
H.-D.
Meyer
and
H.
Wang
, “
On regularizing the MCTDH equations of motion
,”
J. Chem. Phys.
148
,
124105
(
2018
).
64.

Here we follow the nomenclature introduced by Manthe.28 One of us has previously adopted an alternate counting of the layers,30 which differs by one unit and considers MCTDH as a one-layer version of ML-MCTDH.

65.
A. J.
Leggett
,
S.
Chakravarty
,
A. T.
Dorsey
,
M. P. A.
Fisher
,
A.
Garg
, and
W.
Zwerger
, “
Dynamics of the dissipative two-state System
,”
Rev. Mod. Phys.
59
,
1
(
1987
).
66.
U.
Weiss
,
Quantum Dissipative Systems
, 3rd ed. (
World Scientific Publishing
,
Singapore
,
2008
).