Time-dependent restricted-active-space self-consistent-field theory for bosonic many-body systems

We describe the time-dependent restricted-active-space self-consistent-field (TD-RASSCF) method for a system of interacting bosons. We provide the theory of the method and discuss its numerical implementation. The method provides a general wavefunction based approach to solve the time-dependent and time-independent Schr\"odinger equation for a system of bosons. It is based on the time-dependent variational principle to optimize at each instant of time a set of time-dependent coefficients and time-dependent orbitals used to describe the total wavefunction. Including the concept of a restricted-active-space, the exponential growth of the configurational space, resulting from all possible distributions of $N$ bosons in $M$ orbitals, can be controlled trough a specific excitation scheme. We show, by illustrative time-independent and time-dependent examples, that the method provides an accurate description of the system with a substantially smaller configurational space than the one required in the multi-configurational time-dependent Hartree method for bosons (MCTDHB). The TD-RASSCF method can also tackle problems beyond the reach of the MCTDHB method when a large number of orbitals are required.


I. INTRODUCTION
Since the first realizations of Bose-Einstein condensates (BEC) [1][2][3], the experimental and theoretical investigation of trapped cold atoms has attracted much attention. It is nowadays experimentally possible to design and control systems with a specific number of atoms [4,5] trapped in various potential shapes [6,7] and dimensions [8], with tunable inter-particle interactions [9,10], and to provide a controllable transition from a few-to a manyparticle system. Such a detailed control of cold atom systems has opened the possibility to simulate various physical systems [11] from solid-state physics [12] to blackholes analogs [13] through matter-light interaction [14] and electrons dynamics in molecules [15].
Various theoretical models [16] have been used so far to describe static and dynamical properties of many-boson systems, among which, a handful are exactly solvable. One of the most prominent models was introduced by Lieb and Liniger [17,18], to describe a system of spinless bosons interacting through a two-body contact interaction: using the Bethe Ansatz and periodic boundary conditions the resulting Schrödinger equation can be solved exactly for any interaction strength and an arbitrary number of bosons. Unfortunately, this model is exactly solvable only without a trapping potential. In the limit of infinite interaction strength, the Tonks-Girardeau model use the Fermi-Bose mapping to map the wavefunction of bosons into a fermonic wavefunction of noninteracting fermions with frozen parallel spins [19]. This mapping provides the exact solution for the ground-state of the system for arbitrary trapping potentials and re- * camille.leveque@phys.au.dk mains valid also for the excited states, as well as nonequilibrium solutions also for any external potential [20]. In the case of non-interacting bosons or more generally in the Gross-Pitaevskii (GP) limit, i.e., N → ∞ and N λ = constant, with N the number of bosons and λ the interaction strength, the GP equation or its timedependent (TD-GP) analog provides the exact description of the system. In this situation, the exact wavefunction of the system is described by a single product of single-particle functions and the interactions between the particles are correctly described by the mean-field approach. The above models assume that the bosons interact through a pair-wise contact potential. Considering other types of interaction potentials between the particles, other models can be solved exactly with an external potential. One model uses an inverse-harmonic interaction between the particles and can be solved exactly with a harmonic trapping potential [21], while an other model considers a harmonic interaction potential [22,23]. The latter model has the peculiarity that it can be solved exactly numerically also for time-dependent Hamiltonians with a time-dependent trapping potential or a timedependent interaction potential [24].
These exact models, unfortunately, do not cover the large variety of interaction or trapping potentials that are encountered in experiments. Nonetheless, they are of primary interest as they provide a unique way to benchmark numerical methods and approximations. The GP equation can be simplified when the potential and interaction energies are much larger than the kinetic energy, giving rise to the Thomas-Fermi approximation when the kinetic energy is neglected [25]. On an other hand, to overcome the lack of correlation in the GP theory and to take into account a small depletion of the BEC, i.e., to account for atoms which are not in the condensate, a perturbative expansion of the particle number in the condensate leads to arXiv:1612.04419v1 [quant-ph] 13 Dec 2016 Bogoliubov theory [26][27][28]. In the specific case of periodic trapping potentials, such as optical lattices [8], for weak contact interactions and deep lattices the Bose-Hubbard model (BHM) [29] is obtained by expanding the Bose field operator in term of the Wannier functions of the lowest Bloch band and neglecting the tunneling between nonconsecutive sites and interactions between different sites. The BHM and its various extensions have been extensively and successfully used to describe the ground state of trapped atoms in optical lattices and their dynamics [30]. A more general and efficient numerical approach to deal with optical lattices is the density-matrix renormalization group (DMRG) method [31][32][33] based on the matrix product states Ansatz [34]. The method has been used to provide accurate results for ground and exited states of the system, and more recently has been used to investigate time-dependent systems [35][36][37]. The second wide-spread and promising numerical method to study trapped atoms is the quantum Monte Carlo (QMC) approach. It includes, among others, the variational Monte Carlo (VMC) [38] and the diffusion Monte Carlo (DMC) [39][40][41][42] methods, which used a Bijl-Jastrow decomposition of the wavefunction [43,44], but are, however, not applicable to time-dependent systems.
Along with the above theory developments it has been a long standing idea to explore quantum chemistry methodologies to describe a time-independent system of trapped cold atoms. This idea was, to the best of our knowledge, introduced by the work of Ersy [45], applying the mean-field Hartree-Fock (HF) theory and the configuration interaction (CI) method up to double excitations (CISD) to harmonically trapped bosons. The HF method for bosons can be viewed has a variant of the GP theory but has the advantage that it provides a set of optimized virtual orbitals, i.e., non-occupied orbitals, that can be subsequently used in a CI expansion of the wavefunction. The CI expansion corrects the lack of correlation between the particles, not included at the HF level. The CI method is in principle exact but requires a severe truncation of the CI expansion to be numerically tractable. Later, Streltsov et al [46] introduced the multiconfigurational Hartree theory for bosons (MCHB), which is an extension of the multiconfiguration self-consistent field (MCSCF) method introduced for fermions and widely used in electronic-structure calculations in atoms and molecules [47]. The MCHB method uses a CI expansion Ansatz for the many-body wavefunction in which both the coefficients of the expansion and the orbitals are variationally optimized, providing better accuracy with substantially less configurations and orbitals. The coupled-cluster (CC) method was originally introduced in nuclear physics [48,49] and subsequently extended to describe electronic wavefunctions in atoms and molecules [50]. This framework was also extended to bosons up to double excitations (CCSD) by Cederbaum et al, and successfully applied to various particle numbers and interaction strengths [51].
Over the past decade, numerous numerical methods have been developed [52][53][54][55][56][57] to tackle the problem of time-dependent multi-electron dynamics induced by laser pulses that are strong or short or both [58][59][60]. In short, the various successful methods used so far to investigate static properties of atoms and molecules have been extended to solve the time-dependent Schrödinger equation including a time-dependent operator. Among these methods, the multiconfigurational time-dependent Hartree-Fock method [61][62][63][64] variationally optimizes a set of time-dependent orbitals and CI coefficients, following the idea of the multiconfigurational time-dependent Hartree (MCTDH) method [65,66], originally introduced to describe molecular dynamics. The MCTDHF method has been extended to identical bosons, within the framework of the MCTDH for bosons MCTDHB [67], in which the indistinguishability is taken into account using permanents instead of Slater determinants. Further development includes the case of particle mixtures of different type of bosons and fermions [68,69]. The fundamental concept of using a set of time-dependent single-particle functions or orbitals to expand the total wavefunction offers the possibility to use substantially less orbitals than in the case of time-independent orbitals, because the former basis optimally adapts during the evolution of the system. Recently, the framework of the multi-layer (ML) MCTDH method [70][71][72] was extended to systems of bosons and mixtures of them [73,74]. The method uses a ML expansion to reduce the size of the wavefunction in comparison to the MCTDHB method for multi-species or multi-dimensional systems. It is particularly effective for systems which can be subdivided in strongly interacting subsystems while the individual subsystems interact only weakly with each other. In the case of a one-dimensional system consisting of only a single type of particles, the ML-MCTDHB and MCTDHB wavefunctions are identical [75]. The MCTDHB method shed new light on the dynamics of trapped cold atoms, especially when fragmentation occurs and more than one orbital is populated -a situation which can not be describe by the TD-GP theory. Fragmentation occurs in different systems such as during the dynamics at a Josephson junction [76], which is a universal phenomenon [77], and can not be described, even qualitatively, using the TD-GP or BH theories. In double-well trapping potentials, fragmentation of the BEC is also obtained for the ground-state [78] for large barrier height between the two wells and the GP theory fails to describe the variance of position and momentum operators [79]. Multiconfigurational methods are also required to accurately describes the formation and dynamics of fragmented states with repulsive or attractive interactions between the particles [80][81][82] and tunneling of a many boson system to open space [83] or tunneling of trapped vortices [84]. Using time-dependent orbitals reduce the number of orbitals and thus the number of configurations required to describe accurately time-evolving systems in comparison with methods with time-independent orbitals. Nevertheless, simulations using such full-configurational wavefunc-tions remain a difficult task due to the exponential scaling of the configurational space, i.e., the dimensionality determined by the number of ways to arrange N particles in M orbitals, especially for bosons.
This challenge leads us to the quest for a method which maintains the appealing properties of the time-dependent orbitals based methods mentioned above, but is free from the exponential scaling problem. One such method uses the concept of a restricted active-space (RAS), well-known in quantum chemistry, where it has been applied with time-independent molecular orbitals [85]. The RAS based method was successfully extended to time-dependent orbitals in the time-dependent restricted active-space self-consistent-field method (TD-RASSCF) to deal with electron dynamics in atoms [86,87]. Introducing a RAS scheme by fixing the promotion of the electrons between three sets of orbital spaces can considerably reduced the number of configurations. In addition, the theory has the specificity to include, as limiting cases, the TD Hartree-Fock (HF), the TD complete active space self-consistent field (TD-CASSCF) [56] and the MCTDHF frameworks, as a particular RAS schemes are applied to the MCTDHF wavefunction. Successful applications of the TD-RASSCF method include calculations of the ground-states (GS) of atoms, and timedependent dynamics in the presence of strong laser fields to describe, for instance, high-order harmonic generation [86,87]. The aim of this work is to extend the TD-RASSCF method to systems of spinless interacting bosons. To follow the generic naming introduced for the MCTDH methods, we call this method TD-RASSCF-B where the additional B stands for bosons and we will refer to the original TD-RASCSF method for fermions [86][87][88] as TD-RASSCF-F to avoid any confusion concerning the particles considered. As a main finding, we derive the working equations of the TD-RASSCF-B method and we present the general set of working equations for the TD-RASSCF method where the type of particles plays a role in the symmetry of the creation and annihilation operators, only. The applications of the method to compute the GS energy of trapped bosons show that the TD-RASSCF-B theory provides accurate results, in comparison to MCTDHB, while the expansion of the wavefunction is considerably reduced. Moreover, the MCT-DHB accuracy can be overtaken by using large numbers of time-dependent orbitals while the number of configurations remains small thanks to the RAS schemes. The investigation of the breathing dynamic of a BEC illustrates how the TD-GP theory fails to describe the timeevolution of the system, while various examples of the TD-RASSCF-B method qualitatively or quantitatively reproduce the exact dynamics obtained using the MCT-DHB method, depending of the choice of the excitation scheme.
The paper is organized as follows. In Sec. II A we introduce the TD-RASSCF-B Ansatz for the wavefunction and in Sec. II B we derive the equations of motion for the set of coefficients and orbitals. In Sec. III the method is applied and compared to the MCTDHB method to study the static properties of a system consisting of N = 100 bosons trapped in a harmonic potential. The applicability of the method to time-dependent systems is illustrated by two examples of a breathing dynamics following a sudden quenching of the two-body interaction in Sec. IV. Finally, in Sec. V we conclude and provide perspectives to future work. In the Appendices A to D, we provide the key ingredients for the numerical implementation of the method and discuss the numerical effort in comparison to the MCTDHB method.

