The projector-splitting integrator for the Multi-Configuration Time-Dependent Hartree method revisited

# The projector-splitting integrator for the Multi-Configuration Time-Dependent Hartree method revisited

Matteo Bonfanti    Irene Burghardt Institute of Physical and Theoretical Chemistry, Goethe University Frankfurt,
Max-von-Laue-Str. 7, D-60438 Frankfurt/Main, Germany
July 12, 2019
###### Abstract

We revisit a recently introduced algorithm [C. Lubich, Appl. Math. Res. eXpress 2015, 311 (2015), B. Kloss et al., J. Chem. Phys. 146, 174107 (2017)] for the integration of the Multi-Configuration Time-Dependent Hartree (MCTDH) equations for high-dimensional quantum propagation. The new integrator circumvents the direct inversion of reduced density matrices that appears in the standard method, by employing an auxiliary set of non-orthogonal single-particle functions. Here, we review the formulation of this new algorithm, with the aim of providing a bridge between the conventional formulation of MCTDH and the original, more mathematically oriented derivation of the new algorithm within a tensor formalism. We recast the alternative form of the MCTDH equations underlying the algorithm in conventional language, and highlight key features of the integration scheme. Furthermore, we emphasize the derivation of the method in the broader context of the time-dependent variational problem, where the new equations of motion are naturally obtained from a specific splitting of the tangent-space projector.

preprint: APS/123-QED

## I Introduction

A number of recent contributions in the mathematical literature have reviewed the Multi-Configuration Time-Dependent Hartree (MCTDH) methodMeyer, Manthe, and Cederbaum (1990); Manthe, Meyer, and Cederbaum (1992); Beck et al. (2000) and its multi-layer (ML-MCTDH) variantWang and Thoss (2003); Manthe (2008); Vendrell and Meyer (2011) from the viewpoint of low-rank tensor approximation techniques.Koch and Lubich (2010); Lubich et al. (2013); Bachmayr, Schneider, and Uschmajew (2016) Among these developments, LubichLubich (2015) proposed a novel MCTDH integration algorithm, which was later implemented and tested on low-dimensional model systems.Kloss, Burghardt, and Lubich (2017) This algorithm relies on the splitting of the tangent-space projectionLubich (2015) and is, hence, termed projector-splitting integrator. This new integration scheme is the focus of the present work.

Within the standard MCTDH approach, the equations of motion for the single–particle functions (SPFs) – i.e., the time-dependent basis of MCTDH – are highly nonlinear because they involve the inverse of a single-particle density matrix, . The advantage of the projector-splitting integrator lies in the fact that a direct inversion of the density matrix is circumvented and the equations are recast in a linear formLubich (2015) (noting that linearity here refers to the form of the equations, while nonlinearity due to the presence of mean-field potentials remains a feature of the new scheme). Potential (near-)singularities of are dealt with at the level of a QR decomposition, which is known to be robust even in the case of matrices with large condition numbers.Golub and Van Loan (2013)

An extensive number of MCTDH applications,Meyer, Gatti, and Worth (2009); Meyer (2012) spread across all fields of quantum dynamics, show that the method in its original form is generally robust and convergeable and that the regularized inversion that is used for ill–conditioned density matrices rarely affects the quality of the results. However, the hierarchical ML-MCTDH variant was found to be more sensitive to initial conditions and to the regularization parameter.Manthe (2015) Furthermore, numerical analysis has raised some concerns regarding the convergence to the exact solution in cases where ill-conditioned density matrices appear.Conte and Lubich (2010) Indeed, problems were reported in several cases described in the literature, e.g., related to the fermionic variant of MCTDH (i.e., MCTDH-F) where the sensitivity of the results to the regularization parameter is found to be increased for certain classes of systems.Hinz, Bauch, and Bonitz (2016)

In addition to bringing improvements in these specific cases, the development of novel integration algorithms may suggest new strategies to avoid singularities in the general context of variational equations of motion. This is a problem that, e.g., seriously affects methods based on non-orthogonal basis functions such as the Gaussian-based MCTDH (G-MCTDH) methodBurghardt, Meyer, and Cederbaum (1999); Burghardt, Nest, and Worth (2003); Burghardt, Giri, and Worth (2008); Römer, Ruckenbauer, and Burghardt (2013) and its variational multi-configurational Gaussian (vMCG) variant.Worth et al. (2004); Worth and Burghardt (2003); Richings et al. (2015)

Against this background, the purpose of the present paper is twofold. First, we aim to provide a bridge between the conventional formulation of MCTDH and some of the more mathematically oriented developments mentioned above, which are usually formulated in tensor languageBachmayr, Schneider, and Uschmajew (2016) that may not be immediately accessible to the chemical physics community. In particular, we will focus on the derivation and formulation of the novel projector splitting scheme.Lubich (2015); Kloss, Burghardt, and Lubich (2017) Second, we give special emphasis to the general notion of the tangent–space projection in the framework of time–dependent variational problems. To show how these concepts can be applied in practice, a detailed discussion of the tangent space projectorKoch and Lubich (2010) for MCTDH will be given.

