Axel U J Lode

Camille Lévêque

Lars Bojer Madsen

Alexej I Streltsov

Ofir E Alon

Colloquium: Multiconfigurational time-dependent Hartree approaches for indistinguishable particles

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

This Colloquium introduces and discusses the development and capabilities of the Multiconfigurational Time-Dependent Hartree (MCTDH) approaches (Beck et al., 2000;Manthe et al., 1992;Meyer et al., 1990) for the dynamics of indistinguishable particles, MCTDH-F for Fermions (Caillat et al., 2005;[START_REF] Kato | [END_REF]Zanghellini et al., 2003) and MCTDH-B for Bosons (Alon et al., 2008b;Streltsov et al., 2007a) or, together, MCTDH-X (Alon et al., 2007b). To introduce and motivate MCTDH-X, we give an account of the theoretical development that lead to its formulation. We illustrate the insight into many-body physics gained thus far from applications of MCTDH-X in the areas of atomic, molecular, and optical physics with applications to realworld, experimentally-realized examples of the dynamics of atoms in trapped Bose-Einstein condensates and of electrons in atoms and molecules. Finally, theoretical and numerical developments, prospects, and possible future avenues of the MCTDH-X approaches are outlined.

The time-dependent many-body Schrödinger equation (TDSE) for interacting, indistinguishable particles is a cornerstone of many areas of physics. Exactly solvable models are very scarce for, both, the timedependent (Fasshauer and Lode, 2016;[START_REF] Lode | Tunneling Dynamics in Open Ultracold Bosonic Systems[END_REF]Lode et al., 2012a) and the time-independent Schrödinger equation (Calogero, 1969;Dukelsky and Schuck, 2001;Girardeau, 1960;Lieb, 1963;Lieb and Liniger, 1963;McGuire, 1964;Sutherland, 1971;Yukalov and Girardeau, 2005) and could so far not be generalized to real-world problems -a numerical approach to tackle the TDSE is therefore needed. A direct numerical solution becomes, however, impracticable -especially for inhomogeneous systems -already for very few particles, because the Hilbert space in which the solution of the TDSE lives grows exponentially with the number of particles considered. Generally, approximations are therefore needed; the problem has to be represented accurately enough to cover the physical properties of the many-body state while -at the same time -the chosen representation has to be compact enough to be manageable computationally. Since the TDSE is so fundamental, there exist many approximations to its solution, with each new methodology being a step in the quest for an ever more accurate description. We name here as examples the (multi-orbital) mean-field (Alon et al., 2007a;Gross, 1961;McLachlan and Ball, 1964;Pitaevskii, 1961) and the configuration interaction [START_REF] Cramer | Essentials of computational chemistry : theories and models[END_REF]David Sherrill and Schaefer, 1999;[START_REF] Jensen | Introduction to Computational Chemistry[END_REF][START_REF] Szabo | Modern quantum chemistry : introduction to advanced electronic structure theory[END_REF] approaches. Mean-field approaches, however, drop some or all of the correlations from the description of the manybody state and configuration interaction or exact diagonalization is restricted to situations, where the initially chosen, static basis remains suitable for all times (Lode et al., 2012a). For the case that Hubbard models are considered, there exist, for instance, the (time-dependent) density renormalization group (Schollwöck, 2005), matrix product states (Schollwöck, 2011), and time-evolved block decimation methods (Vidal, 2004). These latter methods describe correlated many-body dynamics for Hubbard models, but are not directly applicable other cases. To obtain a method to describe correlations in the dynamics of many-body systems, i.e., for the formulation of the MCTDH for fermions (Caillat et al., 2005;[START_REF] Kato | [END_REF]Zanghellini et al., 2003) and the MCTDH for bosons (Alon et al., 2008b;Streltsov et al., 2007a), or, for short, MCTDH-X (Alon et al., 2007b), two basic ingredients were needed: (i), a unification of the time-independent basis of configuration interaction with the time-adaptive ansatz of the (multiorbital) mean-field (or time-dependent Hartree-Fock or self-consistent field methods) for indistinguishable particles and, (ii), an appropriate time-dependent variational principle [START_REF] Dirac | Proc. R. Soc. Lond. A[END_REF]Frenkel, 1934;[START_REF] Kramer | Geometry of the Time-Dependent Variational Principle in Quantum Mechanics[END_REF]McLachlan, 1964).

MCTDH-X is a general method for the solution of the TDSE for interacting indistinguishable particles that yields a well-controlled error (Fasshauer and Lode, 2016;[START_REF] Lode | Tunneling Dynamics in Open Ultracold Bosonic Systems[END_REF]Lode et al., 2012a); its development, theory, applications, and generalizations constitute the main subject of this Colloquium.

In Sec. II, we provide a unified formulation of the equations of motion (EOM) of MCTDH-X for complete as well as restricted configuration space, i.e., for situations where all or only part of the possible Slater determinants or permanents are included in the description, respectively. For the sake of simplicity and instructivity, we restrict our discussion of MCTDH-X with a restricted space to the so-called restricted active space approach. In Sec. II.C, we conclude our exhibition of the MCTDH-X approaches with benchmarks with an exactly solvable model problem, the harmonic interaction model, that show that the method is in principle exact (Fasshauer and Lode, 2016;[START_REF] Lode | Tunneling Dynamics in Open Ultracold Bosonic Systems[END_REF]Lode et al., 2012a).

In Sec. III, we focus on MCTDH-B applications to the physics of quantum correlations and fluctuations and the variance of operators in Bose-Einstein condensates (BECs). We summarize an illustrative application of MCTDH-B to the dynamics of a BEC subject to time-dependent interparticle interactions where MCTDH-B computations were directly compared to experiment (Nguyen et al., 2019). Moreover, we highlight some insight into the intriguing physics of the variances of observables in the so-called infinite-particle-numberlimit (Alon and Cederbaum, 2018;Klaiman and Alon, 2015) that were obtained with the help of MCTDH-B. In Sec. IV, we discuss some insights that MCTDH-F enabled for the correlated dynamics of electrons in atoms and molecules. We give an account of work with MCTDH-F with a focus on studies of photoionization cross-sections and time delays that were experimentally verified (Haxton et al., 2012;Omiste and Madsen, 2018).