II. THEORETICAL FRAMEWORK
A. Ansatz for the many-body wavefunction For the energy regime of interest, the time evolution of a system composed of N bosons is governed by the time-dependent Schrödinger equation: withĤ(t) the many-body Hamiltonian of the system and |Ψ(t) the N -particle wavefunction. Hereafter we set h = 1, unless explicitly specified. We can approximate the wavefunction using linear combinations of suitably symmetrized sets of products of time-dependent singleparticle functions {φ i (r, t)}. In the following, the singleparticle functions are denoted orbitals. To take into account the indistinguishability of the bosons, the total wavefunction is expressed in terms of permanents. For a given number of bosons and orbitals the multiconfigurational wavefunction is constructed by taking into account all the possible arrangement of the particles in the given orbitals, each arrangements being called a configuration |Φ I (t) , This Ansatz converges to the exact wavefunction when the number of orbitals increases to infinity. The configurational space increases exponentially with respect to the number of orbitals and often makes a numerical treatment impossible, even for a small number of orbitals. In the case of a system of N bosons and M orbitals, the dimension of the full-configurational Fock space V FCI can be evaluated as, In the case of the TD-RASSCF-B method, we introduce two orbital spaces, P 1 and P 2 , such that M 1 + M 2 = M , with M 1 and M 2 the number of orbitals in P 1 and P 2 , respectively (see Fig. 1). The P 1 subspace must include enough orbitals such that it can accommodate all the particles. For bosons one orbital, i.e., M 1 = 1 is the lower bound, and there is no restriction concerning the upper bound. In this subspace all the configurations are used to construct the total wavefunction. Concerning the P 2space, particles are promoted from P 1 to P 2 according to a specific excitation scheme, which is based on the highest number of particles that can be promoted. This number is chosen at will and the restriction of the configurational space provides a way to constrain its size, defining the Ansatz for the TD-RASSCF method as, where the configurations are drawn from the space V subject to restrictions. To evaluate the size of the V, we introduce N max the highest number of bosons in P 2 and consider a RAS scheme allowing all occupations of P 2 from 0 to N max . The dimension of the V of Eq. (4), is then given by The first term is the total number of configurations obtained with the N bosons in the M 1 orbitals and no particle in P 2 . The sum takes into account the configurations resulting from the excitation of k bosons in P 2 , with 1 ≤ k ≤ N max . The total number of configurations with k bosons in P 2 is obtained as a product of the possible arrangements of k bosons in M 2 orbitals and (N − k) bosons in M 1 orbitals, see also Appendix A.
The TD-RASSCF-B Ansatz holds some interesting specificities. First, if only P 1 orbitals are used, i.e., M 1 = M and M 2 = 0, then the TD-RASSCF-B and MCTDHB Ansätze are equivalent with the same number of configurations, as seen be replacing M 1 by M in Eq. (5). Note that this is also true for M 2 = 0 and N max = N . Moreover, if only a single time-dependent orbital is considered, i.e., M 1 = 1 and M 2 = 0, the RAS wavefunction includes a single configuration with all particles in one orbital, which is equivalent the to time-dependent GP wavefunction. Thus the theoretical framework of TD-RASSCF-B is very general and holds, as limiting cases, the GP and MCTDHB theories. The TD-RASSCF-B wavefunction is built from a set of timedependent coefficients {C I (t)} and orbitals {|φ i (t) }. To describe its dynamics, we need a set of equations of motion (EOM), which provides the time-evolution of the coefficients and orbitals through their time-derivatives {Ċ I (t)} and {|φ i (t) }. The P-space is a subset of the total single-particle Hilbert space and we can define its orthogonal complement, Q, collecting the virtual orbitals, as depicted in Fig. 1. While in the case of time-independent orbitals these two subspaces remain fixed, in the case of time-dependent orbitals the P-space is variationally optimized at each time and both P-and Q-space are timedependent. We can defineP andQ, the time-dependent projectors onto the subspaces P and Q, respectively, with the propertyP +Q =1, the identity operator. The role of the Q-space emanates from the time-derivative of the P-space orbitals, that can be written as, with one contribution from the P-space and one contribution from the Q-space. In the following we establish the EOM of the TD-RASSCF-B theory, providing the timederivative of the expansion coefficients in Sec. II B 1 and the time-derivative of the orbitals (Sec. II B 2) through the Q-and P-space contributions in Secs. II B 2 a and II B 2 b, respectively.

B. Derivation of the working equations
The EOM for the TD-RASSCF-F theory have been already established in Refs. [86,87]. In the following we provide the derivation of the EOM in the case of the TD-RASSCF-B method, and highlight the differences with respect to the TD-RASSCF-F theory. Starting from the Lagrangian formulation of the time-dependent Schrödinger equation [89], we define the action functional using the TD-RASSCF-B Ansatz, Eq. (4), as, withK ≡ i∂/∂t −Ĥ and δ ij the Kronecker delta function. The Lagrange multipliers, i j (t), ensure that the orbitals remain orthonormal for all time t. In the following the indexes i, j, k, · · · are used to denote the orbitals of the P-space, the indexes a, b, c, · · · denote the orbitals of the Q-space and p, q, r, · · · are used for either P-or Q-space orbitals, see also Fig. 1. For our purpose, we consider only one-and two-body operators, such that the Hamiltonian can be expressed in the framework of second quantization aŝ withb p (b † p ) the annihilation (creation) operator of a particle in the orbital |φ p (t) [see also Appendix B]. These operators satisfy the commutation relation, for bosons and the anti-commutation relation, {b p ,b † q } =b pb † q +b † qbp = δ qp , for fermions, see for instance Ref. [90]. The matrix elements of the onebody and two-body operators in the basis of the timedependent orbitals, are expressed as and v pr respectively. In the following, the explicit time dependence of the operators, coefficients and orbitals is dropped for brevity. According to the time-dependent variational principle [89,91], the best approximation using the wavefunction Ansatz is obtained by seeking stationarity of the action S, i.e., δS = 0, for any variation of the parameters and with the boundary condition |δΨ(t 1 ) = |δΨ(t 2 ) = 0. The variation of the action gives, where the boundary condition is used to remove the additional term i∂ t Ψ|δΨ resulting from the action ofK on Ψ| instead of |δΨ , see Ref. [91]. The variation of the wavefuntion is explicitly written as [86], We can now proceed with the stationarity condition of the action with respect to the parameters {C I }, {|φ i } and { i j } to obtain the EOM of the TD-RASSCF-B method.