The outline of the remainder of this article is as follows. In section II we briefly explain notational issues, and in section III we review the notion of tangent-space projections. In section IV, we discuss the dynamical equations of MCTDH in the form of the projector-splitting algorithm. In Section V, the integration scheme is commented upon. Appendix A contains a brief key to translation between the tensorial and standard notation, and Appendix B provides details of the derivation of the projector-splitting algorithm.

## Ii Notation

We start by giving a brief description of the notation that will be adopted in this paper. Generally, we will adhere to the standard conventions of the MCTDH literature.Beck et al. (2000)

We seek a solution to the time-dependent Schrödinger Equation (TDSE) for a multidimensional state by approximating the Hilbert space as a tensor product of subspaces of low–dimensional SPFs. The wavefunction is then represented according to the usual MCTDH ansatz,

 Ψ({rκ},t)=n1∑j1=1n2∑j2=1…nf∑jf=1Aj1,j2,…jf(t)f∏κ=1φ(κ)jκ(rκ,t) (1)

where is the -th SPF for mode and is the tensor of the expansion coefficients. The SPFs are defined to be orthogonal at all times, , benefitting from the gauge freedom of the MCTDH ansatz Eq. (1).Beck et al. (2000) This standard gauge also implies that . More generally, the gauge can be defined in terms of constraint operators.Beck et al. (2000)

In the tensor formulation that is adopted in the mathematical literature, the MCTDH expansion of Eq. (1) is equivalently interpreted as a reduction of the dimensionality of the coefficient tensor. This is made evident by projecting the expansion Eq. (1) on a time-independent product basis ,

 Ψ({rκ},t)=N1∑i1=1N2∑i2=1…Nf∑if=1Yi1,i2,…if(t)f∏κ=1χ(κ)iκ(rκ) (2)

with

 Yi1,i2,…if(t)=n1∑j1=1…nf∑jf=1Aj1,j2,…jf(t)f∏κU(κ)iκjκ(t) (3)

where is the representation matrix on the primitive grid of the –mode SPFs. From a tensor algebra perspective, Eq. (3) is known as Tucker decomposition of the tensor into the core tensorBachmayr, Schneider, and Uschmajew (2016) and the set of matrices . As the number of SPFs is obviously smaller than the size of the primitive basis, the Tucker decomposition entails a reduction in dimensionality of the original tensor, taking advantage of its possible sparsity.

Following standard practice,Beck et al. (2000) we make use of multi–indices to cast Eq. (1) in a more compact form (omitting the explicit time and coordinate dependence),

 Ψ=∑JAJΦJ (4)

where represents an -dimensional vector of indices and represents a configuration. Due to the orthonormality of the SPFs, the configurations are orthonormal as well, .

In our discussion of the MCTDH equations of motion, we will make use of two additional conventions for multi–indices.Beck et al. (2000) In situations where a summation is carried out over all indices except one, we introduce a reduced multi-index,

 Jκ=(j1,j2,…jκ−1,jκ+1…jf) (5)

When it is necessary to label a tensor with a multi–index with the -th entry substituted with another integer , we write the modified multi–index as

 Jκl=(j1,j2,…jκ−1,l,jκ+1…jf) (6)

With these two definitions, we can define single-hole functions (SHFs) as

 Ψ(κ)l=∑JκAJκl∏n≠κφ(n)jn (7)

and the wavefunction Eq. (1) can be re-written as a product of SPFs and SHFs,Beck et al. (2000)

 Ψ=∑lκφ(κ)lκΨ(κ)lκ (8)

which is most convenient when equations are defined within a given th subspace. In terms of the SHFs, we can further write the –mode single–particle density matrix as the overlap of SHFs,

 ρ(κ)lκl′κ=⟨Ψ(κ)lκ|Ψ(κ)l′κ⟩ (9)

noting that where is the reduced density operator in the th subspace.

## Iii Tangent-space projection of the Time-Dependent Schrödinger Equation

The projector-splitting schemeLubich (2015); Kloss, Burghardt, and Lubich (2017) is best understood when the equations of motion are derived in terms of a tangent-space projection of the TDSE. This is equivalent to the use of the Dirac-Frenkel Variational Principle (DFVP) to derive the MCTDH equations.Beck et al. (2000) Here, we state the main results and refer to to Ref. [Lubich, 2008] for further background from a mathematical perspective.

### iii.1 Tangent space projection

The conventional formulation of the DFVP states that the best approximation to the time evolving wavefunction at a given time is obtained as the function which satisfiesLubich (2008)

 ⟨δΨ|˙Ψ−1ıℏHΨ⟩=0 (10)

where is an allowed variation of the wavefunction which is compatible with the chosen ansatz (see, e.g., Eq. (1)). In mathematical terms, if is the smooth submanifold of the Hilbert space in which we are seeking an approximation to the time–dependent state, is an element of , the tangent space of at .