We conclude in Sec. V with an overview over current theoretical progress with MCTDH-X as well as possible future avenues of method development. We discuss the key ideas of the multilayer (Cao et al., 2017(Cao et al., , 2013;;Manthe and Weike, 2017;Wang and Thoss, 2009) and other generalizations of MCTDH-B (Alon et al., 2014;Grond et al., 2013b) and MCTDH-F (Lötstedt et al., 2019b;Sato and Ishikawa, 2015;[START_REF] Sawada | [END_REF] as well as orbital adaptive time-dependent coupled cluster theories (Kvaal, 2012(Kvaal, , 2013;;Pedersen and Kvaal, 2019;Sato et al., 2018a,b).

Our Colloquium thus gives an overview of the activities in the community that develops and applies MCTDH-X. Achievements made using MCTDH-X on ultracold atoms in BECs and on the correlated dynamics of electrons in atoms and molecules are illustrated and the state-of-theart development on the theory and its generalizations are introduced.

II. MCTDH-X THEORY

To obtain the MCTDH-X equations, one applies a variational principle to the TDSE with a parametrized ansatz. As Kramer and Saraceno aptly assessed [ [START_REF] Kramer | Geometry of the Time-Dependent Variational Principle in Quantum Mechanics[END_REF], p.6]: "As is well-known, a variational principle is a blind and dumb procedure that always provides an answer, but its accuracy depends crucially on the choice of the trial function." Different types of ansatzes thus lead to approximations with different qualitative behavior. Generally, the MCTDH-X type of ansatz is a time-dependent linear combination of a set of fully symmetrized or antisymmetrized products of time-dependent single-particle states or orbitals, the so-called configurations. So, why is the MCTDH-X-ansatz for the wavefunction a good ansatz? One, the time-dependent configurations in the MCTDH-X ansatz are an in-principle complete basis of N -particle Hilbert space and, two, they are constructed such that they are strictly ortho-normalized at any time. These two properties, in combination with the timedependent variational principle, allow to infer the convergence of the method: if a sufficiently large set of configurations has been included in a computation, i.e., the result remains identical when more configurations are included, one can conclude that the employed ansatz spans a sufficiently large portion of N -particle Hilbert space.

Here, we will discuss the archetypical MCTDH-X theory with an ansatz (Alon et al., 2007b(Alon et al., , 2008b;;Caillat et al., 2005) including all possible configurations of N particles in M orbitals. We will also cover the formulation of MCTDH-X with an ansatz obtained with a further truncation of Hilbert space via the restricted active space (RAS) approach (Olsen et al., 1988) as put forward in Refs. (Lévêque and Madsen, 2017;Miyagi and Madsen, 2013). We note that the RAS approach originates from ground-state quantum chemistry, but -although physical insight into the emergent quantum dynamics may help to choose a sensible RAS scheme -it may not be the best choice for the emergent dynamics of many-body systems. The EOM of MCTDH-X for completely general configuration spaces -of which the RAS is a special case -have been put forward for a single kind of indistinguishable particles in Ref. (Haxton and McCurdy, 2015) and even for multiple species of indistinguishable particles in Ref. (Anzaki et al., 2017). We chose to present the specialized RAS truncation scheme for MCTDH-X in this Colloquium, because applications of them exist for both fermions and bosons. Moreover, as a truncation scheme we find the construction of the RAS instructive, illustrative, and simple, while the obtained EOM hint at some of the changes triggered by the truncation of the Hilbert space in comparison to the standard MCTDH-X with a complete configuration space.

Moreover, as common for ultracold atoms and electron/nuclear dynamics, we focus on Hamiltonians of the form:

Ĥ = N i=1 ĥ(r i ; t) + N i<j W (r i , r j ; t).
(1)

Here, the positions of the k-th particle is denoted by r k , ĥ(r; t) is a general, possibly time-dependent, singleparticle operator and Ŵ (r, r ; t) is a general, possibly time-dependent, two-particle operator.

A. Unified equations of motion

We now discuss the EOM of MCTDH-X and their derivation for the case where all possible configurations of N particles in M time-dependent orbitals are included in the ansatz,

P Q (a) 3,1,5,2,1 > | 2,2,3,3,2 > | 2,1,1,1,1 > | -- 1,2,0,2,1 > | - (
|Ψ F CI = n C n (t)| n; t ; n = (n 1 , ..., n M ) T ; | n; t = N M i=1 b † i (t) ni |vac (2)
Here, the normalization N is

1 √ M i=1 ni! ( 1 √ N ! )
for bosons (fermions), the number of particles N is considered constant, N = i n i , and b † j (t) creates a particle in the single-particle state Φ j (r; t),

Φ j (r, t) = r| b † j (t)|vac . (3)
Here, and in the following, we use the symbol r to summarize the degrees of freedom (spin and space) of the orbitals. The coefficients,

C n (t) = n|Ψ F CI , (4)
are the complex time-dependent weights of each configuration's contribution to the many-body state |Ψ F CI .

Here, and in the following, we drop the dependence on time for notational convenience. For bosons, there are

N +M -1 N
coefficients and for fermions, there are M N coefficients. To obtain the EOM, one can apply the timedependent variational principle [START_REF] Kramer | Geometry of the Time-Dependent Variational Principle in Quantum Mechanics[END_REF] for the TDSE,

Ĥ|Ψ = i∂ t |Ψ , (5)
and use |Ψ F CI as an ansatz. The action reads:

S = dt Ψ F CI | Ĥ -i∂ t |Ψ F CI - kj µ kj (t) [ Φ k |Φ j -δ kj ] . (6)
We add the Lagrange multipliers µ kj (t) to ensure the ortho-normalization of the single-particle states ( Φ k |Φ j = δ kj ) at any time. We demand, independently, the stationarity of S with respect to variations of the orbitals {Φ i (r, t)} and the coefficients {C n (t)},

δS[{Φ i (r; t)}, {C n (t)}] δΦ * i (r; t) ! = 0, δS[{Φ i (r; t)}, {C n (t)}] δC * n (t) ! = 0. ( 7
)
After a straightforward derivation (Alon et al., 2007b(Alon et al., , 2008b;;Caillat et al., 2005) we arrive at a coupled set of non-linear coupled integro-differential EOM for the orbitals,

i∂ t |Φ j = Q ĥ|Φ j + M k,s,q,l=1 {ρ} -1 jk ρ kslq Ŵsl (r; t)|Φ q , Q = 1 - i |Φ i Φ i | (8)
In our derivation of this EOM we have, for convenience, set the gauge that removes the ambiguity in the choice of the orbitals (Alon et al., 2007b(Alon et al., , 2008b;;Manthe et al., 1992;Meyer et al., 1990) to be

η ij = Φ i |∂ t Φ j = 0; ∀i, j ∈ 1, ..., M. (9)
Other choices for η ij are possible (Beck et al., 2000;Caillat et al., 2005) and provide some flexibility in designing the numerical approaches for the time-integration of the EOM, like splitting and regularization methods (Kloss et al., 2017;Koch et al., 2013;Lubich and Oseledets, 2014;Lubich et al., 2018;Meyer and Wang, 2018).

Here, we used the matrix elements of the reduced onebody and two-body density matrices,

ρ kq = Ψ| b † k bq |Ψ , (10)
ρ kslq = Ψ| b † k b † s bq bl |Ψ , (11)
respectively. Since these matrix elements, ρ kq , ρ ksql are functions of the coefficients in the ansatz, Eq. (2), the orbitals' time-evolution is explicitly dependent on the coefficients' time-evolution. The projector Q in the EOM emerges as a result of the elimination of the Lagrange multipliers µ kj in the action S [Eq. ( 6)]; it is therefore a direct consequence of the ortho-normalization of the orbitals Φ j (r; t) at any time. We further defined the local interaction potentials, Ŵsl (r; t) = Φ * s (r ; t) Ŵ (r, r ; t)Φ l (r ; t)dr .

(12)

The EOM for the coefficients [Eq. ( 4)] form a linear set of equations,

i∂ t C n (t) = n n; t| Ĥ| n ; t C n , (13)
which is coupled to the orbital's EOM [Eq. ( 8)] as the expectation value n; t| Ĥ| n ; t is a function of the orbitals, as can easily be understood by expressing the Hamiltonian in second quantized notation:

Ĥ = M k,q=1 h kq b † k bq + M k,s,q,l=1 W ksql b † k b † s bl bq . ( 14
)
Here, we used the matrix elements of the one-and twobody Hamiltonian,

h kq = Φ k | ĥ(r i ; t)|Φ q , (15)
W ksql = dr dr Φ k (r; t)Φ s (r ; t)× (16)
W (r, r ; t)Φ q (r ; t)Φ l (r ; t)

respectively. The Hamiltonian, Eq. ( 14), is a function of h kq , W ksql that are, in turn, functions of the orbitals Φ k (r; t). Therefore, the coefficients' time-evolution, governed by the EOM (8), also directly depends on the orbitals. The EOM of the MCTDH-X method, Eqs. ( 8) and (13), are thus coupled.

B. Restricted spaces

Configurations can be removed from the full set employed in the ansatz |Ψ F CI for the wavefunction that was used in the derivation of the MCTDH-X EOM, Eqs. ( 8) and ( 13). This restriction of the configuration space reduces the numerical effort and may thus enable computations for cases where the number of terms in the ansatz |Ψ F CI is intractably large. Moreover, the changes in the emergent dynamics triggered by the restriction of the configuration space may lead to a physical insight into what parts of Hilbert space are explored by the many-body state.

General restrictions to the configuration space are possible and lead to general MCTDH-X EOM that are discussed, for instance, in Refs. (Anzaki et al., 2017;Haxton and McCurdy, 2015). It is important to stress here that the MacLachlan (McLachlan, 1964) and Lagrangian [START_REF] Kramer | Geometry of the Time-Dependent Variational Principle in Quantum Mechanics[END_REF] variational principles, as well as their union, the Dirac-Frenkel variational principle [START_REF] Dirac | [END_REF]Frenkel, 1934), lead to the same unified MCTDH-X EOM only in the case that the ansatz for the wavefunction contains all possible configurations, i.e., as given in Eq. (2). For general ansatzes with a restricted set of configurations, however, the McLachlan and Lagrangian variational principles are inequivalent (Haxton and McCurdy, 2015).

Here, we focus on the restricted active space (RAS) approach for the restrictions of the configuration space (Olsen et al., 1988) of MCTDH-X, because we find its strategy for the construction of many-body Hilbert space instructive and suitable to illustrate the changes that arise when one deals with a truncated configuration space. Moreover, there are applications of the RAS approach in combination with MCTDH-X for, both, bosons and fermions. These methods are referred to as time-dependent RAS self-consistent-field (TD-RASSCF) for fermions (Miyagi andMadsen, 2013, 2014a,b) (TD-RASSCF-F) and TD-RASSCF-B for bosons (Lévêque and Madsen, 2017;Léveque and Madsen, 2018). For the sake of clarity and coherence of presentation in our Colloquium, we will refer to TD-RASSCF-F and TD-RASSCF-B as RAS-MCTDH-B and RAS-MCTDH-F, respectively, and RAS-MCTDH-X, together.

In the original formulation of the RAS-MCTDH-F (Miyagi andMadsen, 2013, 2014a), three subspaces of adaptive orbitals were considered: P 0 , P 1 , and P 2 with M 0 frozen, M 1 unrestricted, and M 2 restricted occupations, respectively. Since the P 0 space with orbitals with frozen occupations can not be meaningfully assigned for bosons, we limit ourselves here to the case of RAS-MCTDH-X with two active subspaces -P 1 and P 2 -to restrict the number of configurations with a total number M = M 1 + M 2 orbitals. The number of orbitals in the P 1 subspace must be large enough to accommodate all the particles, i.e., one configuration, at least, has no particles in the P 2 subspace. For bosons, the P 1 subspace includes at least one orbital and for fermions M 1 ≥ N holds. The restriction on the configuration space follows from specifying a maximum number of particles, N max , that can occupy the P 2 subspace. The ansatz for the RAS-MCTDH-X method reads,

|Ψ RAS = n∈V C n (t)| n, t , (17)
where the configurations span the space V that is obtained by restricting the total configurational space of Eq. (2) using the RAS determined through the parameters M 1 , M 2 , N max . The RAS-MCTDH-X wavefunction can be seen as a bridge between the mean-field approaches, TD-Hartree-Fock for fermions and TD-Gross-Pitaevskii for bosons, and the MCTDH-X approach, which are all limiting cases of the RAS-MCTDH-X ansatz. The EOM for the set of time-dependent coefficients {C n (t)} and orbitals {|Φ i (t) } M i=1 are derived following the recipes of the MCTDH-X framework, see Sec. II.A, albeit, here, with a real (Lagrangian) action functional (Lévêque and Madsen, 2017;Miyagi andMadsen, 2013, 2014b). A set of equations for the coefficients and the orbitals is obtained:

i∂ t C n = ij (h ij -iη ij ) Φ n | b † i bj |Ψ + 1 2 ijkl W ikjl Φ n | b † i b † k bl bj |Ψ , (18)
and

i Q[∂ t |Φ j ] = Q ĥ|Φ j + M k,s,q,l=1
{ρ} -1 jk ρ kslq Ŵsl (r; t)|Φ q ,

(19) respectively. The set of equations for the orbitals is similar to the one obtained for MCTDH-X, see Eq. ( 8), except that the projector Q appears on both sides of Eq. ( 19) and the set of equations for the coefficients includes an additional term, namely, η ij = Φ i | Φj . This gauge freedom in the MCTDH-X equations, set to zero, cannot be chosen arbitrary to simplify the equations any more, because the P 1 and P 2 orbitals are not equivalent and the transformation of the orbitals from one subspace to another have to be taken into account explicitly. Thus, for each pair of orbitals {i , j }, with i ∈ P 1 and j ∈ P 2 , the matrix element η i j is evaluated via an additional set of equations. The choice of the excitation scheme to promote particles from the P 1 to the P 2 subspace plays an important role to simplify the evaluation of η i j . Here, we present the case of the so-called general excitation scheme (Miyagi and Madsen, 2014b), where all successive occupation numbers 0, ..., N max of the P 2 subspace are considered. The matrix elements η i j are evaluated from,

k l (iη k l -h k l )ζ l j k i = 1 2 klmn W kmln ζ lnj kmi , (20)
where the fourth-and sixth-order tensors are defined by

ζ l j k i = Ψ| b † i bj ( 1 -Π) b † k bl |Ψ (21) ζ lnj kmi = Ψ| b † i bj ( 1 -Π) b † k b † m bn bl |Ψ , (22)
with Π = n∈V | n, t n, t| being the projector onto the RAS configurational space. The time-derivative of the orbitals can be expressed as

∂ t |Φ j = i η ij |Φ i + Q[∂ t |Φ j ].
The η ij |Φ i -term describes the transformation of the P 1 and P 2 orbitals into each other, and the Q[∂ t |Φ j term describes the extension of the orbitals into the space not spanned by the orbitals at time t. From Eqs. (20)

... ... ...

V 0 V 1 V 2 V = ... Bosons Q P 1 P 2 (a) Fermions (b) FIG.
2 Illustration of restricted active space schemes for the restriction of configuration spaces for bosons [(a)] and fermions [(b)]. The space of active orbitals is partitioned into two sets, P1,2 and the space of virtual orbitals is denoted by Q. All possible configurations of N particles in the M1 orbitals of the P1 space are considered in the ansatz of RAS-MCTDH-X. For the P2 space, a maximal occupation of all M2 orbitals together is fixed to be Nmax. Thus, in RAS-MCTDH-X, of all possible configurations of N particles in P1 ⊕P2, those configurations where there are more than Nmax particles in P2 are dropped. The total Hilbert space V spanned by the ansatz of RAS-MCTDH-X ( 17) is a direct sum of the sub-spaces that contain 0, ..., Nmax particles (blue boxes). and ( 19) the time-derivative of the orbitals can be evaluated, and from Eq. ( 18) the time-derivative of the coefficients are available after solving Eq. ( 20) for the matrix elements η i j (t). The restriction of the configuration space thus leads to more complicated EOM, but the (drastic) reduction of the number of configurations enables faster or in some situations more accurate descriptions of many-body systems than MCTDH-X. Note that the EOM for other RAS-excitation-schemes can be found in (Lévêque and Madsen, 2017;Miyagi andMadsen, 2013, 2014b). For the so-called complete active space approach with an additional space hosting orbitals with occupations that are fixed, see (Sato and Ishikawa, 2015).

C. Benchmarks with an exactly solvable model

Since the introduction of MCTDH-B and MCTDH-F, many benchmarks of the predictions of these approaches have been performed. Most of these benchmarks consist in a comparison of the predictions of the MCTDH-X approaches to other theoretical approaches like, for instance, exact diagonalization with a timeindependent one-particle basis set. Example benchmarks against other approaches include the ionization of 4 He (Hochstuhl and Bonitz, 2011) or the photoionization of 9 Be (Haxton et al., 2011) in the case of MCTDH-F or a comparison with the Bose-Hubbard model (Sakmann et al., 2009;Streltsov et al., 2006) in the case of MCTDH-B. We note that the interesting MCTDH-X applications are those cases, where diagonalization is no longer affordable numerically. Here, we focus on available benchmarks of MCTDH-X with exactly solvable models, specifically, on the harmonic interaction model (HIM) (Armstrong et al., 2011;Cohen and Lee, 1985;Gajda, 2006;Yan, 2003;Za luska-Kotur et al., 2000) that describes N indistinguishable harmonically-trapped particles interacting via a harmonic interaction potential that is proportional to the square of their distance. The HIM has the unique feature that it straightforwardly can be generalized to include time-dependence in the harmonic trapping of and the harmonic interactions between particles while remaining exactly solvable (Fasshauer and Lode, 2016;[START_REF] Lode | Tunneling Dynamics in Open Ultracold Bosonic Systems[END_REF]Lode et al., 2012a). This time-dependent HIM (TD-HIM) is a well-suited test case for MCTDH-X, because it represents one of the rare cases where a numerically-exact solution to the TDSE for a correlated problem with a time-dependent Hamiltonian can be obtained via a mapping to a time-dependent one-body Schrödinger equation that can be integrated numerically at any desired level of accuracy.

The Hamiltonian of the TD-HIM reads

Ĥ (t) = N i=1 (- 1 2 ∂ 2 r + 1 2 ω T D (t) 2 r 2 )+K T D (t) N i<j (r i -r j ) 2 , (23)
where the time-dependent trap frequency, ω T D , and the time-dependent interaction strength, K T D , are given by:

ω T D (t) = ω [1 + f (t)] ; K T D (t) = K 1 - ω 2 0 2N K f (t) .
(24) We compare solutions of the TDSE with this Hamiltonian to (RAS-)MCTDH-B ones in Fig. 3.

The convergence of (RAS-)MCTDH-B towards the exact result for an increasing number of variational parameters in the wavefunction is demonstrated by the results in Fig. 3 for N = 10 bosons; for a demonstration with fermions and MCTDH-F see (Fasshauer and Lode, 2016).

III. MCTDH-B AND BOSE-EINSTEIN CONDENSATES

For the sake of brevity, we restrict our discussion here to the quantum dynamics obtained with MCTDH-B modelling an experiment with a quasi-one-dimensional BEC subject to a time-dependent interparticle interaction in Sec. III.B as well as to the appealing many-body physics in the variance of observables in Sec. III.C. Before turning to these applications of MCTDH-B, we introduce the relevant quantities of interest. 23) and ( 24). The time-dependent center-of-mass energy (t) is plotted in comparison to the exact values in the lower panel for different particle and orbital numbers. A convergence with an increasing number of orbitals, i.e., amount of variational parameters in the (RAS-)MCTDH-B wavefunction, is observed. See (Lode et al., 2012a) for details on (t).

A. Analyzing many-body states of bosons

The key insights that MCTDH-B has to offer are due to the fact that it is a wavefunction-based approach: from the approximate solution |Ψ to the TDSE, correlations and coherence can be quantified, for instance, using reduced density matrices and their eigenvalues (Sakmann et al., 2008):

ρ (p) (r 1 , ..., r p , r 1 , ..., r p ; t) = Tr p+1,...,N [|Ψ Ψ|] (25) = N ! (N -p)! dr p+1 • • • dr N Ψ * (r 1 , ..., r p , r p+1 , ..., r N ; t) ×Ψ(r 1 , ..., r p , r p+1 , ..., r N ; t).
The diagonal of the p-th order density matrix, i.e., ρ (p) (r 1 , ..., r p , r 1 , ..., r p ; t) ≡ ρ (p) (r 1 , ..., r p ; t), is the probability to find particles 1, ..., p at positions r 1 , ..., r p , respectively, and is referred to as the p-body density. In the case of p = 1, by convention, one drops the (1)superscript and speaks of just the density ρ(r), i.e., ρ(r; t) ≡ ρ (1) (r; t) ≡ ρ (1) (r, r = r; t) is implied. In this subsection, we present observables like ρ (p) derived using the wavefunction Ψ in position space; the equations are, however, also valid for momentum space analogons of the observables when the wavefunction in mo-mentum space is used and r is replaced by k. The offdiagonal part of the p-th order reduced density matrix, ρ (p) (r 1 = r 1 , ..., r p = r p , r 1 , ..., r p ; t), determines the pth order coherence. To further quantify the p-th order coherence, the p-th order Glauber correlation function,

g (p) (r 1 , ..., r p , r 1 , ..., r p ; t) = ρ (p) (r 1 , ..., r p , r 1 , ..., r p ; t) p k=1 ρ (1) (r k ; t)ρ (1) (r k ; t) , (26)
is a good measure. Essentially, g (p) gives a spatially resolved picture of the representability of the p-th order density matrix by a product of one-body densities:

g (p)
= 1 implies that the p-body density cannot be represented by a product of one-body densities. In this g (p) = 1 case, therefore, the many-body state contains quantum correlations (of p-th order). Such quantum correlations entail fluctuations of observables and can be probed (experimentally) with single-shot images or via the variance of operators (see below).

One important correlation effect that has been discussed in many works applying MCTDH-B is fragmentation (Mueller et al., 2006;Nozières and Saint James, 1982a;Spekkens and Sipe, 1999), i.e., the situation when the reduced one-body density matrix ρ (1) (r, r ; t) of interacting bosons acquires several macroscopic eigenvalues, see for instance [START_REF] Lode | Tunneling Dynamics in Open Ultracold Bosonic Systems[END_REF][START_REF] Lode | [END_REF]Lode and Bruder, 2017;Lode et al., 2012b;[START_REF] Sakmann | Many-Body Schrödinger Dynamics of Bose-Einstein Condensates[END_REF]Sakmann et al., 2009Sakmann et al., , 2010;;Streltsov et al., 2008Streltsov et al., , 2009Streltsov et al., , 2011)). If ρ (1) (r, r ; t) has only one single significant eigenvalue, then the state is referred to as condensed (Penrose and Onsager, 1956).

To discuss fragmentation and condensation, we thus write ρ (1) (r, r ; t) using its eigenvalues n

(1) i and its eigenfunctions Φ

(N O) i (r; t): ρ (1) (r, r ; t) = i n (1) i Φ (N O), * i (r ; t)Φ (N O) i (r; t). (27)
We note that the n

(1) i are nothing but the eigenvalues of the matrix elements ρ kq in Eq. ( 10). In practice, the n

(1) i are therefore computed by straightforwardly diagonalizing the M × M matrix ρ kq . Analogously, the eigenvalues n

(2) i of the two-body density ρ (2) are available via the diagonalization of ρ ksql .

In cold-atom experiments, the standard measurement is absorption images. Such single-shot images correspond to a projective measurement of the many-body state |Ψ (Castin and Dalibard, 1997;Dziarmaga et al., 2003;Javanainen and Yoo, 1996;Sakmann and Kasevich, 2016). In the ideal case, each image contains information about the position and momentum of every particle. Each measurement thus corresponds to a random sample s k of positions that is distributed according to the N -body probability distribution P (r 1 , ..., r N

; t) = ρ (N ) (r 1 , ..., r N ; t) = |Ψ(r 1 , ..., r N ; t)| 2 : s k = {s k 1 , ..., s k N } ∼ |Ψ(r 1 , ..., r N ; t)| 2 (28)
To directly model these images with a wavefunction computed by MCTDH-X, one has to draw random samples from the N -body density, i.e., compute a set of socalled single-shot simulations s k , k = 1, ..., N shots . The numerical difficulty in sampling high-dimensional probability distributions can be overcome by factorizing the N -particle probability into a set of conditional probabilities,

P (r 1 , ..., r N ; t) =P (r 1 ; t)P (r 2 |r 1 ; t) × • • • (29) • • • × P (r N |r 1 , ..., r N -1 ; t).
To obtain a simulation s = (s 1 , ..., s N ) of a single-shot, the first particle's position s 1 is drawn from the one-body density

s 1 ∼ P (r; t) = ρ(r; t) = Ψ| Ψ † (r; t) Ψ(r; t)|Ψ . ( 30
)
Here

Ψ(r) = M j=1 bj Φ j (r; t) [ Ψ † (r) = M j=1
b † j Φ * j (r; t)] is the standard field annihilation [creation] operator. The second particle's position, s 2 , is then sampled from the conditional probability that is computed from a reduced many-body state, Ψ (1) where a particle has been detected at s 1 ,

s 2 ∼ P (r 2 |s 1 ; t) = Ψ (1) | Ψ † (r; t) Ψ(r; t)|Ψ (1) , |Ψ (1) = N Ψ(s 1 )|Ψ , (31)
where N represents the normalization constant. This procedure is continued until all particles have been detected at positions s 1 , ..., s N and the single-shot image, i.e., the vector of positions s = (s 1 , ..., s N ) is obtained.

In principle, all information about the N -body density ρ (N ) (r 1 , ..., r N ; t) can be extracted from single-shot images.

We now discuss the variances of observables that are sums of one-body operators  = N i=1 â(r i ):

1 N ∆ 2 Â = 1 N Â2 -Â 2 (32)
Formally, two-particle operators contribute to the value of this variance, because of the Â2 term in Eq. ( 32),

Â2 = N j=1 â2 (r j ) + N k>j=1 2â(r j )â(r k ). ( 33
)
Using the one-body and two-body reduced density matrices [Eq. ( 25)] to evaluate Eq. ( 32), we obtain

1 N ∆ 2 Â = dr ρ(r) N a(r) 2 -N dr ρ(r) N a(r) 2 + dr 1 dr 2 ρ (2) (r 1 , r 2 , r 1 , r 2 ) N a(r 1 )a(r 2 ). (34)
Evidently, the operator Â2 [Eq. ( 33)] and the variance ∆ 2 Â depend on the coordinates of two particles and are, thereby, two-body operators that can be used to probe many-body physics. Typical choices for Â, which we shall discuss below in Sec. III.C, include the many-body position and momentum operators, X = Faraday waves and "granulation" of a BEC driven with a modulated interparticle interaction strength have been observed in a recent experiment in a quasi-onedimensional setup at Rice University (Nguyen et al., 2019).

Faraday waves result for modulation frequencies on or close to resonance with the transversal trapping (Faraday, 1831) even at rather small-amplitude modulations: Faraday waves are regular, standing, periodic patterns, seen for instance in liquids in a vessel that is shaken. In the experimental realization, the single-shot images of Faraday waves were repeatable [START_REF] Engels | Physical Review Letters[END_REF]Nguyen et al., 2019).

Granulation (Yukalov et al., 2015(Yukalov et al., , 2014) ) results for larger-amplitude modulations with frequencies much lower than the radial confinement: the BEC breaks into "grains" of varying size. The sizes of these grains are broadly distributed, and the grains persist for up to four seconds, i.e., much longer than the modulation time. In the experimental realization, the single-shot images of the granular state -as a direct consequence of quantum correlations -were different, even if all parameters in the experiment were kept fixed (Nguyen et al., 2019).

We stress that the presence of quantum fluctuations and correlations in a many-body state can not be inferred from the density alone. Models like the time-dependent Gross-Pitaevskii mean-field or the time-dependent density functional theory that a priori -by the construction of their ansatz -are aimed at the density may therefore not be able to describe quantum fluctuations and correlations accurately.

A statistical analysis of many observations of the quantum state -i.e., of many (simulated) absorption images in the case of ultracold atoms -is needed in order to study and precisely quantify effects like quantum correlations and fluctuations.

Here, we focus on the case where granulation emerges in the BEC, since the quantum correlations and fluctuations that arise in sync with granulation make this a good example where the application of a wavefunctionbased model like MCTDH-B is crucial, because MCTDH-B (and also MCTDH-F) does incorporate quantum correlations in its ansatz [cf. Eq. ( 2)]. Moreover, the experimental observations in single-shot images can also directly be obtained from the MCTDH-B simulations.

Such a direct comparison of single-shot images simulated from MCTDH-B-computed wavefunctions with the experimental observations on granulation was performed in Ref. (Nguyen et al., 2019). The one-body Hamiltonian used to model the granulation experiment was ĥ

(x) = -1 2 ∂ 2 ∂x 2 + 1
2 Ω 2 x 2 , i.e., a kinetic energy term and a parabolic trap in dimensionless units -the total Hamiltonian was divided by 2 /(mL 2 ), where m is the mass of 7 Li and a length scale L such that Ω ≈ 0.1, see (Nguyen et al., 2019) for details. The time-dependent interaction potential was modelled as

W (x, x ; t) = λ(t)δ(x -x ), ( 35
)
where

λ(t) = λ 0 -β 1 + β1 β2-β3 sin ωt
is the timedependent interaction strength. Here, The MCTDH-B-simulated and the experimental single-shot images do qualitatively agree, see Figs. 4(a),(b).

β 2 = |( B - B ∞ )/∆|, β 1 = -β 2 /(β 2 -1) = 3.
In our present example of the granulation of a BEC, a contrast parameter D that measures discrepancies by more than 20% of experimental and simulated single-shot images from a Thomas-Fermi profile was defined to quantify the amount of fluctuations in the many-body system, see Fig. 4(c).

Since there is no evidence for thermal effects in the experimental realization of granulation, the observed fluctuations are necessarily attributed to quantum correlations. From the contrast parameter [Fig. 4 (c)], we understand that granulation emerges beyond modulation frequencies of ω c ≈ (2π)30Hz and appears side-by-side with quantum correlations, as seen from a significant growth of multiple eigenvalues of the one-body and two-body density matrices [Fig. 4 (d)].

The agreement between the contrast parameter obtained from experimental and simulated single-shot simulations [Fig. 4 (c)] heralds the reliability of the MCTDH-B-prediction for the many-body wavefunction, and the quantum correlations and fluctuations embedded in it.

C. Many-body physics and variances

The inter-connection between mean-field and manybody descriptions of a BEC has attracted considerable attention (Calogero and Degasperis, 1975;Nozières and Saint James, 1982b). Whereas the Gross-Pitaevskii theory has widely been employed in earlier investigations (Burger et al., 1999;Ruprecht et al., 1995), there is nowadays a growing consensus of the need for models that go beyond mean field, as highlighted in Sec. III.B.

Exact and appealing relations between many-body and mean-field descriptions of ultracold bosons are obtained in the so-called infinite-particle-number limit (IPNL), i.e., in the limit where the product of the interaction strength and the number of particles N is kept fixed while the number of particles tends to infinity (Castin and Dum, 1998;Cederbaum, 2017;Erdős et al., 2007a,b;Lieb and Seiringer, 2002;Lieb et al., 2000). In this IPNL, the energy and density per particle, E N and ρ(r) N , respectively, of the BEC computed at the many-body and mean-field levels of theory for N → ∞ are equal; the BEC is 100% condensed.

The Gross-Pitaevskii mean-field theory is obtained as the limiting case when only a single orbital is used with MCTDH-B and computations for a large number of bosons can be done with (RAS-)MCTDH-B, in particular, when the considered state is almost 100% condensed. MCTDH-B is thus very well-suited to investigate the inter-connection between mean-field and many-body descriptions in the IPNL; we will focus on some of the pertinent applications of MCTDH-B in the following dis-cussion.

Even in the IPNL, however, correlations are embedded within a BEC and show in the variance of operators. For the position operator, X = N j=1 xj , where xj is the position of the j-th particle, the effect of correlations can be clearly seen in its variance

1 N ∆ 2 X = 1 N X2 -X 2 ,
see Eq. ( 33) and (Klaiman and Alon, 2015;Klaiman et al., 2016). The reason is that an excitation of as little as a fraction of a particle outside the condensed mode, may interact with a macroscopic number of particles in the condensed mode. Formally, two-particle operators contribute to the evaluation of the variance of oneparticle operators, cf. Eq. ( 33). This is an intriguing result, in particular, because the state is 100% condensed at the IPNL, i.e., the reduced one-particle and two-particle density matrices per particle, ρ (1) N and ρ (2) N (N -1) [Eq. ( 25) for p = 1, 2], respectively, do have only a single macroscopic eigenvalue. In practice, one thus finds a difference when the variance is computed at the many-body and mean-field levels, see Fig. 5(a)-(c) for an example with 1 N ∆ 2 X for bosons in a double well. This difference can be seen as an aspect of the finding that the overlap of the many-body and mean-field wavefunctions can become much smaller than 1 (Klaiman and Cederbaum, 2016). The variance of operators can thus be used to investigate the correlations in BECs that are ignored in mean-field models.

In turn, even at the IPNL the many-body wavefunction is extremely complex and very different from the mean-field one. This difference is caused by only a small amount of bosons outside the condensed mode (Cederbaum, 2017). Since the mean-field and many-body wavefunctions are different, the properties derived from them may also be different. This is particularly true starting from two-body properties, such as the many-particle position variance. When the variance is computed from a mean-field wavefunction it directly relates to the onebody density, because the wavefunction is built as a product of one single-particle state. When the variance, in turn, is computed from a many-body wavefunction it directly relates to the one-body and two-body density, i.e., it contains information about correlations in the wavefunction that is not necessarily built as a product of one single-particle state. The relation between the density of a BEC and the correlations within a BEC can therefore be probed via the variance of operators. The variance can be used as a sensitive diagnostic tool, for the excitations of BECs (Beinke et al., 2018;[START_REF] Theisen | [END_REF], for analyzing the impact of the range of interactions (Haldar and Alon, 2018), and for assessing convergence of numerical approaches like MCTDH-B (Alon and Cederbaum, 2018;Cosme et al., 2016), see Fig. 5(d)-(f) for an example convergence test with the position and momentum space variance, 1 N ∆ 2 X (t) and 1 N ∆ 2 p(t), respectively, in quench dynamics of attractively interacting anharmonically trapped bosons.

The many-body features of the variance of operators in a BEC depend on the strength and sign of the interaction, the geometry of the trap, and the observable under investigation, e.g., the position, momentum, or angular momentum, see, respectively, (Klaiman and Alon, 2015;Klaiman et al., 2016;Sakmann and Schmiedmayer, 2018).

For bosonic systems in two-dimensional ring-shaped traps with a tight transversal confinement where the density exhibits a quasi-one-dimensional behavior and vanishes inside the ring, but, interestingly, the variance remains a manifestly two-dimensional quantity that lives in the inner part of the ring (Alon, 2019).

IV. MCTDH-F AND ELECTRONS IN ATOMS AND MOLECULES

Here we discuss selected applications of MCTDH-F, in some cases with the incorporation of a CAS or RAS scheme, to electron dynamics in atoms and molecules. Before discussing some of the applications of MCTDH-F that contain a comparison with experiment in Sec. IV.B, we introduce the used observables in the following Subsection.

FIG. 5 The position space variance, (a), of N = 1000, 10000, 100000 bosons with contact interactions in a double well as a function of the barrier height on the manybody level (four colored lines atop of each other) drastically differs from the mean-field description (black-dashed line), although the energy per particle (b) and depletion (c) suggest that a mean-field description is applicable; see also Fig. 1 of (Klaiman and Alon, 2015). Convergence of the position variance per particle (d), momentum variance per particle (e), and number of depleted particles (f) in the dynamics of N = 10 attractive bosons in an anharmonic trap, V (x) = 0.05x 4 , following a quench of the strength of the contact interactions [cf. Eq. ( 35) with λ(t) = -0.02 for t ≤ 0 and λ(t) = -0.04 for t > 0], panels (d)-(f) adapted from Fig. 6 in (Alon and Cederbaum, 2018).

A. Extraction of observables

Using (RAS-)MCTDH-F, photoionization cross sections have been calculated using the flux method (Jäckle and Meyer, 1996). The procedure, involving exterior complex scaling, has been described in detail (Haxton et al., 2011) and applications were presented for Be and molecular HF (Haxton et al., 2012).

A direct method is based on expressing the observables of interest in terms of the reduced one-body density [Eq. ( 25) for p = 1], for details see (Madsen et al., 2018;Omiste et al., 2017). To obtain an expression for the photoelectron momentum distribution, the starting point can be the density in coordinate space. The photoelectron distribution can then be obtained by a suitable integral transformation. The density in coordinate space at position r is obtained as the expectation value ρ(r; t) = Ψ| Ψ † (r) Ψ(r)|Ψ , cf. Eq. ( 30). In second quantization, using the orbitals [Eq. ( 3)] and matrix elements of the one-body density [Eq. ( 10)], we obtain the density [see also Eqs. ( 25) and ( 27)]:

ρ(r; t) = kq ρ kq Φ * k (r; t)Φ q (r; t). ( 36
)
To obtain the photoelectron distribution, a projection on an exact scattering state, ψ k (r), should be performed. If this projection is restricted to a region of the simulation volume, beyond an ionization radius, where the effect of potential from the remaining ion is small, the projection can be performed to plane waves; if the long-range Coulomb interaction is still important in that region, the projection may be done to Coulomb scattering waves (Madsen et al., 2007;Omiste et al., 2017). The photoelectron momentum P , is then [cf. Eq. ( 36)]:

dP dk = kq ρ kq Φ * k (k; t) Φq (k; t), (37)
where

Φj (k, t) = drψ * k (r)Φ j (r; t), (38)
and the prime on the integral sign denotes that the integral is only to be evaluated in the outer part of the simulation volume. From the momentum distribution, the energy distribution and the angular distribution can be obtained by integration. The cross section can by obtained from the time-dependent calculation once the ionization probability P 1 is known (Foumouo et al., 2006;Madsen et al., 2000). For example, the photoionization cross section can be extracted by (Foumouo et al., 2006)

σ 1 (Mb) = 1.032 × 10 14 ω 2 P 1 /(n p I 0 ), ( 39
)
where ω is the angular frequency of the laser, I 0 is the peak intensity of the laser pulse in W/cm 2 , and n p is the number of cycles and P 1 is the ionization probability.

Another quantity which we use below and which has received significant interest in strong-field and attosecond physics in recent years is time-delay in photoemission. This field was recently reviewed (Pazourek et al., 2015). The time-delay τ can be extracted in a three-step procedure that we now discuss. (i) From the computed wavefunction, one extracts the expectation value of the radial distance in a given direction, r(t) , as a function of time and the linear momentum of the photoelectron, k that can be evaluated in different ways. For instance k , can be evaluated via integrating only in the outer part of the simulation volume (Omiste et al., 2017;Omiste and Madsen, 2018). (ii) Using r(t) and k the effective ionization time,

t Coul = t - r(t) k = τ EWS + ∆t Coul , (40)
can be evaluated. Here, τ EWS is the Eisenbud-Wigner-Smith (EWS) time-delay, i.e., the time-delay without the interaction with the Coulomb tail of the ion and ∆t Coul = Z k 3 1 -ln(2k 2 t) is the distortion caused by the long-ranged Coulomb potential, where Z is the charge of the ion. (iii) Finally, the time-delay time is evaluated using

τ = τ EWS + τ CLC . ( 41
)
Here, τ EWS can be evaluated from Eq. ( 40) and the

Coulomb-laser-coupling τ CLC = Z k 3 2 -ln( πk 2 ω IR ) is
known, because Z, k and the frequency of the infrared pulse, ω IR are known. Thus the time-delay τ can readily be extracted from the solution of the (RAS-)MCTDH-F EOM (Omiste and Madsen, 2018).

B. Examples involving comparison with experimental results

The processes we will focus on here are in the research area of laser-matter interactions. They are characterized by linear or perturbative interactions, where relatively few photons are exchanged with the external electromagnetic field. This reflects the current challenges with making the MCTDH-F computationally efficient in full dimension and for nonperturbative dynamics where many photons are exchanged. For validation of the MCTDH-F methodology, comparisons with experiments have focused on calculating photoabsorption cross sections (Haxton et al., 2012;Omiste andMadsen, 2018, 2019), where accurate experimental data are available. In addition, XUV transient absorption spectra (Liao et al., 2017) and time delays in photoionization dynamics (Omiste et al., 2017;Omiste and Madsen, 2018) have been considered. Here, we consider cross section and time-delay studies as illustrative examples.

Photoionization cross sections

In the case of photoionization, Fig. 6(a) shows a comparison for atomic neon between the predictions of theory at different levels of approximation and experimental cross section data.

The values of the theoretical cross sections in Fig. 6(a) are obtained by the procedure described in Sec. IV.A. From the agreement between theory and experiment in Fig. 6(a), it can be concluded that it is possible to obtain a precise prediction of the photoionization cross section using an explicitly time-dependent method, the RAS-MCTDH-F, using the procedure discussed in relation to Eq. (39). A second key point to be noticed from the figure is related to the choice of the P 1 and P 2 subspaces and the number orbitals in them. We consider here the RAS-MCTDH-F-D method, cf. Fig. 2 for an illustration of the P 1 and P 2 spaces. The 'D' in the acronym of the method denotes "doubles": only double excitations from the P 1 to the P 2 spaces are allowed. In this example there is no space P 0 with always occupied orbitals like the one used to construct "complete active space" methods (Sato and Ishikawa, 2015). Such a choice of active space and excitation scheme reduces the number of configurations compared with the MCTDH-F method with no restrictions, and as is seen from the figure, can still yield accurate results. The results in the figure show how convergence is obtained by increasing the number of or-FIG. 6 (a) Theoretical total photoionization cross section extracted from a calculation with a 10-cycle linearly polarized pulse with peak intensity 10 14 W/cm 2 as a function of the central angular frequency ω of the laser for several RAS schemes compared to the experimental data by Marr et al. (Marr and West, 1976) and Samson et al. (Samson and Stolte, 2002). (b) Relative time-delay of ionization in Ne, τ2p-2s, as a function of the central frequency of the XUV pulse for a 780 nm IR pulse for (M1, M2) = (5, 0) and (5, 4) together with calculations [START_REF] Dahlström | [END_REF]Feist et al., 2014;Moore et al., 2011) and the measurement (Schultze et al., 2010). Total photoionization cross section (c) and its partitions in the 1π -1 (e), 3σ -1 (d), and 2σ -1 (f) final states. The MCTDH-F computations here used nine orbitals while the complex Kohn ones used eight. The cross sections were computed via the flux into an exterior complex scaling region (Moiseyev, 1998), see (Haxton et al., 2012) for details. The overall agreement between MCTDH-F and experimental results (Brion and Thomson, 1984) is good: for all four depicted cases the salient features are reproduced for the total, 1π -1 , and 3σ -1 cross sections. Reproduced with permission from (Omiste and Madsen, 2018) and (Haxton et al., 2012) with copyright from APS. bitals in P 2 from M 2 = 0 to M 2 = 9. In this manner the accuracy of different approximations from the meanfield TDHF to approaches including more correlation is systematically explored.

Comparisons between theory and experiment for photoionization cross sections have also been performed for atomic beryllium and the HF molecule (Haxton et al., 2012). In these latter cases, full MCTDH-F is considered. Similar to the RAS-MCTDH-F example above, convergence of the MCTDH-F results for the cross sections were obtained with increasing number of active orbitals. We highlight here the good agreement of the photoionization cross sections obtained for HF molecules with the experimental (Brion and Thomson, 1984) and complex Kohn theoretical (Rescigno and Orel, 1991;Schneider and Rescigno, 1988) results, see Fig. 6(c)-(f).

Time delay in photoionization

RAS-MCTDH-F, has been applied to the time-delay in photoionization in neon (Omiste and Madsen, 2018), where experimental data is available (Isinger et al., 2017;Schultze et al., 2010). It is the advent of new light sources for ultrashort light pulses with durations down to the attosecond timescale that has allowed addressing questions like time delay in photoionization in experiments. In Fig. 6(b), the time delay in photoionization between the 2s and the 2p electrons in neon is shown in units of attoseconds (1 as = 10 -18 s). A collection of theoretical results and a measurement point (Schultze et al., 2010) are presented in Fig. 6(b) as a function of photon energy ω (in atomic units, = 1, and for convenience the values in atomic units have been converted to eV, 1 a.u. = 27.21 eV).

The positive value of the time delay can be interpreted as if it takes longer time for the 2p than for the 2s orbital to ionize. Such an interpretation in terms of orbitals, however, assumes a mean-field picture. Theory and experiment have addressed the question about relative time delay between ionization into the two channels

Ne [(1s 2 2s 2 2p 6 ) 1 S e ] → Ne + [(1s 2 2s 2 2p 5 ) 2 P o ] + e -(s, d) (42) Ne [(1s 2 2s 2 2p 6 ) 1 S e ] → Ne + [(1s 2 2s2p 6 ) 2 S e ] + e -(p), (43)
where the dominant configurations have been used to denote the ground state in the neutral as well as the ground and excited state in the ion. Note that dipole selection rules dictate the possible values of the angular momenta in the final channels. From Fig. 6(b), it is seen that all the theories predict a decreasing time delay as a function of the photon energy in the considered energy range. All theoretical values are also smaller than the experimental result. Recently, measurements with an interferometric technique (Isinger et al., 2017) reported a lower value of the time delay in better agreement with the theory results. In the following, we focus on the RAS-MCTDH-F results with (M 1 , M 2 ) at (5, 0) and at (5, 4), see Fig. 6(b). For neon, the (5, 0) results correspond to the TDHF case, i.e., one active orbital for each pair of electrons. The (5, 4) case includes more correlation and has 5 orbitals in P 1 and 4 orbitals in P 2 . The transitions between P 1 and P 2 occur by double excitation. As seen from Fig. 6(b), part of the overall trend of the time delay can be described at the TDHF-level of theory.

Note that there are other cases of interest, where the ionization step can not be captured by TDHF. For example in beryllium, photoionization of the ground state into the channel Be + [(1s 2 2p) 2 P o ]+e -(s or d) changes two orbitals in the dominant configurations by the action of the one-body photoionization operator. Therefore, that process can not be described by TDHF (Omiste et al., 2017).

V. APPLICATIONS, THEORETICAL, AND NUMERICAL DEVELOPMENT

We now discuss theoretical and numerical developments within and beyond (RAS-)MCTDH-X.

A. MCTDH-X-based development

Numerical methods

Since the introduction of MCTDH-F (Caillat et al., 2005;[START_REF] Kato | [END_REF]Zanghellini et al., 2003) and MCTDH-B (Alon et al., 2007b(Alon et al., , 2008a;;Streltsov et al., 2007b) many numerical methods were developed that extend the applicability of MCTDH-X.

The development revolving around the MCTDH-X methods put forward many numerical techniques and theory extensions. For long-ranged interparticle interactions where the interaction potential is a function of the distance of the particles, W (r i , r j ; t) = W (r i -r j ; t), the so-called interaction matrix evaluation via successive transforms (IMEST) has been developed [START_REF] Sakmann | Many-Body Schrödinger Dynamics of Bose-Einstein Condensates[END_REF]. IMEST rewrites the local interaction potentials as a collocation using fast Fourier transforms. IMEST has been applied for solving the TDSE with MCTDH-X, for (time-dependent) harmonic interparticle interactions (Fasshauer and Lode, 2016;[START_REF] Lode | Tunneling Dynamics in Open Ultracold Bosonic Systems[END_REF]Lode et al., 2012a), dipolar interactions (Chatterjee and Lode, 2018;Chatterjee et al., 2019Chatterjee et al., , 2018)), and general longrange interaction potentials (Fischer et al., 2015;Haldar and Alon, 2018) and screened Coulomb interactions (Fasshauer and Lode, 2016).

The development of an implementation of MCTDH-F using a multiresolution Cartesian grid [START_REF] Sawada | [END_REF] holds promise to provide improved adaptive representations for the dynamics of the wavefunction of electrons in atoms and molecules. Moreover, we note the implementation of the infinite-range exterior complex scaling method (Orimo et al., 2018) and the introduction of a space partitioning concept (Miyagi and Madsen, 2017) in combination with RAS-MCTDH-F. We mention that it has been shown that the inclusion of complex absorbing potentials to describe situations like ionization where particles are leaving the region of interest requires one to use a Master equation of Lindblad form for the time-evolution of the density matrix, see Ref. (Selstø and Kvaal, 2010). To solve this Master equation, ρ-MCTDH-F has been formulated in (Kvaal, 2011).

The efficient evaluation of the Coulomb interaction integrals [Eq. ( 16) with Ŵ being the Coulomb interaction] is instrumental to study real-world dynamics of electrons in atoms and molecules in three spatial dimensions. We mention here a sinc-DVR approach that enables an efficient collocation, i.e., Fast-Fourier-transformbased evaluation of the Coulomb interactions in by exploiting the triple-Toeplitz structure of the kinetic energy operator [START_REF] Jones | [END_REF].

We note the recent successful implementation and application of the adaptive removal and addition of configurations, so-called dynamical pruning (Larsson and Tannor, 2017;Wodraszka and Carrington, 2017) for dynamics computed with MCTDH-B (Köhler et al., 2019).

Theoretical progress

The MCTDH-X methodology has also been used to obtain descriptions of the dynamics generated by Hubbard Hamiltonians. In Ref. (Lode and Bruder, 2016), the operators that create/annihilate particles in the time-independent first-band Wannier basis functions of the Hubbard lattice are expressed as effective, creation/annihilation operators that create particles in a time-dependent superposition of all lattice sites. The resulting EOM are identical to the MCTDH-X EOM [Eqs. (8),( 13)], albeit with a special representation of the kinetic and potential energy. In (Sakmann et al., 2011), generalized time-dependent Wannier functions which are a superposition of many bands are proposed, to increase the accuracy of the representation of the many-boson wavefunction beyond the single-band Hubbard model. In (Alon et al., 2014;Grond et al., 2013a), a linear-response framework for the EOM of MCTDH-X, the so-called LR-MCTDH-X, is put forward that allows to obtain highly accurate information about the excitation spectrum of the considered many-body Hamiltonian, as benchmarked in (Beinke et al., 2018(Beinke et al., , 2017)). Recently, the Fourier transform of the auto-correlation function was used to obtain the spectrum for a bosonic many-body system (Lévêque and Madsen, 2019).

For the dynamics of electrons in molecules, an approach termed "multi-configuration electron-nuclear dynamics method" (MCEND) was developed (Nest, 2009) and applied to lithium hydride (Ulusoy and Nest, 2012). This MCEND method represents the total molecular wavefunction as a direct (tensor) product of an MCTDHtype wavefunction for the nuclei with an MCTDH-Ftype wavefunction of the electrons. Other approaches to deal with coupled electronic and nuclear dynamics have been developed and applied for diatomics (Haxton et al., 2011(Haxton et al., , 2015;;Kato and Yamanouchi, 2009;Lötstedt et al., 2019b).

Recent developments of the so-called extended-MCTDH-F in Ref. (Kato and Yamanouchi, 2009) consider coupled electron-nuclear dynamics and molecular wavefunctions and include extensive investigations on H 2 (Ide et al., 2014;[START_REF] Kato | AIP Conference Proceedings[END_REF] and H + 2 in intense laser fields (Lötstedt et al., 2019a,b) as well as a strategy to efficiently partition the CI space of MCTDH-F (Lötstedt et al., 2016).

The multiple active space model put forward in Ref. (Sato and Ishikawa, 2015) introduces a flexible and possibly adaptive approach to construct representations for the N -body Hilbert space with multiconfigurational methods.

We mention here the development, application, and successful benchmark against MCTDH-F predictions for high-harmonic generation of a method that time-evolves the two-body density matrix [cf. Eq. ( 25) for p = 2] without resorting to a wavefunction at all (Lackner et al., 2015).

The unfavorable scaling of the number of coefficients in the MCTDH-X ansatz with the number of orbitals impedes the application of MCTDH-X to systems with many electrons or many bosons with more than a few orbitals. Truncation strategies for the coefficient vector include the RAS approach from quantum chemistry (Olsen et al., 1988) that results in RAS-MCTDH-F (Miyagi andMadsen, 2013, 2014b) and RAS-MCTDH-B (Lévêque and Madsen, 2017;Léveque and Madsen, 2018) theories including a special consideration of single-particle excitations (Miyagi and Madsen, 2014a). The "complete active space" (CAS) truncation approach to limit the number of coefficients was also investigated, see (Sato and Ishikawa, 2013) and, including a generalization to several active spaces (Sato and Ishikawa, 2015). For an MCTDH-F formulation for completely general configuration spaces where different variational principles become inequivalent, see Ref. (Haxton and McCurdy, 2015).

For a review of time-dependent multiconfigurational theories for electronic and nuclear motion in molecules in intense fields see (Kato et al., 2018), for an overview of RAS-MCTDH-X theory see (Madsen et al., 2018).

B. MCTDH-B applications

The archetypical example for the emergence of fragmentation in systems of interacting bosons is the double well potential (Spekkens and Sipe, 1999). Using MCTDH-B for bosons in double-well traps, the reduced density matrices and Glauber correlation functions (Sakmann et al., 2008), the dynamical emergence [START_REF] Sakmann | Many-Body Schrödinger Dynamics of Bose-Einstein Condensates[END_REF]Sakmann et al., 2009Sakmann et al., , 2010;;Streltsov et al., 2007b) and the universality (Sakmann et al., 2014) of fragmentation have been investigated. It is worthwhile to highlight that the works (Sakmann et al., 2009(Sakmann et al., , 2010) ) report converged solutions of the TDSE and demonstrate that the commonly applied Bose-Hubbard model may fail to describe the many-body states for parameter regimes where it was deemed to yield a good approximation to the many-body state. We note that the excitation spectra of interacting bosons in double wells (Grond et al., 2013a;[START_REF] Theisen | [END_REF], in lattices (Beinke et al., 2017), and under rotation (Beinke et al., 2018) have been investigated with LR-MCTDH-B.

Solitons in BEC are thought to be coherent and condensed; several investigations with MCTDH-B (Cosme et al., 2016;Streltsov et al., 2008Streltsov et al., , 2011)), however, have shown that fragmentation and correlations do emerge in their dynamics.

Vortices in ultracold bosonic atoms are conventionally modelled by mean-field approaches (Gross, 1961;Pitaevskii, 1961). Applications of MCTDH-B to interacting bosonic atoms have, however, demonstrated that correlations and fragmentation may emerge as soon as the many-body state contains significant angular momentum (Beinke et al., 2015;Tsatsos and Lode, 2015;Weiner et al., 2017). This emergence of correlations and fragmentation marks the breakdown of the mean-field description and is anticipated from pronounced manybody effects in the excitation spectra of bosonic systems with angular momentum as obtained from LR-MCTDH-B (Beinke et al., 2018).

BECs in high-finesse optical cavities have been used as a quantum simulator for the Dicke model (Baumann et al., 2010;Brennecke et al., 2007). Using MCTDH-B it was shown that the phase diagram of the cold-atom system in the cavity is richer than the phase diagram of the Dicke model and thus the mapping to the Dicke model may break down (Lin et al., 2018;Lode and Bruder, 2017;Lode et al., 2018).

C. MCTDH-F applications

The MCTDH-F was first applied to strong-field ionization of one-dimensional (1D) model molecules with up to eight electrons (Zanghellini et al., 2003), harmonic quantum dots and a 1D model of helium (Zanghellini et al., 2004), and a 1D jellium model (Nest et al., 2007a). Total ionization spectra in strong laser fields were reported for 1D systems with up to six active electrons and strong correlation effects were reported in the shape of photoelectron peaks and the dependence of ionization on molecule size (Caillat et al., 2005). Later, the effect of the reduction in dimensionality from three to one dimension was discussed (Jordan et al., 2006). In the strongfield regime, multielectron and polarization effects have been considered in connection with application to highorder harmonic generation at fixed internuclear distance in model systems (Jordan and Scrinzi, 2008;Miyagi andMadsen, 2013, 2014b;Sukiasyan et al., 2009Sukiasyan et al., , 2010)), in carbon monoxide (Ohmura et al., 2018) as well as helium, beryllium, and neon [START_REF] Sawada | [END_REF].

In molecules, MCTDH-F was applied to H 2 at fixed internuclear distance (Kato andKono, 2004, 2008). The MCTDH-F results reported for molecules include calculations of vertical excitation energies, transition dipole moments, and oscillator strengths for lithium hydride and methane (Nest et al., 2007b), as well as considerations of the response of lithium hydride to few-cycle intense pump fields followed by a probe pulse (Nest et al., 2008). Work on characterizing multielectron dynamics by considering energies and amplitudes was reported (Ohmura et al., 2014). The inclusion of nuclear motion has also been considered (Anzaki et al., 2017;Haxton et al., 2011;Kato and Yamanouchi, 2009;Nest, 2009).

Concerning few-photon processes, MCTDH-F has been applied to the simulation of the two-photon ionization of helium including a comparison with the time-dependent configuration interaction method (Hochstuhl and Bonitz, 2011). The population transfer between two valence states of the lithium atom with a Raman process via intermediate autoionizing states well above the ionization threshold was investigated (Li et al., 2014). A two-color core-hole stimulated Raman process was studied in nitric oxide (Haxton and McCurdy, 2014) and Raman excitations of atoms through continuum levels were considered for neon (Greenman et al., 2017). Moreover, a procedure was suggested for using transient absorption spectroscopy above the ionization threshold to measure the polarization of the continuum induced by an intense optical pulse (Li et al., 2016). Recently, a comparison of MCTDH-F and experimental results was reported in a study using XUV transient absorption spectroscopy to study autoionizing Rydberg states of oxygen (Liao et al., 2017). RAS-MCTDH-F was applied to study electron correlation and time delay in beryllium (Omiste et al., 2017), neon (Omiste and Madsen, 2018), and effects of performing calculations with or without a filled core space (Omiste and Madsen, 2019).

D. Multilayer generalizations of MCTDH-B and MCTDH-F

Multilayering approaches (Vendrell and Meyer, 2011;Wang and Thoss, 2003) provide a powerful and promising generalization of the standard MCTDH. In the multilayer (ML) strategy, the MCTDH is applied recursively: first, the wavefunction is represented as a sum of products of "single-particle" functions that still may be highdimensional functions (first layer); second, the "singleparticle" functions of the first layer are again represented by an MCTDH-type wavefunction, i.e., a sum of products of (second-layer) "single-particle" functions, and so on. In the bottom or last layer, the single-particle functions are expanded on a primitive time-independent basis.

This multilayering strategy has recently also been applied for MCTDH for indistinguishable particles with two different approaches: one, the so-called ML-MCTDH in (optimized) second quantized representation [ML-MCTDH-(o)SQR], where the commutation relations of indistinguishable bosons or fermions are included in the underlying time-independent or an optimized time-dependent "primitive" basis, see Refs. (Wang and Thoss, 2009) and (Manthe and Weike, 2017), respectively. The ML-MCTDH-SQR approach has, for instance, been successfully applied to the dynamics of the spin-boson model (Wang and Shao, 2019;Wang andThoss, 2009, 2016) and transport in the Anderson impurity model (Wang and Thoss, 2018).

A different approach, the ML-MCTDH-X, is to use an MCTDH-type representation for the combined wavefunction of multiple species or the single-particle functions of indistinguishable bosons or fermions, which are, in turn, represented by an MCTDH-X-type expansion (Cao et al., 2017(Cao et al., , 2013;;Krönke et al., 2013). This approach has been applied successfully to mixtures of ultracold bosons and fermions [START_REF] Erdmann | [END_REF]Mistakidis et al., 2018;Siegl et al., 2018) and bosons in more than one spatial dimensions (Bolsinger et al., 2017a,b).

We mention here a fundamental relation between density matrix renormalization group methods (Schollwöck, 2005(Schollwöck, , 2011) ) and ML-MCTDH-X/ML-MCTDH-SQR: mathematically, both methods fall into the class of socalled hierarchical low-rank tensor approximations, a concept which has enabled progress in devising new efficient time integration schemes (Falcó et al., 2019;Lubich et al., 2018) that are applicable for (RAS-)MCTDH-X.

E. Self-consistent field coupled-cluster

To reduce the numerical effort in solving the TDSE to become polynomial, the so-called coupled cluster method (CC) ( Čížek, 1966, 1969;Čížek and Paldus, 1971;Coester and Kümmel, 1960) can be employed. Although CC uses a different type of ansatz than MCTDH-X, we mention it here, because recent developments include approaches with a time-dependent, variationally optimized basis and are thus related to MCTDH-X and RAS-MCTDH-X.

The conventional CC uses time-dependent excitation amplitudes, but does not use a set of time-dependent orbitals in the representation of the wavefunction. The standard CC's ansatz can be generalized to include timedependent amplitudes and orbitals. This generalization of the ansatz in combination with a generalized, so-called bivariational principle, leads to the equations-of-motion of the orbital-adapted time-dependent coupled cluster theory (Kvaal, 2012(Kvaal, , 2013;;Pedersen and Kvaal, 2019). We identify the application of the bivariational principle for the derivation of the MCTDH-X EOM for ansatzes with restricted configuration spaces [like in Eq. ( 17)] as an open question.

When a real-valued variational principle is used, the fully time-dependent coupled cluster ansatz yields the EOM of the time-dependent optimized CC (Sato et al., 2018a,b). The latter theory allows the self-consistent computation of eigenstates via imaginary time propagation and has been applied to single-and double ionization as well as high-harmonic-generation in argon (Sato et al., 2018a).

VI. CONCLUSIONS AND FRONTIERS

In this Colloquium, we introduced the MCTDH-B and the MCTDH-F methods for full and for restricted configuration spaces. We highlighted the use and versatility of MCTDH-X with benchmarks against exactly solvable models as well as direct comparisons with experimental applications.

The development of methods for the time-dependent many-body Schrödinger equation in the field of MCTDH-X and beyond, that we have portrayed in our present Colloquium, has yielded highly efficient and flexible numerical approaches. This flexibility, however, comes with an increasing number of parameters to tune the performance and accuracy of the given approach -we name here as examples the tree structure in multilayering approaches (Cao et al., 2017;Manthe and Weike, 2017;Wang and Thoss, 2009), and the partitioning of Hilbert space into multiple occupation-restricted active spaces (Sato and Ishikawa, 2015) or into P 1,2 (Fig. 2) in the RAS-MCTDH-X approach (Lévêque and Madsen, 2017;Léveque and Madsen, 2018;Miyagi andMadsen, 2013, 2014b). We thus observe that the recent methodological developments demand an ever larger and more complicated set of parameters -like the mulilayering tree or the partitioning of Hilbert space into active spacesto be configured by their users. Such a development towards higher complexity in the application of methods is not desirable, because it makes applications ever more tedious. The trend towards more complexity could possibly be overcome by introducing additional adaptivity.

We mention here the recent fascinating developments with adaptive tensor representations (Ballani and Grasedyck, 2014;Grasedyck et al., 2013), an adaptive number of configurations (Haxton and McCurdy, 2015;Larsson and Tannor, 2017;Lévêque and Madsen, 2017;Miyagi andMadsen, 2013, 2014b;Wodraszka and Carrington, 2017), adaptive grids [START_REF] Sawada | [END_REF], and an adaptive construction of many-particle Hilbert space (Sato and Ishikawa, 2015). We thus envision a flexible implementation that combines multiple of the above multiconfigurational methods in an adaptive framework to solve the many-particle Schrödinger equation: according to a simple/single input -for instance an error threshold -Hilbert space is automatically and adaptively partitioned and represented while for each portion of it (an adaptive version of) the best-suited of the multiconfiguration methods is used.

Interestingly, the extended-MCTDH-F and multiconfiguration electron-nuclear dynamics method (MCEND) ansatzes, proposed in (Kato and Yamanouchi, 2009) and (Nest, 2009), respectively, represent the total wavefunction as a (tensor) product of wavefunctions of different species of particles. In the case of extended-MCTDH-F, the wavefunction is a product of two MCTDH-F-type wavefunctions and in the case of MCEND, the wavefunction is a product of an MCTDH-F-type wavefunction with an MCTDH-type wavefunction for distinguishable particles.

Such a multi-species wavefunction -as well as bulk of the multiconfigurational methods developed for restricted, multiple, and general active spaces -is amenable to multilayering approaches. The combination of truncation methods for the configuration space, including the so-called dynamical pruning approaches (Köhler et al., 2019;Larsson and Tannor, 2017;Wodraszka and Carrington, 2017), with ML-MCTDH-X/ML-MCTDH-(o)SQR is one of the frontiers that we see in the further development with MCTDH-X approaches.

FIG. 1

 1 FIG. 1 Illustration of the configuration space of MCTDH-B [(a)] and of MCTDH-F [(b)]. The space spanned by the timedependent single-particle basis for which all configurations are considered is denoted by P and its complement is denoted by Q. For bosonic particles, (a), the occupation numbers nj are unrestricted, cf. the given five-orbital configuration vectors |n1, ..., n5 . For spin-1 2 -fermions, (b), the Pauli exclusion limits the occupations to be at most two electrons per spinorbital, nj ≤ 2, see the given configurations ( 1 [1] indicates a spin-down [-up] fermion).

FIG. 3

 3 FIG.3Benchmark of RAS-MCTDH-B against exact TD-HIM results for N = 10 bosons. Here, we use f (t) = sin(t) cos(2t) sin(0.5t) sin(0.4t) (upper panel) and K = 0.5 in Eqs. (23) and (24). The time-dependent center-of-mass energy (t) is plotted in comparison to the exact values in the lower panel for different particle and orbital numbers. A convergence with an increasing number of orbitals, i.e., amount of variational parameters in the (RAS-)MCTDH-B wavefunction, is observed. See(Lode et al., 2012a) for details on (t).

  fluctuations and correlations in systems of ultracold bosons

  10 and β 3 = |∆B/∆| are the parameters of the applied time-dependent magnetic field B(t) = B + ∆B sin ωt, where B ∞ = 736.8G, B = 590.9G, and ∆ = 192.3G. Importantly, the sinusoidal modulation of the magnetic field creates a periodic but non-sinusoidal modulation of the interparticle interaction strength λ(t).

FIG. 4

 4 FIG. 4 Experimental and theoretical single-shot line density profiles. (a) Experimental data and (b) many-body simulations for different modulation frequencies. (a) The rows show data for three independent experimental images (single shots) for the indicated ω, where ω = 0 corresponds to no modulation. The interaction between particles was modulated for tm = 250ms around an average value of 8a0 with a maximum of 20a0 and a minimum of 0.7a0; subsequently, the interactions are held constant for another t h = tm = 250ms. (b) The first column shows the density ρ(z, t) [Eq. (25) for p = 1, r1 = r 1 = z] as calculated from the one-dimensional MCTDH-B computations, while the second and third columns display two simulated single shot images [Eq. (28)]. We observe that granulation is present in single-shot images, but absent in the average, ρ(z, t). Quantum fluctuations characterize the emergence of granulation: (c) Comparison of the deviations from a Thomas-Fermi distribution as quantified by the contrast parameter D = D(ω) [see (Nguyen et al., 2019) for details about D] for single shots simulated from wavefunctions computed with MCTDH-B and single shots taken in experiment (EXP). MCTDH-B predicts the threshold value, ωc ≈ (2π)30Hz, where deviations become large and grains form. Each symbol and its error bar are the mean and standard error of the mean of at least 4 experimental measurements of D, while 100 single shots at each ω have been simulated from the MCTDH-B wavefunctions. (d) Eigenvalues of the first-and second-order reduced density matrices. A growth of n (1) 2 , n (2) 2 , and n (2) 3 are observed to occur for ω > ωc, indicating the emergence of correlations and fragmentation. The growth of both n (1) 2

ACKNOWLEDGMENTS

Financial support by the Austrian Science Foundation (FWF) under grant No. P-32033 and M-2653, the Israel Science Foundation (Grant No. 600/15), and by the VKR Center of Excellence, QUSCOPE, and computation time on the HazelHen Cray computer at the HLRS Stuttgart is gratefully acknowledged.