Equations of motion for the coefficients
The variation w.r.t. the Lagrange multipliers leads to the conservation of the orthonormality of the orbitals. We then consider the variation of the action functional with respect to the expansion coefficients. The action S depends on the expansion coefficient C * I only through the bra δΨ| of the first expectation value in Eq. (11). Thus, the stationarity condition, δS/δC * I = 0, readily leads to Φ I |i∂ t −Ĥ|Ψ = 0. Moreover, the derivative of the wavefunction with respect to time reads, with η p q ≡ φ p |φ q , which results from the time-derivative of the orbitals used to build the configurations in |Ψ . Hereafter, the operator in bracket in Eq. (13), pq η p qb † pbq , will be calledD, for brevity. We can now rewrite the stationary condition δS/δC * I = 0 using the explicit form of ∂ t |Ψ , Eq. (13), as or equivalently, using the expressions ofĤ andD, The indexes in the summations are now restricted to the P-space. It is clear that if either annihilation or creation operators act on an orbital of Q, the inner product with all RAS configurations Φ I | vanishes. The EOM for the expansion coefficients, the amplitude equations (15), are identical to those obtained for fermions in the TD-RASSCF-F theory, see Refs. [86,87], and those of the MCTDHB [67] and MCTDHF [92] theories. It is worthwhile to keep in mind that the action of the creation and annihilation operators differs for fermions and bosons. The η i j matrix elements in the amplitude equations [Eq. (15)] describe the rotation of the orbitals into one another and are also present in the EOM of the MCT-DHB/F methods. In these latter cases, besides to be elements of an anti-Hermitian matrix, there are no constraints on the η i j and their values are usually set to zero. The same is true in the TD-RASSCF-B/F methods for equivalent orbitals, i.e., for pairs of orbitals which belong to the same P i -space (i=1, 2). For orbitals which do not belong to the same P i -space, the η i j matrix elements must be evaluated, as discussed in Refs. [86,87] and in Sec. II B 2 c and II B 2 d.

Equations of motion for the orbitals
Seeking stationarity of S with respect to a variation of an orbital φ i |, i.e., δS/δ φ i | = 0, gives with Ψ q i | ≡ Ψ|b † ib q . The index q in the above equation runs over all the orbitals, i.e., the orbitals of the P-space and the Q-space, see Fig. 1. The EOM for the orbitals of the P-and Q-space are obtained by projecting Eq. (16) on either an orbital of the P-space, φ j |, or of the Q-space, φ a |, as done in the following.

a. Equations of motion for the Q-space orbitals
Starting with the EOM for the Q-space orbitals, we multiply Eq. (16) from the left with an orbital φ a | belonging to the Q-space and obtain, where we used the orthogonality between the orbitals of the P and Q spaces to get rid of the Lagrange multipliers. Moreover the inner product Ψ a i |Φ I = Ψ|b † ib a |Φ I vanishes because in all configurations |Φ I the orbital |φ a is unoccupied. Using the explicit expression of the Hamiltonian, Eq. (8), and for the operatorD, Eq. (13), we obtain Using the commutation relation for the creation/annihilation operators for bosons (fermions), we can reestablish the normal ordering of the chains of operators, with the upper sign holding for bosons and the lower for fermions. The chain of four operators in Eq. (19) and six operators in Eq. (20) both annihilate a particle in orbital |φ a , from the Q-space, which is not include in |Ψ and thus vanish. The l.h.s. of Eq. (18) now reads, where we restrict the summation over j ∈ P, the summation over the Q-space orbitals being zero. In the same way, inserting Eq. (20) in the r.h.s. of Eq. (18) simplifies its expression to (10)]. Interestingly, this equation is exactly the same for bosons and fermions. The time-derivative of the orbitals, included in the term η a j , requires the explicit consideration of the Q-space orbitals. This issue is circumvented by using the projector Q onto the subspace spanned by the Q-space orbitals, withP the projector onto the P-space. Introducing the one-body density matrix, ρ j i = Ψ|b † ib j |Ψ and the twobody density matrix ρ jl ik = Ψ|b † ib † kb lbj |Ψ , we obtain, the mean-field operator, which describes the interaction between the particles. The role of the Q-space appears in the time-derivative of the orbitals of the P-space through the termQ|φ i , see Eq. (6). We rearrange Eq. (24), see appendix C, to uncouple the contribution of eacĥ Q|φ i , ∀i ∈ P, and obtain with ρ −1 the inverse of the one-body density matrix. The MCTHB theory leads also to Eq. (26), see Ref. [67], but the l.h.s. is subsequently simplified thanks to the choice of the matrix elements η i j = 0 and usingQ =1 −P , see Eq. (30) below. As discussed in II B 1, such a fixed choice of η i j is not possible in the TD-RASSCF theory. The derivation of Q-space EOM differ slightly for bosons and fermions, see Eqs. (19) and (20), but the final result, Eq. (26), is the same for both types of particles.

b. Equations of motion for the P-space orbitals
Going back to the stationary condition for the variation of the action functional with respect to an orbital, Eq. (16), we multiply this latter on the left by an orbital of the P-space, φ j |, leading to, This equation still contains the Lagrange multiplier i j . A variation of S with respect to the orbital |φ j and its projection onto the orbital φ i |, leads to an equation containing the same Lagrange multiplier, and subtracting Eq. (27) and Eq. (28) gives the EOM for the P-space orbitals, i.e., The P-space EOM provide the contribution of the P-space in the time-derivative of the orbitals, see Eq. (6),P through the evaluation of the matrix elements η j i included in the operatorD. Nonetheless, solving Eq. (29) is not a trivial task because of the presence ofρ j i , which couples the amplitude and P-space orbitals equations. In the case of the wavefunction based on the RAS Ansatz, a freedom in the choice of the elements η j i is still possible for pairs of orbitals which belong to the same P i -space (i = 1, 2), and we use η j i = 0, ∀{i, j} ∈ P 1 or P 2 . The P-space equation [Eq. (29)] remains to be solved only for pairs of orbitals {i , j }, which belong to different P i -spaces,  (6) can be evaluated. In the derivation of the TD-RASSCF-F method, a way to circumvent the difficulty of solving Eq. (29) was proposed [86,87]. This approach will be used in the following also for bosons.

c. Even excitation RAS scheme
First we suggest to consider the case in which only an even number of particles is promoted from P 1 to P 2 , see Fig. 2a. In this case,ρ j i explicitly reads, The action ofb † i b j on the wavefunction |Ψ annihilates one particle in P 2 and creates one in P 1 . Since only an even number of particles is present in P 2 ,b † i b j |Ψ would contain only configurations with an odd number of particles in P 2 , which makes the inner product with Φ I |, ∀I ∈ RAS vanish. In the same mannerb † i b j acting on |Φ I is either zero, if |φ j is unoccupied in the configuration |Φ I , or gives an odd number of particles in P 2 . In this specific excitation scheme,ρ j i = 0, for all pairs of orbitals {i , j }, leaving the amplitudes and the P-space orbitals equations uncoupled. Using the explicit expressions of the Hamiltonian [Eq. (8)] and the operatorD, We can simplify this expression, starting with Now we turn to the chains of six operators in the last term in Eq. (33). The first product of operators is expressed and the second product of operators as The sum of Eqs. (35) and (36) enters Eq. (33), and we see that only chains of four operators remain. Using the fact that v ik jl = v ki lj [Eq. (10)] and that η p q must be evaluated for orbitals which belong to different P i space, {l , k }, Eq. (33) can be rewritten,