In Ref. [Broeckhove et al., 1988], it was shown that the DFVP is equivalent to the McLachlan Variational Principle (MLVP),McLachlan (1964) provided that the tangent space is a complex linear space. According to the MLVP, the best solution to the time–dependent problem is obtained when we approximate the exact derivative with a vector that has the minimal distance – with the metric induced by the scalar product – from the exact derivative of the wavefunction,Kucar, Meyer, and Cederbaum (1987) which is given by the TDSE as . That is,

 δ∥˙Ψ−1ıℏHΨ∥2=0 (11)

This “geometrical” condition can be solved by introducing tangent–space projectorsRaab (2000); Lubich (2004, 2008) such that, at any time, the “best” approximate derivative is constructed as the orthogonal projection of the full derivative onto the tangent space . Hence,Lubich (2008)

 ˙Ψ=P(Ψ)1ıℏHΨ (12)

where is the orthogonal projector onto the tangent space . Given that is a vector of by construction, it belongs to the range of the projector . Thus Eq. (12) can be rearranged as a projected TDSE:Lubich (2008)

 P(Ψ)[˙Ψ−1ıℏHΨ]=0 (13)

In practice, Eq. (13) is the most convenient form to derive the equations of motion from the MLVP, once an explicit formula for the projector is known.

The notion of the tangent–space projection, pictorially represented in Fig. 1, gives an illuminating and immediate understanding of the variational principle and its implications. When the shape of the wavefunction is defined according to a chosen ansatz, we are constraining the time evolution to a submanifold of the full Hilbert space. The dynamics that is returned by the DFVP is such that at any given time the derivative of the approximate wavefunction is optimal, in the sense that it is closest to the exact value of the derivative at that point . However, is obviously constrained to reside within the tangent space since the evolving wavefunction cannot “escape” from . No global condition is given for the dynamical propagation; instead, the approximation is chosen such as to guarantee that at any time the wavefunction evolution diverges the least possible from the exact dynamics.

As a consequence, one cannot exclude that small errors that are incurred at each instant of the propagation may add up to a large deviation of the dynamics from the exact evolution at longer times. Of course, as the submanifold approaches the full size of the space, the projector approaches unity and the approximated evolution tends to the prediction of the TDSE.

### iii.2 Tangent space projection for the MCTDH ansatz

From a practical point of view, the construction of the tangent–space projection for a specific wavefunction ansatz can be obtained from the different linearly independent components that constitute the first–order variation of the wavefunction.

To illustrate the construction of the tangent–space projector, we now turn to the MCTDH case as an example. In this case, the tangent space is a complex linear space, since both the vector and the SPFs are assumed to be complex-valued. As a consequence, the DFVP and the MPVP are equivalent, and both are in turn equivalent to a least–action principle.Broeckhove et al. (1988)

Referring to Eq. (1), the first-order variation of the MCTDH ansatz is given by

 δΨ=∑JδAJΦJ+∑κ(∑lκδφ(κ)lκΨ(κ)lκ) (14)

i.e., a sum of terms relating to the variation of the coefficients and the SPFs in the th subspaces. As a remark about notation, indicates the variation of the corresponding quantity, specifically is the variation of the –th vector component, and is the variation of the –th SPF for mode . Eq. (14) then represents a generic vector of . In contrast, we will use dotted quantities , and to indicate the specific tangent-space vectors resulting from the DFVP.

The tangent-space projector of Eq. (12) directly relates to the first order wavefunction variation of Eq. (14), such that , i.e., the first-order variation lies in the tangent space by construction. Hence, we naturally aim to construct with a similar partitioning as Eq. (14),Koch and Lubich (2010)

 P(Ψ)=P0(Ψ)+∑κPκ(Ψ) (15)

While it is tempting to construct a one-to-one correspondence between the components of in Eq. (15) and the components of in Eq. (14), we will need to make sure that and refer to orthogonal projections. In the following, we use the notation for the parts of the linear variation that relate to the subprojections of Eq. (15), and we will show below how these connect to the r.h.s. of Eq. (14).

To start with, is chosen as the projector onto the configurations ,

 P0(Ψ)=∑J|ΦJ⟩⟨ΦJ| (16)

which, when acting on a generic state, expands this state as a linear combination of configurations . The associated portion of the linear variation Eq. (14), here denoted , reads

 δΨ0 = P0(Ψ)δΨ = ∑JδAJΦJ+P0(Ψ)∑κ,lκδφ(κ)lκΨ(κ)lκ (17)

Since not only relates to the coefficient variation but also acts on the second term corresponding to SPF variations, the definition of the remaining subspace projections of Eq. (15) needs to be chosen such as to project onto the complementary subspace of the range of :

 ∑κδΨκ = ∑κPκ(Ψ)δΨ = (1−P0(Ψ))∑κ,lδφ(κ)lΨ(κ)l (18)

With a few lines of algebra, using again the orthogonality of the SPFs, we find that

 ∑κδΨκ = ∑κ,lΨ(κ)l(1−P(κ))δφ(κ)l (19)

where

 P(κ)=∑i|φ(κ)i⟩⟨φ(κ)i| (20)

