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

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 fermonic 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 ={\rm{constant}}$, with N the number of bosons and λ 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.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:

$\begin{eqnarray}&&{\rm{i}}{\hslash }\displaystyle \frac{\partial }{\partial t}| {\rm{\Psi }}(t)\rangle =\hat{H}(t)| {\rm{\Psi }}(t)\rangle ,\end{eqnarray} \tag{ 1 }$

with $\hat{H}(t)$ the many-body Hamiltonian of the system and $| {\rm{\Psi }}(t)\rangle$ the N-particle wavefunction. Hereafter we set ${\hslash }=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}({\bf{r}},{\rm{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 $| {{\rm{\Phi }}}_{I}(t)\rangle$,

$\begin{eqnarray}&&| {\rm{\Psi }}(t)\rangle =\displaystyle \sum _{I\in {{ \mathcal V }}_{\mathrm{FCI}}}{C}_{I}(t)| {{\rm{\Phi }}}_{I}(t)\rangle .\end{eqnarray} \tag{ 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 ${\rm{\dim }}({{ \mathcal V }}_{\mathrm{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 }}_{\mathrm{FCI}}$ can be evaluated as,

$\begin{eqnarray}&&{\rm{\dim }}({{ \mathcal V }}_{\mathrm{FCI}})=\left(\begin{array}{c}N+M-1\\ N\end{array}\right)=\displaystyle \frac{(N+M-1)!}{N!(M-1)!}.\end{eqnarray} \tag{ 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₁ and M₂ 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,

$\begin{eqnarray}&&| {\rm{\Psi }}(t)\rangle =\displaystyle \sum _{I\in { \mathcal V }}{C}_{I}(t)| {{\rm{\Phi }}}_{I}(t)\rangle ,\end{eqnarray} \tag{ 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}_{\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}_{\max }$. The dimension of the space ${ \mathcal V }$ of equation (4), is then given by

$\begin{eqnarray}&&{\rm{\dim }}({ \mathcal V })=\left(\begin{array}{c}N+{M}_{1}-1\\ N\end{array}\right)+\displaystyle \sum _{k=1}^{{N}_{\max }}\left(\begin{array}{c}k+{M}_{2}-1\\ k\end{array}\right)\left(\begin{array}{c}(N-k)+{M}_{1}-1\\ N-k\end{array}\right).\end{eqnarray} \tag{ 5 }$

Here, the first term is the total number of configurations obtained with the N bosons in the M₁ 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\leqslant k\leqslant {N}_{\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₂ orbitals and $(N-k)$ bosons in M₁ orbitals, see also appendix A.

Zoom In Zoom Out Reset image size

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₁ by M in equation (5). Note that this is also true for ${M}_{2}\ne 0$ and ${N}_{\max }=N$. Moreover, if only a single time-dependent orbital is considered, i.e., ${M}_{1}=1$ and ${M}_{2}=0$, the RAS wavefunction includes a single configuration with all particles in one orbital, which is equivalent 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 ${ \mathcal P }$-space is a subset of the total single-particle Hilbert space and we can define its orthogonal complement, ${ \mathcal 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 ${ \mathcal P }$-space is variationally optimized at each time and both ${ \mathcal P }$- and ${ \mathcal Q }$-space are time-dependent. We can define $\hat{P}$ and $\hat{Q}$, the time-dependent projectors onto the subspaces ${ \mathcal P }$ and ${ \mathcal Q }$, respectively, with the property $\hat{P}+\hat{Q}=\hat{1}$, the identity operator. The role of the ${ \mathcal Q }$-space emanates from the time-derivative of the ${ \mathcal P }$-space orbitals, that can be written as,

$\begin{eqnarray}&&| {\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 ,\end{eqnarray} \tag{ 6 }$

with one contribution from the ${ \mathcal P }$-space and one contribution from the ${ \mathcal 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 ${ \mathcal Q }$- and ${ \mathcal 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,

$\begin{eqnarray}&&S[\{{C}_{I}(t)\},\{| {\phi }_{i}(t)\rangle \},\{{\epsilon }_{j}^{i}(t)\}]={\int }_{{t}_{1}}^{{t}_{2}}\left[\langle {\rm{\Psi }}(t)| \hat{K}| {\rm{\Psi }}(t)\rangle +\displaystyle \sum _{{ij}}{\epsilon }_{j}^{i}(t)(\langle {\phi }_{i}(t)| {\phi }_{j}(t)\rangle -{\delta }_{{ij}})\right]{\rm{d}}{t},\end{eqnarray} \tag{ 7 }$

with $\hat{K}\equiv {\rm{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,k,\ldots$ are used to denote the orbitals of the ${ \mathcal P }$-space, the indexes $a,b,c,\ldots$ denote the orbitals of the ${ \mathcal Q }$-space and $p,q,r,\ldots \,$ are used for either ${ \mathcal P }$- or ${ \mathcal 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

$\begin{eqnarray}&&\hat{H}(t)=\displaystyle \sum _{{pq}}{h}_{q}^{p}(t){\hat{b}}_{p}^{\dagger }{\hat{b}}_{q}+\displaystyle \frac{1}{2}\displaystyle \sum _{{pqrs}}{v}_{{qs}}^{{pr}}(t){\hat{b}}_{p}^{\dagger }{\hat{b}}_{r}^{\dagger }{\hat{b}}_{s}{\hat{b}}_{q},\end{eqnarray} \tag{ 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

$\begin{eqnarray}&&{h}_{q}^{p}(t)=\int {\phi }_{p}^{* }({\bf{r}},t)h({\bf{r}},t){\phi }_{q}({\bf{r}},t){\rm{d}}{\bf{r}},\end{eqnarray} \tag{ 9 }$

and

$\begin{eqnarray}&&{v}_{{qs}}^{{pr}}(t)=\int \int {\phi }_{p}^{* }({\bf{r}},t){\phi }_{r}^{* }({\bf{r}}^{\prime} ,t)W({\bf{r}},{\bf{r}}^{\prime} ,t){\phi }_{q}({\bf{r}},t){\phi }_{s}({\bf{r}}^{\prime} ,t){\rm{d}}{\bf{r}}{\rm{d}}{\bf{r}}^{\prime} ,\end{eqnarray} \tag{ 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 {\rm{\Psi }}({t}_{1})\rangle =| \delta {\rm{\Psi }}({t}_{2})\rangle =0$. The variation of the action gives,

$\begin{eqnarray}&&\delta S={\int }_{{t}_{1}}^{{t}_{2}}\left[\langle \delta {\rm{\Psi }}| \hat{K}{\rm{\Psi }}\rangle +\langle \hat{K}{\rm{\Psi }}| \delta {\rm{\Psi }}\rangle +\displaystyle \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]{\rm{d}}{t},\end{eqnarray} \tag{ 11 }$

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

$\begin{eqnarray}&&| \delta {\rm{\Psi }}\rangle =\displaystyle \sum _{I\in { \mathcal V }}\delta {C}_{I}| {{\rm{\Phi }}}_{I}\rangle +\displaystyle \sum _{{pq}}{\hat{b}}_{p}^{\dagger }{\hat{b}}_{q}| {\rm{\Psi }}\rangle \langle {\phi }_{p}| \delta {\phi }_{q}\rangle .\end{eqnarray} \tag{ 12 }$

We can now proceed with the stationarity condition of the action with respect to the parameters $\{{C}_{I}\}$, $\{| {\phi }_{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 {\rm{\Psi }}|$ of the first expectation value in equation (11). Thus, the stationarity condition, $\delta S/\delta {C}_{I}^{* }=0$, readily leads to $\langle {{\rm{\Phi }}}_{I}| {\rm{i}}{\partial }_{t}-\hat{H}| {\rm{\Psi }}\rangle =0$. Moreover, the derivative of the wavefunction with respect to time reads,

$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}| {\rm{\Psi }}\rangle =\displaystyle \sum _{I\in { \mathcal V }}{\dot{C}}_{I}| {{\rm{\Phi }}}_{I}\rangle +\left[\displaystyle \sum _{{pq}}{\eta }_{q}^{p}{\hat{b}}_{p}^{\dagger }{\hat{b}}_{q}\right]| {\rm{\Psi }}\rangle ,\end{eqnarray} \tag{ 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 $| {\rm{\Psi }}\rangle$. Hereafter, the operator in bracket in equation (13) will be called $\hat{D}$, i.e.,

$\begin{eqnarray}&&\hat{D}\equiv \displaystyle \sum _{{pq}}{\eta }_{q}^{p}{\hat{b}}_{p}^{\dagger }{\hat{b}}_{q},\end{eqnarray} \tag{ 14 }$

for brevity. We can now rewrite the stationary condition $\delta S/\delta {C}_{I}^{* }=0$ using the explicit form of ${\partial }_{t}| {\rm{\Psi }}\rangle$, equation (13), as

$\begin{eqnarray}&&{\rm{i}}{\dot{C}}_{I}+\langle {{\rm{\Phi }}}_{I}| ({\rm{i}}\hat{D}-\hat{H})| {\rm{\Psi }}\rangle =0,\forall I\in { \mathcal V },\end{eqnarray} \tag{ 15 }$

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

$\begin{eqnarray}&&{\rm{i}}{\dot{C}}_{I}=\displaystyle \sum _{{ij}}({h}_{j}^{i}-{\rm{i}}{\eta }_{j}^{i})\langle {{\rm{\Phi }}}_{I}| {\hat{b}}_{i}^{\dagger }{\hat{b}}_{j}| {\rm{\Psi }}\rangle +\displaystyle \frac{1}{2}\displaystyle \sum _{{ijkl}}{v}_{{jl}}^{{ik}}\langle {{\rm{\Phi }}}_{I}| {\hat{b}}_{i}^{\dagger }{\hat{b}}_{k}^{\dagger }{\hat{b}}_{l}{\hat{b}}_{j}| {\rm{\Psi }}\rangle .\end{eqnarray} \tag{ 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 {{\rm{\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

$\begin{eqnarray}&&\displaystyle \sum _{q}| {\phi }_{q}\rangle \langle {{\rm{\Psi }}}_{i}^{q}| \left[\displaystyle \sum _{I\in { \mathcal V }}{\rm{i}}{\dot{C}}_{I}| {{\rm{\Phi }}}_{I}\rangle +({\rm{i}}\hat{D}-\hat{H})| {\rm{\Psi }}\rangle \right]+\displaystyle \sum _{j}{\epsilon }_{j}^{i}| {\phi }_{j}\rangle =0,\end{eqnarray} \tag{ 17 }$

with $\langle {{\rm{\Psi }}}_{i}^{q}| \,\equiv \,\langle {\rm{\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 ${ \mathcal P }$-space and the ${ \mathcal Q }$-space, see figure 1. The EOM for the orbitals of the ${ \mathcal P }$- and ${ \mathcal Q }$-space are obtained by projecting equation (17) on either an orbital of the ${ \mathcal P }$-space, $\langle {\phi }_{j}|$, or of the ${ \mathcal Q }$-space, $\langle {\phi }_{a}|$, as done in the following.

2.2.2.1. Equations of motion for the ${ \mathcal Q }$-space orbitals

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

$\begin{eqnarray}&&\displaystyle \sum _{I\in { \mathcal V }}{\rm{i}}{\dot{C}}_{I}\langle {{\rm{\Psi }}}_{i}^{a}| {{\rm{\Phi }}}_{I}\rangle +\langle {{\rm{\Psi }}}_{i}^{a}| ({\rm{i}}\hat{D}-\hat{H})| {\rm{\Psi }}\rangle =0,\end{eqnarray} \tag{ 18 }$

where we used the orthogonality between the orbitals of the ${ \mathcal P }$ and ${ \mathcal Q }$ spaces to get rid of the Lagrange multipliers. Moreover the inner product $\langle {{\rm{\Psi }}}_{i}^{a}| {{\rm{\Phi }}}_{I}\rangle =\langle {\rm{\Psi }}| {\hat{b}}_{i}^{\dagger }{\hat{b}}_{a}| {{\rm{\Phi }}}_{I}\rangle$ vanishes because in all configurations $| {{\rm{\Phi }}}_{I}\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

$\begin{eqnarray}&&\displaystyle \sum _{{pq}}({\rm{i}}{\eta }_{q}^{p}-{h}_{q}^{p})\langle {\rm{\Psi }}| {\hat{b}}_{i}^{\dagger }{\hat{b}}_{a}{\hat{b}}_{p}^{\dagger }{\hat{b}}_{q}| {\rm{\Psi }}\rangle =\displaystyle \frac{1}{2}\displaystyle \sum _{{pqrs}}{v}_{{qs}}^{{pr}}\langle {\rm{\Psi }}| {\hat{b}}_{i}^{\dagger }{\hat{b}}_{a}{\hat{b}}_{p}^{\dagger }{\hat{b}}_{r}^{\dagger }{\hat{b}}_{s}{\hat{b}}_{q}| {\rm{\Psi }}\rangle .\end{eqnarray} \tag{ 19 }$

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

$\begin{eqnarray}&&{\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},\end{eqnarray} \tag{ 20 }$

$\begin{eqnarray}&&{\hat{b}}_{i}^{\dagger }{\hat{b}}_{a}{\hat{b}}_{p}^{\dagger }{\hat{b}}_{r}^{\dagger }{\hat{b}}_{s}{\hat{b}}_{q}={\hat{b}}_{i}^{\dagger }{\hat{b}}_{r}^{\dagger }{\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}^{\dagger }{\hat{b}}_{s}{\hat{b}}_{q}{\hat{b}}_{a},\end{eqnarray} \tag{ 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 ${ \mathcal Q }$-space, which is not included in $| {\rm{\Psi }}\rangle$ and thus vanish. The lhs of equation (19) now reads,

$\begin{eqnarray}&&\displaystyle \sum _{q}({\rm{i}}{\eta }_{q}^{a}-{h}_{q}^{a})\langle {\rm{\Psi }}| {\hat{b}}_{i}^{\dagger }{\hat{b}}_{q}| {\rm{\Psi }}\rangle =\displaystyle \sum _{j}({\rm{i}}{\eta }_{j}^{a}-{h}_{j}^{a})\langle {\rm{\Psi }}| {\hat{b}}_{i}^{\dagger }{\hat{b}}_{j}| {\rm{\Psi }}\rangle ,\end{eqnarray} \tag{ 22 }$

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

$\begin{eqnarray}&&\displaystyle \sum _{j}({\rm{i}}{\eta }_{j}^{a}-{h}_{j}^{a})\langle {\rm{\Psi }}| {\hat{b}}_{i}^{\dagger }{\hat{b}}_{j}| {\rm{\Psi }}\rangle =\displaystyle \sum _{{jkl}}{v}_{{jl}}^{{ak}}\langle {\rm{\Psi }}| {\hat{b}}_{i}^{\dagger }{\hat{b}}_{k}^{\dagger }{\hat{b}}_{l}{\hat{b}}_{j}| {\rm{\Psi }}\rangle .\end{eqnarray} \tag{ 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 ${ \mathcal Q }$-space orbitals. This issue is circumvented by using the projector $\hat{Q}$ onto the subspace spanned by the ${ \mathcal Q }$-space orbitals,

$\begin{eqnarray}\hat{Q} & = & \displaystyle \sum _{a}| {\phi }_{a}\rangle \langle {\phi }_{a}| \\ & = & \hat{1}-\displaystyle \sum _{i}| {\phi }_{i}\rangle \langle {\phi }_{i}| \\ & = & \hat{1}-\hat{P},\end{eqnarray} \tag{ 24 }$

with $\hat{P}$ the projector onto the ${ \mathcal P }$-space. Introducing the one-body density matrix, ${\rho }_{i}^{j}=\langle {\rm{\Psi }}| {\hat{b}}_{i}^{\dagger }{\hat{b}}_{j}| {\rm{\Psi }}\rangle$ and the two-body density matrix ${\rho }_{{ik}}^{{jl}}=\langle {\rm{\Psi }}| {\hat{b}}_{i}^{\dagger }{\hat{b}}_{k}^{\dagger }{\hat{b}}_{l}{\hat{b}}_{j}| {\rm{\Psi }}\rangle$, we obtain,

$\begin{eqnarray}&&i\displaystyle \sum _{j}\hat{Q}| {\dot{\phi }}_{j}\rangle {\rho }_{i}^{j}=\hat{Q}\left[\displaystyle \sum _{j}\mathop{h}\limits^{}| {\phi }_{j}\rangle {\rho }_{i}^{j}+\displaystyle \sum _{{jkl}}{\mathop{W}\limits^{}}_{l}^{k}| {\phi }_{j}\rangle {\rho }_{{ik}}^{{jl}}\right],\end{eqnarray} \tag{ 25 }$

with,

$\begin{eqnarray}&&{\mathop{W}\limits^{}}_{l}^{k}({\bf{r}})=\int {\phi }_{k}^{* }({\bf{r}}^{\prime} )W({\bf{r}},{\bf{r}}^{\prime} ){\phi }_{l}({\bf{r}}^{\prime} ){\rm{d}}{\bf{r}}^{\prime} ,\end{eqnarray} \tag{ 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 ${ \mathcal Q }$-space appears in the time-derivative of the orbitals of the ${ \mathcal 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

$\begin{eqnarray}&&\hat{Q}| {\dot{\phi }}_{i}\rangle =-{\rm{i}}(\hat{1}-\hat{P})\left[\mathop{h}\limits^{}| {\phi }_{i}\rangle +\displaystyle \sum _{{jklm}}({\mathop{\mathop{{\boldsymbol{\rho }}}\limits_{\_}}\limits_{\_}}^{-1}{)}_{i}^{m}{\rho }_{{mk}}^{{jl}}{\mathop{W}\limits^{}}_{l}^{k}| {\phi }_{j}\rangle \right],\end{eqnarray} \tag{ 27 }$

with ${\underline{\underline{{\boldsymbol{\rho }}}}}^{-1}$ the inverse of the one-body density matrix. The MCTHB theory leads also to equation (27), see [70], but the lhs is subsequently simplified thanks to the choice of the matrix elements ${\eta }_{j}^{i}=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 }_{j}^{i}$ is not possible in the TD-RASSCF theory. The derivations of the ${ \mathcal 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 ${ \mathcal 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 ${ \mathcal P }$-space, $\langle {\phi }_{j}|$, leading to,

$\begin{eqnarray}&&\displaystyle \sum _{I\in { \mathcal V }}{\rm{i}}{\dot{C}}_{I}\langle {{\rm{\Psi }}}_{i}^{j}| {{\rm{\Phi }}}_{I}\rangle +\langle {{\rm{\Psi }}}_{i}^{j}| ({\rm{i}}\hat{D}-\hat{H})| {\rm{\Psi }}\rangle +{\epsilon }_{j}^{i}=0,\end{eqnarray} \tag{ 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 }_{i}|$, leads to an equation containing the same Lagrange multiplier,

$\begin{eqnarray}&&\displaystyle \sum _{I\in { \mathcal V }}(-{\rm{i}}{\dot{C}}_{I}^{* }\langle {{\rm{\Phi }}}_{I}| {{\rm{\Psi }}}_{j}^{i}\rangle )+\langle {\rm{\Psi }}| ({\rm{i}}\hat{D}-\hat{H})| {{\rm{\Psi }}}_{j}^{i}\rangle +{\epsilon }_{j}^{i}=0,\end{eqnarray} \tag{ 29 }$

and subtracting equations (28) and (29) gives the EOM for the ${ \mathcal P }$-space orbitals, i.e.,

$\begin{eqnarray}&&\langle {\rm{\Psi }}| ({\rm{i}}\hat{D}-\hat{H})| {{\rm{\Psi }}}_{j}^{i}\rangle -\langle {{\rm{\Psi }}}_{i}^{j}| ({\rm{i}}\hat{D}-\hat{H})| {\rm{\Psi }}\rangle ={\rm{i}}{\dot{\rho }}_{i}^{j},\end{eqnarray} \tag{ 30 }$

where we have introduced ${\dot{\rho }}_{i}^{j}\equiv {\sum }_{I\in { \mathcal V }}({\dot{C}}_{I}^{* }\langle {{\rm{\Phi }}}_{I}| {{\rm{\Psi }}}_{j}^{i}\rangle +\langle {{\rm{\Psi }}}_{i}^{j}| {{\rm{\Phi }}}_{I}\rangle {\dot{C}}_{I})$, with $| {{\rm{\Psi }}}_{j}^{i}\rangle \equiv {\hat{b}}_{i}^{\dagger }{\hat{b}}_{j}| {\rm{\Psi }}\rangle$. The ${ \mathcal P }$-space EOM provide the contribution of the ${ \mathcal P }$-space in the time-derivative of the orbitals, see equation (6),

$\begin{eqnarray}&&\hat{P}| {\dot{\phi }}_{i}\rangle =\displaystyle \sum _{j}| {\phi }_{j}\rangle {\eta }_{i}^{j},\end{eqnarray} \tag{ 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 ${ \mathcal 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}\ \mathrm{or}\ {{ \mathcal P }}_{2}$. The ${ \mathcal P }$-space equation (equation (30)) remains to be solved only for pairs of orbitals $\{i^{\prime} ,j^{\prime\prime} \}$, which belong to different ${{ \mathcal P }}_{i}$-spaces,

$\begin{eqnarray}&&\langle {\rm{\Psi }}| ({\rm{i}}\hat{D}-\hat{H})| {{\rm{\Psi }}}_{j^{\prime\prime} }^{i^{\prime} }\rangle -\langle {{\rm{\Psi }}}_{i^{\prime} }^{j^{\prime\prime} }| ({\rm{i}}\hat{D}-\hat{H})| {\rm{\Psi }}\rangle ={\rm{i}}{\mathop{\rho }\limits^{}}_{i^{\prime} }^{j^{\prime\prime} },\end{eqnarray} \tag{ 32 }$

but remains coupled to the amplitude equations through ${\dot{\rho }}_{i^{\prime} }^{j^{\prime\prime} }$. 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^{\prime} }^{j^{\prime\prime} }$ explicitly reads,

$\begin{eqnarray}&&{\dot{\rho }}_{i^{\prime} }^{j^{\prime\prime} }=\displaystyle \sum _{I\in { \mathcal V }}({\dot{C}}_{I}^{* }\langle {{\rm{\Phi }}}_{I}| {\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{j^{\prime\prime} }| {\rm{\Psi }}\rangle +\langle {\rm{\Psi }}| {\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{j^{\prime\prime} }| {{\rm{\Phi }}}_{I}\rangle {\dot{C}}_{I}).\end{eqnarray} \tag{ 33 }$

The action of ${\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{j^{\prime\prime} }$ on the wavefunction $| {\rm{\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^{\prime} }^{\dagger }{\hat{b}}_{j^{\prime\prime} }| {\rm{\Psi }}\rangle$ would contain only configurations with an odd number of particles in ${{ \mathcal P }}_{2}$, which makes the inner product with $\langle {{\rm{\Phi }}}_{I}| ,\forall I\in { \mathcal V }$ vanish. In the same manner ${\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{j^{\prime\prime} }$ acting on $| {{\rm{\Phi }}}_{I}\rangle$ is either zero, if $| {\phi }_{j^{\prime\prime} }\rangle$ is unoccupied in the configuration $| {{\rm{\Phi }}}_{I}\rangle$, or gives an odd number of particles in ${{ \mathcal P }}_{2}$. In this specific excitation scheme, ${\dot{\rho }}_{i^{\prime} }^{j^{\prime\prime} }=0$, for all pairs of orbitals $\{i^{\prime} ,j^{\prime\prime} \}$, leaving the amplitudes and the ${ \mathcal P }$-space orbitals equations uncoupled. Using the explicit expressions of the Hamiltonian (equation (8)) and the operator $\hat{D}$ (equation (14)), equation (32) reads,

$\begin{eqnarray}&&\displaystyle \sum _{{pq}}({h}_{q}^{p}-{\rm{i}}{\eta }_{q}^{p})[\langle {\rm{\Psi }}| {\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{j^{\prime\prime} }{\hat{b}}_{p}^{\dagger }{\hat{b}}_{q}-{\hat{b}}_{p}^{\dagger }{\hat{b}}_{q}{\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{j^{\prime\prime} }| {\rm{\Psi }}\rangle ]+\displaystyle \frac{1}{2}\displaystyle \sum _{{pqrs}}{v}_{{qs}}^{{pr}}[\langle {\rm{\Psi }}| {\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{j^{\prime\prime} }{\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^{\prime} }^{\dagger }{\hat{b}}_{j^{\prime\prime} }| {\rm{\Psi }}\rangle ]=0.\end{eqnarray} \tag{ 34 }$

We can simplify this expression, starting with

$\begin{eqnarray} & & \langle {\rm{\Psi }}| {\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{j^{\prime\prime} }{\hat{b}}_{p}^{\dagger }{\hat{b}}_{q}-{\hat{b}}_{p}^{\dagger }{\hat{b}}_{q}{\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{j^{\prime\prime} }| {\rm{\Psi }}\rangle \\ \quad & = & \langle {\rm{\Psi }}| {\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{q}{\delta }_{{pj}^{\prime\prime} }\pm {\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{p}^{\dagger }{\hat{b}}_{j^{\prime\prime} }{\hat{b}}_{q}-{\hat{b}}_{p}^{\dagger }{\hat{b}}_{j^{\prime\prime} }{\delta }_{i^{\prime} q}\mp {\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{p}^{\dagger }{\hat{b}}_{j^{\prime\prime} }{\hat{b}}_{q}| {\rm{\Psi }}\rangle \\ \quad & = & \langle {\rm{\Psi }}| {\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{q}{\delta }_{{pj}^{\prime\prime} }-{\hat{b}}_{p}^{\dagger }{\hat{b}}_{j^{\prime\prime} }{\delta }_{i^{\prime} q}| {\rm{\Psi }}\rangle \\ \quad & = & {\rho }_{i^{\prime} }^{q}{\delta }_{{pj}^{\prime\prime} }-{\rho }_{p}^{j^{\prime\prime} }{\delta }_{i^{\prime} q}\\ \quad & \equiv & {A}_{{pi}^{\prime} }^{{qj}^{\prime\prime} }.\end{eqnarray} \tag{ 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

$\begin{eqnarray}&&{\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{j^{\prime\prime} }{\hat{b}}_{p}^{\dagger }{\hat{b}}_{r}^{\dagger }{\hat{b}}_{s}{\hat{b}}_{q}={\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{r}^{\dagger }{\hat{b}}_{s}{\hat{b}}_{q}{\delta }_{{pj}^{\prime\prime} }\pm {\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{p}^{\dagger }{\hat{b}}_{s}{\hat{b}}_{q}{\delta }_{{rj}^{\prime\prime} }+{\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{p}^{\dagger }{\hat{b}}_{r}^{\dagger }{\hat{b}}_{j^{\prime\prime} }{\hat{b}}_{s}{\hat{b}}_{q}\end{eqnarray} \tag{ 36 }$

and the second product of operators as

$\begin{eqnarray}&&-{\hat{b}}_{p}^{\dagger }{\hat{b}}_{r}^{\dagger }{\hat{b}}_{s}{\hat{b}}_{q}{\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{j^{\prime\prime} }=-{\hat{b}}_{p}^{\dagger }{\hat{b}}_{r}^{\dagger }{\hat{b}}_{s}{\hat{b}}_{j^{\prime\prime} }{\delta }_{i^{\prime} q}\mp {\hat{b}}_{p}^{\dagger }{\hat{b}}_{r}^{\dagger }{\hat{b}}_{q}{\hat{b}}_{j^{\prime\prime} }{\delta }_{i^{\prime} s}-{\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{p}^{\dagger }{\hat{b}}_{r}^{\dagger }{\hat{b}}_{j^{\prime\prime} }{\hat{b}}_{s}{\hat{b}}_{q}.\end{eqnarray} \tag{ 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}_{{jl}}^{{ik}}={v}_{{lj}}^{{ki}}$ (equation (10)) and that ${\eta }_{q}^{p}$ must be evaluated for orbitals which belong to different ${{ \mathcal P }}_{i}$ space, $\{l^{\prime} ,k^{\prime\prime} \}$, equation (34) can be rewritten,

$\begin{eqnarray}&&\displaystyle \sum _{k^{\prime\prime} l^{\prime} }({h}_{l^{\prime} }^{k^{\prime\prime} }-{\rm{i}}{\eta }_{l^{\prime} }^{k^{\prime\prime} }){A}_{k^{\prime\prime} i^{\prime} }^{l^{\prime} j^{\prime\prime} }+\displaystyle \sum _{{klm}}({v}_{{kl}}^{j^{\prime\prime} m}{\rho }_{i^{\prime} m}^{{kl}}-{v}_{i^{\prime} m}^{{kl}}{\rho }_{{kl}}^{j^{\prime\prime} m})=0.\end{eqnarray} \tag{ 38 }$

This equation, used to determine ${\eta }_{l^{\prime} }^{k^{\prime\prime} }$, 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^{\prime} }^{k^{\prime\prime} }$, 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{\tilde{H}}={\rm{i}}\partial /\partial t-\hat{H}$. It means that the single excitations of the wavefunction, i.e., ${\hat{b}}_{i}^{\dagger }{\hat{b}}_{j}| {\rm{\Psi }}\rangle \equiv | {{\rm{\Psi }}}_{j}^{i}\rangle$, are not coupled to the wavefunction $| {\rm{\Psi }}\rangle$ through the action of $\hat{\tilde{H}}$, $\langle {\rm{\Psi }}| \hat{\tilde{H}}| {{\rm{\Psi }}}_{j}^{i}\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.

Zoom In Zoom Out Reset image size

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 }\leqslant 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,

$\begin{eqnarray}&&{ \mathcal V }={{ \mathcal V }}_{0}\oplus {{ \mathcal V }}_{1}\oplus \cdots \,\oplus \,{{ \mathcal V }}_{{N}_{\max }}.\end{eqnarray} \tag{ 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,

$\begin{eqnarray}&&{\rm{i}}{\dot{\rho }}_{i^{\prime} }^{j^{\prime\prime} }=\displaystyle \sum _{I\in { \mathcal V }}(\langle {\rm{\Psi }}| ({\rm{i}}\hat{D}-\hat{H})| {{\rm{\Phi }}}_{I}\rangle \langle {{\rm{\Phi }}}_{I}| {{\rm{\Psi }}}_{j^{\prime\prime} }^{i^{\prime} }\rangle -\langle {{\rm{\Psi }}}_{i^{\prime} }^{j^{\prime\prime} }| {{\rm{\Phi }}}_{I}\rangle \langle {{\rm{\Phi }}}_{I}| ({\rm{i}}\hat{D}-\hat{H})| {\rm{\Psi }}\rangle ).\end{eqnarray} \tag{ 40 }$

We introduce the projector onto the RAS space ${ \mathcal V }$ as $\hat{{\rm{\Pi }}}={\sum }_{I\in { \mathcal V }}| {{\rm{\Phi }}}_{I}\rangle \langle {{\rm{\Phi }}}_{I}|$. Using the above expression of ${\dot{\rho }}_{i}^{j}$, we obtain a new formulation of the ${ \mathcal P }$-space orbital equation,

$\begin{eqnarray}&&\langle {\rm{\Psi }}| ({\rm{i}}\hat{D}-\hat{H})(\hat{1}-\hat{{\rm{\Pi }}})| {{\rm{\Psi }}}_{j^{\prime\prime} }^{i^{\prime} }\rangle -\langle {{\rm{\Psi }}}_{i^{\prime} }^{j^{\prime\prime} }| (\hat{1}-\hat{{\rm{\Pi }}})({\rm{i}}\hat{D}-\hat{H})| {\rm{\Psi }}\rangle =0.\end{eqnarray} \tag{ 41 }$

For $| {\phi }_{i^{\prime} }\rangle \in {P}_{1}$ and $| {\phi }_{j^{\prime\prime} }\rangle \in {P}_{2}$, we note that $| {{\rm{\Psi }}}_{j^{\prime\prime} }^{i^{\prime} }\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{{\rm{\Pi }}})| {{\rm{\Psi }}}_{j^{\prime\prime} }^{i^{\prime} }\rangle =0$. On the other hand, $\langle {{\rm{\Psi }}}_{i^{\prime} }^{j^{\prime\prime} }| \text{}=\text{}\langle {\rm{\Psi }}| {\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{j^{\prime\prime} }$, 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 {{\rm{\Psi }}}_{i^{\prime} }^{j^{\prime\prime} }| (\hat{1}-\hat{{\rm{\Pi }}})\ne 0$ and equation (41) simplifies to,

$\begin{eqnarray}&&\langle {{\rm{\Psi }}}_{i^{\prime} }^{j^{\prime\prime} }| (\hat{1}-\hat{{\rm{\Pi }}})({\rm{i}}\hat{D}-\hat{H})| {\rm{\Psi }}\rangle =0,{\rm{with}}\langle {{\rm{\Psi }}}_{i^{\prime} }^{j^{\prime\prime} }| (\hat{1}-\hat{{\rm{\Pi }}})=\displaystyle \sum _{I\in {{ \mathcal V }}_{{N}_{\max }}}{C}_{I}^{* }\langle {{\rm{\Phi }}}_{I}| {\hat{b}}_{i^{\prime} }^{\dagger }{\hat{b}}_{j^{\prime\prime} }.\end{eqnarray} \tag{ 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^{\prime} ,k^{\prime\prime} \}$ which belong to different ${{ \mathcal P }}_{i}$ subspaces, equation (42) is equivalent to

$\begin{eqnarray}&&\displaystyle \sum _{k^{\prime\prime} l^{\prime} }({\rm{i}}{\eta }_{l^{\prime} }^{k^{\prime\prime} }-{h}_{l^{\prime} }^{k^{\prime\prime} }){\zeta }_{k^{\prime\prime} i^{\prime} }^{l^{\prime} j^{\prime\prime} }=\displaystyle \frac{1}{2}\displaystyle \sum _{{klmn}}{v}_{{ln}}^{{km}}{\zeta }_{{kmi}^{\prime} }^{{lnj}^{\prime\prime} },\end{eqnarray} \tag{ 43 }$

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

$\begin{eqnarray}&&{\zeta }_{k^{\prime\prime} i^{\prime} }^{l^{\prime} j^{\prime\prime} }=\langle {{\rm{\Psi }}}_{i^{\prime} }^{j^{\prime\prime} }| (\hat{1}-\hat{{\rm{\Pi }}}){\hat{b}}_{k^{\prime\prime} }^{\dagger }{\hat{b}}_{l^{\prime} }| {\rm{\Psi }}\rangle ,\end{eqnarray} \tag{ 44 }$

$\begin{eqnarray}&&{\zeta }_{{kmi}^{\prime} }^{{lnj}^{\prime\prime} }=\langle {{\rm{\Psi }}}_{i^{\prime} }^{j^{\prime\prime} }| (\hat{1}-\hat{{\rm{\Pi }}}){\hat{b}}_{k}^{\dagger }{\hat{b}}_{m}^{\dagger }{\hat{b}}_{n}{\hat{b}}_{l}| {\rm{\Psi }}\rangle .\end{eqnarray} \tag{ 45 }$

Here again, the ${ \mathcal P }$-space EOM for the determination of the ${\eta }_{l^{\prime} }^{k^{\prime\prime} }$, equation (43), are identical for bosons and fermions [91, 92]. These equations are solved to determine the ${\eta }_{l^{\prime} }^{k^{\prime\prime} }$ for each pairs of orbitals belonging to different ${{ \mathcal P }}_{i}$-spaces. The value of ${\eta }_{l^{\prime} }^{k^{\prime\prime} }$ 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 ${ \mathcal 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 }\gt 38$ with N = 50 or ${N}_{\max }\gt 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 MCTDH theory for both bosons and fermions was derived.

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,

$\begin{eqnarray}&&\hat{H}=\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{N}\left(-\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}^{2}}+{x}_{i}^{2}\right)+\lambda \displaystyle \sum _{i\lt j}\delta ({x}_{i}-{x}_{j}),\end{eqnarray} \tag{ 46 }$

using the unit of length ${l}_{0}=\sqrt{{\hslash }{(m{\omega }_{x})}^{-1}}$ and the unit of energy ${E}_{0}={\hslash }{\omega }_{x}$, with m the mass of the particles. Assuming no confinement induced resonances [107], the interaction strength, λ, is related to the 3D s-wave scattering length of the particles, a_s, through $\lambda =2{a}_{s}{l}_{0}{l}_{\perp }^{-2}$, with l_⊥ the transversal harmonic oscillator length. Experimentally, the 1D interaction strength can be tuned either by controlling the longitudinal and transversal frequencies or using an external magnetic field [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 \ne 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 ${{ \mathcal 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,

$\begin{eqnarray}&&{{ \mathcal E }}_{\mathrm{corr}}={E}_{\mathrm{GP}}-{E}_{\mathrm{method}},\end{eqnarray} \tag{ 47 }$

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

$\begin{eqnarray}&&{{ \mathcal E }}_{\mathrm{ref}}={E}_{\mathrm{GP}}-{E}_{\mathrm{MCTDHB}}^{5},\end{eqnarray} \tag{ 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 }}_{\mathrm{ref}}$. Secondly, we define the relative correlation energy, ${{ \mathcal E }}_{\mathrm{rel}}^{X}$, as the difference between the energy obtained from a MCTDHB calculation with X orbitals and the GP energy, i.e.,

$\begin{eqnarray}&&{{ \mathcal E }}_{{\rm{rel}}}^{X}={E}_{\mathrm{GP}}-{E}_{\mathrm{MCTDHB}}^{X}.\end{eqnarray} \tag{ 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 MCTDHB wavefunction. Thus, when the correlation energy of a RAS scheme with X orbitals is equal to ${{ \mathcal E }}_{\mathrm{rel}}^{X}$ the calculation is converged.

In table 1, we report the results for the GS energy (in units of E₀) obtained with $\lambda =0.01$. The reference for the correlation energy is ${{ \mathcal E }}_{\mathrm{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 }}_{\mathrm{ref}}$. Specifically we obtain 51%, 78% and 91% of ${{ \mathcal E }}_{\mathrm{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⁻³, 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 ${{ \mathcal P }}_{1}$ and ${{ \mathcal P }}_{2}$ spaces. In table 1 we report the results obtained for M = 2 to 8 orbitals with a single ${{ \mathcal P }}_{1}$ orbital, i.e., ${M}_{1}=1$ and ${M}_{2}=M-{M}_{1}$ ${{ \mathcal 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 ${{ \mathcal P }}_{1}$ to ${{ \mathcal P }}_{2}$. We follow this notation up to -SDTQ56789 and for larger excitations, we just indicate the value of ${N}_{\max }$ (e.g., ''-10'' means that all excitations from ${{ \mathcal P }}_{1}$ to ${{ \mathcal 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 }}_{\mathrm{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 }}_{\mathrm{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 }}_{\mathrm{rel}}^{X=3,4,5}$ correlation energies. Thus, using the -SD scheme with 5 orbitals $98.8 \%$ of ${{ \mathcal E }}_{\mathrm{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 -SDTQ56 schemes, we find that they include $99.99 \%$ and $99.9999 \%$ of ${{ \mathcal E }}_{\mathrm{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 }}_{\mathrm{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, whoes practical applicability is limited by the exponential growth of the number of configurations.

We also consider RAS schemes with only even excitations (see table 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 ${ \mathcal 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}_{\max }\geqslant 6$ and $M\geqslant 3$, the energy converges with respect to ${N}_{\max }$, as the energy does not change by increasing ${N}_{\max }$ further, but with energy slightly larger than the MCTDHB ones. The energy difference between the -DQ68 scheme and the MCTDHB calculation, with 5 orbitals for both methods, is $\sim 1.1\times {10}^{-5}$ and includes $99.96 \%$ of ${{ \mathcal E }}_{\mathrm{ref}}$. In the case of only even excitations, we also find that for $M\gt 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\gt 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 λ. As a limiting case, we note that in the Tonks-Girardeau model [20], obtained for an infinite value of λ, 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.

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,

$\begin{eqnarray}&&{\hat{H}}_{0}=\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{N}\left(-\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}^{2}}+{x}_{i}^{2}\right),\end{eqnarray} \tag{ 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,

$\begin{eqnarray}&&\hat{H}=\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{N}\left(-\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}^{2}}+{x}_{i}^{2}\right)+\displaystyle \sum _{i\lt j}\lambda {({x}_{i}-{x}_{j})}^{2},\end{eqnarray} \tag{ 51 }$

with λ the strength of the two-body interaction. This sudden change in the interaction between the bosons leads to a breathing dynamics of the BEC with frequencies ${{\rm{\Omega }}}_{n}=2n\sqrt{{\omega }_{x}^{2}+2N\lambda }$, with ${\omega }_{x}$ the frequency of the harmonic trap, see [25]. For positive values of λ the two-body interaction is attractive, while for negative values the interaction is repulsive and leads to unbound dynamics for $\lambda \lt -{\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\leqslant t\leqslant 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 }_{\mathrm{MCTDHB}}=3.46$. This value agrees perfectly with the analytical prediction ${{\rm{\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 ({{\rm{\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\lt 1$) as depicted in figure 3(a) and we obtain a lower frequency ${\omega }_{\mathrm{GP}}=3.35$. First, we perform TD-RASSCF-B simulations using a single ${{ \mathcal P }}_{1}$ orbital (${M}_{1}=1$) and ${M}_{2}=3$ ${{ \mathcal P }}_{2}$ orbitals for different RAS schemes reported in figure 3(a). For a short time, i.e., $0\leqslant t\leqslant 5$, all RAS schemes describe accurately the dynamics of the system, in contrast to the GP result. On the scale of the figure, the convergence to the 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\leqslant t\leqslant 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\gt 11.5$ and convergence is achieved for the -SDTQ56 excitation scheme with 84 configurations. We also investigate the role of the ${{ \mathcal 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 ${{ \mathcal P }}_{1}$ orbital, which differs from the MCTDHB results for $t\gt 11.5$ (not shown) while the -SD scheme with 2 ${{ \mathcal P }}_{1}$ orbitals differs for $t\gt 13.5$, see figure 3(b). We obtain converged results for the -SDT RAS scheme with 90 configurations.

Zoom In Zoom Out Reset image size

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\gt 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 ${{ \mathcal 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\gt 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\leqslant t\leqslant 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 }_{\mathrm{MCTDHB}}=6.63$, in perfect agreement with the analytical frequency ${{\rm{\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 ({{\rm{\Omega }}}_{1}t)$, is obtained. This is not surprising since the value of λ is now 5 times larger than the one of the previous example. Along with the MCTDHB result we report, for comparison, the result obtained using the mean-field GP theory, see figure 4(a), which strongly deviates from the MCTDHB result for $t\gt 0.3$ with a larger amplitude in the oscillations and gives a lower frequency for the oscillations, ${\omega }_{\mathrm{GP}}=6.32$. We start the TD-RASSCF-B simulations using M = 8 orbitals with one single ${{ \mathcal P }}_{1}$ orbital and ${M}_{2}=7$ ${{ \mathcal 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 ${{ \mathcal 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 ${{ \mathcal 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\leqslant t\leqslant 15$. For lower excitation schemes, the results are accurate for $0\leqslant t\leqslant 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.

Zoom In Zoom Out Reset image size

Using different numbers of ${{ \mathcal 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\gt 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 ${{ \mathcal P }}_{1}$-space and a factor ∼10 in comparison to the MCTDHB configurational space. The reduction of the configurational space impacts strongly the required CPU times of the simulations. For instance, the -SDTQ5678 RAS scheme with ${M}_{1}=1$ took ∼18.7 CPU hours, while the -SDTQ56 scheme with ${M}_{1}=2$, the -SDT scheme with ${M}_{1}=3$ and the -SD scheme with ${M}_{1}=4$ took ∼19.8, ∼6.9 and ∼11.1 CPU hours, respectively, on a 2.4 GHz Intel E5-2680 CPU. Using these RAS schemes, the CPU time is drastically reduced in comparison to the ∼113.4 CPU hours on a 2.5 GHz Intel E5-2680 CPU needed to perform the MCTDHB simulation, but the dynamics is accurately described. Moreover, the -SDTQ5 RAS scheme with ${M}_{1}=3$ converged to the exact solution with ∼4 times less CPU time. 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 ${{ \mathcal 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}_{1}=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 (-{{\rm{i}}E}_{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 asses 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).

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 ${ \mathcal 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.

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

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 $\{| {{\rm{\Phi }}}_{I}\rangle \}$. Moreover, additional illustrative examples are presented and discussed.

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}_{1},{n}_{2},\cdots ,{n}_{M}\rangle$, represents a configuration with n₁ bosons in the orbital $| {\phi }_{1}\rangle$, n₂ bosons in the orbital $| {\phi }_{2}\rangle$, etc. Such a configuration is indexed by an unique integer J defined as,

$\begin{eqnarray}&&J=1+\displaystyle \sum _{k=1}^{M}\left(\begin{array}{c}N+M-1-k-\displaystyle \sum _{l=1}^{k}{n}_{l}\\ M-k\end{array}\right).\end{eqnarray} \tag{ A1 }$

Thus, for each configuration we store its complex coefficient C_J in an array according to the index J provided by the above mapping. In [119], this mapping was employed to avoid the storage of the configuration vectors $| {n}_{1},{n}_{2},\cdots ,{n}_{M}\rangle$, 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_i 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 ${{ \mathcal P }}_{1}$ and ${{ \mathcal P }}_{2}$ spaces, respectively. For ${{ \mathcal V }}_{0}$, i.e., configurations with particles only in ${{ \mathcal P }}_{1}$, 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 ${{ \mathcal P }}_{2}$, we introduce the excitation number n_exc which is equivalent to the number of particles in ${{ \mathcal P }}_{2}$. The number of remaining particles in ${{ \mathcal P }}_{1}$ is $N-{n}_{{\rm{exc}}}$. For each excitation we enumerate the configurations of $N-{n}_{{\rm{exc}}}$ particles in the ${{ \mathcal P }}_{1}$ orbitals and compute their index. The same is performed for the configurations with n_exc in ${{ \mathcal P }}_{2}$ orbitals. Then the permanents for the total number of bosons are obtained by combining them,

$\begin{eqnarray}&&| {{\bf{n}}}^{{n}_{{\rm{exc}}}}\rangle =\left(\displaystyle \sum _{i=1}^{{N}_{1}^{{n}_{{\rm{exc}}}}}| {n}_{1},{n}_{2},\cdots ,{n}_{{M}_{1}}\rangle {}_{i}\right)\otimes \left(\displaystyle \sum _{i=1}^{{N}_{2}^{{n}_{{\rm{exc}}}}}| {n}_{{M}_{1}+1},\cdots ,{n}_{M}\rangle {}_{i}\right),\end{eqnarray} \tag{ A2 }$

where $| {{\bf{n}}}^{{n}_{{\rm{exc}}}}\rangle$ is an array of dimension $M\times ({N}_{1}^{{n}_{{\rm{exc}}}}\times {N}_{2}^{{n}_{{\rm{exc}}}})$, with ${N}_{1}^{{n}_{{\rm{exc}}}}$ (${N}_{2}^{{n}_{{\rm{exc}}}}$) the total number of configurations obtained from arranging the ${N}_{}-{n}_{{\rm{exc}}}$ (n_exc) particles in the ${{ \mathcal P }}_{1}$ (${{ \mathcal P }}_{2}$) orbitals. We evaluate the indexes of the configurations resulting from the ${{ \mathcal P }}_{1}$ subsystem, ${J}_{{P}_{1}}^{{n}_{{\rm{exc}}}}$, and from the ${{ \mathcal P }}_{2}$ subsystem, ${J}_{{P}_{2}}^{{n}_{{\rm{exc}}}}$, 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}_{{P}_{1}}^{{n}_{{\rm{exc}}}},{J}_{{P}_{2}}^{{n}_{{\rm{exc}}}},{n}_{{\rm{exc}}}\}$, which stores the position of the configuration in the configurational vector. This scheme is thus applied for all excitations, n_exc, included in the RAS scheme. With this storage or construction of the wavefunction, for a given configuration we have access to its coefficient in the following way. (i) We evaluate the index ${J}_{{P}_{1}}^{{n}_{{\rm{exc}}}}$ for the ${{ \mathcal P }}_{1}$ subsystem, (ii) we evaluate the index ${J}_{{P}_{2}}^{{n}_{{\rm{exc}}}}$ for the ${{ \mathcal P }}_{1}$ subsystem, (iii) we know or evaluate the excitation, n_exc, of the configuration (iv) we access to the index of the coefficient which is stored and the three component array at position $\{{J}_{{P}_{1}}^{{n}_{{\rm{exc}}}},{J}_{{P}_{2}}^{{n}_{{\rm{exc}}}},{n}_{{\rm{exc}}}\}$. This scheme has, as a draw back, the requirement to store the list of occupation numbers and the indexes to be efficient for numerical evaluation. But the idea behind the TD-RASSCF-B method is to reduce the size of the configurational space, which made such a storage manageable for the applications done so far.

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,

$\begin{eqnarray}&&{b}_{i}^{\dagger }{b}_{j}| {n}_{1},\cdots ,{n}_{j},\cdots ,{n}_{i},\cdots ,{n}_{M}\rangle =\sqrt{{n}_{j}}\sqrt{{n}_{i}+1}| {n}_{1},\cdots ,{n}_{j}-1,\cdots ,{n}_{i}+1,\cdots ,{n}_{M}\rangle .\end{eqnarray} \tag{ 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,

$\begin{eqnarray}&&{b}_{i}^{\dagger }{b}_{j}\left[\displaystyle \sum _{I\in { \mathcal V }}{C}_{I}| {{\rm{\Phi }}}_{I}\rangle \right]=\displaystyle \sum _{I\in { \mathcal V }}{C}_{I}\sqrt{{n}_{j}}\sqrt{{n}_{i}+1}| {{\rm{\Phi }}}_{I^{\prime} }\rangle ,\end{eqnarray} \tag{ B2 }$

with $| {{\rm{\Phi }}}_{I^{\prime} }\rangle$ the new configuration resulting from the action of ${b}_{i}^{\dagger }{b}_{j}$ on the initial configuration $| {{\rm{\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.,

$\begin{eqnarray}&&{b}_{i}^{\dagger }{b}_{j}\left[\displaystyle \sum _{I\in { \mathcal V }}{C}_{I}| {{\rm{\Phi }}}_{I}\rangle \right]=\displaystyle \sum _{I\in { \mathcal V }}{C^{\prime} }_{I}| {{\rm{\Phi }}}_{I}\rangle .\end{eqnarray} \tag{ B3 }$

In the basis of the configurational states, $\{| {{\rm{\Phi }}}_{I}\rangle \}$, the wavefunction is characterized by its coefficients only, and is stored as a vector. The new set of coefficients, $\{C{^{\prime} }_{I}\}$, resulting from the action of ${b}_{i}^{\dagger }{b}_{j}$, is obtained as,

$\begin{eqnarray}&&C{^{\prime} }_{I}=\langle {{\rm{\Phi }}}_{I}| {b}_{i}^{\dagger }{b}_{j}\left[\displaystyle \sum _{J\in { \mathcal V }}{C}_{J}| {{\rm{\Phi }}}_{J}\rangle \right].\end{eqnarray} \tag{ B4 }$

In practice, we apply ${b}_{i}^{\dagger }{b}_{j}$ on the bra $\langle {{\rm{\Phi }}}_{I}|$, which provides a new configurational state $\langle {{\rm{\Phi }}}_{J}|$. The configurational vectors are orthonormal, and only the configuration $| {{\rm{\Phi }}}_{J}\rangle$ in the sum remains from the projection. Thus, we evaluate the index of $\langle {{\rm{\Phi }}}_{J}|$ to directly access the coefficient C_J, i.e.,

$\begin{eqnarray}&&C{^{\prime} }_{I}={C}_{J}\sqrt{{n}_{i}}\sqrt{{n}_{j+1}}.\end{eqnarray} \tag{ 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 ${{ \mathcal P }}_{1}$ only $\{i^{\prime} ,j^{\prime} \}$, (ii) belong to ${{ \mathcal P }}_{2}$ only $\{i^{\prime\prime} ,j^{\prime\prime} \}$, or belong to ${{ \mathcal P }}_{1}$ and ${{ \mathcal P }}_{2}$, (iii) $\{i^{\prime} ,j^{\prime\prime} \}$ and (iv) $\{i^{\prime\prime} ,j^{\prime} \}$. For (i) and (ii) the excitation, or the number of particles in ${{ \mathcal P }}_{2}$, do not change in the resulting configuration. For the situation (iii) one particle is removed from ${{ \mathcal P }}_{2}$ and added in ${{ \mathcal P }}_{1}$ and the opposite happens for (iv). The excitation of the final configuration is thus known without counting the number of particles in ${{ \mathcal 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ⁱ (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}$.

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^{\prime} }^{j^{\prime\prime} }$. 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 $\underline{\underline{{\boldsymbol{\eta }}}}$, with elements ${\eta }_{i^{\prime} }^{j^{\prime\prime} }$, is anti-hermitian and thus ${\eta }_{j^{\prime\prime} }^{i^{\prime} }=-{({\eta }_{i^{\prime} }^{j^{\prime\prime} })}^{* }$, here $i^{\prime}$ and $j^{\prime\prime}$ correspond to orbitals of the ${{ \mathcal P }}_{1}$ and ${{ \mathcal P }}_{2}$ subspace, respectively. Writing down equations (38) or (43) for any set of orbitals $\{i^{\prime} ,j^{\prime\prime} \}$, 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^{\prime} ,j^{\prime\prime} \}$, this system can be written in a matrix form,

$\begin{eqnarray}&&\underline{\underline{{\bf{A}}}}\cdot \underline{{\bf{X}}}=\underline{{\bf{B}}},\end{eqnarray} \tag{ C1 }$

where the matrix $\underline{\underline{{\bf{A}}}}$, of dimension (${M}_{1}\times {M}_{2}$, ${M}_{1}\times {M}_{2}$), contains the values of ${A}_{k^{\prime\prime} i^{\prime} }^{l^{\prime} j^{\prime\prime} }$ or ${\zeta }_{k^{\prime\prime} i^{\prime} }^{l^{\prime} j^{\prime\prime} }$, for all set of $\{i^{\prime} ,j^{\prime\prime} \}$ and $\{l^{\prime} ,k^{\prime\prime} \}$. The vector $\underline{{\bf{B}}}$, of dimension (${M}_{1}\times {M}_{2}$), contains the rhs of equation (38) or (43) for each set of $\{i^{\prime} ,j^{\prime\prime} \}$ and the vector $\underline{{\bf{X}}}$ with the same dimension as $\underline{{\bf{B}}}$ contains the unknown values of ${\eta }_{i^{\prime} }^{j^{\prime\prime} }$ 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^{\prime} }^{j^{\prime\prime} }$ are trivially obtained from the elements $X\{i^{\prime} ,j^{\prime\prime} \}$ of the vector $\underline{{\bf{X}}}$,

$\begin{eqnarray}{\eta }_{i^{\prime} }^{j^{\prime\prime} }=\left\{\begin{array}{rl}{\rm{i}}(X\{i^{\prime} ,j^{\prime\prime} \}-{h}_{i^{\prime} }^{j^{\prime\prime} }) & \mathrm{for}\ \mathrm{equation}(38),\\ -{\rm{i}}(X\{i^{\prime} ,j^{\prime\prime} \}+{h}_{i^{\prime} }^{j^{\prime\prime} }) & \mathrm{for}\ \mathrm{equation}(43).\end{array}\right.\end{eqnarray} \tag{ C2 }$

After evaluating the matrix elements ${\eta }_{i^{\prime} }^{j^{\prime\prime} }$, the time derivative of the coefficients $\{{\dot{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,

$\begin{eqnarray}&&\hat{P}| \dot{{\phi }_{i}}\rangle =\displaystyle \sum _{j}^{M}| {\phi }_{j}\rangle {\eta }_{i}^{j}.\end{eqnarray} \tag{ C3 }$

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

$\begin{eqnarray}&&i\hat{Q}\underline{\underline{{\boldsymbol{\rho }}}}\underline{\dot{{\boldsymbol{X}}}}=\hat{Q}[\underline{\underline{{\boldsymbol{\rho }}}}\underline{\tilde{{\boldsymbol{h}}}}+\underline{\tilde{{\boldsymbol{W}}}}],\end{eqnarray} \tag{ C4 }$

with $\underline{\underline{{\boldsymbol{\rho }}}}$ the one-body reduced density matrix, $\underline{\dot{{\boldsymbol{X}}}}$ a vector collecting the time derivative of the orbitals, $\underline{\tilde{{\boldsymbol{h}}}}$ and $\underline{\tilde{{\boldsymbol{W}}}}$ are both vectors with elements $h({\boldsymbol{r}},t)| {\phi }_{i}\rangle$ and ${\sum }_{{jlk}}{\mathop{W}\limits^{}}_{l}^{k}| {\phi }_{j}\rangle {\rho }_{{ik}}^{{jl}}$, respectively. To obtain equation (C4), we used the fact that $\underline{\underline{{\boldsymbol{\rho }}}}$ commutes with the projector $\hat{Q}$, as easily seen from the equality $\hat{Q}=\hat{1}-\hat{P}$. Multiplying on the left by the inverse of the one-body density matrix, ${\underline{\underline{{\boldsymbol{\rho }}}}}^{-1}$, we have, for the orbital $| {\phi }_{i}\rangle$,

$\begin{eqnarray}\hat{Q}| {\dot{\phi }}_{i}\rangle & = & -{\rm{i}}\hat{Q}\left[\mathop{h}\limits^{}| {\phi }_{i}\rangle +\displaystyle \sum _{{jklm}}{({\mathop{\mathop{{\boldsymbol{\rho }}}\limits_{\_}}\limits_{\_}}^{-1})}_{i}^{m}{\rho }_{{mk}}^{{jl}}{\mathop{W}\limits^{}}_{l}^{k}| {\phi }_{j}\rangle \right]\\ & = & -{\rm{i}}(\hat{1}-\hat{P})\left[\mathop{h}\limits^{}| {\phi }_{i}\rangle +\displaystyle \sum _{{jklm}}{({\mathop{\mathop{{\boldsymbol{\rho }}}\limits_{\_}}\limits_{\_}}^{-1})}_{i}^{m}{\rho }_{{mk}}^{{jl}}{\mathop{W}\limits^{}}_{l}^{k}| {\phi }_{j}\rangle \right].\end{eqnarray} \tag{ 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{\underline{{\boldsymbol{A}}}}$, we follow the numerical implementation used for the MCTDH method [121] . The method uses a regularized one-body density matrix defined as [121],

$\begin{eqnarray}&&{\underline{\underline{{\boldsymbol{\rho }}}}}^{\mathrm{reg}}=\underline{\underline{{\boldsymbol{\rho }}}}+\epsilon {e}^{-\underline{\underline{{\boldsymbol{\rho }}}}/\epsilon },\end{eqnarray} \tag{ C6 }$

with a small real number (see appendix E.1). The same regularization is used for the matrix $\underline{\underline{{\boldsymbol{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{P}$ is numerically expressed by,

$\begin{eqnarray}&&\hat{P}=\displaystyle \sum _{ij}^{M}| {\phi }_{i}\rangle {({ \mathcal O })}_{{ij}}^{-1}\langle {\phi }_{j}| ,\end{eqnarray} \tag{ C7 }$

where ${{ \mathcal O }}_{{ij}}=\langle {\phi }_{i}| {\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.

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_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}}^{{jl}}$ (see text above equation (25)) are computed. These updates require ${M}^{4}{N}_{{\rm{grid}}}^{2}$ and ${M}^{4}{\dim ({ \mathcal V }}_{\mathrm{FCI}})$ operations, respectively, with ${{ \mathcal V }}_{\mathrm{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 }}_{\mathrm{FCI}})$ and ${M}^{4}{N}_{{\rm{grid}}}^{2}$ operations, respectively. The total cost is thus, approximately, $2{M}^{4}({N}_{{\rm{grid}}}^{2}+{\dim ({ \mathcal V }}_{\mathrm{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 ${ \mathcal 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}_{{\rm{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}{\rm{d}}{\rm{i}}{\rm{m}}({ \mathcal V })$, with ${ \mathcal V }$ the configurational space for the considered RAS scheme. In addition, we also need to solve the ${ \mathcal P }$-space equations, which requires on the order of M⁴ operations for excitation schemes with only even excitations (see equation (38)) and ${M}^{6}{\dim ({ \mathcal V }}_{{N}_{\max }})$ for the general RAS scheme (see equation (43)), with ${\dim ({ \mathcal V }}_{{N}_{\max }})$ the number of configurations including ${N}_{\max }$ particles in the ${{ \mathcal 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}_{\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,

$\begin{eqnarray}&&{\rm{\dim }}({ \mathcal V })=\left(\begin{array}{c}N+{M}_{1}-1\\ N\end{array}\right)+\displaystyle \sum _{k=1}^{{N}_{\max }/2}\left(\begin{array}{c}2k+{M}_{2}-1\\ 2k\end{array}\right)\left(\begin{array}{c}(N-2k)+{M}_{1}-1\\ N-2k\end{array}\right).\end{eqnarray} \tag{ 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 $2{M}^{4}({N}_{{\rm{grid}}}^{2}+\dim ({ \mathcal V })+{M}^{2}{\dim ({ \mathcal V }}_{{{N}}_{\max }}))$ and in the case of only even excitations it scales as $2{M}^{4}({N}_{{\rm{g}}{\rm{r}}{\rm{i}}{\rm{d}}}^{2}+\dim ({ \mathcal V })+1/2)$. To compare the numerical cost between the MCTDHB and TD-RASSCF-B methods, we can introduce ${\rm{\Delta }}({Op})$, the difference between the MCTDHB and TD-RASSCF-B operations to remove the constant number of operation resulting from N_grid,

$\begin{eqnarray}{\rm{\Delta }}(Op)=\phantom{\rule{}{2ex}}\{\begin{array}{cc}2{M}^{4}(\dim ({{ \mathcal V }}_{{\rm{F}}{\rm{C}}{\rm{I}}})-\dim ({ \mathcal V })-1/2) & {\rm{e}}{\rm{v}}{\rm{e}}{\rm{n}}\,{\rm{e}}{\rm{x}}{\rm{c}}{\rm{i}}{\rm{t}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}},\\ 2{M}^{4}(\dim ({{ \mathcal V }}_{{\rm{F}}{\rm{C}}{\rm{I}}})-\dim ({ \mathcal V })-{M}^{2}\dim ({{ \mathcal V }}_{{N}_{\max}})/2) & {\rm{g}}{\rm{e}}{\rm{n}}{\rm{e}}{\rm{r}}{\rm{a}}{\rm{l}}\,{\rm{s}}{\rm{c}}{\rm{h}}{\rm{e}}{\rm{m}}{\rm{e}}.\end{array}\end{eqnarray} \tag{ D2 }$

From the expression of ${\rm{\Delta }}({Op})$, a positive value represents a computational gain with the TD-RASSCF-B in comparison to the MCTDHB method, while a negative value is obtained when the MCTDHB method is more efficient. In the case of a scheme with only even excitations, ${\rm{\Delta }}({Op})$ is proportional to the size difference of the MCTDHB and TD-RASSCF-B configurational spaces and is always positive, which means that the TD-RASSCF-B method is always more efficient. In the case of the general excitation scheme the sixth-order tensor of the ${ \mathcal P }$-space equations, equation (45), can provide an overhead for the computation. To illustrate the computational efficiency we evaluate ${\rm{\Delta }}({Op})$ 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 ${{ \mathcal P }}_{1}$ orbital and $M-1$ orbitals in ${{ \mathcal 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 ${ \mathcal P }$-space equation. In the case of the general RAS scheme, there is always a value of ${N}_{\max }$ which leads to more operations in the TD-RASSCF-B than in the MCTDHB method, due to the evaluation of the six-order tensor in the ${ \mathcal P }$-space equation. But as shown in figure D1, this value is rather large, i.e. ${N}_{\max }=6$ for 10, ${N}_{\max }=40$ for 50 particles and ${N}_{\max }=90$ for 100 particles for the schemes depicted in figures D1(a), (b) and (c), respectively.

Zoom In Zoom Out Reset image size

E.1. Numerical parameters

E.2. More illustrative examples