d. General RAS scheme
Considering only even excitations provides an efficient and simple way to uncouple the equations of the TD-RASSCF-B method. Nonetheless, it is also possible to consider both even and odd excitations in the configurational space. In the following, we specifically consider a RAS scheme with all successive numbers of particles occupying P 2 from 0 to N max , where N max , defined in Sec. II A, is the highest number of particles allowed in P 2 , see Fig. 2b. Note that N max must fulfill the condition N max ≤ N . For instance, taking N max = 4, we consider the promotion of 0,1,2,3 and 4 particles from P 1 to P 2 . In this way, the configurational space is span by the direct sum of N max + 1 subspaces, Using the expression of the time derivative of the {C I } coefficients, Eq. (14), the time derivative of the one-body density matrix, present in Eq. (29), can be expressed as, We introduce the projector onto the RAS space V asΠ = I∈V |Φ I Φ I |. Using the above expression ofρ j i , we obtain a new formulation of the P-space orbital equation, For |φ i ∈ P 1 and |φ j ∈ P 2 , we note that |Ψ i j belongs to V, with one particle from P 2 being annihilated and one particle in provides configurations with a creation of an additional particle in P 2 , which may lie in V Nmax+1 , not included in V. In this case, Ψ j i |(1 − Π) = 0 and Eq. (40) simplifies to, Using the expression of the Hamiltonian, [Eq. (8)], of the operatorD, [Eq. (13)], and keeping in mind that η k l has only to be determined for pairs of orbitals {l , k } which belong to different P i subspaces, Eq. (41) is equivalent to where the fourth-and sixth-order tensors are defined by Here again, the P-space EOM for the determination of the η k l , Eq. (42), are identical for bosons and fermions [86,87]. These equations are solved to determine the η k l for each pairs of orbitals belonging to different P ispace. The value of η k l is subsequently used to solve the amplitudes equations [Eq. (15)] and to evaluate the time-derivative of the P-space orbitals from Eqs. (30) and (26), as for the case of the even excitation scheme.
The general excitations scheme and the only even excitations schemes were originally introduced in the case of fermions in Refs. [86,87]. We mention that recently Haxton et al. [93] derived a general RAS scheme for fermions, in the sense that the configurational space can be build from of any arbitrary configurations. For both excitation schemes presented in this work, the time-derivative of the coefficients and orbitals are obtained by solving the amplitudes equations, Eq. (15), the Q-space equations, Eq. (26) and the P-space equations Eq. (42) and Eq. (37) for the general RAS scheme and the only even excitation scheme, respectively. In the case of the MCT-DHB method, the amplitude [Eq. (15)] and the Q-space equations [Eq. (26)] are also solved to obtain the timederivative of the wavefunction, see Appendix C. The numerical efficiency to solve these equations scale differently with the number of configurations and the number of orbitals, as detailed in Appendix D. For a given number of orbitals, the TD-RASSCF-B method is more efficient to solve Eqs. (15) and (26), irrespectively of the excitation scheme used. Nonetheless, in the TD-RASSCF-B framework one additional system of equations needs to be solved, namely the P-space equations [Eq. (37) or (42)]. For only even excitations, the number of operations required to obtain the time-derivative of the wavefunction is always smaller in the case of the TD-RASSCF-B method than in the MCTDHB method. In the case of the general excitation scheme, the evaluation of the sixth-order tensor, Eq. (44), requires a significantly large number of operations. Thus, the TD-RASSCF-B method may require more operations than MCTDHB for large values of N max and large numbers of orbitals. As shown in Appendix D, this happens only for large values of N max , for instance for N max > 38 with N = 50 or N max > 909 for N = 1000 bosons. Except for these high excitation schemes, the TD-RASSCF-B method is numerically more efficient than the MCTDHB method, but more importantly the exponential grows of the configurational space with respect to the number of orbitals can be controlled thanks to the RAS Ansatz. In addition, we have shown that the TD-RASSCF equations of motion are the same for bosons and fermions, which means that the TD-RASSCF theory is a general framework including as limiting cases the TD-GP (TD-HF) and the MCTDHB (MCTDHF) theories for bosons (fermions). This result is reminiscent to the work of Alon et al [94] where a unified set of EOM for the MCTDH theory for both bosons and fermions was derived.