is the projector onto the space spanned by the -mode SPFs.Beck et al. (2000) This is the projector appearing in the conventional MCTDH equations. To underline the difference between the different types of projectors, we indicate the tangent–space projector and its components with the calligraphic letter while the subspace projector is given in Roman type .

Eq. (19) shows that the second part of the first–order variation of the wavefunction is spanned by products between the SHFs and SPF variations, with the latter being constrained to the orthogonal complement of the SPFs. In the final equations of motion, this condition will guarantee that the SPF propagation does not involve changes which are already represented by the time evolution of the coefficients. From Eq. (19) we further infer that variations corresponding to different modes are orthogonal, namely

 ⟨δΨκ|δΨκ′⟩=0 (21)

as can be seen, e.g., by letting the projector act on the ket .

From the above, we now identify the projectors as the tensor product of two subspace projectors,

 Pκ(Ψ)=(1−P(κ))⊗¯P(κ) (22)

where the second projector on the r.h.s. refers to the space spanned by the -mode SHFs,

 ¯P(κ)=∑l,l′|Ψ(κ)l′⟩(ρ(κ))−1l′l⟨Ψ(κ)l| (23)

Importantly, represents a projector onto a non–orthogonal basis, which includes the inverse of the overlap matrix. In the above, we used the definition Eq. (9), i.e., the SHF overlap coincides with the single-particle density matrix .

To summarize, Eqs. (15), (16), and (22) fully define the tangent space projection for the MCTDH ansatz.

Anticipating the discussion of Sec. V, the novel projector-splitting algorithm will be shown to rely on a partitioning of the projector of Eq. (22) into two components, , see Eq. (41) below, permitting a new partitioning of the equations of motion.

### iii.3 MCTDH equations of motion

The terms constituting in Eq. (15) give rise to the equations of the MCTDH standard formulation, when applied to the TDSE as partial projections according to Eq. (13). Notably, we will consider the conditions

 P0(Ψ)[˙Ψ−1ıℏHΨ]=0;Pκ(Ψ)[˙Ψ−1ıℏHΨ]=0

where

 ˙Ψ−1ıℏHΨ=∑J˙AJΦJ+∑κ,l˙φ(κ)lΨ(κ)l−1ıℏ∑JAJHΦJ (25)

assuming the standard MCTDH gauge that keeps the SPFs orthonormal during the propagation.

The projection then yields the standard differential equation for the vector,

 ˙AI=1ıℏ∑J⟨ΦI|H|ΦJ⟩AJ (26)

whereas the projection along returns the differential equation for the SPFs of mode ,

 ˙φ(κ)=1ıℏ(1−P(κ))(ρ(κ))−1⟨H⟩(κ)φ(κ) (27)

where is the vector composed by the SPFs for mode and is the the usual mean-field potential matrix, given by

 ⟨H⟩(κ)jk=⟨Ψ(κ)j|H|Ψ(κ)k⟩ (28)

Next, we turn to the reformulation of the MCTDH equations according to Ref. [Lubich, 2015].

## Iv The Projector–Splitting Equations of Motion

From the definition of the tangent-space projector Eq. (15) and the resulting MCTDH equations, it is clear that the SHF projector Eq. (23) is at the origin of the inverse of the density matrix appearing in Eq. (27). Hence, one can envisage an orthogonalizing transformation in the SHF space,Lubich (2015) such that the projector of Eq. (23) takes the alternative form,

 ¯P(κ)=∑l|~Ψ(κ)l⟩⟨~Ψ(κ)l| (29)

This concept is a key ingredient of the projector-splitting algorithm. As a trade-off for the resulting simplification of the equations of motion of the SPFs, the time-dependent transformation between non-orthogonal and orthogonalized SHFs has to be taken into account. As will be shown below, this can be conveniently achieved in terms of the splitting of the projectors appearing in Eq. (22).

### iv.1 SHF orthogonalization

First, we focus on the SHF orthogonalization and give a detailed description of the new quantities introduced by this transformation. By construction, the SHFs within the th subspace, , are non-orthogonal, and their overlap is given in terms of the reduced density matrix , see Eq. (9). For the purpose of the present discussion, we re-write the latter as follows,

 ρ(κ)ll′=∑JκA⋆JκlAJκl′ (30)

where the vectors were defined in Eq. (7). In the following, we will interpret as a matrix composed of column vectors of length obtained by fixing the -th index of the tensor. (In the tensor formulation, these vectors constitute a matrix which is defined as the -mode matricisation of the tensor , see Appendix A.)

Assuming that the column vectors of are linearly independent, we define a linear transformation that brings them in orthonormal form, which we conveniently write as

 AJκl=∑l′S(κ)ll′Q(κ)Jκl′ (31)

where is a lower triangular matrix and the tensor is composed of orthogonal vectors in the sense that was discussed above for the tensor. If we interpret both and as matrices, with indices contained in a single multi-index , we see that Eq. (31) corresponds to a QR decomposition, i.e., the decomposition of the matrix as a product of an orthogonal matrix times a triangular matrix,

 A(κ)=Q(κ)S(κ)T (32)

