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

Camille Lévêque, Lars Bojer Madsen

► To cite this version:

Camille Lévêque, Lars Bojer Madsen. Time-dependent restricted-active-space self-consistent-field theory for bosonic many-body systems. New Journal of Physics, 2017, 19 (4), pp.043007. 10.1088/1367-2630/aa6319. hal-03849394

HAL Id: hal-03849394

https://hal.science/hal-03849394v1

Submitted on 13 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

22 December 2016

REVISED

8 February 2017

ACCEPTED FOR PUBLICATION

27 February 2017

PUBLISHED

6 April 2017

Original content from this
work may be used under
the terms of the Creative
Commons Attribution 3.0
licence.

Any further distribution of
this work must maintain
attribution to the
author(s) and the title of
the work, journal citation
and DOI.

<figure><img src="image_3.png" /></figure>

PAPER

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

Camille Lévêque and Lars Bojer Madsen

Department of Physics and Astronomy, Aarhus University, DK-8000 Aarhus C, Denmark

E-mail: bojer@phys.au.dk

**Keywords:** many-body physics, ultra-cold atoms, Bose–Einstein condensate, time-dependent wavefunction based many-body theory

Abstract

We develop an *ab initio* time-dependent wavefunction based theory for the description of a many-body system of cold interacting bosons. Like the multi-configurational time-dependent Hartree method for bosons (MCTDHB), the theory is based on a configurational interaction *Ansatz* for the many-body wavefunction with time-dependent self-consistent-field orbitals. The theory generalizes the MCTDHB method by incorporating restrictions on the active space of the orbital excitations. The restrictions are specified based on the physical situation at hand. The equations of motion of this time-dependent restricted-active-space self-consistent-field (TD-RASSCF) theory are derived. The similarity between the formal development of the theory for bosons and fermions is discussed. The restrictions on the active space allow the theory to be evaluated under conditions where other wavefunction based methods due to exponential scaling in the numerical effort cannot, and to clearly identify the excitations that are important for an accurate description, significantly beyond the mean-field approach. For ground state calculations we find it to be important to allow a few particles to have the freedom to move in many orbitals, an insight facilitated by the flexibility of the restricted-active-space *Ansatz*. Moreover, we find that a high accuracy can be obtained by including only even excitations in the many-body self-consistent-field wavefunction. Time-dependent simulations of harmonically trapped bosons subject to a quenching of their noncontact interaction, show failure of the mean-field Gross-Pitaevskii approach within a fraction of a harmonic oscillation period. The TD-RASSCF theory remains accurate at much reduced computational cost compared to the MCTDHB method. Exploring the effect of changes of the restricted-active-space allows us to identify that even self-consistent-field excitations are mainly responsible for the accuracy of the method.

**1. Introduction**

Since the first realizations of Bose–Einstein condensates (BEC) [1–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–6] trapped in various potential shapes [7, 8] and dimensions [9], with tunable inter-particle interactions [10, 11], and to provide a controllable transition from a few- to a many-particle system. Such a detailed control of cold atom systems has opened the possibility to simulate various physical systems [12] from solid-state physics [13] to black-holes analogs [14] through matter–light interaction [15] and electrons dynamics in molecules [16].

Various theoretical models [17] 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 [18, 19], 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 one spatial dimension and in the limit of infinite interaction strength, the Tonks-Girardeau model use the Fermi-Bose mapping to map the wavefunction of
bosons into a fermionic wavefunction of non-interacting fermions with frozen parallel spins [20]. This mapping provides the exact solution for the ground-state of the system for arbitrary trapping potentials and remains valid also for the excited states, as well as non-equilibrium solutions also for any external potential [21]. In the case of non-interacting bosons or more generally in the Gross-Pitaevskii (GP) limit, i.e., $N \to \infty$ and $N\lambda = \text{constant}$, with $N$ the number of bosons and $\lambda$ the interaction strength, the GP equation or its time-dependent (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 [22], while an other model considers a harmonic interaction potential [23, 24]. 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 time-dependent interaction potential for bosons [25, 26] and fermions [27].

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 [28]. 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 Bogoliubov theory [29–31]. In the specific case of periodic trapping potentials, such as optical lattices [9], for weak contact interactions and deep lattices the Bose-Hubbard model (BHM) [32] 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 [33]. A more general and efficient numerical approach than the BHM to deal with optical lattices is the density-matrix renormalization group method [34–36] based on the matrix product states *Ansatz* [37]. 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 [38–40]. 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) [41] and the diffusion Monte Carlo (DMC) [42–45] methods, which used a Bijl-Jastrow decomposition of the wavefunction [46, 47], 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 Esry [48], 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* [49] 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 [50]. 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 [51, 52] and subsequently extended to describe electronic wavefunctions in atoms and molecules [53]. 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 [54].

Over the past decade, numerous numerical methods were developed [55–60] to tackle the problem of *time-dependent* multi-electron dynamics induced by laser pulses that are strong or short or both [61–63]. 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 [64–67] variationally optimizes a set of time-dependent orbitals and CI coefficients, following the idea of the multi-configurational time-dependent Hartree (MCTDH) method [68, 69], originally introduced to describe molecular dynamics. The MCTDHF method was extended to identical bosons, within the framework of the MCTDH for bosons MCTDHB [70], in which the indistinguishability is taken into account using permanents instead of Slater determinants. Further
developments include the case of particle mixtures of different type of bosons and fermions [71, 72], dynamics described by Hubbard Hamiltonians [73] and bosons with internal degrees of freedom [74]. 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 [75–77] was extended to systems of bosons and mixtures of them [78, 79]. 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 [80].

The MCTDHB method sheds new light on the dynamics of trapped cold atoms, especially when fragmentation occurs and more than one orbital is populated—a situation which cannot be described by the TD-GP theory. Fragmentation occurs in different systems such as during the dynamics at a Josephson junction [81], which is a universal phenomenon [82], and cannot 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 [83] for large barrier height between the two wells and the GP theory fails to describe the variance of position and momentum operators [84]. Multiconfigurational methods are also required to accurately describe the formation and dynamics of fragmented states with repulsive or attractive interactions between the particles [85–87] and tunneling of a many boson system to open space [88] or tunneling of trapped vortices [89]. 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 wavefunctions 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 theory which maintains the appealing properties of the time-
dependent orbitals based methods mentioned above, but is free from the exponential scaling problem. One such
theory uses the concept of a restricted active-space (RAS), well-known in quantum chemistry, where it was
applied with time-independent molecular orbitals [90]. The RAS based method was successfully extended to
time-dependent orbitals in the time-dependent RAS self-consistent-field method (TD-RASSCF) to deal with
electron dynamics in atoms [91, 92]. Introducing a RAS scheme by fixing the promotion of the electrons
between three sets of orbital spaces can considerably reduce the number of configurations. In addition, the
theory includes, as limiting cases, the TD Hartree-Fock (HF), the TD complete active space self-consistent field
(TD-CASSCF) [59] and the MCTDHF frameworks. Successful applications of the TD-RASSCF method include
calculations of the ground-states (GS) of atoms, and time-dependent dynamics in the presence of strong laser
fields to describe, for instance, high-order harmonic generation [91–93]. The aim of this work is to extend the
TD-RASSCF theory 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 [91, 92, 94] as TD-RASSCF-F to avoid any confusion
concerning the particles considered. 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, i.e., bosons
or fermions, 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 very
accurate results, in comparison to MCTDHB, while the expansion of the wavefunction is considerably reduced.
Moreover, the MCTDHB accuracy can be overtaken by using large numbers of time-dependent orbitals while
the number of configurations remains small thanks to the considered RAS. The TD-RASSCF-B theory can be
evaluated in situations where other accurate wavefunction based methods cannot. In addition, and very
importantly, the theory allows an identification of the types of excitations that are important in obtaining
accuracy. For the GS calculations, it is hence important to allow a few particles to move in many orbitals,
representing physically a situation where a small number of particles are outside the condensate, but move
relatively fast. For the GS we find that accurate results can also be obtained by including only even excitations in
the RAS scheme. Insights, that may seem natural, and are now supported by ab initio theory. The investigation of
the breathing dynamic of a BEC confined to a harmonic trap illustrates how the TD-GP theory fails to describe
the time-evolution of the system within a fraction of a harmonic oscillation period, while various examples of
the TD-RASSCF-B method either qualitatively or quantitatively reproduce the exact dynamics obtained using
the MCTDHB method, depending of the choice of the excitation scheme. In particular, we find that very
accurate results can be obtained retaining only even excitations in the RAS scheme, again an example of physics
insight obtained by the flexibility of the TD-RASSCF-B theory.
The paper is organized as follows. In section 2.1 we introduce the TD-RASSCF-B Ansatz for the wavefunction and in section 2.2 we derive the equations of motion for the set of coefficients and orbitals, and highlight the similarity of the formulation for bosons and fermions. In section 3 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 breathing dynamics following a sudden quench of the two-body interaction in section 4. Finally, in section 5 we conclude and provide perspectives to future work. In the appendices A to E, we provide the key ingredients for the numerical implementation of the method, discuss the numerical effort in comparison to the MCTDHB method and provide additional illustrative examples.

## 2. Theoretical framework

### 2.1. 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:

$$i\hbar \frac{\partial}{\partial t} |\Psi(t)\rangle = \hat{H}(t) |\Psi(t)\rangle, \qquad (1)$$

with $\hat{H}(t)$ the many-body Hamiltonian of the system and $|\Psi(t)\rangle$ the $N$-particle wavefunction. Hereafter we set $\hbar = 1$, unless explicitly specified. We can approximate the wavefunction using linear combinations of suitably symmetrized sets of products of time-dependent single-particle functions $\{\phi_i(\mathbf{r}, t)\}$. In the following, the single-particle 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 multi-configurational wavefunction is constructed by taking into account all the possible arrangement of the particles in the given orbitals, each arrangement being called a configuration $|\Phi_I(t)\rangle$,

$$|\Psi(t)\rangle = \sum_{I \in \mathcal{V}_{\text{FCI}}} C_I(t) |\Phi_I(t)\rangle. \qquad (2)$$

This Ansatz converges to the exact wavefunction when the number of orbitals increases to infinity. The configurational space, i.e., the number of coefficients $\dim(\mathcal{V}_{\text{FCI}})$, 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 $\mathcal{V}_{\text{FCI}}$ can be evaluated as,

$$\dim(\mathcal{V}_{\text{FCI}}) = \binom{N + M - 1}{N} = \frac{(N + M - 1)!}{N!(M - 1)!}. \qquad (3)$$