III. APPLICATION TO A TIME-INDEPENDENT SYSTEM: GROUND STATE ENERGY
In this section, we consider a system of N = 100 bosons trapped in a 1-dimensional (1D) harmonic potential. Experimentally, quasi-1D systems have been obtained by using a tight confinement in the transversal coordinates, freezing in that way the transversal dynamics of the system [95][96][97][98][99]. In the following, we consider an anisotropic harmonic trap such that the longitudinal frequency (ω x ) is much smaller than the transversal frequency (ω ⊥ ), i.e., ω ⊥ ω x , such that the transverse part of the wavefunction can be assumed to be energetically frozen to the ground state and be integrated out. The resulting 1D Hamiltonian for the N boson system reads, using the unit of length l 0 = h(mω x ) −1 and the unit of energy E 0 =hω x , with m the mass of the particles. Assuming no confinement induced resonances [100], the interaction strength, λ, is related to the 3D s-wave scattering length of the particles, a s , through λ = 2a s l 0 l −2 ⊥ , with l ⊥ the transversal harmonic oscillator length. Experimentally, the 1D interaction strength can be tuned either by controlling the longitudinal and transversal frequencies or using an external magnetic field [9,10].
To solve numerically the EOM of the MCTDH and TD-RASSCF-B theories, the time-dependent orbitals are expanded on a time-independent basis or primitive basis, which consists of a sine discrete variable representation (DVR), see Ref. [66]. We use 101 basis functions in a box [−8, 8] and compare the results with larger basis sets to ensure the convergence of the energies presented in Tables I and II. We numerically integrate the EOM using different integration algorithms, namely the 4th order runge-kutta (RK), the adaptive time-step 5th order RK [101] and the Adams-Bashforth-Moulton (ABM) predictor-corrector as implemented in the Heidelberg MCTDH package [102]. The different integration schemes were tested against each other and we report the results obtained using the ABM integrator to the 7th order, as it is the most efficient. We calculate the GS energy using imaginary time propagation [103] of the EOM and give its energy in units of E 0 .
To assess the accuracy of the GP, MCTDHB and TD-RASSCF-B methods we compare the GS energies and by virtue of the variational principle (see for instance Ref. [104]), the lower the energy the higher the accuracy. First, as a general remark, for any value of λ = 0 the energy obtained with the MCTDHB method systematically decreases with increasing number of orbitals and subsequently for increasing number of configurations, see for instance the first line of Table I where the numbers of configurations are indicated in parentheses. Concerning the TD-RASSCF-B method, for a given excitation scheme the energy also decreases when the number of orbitals is increased. In addition, for a given number of orbitals the energy decreases when we increase the highest number of allowed particles in P 2 , N max . To simplify the following discussion, we introduce some quantities to help the comparison between the MCTDHB and TD-RASSCF-B methods. Firstly, we define the correlation energy as the difference between the energy obtained with a given method and the mean-field GP energy, where E method designates the energy obtained with a given method. By definition, the GP correlation energy is equal to zero and is considered as uncorrelated. We use as a reference, E ref , the correlation energy obtained for the MCTDHB method with 5 orbitals, i.e., where the superscript 5 denotes for the number of orbitals. Using this reference, we can easily compare the results obtained for different numbers of orbitals by expressing the correlation energy in percent of E ref . Secondly, we define the relative correlation energy, E X rel , as the difference between the energy obtained from a MCT-DHB calculation with X orbitals and the GP energy, i.e., This quantity is particularly useful to compare the results of different RAS schemes within a given number of orbitals. Indeed, the TD-RASSCF-B Ansatz, with restrictions on the active space, is an approximation to the MCTDHB wavefunction. Thus, when the correlation energy of a RAS scheme with X orbitals is equal to E X rel the calculation is converged.
We first focus on the results obtained with λ = 0.01, the weakest interaction strength considered and we report the results in Table I. The reference for the correlation energy is E ref = −2.9 × 10 −2 (the GP result is obtained from MCTDHB with a single orbital). Increasing the number of orbitals in the MCTDHB calculations from 2 to 5 allows us to account for more and more of E ref . Specifically we obtain 51%, 78% and 91% of E ref for 2, 3, and 4 orbitals, respectively. The 10% variation of the correlation energy between 4 and 5 orbitals indicates that the results are not fully converged with respect to the number of orbitals, and more than 5 orbitals are required to converge the energy below 10 −3 , see Table I. Unfortunately, the MCTDHB wavefunction with 5 orbitals includes already ∼ 4.6 × 10 6 configurations and using 6 (8) orbitals leads to ∼ 96×10 6 (∼ 26×10 9 ) configurations, far beyond the scope of any practical numerical implementation.
The TD-RASSCF-B method provides more flexibility to describe the wavefunction in the sense that we can choose different RAS schemes, different numbers of orbitals and their partitions into P 1 and P 2 spaces. In Table I we report the results obtained for M = 2 to 8 orbitals with a single P 1 orbital, i.e., M 1 = 1 and M 2 = M − M 1 P 2 orbitals, and a few specific cases of the general RAS scheme. We indicate the excitation schemes with the usual notation. For example -SD denotes that single and double excitations are allowed from P 1 to P 2 . We follow this notation up to -SDTQ56789 and for larger excitations, we just indicate the value of N max (e.g., "-10" means that all excitations from P 1 to P 2 up to 10 are included). For each number of orbitals, when we increase the excitation scheme the energy becomes closer to the MCTDHB result and converges to this latter for the -SDTQ5678 RAS scheme, as indicated by the underlined digits in Table I. Thus, E ref is recovered for the -SDTQ5678 scheme with 5 orbitals, but the expansion of wavefunction includes only 495 configurations, i.e., 9×10 3 times fewer configurations than the MCTDHB expansion for 5 orbitals. It is worthwhile to note that this RAS scheme converged for all number of orbitals, and always for much fewer configurations than with the MCTDHB. The least accurate TD-RASSCF-B calculation, presented in Table I, consists of 2 orbitals and the -SD RAS scheme. The correlation energy includes 98.7% of E 2 rel and interestingly when we increase the number of orbitals from 3 to 5, a similar amount of correlation is obtained (98.8% for all values) in comparison to the respective E X=3,4,5 rel correlation energies. Thus, using the -SD scheme with 5 orbitals 98.8% of E ref is obtained but the TD-RASSCF-B wavefunction includes only 15 configurations while the MCTDHB wavefunction includes more than 4.5 × 10 6 configurations, i.e., ∼ 3 × 10 5 times more configurations. Moreover, the energy difference between the -SD scheme and MCDTHB method is systematically below 5 × 10 −4 , lower than the convergence obtained with respect to the number of orbitals. Concerning the correlation energy of the -SDTQ and -SDTQ56 schemes, we find that they include 99.99% and 99.9999% of E X rel , with X = 2 to 5. It is remarkable that the correlation energy remains almost constant while the difference between the number of configurations between the MCTDHB and TD-RASSCF-B increases exponentially with the number of orbitals. These results show that the correlation energy does not strongly depend on the number of configurations used in the wavefunction expansion, as the configurational space of the -SD RAS scheme increases only from 3 to 15 configurations for M = 2 to M = 5 but captures 98.8% of E X rel , with X = 2 to 5. Thus, the correlation depends more critically on the number of orbitals than the number of configurations used in the calculation. To illustrate this point, we compute the GS energy with 6 to 8 orbitals with the TD-RASSCF-B method, see Table I, and we obtain energies lower than the energy of the MCTDHB with 5 orbitals for all excitation schemes used here. It means that the -SD scheme with 8 orbitals and 36 configurations is more accurate than the MCTDHB method with 5 orbitals and ∼ 4.6 × 10 6 configurations. Moreover, comparing the energies obtained for the -SDTQ5678 and -10 excitation schemes, we can conclude that the GS energy has converged with respect to the number of excitations. Thus, the TD-RASSCF-B method, thanks to the restriction imposed on the configurational space, can provide more accurate results than the MCTDHB method, whoes practical applicability is limited by the exponential growth of the number of configurations.
We also consider RAS schemes with only even excitations (see Table II), for which the numerical effort is always reduced in comparison to the MCTDHB method, see Appendix D. The energy difference between the -D and the -SD schemes is below 1.3 × 10 −6 for all numbers of orbitals, indicating that more than 98.7% of the relative correlation energy is obtained with slightly fewer configurations. The same conclusion holds for the com-parison of the -DQ and -SDTQ schemes with an energy difference below 1.5 × 10 −5 , including more than 99.95% of the relative correlation energy. The number of configurations is slightly smaller in the case of the RAS schemes with only even excitations but the numerical efficiency is better as the P-space EOM, Eq. (37), does not require the update of a sixth-order tensor at each time-step as it is the case of the general RAS schemes, see Eq. (44). For values of N max ≥ 6 and M ≥ 3, the energy converges with respect to N max , as the energy does not change by increasing N max further, but with energy slightly larger than the MCTDHB ones. The energy difference between the -DQ68 scheme and the MCTDHB calculation, with 5 orbitals for both methods, is ∼ 1.1 × 10 −5 and includes 99.96% of E ref . In the case of only even excitations, we also find that for M > 5 the energy for all schemes is below the energy of the best MCTDHB calculation performed. The comparison of the converged -DQ68 and -SDTQ5678 RAS schemes, show that the energy difference remains below 1.5 × 10 −5 for M > 5, which is two orders of magnitude small then the convergence obtain with respect to the number orbitals ∼ 10 −3 .
We perform the same analysis for an interaction strength λ = 0.1 and we obtain, as a reference for the correlation energy, E ref = −1.34, see Table I. This value is much larger than the one obtained previously and can be explained by the stronger interaction between the particles. Indeed, for a stronger interaction strength, the energy of the system is lowered by allowing the particles to occupied higher orbitals, i.e., orbitals leading to higher kinetic and potential energies, such that the interaction energy is reduced. In the mean-field GP theory, this possibility is not possible as only one orbital is used to describe the wavefunction. The orbitals that diagonalized the reduced density matrix and their respective eigenvalues, or population, can be used to characterized the system. If the largest eigenvalue is of the same order as N , the system is condensed [105]. As the GP wavefunction includes a single orbital, it can only describe condensed systems. If more than one eigenvalue is of the order of N , then the system is fragmented (see Ref. [106] and the discussion in Ref. [78]). We find that, indeed, the occupation of the lowest natural orbital in the MCTDHB calculation using 5 orbitals decreases from a population of 99.987% for λ = 0.01 to a population of 99.465% for λ = 0.1. This slightly larger depletion of the condensate has a strong impact on the correlation energy, as the mean-field GP theory provides a less accurate description of the system. We point out that increasing the number of orbitals from 4 to 5 in the MCTDHB calculations gives an energy difference ∼ 10 −1 , see Table I, which means that the energy does not converge below this value. The relative correlation energy E 2 rel , E 3 rel and E 4 rel include 40.1%, 68.8% and 86.7% of E ref , respectively. Starting the discussion with the general RAS schemes, see Table I, we note that to converge to the MCTDHB energies and thus include 100% of the relative correlation energies, E X rel with X = 2 to 5, large values of N max are required. For the -SDTQ5 RAS scheme we find that the correlation energy includes ∼ 97% of E X rel , the -10 RAS scheme includes ∼ 99.9% of E X rel , the -15 RAS scheme includes ∼ 99.997% of E X rel and the -20 RAS scheme is converged with more than 99.9999% of E X rel . These results are obtained irrespectively of the number of orbitals, i.e., X = 2 to 5. Even if large values of N max are used, the expansion of the wavefunction using the -20 RAS scheme includes ∼ 10 4 configurations for 5 orbitals while the MCTDHB wavefunction includes ∼ 4.6 × 10 6 configurations. We also use RAS schemes with only even excitations and we report the results in Table II. Similarly to the case with λ = 0.01, we find that, except for 2 orbitals, the energy does not converge to the MCTDHB energy, irrespectively to the value of N max used. Thus the energies obtained with 3, 4 and 5 orbitals include 99.7%, 99.0% and 98.8% of the respective E X rel with the -20 RAS scheme, which is converged. For similar numbers of configurations the general RAS scheme provides more accurate results but remains more demanding in term of computation, see Appendix D. It is important to keep in mind that the convergence with respect to the number of orbitals is ∼ 10 −1 and a convergence one order of magnitude below is achieved with the -SDTQ5 excitation scheme (Table I) and the -DQ excitation scheme (Table  II). As previously, the configurational space of the MCT-DHB wavefunction becomes unworkable for more than 5 orbitals, but the TD-RASSCF-B method can include more orbitals. At the level of the -SDT scheme and 8 orbitals and for higher excitation schemes with 6 to 8 orbitals, the GS energy is always below the energy obtained with the MCTDHB method with 5 orbitals, see Table I. In the same way, using only even excitations provides more accurate results for excitation schemes higher than -D. As a remark, we obtain an energy ∼ 0.3 below the energy of the MCTDHB method with 5 orbitals by using the -15 RAS scheme with 8 orbitals, see Table I.
These preceding examples show that the TD-RASSCF-B method provides an efficient approach for computing the GS energy of trapped cold atoms. This wavefunction based approach gives access to quantities of interest such as the one and two-body reduced densities and the fragmentation using the population analysis of the natural orbitals. The accuracy was compared with the MCTDHB results and we showed that the TD-RASSCF-B method converges for relatively low excitation schemes. Moreover, the possibility to constrain the growth of the configurational space gives the possibility to use more orbitals than in the MCTDHB calculations and better results were systematically obtained using the TD-RASSCF-B method. This result can be understand as follows, the MCTDHB wavefunction for a small number of orbitals generates a large number of configurations, as all orbitals can be equally populated. The main part of these configurations, however, do not contribute to lower the energy of the system as they describe states with many particles occupying the same spatial orbital, which induced a large interaction energy for a large value of λ. As a limiting case, we know that in the Tonks-Girardeau model [19], obtained for an infinite value of λ, each boson occupies a different orbital. Thus, using a larger number of orbitals in the TD-RASSCF-B method introduces configurations for which a small number of particles occupy a larger number of different orbitals, and thus describes more efficiently the system. This flexibility of the TD-RASSCF-B method of choosing more orbitals opens a new possibility to explore the static properties of trapped cold atoms in systems with hundreds of particles and large numbers of orbitals, which are for the moment beyond the possibility of the MCTDHB method.