where is of dimension and is of dimension .

From a numerical viewpoint, the QR decomposition is generally stable and robust, even in cases where the orthogonalization is nontrivial because of near linear dependencies of the vectors extracted from (which is equivalent to the ill-conditioning of the density matrix, in light of Eq. (30)).

Substituting the QR decomposition of Eq. (31) into the MCTDH ansatz, we can see that the matrix effectively gives rise to a non-unitary transformation of the SPFs of the -th mode,

 Ψ=∑JQ(κ)Jφ(1)j1φ(2)j2…[∑lφ(κ)lS(κ)ljκ]… (33)

In other words, is the tensor of the coefficients of an equivalent MCTDH expansion in which the -mode SPFs are substituted with non-orthogonal functions defined by the triangular matrix transformation ,

 ~φ(κ)=S(κ)Tφ(κ) (34)

This transformation can also be understood as a QR decomposition if the representation of the SPFs in the primitive representation is considered as in Eq. (3).

Within the transformed representation, the new SHFs are defined in accordance with the standard definition, namely

 ~Ψ(κ)l=∑JκQ(κ)Jκlf∏n≠κφ(n)jn (35)

Importantly, the SHFs are now orthonormal by construction, because of the orthonormality of the “vector cuts” along of ,

 ⟨~Ψ(κ)l|~Ψ(κ)l′⟩=∑JκQ⋆JκlQJκl′=δll′ (36)

In conclusion, an alternative SPF-SHF decomposition of the MCTDH wavefunction has been constructed,

 Ψ=∑l~φ(κ)l~Ψ(κ)l (37)

which is analogous to the original representation (see Eq. 8) with the difference that the SHFs are now orthogonal and the SPFs are not.

Finally, substituting the QR decomposition of Eq. (31) into the expression of the reduced density matrix of Eq. (30) we obtain

 ρ(κ)=S(κ)S(κ)† (38)

i.e., the reduced density matrix factorizes in the form of a Cholesky decomposition. From this, we can better understand that the uniqueness of the matrix is closely connected to the invertibility of , as the Cholesky decomposition is unique for strictly positive-definite matrices.

### iv.2 Splitting of subspace projections

Using the new SHFs of Eq. (35), we can now express the operator as

 Pκ(Ψ)=(1−∑i|φ(κ)i⟩⟨φ(κ)i|)⊗∑l|~Ψ(κ)l⟩⟨~Ψ(κ)l| (39)

where the orthogonality of the SHFs is made evident by the disappearance of the overlap matrix from the projector. Eq. (39) is a hybrid representation where we keep the SPFs in their orthogonal form. As will become clear in the following, this is motivated by the fact that we will construct a suitable subprojection that singles out the time derivative of .

We further divide each of the projectors of Eq. (15) into two new projectors, splitting into two components:

 P+κ(Ψ)= 1 ⊗∑l|~Ψ(κ)l⟩⟨~Ψ(κ)l| (40a) P−κ(Ψ)= ∑i|φ(κ)i⟩⟨φ(κ)i| ⊗∑l|~Ψ(κ)l⟩⟨~Ψ(κ)l| (40b)

such that the overall projector now reads

 P(Ψ)=P0(Ψ)+f∑κ=1(P+κ(Ψ)−P−κ(Ψ)) (41)

These newly defined projection operators give rise to a different formulation of the differential equations in the -subspaces, which is equivalent to the original MCTDH formulation but makes direct use of the new SPF-SHF decomposition. Meanwhile, Eq. (26) remains unchanged, since the projector is left unchanged by the projector splitting.

From the action of on the time-dependent Schrödinger equation, we now obtain the following expression for the propagation of the non-orthogonal SPFs,

 ˙~φ(κ)=1ıℏ⟨~H⟩(κ)~φ(κ) (42)

see App. B for a detailed derivation. In Eq. (42), the new mean-field potential is defined by integrating over the orthogonalized SHFs,

 ⟨~H⟩(κ)ll′=⟨~Ψ(κ)l|H|~Ψ(κ)l′⟩ (43)

The advantage of Eq. (42) over Eq. (27) is evident: the reduced density matrix has been incorporated in the expression for the SPFs and the evaluation of the expression no longer requires the inversion of a potentially singular matrix.

The price to pay for this transformation is that we are now dealing with non–orthogonal SPFs, and an additional differential equation appears which is generated by the projector . Notably, we obtain

 ˙S(κ)ij=1ıℏ⟨φ(κ)i|∑k(⟨~H⟩(κ)S(κ)T)jk|φ(κ)k⟩ (44)

where we again refer to App. B for details of the derivation. The above expression involves matrix elements of the mean-field operators of Eq. (43) multiplied by the transformation matrix.

The combination of Eq. (26), Eq. (42), and Eq. (44), which define the new equations of motion, necessitate toggling between the two SPF-SHF representations of Eq. (8) and Eq. (37). As will be further discussed below, this is achieved by the QR decomposition steps of Eq. (32) and Eq. (34).