In the case of the TD-RASSCF-B method, we introduce two orbital spaces, $\mathcal{P}_1$ and $\mathcal{P}_2$, such that $M_1 + M_2 = M$, with $M_1$ and $M_2$ the number of orbitals in $\mathcal{P}_1$ and $\mathcal{P}_2$, respectively (see figure 1). The $\mathcal{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 $\mathcal{P}_2$-space, particles are promoted from $\mathcal{P}_1$ to $\mathcal{P}_2$ according to a specific excitation scheme, which restricts the maximal number of particles promoted from $\mathcal{P}_1$ to $\mathcal{P}_2$. 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,

$$|\Psi(t)\rangle = \sum_{I \in \mathcal{V}} C_I(t) |\Phi_I(t)\rangle, \qquad (4)$$

where the configurations are drawn from the space $\mathcal{V}$ subject to restrictions. To evaluate the size of the space $\mathcal{V}$, we introduce $N_{\text{max}}$ the highest number of bosons in $\mathcal{P}_2$ and consider a RAS scheme allowing all occupations of $\mathcal{P}_2$ from 0 to $N_{\text{max}}$. The dimension of the space $\mathcal{V}$ of equation (4), is then given by

$$\dim(\mathcal{V}) = \binom{N + M_1 - 1}{N} + \sum_{k=1}^{N_{\text{max}}} \binom{k + M_2 - 1}{k} \binom{(N - k) + M_1 - 1}{N - k}. \qquad (5)$$

Here, the first term is the total number of configurations obtained with the $N$ bosons in the $M_1$ orbitals and no particle in $\mathcal{P}_2$. The second term in equation (5) takes into account the configurations resulting from the excitation of $k$ bosons in $\mathcal{P}_2$, with $1 \le k \le N_{\text{max}}$. The total number of configurations with $k$ bosons in $\mathcal{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 has some interesting properties. First, if only $\mathcal{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 equation (5). Note that this is also true for $M_2 \ne 0$ and
<figure><img src="image_4.png" /><figcaption>**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 *P*₁ and *P*₂ 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.</figcaption></figure>

$N_{\text{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 to the 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 time-dependent coefficients {$C_i(t)$} and orbitals {$|\phi_i(t)\rangle$}. 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 {$\dot{C}_i(t)$} and {$|\dot{\phi}_i(t)\rangle$}. 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 figure 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 time-dependent. We can define $\hat{P}$ and $\hat{Q}$, the time-dependent projectors onto the subspaces *P* and *Q*, respectively, with the property $\hat{P} + \hat{Q} = \hat{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,

$$|\dot{\phi}_i(t)\rangle = (\hat{P} + \hat{Q})|\dot{\phi}_i(t)\rangle = \hat{P}|\dot{\phi}_i(t)\rangle + \hat{Q}|\dot{\phi}_i(t)\rangle, \quad (6)$$

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 time-derivative of the expansion coefficients in section 2.2.1 and the time-derivative of the orbitals (section 2.2.2) through the *Q*- and *P*-space contributions given in sections 2.2.2.1 and 2.2.2.2, respectively.

## 2.2. Derivation of the working equations

The EOM for the TD-RASSCF-F theory have been already established in [91, 92]. In the following we provide the derivation of the EOM in the case of the TD-RASSCF-B method, and highlight the similarities and differences with respect to the TD-RASSCF-F theory. Starting from the Lagrangian formulation of the time-dependent Schrödinger equation [95], we define the action functional using the TD-RASSCF-B Ansatz, equation (4), as,

$$S[\{C_i(t)\}, \{|\phi_i(t)\rangle\}, \{\epsilon_j^i(t)\}] = \int_{t_1}^{t_2} \left[ \langle\Psi(t)|\hat{K}|\Psi(t)\rangle + \sum_{ij} \epsilon_j^i(t) (\langle\phi_i(t)|\phi_j(t)\rangle - \delta_{ij}) \right] dt, \quad (7)$$

with $\hat{K} \equiv i\partial/\partial t - \hat{H}$ and $\delta_{ij}$ the Kronecker delta function. The Lagrange multipliers, $\epsilon_j^i(t)$, ensure that the orbitals remain orthonormal for all time $t$. In the following the indexes i, j, ... 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 figure 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 as

$$\hat{H}(t) = \sum_{pq} h_q^p(t) \hat{b}_p^\dagger \hat{b}_q + \frac{1}{2} \sum_{pqrs} v_{qs}^{pr}(t) \hat{b}_p^\dagger \hat{b}_r^\dagger \hat{b}_s \hat{b}_q, \quad (8)$$

with $\hat{b}_p(\hat{b}_p^\dagger)$ the annihilation (creation) operator of a particle in the orbital $|\phi_p(t)\rangle$ (see also appendix B). These operators satisfy the commutation relation, $[\hat{b}_p, \hat{b}_q^\dagger] = \hat{b}_p \hat{b}_q^\dagger - \hat{b}_q^\dagger \hat{b}_p = \delta_{qp}$, for bosons and the anti-commutation
relation, {$\hat{b}_p, \hat{b}_q^\dagger$} = $\hat{b}_p \hat{b}_q^\dagger + \hat{b}_q^\dagger \hat{b}_p = \delta_{qp}$, for fermions, see for instance [96]. The matrix elements of the one-body and two-body operators in the basis of the time-dependent orbitals, are expressed as

$$h_q^p(t) = \int \phi_p^*(\mathbf{r}, t) h(\mathbf{r}, t) \phi_q(\mathbf{r}, t) \mathrm{d}\mathbf{r}, \qquad (9)$$

and

$$v_{qs}^{pr}(t) = \iint \phi_p^*(\mathbf{r}, t) \phi_r^*(\mathbf{r}', t) W(\mathbf{r}, \mathbf{r}', t) \phi_q(\mathbf{r}, t) \phi_s(\mathbf{r}', t) \mathrm{d}\mathbf{r}\mathrm{d}\mathbf{r}', \qquad (10)$$

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 [95, 97], the best approximation using the wavefunction Ansatz is obtained by seeking stationarity of the action S, i.e., $\delta S = 0$, for any variation of the parameters and with the boundary condition $|\delta\Psi(t_1)\rangle = |\delta\Psi(t_2)\rangle = 0$. The variation of the action gives,

$$\delta S = \int_{t_1}^{t_2} \left[ \langle \delta\Psi | \hat{K} \Psi \rangle + \langle \hat{K} \Psi | \delta\Psi \rangle + \sum_{ij} [\epsilon_j^i (\langle \delta\phi_i | \phi_j \rangle + \langle \phi_i | \delta\phi_j \rangle) + \delta\epsilon_j^i (\langle \phi_i | \phi_j \rangle - \delta_{ij})] \right] \mathrm{d}t, \qquad (11)$$

where the boundary condition is used to remove the additional term $i\partial_t \langle \Psi | \delta\Psi \rangle$ resulting from the action of $\hat{K}$ on $\langle \Psi |$ instead of $|\delta\Psi\rangle$, see [97]. The variation of the wavefunction is explicitly written as [91],

$$|\delta\Psi\rangle = \sum_{I \in \mathcal{V}} \delta C_I |\Phi_I\rangle + \sum_{pq} \hat{b}_p^\dagger \hat{b}_q |\Psi\rangle \langle \phi_p | \delta\phi_q \rangle. \qquad (12)$$

We can now proceed with the stationarity condition of the action with respect to the parameters {$C_I$}, $\{|\phi_i^i\rangle\}$ and $\{\epsilon_j^i\}$ to obtain the EOM of the TD-RASSCF-B method.

### 2.2.1. 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 $\langle \delta\Psi |$ of the first expectation value in equation (11). Thus, the stationarity condition, $\delta S/\delta C_I^* = 0$, readily leads to $\langle \Phi_I | i\partial_t - \hat{H} | \Psi \rangle = 0$. Moreover, the derivative of the wavefunction with respect to time reads,

$$\frac{\partial}{\partial t} |\Psi\rangle = \sum_{I \in \mathcal{V}} \dot{C}_I |\Phi_I\rangle + \left[ \sum_{pq} \eta_q^p \hat{b}_p^\dagger \hat{b}_q \right] |\Psi\rangle, \qquad (13)$$

with $\eta_q^p \equiv \langle \phi_p | \dot{\phi}_q \rangle$, which results from the time-derivative of the orbitals used to build the configurations in $|\Psi\rangle$. Hereafter, the operator in bracket in equation (13) will be called $\hat{D}$, i.e.,

$$\hat{D} = \sum_{pq} \eta_q^p \hat{b}_p^\dagger \hat{b}_q, \qquad (14)$$

for brevity. We can now rewrite the stationary condition $\delta S/\delta C_I^* = 0$ using the explicit form of $\partial_t |\Psi\rangle$, equation (13), as

$$i\dot{C}_I + \langle \Phi_I | (i\hat{D} - \hat{H}) | \Psi \rangle = 0, \forall I \in \mathcal{V}, \qquad (15)$$

or equivalently, using the expressions of $\hat{H}$ and $\hat{D}$,

$$i\dot{C}_I = \sum_{ij} (h_j^i - i\eta_j^i) \langle \Phi_I | \hat{b}_i^\dagger \hat{b}_j | \Psi \rangle + \frac{1}{2} \sum_{ijkl} v_{jl}^{ik} \langle \Phi_I | \hat{b}_i^\dagger \hat{b}_k^\dagger \hat{b}_l \hat{b}_j | \Psi \rangle. \qquad (16)$$

The indexes in the summations are now restricted to the $\mathcal{P}$-space. It is clear that if either annihilation or creation operators act on an orbital of $\mathcal{Q}$, the inner product with all RAS configurations $\langle \Phi_I |$ vanishes. The EOM for the expansion coefficients, the amplitude equation (16), are identical to those obtained for fermions in the TD-RASSCF-F theory, see [91, 92], and those of the MCTDHB [70] and MCTDHF [98] theories. It is worthwhile to keep in mind that the action of the creation and annihilation operators differs for fermions and bosons. The $\eta_j^i$ matrix elements in the amplitude equations (equation (16)) describe the rotation of the orbitals into one another and are also present in the EOM of the MCTDHB/F methods. In these latter cases, besides to be elements of an anti-Hermitian matrix, there are no constraints on the $\eta_j^i$ 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 $\mathcal{P}_i$-space ($i = 1, 2$). For orbitals which do not belong to the same $\mathcal{P}_i$-space, the $\eta_j^i$ matrix elements must be evaluated, as discussed in [91, 92] and in sections 2.2.2.3 and 2.2.2.4.
### 2.2.2. Equations of motion for the orbitals

Seeking stationarity of S with respect to a variation of an orbital $\langle\phi_i|$, i.e., $\delta S/\delta \langle\phi_i| = 0$, gives

$$ \sum_q |\phi_q\rangle \langle\Psi_i^q| \left[ \sum_{l \in V} i \hat{C}_l |\Phi_l\rangle + (i\hat{D} - \hat{H}) |\Psi\rangle \right] + \sum_j \epsilon_j^i |\phi_j\rangle = 0, \quad (17) $$

with $\langle\Psi_i^q| \equiv \langle\Psi|\hat{b}_i^\dagger \hat{b}_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 figure 1. The EOM for the orbitals of the P- and Q-space are obtained by projecting equation (17) on either an orbital of the P-space, $\langle\phi_i|$, or of the Q-space, $\langle\phi_q|$, as done in the following.

#### 2.2.2.1. Equations of motion for the Q-space orbitals

Starting with the EOM for the Q-space orbitals, we multiply equation (17) from the left with an orbital $\langle\phi_a|$ belonging to the Q-space and obtain,

$$ \sum_{l \in V} i \hat{C}_l \langle \Psi_i^a | \Phi_l \rangle + \langle \Psi_i^a | (i \hat{D} - \hat{H}) | \Psi \rangle = 0, \quad (18) $$

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 $\langle\Psi_i^q|\Phi_l\rangle = \langle\Psi|\hat{b}_i^\dagger \hat{b}_a|\Phi_l\rangle$ vanishes because in all configurations $|\Phi_l\rangle$ the orbital $|\phi_a\rangle$ is unoccupied. Using the explicit expression of the Hamiltonian, equation (8), and for the operator $\hat{D}$, equation (13), we obtain

$$ \sum_{pq} (i\eta_q^p - h_q^p) \langle\Psi|\hat{b}_i^\dagger \hat{b}_a \hat{b}_p^\dagger \hat{b}_q|\Psi\rangle = \frac{1}{2} \sum_{pqrs} v_{qs}^{pr} \langle\Psi|\hat{b}_i^\dagger \hat{b}_a \hat{b}_p^\dagger \hat{b}_r^\dagger \hat{b}_s \hat{b}_q|\Psi\rangle. \quad (19) $$

Using the commutation relation for the creation/annihilation operators for bosons (fermions), we can reestablish the normal ordering of the chains of operators,

$$ \hat{b}_i^\dagger \hat{b}_a \hat{b}_p^\dagger \hat{b}_q = \hat{b}_i^\dagger \hat{b}_q \delta_{pa} + \hat{b}_i^\dagger \hat{b}_p^\dagger \hat{b}_q \hat{b}_a, \quad (20) $$

$$ \hat{b}_i^\dagger \hat{b}_a \hat{b}_p^\dagger \hat{b}_r \hat{b}_s \hat{b}_q = \hat{b}_i^\dagger \hat{b}_r \hat{b}_s \hat{b}_q \delta_{pa} \pm \hat{b}_i^\dagger \hat{b}_p^\dagger \hat{b}_s \hat{b}_q \delta_{ra} + \hat{b}_i^\dagger \hat{b}_p^\dagger \hat{b}_r \hat{b}_s \hat{b}_q \hat{b}_a, \quad (21) $$

with the upper sign holding for bosons and the lower for fermions. The chain of four operators in equation (20) and six operators in equation (21) both annihilate a particle in orbital $|\phi_a\rangle$, from the Q-space, which is not included in $|\Psi\rangle$ and thus vanish. The rhs of equation (19) now reads,

$$ \sum_q (i\eta_q^a - h_q^a) \langle\Psi|\hat{b}_i^\dagger \hat{b}_q|\Psi\rangle = \sum_j (i\eta_j^a - h_j^a) \langle\Psi|\hat{b}_i^\dagger \hat{b}_j|\Psi\rangle, \quad (22) $$

where we restrict the summation over $j \in P$, the summation over the Q-space orbitals being zero. In the same way, inserting equation (21) in the rhs of equation (19) simplifies its expression to

$$ \sum_j (i\eta_j^a - h_j^a) \langle\Psi|\hat{b}_i^\dagger \hat{b}_j|\Psi\rangle = \sum_{jkl} v_{jl}^{ak} \langle\Psi|\hat{b}_i^\dagger \hat{b}_k^\dagger \hat{b}_l \hat{b}_j|\Psi\rangle. \quad (23) $$

Here we used that $v_{jl}^{ak} = v_{lj}^{ka}$ (equation (10)). Interestingly, this equation is exactly the same for bosons and fermions. The time-derivative of the orbitals, included in the term $\eta_j^a$, requires the explicit consideration of the Q-space orbitals. This issue is circumvented by using the projector $\hat{Q}$ onto the subspace spanned by the Q-space orbitals,

$$ \begin{aligned} \hat{Q} &= \sum_a |\phi_a\rangle \langle\phi_a| \\ &= \hat{1} - \sum_i |\phi_i\rangle \langle\phi_i| \\ &= \hat{1} - \hat{P}, \end{aligned} \quad (24) $$

with $\hat{P}$ the projector onto the P-space. Introducing the one-body density matrix, $\rho_i^j = \langle\Psi|\hat{b}_i^\dagger \hat{b}_j|\Psi\rangle$ and the two-body density matrix $\rho_{ik}^j = \langle\Psi|\hat{b}_i^\dagger \hat{b}_k^\dagger \hat{b}_l \hat{b}_j|\Psi\rangle$, we obtain,

$$ i \sum_j \hat{Q} |\dot{\phi}_j\rangle \rho_i^j = \hat{Q} \left[ \sum_j h |\phi_j\rangle \rho_i^j + \sum_{jkl} W_l^k |\phi_j\rangle \rho_{ik}^j \right], \quad (25) $$

with,

$$ W_l^k(\mathbf{r}) = \int \phi_k^*(\mathbf{r}') W(\mathbf{r}, \mathbf{r}') \phi_l(\mathbf{r}') d\mathbf{r}', \quad (26) $$

the mean-field operator, which describes the interaction between the particles and h the one-body Hamiltonian from the rhs of equation (9). The role of the Q-space appears in the time-derivative of the orbitals of the P-space
through the term $\hat{Q}|\dot{\phi}_i\rangle$, see equation (6). We rearrange equation (25), see appendix C, to uncouple the contribution of each $\hat{Q}|\dot{\phi}_i\rangle$, $\forall i \in \mathcal{P}$, and obtain

$$ \hat{Q}|\dot{\phi}_i\rangle = -i(\hat{1} - \hat{P}) \left[ h|\phi_i\rangle + \sum_{jklm} (\underline{\rho}^{-1})^m_i \rho_{mk}^{jl} W_l^k |\phi_j\rangle \right], \quad (27) $$

with $\underline{\rho}^{-1}$ the inverse of the one-body density matrix. The MCTHB theory leads also to equation (27), see [70], but the rhs is subsequently simplified thanks to the choice of the matrix elements $\eta_i^j = 0$ and using $\hat{Q} = \hat{1} - \hat{P}$, see equation (31) below. As discussed in 2.2.1, such a fixed choice of $\eta_i^j$ is not possible in the TD-RASSCF theory. The derivations of the Q-space EOMs differ slightly for bosons and fermions, see equations (20) and (21), but the final result, equation (27), is the same for both types of particles.

### 2.2.2.2. 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, equation (17), we multiply this latter on the left by an orbital of the P-space, $\langle\phi_j|$, leading to,

$$ \sum_{I \in \mathcal{V}} i \dot{C}_I \langle \Psi_i^I | \Phi_I \rangle + \langle \Psi_i^I | (i \hat{D} - \hat{H}) | \Psi \rangle + \epsilon_j^I = 0, \quad (28) $$

with $\hat{D}$ defined in equation (14). This equation still contains the Lagrange multiplier $\epsilon_j^I$. A variation of S with respect to the orbital $|\phi_j\rangle$ and its projection onto the orbital $\langle\phi_j|$, leads to an equation containing the same Lagrange multiplier,

$$ \sum_{I \in \mathcal{V}} (-i \dot{C}_I^* \langle \Phi_I | \Psi_j^I \rangle) + \langle \Psi | (i \hat{D} - \hat{H}) | \Psi_j^I \rangle + \epsilon_j^I = 0, \quad (29) $$

and subtracting equations (28) and (29) gives the EOM for the P-space orbitals, i.e.,

$$ \langle \Psi | (i \hat{D} - \hat{H}) | \Psi_j^I \rangle - \langle \Psi_i^I | (i \hat{D} - \hat{H}) | \Psi \rangle = i \dot{\rho}_I^j, \quad (30) $$

where we have introduced $\dot{\rho}_I^j = \sum_{I \in \mathcal{V}} (\dot{C}_I^* \langle \Phi_I | \Psi_j^I \rangle + \langle \Psi_i^I | \Phi_I \rangle \dot{C}_I)$, with $|\Psi_i^I\rangle \equiv \hat{b}_i^\dagger \hat{b}_i |\Psi\rangle$. The P-space EOM provide the contribution of the P-space in the time-derivative of the orbitals, see equation (6),

$$ \hat{P} |\dot{\phi}_i\rangle = \sum_j |\phi_j\rangle \eta_i^j, \quad (31) $$

through the evaluation of the matrix elements $\eta_i^j$ included in the operator $\hat{D}$, see equation (14). Nonetheless, solving equation (30) is not a trivial task because of the presence of $\dot{\rho}_I^j$, 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 $\eta_i^j$ is still possible for pairs of orbitals which belong to the same $\mathcal{P}_I$-space ($i = 1, 2$), and we use $\eta_i^j = 0, \forall \{i, j\} \in \mathcal{P}_1$ or $\mathcal{P}_2$. The P-space equation (equation (30)) remains to be solved only for pairs of orbitals $\{i', j''\}$, which belong to different $\mathcal{P}_I$-spaces,

$$ \langle \Psi | (i \hat{D} - \hat{H}) | \Psi_j^{j''} \rangle - \langle \Psi_i^{j''} | (i \hat{D} - \hat{H}) | \Psi \rangle = i \rho_i^{j''}, \quad (32) $$

but remains coupled to the amplitude equations through $\dot{\rho}_I^{j''}$. In the meantime it is noted that if equation (32) is solved for $\eta_i^j$, the rhs of equation (31) can be constructed. Moreover equation (27) can be solved, and hence $|\dot{\phi}_i(t)\rangle$ of equation (6) can be evaluated. In the derivation of the TD-RASSCF-F method [91, 92], two different approaches of RAS schemes were proposed to circumvent the difficulty of solving equation (30). These approaches are also used for bosons in the following sections (2.2.2.3) and (2.2.2.4).

### 2.2.2.3. Even excitation RAS scheme

First we suggest to consider the case in which only an even number of particles is promoted from $\mathcal{P}_1$ to $\mathcal{P}_2$, see figure 2(a). In this case, $\dot{\rho}_I^{j''}$ explicitly reads,

$$ \dot{\rho}_I^{j''} = \sum_{I \in \mathcal{V}} (\dot{C}_I^* \langle \Phi_I | \hat{b}_{i'}^\dagger \hat{b}_{j''} | \Psi \rangle + \langle \Psi | \hat{b}_{i'}^\dagger \hat{b}_{j''} | \Phi_I \rangle \dot{C}_I). \quad (33) $$

The action of $\hat{b}_{i'}^\dagger \hat{b}_{j''}$ on the wavefunction $|\Psi\rangle$ annihilates one particle in $\mathcal{P}_2$ and creates one in $\mathcal{P}_1$. Since only an even number of particles is present in $\mathcal{P}_2$, $\hat{b}_{i'}^\dagger \hat{b}_{j''}|\Psi\rangle$ would contain only configurations with an odd number of particles in $\mathcal{P}_2$, which makes the inner product with $\langle\Phi_I|, \forall I \in \mathcal{V}$ vanish. In the same manner $\hat{b}_{i'}^\dagger \hat{b}_{j''}$ acting on $|\Phi_I\rangle$ is either zero, if $|\phi_{j''}\rangle$ is unoccupied in the configuration $|\Phi_I\rangle$, or gives an odd number of particles in $\mathcal{P}_2$. In this specific excitation scheme, $\dot{\rho}_I^{j''} = 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 (equation (8)) and the operator $\hat{D}$ (equation (14)), equation (32) reads,
<figure><img src="image_5.png" /><figcaption>Figure 2. Illustration of the decomposition of the Fock-space used in the TD-RASSCF-B method with N = 8 bosons, a single $\mathcal{P}_1$ orbital, $M_1 = 1$ (thick (green) line) and $M_2 = 5$ $\mathcal{P}_2$ orbitals (thin (red) lines). (a) For a RAS scheme including only even excitations, the RAS Fock-space $\mathcal{V}$ is decomposed into the direct sum of subspaces $V_n$ with $n = 0, 2, 4, \dots, N_{\max}$. Each subspace $V_n$ includes all the configurations with $n$ bosons in $\mathcal{P}_2$ and $N-n$ bosons in $\mathcal{P}_1$. (b) The decomposition of the RAS Fock-space for the general excitation scheme is also written as a direct sum of subspaces $V_n$ except that in that case $n = 0, 1, 2, \dots, N_{\max}$.</figcaption></figure>

$$ \sum_{pq} (h_q^P - i\eta_q^P) [\langle\Psi|\hat{b}_{i'}^\dagger \hat{b}_{j''}^\dagger \hat{b}_p^\dagger \hat{b}_q - \hat{b}_p^\dagger \hat{b}_q \hat{b}_{i'}^\dagger \hat{b}_{j''}|\Psi\rangle + \frac{1}{2} \sum_{pqrs} \nu_{qs}^{pr} [\langle\Psi|\hat{b}_{i'}^\dagger \hat{b}_{j''}^\dagger \hat{b}_p^\dagger \hat{b}_r^\dagger \hat{b}_s \hat{b}_q - \hat{b}_p^\dagger \hat{b}_r^\dagger \hat{b}_s \hat{b}_q \hat{b}_{i'}^\dagger \hat{b}_{j''}^\dagger |\Psi\rangle] = 0. \quad (34) $$

We can simplify this expression, starting with

$$ \begin{aligned} & \langle\Psi|\hat{b}_{i'}^\dagger \hat{b}_{j''}^\dagger \hat{b}_p^\dagger \hat{b}_q - \hat{b}_p^\dagger \hat{b}_q \hat{b}_{i'}^\dagger \hat{b}_{j''}^\dagger |\Psi\rangle \\ &= \langle\Psi|\hat{b}_{i'}^\dagger \hat{b}_q \delta_{pj''} \pm \hat{b}_{i'}^\dagger \hat{b}_p^\dagger \hat{b}_{j''}^\dagger \hat{b}_q - \hat{b}_p^\dagger \hat{b}_{j''}^\dagger \delta_{i'q} \mp \hat{b}_{i'}^\dagger \hat{b}_p^\dagger \hat{b}_{j''}^\dagger \hat{b}_q |\Psi\rangle \\ &= \langle\Psi|\hat{b}_{i'}^\dagger \hat{b}_q \delta_{pj''} - \hat{b}_p^\dagger \hat{b}_{j''}^\dagger \delta_{i'q} |\Psi\rangle \\ &= \rho_i^q \delta_{pj''} - \rho_p^{j''} \delta_{i'q} \\ &= A_{pi'}^{qj''}. \end{aligned} \quad (35) $$

Now we turn to the chains of six operators in the last term in equation (34). The first product of operators is expressed as

$$ \hat{b}_{i'}^\dagger \hat{b}_{j''}^\dagger \hat{b}_p^\dagger \hat{b}_s \hat{b}_q = \hat{b}_{i'}^\dagger \hat{b}_r^\dagger \hat{b}_s \hat{b}_q \delta_{pj''} \pm \hat{b}_{i'}^\dagger \hat{b}_p^\dagger \hat{b}_s \hat{b}_q \delta_{ri''} + \hat{b}_{i'}^\dagger \hat{b}_p^\dagger \hat{b}_r^\dagger \hat{b}_{j''}^\dagger \hat{b}_s \hat{b}_q \quad (36) $$

and the second product of operators as

$$ -\hat{b}_p^\dagger \hat{b}_r^\dagger \hat{b}_s \hat{b}_q \hat{b}_{i'}^\dagger \hat{b}_{j''}^\dagger = -\hat{b}_p^\dagger \hat{b}_r^\dagger \hat{b}_s \hat{b}_{j''}^\dagger \delta_{i'q} \mp \hat{b}_p^\dagger \hat{b}_r^\dagger \hat{b}_q \hat{b}_{j''}^\dagger \delta_{i's} - \hat{b}_{i'}^\dagger \hat{b}_p^\dagger \hat{b}_r^\dagger \hat{b}_{j''}^\dagger \hat{b}_s \hat{b}_q. \quad (37) $$

The sum of equations (36) and (37) enters equation (34), and we see that the chains of six operators cancel each other and only the chains of four operators remain. Using the fact that $v_{ji}^{ik} = v_{ij}^{ki}$ (equation (10)) and that $\eta_q^P$ must be evaluated for orbitals which belong to different $\mathcal{P}_i$ space, $\{l', k''\}$, equation (34) can be rewritten,

$$ \sum_{k''l'} (h_{l'}^{k''} - i\eta_{l'}^{k''}) A_{k''l'}^{l''j''} + \sum_{klm} (v_{kl}^{j''m} \rho_{i'm}^{kl} - v_{i'm}^{kl} \rho_{kl}^{j''m}) = 0. \quad (38) $$

This equation, used to determine $\eta_{l'}^{k''}$, is identical for fermions and bosons, only the evaluation of the one- and two-body reduced density matrices depends on the kind of particles. The coefficients and the $\mathcal{P}$-space orbitals equations are separable and can now be solved (see appendix C). The $\eta_{l'}^{k''}$, obtained from equation (38), are used to determine the time-derivative of the coefficients from equation (16). The time-derivative of the $\mathcal{P}$-space orbitals, equation (6), is obtained from equation (31) in addition to the $\mathcal{Q}$-space equations (equation (27)). Interestingly, for the case of only even excitations where $\dot{\rho}_i^j = 0$, equation (30) is the analog of the generalized
Brillouin's theorem, or Levy-Berthier-Brillouin theorem [99, 100], for the time-dependent Schrödinger equation with $\hat{H} = i\partial/\partial t - \hat{H}$. It means that the single excitations of the wavefunction, i.e., $\hat{b}_i^\dagger \hat{b}_j |\Psi\rangle \equiv |\Psi_i^j\rangle$, are not coupled to the wavefunction $|\Psi\rangle$ through the action of $\hat{H}$, $\langle\Psi|\hat{H}|\Psi_i^j\rangle = 0$. Thus, the TD-RASSCF-B wavefunction with only even excitations always includes the subspace of single excitations, see also [59, 91]. In other words, projecting the TD-RASSCF-B wavefunction onto a time-independent basis would include non-vanishing projection in the singles space spanned by the time-independent basis functions.

### 2.2.2.4. 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 $\mathcal{P}_2$ from 0 to $N_{max}$, where $N_{max}$, defined in section 2.1, is the highest number of particles allowed in $\mathcal{P}_2$, see figure 2(b). Note that $N_{max}$ must fulfill the condition $N_{max} \le N$. For instance, taking $N_{max} = 4$, we consider the promotion of 0, 1, 2, 3 and 4 particles from $\mathcal{P}_1$ to $\mathcal{P}_2$. In this way, the configurational space is the direct sum of $N_{max} + 1$ subspaces,

$$ \mathcal{V} = \mathcal{V}_0 \oplus \mathcal{V}_1 \oplus \cdots \oplus \mathcal{V}_{N_{max}}. \qquad (39) $$

Using the expression of the time derivative of the $\{C_i\}$ coefficients, equation (15), the time derivative of the one-body density matrix, in equation (30), can be expressed as,

$$ i\rho_i^{j''} = \sum_{I \in \mathcal{V}} (\langle \Psi | (i\hat{D} - \hat{H}) | \Phi_I \rangle \langle \Phi_I | \Psi_{j''}^I \rangle - \langle \Psi_{j''}^I | \Phi_I \rangle \langle \Phi_I | (i\hat{D} - \hat{H}) | \Psi \rangle). \qquad (40) $$

We introduce the projector onto the RAS space $\mathcal{V}$ as $\hat{\Pi} = \sum_{I \in \mathcal{V}} |\Phi_I\rangle\langle\Phi_I|$. Using the above expression of $\rho_i^j$, we obtain a new formulation of the $\mathcal{P}$-space orbital equation,

$$ \langle \Psi | (i\hat{D} - \hat{H}) (\hat{1} - \hat{\Pi}) | \Psi_{j''}^I \rangle - \langle \Psi_{j''}^I | (\hat{1} - \hat{\Pi}) (i\hat{D} - \hat{H}) | \Psi \rangle = 0. \qquad (41) $$

For $|\phi_i^j\rangle \in P_1$ and $|\phi_{j''}^k\rangle \in P_2$, we note that $|\Psi_{j''}^I\rangle$ belongs to $\mathcal{V}$, with one particle from $\mathcal{P}_2$ being annihilated and one particle in $\mathcal{P}_1$ created, leading to $(\hat{1} - \hat{\Pi})|\Psi_{j''}^I\rangle = 0$. On the other hand, $\langle \Psi_{j''}^I | = \langle \Psi | \hat{b}_{j''}^\dagger \hat{b}_{j''}^I \rangle$, provides configurations with a creation of an additional particle in $\mathcal{P}_2$, which may lie in $\mathcal{V}_{N_{max}+1}$, not included in $\mathcal{V}$. In this case, $\langle \Psi_{j''}^I | (\hat{1} - \hat{\Pi}) \ne 0$ and equation (41) simplifies to,

$$ \langle \Psi_{j''}^I | (\hat{1} - \hat{\Pi}) (i\hat{D} - \hat{H}) | \Psi \rangle = 0, \quad \text{with} \quad \langle \Psi_{j''}^I | (\hat{1} - \hat{\Pi}) = \sum_{I \in \mathcal{V}_{N_{max}}} C_I^* \langle \Phi_I | \hat{b}_{j''}^I \hat{b}_{j''}^I. \qquad (42) $$

Using the expression of the Hamiltonian (equation (8)) and of the operator $\hat{D}$ (equation (14)) while keeping in mind that $\eta_{l'}^k$ has only to be determined for pairs of orbitals $\{l', k''\}$ which belong to different $\mathcal{P}_i$ subspaces, equation (42) is equivalent to

$$ \sum_{k''l'} (i\eta_{l'}^{k''} - h_{l'}^{k''}) \zeta_{k''l'}^{l'j''} = \frac{1}{2} \sum_{klmn} v_{lm}^{km} \zeta_{kmi'}^{lnj''}, \qquad (43) $$

where the fourth- and sixth-order tensors are defined by

$$ \zeta_{k'mi'}^{l'j''} = \langle \Psi_{i'}^j | (\hat{1} - \hat{\Pi}) \hat{b}_{k'm}^\dagger \hat{b}_{l'} | \Psi \rangle, \qquad (44) $$

$$ \zeta_{kmi'}^{lnj''} = \langle \Psi_{i'}^j | (\hat{1} - \hat{\Pi}) \hat{b}_k^\dagger \hat{b}_m^\dagger \hat{b}_n \hat{b}_l | \Psi \rangle. \qquad (45) $$

Here again, the $\mathcal{P}$-space EOM for the determination of the $\eta_{l'}^k$, equation (43), are identical for bosons and fermions [91, 92]. These equations are solved to determine the $\eta_{l'}^k$ for each pairs of orbitals belonging to different $\mathcal{P}_i$-spaces. The value of $\eta_{l'}^k$ is subsequently used to solve the amplitudes equations (equation (16)) and to evaluate the time-derivative of the $\mathcal{P}$-space orbitals from equations (31) and (27), 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 [91, 92]. 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, equation (16), the $\mathcal{Q}$-space equations, equation (27) and the $\mathcal{P}$-space equations equation (43) and equation (38) for the general RAS scheme and the only even excitation scheme, respectively. In the case of the MCTDHB method, the amplitude (equation (16)) and the $\mathcal{Q}$-space equations (equation (27)) are also solved to obtain the time-derivative 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 than the MCTDHB
method to solve equations (16) and (27), 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 (equation (38) or (43)). 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, equation (45), 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 [101] where a unified set of EOM for the MCTDHB theory for both bosons and fermions was derived.

### 3. Application to a time-independent system: ground state energy

In this section, we consider a system of $N = 100$ bosons trapped in a one-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 [102–106]. In the following, we consider an anisotropic harmonic trap such that the longitudinal frequency ($\omega_x$) is much smaller than the transversal frequency ($\omega_\perp$), i.e., $\omega_\perp \gg \omega_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,

$$ \hat{H} = \frac{1}{2} \sum_{i=1}^{N} \left( -\frac{\partial^2}{\partial x_i^2} + x_i^2 \right) + \lambda \sum_{i<j} \delta(x_i - x_j), \quad (46) $$

using the unit of length $l_0 = \sqrt{\hbar (m\omega_x)^{-1}}$ and the unit of energy $E_0 = \hbar\omega_x$, with m the mass of the particles. Assuming no confinement induced resonances [107], the interaction strength, $\lambda$, is related to the 3D s-wave scattering length of the particles, $a_s$, through $\lambda = 2a_s l_0^{-2}$, with $l_\perp$ 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 [10, 11]. The different parameters used to perform the numerical calculations are reported in appendix E.

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 [108]), the lower the energy, the higher the accuracy. First, as a general remark, for any value of $\lambda \neq 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 1 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,

$$ \mathcal{E}_{\text{corr}} = E_{\text{GP}} - E_{\text{method}}, \quad (47) $$

where $E_{\text{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, $\mathcal{E}_{\text{ref}}$, the correlation energy obtained for the MCTDHB method with 5 orbitals, i.e.,

$$ \mathcal{E}_{\text{ref}} = E_{\text{GP}} - E_{\text{MCTDHB}}^{5}, \quad (48) $$

where the superscript 5 denotes 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 $\mathcal{E}_{\text{ref}}$. Secondly, we define the relative correlation energy, $\mathcal{E}_{\text{rel}}^X$, as the difference between the energy obtained from a MCTDHB calculation with X orbitals and the GP energy, i.e.,

$$ \mathcal{E}_{\text{rel}}^{X} = E_{\text{GP}} - E_{\text{MCTDHB}}^{X}. \quad (49) $$

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
<table><caption>**Table 1.** Ground-state energy (in units of $E_0$, see text after equation (46)) of 100 bosons trapped in a 1D harmonic potential interacting through a contact potential with a strength $\lambda = 0.01$. The TD-RASSCF-B calculations were performed with a single $P_1$ orbital, $M_1 = 1$ and $M_2 = M - 1$ $P_2$ orbitals, with $M$ the total number of orbitals. The results were obtained using either the general RAS scheme, labeled 'All excitations', or only even excitations, labeled 'Even excitations'. The excitation schemes are indicated with the usual notations -S, -SD, ..., up to $N_{max} = 9$ and the RAS schemes are labeled by the value of $N_{max}$ for larger excitations, e.g. -10. In addition MCTDHB calculations were carried out to compare the accuracy and efficiency of the TD-RASSCF-B method. The number of configurations used in the wavefunction expansion are indicated in parentheses under each energy. The result obtained with a single orbital is equivalent to the GF method. To highlight the difference between the TD-RASSCF-B and MCTDHB results, the digits that differ are underlined.</caption><thead><tr><th rowspan="2">Method</th><th colspan="10">Orbitals</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th></tr></thead><tbody><tr><td><u>All excitations</u><br>MCTDHB</td><td>68.76816487<br>(1)</td><td>68.75335446<br>(101)</td><td>68.74538390<br>(5151)</td><td>68.74152088<br>(176851)</td><td>68.73891122<br>(4598126)</td><td>—<br>(96560646)</td><td>—<br>(1705904746)</td><td>—<br>(26075972546)</td></tr><tr><td>-SD</td><td>—</td><td>68.75355024<br>(3)</td><td>68.74565678<br>(6)</td><td>68.74184781<br>(10)</td><td>68.73926413<br>(15)</td><td>68.73761231<br>(21)</td><td>68.736360917<br>(28)</td><td>68.73545355<br>(36)</td></tr><tr><td>-SDTQ</td><td>—</td><td>68.75335660<br>(5)</td><td>68.74538672<br>(15)</td><td>68.74152449<br>(35)</td><td>68.73891508<br>(70)</td><td>68.73724366<br>(126)</td><td>68.73598073<br>(210)</td><td>68.73506372<br>(330)</td></tr><tr><td>-SDTQ56</td><td>—</td><td>68.75335448<br>(7)</td><td>68.74538393<br>(28)</td><td>68.74152092<br>(84)</td><td>68.73891126<br>(210)</td><td>68.73723959<br>(462)</td><td>68.73597655<br>(924)</td><td>68.73505943<br>(1716)</td></tr><tr><td>-SDTQ5678</td><td>—</td><td>68.75335446<br>(9)</td><td>68.74538390<br>(45)</td><td>68.74152088<br>(165)</td><td>68.73891122<br>(495)</td><td>68.73723955<br>(1287)</td><td>68.73597651<br>(3003)</td><td>68.73505938<br>(6435)</td></tr><tr><td>-10</td><td>—</td><td>68.75335446<br>(11)</td><td>68.74538390<br>(66)</td><td>68.74152088<br>(286)</td><td>68.73891122<br>(1001)</td><td>68.73723955<br>(3003)</td><td>68.73597651<br>(8008)</td><td>68.73505938<br>(19448)</td></tr><tr><td><u>Even excitations</u></td><td>—</td><td>68.75355024<br>(2)</td><td>68.74565780<br>(4)</td><td>68.74184905<br>(7)</td><td>68.7392653<br>(11)</td><td>68.73761353<br>(16)</td><td>68.73636218<br>(22)</td><td>68.73545481<br>(29)</td></tr><tr><td>-DQ</td><td>—</td><td>68.75335660<br>(3)</td><td>68.74538974<br>(9)</td><td>68.74153257<br>(22)</td><td>68.73892541<br>(46)</td><td>68.73725598<br>(86)</td><td>68.73599408<br>(148)</td><td>68.73507803<br>(239)</td></tr><tr><td>-DQ6</td><td>—</td><td>68.75335448<br>(4)</td><td>68.74538700<br>(16)</td><td>68.74152917<br>(50)</td><td>68.73892180<br>(130)</td><td>68.73725217<br>(296)</td><td>68.73599018<br>(610)</td><td>68.73507404<br>(1163)</td></tr><tr><td>-DQ68</td><td>—</td><td>68.75335446<br>(5)</td><td>68.74538697<br>(25)</td><td>68.74152914<br>(95)</td><td>68.73892177<br>(295)</td><td>68.73725213<br>(791)</td><td>68.73599014<br>(1897)</td><td>68.73507400<br>(4166)</td></tr><tr><td>-10</td><td>—</td><td>68.75335446<br>(6)</td><td>68.74538697<br>(36)</td><td>68.74152914<br>(161)</td><td>68.73892177<br>(581)</td><td>68.73725213<br>(792)</td><td>68.73599014<br>(4900)</td><td>68.73507400<br>(12174)</td></tr></tbody></table>
MCTDHB wavefunction. Thus, when the correlation energy of a RAS scheme with $X$ orbitals is equal to $\mathcal{E}_{\text{rel}}^X$ the calculation is converged.

In table 1, we report the results for the GS energy (in units of $E_0$) obtained with $\lambda = 0.01$. The reference for the correlation energy is $\mathcal{E}_{\text{ref}} = -2.9 \times 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 $\mathcal{E}_{\text{ref}}$. Specifically we obtain 51%, 78% and 91% of $\mathcal{E}_{\text{ref}}$ for 2, 3, and 4 orbitals, respectively. The $\simeq 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 1. Unfortunately, the MCTDHB wavefunction with 5 orbitals includes already $\sim 4.6 \times 10^6$ configurations and using 6 (8) orbitals leads to $\sim 96 \times 10^6$ ($\sim 26 \times 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 1 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_{\text{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 1. Thus, $\mathcal{E}_{\text{ref}}$ is recovered for the -SDTQ5678 scheme with 5 orbitals, but the expansion of the wavefunction includes only 495 configurations, i.e., $9 \times 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 1, consists of 2 orbitals and the -SD RAS scheme. The correlation energy includes 98.7% of $\mathcal{E}_{\text{rel}}^2$ 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 $\mathcal{E}_{\text{rel}}^{X=3,4,5}$ correlation energies. Thus, using the -SD scheme with 5 orbitals 98.8% of $\mathcal{E}_{\text{ref}}$ is obtained but the TD-RASSCF-B wavefunction includes only 15 configurations while the MCTDHB wavefunction includes more than $4.5 \times 10^6$ configurations, i.e., $\sim 3 \times 10^5$ times more configurations. Moreover, the energy difference between the -SD scheme and the MCDTHB method is below $5 \times 10^{-4}$, i.e., systematically lower than the convergence obtained with respect to the number of orbitals. Concerning the correlation energy of the -SDTQ and -SDTQ568 schemes, we find that they include 99.99% and 99.9999% of $\mathcal{E}_{\text{rel}}^X$, 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 $\mathcal{E}_{\text{rel}}^X$, 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 1, 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 $\sim 4.6 \times 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, whose practical applicability is limited by the exponential growth of the number of configurations.

We also consider RAS schemes with only even excitations (see table 1), 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 \times 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 comparison of the -DQ and -SDTQ schemes with an energy difference below $1.5 \times 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 since the $P$-space EOM, equation (38), does not require the update of a sixth-order tensor at each time-step as it is the case for the general RAS schemes, see equation (45). For values of $N_{\text{max}} \ge 6$ and $M \ge 3$, the energy converges with respect to $N_{\text{max}}$, as the energy does not change by increasing $N_{\text{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 $\sim 1.1 \times 10^{-5}$ and includes 99.96% of $\mathcal{E}_{\text{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 \times 10^{-5}$ for $M > 5$, which is two orders of magnitude smaller than the convergence obtained with respect to the number orbitals, $\sim 10^{-3}$.

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 analysis of the natural occupations of the natural orbitals [109]. 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 (lower energies) were systematically obtained using the TD-RASSCF-B method. This result can be understood as follows. The MCTDHB wavefunction for a small number of orbitals generates a large number of configurations, as all orbitals can be equally populated. Most of these configurations, however, do not contribute in lowering 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 $\lambda$. As a limiting case, we note that in the Tonks-Girardeau model [20], obtained for an infinite value of $\lambda$, each boson occupies a different orbital. 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 a large number of orbitals, which is beyond the possibility of the MCTDHB method. The above conclusions are supported by an analysis of the ground state properties for an interaction strength 10 times larger ($\lambda = 0.1$) presented in appendix E.

## 4. 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 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 $\hat{H}_0$ reads,

$$ \hat{H}_0 = \frac{1}{2} \sum_{i=1}^{N} \left( -\frac{\partial^2}{\partial x_i^2} + x_i^2 \right), \qquad (50) $$

where we use the units described in section 3 and the time is expressed in units of $t_0 = \omega_x^{-1}$. 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 [25, 26], leading to the evolution of the system under the action of the following Hamiltonian,

$$ \hat{H} = \frac{1}{2} \sum_{i=1}^{N} \left( -\frac{\partial^2}{\partial x_i^2} + x_i^2 \right) + \sum_{i<j} \lambda (x_i - x_j)^2, \qquad (51) $$

with $\lambda$ 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 $\Omega_n = 2n\sqrt{\omega_x^2 + 2N\lambda}$, with $\omega_x$ the frequency of the harmonic trap, see [25]. For positive values of $\lambda$ the two-body interaction is attractive, while for negative values the interaction is repulsive and leads to unbound dynamics for $\lambda < -\omega_x/2N$. The different parameters used to perform the numerical calculations are reported in appendix E.

### 4.1. Breathing dynamics with $\lambda = 0.1$

We first consider the dynamics following a quench of the interaction strength from $\lambda = 0$ to $\lambda = 0.1$ (equation (51)). 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 \le t \le 15$, see figure 3. The time evolution of the system is characterized by the one-particle density, $\rho(x = 0, t)$, at the center of the trap $x = 0$ and exhibits a periodic evolution with a frequency $\omega_{MCTDHB} = 3.46$. This value agrees perfectly with the analytical prediction $\Omega_{n=1} = 3.46$ meaning that the first excited state is mainly responsible for the dynamics. Nonetheless, the discrepancy with a pure cosine function $\sim \cos(\Omega_1 t)$ indicates the role of higher excited states with higher harmonic frequencies [25]. The mean-field GP fails to describe the system evolution, even at short time ($t < 1$) as depicted in figure 3(a) and we obtain a lower frequency $\omega_{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 figure 3(a). For a short time, i.e., $0 \le t \le 5$, all RAS schemes describe accurately the dynamics of the system, in
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Time evolution of the one-body density at the center of the harmonic trap, $\rho(x=0, t)$ after an instantaneous quench of the interparticle interaction strength from $\lambda=0$ to $\lambda=0.1$ (equations (50) and (51)). The result obtained using the mean-field GP theory (thin (black) line) strongly disagrees with the numerically exact result (open (red) circles) obtained using MCTDHB with 4 orbitals in all panels. Note that the results are plotted for $0 \le t \le 5$ on the left and $10 \le t \le 15$ on the right. (a) Results of the TD-RASSCF-B method using a single $P_1$ orbital, $M_1=1$, and $M_2=3$ $P_2$ orbitals for different RAS schemes, namely -SD (dashed (green)), -SDTQ (dash-dot (purple)) and -SDTQ56 (thick (blue) line). For short time, the lines are on top of the MCTDHB result. (b) Results of the TD-RASSCF-B method using $M_1=2$ and $M_2=2$ keeping the total number of orbitals constant for the -SD (dashed (green)) and -SDT (thick (blue) line) schemes. The results are on top of the MCTDHB curve, which show the good accuracy of the method. (c) Results using a single $P_1$ orbital, $M_1=1$, and $M_2=3$ $P_2$ orbitals for only even excitations, i.e., -D (dashed (green)) and -DQ (dash-dot (purple)) RAS schemes. (d) Results of the TD-RASSCF-B method using $M_1=2$ and $M_2=2$ keeping the total number of orbitals constant.</figcaption></figure>

contrast to the GP result. On the scale of the figure, the convergence to the MCTDHB result is obtained by using the -SDTQ excitation scheme including 35 configurations, a reduction by a factor of ~8. For a longer time, $10 \le t \le 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 figure 3(b), and provides a better description of long time dynamics than the -SDTQ scheme used previously (figure 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 figure 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 applied the -D and -DQ schemes with $M_1 = 1$ and $M_1 = 2$ figures 3(c) and (d), respectively, and keep $M = 4$. For $M_1 = 1$, figure 3(c), 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 figure 3(d). 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 fails 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.

## 4.2. Breathing dynamics with $\lambda = 0.5$

We continue the illustration of the TD-RASSCF-B method by considering a quenching from $\lambda = 0$ to $\lambda = 0.5$ (equation (51)). This interaction strength was used in [25] to benchmark the MCTDHB method. For the time interval considered here, $0 \le t \le 15$, we find that the result obtained with the MCTDHB method using 8 orbitals and 19448 configurations is numerically exact, in agreement with [25]. As previously, the one-body density exhibits oscillations as a function of the time with a period $\omega_{MCTDHB} = 6.63$, in perfect agreement with the analytical frequency $\Omega_{n=1} = 6.63$. In comparison to the previous results in section 4.1, the shape of the oscillations indicates that the role of higher excited states is stronger as a large deviation from a simple cosine function, $\sim \cos(\Omega_1 t)$, is obtained. This is not surprising since the value of $\lambda$ 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 figure 4(a), which strongly deviates from the MCTDHB result for $t > 0.3$ with a larger amplitude in the oscillations and gives a lower frequency for the oscillations, $\omega_{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, figure 4(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 MCTDHB. 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_1 + M_2 = M = 8$. In figures 4(b)-(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 \le t \le 15$. For lower excitation schemes, the results are accurate for $0 \le t \le 3$ and using the -SDTQ scheme, the frequency of the oscillation at a longer time are correctly obtained, see figure 4(b). For $M_1 = 3$, the results in figure 4(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 albeit with a smaller amplitude. Finally, figure 4(d) reports converged results for $M_1 = 4$ using the -SDT schemes, which includes 5216 configurations. Already, the -SD scheme with 2816 configurations is quite 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 changes 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
<figure><img src="image_7.png" /></figure>

**Figure 4.** Time evolution of the one-body density at the center of the harmonic trap, $\rho(x=0, t)$ after an instantaneous quench of the interparticle interaction strength from $\lambda=0$ to $\lambda=0.5$ (equations (50) and (51)). The result obtained using the mean-field GP theory (thin black line) strongly disagrees with the numerically exact result (open red circles) obtained using MCTDHB with 8 orbitals in each panel. Note that the results are plotted for $0 \le t \le 3$ on the left and $12 \le t \le 15$ on the right. (a) Results of the TD-RASSCF-B method using a single $P_1$ orbital, $M_1=1$, and $M_2=7$ $P_2$ orbitals for different RAS schemes, namely -SD (dashed (green)), -SDTQ (dash-dot (purple)), -SDTQ56 (thick (blue) line) and -SDTQ5678 (dash-dot (black)). (b) Results using $M_1=2$ and $M_2=6$ keeping constant the total number of orbitals. For short time, $0 \le t \le 3$, all RAS schemes are almost completely on top of the MCTDHB result. (c) TD-RASSCF-B simulations using $M_1=3$ and $M_2=5$ orbitals, using the -SDT (solid (purple)) and -SDTQ5 (dash-dot (black)) RAS schemes. A very good agreement is obtained in comparison to the MCTDHB result. (d) Results of the TD-RASSCF-B method using $M_1=4$ and $M_2=4$ orbitals using the -SD and -SDT schemes. This latter converges to the MCTDHB results but includes ~4 times fewer configurations than the MCTDHB wavefunction.

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. Due to the drastic reduction in the size of configuration space, the CPU time needed to perform the TD-RASSCF-B calculations is 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 $\lambda=0.1$, the MCTDHB theory converged to the numerically exact result for 4 orbitals. The TD-RASSCF-B method using different sizes of the $P_1$-space and different RAS schemes can accurately describe the dynamics of the system characterized by $\rho(x=0, t)$ with ~4 times fewer configurations. Moreover, converged TD-RASSCF-B results were obtained with substantially less configurations, for instance with $M_1=2$ and the -SDT excitation scheme. In the case of
$\lambda = 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_f = 4$ and the -SDT RAS scheme. Moreover, all calculations performed with the TD-RASSCF-B method were better than the mean-field GP theory, which failed to describe both scenarios. Physically, the sudden quench of the interaction between the particles projects the uncorrelated GS of the Hamiltonian of equation (50) onto the correlated many-body eigenstates of the Hamiltonian of equation (51). The time evolution of these latter eigenstates is trivial. They accumulate phase according to $\exp(-iE_j t)$, with $E_j$ the energy of the jth eigenstate, and the dynamics results from the coherent superposition of these eigenstates. Thus, the discrepancies that can be observed in figures 3 and 4 between the exact MCTDHB results and the different TD-RASSCF-B results assures the error of the different methods to correctly describe the relevant excited states for the dynamics. Moreover, the fact that the RAS schemes with only even excitations (figure 3(c)) provide an accurate description of the dynamics highlights the role of the correlation as the configurations are only coupled through the two-body operator, see equations (8) and (10).

## 5. Conclusion and outlook

In this work, we presented a general theory for the TD-RASSCF method for indistinguishable particles, which includes the derivation obtained previously for fermions (TD-RASSCF-F) [91, 92, 94] 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. The TD-RASSCF-B provides a better accuracy by including more orbitals while constraining the number of configurations and reducing considerably, in that way, the size of the configurational space. In the ground-state calculations, the flexibility of the RAS scheme allowed us in this sense to identify clearly the physics that needs to be accurately represented: (i) Even (two-body) excitations are more important than single excitations. (ii) While most particles remain in the motional ground-state (SCF) orbital, it is important to account for the relatively fast motion of a few particles by including many orbitals. 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 is possible with the MCTDHB method. We also provided a comparison between the TD-GP, the MCTDHB and the TD-RASSCF-B methods in the case of breathing dynamics induced by a sudden quenching of the interaction strength with two different initial conditions. As the MCTDHB method, the TD-RASSCF-B method does not have any restriction on the choice of the two-body interaction potential used. This was illustrated by the use of a non-contact harmonic interaction between the bosons. We showed that for the studied interaction quench, often used in cold atom experiments, the TD-GP approach becomes inaccurate within a fraction of an oscillation period of the harmonic trapping potential. 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 with substantially less configurations and thus less CPU time than with the MCTDHB method. Finally, as was the case for the ground-state calculations, the RAS scheme allowed us to identify the nature of the excitations that are important for this kind of dynamics. In this case accurate results can be obtained using only even SCF excitations, highlighting the importance of the two-body interaction. In the future the accuracy and numerical efficiency of the TD-RASSCF-B method can be exploited further to solve the TDSE beyond the mean-field approach, and for situations out of reach by other time-dependent wave function based methods.
Dynamical effects such as the four wave mixing process [110] used to produce correlated atom beams [111] or the dynamics of bright [112, 113] and dark [114, 115] 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 [116, 117]. In all of this, the RAS concept will be highly operational in identifying the most important excitations, and hence the physical nature of a given process, as was the case in the illustrative applications of the present work.
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, QU SCOPE.

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. Equations (15), (25) and (38) or (43), only depend on the type of particles through the creation and annihilation operators and the set of configurations {|Φᵢ⟩}. Moreover, additional illustrative examples are presented and discussed.

Appendix A. Compact representation of the wavefunction

Our implementation is based on the general mapping of bosonic operators in Fock space introduced in [118] and implemented, for instance, for bosons [119] and fermions [27] 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₁, n₂, ..., nₘ⟩, represents a configuration with n₁ bosons in the orbital |φ₁⟩, n₂ bosons in the orbital |φ₂⟩, etc. Such a configuration is indexed by a unique integer J defined as,

$$J = 1 + \sum_{k=1}^{M} \left( N + M - 1 - k - \sum_{l=1}^{k} n_l \right) \quad (A1)$$

Thus, for each configuration we store its complex coefficient Cⱼ in an array according to the index J provided by the above mapping. In [119], this mapping was employed to avoid the storage of the configuration vectors |n₁, n₂, ..., nₘ⟩, which can be prohibitively memory consuming in the case of the MCTDHB method. Thus, to access the coefficient of the configurations, a set of M-nested loops over the occupation number nᵢ is used to span the full configurational space and to compute, using equation (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 [119] 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₁ and M₂ the number of orbitals in the P₁ and P₂ spaces, respectively. For V₀, i.e., configurations with particles only in P₁, see equation (39), we enumerate all possible configurations, evaluate their index, using equation (A1), and store them. Then for the excited configurations, i.e., configurations with one or more particles in P₂, we introduce the excitation number nₐₐₓ which is equivalent to the number of particles in P₂. The number of remaining particles in P₁ is N - nₐₐₓ. For each excitation we enumerate the configurations of N - nₐₐₓ particles in the P₁ orbitals and compute their index. The same is performed for the configurations with nₐₐₓ in P₂ orbitals. Then the permanents for the total number of bosons are obtained by combining them,

$$|\mathbf{n}^{n_{\text{exc}}}\rangle = \left( \sum_{i=1}^{N_1^{n_{\text{exc}}}} |n_1, n_2, \dots, n_{M_1}\rangle_i \right) \otimes \left( \sum_{i=1}^{N_2^{n_{\text{exc}}}} |n_{M_1+1}, \dots, n_M\rangle_i \right), \quad (A2)$$

where |nₐₐₓ⟩ is an array of dimension M × (N₁ⁿₐₐₓ × N₂ⁿₐₐₓ), with N₁ⁿₐₐₓ (N₂ⁿₐₐₓ) the total number of configurations obtained from arranging the N - nₐₐₓ (nₐₐₓ) particles in the P₁ (P₂) orbitals. We evaluate the indexes of the configurations resulting from the P₁ subsystem, Jₚ₁ⁿₐₐₓ, and from the P₂ subsystem, Jₚ₂ⁿₐₐₓ, applying equation (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ₚ₁ⁿₐₐₓ, Jₚ₂ⁿₐₐₓ, nₐₐₓ}, which stores the position of the configuration in the configurational vector. This scheme is thus applied for all excitations, nₐₐₓ, 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ₚ₁ⁿₐₐₓ for the P₁ subsystem, (ii) we evaluate the index Jₚ₁ⁿₐₐₓ for the P₁ subsystem, (iii) we know or evaluate the excitation, nₐₐₓ, of the configuration (iv) we access to the index of the coefficient which is stored and the three component array at position {Jₚ₁ⁿₐₐₓ, Jₚ₂ⁿₐₐₓ, nₐₐₓ}. 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 reduce the size of the configurational space, which made such a storage manageable for the applications done so far.
## Appendix B. Applying operators in second quantization

In second quantization, the action of the Hamiltonian of equation (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 [118, 119]. Concerning the one-body term, we have,

$$b_i^\dagger b_j |n_1, \dots, n_j, \dots, n_i, \dots, n_M\rangle = \sqrt{n_j} \sqrt{n_i+1} |n_1, \dots, n_j-1, \dots, n_i+1, \dots, n_M\rangle. \quad (B1)$$

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,

$$b_i^\dagger b_j \left[ \sum_{I \in \mathcal{V}} C_I |\Phi_I\rangle \right] = \sum_{I \in \mathcal{V}} C_I \sqrt{n_j} \sqrt{n_i+1} |\Phi_I\rangle, \quad (B2)$$

with $|\Phi_I'\rangle$ the new configuration resulting from the action of $b_i^\dagger b_j$ on the initial configuration $|\Phi_I\rangle$. The result can be interpreted as a reordering of the configuration in the wavefunction, as in equation (B2) or inversely to a reordering of the coefficients with a new factor $(\sqrt{n_j} \sqrt{n_i+1})$ if the configuration are reorganized in the initial order, i.e.,

$$b_i^\dagger b_j \left[ \sum_{I \in \mathcal{V}} C_I |\Phi_I\rangle \right] = \sum_{I \in \mathcal{V}} C_I' |\Phi_I\rangle. \quad (B3)$$

In the basis of the configurational states, $\{|\Phi_I\rangle\}$, 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^\dagger b_j$, is obtained as,

$$C_I' = \langle \Phi_I | b_i^\dagger b_j \left[ \sum_{J \in \mathcal{V}} C_J | \Phi_J \rangle \right]. \quad (B4)$$

In practice, we apply $b_i^\dagger b_j$ on the bra $\langle\Phi_I|$, which provides a new configurational state $\langle\Phi_J|$. The configurational vectors are orthonormal, and only the configuration $|\Phi_J\rangle$ in the sum remains from the projection. Thus, we evaluate the index of $\langle\Phi_J|$ to directly access the coefficient $C_J$, i.e.,

$$C_I' = C_J \sqrt{n_i} \sqrt{n_{j+1}}. \quad (B5)$$

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 explained above can be applied. In addition, a test to check if the resulting configuration remains in the RAS space has to be performed. 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^\dagger 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 particles in $P_2$, 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 maximal number of excited particles allowed for the general RAS scheme (section 2.2.2.4), and both (iii) and (iv) are not occuring for the scheme with only even excitations (section 2.2.2.3). The action of the one-body operator of the Hamiltonian is now straightforward, the coefficient vectors obtained by applying the $b_i^\dagger b_j$ operators are multiplied by the corresponding matrix element $h_j^I$ (equation (9)) and summed for each couple of $\{i, j\}$, with the restriction mentioned above for the RAS wavefunction. The two-body operator, see equation (10), included in the Hamiltonian of equation (8) and the four- and six-order tensors specific to the RAS schemes (equations (44) and (45)) 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 $b_i^\dagger b_j^\dagger b_k b_l$, we have the equalities $b_i^\dagger b_j^\dagger b_k b_l = b_i^\dagger b_j^\dagger b_l b_k = b_j^\dagger b_i^\dagger b_k b_l = b_i^\dagger b_j^\dagger b_l b_k$.

## Appendix C. Numerical implementation for the TD-RAS equations

The EOM for the TD-RASSCF-B and F methods, equations (15), (25) and (38) or (43), 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 $\eta_{i''}^{ij}$. These elements are evaluated from equations (38) or (43) depending on the RAS scheme used, but both are solved in the same way. We recall that the matrix $\eta$, with elements $\eta_{i''}^{ij}$, is anti-hermitian and thus $\eta_{i''}^{ij} = -(\eta_{i'}^{ij})^*$, here $i'$ and $j''$ correspond to orbitals of the $P_1$ and $P_2$
subspace, respectively. Writing down equations (38) or (43) for any set of orbitals {$i', j''$}, provides a system of $M_1 \times M_2$ linear equations with $M_1 \times M_2$ unknowns. Introducing a composite index for the couple of {$i', j''$}, this system can be written in a matrix form,

$$ \underline{\mathbf{A}} \cdot \underline{\mathbf{X}} = \underline{\mathbf{B}}, \qquad (C1) $$

where the matrix $\underline{\mathbf{A}}$, of dimension ($M_1 \times M_2, M_1 \times M_2$), contains the values of $A_{k'k'i'}^{l'j''}$ or $\zeta_{k'k'i'}^{l'j''}$, for all set of {$i', j''$} and {$l', k''$}. The vector $\underline{\mathbf{B}}$, of dimension ($M_1 \times M_2$), contains the rhs of equation (38) or (43) for each set of {$i', j''$} and the vector $\underline{\mathbf{X}}$ with the same dimension as $\underline{\mathbf{B}}$ contains the unknown values of $\eta_{i'}^{j''}$ and the matrix elements of the one-body operator. The system of linear equations, equation (C1), can be solved using standard numerical routines included, for instance, in the LAPACK library [120]. The values for $\eta_{i'}^{j''}$ are trivially obtained from the elements $\underline{\mathbf{X}}$ {$i', j''$} of the vector $\underline{\mathbf{X}}$,

$$ \eta_{i'}^{j''} = \begin{cases} i(\underline{\mathbf{X}}\{i', j''\} - \underline{\mathbf{h}}_{i'}^{j''}) & \text{for equation (38),} \\ -i(\underline{\mathbf{X}}\{i', j''\} + \underline{\mathbf{h}}_{i'}^{j''}) & \text{for equation (43).} \end{cases} \qquad (C2) $$

After evaluating the matrix elements $\eta_{i'}^{j''}$, the time derivative of the coefficients {$\dot{\underline{\mathbf{C}}}_i$} can be computed from equation (16) and the contribution of the $\mathcal{P}$-space orbitals to the time derivative of the orbitals is obtained from,

$$ \hat{\underline{\mathbf{P}}} |\dot{\underline{\mathbf{\phi}}}_i\rangle = \sum_j^M |\underline{\mathbf{\phi}}_j\rangle \eta_i^j. \qquad (C3) $$

It remains to evaluate the contribution from the orbitals of the $\mathcal{Q}$-space, i.e. $\hat{\underline{\mathbf{Q}}} |\dot{\underline{\mathbf{\phi}}}_i\rangle$ from equation (25). This latter can be expressed in a matrix form,

$$ i \hat{\underline{\mathbf{Q}}} \underline{\mathbf{\rho}} \dot{\underline{\mathbf{X}}} = \hat{\underline{\mathbf{Q}}} [\underline{\mathbf{\rho}} \tilde{\underline{\mathbf{h}}} + \tilde{\underline{\mathbf{W}}}], \qquad (C4) $$

with $\underline{\mathbf{\rho}}$ the one-body reduced density matrix, $\dot{\underline{\mathbf{X}}}$ a vector collecting the time derivative of the orbitals, $\tilde{\underline{\mathbf{h}}}$ and $\tilde{\underline{\mathbf{W}}}$ are both vectors with elements $\underline{\mathbf{h}}(\underline{\mathbf{r}}, t)|\underline{\mathbf{\phi}}_i\rangle$ and $\sum_{jlk} \underline{\mathbf{W}}_l^k |\underline{\mathbf{\phi}}_j\rangle \underline{\mathbf{\rho}}_{ik}^{jl}$, respectively. To obtain equation (C4), we used the fact that $\underline{\mathbf{\rho}}$ commutes with the projector $\hat{\underline{\mathbf{Q}}}$, as easily seen from the equality $\hat{\underline{\mathbf{Q}}} = \hat{\mathbf{1}} - \hat{\underline{\mathbf{P}}}$.

$$ \begin{aligned} \hat{\underline{\mathbf{Q}}} |\dot{\underline{\mathbf{\phi}}}_i\rangle &= -i\hat{\underline{\mathbf{Q}}} \left[ \underline{\mathbf{h}} |\underline{\mathbf{\phi}}_i\rangle + \sum_{jklm} (\underline{\mathbf{\rho}}^{-1})_i^m \rho_{mk}^{jl} \underline{\mathbf{W}}_l^k |\underline{\mathbf{\phi}}_j\rangle \right] \\ &= -i(\hat{\mathbf{1}} - \hat{\underline{\mathbf{P}}}) \left[ \underline{\mathbf{h}} |\underline{\mathbf{\phi}}_i\rangle + \sum_{jklm} (\underline{\mathbf{\rho}}^{-1})_i^m \rho_{mk}^{jl} \underline{\mathbf{W}}_l^k |\underline{\mathbf{\phi}}_j\rangle \right]. \end{aligned} \qquad (C5) $$

The right hand side of the above equation is similar to the one that is solved in the MCTDH-based methods [67, 68, 70] and to avoid singularities in the inverse of the one-body reduced density matrix and in equation (C1) for the matrix $\underline{\mathbf{A}}$, we follow the numerical implementation used for the MCTDH method [121]. The method uses a regularized one-body density matrix defined as [121],

$$ \underline{\mathbf{\rho}}^{\text{reg}} = \underline{\mathbf{\rho}} + \epsilon e^{-\underline{\mathbf{\rho}}/\epsilon}, \qquad (C6) $$

with $\epsilon$ a small real number (see appendix E.1). The same regularization is used for the matrix $\underline{\mathbf{A}}$. The numerical implementation of the projector onto the $\mathcal{P}$-space orbitals also follows [121] to ensure stability of the orthonormality of the orbital during the propagation. The projector $\hat{\underline{\mathbf{P}}}$ is numerically expressed by,

$$ \hat{\underline{\mathbf{P}}} = \sum_{ij}^{M} |\underline{\mathbf{\phi}}_i\rangle (\underline{\mathcal{O}})_{ij}^{-1} \langle \underline{\mathbf{\phi}}_j|, \qquad (C7) $$

where $\underline{\mathcal{O}}_{ij} = \langle \underline{\mathbf{\phi}}_i | \underline{\mathbf{\phi}}_j \rangle$ is the overlap matrix of the time-dependent orbitals. Thus, the projector defined in equation (C7) remains orthonormal even if the orthonormality of the orbitals is lost as long as they remain linearly independent.

## Appendix D. Numerical efficiency of the method

Comparing the efficiency between different methods is a difficult task as it depends on 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_{\text{grid}}$ the number of grid points that are used to describe the time-dependent orbital in the time-independent basis, usually a DVR [122], 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_{kl}^{ij}$ (equation (10)) and the two-body reduced density matrix $\rho_{ik}^{ij}$ (see text above equation (25)) are computed. These updates require $M^4 N_{\text{grid}}^2$ and $M^4 \dim(\mathcal{V}_{\text{FCI}})$ operations, respectively, with $\mathcal{V}_{\text{FCI}}$ the configurational space of the MCTDHB wavefunction evaluated from equation (3). Then computing the time derivative of the coefficients and the orbitals require on the order of $M^4 \dim(\mathcal{V}_{\text{FCI}})$ and $M^4 N_{\text{grid}}^2$ operations, respectively. The total cost is thus, approximately, $2M^4 (N_{\text{grid}}^2 + \dim(\mathcal{V}_{\text{FCI}}))$. Considering now the case of the TD-RASSCF-B method. The evaluation of the matrix elements of the 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., on the order of $M^4 N_{\text{grid}}^2$ operations, respectively. The evaluation of the time derivative of the coefficients and the matrix elements of the two-body reduced density matrix scale as $M^4 \dim(\mathcal{V})$, with $\mathcal{V}$ the configurational space for the considered RAS scheme. In addition, we also need to solve the P-space equations, which requires on the order of $M^4$ operations for excitation schemes with only even excitations (see equation (38)) and $M^6 \dim(\mathcal{V}_{N_{\text{max}}})$ for the general RAS scheme (see equation (43)), with $\dim(\mathcal{V}_{N_{\text{max}}})$ the number of configurations including $N_{\text{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 equation (5) and $\dim(\mathcal{V}_{N_{\text{max}}})$ is the last term of the summation. The dimension of the configurational space, $\mathcal{V}$, including only even excitations can be evaluated in a similar way,

$$ \dim(\mathcal{V}) = \binom{N + M_1 - 1}{N} + \sum_{k=1}^{N_{\text{max}}/2} \binom{2k + M_2 - 1}{2k} \binom{(N - 2k) + M_1 - 1}{N - 2k} \quad (\text{D1}) $$

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_{\text{grid}}^2 + \dim(\mathcal{V}) + M^2 \dim(\mathcal{V}_{N_{\text{max}}}))$ and in the case of only even excitations it scales as $2M^4 (N_{\text{grid}}^2 + \dim(\mathcal{V}) + 1/2)$. To compare the numerical cost between the MCTDHB and TD-RASSCF-B methods, we can introduce $\Delta(\mathcal{O}_p)$, the difference between the MCTDHB and TD-RASSCF-B operations to remove the constant number of operation resulting from $N_{\text{grid}}$,

$$ \Delta(\mathcal{O}_p) = \begin{cases} 2M^4 (\dim(\mathcal{V}_{\text{FCI}}) - \dim(\mathcal{V}) - 1/2) & \text{even excitation,} \\ 2M^4 (\dim(\mathcal{V}_{\text{FCI}}) - \dim(\mathcal{V}) - M^2 \dim(\mathcal{V}_{N_{\text{max}}})/2) & \text{general scheme.} \end{cases} \quad (\text{D2}) $$

From the expression of $\Delta(\mathcal{O}_p)$, 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, $\Delta(\mathcal{O}_p)$ 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 sixth-order tensor of the P-space equations, equation (45), can provide an overhead for the computation. To illustrate the computational efficiency we evaluate $\Delta(\mathcal{O}_p)$ for 10, 50 and 100 bosons in $M = 2$ to 8 orbitals, see figure D1. For the TD-RASSCF-B, we consider the case of a single $P_1$ orbital and $M = 1$ orbitals in $P_2$. From the figure, 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_{\text{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 figure D1, this value is rather large, i.e. $N_{\text{max}} = 6$ for 10, $N_{\text{max}} = 40$ for 50 particles and $N_{\text{max}} = 90$ for 100 particles for the schemes depicted in figures D1(a), (b) and (c), respectively.

# Appendix E. Numerical parameters and more illustrative examples

## E.1. Numerical parameters

To solve numerically the EOM of the MCTDHB 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 [69]. 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 table 1 (see section 3) and II (see appendix E.2). We numerically integrate the EOM using different integration algorithms, namely the 4th order runge-kutta (RK) [123], the adaptive time-step 5th order RK and the Adams-Bashforth-Moulton (ABM) predictor-corrector as also implemented in the Heidelberg MCTDH package [124]. Each simulation always used the same integration routine for both coefficients and orbitals. 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. The error tolerance for the integration has been fixed to $10^{-8}$, and we controlled the convergence of our results with respect to this parameter by performing simulations with higher and lower error tolerances. During a calculation, the inverse of the reduced one-body density matrix can be singular when one or more of the natural orbitals are unoccupied. To prevent the singularity we used a regularization scheme, see appendix C, with a
<figure><img src="image_8.png" /></figure>

**Figure D1.** Difference between the number of operations used to evaluate the time derivative of the wavefunction with the MCTDHB and TD-RASSCF-B method. The difference $\Delta(\text{Op})$ is defined by equation (D2) in appendix D. We report the value of $\Delta(\text{Op})$ as a function of the number of orbitals for a system consisting of (a) $N = 10$ particles, (b) $N = 50$ particles and (c) $N = 100$ particles. The circles indicate the results for general RAS schemes (left axis) and the diamond symbols depict the results for the case of RAS schemes with only even excitations ((red) right axis). For the RAS schemes with only even excitations, the number of operations is always fewer than for the MCTDHB method. For the general RAS schemes, for high excitation schemes and large numbers of orbitals the number of operations can be larger than the one of the MCTDHB method. For instance, considering $N = 10$ particles (panel (a)) and the -SDTQ56 scheme (blue circle) the number of operations with the TD-RASSCF-B method is higher than with the MCTDHB method for $M \ge 4$ orbitals. The same happens for $N = 50$ particles and the $-40$ scheme with $M \ge 6$ orbitals and for $N = 100$ particles and the scheme $-90$ with $M \ge 5$ orbitals.

regularization parameter $\epsilon = 10^{-8}$. Again, the convergence with respect to the value of $\epsilon$ was ensured by performing calculations with $10^{-7} \le \epsilon \le 10^{-12}$ and no significant difference in the results was found. For instance, the reduced one-body density matrix is singular when computing the ground state of the Hamiltonian of equation (50), but the energy is known analytically ($5E_0$) as the particles do not interact with each other. The GS energies obtained numerically with $\epsilon = 10^{-8}$ are all accurate up to $10^{-13}$ for any number of orbitals and excitations. We calculate the GS energy using imaginary time propagation [125] of the EOM and give all energies in units of $E_0$.

## E.2. More illustrative examples

In this Appendix, we consider the same system as the one presented in section 3 except that the interaction strength is now 10 times larger with $\lambda = 0.1$, see equation (46). We perform an analysis of the GS energy obtained with the different methods, similarly to section 3, and the results are reported in table 2. The reference for the correlation energy (see equation (48)) is $E_{\text{ref}} = -1.34$, see table 2. This value is much larger than the one obtained in section 3 with $\lambda = 0.01$ 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 occupy 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 is not possible as only one orbital is used to describe the wavefunction. The orbitals that diagonalize the reduced density matrix and their respective eigenvalues, i.e., their natural
<table><caption>**Table 2.** Same as table 1 but for ground-state energies of 100 bosons trapped in a 1D harmonic potential interacting through a contact potential with a strength $\lambda = 0.1$.</caption><thead><tr><th rowspan="2">Method</th><th colspan="10">Orbitals</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th></tr></thead><tbody><tr><td><strong>All excitations</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>MCTDHB</td><td>193.5509587<br>(1)</td><td>193.0154216<br>(101)</td><td>192.6308389<br>(5151)</td><td>192.3920265<br>(176851)</td><td>192.2138048<br>(4598126)</td><td>—<br>(96560646)</td><td>—<br>(1705904746)</td><td>—<br>(26075972546)</td></tr><tr><td>-SDT</td><td>—</td><td>193.0783470<br>(4)</td><td>192.7594351<br>(10)</td><td>192.5665491<br>(20)</td><td>192.4187690<br>(35)</td><td>192.3153377<br>(56)</td><td>192.2315684<br>(84)</td><td>192.1681159<br>(120)</td></tr><tr><td>-SDTQ5</td><td>—</td><td>193.0310688<br>(6)</td><td>192.6591396<br>(21)</td><td>192.4322339<br>(56)</td><td>192.2605383<br>(126)</td><td>192.1396169<br>(252)</td><td>192.0434375<br>(462)</td><td>191.9701535<br>(792)</td></tr><tr><td>-SDTQ567</td><td>—</td><td>193.0191782<br>(8)</td><td>192.6370892<br>(36)</td><td>192.4017271<br>(120)</td><td>192.2250102<br>(330)</td><td>192.1000115<br>(792)</td><td>192.0013608<br>(1716)</td><td>191.9259461<br>(3432)</td></tr><tr><td>-SDTQ56789</td><td>—</td><td>193.0162747<br>(10)</td><td>192.6321926<br>(55)</td><td>192.3943825<br>(220)</td><td>192.2165110<br>(715)</td><td>192.0903987<br>(2002)</td><td>191.9911974<br>(5005)</td><td>191.9152346<br>(11440)</td></tr><tr><td>-10</td><td>—</td><td>193.0158613<br>(11)</td><td>192.6315747<br>(66)</td><td>192.3933448<br>(286)</td><td>192.2153160<br>(1001)</td><td>192.0890211<br>(3003)</td><td>191.9894655<br>(8008)</td><td>191.9133983<br>(19448)</td></tr><tr><td>-15</td><td>—</td><td>193.0154291<br>(16)</td><td>192.6308506<br>(136)</td><td>192.3920555<br>(816)</td><td>192.2138378<br>(3876)</td><td>192.0872814<br>(15504)</td><td>191.9879161<br>(54264)</td><td>191.9117429<br>(170544)</td></tr><tr><td>-20</td><td>—</td><td>193.0154220<br>(21)</td><td>192.6308392<br>(231)</td><td>192.3920271<br>(1771)</td><td>192.2138055<br>(10626)</td><td>192.0872403<br>(53130)</td><td>191.9878729<br>(230230)</td><td>—<br>(888030)</td></tr><tr><td>-23</td><td>—</td><td>193.0154216<br>(24)</td><td>192.6308389<br>(300)</td><td>192.3920265<br>(2600)</td><td>192.2138049<br>(17550)</td><td>192.0872393<br>(98280)</td><td>191.9878720<br>(475020)</td><td>—<br>(2035800)</td></tr><tr><td>-25</td><td>—</td><td>193.0154216<br>(26)</td><td>192.6308389<br>(351)</td><td>192.3920265<br>(3276)</td><td>192.2138048<br>(23751)</td><td>192.0872393<br>(142506)</td><td>191.9878719<br>(736281)</td><td>—<br>(3365856)</td></tr><tr><td><strong>Even excitations</strong></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>-D</td><td>—</td><td>193.1320396<br>(2)</td><td>192.8258797<br>(4)</td><td>192.6453389<br>(7)</td><td>192.5031705<br>(11)</td><td>192.4050140<br>(16)</td><td>192.3243169<br>(22)</td><td>192.2636708<br>(29)</td></tr><tr><td>-DQ</td><td>—</td><td>193.0449780<br>(3)</td><td>192.6788633<br>(9)</td><td>192.4619220<br>(22)</td><td>192.2951569<br>(46)</td><td>192.1793448<br>(86)</td><td>192.0861393<br>(148)</td><td>192.0157234<br>(239)</td></tr><tr><td>-DQ6</td><td>—</td><td>193.0227865<br>(4)</td><td>192.6439550<br>(16)</td><td>192.4175380<br>(50)</td><td>192.2457797<br>(130)</td><td>192.1259116<br>(296)</td><td>192.0303689<br>(610)</td><td>191.9579281<br>(1163)</td></tr><tr><td>-DQ68</td><td>—</td><td>193.0171558<br>(5)</td><td>192.6356968<br>(25)</td><td>192.4065816<br>(95)</td><td>192.2338183<br>(295)</td><td>192.1129247<br>(791)</td><td>192.0169500<br>(1897)</td><td>191.9440484<br>(4166)</td></tr><tr><td>-10</td><td>—</td><td>193.0158072<br>(6)</td><td>192.6338334<br>(36)</td><td>192.4039516<br>(161)</td><td>192.2310023<br>(581)</td><td>192.1098389<br>(1792)</td><td>192.0137924<br>(4900)</td><td>191.9407800<br>(12174)</td></tr><tr><td>-20</td><td>—</td><td>193.0154216<br>(11)</td><td>192.6333324<br>(121)</td><td>192.4031759<br>(946)</td><td>192.2301918<br>(5786)</td><td>192.1089342<br>(29458)</td><td>192.0128774<br>(129844)</td><td>191.9398292<br>(508937)</td></tr><tr><td>-30</td><td>—</td><td>193.0154216<br>(16)</td><td>192.6333323<br>(256)</td><td>192.4031756<br>(2856)</td><td>192.2301915<br>(24616)</td><td>192.1089338<br>(174624)</td><td>192.0128771<br>(1061208)</td><td>—<br>(5678340)</td></tr></tbody></table>
occupations, can be used to characterize the system. If the largest eigenvalue is of the same order as N, the system is condensed [109]. 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 [126] and the discussion in [83]). We find that, indeed, the natural occupation of the lowest natural orbital in the MCTDHB calculation using 5 orbitals decreases from an occupation of 99.987% for $\lambda = 0.01$ to an occupation of 99.465% for $\lambda = 0.1$. This slightly larger depletion of the condensate has a strong impact on the correlation energy, and 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 $\sim 10^{-1}$, see table 2, which means that the energy does not converge below this value. The relative correlation energy, see equation (49), $\mathcal{E}_{\text{rel}}^2$, $\mathcal{E}_{\text{rel}}^3$ and $\mathcal{E}_{\text{rel}}^4$ include 40.1%, 68.8% and 86.7% of $\mathcal{E}_{\text{ref}}$, respectively.

Starting the discussion with the general RAS schemes, see table 2, we note that to converge to the MCTDHB energies and thus include 100% of the relative correlation energies, $\mathcal{E}_{\text{rel}}^X$ with $X = 2$ to 5, large values of $N_{\text{max}}$ are required. For the -SDTQ5 RAS scheme we find that the correlation energy includes $\sim 97\%$ of $\mathcal{E}_{\text{rel}}^X$, the -10 RAS scheme includes $\sim 99.9\%$ of $\mathcal{E}_{\text{rel}}^X$, the -15 RAS scheme includes $\sim 99.997\%$ of $\mathcal{E}_{\text{rel}}^X$ and the -20 RAS scheme is converged with more than 99.9999% of $\mathcal{E}_{\text{rel}}^X$. These results are obtained irrespectively of the number of orbitals, i.e., $X = 2$ to 5. Even if large values of $N_{\text{max}}$ are used, the expansion of the wavefunction using the -20 RAS scheme includes $\sim 10^4$ configurations for 5 orbitals while the MCTDHB wavefunction includes $\sim 4.6 \times 10^6$ configurations. We also use RAS schemes with only even excitations and we report the results in table 2. Similarly to the case with $\lambda = 0.01$, see section 3, we find that, except for 2 orbitals, the energy does not converge to the MCTDHB energy, irrespective of the value of $N_{\text{max}}$ used. Thus the energies obtained with 3, 4 and 5 orbitals include 99.7%, 99.0% and 98.8% of the respective $\mathcal{E}_{\text{rel}}^X$ 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 $\sim 10^{-1}$ and a convergence to $\sim 10^{-2}$ is achieved with the -SDTQ5 excitation scheme and the -DQ excitation scheme (table 2). As previously, the configurational space of the MCTDHB 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 2. 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 $\sim 0.3$ below the energy of the MCTDHB method with 5 orbitals by using the -15 RAS scheme with 8 orbitals, see table 2.

## References

[1] Anderson M H, Ensher J R, Matthews M R, Wieman C E and Cornell E A 1995 *Science* **269** 198

[2] Bradley C C, Sackett C A, Tollett J J and Hulet R G 1995 *Phys. Rev. Lett.* **75** 1687

[3] Davis K B, Mewes M O, Andrews M R, van Druten N J, Durfee D S, Kurn D M and Ketterle W 1995 *Phys. Rev. Lett.* **75** 3969

[4] Serwane F, Zürn G, Lompe T, Ottenstein T B, Wenz A N and Jochim S 2011 *Science* **332** 336

[5] Kaufman A M, Lester B J, Reynolds C M, Wall M L, Foss-Feig M, Hazzard K R A, Rey A M and Regal C A 2014 *Science* **345** 306

[6] Preiss P M, Ma R, Tai M E, Lukin A, Rispoli M, Zupancic P, Lahini Y, Islam R and Greiner M 2015 *Science* **347** 1229

[7] Bloch I 2005 *Nat. Phys.* **1** 23

[8] Jotzu G, Messer M, Desbuquois R, Lebrat M, Uehlinger T, Greif D and Esslinger T 2014 *Nature* **515** 237

[9] Greiner M, Mandel O, Esslinger T, Hansch T W and Bloch I 2002 *Nature* **415** 39

[10] Courteille P, Freeland R S, Heinzen D J, van Abeelen F A and Verhaar B J 1998 *Phys. Rev. Lett.* **81** 69

[11] Inouye S, Andrews M R, Stenger J, Miesner H-J, Stamper-Kurn D M and Ketterle W 1998 *Nature* **392** 151

[12] Bloch I, Dalibard J and Nascimbene S 2012 *Nat. Phys.* **8** 267

[13] Anderson B P and Kasevich M A 1998 *Science* **282** 1686

[14] Steinhauer J 2016 *Nat. Phys.* **12** 959

[15] Sala S, Förster J and Saenz A 2017 *Phys. Rev. A* **95** 011403

[16] Lühmann D-S, Weitenberg C and Sengstock K 2015 *Phys. Rev. X* **5** 031016

[17] Bloch I, Dalibard J and Zwerger W 2008 *Rev. Mod. Phys.* **80** 885

[18] Lieb E H and Liniger W 1963 *Phys. Rev.* **130** 1605

[19] Lieb E H 1963 *Phys. Rev.* **130** 1616

[20] Girardeau M 1960 *J. Math. Phys.* **1** 516

[21] Yukalov V I and Girardeau M D 2005 *Laser Phys. Lett.* **2** 375

[22] Calogero F 1969 *J. Math. Phys.* **10** 2191

[23] Cohen L and Lee C 1985 *J. Math. Phys.* **26** 3105

[24] Yan J 2003 *J. Stat. Phys.* **113** 623

[25] Lode A U J, Sakmann K, Alon O E, Cederbaum L S and Streltsov A I 2012 *Phys. Rev. A* **86** 063606

[26] Lode A U J 2015 * Tunneling Dynamics in Open Ultracold Bosonic Systems * 1st edn (Berlin: Springer)

[27] Fasshauer E and Lode A U J 2016 *Phys. Rev. A* **93** 033635

[28] Baym G and Pethick C J 1996 *Phys. Rev. Lett.* **76** 6

[29] Bogoliubov N N 1947 *J. Phys.* **11** 23

[30] Lee T D and Yang C N 1957 *Phys. Rev.* **105** 1119

[31] Lee T D, Huang K and Yang C N 1957 *Phys. Rev.* **106** 1135
[32] Fisher M P A, Weichman P B, Grinstein G and Fisher D S 1989 *Phys. Rev.* B **40** 546

[33] Dutta O, Gajda M, Hauke P, Lewenstein M, Lühmann D-S, Malomed B A, Sowiński T and Zakrzewski J 2015 *Rep. Prog. Phys.* **78** 066001

[34] White S R and Noack R M 1992 *Phys. Rev. Lett.* **68** 3487

[35] White S R 1992 *Phys. Rev. Lett.* **69** 2863

[36] White S R 1993 *Phys. Rev.* B **48** 10345

[37] Schollwöck U 2005 *Rev. Mod. Phys.* **77** 259

[38] Vidal G 2004 *Phys. Rev. Lett.* **93** 040502

[39] Daley A J, Kollath C, Schollwöck U and Vidal G 2004 *J. Stat. Mech.* P04005

[40] White S R and Feiguin A E 2004 *Phys. Rev. Lett.* **93** 076401

[41] McMillan W L 1965 *Phys. Rev.* **138** A442

[42] Anderson J B 1975 *J. Chem. Phys.* **63** 1499

[43] Reynolds P J, Ceperley D M, Alder B J and Lester W A 1982 *J. Chem. Phys.* **77** 5593

[44] Blume D and Greene C H 2001 *Phys. Rev.* A **63** 063601

[45] Astrakharchik G E, Blume D, Giorgini S and Granger B E 2004 *J. Phys.* B **37** S205

[46] Bijl A 1940 *Physica* **7** 869

[47] Jastrow R 1955 *Phys. Rev.* **98** 1479

[48] Esry B D 1997 *Phys. Rev.* A **55** 1147

[49] Streltsov A I, Alon O E and Cederbaum L S 2006 *Phys. Rev.* A **73** 063626

[50] Löwdin P-O 1955 *Phys. Rev.* **97** 1474

[51] Coester F 1958 *Nucl. Phys.* **7** 421

[52] Coester F and Kümmel H 1960 *Nucl. Phys.* **17** 477

[53] Čížek J 1966 *J. Chem. Phys.* **45** 4256

[54] Cederbaum L S, Alon O E and Streltsov A I 2006 *Phys. Rev.* A **73** 043609

[55] Lysaght M A, van der Hart H W and Burke P G 2009 *Phys. Rev.* A **79** 053411

[56] Hochstuhl D and Bonitz M 2012 *Phys. Rev.* A **86** 053424

[57] Pabst S, Greenman L, Mazziotti D A and Santra R 2012 *Phys. Rev.* A **85** 023411

[58] Kvaal S 2012 *J. Chem. Phys.* **136** 194109

[59] Sato T and Ishikawa K L 2013 *Phys. Rev.* A **88** 023402

[60] Bauch S, Sørensen L K and Madsen L B 2014 *Phys. Rev.* A **90** 062508

[61] Popmintchev T, Chen M-C, Arpin P, Murnane M M and Kapteyn H C 2010 *Nat. Photon.* **4** 822

[62] Calegari F et al 2014 *Science* **346** 336

[63] Kraus P M et al 2015 *Science* **350** 790

[64] Zanghellini J, Kitzler M, Fabian C, Brabec T and Scrinzi A 2003 *Laser Phys.* **13** 1064

[65] Kato T and Kono H 2004 *Chem. Phys. Lett.* **392** 533

[66] Nest M, Klamroth T and Saalfrank P 2005 *J. Chem. Phys.* **122** 124102

[67] Haxton D J, Lawler K V and McCurdy C W 2011 *Phys. Rev.* A **83** 063416

[68] Meyer H-D, Manthe U and Cederbaum L S 1990 *Chem. Phys. Lett.* **165** 73

[69] Beck M, Jäckle A, Worth G and Meyer H-D 2000 *Phys. Rep.* **324** 1

[70] Alon O E, Streltsov A I and Cederbaum L S 2008 *Phys. Rev.* A **77** 033613

[71] Alon O E, Streltsov A I and Cederbaum L S 2007 *Phys. Rev.* A **76** 062501

[72] Alon O E, Streltsov A I, Sakmann K, Lode A U, Grond J and Cederbaum L S 2012 *Chem. Phys.* **401** 2

[73] Lode A U J and Bruder C 2016 *Phys. Rev.* A **94** 013616

[74] Lode A U J 2016 *Phys. Rev.* A **93** 063601

[75] Wang H and Thoss M 2003 *J. Chem. Phys.* **119** 1289

[76] Manthe U 2008 *J. Chem. Phys.* **128** 164116

[77] Vendrell O and Meyer H-D 2011 *J. Chem. Phys.* **134** 044135

[78] Krönke S, Cao L, Vendrell O and Schmelcher P 2013 *New J. Phys.* **15** 063018

[79] Cao L, Krönke S, Vendrell O and Schmelcher P 2013 *J. Chem. Phys.* **139** 134103

[80] Schmitz R, Krönke S, Cao L and Schmelcher P 2013 *Phys. Rev.* A **88** 043601

[81] Sakmann K, Streltsov A I, Alon O E and Cederbaum L S 2009 *Phys. Rev. Lett.* **103** 220601

[82] Sakmann K, Streltsov A I, Alon O E and Cederbaum L S 2014 *Phys. Rev.* A **89** 023602

[83] Sakmann K, Streltsov A I, Alon O E and Cederbaum L S 2008 *Phys. Rev.* A **78** 023615

[84] Klaiman S and Alon O E 2015 *Phys. Rev.* A **91** 063613

[85] Streltsov A I, Alon O E and Cederbaum L S 2008 *Phys. Rev. Lett.* **100** 130401

[86] Streltsov A I, Alon O E and Cederbaum L S 2009 *Phys. Rev.* A **80** 043616

[87] Streltsova O I, Alon O E, Cederbaum L S and Streltsov A I 2014 *Phys. Rev.* A **89** 061602

[88] Lode A U, Streltsov A I, Sakmann K, Alon O E and Cederbaum L S 2012 *Proc. Natl Acad. Sci. USA* **109** 13521

[89] Beinke R, Klaiman S, Cederbaum L S, Streltsov A I and Alon O E 2015 *Phys. Rev.* A **92** 043627

[90] Olsen J, Roos B O, Jørgensen P and Jensen H J A 1988 *J. Chem. Phys.* **89** 2185

[91] Miyagi H and Madsen L B 2013 *Phys. Rev.* A **87** 062511

[92] Miyagi H and Madsen L B 2014 *Phys. Rev.* A **89** 063416

[93] Haxton D J and McCurdy C W 2015 *Phys. Rev.* A **91** 012509

[94] Miyagi H and Madsen L B 2014 *J. Chem. Phys.* **140** 164309

[95] Kramer P and Saraceno M 1981 *Geometry of the Time-Dependent Variational Principle in Quantum Mechanics* (Berlin: Springer)

[96] Cederbaum L S 2016 *J. Phys. Chem.* A **120** 3009

[97] Broeckhove J, Lathouwers L, Kesteloot E and Leuven P V 1988 *Chem. Phys. Lett.* **149** 547

[98] Caillat J, Zanghellini J, Kitzler M, Koch O, Kreuzer W and Scrinzi A 2005 *Phys. Rev.* A **71** 012712

[99] Brillouin L 1933 *J. Phys. Radium* **4** 333

[100] Levy B and Berthier G 1968 *Int. J. Quantum Chem.* **2** 307–19

[101] Alon O E, Streltsov A I and Cederbaum L S 2007 *J. Chem. Phys.* **127** 154103

[102] Görlitz A et al 2001 *Phys. Rev. Lett.* **87** 130402

[103] Moritz H, Stöferle T, Köhl M and Esslinger T 2003 *Phys. Rev. Lett.* **91** 250402
[104] Tolra B L, O'Hara K M, Huckans J H, Phillips W D, Rolston S L and Porto J V 2004 *Phys. Rev. Lett.* **92** 190401

[105] Kinoshita T, Wenger T and Weiss D S 2005 *Phys. Rev. Lett.* **95** 190406

[106] Hofferberth S, Lesanovsky I, Fischer B, Schumm T and Schmiedmayer J 2007 *Nature* **449** 324

[107] Olshanii M 1998 *Phys. Rev. Lett.* **81** 938

[108] Szabo A and Ostlund N 1996 *Modern Quantum Chemistry* (Mineola, NY: Dover)

[109] Penrose O and Onsager L 1956 *Phys. Rev.* **104** 576

[110] Hilligsoe K M and Mølmer K 2005 *Phys. Rev. A* **71** 041602

[111] Bonneau M, Ruaudel J, Lopes R, Jaskula J-C, Aspect A, Boiron D and Westbrook C I 2013 *Phys. Rev. A* **87** 061603

[112] Khaykovich L, Schreck F, Ferrari G, Bourdel T, Cubizolles J, Carr L D, Castin Y and Salomon C 2002 *Science* **296** 1290

[113] Strecker K E, Partridge G B, Truscott A G and Hulet R G 2002 *Nature* **417** 150

[114] Burger S, Bongs K, Dettmer S, Ertmer W, Sengstock K, Sanpera A, Shlyapnikov G V and Lewenstein M 1999 *Phys. Rev. Lett.* **83** 5198

[115] Denschlag J *et al* 2000 *Science* **287** 97

[116] Aidelsburger M, Atala M, Lohse M, Barreiro J T, Paredes B and Bloch I 2013 *Phys. Rev. Lett.* **111** 185301

[117] Goldman N and Dalibard J 2014 *Phys. Rev. X* **4** 031027

[118] Streltsov A I, Alon O E and Cederbaum L S 2010 *Phys. Rev. A* **81** 022124

[119] Streltsov A I, Sakmann K, Alon O E and Cederbaum L S 2011 *Phys. Rev. A* **83** 043604

[120] Anderson *et al* 1999 *LAPACK Users' Guide* 3rd edn (Philadelphia, PA: Society for Industrial and Applied Mathematics)

[121] Beck M and Meyer H-D 1997 *Z. Phys. D* **42** 113

[122] Heather R W and Light J C 1983 *J. Chem. Phys.* **79** 147

[123] Press W H, Teukolsky S A, Vetterling W T and Flannery B P 1993 *Numerical Recipes in FORTRAN; The Art of Scientific Computing* 2nd edn (New York: Cambridge University Press)

[124] Worth G A, Beck M H, Jäckle A and Meyer H-D 2007 the MCTDH package, version 8.4, University of Heidelberg, Germany http://mctdh.uni-hd.de

[125] Kosloff R and Tal-Ezer H 1986 *Chem. Phys. Lett.* **127** 223

[126] Nozières P and James D S 1982 *J. Phys. France* **43** 1133