IV. APPLICATION TO A TIME-DEPENDENT SYSTEM: DYNAMICS OF BOSONS WITH HARMONIC INTERACTION
As an illustration of an application of the TD-RASSCF-B method to a truly time-dependent problem, we simulate the dynamics of an ensemble of N = 10 bosons in a 1D harmonic trap interacting through a harmonic interaction potential. We consider an initial system of non-interacting bosons, for which the Hamiltonian H 0 reads,Ĥ where we use the units described in Sec. III and the time is expressed in units of t 0 = ω −1 x . The analytical ground state wavefunction and energy are used to ensure the convergence of the imaginary time propagation and, as expected, are the same for all methods. The dynamics is initiated at t = 0 by quenching instantaneously the strength of the two-body interaction, as performed in Ref. [24], leading to the evolution of the system under the action of the following Hamiltonian, with λ the strength of the two-body interaction. This sudden change in the interaction between the bosons leads to a breathing dynamics of the BEC with frequencies Ω n = 2n ω 2 x + 2N λ, with ω x the frequency of the harmonic trap, see Ref. [24]. For positive values of λ the two-body interaction is attractive, while for negative values the interaction is repulsive and leads to unbound dynamics for λ < −ω x /2N . Note that we use the same parameters than in Sec. III for the numerical resolution of the EOM.
A. Breathing dynamics with λ = 0. 1 We first consider the dynamics following a quenching of the interaction strength from λ = 0 to λ = 0.1 [Eq. (50)]. We find that the MCTDHB method with M = 4 orbitals and 286 configurations is numerically exact for the propagation time considered here, i.e., 0 ≤ t ≤ 15, see Fig. 3. The time evolution of the system is characterized by the one-particle density, ρ(x = 0, t), at the center of the trap x = 0 and exhibits a periodic evolution with a frequency ω MCTDHB = 3.46. This value agrees perfectly with the analytical prediction Ω n=1 = 3.46 meaning that the first excited state is mainly responsible for the dynamics. Nonetheless, the discrepancy with a pure cosine function ∼ cos(Ω 1 t) indicates the role of higher excited states with higher harmonic frequencies [24]. The meanfield GP fails to describe the system evolution, even at short time (t < 1) as depicted in Fig. 3 (a) and we obtain a lower frequency ω GP = 3.35. First, we perform TD-RASSCF-B simulations using a single P 1 orbital (M 1 = 1) and M 2 = 3 P 2 orbitals for different RAS schemes reported in Fig. 3 (a). For a short time, i.e., 0 ≤ t ≤ 5, all RAS schemes describe accurately the dynamics of the system, in contrast to the GP result. On the scale of the figure, the convergence to the MCT-DHB result is obtained by using the -SDTQ excitation scheme including 35 configurations, a reduction by a factor of ∼ 8. For a longer time, 10 ≤ t ≤ 15, the -SD RAS scheme substantially differs from the MCTDHB result with a shift in the frequency and a smaller amplitude of the oscillations. The -SDTQ scheme only slightly differs by a smaller amplitude for t > 11.5 and convergence is achieved for the -SDTQ56 excitation scheme with 84 configurations. We also investigate the role of the P 1 orbitals on the accuracy of the computations by considering M 1 = 2 and M 2 = 2, such as the total number of orbitals, M = 4, remains unchanged. The -SD excitation scheme converged for short time, see Fig. 3 (b), and provides a better description of long time dynamics than the -SDTQ scheme used previously [ Fig. 3 (a)] but includes 58 configurations. This number of configurations is similar to the 56 configurations obtained by using the -SDTQ5 RAS scheme with a single P 1 orbital, which differs from the MCTDHB results for t > 11.5 (not shown) while the -SD scheme with 2 P 1 orbitals differs for t > 13.5, see Fig. 3 (b). We obtain converged results for the -SDT RAS scheme with 90 configurations. In the previous section, we showed that the RAS schemes with only even excitations provide accurate results for GS energy and reduce the numerical effort (see Appendix D). We apply the -D and -DQ schemes with M 1 = 1 and M 1 = 2 Fig. 4 (a) and (b), respectively, and keep M = 4. For both M 1 = 1 and M 1 = 2, the results do not converge to the MCTDHB results and do not significantly improve for larger excitation schemes. For M 1 = 1, Fig. 4 (a), the results obtained with the -D and the -DQ schemes are in very good agreement with the MCTDHB results concerning both the frequency and the amplitude and start to deviate only for t > 11. In both cases, the number of configurations used to expand the wavefunction is substantially smaller than the expansion of the MCTDHB wavefunction with 7 and 22 configurations, respectively.
When we use 2 orbitals in the P 1 -space both -D and -DQ provide the same results for the dynamics, see Fig.  4 (b). For short time dynamics, the results are similar to the ones obtained previously with M 1 = 1, but for a longer time, the oscillations remain in phase with the MCTDHB result, only the amplitude deviates for t > 13. The wavefunction of the -D scheme includes 38 configurations. Thus, the TD-RASSCF-B method provides an access to describe accurately the dynamics of the interacting system, while the mean-field GP theory failed even for short time. We obtain a good agreement in comparison to the MCTDHB method with a substantial reduction of the configurational space, using few tens instead of the few hundreds of configurations with the MCTDHB method. Moreover, the different parameters of TD-RASSCF-B method which define the wavefunction can be used to converge the results to the MCTDHB calculations. The implications of this reduction on the CPU time are discussed at the end of this Section.
B. Breathing dynamics with λ = 0. 5 We pursue the illustration of the TD-RASSCF-B method by considering a quenching from λ = 0 to λ = 0.5 [Eq. (50)]. This interaction strength was used in Ref. [24] to benchmark the MCTDHB method. For the time interval considered here, 0 ≤ t ≤ 15, we find that the result obtained with the MCTDHB method using 8 orbitals and 19448 configurations is numerically exact, in agreement with Ref. [24]. As previously, the one-body density exhibits oscillations as a function of the time with a period ω MCTDHB = 6.63, in perfect agreement with the analytical frequency Ω n=1 = 6.63. In comparison to the previous results, the shape of the oscillations indicates that the role of higher excited states is stronger as a large deviation from a simple cosine function, ∼ cos(Ω 1 t), is obtained. This is not surprising since the value of λ is now 5 times larger than the one of the previous example. Along with the MCTDHB result we report, for comparison, the result obtained using the mean-field GP theory, see Fig. 5 (a), which strongly deviates from the MCT-DHB result for t > 0.3 with a larger amplitude in the oscillations and gives a lower frequency for the oscillations, ω GP = 6.32. We start the TD-RASSCF-B simulations using M = 8 orbitals with one single P 1 orbital and M 2 = 7 P 2 orbitals, Fig. 5 (a). For times between 0 and 2.5, the excitation schemes larger or equal to -SDTQ provide an accurate description of the dynamics, while the -SDTQ scheme includes 330 configurations in the wavefunction, i.e., a factor of ∼ 60 less than the MCT-DHB. For longer times, the -SDTQ5678 RAS scheme with 6435 configurations is required to accurately describe the MCTDHB results, the lower excitation schemes give substantially different results. To improve the results of the TD-RASSCF-B method, we increase the number of orbital in P 1 and keep constant the total number of orbitals, M = 8. In Fig. 5 (b), (c) and (d) we report the results with M 1 = 2, 3 and 4 P 1 orbitals, respectively. In the case of M 1 = 2, the -SDTQ56 RAS scheme, with 5412 configurations, provides a very accurate description of the dynamics for 0 ≤ t ≤ 15. For lower excitation schemes, the results are accurate for 0 ≤ t ≤ 3 and using the -SDTQ scheme, the frequency of the oscillation at a longer time are correctly obtained, see Fig. 5 (b). For M 1 = 3, the results in Fig. 5 (c) show that the -SDTQ5 RAS scheme including 6882 configurations converged to the MCTDHB results, while the -SDT scheme with 2276 configurations gives the correct frequency for the oscillation but with a smaller amplitude. Finally, Fig.  5 (d) report converged results for M 1 = 4 using the -SDT schemes, which includes 5216 configurations. Nonetheless, the -SD scheme with 2816 configurations is accurate for the considered time of propagation in comparison to the MCTDHB result.
Using different numbers of P 1 orbitals and different RAS schemes points out that the number of configurations required to accurately describe the evolution of the system change substantially. For instance, a similar accuracy is achieved for the -SDTQ5678 scheme with M 1 = 1, the -SDTQ56 scheme with M 1 = 2, the -SDT scheme with M 1 = 3 and the -SD scheme with M 1 = 4. For each case, a smaller amplitude of the oscillations is observed in comparison to the MCTDHB for t > 12. But the numbers of configurations used in the wavefunction expansions are 6435, 5412, 2276 and 2816, respectively. Thus, for a comparable accuracy, the number of configurations can be divided by a factor ∼ 2 by choosing adequately the size of the P 1 -space and a factor ∼ 10 in comparison to the MCTDHB configurational space. The reduction of the configurational space impacts strongly the required CPU times of the simulations. For instance, the -SDTQ5678 RAS scheme with M 1 = 1 took ∼ 18.7 CPU hours, while the -SDTQ56 scheme with M 1 = 2, the -SDT scheme with M 1 = 3 and the -SD scheme with M 1 = 4 took ∼ 19.8, ∼ 6.9 and ∼ 11.1 CPU hours, respectively, on a 2.4 GHz Intel E5-2680 CPU. Using these RAS schemes, the CPU time is drastically reduced in comparison to the ∼ 113.4 CPU hours on a 2.5 GHz Intel E5-2680 CPU needed to perform the MCTDHB simulation, but the dynamics is accurately described. Moreover, the -SDTQ5 RAS scheme with M 1 = 3 converged to the exact solution with ∼ 4 times less CPU time. Albeit the drastic reduction in the size of configuration space, the CPU time needed to perform the TD-RASSCF-B calculations is also substantially reduced.
We briefly summarize the findings concerning the dynamical evolution of trapped atoms after a quenching of the interaction strength of an attractive harmonic interaction. In the case of λ = 0.1, the MCTDHB theory converged to the numerically exact result for 4 orbitals. The TD-RASSCF-B method using different size of the P 1space and different RAS schemes can accurately describe the dynamics of the system characterized by ρ(x = 0, t) with ∼ 4 times fewer configurations. Moreover, converged TD-RASSCF-B results were obtained with sub-stantially less configurations, for instance with M 1 = 2 and the -SDT excitation scheme. In the case of λ = 0.5, the exact solution was obtained with 8 orbitals using the MCTDHB method leading to 19448 configurations. A larger number of orbitals is needed for the stronger interaction between the particles. Accurate results were obtained with the TD-RASSCF-B method, reducing by a factor ∼ 10 the size of the configurational space and converged results were obtained with 4 times less configurations, for instance considering M 1 = 4 and the -SDT RAS scheme. Moreover, all calculations performed with the TD-RASSCF-B method were better than the meanfield GP theory, which failed to describe both scenarios.