The above equations have been obtained with the standard gauge condition for the original SPFs,

 ⟨˙φ(κ)i|φ(κ)j⟩=0∀i,j (45)

along with an additional gauge condition for the orthonormalized SHFs,

 ⟨˙~Ψ(κ)k|~Ψ(κ)l⟩=0∀k,l (46)

This additional gauge is equivalent to the condition , as follows from Eq. (36), and guarantees that the new SHFs remain orthonormal during the propagation. (Alternatively, the presence of an additional gauge condition can be taken to arise because of the QR decomposition of the coefficients according to Eqs. (31) and (32), which necessitates an additional gauge.HDM ())

In the following, we outline the implementation of the above equations Eq. (26), Eq. (42), and Eq. (44), as described in Refs. [Lubich, 2015; Kloss, Burghardt, and Lubich, 2017].

## V The Projector–Splitting Integration Scheme

From a general perspective, the projector–splitting integrator follows the idea of a second-order scheme which is known as Strang splitting Lubich (2008) in the mathematical literature. Well-known examples of this type of integrators in the physical sciences are the popular velocity-Verlet method in classical molecular dynamics Tuckerman (2010) or the second-order split-operator method in quantum dynamics.Feit, Fleck, and Steiger (1982) These algorithms, each in its own appropriate formalism, share the use of the symmetric Trotter expansion of the exponential, i.e., the approximation

 e(A+B)δ∼e12AδeBδ%e12Aδ (47)

These type of integrators have attractive general properties, including unitarity and the preservation of the underlying symplectic structure of the space in which the solution evolves. In our case, the approximated exponential is the formal solution of Eq. (12) over a short time interval which reads as the propagator in the MCTDH tensor–product space

 U(t)=exp(−ıℏP(Ψ)HΨt) (48)

and is approximated according to the chosen splitting scheme of the projector.

In light of the above, the algorithm described in Refs. [Lubich, 2015; Kloss, Burghardt, and Lubich, 2017] comprises three different steps:

1. a half-step for each of the modes, with the propagation of the SPFs and of the matrix (Eq. (44) and Eq. (42)),

2. a full step for the vector (Eq. 26),

3. another half-step for the SPFs and the matrix.

Before addressing these steps in further detail, we point out two key issues that need to be considered in the implementation of the algorithm:

First, in each integration step the regular SPFs need to be reconstructed from the propagated non-orthogonal SPFs and the transformation matrix . If this was done by an inversion of the matrix, according to Eq. (34), issues about ill–conditioning would arise, in exactly the same way as for the reduced density matrix. The present algorithm circumvents this problem by a QR decomposition of the propagated SPFs according to Eq. (34), . By definition, the resulting regular SPFs are orthogonal.

Second, the propagation of the coefficients and the SPFs has to be consistent, in the sense that the time evolution described by the coefficients is not “repeated” by the SPFs and vice versa. This property, which is visible in the standard MCTDH equations of motion, Eq. (27), in terms of the projector, is now encoded in the evolution. However, a complication arises from the fact that is not only updated according to Eq. (44), but also by QR decomposition of the coefficients, Eq. (31) and of the transformed SPFs Eq. (34). In particular, if the updated SPFs are generated by QR decomposition following propagation of the , using as explained above, the updated is not identical to as obtained by time propagation according to Eq. (44). This needs to be corrected for by additional (back-)propagation steps of serving as a “gauge correction”, as detailed below.

These considerations lead to an algorithm Lubich (2015); Kloss, Burghardt, and Lubich (2017) involving a simultaneous backward and forward-in-time propagation accompanied by two QR decompositions per step. The algorithm is schematically depicted in Fig. 2.

In detail, the first step is carried out as follows:

• QR decomposition of the vector and construction of the corresponding non–orthogonal SPFs (Eq. (34)) and mean–field potentials (Eq. (43)) at time .

• forward integration of Eq. (42) for a half time–step , to obtain new values of the non–orthogonal SPFs at time :

 Δ~φ(κ)=∫t0+δ/2t0dt′fφ(~φ(κ)(t′) | ⟨~H⟩(κ)(t0))

with referring to the r.h.s. of Eq. (42), where the left argument specifies the quantity which is integrated, while the right argument specifies the quantities that are held constant at a certain time (cf. Ref. [Beck and Meyer, 1997] for a similar notation in the context of the constant mean field (CMF) integrator of MCTDH).

• QR decomposition of the updated non–orthogonal SPFs (Eq. (34)) to obtain values of the orthonormal SPFs and of the matrix at time .

• backward integration of Eq. (44) for a half time–step to obtain a value of :

 ΔS(κ)=∫t0t0+δ/2dt′fS(S(κ)(t′) | ⟨~H⟩(κ)(t0),φ(κ)(t0+δ2))

where the same convention for the finite-time integration step was used as above, i.e., refers to the r.h.s. of Eq. (44). This is the “gauge correction” step referred to above: As a result of this step, a new initial value is generated, whose forward propagation according to Eq. (44) matches the combination obtained in the previous step. This value defines the correct initial condition, in the new gauge, for the coefficient in the next step.