V. CONCLUSION AND OUTLOOK
In this work, we presented a general formalism for the time-dependent restricted active-space self-consistent field (TD-RASSCF) method, which includes the first derivation obtained for fermions (TD-RASSCF-F) [86][87][88] and extended it for systems of spinless bosons (TD-RASSCF-B). This TD-RASSCF-B method includes, as limiting cases, the (TD)-GP and the MCTDHB theories and provides a way to tackle the exponential growth of the configurational space in the MCTDHB method. The EOM were derived for two families of RAS schemes, which give the possibility to restrict the full-configurational description of the MCTDHB wavefunction. Through a set of numerical examples, we have shown that the method can provide an accurate description of the static properties of the ground-state of the system. In the case of hundreds of particles, the method can lead to results beyond the reach of the MCTDHB method, providing a better accuracy by including more orbitals while constraining the number of configurations. In this sense, the TD-RASSCF-B method paves the way for numerical investigation of intermediate system sizes with a few tens to hundreds of bosons with a better accuracy than what was possible with the MCTDHB method. We also provided a comparison between the MCTDHB and TD-RASSCF-B method in the case of breathing dynamics induced by a sudden quenching of the interaction strength with two different initial conditions. As for the MCTDHB method, the TD-RASSCF-B method does not have any restriction on the choice of the two-body interaction potential used, as was illustrated by the use of a non-contact harmonic interaction between the bosons. Using as a reference the numerically exact result obtained from the MCTDHB method, we showed that the TD-RASSCF-B method is always more accurate than the mean-field TD-GP theory to describe the dynamics of the system. Moreover, using different RAS schemes and partitions of the P-space we obtained very accurate results for substantially less configurations and thus less CPU time than with the MCTDHB method. This reduction of the configurational space can be efficiently exploited to solve numerically the TDSE beyond the mean-field approach. Dynamical effects such as the four wave mixing (FWM) process [107] used to produce correlated atoms beams [108] or the dynamics of bright [109,110] and dark [111,112] solitons can be investigated ab-initio beyond the mean-field TD-GP theory. In addition, the dynamics induced by a time-dependent Hamiltonian can also be explored using the TD-RASSCF-B method, such as in the case of periodically driven optical lattices [113,114].

ACKNOWLEDGMENTS
The authors are indebted to Dr. Haruhide Miyagi for useful discussions. This work was supported by the ERC-StG (Project No. 277767-TDMET), and the VKR center of excellence, QUSCOPE.

APPENDICES
In these Appendices we provide a brief description of the implementation the TD-RASSCF-B method. The implementation is rather similar for bosons and fermions in the sense that the set of equations that we have to solve, i.e. Eqs. (14), (24) and (37) or (42), only depend on the type of particles trough the creation and annihilation operators and the set of configurations {|Φ I }.

Appendix A: Compact representation of the wavefunction
Our implementation is based on the general mapping of bosonic operators in Fock space introduced in Ref. [115] and implemented, for instance, for bosons [116] and fermions [117] in the framework of multi-configurational TD methods. Assuming M orbitals and N bosons, the configurations are expressed using the occupation number formalism, where |n 1 , n 2 , · · · , n M , represents a configuration with n 1 bosons in the orbital |φ 1 , n 2 bosons in the orbital |φ 2 , etc. Such a configuration is indexed by an unique integer J defined as, Thus, for each configuration we store its complex coefficient C J in an array according to the index J provided by the above mapping. In Ref. [116], this mapping was employed to avoid the storage of the configuration vectors |n 1 , n 2 , · · · , n M , which can be prohibitively memory consuming in the case of the MCTDHB method. Thus, to access the coefficient of the configurations, a set of Mnested loops over the occupation number n i is used to span the full configurational space and to compute, using Eq. (A1), the respective indexes. In the case of the TD-RASSCF-B method only a selected number of configurations are used to expand the wavefunction, and the scheme of Ref. [116] is not readily applicable in the sense that we want to avoid explicit use of the full configurational space. Instead we follow a different strategy for indexing the RAS configurations. We introduce M 1 and M 2 the number of orbitals in the P 1 and P 2 spaces, respectively. For V 0 , i.e., configuration with particles only in P 1 , see Eq. (38), we enumerate all possible configurations, evaluate their index, using Eq. (A1), and store them. Then for the excited configurations, i.e., configurations with one or more particles in P 2 , we introduce the excitation number n exc which is equivalent to the number of particles in P 2 . The number of remaining particles in P 1 is N − n exc . For each excitation we enumerate the configurations of N − n exc particles in the P 1 orbitals and compute their index. The same is performed for the configurations with n exc in P 2 orbitals. Then the permanents for the total number of bosons are obtained by combining them, where |n nexc is an array of dimension M × (N nexc ) the total number of configurations obtained from arranging the N b − n exc (n exc ) particles in the P 1 (P 2 ) orbitals. We evaluate the indexes of the configurations resulting from the P 1 subsystem, J nexc P1 , and from the P 2 subsystem, J nexc P2 , applying Eq. (A1) for the subsystems, separately. The index of the configuration with the total number of bosons is build as a three components array defined by {J nexc P1 , J nexc P2 , n exc }, which stores the position of the configuration in the configurational vector. This scheme is thus applied for all excitations, n exc , included in the RAS scheme. With this storage or construction of the wavefunction, for a given configuration we have access to its coefficient in the following way. (i) We evaluate the index J nexc P1 for the P 1 subsystem, (ii) we evaluate the index J nexc P2 for the P 1 subsystem, (iii) we know or evaluate the excitation, n exc , of the configuration (iv) we access to the index of the coefficient which is stored and the three component array at position {J nexc P1 , J nexc P2 , n exc }. This scheme has, as a draw back, the requirement to store the list of occupation numbers and the indexes to be efficient for numerical evaluation. But the idea behind the TD-RASSCF-B method is to reduced the size of the configurational space, which makes such a storage manageable for applications done so far.

Appendix B: Applying operators in second quantization
In second quantization, the action of the Hamiltonian of Eq. (8) on the wavefunction requires the application of annihilation and creation operators and multiplication by the matrix elements of the one-and two-body operators. To know the action of the Hamiltonian, we first need to know the action of the creation-annihilation operators [116]. Concerning the one-body term, we have, b † i b j |n 1 , · · · , n j , · · · , n i , · · · , n M = √ n j √ n i + 1|n 1 , · · · , n j − 1, · · · , n i + 1, · · · , n M .
For a full-configurational wavefunction, the resulting configuration belongs to the configurational space and its index can be determined as described in Appendix A. Now, considering the action on the wavefunction, with |Φ I the new configuration resulting from the action of b † i b j on the initial configuration |Φ I . The result can be interpreted as a reordering of the configuration in the wavefunction, as in Eq. (B2) or inversely to a reordering of the coefficients with a new factor ( √ n j √ n i + 1) if the configuration are reorganized in the initial order, i.e., In the basis of the configurational states, {|Φ I }, the wavefunction is characterized by its coefficients only, and is stored as a vector. The new set of coefficients, {C I }, resulting from the action of b † i b j , is obtained as, In practice, we apply b † i b j on the bra Φ I |, which provides a new configurational state Φ J |. The configurational vectors are orthonormal, and only the configuration |Φ J in the sum remains from the projection. Thus, we evaluate the index of Φ J | to directly access the coefficient C J , i.e., The action on the total wavefunction is obtained by repeating these operations for each configuration, providing a new coefficient vector. In the case of the RAS wavefunction, the configuration obtained from the successive application of the annihilation and creation operators may not belong to the configurational space. Nonetheless, the scheme applied above can be applied with, in addition, a test to check if the resulting configuration remains in the RAS space. This naive approach can be easily improved thanks to the representation of the wavefunction used and its indexing (see Appendix A). The orbitals {i, j}, on which the operators b † i b j act, can (i) belong to P 1 only {i , j }, (ii) belong to P 2 only {i , j }, or belong to P 1 and P 2 , (iii) {i , j } and (iv) {i , j }. For (i) and (ii) the excitation, or the number of particle in P 2 orbital, do not change in the resulting configuration. For the situation (iii) one particle is removed from P 2 and added in P 1 and the opposite happens for (iv). The excitation of the final configuration is thus known without counting the number of particles in P 2 , which is required to determine the index of the configuration. Moreover, to always remain in the RAS configurational space, the case (iv) is never applied to the configuration with the maximum excitation allowed for the general RAS scheme (Sec. II B 2 d), and both (iii) and (iv) are not used for the scheme with only even excitations (Sec. II B 2 c). The action of the one-body operator of the Hamiltonian is now straightforward, the coefficient vectors obtained by applying the b † i b j operators are multiplied by the corresponding matrix element h i j [Eq. (9)] and summed for each couple of {i, j}, with the restriction mentioned above for the RAS wavefunction. The two-body operator, see Eq. (10), included in the Hamiltonian of Eq. (8) and the fourand six-order tensors specific to the RAS schemes [Eqs. (43) and (44)] can be evaluated using the same strategy as the one detailed for the one-body operator. We mention that using the commutation relation for bosonic creation and annihilation operators can substantially reduced the numerical cost. For instance, if we consider the chain of operators

Appendix C: Numerical implementation for the TD-RAS equations
The EOM for the TD-RASSCF-B and F methods, Eqs. (14), (24) and (37) or (42), are solved to obtain the time derivative of the coefficients and orbitals. The main difference with the MCTDHB and F methods results from the evaluation of the matrix elements η j i . These elements are evaluated from Eqs. (37) or (42) depending of the RAS scheme used, but both are solved in the same way. We recall that the matrix η, with elements η j i , is anti-hermitian and thus η i j = −(η j i ) * , and i and j hold for orbitals of the P 1 and P 2 subspace, respectively. Writing down Eqs. (37) or (42) for any set of orbitals {i , j }, provides a system of M 1 × M 2 linear equations with M 1 × M 2 unknowns. Introducing a composite index for the couple of {i , j }, this system can be written in a matrix form, where the matrix A, of dimension ( contains the values of A l j k i or ζ l j k i , for all set of {i , j } and {l , k }. The vector B, of dimension (M 1 ×M 2 ), contains the r.h.s. of Eq. (37) or (42) for each set of {i , j } and the vector X with the same dimension as B contains the unknown values of η j i and the matrix elements of the one-body operator. The system of linear equations, Eq. (C1), can be solved using a standard numerical routine included, for instance, in the LAPACK library [118]. The values for η j i are trivially obtained from the elements of the vector X, After evaluating the matrix elements η j i , the time derivative of the coefficients {Ċ I } can be computed from Eq. (15) and the contribution of the P-space orbitals to the time derivative of the orbitals is obtained from, It remains to evaluate the contribution from the orbitals of the Q-space, i.e.Q|φ i from Eq. (24). This latter can be expressed in a matrix form, with ρ the one-body reduced density matrix,Ẋ a vector collecting the time derivative of the orbitals,h andW are both vectors with elements h(r, t)|φ i and jlkŴ k l |φ j ρ jl ik , respectively. To obtain Eq. (C4), we used the fact that ρ commutes with the projectorQ, as easily seen from the equalityQ =1 −P . Multiplying on the left by the inverse of the one-body density matrix, ρ −1 , we have, for the orbital |φ i , The right hand side of the above equation is similar to the one that is solved in the MCTDH-based methods [64,65,67] and we follow the numerical implementation used for the MCTDH method [119] to avoid singularities in the inverse of the one-body reduced density matrix and in Eq. (C1) for the matrix A, as well as for the projector onto the P-space orbitals.
Appendix D: Numerical efficiency of the method Comparing the efficiency between different methods is a difficult task as it depends of the specific implementation and integration schemes used. Nonetheless, we can roughly estimate the number of operations required to evaluate the time derivative of the orbitals and coefficients and compare the MCTDHB and TD-RASSCF-B methods in this way. We denote by N grid the number of grid points that are used to describe the timedependent orbital in the time-independent basis, usually a DVR [120], which is the same for both methods. Starting with the MCTDHB method, at each evaluation of the time derivative the matrix elements of the two-body operator v ij kl [Eq. (10)] and the two-body reduced density matrix ρ jl ik [see text above Eq. operations, respectively. The total cost is thus, approximatively, 2M 4 (N 2 grid +V FCI ). Considering now the case of the TD-RASSCF-B method. The evaluation of the matrix elements of two-body operator and the calculation of the Q-space equations for the time derivative of the orbitals require the same number of operations as with the MCTDHB method, i.e., M 4 N 2 grid operations for each. The evaluation of the time derivative of the coefficients and the matrix elements of the two-body reduced density matrix scale as M 4 V, with V the size of the configurational RAS space. In addition, we also need to solve the P-space equations, which requires M 4 operations for excitation schemes with only even excitations and M 4 V Nmax for the general RAS scheme, with V Nmax the number of configuration including N max particles in the P 2 -space. The total number of configurations included in the RAS wavefunction for the general excitation scheme is evaluated using Eq. (5) and V Nmax is the last term of the summation. The dimension of the configurational space including only even excitations can be evaluated in a similar way, Combining the results for the general RAS scheme, the number of operations required to evaluate the time derivative of the coefficients and orbitals scales as 2M 4 (N 2 grid + V + M 2 V Nmax ) and in the case of only even excitations it scales as 2M 4 (N 2 grid + V + 1/2). To compare the numerical cost between the MCTDHB and TD-RASSCF-B methods, we can introduce ∆(Op), the difference between the MCTDHB and TD-RASSCF-B op-erations to remove the constant number of operation resulting from N grid , ∆(Op) = 2M 4 (VFCI − V − 1/2) even excitation, From the expression of ∆(Op), a positive value represents a computational gain with the TD-RASSCF-B in comparison to the MCTDHB method, while a negative value is obtained when the MCTDHB method is more efficient. In the case of a scheme with only even excitations, ∆(Op) is proportional to the size difference of the MCTDHB and TD-RASSCF-B configurational spaces and is always positive, which means that the TD-RASSCF-B method is always more efficient. In the case of the general excitation scheme the six-order tensor of the P-space equations, Eq. (44), can provide an overhead for the computation. To illustrate the computational efficiency we evaluate ∆(Op) for 10, 50 and 100 bosons in M = 2 to 8 orbitals, see Fig. 6. For the TD-RASSCF-B, we consider the case of a single P 1 orbital and M − 1 orbitals in P 2 . The case with only even excitations reduces the computational cost almost exponentially for increasing number of orbitals for any number of particles, which results from the efficiency of solving the P-space equation. In the case of the general RAS scheme, there is always a value of N max which leads to more operations in the TD-RASSCF-B than in the MCTDHB method, due to the evaluation of the six-order tensor in the P-space equation. But as shown in Fig. 6, this value is rather large, i.e. N max = 6 for 10, N max = 40 for 50 particles and N max = 90 for 100 particles for the schemes depicted in Figs. 6 (a), (b) and (c), respectively.

Q-space
(Virtual orbitals) a, b, c, ... P 1 -space P 2 -space excitation (P + Q ) -space p, q, r, ... Figure 1. Division of the single-particle Hilbert space within the TD-RASSCF-B framework. The P-space orbitals, used to expand the wavefunction, are divided in a P1 and P2 space between which the particles can be excited through a specific scheme chosen at will. The orthogonal complement of virtual (i.e., unoccupied) orbitals is referred to as the Q-space. The indexes i, j, k, · · · are used to label the orbitals of the P-space, the indexes a, b, c, · · · to label the orbitals of the Q-space and the indexes p, q, r, · · · are used for orbitals in either the P-or Q-space.       particles. The diamond symbols depict the results for the case of RAS schemes with only even excitations and the number of operations is always fewer than for the MCTDHB method. The circles indicate the results for general RAS schemes, and for high excitation schemes and large numbers of orbitals number of operations can be larger than the one of the MCTDHB method. Note that the axis of the general RAS scheme is on the left while the axis of the RAS scheme with only even excitations in on the right.