• update of the vector at time according to Eq. (32), i.e., , where results from the preceding integration half-step while was held fixed.

This sequence is iterated for each mode, with .

In the second step of the algorithm, the updated vector is integrated for a full time step:

 ΔA=∫t0+δt0dt′fA(A(t′) | φ(κ)(t0+δ2))

where refers to the r.h.s. of Eq. (26).

The third step of the algorithm is logically equivalent to the first one, with the difference that the initial QR decomposition acts on the vector at time . For this reason, the order of the integration of the SPFs and of is reversed:

• QR decomposition of the vector (Eq. (31)) and construction of the corresponding mean–field potential (Eq. (43)) at time .

• backward integration of Eq. (44) for a half time–step towards a () value of (“gauge correction”):

 ΔS(κ) = ∫t0+δ2t0+δdt′fS(S(κ)(t′) | ⟨~H⟩(κ)(t0+δ)φ(κ)(t0+δ2))

The updated is subsequently used to construct non–orthogonal SPFs .

• forward integration of Eq. (42) for a half time–step, to obtain new values of the non–orthogonal SPFs at time :

 Δ~φ(κ)=∫t0+δt0+δ/2dt′fφ(~φ(κ)(t′) | ⟨~H⟩(κ)(t0+δ))
• QR decomposition of the updated non–orthogonal SPFs (Eq. (34)) to obtain values of the orthonormal SPFs and of the matrix at time , which is then used to update the vector (Eq. (31)).

Inspecting Fig. 2, it is evident that the first and second half steps have a symmetric structure. In fact, one can be obtained from the other by inverting the sequence of the operations. The direction of time integration, however, needs to be maintained unaltered so that both time steps globally result in a forward propagation of the SPFs.

The integration scheme bears analogies to the abovementioned CMF integrator of standard MCTDH,Beck et al. (2000); Beck and Meyer (1997) but obeys a time stepping scheme that is directly derived from the projector splitting of Eq. (41).

## Vi Summary and Conclusion

The aim of this article is to make some recent results obtained in the mathematics community more accessible to a chemical physics audience, specifically in the context of the new projector-splitting integrator for MCTDH developed by Lubich Lubich (2015) and recently implemented by Kloss et al.Kloss, Burghardt, and Lubich (2017) To this end, three aspects have been highlighted in the present work: First, the derivation of the modified MCTDH equations presented in Ref. [Lubich, 2015], from the perspective of a suitable splitting of the tangent space projection for MCTDH. Second, the concept of orthogonalizing the SHFs to formally eliminate the inverse of the single-particle density matrix from the equations of motion. Third, the structure of the algorithm designed by Lubich Lubich (2015); Kloss, Burghardt, and Lubich (2017) to make the new propagation scheme efficient.

The first two aspects can be deduced immediately from the tangent space projection of the MCTDH wavefunction Koch and Lubich (2010) (Sec. III), which is, however, not a common tool so far in the derivation of variational equations from the viewpoint of the chemical physics community. We believe that this perspective can be most useful in understanding the structure of the variational equations and designing new approximation schemes.

The third aspect emphasizes how the transformation is handled as an auxiliary quantity that is continuously updated during the algorithm, using QR decomposition steps to toggle between the different representations. As underscored in Ref. [Lubich, 2015], the algorithm does not use any pre-determined gauge as in the standard MCTDH formulation, but adapts the gauge via QR decompositions (or, alternatively, singular value decompositionsLubich (2015)). The back propagation steps for the matrix described above serve as gauge corrections, to ensure that the time evolution of the coefficients are SPFs are matched, as guaranteed by the projector in the original MCTDH equations.

The algorithm as described above generally leads to a robust propagation,Kloss, Burghardt, and Lubich (2017) and circumvents the regularization procedure of standard MCTDH. However, the QR decomposition steps are also affected by singularities of the density matrix, since the QR decomposition is non-unique for rank-deficient matrices and, hence, , see Eq. (38). Therefore, the propagation does have some dependence on how the QR algorithm treats the rank-deficient caseGolub and Van Loan (2013), see the discussion of Ref. [Kloss, Burghardt, and Lubich, 2017]. Whether or not the algorithm handles the initially unoccupied SPFs in an advantageous way, especially as compared with the construction of optimal unoccupied SPFs,Manthe (2015) is currently a matter of debateMey () and needs to be further investigated in numerical studies.

## Vii Acknowledgments

We thank the German-Israeli Foundation for Scientific Research and Development for support of this project under grant number GIF I-1337-302.5/2016. M.B. gratefully acknowledges fellowship support by the Alexander von Humboldt Foundation. We thank H.-D. Meyer for helpful discussions and suggestions in the context of the integration scheme.

## Appendix A Basic notions of tensor algebra

In tensor language, the MCTDH ansatz is known as Tucker formatKoch and Lubich (2010); Lubich et al. (2013); Bachmayr, Schneider, and Uschmajew (2016), in which is decomposed in terms of a core–tensor (the –vector of conventional MCTDH) and a rectangular matrix per physical dimension (the matrix of Eq. (3)) representing the change of basis from the primitive grid to the SPFs. Other alternative decompositions have been explored both in the chemical physics and mathematical literature, notably relating to the ML-MCTDH scheme that corresponds to a hierarchical Tucker decomposition. We refer to Ref. [Bachmayr, Schneider, and Uschmajew, 2016] for a topical survey of the field.

To write the standard tensor operations in concise form, the formalism which is used in Refs. [Lubich, 2015; Kloss, Burghardt, and Lubich, 2017] makes extensive use of the concepts of matricisation and tensorisation, which are a key point in translating the equations to common MCTDH notation.

We illustrate these concepts for operators in a sum-of-products (SOP) form, where the application of the operator can be split into the sequential application of smaller matrices, as has been recognized early on in the context of MCTDH Manthe, Meyer, and Cederbaum (1992); Beck et al. (2000).

In the tensor formalism, this type of operation is explicitly written by means of the matricisation of the tensor. When an operator with matrix representation acting specifically on mode is applied to a tensor , it is convenient to recast the tensor in matrix form such that its rows are labeled by the index of the primitive basis of mode . In formulas, the -mode matricisation is defined as

 A(κ)jκJκ=AJ (49)

with the definition of multi–indices as described in the main text. With , the action of the operator can be computed using conventional matrix multiplication,

 A(κ)′=O(κ)A(κ) (50)

The tensor corresponding to the resulting matrix can be reconstructed by tensorization, i.e. by reordering the components and labeling them by the usual multi-index .

The sequential operations consisting in (i) -mode matricisation (ii) matrix multiplication with a -mode matrix and (iii) tensorization are written concisely as

 A′=A×κO(κ) (51)

or, with a full specification of indices,

 A(κ)′Jκl=∑jκO(κ)l,jκAJ (52)

By repeated application of this definition, the action of an operator in product form, , is written as

 A′=A×fκ=1O=A×1O(1)×2O(2)⋯×fO(f) (53)

which again can be explicitly written with a full specification of the indices as

 A′L=∑j1…∑jfO(1)l1,j1…O(f)lf,jfAJ (54)

The reader should note that the operations defined above can be used to define not only the operator matrix elements, but also the tensor decomposition in Tucker form. In fact, by using the definition of Eq. (53), the MCTDH ansatz in tensor form, i.e., Eq. (3) can be immediately recognized as

 Ψ=A×fκ=1U(κ) (55)

where is the core–tensor and the representations of the SPFs in the primitive basis.

By a full specification of the indices and by removing the primitive basis projection, all equations of Refs. [Lubich, 2015; Kloss, Burghardt, and Lubich, 2017] can be cast in standard MCTDH form. We consider as another instance the definition of the SHFs

 Ψ(κ)=A(κ)f⨂n≠κU(κ)T (56)

The matrix has dimensions ( being the dimension of the -th primitive grid) and its elements are given by

 Ψ(κ)lIκ=∑JκA(κ)lJnf∏n≠κU(n)injn (57)

where we have substituted the expression for the Kronecker matrix product. Again, by substituting the expression for the SPF expansion on the primitive basis we obtain

 Ψ(κ)lIκ=∑JκAJκl∏n≠κ⟨χ(n)in|φ(n)jn⟩ (58)

which is recognized as the primitive basis representation of the -th SHF for mode

 Ψ(κ)lIκ=⟨χ(1)i1…χ(κ−1)iκ−1χ(κ+1)iκ+1…|ψ(κ)l⟩ (59)

## Appendix B Projector splitting equations of motion

Here, we provide details of the derivation of the equations of motion Eq. (43) and Eq. (44) for the quantities and . By analogy with Sec. III.C, we consider the following equations for the projectors :

 P±κ(Ψ)[˙Ψ−1ıℏHΨ]=0 (60)

where is conveniently expressed in the form Eq. (37) such that

 ˙Ψ−1ıℏHΨ=∑l~φ(κ)l˙~Ψ(κ)l+∑l˙~φ(κ)l~Ψ(κ)l−1ıℏ∑lH~φ(κ)l~Ψ(κ)l (61)

From the explicit form of the projector in Eq. (40a) in conjunction with the additional gauge condition Eq. (46) for the orthogonalized SHFs, we obtain

 P+κ(Ψ) ∑l(~φ(κ)l˙~Ψ(κ)l+˙~φ(κ)l~Ψ(κ)l−1ıℏH~φ(κ)l~Ψ(κ)l) = ∑j~Ψ(κ)j(˙~φ(κ)j−1ıℏ∑l⟨~H⟩(κ)jl~φ(κ)l)=0

with the mean-field Hamiltonian matrix element . Eq. (B) immediately yields Eq. (43).

Turning to the projector in Eq. (40b), we obtain

 P−κ(Ψ) ∑l(~φ(κ)l˙~Ψ(κ)l+˙~φ(κ)l~Ψ(κ)l−1ıℏH<